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