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