PWGHF/hfe/AliAnalysisTaskHFE.cxx

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