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