PWGHF/hfe/AliAnalysisTaskHFE.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 // Filling an AliCFContainer with the quantities pt, eta and phi
  18 // for tracks which survivied the particle cuts (MC resp. ESD tracks)
  19 // Track selection is done using the AliHFE package
  20 //
  21 // Author:
  22 //  Raphaelle Bailhache <R.Bailhache@gsi.de>
  23 //  Markus Fasel <M.Fasel@gsi.de>
  24 //  Matus Kalisky <matus.kalisky@cern.ch>
  25 //  MinJung Kweon <minjung@physi.uni-heidelberg.de>
  26 //
  27 #include <TAxis.h>
  28 #include <TBits.h>
  29 #include <TCanvas.h>
  30 #include <TChain.h>
  31 #include <TDirectory.h>
  32 #include <TFile.h>
  33 #include <TH3D.h>
  34 #include <TIterator.h>
  35 #include <TList.h>
  36 #include <TLegend.h>
  37 #include <TMath.h>
  38 #include <TObjArray.h>
  39 #include <TObjString.h>
  40 #include <TParticle.h>
  41 #include <TProfile.h>
  42 #include <TString.h>
  43 #include <TF1.h>
  44 #include <TTree.h>
  45
  46 #include "AliAnalysisManager.h"
  47 #include "AliAODInputHandler.h"
  48 #include "AliAODMCParticle.h"
  49 #include "AliAODTrack.h"
  50 #include "AliAODVertex.h"
  51 #include "AliCentrality.h"
  52 #include "AliCFContainer.h"
  53 #include "AliCFManager.h"
  54 #include "AliESDEvent.h"
  55 #include "AliESDInputHandler.h"
  56 #include "AliESDtrack.h"
  57 #include "AliLog.h"
  58 #include "AliMCEvent.h"
  59 #include "AliMCEventHandler.h"
  60 #include "AliMCParticle.h"
  61 #include "AliMultiplicity.h"
  62 #include "AliPID.h"
  63 #include "AliPIDResponse.h"
  64 #include "AliOADBContainer.h"
  65 #include "AliStack.h"
  66 #include "AliTriggerAnalysis.h"
  67 #include "AliVVertex.h"
  68
  69 #include "AliHFEcollection.h"
  70 #include "AliHFEcontainer.h"
  71 #include "AliHFEcuts.h"
  72 #include "AliHFEelecbackground.h"
  73 #include "AliHFEmcQA.h"
  74 #include "AliHFEpairs.h"
  75 #include "AliHFEpid.h"
  76 #include "AliHFEpidQAmanager.h"
  77 #include "AliHFEpostAnalysis.h"
  78 #include "AliHFEsecVtxs.h"
  79 #include "AliHFEsecVtx.h"
  80 #include "AliHFEsignalCuts.h"
  81 #include "AliHFEtaggedTrackAnalysis.h"
  82 #include "AliHFEtools.h"
  83 #include "AliHFEvarManager.h"
  84 #include "AliAnalysisTaskHFE.h"
  85
  86 ClassImp(AliAnalysisTaskHFE)
  87
  88 //____________________________________________________________
  89 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  90   AliAnalysisTaskSE("PID efficiency Analysis")
  91   , fQAlevel(0)
  92   , fPlugins(0)
  93   , fFillSignalOnly(kTRUE)
  94   , fFillNoCuts(kFALSE)
  95   , fUseFlagAOD(kTRUE)
  96   , fApplyCutAOD(kFALSE)
  97   , fFlags(1<<4)
  98   , fBackGroundFactorApply(kFALSE)
  99   , fRemovePileUp(kFALSE)
 100   , fIdentifiedAsPileUp(kFALSE)
 101   , fIdentifiedAsOutInz(kFALSE)
 102   , fPassTheEventCut(kFALSE)
 103   , fRejectKinkMother(kTRUE)
 104   , fisppMultiBin(kFALSE)
 105   , fPbPbUserCentralityBinning(kFALSE)
 106   , fisNonHFEsystematics(kFALSE)
 107   , fSpecialTrigger(NULL)
 108   , fCentralityF(-1)
 109   , fCentralityPercent(-1)
 110   , fContributors(0.5)
 111   , fWeightBackGround(0.)
 112   , fVz(0.0)
 113   , fHadronBackgroundOADB(NULL)
 114   , fContainer(NULL)
 115   , fVarManager(NULL)
 116   , fSignalCuts(NULL)
 117   , fCFM(NULL)
 118   , fTriggerAnalysis(NULL)
 119   , fPID(NULL)
 120   , fPIDqa(NULL)
 121   , fPIDpreselect(NULL)
 122   , fCuts(NULL)
 123   , fTaggedTrackCuts(NULL)
 124   , fCleanTaggedTrack(kFALSE)
 125   , fVariablesTRDTaggedTrack(kFALSE)
 126   , fCutspreselect(NULL)
 127   , fSecVtx(NULL)
 128   , fElecBackGround(NULL)
 129   , fMCQA(NULL)
 130   , fTaggedTrackAnalysis(NULL)
 131   , fExtraCuts(NULL)
 132   , fQA(NULL)
 133   , fOutput(NULL)
 134   , fHistMCQA(NULL)
 135   , fHistSECVTX(NULL)
 136   , fHistELECBACKGROUND(NULL)
 137   , fQACollection(NULL)
 138 {
 139   //
 140   // Dummy constructor
 141   //
 142   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 143   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
 144   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 145   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
 146   memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
 147
 148
 149 }
 150
 151 //____________________________________________________________
 152 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
 153   AliAnalysisTaskSE(name)
 154   , fQAlevel(0)
 155   , fPlugins(0)
 156   , fFillSignalOnly(kTRUE)
 157   , fFillNoCuts(kFALSE)
 158   , fUseFlagAOD(kTRUE)
 159   , fApplyCutAOD(kFALSE)
 160   , fFlags(1<<4)
 161   , fBackGroundFactorApply(kFALSE)
 162   , fRemovePileUp(kFALSE)
 163   , fIdentifiedAsPileUp(kFALSE)
 164   , fIdentifiedAsOutInz(kFALSE)
 165   , fPassTheEventCut(kFALSE)
 166   , fRejectKinkMother(kTRUE)
 167   , fisppMultiBin(kFALSE)
 168   , fPbPbUserCentralityBinning(kFALSE)
 169   , fisNonHFEsystematics(kFALSE)
 170   , fSpecialTrigger(NULL)
 171   , fCentralityF(-1)
 172   , fCentralityPercent(-1)
 173   , fContributors(0.5)
 174   , fWeightBackGround(0.)
 175   , fVz(0.0)
 176   , fHadronBackgroundOADB(NULL)
 177   , fContainer(NULL)
 178   , fVarManager(NULL)
 179   , fSignalCuts(NULL)
 180   , fCFM(NULL)
 181   , fTriggerAnalysis(NULL)
 182   , fPID(NULL)
 183   , fPIDqa(NULL)
 184   , fPIDpreselect(NULL)
 185   , fCuts(NULL)
 186   , fTaggedTrackCuts(NULL)
 187   , fCleanTaggedTrack(kFALSE)
 188   , fVariablesTRDTaggedTrack(kFALSE)
 189   , fCutspreselect(NULL)
 190   , fSecVtx(NULL)
 191   , fElecBackGround(NULL)
 192   , fMCQA(NULL)
 193   , fTaggedTrackAnalysis(NULL)
 194   , fExtraCuts(NULL)
 195   , fQA(NULL)
 196   , fOutput(NULL)
 197   , fHistMCQA(NULL)
 198   , fHistSECVTX(NULL)
 199   , fHistELECBACKGROUND(NULL)
 200   , fQACollection(0x0)
 201 {
 202   //
 203   // Default constructor
 204   //
 205   DefineOutput(1, TList::Class());
 206   DefineOutput(2, TList::Class());
 207
 208   fPID = new AliHFEpid("hfePid");
 209   fPIDqa = new AliHFEpidQAmanager;
 210   fVarManager = new AliHFEvarManager("hfeVarManager");
 211
 212   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins * kCentBins * kBgLevels);
 213   memset(fkBackGroundFactorArray, 0, sizeof(TF1 *) * 12);
 214   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 215   memset(&fisppMultiBin, kFALSE, sizeof(fisppMultiBin));
 216   memset(fCentralityLimits, 0, sizeof(Float_t) * 12);
 217
 218 }
 219
 220 //____________________________________________________________
 221 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
 222   AliAnalysisTaskSE(ref)
 223   , fQAlevel(0)
 224   , fPlugins(0)
 225   , fFillSignalOnly(ref.fFillSignalOnly)
 226   , fFillNoCuts(ref.fFillNoCuts)
 227   , fUseFlagAOD(ref.fUseFlagAOD)
 228   , fApplyCutAOD(ref.fApplyCutAOD)
 229   , fFlags(ref.fFlags)
 230   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
 231   , fRemovePileUp(ref.fRemovePileUp)
 232   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
 233   , fIdentifiedAsOutInz(ref.fIdentifiedAsOutInz)
 234   , fPassTheEventCut(ref.fPassTheEventCut)
 235   , fRejectKinkMother(ref.fRejectKinkMother)
 236   , fisppMultiBin(ref.fisppMultiBin)
 237   , fPbPbUserCentralityBinning(ref.fPbPbUserCentralityBinning)
 238   , fisNonHFEsystematics(ref.fisNonHFEsystematics)
 239   , fSpecialTrigger(ref.fSpecialTrigger)
 240   , fCentralityF(ref.fCentralityF)
 241   , fCentralityPercent(ref.fCentralityPercent)
 242   , fContributors(ref.fContributors)
 243   , fWeightBackGround(ref.fWeightBackGround)
 244   , fVz(ref.fVz)
 245   , fHadronBackgroundOADB(ref.fHadronBackgroundOADB)
 246   , fContainer(NULL)
 247   , fVarManager(NULL)
 248   , fSignalCuts(NULL)
 249   , fCFM(NULL)
 250   , fTriggerAnalysis(NULL)
 251   , fPID(NULL)
 252   , fPIDqa(NULL)
 253   , fPIDpreselect(NULL)
 254   , fCuts(NULL)
 255   , fTaggedTrackCuts(NULL)
 256   , fCleanTaggedTrack(ref.fCleanTaggedTrack)
 257   , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
 258   , fCutspreselect(NULL)
 259   , fSecVtx(NULL)
 260   , fElecBackGround(NULL)
 261   , fMCQA(NULL)
 262   , fTaggedTrackAnalysis(NULL)
 263   , fExtraCuts(NULL)
 264   , fQA(NULL)
 265   , fOutput(NULL)
 266   , fHistMCQA(NULL)
 267   , fHistSECVTX(NULL)
 268   , fHistELECBACKGROUND(NULL)
 269   , fQACollection(NULL)
 270 {
 271   //
 272   // Copy Constructor
 273   //
 274   ref.Copy(*this);
 275 }
 276
 277 //____________________________________________________________
 278 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 279   //
 280   // Assignment operator
 281   //
 282   if(this == &ref)
 283     ref.Copy(*this);
 284   return *this;
 285 }
 286
 287 //____________________________________________________________
 288 void AliAnalysisTaskHFE::Copy(TObject &o) const {
 289   //
 290   // Copy into object o
 291   //
 292   AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
 293   target.fQAlevel = fQAlevel;
 294   target.fPlugins = fPlugins;
 295   target.fFillSignalOnly = fFillSignalOnly;
 296   target.fFillNoCuts = fFillNoCuts;
 297   target.fUseFlagAOD = fUseFlagAOD;
 298   target.fApplyCutAOD = fApplyCutAOD;
 299   target.fFlags = fFlags;
 300   target.fBackGroundFactorApply = fBackGroundFactorApply;
 301   target.fRemovePileUp = fRemovePileUp;
 302   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
 303   target.fIdentifiedAsOutInz = fIdentifiedAsOutInz;
 304   target.fPassTheEventCut = fPassTheEventCut;
 305   target.fRejectKinkMother = fRejectKinkMother;
 306   target.fisppMultiBin =   fisppMultiBin;
 307   target.fPbPbUserCentralityBinning = fPbPbUserCentralityBinning;
 308   target.fisNonHFEsystematics = fisNonHFEsystematics;
 309   target.fSpecialTrigger = fSpecialTrigger;
 310   target.fCentralityF = fCentralityF;
 311   target.fCentralityPercent = fCentralityPercent;
 312   target.fContributors = fContributors;
 313   target.fWeightBackGround = fWeightBackGround;
 314   target.fVz = fVz;
 315   target.fHadronBackgroundOADB = fHadronBackgroundOADB;
 316   target.fContainer = fContainer;
 317   target.fVarManager = fVarManager;
 318   target.fSignalCuts = fSignalCuts;
 319   target.fCFM = fCFM;
 320   target.fTriggerAnalysis = fTriggerAnalysis;
 321   target.fPID = fPID;
 322   target.fPIDqa = fPIDqa;
 323   target.fPIDpreselect = fPIDpreselect;
 324   target.fCuts = fCuts;
 325   target.fTaggedTrackCuts = fTaggedTrackCuts;
 326   target.fCleanTaggedTrack = fCleanTaggedTrack;
 327   target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
 328   target.fCutspreselect = fCutspreselect;
 329   target.fSecVtx = fSecVtx;
 330   target.fElecBackGround = fElecBackGround;
 331   target.fMCQA = fMCQA;
 332   target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
 333   target.fExtraCuts = fExtraCuts;
 334   target.fQA = fQA;
 335   target.fOutput = fOutput;
 336   target.fHistMCQA = fHistMCQA;
 337   target.fHistSECVTX = fHistSECVTX;
 338   target.fHistELECBACKGROUND = fHistELECBACKGROUND;
 339   target.fQACollection = fQACollection;
 340 }
 341
 342 //____________________________________________________________
 343 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 344   //
 345   // Destructor
 346   //
 347   if(fPID) delete fPID;
 348   if(fPIDpreselect) delete fPIDpreselect;
 349   if(fVarManager) delete fVarManager;
 350   if(fCFM) delete fCFM;
 351   if(fTriggerAnalysis) delete fTriggerAnalysis;
 352   if(fSignalCuts) delete fSignalCuts;
 353   if(fSecVtx) delete fSecVtx;
 354   if(fMCQA) delete fMCQA;
 355   if(fElecBackGround) delete fElecBackGround;
 356   if(fSpecialTrigger) delete fSpecialTrigger;
 357   // Delete output objects only if we are not running in PROOF mode because otherwise this produces a crash during merging
 358   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 359   if(mgr && mgr->GetAnalysisType() != AliAnalysisManager::kProofAnalysis){
 360     if(fPIDqa) delete fPIDqa;
 361     if(fOutput) delete fOutput;
 362     if(fQA) delete fQA;
 363   }
 364 }
 365
 366 //____________________________________________________________
 367 void AliAnalysisTaskHFE::UserCreateOutputObjects(){
 368   //
 369   // Creating output container and output objects
 370   // Here we also Initialize the correction framework container and
 371   // the objects for
 372   // - PID
 373   // - MC QA
 374   // - SecVtx
 375   // QA histograms are created if requested
 376   // Called once per worker
 377   //
 378   AliDebug(3, "Creating Output Objects");
 379   // Automatic determination of the analysis mode
 380   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 381   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 382     SetAODAnalysis();
 383   } else {
 384     SetESDAnalysis();
 385     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 386       SetHasMCData();
 387   }
 388   printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
 389   printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
 390
 391   // Enable Trigger Analysis
 392   fTriggerAnalysis = new AliTriggerAnalysis;
 393   fTriggerAnalysis->EnableHistograms();
 394   fTriggerAnalysis->SetAnalyzeMC(HasMCData());
 395
 396
 397   // Make lists for Output
 398   if(!fQA) fQA = new TList;
 399   fQA->SetOwner();
 400   if(!fOutput) fOutput = new TList;
 401   fOutput->SetOwner();
 402
 403   // First Part: Make QA histograms
 404   fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
 405   fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 406   fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
 407   fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
 408   fQACollection->CreateTH2F("TPCdEdxBeforePID", "TPC dE/dx; p (GeV/c); TPC dE/dx (a.u.)", 1000, 0., 10., 200, 0., 200.);
 409   fQACollection->CreateTH2F("TPCnSigmaBeforePID", "TPC dE/dx; p (GeV/c); TPC dE/dx - <TPC dE/dx>|_{el} (#sigma)", 1000, 0., 10., 1000, -10., 10.);
 410
 411   InitPIDperformanceQA();
 412   InitContaminationQA();
 413   InitHistoITScluster();
 414   fQA->Add(fQACollection);
 415
 416   // Initialize PID
 417   fPID->SetHasMCData(HasMCData());
 418   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
 419   if(IsQAOn(kPIDqa)){
 420     AliInfo("PID QA switched on");
 421     fPIDqa->Initialize(fPID);
 422     fQA->Add(fPIDqa->MakeList("HFEpidQA"));
 423   }
 424   fPID->SortDetectors();
 425
 426   // Initialize correction Framework and Cuts
 427   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
 428   fCFM = new AliCFManager;
 429   fCFM->SetNStepParticle(kNcutSteps);
 430   MakeParticleContainer();
 431   MakeEventContainer();
 432   // Temporary fix: Initialize particle cuts with NULL
 433   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 434     fCFM->SetParticleCutsList(istep, NULL);
 435   if(!fCuts){
 436     AliWarning("Cuts not available. Default cuts will be used");
 437     fCuts = new AliHFEcuts;
 438     fCuts->CreateStandardCuts();
 439   }
 440   if(IsAODanalysis()) fCuts->SetAOD();
 441   // Make clone for V0 tagging step
 442   fCuts->Initialize(fCFM);
 443   if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
 444   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
 445   fVarManager->SetSignalCuts(fSignalCuts);
 446
 447   // add output objects to the List
 448   fOutput->AddAt(fContainer, 0);
 449   fOutput->AddAt(fCFM->GetEventContainer(), 1);
 450
 451   // mcQA----------------------------------
 452   if (HasMCData() && IsQAOn(kMCqa)) {
 453     AliInfo("MC QA on");
 454     if(!fMCQA) fMCQA = new AliHFEmcQA;
 455     if(!fHistMCQA) fHistMCQA = new TList();
 456     fHistMCQA->SetOwner();
 457     if(IsPbPb()) fMCQA->SetPbPb();
 458     if(fisppMultiBin) fMCQA->SetPPMultiBin();
 459     if(TestBit(kTreeStream)){
 460       fMCQA->EnableDebugStreamer();
 461     }
 462     fMCQA->CreatDefaultHistograms(fHistMCQA);
 463     fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0][0][0],fBinLimit);
 464     fQA->Add(fHistMCQA);
 465   }
 466
 467   // secvtx----------------------------------
 468   if (GetPlugin(kSecVtx)) {
 469     AliInfo("Secondary Vertex Analysis on");
 470     if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
 471     fSecVtx->SetHasMCData(HasMCData());
 472
 473     if(!fHistSECVTX) fHistSECVTX = new TList();
 474     fHistSECVTX->SetOwner();
 475     fSecVtx->CreateHistograms(fHistSECVTX);
 476     fOutput->Add(fHistSECVTX);
 477   }
 478
 479   // background----------------------------------
 480   if (GetPlugin(kIsElecBackGround)) {
 481     AliInfo("Electron BackGround Analysis on");
 482     if(!fElecBackGround){
 483       AliWarning("ElecBackGround not available. Default elecbackground will be used");
 484       fElecBackGround = new AliHFEelecbackground;
 485     }
 486     fElecBackGround->SetHasMCData(HasMCData());
 487
 488     if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
 489     fHistELECBACKGROUND->SetOwner();
 490     fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
 491     fOutput->Add(fHistELECBACKGROUND);
 492   }
 493
 494   // tagged tracks
 495   if(GetPlugin(kTaggedTrackAnalysis)){
 496     AliInfo("Analysis on V0-tagged tracks enabled");
 497     fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis(Form("taggedTrackAnalysis%s", GetName()));
 498     fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
 499     fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
 500     AliHFEvarManager *varManager = fTaggedTrackAnalysis->GetVarManager();
 501     TObjArray *array = fVarManager->GetVariables();
 502     Int_t nvars = array->GetEntriesFast();
 503     TString namee;
 504     for(Int_t v = 0; v < nvars; v++) {
 505       AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
 506       if(!variable) continue;
 507       TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
 508       if(!name.CompareTo("source")) namee = TString("species");
 509       else namee = TString(name);
 510       Int_t nbins = variable->GetNumberOfBins();
 511       if(variable->HasUserDefinedBinning()){
 512         varManager->AddVariable(namee, nbins, variable->GetBinning());
 513       } else {
 514         varManager->AddVariable(namee, nbins, variable->GetMinimum(), variable->GetMaximum());
 515       }
 516       //printf("For AliTaggedTrackAnalysis, had the variable %s and the one used %s\n",(const char*)variable->GetName(),(const char*) namee);
 517     }
 518     if(fPIDqa->HasHighResolutionHistos())
 519       fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
 520     fTaggedTrackAnalysis->SetPID(fPID);
 521     fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
 522     fTaggedTrackAnalysis->InitContainer();
 523     fOutput->Add(fTaggedTrackAnalysis->GetContainer());
 524     fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
 525     fQA->Add(fTaggedTrackAnalysis->GetCutQA());
 526     fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
 527   }
 528
 529   PrintStatus();
 530   // Done!!!
 531   PostData(1, fOutput);
 532   PostData(2, fQA);
 533 }
 534
 535 //____________________________________________________________
 536 void AliAnalysisTaskHFE::UserExec(Option_t *){
 537   //
 538   // Run the analysis
 539   //
 540   AliDebug(3, "Starting Single Event Analysis");
 541   if(!fInputEvent){
 542     AliError("Reconstructed Event not available");
 543     return;
 544   }
 545   if(HasMCData()){
 546     AliDebug(4, Form("MC Event: %p", fMCEvent));
 547     if(!fMCEvent){
 548       AliError("No MC Event, but MC Data required");
 549       return;
 550     }
 551   }
 552   if(!fCuts){
 553     AliError("HFE cuts not available");
 554     return;
 555   }
 556   if(!fPID->IsInitialized()){
 557     // Initialize PID with the given run number
 558     fPID->InitializePID(fInputEvent->GetRunNumber());
 559   }
 560
 561   // Initialize hadronic background from OADB Container
 562   AliDebug(2, Form("Apply background factors: %s, OADB Container %p", fBackGroundFactorApply ? "Yes" : "No", fHadronBackgroundOADB));
 563   if(fBackGroundFactorApply && !TestBit(kBackgroundInitialized)){
 564     AliDebug(2, "Initializing Background from OADB");
 565     if(!InitializeHadronBackground(fInputEvent->GetRunNumber())) AliError("Failed initializing hadronic background parameterization from OADB");
 566     else AliDebug(2, "Successfully loaded Background from OADB");
 567     SetBit(kBackgroundInitialized);
 568   }
 569
 570   if(IsESDanalysis() && HasMCData()){
 571     // Protect against missing MC trees
 572     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 573     if(!mcH){
 574       AliError("No MC Event Handler available");
 575       return;
 576     }
 577     if(!mcH->InitOk()) return;
 578     if(!mcH->TreeK()) return;
 579     if(!mcH->TreeTR()) return;
 580   }
 581
 582   // need the centrality for everything (MC also)
 583   fCentralityF = -1;
 584   if(!ReadCentrality()) fCentralityF = -1;
 585   //printf("pass centrality\n");
 586   //printf("Reading fCentralityF %f\n",fCentralityF);
 587
 588   // See if pile up and z in the range
 589   RejectionPileUpVertexRangeEventCut();
 590
 591   // Protect agains missing
 592   if(HasMCData()){
 593     //printf("Has MC data\n");
 594     fSignalCuts->SetMCEvent(fMCEvent);
 595     ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
 596   }
 597
 598   AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse();
 599   if(!pidResponse){
 600     AliDebug(1, "Using default PID Response");
 601     pidResponse = AliHFEtools::GetDefaultPID(HasMCData(), fInputEvent->IsA() == AliAODEvent::Class());
 602   }
 603   fPID->SetPIDResponse(pidResponse);
 604   if(fPIDpreselect) fPIDpreselect->SetPIDResponse(pidResponse);
 605
 606   // Event loop
 607   if(IsAODanalysis()){
 608     //printf("test\n");
 609     ProcessAOD();
 610   } else {
 611     const char *specialTrigger = GetSpecialTrigger(fInputEvent->GetRunNumber());
 612     // Check Trigger selection
 613     if(specialTrigger){
 614       AliDebug(2, Form("Special Trigger requested: %s", specialTrigger));
 615       AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
 616       if(!(ev && ev->IsTriggerClassFired(specialTrigger))){
 617         AliDebug(2, "Event not selected");
 618         return;
 619       } else AliDebug(2, "Event Selected");
 620     } else AliDebug(2, "No Special Trigger requested");
 621
 622     ProcessESD();
 623   }
 624   // Done!!!
 625   PostData(1, fOutput);
 626   PostData(2, fQA);
 627 }
 628
 629 //____________________________________________________________
 630 void AliAnalysisTaskHFE::Terminate(Option_t *){
 631   //
 632   // Terminate not implemented at the moment
 633   //
 634   if(GetPlugin(kPostProcess)){
 635     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 636     fQA = dynamic_cast<TList *>(GetOutputData(2));
 637     if(!fOutput){
 638       AliError("Results not available");
 639       return;
 640     }
 641     if(!fQA){
 642       AliError("QA output not available");
 643       return;
 644     }
 645     fContainer = dynamic_cast<AliHFEcontainer *>(fOutput->FindObject("trackContainer"));
 646     if(!fContainer){
 647       AliError("Track container not found");
 648       return;
 649     }
 650     AliHFEpostAnalysis postanalysis;
 651     postanalysis.SetTaskResults(fContainer);
 652     TList *qalist = dynamic_cast<TList *>(fQA->FindObject("list_TaskQA"));
 653     if(!qalist){
 654       AliError("QA List not found");
 655       return;
 656     }
 657     postanalysis.SetTaskQA(qalist);
 658     printf("Running post analysis\n");
 659     //if(HasMCData())
 660     postanalysis.DrawMCSignal2Background();
 661     postanalysis.DrawEfficiency();
 662     postanalysis.DrawPIDperformance();
 663     postanalysis.DrawCutEfficiency();
 664
 665     if (GetPlugin(kIsElecBackGround)) {
 666       AliHFEelecbackground elecBackGround;
 667       TList *oe = 0x0;
 668       if(!(oe = (TList*)dynamic_cast<TList *>(fOutput->FindObject("HFEelecbackground")))){
 669         return;
 670       }
 671       elecBackGround.Load(oe);
 672       elecBackGround.Plot();
 673       elecBackGround.PostProcess();
 674     }
 675   }
 676 }
 677 //_______________________________________________________________
 678 Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
 679   //
 680   //
 681   //
 682
 683   //printf("test in IsEventInBinZero\n");
 684   if(!fInputEvent){
 685     AliError("Reconstructed Event not available");
 686     return kFALSE;
 687   }
 688
 689   // check vertex
 690   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
 691   if(!vertex) return kTRUE;
 692   //if(vertex) return kTRUE;
 693
 694   // check tracks
 695   if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
 696   //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
 697
 698
 699   return kFALSE;
 700
 701 }
 702 //____________________________________________________________
 703 void AliAnalysisTaskHFE::ProcessMC(){
 704   //
 705   // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
 706   // In case MC QA is on also MC QA loop is done
 707   //
 708   AliDebug(3, "Processing MC Information");
 709   Double_t eventContainer [4];
 710   eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
 711   eventContainer[2] = fCentralityF;
 712   eventContainer[3] = fContributors;
 713   fVz = eventContainer[0];
 714   //printf("z position is %f\n",eventContainer[0]);
 715   //if(fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent))
 716   fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
 717   Int_t nElectrons = 0;
 718   if(IsESDanalysis()){
 719    if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
 720     if (HasMCData() && IsQAOn(kMCqa)) {
 721       AliDebug(2, "Running MC QA");
 722
 723       if(fMCEvent->Stack()){
 724         fMCQA->SetMCEvent(fMCEvent);
 725         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
 726         fMCQA->SetCentrality(fCentralityF);
 727         fMCQA->SetPercentrality(fCentralityPercent);
 728         if(IsPbPb()) { fMCQA->SetPbPb();}
 729         else
 730         {
 731             if(fisppMultiBin) fMCQA->SetPPMultiBin();
 732             else fMCQA->SetPP();
 733         }
 734         fMCQA->Init();
 735
 736         fMCQA->GetMesonKine();
 737
 738         // loop over all tracks for decayed electrons
 739         for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
 740           TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
 741           if(!mcpart) continue;
 742           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
 743           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
 744           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
 745           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
 746           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG); // no accept cut
 747           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG); // no accept cut
 748           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG); // no accept cut
 749         }
 750         //fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 751         //fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 752       }
 753
 754     } // end of MC QA loop
 755    }
 756    // -----------------------------------------------------------------
 757    fCFM->SetMCEventInfo(fMCEvent);
 758    // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 759   } else {
 760     fCFM->SetMCEventInfo(fInputEvent);
 761   }
 762   // Run MC loop
 763   AliVParticle *mctrack = NULL;
 764   AliDebug(3, Form("Number of Tracks: %d", fMCEvent->GetNumberOfTracks()));
 765   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfTracks(); imc++){
 766     if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
 767     AliDebug(4, "Next MC Track");
 768     if(ProcessMCtrack(mctrack)) nElectrons++;
 769   }
 770
 771   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 772   fQACollection->Fill("nElectron", nElectrons);
 773 }
 774
 775 //____________________________________________________________
 776 void AliAnalysisTaskHFE::ProcessESD(){
 777   //
 778   // Run Analysis of reconstructed event in ESD Mode
 779   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
 780   //
 781   AliDebug(3, "Processing ESD Event");
 782   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 783   if(!fESD){
 784     AliError("ESD Event required for ESD Analysis");
 785     return;
 786   }
 787
 788   // Set magnetic field if V0 task on
 789   if(fTaggedTrackAnalysis) {
 790     fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
 791     fTaggedTrackAnalysis->SetCentrality(fCentralityF);
 792     if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
 793     else fTaggedTrackAnalysis->SetPP();
 794   }
 795
 796   // Do event Normalization
 797   Double_t eventContainer[4];
 798   eventContainer[0] = 0.;
 799   if(HasMCData()) eventContainer[0] = fVz;
 800   else {
 801     if(fESD->GetPrimaryVertexSPD()) eventContainer[0] = fESD->GetPrimaryVertexSPD()->GetZ();
 802   }
 803   eventContainer[1] = 0.;
 804   eventContainer[2] = fCentralityF;
 805   eventContainer[3] = fContributors;
 806   if(fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0AND))
 807     eventContainer[1] = 1.;
 808
 809   //
 810   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
 811
 812   //
 813   if(fIdentifiedAsPileUp) return;
 814   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
 815
 816   //
 817   if(TMath::Abs(fCentralityF) < 0) return;
 818   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
 819   //printf("In ProcessESD %f\n",fCentralityF);
 820
 821   //
 822   if(fIdentifiedAsOutInz) return;
 823   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
 824
 825   //
 826   if(!fPassTheEventCut) return;
 827   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
 828
 829
 830
 831   fContainer->NewEvent();
 832
 833   if (GetPlugin(kIsElecBackGround)) {
 834     fElecBackGround->SetEvent(fESD);
 835   }
 836   if (GetPlugin(kSecVtx)) {
 837     fSecVtx->SetEvent(fESD);
 838     fSecVtx->GetPrimaryCondition();
 839   }
 840
 841   if(HasMCData()){
 842     if (GetPlugin(kSecVtx)) {
 843       fSecVtx->SetMCEvent(fMCEvent);
 844       fSecVtx->SetMCQA(fMCQA);
 845     }
 846     if (GetPlugin(kIsElecBackGround)) {
 847       fElecBackGround->SetMCEvent(fMCEvent);
 848     }
 849   }
 850
 851   Double_t container[10];
 852   memset(container, 0, sizeof(Double_t) * 10);
 853   // container for the output THnSparse
 854   Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B']
 855   Double_t dataDca[8]; // [source, pT, dca, dcaSig, centrality]
 856   Double_t eDca[6]; // [source, pT, dca, dcaSig, centrality]
 857   Int_t nElectronCandidates = 0;
 858   AliESDtrack *track = NULL, *htrack = NULL;
 859   AliMCParticle *mctrack = NULL;
 860   AliMCParticle *mctrackmother = NULL;
 861   Int_t pid = 0;
 862
 863   Bool_t signal = kTRUE;
 864
 865   fCFM->SetRecEventInfo(fESD);
 866
 867   // minjung for IP QA(temporary ~ 2weeks)
 868   if(!fExtraCuts){
 869     fExtraCuts = new AliHFEextraCuts("hfetmpCuts","HFE tmp Cuts");
 870   }
 871   fExtraCuts->SetRecEventInfo(fESD);
 872
 873   // Electron background analysis
 874   if (GetPlugin(kIsElecBackGround)) {
 875
 876     AliDebug(2, "Running BackGround Analysis");
 877
 878     fElecBackGround->Reset();
 879
 880   } // end of electron background analysis
 881   //
 882   // Loop ESD
 883   //
 884   AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
 885   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 886     AliDebug(4, "New ESD track");
 887     track = fESD->GetTrack(itrack);
 888     track->SetESDEvent(fESD);
 889
 890     // fill counts of v0-identified particles
 891     Int_t v0pid = -1;
 892     if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
 893     else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
 894     else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
 895     // here the tagged track analysis will run
 896     if(fTaggedTrackAnalysis && v0pid > -1){
 897       AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
 898       fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
 899     }
 900
 901     AliDebug(3, Form("Doing track %d, %p", itrack, track));
 902
 903     //////////////////////////////////////
 904     // preselect
 905     /////////////////////////////////////
 906     if(fPIDpreselect && fCutspreselect) {
 907       if(!PreSelectTrack(track)) continue;
 908     }
 909
 910     signal = kTRUE;
 911
 912     // Fill step without any cut
 913
 914     if(HasMCData()){
 915       // Check if it is electrons near the vertex
 916       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
 917
 918       if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 919       else AliDebug(3, "Signal Electron");
 920
 921       // Fill K pt for Ke3 contributions
 922       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==321)) fQACollection->Fill("Kptspectra",mctrack->Pt());
 923       else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode())==130)) fQACollection->Fill("K0Lptspectra",mctrack->Pt());
 924     }
 925     // Cache new Track information inside the var manager
 926     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
 927
 928     if(fFillNoCuts) {
 929       if(signal || !fFillSignalOnly){
 930         fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
 931         fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
 932       }
 933     }
 934
 935     // RecKine: ITSTPC cuts
 936     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 937
 938
 939     // RecPrim
 940     if(fRejectKinkMother) {
 941       if(track->GetKinkIndex(0) != 0) continue; } // Quick and dirty fix to reject both kink mothers and daughters
 942     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 943
 944     // HFEcuts: ITS layers cuts
 945     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 946
 947     // HFE cuts: TOF PID and mismatch flag
 948     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
 949
 950     // HFE cuts: TPC PID cleanup
 951     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
 952
 953     // HFEcuts: Nb of tracklets TRD0
 954     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
 955
 956     // Fill correlation maps before PID
 957     if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
 958       //printf("Fill correlation maps before PID\n");
 959       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
 960     }
 961
 962     if(HasMCData()){
 963       FillProductionVertex(track);
 964
 965       if(fMCQA && signal){
 966         fMCQA->SetCentrality(fCentralityF);
 967         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){
 968          Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.};
 969          Double_t hfeimpactRtmp=0., hfeimpactnsigmaRtmp=0.;
 970          fExtraCuts->GetHFEImpactParameters(track, hfeimpactRtmp, hfeimpactnsigmaRtmp);
 971          UChar_t itsPixel = track->GetITSClusterMap();
 972          Double_t ilyrhit=0, ilyrstat=0;
 973          for(Int_t ilyr=0; ilyr<6; ilyr++){
 974            if(TESTBIT(itsPixel, ilyr)) ilyrhit += TMath::Power(2,ilyr);
 975            if(fExtraCuts->CheckITSstatus(fExtraCuts->GetITSstatus(track,ilyr))) ilyrstat += TMath::Power(2,ilyr);
 976          }
 977          fMCQA->SetITSInfo(ilyrhit,ilyrstat);
 978          fMCQA->SetHFEImpactParameters(hfeimpactRtmp, hfeimpactnsigmaRtmp);
 979          fMCQA->SetTrkKine(track->Pt(),track->Eta(), track->Phi());
 980          fMCQA->SetContainerStep(3);
 981          for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
 982            weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
 983            if(!fisNonHFEsystematics)break;
 984          }
 985
 986          if(fisNonHFEsystematics){
 987            //Fill additional containers for electron source distinction
 988            Int_t elecSource = 0;
 989            elecSource = fMCQA->GetElecSource(mctrack->Particle());
 990            const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
 991            const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
 992            Int_t iName = 0;
 993            for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
 994              if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
 995              if(elecSource == iSource){
 996                for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
 997                  if(weightElecBgV0[iLevel]>0){
 998                    fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, weightElecBgV0[iLevel]);
 999                  }
1000                  else if(weightElecBgV0[iLevel]<0){
1001                    fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 3, kFALSE, -1*weightElecBgV0[iLevel]);
1002                  }
1003                }
1004                break;
1005              }
1006              iName++;
1007              if(iName == kElecBgSpecies)iName = 0;
1008            }
1009          }
1010          //else{
1011            if(weightElecBgV0[0]>0) {
1012              fVarManager->FillContainer(fContainer, "conversionElecs", 3, kFALSE, weightElecBgV0[0]);
1013              fVarManager->FillContainer(fContainer, "conversionElecs", 4, kTRUE, weightElecBgV0[0]);
1014            }
1015            else if(weightElecBgV0[0]<0) {
1016              fVarManager->FillContainer(fContainer, "mesonElecs", 3, kFALSE, -1*weightElecBgV0[0]);
1017              fVarManager->FillContainer(fContainer, "mesonElecs", 4, kTRUE, -1*weightElecBgV0[0]);
1018            }
1019            //}
1020         }
1021       }
1022
1023       // check if it is the proton from the lambda decay, and yes fill the dca info
1024       Double_t hfeimpactR4p=0., hfeimpactnsigmaR4p=0.;
1025       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 211)){
1026         if(track->Pt()>4.){
1027           fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4p, hfeimpactnsigmaR4p);
1028           fQACollection->Fill("pionDcaSig",track->Pt(),hfeimpactnsigmaR4p);
1029         }
1030       }
1031       else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 2212)){
1032         fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4p, hfeimpactnsigmaR4p);
1033         fQACollection->Fill("protonDcaSig",track->Pt(),hfeimpactnsigmaR4p);
1034         Int_t glabel=TMath::Abs(mctrack->GetMother());
1035         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1036           if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==3122){
1037             fQACollection->Fill("secondaryprotonDcaSig",track->Pt(),hfeimpactnsigmaR4p);
1038           }
1039         }
1040       }
1041       else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) == 11)){ // to increas statistics for Martin
1042         if(signal){
1043           fExtraCuts->GetHFEImpactParameters(track, hfeimpactR4p, hfeimpactnsigmaR4p);
1044           Int_t sourceDca =-1;
1045           if(fSignalCuts->IsCharmElectron(track)){
1046             sourceDca=1;
1047           }
1048           else if(fSignalCuts->IsBeautyElectron(track)){
1049             sourceDca=2;
1050           }
1051           else if(fSignalCuts->IsGammaElectron(track)){
1052             sourceDca=3;
1053           }
1054           else if(fSignalCuts->IsNonHFElectron(track)){
1055             sourceDca=4;
1056           }
1057           else if(fSignalCuts->IsJpsiElectron(track)){
1058             sourceDca=5;
1059           }
1060           else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1061             sourceDca=6;
1062             fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
1063             fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
1064           }
1065           else {
1066             sourceDca=7;
1067           }
1068           eDca[0]=sourceDca;
1069           eDca[1]=track->Pt();
1070           eDca[2]=hfeimpactR4p;
1071           eDca[3]=hfeimpactnsigmaR4p;
1072           eDca[4]=fCentralityF;
1073           eDca[5] = track->Charge()/3;
1074
1075           fQACollection->Fill("eDca", eDca);
1076         }
1077       }
1078     }
1079
1080     if(TMath::Abs(track->Eta()) < 0.5){
1081       if(track->GetInnerParam())
1082         fQACollection->Fill("TPCdEdxBeforePID", track->GetInnerParam()->P(), track->GetTPCsignal());
1083       fQACollection->Fill("TPCnSigmaBeforePID", track->P(), fInputHandler->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron));
1084     }
1085
1086     AliHFEpidObject hfetrack;
1087     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1088     hfetrack.SetRecTrack(track);
1089     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1090     hfetrack.SetCentrality(fCentralityF);
1091     if(IsPbPb()) hfetrack.SetPbPb();
1092     else hfetrack.SetPP();
1093     fPID->SetVarManager(fVarManager);
1094     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
1095     nElectronCandidates++;
1096
1097     // Temporary histogram for chi2/ITS cluster
1098     if(IsPbPb()) {
1099             TBits shared = track->GetTPCSharedMap();
1100             Int_t sharebit=0;
1101             if(shared.CountBits() >= 2) sharebit=1;
1102
1103             Double_t itschi2percluster = 0.0;
1104             Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
1105             if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
1106
1107             Double_t itsChi2[7] = {track->Pt(),track->Eta(), track->Phi(),
1108                                    fCentralityF,track->GetTPCsignalN(), sharebit,itschi2percluster};
1109             fQACollection->Fill("fChi2perITScluster", itsChi2);
1110     }
1111     else{
1112
1113       Double_t itschi2percluster = 0.0;
1114       Double_t itsnbcls = static_cast<Double_t>(track->GetNcls(0));
1115       if(itsnbcls > 0) itschi2percluster = track->GetITSchi2()/itsnbcls;
1116
1117       Double_t itsChi2[3] = {track->Pt(), fCentralityF, itschi2percluster};
1118       fQACollection->Fill("fChi2perITScluster", itsChi2);
1119     }
1120
1121     // Fill Histogram for Hadronic Background
1122     if(HasMCData()){
1123       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
1124         fVarManager->FillContainer(fContainer, "hadronicBackground", UInt_t(0), kFALSE);
1125       else if(mctrack){
1126         // Fill Ke3 contributions
1127         Int_t glabel=TMath::Abs(mctrack->GetMother());
1128         if((mctrackmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(glabel)))){
1129           if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==321)
1130             fQACollection->Fill("Ke3Kecorr",mctrack->Pt(),mctrackmother->Pt());
1131           else if(TMath::Abs(mctrackmother->Particle()->GetPdgCode())==130)
1132             fQACollection->Fill("Ke3K0Lecorr",mctrack->Pt(),mctrackmother->Pt());
1133         }
1134       }
1135     }
1136
1137     // Fill Containers
1138     if(signal) {
1139       // Apply weight for background contamination
1140             if(fBackGroundFactorApply) {
1141               if(IsPbPb() && fCentralityF >= 0) fWeightBackGround =  fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
1142               else    fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
1143
1144               if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1145               else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1146         // weightBackGround as special weight
1147         fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
1148       }
1149       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
1150     }
1151
1152     Bool_t bTagged=kFALSE;
1153     if(GetPlugin(kSecVtx)) {
1154       AliDebug(2, "Running Secondary Vertex Analysis");
1155       if(fSecVtx->Process(track) && signal) {
1156         fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
1157         fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
1158         bTagged=kTRUE;
1159       }
1160     }
1161
1162     if(HasMCData()){
1163       dataE[0] = track->Pt();
1164       dataE[1] = track->Eta();
1165       dataE[2] = track->Phi();
1166       dataE[3] = track->Charge();
1167       dataE[4] = -1.;
1168       dataE[5] = -1.;
1169
1170       // Track selected: distinguish between true and fake
1171       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
1172       if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
1173         Int_t type = 0;
1174         if(fSignalCuts->IsCharmElectron(track))
1175           type = 1;
1176         else if(fSignalCuts->IsBeautyElectron(track))
1177           type = 2;
1178         AliDebug(1, Form("Type: %d\n", type));
1179         if(type){
1180                 dataE[5] = type; // beauty[1] or charm[2]
1181                 dataE[4] = 2;  // signal electron
1182         }
1183         else{
1184                 dataE[4] = 1; // not a signal electron
1185                 dataE[5] = 0;
1186         }
1187       }
1188       else {
1189         // Fill THnSparse with the information for Fake Electrons
1190         dataE[4] = 0;
1191         dataE[5] = 0;
1192       }
1193       // fill the performance THnSparse, if the mc origin could be defined
1194       if(dataE[4] > -1){
1195         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1196         fQACollection->Fill("PIDperformance", dataE);
1197       }
1198     }
1199
1200     // Electron background analysis
1201     if (GetPlugin(kIsElecBackGround)) {
1202
1203       AliDebug(2, "Running BackGround Analysis");
1204
1205       for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
1206         htrack = fESD->GetTrack(jtrack);
1207         if ( itrack == jtrack ) continue;
1208         fElecBackGround->PairAnalysis(track, htrack);
1209       }
1210     } // end of electron background analysis
1211
1212
1213     Int_t sourceDca =-1;
1214     if (GetPlugin(kDEstep)) {
1215       Double_t weightElecBgV0[kBgLevels] = {0.,0.,0.,};
1216       Int_t elecSource = 0;
1217       // minjung for IP QA(temporary ~ 2weeks)
1218       Double_t hfeimpactR=0., hfeimpactnsigmaR=0.;
1219       fExtraCuts->GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
1220       sourceDca=0;
1221       if(HasMCData())
1222       {
1223         if(fSignalCuts->IsCharmElectron(track)){
1224           sourceDca=1;
1225         }
1226         else if(fSignalCuts->IsBeautyElectron(track)){
1227           sourceDca=2;
1228         }
1229         else if(fSignalCuts->IsGammaElectron(track)){
1230           sourceDca=3;
1231         }
1232         else if(fSignalCuts->IsNonHFElectron(track)){
1233           sourceDca=4;
1234         }
1235         else if(fSignalCuts->IsJpsiElectron(track)){
1236           sourceDca=5;
1237         }
1238         else if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1239           sourceDca=6;
1240           fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
1241           fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
1242         }
1243         else {
1244           sourceDca=7;
1245         }
1246
1247         if(fMCQA && signal) {
1248
1249           fMCQA->SetContainerStep(0);
1250           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1251             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1252             if(!fisNonHFEsystematics)break;
1253           }
1254
1255           if(fisNonHFEsystematics){
1256             //Fill additional containers for electron source distinction
1257             elecSource = fMCQA->GetElecSource(mctrack->Particle());
1258             const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1259             const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1260             Int_t iName = 0;
1261             for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
1262               if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1263               if(elecSource == iSource){
1264                 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1265                   if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, weightElecBgV0[iLevel]);
1266                   else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 0, kFALSE, -1*weightElecBgV0[iLevel]);
1267                 }
1268                 break;
1269               }
1270               iName++;
1271               if(iName == kElecBgSpecies)iName = 0;
1272             }
1273           }
1274           //else{
1275           if(weightElecBgV0[0]>0) {
1276             fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBgV0[0]);
1277             fVarManager->FillContainer(fContainer, "conversionElecs", 5, kTRUE, weightElecBgV0[0]);
1278           }
1279           else if(weightElecBgV0[0]<0) {
1280             fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBgV0[0]);
1281             fVarManager->FillContainer(fContainer, "mesonElecs", 5, kTRUE, -1*weightElecBgV0[0]);
1282           }
1283           //}
1284           if(bTagged){ // bg estimation for the secondary vertex tagged signals
1285             if(weightElecBgV0[0]>0) fVarManager->FillContainer(fContainer, "conversionElecs", 2, kFALSE, weightElecBgV0[0]);
1286             else if(weightElecBgV0[0]<0) fVarManager->FillContainer(fContainer, "mesonElecs", 2, kFALSE, -1*weightElecBgV0[0]);
1287           }
1288         }
1289       } // end of MC
1290
1291       dataDca[0]=sourceDca;
1292       dataDca[1]=track->Pt();
1293       dataDca[2]=hfeimpactR;
1294       dataDca[3]=hfeimpactnsigmaR;
1295       dataDca[4]=fCentralityF;
1296       dataDca[5] = 49;
1297       Double_t xr[3]={49,49,49};
1298       if(HasMCData()) {
1299         mctrack->XvYvZv(xr);
1300         dataDca[5] = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
1301       }
1302       dataDca[6] = v0pid;
1303       dataDca[7] = track->Charge()/3;
1304
1305  //     printf("Entries dca: [%.3f|%.3f|%.3f|%f|%f]\n", dataDca[0],dataDca[1],dataDca[2],dataDca[3],dataDca[4]);
1306       if (!HasMCData()) fQACollection->Fill("Dca", dataDca);
1307       else if(signal) fQACollection->Fill("Dca", dataDca);
1308
1309
1310       // Fill Containers for impact parameter analysis
1311       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
1312       if(HasMCData()){
1313         if(fMCQA && signal) {
1314           fMCQA->SetContainerStep(1);
1315           for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1316             weightElecBgV0[iLevel] = fMCQA->GetWeightFactor(mctrack, iLevel); // positive:conversion e, negative: nonHFE
1317             if(!fisNonHFEsystematics)break;
1318           }
1319           if(fisNonHFEsystematics){
1320             //Fill additional containers for electron source distinction
1321             elecSource = fMCQA->GetElecSource(mctrack->Particle());
1322             const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1323             const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1324             Int_t iName = 0;
1325             for(Int_t iSource = AliHFEmcQA::kPi0; iSource <=  AliHFEmcQA::kGammaRho0; iSource++){
1326               if((iSource == AliHFEmcQA::kElse)||(iSource == AliHFEmcQA::kMisID)) continue;
1327               if(elecSource == iSource){
1328                 for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1329                   if(weightElecBgV0[iLevel]>0) fVarManager->FillContainer(fContainer, Form("conversionElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, weightElecBgV0[iLevel]);
1330                   else if(weightElecBgV0[iLevel]<0) fVarManager->FillContainer(fContainer, Form("mesonElecs%s%s",sourceName[iName], levelName[iLevel]), 1, kFALSE, -1*weightElecBgV0[iLevel]);
1331                 }
1332                 break;
1333               }
1334               iName++;
1335               if(iName == kElecBgSpecies)iName = 0;
1336             }
1337           }
1338           // else{
1339             if(weightElecBgV0[0]>0) {
1340               fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBgV0[0]);
1341               fVarManager->FillContainer(fContainer, "conversionElecs", 6, kTRUE, weightElecBgV0[0]);
1342             }
1343             else if(weightElecBgV0[0]<0) {
1344               fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBgV0[0]);
1345               fVarManager->FillContainer(fContainer, "mesonElecs", 6, kTRUE, -1*weightElecBgV0[0]);
1346             }
1347             //}
1348         }
1349       }
1350       if(signal) {
1351         fVarManager->FillContainer(fContainer, "recTrackContDEReco", AliHFEcuts::kStepHFEcutsDca, kFALSE);
1352         fVarManager->FillContainer(fContainer, "recTrackContDEMC", AliHFEcuts::kStepHFEcutsDca, kTRUE);
1353         fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterDE"));
1354       }
1355       if(HasMCData()){
1356         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1357           fQACollection->Fill("hadronsAfterIPcut",track->Pt());
1358           fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
1359         }
1360       }
1361     }
1362
1363   }
1364   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1365 }
1366
1367 //____________________________________________________________
1368 void AliAnalysisTaskHFE::ProcessAOD(){
1369   //
1370   // Run Analysis in AOD Mode
1371   // Function is still in development
1372   //
1373   //printf("Process AOD\n");
1374   AliDebug(3, "Processing AOD Event");
1375   Double_t eventContainer[4];
1376   if(HasMCData()) eventContainer[0] = fVz;
1377   else {
1378     eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
1379   }
1380   eventContainer[1] = 1.; // No Information available in AOD analysis, assume all events have V0AND
1381   eventContainer[2] = fCentralityF;
1382   eventContainer[3] = fContributors;
1383
1384   //printf("value event container %f, %f, %f, %f\n",eventContainer[0],eventContainer[1],eventContainer[2],eventContainer[3]);
1385
1386   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1387   if(!fAOD){
1388     AliError("AOD Event required for AOD Analysis");
1389       return;
1390   }
1391
1392   //printf("Will fill\n");
1393   //
1394   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
1395   //printf("Fill\n");
1396   //
1397   if(fIdentifiedAsPileUp) return;
1398   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
1399
1400   //
1401   if(fIdentifiedAsOutInz) return;
1402   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
1403
1404   //
1405   if(!fPassTheEventCut) return;
1406   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
1407   //printf("pass\n");
1408
1409   fContainer->NewEvent();
1410
1411   AliAODTrack *track = NULL;
1412   AliAODMCParticle *mctrack = NULL;
1413   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
1414   Int_t nElectronCandidates = 0;
1415   Int_t pid;
1416   Bool_t signal;
1417   //printf("Number of track %d\n",(Int_t) fAOD->GetNumberOfTracks());
1418   for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
1419     track = fAOD->GetTrack(itrack); mctrack = NULL;
1420     if(!track) continue;
1421     if(fUseFlagAOD){
1422       if(track->GetFlags() != fFlags) continue;  // Only process AOD tracks where the HFE is set
1423     }
1424     //printf("Pass the flag\n");
1425
1426     signal = kTRUE;
1427     if(HasMCData()){
1428
1429       Int_t label = TMath::Abs(track->GetLabel());
1430       if(label)
1431         mctrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(label));
1432         if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
1433     }
1434     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, kTRUE);
1435
1436     if(fFillNoCuts) {
1437       if(signal || !fFillSignalOnly){
1438         fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
1439         fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
1440       }
1441     }
1442
1443     if(fApplyCutAOD) {
1444       // RecKine: ITSTPC cuts
1445       if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
1446
1447       // RecPrim
1448       if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
1449
1450       // HFEcuts: ITS layers cuts
1451       if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
1452
1453       // HFE cuts: TOF PID and mismatch flag
1454       if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
1455
1456       // HFE cuts: TPC PID cleanup
1457       if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
1458
1459       // HFEcuts: Nb of tracklets TRD0
1460       if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
1461     }
1462
1463     // Fill correlation maps before PID
1464     if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
1465       //printf("Fill correlation maps before PID\n");
1466       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
1467     }
1468
1469     //printf("Will process to PID\n");
1470
1471     // track accepted, do PID
1472     AliHFEpidObject hfetrack;
1473     hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1474     hfetrack.SetRecTrack(track);
1475     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1476     hfetrack.SetCentrality(fCentralityF);
1477     if(IsPbPb()) hfetrack.SetPbPb();
1478     else hfetrack.SetPP();
1479     fPID->SetVarManager(fVarManager);
1480     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;    // we will do PID here as soon as possible
1481
1482
1483     // Apply weight for background contamination
1484     //Double_t weightBackGround = 1.0;
1485     if(signal) {
1486       // Apply weight for background contamination
1487       if(fBackGroundFactorApply) {
1488         if(IsPbPb() && fCentralityF >= 0) fWeightBackGround =  fkBackGroundFactorArray[fCentralityF >= 0 ? fCentralityF : 0]->Eval(TMath::Abs(track->P()));
1489         else    fWeightBackGround =  fkBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
1490
1491         if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
1492         else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
1493         // weightBackGround as special weight
1494         fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
1495       }
1496       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
1497     }
1498
1499     nElectronCandidates++;
1500     if(HasMCData()){
1501       dataE[0] = track->Pt();
1502       dataE[1] = track->Eta();
1503       dataE[2] = track->Phi();
1504       dataE[3] = track->Charge();
1505       dataE[4] = -1;
1506       dataE[5] = -1;
1507       // Track selected: distinguish between true and fake
1508       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack ? mctrack->GetPdgCode(): -1));
1509       if(mctrack && ((pid = TMath::Abs(mctrack->GetPdgCode())) == 11)){
1510
1511         Int_t type = 0;
1512         if(fSignalCuts->IsCharmElectron(track))
1513           type = 1;
1514         else if(fSignalCuts->IsBeautyElectron(track))
1515           type = 2;
1516         AliDebug(1, Form("Type: %d\n", type));
1517         if(type){
1518                 dataE[5] = type; // beauty[1] or charm[2]
1519                 dataE[4] = 2;  // signal electron
1520         }
1521         else{
1522                 dataE[4] = 1; // not a signal electron
1523                 dataE[5] = 0;
1524         }
1525       }
1526       else {
1527         // Fill THnSparse with the information for Fake Electrons
1528         dataE[4] = 0;
1529         dataE[5] = 0;
1530       }
1531       // fill the performance THnSparse, if the mc origin could be defined
1532       if(dataE[4] > -1){
1533         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1534         fQACollection->Fill("PIDperformance", dataE);
1535       }
1536     }
1537   }
1538   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1539 }
1540
1541 //____________________________________________________________
1542 Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
1543   //
1544   // Filter the Monte Carlo Track
1545   // Additionally Fill a THnSparse for Signal To Background Studies
1546   // Works for AOD and MC analysis Type
1547   //
1548   fVarManager->NewTrack(track, NULL, fCentralityF, -1, kTRUE);
1549   Double_t signalContainer[6];
1550
1551   signalContainer[0] = track->Pt();
1552   signalContainer[1] = track->Eta();
1553   signalContainer[2] = track->Phi();
1554   signalContainer[3] = track->Charge()/3;
1555
1556   Double_t vertex[3] = {0.,0.,0.}; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
1557   if(IsESDanalysis()){
1558     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
1559     if(mctrack){
1560       vertex[0] = mctrack->Particle()->Vx();
1561       vertex[1] = mctrack->Particle()->Vy();
1562     }
1563   } else {
1564     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
1565     if(aodmctrack) aodmctrack->XvYvZv(vertex);
1566   }
1567
1568   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
1569   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
1570   signalContainer[4] = 0;
1571   if(fSignalCuts->IsSelected(track)){
1572     //fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCsignal, kFALSE);
1573     // Filling of the Signal/Background histogram using the
1574     // definition of the codes for charm and beauty as below in
1575     // th crearion of the histogram
1576     if(fSignalCuts->IsCharmElectron(track))
1577       signalContainer[4] = 1;
1578     else
1579       signalContainer[4] = 2;
1580   } else {
1581     signalContainer[4] = 0; // (and other background)
1582   }
1583   signalContainer[5] = 0;
1584   // apply cut on the sqrt of the production vertex
1585   Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
1586   if(radVertex < 3.5){
1587     // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
1588     signalContainer[5] = 1;
1589   } else if (radVertex < 7.5){
1590     signalContainer[5] = 2;
1591   }
1592   fQACollection->Fill("SignalToBackgroundMC", signalContainer);
1593
1594   // Step GeneratedZOutNoPileUp
1595   if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0)) return kFALSE;
1596   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
1597   //printf("In ProcessMCtrack %f\n",fCentralityF);
1598
1599   // Step Generated Event Cut
1600   if(!fPassTheEventCut) return kFALSE;
1601   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedEventCut, kFALSE);
1602
1603   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
1604   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCInAcceptance, kFALSE);
1605   return kTRUE;
1606 }
1607
1608 //____________________________________________________________
1609 Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
1610   //
1611   // Preselect tracks
1612   //
1613
1614
1615   Bool_t survived = kTRUE;
1616
1617   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) {
1618     survived = kFALSE;
1619     //printf("Did not pass AliHFEcuts::kStepRecKineITSTPC\n");
1620   }
1621   //else printf("Pass AliHFEcuts::kStepRecKineITSTPC\n");
1622   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) {
1623     survived = kFALSE;
1624     //printf("Did not pass AliHFEcuts::kStepRecPrim\n");
1625   }
1626   //else printf("Pass AliHFEcuts::kStepRecPrim\n");
1627   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) {
1628     survived = kFALSE;
1629     //printf("Did not pass AliHFEcuts::kStepHFEcutsITS\n");
1630   }
1631   //else printf("Pass AliHFEcuts::kStepHFEcutsITS\n");
1632   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTOF, track)) {
1633     survived = kFALSE;
1634     //printf("Did not pass AliHFEcuts::kStepHFEcutsTOF\n");
1635   }
1636   //else printf("Pass AliHFEcuts::kStepHFEcutsTOF\n");
1637   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) {
1638     survived = kFALSE;
1639     //printf("Did not pass AliHFEcuts::kStepHFEcutsTRD\n");
1640   }
1641   //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
1642
1643   if(survived){
1644     // Apply PID
1645     AliHFEpidObject hfetrack;
1646     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1647     hfetrack.SetRecTrack(track);
1648     if(!fPIDpreselect->IsSelected(&hfetrack)) {
1649       //printf("Did not pass AliHFEcuts::kPID\n");
1650       survived = kFALSE;
1651     }
1652     //else printf("Pass AliHFEcuts::kPID\n");
1653   }
1654
1655   return survived;
1656
1657 }
1658 //____________________________________________________________
1659 void AliAnalysisTaskHFE::MakeEventContainer(){
1660   //
1661   // Create the event container for the correction framework and link it
1662   // 1st bin: Vertex z-position
1663   // 2nd bin: V0AND decision (normalization to sigma_inel)
1664   // 3rd bin: Centrality class (for pp defined as number of contributors in vertex.)
1665   // 4th bin: Number of contributors > 0
1666   //
1667
1668   const Int_t kNvar = 4;  // number of variables on the grid:
1669   Int_t nBins[kNvar] = {120, 2, 11, 2};
1670   Double_t binMin[kNvar] = {-30. , 0., 0.0, 0.};
1671   Double_t binMax[kNvar] = {30., 2., 11.0, 2.};
1672
1673   AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1674
1675   Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1676   Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1677   Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1678   Double_t *contributorsBins = AliHFEtools::MakeLinearBinning(nBins[3], binMin[3], binMax[3]);
1679   evCont->SetBinLimits(0, vertexBins);
1680   evCont->SetBinLimits(1, v0andBins);
1681   evCont->SetBinLimits(2, centralityBins);
1682   evCont->SetBinLimits(3, contributorsBins);
1683   delete[] vertexBins; delete[] v0andBins; delete[] centralityBins; delete[] contributorsBins;
1684
1685   fCFM->SetEventContainer(evCont);
1686 }
1687
1688 //____________________________________________________________
1689 void AliAnalysisTaskHFE::MakeParticleContainer(){
1690   //
1691   // Create the particle container for the correction framework manager and
1692   // link it
1693   //
1694   if(!fContainer) fContainer = new AliHFEcontainer("trackContainer");
1695   fVarManager->DefineVariables(fContainer);
1696
1697   // Create Correction Framework containers
1698   fContainer->CreateContainer("MCTrackCont", "Track Container filled with MC information", AliHFEcuts::kNcutStepsMCTrack);
1699   fContainer->CreateContainer("recTrackContReco", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1700   fContainer->CreateContainer("recTrackContMC", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1701
1702   fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 2);
1703   fContainer->CreateContainer("recTrackContDEReco", "Container for displaced electron analysis with Reco information", 1);
1704   fContainer->CreateContainer("recTrackContDEMC", "Container for displaced electron analysis with MC information", 1);
1705   fContainer->CreateContainer("recTrackContSecvtxReco", "Container for secondary vertexing analysis with Reco information", 1);
1706   fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
1707
1708   if(HasMCData()){
1709     fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",7);
1710     fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from meson decays",7);
1711     fContainer->Sumw2("conversionElecs");
1712     fContainer->Sumw2("mesonElecs");
1713
1714     if(fisNonHFEsystematics){
1715       const Char_t *sourceName[kElecBgSpecies]={"Pion","Eta","Omega","Phi","EtaPrime","Rho"};
1716       const Char_t *levelName[kBgLevels]={"Best","Lower","Upper"};
1717       for(Int_t iSource = 0; iSource < kElecBgSpecies; iSource++){
1718         for(Int_t iLevel = 0; iLevel < kBgLevels; iLevel++){
1719           fContainer->CreateContainer(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted conversion electrons from %s grandm., %s level",sourceName[iSource],levelName[iLevel]),4);
1720           fContainer->CreateContainer(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]), Form("Container for weighted electrons from %s decays, %s level",sourceName[iSource],levelName[iLevel]),4);
1721           fContainer->Sumw2(Form("conversionElecs%s%s",sourceName[iSource],levelName[iLevel]));
1722           fContainer->Sumw2(Form("mesonElecs%s%s",sourceName[iSource],levelName[iLevel]));
1723         }
1724       }
1725     }
1726     //fContainer->CreateContainer("charmElecs", "Container for weighted charm electrons",2);
1727   }
1728
1729   fContainer->CreateCorrelationMatrix("correlationstepafterPID","THnSparse with correlations");
1730   fContainer->CreateCorrelationMatrix("correlationstepafterDE","THnSparse with correlations");
1731   if(!fVarManager->IsVariableDefined("centrality")) {
1732     //printf("Create the two other correlation maps\n");
1733     fContainer->CreateCorrelationMatrix("correlationstepbeforePID","THnSparse with correlations");
1734     fContainer->CreateCorrelationMatrix("correlationstepafterTOF","THnSparse with correlations");
1735   }
1736
1737   // Define the step names
1738   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsMCTrack; istep++){
1739     fContainer->SetStepTitle("MCTrackCont", AliHFEcuts::MCCutName(istep), istep);
1740   }
1741   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsRecTrack; istep++){
1742     fContainer->SetStepTitle("recTrackContReco", AliHFEcuts::RecoCutName(istep), istep);
1743     fContainer->SetStepTitle("recTrackContMC", AliHFEcuts::RecoCutName(istep), istep);
1744   }
1745   for(UInt_t ipid = 0; ipid < fPID->GetNumberOfPIDdetectors(); ipid++){
1746     fContainer->SetStepTitle("recTrackContReco", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1747     fContainer->SetStepTitle("recTrackContMC", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1748   }
1749 }
1750
1751 //____________________________________________________________
1752 void AliAnalysisTaskHFE::InitPIDperformanceQA(){
1753   // Add a histogram for Fake electrons
1754   const Int_t nDim=6;
1755   Int_t nBin[nDim] = {40, 8, 18, 2, 3, 3};
1756   //number of variables on the grid:pt,eta,phi,charge,
1757   const Double_t kPtbound[2] = {0.1, 20.};
1758   const Double_t kEtabound[2] = {-0.8, 0.8};
1759   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
1760   const Double_t kChargebound[2] = {-1.1, 1.1};
1761   const Double_t kAddInf1bound[2] = {0., 3.};
1762   const Double_t kAddInf2bound[2] = {0., 3.};
1763   Double_t minima[nDim] = {kPtbound[0], kEtabound[0], kPhibound[0], kChargebound[0], kAddInf1bound[0], kAddInf2bound[0]};
1764   Double_t maxima[nDim] = {kPtbound[1], kEtabound[1], kPhibound[1], kChargebound[1], kAddInf1bound[1], kAddInf2bound[1]};
1765
1766   fQACollection->CreateTHnSparse("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin, minima, maxima);
1767   fQACollection->CreateTHnSparse("SignalToBackgroundMC", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; flavor (0 - no, 1 - charm, 2 - bottom); ITS Cluster (0 - no, 1 - first (and maybe second), 2 - second)", nDim, nBin, minima, maxima);
1768
1769   fQACollection->BinLogAxis("PIDperformance", 0);
1770   fQACollection->BinLogAxis("SignalToBackgroundMC", 0);
1771   fQACollection->Sumw2("PIDperformance");
1772   fQACollection->Sumw2("SignalToBackgroundMC");
1773 }
1774
1775 //____________________________________________________________
1776 void AliAnalysisTaskHFE::InitContaminationQA(){
1777   //
1778   // Add QA for Impact Parameter cut
1779   //
1780
1781   TObjArray *array = fVarManager->GetVariables();
1782   Int_t nvars = array->GetEntriesFast();
1783   for(Int_t v = 0; v < nvars; v++) {
1784     AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
1785     if(!variable) continue;
1786     TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
1787     if(!name.CompareTo("pt")) {
1788       const Int_t nBinPt  = variable->GetNumberOfBins();
1789       const Double_t *kPtRange = variable->GetBinning();
1790
1791       fQACollection->CreateTH1Farray("hadronsBeforeIPcut", "Hadrons before IP cut", nBinPt, kPtRange);
1792       fQACollection->CreateTH1Farray("hadronsAfterIPcut", "Hadrons after IP cut", nBinPt, kPtRange);
1793       fQACollection->CreateTH1Farray("hadronsBeforeIPcutMC", "Hadrons before IP cut; MC p_{t}", nBinPt, kPtRange);
1794       fQACollection->CreateTH1Farray("hadronsAfterIPcutMC", "Hadrons after IP cut; MC p_{t} ", nBinPt, kPtRange);
1795
1796       fQACollection->CreateTH2Farray("Ke3Kecorr", "Ke3 decay e and K correlation; Ke3K p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
1797       fQACollection->CreateTH2Farray("Ke3K0Lecorr", "Ke3 decay e and K0L correlation; Ke3K0L p_{t}; Ke3e p_{t}; ", nBinPt, kPtRange, 20,0.,20.);
1798       fQACollection->CreateTH1Farray("Kptspectra", "Charged Kaons: MC p_{t} ", nBinPt, kPtRange);
1799       fQACollection->CreateTH1Farray("K0Lptspectra", "K0L: MC p_{t} ", nBinPt, kPtRange);
1800
1801       const Double_t kDCAbound[2] = {-5., 5.};
1802       const Double_t kDCAsigbound[2] = {-50., 50.};
1803
1804       const Int_t nDimDca=8;
1805       const Int_t nBinDca[nDimDca] = { 9, nBinPt, 2000, 2000, 12, 500, 6, 2};
1806       const Int_t nDimeDca=6;
1807       const Int_t nBineDca[nDimeDca] = { 9, nBinPt, 2000, 2000, 12, 2};
1808       Double_t minimaDca[nDimDca]  = { -1., 0., kDCAbound[0], kDCAsigbound[0], -1., 0, -1, -1.1};
1809       Double_t maximaDca[nDimDca]  = { 8., 20., kDCAbound[1], kDCAsigbound[1], 11., 50, 5, 1.1};
1810
1811       Double_t *sourceBins = AliHFEtools::MakeLinearBinning(nBinDca[0], minimaDca[0], maximaDca[0]);
1812       Double_t *dcaBins = AliHFEtools::MakeLinearBinning(nBinDca[2], minimaDca[2], maximaDca[2]);
1813       Double_t *dcaSigBins = AliHFEtools::MakeLinearBinning(nBinDca[3], minimaDca[3], maximaDca[3]);
1814       Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBinDca[4], minimaDca[4], maximaDca[4]);
1815       Double_t *eProdRBins = AliHFEtools::MakeLinearBinning(nBinDca[5], minimaDca[5], maximaDca[5]);
1816       Double_t *v0PIDBins = AliHFEtools::MakeLinearBinning(nBinDca[6], minimaDca[6], maximaDca[6]);
1817       Double_t *chargeBins = AliHFEtools::MakeLinearBinning(nBinDca[7], minimaDca[7], maximaDca[7]);
1818
1819       fQACollection->CreateTHnSparseNoLimits("Dca", "Dca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; dcasig; centrality bin; eProdR; v0pid; charge", nDimDca, nBinDca);
1820       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(0, sourceBins);
1821       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(1, kPtRange);
1822       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(2, dcaBins);
1823       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(3, dcaSigBins);
1824       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(4, centralityBins);
1825       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(5, eProdRBins);
1826       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(6, v0PIDBins);
1827       ((THnSparse*)(fQACollection->Get("Dca")))->SetBinEdges(7, chargeBins);
1828
1829       fQACollection->CreateTHnSparseNoLimits("eDca", "eDca; source (0-all, 1-charm,etc); pT [GeV/c]; dca; dcasig; centrality bin; charge", nDimeDca, nBineDca);
1830       ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(0, sourceBins);
1831       ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(1, kPtRange);
1832       ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(2, dcaBins);
1833       ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(3, dcaSigBins);
1834       ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(4, centralityBins);
1835       ((THnSparse*)(fQACollection->Get("eDca")))->SetBinEdges(5, chargeBins);
1836
1837       fQACollection->CreateTH2Farray("secondaryprotonDcaSig", "secondary proton dca significance: dca sig ", nBinPt, kPtRange, 2000, kDCAsigbound[0], kDCAsigbound[1]);
1838       fQACollection->CreateTH2Farray("protonDcaSig", "proton dca significance: dca sig ", nBinPt, kPtRange, 2000, kDCAsigbound[0], kDCAsigbound[1]);
1839       fQACollection->CreateTH2Farray("pionDcaSig", "pion dca significance: dca sig ", nBinPt, kPtRange, 2000, kDCAsigbound[0], kDCAsigbound[1]);
1840
1841       break;
1842     }
1843   }
1844
1845 }
1846
1847 //____________________________________________________________
1848 void AliAnalysisTaskHFE::InitHistoITScluster(){
1849   //
1850     // Initialize a temporary histogram to monitor the chi2/ITS cluster
1851     if(IsPbPb()) {
1852         const Int_t kNDim = 7;
1853         const Int_t kNBins[kNDim] = {88, 20,90,11, 160, 2, 1000};
1854         const Double_t kMin[kNDim] = {0.1, -1,0,  0.,0., 0,  0.};
1855         const Double_t kMax[kNDim] = {20., 1, 2.*TMath::Pi(), 11.,160, 2, 100.};
1856         fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c);eta;phi; centrality class;nclus;sharebit; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
1857         fQACollection->BinLogAxis("fChi2perITScluster", 0);
1858     }
1859     else
1860     {
1861         const Int_t kNDim = 3;
1862         const Int_t kNBins[kNDim] = {44, 11, 1000};
1863         const Double_t kMin[kNDim] = {0.1, 0., 0.};
1864         const Double_t kMax[kNDim] = {20., 11., 100.};
1865         fQACollection->CreateTHnSparse("fChi2perITScluster", "chi2/ITS cluster; p_{T} (GeV/c); centrality class; #chi^{2}/ITS cluster", kNDim, kNBins, kMin, kMax);
1866         fQACollection->BinLogAxis("fChi2perITScluster", 0);
1867     }
1868 }
1869
1870 //____________________________________________________________
1871 void AliAnalysisTaskHFE::SelectSpecialTrigger(const Char_t *trgclust, Int_t runMin, Int_t runMax){
1872   //
1873   // Select only events triggered by a special trigeer cluster
1874   //
1875   if(!fSpecialTrigger) fSpecialTrigger = new AliOADBContainer("SpecialTrigger");
1876   fSpecialTrigger->AppendObject(new TObjString(trgclust), runMin, runMax);
1877 }
1878
1879 //____________________________________________________________
1880 const Char_t * AliAnalysisTaskHFE::GetSpecialTrigger(Int_t run){
1881   //
1882   // Derive selected trigger string for given run
1883   //
1884   if(!fSpecialTrigger) return NULL;
1885   TObjString *trg = dynamic_cast<TObjString *>(fSpecialTrigger->GetObject(run));
1886   if(!trg) return NULL;
1887   return trg->String().Data();
1888 }
1889
1890 //____________________________________________________________
1891 void AliAnalysisTaskHFE::PrintStatus() const {
1892   //
1893   // Print Analysis status
1894   //
1895   printf("\n\tAnalysis Settings\n\t========================================\n\n");
1896   printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
1897   printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
1898   printf("\tDisplaced electron analysis step: %s\n", GetPlugin(kDEstep) ? "YES" : "NO");
1899   printf("\tTagged Track Analysis: %s\n", GetPlugin(kTaggedTrackAnalysis) ? "YES" : "NO");
1900   printf("\n");
1901   printf("\tParticle Identification Detectors:\n");
1902   fPID->PrintStatus();
1903   printf("\n");
1904   printf("\tQA: \n");
1905   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
1906   printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
1907   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
1908   printf("\n");
1909 }
1910
1911 //____________________________________________________________
1912 Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
1913   //
1914   // Find the production vertex of the associated MC track
1915   //
1916   if(!fMCEvent) return kFALSE;
1917   const AliVParticle *mctrack = NULL;
1918   TString objectType = track->IsA()->GetName();
1919   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1920     // Reconstructed track
1921     mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
1922   } else {
1923     // MCParticle
1924     mctrack = track;
1925   }
1926
1927   if(!mctrack) return kFALSE;
1928
1929   Double_t xv = 0.0;
1930   Double_t yv = 0.0;
1931
1932   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1933     // case MCParticle
1934     const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(mctrack);
1935     if(mcpart){
1936       xv =  mcpart->Xv();
1937       yv =  mcpart->Yv();
1938     }
1939   } else {
1940     // case AODMCParticle
1941     const AliAODMCParticle *mcpart = dynamic_cast<const AliAODMCParticle *>(mctrack);
1942     if(mcpart){
1943       xv =  mcpart->Xv();
1944       yv =  mcpart->Yv();
1945     }
1946   }
1947
1948   //printf("xv %f, yv %f\n",xv,yv);
1949   fQACollection->Fill("radius", TMath::Abs(xv),TMath::Abs(yv));
1950
1951   return kTRUE;
1952
1953 }
1954 //__________________________________________
1955 void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
1956   //
1957   // Switch on Plugin
1958   // Available:
1959   //  - Primary vertex studies
1960   //  - Secondary vertex Studies
1961   //  - Post Processing
1962   //
1963   switch(plug){
1964     case kPriVtx: SETBIT(fPlugins, plug); break;
1965     case kSecVtx: SETBIT(fPlugins, plug); break;
1966     case kIsElecBackGround: SETBIT(fPlugins, plug); break;
1967     case kPostProcess: SETBIT(fPlugins, plug); break;
1968     case kDEstep: SETBIT(fPlugins, plug); break;
1969     case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
1970     default: AliError("Unknown Plugin");
1971   };
1972 }
1973 //__________________________________________
1974 Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track){
1975   //
1976   // Check single track cuts for a given cut step
1977   // Fill the particle container
1978   //
1979   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1980   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1981   if(fVarManager->IsSignalTrack()) {
1982     fVarManager->FillContainer(fContainer, "recTrackContReco", cutStep, kFALSE);
1983     fVarManager->FillContainer(fContainer, "recTrackContMC", cutStep, kTRUE);
1984   }
1985   return kTRUE;
1986 }
1987 //___________________________________________________
1988 Bool_t AliAnalysisTaskHFE::ReadCentrality() {
1989   //
1990   // Recover the centrality of the event from ESD or AOD
1991   //
1992
1993   Float_t fCentralityLimitstemp[12];
1994   Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
1995   if(!fPbPbUserCentralityBinning) memcpy(fCentralityLimitstemp,fCentralityLimitsdefault,sizeof(fCentralityLimitsdefault));
1996   else memcpy(fCentralityLimitstemp,fCentralityLimits,sizeof(fCentralityLimitsdefault));
1997
1998
1999   Int_t bin = -1;
2000   if(IsPbPb()) {
2001     // Centrality
2002     AliCentrality *centrality = fInputEvent->GetCentrality();
2003     fCentralityPercent = centrality->GetCentralityPercentile("V0M");
2004     //printf("centrality %f\n",fCentralityPercent);
2005
2006     for(Int_t ibin = 0; ibin < 11; ibin++){
2007       if(fCentralityPercent >= fCentralityLimitstemp[ibin] && fCentralityPercent < fCentralityLimitstemp[ibin+1]){
2008         bin = ibin;
2009         //printf("test bin %f, low %f, high %f, %d\n",fCentralityPercent,fCentralityLimitstemp[ibin],fCentralityLimitstemp[ibin+1],ibin);
2010         break;
2011       }
2012     }
2013
2014     if(bin == -1) bin = 11; // Overflow
2015   } else {
2016     // PP: Tracklet multiplicity, use common definition
2017     Int_t itsMultiplicity = GetITSMultiplicity(fInputEvent);
2018     Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
2019     for(Int_t ibin = 0; ibin < 7; ibin++){
2020       if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
2021         bin = ibin;
2022         break;
2023       }
2024     }
2025     if(bin == -1) bin = 7;  // Overflow
2026   }
2027   fCentralityF = bin;
2028   AliDebug(2, Form("Centrality class %d\n", fCentralityF));
2029
2030
2031   // contributors, to be outsourced
2032   const AliVVertex *vtx;
2033   if(IsAODanalysis()){
2034     AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
2035     if(!fAOD){
2036       AliError("AOD Event required for AOD Analysis");
2037       return kFALSE;
2038     }
2039     vtx = fAOD->GetPrimaryVertex();
2040   } else {
2041     AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
2042     if(!fESD){
2043       AliError("ESD Event required for ESD Analysis");
2044       return kFALSE;
2045     }
2046     vtx = fESD->GetPrimaryVertexSPD();
2047   }
2048   if(!vtx){
2049     fContributors = 0.5;
2050     return kFALSE;
2051   }
2052   else {
2053     Int_t contributorstemp = vtx->GetNContributors();
2054     if( contributorstemp <=  0) fContributors =  0.5;
2055     else fContributors = 1.5;
2056     //printf("Number of contributors %d\n",contributorstemp);
2057   }
2058   return kTRUE;
2059 }
2060
2061 //___________________________________________________
2062 Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
2063   //
2064   // Definition of the Multiplicity according to the JPSI group (F. Kramer)
2065   //
2066   Int_t nTracklets = 0;
2067   Int_t nAcc = 0;
2068   Double_t etaRange = 1.6;
2069
2070   if (ev->IsA() == AliAODEvent::Class()) {
2071     AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
2072     nTracklets = tracklets->GetNumberOfTracklets();
2073     for (Int_t nn = 0; nn < nTracklets; nn++) {
2074       Double_t theta = tracklets->GetTheta(nn);
2075       Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
2076       if (TMath::Abs(eta) < etaRange) nAcc++;
2077     }
2078   } else if (ev->IsA() == AliESDEvent::Class()) {
2079     nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
2080     for (Int_t nn = 0; nn < nTracklets; nn++) {
2081        Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
2082       if (TMath::Abs(eta) < etaRange) nAcc++;
2083     }
2084   } else return -1;
2085
2086   return nAcc;
2087 }
2088
2089 //___________________________________________________
2090 Bool_t AliAnalysisTaskHFE::InitializeHadronBackground(Int_t run){
2091   //
2092   // Initialize background factors array from the AliOADBContainer
2093   // The container is expected to provide a TObjArray with the name
2094   // "hadronBackground" and the size 12 for the given run number
2095   //
2096   AliDebug(1, "Deriving hadronic background parameterization from OADB container");
2097   TObjArray *functions = dynamic_cast<TObjArray *>(fHadronBackgroundOADB->GetObject(run, "hadronBackground"));
2098   if(!functions){
2099     AliDebug(1, "Content in the OADB Container is not a TObjArray");
2100     fBackGroundFactorApply = kFALSE;
2101     return kFALSE;
2102   }
2103   if(functions->GetSize() < 12){
2104     AliDebug(1, Form("Size not matching: 12 expected, %d provided", functions->GetSize()));
2105     fBackGroundFactorApply = kFALSE;
2106     return kFALSE;
2107   }
2108   for(Int_t icent = 0; icent < 12; icent++) fkBackGroundFactorArray[icent] = dynamic_cast<const TF1 *>(functions->UncheckedAt(icent));
2109   return kTRUE;
2110 }
2111
2112 //___________________________________________________
2113 void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
2114   //
2115   // Recover the centrality of the event from ESD or AOD
2116   //
2117  if(IsAODanalysis()){
2118
2119    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
2120    if(!fAOD){
2121      AliError("AOD Event required for AOD Analysis");
2122        return;
2123    }
2124    // PileUp
2125    fIdentifiedAsPileUp = kFALSE;
2126    if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
2127    // Z vertex
2128    fIdentifiedAsOutInz = kFALSE;
2129    //printf("Z vertex %f and out %f\n",fAOD->GetPrimaryVertex()->GetZ(),fCuts->GetVertexRange());
2130    if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2131    // Event Cut
2132    fPassTheEventCut = kTRUE;
2133    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) fPassTheEventCut = kFALSE;
2134
2135
2136  } else {
2137
2138    AliDebug(3, "Processing ESD Centrality");
2139    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
2140    if(!fESD){
2141      AliError("ESD Event required for ESD Analysis");
2142        return;
2143    }
2144    // PileUp
2145    fIdentifiedAsPileUp = kFALSE;
2146    if(fRemovePileUp && fESD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
2147    // Z vertex
2148    fIdentifiedAsOutInz = kFALSE;
2149    if (IsPbPb()) {
2150      //printf("PbPb\n");
2151      if(fESD->GetPrimaryVertex()){
2152        if(TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2153      }
2154    }
2155    else {
2156      //printf("pp\n");
2157      if(fESD->GetPrimaryVertexTracks()){
2158        if(TMath::Abs(fESD->GetPrimaryVertexTracks()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
2159      }
2160    }
2161    //Event Cut
2162    fPassTheEventCut = kTRUE;
2163    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;
2164
2165
2166  }
2167
2168 }
2169