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