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