Added class AliForwarddNdetaTask to do the dN/deta
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicityTask.cxx
1 // 
2 // Calculate the multiplicity in the forward regions event-by-event 
3 // 
4 // Inputs: 
5 //   - AliESDEvent 
6 //
7 // Outputs: 
8 //   - AliAODForwardMult 
9 // 
10 // Histograms 
11 //   
12 // Corrections used 
13 //
14 #include "AliForwardMultiplicityTask.h"
15 #include "AliTriggerAnalysis.h"
16 #include "AliPhysicsSelection.h"
17 #include "AliLog.h"
18 #include "AliESDEvent.h"
19 #include "AliAODHandler.h"
20 #include "AliMultiplicity.h"
21 #include "AliInputEventHandler.h"
22 #include "AliForwardCorrectionManager.h"
23 #include "AliAnalysisManager.h"
24 #include <TH1.h>
25 #include <TDirectory.h>
26 #include <TTree.h>
27 #include <TROOT.h>
28
29
30 //====================================================================
31 AliForwardMultiplicityTask::AliForwardMultiplicityTask()
32   : AliForwardMultiplicityBase(),
33     fHData(0),
34     fESDFMD(),
35     fHistos(),
36     fAODFMD(),
37     fEventInspector(),
38     fEnergyFitter(),
39     fSharingFilter(),
40     fDensityCalculator(),
41     fCorrections(),
42     fHistCollector(),
43     fList(0)
44 {
45   // 
46   // Constructor
47   //
48 }
49
50 //____________________________________________________________________
51 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const char* name)
52   : AliForwardMultiplicityBase(name),
53     fHData(0),
54     fESDFMD(),
55     fHistos(),
56     fAODFMD(false),
57     fEventInspector("event"),
58     fEnergyFitter("energy"),
59     fSharingFilter("sharing"), 
60     fDensityCalculator("density"),
61     fCorrections("corrections"),
62     fHistCollector("collector"),
63     fList(0)
64 {
65   // 
66   // Constructor 
67   // 
68   // Parameters:
69   //    name Name of task 
70   //
71   DefineOutput(1, TList::Class());
72 }
73
74 //____________________________________________________________________
75 AliForwardMultiplicityTask::AliForwardMultiplicityTask(const AliForwardMultiplicityTask& o)
76   : AliForwardMultiplicityBase(o),
77     fHData(o.fHData),
78     fESDFMD(o.fESDFMD),
79     fHistos(o.fHistos),
80     fAODFMD(o.fAODFMD),
81     fEventInspector(o.fEventInspector),
82     fEnergyFitter(o.fEnergyFitter),
83     fSharingFilter(o.fSharingFilter),
84     fDensityCalculator(o.fDensityCalculator),
85     fCorrections(o.fCorrections),
86     fHistCollector(o.fHistCollector),
87     fList(o.fList) 
88 {
89   // 
90   // Copy constructor 
91   // 
92   // Parameters:
93   //    o Object to copy from 
94   //
95   DefineOutput(1, TList::Class());
96 }
97
98 //____________________________________________________________________
99 AliForwardMultiplicityTask&
100 AliForwardMultiplicityTask::operator=(const AliForwardMultiplicityTask& o)
101 {
102   // 
103   // Assignment operator 
104   // 
105   // Parameters:
106   //    o Object to assign from 
107   // 
108   // Return:
109   //    Reference to this object 
110   //
111   AliForwardMultiplicityBase::operator=(o);
112
113   fHData             = o.fHData;
114   fEventInspector    = o.fEventInspector;
115   fEnergyFitter      = o.fEnergyFitter;
116   fSharingFilter     = o.fSharingFilter;
117   fDensityCalculator = o.fDensityCalculator;
118   fCorrections       = o.fCorrections;
119   fHistCollector     = o.fHistCollector;
120   fHistos            = o.fHistos;
121   fAODFMD            = o.fAODFMD;
122   fList              = o.fList;
123
124   return *this;
125 }
126
127 //____________________________________________________________________
128 void
129 AliForwardMultiplicityTask::SetDebug(Int_t dbg)
130 {
131   // 
132   // Set debug level 
133   // 
134   // Parameters:
135   //    dbg Debug level
136   //
137   fEventInspector.SetDebug(dbg);
138   fEnergyFitter.SetDebug(dbg);
139   fSharingFilter.SetDebug(dbg);
140   fDensityCalculator.SetDebug(dbg);
141   fCorrections.SetDebug(dbg);
142   fHistCollector.SetDebug(dbg);
143 }
144
145 //____________________________________________________________________
146 void
147 AliForwardMultiplicityTask::InitializeSubs()
148 {
149   // 
150   // Initialise the sub objects and stuff.  Called on first event 
151   // 
152   //
153   const TAxis* pe = 0;
154   const TAxis* pv = 0;
155
156   if (!ReadCorrections(pe,pv)) return;
157
158   fHistos.Init(*pe);
159   fAODFMD.Init(*pe);
160
161   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
162   fHData->SetStats(0);
163   fHData->SetDirectory(0);
164   fList->Add(fHData);
165
166   fEnergyFitter.Init(*pe);
167   fEventInspector.Init(*pv);
168   fDensityCalculator.Init(*pe);
169   fCorrections.Init(*pe);
170   fHistCollector.Init(*pv);
171
172   this->Print();
173 }
174
175 //____________________________________________________________________
176 void
177 AliForwardMultiplicityTask::UserCreateOutputObjects()
178 {
179   // 
180   // Create output objects 
181   // 
182   //
183   fList = new TList;
184
185   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
186   AliAODHandler*      ah = 
187     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
188   if (!ah) AliFatal("No AOD output handler set in analysis manager");
189     
190     
191   TObject* obj = &fAODFMD;
192   ah->AddBranch("AliAODForwardMult", &obj);
193
194   fEventInspector.DefineOutput(fList);
195   fEnergyFitter.DefineOutput(fList);
196   fSharingFilter.DefineOutput(fList);
197   fDensityCalculator.DefineOutput(fList);
198   fCorrections.DefineOutput(fList);
199
200   PostData(1, fList);
201 }
202 //____________________________________________________________________
203 void
204 AliForwardMultiplicityTask::UserExec(Option_t*)
205 {
206   // 
207   // Process each event 
208   // 
209   // Parameters:
210   //    option Not used
211   //  
212
213   // static Int_t cnt = 0;
214   // cnt++;
215   // Get the input data 
216   AliESDEvent* esd = GetESDEvent();
217
218   // Clear stuff 
219   fHistos.Clear();
220   fESDFMD.Clear();
221   fAODFMD.Clear();
222
223   Bool_t   lowFlux  = kFALSE;
224   UInt_t   triggers = 0;
225   UShort_t ivz      = 0;
226   Double_t vz       = 0;
227   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
228  
229   if (found & AliFMDEventInspector::kNoEvent)    return;
230   if (found & AliFMDEventInspector::kNoTriggers) return;
231
232   // Set trigger bits, and mark this event for storage 
233   fAODFMD.SetTriggerBits(triggers);
234   fAODFMD.SetSNN(fEventInspector.GetEnergy());
235   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
236   MarkEventForStore();
237
238   if (found & AliFMDEventInspector::kNoSPD)      return;
239   if (found & AliFMDEventInspector::kNoFMD)      return;
240   if (found & AliFMDEventInspector::kNoVertex)   return;
241   fAODFMD.SetIpZ(vz);
242
243   if (found & AliFMDEventInspector::kBadVertex) return;
244
245   // We we do not want to use low flux specific code, we disable it here. 
246   if (!fEnableLowFlux) lowFlux = false;
247
248   // Get FMD data 
249   AliESDFMD* esdFMD = esd->GetFMDData();
250   // Apply the sharing filter (or hit merging or clustering if you like)
251   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
252     AliWarning("Sharing filter failed!");
253     return;
254   }
255
256   // Do the energy stuff 
257   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
258     AliWarning("Energy fitter failed");
259     return;
260   }
261
262   // Calculate the inclusive charged particle density 
263   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
264     AliWarning("Density calculator failed!");
265     return;
266   }
267   
268   // Do the secondary and other corrections. 
269   if (!fCorrections.Correct(fHistos, ivz)) { 
270     AliWarning("Corrections failed");
271     return;
272   }
273
274   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
275     AliWarning("Histogram collector failed");
276     return;
277   }
278
279   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
280     fHData->Add(&(fAODFMD.GetHistogram()));
281
282   PostData(1, fList);
283 }
284
285 //____________________________________________________________________
286 void
287 AliForwardMultiplicityTask::Terminate(Option_t*)
288 {
289   // 
290   // End of job
291   // 
292   // Parameters:
293   //    option Not used 
294   //
295
296   TList* list = dynamic_cast<TList*>(GetOutputData(1));
297   if (!list) {
298     AliError(Form("No output list defined (%p)", GetOutputData(1)));
299     if (GetOutputData(1)) GetOutputData(1)->Print();
300     return;
301   }
302   
303   // Get our histograms from the container 
304   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
305   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
306   TH1I* hTriggers    = 0;
307   if (!fEventInspector.FetchHistograms(list, hEventsTr, 
308                                        hEventsTrVtx, hTriggers)) { 
309     AliError(Form("Didn't get histograms from event selector "
310                   "(hEventsTr=%p,hEventsTrVtx=%p)", 
311                   hEventsTr, hEventsTrVtx));
312     list->ls();
313     return;
314   }
315
316   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
317   if (!hData) { 
318     AliError(Form("Couldn't get our summed histogram from output "
319                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
320     list->ls();
321     return;
322   }
323   
324   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
325   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
326   TH1D* norm   = hData->ProjectionX("norm",   0,  1,  "");
327   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
328   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
329   dNdeta->Divide(norm);
330   dNdeta->SetStats(0);
331   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
332                 "width");
333   list->Add(dNdeta);
334   list->Add(norm);
335
336   fEnergyFitter.Fit(list);
337   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
338   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
339   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
340 }
341
342 //
343 // EOF
344 //