PWG4/JetTasks/AliPWG4HighPtTrackQA.cxx

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