PWG4/JetTasks/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 // Example of task (running locally, on AliEn and CAF),
  18 // which provides standard way of calculating acceptance and efficiency
  19 // between different steps of the procedure.
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies
  21 // can be calculated
  22 //-----------------------------------------------------------------------
  23 // Author : Marta Verweij - UU
  24 //-----------------------------------------------------------------------
  25
  26
  27 #ifndef ALIPWG4HIGHPTSPECTRA_CXX
  28 #define ALIPWG4HIGHPTSPECTRA_CXX
  29
  30 #include "AliPWG4HighPtSpectra.h"
  31
  32 #include "TVector3.h"
  33 #include <iostream>
  34 #include "TH1.h"
  35 #include "TH2.h"
  36 #include "TH3.h"
  37 #include "TList.h"
  38 #include "TChain.h"
  39
  40 #include "AliAnalysisManager.h"
  41 #include "AliESDInputHandler.h"
  42 #include "AliESDtrack.h"
  43 #include "AliESDtrackCuts.h"
  44 #include "AliExternalTrackParam.h"
  45
  46 #include "AliLog.h"
  47
  48 #include "AliStack.h"
  49 #include "TParticle.h"
  50 #include "AliMCEvent.h"
  51 #include "AliMCEventHandler.h"
  52 #include "AliCFContainer.h"
  53
  54 //#include "$ALICE_ROOT/PWG4/JetTasks/AliAnalysisHelperJetTasks.h"
  55
  56 //#include <iostream>
  57 using namespace std; //required for resolving the 'cout' symbol
  58
  59 ClassImp(AliPWG4HighPtSpectra)
  60
  61 //__________________________________________________________________________
  62 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra() : AliAnalysisTask("AliPWG4HighPtSpectra", ""),
  63   fReadAODData(0),
  64   fCFManagerPos(0x0),
  65   fCFManagerNeg(0x0),
  66   fESD(0),
  67   fTrackCuts(0),
  68   fTrackCutsTPConly(0),
  69   fHistList(0),
  70   fNEventAll(0),
  71   fNEventSel(0)
  72 {
  73   //
  74   //Default ctor
  75   //
  76 }
  77 //___________________________________________________________________________
  78 AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const Char_t* name) :
  79   AliAnalysisTask(name,""),
  80   fReadAODData(0),
  81   fCFManagerPos(0x0),
  82   fCFManagerNeg(0x0),
  83   fESD(0),
  84   fTrackCuts(),
  85   fTrackCutsTPConly(0),
  86   fHistList(0),
  87   fNEventAll(0),
  88   fNEventSel(0)
  89 {
  90   //
  91   // Constructor. Initialization of Inputs and Outputs
  92   //
  93   AliDebug(2,Form("AliPWG4HighPtSpectra Calling Constructor"));
  94   // Input slot #0 works with a TChain ESD
  95   DefineInput(0, TChain::Class());
  96   // Output slot #0 writes into a TList
  97   DefineOutput(0,TList::Class());
  98   // Output slot #1, #2 writes into a AliCFContainer
  99   DefineOutput(1,AliCFContainer::Class());
 100   DefineOutput(2,AliCFContainer::Class());
 101   // Output slot #3 writes into a AliESDtrackCuts
 102   DefineOutput(3, AliESDtrackCuts::Class());
 103   DefineOutput(4, AliESDtrackCuts::Class());
 104 }
 105
 106 //________________________________________________________________________
 107 void AliPWG4HighPtSpectra::LocalInit()
 108 {
 109   //
 110   // Only called once at beginning
 111   //
 112   PostData(3,fTrackCuts);
 113   PostData(4,fTrackCutsTPConly);
 114 }
 115
 116 //________________________________________________________________________
 117 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
 118 {
 119   // Connect ESD here
 120   // Called once
 121   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 122   //  cout << "cout >> AliPWG4HighPtSpectra::ConnectInputData" << endl;
 123   printf(">> AliPWG4HighPtSpectra::ConnectInputData \n");
 124
 125   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 126   if (!tree) {
 127     AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
 128   } else {
 129
 130     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 131
 132     if (!esdH) {
 133       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 134     } else {
 135       fESD = esdH->GetEvent();
 136     }
 137   }
 138
 139 }
 140 //_________________________________________________
 141 void AliPWG4HighPtSpectra::Exec(Option_t *)
 142 {
 143   //
 144   // Main loop function
 145   //
 146   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
 147
 148   // All events without selection
 149   fNEventAll->Fill(0.);
 150
 151   if (!fESD) {
 152     AliDebug(2,Form("ERROR: fESD not available"));
 153     PostData(0,fHistList);
 154     PostData(1,fCFManagerPos->GetParticleContainer());
 155     PostData(2,fCFManagerNeg->GetParticleContainer());
 156     return;
 157   }
 158
 159   UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 160   if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
 161     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 162     PostData(0,fHistList);
 163     PostData(1,fCFManagerPos->GetParticleContainer());
 164     PostData(2,fCFManagerNeg->GetParticleContainer());
 165     return;
 166   }
 167
 168   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 169   // This handler can return the current MC event
 170
 171   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 172
 173   AliStack* stack = 0x0;
 174   AliMCEvent* mcEvent = 0x0;
 175
 176   if(eventHandler) {
 177     mcEvent = eventHandler->MCEvent();
 178     if (!mcEvent) {
 179       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
 180       PostData(0,fHistList);
 181       PostData(1,fCFManagerPos->GetParticleContainer());
 182       PostData(2,fCFManagerNeg->GetParticleContainer());
 183       return;
 184     }
 185
 186     AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
 187
 188     stack = mcEvent->Stack();                //Particles Stack
 189
 190     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
 191   }
 192
 193   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 194   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 195   // Need vertex cut
 196   TString vtxName(vtx->GetName());
 197   if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 198     // SPD vertex
 199     vtx = fESD->GetPrimaryVertexSPD();
 200     if(vtx->GetNContributors()<2) {
 201       vtx = 0x0;
 202       // Post output data
 203       PostData(0,fHistList);
 204       PostData(1,fCFManagerPos->GetParticleContainer());
 205       PostData(2,fCFManagerNeg->GetParticleContainer());
 206       return;
 207     }
 208   }
 209
 210   double primVtx[3];
 211   vtx->GetXYZ(primVtx);
 212   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 213     PostData(0,fHistList);
 214     PostData(1,fCFManagerPos->GetParticleContainer());
 215     PostData(2,fCFManagerNeg->GetParticleContainer());
 216     return;
 217   }
 218
 219   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2){
 220     // Post output data
 221     PostData(0,fHistList);
 222     PostData(1,fCFManagerPos->GetParticleContainer());
 223     PostData(2,fCFManagerNeg->GetParticleContainer());
 224     return;
 225   }
 226   Int_t nTracks = fESD->GetNumberOfTracks();
 227   AliDebug(2,Form("nTracks %d", nTracks));
 228
 229   if(!fTrackCuts) {
 230     // Post output data
 231     PostData(0,fHistList);
 232     PostData(1,fCFManagerPos->GetParticleContainer());
 233     PostData(2,fCFManagerNeg->GetParticleContainer());
 234     return;
 235   }
 236
 237   // Selected events for analysis
 238   fNEventSel->Fill(0.);
 239
 240
 241   Double_t containerInputRec[3]       = {0.,0.,0.};
 242   Double_t containerInputTPConly[3]   = {0.,0.,0.};
 243   Double_t containerInputMC[3]        = {0.,0.,0.};
 244   Double_t containerInputRecMC[3]     = {0.,0.,0.};
 245   Double_t containerInputTPConlyMC[3] = {0.,0.,0.};
 246
 247   //Now go to rec level
 248   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 249     {
 250       if(!fESD->GetTrack(iTrack) ) continue;
 251       AliESDtrack* track = fESD->GetTrack(iTrack);
 252       if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
 253       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
 254       if(!track || !trackTPC) continue;
 255
 256       //fill the container
 257       containerInputRec[0] = track->Pt();
 258       containerInputRec[1] = track->Phi();
 259       containerInputRec[2] = track->Eta();
 260
 261       //Store TPC Inner Params for TPConly tracks
 262       containerInputTPConly[0] = trackTPC->Pt();
 263       containerInputTPConly[1] = trackTPC->Phi();
 264       containerInputTPConly[2] = trackTPC->Eta();
 265
 266       AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
 267       if(trackTPCESD) {
 268         if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
 269           if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 270           if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 271
 272           //Only fill the MC containers if MC information is available
 273           if(eventHandler) {
 274             Int_t label = TMath::Abs(track->GetLabel());
 275             TParticle *particle = stack->Particle(label) ;
 276             if(!particle) continue;
 277
 278             containerInputTPConlyMC[0] = particle->Pt();
 279             containerInputTPConlyMC[1] = particle->Phi();
 280             containerInputTPConlyMC[2] = particle->Eta();
 281
 282             //Container with primaries
 283             if(stack->IsPhysicalPrimary(label)) {
 284               if(particle->GetPDG()->Charge()>0.) {
 285                 fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedTPCOnlyMC);
 286               }
 287               if(particle->GetPDG()->Charge()<0.) {
 288                 fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedTPCOnlyMC);
 289               }
 290             }
 291           }
 292
 293         }
 294       }
 295
 296       if (fTrackCuts->AcceptTrack(track)) {
 297         if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 298         if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 299
 300
 301         //Only fill the MC containers if MC information is available
 302         if(eventHandler) {
 303           Int_t label = TMath::Abs(track->GetLabel());
 304           TParticle *particle = stack->Particle(label) ;
 305           if(!particle) continue;
 306
 307           containerInputRecMC[0] = particle->Pt();
 308           containerInputRecMC[1] = particle->Phi();
 309           containerInputRecMC[2] = particle->Eta();
 310
 311           //Container with primaries
 312           if(stack->IsPhysicalPrimary(label)) {
 313             if(particle->GetPDG()->Charge()>0.) {
 314               fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 315             }
 316             if(particle->GetPDG()->Charge()<0.) {
 317               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
 318             }
 319           }
 320
 321           //Container with secondaries
 322           if (!stack->IsPhysicalPrimary(label) ) {
 323             if(particle->GetPDG()->Charge()>0.) {
 324               fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 325             }
 326             if(particle->GetPDG()->Charge()<0.) {
 327               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 328             }
 329           }
 330         }
 331
 332       }//trackCuts
 333
 334       delete trackTPCESD;
 335     }//track loop
 336
 337
 338   //Fill MC containters if particles are findable
 339   if(eventHandler) {
 340     for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
 341       {
 342         AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
 343         if(!mcPart) continue;
 344         //fill the container
 345         containerInputMC[0] = mcPart->Pt();
 346         containerInputMC[1] = mcPart->Phi();
 347         containerInputMC[2] = mcPart->Eta();
 348
 349         if(stack->IsPhysicalPrimary(iPart)) {
 350           if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 351           if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 352         }
 353       }
 354   }
 355
 356   PostData(0,fHistList);
 357   PostData(1,fCFManagerPos->GetParticleContainer());
 358   PostData(2,fCFManagerNeg->GetParticleContainer());
 359
 360 }
 361
 362
 363 //___________________________________________________________________________
 364 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 365 {
 366   // The Terminate() function is the last function to be called during
 367   // a query. It always runs on the client, it can be used to present
 368   // the results graphically or save the results to file.
 369
 370 }
 371
 372 //___________________________________________________________________________
 373 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 374   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 375   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 376   //
 377   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 378
 379   Bool_t oldStatus = TH1::AddDirectoryStatus();
 380   TH1::AddDirectory(kFALSE);
 381
 382   //slot #1
 383   OpenFile(0);
 384   fHistList = new TList();
 385   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 386   fHistList->Add(fNEventAll);
 387   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 388   fHistList->Add(fNEventSel);
 389
 390   TH1::AddDirectory(oldStatus);
 391
 392 }
 393
 394 #endif