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 if(fTrackType==7) {
 955       //use global constrained track
 956       track = esdtrack;
 957       //     track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 958     }
 959     else
 960       track = esdtrack;
 961
 962     if(!track) {
 963       continue;
 964     }
 965
 966     if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
 967       //Cut on chi2 of constrained fit
 968       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
 969         fh1NTracksReject->Fill("chi2",1);
 970         if(track) delete track;
 971         continue;
 972       }
 973     }
 974
 975     if(fTrackType!=4)
 976       FillSystematicCutHist(track);
 977
 978     fPtAll->Fill(track->Pt());
 979
 980     if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
 981       fh1NTracksReject->Fill("trackCuts",1);
 982       if(fTrackType==1 || fTrackType==2) {
 983         if(track) delete track;
 984       }
 985       continue;
 986     }
 987
 988     if(fTrackType==7) {
 989       if(fTrackCutsITSLoose ) {
 990         if(fTrackCutsITSLoose->AcceptTrack(track) )
 991           continue;
 992       }
 993
 994       if(esdtrack->GetConstrainedParam())
 995         track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 996     }
 997
 998     if(!track) continue;
 999
1000     fh1NTracksSel->Fill(0.);
1001
1002     fVariables->Reset(0.);
1003
1004     fVariables->SetAt(track->Pt(),0);
1005     fVariables->SetAt(track->Phi(),1);
1006     fVariables->SetAt(track->Eta(),2);
1007
1008     Float_t dca2D = 0.;
1009     Float_t dcaz  = 0.;
1010
1011     if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
1012       track->GetImpactParametersTPC(dca2D,dcaz); //TPConly
1013     }
1014     else
1015       track->GetImpactParameters(dca2D,dcaz);    //Global
1016
1017     fVariables->SetAt(dca2D,3);
1018     fVariables->SetAt(dcaz,5);
1019
1020     fVariables->SetAt((float)track->GetTPCNcls(),5);
1021
1022     Int_t nPointITS = 0;
1023     UChar_t itsMap = track->GetITSClusterMap();
1024     for (Int_t i=0; i < 6; i++) {
1025       if (itsMap & (1 << i))
1026         nPointITS ++;
1027     }
1028     fVariables->SetAt((float)nPointITS,6);
1029     Float_t chi2C = (float)track->GetConstrainedChi2();
1030     if(fTrackType==1 || fTrackType==2 || fTrackType==4)
1031       chi2C = (float)track->GetConstrainedChi2TPC();
1032     fVariables->SetAt(chi2C,7);
1033     fVariables->SetAt(fTrackCuts->GetSigmaToVertex(track),8);// Calculates the number of sigma to the vertex for a track.
1034
1035     fVariables->SetAt(GetTrackLengthTPC(track),9);
1036
1037     if(fVariables->At(5)>0.) fVariables->SetAt(track->GetTPCchi2()/fVariables->At(5),10);
1038
1039     fVariables->SetAt(track->GetTPCClusterInfo(2,1),11); //#crossed rows
1040     Float_t crossedRowsTPCNClsF = 1.;//track->GetTPCClusterInfo(2,0);
1041     if(track->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/track->GetTPCNclsF();
1042     fVariables->SetAt(crossedRowsTPCNClsF,12);//(#crossed rows)/(#findable clusters)
1043     fVariables->SetAt(track->GetSigmaY2(),13);
1044     fVariables->SetAt(track->GetSigmaZ2(),14);
1045     fVariables->SetAt(track->GetSigmaSnp2(),15);
1046     fVariables->SetAt(track->GetSigmaTgl2(),16);
1047     fVariables->SetAt(track->GetSigma1Pt2(),17);
1048
1049     fVariables->SetAt(track->GetTPCNclsIter1(),18);
1050     fVariables->SetAt(track->GetTPCchi2Iter1(),19);
1051
1052     FillHistograms();
1053
1054     //      int mult = fTrackCuts->CountAcceptedTracks(fESD);
1055
1056     if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5) {
1057       if(track) delete track;
1058     }
1059
1060   }//track loop
1061
1062 }
1063
1064 //________________________________________________________________________
1065 void AliPWG4HighPtTrackQA::DoAnalysisAOD() {
1066
1067   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(fEvent);
1068   if(!aod)return;
1069   AliExternalTrackParam *exParam = new  AliExternalTrackParam();
1070   for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
1071
1072     AliAODTrack *aodtrack = aod->GetTrack(iTrack);
1073     if( !aodtrack->TestFilterMask(fFilterMask) ) continue;
1074
1075     fVariables->Reset(0.);
1076
1077     fVariables->SetAt(aodtrack->Pt(),0);
1078     fVariables->SetAt(aodtrack->Phi(),1);
1079     fVariables->SetAt(aodtrack->Eta(),2);
1080
1081     Double_t dca[2] = {1e6,1e6};
1082     Double_t covar[3] = {1e6,1e6,1e6};
1083     if(aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),100.,dca,covar)) {
1084       fVariables->SetAt(dca[0],3);
1085       fVariables->SetAt(dca[1],4);
1086     }
1087
1088     fVariables->SetAt((float)aodtrack->GetTPCNcls(),5);
1089     fVariables->SetAt((float)aodtrack->GetITSNcls(),6);
1090     fVariables->SetAt(aodtrack->Chi2perNDF(),7);
1091     fVariables->SetAt(0.,8);
1092     fVariables->SetAt(GetTrackLengthTPC(aodtrack),9);
1093     fVariables->SetAt(aodtrack->Chi2perNDF(),10);
1094     fVariables->SetAt(GetTPCClusterInfo(aodtrack,2),11);
1095     Float_t crossedRowsTPCNClsF = 0.;
1096     if(aodtrack->GetTPCNclsF()>0.) crossedRowsTPCNClsF = fVariables->At(11)/aodtrack->GetTPCNclsF();
1097     fVariables->SetAt(crossedRowsTPCNClsF,12);
1098
1099     //get covariance matrix
1100     Double_t cov[21] = {0,};
1101     aodtrack->GetCovMatrix(cov);
1102     Double_t pxpypz[3] = {0,};
1103     aodtrack->PxPyPz(pxpypz);
1104     Double_t xyz[3] = {0,};
1105     aodtrack->GetXYZ(xyz);
1106     Short_t sign = aodtrack->Charge();
1107     exParam->Set(xyz,pxpypz,cov,sign);
1108
1109     fVariables->SetAt(exParam->GetSigmaY2(),13);
1110     fVariables->SetAt(exParam->GetSigmaZ2(),14);
1111     fVariables->SetAt(exParam->GetSigmaSnp2(),15);
1112     fVariables->SetAt(exParam->GetSigmaTgl2(),16);
1113     fVariables->SetAt(exParam->GetSigma1Pt2(),17);
1114
1115     fVariables->SetAt(0.,18);
1116     fVariables->SetAt(0.,19);
1117
1118     fPtAll->Fill(fVariables->At(0));
1119
1120     FillHistograms();
1121
1122   }
1123
1124 }
1125
1126 //________________________________________________________________________
1127 void AliPWG4HighPtTrackQA::FillHistograms() {
1128
1129   fPtSel->Fill(fVariables->At(0));
1130   fPtPhi->Fill(fVariables->At(0),fVariables->At(1));
1131   fPtEta->Fill(fVariables->At(0),fVariables->At(2));
1132   fPtDCA2D->Fill(fVariables->At(0),fVariables->At(3));
1133   fPtDCAZ->Fill(fVariables->At(0),fVariables->At(4));
1134   fPtNClustersTPC->Fill(fVariables->At(0),fVariables->At(5));
1135   fPtNPointITS->Fill(fVariables->At(0),fVariables->At(6));
1136
1137
1138   fPtNClustersTPCIter1->Fill(fVariables->At(0),fVariables->At(18));
1139   if(fVariables->At(18)>0.)
1140     fPtChi2PerClusterTPCIter1->Fill(fVariables->At(0),fVariables->At(19)/fVariables->At(18));
1141
1142   fPtChi2C->Fill(fVariables->At(0),fVariables->At(7));
1143   fPtNSigmaToVertex->Fill(fVariables->At(0),fVariables->At(8));
1144   fPtRelUncertainty1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1145   fPtRelUncertainty1PtNClus->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(5));
1146   fPtRelUncertainty1PtNClusIter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(18));
1147   fPtRelUncertainty1PtChi2->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(10));
1148   if(fVariables->At(18)>0.)
1149     fPtRelUncertainty1PtChi2Iter1->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(19)/fVariables->At(18));
1150   fPtRelUncertainty1PtPhi->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(1));
1151   fPtRelUncertainty1PtTrkLength->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)),fVariables->At(9));
1152
1153   fPtUncertainty1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1154   fPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1155   fPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1156   fPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1157   fPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1158   fPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1159
1160   fProfPtSigmaY2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(13)));
1161   fProfPtSigmaZ2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(14)));
1162   fProfPtSigmaSnp2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(15)));
1163   fProfPtSigmaTgl2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(16)));
1164   fProfPtSigma1Pt2->Fill(1./fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1165   fProfPtSigma1Pt->Fill(fVariables->At(0),TMath::Sqrt(fVariables->At(17)));
1166   fProfPtPtSigma1Pt->Fill(fVariables->At(0),fVariables->At(0)*TMath::Sqrt(fVariables->At(17)));
1167
1168   fPtChi2PerClusterTPC->Fill(fVariables->At(0),fVariables->At(10));
1169   fPtNCrossedRows->Fill(fVariables->At(0),fVariables->At(11));
1170   fPtNCrossedRowsNClusF->Fill(fVariables->At(0),fVariables->At(12));
1171   fPtNCrRNCrRNClusF->Fill(fVariables->At(0),fVariables->At(11),fVariables->At(12));
1172 }
1173
1174 //________________________________________________________________________
1175 Bool_t AliPWG4HighPtTrackQA::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
1176   //
1177   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
1178   // This is to called in Notify and should provide the path to the AOD/ESD file
1179   // Copied from AliAnalysisTaskJetSpectrum2
1180   //
1181
1182   TString file(currFile);
1183   fXsec = 0;
1184   fTrials = 1;
1185
1186   if(file.Contains("root_archive.zip#")){
1187     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
1188     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
1189     file.Replace(pos+1,20,"");
1190   }
1191   else {
1192     // not an archive take the basename....
1193     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
1194   }
1195   //  Printf("%s",file.Data());
1196
1197
1198   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
1199   if(!fxsec){
1200     // next trial fetch the histgram file
1201     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
1202     if(!fxsec){
1203         // not a severe condition but inciate that we have no information
1204       return kFALSE;
1205     }
1206     else{
1207       // find the tlist we want to be independtent of the name so use the Tkey
1208       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
1209       if(!key){
1210         fxsec->Close();
1211         return kFALSE;
1212       }
1213       TList *list = dynamic_cast<TList*>(key->ReadObj());
1214       if(!list){
1215         fxsec->Close();
1216         return kFALSE;
1217       }
1218       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1219       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1220       fxsec->Close();
1221     }
1222   } // no tree pyxsec.root
1223   else {
1224     TTree *xtree = (TTree*)fxsec->Get("Xsection");
1225     if(!xtree){
1226       fxsec->Close();
1227       return kFALSE;
1228     }
1229     UInt_t   ntrials  = 0;
1230     Double_t  xsection  = 0;
1231     xtree->SetBranchAddress("xsection",&xsection);
1232     xtree->SetBranchAddress("ntrials",&ntrials);
1233     xtree->GetEntry(0);
1234     fTrials = ntrials;
1235     fXsec = xsection;
1236     fxsec->Close();
1237   }
1238   return kTRUE;
1239 }
1240 //________________________________________________________________________
1241 Bool_t AliPWG4HighPtTrackQA::Notify()
1242 {
1243   //
1244   // Implemented Notify() to read the cross sections
1245   // and number of trials from pyxsec.root
1246   // Copied from AliAnalysisTaskJetSpectrum2
1247   //
1248
1249   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1250   Float_t xsection = 0;
1251   Float_t ftrials  = 1;
1252
1253   fAvgTrials = 1;
1254   if(tree){
1255     TFile *curfile = tree->GetCurrentFile();
1256     if (!curfile) {
1257       Error("Notify","No current file");
1258       return kFALSE;
1259     }
1260     if(!fh1Xsec||!fh1Trials){
1261       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1262       return kFALSE;
1263     }
1264     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
1265     fh1Xsec->Fill("<#sigma>",xsection);
1266     // construct a poor man average trials
1267     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
1268     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
1269   }
1270   return kTRUE;
1271 }
1272
1273 //________________________________________________________________________
1274 AliGenPythiaEventHeader*  AliPWG4HighPtTrackQA::GetPythiaEventHeader(AliMCEvent *mcEvent){
1275
1276   if(!mcEvent)return 0;
1277   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
1278   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
1279   if(!pythiaGenHeader){
1280     // cocktail ??
1281     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
1282
1283     if (!genCocktailHeader) {
1284       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
1285       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
1286       return 0;
1287     }
1288     TList* headerList = genCocktailHeader->GetHeaders();
1289     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1290       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1291       if (pythiaGenHeader)
1292         break;
1293     }
1294     if(!pythiaGenHeader){
1295       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
1296       return 0;
1297     }
1298   }
1299   return pythiaGenHeader;
1300
1301 }
1302
1303 //_______________________________________________________________________
1304 Float_t AliPWG4HighPtTrackQA::GetTPCClusterInfo(AliAODTrack *tr,Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1305 {
1306   //MV: copied from AliESDtrack since method is not available in AliAODTrack
1307
1308   //
1309   // TPC cluster information
1310   // type 0: get fraction of found/findable clusters with neighbourhood definition
1311   //      1: findable clusters with neighbourhood definition
1312   //      2: found clusters
1313   //
1314   // definition of findable clusters:
1315   //            a cluster is defined as findable if there is another cluster
1316   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1317   //           effects with a very simple algorithm.
1318   //
1319
1320   TBits fTPCClusterMap = tr->GetTPCClusterMap();
1321   if (type==2) return fTPCClusterMap.CountBits();
1322
1323   Int_t found=0;
1324   Int_t findable=0;
1325   Int_t last=-nNeighbours;
1326
1327   for (Int_t i=row0; i<row1; ++i){
1328     //look to current row
1329     if (fTPCClusterMap[i]) {
1330       last=i;
1331       ++found;
1332       ++findable;
1333       continue;
1334     }
1335     //look to nNeighbours before
1336     if ((i-last)<=nNeighbours) {
1337       ++findable;
1338       continue;
1339     }
1340     //look to nNeighbours after
1341     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1342       if (fTPCClusterMap[j]){
1343         ++findable;
1344         break;
1345       }
1346     }
1347   }
1348   if (type==1) return findable;
1349
1350   if (type==0){
1351     Float_t fraction=0;
1352     if (findable>0)
1353       fraction=(Float_t)found/(Float_t)findable;
1354     else
1355       fraction=0;
1356     return fraction;
1357   }
1358   return 0;  // undefined type - default value
1359 }
1360
1361 //_______________________________________________________________________
1362 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(AliESDtrack *track) {
1363   //
1364   // returns distance between 1st and last hit in TPC
1365   // distance given in number of padrows
1366   //
1367
1368   TBits fTPCClusterMap = track->GetTPCClusterMap();
1369   int firstHit = 0;
1370   int lastHit = 0;
1371
1372   for(int i=0; i<=159; i++) {
1373     if(fTPCClusterMap[i]>0) firstHit = i;
1374   }
1375   for(int i=159; i>=0; i--) {
1376     if(fTPCClusterMap[i]>0) lastHit = i;
1377   }
1378
1379   Int_t trackLength = lastHit - firstHit;
1380
1381   return trackLength;
1382 }
1383
1384 //_______________________________________________________________________
1385 Int_t AliPWG4HighPtTrackQA::GetTrackLengthTPC(AliAODTrack *track) {
1386   //
1387   // returns distance between 1st and last hit in TPC
1388   // distance given in number of padrows
1389   //
1390
1391   TBits fTPCClusterMap = track->GetTPCClusterMap();
1392   int firstHit = 0;
1393   int lastHit = 0;
1394
1395   for(int i=0; i<=159; i++) {
1396     if(fTPCClusterMap[i]>0) firstHit = i;
1397   }
1398   for(int i=159; i>=0; i--) {
1399     if(fTPCClusterMap[i]>0) lastHit = i;
1400   }
1401
1402   Int_t trackLength = lastHit - firstHit;
1403
1404   return trackLength;
1405 }
1406
1407 //_______________________________________________________________________
1408 void AliPWG4HighPtTrackQA::FillSystematicCutHist(AliESDtrack *track) {
1409
1410   fSystTrackCuts->Fill("noCut",1);
1411   if(TMath::Abs(track->Eta())>0.9)
1412     fSystTrackCuts->Fill("eta",1);
1413   if(track->Pt()<0.15 || track->Pt()>1e10)
1414     fSystTrackCuts->Fill("0.15<pT<1e10",1);
1415   if(track->GetKinkIndex(0)>0)
1416     fSystTrackCuts->Fill("kink",1);
1417   if(track->GetTPCclusters(0)<70)
1418     fSystTrackCuts->Fill("NClusterTPC",1);
1419   if(track->GetTPCclusters(0)>0.) {
1420     if(track->GetTPCchi2()/Float_t(track->GetTPCclusters(0))>4.)
1421       fSystTrackCuts->Fill("Chi2PerNClusTPC",1);
1422   }
1423
1424   Float_t dcaToVertexXY = 0.;
1425   Float_t dcaToVertexZ  = 0.;
1426   track->GetImpactParameters(dcaToVertexXY,dcaToVertexZ);
1427   Float_t fCutMaxDCAToVertexXY = 2.4;
1428   Float_t fCutMaxDCAToVertexZ  = 3.2;
1429   Float_t dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
1430   if(dcaToVertex>1.)
1431     fSystTrackCuts->Fill("DCA2D",1);
1432
1433 }
1434
1435 //________________________________________________________________________
1436 void AliPWG4HighPtTrackQA::Terminate(Option_t *)
1437 {
1438   // The Terminate() function is the last function to be called during
1439   // a query. It always runs on the client, it can be used to present
1440   // the results graphically or save the results to file.
1441
1442 }
1443
1444 #endif