Added some more scripts
[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     fHData(0),
18     fFirstEvent(true),
19     fESDFMD(),
20     fHistos(),
21     fAODFMD(),
22     fEventInspector(),
23     fEnergyFitter(),
24     fSharingFilter(),
25     fDensityCalculator(),
26     fCorrections(),
27     fHistCollector(),
28     fList(0)
29 {
30 }
31
32 //____________________________________________________________________
33 AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
34   : AliAnalysisTaskSE(name), 
35     fHData(0),
36     fFirstEvent(true),
37     fESDFMD(),
38     fHistos(),
39     fAODFMD(kTRUE),
40     fEventInspector("event"),
41     fEnergyFitter("energy"),
42     fSharingFilter("sharing"), 
43     fDensityCalculator("density"),
44     fCorrections("corrections"),
45     fHistCollector("collector"),
46     fList(0)
47 {
48   DefineOutput(1, TList::Class());
49 }
50
51 //____________________________________________________________________
52 AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
53   : AliAnalysisTaskSE(o),
54     fHData(o.fHData),
55     fFirstEvent(true),
56     fESDFMD(o.fESDFMD),
57     fHistos(o.fHistos),
58     fAODFMD(o.fAODFMD),
59     fEventInspector(o.fEventInspector),
60     fEnergyFitter(o.fEnergyFitter),
61     fSharingFilter(o.fSharingFilter),
62     fDensityCalculator(o.fDensityCalculator),
63     fCorrections(o.fCorrections),
64     fHistCollector(o.fHistCollector),
65     fList(o.fList) 
66 {
67 }
68
69 //____________________________________________________________________
70 AliForwardMultiplicity&
71 AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
72 {
73   AliAnalysisTaskSE::operator=(o);
74
75   fHData             = o.fHData;
76   fFirstEvent        = o.fFirstEvent;
77   fEventInspector    = o.fEventInspector;
78   fEnergyFitter      = o.fEnergyFitter;
79   fSharingFilter     = o.fSharingFilter;
80   fDensityCalculator = o.fDensityCalculator;
81   fCorrections       = o.fCorrections;
82   fHistCollector     = o.fHistCollector;
83   fHistos            = o.fHistos;
84   fAODFMD            = o.fAODFMD;
85   fList              = o.fList;
86
87   return *this;
88 }
89
90 //____________________________________________________________________
91 void
92 AliForwardMultiplicity::SetDebug(Int_t dbg)
93 {
94   fEventInspector.SetDebug(dbg);
95   fEnergyFitter.SetDebug(dbg);
96   fSharingFilter.SetDebug(dbg);
97   fDensityCalculator.SetDebug(dbg);
98   fCorrections.SetDebug(dbg);
99   fHistCollector.SetDebug(dbg);
100 }
101
102 //____________________________________________________________________
103 void
104 AliForwardMultiplicity::Init()
105 {
106   fFirstEvent = true;
107 }
108
109 //____________________________________________________________________
110 void
111 AliForwardMultiplicity::InitializeSubs()
112 {
113   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
114   pars->Init(kTRUE);
115
116
117   TAxis e(pars->GetNetaBins(),  pars->GetEtaMin(),  pars->GetEtaMax());
118   TAxis v(pars->GetNvtxBins(), -pars->GetVtxCutZ(), pars->GetVtxCutZ());
119                         
120   fHistos.Init(e);
121   fAODFMD.Init(e);
122
123   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
124   fHData->SetStats(0);
125   fHData->SetDirectory(0);
126   fList->Add(fHData);
127
128   fEnergyFitter.Init(e);
129   fEventInspector.Init(v);
130   fHistCollector.Init(v);
131 }
132
133 //____________________________________________________________________
134 void
135 AliForwardMultiplicity::UserCreateOutputObjects()
136 {
137   fList = new TList;
138
139   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
140   AliAODHandler*      ah = dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
141   if (!ah) AliFatal("No AOD output handler set in analysis manager");
142     
143     
144   TObject* obj = &fAODFMD;
145   ah->AddBranch("AliAODForwardMult", &obj);
146
147   fEventInspector.DefineOutput(fList);
148   fEnergyFitter.DefineOutput(fList);
149   fSharingFilter.DefineOutput(fList);
150   fDensityCalculator.DefineOutput(fList);
151   fCorrections.DefineOutput(fList);
152 }
153 //____________________________________________________________________
154 void
155 AliForwardMultiplicity::UserExec(Option_t*)
156 {
157   // Get the input data 
158   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
159   if (!esd) { 
160     AliWarning("No ESD event found for input event");
161     return;
162   }
163
164   // On the first event, initialize the parameters 
165   if (fFirstEvent && esd->GetESDRun()) { 
166     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
167     AliInfo(Form("Initializing with parameters from the ESD:\n"
168                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
169                  "         AliESDEvent::GetBeamType()     ->%s\n"
170                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
171                  "         AliESDEvent::GetMagneticField()->%f\n"
172                  "         AliESDEvent::GetRunNumber()    ->%d\n",
173                  esd->GetBeamEnergy(), 
174                  esd->GetBeamType(),
175                  esd->GetCurrentL3(), 
176                  esd->GetMagneticField(),
177                  esd->GetRunNumber()));
178     pars->SetParametersFromESD(esd);
179     pars->PrintStatus();
180     fFirstEvent = false;
181
182     InitializeSubs();
183   }
184   // Clear stuff 
185   fHistos.Clear();
186   fESDFMD.Clear();
187   fAODFMD.Clear();
188
189   Bool_t   lowFlux  = kFALSE;
190   UInt_t   triggers = 0;
191   Int_t    ivz      = -1;
192   Double_t vz       = 0;
193   UInt_t   found    = fEventInspector.Process(esd, triggers, lowFlux, ivz, vz);
194   if (found & AliFMDEventInspector::kNoEvent)    return;
195   if (found & AliFMDEventInspector::kNoTriggers) return;
196  
197   // Set trigger bits, and mark this event for storage 
198   fAODFMD.SetTriggerBits(triggers);
199   MarkEventForStore();
200
201   if (found & AliFMDEventInspector::kNoSPD)     return;
202   if (found & AliFMDEventInspector::kNoFMD)     return;
203   if (found & AliFMDEventInspector::kNoVertex)  return;
204   fAODFMD.SetIpZ(vz);
205
206   if (found & AliFMDEventInspector::kBadVertex) return;
207
208   // Get FMD data 
209   AliESDFMD* esdFMD = esd->GetFMDData();
210   // Apply the sharing filter (or hit merging or clustering if you like)
211   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
212     AliWarning("Sharing filter failed!");
213     return;
214   }
215
216   // Do the energy stuff 
217   if (!fEnergyFitter.Accumulate(*esdFMD, triggers & AliAODForwardMult::kEmpty)){
218     AliWarning("Energy fitter failed");
219     return;
220   }
221
222   // Calculate the inclusive charged particle density 
223   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
224     AliWarning("Density calculator failed!");
225     return;
226   }
227   
228   // Do the secondary and other corrections. 
229   if (!fCorrections.Correct(fHistos, ivz)) { 
230     AliWarning("Corrections failed");
231     return;
232   }
233
234   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
235     AliWarning("Histogram collector failed");
236     return;
237   }
238
239   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
240     fHData->Add(&(fAODFMD.GetHistogram()));
241
242   PostData(1, fList);
243 }
244
245 //____________________________________________________________________
246 void
247 AliForwardMultiplicity::Terminate(Option_t*)
248 {
249   TList* list = dynamic_cast<TList*>(GetOutputData(1));
250   if (!list) {
251     AliError(Form("No output list defined (%p)", GetOutputData(1)));
252     if (GetOutputData(1)) GetOutputData(1)->Print();
253     return;
254   }
255   
256   // Get our histograms from the container 
257   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
258   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
259   TH1I* hTriggers    = 0;
260   if (!fEventInspector.FetchHistograms(list, hEventsTr, 
261                                        hEventsTrVtx, hTriggers)) { 
262     AliError(Form("Didn't get histograms from event selector "
263                   "(hEventsTr=%p,hEventsTrVtx=%p)", 
264                   hEventsTr, hEventsTrVtx));
265     list->ls();
266     return;
267   }
268
269   TH2D* hData        = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
270   if (!hData) { 
271     AliError(Form("Couldn't get our summed histogram from output "
272                   "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
273     list->ls();
274     return;
275   }
276   
277   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
278   TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
279   TH1D* norm   = hData->ProjectionX("dNdeta", 0, 1,  "");
280   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
281   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
282   dNdeta->Divide(norm);
283   dNdeta->SetStats(0);
284   dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
285                 "width");
286   list->Add(dNdeta);
287
288   fEnergyFitter.Fit(list);
289   fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
290   fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
291   fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
292 }
293
294 //____________________________________________________________________
295 void
296 AliForwardMultiplicity::MarkEventForStore() const
297 {
298   // Make sure the AOD tree is filled 
299   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
300   AliAODHandler*      ah = 
301     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
302   if (!ah)  
303     AliFatal("No AOD output handler set in analysis manager");
304
305   ah->SetFillAOD(kTRUE);
306 }
307
308 //____________________________________________________________________
309 void
310 AliForwardMultiplicity::Print(Option_t*) const
311 {
312 }
313
314 //
315 // EOF
316 //