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       //fill the container
 487       containerInputMC[0] = mcPart->Pt();
 488       containerInputMC[1] = mcPart->Phi();
 489       containerInputMC[2] = mcPart->Eta();
 490       //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
 491       containerInputMC[3] = 159.;
 492
 493       if(fStack->IsPhysicalPrimary(iPart)) {
 494         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 495         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 496       }
 497     }
 498   }
 499
 500   PostData(0,fHistList);
 501   PostData(1,fCFManagerPos->GetParticleContainer());
 502   PostData(2,fCFManagerNeg->GetParticleContainer());
 503
 504 }
 505 //________________________________________________________________________
 506 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials){
 507   //
 508   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 509   // This is to called in Notify and should provide the path to the AOD/ESD file
 510   // Copied from AliAnalysisTaskJetSpectrum2
 511   //
 512
 513   TString file(currFile);
 514   fXsec = 0;
 515   fTrials = 1;
 516
 517   if(file.Contains("root_archive.zip#")){
 518     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 519     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 520     file.Replace(pos+1,20,"");
 521   }
 522   else {
 523     // not an archive take the basename....
 524     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 525   }
 526   //  Printf("%s",file.Data());
 527
 528   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
 529   if(!fxsec){
 530     // next trial fetch the histgram file
 531     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 532     if(!fxsec){
 533         // not a severe condition but inciate that we have no information
 534       return kFALSE;
 535     }
 536     else{
 537       // find the tlist we want to be independtent of the name so use the Tkey
 538       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 539       if(!key){
 540         fxsec->Close();
 541         return kFALSE;
 542       }
 543       TList *list = dynamic_cast<TList*>(key->ReadObj());
 544       if(!list){
 545         fxsec->Close();
 546         return kFALSE;
 547       }
 548       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 549       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 550       fxsec->Close();
 551     }
 552   } // no tree pyxsec.root
 553   else {
 554     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 555     if(!xtree){
 556       fxsec->Close();
 557       return kFALSE;
 558     }
 559     UInt_t   ntrials  = 0;
 560     Double_t  xsection  = 0;
 561     xtree->SetBranchAddress("xsection",&xsection);
 562     xtree->SetBranchAddress("ntrials",&ntrials);
 563     xtree->GetEntry(0);
 564     fTrials = ntrials;
 565     fXsec = xsection;
 566     fxsec->Close();
 567   }
 568   return kTRUE;
 569 }
 570 //________________________________________________________________________
 571 Bool_t AliPWG4HighPtSpectra::Notify()
 572 {
 573   //
 574   // Implemented Notify() to read the cross sections
 575   // and number of trials from pyxsec.root
 576   // Copied from AliAnalysisTaskJetSpectrum2
 577   //
 578
 579   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 580   Float_t xsection = 0;
 581   Float_t ftrials  = 1;
 582
 583   fAvgTrials = 1;
 584   if(tree){
 585     TFile *curfile = tree->GetCurrentFile();
 586     if (!curfile) {
 587       Error("Notify","No current file");
 588       return kFALSE;
 589     }
 590     if(!fh1Xsec||!fh1Trials){
 591       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 592       return kFALSE;
 593     }
 594      PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 595     fh1Xsec->Fill("<#sigma>",xsection);
 596     // construct a poor man average trials
 597     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 598     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
 599   }
 600   return kTRUE;
 601 }
 602
 603 //________________________________________________________________________
 604 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent){
 605
 606   if(!mcEvent)return 0;
 607   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 608   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 609   if(!pythiaGenHeader){
 610     // cocktail ??
 611     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 612
 613     if (!genCocktailHeader) {
 614       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 615       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 616       return 0;
 617     }
 618     TList* headerList = genCocktailHeader->GetHeaders();
 619     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 620       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 621       if (pythiaGenHeader)
 622         break;
 623     }
 624     if(!pythiaGenHeader){
 625       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 626       return 0;
 627     }
 628   }
 629   return pythiaGenHeader;
 630
 631 }
 632
 633
 634 //___________________________________________________________________________
 635 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 636 {
 637   // The Terminate() function is the last function to be called during
 638   // a query. It always runs on the client, it can be used to present
 639   // the results graphically or save the results to file.
 640
 641 }
 642
 643 //___________________________________________________________________________
 644 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 645   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 646   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 647   //
 648   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 649
 650   Bool_t oldStatus = TH1::AddDirectoryStatus();
 651   TH1::AddDirectory(kFALSE);
 652
 653   //slot #1
 654   OpenFile(0);
 655   fHistList = new TList();
 656   fHistList->SetOwner(kTRUE);
 657   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 658   fHistList->Add(fNEventAll);
 659   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 660   fHistList->Add(fNEventSel);
 661
 662   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 663   //Set labels
 664   fNEventReject->Fill("noESD",0);
 665   fNEventReject->Fill("Trigger",0);
 666   fNEventReject->Fill("NTracks<2",0);
 667   fNEventReject->Fill("noVTX",0);
 668   fNEventReject->Fill("VtxStatus",0);
 669   fNEventReject->Fill("NCont<2",0);
 670   fNEventReject->Fill("ZVTX>10",0);
 671   fNEventReject->Fill("cent",0);
 672   fNEventReject->Fill("cent>90",0);
 673   fHistList->Add(fNEventReject);
 674
 675   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
 676   fHistList->Add(fh1Centrality);
 677
 678   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 679   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 680   fHistList->Add(fh1Xsec);
 681
 682   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 683   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 684   fHistList->Add(fh1Trials);
 685
 686   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 687   fHistList->Add(fh1PtHard);
 688   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 689   fHistList->Add(fh1PtHardTrials);
 690
 691   Int_t fgkNPtBins = 100;
 692   Float_t kMinPt   = 0.;
 693   Float_t kMaxPt   = 100.;
 694   Double_t *binsPt = new Double_t[fgkNPtBins+1];
 695   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
 696
 697   Int_t fgkNRel1PtUncertaintyBins=50;
 698   Float_t fgkRel1PtUncertaintyMin = 0.;
 699   Float_t fgkRel1PtUncertaintyMax = 1.;
 700
 701   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
 702   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
 703
 704   fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 705   fHistList->Add(fPtRelUncertainty1PtPrim);
 706
 707   fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 708   fHistList->Add(fPtRelUncertainty1PtSec);
 709
 710   TH1::AddDirectory(oldStatus);
 711
 712   PostData(0,fHistList);
 713   PostData(1,fCFManagerPos->GetParticleContainer());
 714   PostData(2,fCFManagerNeg->GetParticleContainer());
 715
 716   if(binsPt)                delete [] binsPt;
 717   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 718
 719 }
 720
 721 #endif