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