Various fixes to deal with centrality
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMCMultiplicityTask.cxx
1 // 
2 // Calculate the multiplicity in the forward regions event-by-event 
3 // 
4 // Inputs: 
5 //   - AliESDEvent 
6 //   - Kinematics
7 //   - Track references
8 //
9 // Outputs: 
10 //   - AliAODForwardMult 
11 // 
12 // Histograms 
13 //   
14 // Corrections used 
15 // 
16 #include "AliForwardMCMultiplicityTask.h"
17 #include "AliTriggerAnalysis.h"
18 #include "AliPhysicsSelection.h"
19 #include "AliLog.h"
20 #include "AliESDEvent.h"
21 #include "AliAODHandler.h"
22 #include "AliMultiplicity.h"
23 #include "AliInputEventHandler.h"
24 #include "AliForwardCorrectionManager.h"
25 #include "AliAnalysisManager.h"
26 #include <TH1.h>
27 #include <TDirectory.h>
28 #include <TTree.h>
29 #include <TROOT.h>
30
31 //====================================================================
32 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
33   : AliForwardMultiplicityBase(),
34     fHData(0),
35     fESDFMD(),
36     fHistos(),
37     fAODFMD(),
38     fMCESDFMD(),
39     fMCHistos(),
40     fMCAODFMD(),
41     fPrimary(0),
42     fEventInspector(),
43     fEnergyFitter(),
44     fSharingFilter(),
45     fDensityCalculator(),
46     fCorrections(),
47     fHistCollector(),
48     fList(0)
49 {
50   // 
51   // Constructor
52   //
53 }
54
55 //____________________________________________________________________
56 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
57   : AliForwardMultiplicityBase(name), 
58     fHData(0),
59     fESDFMD(),
60     fHistos(),
61     fAODFMD(kFALSE),
62     fMCESDFMD(),
63     fMCHistos(),
64     fMCAODFMD(kTRUE),
65     fPrimary(0),
66     fEventInspector("event"),
67     fEnergyFitter("energy"),
68     fSharingFilter("sharing"), 
69     fDensityCalculator("density"),
70     fCorrections("corrections"),
71     fHistCollector("collector"),
72     fList(0)
73 {
74   // 
75   // Constructor 
76   // 
77   // Parameters:
78   //    name Name of task 
79   //
80   DefineOutput(1, TList::Class());
81 }
82
83 //____________________________________________________________________
84 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
85   : AliForwardMultiplicityBase(o),
86     fHData(o.fHData),
87     fESDFMD(o.fESDFMD),
88     fHistos(o.fHistos),
89     fAODFMD(o.fAODFMD),
90     fMCESDFMD(o.fMCESDFMD),
91     fMCHistos(o.fMCHistos),
92     fMCAODFMD(o.fMCAODFMD),
93     fPrimary(o.fPrimary),
94     fEventInspector(o.fEventInspector),
95     fEnergyFitter(o.fEnergyFitter),
96     fSharingFilter(o.fSharingFilter),
97     fDensityCalculator(o.fDensityCalculator),
98     fCorrections(o.fCorrections),
99     fHistCollector(o.fHistCollector),
100     fList(o.fList) 
101 {
102   // 
103   // Copy constructor 
104   // 
105   // Parameters:
106   //    o Object to copy from 
107   //
108   DefineOutput(1, TList::Class());
109 }
110
111 //____________________________________________________________________
112 AliForwardMCMultiplicityTask&
113 AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
114 {
115   // 
116   // Assignment operator 
117   // 
118   // Parameters:
119   //    o Object to assign from 
120   // 
121   // Return:
122   //    Reference to this object 
123   //
124   AliForwardMultiplicityBase::operator=(o);
125
126   fHData             = o.fHData;
127   fEventInspector    = o.fEventInspector;
128   fEnergyFitter      = o.fEnergyFitter;
129   fSharingFilter     = o.fSharingFilter;
130   fDensityCalculator = o.fDensityCalculator;
131   fCorrections       = o.fCorrections;
132   fHistCollector     = o.fHistCollector;
133   fHistos            = o.fHistos;
134   fAODFMD            = o.fAODFMD;
135   fMCHistos          = o.fMCHistos;
136   fMCAODFMD          = o.fMCAODFMD;
137   fPrimary           = o.fPrimary;
138   fList              = o.fList;
139
140   return *this;
141 }
142
143 //____________________________________________________________________
144 void
145 AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
146 {
147   // 
148   // Set debug level 
149   // 
150   // Parameters:
151   //    dbg debug level
152   //
153   fEventInspector.SetDebug(dbg);
154   fEnergyFitter.SetDebug(dbg);
155   fSharingFilter.SetDebug(dbg);
156   fDensityCalculator.SetDebug(dbg);
157   fCorrections.SetDebug(dbg);
158   fHistCollector.SetDebug(dbg);
159 }
160 //____________________________________________________________________
161 void
162 AliForwardMCMultiplicityTask::InitializeSubs()
163 {
164   // 
165   // Initialise the sub objects and stuff.  Called on first event 
166   // 
167   //
168   const TAxis* pe = 0;
169   const TAxis* pv = 0;
170
171   if (!ReadCorrections(pe,pv,true)) return;
172
173   fHistos.Init(*pe);
174   fAODFMD.Init(*pe);
175   fMCHistos.Init(*pe);
176   fMCAODFMD.Init(*pe);
177
178   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
179   fHData->SetStats(0);
180   fHData->SetDirectory(0);
181
182   fList->Add(fHData);
183
184   fEnergyFitter.Init(*pe);
185   fEventInspector.Init(*pv);
186   fDensityCalculator.Init(*pe);
187   fCorrections.Init(*pe);
188   fHistCollector.Init(*pv,*pe);
189
190   this->Print();
191 }
192
193 //____________________________________________________________________
194 void
195 AliForwardMCMultiplicityTask::UserCreateOutputObjects()
196 {
197   // 
198   // Create output objects 
199   // 
200   //
201   fList = new TList;
202   fList->SetOwner();
203
204   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
205   AliAODHandler*      ah = 
206     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
207   if (!ah) AliFatal("No AOD output handler set in analysis manager");
208     
209     
210   TObject* obj = &fAODFMD;
211   ah->AddBranch("AliAODForwardMult", &obj);
212
213   TObject* mcobj = &fMCAODFMD;
214   ah->AddBranch("AliAODForwardMult", &mcobj);
215
216   fPrimary = new TH2D("primary", "MC Primaries", 
217                       200, -4, 6, 20, 0, 2*TMath::Pi());
218   fPrimary->SetXTitle("#eta");
219   fPrimary->SetYTitle("#varphi [radians]");
220   fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
221   fPrimary->Sumw2();
222   fPrimary->SetStats(0);
223   fPrimary->SetDirectory(0);
224   ah->AddBranch("TH2D", &fPrimary);
225
226   fEventInspector.DefineOutput(fList);
227   fEnergyFitter.DefineOutput(fList);
228   fSharingFilter.DefineOutput(fList);
229   fDensityCalculator.DefineOutput(fList);
230   fCorrections.DefineOutput(fList);
231   fHistCollector.DefineOutput(fList);
232
233   PostData(1, fList);
234 }
235 //____________________________________________________________________
236 void
237 AliForwardMCMultiplicityTask::UserExec(Option_t*)
238 {
239   // 
240   // Process each event 
241   // 
242   // Parameters:
243   //    option Not used
244   //  
245
246   // Get the input data 
247   AliESDEvent* esd     = GetESDEvent();
248   AliMCEvent*  mcEvent = MCEvent();
249
250   // Clear stuff 
251   fHistos.Clear();
252   fESDFMD.Clear();
253   fAODFMD.Clear();
254   fMCHistos.Clear();
255   fMCESDFMD.Clear();
256   fMCAODFMD.Clear();
257   fPrimary->Reset();
258
259   Bool_t   lowFlux  = kFALSE;
260   UInt_t   triggers = 0;
261   UShort_t ivz      = 0;
262   Double_t vz       = 0;
263   Double_t cent     = 0;
264   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, 
265                                               ivz, vz, cent);
266   UShort_t ivzMC    = 0;
267   Double_t vzMC     = 0;
268   Double_t phiR     = 0;
269   Double_t b        = 0;
270   Int_t    npart    = 0;
271   Int_t    nbin     = 0;
272   // UInt_t   foundMC  = 
273   fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, 
274                             npart, nbin, phiR);
275   fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
276   
277   //Store all events
278   MarkEventForStore();
279   
280   Bool_t isAccepted = true;
281   if (found & AliFMDEventInspector::kNoEvent)    isAccepted = false; // return;
282   if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
283   //MarkEventForStore();
284   // Set trigger bits, and mark this event for storage 
285   fAODFMD.SetTriggerBits(triggers);
286   fAODFMD.SetSNN(fEventInspector.GetEnergy());
287   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
288   fAODFMD.SetCentrality(cent);
289
290   fMCAODFMD.SetTriggerBits(triggers);
291   fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
292   fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
293   fMCAODFMD.SetCentrality(cent);
294   
295   //All events should be stored - HHD
296   //if (isAccepted) MarkEventForStore();
297
298   // Disable this check on SPD - will bias data 
299   // if (found & AliFMDEventInspector::kNoSPD)  isAccepted = false; // return;
300   if (found & AliFMDEventInspector::kNoFMD)     isAccepted = false; // return;
301   if (found & AliFMDEventInspector::kNoVertex)  isAccepted = false; // return;
302
303   if (isAccepted) {
304     fAODFMD.SetIpZ(vz);
305     fMCAODFMD.SetIpZ(vz);
306   }
307   if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
308
309   // We we do not want to use low flux specific code, we disable it here. 
310   if (!fEnableLowFlux) lowFlux = false;
311   
312   
313
314   // Get FMD data 
315   AliESDFMD*  esdFMD  = esd->GetFMDData();
316
317   // Apply the sharing filter (or hit merging or clustering if you like)
318   if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
319     AliWarning("Sharing filter failed!");
320     return;
321   }
322   if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) { 
323     AliWarning("MC Sharing filter failed!");
324     return;
325   }
326   if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
327   // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
328   
329   //MarkEventForStore();
330   fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
331
332   // Do the energy stuff 
333   if (!fEnergyFitter.Accumulate(*esdFMD, cent, 
334                                 triggers & AliAODForwardMult::kEmpty)){
335     AliWarning("Energy fitter failed");
336     return;
337   }
338
339   // Calculate the inclusive charged particle density 
340   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
341     AliWarning("Density calculator failed!");
342     return;
343   }
344   if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) { 
345     AliWarning("MC Density calculator failed!");
346     return;
347   }
348   fDensityCalculator.CompareResults(fHistos, fMCHistos);
349   
350   // Do the secondary and other corrections. 
351   if (!fCorrections.Correct(fHistos, ivz)) { 
352     AliWarning("Corrections failed");
353     return;
354   }
355   if (!fCorrections.CorrectMC(fMCHistos, ivz)) { 
356     AliWarning("MC Corrections failed");
357     return;
358   }
359   fCorrections.CompareResults(fHistos, fMCHistos);
360     
361   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
362     AliWarning("Histogram collector failed");
363     return;
364   }
365   if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
366     AliWarning("MC Histogram collector failed");
367     return;
368   }
369
370   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
371     fHData->Add(&(fAODFMD.GetHistogram()));
372
373   PostData(1, fList);
374 }
375
376 //____________________________________________________________________
377 void
378 AliForwardMCMultiplicityTask::Terminate(Option_t*)
379 {
380   // 
381   // End of job
382   // 
383   // Parameters:
384   //    option Not used 
385   //
386   TList* list = dynamic_cast<TList*>(GetOutputData(1));
387   if (!list) {
388     AliError(Form("No output list defined (%p)", GetOutputData(1)));
389     if (GetOutputData(1)) GetOutputData(1)->Print();
390     return;
391   }
392   
393   // Get our histograms from the container 
394   TH1I* hEventsTr    = 0;
395   TH1I* hEventsTrVtx = 0;
396   TH1I* hTriggers    = 0;
397   if (!fEventInspector.FetchHistograms(list, hEventsTr, 
398                                        hEventsTrVtx, hTriggers)) { 
399     AliError(Form("Didn't get histograms from event selector "
400                   "(hEventsTr=%p,hEventsTrVtx=%p)", 
401                   hEventsTr, hEventsTrVtx));
402     list->ls();
403     return;
404   }
405
406   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
407   if (!hData) { 
408     AliError(Form("Couldn't get our summed histogram from output "
409                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
410     list->ls();
411     return;
412   }
413   
414   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
415   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
416   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
417   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
418   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
419   dNdeta->Divide(norm);
420   dNdeta->SetStats(0);
421   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
422                 "width");
423   list->Add(dNdeta);
424   list->Add(norm);
425
426   fEnergyFitter.Fit(list);
427   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
428   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
429   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
430 }
431
432
433 //
434 // EOF
435 //