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