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