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