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