ca934d6fb073bf50feabcc3693833303fce9ac03
[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     fRingSums(),
38     fEventInspector(),
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     fRingSums(),
58     fEventInspector("event"),
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     fRingSums(o.fRingSums),
82     fEventInspector(o.fEventInspector),
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   fSharingFilter     = o.fSharingFilter;
116   fDensityCalculator = o.fDensityCalculator;
117   fCorrections       = o.fCorrections;
118   fHistCollector     = o.fHistCollector;
119   fHistos            = o.fHistos;
120   fAODFMD            = o.fAODFMD;
121   fRingSums          = o.fRingSums;
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   fSharingFilter.SetDebug(dbg);
139   fDensityCalculator.SetDebug(dbg);
140   fCorrections.SetDebug(dbg);
141   fHistCollector.SetDebug(dbg);
142 }
143
144 //____________________________________________________________________
145 void
146 AliForwardMultiplicityTask::InitializeSubs()
147 {
148   // 
149   // Initialise the sub objects and stuff.  Called on first event 
150   // 
151   //
152   const TAxis* pe = 0;
153   const TAxis* pv = 0;
154
155   if (!ReadCorrections(pe,pv)) return;
156
157   fHistos.Init(*pe);
158   fAODFMD.Init(*pe);
159   fRingSums.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   TList* rings = new TList;
167   rings->SetName("ringSums");
168   rings->SetOwner();
169   fList->Add(rings);
170
171   rings->Add(fRingSums.Get(1, 'I'));
172   rings->Add(fRingSums.Get(2, 'I'));
173   rings->Add(fRingSums.Get(2, 'O'));
174   rings->Add(fRingSums.Get(3, 'I'));
175   rings->Add(fRingSums.Get(3, 'O'));
176   fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
177   fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
178   fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
179   fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
180   fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
181
182   fEventInspector.Init(*pv);
183   fSharingFilter.Init();
184   fDensityCalculator.Init(*pe);
185   fCorrections.Init(*pe);
186   fHistCollector.Init(*pv,*pe);
187
188   this->Print();
189 }
190
191 //____________________________________________________________________
192 void
193 AliForwardMultiplicityTask::UserCreateOutputObjects()
194 {
195   // 
196   // Create output objects 
197   // 
198   //
199   fList = new TList;
200   fList->SetOwner();
201
202   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
203   AliAODHandler*      ah = 
204     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
205   if (!ah) AliFatal("No AOD output handler set in analysis manager");
206     
207     
208   TObject* obj = &fAODFMD;
209   ah->AddBranch("AliAODForwardMult", &obj);
210
211   fEventInspector.DefineOutput(fList);
212   fSharingFilter.DefineOutput(fList);
213   fDensityCalculator.DefineOutput(fList);
214   fCorrections.DefineOutput(fList);
215   fHistCollector.DefineOutput(fList);
216
217   PostData(1, fList);
218 }
219 //____________________________________________________________________
220 void
221 AliForwardMultiplicityTask::UserExec(Option_t*)
222 {
223   // 
224   // Process each event 
225   // 
226   // Parameters:
227   //    option Not used
228   //  
229
230   // static Int_t cnt = 0;
231   // cnt++;
232   // Get the input data 
233   AliESDEvent* esd = GetESDEvent();
234
235   // Clear stuff 
236   fHistos.Clear();
237   fESDFMD.Clear();
238   fAODFMD.Clear();
239   
240   Bool_t   lowFlux   = kFALSE;
241   UInt_t   triggers  = 0;
242   UShort_t ivz       = 0;
243   Double_t vz        = 0;
244   Double_t cent      = -1;
245   UShort_t nClusters = 0;
246   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
247                                                ivz, vz, cent, nClusters);
248   
249   if (found & AliFMDEventInspector::kNoEvent)    return;
250   if (found & AliFMDEventInspector::kNoTriggers) return;
251
252   // Set trigger bits, and mark this event for storage 
253   fAODFMD.SetTriggerBits(triggers);
254   fAODFMD.SetSNN(fEventInspector.GetEnergy());
255   fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
256   fAODFMD.SetCentrality(cent);
257   fAODFMD.SetNClusters(nClusters);
258   MarkEventForStore();
259   
260   if (found & AliFMDEventInspector::kNoSPD)      return;
261   if (found & AliFMDEventInspector::kNoFMD)      return;
262   if (found & AliFMDEventInspector::kNoVertex)   return;
263   
264   if (triggers & AliAODForwardMult::kPileUp) return;
265   
266   fAODFMD.SetIpZ(vz);
267
268   if (found & AliFMDEventInspector::kBadVertex) return;
269
270   // We we do not want to use low flux specific code, we disable it here. 
271   if (!fEnableLowFlux) lowFlux = false;
272
273   // Get FMD data 
274   AliESDFMD* esdFMD = esd->GetFMDData();
275   //  // Apply the sharing filter (or hit merging or clustering if you like)
276   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
277     AliWarning("Sharing filter failed!");
278     return;
279   }
280
281   // Calculate the inclusive charged particle density 
282   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent)) { 
283     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { 
284     AliWarning("Density calculator failed!");
285     return;
286   }
287   
288   // Do the secondary and other corrections. 
289   if (!fCorrections.Correct(fHistos, ivz)) { 
290     AliWarning("Corrections failed");
291     return;
292   }
293
294   if (!fHistCollector.Collect(fHistos, fRingSums, 
295                               ivz, fAODFMD.GetHistogram())) {
296     AliWarning("Histogram collector failed");
297     return;
298   }
299
300   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
301     fHData->Add(&(fAODFMD.GetHistogram()));
302
303   PostData(1, fList);
304 }
305
306 //____________________________________________________________________
307 void
308 AliForwardMultiplicityTask::Terminate(Option_t*)
309 {
310   // 
311   // End of job
312   // 
313   // Parameters:
314   //    option Not used 
315   //
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,  0,  "");
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   MakeRingdNdeta(list, "ringSums", list, "ringResults");
358
359   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
360   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
361   fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
362 }
363
364 //
365 // EOF
366 //