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         continue;
 427       else {
 428         track = static_cast<AliVTrack*>(aodtrack);
 429       }
 430       //fill the container
 431       containerInputRec[0] = track->Pt();
 432       containerInputRec[1] = track->Phi();
 433       containerInputRec[2] = track->Eta();
 434       containerInputRec[3] = track->GetTPCNcls();
 435
 436       if(track->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 437       if(track->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 438
 439       if(fArrayMCAOD) {
 440         Int_t label = TMath::Abs(track->GetLabel());
 441         if(label>fArrayMCAOD->GetEntries()) {
 442           continue;
 443         }
 444         AliAODMCParticle *particle = (AliAODMCParticle*) fArrayMCAOD->At(label);
 445         if(!particle) {
 446           continue;
 447         }
 448         //Only select particles generated by HIJING if requested
 449         if(fbSelectHIJING) {
 450           if(!IsHIJINGParticle(label)) {
 451             continue;
 452           }
 453         }
 454         containerInputRecMC[0] = particle->Pt();
 455         containerInputRecMC[1] = particle->Phi();
 456         containerInputRecMC[2] = particle->Eta();
 457         containerInputRecMC[3] = track->GetTPCNcls();
 458
 459         //Container with primaries
 460         if(particle->IsPhysicalPrimary()) {
 461           if(particle->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepReconstructedMC,particle)) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 462           if(particle->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepReconstructedMC,particle)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 463           //Fill pT resolution plots for primaries
 464           //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
 465         }
 466
 467         //Container with secondaries
 468         if (!particle->IsPhysicalPrimary() ) {
 469           if(particle->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
 470           if(particle->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepSecondaries);
 471           //Fill pT resolution plots for primaries
 472           //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
 473         }
 474       }
 475     }//track loop
 476
 477     //Fill MC containers if particles are findable
 478     if(fArrayMCAOD) {
 479       int noPart = fArrayMCAOD->GetEntriesFast();
 480
 481       for(int iPart = 1; iPart<noPart; iPart++) {
 482         AliAODMCParticle *mcPart = (AliAODMCParticle*) fArrayMCAOD->At(iPart);
 483         if(!mcPart) continue;
 484
 485         //Only select particles generated by HIJING if requested
 486         if(fbSelectHIJING) {
 487           if(!IsHIJINGParticle(iPart))
 488             continue;
 489         }
 490
 491         Int_t pdg = mcPart->PdgCode();
 492
 493         // select charged pions, protons, kaons , electrons, muons
 494         if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
 495            321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
 496
 497           //fill the container
 498           containerInputMC[0] = mcPart->Pt();
 499           containerInputMC[1] = mcPart->Phi();
 500           containerInputMC[2] = mcPart->Eta();
 501           //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
 502           containerInputMC[3] = 159.;
 503
 504           if(mcPart->IsPhysicalPrimary()) {
 505             if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 506             if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 507           }
 508         }
 509       }
 510     }
 511
 512   }//end of AOD analysis
 513   else {//ESD analysis
 514     for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 515     {
 516       //Get track for analysis
 517       AliESDtrack *track = 0x0;
 518       AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
 519       if(!esdtrack)
 520         continue;
 521
 522
 523       if(fTrackType==4) {
 524         if (!(fTrackCuts->AcceptTrack(esdtrack)))
 525           continue;
 526       }
 527
 528       if(fTrackType==1)
 529         track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 530       else if(fTrackType==2 || fTrackType==4) {
 531         track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());
 532         if(!track)
 533           continue;
 534
 535         AliExternalTrackParam exParam;
 536         Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
 537         if( !relate ) {
 538           if(track) delete track;
 539           continue;
 540         }
 541         track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 542       }
 543       else if(fTrackType==5 || fTrackType==6) {
 544         if(fTrackCuts->AcceptTrack(esdtrack)) {
 545           continue;
 546         }
 547         else {
 548           if( !(fTrackCutsReject->AcceptTrack(esdtrack)) && fTrackCuts->AcceptTrack(esdtrack) ) {
 549
 550             if(fTrackType==5) {
 551               //use TPConly constrained track
 552               track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
 553               if(!track)
 554                 continue;
 555
 556               AliExternalTrackParam exParam;
 557               Bool_t relate = track->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam);
 558               if( !relate ) {
 559                 if(track) delete track;
 560                 continue;
 561               }
 562               track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 563             }
 564             else if(fTrackType==6) {
 565               //use global constrained track
 566               track = new AliESDtrack(*esdtrack);
 567               track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 568             }
 569           }
 570         }
 571       }
 572       else if(fTrackType==7) {
 573         //use global constrained track
 574         track = new AliESDtrack(*esdtrack);
 575       }
 576       else
 577         track = esdtrack;
 578
 579       if(!track)
 580         continue;
 581
 582
 583       if(fTrackType==2 || fTrackType==4 || fTrackType==5) {
 584         //Cut on chi2 of constrained fit
 585         if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax && fSigmaConstrainedMax>0.) {
 586           if(track) delete track;
 587           continue;
 588         }
 589       }
 590       if (!(fTrackCuts->AcceptTrack(track)) && fTrackType!=4 && fTrackType!=5 && fTrackType!=6) {
 591         if(fTrackType==1 || fTrackType==2 || fTrackType==7) {
 592           if(track) delete track;
 593         }
 594         continue;
 595       }
 596
 597       if(fTrackType==7) {
 598         if(fTrackCutsReject ) {
 599           if(fTrackCutsReject->AcceptTrack(track) ) {
 600             if(track) delete track;
 601             continue;
 602           }
 603         }
 604         if(esdtrack->GetConstrainedParam())
 605           track->Set(esdtrack->GetConstrainedParam()->GetX(),esdtrack->GetConstrainedParam()->GetAlpha(),esdtrack->GetConstrainedParam()->GetParameter(),esdtrack->GetConstrainedParam()->GetCovariance());
 606       }
 607
 608       if(!track) {
 609         if(fTrackType==1 || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
 610           if(track) delete track;
 611         }
 612         continue;
 613       }
 614
 615       //fill the container
 616       containerInputRec[0] = track->Pt();
 617       containerInputRec[1] = track->Phi();
 618       containerInputRec[2] = track->Eta();
 619       containerInputRec[3] = track->GetTPCNclsIter1();
 620
 621       if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 622       if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 623
 624
 625       //Only fill the MC containers if MC information is available
 626       if(fMC) {
 627         Int_t label = TMath::Abs(track->GetLabel());
 628         if(label>fStack->GetNtrack()) {
 629           if(fTrackType==1 || fTrackType==2 || fTrackType==7)
 630             delete track;
 631           continue;
 632         }
 633         //Only select particles generated by HIJING if requested
 634         if(fbSelectHIJING) {
 635           if(!IsHIJINGParticle(label)) {
 636             if(fTrackType==1 || fTrackType==2 || fTrackType==7)
 637               delete track;
 638             continue;
 639           }
 640         }
 641         TParticle *particle = fStack->Particle(label) ;
 642         if(!particle) {
 643           if(fTrackType==1 || fTrackType==2 || fTrackType==7)
 644             delete track;
 645           continue;
 646         }
 647         containerInputRecMC[0] = particle->Pt();
 648         containerInputRecMC[1] = particle->Phi();
 649         containerInputRecMC[2] = particle->Eta();
 650         containerInputRecMC[3] = track->GetTPCNclsIter1();
 651
 652         //Container with primaries
 653         if(fStack->IsPhysicalPrimary(label)) {
 654           if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 655           if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 656           //Fill pT resolution plots for primaries
 657           fPtRelUncertainty1PtPrim->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
 658         }
 659
 660         //Container with secondaries
 661         if (!fStack->IsPhysicalPrimary(label) ) {
 662           if(particle->GetPDG()->Charge()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
 663           if(particle->GetPDG()->Charge()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
 664           //Fill pT resolution plots for primaries
 665           fPtRelUncertainty1PtSec->Fill(containerInputRec[0],containerInputRec[0]*TMath::Sqrt(track->GetSigma1Pt2()));
 666         }
 667       }
 668
 669       if(fTrackType==1  || fTrackType==2 || fTrackType==4 || fTrackType==5 || fTrackType==6 || fTrackType==7) {
 670         if(track) delete track;
 671       }
 672     }//track loop
 673
 674     //Fill MC containers if particles are findable
 675     if(fMC) {
 676       for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
 677         AliMCParticle *mcPart  = (AliMCParticle*)fMC->GetTrack(iPart);
 678         if(!mcPart) continue;
 679
 680         //Only select particles generated by HIJING if requested
 681         if(fbSelectHIJING) {
 682           if(!IsHIJINGParticle(iPart))
 683             continue;
 684         }
 685
 686         Int_t pdg = mcPart->PdgCode();
 687
 688         // select charged pions, protons, kaons , electrons, muons
 689         if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) ==
 690            321 || TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
 691
 692           //fill the container
 693           containerInputMC[0] = mcPart->Pt();
 694           containerInputMC[1] = mcPart->Phi();
 695           containerInputMC[2] = mcPart->Eta();
 696           //      AliESDtrack *esdtrack = fESD->GetTrack(mcPart->GetLabel());
 697           containerInputMC[3] = 159.;
 698
 699           if(fStack->IsPhysicalPrimary(iPart)) {
 700             if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 701             if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 702           }
 703         }
 704       }
 705     }
 706   }//end of ESD analysis
 707
 708   PostData(0,fHistList);
 709   PostData(1,fCFManagerPos->GetParticleContainer());
 710   PostData(2,fCFManagerNeg->GetParticleContainer());
 711
 712 }
 713
 714 //________________________________________________________________________
 715 Bool_t AliPWG4HighPtSpectra::PythiaInfoFromFile(const char* currFile,Float_t &fXsec,Float_t &fTrials)
 716 {
 717   //
 718   // get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 719   // This is to called in Notify and should provide the path to the AOD/ESD file
 720   // Copied from AliAnalysisTaskJetSpectrum2
 721   //
 722
 723   TString file(currFile);
 724   fXsec = 0;
 725   fTrials = 1;
 726
 727   if(file.Contains("root_archive.zip#")){
 728     Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 729     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 730     file.Replace(pos+1,20,"");
 731   }
 732   else {
 733     // not an archive take the basename....
 734     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 735   }
 736   //  Printf("%s",file.Data());
 737
 738   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
 739   if(!fxsec){
 740     // next trial fetch the histgram file
 741     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 742     if(!fxsec){
 743       // not a severe condition but indicates that we have no information
 744       return kFALSE;
 745     }
 746     else{
 747       // find the tlist we want to be independtent of the name so use the Tkey
 748       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 749       if(!key){
 750         fxsec->Close();
 751         return kFALSE;
 752       }
 753       TList *list = dynamic_cast<TList*>(key->ReadObj());
 754       if(!list){
 755         fxsec->Close();
 756         return kFALSE;
 757       }
 758       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 759       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 760       fxsec->Close();
 761     }
 762   } // no tree pyxsec.root
 763   else {
 764     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 765     if(!xtree){
 766       fxsec->Close();
 767       return kFALSE;
 768     }
 769     UInt_t   ntrials  = 0;
 770     Double_t  xsection  = 0;
 771     xtree->SetBranchAddress("xsection",&xsection);
 772     xtree->SetBranchAddress("ntrials",&ntrials);
 773     xtree->GetEntry(0);
 774     fTrials = ntrials;
 775     fXsec = xsection;
 776     fxsec->Close();
 777   }
 778   return kTRUE;
 779 }
 780
 781 //________________________________________________________________________
 782 Bool_t AliPWG4HighPtSpectra::Notify()
 783 {
 784   //
 785   // Implemented Notify() to read the cross sections
 786   // and number of trials from pyxsec.root
 787   // Copied from AliAnalysisTaskJetSpectrum2
 788   //
 789
 790   if(fNoPythiaInfo)
 791     return kTRUE;
 792
 793   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 794   Float_t xsection = 0;
 795   Float_t ftrials  = 1;
 796
 797   if(tree){
 798     TFile *curfile = tree->GetCurrentFile();
 799     if (!curfile) {
 800       Error("Notify","No current file");
 801       return kFALSE;
 802     }
 803     if(!fh1Xsec||!fh1Trials){
 804       //      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 805       return kFALSE;
 806     }
 807     PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 808     fh1Xsec->Fill("<#sigma>",xsection);
 809     fh1Trials->Fill("#sum{ntrials}",ftrials);
 810   }
 811   return kTRUE;
 812 }
 813
 814 //________________________________________________________________________
 815 AliGenPythiaEventHeader*  AliPWG4HighPtSpectra::GetPythiaEventHeader(AliMCEvent *mcEvent)
 816 {
 817
 818   if(!mcEvent)return 0;
 819   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 820   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 821   if(!pythiaGenHeader){
 822     // cocktail ??
 823     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 824
 825     if (!genCocktailHeader) {
 826       //      AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 827       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 828       return 0;
 829     }
 830     TList* headerList = genCocktailHeader->GetHeaders();
 831     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 832       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 833       if (pythiaGenHeader)
 834         break;
 835     }
 836     if(!pythiaGenHeader){
 837       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Pythia event header not found");
 838       return 0;
 839     }
 840   }
 841   return pythiaGenHeader;
 842
 843 }
 844
 845 //________________________________________________________________________
 846 AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeader(AliMCEvent *mcEvent)
 847 {
 848
 849   if(!mcEvent)return 0;
 850   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
 851   AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 852   if(!hijingGenHeader){
 853     // cocktail ??
 854     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 855
 856     if (!genCocktailHeader) {
 857       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Unknown header type (not Pythia or Cocktail)");
 858       //      AliWarning(Form("%s %d: Unknown header type (not Pythia or Cocktail)",(char*)__FILE__,__LINE__));
 859       return 0;
 860     }
 861     TList* headerList = genCocktailHeader->GetHeaders();
 862     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 863       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
 864       if (hijingGenHeader)
 865         break;
 866     }
 867     if(!hijingGenHeader){
 868       AliWarningGeneral(Form(" %s:%d",(char*)__FILE__,__LINE__),"Hijing event header not found");
 869       return 0;
 870     }
 871   }
 872   return hijingGenHeader;
 873
 874 }
 875
 876 //________________________________________________________________________
 877 AliGenHijingEventHeader*  AliPWG4HighPtSpectra::GetHijingEventHeaderAOD(AliAODEvent *aodEvent)
 878 {
 879   AliGenHijingEventHeader* hijingGenHeader = 0x0;
 880   if(aodEvent) {
 881     AliAODMCHeader* mcHeader = (AliAODMCHeader*) aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 882     TList* headerList = mcHeader->GetCocktailHeaders();
 883     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 884       hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
 885       if (hijingGenHeader)
 886         break;
 887     }
 888   }
 889   return hijingGenHeader;
 890 }
 891
 892 //___________________________________________________________________________
 893 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 894 {
 895   // The Terminate() function is the last function to be called during
 896   // a query. It always runs on the client, it can be used to present
 897   // the results graphically or save the results to file.
 898
 899 }
 900
 901 //___________________________________________________________________________
 902 void AliPWG4HighPtSpectra::CreateOutputObjects()
 903 {
 904   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 905   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 906   //
 907   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 908
 909   Bool_t oldStatus = TH1::AddDirectoryStatus();
 910   TH1::AddDirectory(kFALSE);
 911
 912   //slot #1
 913   OpenFile(0);
 914   fHistList = new TList();
 915   fHistList->SetOwner(kTRUE);
 916   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 917   fHistList->Add(fNEventAll);
 918   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 919   fHistList->Add(fNEventSel);
 920
 921   fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
 922   //Set labels
 923   fNEventReject->Fill("noESD/AOD",0);
 924   fNEventReject->Fill("Trigger",0);
 925   fNEventReject->Fill("NTracks<2",0);
 926   fNEventReject->Fill("noVTX",0);
 927   fNEventReject->Fill("VtxStatus",0);
 928   fNEventReject->Fill("NCont<2",0);
 929   fNEventReject->Fill("ZVTX>10",0);
 930   fNEventReject->Fill("cent",0);
 931   fNEventReject->Fill("cent>90",0);
 932   fHistList->Add(fNEventReject);
 933
 934   fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
 935   fHistList->Add(fh1Centrality);
 936
 937   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 938   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 939   fHistList->Add(fh1Xsec);
 940
 941   fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 942   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 943   fHistList->Add(fh1Trials);
 944
 945   fh1PtHard       = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
 946   fHistList->Add(fh1PtHard);
 947   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
 948   fHistList->Add(fh1PtHardTrials);
 949
 950   Int_t fgkNPtBins = 100;
 951   Float_t kMinPt   = 0.;
 952   Float_t kMaxPt   = 100.;
 953   Double_t *binsPt = new Double_t[fgkNPtBins+1];
 954   for(Int_t i=0; i<=fgkNPtBins; i++) binsPt[i]=(Double_t)kMinPt + (kMaxPt-kMinPt)/fgkNPtBins*(Double_t)i ;
 955
 956   Int_t fgkNRel1PtUncertaintyBins=50;
 957   Float_t fgkRel1PtUncertaintyMin = 0.;
 958   Float_t fgkRel1PtUncertaintyMax = 1.;
 959
 960   Double_t *binsRel1PtUncertainty=new Double_t[fgkNRel1PtUncertaintyBins+1];
 961   for(Int_t i=0; i<=fgkNRel1PtUncertaintyBins; i++) binsRel1PtUncertainty[i]=(Double_t)fgkRel1PtUncertaintyMin + (fgkRel1PtUncertaintyMax-fgkRel1PtUncertaintyMin)/fgkNRel1PtUncertaintyBins*(Double_t)i ;
 962
 963   fPtRelUncertainty1PtPrim = new TH2F("fPtRelUncertainty1PtPrim","fPtRelUncertainty1PtPrim",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 964   fHistList->Add(fPtRelUncertainty1PtPrim);
 965
 966   fPtRelUncertainty1PtSec = new TH2F("fPtRelUncertainty1PtSec","fPtRelUncertainty1PtSec",fgkNPtBins,binsPt,fgkNRel1PtUncertaintyBins,binsRel1PtUncertainty);
 967   fHistList->Add(fPtRelUncertainty1PtSec);
 968
 969   TH1::AddDirectory(oldStatus);
 970
 971   OpenFile(1);
 972   OpenFile(2);
 973
 974   PostData(0,fHistList);
 975   PostData(1,fCFManagerPos->GetParticleContainer());
 976   PostData(2,fCFManagerNeg->GetParticleContainer());
 977
 978   if(binsPt)                delete [] binsPt;
 979   if(binsRel1PtUncertainty) delete [] binsRel1PtUncertainty;
 980
 981 }
 982
 983 //________________________________________________________________________
 984 Bool_t AliPWG4HighPtSpectra::IsHIJINGParticle(Int_t label)
 985 {
 986   //
 987   // Return kTRUE in case particle is from HIJING event
 988   //
 989
 990   AliGenHijingEventHeader* hijingHeader;
 991   if(fReadAODData)
 992     hijingHeader = GetHijingEventHeaderAOD(fAOD);
 993   else
 994     hijingHeader = GetHijingEventHeader(fMC);
 995
 996   Int_t nproduced = hijingHeader->NProduced();
 997
 998   if(fReadAODData) {
 999     AliAODMCParticle* mom = (AliAODMCParticle*) fArrayMCAOD->At(label);
1000     Int_t iMom = label;
1001     Int_t iParent = mom->GetMother();
1002     while(iParent!=-1){
1003       iMom = iParent;
1004       mom = (AliAODMCParticle*) fArrayMCAOD->At(iMom);
1005       iParent = mom->GetMother();
1006     }
1007
1008     if(iMom<nproduced) {
1009       return kTRUE;
1010     } else {
1011       return kFALSE;
1012     }
1013   }
1014   else { //ESD analysis
1015     TParticle * mom = fStack->Particle(label);
1016     Int_t iMom = label;
1017     Int_t iParent = mom->GetFirstMother();
1018     while(iParent!=-1){
1019       iMom = iParent;
1020       mom = fStack->Particle(iMom);
1021       iParent = mom->GetFirstMother();
1022     }
1023
1024     if(iMom<nproduced) {
1025       return kTRUE;
1026     } else {
1027       return kFALSE;
1028     }
1029   }
1030 }
1031
1032 #endif