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