Mega commit of many changes to PWGLFforward
[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 #include "AliVVZERO.h"
42
43 //====================================================================
44 const char* AliFMDEventInspector::fgkFolderName = "fmdEventInspector";
45
46 //____________________________________________________________________
47 AliFMDEventInspector::AliFMDEventInspector()
48   : TNamed(),
49     fHEventsTr(0), 
50     fHEventsTrVtx(0),
51     fHEventsAccepted(0),
52     fHEventsAcceptedXY(0),
53     fHTriggers(0),
54     fHTriggerCorr(0),
55     fHType(0),
56     fHWords(0),
57     fHCent(0),
58     fHCentVsQual(0),
59     fHStatus(0),
60     fHVtxStatus(0),
61     fHTrgStatus(0),
62     fLowFluxCut(1000),
63     fMaxVzErr(0.2),
64     fList(0),
65     fEnergy(0),
66     fField(999), 
67     fCollisionSystem(kUnknown),
68     fDebug(0),
69     fCentAxis(0),
70     fVtxAxis(10,-10,10),
71     fUseFirstPhysicsVertex(false),
72     fUseV0AND(false),
73     fMinPileupContrib(3), 
74   fMinPileupDistance(0.8),
75   fUseDisplacedVertices(false),
76   fDisplacedVertex(),
77   fCollWords(),
78   fBgWords(),
79   fCentMethod("V0M"),
80   fMinCent(-1.0),
81   fMaxCent(-1.0),
82   fUsepA2012Vertex(false),
83   fRunNumber(0)
84 {
85   // 
86   // Constructor 
87   //
88   DGUARD(fDebug,1,"Default CTOR of AliFMDEventInspector");
89 }
90
91 //____________________________________________________________________
92 AliFMDEventInspector::AliFMDEventInspector(const char* name)
93   : TNamed(fgkFolderName, name),
94     fHEventsTr(0), 
95     fHEventsTrVtx(0), 
96     fHEventsAccepted(0),
97     fHEventsAcceptedXY(0),
98     fHTriggers(0),
99     fHTriggerCorr(0),
100     fHType(0),
101     fHWords(0),
102     fHCent(0),
103     fHCentVsQual(0),
104     fHStatus(0),
105     fHVtxStatus(0),
106     fHTrgStatus(0),
107     fLowFluxCut(1000),
108     fMaxVzErr(0.2),
109     fList(0),
110     fEnergy(0),
111     fField(999), 
112     fCollisionSystem(kUnknown),
113     fDebug(0),
114     fCentAxis(0),
115     fVtxAxis(10,-10,10),
116     fUseFirstPhysicsVertex(false),
117     fUseV0AND(false),
118     fMinPileupContrib(3), 
119     fMinPileupDistance(0.8),
120     fUseDisplacedVertices(false),
121   fDisplacedVertex(),
122   fCollWords(),
123   fBgWords(),
124   fCentMethod("V0M"),
125   fMinCent(-1.0),
126   fMaxCent(-1.0),
127   fUsepA2012Vertex(false),
128   fRunNumber(0)
129 {
130   // 
131   // Constructor 
132   // 
133   // Parameters:
134   //   name Name of object
135   //
136   DGUARD(fDebug,1,"Named CTOR of AliFMDEventInspector: %s", name);
137 }
138
139 //____________________________________________________________________
140 AliFMDEventInspector::AliFMDEventInspector(const AliFMDEventInspector& o)
141   : TNamed(o), 
142     fHEventsTr(o.fHEventsTr), 
143     fHEventsTrVtx(o.fHEventsTrVtx), 
144     fHEventsAccepted(o.fHEventsAccepted),
145     fHEventsAcceptedXY(o.fHEventsAcceptedXY),
146     fHTriggers(o.fHTriggers),
147     fHTriggerCorr(o.fHTriggerCorr),
148     fHType(o.fHType),
149     fHWords(o.fHWords),
150     fHCent(o.fHCent),
151     fHCentVsQual(o.fHCentVsQual),
152     fHStatus(o.fHStatus),
153     fHVtxStatus(o.fHVtxStatus),
154     fHTrgStatus(o.fHTrgStatus),
155     fLowFluxCut(o.fLowFluxCut),
156     fMaxVzErr(o.fMaxVzErr),
157     fList(o.fList),
158     fEnergy(o.fEnergy),
159     fField(o.fField), 
160     fCollisionSystem(o.fCollisionSystem),
161     fDebug(0),
162     fCentAxis(0),
163     fVtxAxis(o.fVtxAxis),
164     fUseFirstPhysicsVertex(o.fUseFirstPhysicsVertex),
165     fUseV0AND(o.fUseV0AND),
166     fMinPileupContrib(o.fMinPileupContrib), 
167     fMinPileupDistance(o.fMinPileupDistance),
168     fUseDisplacedVertices(o.fUseDisplacedVertices),
169     fDisplacedVertex(o.fDisplacedVertex),
170   fCollWords(),
171   fBgWords(),
172   fCentMethod(o.fCentMethod),
173   fMinCent(o.fMinCent),
174   fMaxCent(o.fMaxCent),
175   fUsepA2012Vertex(o.fUsepA2012Vertex),
176   fRunNumber(o.fRunNumber)
177         
178 {
179   // 
180   // Copy constructor 
181   // 
182   // Parameters:
183   //   o Object to copy from 
184   //
185   DGUARD(fDebug,1,"Copy CTOR of AliFMDEventInspector");
186 }
187
188 //____________________________________________________________________
189 AliFMDEventInspector::~AliFMDEventInspector()
190 {
191   // 
192   // Destructor 
193   //
194   DGUARD(fDebug,1,"DTOR of AliFMDEventInspector");
195   // if (fList)         delete fList;
196 }
197 //____________________________________________________________________
198 AliFMDEventInspector&
199 AliFMDEventInspector::operator=(const AliFMDEventInspector& o)
200 {
201   // 
202   // Assignement operator
203   // 
204   // Parameters:
205   //   o Object to assign from 
206   // 
207   // Return:
208   //    Reference to this object
209   //
210   DGUARD(fDebug,3,"Assignment of AliFMDEventInspector");
211   if (&o == this) return *this; 
212   TNamed::operator=(o);
213   fHEventsTr         = o.fHEventsTr;
214   fHEventsTrVtx      = o.fHEventsTrVtx;
215   fHEventsAccepted   = o.fHEventsAccepted;
216   fHEventsAcceptedXY = o.fHEventsAcceptedXY;
217   fHTriggers         = o.fHTriggers;
218   fHTriggerCorr      = o.fHTriggerCorr;
219   fHType             = o.fHType;
220   fHWords            = o.fHWords;
221   fHCent             = o.fHCent;
222   fHCentVsQual       = o.fHCentVsQual;
223   fHStatus           = o.fHStatus;
224   fHVtxStatus        = o.fHVtxStatus;
225   fHTrgStatus        = o.fHTrgStatus;
226   fLowFluxCut        = o.fLowFluxCut;
227   fMaxVzErr          = o.fMaxVzErr;
228   fDebug             = o.fDebug;
229   fList              = (o.fList ? new TList : 0);
230   fEnergy            = o.fEnergy;
231   fField             = o.fField;
232   fCollisionSystem   = o.fCollisionSystem;
233   fVtxAxis.Set(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), 
234                o.fVtxAxis.GetXmax());
235   
236   fUseFirstPhysicsVertex = o.fUseFirstPhysicsVertex;
237   fUseV0AND              = o.fUseV0AND;
238   fMinPileupContrib      = o.fMinPileupContrib;
239   fMinPileupDistance     = o.fMinPileupDistance;
240   fUseDisplacedVertices  = o.fUseDisplacedVertices;
241   fDisplacedVertex       = o.fDisplacedVertex;
242   fCentMethod            = o.fCentMethod;
243   fMinCent               = o.fMinCent;
244   fMaxCent               = o.fMaxCent; 
245   fUsepA2012Vertex       = o.fUsepA2012Vertex;
246   fRunNumber             = o.fRunNumber;
247
248   if (fList) { 
249     fList->SetName(GetName());
250     if (fHEventsTr)    fList->Add(fHEventsTr);
251     if (fHEventsTrVtx) fList->Add(fHEventsTrVtx);
252     if (fHTriggers)    fList->Add(fHTriggers);
253     if (fHTriggerCorr) fList->Add(fHTriggerCorr);
254     if (fHType)        fList->Add(fHType);
255     if (fHWords)       fList->Add(fHWords);
256     if (fHCent)        fList->Add(fHCent);
257     if (fHCentVsQual)  fList->Add(fHCentVsQual);
258     if (fHStatus)      fList->Add(fHStatus);
259     if (fHVtxStatus)   fList->Add(fHVtxStatus);
260     if (fHTrgStatus)   fList->Add(fHTrgStatus);
261   }
262   return *this;
263 }
264
265 //____________________________________________________________________
266 void 
267 AliFMDEventInspector::SetCentralityMethod(ECentMethod m)
268 {
269   switch (m) { 
270   case kV0Multiplicity: fCentMethod = "VOM"; break; // VZERO multiplicity 
271   case kV0Amplitude:    fCentMethod = "V0A"; break; // VZERO amplitude    
272   case kV0Charge:       fCentMethod = "V0C"; break; // VZERO charge          
273   case kFMDRough:       fCentMethod = "FMD"; break; // FMD scaled energy l
274   case kNTracks:        fCentMethod = "TRK"; break; // Number of tracks   
275   case kLTracks:        fCentMethod = "TKL"; break; // Number of tracks   
276   case kCL0:            fCentMethod = "CL0"; break; //               
277   case kCL1:            fCentMethod = "CL1"; break; //               
278   case kCND:            fCentMethod = "CND"; break; //               
279   case kNParticles:     fCentMethod = "NPA"; break; // Neutral particles  
280   case kNeutrons:       fCentMethod = "ZNA"; break; // ZDC neutron amplitu
281   case kV0vsFMD:        fCentMethod = "V0MvsFMD"; break; // VZERO versus FMD   
282   case kV0vsNTracks:    fCentMethod = "TKLvsVOM"; break; // Tracks versus VZERO
283   case kZEMvsZDC:       fCentMethod = "ZEMvsZDC"; break; // ZDC              
284   default:              fCentMethod = "V0M"; break;
285   }
286 }
287
288 //____________________________________________________________________
289 void 
290 AliFMDEventInspector::SetMinCentrality(Double_t minCent)
291 {
292   AliWarning("\n"
293              "*******************************************************\n"
294              "* Setting centrality cuts in this stage is deprecated *\n"
295              "*******************************************************");
296   fMinCent = minCent;
297 }
298 //____________________________________________________________________
299 void 
300 AliFMDEventInspector::SetMaxCentrality(Double_t maxCent)
301 {
302   AliWarning("\n"
303              "*******************************************************\n"
304              "* Setting centrality cuts in this stage is deprecated *\n"
305              "*******************************************************");
306   fMaxCent = maxCent;
307 }
308
309 //____________________________________________________________________
310 Bool_t
311 AliFMDEventInspector::FetchHistograms(const TList* d, 
312                                       TH1I*& hEventsTr, 
313                                       TH1I*& hEventsTrVtx, 
314                                       TH1I*& hEventsAcc,
315                                       TH1I*& hTriggers) const
316 {
317   // 
318   // Fetch our histograms from the passed list 
319   // 
320   // Parameters:
321   //   d             Input
322   //   hEventsTr     On return, pointer to histogram, or null
323   //   hEventsTrVtx  On return, pointer to histogram, or null
324   //   hTriggers     On return, pointer to histogram, or null
325   // 
326   // Return:
327   //    true on success, false otherwise 
328   //
329   DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector");
330   hEventsTr    = 0;
331   hEventsTrVtx = 0;
332   hEventsAcc   = 0;
333   hTriggers    = 0;  
334   TList* dd    = dynamic_cast<TList*>(d->FindObject(GetName()));
335   if (!dd) return kFALSE;
336   
337   hEventsTr    = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
338   hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
339   hEventsAcc   = dynamic_cast<TH1I*>(dd->FindObject("nEventsAccepted"));
340   hTriggers    = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
341
342   if (!hEventsTr    || 
343       !hEventsTrVtx || 
344       !hEventsAcc   ||
345       !hTriggers) return kFALSE;
346   return kTRUE;
347 }
348 //____________________________________________________________________
349 void
350 AliFMDEventInspector::CacheConfiguredTriggerClasses(TList& cache, 
351                                                     const TList* classes,
352                                                     AliOADBPhysicsSelection* o)
353 {
354   TIter nextClass(classes);
355   TObjString* trigClass = 0;
356   // Loop over all trigger classes.  Trigger classes have the format 
357   //
358   //   class          := positive_words SPACE(s) negative_words 
359   //   positive_words := 
360   //                  |  '+' words
361   //   negative_words := 
362   //                  |  '-' words
363   //   words          := word 
364   //                  |  word ',' words 
365   //   
366   while ((trigClass = static_cast<TObjString*>(nextClass()))) {
367     // Tokenize on space to get positive and negative parts 
368     TString     side   = o->GetBeamSide(trigClass->String());
369     TObjArray*  parts  = trigClass->String().Tokenize(" ");
370     TObjString* part   = 0;
371     TIter       nextPart(parts);
372     while ((part = static_cast<TObjString*>(nextPart()))) {
373       // We only care about the positive ones 
374       if (part->GetName()[0] != '+') continue;
375       part->String().Remove(0,1);
376         
377       // Tokenize on a comma to get the words 
378       TObjArray*  words = part->String().Tokenize(",");
379       TObjString* word  = 0;
380       TIter       nextWord(words);
381       while ((word = static_cast<TObjString*>(nextWord()))) {
382         TNamed* store = new TNamed(word->String(), side);
383         cache.Add(store);
384         DMSG(fDebug,3,"Caching %s trigger word %s", 
385              store->GetTitle(), store->GetName());
386       } // while (word)
387       delete words;
388     }
389     delete parts;
390   }
391 }
392
393 //____________________________________________________________________
394 void
395 AliFMDEventInspector::SetupForData(const TAxis& vtxAxis)
396 {
397   // 
398   // Initialize the object - this is called on the first seen event. 
399   // 
400   // Parameters:
401   //   vtxAxis Vertex axis in use 
402   //
403   DGUARD(fDebug,1,"Initialize in AliFMDEventInspector");
404
405   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
406
407   // Get the input handler - should always be there 
408   AliInputEventHandler* ih = 
409     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
410   if (!ih) { 
411     AliWarning("No input handler");
412     return;
413   }
414   // Get the physics selection - should always be there 
415   AliPhysicsSelection* ps = 
416     static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
417   if (!ps) {
418     AliWarning("No physics selection");
419     return;
420   }
421   // Get the configured triggers
422   AliOADBPhysicsSelection* oadb = 
423     const_cast<AliOADBPhysicsSelection*>(ps->GetOADBPhysicsSelection());
424   if (!oadb) {
425     AliWarning("No OADB physics selection object");
426     return;
427   }
428   // Get the configured trigger words from the physics selection 
429   const TList* collTriggClasses = ps->GetCollisionTriggerClasses();
430   const TList* bgTriggClasses   = ps->GetBGTriggerClasses();
431   if (!collTriggClasses) { 
432     AliWarning("No configured collision trigger classes");
433     return;
434   }
435   if (!bgTriggClasses) { 
436     AliWarning("No configured background trigger classes");
437     return;
438   }
439   CacheConfiguredTriggerClasses(fCollWords, collTriggClasses, oadb);
440   CacheConfiguredTriggerClasses(fBgWords,   bgTriggClasses,   oadb);
441   // fCollWords.ls();
442   // fBgWords.ls();
443   
444   
445   TArrayD limits;
446   if ((fMinCent < 0 && fMaxCent < 0) || fMaxCent <= fMinCent) {
447     // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
448     // ----- 92 number --------- ---- 1 ---
449     limits.Set(93);
450     for (Int_t i = 0; i < 92; i++) limits[i] = -1.5 + i;
451     limits[92] = 100.5;
452   }
453   else {
454     Int_t n = fMaxCent-fMinCent+2;
455     limits.Set(n);
456     for (Int_t i = 0; i < n; i++) { 
457       limits[i] = fMinCent + i - .5;
458     }
459   }
460       
461   fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
462   
463   fCentAxis  = new TAxis(limits.GetSize()-1, limits.GetArray());
464   fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger", 
465                         4*vtxAxis.GetNbins(), 
466                         2*vtxAxis.GetXmin(), 
467                         2*vtxAxis.GetXmax());
468   fHEventsTr->SetXTitle("v_{z} [cm]");
469   fHEventsTr->SetYTitle("# of events");
470   fHEventsTr->SetFillColor(kRed+1);
471   fHEventsTr->SetFillStyle(3001);
472   fHEventsTr->SetDirectory(0);
473   // fHEventsTr->Sumw2();
474   fList->Add(fHEventsTr);
475
476   fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
477   fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex"); 
478   fHEventsTrVtx->SetFillColor(kBlue+1);
479   fHEventsTrVtx->SetDirectory(0);
480   // fHEventsTrVtx->Sumw2();
481   fList->Add(fHEventsTrVtx);
482
483   fHEventsAccepted = new TH1I("nEventsAccepted", 
484                               "Number of events  w/trigger and vertex in range",
485                               2*vtxAxis.GetNbins(), 
486                               2*vtxAxis.GetXmin(), 
487                               2*vtxAxis.GetXmax());
488   fHEventsAccepted->SetXTitle("v_{z} [cm]");
489   fHEventsAccepted->SetYTitle("# of events");
490   fHEventsAccepted->SetFillColor(kGreen+1);
491   fHEventsAccepted->SetFillStyle(3001);
492   fHEventsAccepted->SetDirectory(0);
493   // fHEventsAccepted->Sumw2();
494   fList->Add(fHEventsAccepted);
495                               
496   fHEventsAcceptedXY = new TH2D("nEventsAcceptedXY", 
497                               "XY vertex w/trigger and Z vertex in range",
498                                 1000,-1,1,1000,-1,1);
499   
500   fHEventsAcceptedXY->SetXTitle("v_{x} [cm]");
501   fHEventsAcceptedXY->SetYTitle("v_{y} [cm]");
502   fHEventsAcceptedXY->SetDirectory(0);
503   // fHEventsAccepted->Sumw2();
504   fList->Add(fHEventsAcceptedXY);
505
506   
507   fHTriggers = new TH1I("triggers", "Triggers", kOffline+1, 0, kOffline+1);
508   fHTriggers->SetFillColor(kRed+1);
509   fHTriggers->SetFillStyle(3001);
510   fHTriggers->SetStats(0);
511   fHTriggers->SetDirectory(0);
512
513   fHTriggerCorr = new TH2I("triggerCorr", "Trigger correlation", 
514                            kOffline+1, 0, kOffline+1, 
515                            kOffline+1, 0, kOffline+1);
516   fHTriggerCorr->SetStats(0);
517   fHTriggerCorr->SetDirectory(0);
518   fHTriggerCorr->SetXTitle("Requirement");
519   fHTriggerCorr->SetYTitle("Companion");
520
521   Int_t binNum[] = { kInel   +1,
522                      kInelGt0+1,
523                      kNSD    +1,
524                      kV0AND  +1,
525                      kEmpty  +1,
526                      kA      +1,
527                      kB      +1,
528                      kC      +1,
529                      kE      +1,
530                      kPileUp +1,
531                      kMCNSD  +1,
532                      kSatellite+1,
533                      kOffline+1 };
534   const char* binLbl[] = { "INEL",       
535                            "INEL>0",
536                            "NSD",        
537                            "VOAND",
538                            "Empty",      
539                            "A",  
540                            "B",  
541                            "C",  
542                            "E",  
543                            "Pileup",
544                            "NSD_{MC}", 
545                            "Satellite",
546                            "Offline" };
547   for (Int_t i = 0; i < kOffline+1; i++) {
548     fHTriggers->GetXaxis()->SetBinLabel(binNum[i], binLbl[i]);
549     fHTriggerCorr->GetXaxis()->SetBinLabel(binNum[i], binLbl[i]);
550     fHTriggerCorr->GetYaxis()->SetBinLabel(binNum[i], binLbl[i]);
551   }
552   fList->Add(fHTriggers);
553   fList->Add(fHTriggerCorr);
554   
555
556   fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)", 
557                                  fLowFluxCut), 2, -.5, 1.5);
558   fHType->SetFillColor(kRed+1);
559   fHType->SetFillStyle(3001);
560   fHType->SetStats(0);
561   fHType->SetDirectory(0);
562   fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
563   fHType->GetXaxis()->SetBinLabel(2,"High-flux");
564   fList->Add(fHType);
565
566 #if 0 
567   // This histogram disabled as it causes problems in the merge 
568   fHWords = new TH1I("words", "Trigger words seen", 1, 0, 0); 
569   fHWords->SetFillColor(kBlue+1);
570   fHWords->SetFillStyle(3001);
571   fHWords->SetStats(0);
572   fHWords->SetDirectory(0);
573   fHWords->SetBit(TH1::kCanRebin);
574   fList->Add(fHWords);
575 #endif
576
577   fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
578   fHCent->SetFillColor(kBlue+1);
579   fHCent->SetFillStyle(3001);
580   fHCent->SetStats(0);
581   fHCent->SetDirectory(0);
582   fHCent->SetXTitle("Centrality [%]");
583   fHCent->SetYTitle("Events");
584   fList->Add(fHCent);
585
586   fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality", 
587                           5, 0, 5, limits.GetSize()-1, limits.GetArray());
588   fHCentVsQual->SetXTitle("Quality");
589   fHCentVsQual->SetYTitle("Centrality [%]");
590   fHCentVsQual->SetZTitle("Events");
591   fHCentVsQual->GetXaxis()->SetBinLabel(1, "OK");
592   fHCentVsQual->GetXaxis()->SetBinLabel(2, "Outside v_{z} cut");
593   fHCentVsQual->GetXaxis()->SetBinLabel(3, "V0 vs SPD outlier");
594   fHCentVsQual->GetXaxis()->SetBinLabel(4, "V0 vs TPC outlier");
595   fHCentVsQual->GetXaxis()->SetBinLabel(5, "V0 vs ZDC outlier");
596   fHCentVsQual->SetDirectory(0);
597   fList->Add(fHCentVsQual);
598
599   fHStatus = new TH1I("status", "Status", 7, 1, 8);
600   fHStatus->SetFillColor(kBlue+1);
601   fHStatus->SetFillStyle(3001);
602   fHStatus->SetStats(0);
603   fHStatus->SetDirectory(0);
604   TAxis* xAxis = fHStatus->GetXaxis();
605   xAxis->SetBinLabel(1, "OK");
606   xAxis->SetBinLabel(2, "No event");
607   xAxis->SetBinLabel(3, "No triggers");
608   xAxis->SetBinLabel(4, "No SPD");
609   xAxis->SetBinLabel(5, "No FMD");
610   xAxis->SetBinLabel(6, "No vertex");
611   xAxis->SetBinLabel(7, "Bad vertex");
612   fList->Add(fHStatus);
613
614   fHVtxStatus = new TH1I("vtxStatus","Vertex Status",
615                          kNotVtxZ,kVtxOK,kNotVtxZ+1);
616   fHVtxStatus->SetFillColor(kGreen+1);
617   fHVtxStatus->SetFillStyle(3001);
618   fHVtxStatus->SetStats(0);
619   fHVtxStatus->SetDirectory(0);
620   xAxis = fHVtxStatus->GetXaxis();
621   xAxis->SetBinLabel(kVtxOK,      "OK");
622   xAxis->SetBinLabel(kNoVtx,      "None/bad status");
623   xAxis->SetBinLabel(kNoSPDVtx,   "No SPD/bad status");
624   xAxis->SetBinLabel(kFewContrib, "N_{contrib} <= 0");
625   xAxis->SetBinLabel(kUncertain,  Form("#delta z > %4.2f", fMaxVzErr));
626   xAxis->SetBinLabel(kNotVtxZ,    "Not Z vertexer");
627   fList->Add(fHVtxStatus);
628
629   fHTrgStatus = new TH1I("trgStatus", "Trigger Status", 
630                          kOther, kNoTrgWords, kOther+1);
631   fHTrgStatus->SetFillColor(kMagenta+1);
632   fHTrgStatus->SetFillStyle(3001);
633   fHTrgStatus->SetStats(0);
634   fHTrgStatus->SetDirectory(0);
635   xAxis = fHTrgStatus->GetXaxis();
636   xAxis->SetBinLabel(kNoTrgWords,       "No words");
637   xAxis->SetBinLabel(kPP2760Fast,       "FAST in pp@#sqrt{s}=2.76TeV"); 
638   xAxis->SetBinLabel(kMUON,             "Muon trigger");
639   xAxis->SetBinLabel(kTriggered,        "Triggered");
640   xAxis->SetBinLabel(kMinBias,          "CINT1 (V0A||V0C||FASTOR)");
641   xAxis->SetBinLabel(kMinBiasNoSPD,     "CINT5 (V0A||V0C)");
642   xAxis->SetBinLabel(kV0AndTrg,         "CINT7 (V0A&&V0C)");
643   xAxis->SetBinLabel(kHighMult,         "N>>0");
644   xAxis->SetBinLabel(kCentral,          "Central"); 
645   xAxis->SetBinLabel(kSemiCentral,      "Semi-central"); 
646   xAxis->SetBinLabel(kDiffractive,      "Diffractive");
647   xAxis->SetBinLabel(kUser,             "User");
648   xAxis->SetBinLabel(kOther,            "Other");
649   fList->Add(fHTrgStatus);
650
651   if (fUseDisplacedVertices) fDisplacedVertex.SetupForData(fList);
652 }
653
654 //____________________________________________________________________
655 void
656 AliFMDEventInspector::StoreInformation()
657 {
658   // Write TNamed objects to output list containing information about
659   // the running conditions 
660   DGUARD(fDebug,2,"Store information from AliFMDEventInspector");
661   if (!fList) return;
662
663   
664   fList->Add(AliForwardUtil::MakeParameter("sys", fCollisionSystem));
665   fList->Add(AliForwardUtil::MakeParameter("sNN", fEnergy));
666   fList->Add(AliForwardUtil::MakeParameter("field", fField));
667   fList->Add(AliForwardUtil::MakeParameter("runNo", fRunNumber));
668   fList->Add(AliForwardUtil::MakeParameter("lowFlux", fLowFluxCut));
669   fList->Add(AliForwardUtil::MakeParameter("fpVtx",fUseFirstPhysicsVertex));
670   fList->Add(AliForwardUtil::MakeParameter("v0and",fUseV0AND));
671   fList->Add(AliForwardUtil::MakeParameter("nPileUp", fMinPileupContrib));
672   fList->Add(AliForwardUtil::MakeParameter("dPileup", fMinPileupDistance));
673   fList->Add(AliForwardUtil::MakeParameter("satellite", fUseDisplacedVertices));
674   fList->Add(AliForwardUtil::MakeParameter("alirootRev", 
675                                            AliForwardUtil::AliROOTRevision()));
676   fList->Add(AliForwardUtil::MakeParameter("alirootBranch", 
677                                            AliForwardUtil::AliROOTBranch()));
678
679 }
680
681 //____________________________________________________________________
682 void
683 AliFMDEventInspector::CreateOutputObjects(TList* dir)
684 {
685   // 
686   // Define the output histograms.  These are put in a sub list of the
687   // passed list.   The histograms are merged before the parent task calls 
688   // AliAnalysisTaskSE::Terminate 
689   // 
690   //   dir Directory to add to 
691   //
692   DGUARD(fDebug,1,"Define output from AliFMDEventInspector");
693   fList = new TList;
694   fList->SetName(GetName());
695   fList->SetOwner();
696   dir->Add(fList);
697 }
698
699 //____________________________________________________________________
700 UInt_t
701 AliFMDEventInspector::Process(const AliESDEvent* event, 
702                               UInt_t&            triggers,
703                               Bool_t&            lowFlux,
704                               UShort_t&          ivz, 
705                               TVector3&          ip,
706                               Double_t&          cent,
707                               UShort_t&          nClusters)
708 {
709   // 
710   // Process the event 
711   // 
712   // Parameters:
713   //   event     Input event 
714   //   triggers  On return, the triggers fired 
715   //   lowFlux   On return, true if the event is considered a low-flux 
716   //                  event (according to the setting of fLowFluxCut) 
717   //   ivz       On return, the found vertex bin (1-based).  A zero
718   //                  means outside of the defined vertex range
719   //   vz        On return, the z position of the interaction
720   //   cent      On return, the centrality - if not available < 0
721   // 
722   // Return:
723   //    0 (or kOk) on success, otherwise a bit mask of error codes 
724   //
725   DGUARD(fDebug,1,"Process event in AliFMDEventInspector"); 
726   // --- Check that we have an event ---------------------------------
727   if (!event) { 
728     AliWarning("No ESD event found for input event");
729     fHStatus->Fill(2);
730     return kNoEvent;
731   }
732
733   // --- Read trigger information from the ESD and store in AOD object
734   if (!ReadTriggers(*event, triggers, nClusters)) { 
735     if (fDebug > 2) {
736       AliWarning("Failed to read triggers from ESD"); }
737     fHStatus->Fill(3);
738     return kNoTriggers;
739   }
740
741   // --- Check if this is a high-flux event --------------------------
742   const AliMultiplicity* testmult = event->GetMultiplicity();
743   if (!testmult) {
744     if (fDebug > 3) {
745       AliWarning("No central multiplicity object found"); }
746   }
747   else 
748     lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
749
750   fHType->Fill(lowFlux ? 0 : 1);
751
752   // --- Process satellite event information is requested ------------
753   if (fUseDisplacedVertices) { 
754     if (!fDisplacedVertex.Process(event)) 
755       AliWarning("Failed to process satellite event");
756   }
757   
758  
759   // --- Get the interaction point -----------------------------------
760   Bool_t vzOk = ReadVertex(*event, ip);
761   fHEventsTr->Fill(ip.Z());
762   if (!vzOk) { 
763     if (fDebug > 3) {
764       AliWarning("Failed to read vertex from ESD"); }
765     fHStatus->Fill(6);
766     return kNoVertex;
767   }
768
769    // --- Read centrality information 
770   cent          = -10;
771   UShort_t qual = 0;
772   if (!ReadCentrality(*event, cent, qual)) {
773     if (fDebug > 3) 
774       AliWarning("Failed to get centrality");
775   }
776   // --- check centrality cut
777  
778   if(fMinCent > -0.0001 && cent < fMinCent) return  kNoEvent; 
779   if(fMaxCent > -0.0001 && cent > fMaxCent) return  kNoEvent; 
780   fHCent->Fill(cent);
781   if (qual == 0) fHCentVsQual->Fill(0., cent);
782   else { 
783     for (UShort_t i = 0; i < 4; i++) 
784       if (qual & (1 << i)) fHCentVsQual->Fill(Double_t(i+1), cent);
785   }
786         
787
788   fHEventsTrVtx->Fill(ip.Z());
789   
790   // --- Get the vertex bin ------------------------------------------
791   ivz = fVtxAxis.FindBin(ip.Z());
792   if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) { 
793     if (fDebug > 3) {
794       AliWarning(Form("Vertex @ %f outside of range [%f,%f]", 
795                       ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); 
796     }
797     ivz = 0;
798     fHStatus->Fill(7);
799     return kBadVertex;
800   }
801   fHEventsAccepted->Fill(ip.Z());
802   fHEventsAcceptedXY->Fill(ip.X(),ip.Y());
803   
804   // --- Check the FMD ESD data --------------------------------------
805   if (!event->GetFMDData()) { 
806     if (fDebug > 3) {
807       AliWarning("No FMD data found in ESD"); }
808     fHStatus->Fill(5);
809     return kNoFMD;
810   }
811
812   fHStatus->Fill(1);
813   return kOk;
814 }
815
816 //____________________________________________________________________
817 Bool_t
818 AliFMDEventInspector::ReadCentrality(const AliESDEvent& esd, 
819                                      Double_t& cent, 
820                                      UShort_t& qual) const
821 {
822   // 
823   // Read centrality from event 
824   // 
825   // Parameters:
826   //    esd  Event 
827   //    cent On return, the centrality or negative if not found
828   // 
829   // Return:
830   //    False on error, true otherwise 
831   //
832   DGUARD(fDebug,2,"Read the centrality in AliFMDEventInspector");
833
834   cent = -1;
835   qual = 1;
836   AliCentrality* centObj = const_cast<AliESDEvent&>(esd).GetCentrality();
837   if (centObj) {
838     cent = centObj->GetCentralityPercentile(fCentMethod);  
839     qual = centObj->GetQuality();
840   }
841
842   if (qual > 0 && fUseDisplacedVertices && fDisplacedVertex.IsSatellite()) {
843     cent = fDisplacedVertex.GetCentralityPercentile();
844     qual = 0;
845   }
846
847   return true;
848 }
849
850 //____________________________________________________________________
851 Bool_t
852 AliFMDEventInspector::CheckpAExtraV0(const AliESDEvent& esd) const
853 {
854   if (fCollisionSystem != AliForwardUtil::kPPb) return true;
855
856    AliVVZERO* esdV0 = esd.GetVZEROData();
857    if ((esdV0->GetV0ADecision()!=1) || (esdV0->GetV0CDecision()!=1)) 
858      return false;
859    return true;
860 }
861
862 //____________________________________________________________________
863 Bool_t
864 AliFMDEventInspector::ReadTriggers(const AliESDEvent& esd, UInt_t& triggers,
865                                    UShort_t& nClusters)
866 {
867   // 
868   // Read the trigger information from the ESD event 
869   // 
870   // Parameters:
871   //   esd        ESD event 
872   //   triggers   On return, contains the trigger bits 
873   // 
874   // Return:
875   //    @c true on success, @c false otherwise 
876   //
877   DGUARD(fDebug,2,"Read the triggers in AliFMDEventInspector");
878   triggers = 0;
879
880   // Get the analysis manager - should always be there 
881   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
882   DMSG(fDebug,10,"Got analysis manager %p", am);
883   if (!am) { 
884     AliWarning("No analysis manager defined!");
885     return kFALSE;
886   }
887
888   // Get the input handler - should always be there 
889   AliInputEventHandler* ih = 
890     static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
891   DMSG(fDebug,10,"Got input handler %p", ih);
892   if (!ih) { 
893     AliWarning("No input handler");
894     return kFALSE;
895   }
896
897   // Check if this is a collision candidate (MB)
898   ///
899   // Historic remark: Note, that we should use the value cached in the
900   //   input handler rather than calling IsCollisionCandiate directly
901   //   on the AliPhysicsSelection obejct.  If we called the latter
902   //   then the AliPhysicsSelection object would overcount by a factor
903   //   of 2! :-(
904   UInt_t  trgMask  = ih->IsEventSelected();
905   Bool_t  offline  = trgMask;
906   Bool_t  fastonly = (trgMask & AliVEvent::kFastOnly);
907   TString trigStr  = esd.GetFiredTriggerClasses();
908
909   if (trigStr.IsNull()) fHTrgStatus->Fill(kNoTrgWords);
910   if (fHWords) fHWords->Fill(trigStr.Data(), 1);
911   
912   if(fUseDisplacedVertices) {
913     DMSG(fDebug,3,"Using displaced vertex stuff");
914     // if (TMath::Abs(fDisplacedVertex.GetVertexZ()) >= 999) offline = false;
915     if (fDisplacedVertex.IsSatellite()) 
916       triggers |= AliAODForwardMult::kSatellite;
917   }
918   
919   if (CheckFastPartition(fastonly)) {
920     fHTrgStatus->Fill(kPP2760Fast);
921     offline = false;
922   }
923   
924   if (offline && CheckCosmics(trigStr)) {
925     fHTrgStatus->Fill(kMUON);
926     offline = false;
927   }
928   if (offline) fHTrgStatus->Fill(kTriggered);
929   Int_t f = 0;
930   if (trgMask & AliVEvent::kMB)          f += fHTrgStatus->Fill(kMinBias);
931   if (trgMask & AliVEvent::kCINT5)       f += fHTrgStatus->Fill(kMinBiasNoSPD);
932   if (trgMask & AliVEvent::kINT7)        f += fHTrgStatus->Fill(kV0AndTrg);
933   if (trgMask & AliVEvent::kHighMult)    f += fHTrgStatus->Fill(kHighMult);
934   if (trgMask & AliVEvent::kCentral)     f += fHTrgStatus->Fill(kCentral);
935   if (trgMask & AliVEvent::kSemiCentral) f += fHTrgStatus->Fill(kSemiCentral);
936   if (trgMask & AliVEvent::kDG5)         f += fHTrgStatus->Fill(kDiffractive);
937   if (trgMask & AliVEvent::kUserDefined) f += fHTrgStatus->Fill(kUser);
938   if (f <= 0) fHTrgStatus->Fill(kOther);
939   
940   // if (!CheckpAExtraV0(esd))             offline = false;
941
942   DMSG(fDebug,2,"Event is %striggered by off-line", offline ? "" : "NOT ");
943
944   if (offline) {
945     triggers |= AliAODForwardMult::kOffline;
946     triggers |= AliAODForwardMult::kInel;
947     CheckINELGT0(esd, nClusters, triggers);
948   }
949   
950   CheckNSD(esd,triggers);
951   CheckPileup(esd, triggers);
952   CheckEmpty(trigStr, triggers);
953   // if (CheckPileup(esd, triggers)) fHTriggers->Fill(kPileUp+.5);
954   // if (CheckEmpty(trigStr, triggers)) fHTriggers->Fill(kEmpty+.5);
955
956   CheckWords(esd, triggers); 
957
958 #define TEST_TRIG_BIN(RET,BIN,TRIGGERS) \
959   do { switch (BIN) { \
960     case kInel:     RET = triggers & AliAODForwardMult::kInel;      break; \
961     case kInelGt0:  RET = triggers & AliAODForwardMult::kInelGt0;   break; \
962     case kNSD:      RET = triggers & AliAODForwardMult::kNSD;       break; \
963     case kV0AND:    RET = triggers & AliAODForwardMult::kV0AND;     break; \
964     case kEmpty:    RET = triggers & AliAODForwardMult::kEmpty;     break; \
965     case kA:        RET = triggers & AliAODForwardMult::kA;         break; \
966     case kB:        RET = triggers & AliAODForwardMult::kB;         break; \
967     case kC:        RET = triggers & AliAODForwardMult::kC;         break; \
968     case kE:        RET = triggers & AliAODForwardMult::kE;         break; \
969     case kPileUp:   RET = triggers & AliAODForwardMult::kPileUp;    break; \
970     case kMCNSD:    RET = triggers & AliAODForwardMult::kMCNSD;     break; \
971     case kSatellite:RET = triggers & AliAODForwardMult::kSatellite; break; \
972     case kOffline:  RET = triggers & AliAODForwardMult::kOffline;   break; \
973     default:        RET = false; } } while(false)
974       
975   if (!fHTriggers) { 
976     AliWarning("Histogram of triggers not defined - has init been called");
977     return false;
978   }
979   
980   for (Int_t i = 0; i < kOffline+1; i++) { 
981     Bool_t hasX = false;
982     TEST_TRIG_BIN(hasX, i, triggers);
983     if (!hasX) continue;
984     fHTriggers->Fill(i+.5);
985     for (Int_t j = 0; j < kOffline+1; j++) { 
986       Bool_t hasY = false;
987       TEST_TRIG_BIN(hasY, j, triggers);
988       if (!hasY) continue;
989       
990       fHTriggerCorr->Fill(i+.5, j+.5);
991     }
992   }
993   return kTRUE;
994 }
995
996 //____________________________________________________________________
997 Bool_t
998 AliFMDEventInspector::CheckFastPartition(bool fastonly) const
999 {
1000   // For the 2.76 TeV p+p run, the FMD ran in the slow partition 
1001   // so it received no triggers from the fast partition. Therefore
1002   // the fast triggers are removed here but not for MC where all 
1003   // triggers are fast.
1004   if (TMath::Abs(fEnergy - 2750.) > 20) return false;
1005   if (fCollisionSystem != AliForwardUtil::kPP) return false;
1006   if (fastonly)
1007     DMSG(fDebug,1,"Fast trigger in pp @ sqrt(s)=2.76TeV removed");
1008
1009   return fastonly;
1010 }
1011
1012 //____________________________________________________________________
1013 Bool_t
1014 AliFMDEventInspector::CheckCosmics(const TString& trigStr) const
1015 {
1016   // MUON triggers are not strictly minimum bias (MB) so they are
1017   // removed (HHD)
1018   if(trigStr.Contains("CMUS1")) {
1019     DMSG(fDebug,1,"Cosmic trigger ins't min-bias, removed");
1020     return true;
1021   }
1022   return false;
1023 }
1024
1025 //____________________________________________________________________
1026 Bool_t
1027 AliFMDEventInspector::CheckINELGT0(const AliESDEvent& esd, 
1028                                    UShort_t& nClusters, 
1029                                    UInt_t& triggers) const
1030 {
1031   nClusters = 0;
1032
1033   // If this is inel, see if we have a tracklet 
1034   const AliMultiplicity* spdmult = esd.GetMultiplicity();
1035   if (!spdmult) {
1036     AliWarning("No SPD multiplicity");
1037     return false;
1038   }
1039
1040   // Check if we have one or more tracklets 
1041   // in the range -1 < eta < 1 to set the INEL>0 
1042   // trigger flag. 
1043   // 
1044   // Also count tracklets as a single cluster 
1045   Int_t n = spdmult->GetNumberOfTracklets();
1046   for (Int_t j = 0; j < n; j++) { 
1047     if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
1048       triggers |= AliAODForwardMult::kInelGt0;
1049       nClusters++;
1050     }
1051   }
1052   n = spdmult->GetNumberOfSingleClusters();
1053   for (Int_t j = 0; j < n; j++) { 
1054     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
1055     if (TMath::Abs(eta) < 1) nClusters++;
1056   }
1057   if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
1058
1059   return triggers & AliAODForwardMult::kNClusterGt0;
1060 }
1061
1062 //____________________________________________________________________
1063 Bool_t
1064 AliFMDEventInspector::CheckNSD(const AliESDEvent& esd, UInt_t& triggers) const
1065 {
1066   // Analyse some trigger stuff 
1067   AliTriggerAnalysis ta;
1068   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kV0AND)) {
1069     triggers |= AliAODForwardMult::kV0AND;
1070     if (fUseV0AND) 
1071       triggers |= AliAODForwardMult::kNSD;
1072   }
1073   if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kNSD1)) 
1074     triggers |= AliAODForwardMult::kNSD;
1075   return triggers & AliAODForwardMult::kNSD;
1076 }
1077 //____________________________________________________________________
1078 Bool_t
1079 AliFMDEventInspector::CheckPileup(const AliESDEvent& esd, 
1080                                   UInt_t& triggers) const
1081 {
1082   // Check for multiple vertices (pile-up) with at least 3
1083   // contributors and at least 0.8cm from the primary vertex
1084   if(fCollisionSystem != AliForwardUtil::kPP) return false;
1085
1086   Bool_t pileup =  esd.IsPileupFromSPD(fMinPileupContrib,fMinPileupDistance);
1087   if (pileup) triggers |= AliAODForwardMult::kPileUp;
1088   return pileup;
1089 }
1090
1091 //____________________________________________________________________
1092 Bool_t
1093 AliFMDEventInspector::CheckEmpty(const TString& trigStr, UInt_t& triggers) const
1094 {
1095   if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
1096     triggers |= AliAODForwardMult::kEmpty;
1097     return true;
1098   }
1099   return false;
1100 }
1101 //____________________________________________________________________
1102 Bool_t
1103 AliFMDEventInspector::CheckWords(const AliESDEvent& esd, UInt_t& triggers) const
1104 {
1105   TObject* word = 0;
1106   TIter nextColl(&fCollWords);
1107   while ((word = nextColl())) {
1108     DMSG(fDebug,10,"Checking if %s trigger %s is fired", 
1109          word->GetTitle(), word->GetName());
1110     if (!esd.IsTriggerClassFired(word->GetName())) continue;
1111
1112     TString beamSide = word->GetTitle();
1113     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1114
1115     if (!beamSide.EqualTo("B")) continue;
1116     triggers |= AliAODForwardMult::kB;
1117     break; // No more to do here
1118   }
1119   TIter nextBg(&fBgWords);
1120   UInt_t all = (AliAODForwardMult::kA | 
1121                 AliAODForwardMult::kC | 
1122                 AliAODForwardMult::kE);
1123   while ((word = nextBg())) {
1124     DMSG(fDebug,10,"Checking if %s trigger %s is fired", 
1125          word->GetTitle(), word->GetName());
1126     if (!esd.IsTriggerClassFired(word->GetName())) continue;
1127     
1128     TString beamSide = word->GetTitle();
1129     DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1130
1131     if (beamSide.Contains("A")) triggers |= AliAODForwardMult::kA;
1132     if (beamSide.Contains("C")) triggers |= AliAODForwardMult::kC;
1133     if (beamSide.Contains("E")) triggers |= AliAODForwardMult::kE;
1134
1135     if ((triggers & all) == all) break; // No more to do
1136   }
1137   return true;
1138 }
1139
1140
1141 //____________________________________________________________________
1142 Bool_t
1143 AliFMDEventInspector::ReadVertex(const AliESDEvent& esd, TVector3& ip)
1144 {
1145   // 
1146   // Read the vertex information from the ESD event 
1147   // 
1148   // Parameters:
1149   //   esd  ESD event 
1150   //   vz   On return, the vertex Z position 
1151   // 
1152   // Return:
1153   //    @c true on success, @c false otherwise 
1154   //
1155   DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
1156   ip.SetXYZ(1024, 1024, 0);
1157   
1158   EVtxStatus s = kNoVtx;
1159   if (fUseFirstPhysicsVertex) 
1160     s = CheckPWGUDVertex(esd, ip);
1161   else if (fUsepA2012Vertex) 
1162     s = CheckpA2012Vertex(esd,ip);      
1163   else 
1164     s = CheckVertex(esd, ip);
1165   
1166   if (s != kVtxOK && fUseDisplacedVertices && fDisplacedVertex.IsSatellite()) {
1167     s = kVtxOK;
1168     ip.SetZ(fDisplacedVertex.GetVertexZ());
1169   }
1170
1171   fHVtxStatus->Fill(s);
1172
1173   return s == kVtxOK;
1174 }
1175
1176 //____________________________________________________________________
1177 AliFMDEventInspector::EVtxStatus
1178 AliFMDEventInspector::CheckPWGUDVertex(const AliESDEvent& esd, 
1179                                        TVector3& ip)  const
1180 {
1181   // This is the code used by the 1st physics people 
1182   const AliESDVertex* vertex    = esd.GetPrimaryVertex();
1183   if (!vertex  || !vertex->GetStatus()) {
1184     DMSG(fDebug,2,"No primary vertex (%p) or bad status %d", 
1185          vertex, (vertex ? vertex->GetStatus() : -1));
1186     return kNoVtx;
1187   }
1188   const AliESDVertex* vertexSPD = esd.GetPrimaryVertexSPD();
1189   if (!vertexSPD || !vertexSPD->GetStatus()) {
1190     DMSG(fDebug,2,"No primary SPD vertex (%p) or bad status %d", 
1191          vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1));
1192     return kNoSPDVtx;
1193   }
1194     
1195   // if vertex is from SPD vertexZ, require more stringent cuts 
1196   if (vertex->IsFromVertexerZ()) { 
1197     if (vertex->GetDispersion() > fMaxVzErr || 
1198         vertex->GetZRes() > 1.25 * fMaxVzErr) {
1199       DMSG(fDebug,2,"Dispersion %f > %f or resolution %f > %f",
1200            vertex->GetDispersion(), fMaxVzErr,
1201            vertex->GetZRes(), 1.25 * fMaxVzErr);
1202       return kUncertain;
1203     }
1204   }
1205   ip.SetZ(vertex->GetZ());
1206   
1207   if(!vertex->IsFromVertexerZ()) {
1208     ip.SetX(vertex->GetX());
1209     ip.SetY(vertex->GetY());
1210   }
1211   return kVtxOK;
1212 }
1213 //____________________________________________________________________
1214 AliFMDEventInspector::EVtxStatus
1215 AliFMDEventInspector::CheckpA2012Vertex(const AliESDEvent& esd, 
1216                                         TVector3& ip)  const
1217 {      
1218   const AliESDVertex *vertex = esd.GetPrimaryVertexSPD();
1219   if (!vertex) return kNoSPDVtx;
1220   if (vertex->GetNContributors() <= 0) return kFewContrib;
1221   
1222   TString vtxTyp = vertex->GetTitle();
1223   if (vtxTyp.Contains("vertexer: Z")) return kNotVtxZ;
1224
1225   if (vertex->GetDispersion() >= 0.04 || vertex->GetZRes()>=0.25) 
1226     return kUncertain;
1227
1228   ip.SetX(vertex->GetX());
1229   ip.SetY(vertex->GetY());
1230   ip.SetZ(vertex->GetZ());              
1231
1232   return kVtxOK;
1233 }
1234
1235 //____________________________________________________________________
1236 AliFMDEventInspector::EVtxStatus
1237 AliFMDEventInspector::CheckVertex(const AliESDEvent& esd, 
1238                                   TVector3& ip) const
1239 {
1240   // Use standard SPD vertex (perhaps preferable for Pb+Pb)
1241   // Get the vertex 
1242   const AliESDVertex* vertex = esd.GetPrimaryVertexSPD();
1243   if (!vertex) { 
1244     if (fDebug > 2) {
1245       AliWarning("No SPD vertex found in ESD"); }
1246     return kNoSPDVtx;
1247   }
1248     
1249   // #if 0 // Check disabled - seem to kill a lot of PbPb events
1250   // Check that enough tracklets contributed 
1251   if(vertex->GetNContributors() <= 0) {
1252     DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
1253          vertex->GetNContributors());
1254     ip.SetZ(0);
1255     return kFewContrib;
1256   } 
1257   // #endif
1258
1259   // Check that the uncertainty isn't too large 
1260   if (vertex->GetZRes() > fMaxVzErr) { 
1261     DMSG(fDebug,2,"Uncertaintity in Z of vertex is too large %f > %f", 
1262          vertex->GetZRes(), fMaxVzErr);
1263     return kUncertain;
1264   }
1265     
1266   // Get the z coordiante 
1267   ip.SetZ(vertex->GetZ());
1268   const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
1269     
1270   
1271   if(!vertexXY->IsFromVertexerZ()) {
1272     ip.SetX(vertexXY->GetX());
1273     ip.SetY(vertexXY->GetY());
1274   }
1275   return kVtxOK;
1276 }
1277
1278 //____________________________________________________________________
1279 Bool_t
1280 AliFMDEventInspector::ReadRunDetails(const AliESDEvent* esd)
1281 {
1282   // 
1283   // Read the collision system, collision energy, and L3 field setting
1284   // from the ESD
1285   // 
1286   // Parameters:
1287   //   esd ESD to get information from 
1288   // 
1289   // Return:
1290   //    true on success, false 
1291   //
1292   // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
1293   //           esd->GetBeamType(), 2*esd->GetBeamEnergy(), 
1294   //           esd->GetMagneticField()));
1295   DGUARD(fDebug,2,"Read the run details in AliFMDEventInspector");
1296   const char* sys  = esd->GetBeamType();
1297   Float_t     cms  = 2 * esd->GetBeamEnergy();
1298   Float_t     fld  = esd->GetMagneticField();
1299   fCollisionSystem = AliForwardUtil::ParseCollisionSystem(sys);
1300   fEnergy          = AliForwardUtil::ParseCenterOfMassEnergy(fCollisionSystem, 
1301                                                              cms);
1302   fField           = AliForwardUtil::ParseMagneticField(fld);
1303   fRunNumber       = esd->GetRunNumber();
1304   StoreInformation();
1305   if (fCollisionSystem   == AliForwardUtil::kUnknown) { 
1306     AliWarningF("Unknown collision system: %s - please check", sys);
1307     return false;
1308   }
1309   if (fEnergy            <= 0) {
1310     AliWarningF("Unknown CMS energy: %f (%d) - please check", cms, fEnergy);
1311     return false;
1312   }
1313   if (TMath::Abs(fField) >  10) {
1314     AliWarningF("Unknown L3 field setting: %f (%d) - please check", fld,fField);
1315     return false;
1316   }
1317
1318   return true;
1319 }
1320
1321
1322 //____________________________________________________________________
1323 const Char_t*
1324 AliFMDEventInspector::CodeString(UInt_t code)
1325 {
1326   static TString s;
1327   s = "";
1328   if (code & kNoEvent)    s.Append("NOEVENT ");
1329   if (code & kNoTriggers) s.Append("NOTRIGGERS ");
1330   if (code & kNoSPD)      s.Append("NOSPD ");
1331   if (code & kNoFMD)      s.Append("NOFMD ");
1332   if (code & kNoVertex)   s.Append("NOVERTEX ");
1333   if (code & kBadVertex)  s.Append("BADVERTEX ");
1334   return s.Data();
1335 }
1336 //____________________________________________________________________
1337 void
1338 AliFMDEventInspector::Print(Option_t*) const
1339 {
1340   // 
1341   // Print information
1342   // 
1343   //   option Not used 
1344   //
1345   char ind[gROOT->GetDirLevel()+1];
1346   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1347   ind[gROOT->GetDirLevel()] = '\0';
1348   TString sNN(AliForwardUtil::CenterOfMassEnergyString(fEnergy));
1349   sNN.Strip(TString::kBoth, '0');
1350   sNN.ReplaceAll("GeV", " GeV");
1351   TString field(AliForwardUtil::MagneticFieldString(fField));
1352   field.ReplaceAll("p",  "+");
1353   field.ReplaceAll("m",  "-");
1354   field.ReplaceAll("kG", " kG");
1355   
1356   std::cout << std::boolalpha 
1357             << ind << ClassName() << ": " << GetName() << '\n'
1358             << ind << " Vertex bins:            " << fVtxAxis.GetNbins() << '\n'
1359             << ind << " Vertex range:           [" << fVtxAxis.GetXmin() 
1360             << "," << fVtxAxis.GetXmax() << "]\n"
1361             << ind << " Low flux cut:           " << fLowFluxCut << '\n'
1362             << ind << " Max(delta v_z):         " << fMaxVzErr << " cm\n"
1363             << ind << " Min(nContrib_pileup):   " << fMinPileupContrib << '\n'
1364             << ind << " Min(v-pileup):          " << fMinPileupDistance << '\n'
1365             << ind << " System:                 " 
1366             << AliForwardUtil::CollisionSystemString(fCollisionSystem) << '\n'
1367             << ind << " CMS energy per nucleon: " << sNN << '\n'
1368             << ind << " Field:                  " <<  field << '\n'
1369             << ind << " Satellite events:       " << fUseDisplacedVertices<<'\n'
1370             << ind << " Centrality method:      " << fCentMethod << '\n'
1371             << std::noboolalpha;
1372   if (!fCentAxis) { std::cout << std::flush; return; }
1373   Int_t nBin = fCentAxis->GetNbins();
1374   std::cout << ind << " Centrality axis:        " << nBin << " bins"
1375             << std::flush;
1376   for (Int_t i = 0; i < nBin; i++) { 
1377     if ((i % 10) == 0) std::cout << '\n' << ind << "  ";
1378     std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
1379   }
1380   std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
1381 }
1382
1383   
1384 //
1385 // EOF
1386 //
1387