PWG4/JetTasks/AliPWG4HighPtSpectra.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 // Efficiency between different steps of the procedure.
  18 // The ouptut of the task is a AliCFContainer from which the efficiencies
  19 // can be calculated
  20 //-----------------------------------------------------------------------
  21 // Author : Marta Verweij - UU
  22 //-----------------------------------------------------------------------
  23
  24
  25 #ifndef ALIPWG4HIGHPTSPECTRA_CXX
  26 #define ALIPWG4HIGHPTSPECTRA_CXX
  27
  28 #include "AliPWG4HighPtSpectra.h"
  29
  30 #include "TVector3.h"
  31 #include <iostream>
  32 #include "TH1.h"
  33 #include "TH2.h"
  34 #include "TH3.h"
  35 #include "TProfile.h"
  36 #include "TList.h"
  37 #include "TChain.h"
  38 #include "TKey.h"
  39 #include "TSystem.h"
  40 #include "TFile.h"
  41
  42 #include "AliAnalysisManager.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliESDtrack.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliExternalTrackParam.h"
  47 #include "AliCentrality.h"
  48
  49 #include "AliLog.h"
  50
  51 #include "AliStack.h"
  52 #include "TParticle.h"
  53 #include "AliMCEvent.h"
  54 #include "AliMCEventHandler.h"
  55 #include "AliCFContainer.h"
  56 #include "AliGenPythiaEventHeader.h"
  57 #include "AliGenCocktailEventHeader.h"
  58 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
  59
  60 //#include <iostream>
  61 using namespace std; //required for resolving the 'cout' symbol
  62
  63 ClassImp(AliPWG4HighPtSpectra)
  64
  65 //__________________________________________________________________________
  66 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
  67   fReadAODData(0),
  68   fCFManagerPos(0x0),
  69   fCFManagerNeg(0x0),
  70   fESD(0x0),
  71   fMC(0x0),
  72   fStack(0x0),
  73   fVtx(0x0),
  74   fIsPbPb(0),
  75   fCentClass(10),
  76   fTrackType(0),
  77   fTrackCuts(0x0),
  78   fTrackCutsReject(0x0),
  79   fSigmaConstrainedMax(5.),
  80   fAvgTrials(1),
  81   fHistList(0),
  82   fNEventAll(0),
  83   fNEventSel(0),
  84   fNEventReject(0),
  85   fh1Centrality(0x0),
  86   fh1Xsec(0),
  87   fh1Trials(0),
  88   fh1PtHard(0),
  89   fh1PtHardTrials(0)
  90 {
  91   //
  92   //Default ctor
  93   //
  94 }
  95 //___________________________________________________________________________
  96 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
  97   AliAnalysisTask(name,""),
  98   fReadAODData(0),
  99   fCFManagerPos(0x0),
 100   fCFManagerNeg(0x0),
 101   fESD(0x0),
 102   fMC(0x0),
 103   fStack(0x0),
 104   fVtx(0x0),
 105   fIsPbPb(0),
 106   fCentClass(10),
 107   fTrackType(0),
 108   fTrackCuts(0x0),
 109   fTrackCutsReject(0x0),
 110   fSigmaConstrainedMax(5.),
 111   fAvgTrials(1),
 112   fHistList(0),
 113   fNEventAll(0),
 114   fNEventSel(0),
 115   fNEventReject(0),
 116   fh1Centrality(0x0),
 117   fh1Xsec(0),
 118   fh1Trials(0),
 119   fh1PtHard(0),
 120   fh1PtHardTrials(0)
 121 {
 122   //
 123   // Constructor. Initialization of Inputs and Outputs
 124   //
 125   AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
 126   // Input slot #0 works with a TChain ESD
 127   DefineInput(0, TChain::Class());
 128   // Output slot #0 writes into a TList
 129   DefineOutput(0,TList::Class());
 130   // Output slot #1, #2 writes into a AliCFContainer
 131   DefineOutput(1,AliCFContainer::Class());
 132   DefineOutput(2,AliCFContainer::Class());
 133   // Output slot #3 writes into a AliESDtrackCuts
 134   DefineOutput(3, AliESDtrackCuts::Class());
 135   DefineOutput(4, AliESDtrackCuts::Class());
 136 }
 137
 138 //________________________________________________________________________
 139 void AliPWG4HighPtSpectra::LocalInit()
 140 {
 141   //
 142   // Only called once at beginning
 143   //
 144   PostData(3,fTrackCuts);
 145 }
 146
 147 //________________________________________________________________________
 148 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
 149 {
 150   // Connect ESD here
 151   // Called once
 152   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 153
 154   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 155   if (!tree) {
 156     AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
 157     return;
 158   }
 159
 160   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 161
 162   if (!esdH) {
 163     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 164     return;
 165   } else
 166     fESD = esdH->GetEvent();
 167
 168   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 169   if (!eventHandler) {
 170     AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
 171   }
 172   else
 173     fMC = eventHandler->MCEvent();
 174
 175 }
 176
 177 //________________________________________________________________________
 178 Bool_t AliPWG4HighPtSpectra::SelectEvent() {
 179   //
 180   // Decide if event should be selected for analysis
 181   //
 182
 183   // Checks following requirements:
 184   // - fESD available
 185   // - trigger info from AliPhysicsSelection
 186   // - number of reconstructed tracks > 1
 187   // - primary vertex reconstructed
 188   // - z-vertex < 10 cm
 189
 190   Bool_t selectEvent = kTRUE;
 191
 192   //fESD object available?
 193   if (!fESD) {
 194     AliDebug(2,Form("ERROR: fInputEvent not available\n"));
 195     fNEventReject->Fill("noESD",1);
 196     selectEvent = kFALSE;
 197     return selectEvent;
 198   }
 199
 200   //Trigger
 201   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 202   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 203     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 204     fNEventReject->Fill("Trigger",1);
 205     selectEvent = kFALSE;
 206     return selectEvent;
 207   }
 208
 209   //Check if number of reconstructed tracks is larger than 1
 210   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2)  {
 211     fNEventReject->Fill("NTracks<2",1);
 212     selectEvent = kFALSE;
 213     return selectEvent;
 214   }
 215
 216   //Check if vertex is reconstructed
 217   fVtx = fESD->GetPrimaryVertexSPD();
 218
 219   if(!fVtx) {
 220     fNEventReject->Fill("noVTX",1);
 221     selectEvent = kFALSE;
 222     return selectEvent;
 223   }
 224
 225   if(!fVtx->GetStatus()) {
 226     fNEventReject->Fill("VtxStatus",1);
 227     selectEvent = kFALSE;
 228     return selectEvent;
 229   }
 230
 231   // Need vertex cut
 232   //  TString vtxName(fVtx->GetName());
 233   if(fVtx->GetNContributors()<2) {
 234     fNEventReject->Fill("NCont<2",1);
 235     selectEvent = kFALSE;
 236     return selectEvent;
 237   }
 238
 239   //Check if z-vertex < 10 cm
 240   double primVtx[3];
 241   fVtx->GetXYZ(primVtx);
 242   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 243     fNEventReject->Fill("ZVTX>10",1);
 244     selectEvent = kFALSE;
 245     return selectEvent;
 246   }
 247
 248   //Centrality selection should only be done in case of PbPb
 249   if(IsPbPb()) {
 250     Float_t cent = 0.;
 251     if(fCentClass!=CalculateCentrality(fESD) && fCentClass!=10) {
 252       fNEventReject->Fill("cent",1);
 253       selectEvent = kFALSE;
 254       return selectEvent;
 255     }
 256     else {
 257       if(dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()) {
 258         cent = dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()->GetCentralityPercentile("V0M");
 259       }
 260       if(cent>90.) {
 261         fNEventReject->Fill("cent>90",1);
 262         selectEvent = kFALSE;
 263         return selectEvent;
 264       }
 265       fh1Centrality->Fill(cent);
 266     }
 267   }
 268
 269   return selectEvent;
 270
 271 }
 272
 273 //________________________________________________________________________
 274 Int_t AliPWG4HighPtSpectra::CalculateCentrality(AliESDEvent *esd){
 275
 276
 277   Float_t cent = 999;
 278
 279   if(esd){
 280     if(esd->GetCentrality()){
 281       cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
 282     }
 283   }
 284
 285   if(cent<0)  return 5;
 286   if(cent>80)return 4;
 287   if(cent>50)return 3;
 288   if(cent>30)return 2;
 289   if(cent>10)return 1;
 290   return 0;
 291
 292 }
 293
 294 //_________________________________________________
 295 void AliPWG4HighPtSpectra::Exec(Option_t *)
 296 {
 297   //
 298   // Main loop function
 299   //
 300   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
 301
 302   // All events without selection
 303   fNEventAll->Fill(0.);
 304
 305   if(!SelectEvent()) {
 306     fNEventReject->Fill("NTracks<2",1);
 307     // Post output data
 308     PostData(0,fHistList);
 309     PostData(1,fCFManagerPos->GetParticleContainer());
 310     PostData(2,fCFManagerNeg->GetParticleContainer());
 311     return;
 312   }
 313
 314   //MCEvent available?
 315   //if yes: get stack
 316   if(fMC) {
 317     AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
 318     fStack = fMC->Stack();                //Particles Stack
 319     AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
 320   }
 321
 322   // ---- Get MC Header information (for MC productions in pThard bins) ----
 323   Double_t ptHard = 0.;
 324   Double_t nTrials = 1; // trials for MC trigger weight for real data
 325
 326   if(fMC){
 327     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
 328      if(pythiaGenHeader){
 329        nTrials = pythiaGenHeader->Trials();
 330        ptHard  = pythiaGenHeader->GetPtHard();
 331
 332        fh1PtHard->Fill(ptHard);
 333        fh1PtHardTrials->Fill(ptHard,nTrials);
 334
 335        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 336      }
 337   }
 338
 339   Int_t nTracks = fESD->GetNumberOfTracks();
 340   AliDebug(2,Form("nTracks %d", nTracks));
 341
 342   if(!fTrackCuts) {
 343     fNEventReject->Fill("noTrackCuts",1);
 344     // Post output data
 345     PostData(0,fHistList);
 346     PostData(1,fCFManagerPos->GetParticleContainer());
 347     PostData(2,fCFManagerNeg->GetParticleContainer());
 348     return;
 349   }
 350
 351   // Selected events for analysis
 352   fNEventSel->Fill(0.);
 353
 354
 355   Double_t containerInputRec[3]       = {0.,0.,0.};
 356   Double_t containerInputMC[3]        = {0.,0.,0.};
 357   Double_t containerInputRecMC[3]     = {0.,0.,0.}; //reconstructed yield as function of MC variable
 358
 359   //Now go to rec level
 360   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 361     {
 362       //Get track for analysis
 363       AliESDtrack *track = 0x0;
 364       AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 365       if(!esdtrack) continue;
 366
 367       if(fTrackType==1) {
 368         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 369         if(!track) continue;
 370       }
 371       else if(fTrackType==2) {
 372         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 373         if(!track) continue;
 374
 375         AliExternalTrackParam exParam;
 376         Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
 377         if( !relate ) {
 378           delete track;
 379           continue;
 380         }
 381         track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 382       }
 383       else if(fTrackType==7) {
 384         //use global constrained track
 385         track = esdtrack;
 386         //      track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 387
 388       }
 389       else
 390         track = esdtrack;
 391
 392
 393       if(fTrackType==2) {
 394         //Cut on chi2 of constrained fit
 395         if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
 396           delete track;
 397           continue;
 398         }
 399       }
 400
 401       if (fTrackCuts->AcceptTrack(track)) {
 402
 403         if(fTrackType==7) {
 404           if(fTrackCutsReject ) {
 405             if(fTrackCutsReject->AcceptTrack(track) )
 406               continue;
 407           }
 408
 409           if(esdtrack->GetConstrainedParam())
 410             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 411         }
 412
 413         //fill the container
 414         containerInputRec[0] = track->Pt();
 415         containerInputRec[1] = track->Phi();
 416         containerInputRec[2] = track->Eta();
 417
 418         if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 419         if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 420
 421
 422         //Only fill the MC containers if MC information is available
 423         if(fMC) {
 424           Int_t label = TMath::Abs(track->GetLabel());
 425           TParticle *particle = fStack->Particle(label) ;
 426           if(!particle) {
 427             if(fTrackType==1 || fTrackType==2)
 428               delete track;
 429             continue;
 430           }
 431           containerInputRecMC[0] = particle->Pt();
 432           containerInputRecMC[1] = particle->Phi();
 433           containerInputRecMC[2] = particle->Eta();
 434
 435           //Container with primaries
 436           if(fStack->IsPhysicalPrimary(label)) {
 437             if(particle->GetPDG()->Charge()>0.) {
 438               fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 439             }
 440             if(particle->GetPDG()->Charge()<0.) {
 441               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 442             }
 443           }
 444
 445           //Container with secondaries
 446           if (!fStack->IsPhysicalPrimary(label) ) {
 447             if(particle->GetPDG()->Charge()>0.) {
 448               fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
 449             }
 450             if(particle->GetPDG()->Charge()<0.) {
 451               fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
 452             }
 453           }
 454         }
 455
 456       }//trackCuts global tracks
 457
 458       if(fTrackType==1 || fTrackType==2)
 459         delete track;
 460     }//track loop
 461
 462
 463   //Fill MC containters if particles are findable
 464   if(fMC) {
 465     for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
 466       AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
 467       if(!mcPart) continue;
 468       //fill the container
 469       containerInputMC[0] = mcPart->Pt();
 470       containerInputMC[1] = mcPart->Phi();
 471       containerInputMC[2] = mcPart->Eta();
 472
 473       if(fStack->IsPhysicalPrimary(iPart)) {
 474         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 475         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 476       }
 477     }
 478   }
 479
 480   PostData(0,fHistList);
 481   PostData(1,fCFManagerPos->GetParticleContainer());
 482   PostData(2,fCFManagerNeg->GetParticleContainer());
 483
 484 }
 485 //________________________________________________________________________
 486 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 487   //
 488   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 489   // This is to called in Notify and should provide the path to the AOD/ESD file
 490   // Copied from AliAnalysisTaskJetSpectrum2
 491   //
 492
 493   TString file(currFile);
 494   fXsec = 0;
 495   fTrials = 1;
 496
 497   if(file.Contains("root_archive.zip#")){
 498     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 499     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 500     file.Replace(pos+1,20,"");
 501   }
 502   else {
 503     // not an archive take the basename....
 504     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 505   }
 506   //  Printf("%s",file.Data());
 507
 508   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
 509   if(!fxsec){
 510     // next trial fetch the histgram file
 511     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 512     if(!fxsec){
 513         // not a severe condition but inciate that we have no information
 514       return kFALSE;
 515     }
 516     else{
 517       // find the tlist we want to be independtent of the name so use the Tkey
 518       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 519       if(!key){
 520         fxsec->Close();
 521         return kFALSE;
 522       }
 523       TList *list = dynamic_cast<TList*>(key->ReadObj());
 524       if(!list){
 525         fxsec->Close();
 526         return kFALSE;
 527       }
 528       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 529       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 530       fxsec->Close();
 531     }
 532   } // no tree pyxsec.root
 533   else {
 534     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 535     if(!xtree){
 536       fxsec->Close();
 537       return kFALSE;
 538     }
 539     UInt_t   ntrials  = 0;
 540     Double_t  xsection  = 0;
 541     xtree->SetBranchAddress("xsection",&xsection);
 542     xtree->SetBranchAddress("ntrials",&ntrials);
 543     xtree->GetEntry(0);
 544     fTrials = ntrials;
 545     fXsec = xsection;
 546     fxsec->Close();
 547   }
 548   return kTRUE;
 549 }
 550 //________________________________________________________________________
 551 Bool_t AliPWG4HighPtSpectra::Notify()
 552 {
 553   //
 554   // Implemented Notify() to read the cross sections
 555   // and number of trials from pyxsec.root
 556   // Copied from AliAnalysisTaskJetSpectrum2
 557   //
 558
 559   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 560   Float_t xsection = 0;
 561   Float_t ftrials  = 1;
 562
 563   fAvgTrials = 1;
 564   if(tree){
 565     TFile *curfile = tree->GetCurrentFile();
 566     if (!curfile) {
 567       Error("Notify","No current file");
 568       return kFALSE;
 569     }
 570     if(!fh1Xsec||!fh1Trials){
 571       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 572       return kFALSE;
 573     }
 574      PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 575     fh1Xsec->Fill("<#sigma>",xsection);
 576     // construct a poor man average trials
 577     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 578     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 579   }
 580   return kTRUE;
 581 }
 582
 583 //________________________________________________________________________
 584 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
 585
 586   if(!mcEvent)return 0;
 587   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 588   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 589   if(!pythiaGenHeader){
 590     // cocktail ??
 591     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 592
 593     if (!genCocktailHeader) {
 594       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 595       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 596       return 0;
 597     }
 598     TList* headerList = genCocktailHeader->GetHeaders();
 599     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 600       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 601       if (pythiaGenHeader)
 602         break;
 603     }
 604     if(!pythiaGenHeader){
 605       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 606       return 0;
 607     }
 608   }
 609   return pythiaGenHeader;
 610
 611 }
 612
 613
 614 //___________________________________________________________________________
 615 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 616 {
 617   // The Terminate() function is the last function to be called during
 618   // a query. It always runs on the client, it can be used to present
 619   // the results graphically or save the results to file.
 620
 621 }
 622
 623 //___________________________________________________________________________
 624 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 625   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 626   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 627   //
 628   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 629
 630   Bool_t oldStatus = TH1::AddDirectoryStatus();
 631   TH1::AddDirectory(kFALSE);
 632
 633   //slot #1
 634   OpenFile(0);
 635   fHistList = new TList();
 636   fHistList->SetOwner(kTRUE);
 637   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 638   fHistList->Add(fNEventAll);
 639   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 640   fHistList->Add(fNEventSel);
 641
 642   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 643   //Set labels
 644   fNEventReject->Fill("noESD",0);
 645   fNEventReject->Fill("Trigger",0);
 646   fNEventReject->Fill("NTracks<2",0);
 647   fNEventReject->Fill("noVTX",0);
 648   fNEventReject->Fill("VtxStatus",0);
 649   fNEventReject->Fill("NCont<2",0);
 650   fNEventReject->Fill("ZVTX>10",0);
 651   fNEventReject->Fill("cent",0);
 652   fNEventReject->Fill("cent>90",0);
 653   fHistList->Add(fNEventReject);
 654
 655   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
 656   fHistList->Add(fh1Centrality);
 657
 658   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 659   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 660   fHistList->Add(fh1Xsec);
 661
 662   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 663   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 664   fHistList->Add(fh1Trials);
 665
 666   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 667   fHistList->Add(fh1PtHard);
 668   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 669   fHistList->Add(fh1PtHardTrials);
 670
 671   TH1::AddDirectory(oldStatus);
 672
 673   PostData(0,fHistList);
 674   PostData(1,fCFManagerPos->GetParticleContainer());
 675   PostData(2,fCFManagerNeg->GetParticleContainer());
 676
 677 }
 678
 679 #endif