PWG3/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 <TCanvas.h>
  29 #include <TChain.h>
  30 #include <TDirectory.h>
  31 #include <TFile.h>
  32 #include <TH3D.h>
  33 #include <TIterator.h>
  34 #include <TList.h>
  35 #include <TLegend.h>
  36 #include <TMath.h>
  37 #include <TObjArray.h>
  38 #include <TParticle.h>
  39 #include <TProfile.h>
  40 //#include <TString.h>
  41 #include <TF1.h>
  42 #include <TTree.h>
  43
  44 #include "AliAODInputHandler.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODVertex.h"
  48 #include "AliCFContainer.h"
  49 #include "AliCFManager.h"
  50 #include "AliESDEvent.h"
  51 #include "AliESDInputHandler.h"
  52 #include "AliESDpid.h"
  53 #include "AliESDtrack.h"
  54 #include "AliCentrality.h"
  55 #include "AliLog.h"
  56 #include "AliAnalysisManager.h"
  57 #include "AliMCEvent.h"
  58 #include "AliMCEventHandler.h"
  59 #include "AliMCParticle.h"
  60 #include "AliPID.h"
  61 #include "AliStack.h"
  62 #include "AliTriggerAnalysis.h"
  63 #include "AliVVertex.h"
  64
  65 #include "AliHFEcollection.h"
  66 #include "AliHFEcontainer.h"
  67 #include "AliHFEcuts.h"
  68 #include "AliHFEelecbackground.h"
  69 #include "AliHFEmcQA.h"
  70 #include "AliHFEpairs.h"
  71 #include "AliHFEpid.h"
  72 #include "AliHFEpidQAmanager.h"
  73 #include "AliHFEpostAnalysis.h"
  74 #include "AliHFEsecVtxs.h"
  75 #include "AliHFEsecVtx.h"
  76 #include "AliHFEsignalCuts.h"
  77 #include "AliHFEtaggedTrackAnalysis.h"
  78 #include "AliHFEtools.h"
  79 #include "AliHFEvarManager.h"
  80 #include "AliAnalysisTaskHFE.h"
  81
  82 ClassImp(AliAnalysisTaskHFE)
  83
  84 //____________________________________________________________
  85 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  86   AliAnalysisTaskSE("PID efficiency Analysis")
  87   , fQAlevel(0)
  88   , fPlugins(0)
  89   , fFillSignalOnly(kTRUE)
  90   , fBackGroundFactorApply(kFALSE)
  91   , fRemovePileUp(kFALSE)
  92   , fIdentifiedAsPileUp(kFALSE)
  93   , fIdentifiedAsOutInz(kFALSE)
  94   , fPassTheEventCut(kFALSE)
  95   , fHasSpecialTriggerSelection(kFALSE)
  96   , fRejectKinkMother(kTRUE)
  97   , fSpecialTrigger("NONE")
  98   , fCentralityF(99.0)
  99   , fContributors(0.5)
 100   , fWeightBackGround(0.)
 101   , fVz(0.0)
 102   , fContainer(NULL)
 103   , fVarManager(NULL)
 104   , fSignalCuts(NULL)
 105   , fCFM(NULL)
 106   , fTriggerAnalysis(NULL)
 107   , fPID(NULL)
 108   , fPIDqa(NULL)
 109   , fPIDpreselect(NULL)
 110   , fCuts(NULL)
 111   , fTaggedTrackCuts(NULL)
 112   , fCleanTaggedTrack(kFALSE)
 113   , fVariablesTRDTaggedTrack(kFALSE)
 114   , fCutspreselect(NULL)
 115   , fSecVtx(NULL)
 116   , fElecBackGround(NULL)
 117   , fMCQA(NULL)
 118   , fTaggedTrackAnalysis(NULL)
 119   , fQA(NULL)
 120   , fOutput(NULL)
 121   , fHistMCQA(NULL)
 122   , fHistSECVTX(NULL)
 123   , fHistELECBACKGROUND(NULL)
 124   , fQACollection(NULL)
 125 {
 126   //
 127   // Dummy constructor
 128   //
 129 }
 130
 131 //____________________________________________________________
 132 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
 133   AliAnalysisTaskSE(name)
 134   , fQAlevel(0)
 135   , fPlugins(0)
 136   , fFillSignalOnly(kTRUE)
 137   , fBackGroundFactorApply(kFALSE)
 138   , fRemovePileUp(kFALSE)
 139   , fIdentifiedAsPileUp(kFALSE)
 140   , fIdentifiedAsOutInz(kFALSE)
 141   , fPassTheEventCut(kFALSE)
 142   , fHasSpecialTriggerSelection(kFALSE)
 143   , fRejectKinkMother(kTRUE)
 144   , fSpecialTrigger("NONE")
 145   , fCentralityF(99.0)
 146   , fContributors(0.5)
 147   , fWeightBackGround(0.)
 148   , fVz(0.0)
 149   , fContainer(NULL)
 150   , fVarManager(NULL)
 151   , fSignalCuts(NULL)
 152   , fCFM(NULL)
 153   , fTriggerAnalysis(NULL)
 154   , fPID(NULL)
 155   , fPIDqa(NULL)
 156   , fPIDpreselect(NULL)
 157   , fCuts(NULL)
 158   , fTaggedTrackCuts(NULL)
 159   , fCleanTaggedTrack(kFALSE)
 160   , fVariablesTRDTaggedTrack(kFALSE)
 161   , fCutspreselect(NULL)
 162   , fSecVtx(NULL)
 163   , fElecBackGround(NULL)
 164   , fMCQA(NULL)
 165   , fTaggedTrackAnalysis(NULL)
 166   , fQA(NULL)
 167   , fOutput(NULL)
 168   , fHistMCQA(NULL)
 169   , fHistSECVTX(NULL)
 170   , fHistELECBACKGROUND(NULL)
 171   , fQACollection(0x0)
 172 {
 173   //
 174   // Default constructor
 175   //
 176   DefineOutput(1, TList::Class());
 177   DefineOutput(2, TList::Class());
 178
 179   fPID = new AliHFEpid("hfePid");
 180   fPIDqa = new AliHFEpidQAmanager;
 181   fVarManager = new AliHFEvarManager("hfeVarManager");
 182
 183   memset(fElecBackgroundFactor, 0, sizeof(Double_t) * kElecBgSpecies * kBgPtBins);
 184   memset(fBinLimit, 0, sizeof(Double_t) * (kBgPtBins+1));
 185 }
 186
 187 //____________________________________________________________
 188 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
 189   AliAnalysisTaskSE(ref)
 190   , fQAlevel(0)
 191   , fPlugins(0)
 192   , fFillSignalOnly(ref.fFillSignalOnly)
 193   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
 194   , fRemovePileUp(ref.fRemovePileUp)
 195   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
 196   , fIdentifiedAsOutInz(ref.fIdentifiedAsOutInz)
 197   , fPassTheEventCut(ref.fPassTheEventCut)
 198   , fHasSpecialTriggerSelection(ref.fHasSpecialTriggerSelection)
 199   , fRejectKinkMother(ref.fRejectKinkMother)
 200   , fSpecialTrigger(ref.fSpecialTrigger)
 201   , fCentralityF(ref.fCentralityF)
 202   , fContributors(ref.fContributors)
 203   , fWeightBackGround(ref.fWeightBackGround)
 204   , fVz(ref.fVz)
 205   , fContainer(NULL)
 206   , fVarManager(NULL)
 207   , fSignalCuts(NULL)
 208   , fCFM(NULL)
 209   , fTriggerAnalysis(NULL)
 210   , fPID(NULL)
 211   , fPIDqa(NULL)
 212   , fPIDpreselect(NULL)
 213   , fCuts(NULL)
 214   , fTaggedTrackCuts(NULL)
 215   , fCleanTaggedTrack(ref.fCleanTaggedTrack)
 216   , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
 217   , fCutspreselect(NULL)
 218   , fSecVtx(NULL)
 219   , fElecBackGround(NULL)
 220   , fMCQA(NULL)
 221   , fTaggedTrackAnalysis(NULL)
 222   , fQA(NULL)
 223   , fOutput(NULL)
 224   , fHistMCQA(NULL)
 225   , fHistSECVTX(NULL)
 226   , fHistELECBACKGROUND(NULL)
 227   , fQACollection(NULL)
 228 {
 229   //
 230   // Copy Constructor
 231   //
 232   ref.Copy(*this);
 233 }
 234
 235 //____________________________________________________________
 236 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 237   //
 238   // Assignment operator
 239   //
 240   if(this == &ref)
 241     ref.Copy(*this);
 242   return *this;
 243 }
 244
 245 //____________________________________________________________
 246 void AliAnalysisTaskHFE::Copy(TObject &o) const {
 247   //
 248   // Copy into object o
 249   //
 250   AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
 251   target.fQAlevel = fQAlevel;
 252   target.fPlugins = fPlugins;
 253   target.fFillSignalOnly = fFillSignalOnly;
 254   target.fBackGroundFactorApply = fBackGroundFactorApply;
 255   target.fRemovePileUp = fRemovePileUp;
 256   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
 257   target.fIdentifiedAsOutInz = fIdentifiedAsOutInz;
 258   target.fPassTheEventCut = fPassTheEventCut;
 259   target.fHasSpecialTriggerSelection = fHasSpecialTriggerSelection;
 260   target.fRejectKinkMother = fRejectKinkMother;
 261   target.fSpecialTrigger = fSpecialTrigger;
 262   target.fCentralityF = fCentralityF;
 263   target.fContributors = fContributors;
 264   target.fWeightBackGround = fWeightBackGround;
 265   target.fVz = fVz;
 266   target.fContainer = fContainer;
 267   target.fVarManager = fVarManager;
 268   target.fSignalCuts = fSignalCuts;
 269   target.fCFM = fCFM;
 270   target.fTriggerAnalysis = fTriggerAnalysis;
 271   target.fPID = fPID;
 272   target.fPIDqa = fPIDqa;
 273   target.fPIDpreselect = fPIDpreselect;
 274   target.fCuts = fCuts;
 275   target.fTaggedTrackCuts = fTaggedTrackCuts;
 276   target.fCleanTaggedTrack = fCleanTaggedTrack;
 277   target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
 278   target.fCutspreselect = fCutspreselect;
 279   target.fSecVtx = fSecVtx;
 280   target.fElecBackGround = fElecBackGround;
 281   target.fMCQA = fMCQA;
 282   target.fTaggedTrackAnalysis = fTaggedTrackAnalysis;
 283   target.fQA = fQA;
 284   target.fOutput = fOutput;
 285   target.fHistMCQA = fHistMCQA;
 286   target.fHistSECVTX = fHistSECVTX;
 287   target.fHistELECBACKGROUND = fHistELECBACKGROUND;
 288   target.fQACollection = fQACollection;
 289 }
 290
 291 //____________________________________________________________
 292 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 293   //
 294   // Destructor
 295   //
 296   if(fPID) delete fPID;
 297   if(fVarManager) delete fVarManager;
 298   if(fPIDqa) delete fPIDqa;
 299   if(fSignalCuts) delete fSignalCuts;
 300   if(fCFM) delete fCFM;
 301   if(fSecVtx) delete fSecVtx;
 302   if(fMCQA) delete fMCQA;
 303   if(fElecBackGround) delete fElecBackGround;
 304   if(fTriggerAnalysis) delete fTriggerAnalysis;
 305   if(fPIDpreselect) delete fPIDpreselect;
 306   if(fQA) delete fQA;
 307   if(fOutput) delete fOutput;
 308 }
 309
 310 //____________________________________________________________
 311 void AliAnalysisTaskHFE::UserCreateOutputObjects(){
 312   //
 313   // Creating output container and output objects
 314   // Here we also Initialize the correction framework container and
 315   // the objects for
 316   // - PID
 317   // - MC QA
 318   // - SecVtx
 319   // QA histograms are created if requested
 320   // Called once per worker
 321   //
 322   AliDebug(3, "Creating Output Objects");
 323   // Automatic determination of the analysis mode
 324   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 325   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 326     SetAODAnalysis();
 327   } else {
 328     SetESDAnalysis();
 329     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 330       SetHasMCData();
 331   }
 332   printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
 333   printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
 334
 335   // Enable Trigger Analysis
 336   fTriggerAnalysis = new AliTriggerAnalysis;
 337   fTriggerAnalysis->EnableHistograms();
 338   fTriggerAnalysis->SetAnalyzeMC(HasMCData());
 339
 340
 341   // Make lists for Output
 342   if(!fQA) fQA = new TList;
 343   fQA->SetOwner();
 344   if(!fOutput) fOutput = new TList;
 345   fOutput->SetOwner();
 346
 347   // First Part: Make QA histograms
 348   fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
 349   fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 350   fQACollection->CreateProfile("conr", "Electron PID contamination", 20, 0, 20);
 351   fQACollection->CreateTH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi());
 352   fQACollection->CreateTH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi());
 353   fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
 354   fQACollection->CreateProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20);
 355   fQACollection->CreateProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20);
 356   fQACollection->CreateTH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20);
 357   fQACollection->CreateTH1F("mccharge", "MC Charge", 200, -100, 100);
 358   fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
 359   // Temporary histograms for TPC number of clusters for all signal tracks (MC true electrons) and for selected tracks (Markus Fasel)
 360   fQACollection->CreateTH2F("TPCclusters2_1_Signal", "TPCclusterInfo for findable clusters for 2 neighbors for signal tracks", 30, 0.1, 10., 162, 0., 161.);
 361   fQACollection->CreateTH2F("TPCclusters2_0_Signal", "TPCclusterInfo for the ratio for 2 neighbors for signal tracks", 30, 0.1, 10., 100, 0., 1.);
 362   fQACollection->CreateTH2F("TPCclusters2_1_Selected", "TPCclusterInfo for findable clusters for 2 neighbors for selected tracks", 30, 0.1, 10., 162, 0., 161.);
 363   fQACollection->CreateTH2F("TPCclusters2_0_Selected", "TPCclusterInfo for the ratio for 2 neighbors for selected tracks", 30, 0.1, 10., 110, 0., 1.1);
 364   fQACollection->CreateTH2F("TPCncls_Signal", "TPC Number of clusters for signal tracks", 30, 0.1, 10., 162, 0., 161.);
 365   fQACollection->CreateTH2F("TPCclr_Signal", "TPC cluster ratio for signal tracks", 30, 0.1, 10., 110, 0., 1.1);
 366   fQACollection->BinLogAxis("TPCclusters2_1_Signal", 0);
 367   fQACollection->BinLogAxis("TPCclusters2_0_Signal", 0);
 368   fQACollection->BinLogAxis("TPCclusters2_1_Selected", 0);
 369   fQACollection->BinLogAxis("TPCclusters2_0_Selected", 0);
 370   fQACollection->BinLogAxis("TPCncls_Signal", 0);
 371   fQACollection->BinLogAxis("TPCclr_Signal", 0);
 372   // Temporary histograms for TRD efficiency maps (Markus Fasel)
 373   fQACollection->CreateTH2F("TRDmapPosBefore", "Pos. charged tracks before TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
 374   fQACollection->CreateTH2F("TRDmapNegBefore", "Neg. charged tracks before TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
 375   fQACollection->CreateTH2F("TRDmapPosAfter", "Pos. charged tracks after TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
 376   fQACollection->CreateTH2F("TRDmapNegAfter", "Neg. charged tracks after TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
 377
 378   InitPIDperformanceQA();
 379   InitContaminationQA();
 380   fQA->Add(fQACollection->GetList());
 381
 382   // Initialize PID
 383   fPID->SetHasMCData(HasMCData());
 384   if(!fPID->GetNumberOfPIDdetectors()) fPID->AddDetector("TPC", 0);
 385   fPID->InitializePID();
 386   if(IsQAOn(kPIDqa)){
 387     AliInfo("PID QA switched on");
 388     fPIDqa->Initialize(fPID);
 389     fQA->Add(fPIDqa->MakeList("HFEpidQA"));
 390   }
 391
 392   // Initialize correction Framework and Cuts
 393   const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack + AliHFEcuts::kNcutStepsSecvtxTrack;
 394   fCFM = new AliCFManager;
 395   fCFM->SetNStepParticle(kNcutSteps);
 396   MakeParticleContainer();
 397   MakeEventContainer();
 398   // Temporary fix: Initialize particle cuts with NULL
 399   for(Int_t istep = 0; istep < kNcutSteps; istep++)
 400     fCFM->SetParticleCutsList(istep, NULL);
 401   if(!fCuts){
 402     AliWarning("Cuts not available. Default cuts will be used");
 403     fCuts = new AliHFEcuts;
 404     fCuts->CreateStandardCuts();
 405   }
 406   if(IsAODanalysis()) fCuts->SetAOD();
 407   // Make clone for V0 tagging step
 408   fCuts->Initialize(fCFM);
 409   if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
 410   fSignalCuts = new AliHFEsignalCuts("HFEsignalCuts", "HFE MC Signal definition");
 411   fVarManager->SetSignalCuts(fSignalCuts);
 412
 413   // add output objects to the List
 414   fOutput->AddAt(fContainer, 0);
 415   fOutput->AddAt(fCFM->GetEventContainer(), 1);
 416
 417   // mcQA----------------------------------
 418   if (HasMCData() && IsQAOn(kMCqa)) {
 419     AliInfo("MC QA on");
 420     if(!fMCQA) fMCQA = new AliHFEmcQA;
 421     if(!fHistMCQA) fHistMCQA = new TList();
 422     fHistMCQA->SetOwner();
 423     fMCQA->CreatDefaultHistograms(fHistMCQA);
 424     if(!fFillSignalOnly) fMCQA->SetBackgroundWeightFactor(fElecBackgroundFactor[0],fBinLimit);
 425     fQA->Add(fHistMCQA);
 426   }
 427
 428   // secvtx----------------------------------
 429   if (GetPlugin(kSecVtx)) {
 430     AliInfo("Secondary Vertex Analysis on");
 431     if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
 432     fSecVtx->SetHasMCData(HasMCData());
 433
 434     if(!fHistSECVTX) fHistSECVTX = new TList();
 435     fHistSECVTX->SetOwner();
 436     fSecVtx->CreateHistograms(fHistSECVTX);
 437     fOutput->Add(fHistSECVTX);
 438   }
 439
 440   // background----------------------------------
 441   if (GetPlugin(kIsElecBackGround)) {
 442     AliInfo("Electron BackGround Analysis on");
 443     if(!fElecBackGround){
 444       AliWarning("ElecBackGround not available. Default elecbackground will be used");
 445       fElecBackGround = new AliHFEelecbackground;
 446     }
 447     fElecBackGround->SetHasMCData(HasMCData());
 448
 449     if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
 450     fHistELECBACKGROUND->SetOwner();
 451     fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
 452     fOutput->Add(fHistELECBACKGROUND);
 453   }
 454
 455   // tagged tracks
 456   if(GetPlugin(kTaggedTrackAnalysis)){
 457     AliInfo("Analysis on V0-tagged tracks enabled");
 458     fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis(Form("taggedTrackAnalysis%s", GetName()));
 459     fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
 460     fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
 461     AliHFEvarManager *varManager = fTaggedTrackAnalysis->GetVarManager();
 462     TObjArray *array = fVarManager->GetVariables();
 463     Int_t nvars = array->GetEntriesFast();
 464     TString namee;
 465     for(Int_t v = 0; v < nvars; v++) {
 466       AliHFEvarManager::AliHFEvariable *variable = (AliHFEvarManager::AliHFEvariable *) array->At(v);
 467       if(!variable) continue;
 468       TString name(((AliHFEvarManager::AliHFEvariable *)variable)->GetName());
 469       if(!name.CompareTo("source")) namee = TString("species");
 470       else namee = TString(name);
 471       varManager->AddVariable(namee);
 472       //printf("For AliTaggedTrackAnalysis, had the variable %s and the one used %s\n",(const char*)variable->GetName(),(const char*) namee);
 473     }
 474     if(fPIDqa->HasHighResolutionHistos())
 475       fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
 476     fTaggedTrackAnalysis->SetPID(fPID);
 477     fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
 478     fTaggedTrackAnalysis->InitContainer();
 479     fOutput->Add(fTaggedTrackAnalysis->GetContainer());
 480     fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
 481     fQA->Add(fTaggedTrackAnalysis->GetCutQA());
 482     fQA->Add(fTaggedTrackAnalysis->GetQAcollection());
 483   }
 484   PrintStatus();
 485 }
 486
 487 //____________________________________________________________
 488 void AliAnalysisTaskHFE::UserExec(Option_t *){
 489   //
 490   // Run the analysis
 491   //
 492   AliDebug(3, "Starting Single Event Analysis");
 493   if(!fInputEvent){
 494     AliError("Reconstructed Event not available");
 495     return;
 496   }
 497   if(HasMCData()){
 498     AliDebug(4, Form("MC Event: %p", fMCEvent));
 499     if(!fMCEvent){
 500       AliError("No MC Event, but MC Data required");
 501       return;
 502     }
 503   }
 504   if(!fCuts){
 505     AliError("HFE cuts not available");
 506     return;
 507   }
 508
 509
 510   if(IsESDanalysis() && HasMCData()){
 511     // Protect against missing MC trees
 512     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 513     if(!mcH){
 514       AliError("No MC Event Handler available");
 515       return;
 516     }
 517     if(!mcH->InitOk()) return;
 518     if(!mcH->TreeK()) return;
 519     if(!mcH->TreeTR()) return;
 520   }
 521
 522   // need the centrality for everything (MC also)
 523   fCentralityF = -100.0;
 524   if(!ReadCentrality()) fCentralityF = -100.0;
 525   //printf("pass centrality\n");
 526   //printf("Reading fCentralityF %f\n",fCentralityF);
 527
 528   // See if pile up and z in the range
 529   RejectionPileUpVertexRangeEventCut();
 530
 531   // Protect agains missing
 532   if(HasMCData()){
 533     //printf("Has MC data\n");
 534     fSignalCuts->SetMCEvent(fMCEvent);
 535     ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
 536   }
 537
 538   if(IsAODanalysis()){
 539     AliAODpidUtil *aodworkingpid = AliHFEtools::GetDefaultAODPID(HasMCData());
 540     fPID->SetAODpid(aodworkingpid);
 541     ProcessAOD();
 542   } else {
 543     // Check Trigger selection
 544     if(fHasSpecialTriggerSelection){
 545       AliESDEvent *ev = dynamic_cast<AliESDEvent *>(fInputEvent);
 546       if(!(ev && ev->IsTriggerClassFired(fSpecialTrigger.Data()))) return;
 547     }
 548     AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
 549     if(!inH){
 550       AliError("No ESD Input handler available");
 551       return;
 552     }
 553     AliESDpid *workingPID = inH->GetESDpid();
 554     if(!workingPID){
 555       AliDebug(1, "Using default ESD PID");
 556       workingPID = AliHFEtools::GetDefaultPID(HasMCData());
 557     } else {
 558       AliDebug(1, "Using ESD PID from the input handler");
 559     }
 560     fPID->SetESDpid(workingPID);
 561     if(fPIDpreselect) fPIDpreselect->SetESDpid(workingPID);
 562
 563     ProcessESD();
 564   }
 565   // Done!!!
 566   PostData(1, fOutput);
 567   PostData(2, fQA);
 568 }
 569
 570 //____________________________________________________________
 571 void AliAnalysisTaskHFE::Terminate(Option_t *){
 572   //
 573   // Terminate not implemented at the moment
 574   //
 575   if(GetPlugin(kPostProcess)){
 576     fOutput = dynamic_cast<TList *>(GetOutputData(1));
 577     fQA = dynamic_cast<TList *>(GetOutputData(2));
 578     if(!fOutput){
 579       AliError("Results not available");
 580       return;
 581     }
 582     if(!fQA){
 583       AliError("QA output not available");
 584       return;
 585     }
 586     fContainer = dynamic_cast<AliHFEcontainer *>(fOutput->FindObject("trackContainer"));
 587     if(!fContainer){
 588       AliError("Track container not found");
 589       return;
 590     }
 591     AliHFEpostAnalysis postanalysis;
 592     postanalysis.SetTaskResults(fContainer);
 593     TList *qalist = dynamic_cast<TList *>(fQA->FindObject("list_TaskQA"));
 594     if(!qalist){
 595       AliError("QA List not found");
 596       return;
 597     }
 598     postanalysis.SetTaskQA(qalist);
 599     printf("Running post analysis\n");
 600     //if(HasMCData())
 601     postanalysis.DrawMCSignal2Background();
 602     postanalysis.DrawEfficiency();
 603     postanalysis.DrawPIDperformance();
 604     postanalysis.DrawCutEfficiency();
 605
 606     if (GetPlugin(kIsElecBackGround)) {
 607       AliHFEelecbackground elecBackGround;
 608       TList *oe = 0x0;
 609       if(!(oe = (TList*)dynamic_cast<TList *>(fOutput->FindObject("HFEelecbackground")))){
 610         return;
 611       }
 612       elecBackGround.Load(oe);
 613       elecBackGround.Plot();
 614       elecBackGround.PostProcess();
 615     }
 616   }
 617 }
 618 //_______________________________________________________________
 619 Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
 620   //
 621   //
 622   //
 623
 624   //printf("test in IsEventInBinZero\n");
 625   if(!fInputEvent){
 626     AliError("Reconstructed Event not available");
 627     return kFALSE;
 628   }
 629
 630   // check vertex
 631   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
 632   if(!vertex) return kTRUE;
 633   //if(vertex) return kTRUE;
 634
 635   // check tracks
 636   if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
 637   //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
 638
 639
 640   return kFALSE;
 641
 642 }
 643 //____________________________________________________________
 644 void AliAnalysisTaskHFE::ProcessMC(){
 645   //
 646   // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
 647   // In case MC QA is on also MC QA loop is done
 648   //
 649   AliDebug(3, "Processing MC Information");
 650   Double_t eventContainer [4];
 651   eventContainer[0] = fMCEvent->GetPrimaryVertex()->GetZ();
 652   eventContainer[2] = fCentralityF;
 653   eventContainer[3] = fContributors;
 654   fVz = eventContainer[0];
 655   //printf("z position is %f\n",eventContainer[0]);
 656   //if(fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent))
 657   fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
 658   Int_t nElectrons = 0;
 659   if(IsESDanalysis()){
 660    if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
 661     if (HasMCData() && IsQAOn(kMCqa)) {
 662       AliDebug(2, "Running MC QA");
 663
 664       if(fMCEvent->Stack()){
 665               fMCQA->SetMCEvent(fMCEvent);
 666         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
 667         fMCQA->Init();
 668
 669         fMCQA->GetMesonKine();
 670
 671         // loop over all tracks for decayed electrons
 672         for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
 673           TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
 674           if(!mcpart) continue;
 675           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
 676           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
 677           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
 678           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
 679           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
 680           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 681           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 682           if (TMath::Abs(mcpart->Eta()) < 0.9) {
 683             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 684             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 685             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 686           }
 687           if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) {
 688             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 689             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 690             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 691           }
 692         }
 693         fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 694         fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 695       }
 696
 697     } // end of MC QA loop
 698    }
 699    // -----------------------------------------------------------------
 700    fCFM->SetMCEventInfo(fMCEvent);
 701    // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 702   } else {
 703     fCFM->SetMCEventInfo(fInputEvent);
 704   }
 705   // Run MC loop
 706   AliVParticle *mctrack = NULL;
 707   AliDebug(3, Form("Number of Tracks: %d", fMCEvent->GetNumberOfTracks()));
 708   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfTracks(); imc++){
 709     if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
 710     AliDebug(4, "Next MC Track");
 711     if(ProcessMCtrack(mctrack)) nElectrons++;
 712   }
 713
 714   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 715   fQACollection->Fill("nElectron", nElectrons);
 716 }
 717
 718 //____________________________________________________________
 719 void AliAnalysisTaskHFE::ProcessESD(){
 720   //
 721   // Run Analysis of reconstructed event in ESD Mode
 722   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
 723   //
 724   AliDebug(3, "Processing ESD Event");
 725   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 726   if(!fESD){
 727     AliError("ESD Event required for ESD Analysis");
 728     return;
 729   }
 730
 731   // Set magnetic field if V0 task on
 732   if(fTaggedTrackAnalysis) {
 733     fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
 734     fTaggedTrackAnalysis->SetCentrality(fCentralityF);
 735   }
 736
 737   // Do event Normalization
 738   Double_t eventContainer[4];
 739   eventContainer[0] = 0.0;
 740   if(HasMCData()) eventContainer[0] = fVz;
 741   else {
 742     if(IsPbPb()) {
 743       //printf("Z PbPb\n");
 744       if(fESD->GetPrimaryVertexSPD()) eventContainer[0] = fESD->GetPrimaryVertexSPD()->GetZ();
 745     }
 746     else {
 747       //printf("z pp\n");
 748       if(fESD->GetPrimaryVertexTracks()) eventContainer[0] = fESD->GetPrimaryVertexTracks()->GetZ();
 749     }
 750   }
 751   eventContainer[1] = 0.;
 752   eventContainer[2] = fCentralityF;
 753   eventContainer[3] = fContributors;
 754   if(fTriggerAnalysis->IsOfflineTriggerFired(fESD, AliTriggerAnalysis::kV0AND))
 755     eventContainer[1] = 1.;
 756
 757   //
 758   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
 759
 760   //
 761   if(fIdentifiedAsPileUp) return;
 762   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
 763
 764   //
 765   if(TMath::Abs(fCentralityF+100.0) < 0.01) return;
 766   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
 767   //printf("In ProcessESD %f\n",fCentralityF);
 768
 769   //
 770   if(fIdentifiedAsOutInz) return;
 771   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
 772
 773   //
 774   if(!fPassTheEventCut) return;
 775   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
 776
 777
 778
 779   fContainer->NewEvent();
 780
 781   if (GetPlugin(kIsElecBackGround)) {
 782     fElecBackGround->SetEvent(fESD);
 783   }
 784   if (GetPlugin(kSecVtx)) {
 785     fSecVtx->SetEvent(fESD);
 786     fSecVtx->GetPrimaryCondition();
 787   }
 788
 789   if(HasMCData()){
 790     if (GetPlugin(kSecVtx)) {
 791       fSecVtx->SetMCEvent(fMCEvent);
 792       fSecVtx->SetMCQA(fMCQA);
 793     }
 794     if (GetPlugin(kIsElecBackGround)) {
 795       fElecBackGround->SetMCEvent(fMCEvent);
 796     }
 797   }
 798
 799   Double_t container[10];
 800   memset(container, 0, sizeof(Double_t) * 10);
 801   // container for the output THnSparse
 802   Double_t dataE[6]; // [pT, eta, Phi, type, 'C' or 'B']
 803   Int_t nElectronCandidates = 0;
 804   AliESDtrack *track = NULL, *htrack = NULL;
 805   AliMCParticle *mctrack = NULL;
 806   Int_t pid = 0;
 807
 808   Bool_t signal = kTRUE;
 809
 810   fCFM->SetRecEventInfo(fESD);
 811   // Electron background analysis
 812   if (GetPlugin(kIsElecBackGround)) {
 813
 814     AliDebug(2, "Running BackGround Analysis");
 815
 816     fElecBackGround->Reset();
 817
 818   } // end of electron background analysis
 819   //
 820   // Loop ESD
 821   //
 822   AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
 823   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 824     AliDebug(4, "New ESD track");
 825     track = fESD->GetTrack(itrack);
 826     track->SetESDEvent(fESD);
 827
 828     // fill counts of v0-identified particles
 829     Int_t v0pid = -1;
 830     if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
 831     else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
 832     else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
 833     // here the tagged track analysis will run
 834     if(fTaggedTrackAnalysis && v0pid > -1){
 835       AliDebug(1, Form("Track identified as %s", AliPID::ParticleName(v0pid)));
 836       fTaggedTrackAnalysis->ProcessTrack(track, v0pid);
 837     }
 838
 839     AliDebug(3, Form("Doing track %d, %p", itrack, track));
 840
 841     //////////////////////////////////////
 842     // preselect
 843     /////////////////////////////////////
 844     if(fPIDpreselect && fCutspreselect) {
 845       if(!PreSelectTrack(track)) continue;
 846     }
 847
 848     signal = kTRUE;
 849
 850     // Fill step without any cut
 851
 852     if(HasMCData()){
 853       // Check if it is electrons near the vertex
 854       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
 855
 856       if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 857       else AliDebug(3, "Signal Electron");
 858     }
 859     // Cache new Track information inside the var manager
 860     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
 861
 862     if(signal) {
 863       fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
 864       fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
 865       if((track->GetStatus() & AliESDtrack::kTPCout)
 866           && (TMath::Abs(track->Eta()) < 0.8)
 867           && (track->GetKinkIndex(0) == 0)){
 868         fQACollection->Fill("TPCclusters2_1_Signal", track->Pt(), track->GetTPCClusterInfo(2,1));
 869         fQACollection->Fill("TPCclusters2_0_Signal", track->Pt(), track->GetTPCNclsF() > 0 ?  track->GetTPCClusterInfo(2,1)/track->GetTPCNclsF() : 0.);
 870         fQACollection->Fill("TPCncls_Signal", track->Pt(), track->GetTPCNcls());
 871         fQACollection->Fill("TPCclr_Signal", track->Pt(), track->GetTPCNclsF() > 0 ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 0.);
 872       }
 873     }
 874
 875     // RecKine: ITSTPC cuts
 876     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
 877
 878     // Check TRD criterions (outside the correction framework)
 879     if(track->GetTRDncls()){
 880       fQACollection->Fill("chi2TRD", track->GetTRDchi2()/track->GetTRDncls());
 881       fQACollection->Fill("alpha_rec", track->GetAlpha());
 882       fQACollection->Fill("pidquality", container[0], track->GetTRDpidQuality());
 883       fQACollection->Fill("ntrdclusters", container[0], track->GetTRDncls());
 884     }
 885
 886
 887     // RecPrim
 888     if(fRejectKinkMother) {
 889       if(track->GetKinkIndex(0) != 0) continue; } // Quick and dirty fix to reject both kink mothers and daughters
 890     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
 891
 892     // Temporary eta-phi maps for TRD (Markus Fasel)
 893     Int_t minTrackletsTRD = fCuts->GetMinTrackletsTRD();
 894     if(minTrackletsTRD && track->Pt() > 2.){
 895       if(track->Charge() > 0){
 896         // positive charge
 897         fQACollection->Fill("TRDmapPosBefore", track->Eta(), track->Phi());
 898         if(track->GetTRDntrackletsPID() >= minTrackletsTRD) fQACollection->Fill("TRDmapPosAfter",track->Eta(), track->Phi());
 899       } else {
 900         // negative charge
 901         fQACollection->Fill("TRDmapNegBefore", track->Eta(), track->Phi());
 902         if(track->GetTRDntrackletsPID() >= minTrackletsTRD) fQACollection->Fill("TRDmapNegAfter",track->Eta(), track->Phi());
 903       }
 904     }
 905
 906     // HFEcuts: ITS layers cuts
 907     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue;
 908
 909     // HFE cuts: TOF PID and mismatch flag
 910     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) continue;
 911
 912     // HFEcuts: Nb of tracklets TRD0
 913     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
 914
 915     // Fill correlation maps before PID
 916     if(signal && fContainer->GetCorrelationMatrix("correlationstepbeforePID")) {
 917       //printf("Fill correlation maps before PID\n");
 918       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepbeforePID"));
 919     }
 920
 921     if(HasMCData()){
 922       FillProductionVertex(track);
 923     }
 924
 925     // track accepted, do PID
 926     AliHFEpidObject hfetrack;
 927     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
 928     hfetrack.SetRecTrack(track);
 929     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
 930     hfetrack.SetCentrality(fCentralityF);
 931     fPID->SetVarManager(fVarManager);
 932     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
 933     nElectronCandidates++;
 934     fQACollection->Fill("TPCclusters2_1_Selected", track->Pt(), track->GetTPCClusterInfo(2,1));
 935     fQACollection->Fill("TPCclusters2_0_Selected", track->Pt(), track->GetTPCClusterInfo(2,0));
 936
 937     // Fill Histogram for Hadronic Background
 938     if(HasMCData()){
 939       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
 940         fVarManager->FillContainer(fContainer, "hadronicBackground", UInt_t(0), kFALSE);
 941     }
 942
 943     // Fill Containers
 944     if(signal) {
 945       // Apply weight for background contamination
 946             if(fBackGroundFactorApply==kTRUE) {
 947               if(IsPbPb()) fWeightBackGround =  fBackGroundFactorArray[(Int_t)fCentralityF]->Eval(TMath::Abs(track->P()));
 948               else    fWeightBackGround =  fBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
 949
 950               if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
 951               else if(fWeightBackGround > 1.0) fWeightBackGround = 1.0;
 952         // weightBackGround as special weight
 953         fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, fWeightBackGround);
 954       }
 955       fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterPID"));
 956     }
 957
 958     Bool_t bTagged=kFALSE;
 959     if(GetPlugin(kSecVtx)) {
 960       AliDebug(2, "Running Secondary Vertex Analysis");
 961       if(fSecVtx->Process(track) && signal) {
 962         fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
 963         fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
 964         bTagged=kTRUE;
 965       }
 966     }
 967
 968     if(HasMCData()){
 969       dataE[0] = track->Pt();
 970       dataE[1] = track->Eta();
 971       dataE[2] = track->Phi();
 972       dataE[3] = track->Charge();
 973       dataE[4] = -1.;
 974       dataE[5] = -1.;
 975
 976       // Track selected: distinguish between true and fake
 977       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
 978       if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
 979         Int_t type = 0;
 980         if(fSignalCuts->IsCharmElectron(track))
 981           type = 1;
 982         else if(fSignalCuts->IsBeautyElectron(track))
 983           type = 2;
 984         AliDebug(1, Form("Type: %d\n", type));
 985         if(type){
 986                 dataE[5] = type; // beauty[1] or charm[2]
 987                 dataE[4] = 2;  // signal electron
 988         }
 989         else{
 990                 dataE[4] = 1; // not a signal electron
 991                 dataE[5] = 0;
 992         }
 993       }
 994       else {
 995         // Fill THnSparse with the information for Fake Electrons
 996         dataE[4] = 0;
 997         dataE[5] = 0;
 998       }
 999       // fill the performance THnSparse, if the mc origin could be defined
