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