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