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