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