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