Added some more scripts
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicityTask.cxx
1 #include "AliForwardMultiplicityTask.h"
2 #include "AliTriggerAnalysis.h"
3 #include "AliPhysicsSelection.h"
4 #include "AliLog.h"
5 // #include "AliFMDAnaParameters.h"
6 #include "AliESDEvent.h"
7 #include "AliAODHandler.h"
8 #include "AliMultiplicity.h"
9 #include "AliInputEventHandler.h"
10 #include "AliForwardCorrectionManager.h"
11 #include "AliAnalysisManager.h"
12 #include <TH1.h>
13 #include <TDirectory.h>
14 #include <TTree.h>
15 #include <TROOT.h>
16 #include <iostream>
17 #include <iomanip>
18
19 //====================================================================
20 AliForwardMultiplicityTask::AliForwardMultiplicityTask()
21   : AliAnalysisTaskSE(),
22     fEnableLowFlux(true),
23     fHData(0),
24     fFirstEvent(true),
25     fESDFMD(),
26     fHistos(),
27     fAODFMD(),
28     fEventInspector(),
29     fEnergyFitter(),
30     fSharingFilter(),
31     fDensityCalculator(),
32     fCorrections(),
33     fHistCollector(),
34     fList(0)
35 {
36 }
37
38 //____________________________________________________________________
39 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
40   : AliAnalysisTaskSE(name), 
41     fEnableLowFlux(true),
42     fHData(0),
43     fFirstEvent(true),
44     fESDFMD(),
45     fHistos(),
46     fAODFMD(kTRUE),
47     fEventInspector("event"),
48     fEnergyFitter("energy"),
49     fSharingFilter("sharing"), 
50     fDensityCalculator("density"),
51     fCorrections("corrections"),
52     fHistCollector("collector"),
53     fList(0)
54 {
55   DefineOutput(1, TList::Class());
56 }
57
58 //____________________________________________________________________
59 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
60   : AliAnalysisTaskSE(o),
61     fEnableLowFlux(o.fEnableLowFlux),
62     fHData(o.fHData),
63     fFirstEvent(true),
64     fESDFMD(o.fESDFMD),
65     fHistos(o.fHistos),
66     fAODFMD(o.fAODFMD),
67     fEventInspector(o.fEventInspector),
68     fEnergyFitter(o.fEnergyFitter),
69     fSharingFilter(o.fSharingFilter),
70     fDensityCalculator(o.fDensityCalculator),
71     fCorrections(o.fCorrections),
72     fHistCollector(o.fHistCollector),
73     fList(o.fList) 
74 {
75   DefineOutput(1, TList::Class());
76 }
77
78 //____________________________________________________________________
79 AliForwardMultiplicityTask&
80 AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
81 {
82   AliAnalysisTaskSE::operator=(o);
83
84   fEnableLowFlux     = o.fEnableLowFlux;
85   fHData             = o.fHData;
86   fFirstEvent        = o.fFirstEvent;
87   fEventInspector    = o.fEventInspector;
88   fEnergyFitter      = o.fEnergyFitter;
89   fSharingFilter     = o.fSharingFilter;
90   fDensityCalculator = o.fDensityCalculator;
91   fCorrections       = o.fCorrections;
92   fHistCollector     = o.fHistCollector;
93   fHistos            = o.fHistos;
94   fAODFMD            = o.fAODFMD;
95   fList              = o.fList;
96
97   return *this;
98 }
99
100 //____________________________________________________________________
101 void
102 AliForwardMultiplicityTask::SetDebug(Int_t dbg)
103 {
104   fEventInspector.SetDebug(dbg);
105   fEnergyFitter.SetDebug(dbg);
106   fSharingFilter.SetDebug(dbg);
107   fDensityCalculator.SetDebug(dbg);
108   fCorrections.SetDebug(dbg);
109   fHistCollector.SetDebug(dbg);
110 }
111
112 //____________________________________________________________________
113 void
114 AliForwardMultiplicityTask::Init()
115 {
116   fFirstEvent = true;
117 }
118
119 //____________________________________________________________________
120 void
121 AliForwardMultiplicityTask::InitializeSubs()
122 {
123
124   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
125   // pars->Init(kTRUE);
126
127   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
128   fcm.Init(fEventInspector.GetCollisionSystem(), 
129            fEventInspector.GetEnergy(),
130            fEventInspector.GetField());
131   // Check that we have the energy loss fits, needed by 
132   //   AliFMDSharingFilter 
133   //   AliFMDDensityCalculator 
134   if (!fcm.GetELossFit()) { 
135     AliFatal(Form("No energy loss fits"));
136     return;
137   }
138   // Check that we have the double hit correction - (optionally) used by 
139   //  AliFMDDensityCalculator 
140   if (!fcm.GetDoubleHit()) {
141     AliWarning("No double hit corrections"); 
142   }
143   // Check that we have the secondary maps, needed by 
144   //   AliFMDCorrections 
145   //   AliFMDHistCollector
146   if (!fcm.GetSecondaryMap()) {
147     AliFatal("No secondary corrections");
148     return;
149   }
150   // Check that we have the vertex bias correction, needed by 
151   //   AliFMDCorrections 
152   if (!fcm.GetVertexBias()) { 
153     AliFatal("No event vertex bias corrections");
154     return;
155   }
156   // Check that we have the merging efficiencies, optionally used by 
157   //   AliFMDCorrections 
158   if (!fcm.GetMergingEfficiency()) {
159     AliWarning("No merging efficiencies");
160   }
161
162
163   const TAxis* pe = fcm.GetEtaAxis();
164   const TAxis* pv = fcm.GetVertexAxis();
165   if (!pe) AliFatal("No eta axis defined");
166   if (!pv) AliFatal("No vertex axis defined");
167
168   // TAxis e(pars->GetNetaBins(),  pars->GetEtaMin(),  pars->GetEtaMax());
169   // TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
170   // pv->Dump();
171
172   fHistos.Init(*pe);
173   fAODFMD.Init(*pe);
174
175   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
176   fHData->SetStats(0);
177   fHData->SetDirectory(0);
178   fList->Add(fHData);
179
180   fEnergyFitter.Init(*pe);
181   fEventInspector.Init(*pv);
182   fHistCollector.Init(*pv);
183
184   this->Print();
185 }
186
187 //____________________________________________________________________
188 void
189 AliForwardMultiplicityTask::UserCreateOutputObjects()
190 {
191   fList = new TList;
192
193   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
194   AliAODHandler*      ah = 
195     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
196   if (!ah) AliFatal("No AOD output handler set in analysis manager");
197     
198     
199   TObject* obj = &fAODFMD;
200   ah->AddBranch("AliAODForwardMult", &obj);
201
202   fEventInspector.DefineOutput(fList);
203   fEnergyFitter.DefineOutput(fList);
204   fSharingFilter.DefineOutput(fList);
205   fDensityCalculator.DefineOutput(fList);
206   fCorrections.DefineOutput(fList);
207
208   PostData(1, fList);
209 }
210 //____________________________________________________________________
211 void
212 AliForwardMultiplicityTask::UserExec(Option_t*)
213 {
214   // static Int_t cnt = 0;
215   // cnt++;
216   // Get the input data 
217   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
218   // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
219   if (!esd) { 
220     AliWarning("No ESD event found for input event");
221     return;
222   }
223
224   // On the first event, initialize the parameters 
225   if (fFirstEvent && esd->GetESDRun()) { 
226     fEventInspector.ReadRunDetails(esd);
227     
228     AliInfo(Form("Initializing with parameters from the ESD:\n"
229                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
230                  "         AliESDEvent::GetBeamType()     ->%s\n"
231                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
232                  "         AliESDEvent::GetMagneticField()->%f\n"
233                  "         AliESDEvent::GetRunNumber()    ->%d\n",
234                  esd->GetBeamEnergy(), 
235                  esd->GetBeamType(),
236                  esd->GetCurrentL3(), 
237                  esd->GetMagneticField(),
238                  esd->GetRunNumber()));
239
240               
241
242     // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
243     // pars->SetParametersFromESD(esd);
244     // pars->PrintStatus();
245     fFirstEvent = false;
246
247     InitializeSubs();
248   }
249   // Clear stuff 
250   fHistos.Clear();
251   fESDFMD.Clear();
252   fAODFMD.Clear();
253
254   Bool_t   lowFlux  = kFALSE;
255   UInt_t   triggers = 0;
256   UShort_t ivz      = 0;
257   Double_t vz       = 0;
258   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
259   if (found & AliFMDEventInspector::kNoEvent)    return;
260   if (found & AliFMDEventInspector::kNoTriggers) return;
261  
262   // Set trigger bits, and mark this event for storage 
263   fAODFMD.SetTriggerBits(triggers);
264   MarkEventForStore();
265
266   if (found & AliFMDEventInspector::kNoSPD)     return;
267   if (found & AliFMDEventInspector::kNoFMD)     return;
268   if (found & AliFMDEventInspector::kNoVertex)  return;
269   fAODFMD.SetIpZ(vz);
270
271   if (found & AliFMDEventInspector::kBadVertex) return;
272
273   // We we do not want to use low flux specific code, we disable it here. 
274   if (!fEnableLowFlux) lowFlux = false;
275
276   // Get FMD data 
277   AliESDFMD* esdFMD = esd->GetFMDData();
278   // Apply the sharing filter (or hit merging or clustering if you like)
279   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
280     AliWarning("Sharing filter failed!");
281     return;
282   }
283
284   // Do the energy stuff 
285   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
286     AliWarning("Energy fitter failed");
287     return;
288   }
289
290   // Calculate the inclusive charged particle density 
291   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
292     AliWarning("Density calculator failed!");
293     return;
294   }
295   
296   // Do the secondary and other corrections. 
297   if (!fCorrections.Correct(fHistos, ivz)) { 
298     AliWarning("Corrections failed");
299     return;
300   }
301
302   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
303     AliWarning("Histogram collector failed");
304     return;
305   }
306
307   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
308     fHData->Add(&(fAODFMD.GetHistogram()));
309
310   PostData(1, fList);
311 }
312
313 //____________________________________________________________________
314 void
315 AliForwardMultiplicityTask::Terminate(Option_t*)
316 {
317   TList* list = dynamic_cast<TList*>(GetOutputData(1));
318   if (!list) {
319     AliError(Form("No output list defined (%p)", GetOutputData(1)));
320     if (GetOutputData(1)) GetOutputData(1)->Print();
321     return;
322   }
323   
324   // Get our histograms from the container 
325   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
326   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
327   TH1I* hTriggers    = 0;
328   if (!fEventInspector.FetchHistograms(list, hEventsTr, 
329                                        hEventsTrVtx, hTriggers)) { 
330     AliError(Form("Didn't get histograms from event selector "
331                   "(hEventsTr=%p,hEventsTrVtx=%p)", 
332                   hEventsTr, hEventsTrVtx));
333     list->ls();
334     return;
335   }
336
337   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
338   if (!hData) { 
339     AliError(Form("Couldn't get our summed histogram from output "
340                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
341     list->ls();
342     return;
343   }
344   
345   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
346   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
347   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
348   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
349   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
350   dNdeta->Divide(norm);
351   dNdeta->SetStats(0);
352   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
353                 "width");
354   list->Add(dNdeta);
355   list->Add(norm);
356
357   fEnergyFitter.Fit(list);
358   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
359   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
360   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
361 }
362
363 //____________________________________________________________________
364 void
365 AliForwardMultiplicityTask::MarkEventForStore() const
366 {
367   // Make sure the AOD tree is filled 
368   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
369   AliAODHandler*      ah = 
370     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
371   if (!ah)  
372     AliFatal("No AOD output handler set in analysis manager");
373
374   ah->SetFillAOD(kTRUE);
375 }
376
377 //____________________________________________________________________
378 void
379 AliForwardMultiplicityTask::Print(Option_t* option) const
380 {
381   std::cout << "AliForwardMultiplicityTask: " << GetName() << "\n" 
382             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no")
383             << std::endl;
384   gROOT->IncreaseDirLevel();
385   fEventInspector   .Print(option);
386   fEnergyFitter     .Print(option);    
387   fSharingFilter    .Print(option);
388   fDensityCalculator.Print(option);
389   fCorrections      .Print(option);
390   fHistCollector    .Print(option);
391   gROOT->DecreaseDirLevel();
392 }
393
394 //
395 // EOF
396 //