f6393b0ee6e538e5c0518061eecf2b2f44465e79
[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("nEvents", "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   // fHEventsTr->Sumw2();
124
125   fHEventsTrVtx = new TH1I("nEventsTrVtx", 
126                            "Number of events w/trigger and vertex", 
127                            pars->GetNvtxBins(), 
128                            -pars->GetVtxCutZ(), 
129                            pars->GetVtxCutZ());
130   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
131   fHEventsTrVtx->SetYTitle("# of events");
132   fHEventsTrVtx->SetFillColor(kBlue+1);
133   fHEventsTrVtx->SetFillStyle(3001);
134   fHEventsTrVtx->SetDirectory(0);
135   // fHEventsTrVtx->Sumw2();
136
137       
138   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
139   fHTriggers->SetFillColor(kRed+1);
140   fHTriggers->SetFillStyle(3001);
141   fHTriggers->SetStats(0);
142   fHTriggers->SetDirectory(0);
143   fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
144   fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
145   fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
146   fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
147   fHTriggers->GetXaxis()->SetBinLabel(5,"A");
148   fHTriggers->GetXaxis()->SetBinLabel(6,"B");
149   fHTriggers->GetXaxis()->SetBinLabel(7,"C");
150   fHTriggers->GetXaxis()->SetBinLabel(8,"E");
151
152   TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
153   fHistos.Init(e);
154   fAODFMD.Init(e);
155
156   fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
157   fHData->SetStats(0);
158   fHData->SetDirectory(0);
159   fSharingFilter.Init();
160   fHistCollector.Init(*(fHEventsTr->GetXaxis()));
161 }
162
163 //____________________________________________________________________
164 void
165 AliForwardMultiplicity::UserCreateOutputObjects()
166 {
167   fList = new TList;
168
169   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
170   AliAODHandler*      ah = 
171     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
172   if (!ah)  
173     AliFatal("No AOD output handler set in analysis manager");
174     
175     
176   TObject* obj = &fAODFMD;
177   ah->AddBranch("AliAODForwardMult", &obj);
178
179   // fTree = new TTree("T", "T");
180   // fTree->Branch("forward", &fAODFMD);
181
182   PostData(1, fList);
183   // PostData(2, fTree);
184 }
185 //____________________________________________________________________
186 void
187 AliForwardMultiplicity::UserExec(Option_t*)
188 {
189   // Get the input data 
190   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
191   if (!esd) { 
192     AliWarning("No ESD event found for input event");
193     return;
194   }
195
196 #if 0
197   static Int_t nEvents = 0;
198   nEvents++;
199   if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
200 #endif
201
202   // On the first event, initialize the parameters 
203   if (fFirstEvent) { 
204     AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
205     pars->SetParametersFromESD(esd);
206     pars->PrintStatus();
207     fFirstEvent = false;
208
209     InitializeSubs();
210   }
211   // Clear stuff 
212   fHistos.Clear();
213   fESDFMD.Clear();
214   fAODFMD.Clear();
215
216   // Read trigger information from the ESD and store in AOD object
217   if (!ReadTriggers(esd)) { 
218 #ifdef VERBOSE
219     AliWarning("Failed to read triggers from ESD");
220 #endif
221     return;
222   }
223
224   // Mark this event for storage 
225   MarkEventForStore();
226
227   // Check if this is a high-flux event 
228   const AliMultiplicity* testmult = esd->GetMultiplicity();
229   if (!testmult) { 
230 #ifdef VERBOSE
231     AliWarning("No central multiplicity object found");
232 #endif
233     return;
234   }
235   Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
236
237   // Get the FMD ESD data 
238   AliESDFMD* esdFMD = esd->GetFMDData();
239   if (!esdFMD) { 
240 #ifdef VERBOSE
241     AliWarning("No FMD data found in ESD");
242 #endif
243     return;
244   }
245
246   // Get the vertex information 
247   Double_t vz   = 0;
248   Bool_t   vzOk = ReadVertex(esd, vz);
249
250   fHEventsTr->Fill(vz);
251   if (!vzOk) { 
252 #ifdef VERBOSE
253     AliWarning("Failed to read vertex from ESD");
254 #endif
255     return;
256   }
257   fHEventsTrVtx->Fill(vz);
258
259   // Get the vertex bin 
260   Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
261   fAODFMD.SetIpZ(vz);
262   if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) { 
263 #if 0
264     AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
265                     vz, fHEventsTr->GetXaxis()->GetXmin(), 
266                     fHEventsTr->GetXaxis()->GetXmax()));
267 #endif
268     return;
269   }
270
271   // Apply the sharing filter (or hit merging or clustering if you like)
272   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) { 
273 #ifdef VERBOSE
274     AliWarning("Sharing filter failed!");
275 #endif
276     return;
277   }
278
279   // Calculate the inclusive charged particle density 
280   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
281     AliWarning("Density calculator failed!");
282     return;
283   }
284   
285   // Do the secondary and other corrections. 
286   if (!fCorrections.Correct(fHistos, ivz)) { 
287     AliWarning("Corrections failed");
288     return;
289   }
290
291   if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
292     AliWarning("Histogram collector failed");
293     return;
294   }
295
296   if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
297     fHData->Add(&(fAODFMD.GetHistogram()));
298 }
299
300 //____________________________________________________________________
301 void
302 AliForwardMultiplicity::Terminate(Option_t*)
303 {
304   TList* list = dynamic_cast<TList*>(GetOutputData(1));
305   if (!list) {
306     AliError("No output list defined");
307     return;
308   }
309   // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
310   TH1D* dNdeta = fHData->ProjectionX("dNdeta", 1, -1, "e");
311   TH1D* norm   = fHData->ProjectionX("dNdeta", 0, 1,  "");
312   dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
313   dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
314   dNdeta->Divide(norm);
315   dNdeta->SetStats(0);
316   dNdeta->Scale(Double_t(fHEventsTrVtx->GetEntries())/fHEventsTr->GetEntries(),
317                 "width");
318
319   list->Add(fHEventsTr);
320   list->Add(fHEventsTrVtx);
321   list->Add(fHTriggers);
322   list->Add(fHData);
323   list->Add(dNdeta);
324   
325   TList* last = new TList;
326   last->SetName("LastEvent");
327   list->Add(last);
328   last->Add(&fAODFMD.GetHistogram());
329   last->Add(fHistos.fFMD1i);
330   last->Add(fHistos.fFMD2i);
331   last->Add(fHistos.fFMD2o);
332   last->Add(fHistos.fFMD3i);
333   last->Add(fHistos.fFMD3o);
334
335
336   fSharingFilter.ScaleHistograms(fHEventsTr->Integral());
337   fSharingFilter.Output(list);
338
339   fDensityCalculator.ScaleHistograms(fHEventsTrVtx->Integral());
340   fDensityCalculator.Output(list);
341
342   fCorrections.ScaleHistograms(fHEventsTrVtx->Integral());
343   fCorrections.Output(list);
344 }
345
346 //____________________________________________________________________
347 void
348 AliForwardMultiplicity::MarkEventForStore() const
349 {
350   // Make sure the AOD tree is filled 
351   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
352   AliAODHandler*      ah = 
353     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
354   if (!ah)  
355     AliFatal("No AOD output handler set in analysis manager");
356
357   ah->SetFillAOD(kTRUE);
358 }
359 //____________________________________________________________________
360 Bool_t
361 AliForwardMultiplicity::ReadTriggers(AliESDEvent* esd)
362 {
363   // Get the analysis manager - should always be there 
364   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
365   if (!am) { 
366     AliWarning("No analysis manager defined!");
367     return kFALSE;
368   }
369
370   // Get the input handler - should always be there 
371   AliInputEventHandler* ih = 
372     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
373   if (!ih) { 
374     AliWarning("No input handler");
375     return kFALSE;
376   }
377   
378   // Get the physics selection - add that by using the macro 
379   // AddTaskPhysicsSelection.C 
380   AliPhysicsSelection* ps = 
381     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
382   if (!ps) { 
383     AliWarning("No physics selection");
384     return kFALSE;
385   }
386   
387   // Check if this is a collision candidate (INEL)
388   Bool_t inel = ps->IsCollisionCandidate(esd);
389   if (inel) { 
390     fAODFMD.SetTriggerBits(AliAODForwardMult::kInel);
391     fHTriggers->Fill(.5);
392   }
393   
394
395   // IF this is inel, see if we have a tracklet 
396   if (inel) { 
397     const AliMultiplicity* spdmult = esd->GetMultiplicity();
398     if (!spdmult) {
399       AliWarning("No SPD multiplicity");
400     }
401     else { 
402       Int_t n = spdmult->GetNumberOfTracklets();
403       for (Int_t j = 0; j < n; j++) { 
404         if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
405           fAODFMD.SetTriggerBits(AliAODForwardMult::kInelGt0);
406           fHTriggers->Fill(1.5);
407           break;
408         }
409       }
410     }
411   }
412
413   // Analyse some trigger stuff 
414   AliTriggerAnalysis ta;
415   if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
416     fAODFMD.SetTriggerBits(AliAODForwardMult::kNSD);
417     fHTriggers->Fill(2.5);
418   }
419
420   // Get trigger stuff 
421   TString triggers = esd->GetFiredTriggerClasses();
422   if (triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
423     fAODFMD.SetTriggerBits(AliAODForwardMult::kEmpty);
424     fHTriggers->Fill(3.5);
425   }
426
427   if (triggers.Contains("CINT1A-ABCE-NOPF-ALL")) {
428     fAODFMD.SetTriggerBits(AliAODForwardMult::kA);
429     fHTriggers->Fill(4.5);
430   }
431
432   if (triggers.Contains("CINT1B-ABCE-NOPF-ALL")) {
433     fAODFMD.SetTriggerBits(AliAODForwardMult::kB);
434     fHTriggers->Fill(5.5);
435   }
436
437
438   if (triggers.Contains("CINT1C-ABCE-NOPF-ALL")) {
439     fAODFMD.SetTriggerBits(AliAODForwardMult::kC);
440     fHTriggers->Fill(6.5);
441   }
442
443   if (triggers.Contains("CINT1-E-NOPF-ALL")) {
444     fAODFMD.SetTriggerBits(AliAODForwardMult::kE);
445     fHTriggers->Fill(7.5);
446   }
447
448   return kTRUE;
449 }
450 //____________________________________________________________________
451 Bool_t
452 AliForwardMultiplicity::ReadVertex(AliESDEvent* esd, Double_t& vz)
453 {
454   // Get the vertex 
455   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
456   if (!vertex) { 
457 #ifdef VERBOSE
458     AliWarning("No SPD vertex found in ESD");
459 #endif
460     return kFALSE;
461   }
462
463   // Check that enough tracklets contributed 
464   if(vertex->GetNContributors() <= 0) {
465 #ifdef VERBOSE
466     AliWarning(Form("Number of contributors to vertex is %d<=0",
467                     vertex->GetNContributors()));
468 #endif
469     return kFALSE;
470   }
471
472   // Check that the uncertainty isn't too large 
473   if (vertex->GetZRes() > 0.1) { 
474 #ifdef VERBOSE
475     AliWarning(Form("Uncertaintity in Z of vertex is too large %f > 0.1", 
476                     vertex->GetZRes()));
477 #endif
478     return kFALSE;
479   }
480
481   // Get the z coordiante 
482   vz = vertex->GetZ();
483                
484   return kTRUE;
485 }
486
487 //____________________________________________________________________
488 void
489 AliForwardMultiplicity::Print(Option_t*) const
490 {
491 }
492
493 //
494 // EOF
495 //