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