cfd82d5c06a51fd235405d3089f5272cf699eea0
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisTaskNetParticle.cxx
1 //-*- Mode: C++ -*-
2
3 #include "TFile.h"
4 #include "TChain.h"
5 #include "TTree.h"
6 #include "TH1F.h"
7 #include "TH2F.h"
8 #include "TH3F.h"
9 #include "TCanvas.h"
10 #include "TStopwatch.h"
11 #include "TMath.h"
12 #include "THashList.h"
13
14 #include "AliAnalysisTask.h"
15 #include "AliAnalysisManager.h"
16 #include "AliTracker.h" 
17 #include "AliESDEvent.h"
18 #include "AliESDInputHandler.h"
19 #include "AliESDpid.h"
20 #include "AliStack.h"
21 #include "AliMCEvent.h"
22 #include "AliMCEventHandler.h"
23 #include "AliESDtrackCuts.h"
24 #include "AliKineTrackCuts.h"
25 #include "AliMCParticle.h"
26 #include "AliESDVZERO.h"
27 #include "AliAnalysisTaskNetParticle.h"
28 #include "AliGenEventHeader.h"
29 #include "AliCentrality.h"
30 #include "AliAODEvent.h"
31 #include "AliAODInputHandler.h"
32
33 using namespace std;
34
35 /**
36  * Class for NetParticle Distributions
37  * -- AnalysisTask
38  * Authors: Jochen Thaeder <jochen@thaeder.de>
39  *          Michael Weber <m.weber@cern.ch>
40  */
41
42 ClassImp(AliAnalysisTaskNetParticle)
43
44 /*
45  * ---------------------------------------------------------------------------------
46  *                            Constructor / Destructor
47  * ---------------------------------------------------------------------------------
48  */
49
50 //________________________________________________________________________
51 AliAnalysisTaskNetParticle::AliAnalysisTaskNetParticle(const char *name) :
52   AliAnalysisTaskSE(name),
53   fHelper(NULL),
54   fEffCont(NULL),
55   fDCA(NULL),
56   fDist(NULL),
57   fQA(NULL),
58
59   fOutList(NULL),
60   fOutListEff(NULL),
61   fOutListCont(NULL),
62   fOutListDCA(NULL),
63   fOutListQA(NULL),
64
65   fESD(NULL), 
66   fESDHandler(NULL),
67
68   fESDTrackCutsBase(NULL),
69   fESDTrackCuts(NULL),
70   fESDTrackCutsBkg(NULL),
71   fESDTrackCutsEff(NULL),
72
73   fAOD(NULL), 
74   fAODHandler(NULL),
75
76   fIsMC(kFALSE),
77   fIsAOD(kFALSE),
78   fESDTrackCutMode(0),
79   fModeEffCreation(0),
80   fModeDCACreation(0),
81   fModeDistCreation(0),
82   fModeQACreation(0),
83
84   fMCEvent(NULL),
85   fMCStack(NULL),
86
87   fEtaMax(0.8),
88   fEtaMaxEff(0.9),
89   fPtRange(),
90   fPtRangeEff(),
91
92   fAODtrackCutBit(1024) {
93   // Constructor   
94
95   AliLog::SetClassDebugLevel("AliAnalysisTaskNetParticle",10);
96
97   fPtRange[0] = 0.4;
98   fPtRange[1] = 0.8;
99   fPtRangeEff[0] = 0.2;
100   fPtRangeEff[1] = 1.6;
101
102   // -- Output slots
103   // -------------------------------------------------
104   DefineOutput(1, TList::Class());
105   DefineOutput(2, TList::Class());
106   DefineOutput(3, TList::Class());
107   DefineOutput(4, TList::Class());
108   DefineOutput(5, TList::Class());
109 }
110
111 //________________________________________________________________________
112 AliAnalysisTaskNetParticle::~AliAnalysisTaskNetParticle() {
113   // Destructor
114
115   if (fESDTrackCutsBase) delete fESDTrackCutsBase;
116   if (fESDTrackCuts)     delete fESDTrackCuts;
117   if (fESDTrackCutsBkg)  delete fESDTrackCutsBkg;
118   if (fESDTrackCutsEff)  delete fESDTrackCutsEff;
119
120   if (fEffCont)          delete fEffCont;
121   if (fDCA)              delete fDCA;
122   if (fDist)             delete fDist;
123   if (fQA)               delete fQA;
124   if (fHelper)           delete fHelper;
125 }
126
127 /*
128  * ---------------------------------------------------------------------------------
129  *                                 Public Methods
130  * ---------------------------------------------------------------------------------
131  */
132
133 //________________________________________________________________________
134 void AliAnalysisTaskNetParticle::UserCreateOutputObjects() {
135   // Create histograms
136
137   // ------------------------------------------------------------------
138   // -- Create Output Lists
139   // ------------------------------------------------------------------
140   Bool_t oldStatus = TH1::AddDirectoryStatus();
141   TH1::AddDirectory(kFALSE);
142
143   fOutList = new TList;
144   fOutList->SetName(GetName()) ;
145   fOutList->SetOwner(kTRUE);
146  
147   fOutListEff = new TList;
148   fOutListEff->SetName(Form("%s_eff",GetName()));
149   fOutListEff->SetOwner(kTRUE) ;
150
151   fOutListCont = new TList;
152   fOutListCont->SetName(Form("%s_cont",GetName()));
153   fOutListCont->SetOwner(kTRUE) ;
154
155   fOutListDCA = new TList;
156   fOutListDCA->SetName(Form("%s_dca",GetName()));
157   fOutListDCA->SetOwner(kTRUE) ;
158  
159   fOutListQA = new TList;
160   fOutListQA->SetName(Form("%s_qa",GetName()));
161   fOutListQA->SetOwner(kTRUE) ;
162
163   // ------------------------------------------------------------------
164   // -- Initialize all classes
165   // ------------------------------------------------------------------
166   Initialize();
167  
168   // ------------------------------------------------------------------
169   // -- Get event / trigger statistics histograms
170   // ------------------------------------------------------------------
171   fOutList->Add(fHelper->GetHEventStat0());
172   fOutList->Add(fHelper->GetHEventStat1());
173   fOutList->Add(fHelper->GetHTriggerStat());
174   fOutList->Add(fHelper->GetHCentralityStat());
175
176   // ------------------------------------------------------------------
177   // -- Add histograms from efficiency/contamination class
178   // ------------------------------------------------------------------
179   if ((fIsAOD||fIsMC) && fModeEffCreation == 1) {
180     fOutListEff->Add(fEffCont->GetHnEff());
181     fOutListCont->Add(fEffCont->GetHnCont());
182   }
183
184   // ------------------------------------------------------------------
185   // -- Add histograms from DCA class
186   // ------------------------------------------------------------------
187   if (fModeDCACreation == 1)
188     fOutListDCA->Add(fDCA->GetHnDCA());
189
190   // ------------------------------------------------------------------
191   // -- Add histograms from QA class
192   // ------------------------------------------------------------------
193   if (fModeQACreation == 1)
194     fOutListQA->Add(fQA->GetHnQA());
195
196   // ------------------------------------------------------------------
197
198   TH1::AddDirectory(oldStatus);
199
200   PostData(1,fOutList);
201   PostData(2,fOutListEff);
202   PostData(3,fOutListCont);
203   PostData(4,fOutListDCA);
204   PostData(5,fOutListQA);
205
206   return;
207 }
208
209 //________________________________________________________________________
210 void AliAnalysisTaskNetParticle::UserExec(Option_t *) {
211   // Called for each event
212
213   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
214   // -- Setup Event
215   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
216   if (SetupEvent() < 0) {
217     PostData(1,fOutList);
218     PostData(2,fOutListEff);
219     PostData(3,fOutListCont);
220     PostData(4,fOutListDCA);
221     PostData(5,fOutListQA);
222     return;
223   }
224
225   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
226   // -- Process Efficiency / Contamination Determination
227   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
228   if ((fIsMC||fIsAOD) && fModeEffCreation == 1)
229     fEffCont->Process();
230
231   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
232   // -- Process DCA Determination
233   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
234   if (fModeDCACreation == 1)
235     fDCA->Process();
236
237   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
238   // -- Process Distributions 
239   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
240   if (fModeDistCreation == 1)
241     fDist->Process();
242
243   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
244   // -- Fill QA histograms
245   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
246   if (fModeQACreation == 1)
247     fQA->Process();
248
249   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
250   // -- Post output data
251   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
252
253   PostData(1,fOutList);
254   PostData(2,fOutListEff);
255   PostData(3,fOutListCont);
256   PostData(4,fOutListDCA);
257   PostData(5,fOutListQA);
258
259   return;
260 }      
261
262 //________________________________________________________________________
263 void AliAnalysisTaskNetParticle::Terminate(Option_t *){
264   // Terminate
265 }
266
267 /*
268  * ---------------------------------------------------------------------------------
269  *                                 Public Methods
270  * ---------------------------------------------------------------------------------
271  */
272
273 //________________________________________________________________________
274 Int_t AliAnalysisTaskNetParticle::Initialize() {
275   // Initialize event
276
277   // ------------------------------------------------------------------
278   // -- ESD TrackCuts
279   // ------------------------------------------------------------------
280   TString sModeName("");
281   
282   // -- Create ESD track cuts
283   // --------------------------
284   fESDTrackCutsBase = new AliESDtrackCuts;
285   
286   if (fESDTrackCutMode == 0) {
287     fESDTrackCutsBase->SetMinNCrossedRowsTPC(70);                                             // TPC
288     fESDTrackCutsBase->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);                    // TPC
289   }
290   else if (fESDTrackCutMode == 1) {
291     fESDTrackCutsBase->SetMinNClustersTPC(70);                                                // TPC  2010
292   }
293
294   fESDTrackCutsBase->SetMaxChi2PerClusterTPC(4);                                              // TPC  2010
295   fESDTrackCutsBase->SetAcceptKinkDaughters(kFALSE);                                          // TPC  2010
296   fESDTrackCutsBase->SetRequireTPCRefit(kTRUE);                                               // TPC  2010
297
298   if (fESDTrackCutMode == 0) {
299     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff); // ITS
300     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kOff); // ITS
301     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSSD,AliESDtrackCuts::kOff); // ITS
302   } 
303   else if (fESDTrackCutMode == 1) {
304     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); // ITS 2010
305   //  fESDTrackCutsBase->SetMinNClustersITS(4);
306   }
307
308   fESDTrackCutsBase->SetRequireITSRefit(kTRUE);                                               // ITS 2010 
309   fESDTrackCutsBase->SetMaxChi2PerClusterITS(36);                                             // ITS 2010
310
311   fESDTrackCutsBase->SetDCAToVertex2D(kFALSE);                                                // VertexConstrained  2010
312   fESDTrackCutsBase->SetRequireSigmaToVertex(kFALSE);                                         // VertexConstrained  2010
313   fESDTrackCutsBase->SetMaxDCAToVertexZ(2);                                                   // VertexConstrained  2010
314  
315   fESDTrackCutsBase->SetEtaRange(-1.*fEtaMax, fEtaMax);                                       // Acceptance
316   fESDTrackCutsBase->SetPtRange(fPtRange[0],fPtRange[1]);                                     // Acceptance
317
318   // -- Mode : standard cuts
319   if (fESDTrackCutMode == 0) 
320     sModeName = "Std";
321   // -- Mode : for comparison to LF
322   else if (fESDTrackCutMode == 1)
323     sModeName = "LF";
324   // -- Mode : Default
325   else
326     sModeName = "Base";
327   
328   fESDTrackCutsBase->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
329
330   // -- Create ESD track cuts -> Base + DCA
331   // ------------------------------
332   fESDTrackCuts = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
333   fESDTrackCuts->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
334   if (fESDTrackCutMode == 0) 
335     fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");                       // 2010 VertexConstrained  ->  7*(0.0026+0.0050/pt^1.01)
336     //    fESDTrackCuts->SetMaxDCAToVertexXY(0.3);
337   else if (fESDTrackCutMode == 1)
338     fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");                       // 2010 VertexConstrained  ->  7*(0.0026+0.0050/pt^1.01)
339
340   //  fESDTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);                                    // golden cut off
341
342   // -- Create ESD BKG track cuts -> Base + Acceptance(Eff)
343   // ------------------------------
344   fESDTrackCutsBkg = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
345   fESDTrackCutsBkg->SetName(Form("NetParticleCuts2010_%s_Bkg",sModeName.Data()));
346   fESDTrackCutsBkg->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]);                              // Acceptance
347   fESDTrackCutsBkg->SetEtaRange(-1.*fEtaMaxEff, fEtaMaxEff);                                // Acceptance
348   
349   // -- Create ESD Eff track cuts -> Base + DCA + Acceptance(Eff)
350   // ------------------------------
351   fESDTrackCutsEff = static_cast<AliESDtrackCuts*>(fESDTrackCuts->Clone());
352   fESDTrackCutsEff->SetName(Form("NetParticleCuts2010_%s_Eff",sModeName.Data()));
353   fESDTrackCutsEff->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]);                              // Acceptance
354   fESDTrackCutsEff->SetEtaRange(-1.*fEtaMaxEff, fEtaMaxEff);                                // Acceptance
355
356   // ------------------------------------------------------------------
357   // -- Initialize Helper
358   // ------------------------------------------------------------------
359   if (fHelper->Initialize(fESDTrackCutsEff, fIsMC, fAODtrackCutBit, fModeDistCreation))
360     return -1;
361
362   // ------------------------------------------------------------------
363   // -- Create / Initialize Efficiency/Contamination
364   // ------------------------------------------------------------------
365   if ((fIsMC||fIsAOD) && fModeEffCreation == 1) {
366     fEffCont = new AliAnalysisNetParticleEffCont;
367     fEffCont->Initialize(fHelper);
368   }
369
370   // ------------------------------------------------------------------
371   // -- Create / Initialize DCA Determination
372   // ------------------------------------------------------------------
373   if (fModeDCACreation == 1) {
374     fDCA = new AliAnalysisNetParticleDCA;
375     fDCA->SetESDTrackCutsBkg(fESDTrackCutsBkg);
376     fDCA->Initialize(fHelper);
377   }
378
379   // ------------------------------------------------------------------
380   // -- Create / Initialize Distribution Determination
381   // ------------------------------------------------------------------
382   if (fModeDistCreation == 1) {
383     fDist = new AliAnalysisNetParticleDistribution;
384     fDist->SetOutList(fOutList);
385     fDist->Initialize(fHelper, fESDTrackCuts);
386   }
387
388   // ------------------------------------------------------------------
389   // -- Create / Initialize QA Determination
390   // ------------------------------------------------------------------
391   if (fModeQACreation == 1) {
392     fQA = new AliAnalysisNetParticleQA();
393     fQA->Initialize(fHelper);
394   }
395
396   // ------------------------------------------------------------------
397   // -- Reset Event
398   // ------------------------------------------------------------------
399   ResetEvent();
400
401   return 0;
402 }
403
404 /*
405  * ---------------------------------------------------------------------------------
406  *                            Setup/Reset Methods - private
407  * ---------------------------------------------------------------------------------
408  */
409
410 //________________________________________________________________________
411 Int_t AliAnalysisTaskNetParticle::SetupEvent() {
412   // Setup Reading of event
413   // > return 0 for success / accepted event
414   // > return -1 for failed setup
415   // > return -2 for rejected event
416
417   ResetEvent();
418
419   // -- ESD Event
420   // ------------------------------------------------------------------
421   if (!fIsAOD && SetupESDEvent() < 0) {
422     AliError("Setup ESD Event failed");
423     return -1;
424   }
425
426   // -- AOD Event
427   // ------------------------------------------------------------------
428   if (fIsAOD && SetupAODEvent() < 0) {
429     AliError("Setup AOD Event failed");
430     return -1;
431   }
432   
433   // -- Setup MC Event
434   // ------------------------------------------------------------------
435   if (fIsMC && SetupMCEvent() < 0) {
436     AliError("Setup MC Event failed");
437     return -1;
438   }
439
440   // -- Setup Event for Helper / EffCont  / DCA / Dist / QA classes
441   // ------------------------------------------------------------------
442   fHelper->SetupEvent(fESDHandler, fAODHandler, fMCEvent);
443
444   if (fModeEffCreation && (fIsMC || fIsAOD) )
445     fEffCont->SetupEvent(); 
446
447   if (fModeDCACreation == 1)
448     fDCA->SetupEvent();
449
450   if (fModeDistCreation == 1)
451     fDist->SetupEvent(); 
452
453   if (fModeQACreation == 1)
454     fQA->SetupEvent(); 
455
456   // -- Evaluate Event cuts
457   // ------------------------------------------------------------------
458   return fHelper->IsEventRejected() ? -2 : 0;
459 }
460
461 //________________________________________________________________________
462 Int_t AliAnalysisTaskNetParticle::SetupESDEvent() {
463   // -- Setup ESD Event
464   // > return 0 for success 
465   // > return -1 for failed setup
466
467   fESDHandler= dynamic_cast<AliESDInputHandler*> 
468     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
469   if (!fESDHandler) {
470     AliError("Could not get ESD input handler");
471     return -1;
472   } 
473
474   fESD = fESDHandler->GetEvent();
475   if (!fESD) {
476     AliError("Could not get ESD event");
477     return -1;
478   }
479
480   // -- Check PID response
481   // ------------------------------------------------------------------
482   if (!fESDHandler->GetPIDResponse()) {
483     AliError("Could not get PID response");
484     return -1;
485   } 
486
487   // -- Check Vertex
488   // ------------------------------------------------------------------
489   if (!fESD->GetPrimaryVertexTracks()) {
490     AliError("Could not get vertex from tracks");
491     return -1;
492   }
493
494   // -- Check Centrality
495   // ------------------------------------------------------------------
496   if (!fESD->GetCentrality()) {
497     AliError("Could not get centrality");
498     return -1;
499   }
500
501   return 0;
502 }
503
504 //________________________________________________________________________
505 Int_t AliAnalysisTaskNetParticle::SetupAODEvent() {
506   // -- Setup AOD Event
507   // > return 0 for success 
508   // > return -1 for failed setup
509
510   fAODHandler= dynamic_cast<AliAODInputHandler*> 
511     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
512   if (!fAODHandler) {
513     AliError("Could not get AOD input handler");
514     return -1;
515   } 
516
517   fAOD = fAODHandler->GetEvent();
518   if (!fAOD) {
519     AliError("Could not get AOD event");
520     return -1;
521   }
522
523   // -- Check PID response
524   // ------------------------------------------------------------------
525   if (!fAODHandler->GetPIDResponse()) {
526     AliError("Could not get PID response");
527     return -1;
528   } 
529
530   // -- Check Vertex
531   // ------------------------------------------------------------------
532   if (!fAOD->GetPrimaryVertex()) {
533     AliError("Could not get primary vertex");
534     return -1;
535   }
536
537   // -- Check Centrality
538   // ------------------------------------------------------------------
539   if (!fAOD->GetHeader()->GetCentralityP()) {
540     AliError("Could not get centrality");
541     return -1;
542   }
543
544   return 0;
545 }
546
547 //________________________________________________________________________
548 Int_t AliAnalysisTaskNetParticle::SetupMCEvent() {
549   // -- Setup MC Event
550   // > return 0 for success 
551   // > return -1 for failed setup
552
553   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> 
554     (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
555   
556   if (!mcH) {
557     AliError("MC event handler not available");
558     return -1;
559   }
560
561   fMCEvent = mcH->MCEvent();
562   if (!fMCEvent) {
563     AliError("MC event not available");
564     return -1;
565   }
566
567   // -- Get MC header
568   // ------------------------------------------------------------------
569   AliHeader* header = fMCEvent->Header();
570   if (!header) {
571     AliError("MC header not available");
572     return -1;
573   }
574
575   // -- Check Stack
576   // ------------------------------------------------------------------
577   fMCStack = fMCEvent->Stack(); 
578   if (!fMCStack) {
579     AliError("MC stack not available");
580     return -1;
581   }
582     
583   // -- Check GenHeader
584   // ------------------------------------------------------------------
585   if (!header->GenEventHeader()) {
586     AliError("Could not retrieve genHeader from header");
587     return -1;
588   }
589
590   // -- Check primary vertex
591   // ------------------------------------------------------------------
592   if (!fMCEvent->GetPrimaryVertex()){
593     AliError("Could not get MC vertex");
594     return -1;
595   }
596
597   return 0;
598 }
599
600 //________________________________________________________________________
601 void AliAnalysisTaskNetParticle::ResetEvent() {
602   // -- Reset event
603   
604   // -- Reset ESD Event
605   fESD       = NULL;
606
607   // -- Reset AOD Event
608   fAOD       = NULL;
609
610   // -- Reset MC Event
611   if (fIsMC)
612     fMCEvent = NULL;
613
614   // -- Reset Dist Creation 
615   if (fModeDistCreation == 1)
616     fDist->ResetEvent();
617
618   return;
619 }
620
621