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