PWG4/JetTasks/AliPWG4HighPtTrackQA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------------
  17 // This class stores QA variables as function of pT for different type
  18 // of tracks and track selection criteria
  19 // Output: Histograms for different set of cuts
  20 //-----------------------------------------------------------------------
  21 // Author : Marta Verweij - UU
  22 //-----------------------------------------------------------------------
  23
  24 #ifndef ALIPWG4HIGHPTTRACKQA_CXX
  25 #define ALIPWG4HIGHPTTRACKQA_CXX
  26
  27 #include "AliPWG4HighPtTrackQA.h"
  28
  29 #include "TH1.h"
  30 #include "TH2.h"
  31 #include "TH3.h"
  32 #include "TProfile.h"
  33 #include "TList.h"
  34 #include "TFile.h"
  35 #include "TChain.h"
  36 #include "TH3F.h"
  37 #include "TKey.h"
  38 #include "TSystem.h"
  39 #include "TBits.h"
  40
  41 #include "AliAnalysisManager.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliStack.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliExternalTrackParam.h"
  49 #include "AliLog.h"
  50 #include "AliGenPythiaEventHeader.h"
  51 #include "AliGenCocktailEventHeader.h"
  52 #include "AliCentrality.h"
  53 #include "AliAODVertex.h"
  54 #include "AliAODEvent.h"
  55 //#include "AliAnalysisHelperJetTasks.h"
  56
  57 using namespace std; //required for resolving the 'cout' symbol
  58
  59 ClassImp(AliPWG4HighPtTrackQA)
  60
  61 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA()
  62 : AliAnalysisTaskSE(),
  63   fDataType(kESD),
  64   fEvent(0x0),
  65   fESD(0x0),
  66   fVtx(0x0),
  67   fTrackCuts(0x0),
  68   fTrackCutsITSLoose(0x0),
  69   fTrackCutsTPConly(0x0),
  70   fTrackType(0),
  71   fFilterMask(0),
  72   fSigmaConstrainedMax(-1.),
  73   fPtMax(100.),
  74   fIsPbPb(0),
  75   fCentClass(10),
  76   fNVariables(18),
  77   fVariables(0x0),
  78   fAvgTrials(1),
  79   fNEventAll(0),
  80   fNEventSel(0),
  81   fNEventReject(0),
  82   fh1Centrality(0x0),
  83   fh1Xsec(0),
  84   fh1Trials(0),
  85   fh1PtHard(0),
  86   fh1PtHardTrials(0),
  87   fh1NTracksAll(0x0),
  88   fh1NTracksReject(0x0),
  89   fh1NTracksSel(0x0),
  90   fPtAll(0),
  91   fPtSel(0),
  92   fPtPhi(0x0),
  93   fPtEta(0x0),
  94   fPtDCA2D(0x0),
  95   fPtDCAZ(0x0),
  96   fPtNClustersTPC(0x0),
  97   fPtNPointITS(0x0),
  98   fPtChi2C(0x0),
  99   fPtNSigmaToVertex(0x0),
 100   fPtRelUncertainty1Pt(0x0),
 101   fPtUncertainty1Pt(0x0),
 102   fPtChi2PerClusterTPC(0x0),
 103   fPtNCrossedRows(0x0),
 104   fPtNCrossedRowsNClusF(0x0),
 105   fPtNCrRNCrRNClusF(0x0),
 106   fPtSigmaY2(0x0),
 107   fPtSigmaZ2(0x0),
 108   fPtSigmaSnp2(0x0),
 109   fPtSigmaTgl2(0x0),
 110   fPtSigma1Pt2(0x0),
 111   fProfPtSigmaY2(0x0),
 112   fProfPtSigmaZ2(0x0),
 113   fProfPtSigmaSnp2(0x0),
 114   fProfPtSigmaTgl2(0x0),
 115   fProfPtSigma1Pt2(0x0),
 116   fProfPtSigma1Pt(0x0),
 117   fProfPtPtSigma1Pt(0x0),
 118   fSystTrackCuts(0x0),
 119   fHistList(0)
 120 {
 121   SetNVariables(18);
 122
 123   fPtBinEdges[0][0] = 10.;
 124   fPtBinEdges[0][1] = 1.;
 125   fPtBinEdges[1][0] = 20.;
 126   fPtBinEdges[1][1] = 2.;
 127   fPtBinEdges[2][0] = 50.;
 128   fPtBinEdges[2][1] = 5.;
 129
 130 }
 131 //________________________________________________________________________
 132 AliPWG4HighPtTrackQA::AliPWG4HighPtTrackQA(const char *name):
 133   AliAnalysisTaskSE(name),
 134   fDataType(kESD),
 135   fEvent(0x0),
 136   fESD(0x0),
 137   fVtx(0x0),
 138   fTrackCuts(0x0),
 139   fTrackCutsITSLoose(0x0),
 140   fTrackCutsTPConly(0x0),
 141   fTrackType(0),
 142   fFilterMask(0),
 143   fSigmaConstrainedMax(-1.),
 144   fPtMax(100.),
 145   fIsPbPb(0),
 146   fCentClass(10),
 147   fNVariables(18),
 148   fVariables(0x0),
 149   fAvgTrials(1),
 150   fNEventAll(0),
 151   fNEventSel(0),
 152   fNEventReject(0),
 153   fh1Centrality(0x0),
 154   fh1Xsec(0),
 155   fh1Trials(0),
 156   fh1PtHard(0),
 157   fh1PtHardTrials(0),
 158   fh1NTracksAll(0x0),
 159   fh1NTracksReject(0x0),
 160   fh1NTracksSel(0x0),
 161   fPtAll(0),
 162   fPtSel(0),
 163   fPtPhi(0x0),
 164   fPtEta(0x0),
 165   fPtDCA2D(0x0),
 166   fPtDCAZ(0x0),
 167   fPtNClustersTPC(0x0),
 168   fPtNPointITS(0x0),
 169   fPtChi2C(0x0),
 170   fPtNSigmaToVertex(0x0),
 171   fPtRelUncertainty1Pt(0x0),
 172   fPtUncertainty1Pt(0x0),
 173   fPtChi2PerClusterTPC(0x0),
 174   fPtNCrossedRows(0x0),
 175   fPtNCrossedRowsNClusF(0x0),
 176   fPtNCrRNCrRNClusF(0x0),
 177   fPtSigmaY2(0x0),
 178   fPtSigmaZ2(0x0),
 179   fPtSigmaSnp2(0x0),
 180   fPtSigmaTgl2(0x0),
 181   fPtSigma1Pt2(0x0),
 182   fProfPtSigmaY2(0x0),
 183   fProfPtSigmaZ2(0x0),
 184   fProfPtSigmaSnp2(0x0),
 185   fProfPtSigmaTgl2(0x0),
 186   fProfPtSigma1Pt2(0x0),
 187   fProfPtSigma1Pt(0x0),
 188   fProfPtPtSigma1Pt(0x0),
 189   fSystTrackCuts(0x0),
 190   fHistList(0)
 191 {
 192   //
 193   // Constructor. Initialization of Inputs and Outputs
 194   //
 195   AliDebug(2,Form("AliPWG4HighPtTrackQA Calling Constructor"));
 196
 197   SetNVariables(18);
 198
 199   fPtBinEdges[0][0] = 10.;
 200   fPtBinEdges[0][1] = 1.;
 201   fPtBinEdges[1][0] = 20.;
 202   fPtBinEdges[1][1] = 2.;
 203   fPtBinEdges[2][0] = 50.;
 204   fPtBinEdges[2][1] = 5.;
 205
 206   // Input slot #0 works with a TChain ESD
 207   DefineInput(0, TChain::Class());
 208   // Output slot #1 write into a TList
 209   DefineOutput(1, TList::Class());
 210 }
 211
 212 //________________________________________________________________________
 213 void AliPWG4HighPtTrackQA::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
 214
 215   if(region<3) {
 216     fPtBinEdges[region][0] = ptmax;
 217     fPtBinEdges[region][1] = ptBinWidth;
 218   }
 219   else {
 220     AliError("Only 3 regions alowed. Use region 0/1/2\n");
 221     return;
 222   }
 223
 224 }
 225
 226 //________________________________________________________________________
 227 void AliPWG4HighPtTrackQA::UserCreateOutputObjects() {
 228   //Create output objects
 229   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserCreateOutputObjects \n"));
 230
 231   Bool_t oldStatus = TH1::AddDirectoryStatus();
 232   TH1::AddDirectory(kFALSE);
 233
 234   OpenFile(1);
 235   fHistList = new TList();
 236   fHistList->SetOwner(kTRUE);
 237
 238   Float_t fgkPtMin = 0.;
 239   Float_t fgkPtMax = fPtMax;
 240
 241   //  Float_t ptBinEdges[2][2];
 242   // ptBinEdges[0][0] = 10.;
 243   // ptBinEdges[0][1] = 1.;
 244   // ptBinEdges[1][0] = 20.;
 245   // ptBinEdges[1][1] = 2.;
 246   // Float_t binWidth3 = 5.;
 247   // if(fPtMax>100.) {
 248   //   ptBinEdges[0][0] = 100.;
 249   //   ptBinEdges[0][1] = 5.;
 250   //   ptBinEdges[1][0] = 300.;
 251   //   ptBinEdges[1][1] = 10.;
 252   //   binWidth3 = 20.;
 253   // }
 254
 255   //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
 256   const Float_t ptmin1 =  fgkPtMin;
 257   const Float_t ptmax1 =  fPtBinEdges[0][0];
 258   const Float_t ptmin2 =  ptmax1 ;
 259   const Float_t ptmax2 =  fPtBinEdges[1][0];
 260   const Float_t ptmin3 =  ptmax2 ;
 261   const Float_t ptmax3 =  fgkPtMax;
 262   const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
 263   const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
 264   const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
 265   Int_t fgkNPtBins=nbin13;
 266   //Create array with low edges of each bin
 267   Double_t *binsPt=new Double_t[fgkNPtBins+1];
 268   for(Int_t i=0; i<=fgkNPtBins; i++) {
 269     if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
 270     if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 271     if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 272   }
 273
 274   Int_t fgkNPhiBins = 18*6;
 275   Float_t kMinPhi   = 0.;
 276   Float_t kMaxPhi   = 2.*TMath::Pi();
 277   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 278   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 279
 280   Int_t fgkNEtaBins=20;
 281   Float_t fgkEtaMin = -1.;
 282   Float_t fgkEtaMax = 1.;
 283   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 284   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 285
 286   Int_t fgkNNClustersTPCBins=80;
 287   Float_t fgkNClustersTPCMin = 0.5;
 288   Float_t fgkNClustersTPCMax = 160.5;
 289   Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
 290   for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
 291
 292   Int_t fgkNDCA2DBins=80;
 293   Float_t fgkDCA2DMin = -0.2;
 294   Float_t fgkDCA2DMax = 0.2;
 295   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
 296     fgkDCA2DMin = -2.;
 297     fgkDCA2DMax = 2.;
 298   }
 299   Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
 300   for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
 301
 302   Int_t fgkNDCAZBins=80;
 303   Float_t fgkDCAZMin = -2.;
 304   Float_t fgkDCAZMax = 2.;
 305   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
 306     fgkDCAZMin = -5.;
 307     fgkDCAZMax = 5.;
 308   }
 309   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
 310   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
 311
 312   Int_t fgkNNPointITSBins=9;
 313   Float_t fgkNPointITSMin = -0.5;
 314   Float_t fgkNPointITSMax = 8.5;
 315   Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
 316   for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
 317
 318   Int_t fgkNNSigmaToVertexBins=40;
 319   Float_t fgkNSigmaToVertexMin = 0.;
 320   Float_t fgkNSigmaToVertexMax = 8.;
 321   Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
 322   for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
 323
 324   Int_t fgkNChi2CBins=20;
 325   Float_t fgkChi2CMin = 0.;
 326   Float_t fgkChi2CMax = 100.;
 327   Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
 328   for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
 329
 330   Int_t fgkNRel1PtUncertaintyBins=30;
 331   Float_t fgkRel1PtUncertaintyMin = 0.;
 332   Float_t fgkRel1PtUncertaintyMax = 0.3;
 333   if(fTrackType==1 || fTrackType==2 || fTrackType==4)
 334     fgkRel1PtUncertaintyMax = 0.5;
 335   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
 336   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
 337
 338   Int_t fgkNUncertainty1PtBins = 30;
 339   Float_t fgkUncertainty1PtMin = 0.;
 340   Float_t fgkUncertainty1PtMax = 0.1;
 341   if(fTrackType==1 || fTrackType==2 || fTrackType==4)
 342     fgkUncertainty1PtMax = 0.2;
 343   Double_t *binsUncertainty1Pt=new Double_t[fgkNUncertainty1PtBins+1];
 344   for(Int_t i=0; i<=fgkNUncertainty1PtBins; i++) binsUncertainty1Pt[i]=(Double_t)fgkUncertainty1PtMin + (fgkUncertainty1PtMax-fgkUncertainty1PtMin)/fgkNUncertainty1PtBins*(Double_t)i ;
 345
 346   Float_t fgkChi2PerClusMin = 0.;
 347   Float_t fgkChi2PerClusMax = 4.;
 348   Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
 349   Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
 350   for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
 351
 352   Int_t fgkNCrossedRowsNClusFBins  = 50;
 353   Float_t fgkNCrossedRowsNClusFMin = 0.;
 354   Float_t fgkNCrossedRowsNClusFMax = 1.;
 355   Double_t *binsNCrossedRowsNClusF=new Double_t[fgkNCrossedRowsNClusFBins+1];
 356   for(Int_t i=0; i<=fgkNCrossedRowsNClusFBins; i++) binsNCrossedRowsNClusF[i]=(Double_t)fgkNCrossedRowsNClusFMin + (fgkNCrossedRowsNClusFMax-fgkNCrossedRowsNClusFMin)/fgkNCrossedRowsNClusFBins*(Double_t)i ;
 357
 358   Float_t fgk1PtMin = 0.;
 359   Float_t fgk1PtMax = 6.;
 360   Float_t binEdge1Pt1 = 1.;
 361   Float_t binWidth1Pt1 = 0.05;
 362   Int_t fgkN1PtBins1 = (int)((binEdge1Pt1-fgk1PtMin)/binWidth1Pt1);
 363   Float_t binWidth1Pt2 = 0.1;
 364   Int_t fgkN1PtBins2 = (int)((fgk1PtMax-binEdge1Pt1)/binWidth1Pt2);
 365   Int_t fgkN1PtBins = fgkN1PtBins1+fgkN1PtBins2;
 366   Double_t *bins1Pt=new Double_t[fgkN1PtBins+1];
 367
 368   for(Int_t i=0; i<=fgkN1PtBins; i++) {
 369     if(i<=fgkN1PtBins1)
 370       bins1Pt[i]=(Double_t)fgk1PtMin + (Double_t)(binEdge1Pt1-fgk1PtMin)/(Double_t)fgkN1PtBins1*(Double_t)i;
 371     if(i<=fgkN1PtBins && i>fgkN1PtBins1)
 372       bins1Pt[i]=(Double_t)binEdge1Pt1 + (Double_t)(fgk1PtMax-binEdge1Pt1)/(Double_t)fgkN1PtBins2*(Double_t)(i-fgkN1PtBins1);
 373   }
 374
 375   Int_t fgkNSigmaY2Bins = 50;
 376   Float_t fgkSigmaY2Min = 0.;
 377   Float_t fgkSigmaY2Max = 1.;
 378   if(fTrackType==1) fgkSigmaY2Max = 4.;
 379   if(fTrackType==2 || fTrackType==4) fgkSigmaY2Max = 0.1;
 380   Double_t *binsSigmaY2=new Double_t[fgkNSigmaY2Bins+1];
 381   for(Int_t i=0; i<=fgkNSigmaY2Bins; i++) binsSigmaY2[i]=(Double_t)fgkSigmaY2Min + (fgkSigmaY2Max-fgkSigmaY2Min)/fgkNSigmaY2Bins*(Double_t)i ;
 382
 383   Int_t fgkNSigmaZ2Bins = 50;
 384   Float_t fgkSigmaZ2Min = 0.;
 385   Float_t fgkSigmaZ2Max = 0.4;
 386   Double_t *binsSigmaZ2=new Double_t[fgkNSigmaZ2Bins+1];
 387   for(Int_t i=0; i<=fgkNSigmaZ2Bins; i++) binsSigmaZ2[i]=(Double_t)fgkSigmaZ2Min + (fgkSigmaZ2Max-fgkSigmaZ2Min)/fgkNSigmaZ2Bins*(Double_t)i ;
 388
 389   Int_t fgkNSigmaSnp2Bins = 50;
 390   Float_t fgkSigmaSnp2Min = 0.;
 391   Float_t fgkSigmaSnp2Max = 0.05;
 392   if(fTrackType==1) fgkSigmaSnp2Max = 0.2;
 393   if(fTrackType==2 || fTrackType==4) fgkSigmaSnp2Max = 0.1;
 394   Double_t *binsSigmaSnp2=new Double_t[fgkNSigmaSnp2Bins+1];
 395   for(Int_t i=0; i<=fgkNSigmaSnp2Bins; i++) binsSigmaSnp2[i]=(Double_t)fgkSigmaSnp2Min + (fgkSigmaSnp2Max-fgkSigmaSnp2Min)/fgkNSigmaSnp2Bins*(Double_t)i ;
 396
 397   Int_t fgkNSigmaTgl2Bins = 50;
 398   Float_t fgkSigmaTgl2Min = 0.;
 399   Float_t fgkSigmaTgl2Max = 0.1;
 400   if(fTrackType==1) fgkSigmaTgl2Max = 0.2;
 401   if(fTrackType==2 || fTrackType==4) fgkSigmaTgl2Max = 0.1;
 402   Double_t *binsSigmaTgl2=new Double_t[fgkNSigmaTgl2Bins+1];
 403   for(Int_t i=0; i<=fgkNSigmaTgl2Bins; i++) binsSigmaTgl2[i]=(Double_t)fgkSigmaTgl2Min + (fgkSigmaTgl2Max-fgkSigmaTgl2Min)/fgkNSigmaTgl2Bins*(Double_t)i ;
 404
 405   Int_t fgkNSigma1Pt2Bins = 50;
 406   Float_t fgkSigma1Pt2Min = 0.;
 407   Float_t fgkSigma1Pt2Max = 1.;
 408   Double_t *binsSigma1Pt2=new Double_t[fgkNSigma1Pt2Bins+1];
 409   for(Int_t i=0; i<=fgkNSigma1Pt2Bins; i++) binsSigma1Pt2[i]=(Double_t)fgkSigma1Pt2Min + (fgkSigma1Pt2Max-fgkSigma1Pt2Min)/fgkNSigma1Pt2Bins*(Double_t)i ;
 410
 411
 412   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 413   fHistList->Add(fNEventAll);
 414   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 415   fHistList->Add(fNEventSel);
 416   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 417   //Set labels
 418   fNEventReject->Fill("noESD",0);
 419   fNEventReject->Fill("Trigger",0);
 420   fNEventReject->Fill("NTracks<2",0);
 421   fNEventReject->Fill("noVTX",0);
 422   fNEventReject->Fill("VtxStatus",0);
 423   fNEventReject->Fill("NCont<2",0);
 424   fNEventReject->Fill("ZVTX>10",0);
 425   fNEventReject->Fill("cent",0);
 426   fNEventReject->Fill("cent>90",0);
 427   fHistList->Add(fNEventReject);
 428
 429   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
 430   fHistList->Add(fh1Centrality);
 431
 432   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 433   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 434   fHistList->Add(fh1Xsec);
 435
 436   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 437   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 438   fHistList->Add(fh1Trials);
 439
 440   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 441   fHistList->Add(fh1PtHard);
 442   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 443   fHistList->Add(fh1PtHardTrials);
 444
 445   fh1NTracksAll = new TH1F("fh1NTracksAll","fh1NTracksAll",1,-0.5,0.5);
 446   fHistList->Add(fh1NTracksAll);
 447
 448   fh1NTracksReject = new TH1F("fh1NTracksReject","fh1NTracksReject",1,-0.5,0.5);
 449   fh1NTracksReject->Fill("noESDtrack",0);
 450   fh1NTracksReject->Fill("noTPCInner",0);
 451   fh1NTracksReject->Fill("FillTPC",0);
 452   fh1NTracksReject->Fill("noTPConly",0);
 453   fh1NTracksReject->Fill("relate",0);
 454   fh1NTracksReject->Fill("trackCuts",0);
 455   fh1NTracksReject->Fill("laser",0);
 456   fh1NTracksReject->Fill("chi2",0);
 457   fHistList->Add(fh1NTracksReject);
 458
 459   fh1NTracksSel = new TH1F("fh1NTracksSel","fh1NTracksSel",1,-0.5,0.5);
 460   fHistList->Add(fh1NTracksSel);
 461
 462   fSystTrackCuts = new TH1F("fSystTrackCuts","fSystTrackCuts",1,-0.5,0.5);
 463   fSystTrackCuts->Fill("noCut",0);
 464   fSystTrackCuts->Fill("eta",0);
 465   fSystTrackCuts->Fill("0.15<pT<1e10",0);
 466   fSystTrackCuts->Fill("kink",0);
 467   fSystTrackCuts->Fill("NClusterTPC",0);
 468   fSystTrackCuts->Fill("Chi2PerNClusTPC",0);
 469   fSystTrackCuts->Fill("DCA2D",0);
 470   fHistList->Add(fSystTrackCuts);
 471
 472
 473   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
 474   fHistList->Add(fPtAll);
 475   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
 476   fHistList->Add(fPtSel);
 477
 478   fPtPhi = new TH2F("fPtPhi","fPtPhi",fgkNPtBins,binsPt,fgkNPhiBins,binsPhi);
 479   fHistList->Add(fPtPhi);
 480
 481   fPtEta = new TH2F("fPtEta","fPtEta",fgkNPtBins,binsPt,fgkNEtaBins,binsEta);
 482   fHistList->Add(fPtEta);
 483
 484   fPtDCA2D = new TH2F("fPtDCA2D","fPtDCA2D",fgkNPtBins,binsPt,fgkNDCA2DBins,binsDCA2D);
 485   fHistList->Add(fPtDCA2D);
 486
 487   fPtDCAZ = new TH2F("fPtDCAZ","fPtDCAZ",fgkNPtBins,binsPt,fgkNDCAZBins,binsDCAZ);
 488   fHistList->Add(fPtDCAZ);
 489
 490   fPtNClustersTPC = new TH2F("fPtNClustersTPC","fPtNClustersTPC",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 491   fHistList->Add(fPtNClustersTPC);
 492
 493   fPtNPointITS = new TH2F("fPtNPointITS","fPtNPointITS",fgkNPtBins,binsPt,fgkNNPointITSBins,binsNPointITS);
 494   fHistList->Add(fPtNPointITS);
 495
 496   fPtChi2C = new TH2F("fPtChi2C","fPtChi2C",fgkNPtBins,binsPt,fgkNChi2CBins,binsChi2C);
 497   fHistList->Add(fPtChi2C);
 498
 499   fPtNSigmaToVertex = new TH2F("fPtNSigmaToVertex","fPtNSigmaToVertex",fgkNPtBins,binsPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
 500   fHistList->Add(fPtNSigmaToVertex);
 501
 502   fPtRelUncertainty1Pt = new TH2F("fPtRelUncertainty1Pt","fPtRelUncertainty1Pt",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 503   fHistList->Add(fPtRelUncertainty1Pt);
 504
 505   fPtUncertainty1Pt = new TH2F("fPtUncertainty1Pt","fPtUncertainty1Pt",fgkNPtBins,binsPt,fgkNUncertainty1PtBins,binsUncertainty1Pt);
 506   fHistList->Add(fPtUncertainty1Pt);
 507
 508   fPtChi2PerClusterTPC = new TH2F("fPtChi2PerClusterTPC","fPtChi2PerClusterTPC",fgkNPtBins,binsPt,fgkNChi2PerClusBins,binsChi2PerClus);
 509   fHistList->Add(fPtChi2PerClusterTPC);
 510
 511   fPtNCrossedRows = new TH2F("fPtNCrossedRows","fPtNCrossedRows",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC);
 512   fHistList->Add(fPtNCrossedRows);
 513
 514   fPtNCrossedRowsNClusF = new TH2F("fPtNCrossedRowsNClusF","fPtNCrossedRowsNClusF",fgkNPtBins,binsPt,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
 515   fHistList->Add(fPtNCrossedRowsNClusF);
 516
 517   fPtNCrRNCrRNClusF = new TH3F("fPtNCrRNCrRNClusF","fPtNCrRNCrRNClusF",fgkNPtBins,binsPt,fgkNNClustersTPCBins,binsNClustersTPC,fgkNCrossedRowsNClusFBins,binsNCrossedRowsNClusF);
 518   fHistList->Add(fPtNCrRNCrRNClusF);
 519
 520   fPtSigmaY2 = new TH2F("fPtSigmaY2","fPtSigmaY2",fgkN1PtBins,bins1Pt,fgkNSigmaY2Bins,binsSigmaY2);
 521   fHistList->Add(fPtSigmaY2);
 522
 523   fPtSigmaZ2 = new TH2F("fPtSigmaZ2","fPtSigmaZ2",fgkN1PtBins,bins1Pt,fgkNSigmaZ2Bins,binsSigmaZ2);
 524   fHistList->Add(fPtSigmaZ2);
 525
 526   fPtSigmaSnp2 = new TH2F("fPtSigmaSnp2","fPtSigmaSnp2",fgkN1PtBins,bins1Pt,fgkNSigmaSnp2Bins,binsSigmaSnp2);
 527   fHistList->Add(fPtSigmaSnp2);
 528
 529   fPtSigmaTgl2 = new TH2F("fPtSigmaTgl2","fPtSigmaTgl2",fgkN1PtBins,bins1Pt,fgkNSigmaTgl2Bins,binsSigmaTgl2);
 530   fHistList->Add(fPtSigmaTgl2);
 531
 532   fPtSigma1Pt2 = new TH2F("fPtSigma1Pt2","fPtSigma1Pt2",fgkN1PtBins,bins1Pt,fgkNSigma1Pt2Bins,binsSigma1Pt2);
 533   fHistList->Add(fPtSigma1Pt2);
 534
 535   fProfPtSigmaY2 = new TProfile("fProfPtSigmaY2","fProfPtSigmaY2",fgkN1PtBins,bins1Pt);
 536   fHistList->Add(fProfPtSigmaY2);
 537
 538   fProfPtSigmaZ2 = new TProfile("fProfPtSigmaZ2","fProfPtSigmaZ2",fgkN1PtBins,bins1Pt);
 539   fHistList->Add(fProfPtSigmaZ2);
 540
 541   fProfPtSigmaSnp2 = new TProfile("fProfPtSigmaSnp2","fProfPtSigmaSnp2",fgkN1PtBins,bins1Pt);
 542   fHistList->Add(fProfPtSigmaSnp2);
 543
 544   fProfPtSigmaTgl2 = new TProfile("fProfPtSigmaTgl2","fProfPtSigmaTgl2",fgkN1PtBins,bins1Pt);
 545   fHistList->Add(fProfPtSigmaTgl2);
 546
 547   fProfPtSigma1Pt2 = new TProfile("fProfPtSigma1Pt2","fProfPtSigma1Pt2",fgkN1PtBins,bins1Pt);
 548   fHistList->Add(fProfPtSigma1Pt2);
 549
 550   fProfPtSigma1Pt = new TProfile("fProfPtSigma1Pt","fProfPtSigma1Pt;p_{T};#sigma(1/p_{T})",fgkNPtBins,binsPt);
 551   fHistList->Add(fProfPtSigma1Pt);
 552
 553   fProfPtPtSigma1Pt = new TProfile("fProfPtPtSigma1Pt","fProfPtPtSigma1Pt;p_{T};p_{T}#sigma(1/p_{T})",fgkNPtBins,binsPt);
 554   fHistList->Add(fProfPtPtSigma1Pt);
 555
 556   TH1::AddDirectory(oldStatus);
 557
 558   PostData(1, fHistList);
 559
 560   if(binsPhi)               delete [] binsPhi;
 561   if(binsPt)                delete [] binsPt;
 562   if(binsNClustersTPC)      delete [] binsNClustersTPC;
 563   if(binsDCA2D)             delete [] binsDCA2D;
 564   if(binsDCAZ)              delete [] binsDCAZ;
 565   if(binsNPointITS)         delete [] binsNPointITS;
 566   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
 567   if(binsChi2C)             delete [] binsChi2C;
 568   if(binsEta)               delete [] binsEta;
 569   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 570   if(binsUncertainty1Pt)    delete [] binsUncertainty1Pt;
 571   if(binsChi2PerClus)       delete [] binsChi2PerClus;
 572   if(binsChi2PerClus)       delete [] binsNCrossedRowsNClusF;
 573   if(bins1Pt)               delete [] bins1Pt;
 574   if(binsSigmaY2)           delete [] binsSigmaY2;
 575   if(binsSigmaZ2)           delete [] binsSigmaZ2;
 576   if(binsSigmaSnp2)         delete [] binsSigmaSnp2;
 577   if(binsSigmaTgl2)         delete [] binsSigmaTgl2;
 578   if(binsSigma1Pt2)         delete [] binsSigma1Pt2;
 579 }
 580
 581 //________________________________________________________________________
 582 Bool_t AliPWG4HighPtTrackQA::SelectEvent() {
 583   //
 584   // Decide if event should be selected for analysis
 585   //
 586
 587   // Checks following requirements:
 588   // - fEvent available
 589   // - trigger info from AliPhysicsSelection
 590   // - MCevent available
 591   // - number of reconstructed tracks > 1
 592   // - primary vertex reconstructed
 593   // - z-vertex < 10 cm
 594   // - centrality in case of PbPb
 595
 596   Bool_t selectEvent = kTRUE;
 597
 598   //fEvent object available?
 599   if (!fEvent) {
 600     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 601     fNEventReject->Fill("noAliVEvent",1);
 602     selectEvent = kFALSE;
 603     return selectEvent;
 604   }
 605
 606   //Check if number of reconstructed tracks is larger than 1
 607   if(!fEvent->GetNumberOfTracks() || fEvent->GetNumberOfTracks()<2)  {
 608     fNEventReject->Fill("NTracks<2",1);
 609     selectEvent = kFALSE;
 610     return selectEvent;
 611   }
 612
 613   //Check if vertex is reconstructed
 614   if(fDataType==kESD&&dynamic_cast<AliESDEvent*>(fEvent)) {
 615     fVtx = ((AliESDEvent*)fEvent)->GetPrimaryVertexSPD();
 616
 617     if(!fVtx) {
 618       fNEventReject->Fill("noVTX",1);
 619       selectEvent = kFALSE;
 620       return selectEvent;
 621     }
 622
 623     if(!fVtx->GetStatus()) {
 624       fNEventReject->Fill("VtxStatus",1);
 625       selectEvent = kFALSE;
 626       return selectEvent;
 627     }
 628
 629     // Need vertex cut
 630     if(fVtx->GetNContributors()<2) {
 631       fNEventReject->Fill("NCont<2",1);
 632       selectEvent = kFALSE;
 633       return selectEvent;
 634     }
 635
 636     //Check if z-vertex < 10 cm
 637     double primVtx[3];
 638     fVtx->GetXYZ(primVtx);
 639     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 640       fNEventReject->Fill("ZVTX>10",1);
 641       selectEvent = kFALSE;
 642       return selectEvent;
 643     }
 644   }
 645   else if(fDataType==kAOD&&dynamic_cast<AliAODEvent*>(fEvent)) {
 646     const AliAODVertex *vtx = ((AliAODEvent*)fEvent)->GetPrimaryVertexSPD();
 647     if(!vtx) {
 648       fNEventReject->Fill("noVTX",1);
 649       selectEvent = kFALSE;
 650       return selectEvent;
 651     }
 652
 653     // Need vertex cut
 654     if(vtx->GetNContributors()<2) {
 655       fNEventReject->Fill("NCont<2",1);
 656       selectEvent = kFALSE;
 657       return selectEvent;
 658     }
 659
 660     //Check if z-vertex < 10 cm
 661     double primVtx[3];
 662     vtx->GetXYZ(primVtx);
 663     if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 664       fNEventReject->Fill("ZVTX>10",1);
 665       selectEvent = kFALSE;
 666       return selectEvent;
 667     }
 668
 669   }
 670
 671   //Centrality selection should only be done in case of PbPb
 672   if(IsPbPb()) {
 673     Float_t cent = 0.;
 674     if(fCentClass!=CalculateCentrality(fEvent) && fCentClass!=10) {
 675       fNEventReject->Fill("cent",1);
 676       selectEvent = kFALSE;
 677       return selectEvent;
 678     }
 679     else {
 680       if(fDataType==kESD) {
 681         if(dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()) {
 682           cent = dynamic_cast<AliESDEvent*>(fEvent)->GetCentrality()->GetCentralityPercentile("V0M");
 683         }
 684       }
 685       else if(fDataType==kAOD) {
 686         if(dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality())
 687           cent = dynamic_cast<AliAODEvent*>(fEvent)->GetHeader()->GetCentrality();
 688        }
 689       if(cent>90.) {
 690         fNEventReject->Fill("cent>90",1);
 691         selectEvent = kFALSE;
 692         return selectEvent;
 693       }
 694       fh1Centrality->Fill(cent);
 695     }
 696   }
 697
 698   return selectEvent;
 699
 700 }
 701
 702 //________________________________________________________________________
 703 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliVEvent *ev){
 704   //
 705   // Get centrality from ESD or AOD
 706   //
 707
 708   if(fDataType==kESD)
 709     return CalculateCentrality(dynamic_cast<AliESDEvent*>(ev));
 710   else if(fDataType==kAOD)
 711     return CalculateCentrality(dynamic_cast<AliAODEvent*>(ev));
 712   else
 713     return 5;
 714 }
 715
 716 //________________________________________________________________________
 717 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliESDEvent *esd){
 718   //
 719   // Get centrality from ESD
 720   //
 721
 722   Float_t cent = -1;
 723
 724   if(esd){
 725     if(esd->GetCentrality()){
 726       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
 727       if(fDebug>3) printf("centrality: %f\n",cent);
 728     }
 729   }
 730
 731   return GetCentralityClass(cent);
 732
 733 }
 734
 735 //________________________________________________________________________
 736 Int_t AliPWG4HighPtTrackQA::CalculateCentrality(AliAODEvent *aod){
 737   //
 738   // Get centrality from AOD
 739   //
 740
 741   if(!aod) return 5;
 742   Float_t cent = aod->GetHeader()->GetCentrality();
 743   if(fDebug>3) printf("centrality: %f\n",cent);
 744
 745   return GetCentralityClass(cent);
 746
 747 }
 748
 749 //________________________________________________________________________
 750 Int_t AliPWG4HighPtTrackQA::GetCentralityClass(Float_t cent) {
 751   //
 752   // Get centrality class
 753   //
 754
 755   if(cent<0)  return 5; // OB - cent sometimes negative
 756   if(cent>80) return 4;
 757   if(cent>50) return 3;
 758   if(cent>30) return 2;
 759   if(cent>10) return 1;
 760   return 0;
 761
 762 }
 763
 764 //________________________________________________________________________
 765 void AliPWG4HighPtTrackQA::UserExec(Option_t *) {
 766   // Main loop
 767   // Called for each event
 768   AliDebug(2,Form(">> AliPWG4HighPtTrackQA::UserExec \n"));
 769
 770   fEvent = InputEvent();
 771   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 772
 773   // All events without selection
 774   fNEventAll->Fill(0.);
 775
 776   if(!SelectEvent()) {
 777     // Post output data
 778     PostData(1, fHistList);
 779     return;
 780   }
 781
 782
 783   //Need to keep track of selected events
 784   fNEventSel->Fill(0.);
 785
 786   fVariables = new TArrayF(fNVariables);
 787
 788   if(fDataType==kESD) DoAnalysisESD();
 789   if(fDataType==kAOD) DoAnalysisAOD();
 790
 791   //Delete old fVariables
 792   if(fVariables) delete fVariables;
 793
 794   // Post output data
 795   PostData(1, fHistList);
 796
 797 }
 798
 799 //________________________________________________________________________
 800 void AliPWG4HighPtTrackQA::DoAnalysisESD() {
 801
 802   if(!fESD) {
 803     PostData(1, fHistList);
 804     return;
 805   }
 806
 807   // ---- Get MC Header information (for MC productions in pThard bins) ----
 808   Double_t ptHard = 0.;
 809   Double_t nTrials = 1; // trials for MC trigger weight for real data
 810
 811   AliMCEventHandler *eventHandlerMC = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 812   if (eventHandlerMC) {
 813
 814     if(eventHandlerMC->MCEvent()){
 815       AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(eventHandlerMC->MCEvent());
 816       if(pythiaGenHeader){
 817         nTrials = pythiaGenHeader->Trials();
 818         ptHard  = pythiaGenHeader->GetPtHard();
 819
 820         fh1PtHard->Fill(ptHard);
 821         fh1PtHardTrials->Fill(ptHard,nTrials);
 822
 823         fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 824       }
 825     }
 826   }
 827
 828   Int_t nTracks = fESD->GetNumberOfTracks();
 829   AliDebug(2,Form("nTracks ESD%d", nTracks));
 830
 831   /*
 832     Variables to be put in fVariables
 833     0: pt
 834     1: phi
 835     2: eta
 836     3: dca2D
 837     4: dcaZ
 838     5: nClustersTPC
 839     6: nPointITS
 840     7: chi2C
 841     8: nSigmaToVertex
 842     9: relUncertainty1Pt
 843     10: chi2PerClusterTPC
 844     11: #crossed rows
 845     12: (#crossed rows)/(#findable clusters)
 846     13: SigmaY2
 847     14: SigmaZ2
 848     15: SigmaSnp2
 849     16: SigmaTgl2
 850     17: Sigma1Pt2
 851   */
 852
 853   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 854     fh1NTracksAll->Fill(0.);
 855
 856     //Get track for analysis
 857     AliESDtrack *track = 0x0;
 858     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 859     if(!esdtrack) {
 860       fh1NTracksReject->Fill("noESDtrack",1);
 861       continue;
 862     }
 863
 864     if(fTrackType==4) {
 865       FillSystematicCutHist(esdtrack);
 866       if (!(fTrackCuts->AcceptTrack(esdtrack))) {
 867         fh1NTracksReject->Fill("trackCuts",1);
 868         continue;
 869       }
 870     }
 871
 872     if(fTrackType==1)
 873       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 874     else if(fTrackType==2 || fTrackType==4) {
 875       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 876       if(!track) {
 877         fh1NTracksReject->Fill("noTPConly",1);
 878         continue;
 879       }
 880       AliExternalTrackParam exParam;
 881       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
 882       if( !relate ) {
 883         fh1NTracksReject->Fill("relate",1);
 884         delete track;
 885         continue;
 886       }
 887       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 888     }
 889     else if(fTrackType==5 || fTrackType==6) {
 890       if(fTrackCuts->AcceptTrack(esdtrack)) {
 891         continue;
 892       }
 893       else {
 894         if( !(fTrackCutsITSLoose->AcceptTrack(esdtrack)) && fTrackCutsTPConly->AcceptTrack(esdtrack) ) {
 895
 896           if(fTrackType==5) {
 897             //use TPConly constrained track
 898             track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 899             AliExternalTrackParam exParam;
 900             Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
 901             if( !relate ) {
 902               fh1NTracksReject->Fill("relate",1);
 903               delete track;
 904               continue;
 905             }
 906             track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 907           }
 908           else if(fTrackType==6) {
 909             //use global constrained track
 910             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 911
 912           }
 913         }
 914       }
 915     }
 916     else
 917       track = esdtrack;
 918
 919     if(!track) {
 920       //      if(fTrackType==1 || fTrackType==2 || fTrackType==4) delete track;
 921       continue;
 922     }
 923
 924     if(fTrackType==2 || fTrackType==4) {
 925       //Cut on chi2 of constrained fit
 926       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
 927         fh1NTracksReject->Fill("chi2",1);
 928         delete track;
 929         continue;
 930       }
 931     }
 932
 933     if(fTrackType!=4)
 934       FillSystematicCutHist(track);
 935
 936     fPtAll->Fill(track->Pt());
 937
 938     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
 939       fh1NTracksReject->Fill("trackCuts",1);
 940       if(fTrackType==1 || fTrackType==2) delete track;
 941       continue;
 942     }
 943
 944     fh1NTracksSel->Fill(0.);
 945
 946     fVariables->Reset(0.);
 947
 948     fVariables->SetAt(track->Pt(),0);
 949     fVariables->SetAt(track->Phi(),1);
 950     fVariables->SetAt(track->Eta(),2);
 951
 952     Float_t dca2D = 0.;
 953     Float_t dcaz  = 0.;
 954     if(fTrackType==0) {       //Global
 955       track->GetImpactParameters(dca2D,dcaz);
 956     }
 957     else if(fTrackType==1 || fTrackType==2 || fTrackType==4) {  //TPConly
 958       track->GetImpactParametersTPC(dca2D,dcaz);
 959     }
 960     fVariables->SetAt(dca2D,3);
 961     fVariables->SetAt(dcaz,5);
 962
 963     fVariables->SetAt((float)track->GetTPCNcls(),5);
 964
 965     Int_t nPointITS = 0;
 966     UChar_t itsMap = track->GetITSClusterMap();
 967     for (Int_t i=0; i < 6; i++) {
 968       if (itsMap & (1 << i))
 969         nPointITS ++;
 970     }
 971     fVariables->SetAt((float)nPointITS,6);
 972     Float_t chi2C = (float)track->GetConstrainedChi2();
 973     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
 974       chi2C = (float)track->GetConstrainedChi2TPC();
 975     fVariables->SetAt(chi2C,7);
 976     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
 977
 978     fVariables->SetAt(TMath::Sqrt(track->GetSigma1Pt2())*fVariables->At(0),9);
 979
 980     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
 981
 982     fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
 983     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
 984     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
 985     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
 986     fVariables->SetAt(track->GetSigmaY2(),13);
 987     fVariables->SetAt(track->GetSigmaZ2(),14);
 988     fVariables->SetAt(track->GetSigmaSnp2(),15);
 989     fVariables->SetAt(track->GetSigmaTgl2(),16);
 990     fVariables->SetAt(track->GetSigma1Pt2(),17);
 991
 992     FillHistograms();
 993
 994     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
 995
 996     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5) delete track;
 997
 998   }//track loop
 999
