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