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