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