Various improvements
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEventInspector.cxx
1 // 
2 // This class inspects the event 
3 //
4 // Input:
5 //   - AliESDFMD object possibly corrected for sharing
6 //
7 // Output:
8 //   - A histogram of v_z of events with triggers. 
9 //   - A histogram of v_z of events with vertex and triggers 
10 //   - A histogram of trigger counters 
11 // 
12 // Note, that these are added to the master output list 
13 //
14 // Corrections used: 
15 //   - None
16 //
17 #include "AliFMDEventInspector.h"
18 #include "AliLog.h"
19 #include "AliESDEvent.h"
20 #include "AliMultiplicity.h"
21 #include "AliAnalysisManager.h"
22 #include "AliInputEventHandler.h"
23 #include "AliTriggerAnalysis.h"
24 #include "AliPhysicsSelection.h"
25 #include "AliAODForwardMult.h"
26 #include "AliForwardUtil.h"
27 #include <TH1.h>
28 #include <TList.h>
29 #include <TDirectory.h>
30 #include <TROOT.h>
31 #include <iostream>
32 #include <iomanip>
33
34 //====================================================================
35 AliFMDEventInspector::AliFMDEventInspector()
36   : TNamed(),
37     fHEventsTr(0), 
38     fHEventsTrVtx(0),
39     fHTriggers(0),
40     fHType(0),
41     fLowFluxCut(1000),
42     fMaxVzErr(0.1),
43     fList(0),
44     fEnergy(0),
45     fField(999), 
46     fCollisionSystem(kUnknown),
47     fDebug(0)
48 {
49   // 
50   // Constructor 
51   //
52 }
53
54 //____________________________________________________________________
55 AliFMDEventInspector::AliFMDEventInspector(const char* name)
56   : TNamed("fmdEventInspector", name),
57     fHEventsTr(0), 
58     fHEventsTrVtx(0), 
59     fHTriggers(0),
60     fHType(0),
61     fLowFluxCut(1000),
62     fMaxVzErr(0.1),
63     fList(0),
64     fEnergy(0),
65     fField(999), 
66     fCollisionSystem(kUnknown),
67     fDebug(0)
68 {
69   // 
70   // Constructor 
71   // 
72   // Parameters:
73   //   name Name of object
74   //
75 }
76
77 //____________________________________________________________________
78 AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
79   : TNamed(o), 
80     fHEventsTr(o.fHEventsTr), 
81     fHEventsTrVtx(o.fHEventsTrVtx), 
82     fHTriggers(o.fHTriggers),
83     fHType(o.fHType),
84     fLowFluxCut(o.fMaxVzErr),
85     fMaxVzErr(o.fMaxVzErr),
86     fList(o.fList),
87     fEnergy(o.fEnergy),
88     fField(o.fField), 
89     fCollisionSystem(o.fCollisionSystem),
90     fDebug(0)
91 {
92   // 
93   // Copy constructor 
94   // 
95   // Parameters:
96   //   o Object to copy from 
97   //
98 }
99
100 //____________________________________________________________________
101 AliFMDEventInspector::~AliFMDEventInspector()
102 {
103   // 
104   // Destructor 
105   //
106   if (fHEventsTr)    delete fHEventsTr;
107   if (fHEventsTrVtx) delete fHEventsTrVtx;
108   if (fHTriggers)    delete fHTriggers;  
109   if (fHType)        delete fHType;
110   if (fList)         delete fList;
111 }
112 //____________________________________________________________________
113 AliFMDEventInspector&
114 AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
115 {
116   // 
117   // Assignement operator
118   // 
119   // Parameters:
120   //   o Object to assign from 
121   // 
122   // Return:
123   //    Reference to this object
124   //
125   TNamed::operator=(o);
126   fHEventsTr         = o.fHEventsTr;
127   fHEventsTrVtx      = o.fHEventsTrVtx;
128   fHTriggers         = o.fHTriggers;
129   fHType             = o.fHType;
130   fLowFluxCut        = o.fLowFluxCut;
131   fMaxVzErr          = o.fMaxVzErr;
132   fDebug             = o.fDebug;
133   fList              = (o.fList ? new TList : 0);
134   fEnergy            = o.fEnergy;
135   fField             = o.fField;
136   fCollisionSystem   = o.fCollisionSystem;
137   if (fList) { 
138     fList->SetName(GetName());
139     if (fHEventsTr)    fList->Add(fHEventsTr);
140     if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
141     if (fHTriggers)    fList->Add(fHTriggers);
142     if (fHType)        fList->Add(fHType);
143   }
144   return *this;
145 }
146
147 //____________________________________________________________________
148 Bool_t
149 AliFMDEventInspector::FetchHistograms(TList* d, 
150                                       TH1I*& hEventsTr, 
151                                       TH1I*& hEventsTrVtx, 
152                                       TH1I*& hTriggers) const
153 {
154   // 
155   // Fetch our histograms from the passed list 
156   // 
157   // Parameters:
158   //   d             Input
159   //   hEventsTr     On return, pointer to histogram, or null
160   //   hEventsTrVtx  On return, pointer to histogram, or null
161   //   hTriggers     On return, pointer to histogram, or null
162   // 
163   // Return:
164   //    true on success, false otherwise 
165   //
166   hEventsTr    = 0;
167   hEventsTrVtx = 0;
168   hTriggers    = 0;
169   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
170   if (!dd) return kFALSE;
171   
172   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
173   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
174   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
175
176   if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE;
177   return kTRUE;
178 }
179 //____________________________________________________________________
180 void
181 AliFMDEventInspector::Init(const TAxis& vtxAxis)
182 {
183   // 
184   // Initialize the object 
185   // 
186   // Parameters:
187   //   vtxAxis Vertex axis in use 
188   //
189   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
190                         vtxAxis.GetNbins(), 
191                         vtxAxis.GetXmin(), 
192                         vtxAxis.GetXmax());
193   fHEventsTr->SetXTitle("v_{z} [cm]");
194   fHEventsTr->SetYTitle("# of events");
195   fHEventsTr->SetFillColor(kRed+1);
196   fHEventsTr->SetFillStyle(3001);
197   fHEventsTr->SetDirectory(0);
198   // fHEventsTr->Sumw2();
199   fList->Add(fHEventsTr);
200
201   fHEventsTrVtx = new TH1I("nEventsTrVtx", 
202                            "Number of events w/trigger and vertex", 
203                            vtxAxis.GetNbins(), 
204                            vtxAxis.GetXmin(), 
205                            vtxAxis.GetXmax());
206   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
207   fHEventsTrVtx->SetYTitle("# of events");
208   fHEventsTrVtx->SetFillColor(kBlue+1);
209   fHEventsTrVtx->SetFillStyle(3001);
210   fHEventsTrVtx->SetDirectory(0);
211   // fHEventsTrVtx->Sumw2();
212   fList->Add(fHEventsTrVtx);
213
214       
215   fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
216   fHTriggers->SetFillColor(kRed+1);
217   fHTriggers->SetFillStyle(3001);
218   fHTriggers->SetStats(0);
219   fHTriggers->SetDirectory(0);
220   fHTriggers->GetXaxis()->SetBinLabel(kInel   +1,"INEL");
221   fHTriggers->GetXaxis()->SetBinLabel(kInelGt0+1,"INEL>0");
222   fHTriggers->GetXaxis()->SetBinLabel(kNSD    +1,"NSD");
223   fHTriggers->GetXaxis()->SetBinLabel(kEmpty  +1,"Empty");
224   fHTriggers->GetXaxis()->SetBinLabel(kA      +1,"A");
225   fHTriggers->GetXaxis()->SetBinLabel(kB      +1,"B");
226   fHTriggers->GetXaxis()->SetBinLabel(kC      +1,"C");
227   fHTriggers->GetXaxis()->SetBinLabel(kE      +1,"E");
228   fHTriggers->GetXaxis()->SetBinLabel(9,         "spare1");
229   fHTriggers->GetXaxis()->SetBinLabel(10,        "spare2");
230   fList->Add(fHTriggers);
231
232   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)", 
233                                  fLowFluxCut), 2, -.5, 1.5);
234   fHType->SetFillColor(kRed+1);
235   fHType->SetFillStyle(3001);
236   fHType->SetStats(0);
237   fHType->SetDirectory(0);
238   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
239   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
240   fList->Add(fHType);
241   
242 }
243
244 //____________________________________________________________________
245 void
246 AliFMDEventInspector::DefineOutput(TList* dir)
247 {
248   // 
249   // Define the output histograms.  These are put in a sub list of the
250   // passed list.   The histograms are merged before the parent task calls 
251   // AliAnalysisTaskSE::Terminate 
252   // 
253   //   dir Directory to add to 
254   //
255   fList = new TList;
256   fList->SetName(GetName());
257   dir->Add(fList);
258 }
259
260 //____________________________________________________________________
261 UInt_t
262 AliFMDEventInspector::Process(const AliESDEvent* event, 
263                               UInt_t&            triggers,
264                               Bool_t&            lowFlux,
265                               UShort_t&          ivz, 
266                               Double_t&          vz)
267 {
268   // 
269   // Process the event 
270   // 
271   // Parameters:
272   //   event     Input event 
273   //   triggers  On return, the triggers fired 
274   //   lowFlux   On return, true if the event is considered a low-flux 
275   //                  event (according to the setting of fLowFluxCut) 
276   //   ivz       On return, the found vertex bin (1-based).  A zero
277   //                  means outside of the defined vertex range
278   //   vz        On return, the z position of the interaction
279   // 
280   // Return:
281   //    0 (or kOk) on success, otherwise a bit mask of error codes 
282   //
283
284   // Check that we have an event 
285   if (!event) { 
286     AliWarning("No ESD event found for input event");
287     return kNoEvent;
288   }
289
290   // Read trigger information from the ESD and store in AOD object
291   if (!ReadTriggers(event, triggers)) { 
292     if (fDebug > 2) {
293       AliWarning("Failed to read triggers from ESD"); }
294     return kNoTriggers;
295   }
296
297   // Check if this is a high-flux event 
298   const AliMultiplicity* testmult = event->GetMultiplicity();
299   if (!testmult) {
300     if (fDebug > 3) {
301       AliWarning("No central multiplicity object found"); }
302     return kNoSPD;
303   }
304   lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
305   
306   fHType->Fill(lowFlux ? 0 : 1);
307
308   // Check the FMD ESD data 
309   if (!event->GetFMDData()) { 
310     if (fDebug > 3) {
311       AliWarning("No FMD data found in ESD"); }
312     return kNoFMD;
313   }
314
315   // Get the vertex information 
316   vz          = 0;
317   Bool_t vzOk = ReadVertex(event, vz);
318
319   fHEventsTr->Fill(vz);
320   if (!vzOk) { 
321     if (fDebug > 3) {
322       AliWarning("Failed to read vertex from ESD"); }
323     return kNoVertex;
324   }
325   fHEventsTrVtx->Fill(vz);
326
327   // Get the vertex bin 
328   ivz = fHEventsTr->GetXaxis()->FindBin(vz);
329   if (ivz <= 0 || ivz > fHEventsTr->GetXaxis()->GetNbins()) { 
330     if (fDebug > 3) {
331       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
332                       vz, fHEventsTr->GetXaxis()->GetXmin(), 
333                       fHEventsTr->GetXaxis()->GetXmax())); }
334     ivz = 0;
335     return kBadVertex;
336   }
337   return kOk;
338 }
339
340 //____________________________________________________________________
341 Bool_t
342 AliFMDEventInspector::ReadTriggers(const AliESDEvent* esd, UInt_t& triggers)
343 {
344   // 
345   // Read the trigger information from the ESD event 
346   // 
347   // Parameters:
348   //   esd        ESD event 
349   //   triggers   On return, contains the trigger bits 
350   // 
351   // Return:
352   //    @c true on success, @c false otherwise 
353   //
354   triggers = 0;
355
356   // Get the analysis manager - should always be there 
357   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
358   if (!am) { 
359     AliWarning("No analysis manager defined!");
360     return kFALSE;
361   }
362
363   // Get the input handler - should always be there 
364   AliInputEventHandler* ih = 
365     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
366   if (!ih) { 
367     AliWarning("No input handler");
368     return kFALSE;
369   }
370
371 #if 1
372   // Check if this is a collision candidate (INEL)
373   // Note, that we should use the value cached in the input 
374   // handler rather than calling IsCollisionCandiate directly 
375   // on the AliPhysicsSelection obejct.  If we called the latter
376   // then the AliPhysicsSelection object would overcount by a 
377   // factor of 2! :-(
378   Bool_t inel = ih->IsEventSelected();
379 #else 
380   // Get the physics selection - add that by using the macro 
381   // AddTaskPhysicsSelection.C 
382   AliPhysicsSelection* ps = 
383     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
384   if (!ps) { 
385     AliWarning("No physics selection");
386     return kFALSE;
387   }
388   
389   // Check if this is a collision candidate (INEL)
390   Bool_t inel = ps->IsCollisionCandidate(esd);
391 #endif
392   if (inel) { 
393     triggers |= AliAODForwardMult::kInel;
394     fHTriggers->Fill(kInel+0.5);
395   }
396
397   // If this is inel, see if we have a tracklet 
398   if (inel) { 
399     const AliMultiplicity* spdmult = esd->GetMultiplicity();
400     if (!spdmult) {
401       AliWarning("No SPD multiplicity");
402     }
403     else { 
404       // Check if we have one or more tracklets 
405       // in the range -1 < eta < 1 to set the INEL>0 
406       // trigger flag. 
407       Int_t n = spdmult->GetNumberOfTracklets();
408       for (Int_t j = 0; j < n; j++) { 
409         if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
410           triggers |= AliAODForwardMult::kInelGt0;
411           fHTriggers->Fill(kInelGt0+.5);
412           break;
413         }
414       }
415     }
416   }
417
418   // Analyse some trigger stuff 
419   AliTriggerAnalysis ta;
420   if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
421     triggers |= AliAODForwardMult::kNSD;
422     fHTriggers->Fill(kNSD+.5);
423   }
424
425   // Get trigger stuff 
426   TString trigStr = esd->GetFiredTriggerClasses();
427   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
428     triggers |= AliAODForwardMult::kEmpty;
429     fHTriggers->Fill(kEmpty+.5);
430   }
431
432   if (trigStr.Contains("CINT1A-ABCE-NOPF-ALL")) {
433     triggers |= AliAODForwardMult::kA;
434     fHTriggers->Fill(kA+.5);
435   }
436
437   if (trigStr.Contains("CINT1B-ABCE-NOPF-ALL")) {
438     triggers |= AliAODForwardMult::kB;
439     fHTriggers->Fill(kB+.5);
440   }
441
442
443   if (trigStr.Contains("CINT1C-ABCE-NOPF-ALL")) {
444     triggers |= AliAODForwardMult::kC;
445     fHTriggers->Fill(kC+.5);
446   }
447
448   if (trigStr.Contains("CINT1-E-NOPF-ALL")) {
449     triggers |= AliAODForwardMult::kE;
450     fHTriggers->Fill(kE+.5);
451   }
452
453   return kTRUE;
454 }
455 //____________________________________________________________________
456 Bool_t
457 AliFMDEventInspector::ReadVertex(const AliESDEvent* esd, Double_t& vz)
458 {
459   // 
460   // Read the vertex information from the ESD event 
461   // 
462   // Parameters:
463   //   esd  ESD event 
464   //   vz   On return, the vertex Z position 
465   // 
466   // Return:
467   //    @c true on success, @c false otherwise 
468   //
469   vz = 0;
470   // Get the vertex 
471   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
472   if (!vertex) { 
473     if (fDebug > 2) {
474       AliWarning("No SPD vertex found in ESD"); }
475     return kFALSE;
476   }
477   
478   // Check that enough tracklets contributed 
479   if(vertex->GetNContributors() <= 0) {
480     if (fDebug > 2) {
481       AliWarning(Form("Number of contributors to vertex is %d<=0",
482                       vertex->GetNContributors())); }
483     vz = 0;
484     return kFALSE;
485   }
486
487   // Check that the uncertainty isn't too large 
488   if (vertex->GetZRes() > fMaxVzErr) { 
489     if (fDebug > 2) {
490       AliWarning(Form("Uncertaintity in Z of vertex is too large %f > %f", 
491                       vertex->GetZRes(), fMaxVzErr)); }
492     return kFALSE;
493   }
494
495   // Get the z coordiante 
496   vz = vertex->GetZ();
497   return kTRUE;
498 }
499
500 //____________________________________________________________________
501 Bool_t
502 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
503 {
504   // 
505   // Read the collision system, collision energy, and L3 field setting
506   // from the ESD
507   // 
508   // Parameters:
509   //   esd ESD to get information from 
510   // 
511   // Return:
512   //    true on success, false 
513   //
514   // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
515   //           esd->GetBeamType(), 2*esd->GetBeamEnergy(), 
516   //           esd->GetMagneticField()));
517   fCollisionSystem = 
518     AliForwardUtil::ParseCollisionSystem(esd->GetBeamType());
519   fEnergy          = 
520     AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem,   
521                                               2 * esd->GetBeamEnergy());
522   fField           = 
523     AliForwardUtil::ParseMagneticField(esd->GetMagneticField());
524   
525   if (fCollisionSystem   == AliForwardUtil::kUnknown || 
526       fEnergy            <= 0                        || 
527       TMath::Abs(fField) >  10) 
528     return kFALSE;
529
530   return kTRUE;
531 }
532
533 //____________________________________________________________________
534 void
535 AliFMDEventInspector::Print(Option_t*) const
536 {
537   // 
538   // Print information
539   // 
540   //   option Not used 
541   //
542   char ind[gROOT->GetDirLevel()+1];
543   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
544   ind[gROOT->GetDirLevel()] = '\0';
545   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
546   sNN.Strip(TString::kBoth, '0');
547   sNN.ReplaceAll("GeV", " GeV");
548   TString field(AliForwardUtil::MagneticFieldString(fField));
549   field.ReplaceAll("p",  "+");
550   field.ReplaceAll("m",  "-");
551   field.ReplaceAll("kG", " kG");
552   
553   std::cout << ind << "AliFMDEventInspector: " << GetName() << '\n'
554             << ind << " Low flux cut:           " << fLowFluxCut << '\n'
555             << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
556             << ind << " System:                 " 
557             << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
558             << ind << " CMS energy per nucleon: " << sNN << '\n'
559             << ind << " Field:                  " <<  field << std::endl;
560 }
561
562   
563 //
564 // EOF
565 //
566