Update of the HFE package
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisTaskDisplacedElectrons.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 // study displaced electrons from beauty and charm 
18 // with cut on impact parameters in various pT bins
19 // 
20 // 
21 // Authors:
22 //  Hongyan Yang <hongyan@physi.uni-heidelberg.de>
23 //  Carlo Bombonati <Carlo.Bombonati@cern.ch>
24 //
25
26 #include <TChain.h>
27 #include <TFile.h>
28 #include <TH1F.h>
29 #include <TH1I.h>
30 #include "AliLog.h"
31
32 #include <TList.h>
33 #include <TMath.h>
34 #include <TObjArray.h>
35 #include <TParticle.h>
36 #include <TString.h>
37 #include <TCanvas.h>
38
39 #include "AliCFManager.h"
40 #include "AliMCEvent.h"
41 #include "AliMCEventHandler.h"
42 #include "AliMCParticle.h"
43 #include "AliStack.h"
44
45 #include "AliESDEvent.h"
46 #include "AliESDInputHandler.h"
47 #include "AliESDtrack.h"
48 #include "AliESDpid.h"
49
50 #include "AliAnalysisManager.h"
51
52 #include "AliHFEpid.h"
53 #include "AliHFEcuts.h"
54 #include "AliHFEtools.h"
55 #include "AliHFEdisplacedElectrons.h"
56 #include "AliAnalysisTaskDisplacedElectrons.h"
57
58
59 //____________________________________________________________
60 AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons():
61   AliAnalysisTaskSE("Task for displaced electron study")
62   , fDeDebugLevel(0)
63   , fNminITSCluster(0)
64   , fNminPrimVtxContrib(0)
65   , fDePIDdetectors("")
66   , fDePIDstrategy(0)
67   , fDePlugins(0)
68   , fDeCuts(0x0)
69   , fDeDefaultPID(0x0)
70   , fDePID(0x0)
71   , fDeCFM(0x0)
72   , fDisplacedElectrons(0x0)
73   , fDeNEvents(0x0)
74   , fElectronsMcPt(0x0)
75   , fElectronsEsdPt(0x0)
76   , fElectronsDataPt(0x0)
77   , fDeCorrection(0x0)
78   , fDeQA(0x0)
79   , fHistDisplacedElectrons(0x0)
80 {
81   //
82   // Dummy constructor
83   //
84 }
85
86 //____________________________________________________________
87 AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons(const char * name):
88   AliAnalysisTaskSE(name)
89   , fDeDebugLevel(0)
90   , fNminITSCluster(0)
91   , fNminPrimVtxContrib(0)
92   , fDePIDdetectors("")
93   , fDePIDstrategy(0)
94   , fDePlugins(0)
95   , fDeCuts(0x0)
96   , fDeDefaultPID(0x0)
97   , fDePID(0x0)
98   , fDeCFM(0x0)
99   , fDisplacedElectrons(0x0)
100   , fDeNEvents(0x0)
101   , fElectronsMcPt(0x0)
102   , fElectronsEsdPt(0x0)
103   , fElectronsDataPt(0x0)
104   , fDeCorrection(0x0)
105   , fDeQA(0x0)
106   , fHistDisplacedElectrons(0x0)
107 {
108   //
109   // Default constructor
110   //
111   DefineInput(0, TChain::Class());
112   DefineOutput(1, TList::Class());
113   DefineOutput(2, TList::Class());
114   DefineOutput(3, TList::Class());
115
116   // Initialize pid
117   fDeDefaultPID = new AliESDpid;
118   fDePID = new AliHFEpid("DEPID");
119
120 }
121
122 //____________________________________________________________
123 AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons(const AliAnalysisTaskDisplacedElectrons &ref):
124   AliAnalysisTaskSE(ref)
125   , fDeDebugLevel(ref.fDeDebugLevel)
126   , fNminITSCluster(ref.fNminITSCluster)
127   , fNminPrimVtxContrib(ref.fNminPrimVtxContrib)
128   , fDePIDdetectors(ref.fDePIDdetectors)
129   , fDePIDstrategy(ref.fDePIDstrategy)
130   , fDePlugins(ref.fDePlugins)
131   , fDeCuts(ref.fDeCuts)
132   , fDeDefaultPID(ref.fDeDefaultPID)
133   , fDePID(ref.fDePID)
134   , fDeCFM(ref.fDeCFM)
135   , fDisplacedElectrons(ref.fDisplacedElectrons)
136   , fDeNEvents(ref.fDeNEvents)
137   , fElectronsMcPt(ref.fElectronsMcPt)
138   , fElectronsEsdPt(ref.fElectronsEsdPt)
139   , fElectronsDataPt(ref.fElectronsDataPt)
140   , fDeCorrection(ref.fDeCorrection)
141   , fDeQA(ref.fDeQA)
142   , fHistDisplacedElectrons(ref.fHistDisplacedElectrons)
143 {
144   //
145   // Copy Constructor
146   //
147 }
148
149 //____________________________________________________________
150 AliAnalysisTaskDisplacedElectrons &AliAnalysisTaskDisplacedElectrons::operator=(const AliAnalysisTaskDisplacedElectrons &ref){
151   //
152   // Assignment operator
153   //
154   if(this == &ref) return *this;
155   AliAnalysisTask::operator=(ref);
156   fDeDebugLevel = ref.fDeDebugLevel;
157   fNminITSCluster = ref.fNminITSCluster;
158   fNminPrimVtxContrib = ref.fNminPrimVtxContrib;
159   fDePIDdetectors = ref.fDePIDdetectors;
160   fDePIDstrategy = ref.fDePIDstrategy;
161   fDePlugins = ref.fDePlugins;
162   fDeDefaultPID = ref.fDeDefaultPID;
163   fDePID = ref.fDePID;
164   fDeCuts = ref.fDeCuts;
165   fDeCFM = ref.fDeCFM;
166   fDisplacedElectrons = ref.fDisplacedElectrons;
167   fDeNEvents = ref.fDeNEvents;
168   fElectronsMcPt = ref.fElectronsMcPt;
169   fElectronsEsdPt = ref.fElectronsEsdPt;
170   fElectronsDataPt = ref.fElectronsDataPt;
171   fDeCorrection = ref.fDeCorrection;
172   fDeQA = ref.fDeQA;
173   fHistDisplacedElectrons = ref.fHistDisplacedElectrons;
174
175   return *this;
176 }
177
178 //____________________________________________________________
179 AliAnalysisTaskDisplacedElectrons::~AliAnalysisTaskDisplacedElectrons(){
180   //
181   // Destructor
182   //
183
184   if(fDeDefaultPID) delete fDeDefaultPID;
185   if(fDePID) delete fDePID;
186   if(fDeCFM) delete fDeCFM;
187   if(fDisplacedElectrons) delete fDisplacedElectrons;  
188   if(fDeNEvents) delete fDeNEvents;
189   if(fElectronsMcPt) delete fElectronsMcPt;
190   if(fElectronsEsdPt) delete fElectronsEsdPt;
191   if(fElectronsDataPt) delete fElectronsDataPt;
192   if(fDeCorrection){
193     fDeCorrection->Clear();
194     delete fDeCorrection;
195   }
196   if(fDeQA){
197     fDeQA->Clear();
198     delete fDeQA;
199   }
200   if(fHistDisplacedElectrons){ 
201     fHistDisplacedElectrons->Clear();  
202     delete fHistDisplacedElectrons;  
203   }
204 }
205
206 //____________________________________________________________
207 void AliAnalysisTaskDisplacedElectrons::UserCreateOutputObjects(){
208   // create output objects
209   // fDeNEvents
210   // MC and Data containers
211
212
213   if(!fDeQA) fDeQA = new TList;
214   fDeQA->SetName("variousQAhistograms");
215   
216   fDeNEvents = new TH1I("nDeEvents", "Number of Events in the DE Analysis", 2, 0, 2); 
217   const Int_t nBins = 14;
218   const Float_t ptBins[nBins] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,4.0,5.0,7.0,9.0,12.0,16.0,20.0};
219   fElectronsMcPt = new TH1F("mcElectronPt", "MC: p_{T} distribution of identified electrons (mcpid);p_{T} (GeV/c);Counts;", nBins-1, ptBins); 
220   fElectronsEsdPt = new TH1F("esdElectronPt", "ESD: p_{T} distribution of identified electrons (hfepid);p_{T} (GeV/c);Counts;", nBins-1, ptBins); 
221   fElectronsDataPt = new TH1F("dataElectronPt", "DATA: p_{T} distribution of identified electrons (hfepid);p_{T} (GeV/c);Counts;", nBins-1, ptBins); 
222
223   fDeQA->AddAt(fDeNEvents,0);
224   if(HasMCData()){
225     fDeQA->AddAt(fElectronsMcPt, 1);
226     fDeQA->AddAt(fElectronsEsdPt, 2); 
227   }
228   else{
229     fDeQA->AddAt(fElectronsDataPt, 1);  
230   }
231   // Initialize correction Framework and Cuts
232   fDeCFM = new AliCFManager;
233   MakeEventContainer();
234   MakeParticleContainer();
235  
236   if(!fDeCorrection) fDeCorrection = new TList();
237   fDeCorrection->SetName("deCorrections");
238   fDeCorrection->AddAt(fDeCFM->GetEventContainer(), 0);
239   fDeCorrection->AddAt(fDeCFM->GetParticleContainer(), 1);
240   fDeCorrection->Print();
241
242   for(Int_t istep = 0; istep < fDeCFM->GetEventContainer()->GetNStep(); istep++)
243     fDeCFM->SetEventCutsList(istep, 0x0);
244   for(Int_t istep = 0; istep < fDeCFM->GetParticleContainer()->GetNStep(); istep++)
245     fDeCFM->SetParticleCutsList(istep, 0x0);
246
247   if(!fDeCuts){
248     AliWarning("Cuts not available. Default cuts will be used");
249     fDeCuts = new AliHFEcuts;
250     fDeCuts->CreateStandardCuts();    
251   }
252   
253   fDeCuts->Initialize(fDeCFM);
254   
255   fDePID->SetHasMCData(HasMCData());
256   fDePID->AddDetector("TPC", 0);
257   fDePID->AddDetector("TOF", 1);
258         fDePID->ConfigureTOF();
259   fDePID->ConfigureTPCdefaultCut();
260   fDePID->InitializePID();     // Only restrictions to TPC allowed   
261
262
263   // displaced electron study----------------------------------
264   if(GetPlugin(kDisplacedElectrons)){
265     
266     fDisplacedElectrons = new AliHFEdisplacedElectrons;
267     fDisplacedElectrons->SetDebugLevel(fDeDebugLevel);
268     fDisplacedElectrons->SetHasMCData(HasMCData());
269     fDisplacedElectrons->SetMinPrimVtxContrib(fNminPrimVtxContrib);
270     fDisplacedElectrons->SetNitsCluster(fNminITSCluster);
271   
272     if(!fHistDisplacedElectrons) fHistDisplacedElectrons = new TList();
273     fDisplacedElectrons->CreateOutputs(fHistDisplacedElectrons);
274   }
275   
276 }
277
278
279
280 //____________________________________________________________
281 void AliAnalysisTaskDisplacedElectrons::UserExec(Option_t *){
282   //
283   // Run the analysis
284   // 
285
286   if(fDeDebugLevel>=10) AliInfo("analyse single event");
287
288   if(!fInputEvent){
289     AliError("Reconstructed Event not available");
290     return;
291   }
292  
293   //  
294   AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
295   
296   if(!inH){
297     AliError("AliESDInputHandler dynamic cast failed!");
298     return;
299   }
300   
301   // from now on, only ESD are analyzed
302   // using HFE pid, using HFE cuts
303   // using CORRFW
304   AliESDpid *workingPID = inH->GetESDpid();
305   if(workingPID){
306     AliDebug(1, "Using ESD PID from the input handler");
307     //printf("\n ESD PID\n");
308     fDePID->SetESDpid(workingPID);
309   } else { 
310     AliDebug(1, "Using default ESD PID");
311     //printf(" DEFAULT PID!\n\n");
312     fDePID->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData()));
313   }
314
315   if(!fDeCuts){
316     AliError("HFE cuts not available");
317     return;
318   }
319
320   // ---- CHOICE OF ANALYSIS ----
321   if(HasMCData()){
322     // Protect against missing MC trees
323     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
324     if(!mcH){
325       AliError("MC handler cuts not available");
326       return;
327     }
328
329     if(!mcH->InitOk()) return;
330     if(!mcH->TreeK()) return;
331     if(!mcH->TreeTR()) return;            
332     
333     AliDebug(4, Form("MC Event: %p", fMCEvent));
334     if(!fMCEvent){
335       AliError("No MC Event, but MC Data required");
336       return;
337     }
338     
339     ProcessMC(); // PURE MC
340     
341     if(IsESDanalysis()) ProcessESD(); // ESD WITH MC
342     
343   }else if(IsESDanalysis()) ProcessData(); // PURE ESD
344   
345
346   fDeNEvents->Fill(1);  
347   PostData(1, fHistDisplacedElectrons);
348   PostData(2, fDeCorrection);
349   PostData(3, fDeQA);
350   
351 }
352
353 //____________________________________________________________
354 void AliAnalysisTaskDisplacedElectrons::ProcessMC(){
355   //
356   // handle pure MC analysis
357   //
358
359   Int_t nMCelectrons = 0;
360   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
361   if(!fESD){
362     AliError("ESD event not available");
363     return;
364   }
365
366   Double_t mcContainer[4];   // container for the output in THnSparse
367   memset(mcContainer, 0, sizeof(Double_t) * 4);
368
369   fDeCFM->SetMCEventInfo(fMCEvent);
370
371   Double_t nContributor[1] = {0};
372   const AliVVertex *mcPrimVtx = fMCEvent->GetPrimaryVertex();
373   if(mcPrimVtx) nContributor[0] = mcPrimVtx->GetNContributors();
374   
375   // 
376   // cut at MC event level
377   //
378
379   if(!fDeCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent)) return;
380   if(GetPlugin(kCorrection)) fDeCFM->GetEventContainer()->Fill(nContributor,AliHFEcuts::kEventStepGenerated);
381   
382   AliStack *stack = 0x0;
383   
384   if(!fMCEvent->Stack())return;
385   stack = fMCEvent->Stack();
386   Int_t nTracks = stack->GetNtrack();
387
388   AliMCParticle *mcTrack = 0x0;
389
390   for(Int_t itrack = 0; itrack<nTracks; itrack++){
391     if(!(stack->Particle(itrack))) continue;
392     if(mcTrack) mcTrack = 0x0;
393     mcTrack = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(itrack));
394     if(!mcTrack) continue;
395     //TParticle *mcPart = stack->Particle(itrack);
396         
397     mcContainer[0] = mcTrack->Pt();
398     mcContainer[1] = mcTrack->Eta();
399     mcContainer[2] = mcTrack->Phi();
400     mcContainer[3] = mcTrack->Charge();
401     
402
403
404     if (!stack->IsPhysicalPrimary(mcTrack->GetLabel())) continue;
405     // no cut but require primary
406     if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(mcContainer, 0);    
407         
408     // all pions for reference
409     if(TMath::Abs(mcTrack->Particle()->GetPdgCode())==AliHFEdisplacedElectrons::kPDGpion && GetPlugin(kCorrection))
410       fDeCFM->GetParticleContainer()->Fill(mcContainer, 1);
411     
412     // cut for signal: all MC electrons
413     if(TMath::Abs(mcTrack->Particle()->GetPdgCode())==AliHFEdisplacedElectrons::kPDGelectron && GetPlugin(kCorrection)) 
414       fDeCFM->GetParticleContainer()->Fill(mcContainer, 2);
415
416     // cut at track level kinematics: pt and eta
417     if(TMath::Abs(mcContainer[1])>=0.8 || mcContainer[0]>20 || mcContainer[0]<0.1) continue;
418
419     if(TMath::Abs(mcTrack->Particle()->GetPdgCode())==AliHFEdisplacedElectrons::kPDGelectron){
420       nMCelectrons++;
421       fElectronsMcPt->Fill(mcContainer[0]);
422     }  
423
424     if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(mcContainer, 3);
425
426     
427     // fill MC THnSparse
428     fDisplacedElectrons->FillMcOutput(fESD, fMCEvent, mcTrack);
429     
430   }  // mc track loop 
431
432   if(fDeDebugLevel>=10) printf("there are %d electrons in this MC event", nMCelectrons);
433   
434 }
435
436 //____________________________________________________________
437 void AliAnalysisTaskDisplacedElectrons::ProcessESD(){
438   
439   // this is to handle ESD tracks with MC information
440   // MC pid is only used when HFE pid is implemented, for comparison
441   // corrections are taken into account
442   
443   // process data: ESD tracks with MC information
444   
445   const Int_t kStepPID = AliHFEcuts::kStepHFEcutsTRD + 1;
446
447   Double_t esdContainer[4];   // container for the output in THnSparse
448   memset(esdContainer, 0, sizeof(Double_t) * 4);
449   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
450   if(!fESD){
451     AliError("No ESD event available");
452     return;
453   }
454
455   fDeCFM->SetRecEventInfo(fESD);
456   Double_t nContrib[1] = {fESD->GetPrimaryVertex()->GetNContributors()};
457
458   Bool_t alreadyseen = kFALSE;
459   AliLabelContainer cont(fESD->GetNumberOfTracks());
460   
461   Int_t nHFEelectrons = 0;  
462   AliESDtrack *track = 0x0;    
463   AliStack *stack = 0x0;
464
465   if(!(stack = fMCEvent->Stack())) return;
466   
467   //
468   // cut at ESD event level
469   //
470   if(!fDeCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
471   if(GetPlugin(kCorrection)) fDeCFM->GetEventContainer()->Fill(nContrib, AliHFEcuts::kEventStepReconstructed);
472   
473   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
474     track = fESD->GetTrack(itrack);
475     
476     if(GetPlugin(kDisplacedElectrons)) {
477       
478       esdContainer[0] = track->Pt();
479       esdContainer[1] = track->Eta();
480       esdContainer[2] = track->Phi();
481       esdContainer[3] = track->Charge();
482       
483       // before any cut
484       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));  
485       cont.Append(TMath::Abs(track->GetLabel()));  // check double counting
486       if(alreadyseen) continue;  // avoid double counting
487       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepRecNoCut + AliHFEcuts::kNcutStepsMCTrack);       
488       
489       // 1st track cut
490       // RecKine: ITSTPC cuts : ITS & TPC refit, covmatrix: (2, 2, 0.5, 0.5, 2); min_tpccls: 50, chi2_tpccls: 3.5
491       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
492       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack);
493       
494       // 2nd track cut
495       // RecPrim: cut on track quality : DCA to vertex max: 3cm and 10cm; reject kink daughters
496       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
497       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack);
498       
499       // 3rd track cut
500       // HFEcuts: ITS layers cuts: ITS pixel layer: kFirst, kSecond, kBoth, kNone or kAny
501       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
502       if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack);
503
504       /*
505       //  4th track cut
506       // TRD: number of tracklets in TRD
507       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
508       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack);
509       */
510       
511       // 5th track cut
512       // track accepted, do PID 
513       // --> only electron candidate will be processed
514       AliHFEpidObject hfetrack;
515       hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
516       hfetrack.SetRecTrack(track);
517       //if(HasMCData())hfetrack.SetMCTrack(mctrack);
518       
519       if(!fDePID->IsSelected(&hfetrack)) continue;
520       else if(fDeDebugLevel>=10)
521         AliInfo("ESD info: this particle is identified as electron by HFEpid method \n");
522       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(esdContainer, 1+kStepPID + AliHFEcuts::kNcutStepsMCTrack);
523     
524       // Fill Containers
525       nHFEelectrons++;
526       fElectronsEsdPt->Fill(esdContainer[0]);
527       fDisplacedElectrons->FillEsdOutput(fESD, track, stack);
528     
529     }  // displaced electron analysis on ESD with MC plugin
530   } // track loop  
531   
532   if(fDeDebugLevel>=10) printf("there are %d HFE electrons in this ESD event", nHFEelectrons);
533   
534 }
535
536
537
538 //____________________________________________________________
539 void AliAnalysisTaskDisplacedElectrons::ProcessData(){
540
541   // this is a track loop over real data 
542   // no MC information at all 
543   // HFE pid is used
544   
545   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
546   //const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
547   const Int_t kStepPID = AliHFEcuts::kStepHFEcutsTRD + 1;
548
549   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
550   if(!fESD){
551     AliError("No ESD event available");
552     return;
553   }
554
555   Double_t dataContainer[4];   // container for the output in THnSparse
556   memset(dataContainer, 0, sizeof(Double_t) * 4);
557   
558   Bool_t alreadyseen = kFALSE;
559   AliLabelContainer cont(fESD->GetNumberOfTracks());
560
561
562   AliESDtrack *track = 0x0;
563   Int_t nHFEelectrons= 0;
564   
565   fDeCFM->SetRecEventInfo(fESD);
566   Double_t nContrib[1] = {fESD->GetPrimaryVertex()->GetNContributors()};
567   if(!fDeCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
568   if(GetPlugin(kCorrection)) fDeCFM->GetEventContainer()->Fill(nContrib, AliHFEcuts::kEventStepReconstructed);
569   
570   
571   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
572     track = fESD->GetTrack(itrack);
573     
574     if(GetPlugin(kDisplacedElectrons)) {
575       
576       dataContainer[0] = track->Pt();
577       dataContainer[1] = track->Eta();
578       dataContainer[2] = track->Phi();
579       dataContainer[3] = track->Charge();
580       
581       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));  // double counted track
582       cont.Append(TMath::Abs(track->GetLabel()));
583       if(alreadyseen) continue;  // avoid double counting
584       if(GetPlugin(kCorrection))fDeCFM->GetParticleContainer()->Fill(&dataContainer[4], 
585                                                                      1+AliHFEcuts::kStepRecNoCut + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
586
587       // 1st track cut
588       // RecKine: ITSTPC cuts : ITS & TPC refit, covmatrix: (2, 2, 0.5, 0.5, 2); min_tpccls: 50, chi2_tpccls: 3.5
589       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
590       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(&dataContainer[4], 
591                                                                       1+AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
592       
593       // 2nd track cut
594       // RecPrim: cut on track quality : DCA to vertex max: 3cm and 10cm; reject kink daughters
595       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track))  continue;
596       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(&dataContainer[4], 
597                                                                       1+AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
598         
599       //  3rd track cut
600       // HFEcuts: ITS layers cuts: ITS pixel layer: kFirst, kSecond, kBoth, kNone or kAny
601       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
602       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(&dataContainer[4],
603                                                                       1+AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
604       
605       /*
606       //  4th track cut
607       // TRD: number of tracklets in TRD0
608       if(!fDeCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
609       if(GetPlugin(kCorrection)) if(HasMCData())fDeCFM->GetParticleContainer()->Fill(&dataContainer[4], 
610                                                                                      1+AliHFEcuts::kStepHFEcutsTRD + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
611       */
612
613
614       // 5th track cut
615       // track accepted, do PID --> only electron candidate will be processed
616
617       AliHFEpidObject hfetrack;
618       hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
619       hfetrack.SetRecTrack(track);
620       
621       if(!fDePID->IsSelected(&hfetrack)) continue;
622       else if(fDeDebugLevel>=10)
623         AliInfo("ESD info: this particle is identified as electron by HFEpid method \n");
624       if(GetPlugin(kCorrection)) fDeCFM->GetParticleContainer()->Fill(dataContainer,  1+kStepPID + AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack);
625
626       nHFEelectrons++;
627       fElectronsDataPt->Fill(dataContainer[0]);
628       fDisplacedElectrons->FillDataOutput(fESD, track);
629     } // analyze displaced electrons plugin switched on
630   } // track loop  
631
632   if(fDeDebugLevel>=10) printf("there are %d HFE electrons in this DATA event", nHFEelectrons);
633 }
634
635
636 //____________________________________________________________
637 void AliAnalysisTaskDisplacedElectrons::Terminate(Option_t *){
638   //
639   // Terminate not implemented at the moment
640   //  
641   
642   fHistDisplacedElectrons = dynamic_cast<TList *>(GetOutputData(1));
643   fDeCorrection = dynamic_cast<TList *>(GetOutputData(2));
644   fDeQA = dynamic_cast<TList *>(GetOutputData(3));
645   if(!fDeCorrection) AliError("correction list not available\n");
646   if(!fHistDisplacedElectrons) AliError("de list not available\n");
647   if(!fDeQA) AliError("qa list is not available\n");
648   
649   fHistDisplacedElectrons->Print();
650   fDeCorrection->Print();
651   fDeQA->Print();
652
653   AliInfo("analysis done!\n");
654   
655 }
656
657 //____________________________________________________________
658 void AliAnalysisTaskDisplacedElectrons::PrintStatus() const {
659   
660   //
661   // Print Analysis status
662   //
663   printf("\n");
664   printf("\t Analysis Settings\n\t========================================\n");
665   printf("\t running over %s\n", HasMCData()?"MC data":"pp collision data");
666   printf("\t displaced electrons' analysis is %s\n", GetPlugin(kDisplacedElectrons)?"ON":"OFF");
667   printf("\t correction container is %s\n", GetPlugin(kCorrection)?"ON":"OFF");
668   printf("\t hfe pid qa is %s\n", GetPlugin(kDePidQA)?"ON":"OFF");
669   printf("\t post processing  is %s\n", GetPlugin(kPostProcess)?"ON":"OFF");
670   printf("\t cuts: %s\n", (fDeCuts != NULL) ? "YES" : "NO");
671   printf("\t ");
672   printf("\n");
673 }
674
675 //__________________________________________                                                  
676 void AliAnalysisTaskDisplacedElectrons::SwitchOnPlugin(Int_t plug){
677   //                                            
678   // Switch on Plugin          
679   // Available:                                  
680   //  - analyze impact parameter
681   //  - Post Processing                                                                      
682   
683   switch(plug)
684     {
685     case kDisplacedElectrons: 
686       SETBIT(fDePlugins, plug); 
687       break;
688     case kCorrection:
689       SETBIT(fDePlugins, plug); 
690       break;
691     case kDePidQA:
692       SETBIT(fDePlugins, plug); 
693       break;
694     case kPostProcess: 
695       SETBIT(fDePlugins, plug); 
696       break;
697     default: 
698       AliError("Unknown Plugin");
699     };
700 }
701
702 //____________________________________________________________
703 void AliAnalysisTaskDisplacedElectrons::MakeParticleContainer(){
704   //
705   // Create the particle container for the correction framework manager and 
706   // link it
707   //
708   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
709   const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
710  
711   const Int_t kNvar   = 4;
712   //number of variables on the grid:pt,eta, phi, charge
713   const Double_t kPtbound[2] = {0.1, 10.};
714   const Double_t kEtabound[2] = {-0.8, 0.8};
715   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
716
717   //arrays for the number of bins in each dimension
718   Int_t iBin[kNvar];
719   iBin[0] = 40; // bins in pt
720   iBin[1] =  8; // bins in eta 
721   iBin[2] = 18; // bins in phi
722   iBin[3] =  2; // bins in charge
723
724   //arrays for lower bounds :
725   Double_t* binEdges[kNvar];
726   binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
727   binEdges[1] = AliHFEtools::MakeLinearBinning(iBin[1], kEtabound[0], kEtabound[1]);
728   binEdges[2] = AliHFEtools::MakeLinearBinning(iBin[2], kPhibound[0], kPhibound[1]);
729   binEdges[3] = AliHFEtools::MakeLinearBinning(iBin[3], -1.1, 1.1); // Numeric precision
730
731   //------------------------------------------------
732   //     one "container" for MC+ESD+Data          
733   //----------pure MC track-------------------------
734   // 0: MC generated
735   // 1: MC pion total  ---- be careful!!!!
736   // 2: MC electrons total
737   // 3: MC electrons in acceptance 
738   //-------ESD track with MC info-------------------
739   // 4: ESD track with MC: no cut
740   // 5: ESD track with MC: cut on kine its tpc
741   // 6: ESD track with MC: rec prim
742   // 7: ESD track with MC: hfe cuts its
743   // 8: ESD track with MC: hfe cuts trd
744   // 9: ESD track with MC: hfe pid 
745   //-----------data track---------------------------
746   // 10: DATA track wo MC: no cut
747   // 11: DATA track wo MC: cut on kine its tpc
748   // 12: DATA track wo MC: rec prim
749   // 13: DATA track wo MC: hfe cuts its
750   // 14: DATA track wo MC: hfe cuts trd
751   // 15: DATA track wo MC: hfe pid 
752   //------------------------------------------------
753
754   AliCFContainer* container = new AliCFContainer("deTrackContainer", "Container for tracks", 
755                                                  (1 + kNcutStepsTrack + kNcutStepsESDtrack), kNvar, iBin);
756   
757   //setting the bin limits
758   for(Int_t ivar = 0; ivar < kNvar; ivar++){
759     container -> SetBinLimits(ivar, binEdges[ivar]);
760   }
761   fDeCFM->SetParticleContainer(container);
762 }
763
764
765 //____________________________________________________________
766 void AliAnalysisTaskDisplacedElectrons::MakeEventContainer(){
767   //
768   // Create the event container for the correction framework and link it
769   //
770   
771   // event container
772   // 0: MC event
773   // 1: ESD event
774
775   const Int_t kNvar = 1;  // number of variables on the grid: number of tracks per event
776   const Double_t kNTrackBound[2] = {-0.5, 200.5};
777   const Int_t kNBins = 201;
778
779   AliCFContainer *evCont = new AliCFContainer("deEventContainer", "Container for DE events", AliHFEcuts::kNcutStepsEvent, kNvar, &kNBins);
780
781   Double_t trackBins[kNBins];
782   for(Int_t ibin = 0; ibin < kNBins; ibin++) trackBins[ibin] = kNTrackBound[0] + static_cast<Double_t>(ibin);
783   evCont->SetBinLimits(0,trackBins);
784
785   fDeCFM->SetEventContainer(evCont);
786
787 }
788
789
790
791 //__________________________________________________________
792 void AliAnalysisTaskDisplacedElectrons::AddPIDdetector(TString detector){
793   //
794   // Adding PID detector to the task
795   //
796   if(!fDePIDdetectors.Length()) 
797     fDePIDdetectors = detector;
798   else
799     fDePIDdetectors += ":" + detector;
800 }
801
802
803
804 //____________________________________________________________
805 AliAnalysisTaskDisplacedElectrons::AliLabelContainer::AliLabelContainer(Int_t capacity):
806   fContainer(NULL),
807   fBegin(NULL),
808   fEnd(NULL),
809   fLast(NULL),
810   fCurrent(NULL)
811 {
812   //
813   // Default constructor
814   //
815   fContainer = new Int_t[capacity];
816   fBegin = &fContainer[0];
817   fEnd = &fContainer[capacity - 1];
818   fLast = fCurrent = fBegin;
819 }
820
821 //____________________________________________________________
822 Bool_t AliAnalysisTaskDisplacedElectrons::AliLabelContainer::Append(Int_t label){
823   //
824   // Add Label to the container
825   //
826   if(fLast > fEnd) return kFALSE;
827   *fLast++ = label;
828   return kTRUE;
829 }
830
831 //____________________________________________________________
832 Bool_t AliAnalysisTaskDisplacedElectrons::AliLabelContainer::Find(Int_t label) const {
833   //
834   // Find track in the list of labels
835   //
836   for(Int_t *entry = fBegin; entry <= fLast; entry++) 
837     if(*entry == label) return kTRUE;
838   return kFALSE;
839 }
840
841 //____________________________________________________________
842 Int_t AliAnalysisTaskDisplacedElectrons::AliLabelContainer::Next() { 
843   //
844   // Mimic iterator
845   //
846   if(fCurrent > fLast) return -1; 
847   fCurrent++;
848   return *fCurrent;
849 }