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