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 <TH1D.h>
  33 #include <TH1F.h>
  34 #include <TH1I.h>
  35 #include <TH2F.h>
  36 #include <TH3D.h>
  37 #include <TIterator.h>
  38 #include <TList.h>
  39 #include <TLegend.h>
  40 #include <TMath.h>
  41 #include <TObjArray.h>
  42 #include <TParticle.h>
  43 #include <TProfile.h>
  44 #include <TString.h>
  45 #include <TF1.h>
  46 #include <TTree.h>
  47
  48 #include "AliAODInputHandler.h"
  49 #include "AliAODMCParticle.h"
  50 #include "AliAODTrack.h"
  51 #include "AliCFContainer.h"
  52 #include "AliCFManager.h"
  53 #include "AliESDEvent.h"
  54 #include "AliESDInputHandler.h"
  55 #include "AliESDpid.h"
  56 #include "AliESDtrack.h"
  57 #include "AliLog.h"
  58 #include "AliAnalysisManager.h"
  59 #include "AliMCEvent.h"
  60 #include "AliMCEventHandler.h"
  61 #include "AliMCParticle.h"
  62 #include "AliPID.h"
  63 #include "AliStack.h"
  64 #include "AliVVertex.h"
  65
  66 #include "AliHFEpid.h"
  67 #include "AliHFEcollection.h"
  68 #include "AliHFEcuts.h"
  69 #include "AliHFEmcQA.h"
  70 #include "AliHFEpairs.h"
  71 #include "AliHFEpostAnalysis.h"
  72 #include "AliHFEsecVtxs.h"
  73 #include "AliHFEsecVtx.h"
  74 #include "AliHFEelecbackground.h"
  75 #include "AliHFEtools.h"
  76 #include "AliAnalysisTaskHFE.h"
  77
  78 ClassImp(AliAnalysisTaskHFE)
  79
  80 //____________________________________________________________
  81 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
  82   AliAnalysisTaskSE("PID efficiency Analysis")
  83   , fQAlevel(0)
  84   , fPIDdetectors("")
  85   , fPIDstrategy(0)
  86   , fPlugins(0)
  87   , fWeighting(kFALSE)
  88   , fWeightFactors(NULL)
  89   , fWeightFactorsFunction(NULL)
  90   , fCFM(NULL)
  91   , fHadronicBackground(NULL)
  92   , fCorrelation(NULL)
  93   , fPIDperformance(NULL)
  94   , fSignalToBackgroundMC(NULL)
  95   , fPID(NULL)
  96   , fCuts(NULL)
  97   , fSecVtx(NULL)
  98   , fElecBackGround(NULL)
  99   , fMCQA(NULL)
 100   , fNEvents(NULL)
 101   , fNElectronTracksEvent(NULL)
 102   , fQA(NULL)
 103   , fOutput(NULL)
 104   , fHistMCQA(NULL)
 105   , fHistSECVTX(NULL)
 106   , fHistELECBACKGROUND(NULL)
 107 //  , fQAcoll(0x0)
 108 {
 109   //
 110   // Dummy constructor
 111   //
 112 }
 113
 114 //____________________________________________________________
 115 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
 116   AliAnalysisTaskSE(name)
 117   , fQAlevel(0)
 118   , fPIDdetectors("")
 119   , fPIDstrategy(0)
 120   , fPlugins(0)
 121   , fWeighting(kFALSE)
 122   , fWeightFactors(NULL)
 123   , fWeightFactorsFunction(NULL)
 124   , fCFM(NULL)
 125   , fHadronicBackground(NULL)
 126   , fCorrelation(NULL)
 127   , fPIDperformance(NULL)
 128   , fSignalToBackgroundMC(NULL)
 129   , fPID(NULL)
 130   , fCuts(NULL)
 131   , fSecVtx(NULL)
 132   , fElecBackGround(NULL)
 133   , fMCQA(NULL)
 134   , fNEvents(NULL)
 135   , fNElectronTracksEvent(NULL)
 136   , fQA(NULL)
 137   , fOutput(NULL)
 138   , fHistMCQA(NULL)
 139   , fHistSECVTX(NULL)
 140   , fHistELECBACKGROUND(NULL)
 141 //  , fQAcoll(0x0)
 142 {
 143   //
 144   // Default constructor
 145   //
 146   DefineOutput(1, TH1I::Class());
 147   DefineOutput(2, TList::Class());
 148   DefineOutput(3, TList::Class());
 149 //  DefineOutput(4, TList::Class());
 150
 151   // Initialize cuts
 152 }
 153
 154 //____________________________________________________________
 155 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
 156   AliAnalysisTaskSE(ref)
 157   , fQAlevel(0)
 158   , fPIDdetectors()
 159   , fPIDstrategy(0)
 160   , fPlugins(0)
 161   , fWeighting(kFALSE)
 162   , fWeightFactors(NULL)
 163   , fWeightFactorsFunction(NULL)
 164   , fCFM(NULL)
 165   , fHadronicBackground(NULL)
 166   , fCorrelation(NULL)
 167   , fPIDperformance(NULL)
 168   , fSignalToBackgroundMC(NULL)
 169   , fPID(NULL)
 170   , fCuts(NULL)
 171   , fSecVtx(NULL)
 172   , fElecBackGround(NULL)
 173   , fMCQA(NULL)
 174   , fNEvents(NULL)
 175   , fNElectronTracksEvent(NULL)
 176   , fQA(NULL)
 177   , fOutput(NULL)
 178   , fHistMCQA(NULL)
 179   , fHistSECVTX(NULL)
 180   , fHistELECBACKGROUND(NULL)
 181 //  , fQAcoll(ref.fQAcoll)
 182 {
 183   //
 184   // Copy Constructor
 185   //
 186   ref.Copy(*this);
 187 }
 188
 189 //____________________________________________________________
 190 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
 191   //
 192   // Assignment operator
 193   //
 194   if(this == &ref)
 195     ref.Copy(*this);
 196   return *this;
 197 }
 198
 199 //____________________________________________________________
 200 void AliAnalysisTaskHFE::Copy(TObject &o) const {
 201   //
 202   // Copy into object o
 203   //
 204   AliAnalysisTaskHFE &target = dynamic_cast<AliAnalysisTaskHFE &>(o);
 205   target.fQAlevel = fQAlevel;
 206   target.fPIDdetectors = fPIDdetectors;
 207   target.fPIDstrategy = fPIDstrategy;
 208   target.fPlugins = fPlugins;
 209   target.fWeighting = fWeighting;
 210   target.fWeightFactors = fWeightFactors;
 211   target.fWeightFactorsFunction = fWeightFactorsFunction;
 212   target.fCFM = fCFM;
 213   target.fHadronicBackground = fHadronicBackground;
 214   target.fCorrelation = fCorrelation;
 215   target.fPIDperformance = fPIDperformance;
 216   target.fSignalToBackgroundMC = fSignalToBackgroundMC;
 217   target.fPID = fPID;
 218   target.fCuts = fCuts;
 219   target.fSecVtx = fSecVtx;
 220   target.fElecBackGround = fElecBackGround;
 221   target.fMCQA = fMCQA;
 222   target.fNEvents = fNEvents;
 223   target.fNElectronTracksEvent = fNElectronTracksEvent;
 224   target.fQA = fQA;
 225   target.fOutput = fOutput;
 226   target.fHistMCQA = fHistMCQA;
 227   target.fHistSECVTX = fHistSECVTX;
 228   target.fHistELECBACKGROUND = fHistELECBACKGROUND;
 229 }
 230
 231 //____________________________________________________________
 232 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
 233   //
 234   // Destructor
 235   //
 236   return;
 237   if(fPID) delete fPID;
 238   if(fQA){
 239     fQA->Clear();
 240     delete fQA;
 241   }
 242   if(fOutput){
 243     fOutput->Clear();
 244     delete fOutput;
 245   }
 246   if(fWeightFactors) delete fWeightFactors;
 247   if(fWeightFactorsFunction) delete fWeightFactorsFunction;
 248   if(fHistMCQA){
 249     fHistMCQA->Clear();
 250     delete fHistMCQA;
 251   }
 252   if(fHistSECVTX){
 253     fHistSECVTX->Clear();
 254     delete fHistSECVTX;
 255   }
 256   if(fHistELECBACKGROUND){
 257     fHistELECBACKGROUND->Clear();
 258     delete fHistELECBACKGROUND;
 259   }
 260   if(fSecVtx) delete fSecVtx;
 261   if(fElecBackGround) delete fElecBackGround;
 262   if(fMCQA) delete fMCQA;
 263   if(fNEvents) delete fNEvents;
 264   if(fCorrelation){
 265     fCorrelation->Clear();
 266     delete fCorrelation;
 267   }
 268   if(fPIDperformance) delete fPIDperformance;
 269   if(fSignalToBackgroundMC) delete fSignalToBackgroundMC;
 270 //  if(fQAcoll) delete fQAcoll;
 271
 272 }
 273
 274 //____________________________________________________________
 275 void AliAnalysisTaskHFE::UserCreateOutputObjects(){
 276   //
 277   // Creating output container and output objects
 278   // Here we also Initialize the correction framework container and
 279   // the objects for
 280   // - PID
 281   // - MC QA
 282   // - SecVtx
 283   // QA histograms are created if requested
 284   // Called once per worker
 285   //
 286   fPID = new AliHFEpid;
 287   AliDebug(3, "Creating Output Objects");
 288   // Automatic determination of the analysis mode
 289   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 290   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
 291     SetAODAnalysis();
 292   } else {
 293     SetESDAnalysis();
 294     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
 295       SetHasMCData();
 296   }
 297   printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
 298   printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
 299
 300   // example how to use the AliHFEcollection
 301   //fQAcoll = new AliHFEcollection("fQAcoll", "QA");
 302   //fQAcoll->CreateTH1F("fNevents", "Number of Events in the Analysis", 2, 0, 2);
 303   //fQAcoll->CreateProfile("fNtrdclusters", "Number of TRD clusters as function of momentum; p[GeV/c]", 20, 0, 20);
 304
 305   fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
 306   fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
 307   // First Step: TRD alone
 308   if(!fQA) fQA = new TList;
 309   fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
 310   fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
 311   fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
 312   fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
 313   fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
 314   fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
 315   fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
 316   fQA->AddAt(new TH1I("mccharge", "MC Charge", 200, -100, 100), 7);
 317   fQA->AddAt(new TH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0), 8);
 318   fQA->AddAt(new TH1F("secvtxept", "pT of tagged e", 500, 0, 50), 9); // mj: will move to another place soon
 319   fQA->AddAt(new TH2F("secvtxeTPCsig", "TPC signal for tagged e",125, 0, 25, 200, 0, 200 ), 10); // mj: will move to another place soon
 320
 321   if(!fOutput) fOutput = new TList;
 322   // Initialize correction Framework and Cuts
 323   fCFM = new AliCFManager;
 324   MakeParticleContainer();
 325   MakeEventContainer();
 326   // Temporary fix: Initialize particle cuts with NULL
 327   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
 328     fCFM->SetParticleCutsList(istep, NULL);
 329   if(!fCuts){
 330     AliWarning("Cuts not available. Default cuts will be used");
 331     fCuts = new AliHFEcuts;
 332     fCuts->CreateStandardCuts();
 333   }
 334   if(IsAODanalysis()) fCuts->SetAOD();
 335   fCuts->Initialize(fCFM);
 336   if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms());
 337
 338   // add output objects to the List
 339   fOutput->AddAt(fCFM->GetParticleContainer(), 0);
 340   fOutput->AddAt(fCFM->GetEventContainer(), 1);
 341   fOutput->AddAt(fCorrelation, 2);
 342   fOutput->AddAt(fPIDperformance, 3);
 343   fOutput->AddAt(fSignalToBackgroundMC, 4);
 344   fOutput->AddAt(fNElectronTracksEvent, 5);
 345   fOutput->AddAt(fHadronicBackground, 6);
 346
 347   // Initialize PID
 348   if(IsQAOn(kPIDqa)){
 349     AliInfo("PID QA switched on");
 350     //fPID->SetDebugLevel(2);
 351     fPID->SetQAOn();
 352     fQA->Add(fPID->GetQAhistograms());
 353   }
 354   fPID->SetHasMCData(HasMCData());
 355   if(!fPIDdetectors.Length() && ! fPIDstrategy) AddPIDdetector("TPC");
 356   if(fPIDstrategy)
 357     fPID->InitializePID(Form("Strategy%d", fPIDstrategy));
 358   else
 359     fPID->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed
 360
 361   // mcQA----------------------------------
 362   if (HasMCData() && IsQAOn(kMCqa)) {
 363     AliInfo("MC QA on");
 364     if(!fMCQA) fMCQA = new AliHFEmcQA;
 365     if(!fHistMCQA) fHistMCQA = new TList();
 366     fHistMCQA->SetName("MCqa");
 367     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
 368     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
 369     fMCQA->CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
 370     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm
 371     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
 372     fMCQA->CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
 373     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm
 374     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
 375     fMCQA->CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
 376     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm
 377     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
 378     fMCQA->CreateHistograms(AliHFEmcQA::kOthers,3,"mcqa_reccut_");       // create histograms for beauty
 379     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm
 380     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
 381     fMCQA->CreateHistograms(AliHFEmcQA::kOthers,4,"mcqa_recpidcut_");    // create histograms for beauty
 382     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,5,"mcqa_secvtxcut_");     // create histograms for charm
 383     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,5,"mcqa_secvtxcut_");    // create histograms for beauty
 384     fMCQA->CreateHistograms(AliHFEmcQA::kOthers,5,"mcqa_secvtxcut_");    // create histograms for beauty
 385     TIter next(gDirectory->GetList());
 386     TObject *obj;
 387     int counter = 0;
 388     TString objname;
 389     while ((obj = next.Next())) {
 390       objname = obj->GetName();
 391       TObjArray *toks = objname.Tokenize("_");
 392       if (toks->GetEntriesFast()){
 393         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
 394         if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++);
 395       }
 396     }
 397     fQA->Add(fHistMCQA);
 398   }
 399
 400   // secvtx----------------------------------
 401   if (GetPlugin(kSecVtx)) {
 402     AliInfo("Secondary Vertex Analysis on");
 403     fSecVtx = new AliHFEsecVtx;
 404     fSecVtx->SetHasMCData(HasMCData());
 405
 406     if(!fHistSECVTX) fHistSECVTX = new TList();
 407     fSecVtx->CreateHistograms(fHistSECVTX);
 408     fOutput->Add(fHistSECVTX);
 409   }
 410
 411   // background----------------------------------
 412   if (GetPlugin(kIsElecBackGround)) {
 413     AliInfo("Electron BackGround Analysis on");
 414     if(!fElecBackGround){
 415       AliWarning("ElecBackGround not available. Default elecbackground will be used");
 416       fElecBackGround = new AliHFEelecbackground;
 417     }
 418     fElecBackGround->SetHasMCData(HasMCData());
 419
 420     if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
 421     fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
 422     fOutput->Add(fHistELECBACKGROUND);
 423   }
 424 }
 425
 426 //____________________________________________________________
 427 void AliAnalysisTaskHFE::UserExec(Option_t *){
 428   //
 429   // Run the analysis
 430   //
 431   AliDebug(3, "Starting Single Event Analysis");
 432   if(!fInputEvent){
 433     AliError("Reconstructed Event not available");
 434     return;
 435   }
 436   if(HasMCData()){
 437     AliDebug(4, Form("MC Event: %p", fMCEvent));
 438     if(!fMCEvent){
 439       AliError("No MC Event, but MC Data required");
 440       return;
 441     }
 442   }
 443   if(!fCuts){
 444     AliError("HFE cuts not available");
 445     return;
 446   }
 447
 448   if(IsESDanalysis() && HasMCData()){
 449     // Protect against missing MC trees
 450     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 451     if(!mcH->InitOk()) return;
 452     if(!mcH->TreeK()) return;
 453     if(!mcH->TreeTR()) return;
 454   }
 455
 456   // Protect agains missing
 457   if(HasMCData()) ProcessMC();  // Run the MC loop + MC QA in case MC Data are available
 458
 459   if(IsAODanalysis()) ProcessAOD();
 460   else{
 461     AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
 462     AliESDpid *workingPID = inH->GetESDpid();
 463     if(workingPID){
 464       AliDebug(1, "Using ESD PID from the input handler");
 465       fPID->SetESDpid(workingPID);
 466     } else {
 467       AliDebug(1, "Using default ESD PID");
 468       fPID->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData()));
 469     }
 470     ProcessESD();
 471   }
 472   // Done!!!
 473   PostData(1, fNEvents);
 474   PostData(2, fOutput);
 475   PostData(3, fQA);
 476 //  PostData(4, fQAcoll->GetList());
 477 }
 478
 479 //____________________________________________________________
 480 void AliAnalysisTaskHFE::Terminate(Option_t *){
 481   //
 482   // Terminate not implemented at the moment
 483   //
 484   if(GetPlugin(kPostProcess)){
 485     fOutput = dynamic_cast<TList *>(GetOutputData(2));
 486     if(!fOutput){
 487       AliError("Results not available");
 488       return;
 489     }
 490     AliHFEpostAnalysis postanalysis;
 491     postanalysis.SetResults(fOutput);
 492     if(HasMCData())postanalysis.DrawMCSignal2Background();
 493     postanalysis.DrawEfficiency();
 494     postanalysis.DrawPIDperformance();
 495     postanalysis.DrawCutEfficiency();
 496
 497     if (GetPlugin(kIsElecBackGround)) {
 498       AliHFEelecbackground elecBackGround;
 499       TList *oe = 0x0;
 500       if(!(oe = (TList*)dynamic_cast<TList *>(fOutput->FindObject("HFEelecbackground")))){
 501         return;
 502       }
 503       elecBackGround.Load(oe);
 504       elecBackGround.Plot();
 505       elecBackGround.PostProcess();
 506     }
 507   }
 508 }
 509 //_______________________________________________________________
 510 Bool_t AliAnalysisTaskHFE::IsEventInBinZero() {
 511   //
 512   //
 513   //
 514
 515   //printf("test in IsEventInBinZero\n");
 516   if(!fInputEvent){
 517     AliError("Reconstructed Event not available");
 518     return kFALSE;
 519   }
 520
 521   // check vertex
 522   const AliVVertex *vertex = fInputEvent->GetPrimaryVertex();
 523   if(!vertex) return kTRUE;
 524   //if(vertex) return kTRUE;
 525
 526   // check tracks
 527   if(fInputEvent->GetNumberOfTracks()<=0) return kTRUE;
 528   //if(fInputEvent->GetNumberOfTracks()>0) return kTRUE;
 529
 530
 531   return kFALSE;
 532
 533 }
 534 //____________________________________________________________
 535 void AliAnalysisTaskHFE::ProcessMC(){
 536   //
 537   // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
 538   // In case MC QA is on also MC QA loop is done
 539   //
 540   AliDebug(3, "Processing MC Information");
 541   Double_t nContrib = 0;
 542   const AliVVertex *pVertex = fMCEvent->GetPrimaryVertex();
 543   if(pVertex) nContrib = pVertex->GetNContributors();
 544   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent)) return;
 545   fCFM->GetEventContainer()->Fill(&nContrib,AliHFEcuts::kEventStepGenerated);
 546   Int_t nElectrons = 0;
 547   if(IsESDanalysis()){
 548     if (HasMCData() && IsQAOn(kMCqa)) {
 549       AliDebug(2, "Running MC QA");
 550
 551       if(fMCEvent->Stack()){
 552         fMCQA->SetMCEvent(fMCEvent);
 553         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
 554         fMCQA->Init();
 555
 556         Int_t nMCTracks = fMCEvent->Stack()->GetNtrack();
 557
 558         // loop over all tracks for decayed electrons
 559         for (Int_t igen = 0; igen < nMCTracks; igen++){
 560           TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
 561           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
 562           fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
 563           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
 564           fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
 565           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0); // no accept cut
 566           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 567           fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 0); // no accept cut
 568           if (TMath::Abs(mcpart->Eta()) < 0.9) {
 569             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 570             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 571             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
 572           }
 573           if (TMath::Abs(AliHFEtools::GetRapidity(mcpart)) < 0.5) {
 574             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 575             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 576             fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kOthers, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
 577           }
 578         }
 579         fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
 580         fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
 581       }
 582
 583     } // end of MC QA loop
 584     // -----------------------------------------------------------------
 585     fCFM->SetMCEventInfo(fMCEvent);
 586     // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 587   } else {
 588     fCFM->SetMCEventInfo(fInputEvent);
 589   }
 590   // Run MC loop
 591   AliVParticle *mctrack = NULL;
 592   AliDebug(3, Form("Number of Tracks: %d", fMCEvent->GetNumberOfTracks()));
 593   for(Int_t imc = 0; imc <fMCEvent->GetNumberOfTracks(); imc++){
 594     if(!(mctrack = fMCEvent->GetTrack(imc))) continue;
 595     AliDebug(4, "Next Track");
 596     if(ProcessMCtrack(mctrack)) nElectrons++;
 597   }
 598
 599   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
 600   (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
 601 }
 602
 603 //____________________________________________________________
 604 void AliAnalysisTaskHFE::ProcessESD(){
 605   //
 606   // Run Analysis of reconstructed event in ESD Mode
 607   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
 608   //
 609   AliDebug(3, "Processing ESD Event");
 610   Double_t nContrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
 611   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fInputEvent)) return;
 612   fCFM->GetEventContainer()->Fill(&nContrib, AliHFEcuts::kEventStepReconstructed);
 613   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
 614   if(!fESD){
 615     AliError("ESD Event required for ESD Analysis")
 616     return;
 617   }
 618   if (GetPlugin(kIsElecBackGround)) {
 619     fElecBackGround->SetEvent(fESD);
 620   }
 621   if (GetPlugin(kSecVtx)) {
 622     fSecVtx->SetEvent(fESD);
 623     fSecVtx->GetPrimaryCondition();
 624   }
 625
 626   if(HasMCData()){
 627     if (GetPlugin(kSecVtx)) {
 628       fSecVtx->SetMCEvent(fMCEvent);
 629     }
 630     if (GetPlugin(kIsElecBackGround)) {
 631       fElecBackGround->SetMCEvent(fMCEvent);
 632     }
 633   }
 634
 635
 636   Double_t container[10];
 637   memset(container, 0, sizeof(Double_t) * 10);
 638   // container for the output THnSparse
 639   Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
 640   Int_t nElectronCandidates = 0;
 641   AliESDtrack *track = NULL, *htrack = NULL;
 642   AliMCParticle *mctrack = NULL;
 643   TParticle* mctrack4QA = NULL;
 644   Int_t pid = 0;
 645   // For double counted tracks
 646   LabelContainer cont(fESD->GetNumberOfTracks());
 647   Bool_t alreadyseen = kFALSE;
 648
 649   Bool_t signal = kTRUE;
 650
 651   fCFM->SetRecEventInfo(fESD);
 652   // Electron background analysis
 653   if (GetPlugin(kIsElecBackGround)) {
 654
 655     AliDebug(2, "Running BackGround Analysis");
 656
 657     fElecBackGround->Reset();
 658
 659   } // end of electron background analysis
 660   //
 661   // Loop ESD
 662   //
 663   AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
 664   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
 665
 666     track = fESD->GetTrack(itrack);
 667
 668     AliDebug(3, Form("Doing track %d, %p", itrack, track));
 669
 670     container[0] = track->Pt();
 671     container[1] = track->Eta();
 672     container[2] = track->Phi();
 673     container[3] = track->Charge();
 674     container[4] = 0;
 675
 676     dataE[0] = track->Pt();
 677     dataE[1] = track->Eta();
 678     dataE[2] = track->Phi();
 679     dataE[3] = track->Charge();
 680     dataE[4] = -1;
 681     dataE[5] = -1;
 682
 683     signal = kTRUE;
 684     Double_t weight = 1.0;
 685
 686     // Fill step without any cut
 687
 688     if(HasMCData()){
 689       container[4] = container[9] = kOther;
 690       // Check if it is electrons near the vertex
 691       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
 692       mctrack4QA = mctrack->Particle();
 693
 694       container[5] = mctrack->Pt();
 695       container[6] = mctrack->Eta();
 696       container[7] = mctrack->Phi();
 697       container[8] = mctrack->Charge()/3.;
 698
 699       if(fWeighting) weight = FindWeight(container[5],container[6],container[7]);
 700
 701       if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
 702       else AliDebug(3, "Signal Electron");
 703
 704       Int_t signalTrack = 0;
 705       if((signalTrack = IsSignalElectron(track))){
 706         AliDebug(3, Form("Signal: Index = %d\n", signalTrack));
 707         switch(signalTrack){
 708           case 1: container[4] = container[9] = kSignalCharm; break;
 709           case 2: container[4] = container[9] = kSignalBeauty; break;
 710           default: container[4] = container[9] = kOther; break;
 711         };
 712       } else if(IsGammaElectron(track)) container[4] = container[9] = kGammaConv;
 713       AliDebug(3, Form("Signal Decision(%f/%f)", container[4], container[9]));
 714     }
 715     AliDebug(3, Form("Weight? %f", weight));
 716     if(signal) {
 717       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));
 718       cont.Append(TMath::Abs(track->GetLabel()));
 719
 720       fCFM->GetParticleContainer()->Fill(&container[5], AliHFEcuts::kStepRecNoCut,weight);
 721       fCFM->GetParticleContainer()->Fill(&container[0], AliHFEcuts::kStepRecNoCut + 2*AliHFEcuts::kNcutStepsESDtrack,weight);
 722       if(alreadyseen) {
 723         fCFM->GetParticleContainer()->Fill(&container[5], AliHFEcuts::kStepRecNoCut + AliHFEcuts::kNcutStepsESDtrack,weight);
 724       }
 725     }
 726
 727     // RecKine: ITSTPC cuts
 728     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track, container, signal, alreadyseen, weight)) continue;
 729
 730     // Check TRD criterions (outside the correction framework)
 731     if(track->GetTRDncls()){
 732       (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
 733       (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());    // Check the acceptance without tight cuts
 734       (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
 735       (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
 736       //fQAcoll->Fill("fNtrdclusters", container[0], track->GetTRDncls());
 737     }
 738
 739
 740     // RecPrim
 741     if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track, container, signal, alreadyseen,weight)) continue;
 742
 743     // HFEcuts: ITS layers cuts
 744     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track, container, signal, alreadyseen,weight)) continue;
 745
 746     // HFEcuts: Nb of tracklets TRD0
 747     if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track, container, signal, alreadyseen,weight)) continue;
 748     if(signal) {
 749       // dimensions 3&4&5 : pt,eta,phi (MC)
 750       ((THnSparseF *)fCorrelation->At(0))->Fill(container);
 751     }
 752
 753     if(HasMCData() && IsQAOn(kMCqa)) {
 754       // mc qa for after the reconstruction cuts
 755       AliDebug(2, "Running MC QA");
 756       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 3);  // charm
 757       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 3); // beauty
 758       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kOthers,  AliHFEmcQA::kElectronPDG, 3); // beauty
 759     }
 760
 761     if(HasMCData()){
 762       FillProductionVertex(track);
 763     }
 764
 765     // track accepted, do PID
 766     AliHFEpidObject hfetrack;
 767     hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis;
 768     hfetrack.fRecTrack = track;
 769     if(HasMCData()) hfetrack.fMCtrack = mctrack;
 770     if(!fPID->IsSelected(&hfetrack)) continue;
 771     nElectronCandidates++;
 772
 773     // Fill Histogram for Hadronic Background
 774     if(HasMCData()){
 775       if(mctrack && (TMath::Abs(mctrack->Particle()->GetPdgCode()) != 11))
 776         fHadronicBackground->Fill(container, 0);
 777     }
 778
 779     if (HasMCData() && IsQAOn(kMCqa)) {
 780       // mc qa for after the reconstruction and pid cuts
 781       AliDebug(2, "Running MC QA");
 782       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 4);  // charm
 783       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 4); // beauty
 784       fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kOthers,  AliHFEmcQA::kElectronPDG, 4); // beauty
 785     }
 786
 787     // Fill Containers
 788     if(signal) {
 789       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack,weight);
 790       fCFM->GetParticleContainer()->Fill(&container[5], AliHFEcuts::kStepPID,weight);
 791       if(alreadyseen) {
 792         fCFM->GetParticleContainer()->Fill(&container[5], (AliHFEcuts::kStepPID + (AliHFEcuts::kNcutStepsESDtrack)),weight);
 793       }
 794       // dimensions 3&4&5 : pt,eta,phi (MC)
 795       ((THnSparseF *)fCorrelation->At(1))->Fill(container);
 796     }
 797
 798     if(GetPlugin(kSecVtx)) {
 799       AliDebug(2, "Running Secondary Vertex Analysis");
 800       if(track->Pt()>2.0 && nContrib > 1){
 801         fSecVtx->InitHFEpairs();
 802         fSecVtx->InitHFEsecvtxs();
 803         for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
 804           htrack = fESD->GetTrack(jtrack);
 805           if ( itrack == jtrack ) continue; // since it is for tagging single electron, don't need additional condition
 806           if (htrack->Pt()<2.0) continue;
 807           if (!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, htrack)) continue;
 808           if (!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, htrack)) continue;
 809           fSecVtx->PairAnalysis(track, htrack, jtrack); // e-h pairing
 810         }
 811         for(int ip=0; ip<fSecVtx->HFEpairs()->GetEntriesFast(); ip++){
 812           //if(HasMCData()){
 813             AliHFEpairs *pair = (AliHFEpairs*) (fSecVtx->HFEpairs()->UncheckedAt(ip));
 814             //if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.))  // apply various cuts
 815             // apply various cuts
 816             if(pair->GetKFChi2()>5.) // only apply vertex chi2 cut for the moment
 817             //if((pair->GetKFChi2()>5.) || !(pair->GetSignedLxy()>0. && pair->GetSignedLxy()<2.))
 818               fSecVtx->HFEpairs()->RemoveAt(ip);
 819           //}
 820         }
 821         fSecVtx->HFEpairs()->Compress();
 822         if(fSecVtx->HFEpairs()->GetEntriesFast()) fSecVtx->RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
 823         for(int ip=0; ip<fSecVtx->HFEsecvtxs()->GetEntriesFast(); ip++){
 824           AliHFEsecVtxs *secvtx=0x0;
 825           secvtx = (AliHFEsecVtxs*) (fSecVtx->HFEsecvtxs()->UncheckedAt(ip));
 826           if(!(secvtx->GetInvmass()>2.0 && secvtx->GetInvmass()<5.2) || !(secvtx->GetSignedLxy2()>0.08 && secvtx->GetSignedLxy2()<1.5) || !(secvtx->GetKFIP2()>-0.1 && secvtx->GetKFIP2()<0.1))
 827             fSecVtx->HFEsecvtxs()->RemoveAt(ip);
 828           // here you apply cuts, then if it doesn't pass the cut, remove it from the fSecVtx->HFEsecvtxs()
 829         }
 830         if(fSecVtx->HFEsecvtxs()->GetEntriesFast()) {
 831           (dynamic_cast<TH1F *>(fQA->At(9)))->Fill(track->Pt());
 832           (dynamic_cast<TH2F *>(fQA->At(10)))->Fill(track->P(),track->GetTPCsignal());
 833           if (HasMCData() && IsQAOn(kMCqa)) {
 834             // mc qa for after the reconstruction and pid cuts
 835             AliDebug(2, "Running MC QA");
 836             fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 5);  // charm
 837             fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty,  AliHFEmcQA::kElectronPDG, 5); // beauty
 838             fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kOthers,  AliHFEmcQA::kElectronPDG, 5); // beauty
 839           }
 840         }
 841         fSecVtx->DeleteHFEpairs();
 842         fSecVtx->DeleteHFEsecvtxs();
 843       }
 844     }
 845
 846     if(HasMCData()){
 847       // Track selected: distinguish between true and fake
 848       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
 849       if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
 850         Int_t type = IsSignalElectron(track);
 851         AliDebug(1, Form("Type: %d\n", type));
 852         if(type){
 853                 dataE[5] = type; // beauty[1] or charm[2]
 854                 dataE[4] = 2;  // signal electron
 855         }
 856         else{
 857                 dataE[4] = 1; // not a signal electron
 858                 dataE[5] = 0;
 859         }
 860       }
 861       else {
 862         // Fill THnSparse with the information for Fake Electrons
 863         dataE[4] = 0;
 864         dataE[5] = 0;
 865       }
 866       // fill the performance THnSparse, if the mc origin could be defined
 867       if(dataE[4] > -1){
 868         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
 869         fPIDperformance->Fill(dataE);
 870       }
 871     }
 872     // Electron background analysis
 873     if (GetPlugin(kIsElecBackGround)) {
 874
 875       AliDebug(2, "Running BackGround Analysis");
 876
 877       for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
 878               htrack = fESD->GetTrack(jtrack);
 879         if ( itrack == jtrack ) continue;
 880         fElecBackGround->PairAnalysis(track, htrack);
 881       }
 882     } // end of electron background analysis
 883
 884   }
 885   fNEvents->Fill(1);
 886   //fQAcoll->Fill("fNevents", 1);
 887   fNElectronTracksEvent->Fill(nElectronCandidates);
 888 }
 889
 890 //____________________________________________________________
 891 void AliAnalysisTaskHFE::ProcessAOD(){
 892   //
 893   // Run Analysis in AOD Mode
 894   // Function is still in development
 895   //
 896   AliDebug(3, "Processing AOD Event");
 897   Double_t nContrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
 898   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fInputEvent)) return;
 899   fCFM->GetEventContainer()->Fill(&nContrib,AliHFEcuts::kEventStepReconstructed);
 900   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
 901   if(!fAOD){
 902     AliError("AOD Event required for AOD Analysis")
 903     return;
 904   }
 905
 906   AliAODTrack *track = NULL;
 907   AliAODMCParticle *mctrack = NULL;
 908   Double_t container[10]; memset(container, 0, sizeof(Double_t) * 10);
 909   Double_t dataE[6]; // [pT, eta, Phi, Charge, type, 'C' or 'B']
 910   Int_t nElectronCandidates = 0;
 911   Int_t pid;
 912   for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
 913     track = fAOD->GetTrack(itrack);
 914     if(!track) continue;
 915     if(track->GetFlags() != 1<<4) continue;  // Only process AOD tracks where the HFE is set
 916
 917     container[0] = track->Pt();
 918     container[1] = track->Eta();
 919     container[2] = track->Phi();
 920     container[3] = track->Charge();
 921
 922     dataE[0] = track->Pt();
 923     dataE[1] = track->Eta();
 924     dataE[2] = track->Phi();
 925     dataE[3] = track->Charge();
 926     dataE[4] = -1;
 927     dataE[5] = -1;
 928
 929     if(HasMCData()){
 930       Int_t signalTrack = 0;
 931       if((signalTrack = IsSignalElectron(track))){
 932         switch(signalTrack){
 933           case 1: container[4] = container[9] = kSignalCharm; break;
 934           case 2: container[4] = container[9] = kSignalBeauty; break;
 935         };
 936       } else if(IsGammaElectron(track))
 937         container[4] = container[9] = kGammaConv;
 938       else container[4] = container[9] = kOther;
 939
 940       Int_t label = TMath::Abs(track->GetLabel());
 941       if(label){
 942         mctrack = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(label));
 943         container[5] = mctrack->Pt();
 944         container[6] = mctrack->Eta();
 945         container[7] = mctrack->Phi();
 946         container[8] = mctrack->Charge();
 947       }
 948     }
 949     // track accepted, do PID
 950     AliHFEpidObject hfetrack;
 951     hfetrack.fAnalysisType = AliHFEpidObject::kAODanalysis;
 952     hfetrack.fRecTrack = track;
 953     if(HasMCData()) hfetrack.fMCtrack = mctrack;
 954     //if(!fPID->IsSelected(&hfetrack)) continue;    // we will do PID here as soon as possible
 955     // Particle identified - Fill CF Container
 956     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack);
 957     nElectronCandidates++;
 958     if(HasMCData()){
 959       // Track selected: distinguish between true and fake
 960       AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->GetPdgCode()));
 961       if((pid = TMath::Abs(mctrack->GetPdgCode())) == 11){
 962         Int_t type = IsSignalElectron(track);
 963         AliDebug(1, Form("Type: %d\n", type));
 964         if(type){
 965                 dataE[5] = type; // beauty[1] or charm[2]
 966                 dataE[4] = 2;  // signal electron
 967         }
 968         else{
 969                 dataE[4] = 1; // not a signal electron
 970                 dataE[5] = 0;
 971         }
 972       }
 973       else {
 974         // Fill THnSparse with the information for Fake Electrons
 975         dataE[4] = 0;
 976         dataE[5] = 0;
 977       }
 978       // fill the performance THnSparse, if the mc origin could be defined
 979       if(dataE[4] > -1){
 980         AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4],dataE[5]));
 981         fPIDperformance->Fill(dataE);
 982       }
 983     }
 984   }
 985   fNEvents->Fill(1);
 986   fNElectronTracksEvent->Fill(nElectronCandidates);
 987 }
 988
 989 //____________________________________________________________
 990 Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
 991   //
 992   // Filter the Monte Carlo Track
 993   // Additionally Fill a THnSparse for Signal To Background Studies
 994   // Works for AOD and MC analysis Type
 995   //
 996   Double_t container[5], signalContainer[6];
 997   Double_t vertex[3]; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
 998   if(IsESDanalysis()){
 999     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
1000     container[0] = mctrack->Pt();
1001     container[1] = mctrack->Eta();
1002     container[2] = mctrack->Phi();
1003     container[3] = mctrack->Charge()/3;
1004
1005     signalContainer[0] = mctrack->Pt();
1006     signalContainer[1] = mctrack->Eta();
1007     signalContainer[2] = mctrack->Phi();
1008     signalContainer[3] = mctrack->Charge()/3;
1009
1010     vertex[0] = mctrack->Particle()->Vx();
1011     vertex[1] = mctrack->Particle()->Vy();
1012   } else {
1013     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
1014     container[0] = aodmctrack->Pt();
1015     container[1] = aodmctrack->Eta();
1016     container[2] = aodmctrack->Phi();
1017     container[3] = aodmctrack->Charge()/3;
1018
1019     signalContainer[0] = aodmctrack->Pt();
1020     signalContainer[1] = aodmctrack->Eta();
1021     signalContainer[2] = aodmctrack->Phi();
1022     signalContainer[3] = aodmctrack->Charge()/3;
1023
1024     aodmctrack->XvYvZv(vertex);
1025   }
1026   Int_t signal = 0;
1027   if((signal = IsSignalElectron(track))){
1028     switch(signal){
1029       case 1: container[4] = kSignalCharm; break;
1030       case 2: container[4] = kSignalBeauty; break;
1031     };
1032   }else if(IsGammaElectron(track)) container[4] = kGammaConv;
1033   else container[4] = kOther;
1034
1035   // weight
1036   Double_t weight = 1.0;
1037   if(fWeighting) weight = FindWeight(container[0],container[1],container[2]);
1038
1039  if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
1040   TH1 *test = dynamic_cast<TH1I*>(fQA->FindObject("mccharge"));
1041   test->Fill(signalContainer[3]);
1042  fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated,weight);
1043   if((signalContainer[4] = static_cast<Double_t >(IsSignalElectron(track))) > 1e-3) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal,weight);
1044   signalContainer[5] = 0;
1045   // apply cut on the sqrt of the production vertex
1046   Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
1047   if(radVertex < 3.5){
1048     // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
1049     signalContainer[5] = 1;
1050   } else if (radVertex < 7.5){
1051     signalContainer[5] = 2;
1052   }
1053   fSignalToBackgroundMC->Fill(signalContainer);
1054   (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(container[2] - TMath::Pi());
1055   //if(IsESDanalysis()){
1056     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
1057     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance,weight);
1058   //}
1059   return kTRUE;
1060 }
1061
1062 //____________________________________________________________
1063 void AliAnalysisTaskHFE::MakeEventContainer(){
1064   //
1065   // Create the event container for the correction framework and link it
1066   //
1067   const Int_t kNvar = 1;  // number of variables on the grid: number of tracks per event
1068   const Double_t kNTrackBound[2] = {-0.5, 200.5};
1069   const Int_t kNBins = 201;
1070
1071   AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, &kNBins);
1072
1073   Double_t *trackBins = AliHFEtools::MakeLinearBinning(kNBins, kNTrackBound[0], kNTrackBound[1]);
1074   evCont->SetBinLimits(0,trackBins);
1075   delete[] trackBins;
1076
1077   fCFM->SetEventContainer(evCont);
1078 }
1079
1080 //____________________________________________________________
1081 void AliAnalysisTaskHFE::MakeParticleContainer(){
1082   //
1083   // Create the particle container for the correction framework manager and
1084   // link it
1085   //
1086   const Int_t kNvar   = 5;
1087   //number of variables on the grid:pt,eta, phi, charge
1088   const Double_t kPtbound[2] = {0.1, 10.};
1089   const Double_t kEtabound[2] = {-0.8, 0.8};
1090   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
1091
1092   //arrays for the number of bins in each dimension
1093   Int_t iBin[kNvar];
1094   iBin[0] = 40; // bins in pt
1095   iBin[1] =  8; // bins in eta
1096   iBin[2] = 18; // bins in phi
1097   iBin[3] =  2; // bins in charge
1098   iBin[4] =  4; // creation process of the electron
1099
1100   //arrays for lower bounds :
1101   Double_t* binEdges[kNvar];
1102   binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
1103   binEdges[1] = AliHFEtools::MakeLinearBinning(iBin[1], kEtabound[0], kEtabound[1]);
1104   binEdges[2] = AliHFEtools::MakeLinearBinning(iBin[2], kPhibound[0], kPhibound[1]);
1105   binEdges[3] = AliHFEtools::MakeLinearBinning(iBin[3], -1.1, 1.1); // Numeric precision
1106   binEdges[4] = AliHFEtools::MakeLinearBinning(iBin[4], 0, iBin[4]); // Numeric precision
1107   //for(Int_t ib = 0; ib <= iBin[4]; ib++) printf("%f\t", binEdges[4][ib]);
1108   //printf("\n");
1109
1110   //one "container" for MC
1111   AliCFContainer* container = new AliCFContainer("trackContainer", "Container for tracks", (AliHFEcuts::kNcutStepsTrack + 2*AliHFEcuts::kNcutStepsESDtrack), kNvar, iBin);
1112   fHadronicBackground = new AliCFContainer("hadronicBackground", "Container for hadronic Background", 1, kNvar, iBin);
1113
1114   //setting the bin limits
1115   for(Int_t ivar = 0; ivar < kNvar; ivar++){
1116     container -> SetBinLimits(ivar, binEdges[ivar]);
1117     fHadronicBackground -> SetBinLimits(ivar, binEdges[ivar]);
1118   }
1119   fCFM->SetParticleContainer(container);
1120
1121   //create correlation matrix for unfolding
1122   Int_t thnDim[2*kNvar];
1123   for (int k=0; k<kNvar; k++) {
1124     //first half  : reconstructed
1125     //second half : MC
1126     thnDim[k]      = iBin[k];
1127     thnDim[k+kNvar] = iBin[k];
1128   }
1129
1130   if(!fCorrelation) fCorrelation = new TList();
1131   fCorrelation->SetName("correlation");
1132
1133   THnSparseF *correlation0 = new THnSparseF("correlationstepbeforePID","THnSparse with correlations",2*kNvar,thnDim);
1134   THnSparseF *correlation1 = new THnSparseF("correlationstepafterPID","THnSparse with correlations",2*kNvar,thnDim);
1135   for (int k=0; k<kNvar; k++) {
1136     correlation0->SetBinEdges(k,binEdges[k]);
1137     correlation0->SetBinEdges(k+kNvar,binEdges[k]);
1138     correlation1->SetBinEdges(k,binEdges[k]);
1139     correlation1->SetBinEdges(k+kNvar,binEdges[k]);
1140   }
1141   correlation0->Sumw2();
1142   correlation1->Sumw2();
1143
1144   fCorrelation->AddAt(correlation0,0);
1145   fCorrelation->AddAt(correlation1,1);
1146
1147   // Add a histogram for Fake electrons
1148   const Int_t nDim=6;
1149   Int_t nBin[nDim] = {40, 8, 18, 2, 3, 3};
1150   Double_t* binEdges2[nDim];
1151
1152   //values for bin lower bounds
1153   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1154     binEdges2[ivar] = binEdges[ivar];
1155   binEdges2[4] = AliHFEtools::MakeLinearBinning(nBin[4], 0, nBin[4]);
1156   binEdges2[5] = AliHFEtools::MakeLinearBinning(nBin[5], 0, nBin[5]);
1157
1158   fPIDperformance = new THnSparseF("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);
1159   fPIDperformance->Sumw2();
1160   fSignalToBackgroundMC = new THnSparseF("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);
1161   fSignalToBackgroundMC->Sumw2();
1162   for(Int_t idim = 0; idim < nDim; idim++){
1163     fPIDperformance->SetBinEdges(idim, binEdges2[idim]);
1164     fSignalToBackgroundMC->SetBinEdges(idim, binEdges2[idim]);
1165   }
1166   for(Int_t ivar = 0; ivar < kNvar; ivar++)
1167     delete binEdges[ivar];
1168   for(Int_t ivar = kNvar; ivar < nDim; ivar++)
1169     delete binEdges2[ivar];
1170 }
1171
1172 //____________________________________________________________
1173 void AliAnalysisTaskHFE::AddPIDdetector(TString detector){
1174   //
1175   // Adding PID detector to the task
1176   //
1177   if(!fPIDdetectors.Length())
1178     fPIDdetectors = detector;
1179   else
1180     fPIDdetectors += ":" + detector;
1181 }
1182
1183 //____________________________________________________________
1184 void AliAnalysisTaskHFE::PrintStatus() const {
1185   //
1186   // Print Analysis status
1187   //
1188   printf("\n\tAnalysis Settings\n\t========================================\n\n");
1189   printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
1190   printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
1191   printf("\n");
1192   printf("\tParticle Identification Detectors:\n");
1193   TObjArray *detectors = fPIDdetectors.Tokenize(":");
1194   for(Int_t idet = 0; idet < detectors->GetEntries(); idet++)
1195     printf("\t\t%s\n", (dynamic_cast<TObjString *>(detectors->At(idet)))->String().Data());
1196   printf("\n");
1197   printf("\tQA: \n");
1198   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
1199   printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsInDebugMode()) ? "YES" : "NO");
1200   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
1201   printf("\n");
1202 }
1203
1204 //____________________________________________________________
1205 AliAnalysisTaskHFE::LabelContainer::LabelContainer(Int_t capacity):
1206   fContainer(NULL),
1207   fBegin(NULL),
1208   fEnd(NULL),
1209   fLast(NULL),
1210   fCurrent(NULL)
1211 {
1212   //
1213   // Default constructor
1214   //
1215   fContainer = new Int_t[capacity];
1216   fBegin = &fContainer[0];
1217   fEnd = &fContainer[capacity - 1];
1218   fLast = fCurrent = fBegin;
1219 }
1220
1221 //____________________________________________________________
1222 Bool_t AliAnalysisTaskHFE::LabelContainer::Append(Int_t label){
1223   //
1224   // Add Label to the container
1225   //
1226   if(fLast > fEnd) return kFALSE;
1227   *fLast++ = label;
1228   return kTRUE;
1229 }
1230
1231 //____________________________________________________________
1232 Bool_t AliAnalysisTaskHFE::LabelContainer::Find(Int_t label) const {
1233   //
1234   // Find track in the list of labels
1235   //
1236   for(Int_t *entry = fBegin; entry <= fLast; entry++)
1237     if(*entry == label) return kTRUE;
1238   return kFALSE;
1239 }
1240
1241 //____________________________________________________________
1242 Int_t AliAnalysisTaskHFE::LabelContainer::Next(){
1243   //
1244   // Mimic iterator
1245   //
1246   if(fCurrent > fLast) return -1;
1247   return *fCurrent++;
1248 }
1249
1250 //____________________________________________________________
1251 Int_t AliAnalysisTaskHFE::IsSignalElectron(const AliVParticle * const track) const{
1252   //
1253   // Checks whether the identified electron track is coming from heavy flavour
1254   // returns 0 in case of no signal, 1 in case of charm and 2 in case of Bottom
1255   //
1256   enum{
1257     kNoSignal = 0,
1258     kCharm = 1,
1259     kBeauty = 2
1260   };
1261
1262   if(!fMCEvent) return kNoSignal;
1263   const AliVParticle *motherParticle = NULL, *mctrack = NULL;
1264   TString objectType = track->IsA()->GetName();
1265   Int_t label = 0;
1266   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1267     // Reconstructed track
1268     if((label = TMath::Abs(track->GetLabel())) && label < fMCEvent->GetNumberOfTracks())
1269       mctrack = fMCEvent->GetTrack(label);
1270   } else {
1271     // MCParticle
1272     mctrack = track;
1273   }
1274
1275   if(!mctrack) return kNoSignal;
1276
1277   Int_t pid = 0;
1278   Int_t daughterPDG = 0, motherLabel = 0;
1279   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1280     // case MC Particle
1281     daughterPDG = TMath::Abs((dynamic_cast<const AliMCParticle *>(mctrack))->Particle()->GetPdgCode());
1282     motherLabel = (dynamic_cast<const AliMCParticle *>(mctrack))->Particle()->GetFirstMother();
1283     if(motherLabel >= 0 && motherLabel < fMCEvent->GetNumberOfTracks())
1284       motherParticle = fMCEvent->GetTrack(motherLabel);
1285     if(motherParticle)
1286       pid = TMath::Abs((dynamic_cast<const AliMCParticle *>(motherParticle))->Particle()->GetPdgCode());
1287   } else {
1288     // case AODMCParticle
1289     daughterPDG = TMath::Abs((dynamic_cast<const AliAODMCParticle *>(mctrack))->GetPdgCode());
1290     motherLabel = (dynamic_cast<const AliAODMCParticle *>(mctrack))->GetMother();
1291     if(motherLabel >= 0 && motherLabel < fMCEvent->GetNumberOfTracks())
1292       motherParticle = fMCEvent->GetTrack(motherLabel);
1293     if(motherParticle)
1294       pid = TMath::Abs((dynamic_cast<const AliAODMCParticle *>(motherParticle))->GetPdgCode());
1295   }
1296   AliDebug(5, Form("Daughter PDG code: %d", daughterPDG));
1297
1298   if(!pid) return kNoSignal;
1299
1300   // From here the two analysis modes go together
1301   AliDebug(5, Form("Mother PDG code: %d", pid));
1302
1303   // identify signal according to Pdg Code - barions higher ranked than mesons
1304   if(pid / 1000 == 4) return kCharm;            // charmed baryon, 4th position in pdg code == 4
1305   if(pid / 1000 == 5) return kBeauty;           // beauty baryon, 4th position in pdg code == 5
1306   if((pid % 1000) / 100 == 4) return kCharm;    // charmed meson, 3rd position in pdg code == 4
1307   if((pid % 1000) / 100 == 5) return kBeauty;   // beauty meson, 3rd position in pdg code == 5
1308   return kNoSignal;
1309 }
1310
1311 //__________________________________________
1312 Bool_t AliAnalysisTaskHFE::IsGammaElectron(const AliVParticle * const track) const {
1313   //
1314   // Check for MC if the electron is coming from Gamma
1315   //
1316   if(!fMCEvent) return kFALSE;
1317   const AliVParticle *motherParticle = NULL, *mctrack = NULL;
1318   TString objectType = track->IsA()->GetName();
1319   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1320     // Reconstructed track
1321     if(track->GetLabel())
1322       mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
1323   } else {
1324     // MCParticle
1325     mctrack = track;
1326   }
1327
1328   if(!mctrack) return kFALSE;
1329
1330   Int_t motherPDG = 0;
1331   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1332     // case MC Particle
1333     motherParticle = fMCEvent->GetTrack((dynamic_cast<const AliMCParticle *>(mctrack)->Particle()->GetFirstMother()));
1334     if(motherParticle)
1335       motherPDG = TMath::Abs((dynamic_cast<const AliMCParticle *>(motherParticle))->Particle()->GetPdgCode());
1336   } else {
1337     // case AODMCParticle
1338     motherParticle = fMCEvent->GetTrack((dynamic_cast<const AliAODMCParticle *>(mctrack))->GetMother());
1339     if(motherParticle)
1340       motherPDG = TMath::Abs((dynamic_cast<const AliAODMCParticle *>(motherParticle))->GetPdgCode());
1341   }
1342   if(motherPDG!=22) return kFALSE;
1343   else return kTRUE;
1344 }
1345 //____________________________________________________________
1346 Bool_t AliAnalysisTaskHFE::FillProductionVertex(const AliVParticle * const track) const{
1347   //
1348   // Find the production vertex of the associated MC track
1349   //
1350   if(!fMCEvent) return kFALSE;
1351   const AliVParticle *mctrack = NULL;
1352   TString objectType = track->IsA()->GetName();
1353   if(objectType.CompareTo("AliESDtrack") == 0 || objectType.CompareTo("AliAODTrack") == 0){
1354     // Reconstructed track
1355     mctrack = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
1356   } else {
1357     // MCParticle
1358     mctrack = track;
1359   }
1360
1361   if(!mctrack) return kFALSE;
1362
1363   Double_t xv = 0.0;
1364   Double_t yv = 0.0;
1365
1366   if(TString(mctrack->IsA()->GetName()).CompareTo("AliMCParticle") == 0){
1367     // case MCParticle
1368     xv =  (dynamic_cast<const AliMCParticle *>(mctrack)->Xv());
1369     yv =  (dynamic_cast<const AliMCParticle *>(mctrack)->Yv());
1370
1371   } else {
1372     // case AODMCParticle
1373     xv =  (dynamic_cast<const AliAODMCParticle *>(mctrack)->Xv());
1374     yv =  (dynamic_cast<const AliAODMCParticle *>(mctrack)->Yv());
1375   }
1376
1377   //printf("xv %f, yv %f\n",xv,yv);
1378
1379   TH2F *test = dynamic_cast<TH2F*>(fQA->FindObject("radius"));
1380   if(!test) return kFALSE;
1381   else {
1382     test->Fill(TMath::Abs(xv),TMath::Abs(yv));
1383   }
1384
1385   return kTRUE;
1386
1387 }
1388 //__________________________________________
1389 void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
1390   //
1391   // Switch on Plugin
1392   // Available:
1393   //  - Primary vertex studies
1394   //  - Secondary vertex Studies
1395   //  - Post Processing
1396   //
1397   switch(plug){
1398     case kPriVtx: SETBIT(fPlugins, plug); break;
1399     case kSecVtx: SETBIT(fPlugins, plug); break;
1400     case kIsElecBackGround: SETBIT(fPlugins, plug); break;
1401     case kPostProcess: SETBIT(fPlugins, plug); break;
1402     default: AliError("Unknown Plugin");
1403   };
1404 }
1405 //_______________________________________________
1406 void AliAnalysisTaskHFE::SetWeightFactors(TH3D * const weightFactors){
1407   //
1408   // Set the histos with the weights for the efficiency maps
1409   //
1410   fWeighting = kTRUE;
1411   fWeightFactors = weightFactors;
1412 }
1413 //_______________________________________________
1414 void AliAnalysisTaskHFE::SetWeightFactorsFunction(TF1 * const weightFactorsFunction){
1415   //
1416   // Set the histos with the weights for the efficiency maps
1417   //
1418   fWeighting = kTRUE;
1419   fWeightFactorsFunction = weightFactorsFunction;
1420   //printf("SetWeightFactors\n");
1421 }
1422 //_______________________________________________
1423 Double_t AliAnalysisTaskHFE::FindWeight(Double_t pt, Double_t eta, Double_t phi) const {
1424   //
1425   // Find the weight corresponding to pt eta and phi in the TH3D
1426   //
1427   Double_t weight = 1.0;
1428   if(fWeightFactors) {
1429
1430     TAxis *ptaxis = fWeightFactors->GetXaxis();
1431     TAxis *etaaxis = fWeightFactors->GetYaxis();
1432     TAxis *phiaxis = fWeightFactors->GetZaxis();
1433
1434     Int_t ptbin = ptaxis->FindBin(pt);
1435     Int_t etabin = etaaxis->FindBin(eta);
1436     Int_t phibin = phiaxis->FindBin(phi);
1437
1438
1439     weight = fWeightFactors->GetBinContent(ptbin,etabin,phibin);
1440   }
1441   else if(fWeightFactorsFunction) {
1442
1443     weight = fWeightFactorsFunction->Eval(pt,eta,phi);
1444     //printf("pt %f and weight %f\n",pt,weight);
1445
1446   }
1447
1448   //printf("pt %f, eta %f, phi %f, weight %f\n",pt,eta,phi,weight);
1449
1450   return weight;
1451
1452 }
1453 //__________________________________________
1454 Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track, Double_t *container, Bool_t signal, Bool_t alreadyseen,Double_t weight){
1455   //
1456   // Check single track cuts for a given cut step
1457   // Fill the particle container
1458   //
1459   if(!fCFM->CheckParticleCuts(cutStep, track)) return kFALSE;
1460   if(signal) {
1461     fCFM->GetParticleContainer()->Fill(container, cutStep + 2*AliHFEcuts::kNcutStepsESDtrack,weight);
1462     fCFM->GetParticleContainer()->Fill(&container[5], cutStep,weight);
1463     if(alreadyseen) {
1464       fCFM->GetParticleContainer()->Fill(&container[5], cutStep + AliHFEcuts::kNcutStepsESDtrack,weight);
1465     }
1466   }
1467   return kTRUE;
1468 }
1469