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   //  AliMCEventHandler* eventHandler = (AliMCEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 173
 174   AliStack* stack = 0x0;
 175   AliMCEvent* mcEvent = 0x0;
 176
 177   if(eventHandler) {
 178     mcEvent = eventHandler->MCEvent();
 179     if (!mcEvent) {
 180       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
 181       PostData(0,fHistList);
 182       PostData(1,fCFManagerPos->GetParticleContainer());
 183       PostData(2,fCFManagerNeg->GetParticleContainer());
 184     return;
 185     }
 186
 187     AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
 188
 189     stack = mcEvent->Stack();                //Particles Stack
 190
 191     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
 192   }
 193
 194   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 195   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 196   // Need vertex cut
 197   TString vtxName(vtx->GetName());
 198   if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
 199     // SPD vertex
 200     vtx = fESD->GetPrimaryVertexSPD();
 201     if(vtx->GetNContributors()<2) {
 202       vtx = 0x0;
 203       // Post output data
 204       PostData(0,fHistList);
 205       PostData(1,fCFManagerPos->GetParticleContainer());
 206       PostData(2,fCFManagerNeg->GetParticleContainer());
 207       return;
 208     }
 209   }
 210
 211   double primVtx[3];
 212   vtx->GetXYZ(primVtx);
 213   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 214     PostData(0,fHistList);
 215     PostData(1,fCFManagerPos->GetParticleContainer());
 216     PostData(2,fCFManagerNeg->GetParticleContainer());
 217     return;
 218   }
 219
 220   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2){
 221     // Post output data
 222     PostData(0,fHistList);
 223     PostData(1,fCFManagerPos->GetParticleContainer());
 224     PostData(2,fCFManagerNeg->GetParticleContainer());
 225     return;
 226   }
 227   Int_t nTracks = fESD->GetNumberOfTracks();
 228   AliDebug(2,Form("nTracks %d", nTracks));
 229
 230   if(!fTrackCuts) {
 231     // Post output data
 232     PostData(0,fHistList);
 233     PostData(1,fCFManagerPos->GetParticleContainer());
 234     PostData(2,fCFManagerNeg->GetParticleContainer());
 235     return;
 236   }
 237
 238   // Selected events for analysis
 239   fNEventSel->Fill(0.);
 240
 241
 242   Double_t containerInputRec[5] ;
 243   Double_t containerInputTPConly[5];
 244   Double_t containerInputMC[5];
 245   //Now go to rec level
 246   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 247     {
 248       if(!fESD->GetTrack(iTrack) ) continue;
 249       AliESDtrack* track = fESD->GetTrack(iTrack);
 250       if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
 251       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
 252       if(!track || !trackTPC) continue;
 253
 254       Float_t dca2D, dcaZ;
 255       track->GetImpactParameters(dca2D,dcaZ);
 256       Float_t dca2DTPC, dcaZTPC;
 257       track->GetImpactParametersTPC(dca2DTPC,dcaZTPC);
 258       Float_t chi2PerClusterTPC = -1.;
 259       Float_t nClustersTPC = track->GetTPCNcls();//track->GetTPCclusters(0);
 260       if(nClustersTPC>0.) chi2PerClusterTPC = track->GetTPCchi2()/(2.*nClustersTPC-5.);
 261       Float_t chi2PerClusterTPCIter1 = -1.;
 262       Float_t nClustersTPCIter1 = track->GetTPCNclsIter1();
 263       if(nClustersTPCIter1>0.) chi2PerClusterTPCIter1 = track->GetTPCchi2Iter1()/(2.*nClustersTPCIter1-5.);
 264
 265       //fill the container
 266       containerInputRec[0] = track->Pt();
 267       containerInputRec[1] = track->Phi();
 268       containerInputRec[2] = track->Eta();
 269       containerInputRec[3] = dca2D;
 270       containerInputRec[4] = chi2PerClusterTPC;
 271
 272       //Store TPC Inner Params for TPConly tracks
 273       containerInputTPConly[0] = trackTPC->Pt();
 274       containerInputTPConly[1] = trackTPC->Phi();
 275       containerInputTPConly[2] = trackTPC->Eta();
 276       containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same container. Should be corrected back when looking at output.
 277       containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
 278
 279       AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
 280       if(trackTPCESD) {
 281         if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
 282           if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 283           if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 284         }
 285       }
 286
 287       if (fTrackCuts->AcceptTrack(track)) {
 288         if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 289         if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 290
 291
 292         //Only fill the MC containers if MC information is available
 293         if(eventHandler) {
 294           Int_t label = TMath::Abs(track->GetLabel());
 295           TParticle *particle = stack->Particle(label) ;
 296           if(!particle) continue;
 297
 298           containerInputMC[0] = particle->Pt();
 299           containerInputMC[1] = particle->Phi();
 300           containerInputMC[2] = particle->Eta();
 301           containerInputMC[3] = 0.0;
 302           containerInputMC[4] = 0.0;
 303
 304           //Container with primaries
 305           if(stack->IsPhysicalPrimary(label)) {
 306             if(particle->GetPDG()->Charge()>0.) {
 307               fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
 308             }
 309             if(particle->GetPDG()->Charge()<0.) {
 310               fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
 311             }
 312           }
 313
 314           //Container with secondaries
 315           if (!stack->IsPhysicalPrimary(label) ) {
 316             if(particle->GetPDG()->Charge()>0.) {
 317               fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 318             }
 319             if(particle->GetPDG()->Charge()<0.) {
 320               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 321             }
 322           }
 323         }
 324
 325       }//trackCuts
 326
 327       delete trackTPCESD;
 328     }//track loop
 329
 330
 331   //Fill MC containters if particles are findable
 332   if(eventHandler) {
 333     for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
 334       {
 335         AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
 336         if(!mcPart) continue;
 337
 338         //fill the container
 339         containerInputMC[0] = mcPart->Pt();
 340         containerInputMC[1] = mcPart->Phi();
 341         containerInputMC[2] = mcPart->Eta();
 342         containerInputMC[3] = 0.0;
 343         containerInputMC[4] = 0.0;
 344
 345         if(stack->IsPhysicalPrimary(iPart)) {
 346           if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 347           if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
 348         }
 349       }
 350   }
 351
 352   PostData(0,fHistList);
 353   PostData(1,fCFManagerPos->GetParticleContainer());
 354   PostData(2,fCFManagerNeg->GetParticleContainer());
 355
 356 }
 357
 358
 359 //___________________________________________________________________________
 360 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 361 {
 362   // The Terminate() function is the last function to be called during
 363   // a query. It always runs on the client, it can be used to present
 364   // the results graphically or save the results to file.
 365
 366 }
 367
 368 //___________________________________________________________________________
 369 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 370   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 371   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 372   //
 373   AliDebug(2,Form("CreateOutputObjects CreateOutputObjects of task %s", GetName()));
 374
 375   Bool_t oldStatus = TH1::AddDirectoryStatus();
 376   TH1::AddDirectory(kFALSE);
 377
 378   //slot #1
 379   OpenFile(0);
 380   fHistList = new TList();
 381   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 382   fHistList->Add(fNEventAll);
 383   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 384   fHistList->Add(fNEventSel);
 385
 386   TH1::AddDirectory(oldStatus);
 387
 388 }
 389
 390 #endif