Make sure that histograms are obtained from output list in
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicity.cxx
1 #include "AliForwardMultiplicity.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 <TH1.h>
11 #include <TDirectory.h>
12 #include <TTree.h>
13
14 //====================================================================
15 AliForwardMultiplicity::AliForwardMultiplicity()
16   : AliAnalysisTaskSE(),
17     fHEventsTr(0), 
18     fHEventsTrVtx(0),
19     fHTriggers(0),
20     fHData(0),
21     fFirstEvent(true),
22     fLowFluxCut(1000),
23     fESDFMD(),
24     fHistos(),
25     fAODFMD(),
26     fSharingFilter(),
27     fDensityCalculator(),
28     fCorrections(),
29     fHistCollector(),
30     fList(0), 
31     fTree(0)
32 {
33 }
34
35 //____________________________________________________________________
36 AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
37   : AliAnalysisTaskSE(name), 
38     fHEventsTr(0), 
39     fHEventsTrVtx(0), 
40     fHTriggers(0),
41     fHData(0),
42     fFirstEvent(true),
43     fLowFluxCut(1000),
44     fESDFMD(),
45     fHistos(),
46     fAODFMD(kTRUE),
47     fSharingFilter("sharing"), 
48     fDensityCalculator("density"),
49     fCorrections("corrections"),
50     fHistCollector("collector"),
51     fList(0), 
52     fTree(0)
53 {
54   DefineOutput(1, TList::Class());
55   // DefineOutput(2, TTree::Class());
56 }
57
58 //____________________________________________________________________
59 AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
60   : AliAnalysisTaskSE(o),
61     fHEventsTr(o.fHEventsTr), 
62     fHEventsTrVtx(o.fHEventsTrVtx), 
63     fHTriggers(o.fHTriggers),
64     fHData(o.fHData),
65     fFirstEvent(true),
66     fLowFluxCut(1000),
67     fESDFMD(o.fESDFMD),
68     fHistos(o.fHistos),
69     fAODFMD(o.fAODFMD),
70     fSharingFilter(o.fSharingFilter),
71     fDensityCalculator(o.fDensityCalculator),
72     fCorrections(o.fCorrections),
73     fHistCollector(o.fHistCollector),
74     fList(o.fList), 
75     fTree(o.fTree)
76 {
77 }
78
79 //____________________________________________________________________
80 AliForwardMultiplicity&
81 AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
82 {
83   fHEventsTr         = o.fHEventsTr;
84   fHEventsTrVtx      = o.fHEventsTrVtx;
85   fHTriggers         = o.fHTriggers;
86   fHData             = o.fHData;
87   fFirstEvent        = o.fFirstEvent;
88   fSharingFilter     = o.fSharingFilter;
89   fDensityCalculator = o.fDensityCalculator;
90   fCorrections       = o.fCorrections;
91   fHistCollector     = o.fHistCollector;
92   fHistos            = o.fHistos;
93   fAODFMD            = o.fAODFMD;
94   fList              = o.fList;
95   fTree              = o.fTree;
96
97   return *this;
98 }
99
100 //____________________________________________________________________
101 void
102 AliForwardMultiplicity::Init()
103 {
104   fFirstEvent = true;
105 }
106
107 //____________________________________________________________________
108 void
109 AliForwardMultiplicity::InitializeSubs()
110 {
111   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
112   pars->Init(kTRUE);
113
114   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
115                       pars->GetNvtxBins(), 
116                       -pars->GetVtxCutZ(), 
117                       pars->GetVtxCutZ());
118   fHEventsTr->SetXTitle("v_{z} [cm]");
119   fHEventsTr->SetYTitle("# of events");
120   fHEventsTr->SetFillColor(kRed+1);
121   fHEventsTr->SetFillStyle(3001);
122   fHEventsTr->SetDirectory(0);
123   fList->Add(fHEventsTr);
124   // fHEventsTr->Sumw2();
125
126   fHEventsTrVtx = new TH1I("nEventsTrVtx", 
127                            "Number of events w/trigger and vertex", 
128                            pars->GetNvtxBins(), 
129                            -pars->GetVtxCutZ(), 
130                            pars->GetVtxCutZ());
131   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
132   fHEventsTrVtx->SetYTitle("# of events");
133   fHEventsTrVtx->SetFillColor(kBlue+1);
134   fHEventsTrVtx->SetFillStyle(3001);
135   fHEventsTrVtx->SetDirectory(0);
136   fList->Add(fHEventsTrVtx);
137   // fHEventsTrVtx->Sumw2();
138
139       
140   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
141   fHTriggers->SetFillColor(kRed+1);
142   fHTriggers->SetFillStyle(3001);
143   fHTriggers->SetStats(0);
144   fHTriggers->SetDirectory(0);
145   fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
146   fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
147   fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
148   fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
149   fHTriggers->GetXaxis()->SetBinLabel(5,"A");
150   fHTriggers->GetXaxis()->SetBinLabel(6,"B");
151   fHTriggers->GetXaxis()->SetBinLabel(7,"C");
152   fHTriggers->GetXaxis()->SetBinLabel(8,"E");
153   fList->Add(fHTriggers);
154
155   TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
156   fHistos.Init(e);
157   fAODFMD.Init(e);
158
159   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
160   fHData->SetStats(0);
161   fHData->SetDirectory(0);
162   fList->Add(fHData);
163
164   fHistCollector.Init(*(fHEventsTr->GetXaxis()));
165 }
166
167 //____________________________________________________________________
168 void
169 AliForwardMultiplicity::UserCreateOutputObjects()
170 {
171   fList = new TList;
172
173   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
174   AliAODHandler*      ah = 
175     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
176   if (!ah)  
177     AliFatal("No AOD output handler set in analysis manager");
178     
179     
180   TObject* obj = &fAODFMD;
181   ah->AddBranch("AliAODForwardMult", &obj);
182
183   fSharingFilter.DefineOutput(fList);
184   fDensityCalculator.DefineOutput(fList);
185   fCorrections.DefineOutput(fList);
186
187   // fTree = new TTree("T", "T");
188   // fTree->Branch("forward", &fAODFMD);
189
190   PostData(1, fList);
191   // PostData(2, fTree);
192 }
193 //____________________________________________________________________
194 void
195 AliForwardMultiplicity::UserExec(Option_t*)
196 {
197   // Get the input data 
198   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
199   if (!esd) { 
200     AliWarning("No ESD event found for input event");
201     return;
202   }
203
204 #if 0
205   static Int_t nEvents = 0;
206   nEvents++;
207   if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
208 #endif
209
210   // On the first event, initialize the parameters 
211   if (fFirstEvent) { 
212     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
213     pars->SetParametersFromESD(esd);
214     pars->PrintStatus();
215     fFirstEvent = false;
216
217     InitializeSubs();
218   }
219   // Clear stuff 
220   fHistos.Clear();
221   fESDFMD.Clear();
222   fAODFMD.Clear();
223
224   // Read trigger information from the ESD and store in AOD object
225   UInt_t triggers = 0;
226   if (!AliForwardUtil::ReadTriggers(esd, triggers, fHTriggers)) { 
227 #ifdef VERBOSE
228     AliWarning("Failed to read triggers from ESD");
229 #endif
230     return;
231   }
232   fAODFMD.SetTriggerBits(triggers);
233
234   // Mark this event for storage 
235   MarkEventForStore();
236
237   // Check if this is a high-flux event 
238   const AliMultiplicity* testmult = esd->GetMultiplicity();
239   if (!testmult) { 
240 #ifdef VERBOSE
241     AliWarning("No central multiplicity object found");
242 #endif
243     return;
244   }
245   Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
246
247   // Get the FMD ESD data 
248   AliESDFMD* esdFMD = esd->GetFMDData();
249   if (!esdFMD) { 
250 #ifdef VERBOSE
251     AliWarning("No FMD data found in ESD");
252 #endif
253     return;
254   }
255
256   // Get the vertex information 
257   Double_t vz   = 0;
258   Bool_t   vzOk = AliForwardUtil::ReadVertex(esd, vz);
259
260   fHEventsTr->Fill(vz);
261   if (!vzOk) { 
262 #ifdef VERBOSE
263     AliWarning("Failed to read vertex from ESD");
264 #endif
265     return;
266   }
267   fHEventsTrVtx->Fill(vz);
268
269   // Get the vertex bin 
270   Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
271   fAODFMD.SetIpZ(vz);
272   if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) { 
273 #if 0
274     AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
275                     vz, fHEventsTr->GetXaxis()->GetXmin(), 
276                     fHEventsTr->GetXaxis()->GetXmax()));
277 #endif
278     return;
279   }
280
281   // Apply the sharing filter (or hit merging or clustering if you like)
282   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { 
283 #ifdef VERBOSE
284     AliWarning("Sharing filter failed!");
285 #endif
286     return;
287   }
288
289   // Calculate the inclusive charged particle density 
290   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
291     AliWarning("Density calculator failed!");
292     return;
293   }
294   
295   // Do the secondary and other corrections. 
296   if (!fCorrections.Correct(fHistos, ivz)) { 
297     AliWarning("Corrections failed");
298     return;
299   }
300
301   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
302     AliWarning("Histogram collector failed");
303     return;
304   }
305
306   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
307     fHData->Add(&(fAODFMD.GetHistogram()));
308 }
309
310 //____________________________________________________________________
311 void
312 AliForwardMultiplicity::Terminate(Option_t*)
313 {
314   TList* list = dynamic_cast<TList*>(GetOutputData(1));
315   if (!list) {
316     AliError("No output list defined");
317     return;
318   }
319   
320   // Get our histograms from the container 
321   TH1I* hEventsTr    = static_cast<TH1I*>(list->FindObject("nEventsTr"));
322   TH1I* hEventsTrVtx = static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
323   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
324   
325   
326   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
327   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
328   TH1D* norm   = hData->ProjectionX("dNdeta", 0, 1,  "");
329   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
330   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
331   dNdeta->Divide(norm);
332   dNdeta->SetStats(0);
333   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
334                 "width");
335   list->Add(dNdeta);
336   
337   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
338   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
339   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
340 }
341
342 //____________________________________________________________________
343 void
344 AliForwardMultiplicity::MarkEventForStore() const
345 {
346   // Make sure the AOD tree is filled 
347   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
348   AliAODHandler*      ah = 
349     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
350   if (!ah)  
351     AliFatal("No AOD output handler set in analysis manager");
352
353   ah->SetFillAOD(kTRUE);
354 }
355
356 //____________________________________________________________________
357 void
358 AliForwardMultiplicity::Print(Option_t*) const
359 {
360 }
361
362 //
363 // EOF
364 //