Since it contains fixes of coding rule violations, all classes are involved. Further...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisTaskHFE.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15 /*
16  * The analysis task:
17  * Filling an AliCFContainer with the quantities pt, eta and phi
18  * for tracks which survivied the particle cuts (MC resp. ESD tracks)
19  * Track selection is done using the AliHFE package
20  * 
21  * Author:
22  *  Raphaelle Bailhache <R.Bailhache@gsi.de>
23  *  Markus Fasel <M.Fasel@gsi.de>
24  *  MinJung Kweon <minjung@physi.uni-heidelberg.de>
25  */
26 #include <TAxis.h>
27 #include <TCanvas.h>
28 #include <TChain.h>
29 #include <TDirectory.h>
30 #include <TH1F.h>
31 #include <TH1I.h>
32 #include <TH2F.h>
33 #include <TIterator.h>
34 #include <TList.h>
35 #include <TLegend.h>
36 #include <TMath.h>
37 #include <TObjArray.h>
38 #include <TParticle.h>
39 #include <TProfile.h>
40 #include <TString.h>
41 #include <TTree.h>
42
43 #include "AliCFContainer.h"
44 #include "AliCFManager.h"
45 #include "AliESDEvent.h"
46 #include "AliESDInputHandler.h"
47 #include "AliESDtrack.h"
48 #include "AliESDtrackCuts.h"
49 #include "AliLog.h"
50 #include "AliAnalysisManager.h"
51 #include "AliMCEvent.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCParticle.h"
54 #include "AliPID.h"
55 #include "AliStack.h"
56
57 #include "AliHFEpid.h"
58 #include "AliHFEcuts.h"
59 #include "AliHFEmcQA.h"
60 #include "AliHFEsecVtx.h"
61 #include "AliAnalysisTaskHFE.h"
62
63 //____________________________________________________________
64 AliAnalysisTaskHFE::AliAnalysisTaskHFE():
65   AliAnalysisTask("PID efficiency Analysis", "")
66   , fQAlevel(0)
67   , fPIDdetectors("")
68   , fESD(0x0)
69   , fMC(0x0)
70   , fCFM(0x0)
71   , fCorrelation(0x0)
72   , fFakeElectrons(0x0)
73   , fPID(0x0)
74   , fCuts(0x0)
75   , fSecVtx(0x0)
76   , fMCQA(0x0)
77   , fNEvents(0x0)
78   , fNElectronTracksEvent(0x0)
79   , fQA(0x0)
80   , fOutput(0x0)
81   , fHistMCQA(0x0)
82   , fHistSECVTX(0x0)
83 {
84   //
85   // Default constructor
86   //
87   DefineInput(0, TChain::Class());
88   DefineOutput(0, TH1I::Class());
89   DefineOutput(1, TList::Class());
90   DefineOutput(2, TList::Class());
91
92   // Initialize cuts
93   fCuts = new AliHFEcuts;
94   fPID = new AliHFEpid;
95 }
96
97 //____________________________________________________________
98 AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
99   AliAnalysisTask(ref)
100   , fQAlevel(ref.fQAlevel)
101   , fPIDdetectors(ref.fPIDdetectors)
102   , fESD(ref.fESD)
103   , fMC(ref.fMC)
104   , fCFM(ref.fCFM)
105   , fCorrelation(ref.fCorrelation)
106   , fFakeElectrons(ref.fFakeElectrons)
107   , fPID(ref.fPID)
108   , fCuts(ref.fCuts)
109   , fSecVtx(ref.fSecVtx)
110   , fMCQA(ref.fMCQA)
111   , fNEvents(ref.fNEvents)
112   , fNElectronTracksEvent(ref.fNElectronTracksEvent)
113   , fQA(ref.fQA)
114   , fOutput(ref.fOutput)
115   , fHistMCQA(ref.fHistMCQA)
116   , fHistSECVTX(ref.fHistSECVTX)
117 {
118   //
119   // Copy Constructor
120   //
121 }
122
123 //____________________________________________________________
124 AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref){
125   //
126   // Assignment operator
127   //
128   if(this == &ref) return *this;
129   AliAnalysisTask::operator=(ref);
130   fQAlevel = ref.fQAlevel;
131   fPIDdetectors = ref.fPIDdetectors;
132   fESD = ref.fESD;
133   fMC = ref.fMC;
134   fCFM = ref.fCFM;
135   fCorrelation = ref.fCorrelation;
136   fFakeElectrons = ref.fFakeElectrons;
137   fPID = ref.fPID;
138   fCuts = ref.fCuts;
139   fSecVtx = ref.fSecVtx;
140   fMCQA = ref.fMCQA;
141   fNEvents = ref.fNEvents;
142   fNElectronTracksEvent = ref.fNElectronTracksEvent;
143   fQA = ref.fQA;
144   fOutput = ref.fOutput;
145   fHistMCQA = ref.fHistMCQA;
146   fHistSECVTX = ref.fHistSECVTX;
147   return *this;
148 }
149
150 //____________________________________________________________
151 AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
152   //
153   // Destructor
154   //
155   if(fESD) delete fESD;
156   if(fMC) delete fMC;
157   if(fPID) delete fPID;
158   if(fQA){
159     fQA->Clear();
160     delete fQA;
161   }
162   if(fOutput){ 
163     fOutput->Clear();
164     delete fOutput;
165   }
166   if(fHistMCQA){
167     fHistMCQA->Clear();
168     delete fHistMCQA;
169   }
170   if(fHistSECVTX){
171     fHistSECVTX->Clear();
172     delete fHistSECVTX;
173   }
174   if(fCuts) delete fCuts;
175   if(fSecVtx) delete fSecVtx;
176   if(fMCQA) delete fMCQA;
177   if(fNEvents) delete fNEvents;
178   if(fCorrelation) delete fCorrelation;
179   if(fFakeElectrons) delete fFakeElectrons;
180 }
181
182 //____________________________________________________________
183 void AliAnalysisTaskHFE::ConnectInputData(Option_t *){
184 /*  TTree *esdchain = dynamic_cast<TChain *>(GetInputData(0));
185   if(!esdchain){
186     AliError("ESD chain empty");
187     return;
188   } else {
189     esdchain->SetBranchStatus("Tracks", 1);
190   }
191 */
192   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
193   if(!esdH){      
194     AliError("No ESD input handler");
195     return;
196   } else {
197     fESD = esdH->GetEvent();
198   }
199   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
200   if(!mcH){       
201     AliError("No MC truth handler");
202     return;
203   } else {
204     fMC = mcH->MCEvent();
205   }
206 }
207
208 //____________________________________________________________
209 void AliAnalysisTaskHFE::CreateOutputObjects(){
210   fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
211   fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
212   // First Step: TRD alone
213   if(!fQA) fQA = new TList;
214   fQA->AddAt(new TProfile("conr", "Electron PID contamination", 20, 0, 20), 0);
215   fQA->AddAt(new TH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi()), 1);
216   fQA->AddAt(new TH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi()), 2);
217   fQA->AddAt(new TH1F("nElectron", "Number of electrons", 100, 0, 100), 3);
218   fQA->AddAt(new TProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20), 4);
219   fQA->AddAt(new TProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20), 5);
220   fQA->AddAt(new TH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20), 6);
221
222   if(!fOutput) fOutput = new TList;
223   // Initialize correction Framework and Cuts
224   fCFM = new AliCFManager;
225   MakeParticleContainer();
226   // Temporary fix: Initialize particle cuts with 0x0
227   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
228     fCFM->SetParticleCutsList(istep, 0x0);
229   if(IsQAOn(kCUTqa)){
230     AliInfo("QA on for Cuts");
231     fCuts->SetDebugMode();
232     fQA->Add(fCuts->GetQAhistograms());
233   }
234   fCuts->CreateStandardCuts();
235   fCuts->SetMinNTrackletsTRD(0);  // Minimal requirement to get a minimum biased electron sample (only TPC and ITS requirements allowed)
236   fCuts->Initialize(fCFM);
237   // add output objects to the List
238   fOutput->AddAt(fCFM->GetParticleContainer(), 0);
239   fOutput->AddAt(fCorrelation, 1);
240   fOutput->AddAt(fFakeElectrons, 2);
241   fOutput->AddAt(fNElectronTracksEvent, 3);
242
243   // Initialize PID
244   if(IsQAOn(kPIDqa)){
245     AliInfo("PID QA switched on");
246     fPID->SetQAOn();
247     fQA->Add(fPID->GetQAhistograms());
248   }
249   fPID->SetHasMCData(kTRUE);
250   if(!fPIDdetectors.Length()) AddPIDdetector("TPC");
251   fPID->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed 
252
253   // mcQA----------------------------------
254   if (IsQAOn(kMCqa)) {
255     AliInfo("MC QA on");
256     if(!fMCQA) fMCQA = new AliHFEmcQA;
257     if(!fHistMCQA) fHistMCQA = new TList();
258     fHistMCQA->SetName("MCqa");
259     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
260     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
261     fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
262     fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
263     TIter next(gDirectory->GetList());
264     TObject *obj;
265     int counter = 0;
266     TString objname;
267     while ((obj = next.Next())) {
268       objname = obj->GetName();
269       TObjArray *toks = objname.Tokenize("_");
270       if (toks->GetEntriesFast()){
271         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
272         if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++);
273       }
274     }
275     fQA->Add(fHistMCQA);
276   } 
277
278   // secvtx----------------------------------
279   if (IsSecVtxOn()) {
280     AliInfo("Secondary Vertex Analysis on");
281     fSecVtx = new AliHFEsecVtx;
282
283     if(!fHistSECVTX) fHistSECVTX = new TList();
284     fHistSECVTX->SetName("SecVtx");
285     fSecVtx->CreateHistograms("secvtx_");
286     TIter next(gDirectory->GetList());
287     TObject *obj;
288     int counter = 0;
289     TString objname;
290     while ((obj = next.Next())) {
291       objname = obj->GetName();
292       TObjArray *toks = objname.Tokenize("_");
293       if (toks->GetEntriesFast()){
294         TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
295         if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++);
296       }
297     }
298     fOutput->Add(fHistSECVTX);
299   } 
300 }
301
302 //____________________________________________________________
303 void AliAnalysisTaskHFE::Exec(Option_t *){
304   //
305   // Run the analysis
306   // 
307   if(!fESD){
308     AliError("No ESD Event");
309     return;
310   }
311   if(!fMC){
312     AliError("No MC Event");
313     return;
314   }
315   fCFM->SetEventInfo(fMC);
316   fPID->SetMCEvent(fMC);
317
318   //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
319   Double_t container[6];
320
321   // Loop over the Monte Carlo tracks to see whether we have overlooked any track
322   AliMCParticle *mctrack = 0x0;
323   Int_t nElectrons = 0;
324
325   if (IsSecVtxOn()) { 
326     fSecVtx->SetEvent(fESD);
327     fSecVtx->SetStack(fMC->Stack());
328   }
329
330   // run mc QA ------------------------------------------------
331   if (IsQAOn(kMCqa)) {
332     AliDebug(2, "Running MC QA");
333
334     fMCQA->SetStack(fMC->Stack());
335     fMCQA->Init();
336
337     Int_t nPrims = fMC->Stack()->GetNprimary();
338     Int_t nMCTracks = fMC->Stack()->GetNtrack();
339
340     // loop over primary particles for quark and heavy hadrons
341     for (Int_t igen = 0; igen < nPrims; igen++){
342       fMCQA->GetQuarkKine(igen, AliHFEmcQA::kCharm);
343       fMCQA->GetQuarkKine(igen, AliHFEmcQA::kBeauty);
344       fMCQA->GetHadronKine(igen, AliHFEmcQA::kCharm);
345       fMCQA->GetHadronKine(igen, AliHFEmcQA::kBeauty);
346     }
347     fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
348     fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
349
350     // loop over all tracks for decayed electrons
351     for (Int_t igen = 0; igen < nMCTracks; igen++){
352       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 0);
353       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0);
354       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kCharm,  AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
355       fMCQA->GetDecayedKine(igen, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1, kTRUE); // barrel
356     }
357
358   } // end of MC QA loop
359   // -----------------------------------------------------------------
360
361   //
362   // Loop MC
363   //
364   for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
365     mctrack = fMC->GetTrack(imc);
366
367     container[0] = mctrack->Pt();
368     container[1] = mctrack->Eta();
369     container[2] = mctrack->Phi();
370
371     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
372     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
373     (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
374     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
375     // find the label in the vector
376     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
377     nElectrons++;
378   }
379   (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
380
381   // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
382         
383
384   //
385   // Loop ESD
386   //
387
388   Int_t nbtrackrec = fESD->GetNumberOfTracks();
389   Int_t nElectronCandidates = 0;
390   AliESDtrack *track = 0x0, *htrack = 0x0;
391   Int_t pid = 0;
392   // For double counted tracks
393   Int_t *indexRecKineTPC = new Int_t[nbtrackrec];
394   Int_t *indexRecKineITS = new Int_t[nbtrackrec];
395   Int_t *indexRecPrim = new Int_t[nbtrackrec];
396   Int_t *indexHFEcutsITS = new Int_t[nbtrackrec];
397   Int_t *indexHFEcutsTPC = new Int_t[nbtrackrec];
398   Int_t *indexHFEcutsTRD = new Int_t[nbtrackrec];
399   Int_t *indexPid = new Int_t[nbtrackrec];
400
401   for(Int_t nb = 0; nb < nbtrackrec; nb++){
402     indexRecKineTPC[nb] = -1;
403     indexRecKineITS[nb] = -1;
404     indexRecPrim[nb]    = -1;
405     indexHFEcutsITS[nb] = -1;
406     indexHFEcutsTPC[nb] = -1;
407     indexHFEcutsTRD[nb] = -1;
408     indexPid[nb]        = -1;
409   } 
410
411   Bool_t alreadyseen = kFALSE;
412
413   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
414     
415     track = fESD->GetTrack(itrack);
416           
417     container[0] = track->Pt();
418     container[1] = track->Eta();
419     container[2] = track->Phi();
420
421     Bool_t signal = kTRUE;
422           
423     // RecKine: TPC cuts  
424     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineTPC, track)) continue;
425     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineTPC);
426     
427
428     // Check if it is signal electrons
429     if(!(mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel())))) continue;
430     
431     container[3] = mctrack->Pt();
432     container[4] = mctrack->Eta();
433     container[5] = mctrack->Phi();
434     
435     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
436     
437     if(signal) {
438       alreadyseen = kFALSE;
439       for(Int_t nb = 0; nb < itrack; nb++){
440         if(indexRecKineTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
441       }
442       indexRecKineTPC[itrack] = TMath::Abs(track->GetLabel());
443       
444       fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineTPC + AliHFEcuts::kNcutESDSteps + 1));
445       fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 2*(AliHFEcuts::kNcutESDSteps + 1)));
446       if(alreadyseen) {
447         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
448       }
449       else {
450         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
451       }
452     }
453
454     
455     // RecKine: ITS cuts
456     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITS, track)) continue;
457     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecKineITS);
458     if(signal) {
459       alreadyseen = kFALSE;
460       for(Int_t nb = 0; nb < itrack; nb++){
461         if(indexRecKineITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
462       }
463       indexRecKineITS[itrack] = TMath::Abs(track->GetLabel());
464       fCFM->GetParticleContainer()->Fill(container, (AliHFEcuts::kStepRecKineITS + AliHFEcuts::kNcutESDSteps + 1));
465       fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 2*(AliHFEcuts::kNcutESDSteps + 1)));
466       if(alreadyseen) {
467         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
468       }
469       else {
470         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepRecKineITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
471       }
472     }
473
474
475     // Check TRD criterions (outside the correction framework)
476     if(track->GetTRDncls()){
477       (dynamic_cast<TH1F *>(fQA->At(6)))->Fill(track->GetTRDchi2()/track->GetTRDncls());
478       (dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha());    // Check the acceptance without tight cuts
479       (dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
480       (dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
481     }
482
483     
484     // RecPrim
485     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
486     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim);
487     if(signal) {
488       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
489       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
490       alreadyseen = kFALSE;
491       for(Int_t nb = 0; nb < itrack; nb++){
492         if(indexRecPrim[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
493       }
494       indexRecPrim[itrack] = TMath::Abs(track->GetLabel());
495       if(alreadyseen) {
496         fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 4*(AliHFEcuts::kNcutESDSteps + 1));
497       }
498       else {
499         fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + 3*(AliHFEcuts::kNcutESDSteps + 1));
500       }
501     }
502
503
504     // HFEcuts: ITS layers cuts
505     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
506     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS);
507     if(signal) {
508       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
509       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
510       alreadyseen = kFALSE;
511       for(Int_t nb = 0; nb < itrack; nb++){
512         if(indexHFEcutsITS[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
513       }
514       indexHFEcutsITS[itrack] = TMath::Abs(track->GetLabel());
515       if(alreadyseen) {
516         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 4*(AliHFEcuts::kNcutESDSteps + 1)));
517       }
518       else {
519         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsITS + 3*(AliHFEcuts::kNcutESDSteps + 1)));
520       }
521     }
522
523     // HFEcuts: TPC ratio clusters
524     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTPC, track)) continue;
525     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC);
526     if(signal) {
527       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTPC + AliHFEcuts::kNcutESDSteps + 1);
528       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTPC + 2*(AliHFEcuts::kNcutESDSteps + 1));
529       alreadyseen = kFALSE;
530       for(Int_t nb = 0; nb < itrack; nb++){
531         if(indexHFEcutsTPC[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
532       }
533       indexHFEcutsTPC[itrack] = TMath::Abs(track->GetLabel());    
534       if(alreadyseen) {
535         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 4*(AliHFEcuts::kNcutESDSteps + 1)));
536       }
537       else {
538         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTPC + 3*(AliHFEcuts::kNcutESDSteps + 1)));
539       }
540     }
541     
542     // HFEcuts: Nb of tracklets TRD
543     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
544     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD);
545     if(signal) {
546       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 1);
547       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
548       alreadyseen = kFALSE;
549       for(Int_t nb = 0; nb < itrack; nb++){
550         if(indexHFEcutsTRD[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
551       }
552       indexHFEcutsTRD[itrack] = TMath::Abs(track->GetLabel());
553       if(alreadyseen) {
554         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 4*(AliHFEcuts::kNcutESDSteps + 1)));
555       }
556       else {
557         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 3*(AliHFEcuts::kNcutESDSteps + 1)));
558       }
559     }
560
561
562     // track accepted, do PID
563     if(!fPID->IsSelected(track)) continue;
564     nElectronCandidates++;
565
566
567     // Fill Containers
568     fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 1);
569     if(signal) {
570       fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutESDSteps + 2);
571       fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1)+1);
572       alreadyseen = kFALSE;
573       for(Int_t nb = 0; nb < itrack; nb++){
574         if(indexPid[nb] == TMath::Abs(track->GetLabel())) alreadyseen = kTRUE;
575       }
576       indexPid[itrack] = TMath::Abs(track->GetLabel());
577       if(alreadyseen) {
578         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)));
579       }
580       else {
581         fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + 1 + 3*(AliHFEcuts::kNcutESDSteps + 1)));
582       }
583       // dimensions 3&4&5 : pt,eta,phi (MC)
584       fCorrelation->Fill(container);
585     }
586
587
588     // Track selected: distinguish between true and fake
589     if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
590       // pair analysis [mj]
591       if (IsSecVtxOn()) {
592         AliDebug(2, "Running Secondary Vertex Analysis");
593         fSecVtx->InitAnaPair();
594         for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
595           htrack = fESD->GetTrack(jtrack);
596           if ( itrack == jtrack ) continue;
597           //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue;
598           if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue;
599           fSecVtx->AnaPair(track, htrack, jtrack);
600         }
601         // based on the partner of e info, you run secandary vertexing function
602         fSecVtx->RunSECVTX(track);
603       } // end of pair analysis
604     } else {
605       // Fill THnSparse with the information for Fake Electrons
606       switch(pid){
607       case 13:    container[3] = AliPID::kMuon;
608       case 211:   container[3] = AliPID::kPion;
609       case 321:   container[3] = AliPID::kKaon;
610       case 2212:  container[3] = AliPID::kProton;
611       }
612       fFakeElectrons->Fill(container);
613     }
614   }
615   fNEvents->Fill(1);
616   fNElectronTracksEvent->Fill(nElectronCandidates);
617
618   // release memory
619   delete[] indexRecKineTPC;
620   delete[] indexRecKineITS;
621   delete[] indexRecPrim;
622   delete[] indexHFEcutsITS;
623   delete[] indexHFEcutsTPC;
624   delete[] indexHFEcutsTRD;
625   delete[] indexPid;
626   
627   // Done!!!
628   PostData(0, fNEvents);
629   PostData(1, fOutput);
630   PostData(2, fQA);
631 }
632
633 //____________________________________________________________
634 void AliAnalysisTaskHFE::Terminate(Option_t *){
635   //
636   // Terminate not implemented at the moment
637   //
638 }
639
640 //____________________________________________________________
641 void AliAnalysisTaskHFE::MakeParticleContainer(){
642   //
643   // Create the particle container for the correction framework manager and 
644   // link it
645   //
646   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
647   const Double_t kPtmin = 0.1, kPtmax = 10.;
648   const Double_t kEtamin = -0.9, kEtamax = 0.9;
649   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
650
651   //arrays for the number of bins in each dimension
652   Int_t iBin[kNvar];
653   iBin[0] = 40; //bins in pt
654   iBin[1] =  8; //bins in eta 
655   iBin[2] = 18; // bins in phi
656
657   //arrays for lower bounds :
658   Double_t* binEdges[kNvar];
659   for(Int_t ivar = 0; ivar < kNvar; ivar++)
660     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
661
662   //values for bin lower bounds
663   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);  
664   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
665   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
666
667   //one "container" for MC
668   AliCFContainer* container = new AliCFContainer("container","container for tracks", (AliHFEcuts::kNcutSteps + 1 + 4*(AliHFEcuts::kNcutESDSteps + 1)), kNvar, iBin);
669
670   //setting the bin limits
671   for(Int_t ivar = 0; ivar < kNvar; ivar++)
672     container -> SetBinLimits(ivar, binEdges[ivar]);
673   fCFM->SetParticleContainer(container);
674
675   //create correlation matrix for unfolding
676   Int_t thnDim[2*kNvar];
677   for (int k=0; k<kNvar; k++) {
678     //first half  : reconstructed 
679     //second half : MC
680     thnDim[k]      = iBin[k];
681     thnDim[k+kNvar] = iBin[k];
682   }
683
684   fCorrelation = new THnSparseF("correlation","THnSparse with correlations",2*kNvar,thnDim);
685   for (int k=0; k<kNvar; k++) {
686     fCorrelation->SetBinEdges(k,binEdges[k]);
687     fCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
688   }
689   fCorrelation->Sumw2();
690
691   // Add a histogram for Fake electrons
692   thnDim[3] = AliPID::kSPECIES;
693   fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
694   for(Int_t idim = 0; idim < kNvar; idim++)
695     fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
696 }
697
698 void AliAnalysisTaskHFE::AddPIDdetector(Char_t *detector){
699   if(!fPIDdetectors.Length()) 
700     fPIDdetectors = detector;
701   else
702     fPIDdetectors += Form(":%s", detector);
703 }
704
705 //____________________________________________________________
706 void AliAnalysisTaskHFE::PrintStatus(){
707   printf("\n\tAnalysis Settings\n\t========================================\n\n");
708   printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO");
709   printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO");
710   printf("\n");
711   printf("\tParticle Identification Detectors:\n");
712   TObjArray *detectors = fPIDdetectors.Tokenize(":");
713   for(Int_t idet = 0; idet < detectors->GetEntries(); idet++)
714     printf("\t\t%s\n", (dynamic_cast<TObjString *>(detectors->At(idet)))->String().Data());
715   printf("\n");
716   printf("\tQA: \n");
717   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
718   printf("\t\tCUTS: %s\n", IsQAOn(kCUTqa) ? "YES" : "NO");
719   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
720   printf("\n");
721 }