1000       if(dataE[4] > -1){
1001         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1002         fQACollection->Fill("PIDperformance", dataE);
1003       }
1004     }
1005     // Electron background analysis
1006     if (GetPlugin(kIsElecBackGround)) {
1007
1008       AliDebug(2, "Running BackGround Analysis");
1009
1010       for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
1011         htrack = fESD->GetTrack(jtrack);
1012         if ( itrack == jtrack ) continue;
1013         fElecBackGround->PairAnalysis(track, htrack);
1014       }
1015     } // end of electron background analysis
1016
1017     if (GetPlugin(kDEstep)) {
1018       Double_t weightElecBg = 0.;
1019       if(HasMCData()){
1020         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1021           fQACollection->Fill("hadronsBeforeIPcut",track->Pt());
1022           fQACollection->Fill("hadronsBeforeIPcutMC",mctrack->Pt());
1023         }
1024         if(fMCQA && !fFillSignalOnly) {
1025            weightElecBg = fMCQA->GetWeightFactor(mctrack); // positive:conversion e, negative: nonHFE
1026            if(weightElecBg>0) fVarManager->FillContainer(fContainer, "conversionElecs", 0, kFALSE, weightElecBg);
1027            else if(weightElecBg<0) fVarManager->FillContainer(fContainer, "mesonElecs", 0, kFALSE, -1*weightElecBg);
1028            if(bTagged){ // bg estimation for the secondary vertex tagged signals
1029              if(weightElecBg>0) fVarManager->FillContainer(fContainer, "conversionElecs", 2, kFALSE, weightElecBg);
1030              else if(weightElecBg<0) fVarManager->FillContainer(fContainer, "mesonElecs", 2, kFALSE, -1*weightElecBg);
1031            }
1032         }
1033       }
1034       // Fill Containers for impact parameter analysis
1035       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsDca + AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack,track)) continue;
1036       if(HasMCData()){
1037         if(fMCQA && !fFillSignalOnly) {
1038            if(weightElecBg>0) fVarManager->FillContainer(fContainer, "conversionElecs", 1, kFALSE, weightElecBg);
1039            else if(weightElecBg<0) fVarManager->FillContainer(fContainer, "mesonElecs", 1, kFALSE, -1*weightElecBg);
1040         }
1041       }
1042       if(signal) {
1043         fVarManager->FillContainer(fContainer, "recTrackContDEReco", AliHFEcuts::kStepHFEcutsDca, kFALSE);
1044         fVarManager->FillContainer(fContainer, "recTrackContDEMC", AliHFEcuts::kStepHFEcutsDca, kTRUE);
1045         fVarManager->FillCorrelationMatrix(fContainer->GetCorrelationMatrix("correlationstepafterDE"));
1046       }
1047       if(HasMCData()){
1048         if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11)){
1049           fQACollection->Fill("hadronsAfterIPcut",track->Pt());
1050           fQACollection->Fill("hadronsAfterIPcutMC",mctrack->Pt());
1051         }
1052       }
1053     }
1054
1055   }
1056   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1057 }
1058
1059 //____________________________________________________________
1060 void AliAnalysisTaskHFE::ProcessAOD(){
1061   //
1062   // Run Analysis in AOD Mode
1063   // Function is still in development
1064   //
1065   AliDebug(3, "Processing AOD Event");
1066   Double_t eventContainer[4];
1067   if(HasMCData()) eventContainer[0] = fVz;
1068   else {
1069     eventContainer[0] = fInputEvent->GetPrimaryVertex()->GetZ();
1070   }
1071   eventContainer[1] = 1.; // No Information available in AOD analysis, assume all events have V0AND
1072   eventContainer[2] = fCentralityF;
1073   eventContainer[3] = fContributors;
1074
1075   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1076   if(!fAOD){
1077     AliError("AOD Event required for AOD Analysis");
1078     return;
1079   }
1080
1081   //
1082   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoCut);
1083
1084   //
1085   if(fIdentifiedAsPileUp) return;
1086   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
1087
1088   //
1089   if(fIdentifiedAsOutInz) return;
1090   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepZRange);
1091
1092   //
1093   if(!fPassTheEventCut) return;
1094   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepReconstructed);
1095
1096   fContainer->NewEvent();
1097
1098   AliAODTrack *track = NULL;
1099   AliAODMCParticle *mctrack = NULL;
1100   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
1101   Int_t nElectronCandidates = 0;
1102   Int_t pid;
1103   Bool_t signal;
1104   for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
1105     track = fAOD->GetTrack(itrack); mctrack = NULL;
1106     if(!track) continue;
1107     if(track->GetFlags() != 1<<4) continue;  // Only process AOD tracks where the HFE is set
1108
1109     signal = kTRUE;
1110     if(HasMCData()){
1111
1112       Int_t label = TMath::Abs(track->GetLabel());
1113       if(label)
1114         mctrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(label));
1115         if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
1116     }
1117     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, kTRUE);
1118     // track accepted, do PID
1119     AliHFEpidObject hfetrack;
1120     hfetrack.SetAnalysisType(AliHFEpidObject::kAODanalysis);
1121     hfetrack.SetRecTrack(track);
1122     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
1123     hfetrack.SetCentrality(fCentralityF);
1124     fPID->SetVarManager(fVarManager);
1125     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;    // we will do PID here as soon as possible
1126     // Apply weight for background contamination
1127     Double_t weightBackGround = 1.0;
1128
1129     // not correct treatment for pp
1130     if(fBackGroundFactorApply==kTRUE) {
1131             if(IsPbPb()) weightBackGround =  fBackGroundFactorArray[(Int_t)fCentralityF]->Eval(TMath::Abs(track->P()));
1132       else weightBackGround =  fBackGroundFactorArray[0]->Eval(TMath::Abs(track->P()));
1133
1134       if(weightBackGround < 0.0) weightBackGround = 0.0;
1135             else if(weightBackGround > 1.0) weightBackGround = 1.0;
1136             fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, weightBackGround);
1137     }
1138
1139
1140
1141     nElectronCandidates++;
1142     if(HasMCData()){
1143       dataE[0] = track->Pt();
1144       dataE[1] = track->Eta();
1145       dataE[2] = track->Phi();
1146       dataE[3] = track->Charge();
1147       dataE[4] = -1;
1148       dataE[5] = -1;
1149       // Track selected: distinguish between true and fake
1150       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack ? mctrack->GetPdgCode(): -1));
1151       if(mctrack && ((pid = TMath::Abs(mctrack->GetPdgCode())) == 11)){
1152
1153         Int_t type = 0;
1154         if(fSignalCuts->IsCharmElectron(track))
1155           type = 1;
1156         else if(fSignalCuts->IsBeautyElectron(track))
1157           type = 2;
1158         AliDebug(1, Form("Type: %d\n", type));
1159         if(type){
1160                 dataE[5] = type; // beauty[1] or charm[2]
1161                 dataE[4] = 2;  // signal electron
1162         }
1163         else{
1164                 dataE[4] = 1; // not a signal electron
1165                 dataE[5] = 0;
1166         }
1167       }
1168       else {
1169         // Fill THnSparse with the information for Fake Electrons
1170         dataE[4] = 0;
1171         dataE[5] = 0;
1172       }
1173       // fill the performance THnSparse, if the mc origin could be defined
1174       if(dataE[4] > -1){
1175         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
1176         fQACollection->Fill("PIDperformance", dataE);
1177       }
1178     }
1179   }
1180   fQACollection->Fill("nElectronTracksEvent", nElectronCandidates);
1181 }
1182
1183 //____________________________________________________________
1184 Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
1185   //
1186   // Filter the Monte Carlo Track
1187   // Additionally Fill a THnSparse for Signal To Background Studies
1188   // Works for AOD and MC analysis Type
1189   //
1190   fVarManager->NewTrack(track, NULL, fCentralityF, -1, kTRUE);
1191   Double_t signalContainer[6];
1192
1193   signalContainer[0] = track->Pt();
1194   signalContainer[1] = track->Eta();
1195   signalContainer[2] = track->Phi();
1196   signalContainer[3] = track->Charge()/3;
1197
1198  Double_t vertex[3]; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
1199   if(IsESDanalysis()){
1200     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
1201     if(mctrack){
1202       vertex[0] = mctrack->Particle()->Vx();
1203       vertex[1] = mctrack->Particle()->Vy();
1204     }
1205   } else {
1206     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
1207     if(aodmctrack) aodmctrack->XvYvZv(vertex);
1208   }
1209
1210   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
1211   fQACollection->Fill("mccharge", signalContainer[3]);
1212   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
1213   signalContainer[4] = 0;
1214   if(fSignalCuts->IsSelected(track)){
1215     //fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCsignal, kFALSE);
1216     // Filling of the Signal/Background histogram using the
1217     // definition of the codes for charm and beauty as below in
1218     // th crearion of the histogram
1219     if(fSignalCuts->IsCharmElectron(track))
1220       signalContainer[4] = 1;
1221     else
1222       signalContainer[4] = 2;
1223   } else {
1224     signalContainer[4] = 0; // (and other background)
1225   }
1226   signalContainer[5] = 0;
1227   // apply cut on the sqrt of the production vertex
1228   Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
1229   if(radVertex < 3.5){
1230     // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
1231     signalContainer[5] = 1;
1232   } else if (radVertex < 7.5){
1233     signalContainer[5] = 2;
1234   }
1235   fQACollection->Fill("SignalToBackgroundMC", signalContainer);
1236   fQACollection->Fill("alpha_sim", track->Phi() - TMath::Pi());
1237
1238   // Step GeneratedZOutNoPileUp
1239   if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01)) return kFALSE;
1240   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
1241   //printf("In ProcessMCtrack %f\n",fCentralityF);
1242
1243   // Step Generated Event Cut
1244   if(!fPassTheEventCut) return kFALSE;
1245   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedEventCut, kFALSE);
1246
1247   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
1248   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCInAcceptance, kFALSE);
1249   return kTRUE;
1250 }
1251
1252 //____________________________________________________________
1253 Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
1254   //
1255   // Preselect tracks
1256   //
1257
1258
1259   Bool_t survived = kTRUE;
1260
1261   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) {
1262     survived = kFALSE;
1263     //printf("Did not pass AliHFEcuts::kStepRecKineITSTPC\n");
1264   }
1265   //else printf("Pass AliHFEcuts::kStepRecKineITSTPC\n");
1266   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) {
1267     survived = kFALSE;
1268     //printf("Did not pass AliHFEcuts::kStepRecPrim\n");
1269   }
1270   //else printf("Pass AliHFEcuts::kStepRecPrim\n");
1271   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) {
1272     survived = kFALSE;
1273     //printf("Did not pass AliHFEcuts::kStepHFEcutsITS\n");
1274   }
1275   //else printf("Pass AliHFEcuts::kStepHFEcutsITS\n");
1276   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTOF, track)) {
1277     survived = kFALSE;
1278     //printf("Did not pass AliHFEcuts::kStepHFEcutsTOF\n");
1279   }
1280   //else printf("Pass AliHFEcuts::kStepHFEcutsTOF\n");
1281   if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) {
1282     survived = kFALSE;
1283     //printf("Did not pass AliHFEcuts::kStepHFEcutsTRD\n");
1284   }
1285   //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
1286
1287   if(survived){
1288     // Apply PID
1289     AliHFEpidObject hfetrack;
1290     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
1291     hfetrack.SetRecTrack(track);
1292     if(!fPIDpreselect->IsSelected(&hfetrack)) {
1293       //printf("Did not pass AliHFEcuts::kPID\n");
1294       survived = kFALSE;
1295     }
1296     //else printf("Pass AliHFEcuts::kPID\n");
1297   }
1298
1299   return survived;
1300
1301 }
1302 //____________________________________________________________
1303 void AliAnalysisTaskHFE::MakeEventContainer(){
1304   //
1305   // Create the event container for the correction framework and link it
1306   // 1st bin: Vertex z-position
1307   // 2nd bin: V0AND decision (normalization to sigma_inel)
1308   // 3rd bin: Centrality class (for pp defined as number of contributors in vertex.)
1309   //
1310
1311   if(IsPbPb()) {
1312
1313     //printf("This is PbPb!!!!!!!!!!!\n");
1314
1315     const Int_t kNvar = 4;  // number of variables on the grid:
1316     Int_t nBins[kNvar] = {120, 2, 11, 2};
1317     Double_t binMin[kNvar] = {-30. , 0., 0.0, 0.};
1318     Double_t binMax[kNvar] = {30., 2., 11.0, 2.};
1319
1320     AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1321
1322     Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1323     Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1324     Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1325     Double_t *contributorsBins = AliHFEtools::MakeLinearBinning(nBins[3], binMin[3], binMax[3]);
1326     evCont->SetBinLimits(0, vertexBins);
1327     evCont->SetBinLimits(1, v0andBins);
1328     evCont->SetBinLimits(2, centralityBins);
1329     evCont->SetBinLimits(3, contributorsBins);
1330     delete[] vertexBins; delete[] v0andBins; delete[] centralityBins; delete[] contributorsBins;
1331
1332     fCFM->SetEventContainer(evCont);
1333   }
1334   else {
1335
1336     //printf("This is pp!!!!!!!!!!!\n");
1337
1338     const Int_t kNvar = 3;  // number of variables on the grid:
1339     Int_t nBins[kNvar] = {120, 2, 11};
1340     Double_t binMin[kNvar] = {-30. , 0., 0.0};
1341     Double_t binMax[kNvar] = {30., 2., 11.0};
1342
1343     AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, nBins);
1344
1345     Double_t *vertexBins = AliHFEtools::MakeLinearBinning(nBins[0], binMin[0], binMax[0]);
1346     Double_t *v0andBins = AliHFEtools::MakeLinearBinning(nBins[1], binMin[1], binMax[1]);
1347     Double_t *centralityBins = AliHFEtools::MakeLinearBinning(nBins[2], binMin[2], binMax[2]);
1348     evCont->SetBinLimits(0, vertexBins);
1349     evCont->SetBinLimits(1, v0andBins);
1350     evCont->SetBinLimits(2, centralityBins);
1351     delete[] vertexBins; delete[] v0andBins; delete[] centralityBins;
1352
1353     fCFM->SetEventContainer(evCont);
1354   }
1355
1356 }
1357
1358 //____________________________________________________________
1359 void AliAnalysisTaskHFE::MakeParticleContainer(){
1360   //
1361   // Create the particle container for the correction framework manager and
1362   // link it
1363   //
1364
1365   if(!fContainer) fContainer = new AliHFEcontainer("trackContainer");
1366   fVarManager->DefineVariables(fContainer);
1367
1368   // Create Correction Framework containers
1369   fContainer->CreateContainer("MCTrackCont", "Track Container filled with MC information", AliHFEcuts::kNcutStepsMCTrack);
1370   fContainer->CreateContainer("recTrackContReco", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1371   fContainer->CreateContainer("recTrackContMC", "Track Container filled with MC information", AliHFEcuts::kNcutStepsRecTrack + fPID->GetNumberOfPIDdetectors());
1372
1373   fContainer->CreateContainer("hadronicBackground", "Container for Hadronic Background", 2);
1374   fContainer->CreateContainer("recTrackContDEReco", "Container for displaced electron analysis with Reco information", 1);
1375   fContainer->CreateContainer("recTrackContDEMC", "Container for displaced electron analysis with MC information", 1);
1376   fContainer->CreateContainer("recTrackContSecvtxReco", "Container for secondary vertexing analysis with Reco information", 1);
1377   fContainer->CreateContainer("recTrackContSecvtxMC", "Container for secondary vertexing analysis with MC information", 1);
1378
1379   if(HasMCData()){
1380     fContainer->CreateContainer("conversionElecs", "Container for weighted conversion electrons",3);
1381     fContainer->CreateContainer("mesonElecs", "Container for weighted electrons from non-HF meson decays",3);
1382     //fContainer->CreateContainer("charmElecs", "Container for weighted charm electrons",2);
1383   }
1384
1385   fContainer->CreateCorrelationMatrix("correlationstepafterPID","THnSparse with correlations");
1386   fContainer->CreateCorrelationMatrix("correlationstepafterDE","THnSparse with correlations");
1387   if(!fVarManager->IsVariableDefined("centrality")) {
1388     //printf("Create the two other correlation maps\n");
1389     fContainer->CreateCorrelationMatrix("correlationstepbeforePID","THnSparse with correlations");
1390     fContainer->CreateCorrelationMatrix("correlationstepafterTOF","THnSparse with correlations");
1391   }
1392
1393   // Define the step names
1394   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsMCTrack; istep++){
1395     fContainer->SetStepTitle("MCTrackCont", AliHFEcuts::MCCutName(istep), istep);
1396   }
1397   for(UInt_t istep = 0; istep < AliHFEcuts::kNcutStepsRecTrack; istep++){
1398     fContainer->SetStepTitle("recTrackContReco", AliHFEcuts::RecoCutName(istep), istep);
1399     fContainer->SetStepTitle("recTrackContMC", AliHFEcuts::RecoCutName(istep), istep);
1400   }
1401   for(UInt_t ipid = 0; ipid < fPID->GetNumberOfPIDdetectors(); ipid++){
1402     fContainer->SetStepTitle("recTrackContReco", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1403     fContainer->SetStepTitle("recTrackContMC", fPID->SortedDetectorName(ipid), AliHFEcuts::kNcutStepsRecTrack + ipid);
1404   }
1405 }
1406
1407 //____________________________________________________________
1408 void AliAnalysisTaskHFE::InitPIDperformanceQA(){
1409   // Add a histogram for Fake electrons
1410   const Int_t nDim=6;
1411   Int_t nBin[nDim] = {40, 8, 18, 2, 3, 3};
1412   //number of variables on the grid:pt,eta,phi,charge,
1413   const Double_t kPtbound[2] = {0.1, 20.};
1414   const Double_t kEtabound[2] = {-0.8, 0.8};
1415   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
1416   const Double_t kChargebound[2] = {-1.1, 1.1};
1417   const Double_t kAddInf1bound[2] = {0., 3.};
1418   const Double_t kAddInf2bound[2] = {0., 3.};
1419   Double_t minima[nDim] = {kPtbound[0], kEtabound[0], kPhibound[0], kChargebound[0], kAddInf1bound[0], kAddInf2bound[0]};
1420   Double_t maxima[nDim] = {kPtbound[1], kEtabound[1], kPhibound[1], kChargebound[1], kAddInf1bound[1], kAddInf2bound[1]};
1421
1422   fQACollection->CreateTHnSparse("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin, minima, maxima);
1423   fQACollection->CreateTHnSparse("SignalToBackgroundMC", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; charge; flavor (0 - no, 1 - charm, 2 - bottom); ITS Cluster (0 - no, 1 - first (and maybe second), 2 - second)", nDim, nBin, minima, maxima);
1424
1425   fQACollection->BinLogAxis("PIDperformance", 0);
1426   fQACollection->BinLogAxis("SignalToBackgroundMC", 0);
1427   fQACollection->Sumw2("PIDperformance");
1428   fQACollection->Sumw2("SignalToBackgroundMC");
1429 }
1430
1431 //____________________________________________________________
1432 void AliAnalysisTaskHFE::InitContaminationQA(){
1433   //
1434   // Add QA for Impact Parameter cut
1435   //
1436   const Double_t kPtbound[2] = {0.1, 20.};
1437   Int_t iBin[1];
1438   iBin[0] = 44; // bins in pt
1439   fQACollection->CreateTH1F("hadronsBeforeIPcut", "Hadrons before IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
1440   fQACollection->CreateTH1F("hadronsAfterIPcut", "Hadrons after IP cut", iBin[0], kPtbound[0], kPtbound[1], 1);
1441   fQACollection->CreateTH1F("hadronsBeforeIPcutMC", "Hadrons before IP cut: MC p_{t}", iBin[0], kPtbound[0], kPtbound[1], 1);
1442   fQACollection->CreateTH1F("hadronsAfterIPcutMC", "Hadrons after IP cut: MC p_{t} ", iBin[0],kPtbound[0], kPtbound[1], 1);
1443 }
1444
1445 //____________________________________________________________
1446 void AliAnalysisTaskHFE::PrintStatus() const {
1447   //
1448   // Print Analysis status
1449   //
1450   printf("\n\tAnalysis Settings\n\t========================================\n\n");
1451   printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
1452   printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
1453   printf("\tDisplaced electron analysis step: %s\n", GetPlugin(kDEstep) ? "YES" : "NO");
1454   printf("\tTagged Track Analysis: %s\n", GetPlugin(kTaggedTrackAnalysis) ? "YES" : "NO");
1455   printf("\n");
1456   printf("\tParticle Identification Detectors:\n");
1457   fPID->PrintStatus();
1458   printf("\n");
1459   printf("\tQA: \n");
1460   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
1461   printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
1462   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
1463   printf("\n");
1464 }
1465
1466 //____________________________________________________________
1467 Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
1468   //
1469   // Find the production vertex of the associated MC track
1470   //
1471   if(!fMCEvent) return kFALSE;
1472   const AliVParticle *mctrack = NULL;
1473   TString objectType = track->IsA()->GetName();
1474   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1475     // Reconstructed track
1476     mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
1477   } else {
1478     // MCParticle
1479     mctrack = track;
1480   }
1481
1482   if(!mctrack) return kFALSE;
1483
1484   Double_t xv = 0.0;
1485   Double_t yv = 0.0;
1486
1487   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1488     // case MCParticle
1489     const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(mctrack);
1490     if(mcpart){
1491       xv =  mcpart->Xv();
1492       yv =  mcpart->Yv();
1493     }
1494   } else {
1495     // case AODMCParticle
1496     const AliAODMCParticle *mcpart = dynamic_cast<const AliAODMCParticle *>(mctrack);
1497     if(mcpart){
1498       xv =  mcpart->Xv();
1499       yv =  mcpart->Yv();
1500     }
1501   }
1502
1503   //printf("xv %f, yv %f\n",xv,yv);
1504   fQACollection->Fill("radius", TMath::Abs(xv),TMath::Abs(yv));
1505
1506   return kTRUE;
1507
1508 }
1509 //__________________________________________
1510 void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
1511   //
1512   // Switch on Plugin
1513   // Available:
1514   //  - Primary vertex studies
1515   //  - Secondary vertex Studies
1516   //  - Post Processing
1517   //
1518   switch(plug){
1519     case kPriVtx: SETBIT(fPlugins, plug); break;
1520     case kSecVtx: SETBIT(fPlugins, plug); break;
1521     case kIsElecBackGround: SETBIT(fPlugins, plug); break;
1522     case kPostProcess: SETBIT(fPlugins, plug); break;
1523     case kDEstep: SETBIT(fPlugins, plug); break;
1524     case kTaggedTrackAnalysis: SETBIT(fPlugins, plug); break;
1525     default: AliError("Unknown Plugin");
1526   };
1527 }
1528 //__________________________________________
1529 Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track){
1530   //
1531   // Check single track cuts for a given cut step
1532   // Fill the particle container
1533   //
1534   const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack;
1535   if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE;
1536   if(fVarManager->IsSignalTrack()) {
1537     fVarManager->FillContainer(fContainer, "recTrackContReco", cutStep, kFALSE);
1538     fVarManager->FillContainer(fContainer, "recTrackContMC", cutStep, kTRUE);
1539   }
1540   return kTRUE;
1541 }
1542 //___________________________________________________
1543 Bool_t AliAnalysisTaskHFE::ReadCentrality() {
1544   //
1545   // Recover the centrality of the event from ESD or AOD
1546   //
1547  if(IsAODanalysis()){
1548
1549    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1550    if(!fAOD){
1551      AliError("AOD Event required for AOD Analysis");
1552      return kFALSE;
1553    }
1554
1555    if(IsPbPb()) {
1556      // Centrality
1557      AliCentrality *aodCentrality = fAOD->GetCentrality();
1558      Float_t fCentralityFtemp = aodCentrality->GetCentralityPercentile("V0M");
1559
1560      if( fCentralityFtemp >=  0. && fCentralityFtemp <  10.) fCentralityF =  0;
1561      else if ( fCentralityFtemp >=  10. && fCentralityFtemp <  20.) fCentralityF =  1;
1562      else if ( fCentralityFtemp >= 20. && fCentralityFtemp <  30.) fCentralityF = 2;
1563      else if ( fCentralityFtemp >= 30. && fCentralityFtemp <  40.) fCentralityF = 3;
1564      else if ( fCentralityFtemp >= 40. && fCentralityFtemp <  50.) fCentralityF = 4;
1565      else if ( fCentralityFtemp >= 50. && fCentralityFtemp <  60.) fCentralityF = 5;
1566      else if ( fCentralityFtemp >= 60. && fCentralityFtemp <  90.) fCentralityF = 6;
1567      else if ( fCentralityFtemp >= 90. && fCentralityFtemp <=  100.) fCentralityF = 7;
1568      //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <  95.) fCentralityF = 8;
1569      //else if ( fCentralityF_temp >= 95. && fCentralityF_temp <  90.) fCentralityF = 9;
1570      //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
1571      else return kFALSE;
1572
1573      // contributors
1574      fContributors = 0.5;
1575      Int_t contributorstemp = 0;
1576      const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
1577      if(vtxAOD) contributorstemp = vtxAOD->GetNContributors();
1578
1579      //printf("PbPb contributors_temp %d\n",contributors_temp);
1580
1581      if( contributorstemp <=  0) fContributors =  0.5;
1582      else fContributors = 1.5;
1583
1584
1585
1586    }
1587    else {
1588      fCentralityF = 0;
1589      Int_t centralityFtemp = 0;
1590      const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
1591      if(vtxAOD) centralityFtemp = vtxAOD->GetNContributors();
1592
1593      //printf("pp centralityF_temp %d\n",centralityF_temp);
1594
1595      if( centralityFtemp <=  0) fCentralityF =  0;
1596      else if ( centralityFtemp >  0 && centralityFtemp <  2) fCentralityF = 1;
1597      else if ( centralityFtemp >=  2 && centralityFtemp <  3) fCentralityF = 2;
1598      else if ( centralityFtemp >= 3 && centralityFtemp <  4) fCentralityF = 3;
1599      else if ( centralityFtemp >= 4 && centralityFtemp <  5) fCentralityF = 4;
1600      else if ( centralityFtemp >= 5 && centralityFtemp <  10) fCentralityF = 5;
1601      else if ( centralityFtemp >= 10 && centralityFtemp <  20) fCentralityF = 6;
1602      else if ( centralityFtemp >= 20 && centralityFtemp <  30) fCentralityF = 7;
1603      else if ( centralityFtemp >= 30 && centralityFtemp <  40) fCentralityF = 8;
1604      else if ( centralityFtemp >= 40 && centralityFtemp <  50) fCentralityF = 9;
1605      else if ( centralityFtemp >= 50) fCentralityF = 10;
1606
1607    }
1608
1609    return kTRUE;
1610
1611
1612  } else {
1613
1614    AliDebug(3, "Processing ESD Centrality");
1615    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
1616    if(!fESD){
1617      AliError("ESD Event required for ESD Analysis");
1618      return kFALSE;
1619    }
1620    const char* type = fESD->GetBeamType();
1621
1622    if (strstr(type,"Pb-Pb")) {
1623
1624    // Centrality
1625    AliCentrality *esdCentrality = fESD->GetCentrality();
1626    Float_t fCentralityFtemp = esdCentrality->GetCentralityPercentile("V0M");
1627    //printf("PbPb fCentralityF_temp %f\n",fCentralityF_temp);
1628
1629    if( fCentralityFtemp >=  0. && fCentralityFtemp <   10.) fCentralityF =  0;
1630    else if ( fCentralityFtemp >= 10. && fCentralityFtemp <  20.) fCentralityF =  1;
1631    else if ( fCentralityFtemp >= 20. && fCentralityFtemp <  30.) fCentralityF = 2;
1632    else if ( fCentralityFtemp >= 30. && fCentralityFtemp <  40.) fCentralityF = 3;
1633    else if ( fCentralityFtemp >= 40. && fCentralityFtemp <  50.) fCentralityF = 4;
1634    else if ( fCentralityFtemp >= 50. && fCentralityFtemp <  60.) fCentralityF = 5;
1635    else if ( fCentralityFtemp >= 60. && fCentralityFtemp <  80.) fCentralityF = 6;
1636    else if ( fCentralityFtemp >= 80. && fCentralityFtemp <  90.) fCentralityF = 7;
1637    else if ( fCentralityFtemp >= 90. && fCentralityFtemp <=  100.) fCentralityF = 8;
1638    //else if ( fCentralityF_temp >= 80. && fCentralityF_temp <  90.) fCentralityF = 9;
1639    //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
1640    else return kFALSE;
1641
1642    //   Float_t fCentralityF_temp10 = esdCentrality->GetCentralityClass10("V0M");
1643    //   printf("PbPb fCentralityF_temp %f %f %f \n",fCentralityF_temp, fCentralityF_temp10, fCentralityF);
1644
1645    // contributors
1646    fContributors = 0.5;
1647    Int_t contributorstemp = 0;
1648    const AliESDVertex *vtxESD = fESD->GetPrimaryVertexSPD();
1649    if(vtxESD) contributorstemp = vtxESD->GetNContributors();
1650
1651    //printf("PbPb contributors_temp %d\n",contributors_temp);
1652
1653    if( contributorstemp <=  0) fContributors =  0.5;
1654    else fContributors = 1.5;
1655
1656    return kTRUE;
1657
1658    }
1659
1660
1661    if (strstr(type,"p-p")) {
1662      fCentralityF = 0;
1663      Int_t centralityFtemp = 0;
1664      const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
1665      if(vtxESD && vtxESD->GetStatus()) centralityFtemp = vtxESD->GetNContributors();
1666
1667      //printf("pp centralityF_temp %d\n",centralityF_temp);
1668
1669      if( centralityFtemp <=  0) fCentralityF =  0;
1670      else if ( centralityFtemp >  0 && centralityFtemp <  2) fCentralityF = 1;
1671      else if ( centralityFtemp >=  2 && centralityFtemp <  3) fCentralityF = 2;
1672      else if ( centralityFtemp >= 3 && centralityFtemp <  4) fCentralityF = 3;
1673      else if ( centralityFtemp >= 4 && centralityFtemp <  5) fCentralityF = 4;
1674      else if ( centralityFtemp >= 5 && centralityFtemp <  10) fCentralityF = 5;
1675      else if ( centralityFtemp >= 10 && centralityFtemp <  20) fCentralityF = 6;
1676      else if ( centralityFtemp >= 20 && centralityFtemp <  30) fCentralityF = 7;
1677      else if ( centralityFtemp >= 30 && centralityFtemp <  40) fCentralityF = 8;
1678      else if ( centralityFtemp >= 40 && centralityFtemp <  50) fCentralityF = 9;
1679      else if ( centralityFtemp >= 50) fCentralityF = 10;
1680
1681      return kTRUE;
1682
1683    }
1684
1685    return kFALSE;
1686
1687   //printf("centrality %f\n",fCentralityF);
1688
1689  }
1690
1691  //printf("centrality %f\n",fCentralityF);
1692
1693 }
1694 //___________________________________________________
1695 void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
1696   //
1697   // Recover the centrality of the event from ESD or AOD
1698   //
1699  if(IsAODanalysis()){
1700
1701    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
1702    if(!fAOD){
1703      AliError("AOD Event required for AOD Analysis");
1704      return;
1705    }
1706    // PileUp
1707    if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
1708    // Z vertex
1709    if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
1710    // Event Cut
1711    fPassTheEventCut = kTRUE;
1712    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fAOD)) fPassTheEventCut = kFALSE;
1713
1714
1715  } else {
1716
1717    AliDebug(3, "Processing ESD Centrality");
1718    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
1719    if(!fESD){
1720      AliError("ESD Event required for ESD Analysis");
1721      return;
1722    }
1723    // PileUp
1724    fIdentifiedAsPileUp = kFALSE;
1725    if(fRemovePileUp && fESD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
1726    // Z vertex
1727    fIdentifiedAsOutInz = kFALSE;
1728    if (IsPbPb()) {
1729      //printf("PbPb\n");
1730      if(fESD->GetPrimaryVertex()){
1731        if(TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
1732      }
1733    }
1734    else {
1735      //printf("pp\n");
1736      if(fESD->GetPrimaryVertexTracks()){
1737        if(TMath::Abs(fESD->GetPrimaryVertexTracks()->GetZ()) > fCuts->GetVertexRange()) fIdentifiedAsOutInz = kTRUE;
1738      }
1739    }
1740    //Event Cut
1741    fPassTheEventCut = kTRUE;
1742    if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) fPassTheEventCut = kFALSE;
1743
1744
1745  }
1746
1747 }
1748