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