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   // ---- Get MC Header information (for MC productions in pThard bins) ----
 327   Double_t ptHard = 0.;
 328   Double_t nTrials = 1; // trials for MC trigger weight for real data
 329
 330   if(fMC){
 331     AliGenPythiaEventHeader*  pythiaGenHeader = GetPythiaEventHeader(fMC);
 332      if(pythiaGenHeader){
 333        nTrials = pythiaGenHeader->Trials();
 334        ptHard  = pythiaGenHeader->GetPtHard();
 335
 336        fh1PtHard->Fill(ptHard);
 337        fh1PtHardTrials->Fill(ptHard,nTrials);
 338
 339        fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
 340      }
 341   }
 342
 343   Int_t nTracks = fESD->GetNumberOfTracks();
 344   AliDebug(2,Form("nTracks %d", nTracks));
 345
 346   if(!fTrackCuts) {
 347     fNEventReject->Fill("noTrackCuts",1);
 348     // Post output data
 349     PostData(0,fHistList);
 350     PostData(1,fCFManagerPos->GetParticleContainer());
 351     PostData(2,fCFManagerNeg->GetParticleContainer());
 352     return;
 353   }
 354
 355   // Selected events for analysis
 356   fNEventSel->Fill(0.);
 357
 358   const Int_t nvar = 4;
 359
 360   Double_t containerInputRec[nvar]       = {0.,0.,0.,0.};
 361   Double_t containerInputMC[nvar]        = {0.,0.,0.,0.};
 362   Double_t containerInputRecMC[nvar]     = {0.,0.,0.,0.}; //reconstructed yield as function of MC variable
 363
 364   //Now go to rec level
 365   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 366     {
 367       //Get track for analysis
 368       AliESDtrack *track = 0x0;
 369       AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 370       if(!esdtrack) continue;
 371
 372       if(fTrackType==1) {
 373         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 374         if(!track) continue;
 375       }
 376       else if(fTrackType==2) {
 377         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 378         if(!track) continue;
 379
 380         AliExternalTrackParam exParam;
 381         Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
 382         if( !relate ) {
 383           delete track;
 384           continue;
 385         }
 386         track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 387       }
 388       else if(fTrackType==7) {
 389         //use global constrained track
 390         track = new AliESDtrack(*esdtrack);
 391         //      track = esdtrack;
 392         //      track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 393
 394       }
 395       else
 396         track = esdtrack;
 397
 398       if(!track) continue;
 399
 400       if(fTrackType==2) {
 401         //Cut on chi2 of constrained fit
 402         if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
 403           delete track;
 404           continue;
 405         }
 406       }
 407
 408       if (fTrackCuts->AcceptTrack(track)) {
 409
 410         if(fTrackType==7) {
 411           if(fTrackCutsReject ) {
 412             if(fTrackCutsReject->AcceptTrack(track) ) {
 413               if(track) delete track;
 414               continue;
 415             }
 416           }
 417
 418           if(esdtrack->GetConstrainedParam())
 419             track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 420         }
 421
 422         //fill the container
 423         containerInputRec[0] = track->Pt();
 424         containerInputRec[1] = track->Phi();
 425         containerInputRec[2] = track->Eta();
 426         containerInputRec[3] = track->GetTPCNclsIter1();
 427
 428         if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 429         if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 430
 431
 432         //Only fill the MC containers if MC information is available
 433         if(fMC) {
 434           Int_t label = TMath::Abs(track->GetLabel());
 435           TParticle *particle = fStack->Particle(label) ;
 436           if(!particle) {
 437             if(fTrackType==1 || fTrackType==2 || fTrackType==7)
 438               delete track;
 439             continue;
 440           }
 441           containerInputRecMC[0] = particle->Pt();
 442           containerInputRecMC[1] = particle->Phi();
 443           containerInputRecMC[2] = particle->Eta();
 444           containerInputRecMC[3] = track->GetTPCNclsIter1();
 445
 446           //Container with primaries
 447           if(fStack->IsPhysicalPrimary(label)) {
 448             if(particle->GetPDG()->Charge()>0.) {
 449               fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 450             }
 451             if(particle->GetPDG()->Charge()<0.) {
 452               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 453             }
 454
 455             //Fill pT resolution plots for primaries
 456             fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
 457
 458           }
 459
 460           //Container with secondaries
 461           if (!fStack->IsPhysicalPrimary(label) ) {
 462             if(particle->GetPDG()->Charge()>0.) {
 463               fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
 464             }
 465             if(particle->GetPDG()->Charge()<0.) {
 466               fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
 467             }
 468             //Fill pT resolution plots for primaries
 469             fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
 470           }
 471         }
 472
 473       }//trackCuts global tracks
 474
 475       if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
 476         if(track) delete track;
 477       }
 478     }//track loop
 479
 480
 481   //Fill MC containers if particles are findable
 482   if(fMC) {
 483     for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
 484       AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
 485       if(!mcPart) continue;
 486
 487       Int_t pdg = mcPart->PdgCode();
 488
 489       // select charged pions, protons, kaons , electrons, muons
 490       if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
 491          321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
 492
 493         //fill the container
 494         containerInputMC[0] = mcPart->Pt();
 495         containerInputMC[1] = mcPart->Phi();
 496         containerInputMC[2] = mcPart->Eta();
 497         //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
 498         containerInputMC[3] = 159.;
 499
 500         if(fStack->IsPhysicalPrimary(iPart)) {
 501           if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 502           if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 503         }
 504       }
 505     }
 506   }
 507
 508   PostData(0,fHistList);
 509   PostData(1,fCFManagerPos->GetParticleContainer());
 510   PostData(2,fCFManagerNeg->GetParticleContainer());
 511
 512 }
 513 //________________________________________________________________________
 514 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 515   //
 516   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 517   // This is to called in Notify and should provide the path to the AOD/ESD file
 518   // Copied from AliAnalysisTaskJetSpectrum2
 519   //
 520
 521   TString file(currFile);
 522   fXsec = 0;
 523   fTrials = 1;
 524
 525   if(file.Contains("root_archive.zip#")){
 526     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 527     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 528     file.Replace(pos+1,20,"");
 529   }
 530   else {
 531     // not an archive take the basename....
 532     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 533   }
 534   //  Printf("%s",file.Data());
 535
 536   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
 537   if(!fxsec){
 538     // next trial fetch the histgram file
 539     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 540     if(!fxsec){
 541         // not a severe condition but inciate that we have no information
 542       return kFALSE;
 543     }
 544     else{
 545       // find the tlist we want to be independtent of the name so use the Tkey
 546       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 547       if(!key){
 548         fxsec->Close();
 549         return kFALSE;
 550       }
 551       TList *list = dynamic_cast<TList*>(key->ReadObj());
 552       if(!list){
 553         fxsec->Close();
 554         return kFALSE;
 555       }
 556       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 557       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 558       fxsec->Close();
 559     }
 560   } // no tree pyxsec.root
 561   else {
 562     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 563     if(!xtree){
 564       fxsec->Close();
 565       return kFALSE;
 566     }
 567     UInt_t   ntrials  = 0;
 568     Double_t  xsection  = 0;
 569     xtree->SetBranchAddress("xsection",&xsection);
 570     xtree->SetBranchAddress("ntrials",&ntrials);
 571     xtree->GetEntry(0);
 572     fTrials = ntrials;
 573     fXsec = xsection;
 574     fxsec->Close();
 575   }
 576   return kTRUE;
 577 }
 578 //________________________________________________________________________
 579 Bool_t AliPWG4HighPtSpectra::Notify()
 580 {
 581   //
 582   // Implemented Notify() to read the cross sections
 583   // and number of trials from pyxsec.root
 584   // Copied from AliAnalysisTaskJetSpectrum2
 585   //
 586
 587   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 588   Float_t xsection = 0;
 589   Float_t ftrials  = 1;
 590
 591   fAvgTrials = 1;
 592   if(tree){
 593     TFile *curfile = tree->GetCurrentFile();
 594     if (!curfile) {
 595       Error("Notify","No current file");
 596       return kFALSE;
 597     }
 598     if(!fh1Xsec||!fh1Trials){
 599       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 600       return kFALSE;
 601     }
 602      PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 603     fh1Xsec->Fill("<#sigma>",xsection);
 604     // construct a poor man average trials
 605     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 606     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 607   }
 608   return kTRUE;
 609 }
 610
 611 //________________________________________________________________________
 612 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
 613
 614   if(!mcEvent)return 0;
 615   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 616   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 617   if(!pythiaGenHeader){
 618     // cocktail ??
 619     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 620
 621     if (!genCocktailHeader) {
 622       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 623       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 624       return 0;
 625     }
 626     TList* headerList = genCocktailHeader->GetHeaders();
 627     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 628       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 629       if (pythiaGenHeader)
 630         break;
 631     }
 632     if(!pythiaGenHeader){
 633       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 634       return 0;
 635     }
 636   }
 637   return pythiaGenHeader;
 638
 639 }
 640
 641
 642 //___________________________________________________________________________
 643 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 644 {
 645   // The Terminate() function is the last function to be called during
 646   // a query. It always runs on the client, it can be used to present
 647   // the results graphically or save the results to file.
 648
 649 }
 650
 651 //___________________________________________________________________________
 652 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 653   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 654   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 655   //
 656   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 657
 658   Bool_t oldStatus = TH1::AddDirectoryStatus();
 659   TH1::AddDirectory(kFALSE);
 660
 661   //slot #1
 662   OpenFile(0);
 663   fHistList = new TList();
 664   fHistList->SetOwner(kTRUE);
 665   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 666   fHistList->Add(fNEventAll);
 667   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 668   fHistList->Add(fNEventSel);
 669
 670   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 671   //Set labels
 672   fNEventReject->Fill("noESD",0);
 673   fNEventReject->Fill("Trigger",0);
 674   fNEventReject->Fill("NTracks<2",0);
 675   fNEventReject->Fill("noVTX",0);
 676   fNEventReject->Fill("VtxStatus",0);
 677   fNEventReject->Fill("NCont<2",0);
 678   fNEventReject->Fill("ZVTX>10",0);
 679   fNEventReject->Fill("cent",0);
 680   fNEventReject->Fill("cent>90",0);
 681   fHistList->Add(fNEventReject);
 682
 683   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
 684   fHistList->Add(fh1Centrality);
 685
 686   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 687   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 688   fHistList->Add(fh1Xsec);
 689
 690   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 691   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 692   fHistList->Add(fh1Trials);
 693
 694   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 695   fHistList->Add(fh1PtHard);
 696   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 697   fHistList->Add(fh1PtHardTrials);
 698
 699   Int_t fgkNPtBins = 100;
 700   Float_t kMinPt   = 0.;
 701   Float_t kMaxPt   = 100.;
 702   Double_t *binsPt = new Double_t[fgkNPtBins+1];
 703   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
 704
 705   Int_t fgkNRel1PtUncertaintyBins=50;
 706   Float_t fgkRel1PtUncertaintyMin = 0.;
 707   Float_t fgkRel1PtUncertaintyMax = 1.;
 708
 709   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
 710   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
 711
 712   fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 713   fHistList->Add(fPtRelUncertainty1PtPrim);
 714
 715   fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 716   fHistList->Add(fPtRelUncertainty1PtSec);
 717
 718   TH1::AddDirectory(oldStatus);
 719
 720   PostData(0,fHistList);
 721   PostData(1,fCFManagerPos->GetParticleContainer());
 722   PostData(2,fCFManagerNeg->GetParticleContainer());
 723
 724   if(binsPt)                delete [] binsPt;
 725   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 726
 727 }
 728
 729 #endif