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   fTrigger(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   fTrigger(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   DefineOutput(0,TList::Class());
  97   DefineOutput(1,AliCFContainer::Class());
  98   DefineOutput(2,AliCFContainer::Class());
  99 }
 100
 101 // //___________________________________________________________________________
 102 // AliPWG4HighPtSpectra& AliPWG4HighPtSpectra::operator=(const AliPWG4HighPtSpectra& c)
 103 // {
 104 //   //
 105 //   // Assignment operator
 106 //   //
 107 //   if (this!=&c) {
 108 //     AliAnalysisTask::operator=(c) ;
 109 //     fReadAODData = c.fReadAODData ;
 110 //     fCFManagerPos  = c.fCFManagerPos;
 111 //     fCFManagerNeg  = c.fCFManagerNeg;
 112 //     fHistList = c.fHistList;
 113 //     fNEventAll = c.fNEventAll;
 114 //     fNEventSel = c.fNEventSel;
 115 //   }
 116 //   return *this;
 117 // }
 118
 119 // //___________________________________________________________________________
 120 // AliPWG4HighPtSpectra::AliPWG4HighPtSpectra(const AliPWG4HighPtSpectra& c) :
 121 //   AliAnalysisTask(c),
 122 //   fReadAODData(c.fReadAODData),
 123 //   fCFManagerPos(c.fCFManagerPos),
 124 //   fCFManagerNeg(c.fCFManagerNeg),
 125 //   fESD(c.fESD),
 126 //   fTrackCuts(c.fTrackCuts),
 127 //   fTrigger(c.fTrigger),
 128 //   fHistList(c.fHistList),
 129 //   fNEventAll(c.fNEventAll),
 130 //   fNEventSel(c.fNEventSel)
 131 // {
 132 //   //
 133 //   // Copy Constructor
 134 //   //
 135 // }
 136
 137 // //___________________________________________________________________________
 138 // AliPWG4HighPtSpectra::~AliPWG4HighPtSpectra() {
 139 //   //
 140 //   //destructor
 141 //   //
 142 //   Info("~AliPWG4HighPtSpectra","Calling Destructor");
 143 //   if (fCFManagerPos)           delete fCFManagerPos ;
 144 //   if (fCFManagerNeg)           delete fCFManagerNeg ;
 145 //   if (fNEventAll)              delete fNEventAll ;
 146 //   if (fNEventSel)              delete fNEventSel ;
 147 // }
 148 //________________________________________________________________________
 149 void AliPWG4HighPtSpectra::ConnectInputData(Option_t *)
 150 {
 151   // Connect ESD here
 152   // Called once
 153   AliDebug(2,Form(">> AliPWG4HighPtSpectra::ConnectInputData \n"));
 154   //  cout << "cout >> AliPWG4HighPtSpectra::ConnectInputData" << endl;
 155   printf(">> AliPWG4HighPtSpectra::ConnectInputData \n");
 156
 157   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 158   if (!tree) {
 159     AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
 160   } else {
 161
 162     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 163
 164     if (!esdH) {
 165       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 166     } else {
 167       fESD = esdH->GetEvent();
 168     }
 169   }
 170
 171 }
 172 //_________________________________________________
 173 void AliPWG4HighPtSpectra::Exec(Option_t *)
 174 {
 175   //
 176   // Main loop function
 177   //
 178   AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
 179
 180   // All events without selection
 181   fNEventAll->Fill(0.);
 182
 183   if (!fESD) {
 184     AliDebug(2,Form("ERROR: fESD not available"));
 185     PostData(0,fHistList);
 186     PostData(1,fCFManagerPos->GetParticleContainer());
 187     PostData(2,fCFManagerNeg->GetParticleContainer());
 188     return;
 189   }
 190
 191   Bool_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 192   if(!isSelected) { //Select collison candidates
 193     AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
 194     PostData(0,fHistList);
 195     PostData(1,fCFManagerPos->GetParticleContainer());
 196     PostData(2,fCFManagerNeg->GetParticleContainer());
 197     return;
 198   }
 199
 200   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 201   // This handler can return the current MC event
 202
 203   AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 204   //  AliMCEventHandler* eventHandler = (AliMCEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 205
 206   AliStack* stack = 0x0;
 207   AliMCEvent* mcEvent = 0x0;
 208
 209   if(eventHandler) {
 210     mcEvent = eventHandler->MCEvent();
 211     if (!mcEvent) {
 212       AliDebug(2,Form("ERROR: Could not retrieve MC event"));
 213       PostData(0,fHistList);
 214       PostData(1,fCFManagerPos->GetParticleContainer());
 215     PostData(2,fCFManagerNeg->GetParticleContainer());
 216     return;
 217     }
 218
 219     AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
 220
 221     stack = mcEvent->Stack();                //Particles Stack
 222
 223     AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
 224   }
 225
 226   const AliESDVertex *vtx = fESD->GetPrimaryVertex();
 227   AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
 228   // Need vertex cut
 229   if (vtx->GetNContributors() < 2) {
 230     PostData(0,fHistList);
 231     PostData(1,fCFManagerPos->GetParticleContainer());
 232     PostData(2,fCFManagerNeg->GetParticleContainer());
 233     return;
 234   }
 235
 236   double primVtx[3];
 237   vtx->GetXYZ(primVtx);
 238   if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
 239     PostData(0,fHistList);
 240     PostData(1,fCFManagerPos->GetParticleContainer());
 241     PostData(2,fCFManagerNeg->GetParticleContainer());
 242     return;
 243   }
 244
 245   if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2){
 246     // Post output data
 247     PostData(0,fHistList);
 248     PostData(1,fCFManagerPos->GetParticleContainer());
 249     PostData(2,fCFManagerNeg->GetParticleContainer());
 250     return;
 251   }
 252   Int_t nTracks = fESD->GetNumberOfTracks();
 253   AliDebug(2,Form("nTracks %d", nTracks));
 254
 255   if(!fTrackCuts) {
 256     // Post output data
 257     PostData(0,fHistList);
 258     PostData(1,fCFManagerPos->GetParticleContainer());
 259     PostData(2,fCFManagerNeg->GetParticleContainer());
 260     return;
 261   }
 262
 263   // Selected events for analysis
 264   fNEventSel->Fill(0.);
 265
 266
 267   Double_t containerInputRec[5] ;
 268   Double_t containerInputTPConly[5];
 269   Double_t containerInputMC[5];
 270   //Now go to rec level
 271   for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
 272     {
 273       if(!fESD->GetTrack(iTrack) ) continue;
 274       AliESDtrack* track = fESD->GetTrack(iTrack);
 275       if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
 276       AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
 277       if(!track || !trackTPC) continue;
 278
 279       Float_t dca2D, dcaZ;
 280       track->GetImpactParameters(dca2D,dcaZ);
 281       Float_t dca2DTPC, dcaZTPC;
 282       track->GetImpactParametersTPC(dca2DTPC,dcaZTPC);
 283       Float_t chi2PerClusterTPC = -1.;
 284       Float_t nClustersTPC = track->GetTPCNcls();//track->GetTPCclusters(0);
 285       if(nClustersTPC>0.) chi2PerClusterTPC = track->GetTPCchi2()/(2.*nClustersTPC-5.);
 286       Float_t chi2PerClusterTPCIter1 = -1.;
 287       Float_t nClustersTPCIter1 = track->GetTPCNclsIter1();
 288       if(nClustersTPCIter1>0.) chi2PerClusterTPCIter1 = track->GetTPCchi2Iter1()/(2.*nClustersTPCIter1-5.);
 289
 290       //fill the container
 291       containerInputRec[0] = track->Pt();
 292       containerInputRec[1] = track->Phi();
 293       containerInputRec[2] = track->Eta();
 294       containerInputRec[3] = dca2D;
 295       containerInputRec[4] = chi2PerClusterTPC;
 296
 297       containerInputTPConly[0] = trackTPC->Pt();
 298       containerInputTPConly[1] = trackTPC->Phi();
 299       containerInputTPConly[2] = trackTPC->Eta();
 300       containerInputTPConly[3] = dca2DTPC/10.; //Divide by 10 in order to store in same containter. Should be corrected back when looking at output.
 301       containerInputTPConly[4] = chi2PerClusterTPCIter1;//TPC;
 302
 303       if (fTrackCuts->AcceptTrack(track)) {
 304         if(track->GetSign()>0.) {
 305           fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 306           fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 307         }
 308         if(track->GetSign()<0.) {
 309           fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
 310           fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
 311         }
 312
 313
 314         //Only fill the secondary particle container if MC information is available
 315         if(eventHandler) {
 316           Int_t label = TMath::Abs(track->GetLabel());
 317           TParticle *particle = stack->Particle(label) ;
 318           if(!particle) continue;
 319
 320           containerInputMC[0] = particle->Pt();
 321           containerInputMC[1] = particle->Phi();
 322           containerInputMC[2] = particle->Eta();
 323           containerInputMC[3] = 0.0;
 324           containerInputMC[4] = 0.0;
 325
 326           if(particle->GetPDG()->Charge()>0.) {
 327             fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
 328           }
 329           if(particle->GetPDG()->Charge()<0.) {
 330             fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
 331           }
 332
 333           if (!stack->IsPhysicalPrimary(label) ) {
 334             if(particle->GetPDG()->Charge()>0.) {
 335               fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 336             }
 337             if(particle->GetPDG()->Charge()<0.) {
 338               fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
 339             }
 340           }
 341         }
 342
 343       }
 344
 345     }
 346
 347
 348   //Fill MC containters if particles are findable
 349   if(eventHandler) {
 350     for(int iPart = 1; iPart<(mcEvent->GetNumberOfTracks()); iPart++)//stack->GetNprimary();
 351       {
 352         AliMCParticle *mcPart  = (AliMCParticle*)mcEvent->GetTrack(iPart);
 353         if(!mcPart) continue;
 354
 355         //fill the container
 356         containerInputMC[0] = mcPart->Pt();
 357         containerInputMC[1] = mcPart->Phi();
 358         containerInputMC[2] = mcPart->Eta();
 359         containerInputMC[3] = 0.0;
 360         containerInputMC[4] = 0.0;
 361
 362         int counter;
 363         Float_t trackLengthTPC = 0.;
 364         if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(3,mcPart)) {
 365           trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
 366           if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
 367         }
 368         if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(3,mcPart)) {
 369           trackLengthTPC = mcPart->GetTPCTrackLength(fESD->GetMagneticField(),0.1,counter,3.0);
 370           if(trackLengthTPC > (float)fTrackCuts->GetMinNClusterTPC()) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCtrackable) ;
 371         }
 372       }
 373   }
 374
 375   PostData(0,fHistList);
 376   PostData(1,fCFManagerPos->GetParticleContainer());
 377   PostData(2,fCFManagerNeg->GetParticleContainer());
 378
 379 }
 380
 381
 382 //___________________________________________________________________________
 383 void AliPWG4HighPtSpectra::Terminate(Option_t*)
 384 {
 385   // The Terminate() function is the last function to be called during
 386   // a query. It always runs on the client, it can be used to present
 387   // the results graphically or save the results to file.
 388 }
 389
 390 //___________________________________________________________________________
 391 void AliPWG4HighPtSpectra::CreateOutputObjects() {
 392   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 393   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 394   //
 395   AliDebug(2,Form("CreateOutputObjects","CreateOutputObjects of task %s", GetName()));
 396
 397   Bool_t oldStatus = TH1::AddDirectoryStatus();
 398   TH1::AddDirectory(kFALSE);
 399
 400   //slot #1
 401   OpenFile(0);
 402   fHistList = new TList();
 403   fNEventAll = new TH1F("fNEventAll","NEventAll",1,-0.5,0.5);
 404   fHistList->Add(fNEventAll);
 405   fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
 406   fHistList->Add(fNEventSel);
 407
 408   TH1::AddDirectory(oldStatus);
 409
 410 }
 411
 412 #endif