1000 }
1001
1002 //________________________________________________________________________
1003 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1004
1005   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1006   if(!aod)return;
1007   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1008
1009     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1010     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1011
1012     fVariables->Reset(0.);
1013
1014     fVariables->SetAt(aodtrack->Pt(),0);
1015     fVariables->SetAt(aodtrack->Phi(),1);
1016     fVariables->SetAt(aodtrack->Eta(),2);
1017
1018     Double_t dca[2] = {1e6,1e6};
1019     Double_t covar[3] = {1e6,1e6,1e6};
1020     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1021       fVariables->SetAt(dca[0],3);
1022       fVariables->SetAt(dca[1],4);
1023     }
1024
1025     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1026     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1027     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1028     fVariables->SetAt(0.,8);
1029     fVariables->SetAt(0.,9);
1030     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1031     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
1032     Float_t crossedRowsTPCNClsF = 0.;
1033     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1034     fVariables->SetAt(crossedRowsTPCNClsF,12);
1035
1036     fPtAll->Fill(fVariables->At(0));
1037
1038     FillHistograms();
1039
1040   }
1041
1042 }
1043
1044 //________________________________________________________________________
1045 void AliPWG4HighPtTrackQA::FillHistograms() {
1046
1047   fPtSel->Fill(fVariables->At(0));
1048   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1049   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1050   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1051   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1052   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1053   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1054   if(fDataType==kESD) {
1055     fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1056     fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1057     fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(9));
1058     fPtUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)/fVariables->At(9));
1059     fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1060     fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1061     fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1062     fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1063     fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1064
1065     fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1066     fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1067     fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1068     fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1069     fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1070     fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1071     fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1072   }
1073   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1074   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1075   fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
1076   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1077 }
1078
1079 //________________________________________________________________________
1080 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1081   //
1082   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1083   // This is to called in Notify and should provide the path to the AOD/ESD file
1084   // Copied from AliAnalysisTaskJetSpectrum2
1085   //
1086
1087   TString file(currFile);
1088   fXsec = 0;
1089   fTrials = 1;
1090
1091   if(file.Contains("root_archive.zip#")){
1092     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1093     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1094     file.Replace(pos+1,20,"");
1095   }
1096   else {
1097     // not an archive take the basename....
1098     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1099   }
1100   //  Printf("%s",file.Data());
1101
1102
1103   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
1104   if(!fxsec){
1105     // next trial fetch the histgram file
1106     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1107     if(!fxsec){
1108         // not a severe condition but inciate that we have no information
1109       return kFALSE;
1110     }
1111     else{
1112       // find the tlist we want to be independtent of the name so use the Tkey
1113       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1114       if(!key){
1115         fxsec->Close();
1116         return kFALSE;
1117       }
1118       TList *list = dynamic_cast<TList*>(key->ReadObj());
1119       if(!list){
1120         fxsec->Close();
1121         return kFALSE;
1122       }
1123       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1124       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1125       fxsec->Close();
1126     }
1127   } // no tree pyxsec.root
1128   else {
1129     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1130     if(!xtree){
1131       fxsec->Close();
1132       return kFALSE;
1133     }
1134     UInt_t   ntrials  = 0;
1135     Double_t  xsection  = 0;
1136     xtree->SetBranchAddress("xsection",&xsection);
1137     xtree->SetBranchAddress("ntrials",&ntrials);
1138     xtree->GetEntry(0);
1139     fTrials = ntrials;
1140     fXsec = xsection;
1141     fxsec->Close();
1142   }
1143   return kTRUE;
1144 }
1145 //________________________________________________________________________
1146 Bool_t AliPWG4HighPtTrackQA::Notify()
1147 {
1148   //
1149   // Implemented Notify() to read the cross sections
1150   // and number of trials from pyxsec.root
1151   // Copied from AliAnalysisTaskJetSpectrum2
1152   //
1153
1154   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1155   Float_t xsection = 0;
1156   Float_t ftrials  = 1;
1157
1158   fAvgTrials = 1;
1159   if(tree){
1160     TFile *curfile = tree->GetCurrentFile();
1161     if (!curfile) {
1162       Error("Notify","No current file");
1163       return kFALSE;
1164     }
1165     if(!fh1Xsec||!fh1Trials){
1166       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1167       return kFALSE;
1168     }
1169     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1170     fh1Xsec->Fill("<#sigma>",xsection);
1171     // construct a poor man average trials
1172     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1173     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1174   }
1175   return kTRUE;
1176 }
1177
1178 //________________________________________________________________________
1179 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1180
1181   if(!mcEvent)return 0;
1182   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1183   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1184   if(!pythiaGenHeader){
1185     // cocktail ??
1186     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1187
1188     if (!genCocktailHeader) {
1189       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1190       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1191       return 0;
1192     }
1193     TList* headerList = genCocktailHeader->GetHeaders();
1194     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1195       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1196       if (pythiaGenHeader)
1197         break;
1198     }
1199     if(!pythiaGenHeader){
1200       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1201       return 0;
1202     }
1203   }
1204   return pythiaGenHeader;
1205
1206 }
1207
1208 //_______________________________________________________________________
1209 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1210 {
1211   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1212
1213   //
1214   // TPC cluster information
1215   // type 0: get fraction of found/findable clusters with neighbourhood definition
1216   //      1: findable clusters with neighbourhood definition
1217   //      2: found clusters
1218   //
1219   // definition of findable clusters:
1220   //            a cluster is defined as findable if there is another cluster
1221   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1222   //           effects with a very simple algorithm.
1223   //
1224
1225   TBits fTPCClusterMap = tr->GetTPCClusterMap();
1226   if (type==2) return fTPCClusterMap.CountBits();
1227
1228   Int_t found=0;
1229   Int_t findable=0;
1230   Int_t last=-nNeighbours;
1231
1232   for (Int_t i=row0; i<row1; ++i){
1233     //look to current row
1234     if (fTPCClusterMap[i]) {
1235       last=i;
1236       ++found;
1237       ++findable;
1238       continue;
1239     }
1240     //look to nNeighbours before
1241     if ((i-last)<=nNeighbours) {
1242       ++findable;
1243       continue;
1244     }
1245     //look to nNeighbours after
1246     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1247       if (fTPCClusterMap[j]){
1248         ++findable;
1249         break;
1250       }
1251     }
1252   }
1253   if (type==1) return findable;
1254
1255   if (type==0){
1256     Float_t fraction=0;
1257     if (findable>0)
1258       fraction=(Float_t)found/(Float_t)findable;
1259     else
1260       fraction=0;
1261     return fraction;
1262   }
1263   return 0;  // undefined type - default value
1264 }
1265
1266 void AliPWG4HighPtTrackQA::FillSystematicCutHist(AliESDtrack *track) {
1267
1268   fSystTrackCuts->Fill("noCut",1);
1269   if(TMath::Abs(track->Eta())>0.9)
1270     fSystTrackCuts->Fill("eta",1);
1271   if(track->Pt()<0.15 || track->Pt()>1e10)
1272     fSystTrackCuts->Fill("0.15<pT<1e10",1);
1273   if(track->GetKinkIndex(0)>0)
1274     fSystTrackCuts->Fill("kink",1);
1275   if(track->GetTPCclusters(0)<70)
1276     fSystTrackCuts->Fill("NClusterTPC",1);
1277   if(track->GetTPCclusters(0)>0.) {
1278     if(track->GetTPCchi2()/Float_t(track->GetTPCclusters(0))>4.)
1279       fSystTrackCuts->Fill("Chi2PerNClusTPC",1);
1280   }
1281
1282   Float_t dcaToVertexXY = 0.;
1283   Float_t dcaToVertexZ  = 0.;
1284   track->GetImpactParameters(dcaToVertexXY,dcaToVertexZ);
1285   Float_t fCutMaxDCAToVertexXY = 2.4;
1286   Float_t fCutMaxDCAToVertexZ  = 3.2;
1287   Float_t dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1288   if(dcaToVertex>1.)
1289     fSystTrackCuts->Fill("DCA2D",1);
1290
1291 }
1292
1293 //________________________________________________________________________
1294 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1295 {
1296   // The Terminate() function is the last function to be called during
1297   // a query. It always runs on the client, it can be used to present
1298   // the results graphically or save the results to file.
1299
1300 }
1301
1302 #endif