PWG4/JetTasks/AliPWG4HighPtQAMC.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 compares the global reconstruction with the MC information
  18 // Momentum resolution is stored as function of track cuts and pt.
  19 // Output: Histograms for different set of cuts
  20 //-----------------------------------------------------------------------
  21 // Author : Marta Verweij - UU
  22 //-----------------------------------------------------------------------
  23
  24 #ifndef ALIPWG4HighPtQAMC_CXX
  25 #define ALIPWG4HighPtQAMC_CXX
  26
  27 #include "AliPWG4HighPtQAMC.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
  40 #include "AliAnalysisTask.h"
  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 "AliAnalysisHelperJetTasks.h"
  53
  54 using namespace std; //required for resolving the 'cout' symbol
  55
  56 ClassImp(AliPWG4HighPtQAMC)
  57
  58 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC()
  59 : AliAnalysisTask("AliPWG4HighPtQAMC", ""),
  60   fESD(0),
  61   fMC(0),
  62   fStack(0),
  63   fVtx(0x0),
  64   fTrackCuts(0),
  65   fTrackCutsReject(),
  66   fTrackType(0),
  67   fSigmaConstrainedMax(1e6),
  68   fPtMax(100.),
  69   fAvgTrials(1),
  70   fNEventAll(0),
  71   fNEventSel(0),
  72   fNEventReject(0),
  73   fh1Xsec(0),
  74   fh1Trials(0),
  75   fh1PtHard(0),
  76   fh1PtHardTrials(0),
  77   fPtAll(0),
  78   fPtSel(0),
  79   fPtSelFakes(0),
  80   fNPointTPCFakes(0),
  81   fPtSelLargeLabel(0),
  82   fMultRec(0),
  83   fNPointTPCMultRec(0),
  84   fDeltaPtMultRec(0),
  85   fPtAllvsPtMC(0),
  86   fPtAllminPtMCvsPtAll(0),
  87   fPtAllvsPtMCvsMult(0),
  88   fPtAllminPtMCvsPtAllvsMult(0),
  89   fPtAllminPtMCvsPtAllNPointTPC(0),
  90   fPtAllminPtMCvsPtAllNPointTPCIter1(0),
  91   fPtAllminPtMCvsPtAllChi2TPC(0),
  92   fPtAllminPtMCvsPtAllChi2TPCIter1(0),
  93   fPtAllminPtMCvsPtAllDCAR(0),
  94   fPtAllminPtMCvsPtAllDCAZ(0),
  95   fPtAllminPtMCvsPtAllPhi(0),
  96   fPtAllminPtMCvsPtAllNPointITS(0),
  97   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
  98   fPtAllminPtMCvsPtAllChi2C(0),
  99   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
 100   fPtAllMC(0),
 101   fPtSelMC(0),
 102   fHistList(0)
 103 {
 104
 105   fPtBinEdges[0][0] = 10.;
 106   fPtBinEdges[0][1] = 1.;
 107   fPtBinEdges[1][0] = 20.;
 108   fPtBinEdges[1][1] = 2.;
 109   fPtBinEdges[2][0] = 100.;
 110   fPtBinEdges[2][1] = 10.;
 111
 112 }
 113 //________________________________________________________________________
 114 AliPWG4HighPtQAMC::AliPWG4HighPtQAMC(const char *name):
 115   AliAnalysisTask(name,""),
 116   fESD(0),
 117   fMC(0),
 118   fStack(0),
 119   fVtx(0x0),
 120   fTrackCuts(),
 121   fTrackCutsReject(),
 122   fTrackType(0),
 123   fSigmaConstrainedMax(1e6),
 124   fPtMax(100.),
 125   fAvgTrials(1),
 126   fNEventAll(0),
 127   fNEventSel(0),
 128   fNEventReject(0),
 129   fh1Xsec(0),
 130   fh1Trials(0),
 131   fh1PtHard(0),
 132   fh1PtHardTrials(0),
 133   fPtAll(0),
 134   fPtSel(0),
 135   fPtSelFakes(0),
 136   fNPointTPCFakes(0),
 137   fPtSelLargeLabel(0),
 138   fMultRec(0),
 139   fNPointTPCMultRec(0),
 140   fDeltaPtMultRec(0),
 141   fPtAllvsPtMC(0),
 142   fPtAllminPtMCvsPtAll(0),
 143   fPtAllvsPtMCvsMult(0),
 144   fPtAllminPtMCvsPtAllvsMult(0),
 145   fPtAllminPtMCvsPtAllNPointTPC(0),
 146   fPtAllminPtMCvsPtAllNPointTPCIter1(0),
 147   fPtAllminPtMCvsPtAllChi2TPC(0),
 148   fPtAllminPtMCvsPtAllChi2TPCIter1(0),
 149   fPtAllminPtMCvsPtAllDCAR(0),
 150   fPtAllminPtMCvsPtAllDCAZ(0),
 151   fPtAllminPtMCvsPtAllPhi(0),
 152   fPtAllminPtMCvsPtAllNPointITS(0),
 153   fPtAllminPtMCvsPtAllNSigmaToVertex(0),
 154   fPtAllminPtMCvsPtAllChi2C(0),
 155   fPtAllminPtMCvsPtAllRel1PtUncertainty(0),
 156   fPtAllMC(0),
 157   fPtSelMC(0),
 158   fHistList(0)
 159 {
 160   //
 161   // Constructor. Initialization of Inputs and Outputs
 162   //
 163   AliDebug(2,Form("AliPWG4HighPtQAMC Calling Constructor"));
 164
 165   fPtBinEdges[0][0] = 10.;
 166   fPtBinEdges[0][1] = 1.;
 167   fPtBinEdges[1][0] = 20.;
 168   fPtBinEdges[1][1] = 2.;
 169   fPtBinEdges[2][0] = 100.;
 170   fPtBinEdges[2][1] = 10.;
 171
 172   // Input slot #0 works with a TChain ESD
 173   DefineInput(0, TChain::Class());
 174   // Output slot #0, #1 write into a TList
 175   DefineOutput(0, TList::Class());
 176 }
 177
 178 //________________________________________________________________________
 179 void AliPWG4HighPtQAMC::SetPtBinEdges(Int_t region, Double_t ptmax, Double_t ptBinWidth) {
 180
 181   if(region<3) {
 182     fPtBinEdges[region][0] = ptmax;
 183     fPtBinEdges[region][1] = ptBinWidth;
 184   }
 185   else {
 186     AliError("Only 3 regions alowed. Use region 0/1/2\n");
 187     return;
 188   }
 189
 190 }
 191
 192 //________________________________________________________________________
 193 void AliPWG4HighPtQAMC::ConnectInputData(Option_t *)
 194 {
 195   // Connect ESD and MC event handler here
 196   // Called once
 197   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 198   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 199   if (!tree) {
 200     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
 201     return;
 202   }
 203
 204   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 205
 206   if (!esdH) {
 207     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 208     return;
 209   } else
 210     fESD = esdH->GetEvent();
 211
 212   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 213   if (!eventHandler) {
 214     AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
 215   }
 216   else
 217     fMC = eventHandler->MCEvent();
 218
 219 }
 220
 221
 222 //________________________________________________________________________
 223 void AliPWG4HighPtQAMC::CreateOutputObjects() {
 224   //Create output objects
 225   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::CreateOutputObjects \n"));
 226
 227   Bool_t oldStatus = TH1::AddDirectoryStatus();
 228   TH1::AddDirectory(kFALSE);
 229
 230   OpenFile(0);
 231   fHistList = new TList();
 232   fHistList->SetOwner(kTRUE);
 233
 234   Float_t fgkPtMin = 0.;
 235
 236   //fPtBinEdges[region][0] = ptmax of region ; fPtBinEdges[region][1] = binWidth of region
 237   const Float_t ptmin1 =  fgkPtMin;
 238   const Float_t ptmax1 =  fPtBinEdges[0][0];
 239   const Float_t ptmin2 =  ptmax1 ;
 240   const Float_t ptmax2 =  fPtBinEdges[1][0];
 241   const Float_t ptmin3 =  ptmax2 ;
 242   const Float_t ptmax3 =  fPtBinEdges[2][0];//fgkPtMax;
 243   const Int_t nbin11 = (int)((ptmax1-ptmin1)/fPtBinEdges[0][1]);
 244   const Int_t nbin12 = (int)((ptmax2-ptmin2)/fPtBinEdges[1][1])+nbin11;
 245   const Int_t nbin13 = (int)((ptmax3-ptmin3)/fPtBinEdges[2][1])+nbin12;
 246   Int_t fgkNPtBins=nbin13;
 247   //Create array with low edges of each bin
 248   Double_t *binsPt=new Double_t[fgkNPtBins+1];
 249   for(Int_t i=0; i<=fgkNPtBins; i++) {
 250     if(i<=nbin11) binsPt[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
 251     if(i<=nbin12 && i>nbin11) binsPt[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 252     if(i<=nbin13 && i>nbin12) binsPt[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 253   }
 254
 255   Int_t fgkNPhiBins = 18*6;
 256   Float_t kMinPhi   = 0.;
 257   Float_t kMaxPhi   = 2.*TMath::Pi();
 258   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 259   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 260
 261   Int_t fgkResPtBins=80;
 262   Float_t kMinResPt = -1.;
 263   Float_t kMaxResPt = 1.;
 264   Double_t *binsResPt = new Double_t[fgkResPtBins+1];
 265   for(Int_t i=0; i<=fgkResPtBins; i++) binsResPt[i]=(Double_t)kMinResPt + (kMaxResPt-kMinResPt)/fgkResPtBins*(Double_t)i ;
 266
 267   Int_t fgkMultBins=50;
 268   Float_t kMinMult = 0.;
 269   Float_t kMaxMult = 3000.;
 270   Double_t *binsMult = new Double_t[fgkMultBins+1];
 271   for(Int_t i=0; i<=fgkMultBins; i++) binsMult[i]=(Double_t)kMinMult + (kMaxMult-kMinMult)/fgkMultBins*(Double_t)i ;
 272
 273   Int_t fgkNEtaBins=20;
 274   Float_t fgkEtaMin = -1.;
 275   Float_t fgkEtaMax = 1.;
 276   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 277   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 278
 279   Int_t fgkNNClustersTPCBins=80;
 280   Float_t fgkNClustersTPCMin = 0.5;
 281   Float_t fgkNClustersTPCMax = 160.5;
 282   Double_t *binsNClustersTPC=new Double_t[fgkNNClustersTPCBins+1];
 283   for(Int_t i=0; i<=fgkNNClustersTPCBins; i++) binsNClustersTPC[i]=(Double_t)fgkNClustersTPCMin + (fgkNClustersTPCMax-fgkNClustersTPCMin)/fgkNNClustersTPCBins*(Double_t)i ;
 284
 285   Int_t fgkNDCA2DBins=80;
 286   Float_t fgkDCA2DMin = -0.2;
 287   Float_t fgkDCA2DMax = 0.2;
 288   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
 289     fgkDCA2DMin = -2.;
 290     fgkDCA2DMax = 2.;
 291   }
 292   Double_t *binsDCA2D=new Double_t[fgkNDCA2DBins+1];
 293   for(Int_t i=0; i<=fgkNDCA2DBins; i++) binsDCA2D[i]=(Double_t)fgkDCA2DMin + (fgkDCA2DMax-fgkDCA2DMin)/fgkNDCA2DBins*(Double_t)i ;
 294
 295   Int_t fgkNDCAZBins=80;
 296   Float_t fgkDCAZMin = -2.;
 297   Float_t fgkDCAZMax = 2.;
 298   if(fTrackType==1 || fTrackType==2 || fTrackType==4) {
 299     fgkDCAZMin = -5.;
 300     fgkDCAZMax = 5.;
 301   }
 302   Double_t *binsDCAZ=new Double_t[fgkNDCAZBins+1];
 303   for(Int_t i=0; i<=fgkNDCAZBins; i++) binsDCAZ[i]=(Double_t)fgkDCAZMin + (fgkDCAZMax-fgkDCAZMin)/fgkNDCAZBins*(Double_t)i ;
 304
 305  Int_t fgkNNPointITSBins=9;
 306   Float_t fgkNPointITSMin = -0.5;
 307   Float_t fgkNPointITSMax = 8.5;
 308   Double_t *binsNPointITS=new Double_t[fgkNNPointITSBins+1];
 309   for(Int_t i=0; i<=fgkNNPointITSBins; i++) binsNPointITS[i]=(Double_t)fgkNPointITSMin + (fgkNPointITSMax-fgkNPointITSMin)/fgkNNPointITSBins*(Double_t)i ;
 310
 311   Int_t fgkNNSigmaToVertexBins=20;
 312   Float_t fgkNSigmaToVertexMin = 0.;
 313   Float_t fgkNSigmaToVertexMax = 8.;
 314   Double_t *binsNSigmaToVertex=new Double_t[fgkNNSigmaToVertexBins+1];
 315   for(Int_t i=0; i<=fgkNNSigmaToVertexBins; i++) binsNSigmaToVertex[i]=(Double_t)fgkNSigmaToVertexMin + (fgkNSigmaToVertexMax-fgkNSigmaToVertexMin)/fgkNNSigmaToVertexBins*(Double_t)i ;
 316
 317   Int_t fgkNChi2CBins=20;
 318   Float_t fgkChi2CMin = 0.;
 319   Float_t fgkChi2CMax = 100.;
 320   Double_t *binsChi2C=new Double_t[fgkNChi2CBins+1];
 321   for(Int_t i=0; i<=fgkNChi2CBins; i++) binsChi2C[i]=(Double_t)fgkChi2CMin + (fgkChi2CMax-fgkChi2CMin)/fgkNChi2CBins*(Double_t)i ;
 322
 323   Int_t fgkNRel1PtUncertaintyBins=50;
 324   Float_t fgkRel1PtUncertaintyMin = 0.;
 325   Float_t fgkRel1PtUncertaintyMax = 1.;
 326   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
 327   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
 328
 329   Float_t fgkChi2PerClusMin = 0.;
 330   Float_t fgkChi2PerClusMax = 4.;
 331   Int_t fgkNChi2PerClusBins = (int)(fgkChi2PerClusMax*10.);
 332   Double_t *binsChi2PerClus=new Double_t[fgkNChi2PerClusBins+1];
 333   for(Int_t i=0; i<=fgkNChi2PerClusBins; i++) binsChi2PerClus[i]=(Double_t)fgkChi2PerClusMin + (fgkChi2PerClusMax-fgkChi2PerClusMin)/fgkNChi2PerClusBins*(Double_t)i ;
 334
 335
 336   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 337   fHistList->Add(fNEventAll);
 338   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 339   fHistList->Add(fNEventSel);
 340   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 341   //Set labels
 342   fNEventReject->Fill("noESD",0);
 343   fNEventReject->Fill("Trigger",0);
 344   fNEventReject->Fill("noMCEvent",0);
 345   fNEventReject->Fill("noStack",0);
 346   fNEventReject->Fill("NTracks<2",0);
 347   fNEventReject->Fill("noVTX",0);
 348   fNEventReject->Fill("VtxStatus",0);
 349   fNEventReject->Fill("NCont<2",0);
 350   fNEventReject->Fill("ZVTX>10",0);
 351   fHistList->Add(fNEventReject);
 352
 353   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 354   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 355   fHistList->Add(fh1Xsec);
 356
 357   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 358   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 359   fHistList->Add(fh1Trials);
 360
 361   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 362   fHistList->Add(fh1PtHard);
 363   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 364   fHistList->Add(fh1PtHardTrials);
 365
 366   fPtAll = new TH1F("fPtAll","PtAll",fgkNPtBins, binsPt);
 367   fHistList->Add(fPtAll);
 368   fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, binsPt);
 369   fHistList->Add(fPtSel);
 370
 371   fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, binsPt);
 372   fHistList->Add(fPtSelFakes);
 373
 374   fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",fgkNNClustersTPCBins, binsNClustersTPC);
 375   fHistList->Add(fNPointTPCFakes);
 376
 377   fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, binsPt);
 378   fHistList->Add(fPtSelLargeLabel);
 379
 380   fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",fgkMultBins, binsMult);
 381   fHistList->Add(fMultRec);
 382
 383   fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",fgkNNClustersTPCBins, binsNClustersTPC);
 384   fHistList->Add(fNPointTPCMultRec);
 385
 386   fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",50,0.,50.,40,-20.,20.);
 387   fHistList->Add(fDeltaPtMultRec);
 388
 389   fPtAllvsPtMC = new TH2F("fPtAllvsPtMC","fPtAllvsPtMC;p_{T,MC};p_{T,rec}",fgkNPtBins, binsPt,fgkNPtBins, binsPt);
 390   fHistList->Add(fPtAllvsPtMC);
 391
 392   fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, binsPt,fgkResPtBins,binsResPt);
 393   fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
 394   fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 395   fHistList->Add(fPtAllminPtMCvsPtAll);
 396
 397   fPtAllvsPtMCvsMult = new TH3F("fPtAllvsPtMCvsMult","fPtAllvsPtMCvsMult;p_{T,MC};p_{T,rec};N_{tracks}",fgkNPtBins, binsPt,fgkNPtBins, binsPt,fgkMultBins,binsMult);
 398   fHistList->Add(fPtAllvsPtMCvsMult);
 399
 400   fPtAllminPtMCvsPtAllvsMult = new TH3F("fPtAllminPtMCvsPtAllvsMult","fPtAllminPtMCvsPtAllvsMult",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkMultBins,binsMult);
 401   fPtAllminPtMCvsPtAllvsMult->SetXTitle("p_{t}^{MC}");
 402   fPtAllminPtMCvsPtAllvsMult->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 403   fPtAllminPtMCvsPtAllvsMult->SetZTitle("N_{tracks}");
 404   fHistList->Add(fPtAllminPtMCvsPtAllvsMult);
 405
 406   fPtAllminPtMCvsPtAllNPointTPC = new TH3F("fPtAllminPtMCvsPtAllNPointTPC","PtAllminPtMCvsPtAllNPointTPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
 407   fPtAllminPtMCvsPtAllNPointTPC->SetXTitle("p_{t}^{MC}");
 408   fPtAllminPtMCvsPtAllNPointTPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 409   fPtAllminPtMCvsPtAllNPointTPC->SetZTitle("N_{cls,TPC}");
 410   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPC);
 411
 412   fPtAllminPtMCvsPtAllNPointTPCIter1 = new TH3F("fPtAllminPtMCvsPtAllNPointTPCIter1","PtAllminPtMCvsPtAllNPointTPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNClustersTPCBins, binsNClustersTPC);
 413   fPtAllminPtMCvsPtAllNPointTPCIter1->SetXTitle("p_{t}^{MC}");
 414   fPtAllminPtMCvsPtAllNPointTPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 415   fPtAllminPtMCvsPtAllNPointTPCIter1->SetZTitle("N_{clsIter1,TPC}");
 416   fHistList->Add(fPtAllminPtMCvsPtAllNPointTPCIter1);
 417
 418   fPtAllminPtMCvsPtAllChi2TPC = new TH3F("fPtAllminPtMCvsPtAllChi2TPC","PtAllminPtMCvsPtAllChi2TPC",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
 419   fPtAllminPtMCvsPtAllChi2TPC->SetXTitle("p_{t}^{MC}");
 420   fPtAllminPtMCvsPtAllChi2TPC->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 421   fPtAllminPtMCvsPtAllChi2TPC->SetZTitle("#chi^{2} TPC");
 422   fHistList->Add(fPtAllminPtMCvsPtAllChi2TPC);
 423
 424   fPtAllminPtMCvsPtAllChi2TPCIter1 = new TH3F("fPtAllminPtMCvsPtAllChi2TPCIter1","PtAllminPtMCvsPtAllChi2TPCIter1",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2PerClusBins, binsChi2PerClus);
 425   fPtAllminPtMCvsPtAllChi2TPCIter1->SetXTitle("p_{t}^{MC}");
 426   fPtAllminPtMCvsPtAllChi2TPCIter1->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 427   fPtAllminPtMCvsPtAllChi2TPCIter1->SetZTitle("#chi^{2} TPC Iter1");
 428   fHistList->Add(fPtAllminPtMCvsPtAllChi2TPCIter1);
 429
 430   fPtAllminPtMCvsPtAllDCAR = new TH3F("fPtAllminPtMCvsPtAllDCAR","PtAllminPtMCvsPtAllDCAR",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCA2DBins,binsDCA2D);
 431   fPtAllminPtMCvsPtAllDCAR->SetXTitle("p_{t}^{MC}");
 432   fPtAllminPtMCvsPtAllDCAR->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 433   fPtAllminPtMCvsPtAllDCAR->SetZTitle("DCA_{R}");
 434   fHistList->Add(fPtAllminPtMCvsPtAllDCAR);
 435
 436   fPtAllminPtMCvsPtAllDCAZ = new TH3F("fPtAllminPtMCvsPtAllDCAZ","PtAllminPtMCvsPtAllDCAZ",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNDCAZBins,binsDCAZ);
 437   fPtAllminPtMCvsPtAllDCAZ->SetXTitle("p_{t}^{MC}");
 438   fPtAllminPtMCvsPtAllDCAZ->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 439   fPtAllminPtMCvsPtAllDCAZ->SetZTitle("DCA_{Z}");
 440   fHistList->Add(fPtAllminPtMCvsPtAllDCAZ);
 441
 442   fPtAllminPtMCvsPtAllPhi = new TH3F("fPtAllminPtMCvsPtAllPhi","PtAllminPtMCvsPtAllPhi",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNPhiBins,binsPhi);
 443   fPtAllminPtMCvsPtAllPhi->SetXTitle("p_{t}^{MC}");
 444   fPtAllminPtMCvsPtAllPhi->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 445   fPtAllminPtMCvsPtAllPhi->SetZTitle("#phi");
 446   fHistList->Add(fPtAllminPtMCvsPtAllPhi);
 447
 448   fPtAllminPtMCvsPtAllNPointITS = new TH3F("fPtAllminPtMCvsPtAllNPointITS","PtAllminPtMCvsPtAllNPointITS",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNPointITSBins,binsNPointITS);
 449   fPtAllminPtMCvsPtAllNPointITS->SetXTitle("p_{t}^{MC}");
 450   fPtAllminPtMCvsPtAllNPointITS->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 451   fPtAllminPtMCvsPtAllNPointITS->SetZTitle("N_{point,ITS}}");
 452   fHistList->Add(fPtAllminPtMCvsPtAllNPointITS);
 453
 454   fPtAllminPtMCvsPtAllNSigmaToVertex = new TH3F("fPtAllminPtMCvsPtAllNSigmaToVertex","PtAllminPtMCvsPtAllNSigmaToVertex",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNNSigmaToVertexBins,binsNSigmaToVertex);
 455   fPtAllminPtMCvsPtAllNSigmaToVertex->SetXTitle("p_{t}^{MC}");
 456   fPtAllminPtMCvsPtAllNSigmaToVertex->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 457   fPtAllminPtMCvsPtAllNSigmaToVertex->SetZTitle("N#sigma to vertex");
 458   fHistList->Add(fPtAllminPtMCvsPtAllNSigmaToVertex);
 459
 460   fPtAllminPtMCvsPtAllChi2C = new TH3F("fPtAllminPtMCvsPtAllChi2C","PtAllminPtMCvsPtAllChi2C",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNChi2CBins,binsChi2C);
 461   fPtAllminPtMCvsPtAllChi2C->SetXTitle("p_{t}^{MC}");
 462   fPtAllminPtMCvsPtAllChi2C->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 463   fPtAllminPtMCvsPtAllChi2C->SetZTitle("Constrained #chi^{2}");
 464   fHistList->Add(fPtAllminPtMCvsPtAllChi2C);
 465
 466   fPtAllminPtMCvsPtAllRel1PtUncertainty = new TH3F("fPtAllminPtMCvsPtAllRel1PtUncertainty","PtAllminPtMCvsPtAllRel1PtUncertainty",fgkNPtBins, binsPt,fgkResPtBins,binsResPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 467   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetXTitle("p_{t}^{MC}");
 468   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
 469   fPtAllminPtMCvsPtAllRel1PtUncertainty->SetZTitle("Rel1PtUncertainty");
 470   fHistList->Add(fPtAllminPtMCvsPtAllRel1PtUncertainty);
 471
 472
 473   fPtAllMC = new TH1F("fPtAllMC","PtAll",fgkNPtBins, binsPt);
 474   fHistList->Add(fPtAllMC);
 475   fPtSelMC = new TH1F("fPtSelMC","PtSel",fgkNPtBins, binsPt);
 476   fHistList->Add(fPtSelMC);
 477
 478   TH1::AddDirectory(oldStatus);
 479
 480   PostData(0, fHistList);
 481
 482   if(binsPt)                delete [] binsPt;
 483   if(binsResPt)             delete [] binsResPt;
 484   if(binsPhi)               delete [] binsPhi;
 485   if(binsNClustersTPC)      delete [] binsNClustersTPC;
 486   if(binsDCA2D)             delete [] binsDCA2D;
 487   if(binsDCAZ)              delete [] binsDCAZ;
 488   if(binsNPointITS)         delete [] binsNPointITS;
 489   if(binsNSigmaToVertex)    delete [] binsNSigmaToVertex;
 490   if(binsChi2C)             delete [] binsChi2C;
 491   if(binsEta)               delete [] binsEta;
 492   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 493   if(binsChi2PerClus)       delete [] binsChi2PerClus;
 494   if(binsMult)              delete [] binsMult;
 495
 496 }
 497
 498 //________________________________________________________________________
 499 Bool_t AliPWG4HighPtQAMC::SelectEvent() {
 500   //
 501   // Decide if event should be selected for analysis
 502   //
 503
 504   // Checks following requirements:
 505   // - fESD available
 506   // - trigger info from AliPhysicsSelection
 507   // - MCevent available
 508   // - number of reconstructed tracks > 1
 509   // - primary vertex reconstructed
 510   // - z-vertex < 10 cm
 511
 512   Bool_t selectEvent = kTRUE;
 513
 514   //fESD object available?
 515   if (!fESD) {
 516     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 517     fNEventReject->Fill("noESD",1);
 518     selectEvent = kFALSE;
 519     return selectEvent;
 520   }
 521
 522   //Trigger
 523   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 524   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 525     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 526     fNEventReject->Fill("Trigger",1);
 527     selectEvent = kFALSE;
 528     return selectEvent;
 529   }
 530
 531   //MCEvent available?
 532   //if yes: get stack
 533   if(fMC) {
 534     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
 535     fStack = fMC->Stack();                //Particles Stack
 536     if(fStack) {
 537       AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
 538     }
 539     else {
 540       AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
 541       fNEventReject->Fill("noStack",1);
 542       selectEvent = kFALSE;
 543       return selectEvent;
 544     }
 545   } else {
 546     AliDebug(2,Form("ERROR: Could not retrieve stack"));
 547     fNEventReject->Fill("noMCEvent",1);
 548     selectEvent = kFALSE;
 549     return selectEvent;
 550   }
 551
 552   //Check if number of reconstructed tracks is larger than 1
 553   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 554     fNEventReject->Fill("NTracks<2",1);
 555     selectEvent = kFALSE;
 556     return selectEvent;
 557   }
 558
 559   //Check if vertex is reconstructed
 560   if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
 561   else              fVtx = fESD->GetPrimaryVertexSPD();
 562
 563   if(!fVtx) {
 564     fNEventReject->Fill("noVTX",1);
 565     selectEvent = kFALSE;
 566     return selectEvent;
 567   }
 568
 569   if(!fVtx->GetStatus()) {
 570     fNEventReject->Fill("VtxStatus",1);
 571     selectEvent = kFALSE;
 572     return selectEvent;
 573   }
 574
 575   // Need vertex cut
 576   //  TString vtxName(fVtx->GetName());
 577   if(fVtx->GetNContributors()<2) {
 578     fNEventReject->Fill("NCont<2",1);
 579     selectEvent = kFALSE;
 580     return selectEvent;
 581   }
 582
 583   //Check if z-vertex < 10 cm
 584   double primVtx[3];
 585   fVtx->GetXYZ(primVtx);
 586   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 587     fNEventReject->Fill("ZVTX>10",1);
 588     selectEvent = kFALSE;
 589     return selectEvent;
 590   }
 591
 592   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
 593
 594   return selectEvent;
 595
 596 }
 597
 598 //________________________________________________________________________
 599 void AliPWG4HighPtQAMC::Exec(Option_t *) {
 600   // Main loop
 601   // Called for each event
 602   AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
 603   // All events without selection
 604   fNEventAll->Fill(0.);
 605
 606   if(!SelectEvent()) {
 607     // Post output data
 608     PostData(0, fHistList);
 609     return;
 610   }
 611
 612   // ---- Get MC Header information (for MC productions in pThard bins) ----
 613   Double_t ptHard = 0.;
 614   Double_t nTrials = 1; // trials for MC trigger weight for real data
 615
 616   if(fMC){
 617     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
 618      if(pythiaGenHeader){
 619        nTrials = pythiaGenHeader->Trials();
 620        ptHard  = pythiaGenHeader->GetPtHard();
 621
 622        fh1PtHard->Fill(ptHard);
 623        fh1PtHardTrials->Fill(ptHard,nTrials);
 624
 625        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 626      }
 627   }
 628
 629   //Need to keep track of selected events
 630   fNEventSel->Fill(0.);
 631
 632   Int_t nTracks = fESD->GetNumberOfTracks();
 633   AliDebug(2,Form("nTracks ESD%d", nTracks));
 634
 635   int nMCtracks = fStack->GetNtrack();
 636
 637   Float_t pt      = 0.;
 638   Float_t ptMC    = 0.;
 639   Float_t phi     = 0.;
 640   Float_t dca2D   = 0.;
 641   Float_t dcaZ    = 0.;
 642   Int_t nPointITS = 0;
 643   Float_t chi2C   = 0.;
 644   Float_t nSigmaToVertex    = 0.;
 645   Float_t relUncertainty1Pt = 0.;
 646
 647   int mult = fTrackCuts->CountAcceptedTracks(fESD);
 648
 649   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 650
 651     AliESDtrack *track = 0;
 652     AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 653     if(!esdtrack) continue;
 654
 655     if(fTrackType==1) {
 656       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 657       if(!track) continue;
 658     }
 659     else if(fTrackType==2) {
 660       track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 661       if(!track) continue;
 662
 663       AliExternalTrackParam exParam;
 664       Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
 665       if( !relate ) {
 666         delete track;
 667         continue;
 668       }
 669       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 670     }
 671     else if(fTrackType==7) {
 672       //use global constrained track
 673       track = new AliESDtrack(*esdtrack);
 674       //      track = esdtrack;
 675       //      track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 676     }
 677     else
 678       track = esdtrack;
 679
 680     if(!track) continue;
 681
 682     if(fTrackType==2) {
 683       //Cut on chi2 of constrained fit
 684       if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
 685         delete track;
 686         continue;
 687       }
 688     }
 689
 690     Int_t label = TMath::Abs(track->GetLabel());
 691     if(label>=nMCtracks) {
 692       if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
 693       if(fTrackType==1 || fTrackType==2 || fTrackType==7) delete track;
 694       continue;
 695     }
 696
 697     TParticle *particle = fStack->Particle(label) ;
 698     if(!particle) {
 699       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
 700         if(track) delete track;
 701       }
 702       continue;
 703     }
 704
 705     ptMC = particle->Pt();
 706
 707     if(fTrackType==1 || fTrackType==2)
 708       track->GetImpactParametersTPC(dca2D,dcaZ);  //TPConly
 709     else
 710       track->GetImpactParameters(dca2D,dcaZ);     //Global
 711
 712     UChar_t itsMap = track->GetITSClusterMap();
 713     for (Int_t i=0; i < 6; i++) {
 714       if (itsMap & (1 << i))
 715         nPointITS ++;
 716     }
 717
 718     //    fPtAll->Fill(pt);
 719     fPtAllMC->Fill(ptMC);
 720
 721     if (fTrackCuts->AcceptTrack(track)) {
 722
 723       if(fTrackType==7) {
 724         if(fTrackCutsReject ) {
 725           if(fTrackCutsReject->AcceptTrack(track) ) {
 726             if(track) delete track;
 727             continue;
 728           }
 729         }
 730
 731         if(esdtrack->GetConstrainedParam())
 732           track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 733       }
 734
 735       pt  = track->Pt();
 736       phi = track->Phi();
 737
 738       nSigmaToVertex = fTrackCuts->GetSigmaToVertex(track);// Calculates the number of sigma to the vertex for a track.
 739       chi2C = track->GetConstrainedChi2();
 740       relUncertainty1Pt = TMath::Sqrt(track->GetSigma1Pt2())*pt;
 741
 742       fPtSel->Fill(pt);
 743       if(track->GetLabel()<0) {
 744         fPtSelFakes->Fill(pt);
 745         fNPointTPCFakes->Fill(track->GetTPCNcls());
 746       }
 747       fPtSelMC->Fill(ptMC);
 748       fPtAllvsPtMC->Fill(ptMC,pt);
 749       fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
 750       fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
 751       fPtAllminPtMCvsPtAllNPointTPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNclsIter1());
 752       if(track->GetTPCNcls()>0.)
 753         fPtAllminPtMCvsPtAllChi2TPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2()/track->GetTPCNcls());
 754       if(track->GetTPCNclsIter1()>0.)
 755         fPtAllminPtMCvsPtAllChi2TPCIter1->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCchi2Iter1()/track->GetTPCNclsIter1());
 756       fPtAllminPtMCvsPtAllDCAR->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dca2D);
 757       fPtAllminPtMCvsPtAllDCAZ->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),dcaZ);
 758       fPtAllminPtMCvsPtAllPhi->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),phi);
 759       fPtAllminPtMCvsPtAllNPointITS->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nPointITS);
 760       fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
 761       fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
 762       fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
 763
 764       fPtAllvsPtMCvsMult->Fill(ptMC,pt,mult);
 765       fPtAllminPtMCvsPtAllvsMult->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC), mult);
 766
 767       //Check if track is reconstructed multiple times
 768       /*
 769       int multCounter = 1;
 770       for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
 771         //   AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
 772         AliESDtrack *track2;
 773         AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
 774         if(!esdtrack2) continue;
 775
 776         if(fTrackType==1)
 777           track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
 778         else if(fTrackType==2) {
 779           track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
 780           if(!track2) {
 781             continue;
 782           }
 783           AliExternalTrackParam exParam2;
 784           Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
 785           if( !relate ) {
 786             delete track2;
 787             continue;
 788           }
 789           track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
 790         }
 791         else
 792           track2 = esdtrack2;
 793         if(!track2) {
 794           continue;
 795         }
 796
 797         if (fTrackCuts->AcceptTrack(track2)) {
 798           Int_t label2 = TMath::Abs(track2->GetLabel());
 799           if(label==label2) {
 800             fNPointTPCMultRec->Fill(track->GetTPCNcls());
 801             fNPointTPCMultRec->Fill(track2->GetTPCNcls());
 802             fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
 803             multCounter++;
 804           }
 805         }
 806         if(fTrackType==1 || fTrackType==2) delete track2;
 807       }//track2 loop
 808       if(multCounter>1) fMultRec->Fill(multCounter);
 809       */
 810
 811     }//fTrackCuts selection
 812
 813
 814     if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
 815       if(track) delete track;
 816     }
 817
 818   }//ESD track loop
 819
 820   // Post output data
 821   PostData(0, fHistList);
 822
 823 }
 824 //________________________________________________________________________
 825 Bool_t AliPWG4HighPtQAMC::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 826   //
 827   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 828   // This is to called in Notify and should provide the path to the AOD/ESD file
 829   // Copied from AliAnalysisTaskJetSpectrum2
 830   //
 831
 832   TString file(currFile);
 833   fXsec = 0;
 834   fTrials = 1;
 835
 836   if(file.Contains("root_archive.zip#")){
 837     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 838     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 839     file.Replace(pos+1,20,"");
 840   }
 841   else {
 842     // not an archive take the basename....
 843     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 844   }
 845   //  Printf("%s",file.Data());
 846
 847
 848   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
 849   if(!fxsec){
 850     // next trial fetch the histgram file
 851     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 852     if(!fxsec){
 853         // not a severe condition but inciate that we have no information
 854       return kFALSE;
 855     }
 856     else{
 857       // find the tlist we want to be independtent of the name so use the Tkey
 858       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 859       if(!key){
 860         fxsec->Close();
 861         return kFALSE;
 862       }
 863       TList *list = dynamic_cast<TList*>(key->ReadObj());
 864       if(!list){
 865         fxsec->Close();
 866         return kFALSE;
 867       }
 868       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 869       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 870       fxsec->Close();
 871     }
 872   } // no tree pyxsec.root
 873   else {
 874     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 875     if(!xtree){
 876       fxsec->Close();
 877       return kFALSE;
 878     }
 879     UInt_t   ntrials  = 0;
 880     Double_t  xsection  = 0;
 881     xtree->SetBranchAddress("xsection",&xsection);
 882     xtree->SetBranchAddress("ntrials",&ntrials);
 883     xtree->GetEntry(0);
 884     fTrials = ntrials;
 885     fXsec = xsection;
 886     fxsec->Close();
 887   }
 888   return kTRUE;
 889 }
 890 //________________________________________________________________________
 891 Bool_t AliPWG4HighPtQAMC::Notify()
 892 {
 893   //
 894   // Implemented Notify() to read the cross sections
 895   // and number of trials from pyxsec.root
 896   // Copied from AliAnalysisTaskJetSpectrum2
 897   //
 898
 899   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 900   Float_t xsection = 0;
 901   Float_t ftrials  = 1;
 902
 903   fAvgTrials = 1;
 904   if(tree){
 905     TFile *curfile = tree->GetCurrentFile();
 906     if (!curfile) {
 907       Error("Notify","No current file");
 908       return kFALSE;
 909     }
 910     if(!fh1Xsec||!fh1Trials){
 911       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 912       return kFALSE;
 913     }
 914     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 915     fh1Xsec->Fill("<#sigma>",xsection);
 916     // construct a poor man average trials
 917     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 918     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 919   }
 920   return kTRUE;
 921 }
 922
 923 //________________________________________________________________________
 924 AliGenPythiaEventHeader*  AliPWG4HighPtQAMC::GetPythiaEventHeader(AliMCEvent *mcEvent){
 925
 926   if(!mcEvent)return 0;
 927   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 928   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 929   if(!pythiaGenHeader){
 930     // cocktail ??
 931     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 932
 933     if (!genCocktailHeader) {
 934       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 935       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 936       return 0;
 937     }
 938     TList* headerList = genCocktailHeader->GetHeaders();
 939     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 940       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 941       if (pythiaGenHeader)
 942         break;
 943     }
 944     if(!pythiaGenHeader){
 945       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 946       return 0;
 947     }
 948   }
 949   return pythiaGenHeader;
 950
 951 }
 952
 953 //________________________________________________________________________
 954 void AliPWG4HighPtQAMC::Terminate(Option_t *)
 955 {
 956   // The Terminate() function is the last function to be called during
 957   // a query. It always runs on the client, it can be used to present
 958   // the results graphically or save the results to file.
 959
 960 }
 961
 962 #endif