PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.cxx

   1 #include "TCanvas.h"
   2 #include "TParticle.h"
   3 #include "TH1I.h"
   4 #include <TDatabasePDG.h>
   5
   6 #include "AliAnalysisDataSlot.h"
   7 #include "AliAnalysisDataContainer.h"
   8 #include "AliAnalysisManager.h"
   9 #include "AliAODEvent.h"
  10 #include "AliAODMCParticle.h"
  11 #include "AliAODTrack.h"
  12 #include "AliAODTracklets.h"
  13 #include "AliMultiplicity.h"
  14 #include "AliCFManager.h"
  15 #include "AliCFCutBase.h"
  16 #include "AliCFContainer.h"
  17 #include "AliESDEvent.h"
  18 #include "AliESDtrack.h"
  19 #include "AliESDtrackCuts.h"
  20 #include "AliESDVertex.h"
  21 #include "AliInputEventHandler.h"
  22 #include "AliMCEvent.h"
  23 #include "AliVEvent.h"
  24 #include "AliAODEvent.h"
  25 #include "AliVVertex.h"
  26 #include "AliLog.h"
  27 #include "AliStack.h"
  28 #include "AliGenEventHeader.h"
  29 #include "AliAODMCHeader.h"
  30 #include "AliCentrality.h"
  31
  32 #include "AliSingleTrackEffCuts.h"
  33 #include "AliCFSingleTrackEfficiencyTask.h"
  34
  35
  36 ClassImp(AliCFSingleTrackEfficiencyTask)
  37
  38 //__________________________________________________________________________
  39 AliCFSingleTrackEfficiencyTask::AliCFSingleTrackEfficiencyTask() :
  40   fReadAODData(0),
  41   fCFManager(0x0),
  42   fQAHistList(0x0),
  43   fTrackCuts(0x0),
  44   fMCCuts(0x0),
  45   fTriggerMask(AliVEvent::kAny),
  46   fSetFilterBit(kFALSE),
  47   fbit(0),
  48   fEvalCentrality(kFALSE),
  49   fCentralityEstimator("V0M"),
  50   fHistEventsProcessed(0x0)
  51 {
  52   //
  53   //Default constructor
  54   //
  55 }
  56
  57 //___________________________________________________________________________
  58 AliCFSingleTrackEfficiencyTask::AliCFSingleTrackEfficiencyTask(const Char_t* name,AliESDtrackCuts *trackcuts, AliSingleTrackEffCuts * mccuts) :
  59   AliAnalysisTaskSE(name),
  60   fReadAODData(0),
  61   fCFManager(0x0),
  62   fQAHistList(0x0),
  63   fTrackCuts(trackcuts),
  64   fMCCuts(mccuts),
  65   fTriggerMask(AliVEvent::kAny),
  66   fSetFilterBit(kFALSE),
  67   fbit(0),
  68   fEvalCentrality(kFALSE),
  69   fCentralityEstimator("V0M"),
  70   fHistEventsProcessed(0x0)
  71 {
  72   //
  73   // Constructor. Initialization of Inputs and Outputs
  74   //
  75   Info("AliCFSingleTrackEfficiencyTask","Calling Constructor");
  76
  77   // Output slot #1 writes into a TList container (nevents histogran)
  78   DefineOutput(1,TH1I::Class());
  79   // Output slot #2 writes into a TList container (distributions)
  80   DefineOutput(2,AliCFContainer::Class());
  81   // Output slot #3 writes the QA list
  82   DefineOutput(3,TList::Class());
  83   // Output slot #3 writes the ESD track cuts used
  84   DefineOutput(4,AliESDtrackCuts::Class());
  85   // Output slot #4 writes the particle and event selection object
  86   DefineOutput(5,AliSingleTrackEffCuts::Class());
  87 }
  88
  89 //_________________________________________________________________________________________________________________
  90 AliCFSingleTrackEfficiencyTask& AliCFSingleTrackEfficiencyTask::operator=(const AliCFSingleTrackEfficiencyTask& c)
  91 {
  92   //
  93   // Assignment operator
  94   //
  95   if (this!=&c) {
  96     AliAnalysisTaskSE::operator=(c) ;
  97
  98     fReadAODData = c.fReadAODData;
  99     fCFManager  = c.fCFManager;
 100     fQAHistList = c.fQAHistList;
 101
 102     if(c.fTrackCuts) { delete fTrackCuts; fTrackCuts = new AliESDtrackCuts(*(c.fTrackCuts)); }
 103     if(c.fMCCuts) { delete fMCCuts; fMCCuts = new AliSingleTrackEffCuts(*(c.fMCCuts)); }
 104     fTriggerMask = c.fTriggerMask;
 105
 106     fSetFilterBit  = c.fSetFilterBit;
 107     fbit = c.fbit;
 108
 109     fEvalCentrality = c.fEvalCentrality;
 110     fCentralityEstimator = c.fCentralityEstimator;
 111
 112     fHistEventsProcessed = c.fHistEventsProcessed;
 113   }
 114   return *this;
 115 }
 116
 117 //________________________________________________________________________________________________________
 118 AliCFSingleTrackEfficiencyTask::AliCFSingleTrackEfficiencyTask(const AliCFSingleTrackEfficiencyTask& c) :
 119   AliAnalysisTaskSE(c),
 120   fReadAODData(c.fReadAODData),
 121   fCFManager(c.fCFManager),
 122   fQAHistList(c.fQAHistList),
 123   fTrackCuts(c.fTrackCuts),
 124   fMCCuts(c.fMCCuts),
 125   fTriggerMask(c.fTriggerMask),
 126   fSetFilterBit(c.fSetFilterBit),
 127   fbit(c.fbit),
 128   fEvalCentrality(c.fEvalCentrality),
 129   fCentralityEstimator(c.fCentralityEstimator),
 130   fHistEventsProcessed(c.fHistEventsProcessed)
 131 {
 132   //
 133   // Copy Constructor
 134   //
 135 }
 136
 137 //___________________________________________________________________________
 138 AliCFSingleTrackEfficiencyTask::~AliCFSingleTrackEfficiencyTask()
 139 {
 140   //
 141   // Destructor
 142   //
 143   Info("~AliCFSingleTrackEfficiencyTask","Calling Destructor");
 144
 145   if (fCFManager)           delete fCFManager;
 146   if (fHistEventsProcessed) delete fHistEventsProcessed;
 147   if (fQAHistList) { fQAHistList->Clear(); delete fQAHistList; }
 148   if (fTrackCuts)           delete fTrackCuts;
 149   if (fMCCuts)              delete fMCCuts;
 150 }
 151
 152 //_________________________________________________
 153 void AliCFSingleTrackEfficiencyTask::Init()
 154 {
 155   //
 156   // Initialization, checks + copy cuts
 157   //
 158   if(!fMCCuts) {
 159     AliFatal(" MC Cuts not defined");
 160     return;
 161   }
 162   if(!fTrackCuts) {
 163     AliFatal(" Track Cuts not defined");
 164     return;
 165   }
 166
 167   AliESDtrackCuts* copyfTrackCuts = new AliESDtrackCuts(*fTrackCuts);
 168   const char* nameoutput = GetOutputSlot(4)->GetContainer()->GetName();
 169   copyfTrackCuts->SetName(nameoutput);
 170   // Post the data
 171   PostData(4,copyfTrackCuts);
 172
 173   AliSingleTrackEffCuts* copyfMCCuts = new AliSingleTrackEffCuts(*fMCCuts);
 174   nameoutput = GetOutputSlot(5)->GetContainer()->GetName();
 175   copyfMCCuts->SetName(nameoutput);
 176   // Post the data
 177   PostData(5,copyfMCCuts);
 178
 179
 180   if(fEvalCentrality) {
 181     Bool_t isCentEstimatorOk = kFALSE;
 182     TString validEstimators[9] = { "V0M", "V0A", "V0C", "TRK", "TKL", "CL1", "ZNA", "ZNC" "ZPA" };
 183     for(Int_t i=0; i<9; i++ ) {
 184       if(fCentralityEstimator==validEstimators[i]) isCentEstimatorOk = kTRUE;
 185     }
 186     if(!isCentEstimatorOk) {
 187       AliFatal(Form("Chosen centrality estimator %s is not valid\n",fCentralityEstimator.Data()));
 188       return;
 189     }
 190   }
 191
 192   return;
 193 }
 194
 195 //_________________________________________________________________
 196 void AliCFSingleTrackEfficiencyTask::UserExec(Option_t *)
 197 {
 198   //
 199   // User Exec
 200   //
 201
 202   Info("UserExec","Start of method") ;
 203
 204   AliVEvent* event = fInputEvent;
 205
 206   if(!fInputEvent) {
 207     AliFatal("NO EVENT FOUND!");
 208     return;
 209   }
 210   if(!fMCEvent) {
 211     AliFatal("NO MC INFO FOUND");
 212     return;
 213   }
 214
 215   fHistEventsProcessed->Fill(0.5); // # of Event proceed
 216   Bool_t IsEventMCSelected = kFALSE;
 217   Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
 218
 219   //
 220   // Step 0: MC Gen Event Selection
 221   //
 222   if(isAOD) {
 223     if(!event && AODEvent() && IsStandardAOD()) {
 224       // In case there is an AOD handler writing a standard AOD, use the AOD
 225       // event in memory rather than the input (ESD) event.
 226       event = dynamic_cast<AliAODEvent*> (AODEvent());
 227     }
 228     IsEventMCSelected = fMCCuts->IsMCEventSelected(event);//AODs
 229   } else {
 230     IsEventMCSelected = fMCCuts->IsMCEventSelected(fMCEvent);//ESDs
 231   }
 232
 233   // pass to the manager the event info to the cuts that need it
 234   if(!isAOD) fCFManager->SetMCEventInfo(fMCEvent);
 235   else fCFManager->SetMCEventInfo(event);
 236   fCFManager->SetRecEventInfo(event);
 237
 238   if(!IsEventMCSelected) {
 239     AliDebug(3,"MC event not passing the quality criteria \n");
 240     PostData(1,fHistEventsProcessed);
 241     PostData(2,fCFManager->GetParticleContainer());
 242     PostData(3,fQAHistList);
 243     return;
 244   }
 245   fHistEventsProcessed->Fill(1.5); // # of Event after passing MC cuts
 246
 247
 248   //
 249   // Step 1-3: Check the MC generated particles
 250   //
 251   if(!isAOD) CheckESDMCParticles();
 252   else CheckAODMCParticles();
 253
 254   //
 255   // Step 4-7: Reconstructed event and track selection
 256   //
 257   Bool_t isRecoEventOk = fMCCuts->IsRecoEventSelected(event);
 258
 259   if(isRecoEventOk) {
 260     fHistEventsProcessed->Fill(2.5); // # of Event after passing all cuts
 261     CheckReconstructedParticles();
 262   }
 263   else AliDebug(3,"Event not passing quality criteria\n");
 264
 265
 266   PostData(1,fHistEventsProcessed);
 267   PostData(2,fCFManager->GetParticleContainer());
 268   PostData(3,fQAHistList);
 269
 270   return;
 271 }
 272
 273
 274 //___________________________________________________________________________
 275 void AliCFSingleTrackEfficiencyTask::Terminate(Option_t*)
 276 {
 277   //
 278   // Terminate
 279   //
 280
 281   Info("Terminate","Start and end of Method");
 282   AliAnalysisTaskSE::Terminate();
 283
 284   /*
 285   //draw some example plots....
 286   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 287
 288   TH1D* h00 =   cont->ShowProjection(5,0) ;
 289   TH1D* h01 =   cont->ShowProjection(5,1) ;
 290   TH1D* h02 =   cont->ShowProjection(5,2) ;
 291   TH1D* h03 =   cont->ShowProjection(5,3) ;
 292   TH1D* h04 =   cont->ShowProjection(5,4) ;
 293   TH1D* h05 =   cont->ShowProjection(5,5) ;
 294   TH1D* h06 =   cont->ShowProjection(5,6) ;
 295   TH1D* h07 =   cont->ShowProjection(5,7) ;
 296
 297   h00->SetMarkerStyle(23) ;
 298   h01->SetMarkerStyle(24) ;
 299   h02->SetMarkerStyle(25) ;
 300   h03->SetMarkerStyle(26) ;
 301   h04->SetMarkerStyle(27) ;
 302   h05->SetMarkerStyle(28) ;
 303   h06->SetMarkerStyle(28) ;
 304
 305
 306
 307   TCanvas * c =new TCanvas("c","",1400,800);
 308   c->Divide(4,2);
 309
 310   c->cd(1);
 311   h00->Draw("p");
 312   c->cd(2);
 313   h01->Draw("p");
 314   c->cd(3);
 315   h02->Draw("p");
 316   c->cd(4);
 317   h03->Draw("p");
 318   c->cd(5);
 319   h04->Draw("p");
 320   c->cd(6);
 321   h05->Draw("p");
 322   c->cd(7);
 323   h06->Draw("p");
 324   c->cd(8);
 325   h07->Draw("p");
 326
 327   c->SaveAs("plots.eps");
 328   */
 329 }
 330
 331
 332 //___________________________________________________________________________
 333 void AliCFSingleTrackEfficiencyTask::UserCreateOutputObjects()
 334 {
 335   //
 336   // UserCreateOutputObjects
 337   //
 338
 339   Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
 340
 341   //slot #1
 342   OpenFile(1);
 343
 344   const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
 345   //       fHistEventsProcessed = new TH1I(nameoutput"fHistEventsProcessed","fHistEventsProcessed",3,0,3) ;
 346   fHistEventsProcessed = new TH1I(nameoutput,"fHistEventsProcessed",3,0,3) ;
 347   fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"All events");
 348   fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Good MC events");
 349   fHistEventsProcessed->GetXaxis()->SetBinLabel(3,"Good Reconstructed events");
 350
 351   PostData(1,fHistEventsProcessed) ;
 352   PostData(2,fCFManager->GetParticleContainer()) ;
 353   PostData(3,fQAHistList) ;
 354
 355   return;
 356 }
 357
 358
 359 //_________________________________________________________________________
 360 void AliCFSingleTrackEfficiencyTask::CheckESDMCParticles()
 361 {
 362   //
 363   // Check ESD generated particles
 364   //
 365   if (!fMCEvent) {
 366     AliFatal("NO MC INFO FOUND");
 367     return;
 368   }
 369
 370   Double_t containerInput[7] = { 0., 0., 0., 0., 0., 0., 0. };
 371
 372   TArrayF vtxPos(3);
 373   AliGenEventHeader *genHeader = NULL;
 374   genHeader = fMCEvent->GenEventHeader();
 375   genHeader->PrimaryVertex(vtxPos);
 376   containerInput[4] = vtxPos[2]; // z-vtx position
 377
 378   Double_t multiplicity = (Double_t)GetNumberOfTrackletsInEtaRange(-1.0,1.0);
 379   containerInput[5] = multiplicity; //reconstructed number of tracklets
 380
 381   Double_t centrality = -1.;
 382   if(fEvalCentrality) centrality = GetCentrality();
 383   containerInput[6] = centrality;
 384
 385   // loop on the MC Generated Particles
 386   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 387
 388     AliMCParticle *mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
 389     containerInput[0] = (Float_t)mcPart->Pt();
 390     containerInput[1] = mcPart->Eta();
 391     containerInput[2] = mcPart->Phi();
 392     containerInput[3] = mcPart->Theta();
 393
 394     // Step 1. Particle passing through Generation criteria and filling
 395     if( !fMCCuts->IsMCParticleGenerated(mcPart) ) {
 396       AliDebug(3,"MC Particle not passing through genetations criteria\n");
 397       continue;
 398     }
 399     fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCGenCut);
 400
 401     // Step 2. Particle passing through Kinematic criteria and filling
 402     if( !fMCCuts->IsMCParticleInKineAcceptance(mcPart) ) {
 403       AliDebug(3,"MC Particle not in the kine acceptance\n");
 404       continue;
 405     }
 406     fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCKineCut);
 407
 408     // Step 3. Particle passing through Track ref criteria and filling
 409     // did leave signal (enough clusters) on the detector
 410     if( !fMCCuts->IsMCParticleInReconstructable(mcPart) ) {
 411       AliDebug(3,"MC Particle not in the reconstructible\n");
 412       continue;
 413     }
 414     fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCAccpCut);
 415
 416   }// end of particle loop
 417
 418   return;
 419 }
 420
 421
 422 //________________________________________________________________________
 423 void AliCFSingleTrackEfficiencyTask::CheckAODMCParticles()
 424 {
 425   //
 426   // Check AOD generated particles
 427   //
 428   if (!fInputEvent) {
 429     AliFatal("NO EVENT FOUND!");
 430     return;
 431   }
 432
 433   AliAODEvent* event = dynamic_cast<AliAODEvent*>(fInputEvent);
 434   if(!event && AODEvent() && IsStandardAOD()) {
 435     // In case there is an AOD handler writing a standard AOD, use the AOD
 436     // event in memory rather than the input (ESD) event.
 437     event = dynamic_cast<AliAODEvent*> (AODEvent());
 438   }
 439   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
 440   if (!mcArray) {
 441     AliError("Could not find Monte-Carlo in AOD");
 442     return;
 443   }
 444   AliAODMCHeader *mcHeader=NULL;
 445   mcHeader = dynamic_cast<AliAODMCHeader*>(event->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 446   if (!mcHeader) {
 447     AliError("Could not find MC Header in AOD");
 448     return;
 449   }
 450
 451   Double_t containerInput[7] = { 0., 0., 0., 0., 0., 0., 0. };
 452   // Set the z-vertex position
 453   containerInput[4]  = mcHeader->GetVtxZ();
 454   // Multiplicity of the event defined as Ntracklets |eta|<1.0
 455   Double_t multiplicity = (Double_t)GetNumberOfTrackletsInEtaRange(-1.0,1.0);
 456   containerInput[5] = multiplicity;
 457   // Determine the event centrality
 458   Double_t centrality = 0.;
 459   if(fEvalCentrality) centrality = GetCentrality();
 460   containerInput[6] = centrality;
 461
 462   for (Int_t ipart=0; ipart<mcArray->GetEntriesFast(); ipart++) {
 463
 464     AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(ipart));
 465     if (!mcPart){
 466       AliError("Failed casting particle from MC array!, Skipping particle");
 467       continue;
 468     }
 469     containerInput[0] = (Float_t)mcPart->Pt();
 470     containerInput[1] = mcPart->Eta();
 471     containerInput[2] = mcPart->Phi();
 472     containerInput[3] = mcPart->Theta();
 473
 474     // Step 1. Particle passing through Generation criteria and filling
 475     if( !fMCCuts->IsMCParticleGenerated(mcPart) ) {
 476       AliDebug(3,"MC Particle not passing quality criteria\n");
 477       continue;
 478     }
 479     fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCGenCut);
 480
 481     // Step 2. Particle passing through Kinematic criteria and filling
 482     if( !fMCCuts->IsMCParticleInKineAcceptance(mcPart) ) {
 483       AliDebug(3,"MC Particle not in the acceptance\n");
 484       continue;
 485     }
 486     fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCKineCut);
 487
 488     // Step 3. Particle passing through Track Ref criteria and filling
 489     // but no info available for Track ref in AOD fillng same as above
 490     fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCAccpCut);
 491   }
 492
 493   return;
 494 }
 495
 496
 497 //_______________________________________________________________________________
 498 AliESDtrack * AliCFSingleTrackEfficiencyTask::ConvertTrack(AliAODTrack *track)
 499 {
 500   //
 501   // Convert an AOD track to an  ESD track to apply ESDtrackCuts
 502   //
 503
 504   AliAODEvent *aodEvent = (AliAODEvent*)fInputEvent;
 505   const AliAODVertex *primary = aodEvent->GetPrimaryVertex();
 506   Double_t pos[3],cov[6];
 507   primary->GetXYZ(pos);
 508   primary->GetCovarianceMatrix(cov);
 509   const AliESDVertex vESD(pos,cov,100.,100);
 510
 511   AliESDtrack *esdTrack =  new AliESDtrack(track);
 512   // set the TPC cluster info
 513   esdTrack->SetTPCClusterMap(track->GetTPCClusterMap());
 514   esdTrack->SetTPCSharedMap(track->GetTPCSharedMap());
 515   esdTrack->SetTPCPointsF(track->GetTPCNclsF());
 516   // needed to calculate the impact parameters
 517   esdTrack->RelateToVertex(&vESD,0.,3.);
 518   //  std::cout << " primary vtx "<< primary << std::endl;
 519   //  std::cout << " esdvtx "<< vESD.GetName() << std::endl;
 520   //  std::cout<< " esdtrack pt "<< esdTrack.Pt() << " and status " << esdTrack.GetStatus() <<endl;
 521   //  std::cout << " aod track "<< track<< " and status " << track->GetStatus() << std::endl;
 522   //  std::cout << " esd track "<< esdTrack<< " and status " << esdTrack->GetStatus() << std::endl;
 523   return esdTrack;
 524 }
 525
 526
 527 //___________________________________________________________________________
 528 void AliCFSingleTrackEfficiencyTask::CheckReconstructedParticles()
 529 {
 530   //
 531   // Check reconstructed particles
 532   //
 533
 534   AliVEvent* event = fInputEvent;
 535   Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
 536   if(!event && AODEvent() && IsStandardAOD()) {
 537    event  = dynamic_cast<AliAODEvent*> (AODEvent());
 538   }
 539   if (!event) return;
 540
 541   Double_t containerInput[7] = { 0, 0, 0, 0, 0, 0, 0};   // contains reconstructed quantities
 542   Double_t containerInputMC[7] = { 0, 0, 0, 0, 0, 0, 0}; // contains generated quantities
 543
 544   const AliVVertex *vertex = event->GetPrimaryVertex();
 545   // set the z-vertex position
 546   containerInput[4] = vertex->GetZ();
 547   containerInputMC[4] = containerInput[4];
 548   // set the event multiplicity as Ntracklets in |eta|<1.0
 549   Double_t multiplicity = (Double_t)GetNumberOfTrackletsInEtaRange(-1.0,1.0);
 550   containerInput[5] = multiplicity;
 551   containerInputMC[5] = multiplicity;
 552   // Determine the event centrality
 553   Double_t centrality = 0.;
 554   if(fEvalCentrality) centrality = GetCentrality();
 555   containerInput[6] = centrality;
 556   containerInputMC[6] = centrality;
 557
 558   // Reco tracks track loop
 559   AliVParticle* track = NULL;
 560   for (Int_t iTrack = 0; iTrack<event->GetNumberOfTracks(); iTrack++) {
 561
 562     track = event->GetTrack(iTrack);
 563     if(!track) {
 564       AliDebug(3,Form("Track %d not found",iTrack));
 565       continue;
 566     }
 567
 568     Double_t mom[3];
 569     track->PxPyPz(mom);
 570     Double_t pt = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
 571     containerInput[0] = pt;
 572     containerInput[1] = track->Eta();
 573     containerInput[2] = track->Phi();
 574     containerInput[3] = track->Theta();
 575
 576     //
 577     // Step 4. Track that are recostructed, filling
 578     //
 579     fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed);
 580
 581     //
 582     // Step 5. Track that are recostructed and pass acceptance cuts, filling
 583     //
 584     if (!fMCCuts->IsRecoParticleKineAcceptance(track)) continue;
 585     fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoKineCuts);
 586
 587     // is track associated to particle ? if yes + implimenting the physical primary..
 588     Int_t label = TMath::Abs(track->GetLabel());
 589     if (label<=0) {
 590       AliDebug(3,"Particle not matching MC label \n");
 591       continue;
 592     }
 593
 594     //
 595     // Step 6-7. Track that are recostructed + Kine + Quality criteria, filling
 596     //
 597
 598     // check particle selections at MC level
 599     AliVParticle *mcPart  = (AliVParticle*)fMCEvent->GetTrack(label);
 600     if(!mcPart) continue;
 601     containerInputMC[0] = (Float_t)mcPart->Pt();
 602     containerInputMC[1] = mcPart->Eta();
 603     containerInputMC[2] = mcPart->Phi();
 604     containerInputMC[3] = mcPart->Theta();
 605
 606     if (!fMCCuts->IsMCParticleGenerated(mcPart)) continue;
 607     //    cout<< "MC matching did work"<<endl;
 608
 609
 610     // for filter bit selection
 611     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
 612     if(isAOD && !aodTrack) continue;
 613     if(isAOD && fSetFilterBit) if (!aodTrack->TestFilterMask(BIT(fbit))) continue;
 614     //    cout<<" Filter bit check passed"<<endl;
 615
 616     Bool_t isESDtrack = track->IsA()->InheritsFrom("AliESDtrack");
 617     AliESDtrack *tmptrack = NULL;
 618     if(isESDtrack) {
 619       tmptrack = dynamic_cast<AliESDtrack*>(track); // ESDs
 620     } else {
 621       if (aodTrack) tmptrack = ConvertTrack(aodTrack); // AODs
 622     }
 623     if (!tmptrack) continue;
 624
 625     // exclude global constrained and TPC only tracks (filter bits 128 and 512)
 626     Int_t id = tmptrack->GetID();
 627     if(isAOD && id<0) {
 628       AliDebug(3,"Track removed bc corresponts to either filter bit 128 or 512 (TPC only tracks)\n");
 629       continue;
 630     }
 631
 632     // Apply ESD track cuts
 633     if( !fTrackCuts->IsSelected(tmptrack) ){
 634       AliDebug(3,"Reconstructed track not passing quality criteria\n");
 635       if(isAOD) delete tmptrack;
 636       continue;
 637     }
 638     //    cout<<" analysis cuts passed"<<endl;
 639     fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
 640     fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoQualityCuts);
 641     if(isAOD) delete tmptrack;
 642     //    cout << " checking particle pid"<<endl;
 643
 644     //
 645     // Step8, PID check
 646     //
 647     if( !fMCCuts->IsRecoParticlePID(track) ){
 648       AliDebug(3,"Reconstructed track not passing PID criteria\n");
 649       continue;
 650     }
 651     fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID);
 652     //    cout << " all steps filled"<<endl;
 653
 654   }
 655   return;
 656 }
 657
 658 //______________________________________________________________________
 659 Int_t AliCFSingleTrackEfficiencyTask::GetNumberOfTrackletsInEtaRange(Double_t mineta, Double_t maxeta)
 660 {
 661   //
 662   // counts tracklets in given eta range
 663   //
 664
 665   AliAODEvent* event = dynamic_cast<AliAODEvent*>(fInputEvent);
 666   Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
 667   if(!event && AODEvent() && IsStandardAOD()) {
 668    event  = dynamic_cast<AliAODEvent*> (AODEvent());
 669   }
 670   if (!event) return -1;
 671   Int_t count=0;
 672
 673   if(isAOD) {
 674     AliAODTracklets* tracklets = event->GetTracklets();
 675     Int_t nTr=tracklets->GetNumberOfTracklets();
 676     for(Int_t iTr=0; iTr<nTr; iTr++){
 677       Double_t theta=tracklets->GetTheta(iTr);
 678       Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
 679       if(eta>mineta && eta<maxeta) count++;
 680     }
 681   } else {
 682     AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
 683     if (!esdEvent) return -1;
 684     const AliMultiplicity *mult = esdEvent->GetMultiplicity();
 685     if (mult) {
 686       if (mult->GetNumberOfTracklets()>0) {
 687         for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
 688           Double_t eta = -TMath::Log( TMath::Tan( mult->GetTheta(n) / 2.0 ) );
 689           if(eta>mineta && eta<maxeta) count++;
 690         }
 691       }
 692     }
 693   }
 694
 695   return count;
 696 }
 697
 698 //______________________________________________________________________
 699 Double_t AliCFSingleTrackEfficiencyTask::GetCentrality()
 700 {
 701   //
 702   // Get centrality
 703   //
 704
 705   Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
 706   Double_t cent = -1;
 707
 708   if(isAOD) {
 709     AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 710     if(!aodEvent) return cent;
 711     AliAODHeader* header = aodEvent->GetHeader();
 712     if(!header) return cent;
 713     AliCentrality *centrality = header->GetCentralityP();
 714     if(!centrality) return cent;
 715     //    cout<<" about to get cent perc AOD"<<endl;
 716     cent = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
 717   } else {
 718     AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
 719     if(!esdEvent) return cent;
 720     AliCentrality *centrality = esdEvent->GetCentrality();
 721     if(!centrality) return cent;
 722     cent = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
 723   }
 724
 725   return cent;
 726 }