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