More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEventInspector.cxx
1 #include "AliFMDEventInspector.h"
2 #include "AliLog.h"
3 #include "AliESDEvent.h"
4 #include "AliMultiplicity.h"
5 #include "AliAnalysisManager.h"
6 #include "AliInputEventHandler.h"
7 #include "AliTriggerAnalysis.h"
8 #include "AliPhysicsSelection.h"
9 #include "AliAODForwardMult.h"
10 #include "AliForwardUtil.h"
11 #include <TH1.h>
12 #include <TList.h>
13 #include <TDirectory.h>
14 #include <TROOT.h>
15 #include <iostream>
16 #include <iomanip>
17
18 //====================================================================
19 AliFMDEventInspector::AliFMDEventInspector()
20   : TNamed(),
21     fHEventsTr(0), 
22     fHEventsTrVtx(0),
23     fHTriggers(0),
24     fHType(0),
25     fLowFluxCut(1000),
26     fMaxVzErr(0.1),
27     fList(0),
28     fEnergy(0),
29     fField(999), 
30     fCollisionSystem(kUnknown),
31     fDebug(0)
32 {
33 }
34
35 //____________________________________________________________________
36 AliFMDEventInspector::AliFMDEventInspector(const char* name)
37   : TNamed("fmdEventInspector", name),
38     fHEventsTr(0), 
39     fHEventsTrVtx(0), 
40     fHTriggers(0),
41     fHType(0),
42     fLowFluxCut(1000),
43     fMaxVzErr(0.1),
44     fList(0),
45     fEnergy(0),
46     fField(999), 
47     fCollisionSystem(kUnknown),
48     fDebug(0)
49 {
50 }
51
52 //____________________________________________________________________
53 AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
54   : TNamed(o), 
55     fHEventsTr(o.fHEventsTr), 
56     fHEventsTrVtx(o.fHEventsTrVtx), 
57     fHTriggers(o.fHTriggers),
58     fHType(o.fHType),
59     fLowFluxCut(o.fMaxVzErr),
60     fMaxVzErr(o.fMaxVzErr),
61     fList(o.fList),
62     fEnergy(o.fEnergy),
63     fField(o.fField), 
64     fCollisionSystem(o.fCollisionSystem),
65     fDebug(0)
66 {
67 }
68
69 //____________________________________________________________________
70 AliFMDEventInspector::~AliFMDEventInspector()
71 {
72   if (fHEventsTr)    delete fHEventsTr;
73   if (fHEventsTrVtx) delete fHEventsTrVtx;
74   if (fHTriggers)    delete fHTriggers;  
75   if (fHType)        delete fHType;
76   if (fList)         delete fList;
77 }
78 //____________________________________________________________________
79 AliFMDEventInspector&
80 AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
81 {
82   TNamed::operator=(o);
83   fHEventsTr         = o.fHEventsTr;
84   fHEventsTrVtx      = o.fHEventsTrVtx;
85   fHTriggers         = o.fHTriggers;
86   fHType             = o.fHType;
87   fLowFluxCut        = o.fLowFluxCut;
88   fMaxVzErr          = o.fMaxVzErr;
89   fDebug             = o.fDebug;
90   fList              = (o.fList ? new TList : 0);
91   fEnergy            = o.fEnergy;
92   fField             = o.fField;
93   fCollisionSystem   = o.fCollisionSystem;
94   if (fList) { 
95     fList->SetName(GetName());
96     if (fHEventsTr)    fList->Add(fHEventsTr);
97     if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
98     if (fHTriggers)    fList->Add(fHTriggers);
99     if (fHType)        fList->Add(fHType);
100   }
101   return *this;
102 }
103
104 //____________________________________________________________________
105 Bool_t
106 AliFMDEventInspector::FetchHistograms(TList* d, 
107                                       TH1I*& hEventsTr, 
108                                       TH1I*& hEventsTrVtx, 
109                                       TH1I*& hTriggers) const
110 {
111   hEventsTr    = 0;
112   hEventsTrVtx = 0;
113   hTriggers    = 0;
114   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
115   if (!dd) return kFALSE;
116   
117   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
118   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
119   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
120
121   if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE;
122   return kTRUE;
123 }
124 //____________________________________________________________________
125 void
126 AliFMDEventInspector::Init(const TAxis& vtxAxis)
127 {
128   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
129                         vtxAxis.GetNbins(), 
130                         vtxAxis.GetXmin(), 
131                         vtxAxis.GetXmax());
132   fHEventsTr->SetXTitle("v_{z} [cm]");
133   fHEventsTr->SetYTitle("# of events");
134   fHEventsTr->SetFillColor(kRed+1);
135   fHEventsTr->SetFillStyle(3001);
136   fHEventsTr->SetDirectory(0);
137   // fHEventsTr->Sumw2();
138   fList->Add(fHEventsTr);
139
140   fHEventsTrVtx = new TH1I("nEventsTrVtx", 
141                            "Number of events w/trigger and vertex", 
142                            vtxAxis.GetNbins(), 
143                            vtxAxis.GetXmin(), 
144                            vtxAxis.GetXmax());
145   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
146   fHEventsTrVtx->SetYTitle("# of events");
147   fHEventsTrVtx->SetFillColor(kBlue+1);
148   fHEventsTrVtx->SetFillStyle(3001);
149   fHEventsTrVtx->SetDirectory(0);
150   // fHEventsTrVtx->Sumw2();
151   fList->Add(fHEventsTrVtx);
152
153       
154   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
155   fHTriggers->SetFillColor(kRed+1);
156   fHTriggers->SetFillStyle(3001);
157   fHTriggers->SetStats(0);
158   fHTriggers->SetDirectory(0);
159   fHTriggers->GetXaxis()->SetBinLabel(kInel   +1,"INEL");
160   fHTriggers->GetXaxis()->SetBinLabel(kInelGt0+1,"INEL>0");
161   fHTriggers->GetXaxis()->SetBinLabel(kNSD    +1,"NSD");
162   fHTriggers->GetXaxis()->SetBinLabel(kEmpty  +1,"Empty");
163   fHTriggers->GetXaxis()->SetBinLabel(kA      +1,"A");
164   fHTriggers->GetXaxis()->SetBinLabel(kB      +1,"B");
165   fHTriggers->GetXaxis()->SetBinLabel(kC      +1,"C");
166   fHTriggers->GetXaxis()->SetBinLabel(kE      +1,"E");
167   fHTriggers->GetXaxis()->SetBinLabel(9,         "spare1");
168   fHTriggers->GetXaxis()->SetBinLabel(10,        "spare2");
169   fList->Add(fHTriggers);
170
171   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)", 
172                                  fLowFluxCut), 2, -.5, 1.5);
173   fHType->SetFillColor(kRed+1);
174   fHType->SetFillStyle(3001);
175   fHType->SetStats(0);
176   fHType->SetDirectory(0);
177   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
178   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
179   fList->Add(fHType);
180   
181 }
182
183 //____________________________________________________________________
184 void
185 AliFMDEventInspector::DefineOutput(TList* dir)
186 {
187   fList = new TList;
188   fList->SetName(GetName());
189   dir->Add(fList);
190 }
191
192 //____________________________________________________________________
193 UInt_t
194 AliFMDEventInspector::Process(const AliESDEvent* event, 
195                               UInt_t&            triggers,
196                               Bool_t&            lowFlux,
197                               UShort_t&          ivz, 
198                               Double_t&          vz)
199 {
200   // Check that we have an event 
201   if (!event) { 
202     AliWarning("No ESD event found for input event");
203     return kNoEvent;
204   }
205
206   // Read trigger information from the ESD and store in AOD object
207   if (!ReadTriggers(event, triggers)) { 
208     if (fDebug > 2) {
209       AliWarning("Failed to read triggers from ESD"); }
210     return kNoTriggers;
211   }
212
213   // Check if this is a high-flux event 
214   const AliMultiplicity* testmult = event->GetMultiplicity();
215   if (!testmult) {
216     if (fDebug > 3) {
217       AliWarning("No central multiplicity object found"); }
218     return kNoSPD;
219   }
220   lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
221   
222   fHType->Fill(lowFlux ? 0 : 1);
223
224   // Check the FMD ESD data 
225   if (!event->GetFMDData()) { 
226     if (fDebug > 3) {
227       AliWarning("No FMD data found in ESD"); }
228     return kNoFMD;
229   }
230
231   // Get the vertex information 
232   vz          = 0;
233   Bool_t vzOk = ReadVertex(event, vz);
234
235   fHEventsTr->Fill(vz);
236   if (!vzOk) { 
237     if (fDebug > 3) {
238       AliWarning("Failed to read vertex from ESD"); }
239     return kNoVertex;
240   }
241   fHEventsTrVtx->Fill(vz);
242
243   // Get the vertex bin 
244   ivz = fHEventsTr->GetXaxis()->FindBin(vz);
245   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
246     if (fDebug > 3) {
247       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
248                       vz, fHEventsTr->GetXaxis()->GetXmin(), 
249                       fHEventsTr->GetXaxis()->GetXmax())); }
250     ivz = 0;
251     return kBadVertex;
252   }
253   return kOk;
254 }
255
256 //____________________________________________________________________
257 Bool_t
258 AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
259 {
260   triggers = 0;
261
262   // Get the analysis manager - should always be there 
263   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
264   if (!am) { 
265     AliWarning("No analysis manager defined!");
266     return kFALSE;
267   }
268
269   // Get the input handler - should always be there 
270   AliInputEventHandler* ih = 
271     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
272   if (!ih) { 
273     AliWarning("No input handler");
274     return kFALSE;
275   }
276   
277   // Get the physics selection - add that by using the macro 
278   // AddTaskPhysicsSelection.C 
279   AliPhysicsSelection* ps = 
280     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
281   if (!ps) { 
282     AliWarning("No physics selection");
283     return kFALSE;
284   }
285   
286   // Check if this is a collision candidate (INEL)
287   Bool_t inel = ps->IsCollisionCandidate(esd);
288   if (inel) { 
289     triggers |= AliAODForwardMult::kInel;
290     fHTriggers->Fill(kInel+0.5);
291   }
292
293   // IF this is inel, see if we have a tracklet 
294   if (inel) { 
295     const AliMultiplicity* spdmult = esd->GetMultiplicity();
296     if (!spdmult) {
297       AliWarning("No SPD multiplicity");
298     }
299     else { 
300       Int_t n = spdmult->GetNumberOfTracklets();
301       for (Int_t j = 0; j < n; j++) { 
302         if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
303           triggers |= AliAODForwardMult::kInelGt0;
304           fHTriggers->Fill(kInelGt0+.5);
305           break;
306         }
307       }
308     }
309   }
310
311   // Analyse some trigger stuff 
312   AliTriggerAnalysis ta;
313   if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
314     triggers |= AliAODForwardMult::kNSD;
315     fHTriggers->Fill(kNSD+.5);
316   }
317
318   // Get trigger stuff 
319   TString trigStr = esd->GetFiredTriggerClasses();
320   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
321     triggers |= AliAODForwardMult::kEmpty;
322     fHTriggers->Fill(kEmpty+.5);
323   }
324
325   if (trigStr.Contains("CINT1A-ABCE-NOPF-ALL")) {
326     triggers |= AliAODForwardMult::kA;
327     fHTriggers->Fill(kA+.5);
328   }
329
330   if (trigStr.Contains("CINT1B-ABCE-NOPF-ALL")) {
331     triggers |= AliAODForwardMult::kB;
332     fHTriggers->Fill(kB+.5);
333   }
334
335
336   if (trigStr.Contains("CINT1C-ABCE-NOPF-ALL")) {
337     triggers |= AliAODForwardMult::kC;
338     fHTriggers->Fill(kC+.5);
339   }
340
341   if (trigStr.Contains("CINT1-E-NOPF-ALL")) {
342     triggers |= AliAODForwardMult::kE;
343     fHTriggers->Fill(kE+.5);
344   }
345
346   return kTRUE;
347 }
348 //____________________________________________________________________
349 Bool_t
350 AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
351 {
352   vz = 0;
353   // Get the vertex 
354   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
355   if (!vertex) { 
356     if (fDebug > 2) {
357       AliWarning("No SPD vertex found in ESD"); }
358     return kFALSE;
359   }
360   
361   // Check that enough tracklets contributed 
362   if(vertex->GetNContributors() <= 0) {
363     if (fDebug > 2) {
364       AliWarning(Form("Number of contributors to vertex is %d<=0",
365                       vertex->GetNContributors())); }
366     vz = 0;
367     return kFALSE;
368   }
369
370   // Check that the uncertainty isn't too large 
371   if (vertex->GetZRes() > fMaxVzErr) { 
372     if (fDebug > 2) {
373       AliWarning(Form("Uncertaintity in Z of vertex is too large %f > %f", 
374                       vertex->GetZRes(), fMaxVzErr)); }
375     return kFALSE;
376   }
377
378   // Get the z coordiante 
379   vz = vertex->GetZ();
380   return kTRUE;
381 }
382
383 //____________________________________________________________________
384 Bool_t
385 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
386 {
387   fCollisionSystem = 
388     AliForwardUtil::ParseCollisionSystem(esd->GetBeamType());
389   fEnergy          = 
390     AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,   
391                                               2 * esd->GetBeamEnergy());
392   fField           = 
393     AliForwardUtil::ParseMagneticField(esd->GetMagneticField());
394   
395   if (fCollisionSystem   == AliForwardUtil::kUnknown || 
396       fEnergy            <= 0                        || 
397       TMath::Abs(fField) >  10) 
398     return kFALSE;
399
400   return kTRUE;
401 }
402
403 //____________________________________________________________________
404 void
405 AliFMDEventInspector::Print(Option_t*) const
406 {
407   char ind[gROOT->GetDirLevel()+1];
408   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
409   ind[gROOT->GetDirLevel()] = '\0';
410   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
411   sNN.Strip(TString::kBoth, '0');
412   sNN.ReplaceAll("GeV", " GeV");
413   TString field(AliForwardUtil::MagneticFieldString(fField));
414   field.ReplaceAll("p",  "+");
415   field.ReplaceAll("m",  "-");
416   field.ReplaceAll("kG", " kG");
417   
418   std::cout << ind << "AliFMDEventInspector: " << GetName() << '\n'
419             << ind << " Low flux cut:           " << fLowFluxCut << '\n'
420             << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
421             << ind << " System:                 " 
422             << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
423             << ind << " CMS energy per nucleon: " << sNN << '\n'
424             << ind << " Field:                  " <<  field << std::endl;
425 }
426
427   
428 //
429 // EOF
430 //
431