93f2de5520fbb9ac872b776a55458a052c32ca5c
[u/mrichter/AliRoot.git] / PWGLF / 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 "AliMCEventHandler.h"
23 #include "AliInputEventHandler.h"
24 #include "AliTriggerAnalysis.h"
25 #include "AliPhysicsSelection.h"
26 #include "AliOADBPhysicsSelection.h"
27 #include "AliAODForwardMult.h"
28 #include "AliForwardUtil.h"
29 #include "AliCentrality.h"
30 #include <TH1.h>
31 #include <TList.h>
32 #include <TDirectory.h>
33 #include <TROOT.h>
34 #include <TParameter.h>
35 #include <iostream>
36 #include <iomanip>
37 #include "AliMCEvent.h"
38 #include "AliHeader.h"
39 #include "AliGenEventHeader.h"
40 #include "AliCollisionGeometry.h"
41
42 //====================================================================
43 AliFMDEventInspector::AliFMDEventInspector()
44   : TNamed(),
45     fHEventsTr(0), 
46     fHEventsTrVtx(0),
47     fHEventsAccepted(0),
48     fHEventsAcceptedXY(0),
49     fHTriggers(0),
50     fHType(0),
51   fHWords(0),
52     fHCent(0),
53   fHCentVsQual(0),
54   fHStatus(0),
55     fLowFluxCut(1000),
56     fMaxVzErr(0.2),
57     fList(0),
58     fEnergy(0),
59     fField(999), 
60     fCollisionSystem(kUnknown),
61     fDebug(0),
62     fCentAxis(0),
63     fVtxAxis(10,-10,10),
64     fUseFirstPhysicsVertex(true),
65     fUseV0AND(false),
66     fMinPileupContrib(3), 
67     fMinPileupDistance(0.8),
68     fUseDisplacedVertices(false),
69   fDisplacedVertex(),
70   fCollWords(),
71   fBgWords()
72 {
73   // 
74   // Constructor 
75   //
76   DGUARD(fDebug,1,"Default CTOR of AliFMDEventInspector");
77 }
78
79 //____________________________________________________________________
80 AliFMDEventInspector::AliFMDEventInspector(const char* name)
81   : TNamed("fmdEventInspector", name),
82     fHEventsTr(0), 
83     fHEventsTrVtx(0), 
84     fHEventsAccepted(0),
85     fHEventsAcceptedXY(0),
86     fHTriggers(0),
87     fHType(0),
88     fHWords(0),
89     fHCent(0),
90     fHCentVsQual(0),
91     fHStatus(0),
92     fLowFluxCut(1000),
93     fMaxVzErr(0.2),
94     fList(0),
95     fEnergy(0),
96     fField(999), 
97     fCollisionSystem(kUnknown),
98     fDebug(0),
99     fCentAxis(0),
100     fVtxAxis(10,-10,10),
101     fUseFirstPhysicsVertex(true),
102     fUseV0AND(false),
103     fMinPileupContrib(3), 
104     fMinPileupDistance(0.8),
105     fUseDisplacedVertices(false),
106   fDisplacedVertex(),
107   fCollWords(),
108   fBgWords()
109 {
110   // 
111   // Constructor 
112   // 
113   // Parameters:
114   //   name Name of object
115   //
116   DGUARD(fDebug,1,"Named CTOR of AliFMDEventInspector: %s", name);
117 }
118
119 //____________________________________________________________________
120 AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
121   : TNamed(o), 
122     fHEventsTr(o.fHEventsTr), 
123     fHEventsTrVtx(o.fHEventsTrVtx), 
124     fHEventsAccepted(o.fHEventsAccepted),
125     fHEventsAcceptedXY(o.fHEventsAcceptedXY),
126     fHTriggers(o.fHTriggers),
127     fHType(o.fHType),
128     fHWords(o.fHWords),
129     fHCent(o.fHCent),
130     fHCentVsQual(o.fHCentVsQual),
131     fHStatus(o.fHStatus),
132     fLowFluxCut(o.fLowFluxCut),
133     fMaxVzErr(o.fMaxVzErr),
134     fList(o.fList),
135     fEnergy(o.fEnergy),
136     fField(o.fField), 
137     fCollisionSystem(o.fCollisionSystem),
138     fDebug(0),
139     fCentAxis(0),
140     fVtxAxis(o.fVtxAxis),
141     fUseFirstPhysicsVertex(o.fUseFirstPhysicsVertex),
142     fUseV0AND(o.fUseV0AND),
143     fMinPileupContrib(o.fMinPileupContrib), 
144     fMinPileupDistance(o.fMinPileupDistance),
145     fUseDisplacedVertices(o.fUseDisplacedVertices),
146     fDisplacedVertex(o.fDisplacedVertex),
147   fCollWords(),
148   fBgWords()
149 {
150   // 
151   // Copy constructor 
152   // 
153   // Parameters:
154   //   o Object to copy from 
155   //
156   DGUARD(fDebug,1,"Copy CTOR of AliFMDEventInspector");
157 }
158
159 //____________________________________________________________________
160 AliFMDEventInspector::~AliFMDEventInspector()
161 {
162   // 
163   // Destructor 
164   //
165   if (fList)         delete fList;
166 }
167 //____________________________________________________________________
168 AliFMDEventInspector&
169 AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
170 {
171   // 
172   // Assignement operator
173   // 
174   // Parameters:
175   //   o Object to assign from 
176   // 
177   // Return:
178   //    Reference to this object
179   //
180   DGUARD(fDebug,3,"Assignment of AliFMDEventInspector");
181   if (&o == this) return *this; 
182   TNamed::operator=(o);
183   fHEventsTr         = o.fHEventsTr;
184   fHEventsTrVtx      = o.fHEventsTrVtx;
185   fHEventsAccepted   = o.fHEventsAccepted;
186   fHEventsAcceptedXY = o.fHEventsAcceptedXY;
187   fHTriggers         = o.fHTriggers;
188   fHType             = o.fHType;
189   fHWords            = o.fHWords;
190   fHCent             = o.fHCent;
191   fHCentVsQual       = o.fHCentVsQual;
192   fHStatus           = o.fHStatus;
193   fLowFluxCut        = o.fLowFluxCut;
194   fMaxVzErr          = o.fMaxVzErr;
195   fDebug             = o.fDebug;
196   fList              = (o.fList ? new TList : 0);
197   fEnergy            = o.fEnergy;
198   fField             = o.fField;
199   fCollisionSystem   = o.fCollisionSystem;
200   fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), 
201                o.fVtxAxis.GetXmax());
202   
203   fUseFirstPhysicsVertex = o.fUseFirstPhysicsVertex;
204   fUseV0AND              = o.fUseV0AND;
205   fMinPileupContrib      = o.fMinPileupContrib;
206   fMinPileupDistance     = o.fMinPileupDistance;
207   fUseDisplacedVertices  = o.fUseDisplacedVertices;
208   fDisplacedVertex       = o.fDisplacedVertex;
209   if (fList) { 
210     fList->SetName(GetName());
211     if (fHEventsTr)    fList->Add(fHEventsTr);
212     if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
213     if (fHTriggers)    fList->Add(fHTriggers);
214     if (fHType)        fList->Add(fHType);
215     if (fHWords)       fList->Add(fHWords);
216     if (fHCent)        fList->Add(fHCent);
217     if (fHCentVsQual)  fList->Add(fHCentVsQual);
218     if (fHStatus)      fList->Add(fHStatus);
219   }
220   return *this;
221 }
222
223 //____________________________________________________________________
224 Bool_t
225 AliFMDEventInspector::FetchHistograms(const TList* d, 
226                                       TH1I*& hEventsTr, 
227                                       TH1I*& hEventsTrVtx, 
228                                       TH1I*& hTriggers) const
229 {
230   // 
231   // Fetch our histograms from the passed list 
232   // 
233   // Parameters:
234   //   d             Input
235   //   hEventsTr     On return, pointer to histogram, or null
236   //   hEventsTrVtx  On return, pointer to histogram, or null
237   //   hTriggers     On return, pointer to histogram, or null
238   // 
239   // Return:
240   //    true on success, false otherwise 
241   //
242   DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector");
243   hEventsTr    = 0;
244   hEventsTrVtx = 0;
245   hTriggers    = 0;
246   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
247   if (!dd) return kFALSE;
248   
249   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
250   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
251   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
252
253   if (!hEventsTr || !hEventsTrVtx || !hTriggers) return kFALSE;
254   return kTRUE;
255 }
256 //____________________________________________________________________
257 void
258 AliFMDEventInspector::CacheConfiguredTriggerClasses(TList& cache, 
259                                                     const TList* classes,
260                                                     AliOADBPhysicsSelection* o)
261 {
262   TIter nextClass(classes);
263   TObjString* trigClass = 0;
264   // Loop over all trigger classes.  Trigger classes have the format 
265   //
266   //   class          := positive_words SPACE(s) negative_words 
267   //   positive_words := 
268   //                  |  '+' words
269   //   negative_words := 
270   //                  |  '-' words
271   //   words          := word 
272   //                  |  word ',' words 
273   //   
274   while ((trigClass = static_cast<TObjString*>(nextClass()))) {
275     // Tokenize on space to get positive and negative parts 
276     TString     side   = o->GetBeamSide(trigClass->String());
277     TObjArray*  parts  = trigClass->String().Tokenize(" ");
278     TObjString* part   = 0;
279     TIter       nextPart(parts);
280     while ((part = static_cast<TObjString*>(nextPart()))) {
281       // We only care about the positive ones 
282       if (part->GetName()[0] != '+') continue;
283       part->String().Remove(0,1);
284         
285       // Tokenize on a comma to get the words 
286       TObjArray*  words = part->String().Tokenize(",");
287       TObjString* word  = 0;
288       TIter       nextWord(words);
289       while ((word = static_cast<TObjString*>(nextWord()))) {
290         TNamed* store = new TNamed(word->String(), side);
291         cache.Add(store);
292         DMSG(fDebug,3,"Caching %s trigger word %s", 
293              store->GetTitle(), store->GetName());
294       } // while (word)
295       delete words;
296     }
297     delete parts;
298   }
299 }
300
301 //____________________________________________________________________
302 void
303 AliFMDEventInspector::Init(const TAxis& vtxAxis)
304 {
305   // 
306   // Initialize the object - this is called on the first seen event. 
307   // 
308   // Parameters:
309   //   vtxAxis Vertex axis in use 
310   //
311   DGUARD(fDebug,1,"Initialize in AliFMDEventInspector");
312
313   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
314
315   // Get the input handler - should always be there 
316   AliInputEventHandler* ih = 
317     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
318   if (!ih) { 
319     AliWarning("No input handler");
320     return;
321   }
322   // Get the physics selection - should always be there 
323   AliPhysicsSelection* ps = 
324     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
325   if (!ps) {
326     AliWarning("No physics selection");
327     return;
328   }
329   // Get the configured triggers
330   AliOADBPhysicsSelection* oadb = 
331     const_cast<AliOADBPhysicsSelection*>(ps->GetOADBPhysicsSelection());
332   if (!oadb) {
333     AliWarning("No OADB physics selection object");
334     return;
335   }
336   // Get the configured trigger words from the physics selection 
337   const TList* collTriggClasses = ps->GetCollisionTriggerClasses();
338   const TList* bgTriggClasses   = ps->GetBGTriggerClasses();
339   if (!collTriggClasses) { 
340     AliWarning("No configured collision trigger classes");
341     return;
342   }
343   if (!bgTriggClasses) { 
344     AliWarning("No configured background trigger classes");
345     return;
346   }
347   CacheConfiguredTriggerClasses(fCollWords, collTriggClasses, oadb);
348   CacheConfiguredTriggerClasses(fBgWords,   bgTriggClasses,   oadb);
349   // fCollWords.ls();
350   // fBgWords.ls();
351   
352   
353   // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
354   // ----- 92 number --------- ---- 1 ---
355   TArrayD limits(93);
356   for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
357   limits[92] = 100.5;
358
359   fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
360   
361   fCentAxis  = new TAxis(limits.GetSize()-1, limits.GetArray());
362   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
363                         4*vtxAxis.GetNbins(), 
364                         2*vtxAxis.GetXmin(), 
365                         2*vtxAxis.GetXmax());
366   fHEventsTr->SetXTitle("v_{z} [cm]");
367   fHEventsTr->SetYTitle("# of events");
368   fHEventsTr->SetFillColor(kRed+1);
369   fHEventsTr->SetFillStyle(3001);
370   fHEventsTr->SetDirectory(0);
371   // fHEventsTr->Sumw2();
372   fList->Add(fHEventsTr);
373
374   fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
375   fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex"); 
376   fHEventsTrVtx->SetFillColor(kBlue+1);
377   fHEventsTrVtx->SetDirectory(0);
378   // fHEventsTrVtx->Sumw2();
379   fList->Add(fHEventsTrVtx);
380
381   fHEventsAccepted = new TH1I("nEventsAccepted", 
382                               "Number of events  w/trigger and vertex in range",
383                               2*vtxAxis.GetNbins(), 
384                               2*vtxAxis.GetXmin(), 
385                               2*vtxAxis.GetXmax());
386   fHEventsAccepted->SetXTitle("v_{z} [cm]");
387   fHEventsAccepted->SetYTitle("# of events");
388   fHEventsAccepted->SetFillColor(kGreen+1);
389   fHEventsAccepted->SetFillStyle(3001);
390   fHEventsAccepted->SetDirectory(0);
391   // fHEventsAccepted->Sumw2();
392   fList->Add(fHEventsAccepted);
393                               
394   fHEventsAcceptedXY = new TH2D("nEventsAcceptedXY", 
395                               "XY vertex w/trigger and Z vertex in range",
396                                 1000,-1,1,1000,-1,1);
397   
398   fHEventsAcceptedXY->SetXTitle("v_{x} [cm]");
399   fHEventsAcceptedXY->SetYTitle("v_{y} [cm]");
400   fHEventsAcceptedXY->SetDirectory(0);
401   // fHEventsAccepted->Sumw2();
402   fList->Add(fHEventsAcceptedXY);
403
404   
405   fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1);
406   fHTriggers->SetFillColor(kRed+1);
407   fHTriggers->SetFillStyle(3001);
408   fHTriggers->SetStats(0);
409   fHTriggers->SetDirectory(0);
410   fHTriggers->GetXaxis()->SetBinLabel(kInel   +1,"INEL");
411   fHTriggers->GetXaxis()->SetBinLabel(kInelGt0+1,"INEL>0");
412   fHTriggers->GetXaxis()->SetBinLabel(kNSD    +1,"NSD");
413   fHTriggers->GetXaxis()->SetBinLabel(kV0AND  +1,"VOAND");
414   fHTriggers->GetXaxis()->SetBinLabel(kEmpty  +1,"Empty");
415   fHTriggers->GetXaxis()->SetBinLabel(kA      +1,"A");
416   fHTriggers->GetXaxis()->SetBinLabel(kB      +1,"B");
417   fHTriggers->GetXaxis()->SetBinLabel(kC      +1,"C");
418   fHTriggers->GetXaxis()->SetBinLabel(kE      +1,"E");
419   fHTriggers->GetXaxis()->SetBinLabel(kPileUp +1,"Pileup");
420   fHTriggers->GetXaxis()->SetBinLabel(kMCNSD  +1,"NSD_{MC}");
421   fHTriggers->GetXaxis()->SetBinLabel(kOffline+1,"Offline");
422   fList->Add(fHTriggers);
423
424   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)", 
425                                  fLowFluxCut), 2, -.5, 1.5);
426   fHType->SetFillColor(kRed+1);
427   fHType->SetFillStyle(3001);
428   fHType->SetStats(0);
429   fHType->SetDirectory(0);
430   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
431   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
432   fList->Add(fHType);
433
434 #if 0 
435   // This histogram disabled as it causes problems in the merge 
436   fHWords = new TH1I("words", "Trigger words seen", 1, 0, 0); 
437   fHWords->SetFillColor(kBlue+1);
438   fHWords->SetFillStyle(3001);
439   fHWords->SetStats(0);
440   fHWords->SetDirectory(0);
441   fHWords->SetBit(TH1::kCanRebin);
442   fList->Add(fHWords);
443 #endif
444
445   fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
446   fHCent->SetFillColor(kBlue+1);
447   fHCent->SetFillStyle(3001);
448   fHCent->SetStats(0);
449   fHCent->SetDirectory(0);
450   fHCent->SetXTitle("Centrality [%]");
451   fHCent->SetYTitle("Events");
452   fList->Add(fHCent);
453
454   fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality", 
455                           5, 0, 5, limits.GetSize()-1, limits.GetArray());
456   fHCentVsQual->SetXTitle("Quality");
457   fHCentVsQual->SetYTitle("Centrality [%]");
458   fHCentVsQual->SetZTitle("Events");
459   fHCentVsQual->GetXaxis()->SetBinLabel(1, "OK");
460   fHCentVsQual->GetXaxis()->SetBinLabel(2, "Outside v_{z} cut");
461   fHCentVsQual->GetXaxis()->SetBinLabel(3, "V0 vs SPD outlier");
462   fHCentVsQual->GetXaxis()->SetBinLabel(4, "V0 vs TPC outlier");
463   fHCentVsQual->GetXaxis()->SetBinLabel(5, "V0 vs ZDC outlier");
464   fList->Add(fHCentVsQual);
465
466   fHStatus = new TH1I("status", "Status", 7, 1, 8);
467   fHStatus->SetFillColor(kRed+1);
468   fHStatus->SetFillStyle(3001);
469   fHStatus->SetStats(0);
470   fHStatus->SetDirectory(0);
471   fHStatus->GetXaxis()->SetBinLabel(1, "OK");
472   fHStatus->GetXaxis()->SetBinLabel(2, "No event");
473   fHStatus->GetXaxis()->SetBinLabel(3, "No triggers");
474   fHStatus->GetXaxis()->SetBinLabel(4, "No SPD");
475   fHStatus->GetXaxis()->SetBinLabel(5, "No FMD");
476   fHStatus->GetXaxis()->SetBinLabel(6, "No vertex");
477   fHStatus->GetXaxis()->SetBinLabel(7, "Bad vertex");
478   fList->Add(fHStatus);
479 }
480
481 //____________________________________________________________________
482 void
483 AliFMDEventInspector::StoreInformation(Int_t runNo)
484 {
485   // Write TNamed objects to output list containing information about
486   // the running conditions 
487   DGUARD(fDebug,2,"Store information from AliFMDEventInspector");
488   if (!fList) return;
489
490   
491   fList->Add(AliForwardUtil::MakeParameter("sys", fCollisionSystem));
492   fList->Add(AliForwardUtil::MakeParameter("sNN", fEnergy));
493   fList->Add(AliForwardUtil::MakeParameter("field", fField));
494   fList->Add(AliForwardUtil::MakeParameter("runNo", runNo));
495   fList->Add(AliForwardUtil::MakeParameter("lowFlux", fLowFluxCut));
496   fList->Add(AliForwardUtil::MakeParameter("fpVtx",fUseFirstPhysicsVertex));
497   fList->Add(AliForwardUtil::MakeParameter("v0and",fUseV0AND));
498   fList->Add(AliForwardUtil::MakeParameter("nPileUp", fMinPileupContrib));
499   fList->Add(AliForwardUtil::MakeParameter("dPileup", fMinPileupDistance));
500 }
501
502 //____________________________________________________________________
503 void
504 AliFMDEventInspector::DefineOutput(TList* dir)
505 {
506   // 
507   // Define the output histograms.  These are put in a sub list of the
508   // passed list.   The histograms are merged before the parent task calls 
509   // AliAnalysisTaskSE::Terminate 
510   // 
511   //   dir Directory to add to 
512   //
513   DGUARD(fDebug,1,"Define output from AliFMDEventInspector");
514   fList = new TList;
515   fList->SetName(GetName());
516   dir->Add(fList);
517 }
518
519 //____________________________________________________________________
520 UInt_t
521 AliFMDEventInspector::Process(const AliESDEvent* event, 
522                               UInt_t&            triggers,
523                               Bool_t&            lowFlux,
524                               UShort_t&          ivz, 
525                               Double_t&          vz,
526                               Double_t&          cent,
527                               UShort_t&          nClusters)
528 {
529   // 
530   // Process the event 
531   // 
532   // Parameters:
533   //   event     Input event 
534   //   triggers  On return, the triggers fired 
535   //   lowFlux   On return, true if the event is considered a low-flux 
536   //                  event (according to the setting of fLowFluxCut) 
537   //   ivz       On return, the found vertex bin (1-based).  A zero
538   //                  means outside of the defined vertex range
539   //   vz        On return, the z position of the interaction
540   //   cent      On return, the centrality - if not available < 0
541   // 
542   // Return:
543   //    0 (or kOk) on success, otherwise a bit mask of error codes 
544   //
545   DGUARD(fDebug,1,"Process event in AliFMDEventInspector");
546
547   // --- Check that we have an event ---------------------------------
548   if (!event) { 
549     AliWarning("No ESD event found for input event");
550     fHStatus->Fill(2);
551     return kNoEvent;
552   }
553
554   // --- Read trigger information from the ESD and store in AOD object
555   if (!ReadTriggers(*event, triggers, nClusters)) { 
556     if (fDebug > 2) {
557       AliWarning("Failed to read triggers from ESD"); }
558     fHStatus->Fill(3);
559     return kNoTriggers;
560   }
561
562   // --- Check if this is a high-flux event --------------------------
563   const AliMultiplicity* testmult = event->GetMultiplicity();
564   if (!testmult) {
565     if (fDebug > 3) {
566       AliWarning("No central multiplicity object found"); }
567   }
568   else 
569     lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
570
571   fHType->Fill(lowFlux ? 0 : 1);
572
573   // --- Process satellite event information is requested ------------
574   if (fUseDisplacedVertices) { 
575     if (!fDisplacedVertex.Process(event)) 
576       AliWarning("Failed to process satellite event");
577   }
578   
579   // --- Read centrality information 
580   cent          = -10;
581   UShort_t qual = 0;
582   if (!ReadCentrality(*event, cent, qual)) {
583     if (fDebug > 3) 
584       AliWarning("Failed to get centrality");
585   }
586   fHCent->Fill(cent);
587   if (qual == 0) fHCentVsQual->Fill(0., cent);
588   else { 
589     for (UShort_t i = 0; i < 4; i++) 
590       if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
591   }
592
593   // --- Get the vertex information ----------------------------------
594   Double_t vx = 0;
595   Double_t vy = 0;
596   vz          = 0;
597   
598   Bool_t vzOk = ReadVertex(*event, vz,vx,vy);
599
600   fHEventsTr->Fill(vz);
601   if (!vzOk) { 
602     if (fDebug > 3) {
603       AliWarning("Failed to read vertex from ESD"); }
604     fHStatus->Fill(6);
605     return kNoVertex;
606   }
607   fHEventsTrVtx->Fill(vz);
608   
609   // --- Get the vertex bin ------------------------------------------
610   ivz = fVtxAxis.FindBin(vz);
611   if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) { 
612     if (fDebug > 3) {
613       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
614                       vz, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); 
615     }
616     ivz = 0;
617     fHStatus->Fill(7);
618     return kBadVertex;
619   }
620   fHEventsAccepted->Fill(vz);
621   fHEventsAcceptedXY->Fill(vx,vy);
622   
623   // --- Check the FMD ESD data --------------------------------------
624   if (!event->GetFMDData()) { 
625     if (fDebug > 3) {
626       AliWarning("No FMD data found in ESD"); }
627     fHStatus->Fill(5);
628     return kNoFMD;
629   }
630
631   fHStatus->Fill(1);
632   return kOk;
633 }
634
635 //____________________________________________________________________
636 Bool_t
637 AliFMDEventInspector::ReadCentrality(const AliESDEvent& esd, 
638                                      Double_t& cent, 
639                                      UShort_t& qual) const
640 {
641   // 
642   // Read centrality from event 
643   // 
644   // Parameters:
645   //    esd  Event 
646   //    cent On return, the centrality or negative if not found
647   // 
648   // Return:
649   //    False on error, true otherwise 
650   //
651   DGUARD(fDebug,2,"Read the centrality in AliFMDEventInspector");
652
653   if(fUseDisplacedVertices) {
654     Double_t zvtx = fDisplacedVertex.GetVertexZ();
655     qual          = 1;
656     if(TMath::Abs(zvtx) < 999) {
657       cent = fDisplacedVertex.GetCentralityPercentile();
658       qual = 0;
659     }
660     return true;
661   }
662   
663   cent = -1;
664   qual = 0;
665   AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
666   if (!centObj)  return true;
667
668   cent = centObj->GetCentralityPercentile("V0M");  
669   qual = centObj->GetQuality();
670
671   return true;
672 }
673
674 //____________________________________________________________________
675 Bool_t
676 AliFMDEventInspector::ReadTriggers(const AliESDEvent& esd, UInt_t& triggers,
677                                    UShort_t& nClusters)
678 {
679   // 
680   // Read the trigger information from the ESD event 
681   // 
682   // Parameters:
683   //   esd        ESD event 
684   //   triggers   On return, contains the trigger bits 
685   // 
686   // Return:
687   //    @c true on success, @c false otherwise 
688   //
689   DGUARD(fDebug,2,"Read the triggers in AliFMDEventInspector");
690   triggers = 0;
691
692   // Get the analysis manager - should always be there 
693   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
694   DMSG(fDebug,10,"Got analysis manager %p", am);
695   if (!am) { 
696     AliWarning("No analysis manager defined!");
697     return kFALSE;
698   }
699
700   // Get the input handler - should always be there 
701   AliInputEventHandler* ih = 
702     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
703   DMSG(fDebug,10,"Got input handler %p", ih);
704   if (!ih) { 
705     AliWarning("No input handler");
706     return kFALSE;
707   }
708
709   // Check if this is a collision candidate (MB)
710   // Note, that we should use the value cached in the input 
711   // handler rather than calling IsCollisionCandiate directly 
712   // on the AliPhysicsSelection obejct.  If we called the latter
713   // then the AliPhysicsSelection object would overcount by a 
714   // factor of 2! :-(
715   Bool_t  offline  = ih->IsEventSelected() ;
716   Bool_t  fastonly = (ih->IsEventSelected() & AliVEvent::kFastOnly);
717   TString trigStr  = esd.GetFiredTriggerClasses();
718
719   if (fHWords) fHWords->Fill(trigStr.Data(), 1);
720   
721   if(fUseDisplacedVertices) {
722     DMSG(fDebug,3,"Using displaced vertex stuff");
723     if (TMath::Abs(fDisplacedVertex.GetVertexZ()) >= 999) offline = false;
724   }
725   
726   if (CheckFastPartition(fastonly))     offline = false;
727   if (offline && CheckCosmics(trigStr)) offline = false;
728
729   DMSG(fDebug,2,"Event is %striggered by off-line", offline ? "" : "NOT ");
730
731   if (offline) {
732     triggers |= AliAODForwardMult::kOffline;
733     triggers |= AliAODForwardMult::kInel;
734     if (!fHTriggers) { 
735       AliWarning("Histogram of triggers not defined - has init been called");
736       return false;
737     }
738     fHTriggers->Fill(kOffline+0.5);
739     
740     CheckINELGT0(esd, nClusters, triggers);
741   }
742   
743   CheckNSD(esd,triggers);
744   if (CheckPileup(esd, triggers)) fHTriggers->Fill(kPileUp+.5);
745   if (CheckEmpty(trigStr, triggers)) fHTriggers->Fill(kEmpty+.5);
746
747   CheckWords(esd, triggers);
748
749   // Now check - if we have a collision - for offline triggers and
750   // fill histogram.
751   if (triggers & AliAODForwardMult::kB) {
752     fHTriggers->Fill(kB+.5);
753     if (triggers & AliAODForwardMult::kInel) 
754       fHTriggers->Fill(kInel);
755     
756     if (triggers & AliAODForwardMult::kInelGt0)
757       fHTriggers->Fill(kInelGt0+.5);
758     
759     if (triggers & AliAODForwardMult::kNSD)
760       fHTriggers->Fill(kNSD+.5);
761
762     if (triggers & AliAODForwardMult::kV0AND)
763       fHTriggers->Fill(kV0AND+.5);
764   }
765   if (triggers & AliAODForwardMult::kA) fHTriggers->Fill(kA+.5);
766   if (triggers & AliAODForwardMult::kC) fHTriggers->Fill(kC+.5);
767   if (triggers & AliAODForwardMult::kE) fHTriggers->Fill(kE+.5);
768   
769   return kTRUE;
770 }
771
772 //____________________________________________________________________
773 Bool_t
774 AliFMDEventInspector::CheckFastPartition(bool fastonly) const
775 {
776   // For the 2.76 TeV p+p run, the FMD ran in the slow partition 
777   // so it received no triggers from the fast partition. Therefore
778   // the fast triggers are removed here but not for MC where all 
779   // triggers are fast.
780   if (TMath::Abs(fEnergy - 2750.) > 20) return false;
781   if (fCollisionSystem != AliForwardUtil::kPP) return false;
782   if (fastonly)
783     DMSG(fDebug,1,"Fast trigger in pp @ sqrt(s)=2.76TeV removed");
784
785   return fastonly;
786 }
787
788 //____________________________________________________________________
789 Bool_t
790 AliFMDEventInspector::CheckCosmics(const TString& trigStr) const
791 {
792   // MUON triggers are not strictly minimum bias (MB) so they are
793   // removed (HHD)
794   if(trigStr.Contains("CMUS1")) {
795     DMSG(fDebug,1,"Cosmic trigger ins't min-bias, removed");
796     return true;
797   }
798   return false;
799 }
800
801 //____________________________________________________________________
802 Bool_t
803 AliFMDEventInspector::CheckINELGT0(const AliESDEvent& esd, 
804                                    UShort_t& nClusters, 
805                                    UInt_t& triggers) const
806 {
807   nClusters = 0;
808
809   // If this is inel, see if we have a tracklet 
810   const AliMultiplicity* spdmult = esd.GetMultiplicity();
811   if (!spdmult) {
812     AliWarning("No SPD multiplicity");
813     return false;
814   }
815
816   // Check if we have one or more tracklets 
817   // in the range -1 < eta < 1 to set the INEL>0 
818   // trigger flag. 
819   // 
820   // Also count tracklets as a single cluster 
821   Int_t n = spdmult->GetNumberOfTracklets();
822   for (Int_t j = 0; j < n; j++) { 
823     if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
824       triggers |= AliAODForwardMult::kInelGt0;
825       nClusters++;
826     }
827   }
828   n = spdmult->GetNumberOfSingleClusters();
829   for (Int_t j = 0; j < n; j++) { 
830     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
831     if (TMath::Abs(eta) < 1) nClusters++;
832   }
833   if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
834
835   return triggers & AliAODForwardMult::kNClusterGt0;
836 }
837
838 //____________________________________________________________________
839 Bool_t
840 AliFMDEventInspector::CheckNSD(const AliESDEvent& esd, UInt_t& triggers) const
841 {
842   // Analyse some trigger stuff 
843   AliTriggerAnalysis ta;
844   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kV0AND)) {
845     triggers |= AliAODForwardMult::kV0AND;
846     if (fUseV0AND) 
847       triggers |= AliAODForwardMult::kNSD;
848   }
849   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kNSD1)) 
850     triggers |= AliAODForwardMult::kNSD;
851   return triggers & AliAODForwardMult::kNSD;
852 }
853 //____________________________________________________________________
854 Bool_t
855 AliFMDEventInspector::CheckPileup(const AliESDEvent& esd, 
856                                   UInt_t& triggers) const
857 {
858   // Check for multiple vertices (pile-up) with at least 3
859   // contributors and at least 0.8cm from the primary vertex
860   if(fCollisionSystem != AliForwardUtil::kPP) return false;
861
862   Bool_t pileup =  esd.IsPileupFromSPD(fMinPileupContrib,fMinPileupDistance);
863   if (pileup) triggers |= AliAODForwardMult::kPileUp;
864   return pileup;
865 }
866
867 //____________________________________________________________________
868 Bool_t
869 AliFMDEventInspector::CheckEmpty(const TString& trigStr, UInt_t& triggers) const
870 {
871   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
872     triggers |= AliAODForwardMult::kEmpty;
873     return true;
874   }
875   return false;
876 }
877 //____________________________________________________________________
878 Bool_t
879 AliFMDEventInspector::CheckWords(const AliESDEvent& esd, UInt_t& triggers) const
880 {
881   TObject* word = 0;
882   TIter nextColl(&fCollWords);
883   while ((word = nextColl())) {
884     DMSG(fDebug,10,"Checking if %s trigger %s is fired", 
885          word->GetTitle(), word->GetName());
886     if (!esd.IsTriggerClassFired(word->GetName())) continue;
887
888     TString beamSide = word->GetTitle();
889     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
890
891     if (!beamSide.EqualTo("B")) continue;
892     triggers |= AliAODForwardMult::kB;
893     break; // No more to do here
894   }
895   TIter nextBg(&fBgWords);
896   UInt_t all = (AliAODForwardMult::kA | 
897                 AliAODForwardMult::kC | 
898                 AliAODForwardMult::kE);
899   while ((word = nextBg())) {
900     DMSG(fDebug,10,"Checking if %s trigger %s is fired", 
901          word->GetTitle(), word->GetName());
902     if (!esd.IsTriggerClassFired(word->GetName())) continue;
903     
904     TString beamSide = word->GetTitle();
905     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
906
907     if (beamSide.Contains("A")) triggers |= AliAODForwardMult::kA;
908     if (beamSide.Contains("C")) triggers |= AliAODForwardMult::kC;
909     if (beamSide.Contains("E")) triggers |= AliAODForwardMult::kE;
910
911     if ((triggers & all) == all) break; // No more to do
912   }
913   return true;
914 }
915
916
917 //____________________________________________________________________
918 Bool_t
919 AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, 
920                                  Double_t& vz, 
921                                  Double_t& vx, 
922                                  Double_t& vy)
923 {
924   // 
925   // Read the vertex information from the ESD event 
926   // 
927   // Parameters:
928   //   esd  ESD event 
929   //   vz   On return, the vertex Z position 
930   // 
931   // Return:
932   //    @c true on success, @c false otherwise 
933   //
934   DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
935   vz = 0;
936   vx = 1024;
937   vy = 1024;
938   
939   if(fUseDisplacedVertices) {
940     Double_t zvtx = fDisplacedVertex.GetVertexZ();
941       
942     if(TMath::Abs(zvtx) < 999) {
943       vz = zvtx;
944       return true;
945     }
946     return false;
947   }
948
949   if(fUseFirstPhysicsVertex) return CheckPWGUDVertex(esd, vz, vx, vy);
950   
951   
952   return CheckVertex(esd, vz, vx, vy);
953 }
954
955 //____________________________________________________________________
956 Bool_t
957 AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd, 
958                                        Double_t& vz, 
959                                        Double_t& vx, 
960                                        Double_t& vy) const
961 {
962   // This is the code used by the 1st physics people 
963   const AliESDVertex* vertex    = esd.GetPrimaryVertex();
964   if (!vertex  || !vertex->GetStatus()) {
965     DMSG(fDebug,2,"No primary vertex (%p) or bad status %d", 
966          vertex, (vertex ? vertex->GetStatus() : -1));
967     return false;
968   }
969   const AliESDVertex* vertexSPD = esd.GetPrimaryVertexSPD();
970   if (!vertexSPD || !vertexSPD->GetStatus()) {
971     DMSG(fDebug,2,"No primary SPD vertex (%p) or bad status %d", 
972          vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1));
973     return false;
974   }
975     
976   // if vertex is from SPD vertexZ, require more stringent cuts 
977   if (vertex->IsFromVertexerZ()) { 
978     if (vertex->GetDispersion() > fMaxVzErr || 
979         vertex->GetZRes() > 1.25 * fMaxVzErr) {
980       DMSG(fDebug,2,"Dispersion %f > %f or resolution %f > %f",
981            vertex->GetDispersion(), fMaxVzErr,
982            vertex->GetZRes(), 1.25 * fMaxVzErr);
983       return false;
984     }
985   }
986   vz = vertex->GetZ();
987   
988   if(!vertex->IsFromVertexerZ()) {
989     vx = vertex->GetX();
990     vy = vertex->GetY();
991   }
992   return true;
993 }
994 //____________________________________________________________________
995 Bool_t
996 AliFMDEventInspector::CheckVertex(const AliESDEvent& esd, 
997                                   Double_t& vz, 
998                                   Double_t& vx, 
999                                   Double_t& vy) const
1000 {
1001   // Use standard SPD vertex (perhaps preferable for Pb+Pb)
1002   // Get the vertex 
1003   const AliESDVertex* vertex = esd.GetPrimaryVertexSPD();
1004   if (!vertex) { 
1005     if (fDebug > 2) {
1006       AliWarning("No SPD vertex found in ESD"); }
1007     return false;
1008   }
1009     
1010   // Check that enough tracklets contributed 
1011   if(vertex->GetNContributors() <= 0) {
1012     DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
1013          vertex->GetNContributors());
1014     vz = 0;
1015     return false;
1016   } 
1017   // Check that the uncertainty isn't too large 
1018   if (vertex->GetZRes() > fMaxVzErr) { 
1019     DMSG(fDebug,2,"Uncertaintity in Z of vertex is too large %f > %f", 
1020          vertex->GetZRes(), fMaxVzErr);
1021     return false;
1022   }
1023     
1024   // Get the z coordiante 
1025   vz = vertex->GetZ();
1026   const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
1027     
1028   if(!vertexXY->IsFromVertexerZ()) {
1029     vx = vertexXY->GetX();
1030     vy = vertexXY->GetY();
1031   }
1032   return true;
1033 }
1034
1035 //____________________________________________________________________
1036 Bool_t
1037 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
1038 {
1039   // 
1040   // Read the collision system, collision energy, and L3 field setting
1041   // from the ESD
1042   // 
1043   // Parameters:
1044   //   esd ESD to get information from 
1045   // 
1046   // Return:
1047   //    true on success, false 
1048   //
1049   // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
1050   //           esd->GetBeamType(), 2*esd->GetBeamEnergy(), 
1051   //           esd->GetMagneticField()));
1052   DGUARD(fDebug,2,"Read the run details in AliFMDEventInspector");
1053   const char* sys  = esd->GetBeamType();
1054   Float_t     cms  = 2 * esd->GetBeamEnergy();
1055   Float_t     fld  = esd->GetMagneticField();
1056   fCollisionSystem = AliForwardUtil::ParseCollisionSystem(sys);
1057   fEnergy          = AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem, 
1058                                                              cms);
1059   fField           = AliForwardUtil::ParseMagneticField(fld);
1060
1061   StoreInformation(esd->GetRunNumber());
1062   if (fCollisionSystem   == AliForwardUtil::kUnknown) { 
1063     AliWarningF("Unknown collision system: %s - please check", sys);
1064     return false;
1065   }
1066   if (fEnergy            <= 0) {
1067     AliWarningF("Unknown CMS energy: %f (%d) - please check", cms, fEnergy);
1068     return false;
1069   }
1070   if (TMath::Abs(fField) >  10) {
1071     AliWarningF("Unknown L3 field setting: %f (%d) - please check", fld,fField);
1072     return false;
1073   }
1074
1075   return true;
1076 }
1077
1078
1079 //____________________________________________________________________
1080 const Char_t*
1081 AliFMDEventInspector::CodeString(UInt_t code)
1082 {
1083   static TString s;
1084   s = "";
1085   if (code & kNoEvent)    s.Append("NOEVENT ");
1086   if (code & kNoTriggers) s.Append("NOTRIGGERS ");
1087   if (code & kNoSPD)      s.Append("NOSPD ");
1088   if (code & kNoFMD)      s.Append("NOFMD ");
1089   if (code & kNoVertex)   s.Append("NOVERTEX ");
1090   if (code & kBadVertex)  s.Append("BADVERTEX ");
1091   return s.Data();
1092 }
1093 //____________________________________________________________________
1094 void
1095 AliFMDEventInspector::Print(Option_t*) const
1096 {
1097   // 
1098   // Print information
1099   // 
1100   //   option Not used 
1101   //
1102   char ind[gROOT->GetDirLevel()+1];
1103   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1104   ind[gROOT->GetDirLevel()] = '\0';
1105   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
1106   sNN.Strip(TString::kBoth, '0');
1107   sNN.ReplaceAll("GeV", " GeV");
1108   TString field(AliForwardUtil::MagneticFieldString(fField));
1109   field.ReplaceAll("p",  "+");
1110   field.ReplaceAll("m",  "-");
1111   field.ReplaceAll("kG", " kG");
1112   
1113   std::cout << std::boolalpha 
1114             << ind << ClassName() << ": " << GetName() << '\n'
1115             << ind << " Vertex bins:            " << fVtxAxis.GetNbins() << '\n'
1116             << ind << " Vertex range:           [" << fVtxAxis.GetXmin() 
1117             << "," << fVtxAxis.GetXmax() << "]\n"
1118             << ind << " Low flux cut:           " << fLowFluxCut << '\n'
1119             << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
1120             << ind << " Min(nContrib_pileup):   " << fMinPileupContrib << '\n'
1121             << ind << " Min(v-pileup):          " << fMinPileupDistance << '\n'
1122             << ind << " System:                 " 
1123             << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
1124             << ind << " CMS energy per nucleon: " << sNN << '\n'
1125             << ind << " Field:                  " <<  field << '\n'
1126             << ind << " Satellite events:       " << fUseDisplacedVertices 
1127             << "\n" << std::noboolalpha;
1128   if (!fCentAxis) { std::cout << std::flush; return; }
1129   Int_t nBin = fCentAxis->GetNbins();
1130   std::cout << ind << " Centrality axis:        " << nBin << " bins"
1131             << std::flush;
1132   for (Int_t i = 0; i < nBin; i++) { 
1133     if ((i % 10) == 0) std::cout << '\n' << ind << "  ";
1134     std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
1135   }
1136   std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
1137 }
1138
1139   
1140 //
1141 // EOF
1142 //
1143