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