PWGHF/hfe/AliAnalysisTaskDisplacedElectrons.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 // The analysis task:
  17 // study displaced electrons from beauty and charm
  18 // with cut on impact parameters in various pT bins
  19 //
  20 //
  21 // Authors:
  22 //  Hongyan Yang <hongyan@physi.uni-heidelberg.de>
  23 //  Carlo Bombonati <Carlo.Bombonati@cern.ch>
  24 //
  25
  26 #include <TChain.h>
  27 #include <TFile.h>
  28 #include <TH1F.h>
  29 #include <TH1I.h>
  30 #include "AliLog.h"
  31
  32 #include <TList.h>
  33 #include <TMath.h>
  34 #include <TObjArray.h>
  35 #include <TParticle.h>
  36 #include <TString.h>
  37 #include <TCanvas.h>
  38
  39 #include "AliAODEvent.h"
  40 #include "AliCFManager.h"
  41 #include "AliMCEvent.h"
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCParticle.h"
  44 #include "AliPIDResponse.h"
  45 #include "AliStack.h"
  46
  47 #include "AliESDEvent.h"
  48 #include "AliESDInputHandler.h"
  49 #include "AliESDtrack.h"
  50
  51 #include "AliAnalysisManager.h"
  52
  53 #include "AliHFEpid.h"
  54 #include "AliHFEcuts.h"
  55 #include "AliHFEtools.h"
  56 #include "AliHFEdisplacedElectrons.h"
  57 #include "AliAnalysisTaskDisplacedElectrons.h"
  58
  59
  60 //____________________________________________________________
  61 AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons():
  62   AliAnalysisTaskSE("Task for displaced electron study")
  63   , fDeDebugLevel(0)
  64   , fNminITSCluster(0)
  65   , fNminPrimVtxContrib(0)
  66   , fDePIDdetectors("")
  67   , fDePIDstrategy(0)
  68   , fDePlugins(0)
  69   , fDeCuts(0x0)
  70   , fDePID(0x0)
  71   , fDeCFM(0x0)
  72   , fDisplacedElectrons(0x0)
  73   , fDeNEvents(0x0)
  74   , fElectronsMcPt(0x0)
  75   , fElectronsEsdPt(0x0)
  76   , fElectronsDataPt(0x0)
  77   , fDeCorrection(0x0)
  78   , fDeQA(0x0)
  79   , fHistDisplacedElectrons(0x0)
  80 {
  81   //
  82   // Dummy constructor
  83   //
  84 }
  85
  86 //____________________________________________________________
  87 AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons(const char * name):
  88   AliAnalysisTaskSE(name)
  89   , fDeDebugLevel(0)
  90   , fNminITSCluster(0)
  91   , fNminPrimVtxContrib(0)
  92   , fDePIDdetectors("")
  93   , fDePIDstrategy(0)
  94   , fDePlugins(0)
  95   , fDeCuts(0x0)
  96   , fDePID(0x0)
  97   , fDeCFM(0x0)
  98   , fDisplacedElectrons(0x0)
  99   , fDeNEvents(0x0)
 100   , fElectronsMcPt(0x0)
 101   , fElectronsEsdPt(0x0)
 102   , fElectronsDataPt(0x0)
 103   , fDeCorrection(0x0)
 104   , fDeQA(0x0)
 105   , fHistDisplacedElectrons(0x0)
 106 {
 107   //
 108   // Default constructor
 109   //
 110   DefineInput(0, TChain::Class());
 111   DefineOutput(1, TList::Class());
 112   DefineOutput(2, TList::Class());
 113   DefineOutput(3, TList::Class());
 114
 115   // Initialize pid
 116   fDePID = new AliHFEpid("DEPID");
 117
 118 }
 119
 120 //____________________________________________________________
 121 AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons(const AliAnalysisTaskDisplacedElectrons &ref):
 122   AliAnalysisTaskSE(ref)
 123   , fDeDebugLevel(ref.fDeDebugLevel)
 124   , fNminITSCluster(ref.fNminITSCluster)
 125   , fNminPrimVtxContrib(ref.fNminPrimVtxContrib)
 126   , fDePIDdetectors(ref.fDePIDdetectors)
 127   , fDePIDstrategy(ref.fDePIDstrategy)
 128   , fDePlugins(ref.fDePlugins)
 129   , fDeCuts(ref.fDeCuts)
 130   , fDePID(ref.fDePID)
 131   , fDeCFM(ref.fDeCFM)
 132   , fDisplacedElectrons(ref.fDisplacedElectrons)
 133   , fDeNEvents(ref.fDeNEvents)
 134   , fElectronsMcPt(ref.fElectronsMcPt)
 135   , fElectronsEsdPt(ref.fElectronsEsdPt)
 136   , fElectronsDataPt(ref.fElectronsDataPt)
 137   , fDeCorrection(ref.fDeCorrection)
 138   , fDeQA(ref.fDeQA)
 139   , fHistDisplacedElectrons(ref.fHistDisplacedElectrons)
 140 {
 141   //
 142   // Copy Constructor
 143   //
 144 }
 145
 146 //____________________________________________________________
 147 AliAnalysisTaskDisplacedElectrons &AliAnalysisTaskDisplacedElectrons::operator=(const AliAnalysisTaskDisplacedElectrons &ref){
 148   //
 149   // Assignment operator
 150   //
 151   if(this == &ref) return *this;
 152   AliAnalysisTask::operator=(ref);
 153   fDeDebugLevel = ref.fDeDebugLevel;
 154   fNminITSCluster = ref.fNminITSCluster;
 155   fNminPrimVtxContrib = ref.fNminPrimVtxContrib;
 156   fDePIDdetectors = ref.fDePIDdetectors;
 157   fDePIDstrategy = ref.fDePIDstrategy;
 158   fDePlugins = ref.fDePlugins;
 159   fDePID = ref.fDePID;
 160   fDeCuts = ref.fDeCuts;
 161   fDeCFM = ref.fDeCFM;
 162   fDisplacedElectrons = ref.fDisplacedElectrons;
 163   fDeNEvents = ref.fDeNEvents;
 164   fElectronsMcPt = ref.fElectronsMcPt;
 165   fElectronsEsdPt = ref.fElectronsEsdPt;
 166   fElectronsDataPt = ref.fElectronsDataPt;
 167   fDeCorrection = ref.fDeCorrection;
 168   fDeQA = ref.fDeQA;
 169   fHistDisplacedElectrons = ref.fHistDisplacedElectrons;
 170
 171   return *this;
 172 }
 173
 174 //____________________________________________________________
 175 AliAnalysisTaskDisplacedElectrons::~AliAnalysisTaskDisplacedElectrons(){
 176   //
 177   // Destructor
 178   //
 179
 180   if(fDePID) delete fDePID;
 181   if(fDeCFM) delete fDeCFM;
 182   if(fDisplacedElectrons) delete fDisplacedElectrons;
 183   if(fDeNEvents) delete fDeNEvents;
 184   if(fElectronsMcPt) delete fElectronsMcPt;
 185   if(fElectronsEsdPt) delete fElectronsEsdPt;
 186   if(fElectronsDataPt) delete fElectronsDataPt;
 187   if(fDeCorrection){
 188     fDeCorrection->Clear();
 189     delete fDeCorrection;
 190   }
 191   if(fDeQA){
 192     fDeQA->Clear();
 193     delete fDeQA;
 194   }
 195   if(fHistDisplacedElectrons){
 196     fHistDisplacedElectrons->Clear();
 197     delete fHistDisplacedElectrons;
 198   }
 199 }
 200
 201 //____________________________________________________________
 202 void AliAnalysisTaskDisplacedElectrons::UserCreateOutputObjects(){
 203   // create output objects
 204   // fDeNEvents
 205   // MC and Data containers
 206
 207
 208   if(!fDeQA) fDeQA = new TList;
 209   fDeQA->SetName("variousQAhistograms");
 210
 211   fDeNEvents = new TH1I("nDeEvents", "Number of Events in the DE Analysis", 2, 0, 2);
 212   const Int_t nBins = 14;
 213   const Float_t ptBins[nBins] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,4.0,5.0,7.0,9.0,12.0,16.0,20.0};
 214   fElectronsMcPt = new TH1F("mcElectronPt", "MC: p_{T} distribution of identified electrons (mcpid);p_{T} (GeV/c);Counts;", nBins-1, ptBins);
 215   fElectronsEsdPt = new TH1F("esdElectronPt", "ESD: p_{T} distribution of identified electrons (hfepid);p_{T} (GeV/c);Counts;", nBins-1, ptBins);
 216   fElectronsDataPt = new TH1F("dataElectronPt", "DATA: p_{T} distribution of identified electrons (hfepid);p_{T} (GeV/c);Counts;", nBins-1, ptBins);
 217
 218   fDeQA->AddAt(fDeNEvents,0);
 219   if(HasMCData()){
 220     fDeQA->AddAt(fElectronsMcPt, 1);
 221     fDeQA->AddAt(fElectronsEsdPt, 2);
 222   }
 223   else{
 224     fDeQA->AddAt(fElectronsDataPt, 1);
 225   }
 226   // Initialize correction Framework and Cuts
 227   fDeCFM = new AliCFManager;
 228   MakeEventContainer();
 229   MakeParticleContainer();
 230
 231   if(!fDeCorrection) fDeCorrection = new TList();
 232   fDeCorrection->SetName("deCorrections");
 233   fDeCorrection->AddAt(fDeCFM->GetEventContainer(), 0);
 234   fDeCorrection->AddAt(fDeCFM->GetParticleContainer(), 1);
 235   fDeCorrection->Print();
 236
 237   for(Int_t istep = 0; istep < fDeCFM->GetEventContainer()->GetNStep(); istep++)
 238     fDeCFM->SetEventCutsList(istep, 0x0);
 239   for(Int_t istep = 0; istep < fDeCFM->GetParticleContainer()->GetNStep(); istep++)
 240     fDeCFM->SetParticleCutsList(istep, 0x0);
 241
 242   if(!fDeCuts){
 243     AliWarning("Cuts not available. Default cuts will be used");
 244     fDeCuts = new AliHFEcuts;
 245     fDeCuts->CreateStandardCuts();
 246   }
 247
 248   fDeCuts->Initialize(fDeCFM);
 249
 250   fDePID->SetHasMCData(HasMCData());
 251   fDePID->AddDetector("TPC", 0);
 252   fDePID->AddDetector("TOF", 1);
 253         fDePID->ConfigureTOF();
 254   fDePID->ConfigureTPCdefaultCut();
 255   fDePID->InitializePID();     // Only restrictions to TPC allowed
 256
 257
 258   // displaced electron study----------------------------------
 259   if(GetPlugin(kDisplacedElectrons)){
 260
 261     fDisplacedElectrons = new AliHFEdisplacedElectrons;
 262     fDisplacedElectrons->SetDebugLevel(fDeDebugLevel);
 263     fDisplacedElectrons->SetHasMCData(HasMCData());
 264     fDisplacedElectrons->SetMinPrimVtxContrib(fNminPrimVtxContrib);
 265     fDisplacedElectrons->SetNitsCluster(fNminITSCluster);
 266
 267     if(!fHistDisplacedElectrons) fHistDisplacedElectrons = new TList();
 268     fDisplacedElectrons->CreateOutputs(fHistDisplacedElectrons);
 269   }
 270
 271 }
 272
 273
 274
 275 //____________________________________________________________
 276 void AliAnalysisTaskDisplacedElectrons::UserExec(Option_t *){
 277   //
 278   // Run the analysis
 279   //
 280
 281   if(fDeDebugLevel>=10) AliInfo("analyse single event");
 282
 283   if(!fInputEvent){
 284     AliError("Reconstructed Event not available");
 285     return;
 286   }
 287
 288   //
 289   AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
 290
 291   if(!inH){
 292     AliError("AliESDInputHandler dynamic cast failed!");
 293     return;
 294   }
 295
 296   // from now on, only ESD are analyzed
 297   // using HFE pid, using HFE cuts
 298   // using CORRFW
 299   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 300   if(!pidResponse){
 301     AliDebug(1, "Using default PID Response");
 302     pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliAODEvent::Class());
 303   }
 304   fDePID->SetPIDResponse(pidResponse);
 305
 306   if(!fDeCuts){
 307     AliError("HFE cuts not available");
 308     return;
 309   }
 310
 311   // ---- CHOICE OF ANALYSIS ----
 312   if(HasMCData()){
 313     // Protect against missing MC trees
 314     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 315     if(!mcH){
 316       AliError("MC handler cuts not available");
 317       return;
 318     }
 319
 320     if(!mcH->InitOk()) return;
 321     if(!mcH->TreeK()) return;
 322     if(!mcH->TreeTR()) return;
 323
 324     AliDebug(4, Form("MC Event: %p", fMCEvent));
 325     if(!fMCEvent){
 326       AliError("No MC Event, but MC Data required");
 327       return;
 328     }
 329
 330     ProcessMC(); // PURE MC
 331
 332     if(IsESDanalysis()) ProcessESD(); // ESD WITH MC
 333
 334   }else if(IsESDanalysis()) ProcessData(); // PURE ESD
 335
 336
 337   fDeNEvents->Fill(1);
 338   PostData(1, fHistDisplacedElectrons);
 339   PostData(2, fDeCorrection);
 340   PostData(3, fDeQA);
 341
 342 }
 343
 344 //____________________________________________________________
 345 void AliAnalysisTaskDisplacedElectrons::ProcessMC(){
 346   //
 347   // handle pure MC analysis
 348   //
 349
 350   Int_t nMCelectrons = 0;
 351   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 352   if(!fESD){
 353     AliError("ESD event not available");
 354     return;
 355   }
 356
 357   Double_t mcContainer[4];   // container for the output in THnSparse
 358   memset(mcContainer, 0, sizeof(Double_t) * 4);
 359
 360   fDeCFM->SetMCEventInfo(fMCEvent);
 361
 362   Double_t nContributor[1] = {0};
 363   const AliVVertex *mcPrimVtx = fMCEvent->GetPrimaryVertex();
 364   if(mcPrimVtx) nContributor[0] = mcPrimVtx->GetNContributors();
 365
 366   //
 367   // cut at MC event level
 368   //
 369
 370   if(!fDeCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent)) return;
 371   if(GetPlugin(kCorrection)) fDeCFM->GetEventContainer()->Fill(nContributor,AliHFEcuts::kEventStepGenerated);
 372
 373   AliStack *stack = 0x0;
 374
 375   if(!fMCEvent->Stack())return;
 376   stack = fMCEvent->Stack();
 377   Int_t nTracks = stack->GetNtrack();
 378
 379   AliMCParticle *mcTrack = 0x0;
 380
 381   for(Int_t itrack = 0; itrack<nTracks; itrack++){
 382     if(!(stack->Particle(itrack))) continue;
 383     if(mcTrack) mcTrack = 0x0;
 384     mcTrack = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(itrack));
 385     if(!mcTrack) continue;
 386     //TParticle *mcPart = stack->Particle(itrack);
 387
 388     mcContainer[0] = mcTrack->Pt();
 389     mcContainer[1] = mcTrack->Eta();
 390     mcContainer[2] = mcTrack->Phi();
 391     mcContainer[3] = mcTrack->Charge();
 392
 393
 394
 395     if (!stack->IsPhysicalPrimary(mcTrack->GetLabel())) continue;
 396     // no cut but require primary
 397     if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(mcContainer, 0);
 398
 399     // all pions for reference
 400     if(TMath::Abs(mcTrack->Particle()->GetPdgCode())==AliHFEdisplacedElectrons::kPDGpion && GetPlugin(kCorrection))
 401       fDeCFM->GetParticleContainer()->Fill(mcContainer, 1);
 402
 403     // cut for signal: all MC electrons
 404     if(TMath::Abs(mcTrack->Particle()->GetPdgCode())==AliHFEdisplacedElectrons::kPDGelectron && GetPlugin(kCorrection))
 405       fDeCFM->GetParticleContainer()->Fill(mcContainer, 2);
 406
 407     // cut at track level kinematics: pt and eta
 408     if(TMath::Abs(mcContainer[1])>=0.8 || mcContainer[0]>20 || mcContainer[0]<0.1) continue;
 409
 410     if(TMath::Abs(mcTrack->Particle()->GetPdgCode())==AliHFEdisplacedElectrons::kPDGelectron){
 411       nMCelectrons++;
 412       fElectronsMcPt->Fill(mcContainer[0]);
 413     }
 414
 415     if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(mcContainer, 3);
 416
 417
 418     // fill MC THnSparse
 419     fDisplacedElectrons->FillMcOutput(fESD, fMCEvent, mcTrack);
 420
 421   }  // mc track loop
 422
 423   if(fDeDebugLevel>=10) printf("there are %d electrons in this MC event", nMCelectrons);
 424
 425 }
 426
 427 //____________________________________________________________
 428 void AliAnalysisTaskDisplacedElectrons::ProcessESD(){
 429
 430   // this is to handle ESD tracks with MC information
 431   // MC pid is only used when HFE pid is implemented, for comparison
 432   // corrections are taken into account
 433
 434   // process data: ESD tracks with MC information
 435
 436   const Int_t kStepPID = AliHFEcuts::kStepHFEcutsTRD + 1;
 437
 438   Double_t esdContainer[4];   // container for the output in THnSparse
 439   memset(esdContainer, 0, sizeof(Double_t) * 4);
 440   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 441   if(!fESD){
 442     AliError("No ESD event available");
 443     return;
 444   }
 445
 446   fDeCFM->SetRecEventInfo(fESD);
 447   Double_t nContrib[1] = {fESD->GetPrimaryVertex()->GetNContributors()};
 448
 449   Bool_t alreadyseen = kFALSE;
 450   AliLabelContainer cont(fESD->GetNumberOfTracks());
 451
 452   Int_t nHFEelectrons = 0;
 453   AliESDtrack *track = 0x0;
 454   AliStack *stack = 0x0;
 455
 456   if(!(stack = fMCEvent->Stack())) return;
 457
 458   //
 459   // cut at ESD event level
 460   //
 461   if(!fDeCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
 462   if(GetPlugin(kCorrection)) fDeCFM->GetEventContainer()->Fill(nContrib, AliHFEcuts::kEventStepReconstructed);
 463
 464   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 465     track = fESD->GetTrack(itrack);
 466
 467     if(GetPlugin(kDisplacedElectrons)) {
 468
 469       esdContainer[0] = track->Pt();
 470       esdContainer[1] = track->Eta();
 471       esdContainer[2] = track->Phi();
 472       esdContainer[3] = track->Charge();
 473
 474       // before any cut
 475       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));
 476       cont.Append(TMath::Abs(track->GetLabel()));  // check double counting
 477       if(alreadyseen) continue;  // avoid double counting
 478       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepRecNoCut + AliHFEcuts::kNcutStepsMCTrack);
 479
 480       // 1st track cut
 481       // RecKine: ITSTPC cuts : ITS & TPC refit, covmatrix: (2, 2, 0.5, 0.5, 2); min_tpccls: 50, chi2_tpccls: 3.5
 482       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 483       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
 484
 485       // 2nd track cut
 486       // RecPrim: cut on track quality : DCA to vertex max: 3cm and 10cm; reject kink daughters
 487       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 488       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack);
 489
 490       // 3rd track cut
 491       // HFEcuts: ITS layers cuts: ITS pixel layer: kFirst, kSecond, kBoth, kNone or kAny
 492       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 493       if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack);
 494
 495       /*
 496       //  4th track cut
 497       // TRD: number of tracklets in TRD
 498       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 499       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
 500       */
 501
 502       // 5th track cut
 503       // track accepted, do PID
 504       // --> only electron candidate will be processed
 505       AliHFEpidObject hfetrack;
 506       hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 507       hfetrack.SetRecTrack(track);
 508       //if(HasMCData())hfetrack.SetMCTrack(mctrack);
 509
 510       if(!fDePID->IsSelected(&hfetrack)) continue;
 511       else if(fDeDebugLevel>=10)
 512         AliInfo("ESD info: this particle is identified as electron by HFEpid method \n");
 513       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+kStepPID + AliHFEcuts::kNcutStepsMCTrack);
 514
 515       // Fill Containers
 516       nHFEelectrons++;
 517       fElectronsEsdPt->Fill(esdContainer[0]);
 518       fDisplacedElectrons->FillEsdOutput(fESD, track, stack);
 519
 520     }  // displaced electron analysis on ESD with MC plugin
 521   } // track loop
 522
 523   if(fDeDebugLevel>=10) printf("there are %d HFE electrons in this ESD event", nHFEelectrons);
 524
 525 }
 526
 527
 528
 529 //____________________________________________________________
 530 void AliAnalysisTaskDisplacedElectrons::ProcessData(){
 531
 532   // this is a track loop over real data
 533   // no MC information at all
 534   // HFE pid is used
 535
 536   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
 537   //const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
 538   const Int_t kStepPID = AliHFEcuts::kStepHFEcutsTRD + 1;
 539
 540   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 541   if(!fESD){
 542     AliError("No ESD event available");
 543     return;
 544   }
 545
 546   Double_t dataContainer[4];   // container for the output in THnSparse
 547   memset(dataContainer, 0, sizeof(Double_t) * 4);
 548
 549   Bool_t alreadyseen = kFALSE;
 550   AliLabelContainer cont(fESD->GetNumberOfTracks());
 551
 552
 553   AliESDtrack *track = 0x0;
 554   Int_t nHFEelectrons= 0;
 555
 556   fDeCFM->SetRecEventInfo(fESD);
 557   Double_t nContrib[1] = {fESD->GetPrimaryVertex()->GetNContributors()};
 558   if(!fDeCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
 559   if(GetPlugin(kCorrection)) fDeCFM->GetEventContainer()->Fill(nContrib, AliHFEcuts::kEventStepReconstructed);
 560
 561
 562   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 563     track = fESD->GetTrack(itrack);
 564
 565     if(GetPlugin(kDisplacedElectrons)) {
 566
 567       dataContainer[0] = track->Pt();
 568       dataContainer[1] = track->Eta();
 569       dataContainer[2] = track->Phi();
 570       dataContainer[3] = track->Charge();
 571
 572       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));  // double counted track
 573       cont.Append(TMath::Abs(track->GetLabel()));
 574       if(alreadyseen) continue;  // avoid double counting
 575       if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(&dataContainer[4],
 576                                                                      1+AliHFEcuts::kStepRecNoCut + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
 577
 578       // 1st track cut
 579       // RecKine: ITSTPC cuts : ITS & TPC refit, covmatrix: (2, 2, 0.5, 0.5, 2); min_tpccls: 50, chi2_tpccls: 3.5
 580       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 581       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(&dataContainer[4],
 582                                                                       1+AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
 583
 584       // 2nd track cut
 585       // RecPrim: cut on track quality : DCA to vertex max: 3cm and 10cm; reject kink daughters
 586       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track))  continue;
 587       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(&dataContainer[4],
 588                                                                       1+AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
 589
 590       //  3rd track cut
 591       // HFEcuts: ITS layers cuts: ITS pixel layer: kFirst, kSecond, kBoth, kNone or kAny
 592       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 593       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(&dataContainer[4],
 594                                                                       1+AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
 595
 596       /*
 597       //  4th track cut
 598       // TRD: number of tracklets in TRD0
 599       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
 600       if(GetPlugin(kCorrection)) if(HasMCData())fDeCFM->GetParticleContainer()->Fill(&dataContainer[4],
 601                                                                                      1+AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
 602       */
 603
 604
 605       // 5th track cut
 606       // track accepted, do PID --> only electron candidate will be processed
 607
 608       AliHFEpidObject hfetrack;
 609       hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 610       hfetrack.SetRecTrack(track);
 611
 612       if(!fDePID->IsSelected(&hfetrack)) continue;
 613       else if(fDeDebugLevel>=10)
 614         AliInfo("ESD info: this particle is identified as electron by HFEpid method \n");
 615       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(dataContainer,  1+kStepPID + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
 616
 617       nHFEelectrons++;
 618       fElectronsDataPt->Fill(dataContainer[0]);
 619       fDisplacedElectrons->FillDataOutput(fESD, track);
 620     } // analyze displaced electrons plugin switched on
 621   } // track loop
 622
 623   if(fDeDebugLevel>=10) printf("there are %d HFE electrons in this DATA event", nHFEelectrons);
 624 }
 625
 626
 627 //____________________________________________________________
 628 void AliAnalysisTaskDisplacedElectrons::Terminate(Option_t *){
 629   //
 630   // Terminate not implemented at the moment
 631   //
 632
 633   fHistDisplacedElectrons = dynamic_cast<TList *>(GetOutputData(1));
 634   fDeCorrection = dynamic_cast<TList *>(GetOutputData(2));
 635   fDeQA = dynamic_cast<TList *>(GetOutputData(3));
 636   if(!fDeCorrection) AliError("correction list not available\n");
 637   if(!fHistDisplacedElectrons) AliError("de list not available\n");
 638   if(!fDeQA) AliError("qa list is not available\n");
 639
 640   fHistDisplacedElectrons->Print();
 641   fDeCorrection->Print();
 642   fDeQA->Print();
 643
 644   AliInfo("analysis done!\n");
 645
 646 }
 647
 648 //____________________________________________________________
 649 void AliAnalysisTaskDisplacedElectrons::PrintStatus() const {
 650
 651   //
 652   // Print Analysis status
 653   //
 654   printf("\n");
 655   printf("\t Analysis Settings\n\t========================================\n");
 656   printf("\t running over %s\n", HasMCData()?"MC data":"pp collision data");
 657   printf("\t displaced electrons' analysis is %s\n", GetPlugin(kDisplacedElectrons)?"ON":"OFF");
 658   printf("\t correction container is %s\n", GetPlugin(kCorrection)?"ON":"OFF");
 659   printf("\t hfe pid qa is %s\n", GetPlugin(kDePidQA)?"ON":"OFF");
 660   printf("\t post processing  is %s\n", GetPlugin(kPostProcess)?"ON":"OFF");
 661   printf("\t cuts: %s\n", (fDeCuts != NULL) ? "YES" : "NO");
 662   printf("\t ");
 663   printf("\n");
 664 }
 665
 666 //__________________________________________
 667 void AliAnalysisTaskDisplacedElectrons::SwitchOnPlugin(Int_t plug){
 668   //
 669   // Switch on Plugin
 670   // Available:
 671   //  - analyze impact parameter
 672   //  - Post Processing
 673
 674   switch(plug)
 675     {
 676     case kDisplacedElectrons:
 677       SETBIT(fDePlugins, plug);
 678       break;
 679     case kCorrection:
 680       SETBIT(fDePlugins, plug);
 681       break;
 682     case kDePidQA:
 683       SETBIT(fDePlugins, plug);
 684       break;
 685     case kPostProcess:
 686       SETBIT(fDePlugins, plug);
 687       break;
 688     default:
 689       AliError("Unknown Plugin");
 690     };
 691 }
 692
 693 //____________________________________________________________
 694 void AliAnalysisTaskDisplacedElectrons::MakeParticleContainer(){
 695   //
 696   // Create the particle container for the correction framework manager and
 697   // link it
 698   //
 699   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
 700   const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
 701
 702   const Int_t kNvar   = 4;
 703   //number of variables on the grid:pt,eta, phi, charge
 704   const Double_t kPtbound[2] = {0.1, 10.};
 705   const Double_t kEtabound[2] = {-0.8, 0.8};
 706   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
 707
 708   //arrays for the number of bins in each dimension
 709   Int_t iBin[kNvar];
 710   iBin[0] = 40; // bins in pt
 711   iBin[1] =  8; // bins in eta
 712   iBin[2] = 18; // bins in phi
 713   iBin[3] =  2; // bins in charge
 714
 715   //arrays for lower bounds :
 716   Double_t* binEdges[kNvar];
 717   binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
 718   binEdges[1] = AliHFEtools::MakeLinearBinning(iBin[1], kEtabound[0], kEtabound[1]);
 719   binEdges[2] = AliHFEtools::MakeLinearBinning(iBin[2], kPhibound[0], kPhibound[1]);
 720   binEdges[3] = AliHFEtools::MakeLinearBinning(iBin[3], -1.1, 1.1); // Numeric precision
 721
 722   //------------------------------------------------
 723   //     one "container" for MC+ESD+Data
 724   //----------pure MC track-------------------------
 725   // 0: MC generated
 726   // 1: MC pion total  ---- be careful!!!!
 727   // 2: MC electrons total
 728   // 3: MC electrons in acceptance
 729   //-------ESD track with MC info-------------------
 730   // 4: ESD track with MC: no cut
 731   // 5: ESD track with MC: cut on kine its tpc
 732   // 6: ESD track with MC: rec prim
 733   // 7: ESD track with MC: hfe cuts its
 734   // 8: ESD track with MC: hfe cuts trd
 735   // 9: ESD track with MC: hfe pid
 736   //-----------data track---------------------------
 737   // 10: DATA track wo MC: no cut
 738   // 11: DATA track wo MC: cut on kine its tpc
 739   // 12: DATA track wo MC: rec prim
 740   // 13: DATA track wo MC: hfe cuts its
 741   // 14: DATA track wo MC: hfe cuts trd
 742   // 15: DATA track wo MC: hfe pid
 743   //------------------------------------------------
 744
 745   AliCFContainer* container = new AliCFContainer("deTrackContainer", "Container for tracks",
 746                                                  (1 + kNcutStepsTrack + kNcutStepsESDtrack), kNvar, iBin);
 747
 748   //setting the bin limits
 749   for(Int_t ivar = 0; ivar < kNvar; ivar++){
 750     container -> SetBinLimits(ivar, binEdges[ivar]);
 751   }
 752   fDeCFM->SetParticleContainer(container);
 753 }
 754
 755
 756 //____________________________________________________________
 757 void AliAnalysisTaskDisplacedElectrons::MakeEventContainer(){
 758   //
 759   // Create the event container for the correction framework and link it
 760   //
 761
 762   // event container
 763   // 0: MC event
 764   // 1: ESD event
 765
 766   const Int_t kNvar = 1;  // number of variables on the grid: number of tracks per event
 767   const Double_t kNTrackBound[2] = {-0.5, 200.5};
 768   const Int_t kNBins = 201;
 769
 770   AliCFContainer *evCont = new AliCFContainer("deEventContainer", "Container for DE events", AliHFEcuts::kNcutStepsEvent, kNvar, &kNBins);
 771
 772   Double_t trackBins[kNBins];
 773   for(Int_t ibin = 0; ibin < kNBins; ibin++) trackBins[ibin] = kNTrackBound[0] + static_cast<Double_t>(ibin);
 774   evCont->SetBinLimits(0,trackBins);
 775
 776   fDeCFM->SetEventContainer(evCont);
 777
 778 }
 779
 780
 781
 782 //__________________________________________________________
 783 void AliAnalysisTaskDisplacedElectrons::AddPIDdetector(TString detector){
 784   //
 785   // Adding PID detector to the task
 786   //
 787   if(!fDePIDdetectors.Length())
 788     fDePIDdetectors = detector;
 789   else
 790     fDePIDdetectors += ":" + detector;
 791 }
 792
 793
 794
 795 //____________________________________________________________
 796 AliAnalysisTaskDisplacedElectrons::AliLabelContainer::AliLabelContainer(Int_t capacity):
 797   fContainer(NULL),
 798   fBegin(NULL),
 799   fEnd(NULL),
 800   fLast(NULL),
 801   fCurrent(NULL)
 802 {
 803   //
 804   // Default constructor
 805   //
 806   fContainer = new Int_t[capacity];
 807   fBegin = &fContainer[0];
 808   fEnd = &fContainer[capacity - 1];
 809   fLast = fCurrent = fBegin;
 810 }
 811
 812 //____________________________________________________________
 813 Bool_t AliAnalysisTaskDisplacedElectrons::AliLabelContainer::Append(Int_t label){
 814   //
 815   // Add Label to the container
 816   //
 817   if(fLast > fEnd) return kFALSE;
 818   *fLast++ = label;
 819   return kTRUE;
 820 }
 821
 822 //____________________________________________________________
 823 Bool_t AliAnalysisTaskDisplacedElectrons::AliLabelContainer::Find(Int_t label) const {
 824   //
 825   // Find track in the list of labels
 826   //
 827   for(Int_t *entry = fBegin; entry <= fLast; entry++)
 828     if(*entry == label) return kTRUE;
 829   return kFALSE;
 830 }
 831
 832 //____________________________________________________________
 833 Int_t AliAnalysisTaskDisplacedElectrons::AliLabelContainer::Next() {
 834   //
 835   // Mimic iterator
 836   //
 837   if(fCurrent > fLast) return -1;
 838   fCurrent++;
 839   return *fCurrent;
 840 }