AliPhysics  db95e02 (db95e02)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFMDEventInspector.cxx
Go to the documentation of this file.
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 "AliProdInfo.h"
19 #include "AliLog.h"
20 #include "AliESDEvent.h"
21 #include "AliMultiplicity.h"
22 #include "AliAnalysisManager.h"
23 #include "AliMCEventHandler.h"
24 #include "AliInputEventHandler.h"
25 #include "AliTriggerAnalysis.h"
26 #include "AliPhysicsSelection.h"
27 #include "AliOADBPhysicsSelection.h"
28 #include "AliAODForwardMult.h"
29 #include "AliForwardUtil.h"
30 #include "AliCentrality.h"
31 #include "AliVVZERO.h"
32 #include <TH1.h>
33 #include <TList.h>
34 #include <TDirectory.h>
35 #include <TROOT.h>
36 #include <TParameter.h>
37 #include <TMatrixDSym.h>
38 #include <TPRegexp.h>
39 #include <iostream>
40 #include <iomanip>
41 #include "AliMCEvent.h"
42 #include "AliHeader.h"
43 #include "AliGenEventHeader.h"
44 #include "AliCollisionGeometry.h"
45 #include "AliVVZERO.h"
46 
47 //====================================================================
48 namespace {
50  // (AliAODForwardMult::kInel |
51  // AliAODForwardMult::kV0AND |
52  // AliAODForwardMult::kADOR |
53  // AliAODForwardMult::kADAND |
54  // AliAODForwardMult::kSatellite |
55  // AliAODForwardMult::kOffline |
56  // AliAODForwardMult::kInclusive);
57 
58 
59  AliPhysicsSelection* GetPhysicsSelection()
60  {
61  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
62  if (!am) {
63  Warning("GetPhysicsSelection", "No analysis manager");
64  return 0;
65  }
66 
67  // Get the input handler - should always be there
68  AliInputEventHandler* ih =
69  dynamic_cast<AliInputEventHandler*>(am->GetInputEventHandler());
70  if (!ih) {
71  Warning("GetPhysicsSelection", "No input handler");
72  return 0;
73  }
74  // Get the physics selection - should always be there
75  AliPhysicsSelection* ps =
76  dynamic_cast<AliPhysicsSelection*>(ih->GetEventSelection());
77  if (!ps) {
78  Warning("GetPhysicsSelection", "No physics selection");
79  return 0;
80  }
81  return ps;
82  }
83 }
84 
85 //====================================================================
86 const char* AliFMDEventInspector::fgkFolderName = "fmdEventInspector";
87 
88 //____________________________________________________________________
90  : TNamed(),
91  fHEventsTr(0),
92  fHEventsTrVtx(0),
93  fHEventsAccepted(0),
94  fHEventsAcceptedXY(0),
95  fHTriggers(0),
96  fHType(0),
97  fHWords(0),
98  fHCent(0),
99  fHCentVsQual(0),
100  fHStatus(0),
101  fHVtxStatus(0),
102  fHTrgStatus(0),
103  fHPileup(0),
104  fLowFluxCut(1000),
105  fMaxVzErr(0.2),
106  fList(0),
107  fEnergy(0),
108  fField(999),
109  fCollisionSystem(kUnknown),
110  fDebug(0),
111  fCentAxis(0),
112  fVtxAxis(10,-10,10),
113  fVtxMethod(kNormal),
114  // fUseFirstPhysicsVertex(false),
115  fUseV0AND(false),
116  fPileupFlags(kSPD|kTracks|kOutOfBunch|kSPDBins),
117  fMinPileupContrib(3),
118  fMinPileupDistance(0.8),
119  // fUseDisplacedVertices(false),
120  fDisplacedVertex(),
121  fCollWords(),
122  fBgWords(),
123  fCentMethod("default"),
124  fMinCent(-1.0),
125  fMaxCent(-1.0),
126  // fUsepA2012Vertex(false),
127  fRunNumber(0),
128  fMC(false),
129  fProdYear(-1),
130  fProdLetter('?'),
131  fProdPass(-1),
132  fProdSVN(-1),
133  fProdMC(false)
134 {
135  //
136  // Constructor
137  //
138  DGUARD(fDebug,1,"Default CTOR of AliFMDEventInspector");
139 }
140 
141 //____________________________________________________________________
143  : TNamed(fgkFolderName, name),
144  fHEventsTr(0),
145  fHEventsTrVtx(0),
146  fHEventsAccepted(0),
147  fHEventsAcceptedXY(0),
148  fHTriggers(0),
149  fHType(0),
150  fHWords(0),
151  fHCent(0),
152  fHCentVsQual(0),
153  fHStatus(0),
154  fHVtxStatus(0),
155  fHTrgStatus(0),
156  fHPileup(0),
157  fLowFluxCut(1000),
158  fMaxVzErr(0.2),
159  fList(0),
160  fEnergy(0),
161  fField(999),
162  fCollisionSystem(kUnknown),
163  fDebug(0),
164  fCentAxis(0),
165  fVtxAxis(10,-10,10),
166  fVtxMethod(kNormal),
167 // fUseFirstPhysicsVertex(false),
168  fUseV0AND(false),
169  fPileupFlags(kSPD|kTracks|kOutOfBunch|kSPDBins),
170  fMinPileupContrib(3),
171  fMinPileupDistance(0.8),
172 // fUseDisplacedVertices(false),
173  fDisplacedVertex(),
174  fCollWords(),
175  fBgWords(),
176  fCentMethod("default"),
177  fMinCent(-1.0),
178  fMaxCent(-1.0),
179 // fUsepA2012Vertex(false),
180  fRunNumber(0),
181  fMC(false),
182  fProdYear(-1),
183  fProdLetter('?'),
184  fProdPass(-1),
185  fProdSVN(-1),
186  fProdMC(false)
187 {
188  //
189  // Constructor
190  //
191  // Parameters:
192  // name Name of object
193  //
194  DGUARD(fDebug,1,"Named CTOR of AliFMDEventInspector: %s", name);
195  if (!GetPhysicsSelection()) {
196  AliWarning("No physics selection defined in manager\n"
197  "=======================================================\n"
198  " The use of AliFMDEventInspector _suggests_ a valid\n"
199  " AliPhysicsSelection registered with the input handler.\n\n"
200  " Please check our train setup\n"
201  "=======================================================\n");
202  }
203 }
204 
205 //____________________________________________________________________
207 {
208  //
209  // Destructor
210  //
211  DGUARD(fDebug,1,"DTOR of AliFMDEventInspector");
212  // if (fList) delete fList;
213 }
214 
215 //____________________________________________________________________
216 void
218 {
219  switch (m) {
220  case kV0Multiplicity: fCentMethod = "VOM"; break; // VZERO multiplicity
221  case kV0Amplitude: fCentMethod = "V0A"; break; // VZERO amplitude
222  case kV0Charge: fCentMethod = "V0C"; break; // VZERO charge
223  case kFMDRough: fCentMethod = "FMD"; break; // FMD scaled energy l
224  case kNTracks: fCentMethod = "TRK"; break; // Number of tracks
225  case kLTracks: fCentMethod = "TKL"; break; // Number of tracks
226  case kCL0: fCentMethod = "CL0"; break; //
227  case kCL1: fCentMethod = "CL1"; break; //
228  case kCND: fCentMethod = "CND"; break; //
229  case kNParticles: fCentMethod = "NPA"; break; // Neutral particles
230  case kNeutrons: fCentMethod = "ZNA"; break; // ZDC neutron amplitu
231  case kV0vsFMD: fCentMethod = "V0MvsFMD"; break; // VZERO versus FMD
232  case kV0vsNTracks: fCentMethod = "TKLvsVOM"; break; // Tracks versus VZERO
233  case kZEMvsZDC: fCentMethod = "ZEMvsZDC"; break; // ZDC
234  default: fCentMethod = "default"; break;
235  }
236 }
237 
238 //____________________________________________________________________
239 void
241 {
242  AliWarning("\n"
243  "*******************************************************\n"
244  "* Setting centrality cuts in this stage is deprecated *\n"
245  "*******************************************************");
246  fMinCent = minCent;
247 }
248 //____________________________________________________________________
249 void
251 {
252  AliWarning("\n"
253  "*******************************************************\n"
254  "* Setting centrality cuts in this stage is deprecated *\n"
255  "*******************************************************");
256  fMaxCent = maxCent;
257 }
258 
259 //____________________________________________________________________
260 Bool_t
262  TH1I*& hEventsTr,
263  TH1I*& hEventsTrVtx,
264  TH1I*& hEventsAcc,
265  TH1I*& hTriggers) const
266 {
267  //
268  // Fetch our histograms from the passed list
269  //
270  // Parameters:
271  // d Input
272  // hEventsTr On return, pointer to histogram, or null
273  // hEventsTrVtx On return, pointer to histogram, or null
274  // hTriggers On return, pointer to histogram, or null
275  //
276  // Return:
277  // true on success, false otherwise
278  //
279  DGUARD(fDebug,3,"Fetch histograms in AliFMDEventInspector");
280  hEventsTr = 0;
281  hEventsTrVtx = 0;
282  hEventsAcc = 0;
283  hTriggers = 0;
284  TList* dd = dynamic_cast<TList*>(d->FindObject(GetName()));
285  if (!dd) return kFALSE;
286 
287  hEventsTr = dynamic_cast<TH1I*>(dd->FindObject("nEventsTr"));
288  hEventsTrVtx = dynamic_cast<TH1I*>(dd->FindObject("nEventsTrVtx"));
289  hEventsAcc = dynamic_cast<TH1I*>(dd->FindObject("nEventsAccepted"));
290  hTriggers = dynamic_cast<TH1I*>(dd->FindObject("triggers"));
291 
292  if (!hEventsTr ||
293  !hEventsTrVtx ||
294  !hEventsAcc ||
295  !hTriggers) return kFALSE;
296  return kTRUE;
297 }
298 //____________________________________________________________________
299 void
301  const TList* classes,
303 {
304  TIter nextClass(classes);
305  TObjString* trigClass = 0;
306  // Loop over all trigger classes. Trigger classes have the format
307  //
308  // class := positive_words SPACE(s) negative_words
309  // positive_words :=
310  // | '+' words
311  // negative_words :=
312  // | '-' words
313  // words := word
314  // | word ',' words
315  //
316  Printf("Inspecting physics selection: %p", o);
317  // o->Print();
318  while ((trigClass = static_cast<TObjString*>(nextClass()))) {
319  // Tokenize on space to get positive and negative parts
320  TString side = o->GetBeamSide(trigClass->String());
321  TObjArray* parts = trigClass->String().Tokenize(" ");
322  TObjString* part = 0;
323  TIter nextPart(parts);
324  while ((part = static_cast<TObjString*>(nextPart()))) {
325  // We only care about the positive ones
326  // if (part->GetName()[0] == '-') continue;
327  // part->String().Remove(0,1);
328 
329  // Tokenize on a comma to get the words
330  TObjArray* words = part->String().Tokenize(",");
331  TObjString* word = 0;
332  TIter nextWord(words);
333  while ((word = static_cast<TObjString*>(nextWord()))) {
334  char first = word->String()[0];
335  switch (first) {
336  case '-': // Negatives
337  case '&': // bunch crossing
338  case '*': // ?
339  continue;
340  break;
341  case '+':
342  word->String().Remove(0,1);
343  break;
344  }
345  if (cache.FindObject(word->String())) {
346  DMSG(fDebug, 3, "Word %s already cached", word->GetName());
347  continue;
348  }
349  TNamed* store = new TNamed(word->String(), side);
350  cache.Add(store);
351  DMSG(fDebug,3,"Caching %s trigger word %s",
352  store->GetTitle(), store->GetName());
353  } // while (word)
354  delete words;
355  }
356  delete parts;
357  }
358 }
359 
360 //____________________________________________________________________
361 void
363 {
364  //
365  // Initialize the object - this is called on the first seen event.
366  //
367  // Parameters:
368  // vtxAxis Vertex axis in use
369  //
370  DGUARD(fDebug,1,"Initialize in AliFMDEventInspector");
371 
372  // Get the physics selection - should always be there
373  AliPhysicsSelection* ps = GetPhysicsSelection();
374  if (!ps) {
375  AliWarning("No physics selection");
376  }
377  else {
378  // Get the configured triggers
379  AliOADBPhysicsSelection* oadb =
380  const_cast<AliOADBPhysicsSelection*>(ps->GetOADBPhysicsSelection());
381  if (!oadb) {
382  AliWarning("No OADB physics selection object");
383  return;
384  }
385  // Get the configured trigger words from the physics selection
386  const TList* collTriggClasses = ps->GetCollisionTriggerClasses();
387  const TList* bgTriggClasses = ps->GetBGTriggerClasses();
388  if (!collTriggClasses) {
389  AliWarning("No configured collision trigger classes");
390  return;
391  }
392  if (!bgTriggClasses) {
393  AliWarning("No configured background trigger classes");
394  return;
395  }
396 #if 1
397  // This histogram disabled as it causes problems in the merge
398  fHWords = new TH1I("words", "Trigger words seen", 1, 0, 0);
399  fHWords->SetFillColor(kBlue-2);
400  fHWords->SetFillStyle(3001);
401  fHWords->SetStats(0);
402  fHWords->SetDirectory(0);
403 #endif
404 
405  CacheConfiguredTriggerClasses(fCollWords, collTriggClasses, oadb);
406  CacheConfiguredTriggerClasses(fBgWords, bgTriggClasses, oadb);
407  fCollWords.Sort();
408  fBgWords.Sort();
409  DMSG(fDebug, 3, "fHWords=%p", fHWords);
410  if (fHWords) {
411  Int_t nColl = fCollWords.GetEntries();
412  Int_t nBg = fBgWords.GetEntries();
413  TAxis* xAxis = fHWords->GetXaxis();
414 
415  DMSG(fDebug,3, "Got a total of %d+%d=%d trigger words",
416  nColl, nBg, nColl+nBg);
417  xAxis->Set(nColl+nBg, 0.5, nColl+nBg+.5);
418  TIter nextC(&fCollWords);
419  TIter nextB(&fBgWords);
420  Int_t bin = 1;
421  TObject* word = 0;
422  while((word = nextC())) {
423  DMSG(fDebug, 3, "Defining bin %d to be %s", bin, word->GetName());
424  word->SetUniqueID(bin);
425  xAxis->SetBinLabel(bin++, word->GetName());
426  }
427  while((word = nextB())) {
428  DMSG(fDebug, 3, "Defining bin %d to be %s", bin, word->GetName());
429  word->SetUniqueID(bin);
430  xAxis->SetBinLabel(bin++, word->GetName());
431  }
432  fHWords->Rebuild();
433  }
434  if (fDebug >= 3) {
435  DMSG(fDebug, 3, "Collision trigger classes");
436  fCollWords.ls();
437  DMSG(fDebug, 3, "Background trigger classes");
438  fBgWords.ls();
439  }
440  }
441 
442  TArrayD limits;
443  if ((fMinCent < 0 && fMaxCent < 0) || fMaxCent <= fMinCent) {
444  // Was:
445  // -1.5 -0.5 0.5 1.5 ... 89.5 ... 100.5
446  // ----- 92 number --------- ---- 1 ---
447  // Now:
448  // -0.5 0.5 1.5 ... 89.5 ... 100.5
449  // ----- 91 number ---- ---- 1 ---
450  limits.Set(92);
451  for (Int_t i = 0; i < 91; i++) limits[i] = -.5 + i;
452  limits[91] = 100.5;
453  }
454  else {
455  Int_t n = Int_t(fMaxCent-fMinCent+2);
456  limits.Set(n);
457  for (Int_t i = 0; i < n; i++) {
458  limits[i] = fMinCent + i - .5;
459  }
460  }
461 
462  if (vtxAxis.GetXbins() && vtxAxis.GetXbins()->GetArray())
463  fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXbins()->GetArray());
464  else
465  fVtxAxis.Set(vtxAxis.GetNbins(), vtxAxis.GetXmin(), vtxAxis.GetXmax());
466 
467  fCentAxis = new TAxis(limits.GetSize()-1, limits.GetArray());
468  fHEventsTr = 0;
469  if (vtxAxis.GetXbins() && vtxAxis.GetXbins()->GetArray())
470  fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
471  vtxAxis.GetNbins(), vtxAxis.GetXbins()->GetArray());
472  else
473  fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
474  4*vtxAxis.GetNbins(),
475  2*vtxAxis.GetXmin(),
476  2*vtxAxis.GetXmax());
477  fHEventsTr->SetXTitle("v_{z} [cm]");
478  fHEventsTr->SetYTitle("# of events");
479  fHEventsTr->SetFillColor(kRed+1);
480  fHEventsTr->SetFillStyle(3001);
481  fHEventsTr->SetDirectory(0);
482  // fHEventsTr->Sumw2();
483  fList->Add(fHEventsTr);
484 
485  fHEventsTrVtx = static_cast<TH1I*>(fHEventsTr->Clone("nEventsTrVtx"));
486  fHEventsTrVtx->SetTitle("Number of events w/trigger and vertex");
487  fHEventsTrVtx->SetFillColor(kBlue+1);
488  fHEventsTrVtx->SetDirectory(0);
489  // fHEventsTrVtx->Sumw2();
490  fList->Add(fHEventsTrVtx);
491 
492  fHEventsAccepted = 0;
493  if (vtxAxis.GetXbins() && vtxAxis.GetXbins()->GetArray())
495  new TH1I("nEventsAccepted",
496  "Number of events w/trigger and vertex in range",
497  vtxAxis.GetNbins(), vtxAxis.GetXbins()->GetArray());
498  else
500  new TH1I("nEventsAccepted",
501  "Number of events w/trigger and vertex in range",
502  2*vtxAxis.GetNbins(),
503  2*vtxAxis.GetXmin(),
504  2*vtxAxis.GetXmax());
505  fHEventsAccepted->SetXTitle("v_{z} [cm]");
506  fHEventsAccepted->SetYTitle("# of events");
507  fHEventsAccepted->SetFillColor(kGreen+1);
508  fHEventsAccepted->SetFillStyle(3001);
509  fHEventsAccepted->SetDirectory(0);
510  // fHEventsAccepted->Sumw2();
511  fList->Add(fHEventsAccepted);
512 
513  fHEventsAcceptedXY = new TH2D("nEventsAcceptedXY",
514  "XY vertex w/trigger and Z vertex in range",
515  1000,-1,1,1000,-1,1);
516 
517  fHEventsAcceptedXY->SetXTitle("v_{x} [cm]");
518  fHEventsAcceptedXY->SetYTitle("v_{y} [cm]");
519  fHEventsAcceptedXY->SetDirectory(0);
520  // fHEventsAccepted->Sumw2();
522 
523  fHTriggers = AliAODForwardMult::MakeTriggerHistogram("triggers",kTriggerMask);
524  fHTriggers->SetStats(0);
525  fHTriggers->SetDirectory(0);
526  fList->Add(fHTriggers);
527 
528 
529  fHType = new TH1I("type", Form("Event type (cut: SPD mult>%d)",
530  fLowFluxCut), 2, -.5, 1.5);
531  fHType->SetFillColor(kRed+1);
532  fHType->SetFillStyle(3001);
533  fHType->SetStats(0);
534  fHType->SetDirectory(0);
535  fHType->GetXaxis()->SetBinLabel(1,"Low-flux");
536  fHType->GetXaxis()->SetBinLabel(2,"High-flux");
537  fList->Add(fHType);
538 
539  if (fHWords) fList->Add(fHWords);
540 
541  fHCent = new TH1F("cent", "Centrality", limits.GetSize()-1,limits.GetArray());
542  fHCent->SetFillColor(kBlue+1);
543  fHCent->SetFillStyle(3001);
544  fHCent->SetStats(0);
545  fHCent->SetDirectory(0);
546  fHCent->SetXTitle("Centrality [%]");
547  fHCent->SetYTitle("Events");
548  fList->Add(fHCent);
549 
550  fHCentVsQual = new TH2F("centVsQuality", "Quality vs Centrality",
551  9, 196.5, 205.5,
552  limits.GetSize()-1, limits.GetArray());
553  fHCentVsQual->SetXTitle("Quality");
554  fHCentVsQual->SetYTitle("Centrality [%]");
555  fHCentVsQual->SetZTitle("Events");
556  TAxis* qAxis = fHCentVsQual->GetXaxis();
557  qAxis->SetBinLabel(1, "OK");
558  qAxis->SetBinLabel(2, "198: beyond anchor");
559  qAxis->SetBinLabel(3, "199: no calib");
560  qAxis->SetBinLabel(4, "200: not desired trigger");
561  qAxis->SetBinLabel(5, "201: not INEL0 with tracklets");
562  qAxis->SetBinLabel(6, "202: vertex Z not within 10cm");
563  qAxis->SetBinLabel(7, "203: tagged as pileup (SPD)");
564  qAxis->SetBinLabel(8, "204: inconsistent SPD/tracking vertex");
565  qAxis->SetBinLabel(9, "205: rejected by tracklets-vs-clusters");
566 
567  fHCentVsQual->SetDirectory(0);
568  fList->Add(fHCentVsQual);
569 
570  fHStatus = new TH1I("status", "Status", 8, 1, 9);
571  fHStatus->SetFillColor(kBlue+1);
572  fHStatus->SetFillStyle(3001);
573  fHStatus->SetStats(0);
574  fHStatus->SetDirectory(0);
575  TAxis* xAxis = fHStatus->GetXaxis();
576  xAxis->SetBinLabel(1, "OK");
577  xAxis->SetBinLabel(2, "No event");
578  xAxis->SetBinLabel(3, "No triggers");
579  xAxis->SetBinLabel(4, "No SPD");
580  xAxis->SetBinLabel(5, "No FMD");
581  xAxis->SetBinLabel(6, "No vertex");
582  xAxis->SetBinLabel(7, "Bad vertex");
583  xAxis->SetBinLabel(8, "No centrality");
584  fList->Add(fHStatus);
585 
586  fHVtxStatus = new TH1I("vtxStatus","Vertex Status",
588  fHVtxStatus->SetFillColor(kGreen+1);
589  fHVtxStatus->SetFillStyle(3001);
590  fHVtxStatus->SetStats(0);
591  fHVtxStatus->SetDirectory(0);
592  xAxis = fHVtxStatus->GetXaxis();
593  xAxis->SetBinLabel(kVtxOK, "OK");
594  xAxis->SetBinLabel(kNoVtx, "None/bad status");
595  xAxis->SetBinLabel(kNoSPDVtx, "No SPD/bad status");
596  xAxis->SetBinLabel(kFewContrib, "N_{contrib} <= 0");
597  xAxis->SetBinLabel(kUncertain, Form("#delta z > %4.2f", fMaxVzErr));
598  xAxis->SetBinLabel(kNotVtxZ, "Not Z vertexer");
599  fList->Add(fHVtxStatus);
600 
601  fHTrgStatus = new TH1I("trgStatus", "Trigger Status",
603  fHTrgStatus->SetFillColor(kMagenta+1);
604  fHTrgStatus->SetFillStyle(3001);
605  fHTrgStatus->SetStats(0);
606  fHTrgStatus->SetDirectory(0);
607  xAxis = fHTrgStatus->GetXaxis();
608  xAxis->SetBinLabel(kNoTrgWords, "No words");
609  xAxis->SetBinLabel(kPP2760Fast, "FAST in pp@#sqrt{s}=2.76TeV");
610  xAxis->SetBinLabel(kMUON, "Muon trigger");
611  xAxis->SetBinLabel(kTriggered, "Triggered");
612  xAxis->SetBinLabel(kMinBias, "CINT1 (V0A||V0C||FASTOR)");
613  xAxis->SetBinLabel(kMinBiasNoSPD, "CINT5 (V0A||V0C)");
614  xAxis->SetBinLabel(kV0AndTrg, "CINT7 (V0A&&V0C)");
615  xAxis->SetBinLabel(kMinBiasAD, "CINT10 (ADA||ADC||V0A||V0C||FASTOR)");
616  xAxis->SetBinLabel(kHighMult, "N>>0");
617  xAxis->SetBinLabel(kCentral, "Central");
618  xAxis->SetBinLabel(kSemiCentral, "Semi-central");
619  xAxis->SetBinLabel(kDiffractive, "Diffractive");
620  xAxis->SetBinLabel(kADOr, "(ADA || ADC)");
621  xAxis->SetBinLabel(kADAnd, "(ADA && ADC)");
622  xAxis->SetBinLabel(kUser, "User");
623  xAxis->SetBinLabel(kOther, "Other");
624  fList->Add(fHTrgStatus);
625 
626  fHPileup = new TH1I("pileupStatus", "Pile-up status", 4, 0, 4);
627  fHPileup->SetFillColor(kYellow+2);
628  fHPileup->SetFillStyle(3001);
629  fHPileup->SetStats(0);
630  fHPileup->SetDirectory(0);
631  xAxis = fHPileup->GetXaxis();
632  xAxis->SetBinLabel(1, "SPD tracklets");
633  xAxis->SetBinLabel(2, "Tracks");
634  xAxis->SetBinLabel(3, "Out-of-bunch");
635  xAxis->SetBinLabel(4, "SPD tracklets (bins)");
636  fList->Add(fHPileup);
637 
639 }
640 
641 //____________________________________________________________________
642 void
644 {
645  // Write TNamed objects to output list containing information about
646  // the running conditions
647  DGUARD(fDebug,2,"Store information from AliFMDEventInspector");
648  if (!fList) return;
649 
650 
657  //fList->Add(AliForwardUtil::MakeParameter("fpVtx",fUseFirstPhysicsVertex));
662  fList->Add(AliForwardUtil::MakeParameter("satellite", AllowDisplaced()));
663  fList->Add(AliForwardUtil::MakeParameter("alirootRev",
665  fList->Add(AliForwardUtil::MakeParameter("alirootBranch",
668 
669 }
670 
671 //____________________________________________________________________
672 void
674 {
675  if (!fList) return;
676 
677  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
678  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
679  if (!inputHandler) {
680  AliWarning("Got no input handler");
681  return;
682  }
683  TList *uiList = inputHandler->GetUserInfo();
684  if (!uiList) {
685  AliWarning("Got no user list from input tree");
686  return;
687  }
688 
689  AliProdInfo p(uiList);
690  p.List();
691  if (p.GetAlirootSvnVersion() <= 0) return;
692 
693  // Make our output list
694  TList* out = new TList;
695  out->SetOwner(true);
696  out->SetName("production");
697  // out->SetTitle("ESD production details");
698  fList->Add(out);
699 
700  TString period = p.GetLHCPeriod();
701  // TString aliROOTVersion = p.GetAlirootVersion();
702  fProdSVN = p.GetAlirootSvnVersion();
703  // TString rootVersion = p.GetRootVersion();
704  Int_t rootSVN = p.GetRootSvnVersion();
705  fProdPass = p.GetRecoPass();
706  fProdMC = p.IsMC();
707 
708  if (period.Length() > 0) {
709  TObjArray* pp = TPRegexp("LHC([0-9]+)([a-z]+)").MatchS(period);
710  fProdYear = static_cast<TObjString*>(pp->At(1))->String().Atoi();
711  fProdLetter = static_cast<TObjString*>(pp->At(2))->String()[0];
712  pp->Delete();
713  }
714 
715  out->Add(AliForwardUtil::MakeParameter("year", fProdYear));
716  out->Add(AliForwardUtil::MakeParameter("letter", Int_t(fProdLetter)));
717  out->Add(AliForwardUtil::MakeParameter("alirootSVN", fProdSVN));
718  out->Add(AliForwardUtil::MakeParameter("rootSVN", rootSVN));
719  out->Add(AliForwardUtil::MakeParameter("pass", fProdPass));
720  out->Add(AliForwardUtil::MakeParameter("mc", fProdMC));
721 }
722 
723 
724 
725 
726 //____________________________________________________________________
727 void
729 {
730  //
731  // Define the output histograms. These are put in a sub list of the
732  // passed list. The histograms are merged before the parent task calls
733  // AliAnalysisTaskSE::Terminate
734  //
735  // dir Directory to add to
736  //
737  DGUARD(fDebug,1,"Define output from AliFMDEventInspector");
738  fList = new TList;
739  fList->SetName(GetName());
740  fList->SetOwner();
741  dir->Add(fList);
742 }
743 
744 //____________________________________________________________________
745 UInt_t
747  UInt_t& triggers,
748  Bool_t& lowFlux,
749  UShort_t& ivz,
750  TVector3& ip,
751  Double_t& cent,
752  UShort_t& nClusters)
753 {
754  //
755  // Process the event
756  //
757  // Parameters:
758  // event Input event
759  // triggers On return, the triggers fired
760  // lowFlux On return, true if the event is considered a low-flux
761  // event (according to the setting of fLowFluxCut)
762  // ivz On return, the found vertex bin (1-based). A zero
763  // means outside of the defined vertex range
764  // vz On return, the z position of the interaction
765  // cent On return, the centrality - if not available < 0
766  //
767  // Return:
768  // 0 (or kOk) on success, otherwise a bit mask of error codes
769  //
770  DGUARD(fDebug,1,"Process event in AliFMDEventInspector");
771  // --- Check that we have an event ---------------------------------
772  if (!event) {
773  AliWarning("No ESD event found for input event");
774  fHStatus->Fill(2);
775  return kNoEvent;
776  }
777 
778  // --- Process satellite event information is requested ------------
779  if (AllowDisplaced()) {
780  if (!fDisplacedVertex.Process(event))
781  AliWarning("Failed to process satellite event");
782  }
783 
784  // --- Read trigger information from the ESD and store in AOD object
785  if (!ReadTriggers(*event, triggers, nClusters)) {
786  if (fDebug > 2) {
787  AliWarning("Failed to read triggers from ESD"); }
788  fHStatus->Fill(3);
789  return kNoTriggers;
790  }
791  Bool_t trig = (triggers & kTriggerMask);
793 
794  // --- Check if this is a high-flux event --------------------------
795  const AliMultiplicity* testmult = event->GetMultiplicity();
796  if (!testmult) {
797  if (fDebug > 3) {
798  AliWarning("No central multiplicity object found"); }
799  }
800  else
801  lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
802 
803  fHType->Fill(lowFlux ? 0 : 1);
804 
805  // --- Get the interaction point -----------------------------------
806  Bool_t vzOk = ReadVertex(*event, ip);
807  fHEventsTr->Fill(ip.Z());
808  if (!vzOk) {
809  if (fDebug > 3) {
810  AliWarning("Failed to read vertex from ESD"); }
811  fHStatus->Fill(6);
812  return kNoVertex;
813  }
815 
816  // --- Read centrality information
817  cent = -10;
818  UShort_t qual = 0;
819  if (!ReadCentrality(*event, cent, qual)) {
820  if (fDebug > 3)
821  AliWarning("Failed to get centrality");
822  }
823  // --- check centrality cut
824  if (fMinCent > -0.0001 && cent < fMinCent) return kNoEvent;
825  if (fMaxCent > -0.0001 && cent > fMaxCent) return kNoEvent;
826 
827  if (cent >= 0) {
828  fHCent->Fill(cent);
829  fHCentVsQual->Fill((qual == 0 ? 197 : qual), cent);
830  }
831  fHEventsTrVtx->Fill(ip.Z());
832 
833  // --- Get the vertex bin ------------------------------------------
834  ivz = fVtxAxis.FindBin(ip.Z());
835  if (ivz <= 0 || ivz > fVtxAxis.GetNbins()) {
836  if (fDebug > 3) {
837  AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
838  ip.Z(), fVtxAxis.GetXmin(), fVtxAxis.GetXmax()));
839  }
840  ivz = 0;
841  fHStatus->Fill(7);
842  return kBadVertex;
843  }
844  fHEventsAccepted->Fill(ip.Z());
845  fHEventsAcceptedXY->Fill(ip.X(),ip.Y());
846  if (trig) fHTriggers->Fill(AliAODForwardMult::kAccepted);
847 
848  // --- Check the FMD ESD data --------------------------------------
849  if (!event->GetFMDData()) {
850  if (fDebug > 3) {
851  AliWarning("No FMD data found in ESD"); }
852  fHStatus->Fill(5);
853  return kNoFMD;
854  }
855 
856  fHStatus->Fill(1);
857  return kOk;
858 }
859 
860 //____________________________________________________________________
861 Bool_t
863  Double_t& cent,
864  UShort_t& qual) const
865 {
866  //
867  // Read centrality from event
868  //
869  // Parameters:
870  // esd Event
871  // cent On return, the centrality or negative if not found
872  //
873  // Return:
874  // False on error, true otherwise
875  //
876  DGUARD(fDebug,2,"Read the centrality in AliFMDEventInspector");
877 
878  Int_t iqual;
879  cent = AliForwardUtil::GetCentrality(esd, "V0M", iqual, (fDebug>1));
880  qual = iqual;
881  // We overwrite with satellite events, so we can be sure to get the
882  // centrality determination from the satellite vertex selection
885  qual = 0;
886  }
887 
888  return true;
889 }
890 
891 //____________________________________________________________________
892 Bool_t
894 {
895  if (fCollisionSystem != AliForwardUtil::kPPb) return true;
896 
897  AliVVZERO* esdV0 = esd.GetVZEROData();
898  if ((esdV0->GetV0ADecision()!=1) || (esdV0->GetV0CDecision()!=1))
899  return false;
900  return true;
901 }
902 
903 //____________________________________________________________________
904 Bool_t
906  UShort_t& nClusters)
907 {
908  //
909  // Read the trigger information from the ESD event
910  //
911  // Parameters:
912  // esd ESD event
913  // triggers On return, contains the trigger bits
914  //
915  // Return:
916  // @c true on success, @c false otherwise
917  //
918  DGUARD(fDebug,2,"Read the triggers in AliFMDEventInspector");
919  triggers = 0;
920 
921  // Get the analysis manager - should always be there
922  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
923  DMSG(fDebug,10,"Got analysis manager %p", am);
924  if (!am) {
925  AliWarning("No analysis manager defined!");
926  return kFALSE;
927  }
928 
929  // Get the input handler - should always be there
930  AliInputEventHandler* ih =
931  static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
932  DMSG(fDebug,10,"Got input handler %p", ih);
933  if (!ih) {
934  AliWarning("No input handler");
935  return kFALSE;
936  }
937  AliPhysicsSelection* ps =
938  dynamic_cast<AliPhysicsSelection*>(ih->GetEventSelection());
939  UInt_t trgMask = 0;
940  AliTriggerAnalysis ta;
941  if (ps) {
942  trgMask = ih->IsEventSelected();
943  }
944  else {
945  // Replay HW triggers
946  Int_t v0A = ta.V0Trigger(&esd, AliTriggerAnalysis::kASide, true, false);
947  Int_t v0C = ta.V0Trigger(&esd, AliTriggerAnalysis::kCSide, true, false);
948  Int_t fso = ta.SPDFiredChips(&esd, 1, false, 0);
949  Bool_t v0Abg = (v0A == 2);
950  Bool_t v0Cbg = (v0C == 2);
951  Bool_t v0Abb = (v0A == 1 || v0A == 3); // Include fake
952  Bool_t v0Cbb = (v0A == 1 || v0A == 3); // Include fake
953  Bool_t isBg = (v0Abg || v0Cbg);
954  if (!isBg && (v0Abb || v0Cbb || fso > 0)) trgMask |= AliVEvent::kMB;
955  if (!isBg && (v0Abb || v0Cbb)) trgMask |= AliVEvent::kCINT5;
956  if (!isBg && (v0Abb && v0Cbb)) trgMask |= AliVEvent::kINT7;
957  DMSG(fDebug, 1, "HW replay: v0A=%d v0C=%d gfo=%d mask=0x%x",
958  v0A, v0C, fso, trgMask);
959  }
960 
961  // Check if this is a collision candidate (MB)
963  // Historic remark: Note, that we should use the value cached in the
964  // input handler rather than calling IsCollisionCandiate directly
965  // on the AliPhysicsSelection obejct. If we called the latter
966  // then the AliPhysicsSelection object would overcount by a factor
967  // of 2! :-(
968  // UInt_t trgMask = ih->IsEventSelected();
969  Bool_t offline = trgMask;
970  Bool_t fastonly = (trgMask & AliVEvent::kFastOnly);
971  TString trigStr = esd.GetFiredTriggerClasses();
972 
973  if (trigStr.IsNull()) fHTrgStatus->Fill(kNoTrgWords);
974 
975  if(AllowDisplaced()) {
976  DMSG(fDebug,3,"Using displaced vertex stuff");
977  // if (TMath::Abs(fDisplacedVertex.GetVertexZ()) >= 999) offline = false;
979  triggers |= AliAODForwardMult::kSatellite;
980  }
981 
982  if (CheckFastPartition(fastonly)) {
983  fHTrgStatus->Fill(kPP2760Fast);
984  offline = false;
985  }
986 
987  if (offline && CheckCosmics(trigStr)) {
988  fHTrgStatus->Fill(kMUON);
989  offline = false;
990  }
991  if (offline) fHTrgStatus->Fill(kTriggered);
992  Int_t f = 0;
993  if (trgMask & AliVEvent::kMB) f += fHTrgStatus->Fill(kMinBias);
994  if (trgMask & AliVEvent::kCINT5) f += fHTrgStatus->Fill(kMinBiasNoSPD);
995  if (trgMask & AliVEvent::kINT7) f += fHTrgStatus->Fill(kV0AndTrg);
996  if (trgMask & AliVEvent::kHighMult) f += fHTrgStatus->Fill(kHighMult);
997  if (trgMask & AliVEvent::kCentral) f += fHTrgStatus->Fill(kCentral);
998  if (trgMask & AliVEvent::kSemiCentral) f += fHTrgStatus->Fill(kSemiCentral);
999  if (trgMask & AliVEvent::kDG5) f += fHTrgStatus->Fill(kDiffractive);
1000  if (trgMask & AliVEvent::kUserDefined) f += fHTrgStatus->Fill(kUser);
1001  if (f <= 0) fHTrgStatus->Fill(kOther);
1002 
1003  // if (!CheckpAExtraV0(esd)) offline = false;
1004 
1005  DMSG(fDebug,2,"Event is %striggered by off-line", offline ? "" : "NOT ");
1006 
1007  if (offline) triggers |= AliAODForwardMult::kOffline;
1008  // Only flag as in-elastic if a min-bias trigger is here
1009  if (trgMask & AliVEvent::kMB) {
1010  triggers |= AliAODForwardMult::kInel;
1011  CheckINELGT0(esd, nClusters, triggers);
1012  }
1013  // Flag as V0-AND if kINT7 is there
1014  if (trgMask & AliVEvent::kINT7) {
1015  triggers |= AliAODForwardMult::kV0AND;
1016  if (fUseV0AND)
1017  // If we want to use V0-AND for NSD, flag as NSD
1018  triggers |= AliAODForwardMult::kNSD;
1019  }
1020  // Check if NSD trigger fired off-line
1021  if (ta.IsOfflineTriggerFired(&esd, AliTriggerAnalysis::kNSD1))
1022  triggers |= AliAODForwardMult::kNSD;
1023 
1024  // AliTriggerAnalysis ta;
1025  if (ta.IsSPDClusterVsTrackletBG(&esd, false))
1026  // SPD outlier
1027  triggers |= AliAODForwardMult::kSPDOutlier;
1028 
1029  CheckPileup(esd, triggers);
1030  CheckEmpty(trigStr, triggers);
1031  // if (CheckPileup(esd, triggers)) fHTriggers->Fill(kPileUp+.5);
1032  // if (CheckEmpty(trigStr, triggers)) fHTriggers->Fill(kEmpty+.5);
1033 
1034  CheckWords(esd, triggers);
1035 
1036  AliAODForwardMult::FillTriggerHistogram(kTriggerMask, triggers, fHTriggers);
1037 #if 0
1038 #define TEST_TRIG_BIN(RET,BIN,TRIGGERS) \
1039  do { switch (BIN) { \
1040  case kInel: RET = triggers & AliAODForwardMult::kInel; break; \
1041  case kInelGt0: RET = triggers & AliAODForwardMult::kInelGt0; break; \
1042  case kNSD: RET = triggers & AliAODForwardMult::kNSD; break; \
1043  case kV0AND: RET = triggers & AliAODForwardMult::kV0AND; break; \
1044  case kEmpty: RET = triggers & AliAODForwardMult::kEmpty; break; \
1045  case kA: RET = triggers & AliAODForwardMult::kA; break; \
1046  case kB: RET = triggers & AliAODForwardMult::kB; break; \
1047  case kC: RET = triggers & AliAODForwardMult::kC; break; \
1048  case kE: RET = triggers & AliAODForwardMult::kE; break; \
1049  case kPileUp: RET = triggers & AliAODForwardMult::kPileUp; break; \
1050  case kMCNSD: RET = triggers & AliAODForwardMult::kMCNSD; break; \
1051  case kSatellite: RET = triggers & AliAODForwardMult::kSatellite; break; \
1052  case kSpdOutlier:RET = triggers & AliAODForwardMult::kSPDOutlier; break; \
1053  case kOffline: RET = triggers & AliAODForwardMult::kOffline; break; \
1054  default: RET = false; } } while(false)
1055 
1056  if (!fHTriggers) {
1057  AliWarning("Histogram of triggers not defined - has init been called");
1058  return false;
1059  }
1060 
1061  for (Int_t i = 0; i < kOffline+1; i++) {
1062  Bool_t hasX = false;
1063  TEST_TRIG_BIN(hasX, i, triggers);
1064  if (!hasX) continue;
1065  fHTriggers->Fill(i+.5);
1066  for (Int_t j = 0; j < kOffline+1; j++) {
1067  Bool_t hasY = false;
1068  TEST_TRIG_BIN(hasY, j, triggers);
1069  if (!hasY) continue;
1070 
1071  fHTriggerCorr->Fill(i+.5, j+.5);
1072  }
1073  }
1074 #endif
1075  return kTRUE;
1076 }
1077 
1078 //____________________________________________________________________
1079 Bool_t
1081 {
1082  // For the 2.76 TeV p+p run (LHC11a), the FMD ran in the slow partition
1083  // so it received no triggers from the fast partition. Therefore
1084  // the fast triggers are removed here but not for MC where all
1085  // triggers are fast.
1086  if (TMath::Abs(fEnergy - 2750.) > 20) return false;
1087  if (fCollisionSystem != AliForwardUtil::kPP) return false;
1088  if (fMC) return false; // All MC events for pp @ 2.76TeV are `fast'
1089  // if (fastonly)
1090  // DMSG(fDebug,1,"Fast trigger in pp @ sqrt(s)=2.76TeV removed");
1091 
1092  return fastonly;
1093 }
1094 
1095 //____________________________________________________________________
1096 Bool_t
1098 {
1099  // MUON triggers are not strictly minimum bias (MB) so they are
1100  // removed (HHD)
1101  if(trigStr.Contains("CMUS1")) {
1102  DMSG(fDebug,1,"Cosmic trigger ins't min-bias, removed");
1103  return true;
1104  }
1105  return false;
1106 }
1107 
1108 //____________________________________________________________________
1109 Bool_t
1111  UShort_t& nClusters,
1112  UInt_t& triggers) const
1113 {
1114  Int_t nTracklets = 0;
1115  nClusters = 0;
1116 
1117  // If this is inel, see if we have a tracklet
1118  const AliMultiplicity* spdmult = esd.GetMultiplicity();
1119  if (!spdmult) {
1120  AliWarning("No SPD multiplicity");
1121  return false;
1122  }
1123 
1124  // Check if we have one or more tracklets
1125  // in the range -1 < eta < 1 to set the INEL>0
1126  // trigger flag.
1127  //
1128  // Also count tracklets as a single cluster
1129  Int_t n = spdmult->GetNumberOfTracklets();
1130  for (Int_t j = 0; j < n; j++) {
1131  if(TMath::Abs(spdmult->GetEta(j)) < 1) nTracklets++;
1132  }
1133  // If we at this point had non-zero nClusters, it's INEL>0
1134  if (nTracklets > 0) triggers |= AliAODForwardMult::kInelGt0;
1135  nClusters = nTracklets;
1136 
1137  // Loop over single clusters
1138  n = spdmult->GetNumberOfSingleClusters();
1139  for (Int_t j = 0; j < n; j++) {
1140  Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
1141  if (TMath::Abs(eta) < 1) nClusters++;
1142  }
1143  if (nClusters > 0) triggers |= AliAODForwardMult::kNClusterGt0;
1144 
1145  return triggers & AliAODForwardMult::kNClusterGt0;
1146 }
1147 
1148 //____________________________________________________________________
1149 Bool_t
1151  UInt_t& triggers) const
1152 {
1153  // Check for multiple vertices (pile-up) with at least 3
1154  // contributors and at least 0.8cm from the primary vertex
1155  // if(fCollisionSystem != AliForwardUtil::kPP) return false;
1156 
1157  UInt_t locTrig = 0;
1158  // Check for standard SPD pile-up
1159  if ((fPileupFlags & kSPD) &&
1160  esd.IsPileupFromSPD(fMinPileupContrib,fMinPileupDistance)) {
1161  fHPileup->Fill(0.5, 1);
1162  locTrig |= AliAODForwardMult::kPileupSPD;
1163  }
1164 
1165  // Check for multi-vertex pileup
1166  if ((fPileupFlags & kTracks) && CheckMultiVertex(esd)) {
1167  fHPileup->Fill(1.5, 1);
1169  }
1170 
1171  // Check for out-of-bunch pileup
1172  if ((fPileupFlags & kOutOfBunch) &&
1173  (esd.GetHeader()->GetIRInt2ClosestInteractionMap() != 0 ||
1174  esd.GetHeader()->GetIRInt1ClosestInteractionMap() != 0)) {
1175  fHPileup->Fill(2.5, 1);
1176  locTrig |= AliAODForwardMult::kPileupBC;
1177  }
1178  // check SPD by multiplicity bin
1179  if ((fPileupFlags & kSPDBins) &&
1180  (esd.IsPileupFromSPDInMultBins())) {
1181  // Int_t nTracklets=GetMultiplicity()->GetNumberOfTracklets();
1182  // if (nTracklets<20) return IsPileupFromSPD(3,0.8);
1183  // else if (nTracklets<50) return IsPileupFromSPD(4,0.8);
1184  // else return IsPileupFromSPD(5,0.8);
1185 
1186  locTrig |= AliAODForwardMult::kPileupBins;
1187  fHPileup->Fill(3.5, 1);
1188  }
1189 
1190 
1191  // Our result
1192  triggers |= locTrig;
1193  Bool_t pileup = locTrig != 0;
1194  if (pileup) triggers |= AliAODForwardMult::kPileUp;
1195  return pileup;
1196 }
1197 
1198 //____________________________________________________________________
1199 Bool_t
1201  Bool_t checkOtherBC) const
1202 {
1203  // Adapted from AliAnalysisUtils
1204  //
1205  // Parameters
1206  const Double_t maxChi2nu = 5;
1207 
1208  // Number of vertices
1209  Int_t n = esd.GetNumberOfPileupVerticesTracks();
1210  if (n < 1)
1211  // No additional vertices from tracks -> no pileup
1212  return false;
1213 
1214  // Primary vertex
1215  const AliESDVertex* primary = esd.GetPrimaryVertexTracks();
1216  if (primary->GetStatus() != 1)
1217  // No good primary vertex, but many candidates -> pileup
1218  return true;
1219 
1220  // Bunch crossing number
1221  Int_t bcPrimary = primary->GetBC();
1222 
1223  for (Int_t i = 0; i < n; i++) {
1224  const AliESDVertex* other = esd.GetPileupVertexTracks(i);
1225 
1226  if (other->GetNContributors() < fMinPileupContrib)
1227  // Not enough contributors to this vertex
1228  continue;
1229 
1230  if (other->GetChi2perNDF() > maxChi2nu)
1231  // Poorly determined vertex
1232  continue;
1233  if (checkOtherBC) {
1234  Int_t bcOther = other->GetBC();
1235  if (bcOther != AliVTrack::kTOFBCNA && TMath::Abs(bcOther-bcPrimary) > 2)
1236  // Pile-up from other BC
1237  return true;
1238  }
1239  // Calculate the distance
1240  Double_t d = -1;
1241  Double_t dx = primary->GetX() - other->GetX();
1242  Double_t dy = primary->GetY() - other->GetY();
1243  Double_t dz = primary->GetZ() - other->GetZ();
1244  Double_t covPrimary[6], covOther[6];
1245  primary->GetCovarianceMatrix(covPrimary);
1246  other->GetCovarianceMatrix(covOther);
1247  TMatrixDSym v(3);
1248  v(0,0) = covPrimary[0] + covOther[0]; // diagonal
1249  v(1,1) = covPrimary[2] + covOther[2]; // diagonal
1250  v(2,2) = covPrimary[5] + covOther[5]; // diagonal
1251  v(1,0) = v(0,1) = covPrimary[1]+covOther[1]; // Band
1252  v(0,2) = v(1,2) = v(2,0) = v(2,1) = 0; // Off-diagonal+band
1253  // Try to invert
1254  v.InvertFast();
1255  if (!v.IsValid())
1256  // Question if kStatus bit is every set after InvertFast?
1257  continue;
1258  d = (v(0,0) * dx * dx + v(1,1) * dy * dy + v(2,2) * dz * dz +
1259  2 * (v(0,1) * dx * dy + v(0,2) * dx * dz + v(1,2) * dy * dz));
1260 
1261  if (d < 0 || TMath::Sqrt(d) < fMinPileupDistance)
1262  // Bad distance, or not fare enough from each other
1263  continue;
1264 
1265  // Well separated vertices -> pileup
1266  return true;
1267  }
1268  return false;
1269 }
1270 
1271 //____________________________________________________________________
1272 Bool_t
1273 AliFMDEventInspector::CheckEmpty(const TString& trigStr, UInt_t& triggers) const
1274 {
1275  if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")||
1276  trigStr.Contains("CBEAMB-B-NOPF-ALLNOTRD")) {
1277  triggers |= AliAODForwardMult::kEmpty;
1278  return true;
1279  }
1280  return false;
1281 }
1282 //____________________________________________________________________
1283 Bool_t
1285 {
1286  TObject* word = 0;
1287  TIter nextColl(&fCollWords);
1288  while ((word = nextColl())) {
1289  DMSG(fDebug,10,"Checking if %s (Coll) trigger %s is fired (%d)",
1290  word->GetTitle(), word->GetName(), word->GetUniqueID());
1291  if (!esd.IsTriggerClassFired(word->GetName())) continue;
1292  // if (fHWords) fHWords->Fill(word->GetName(), 1.);
1293  if (fHWords) fHWords->Fill(word->GetUniqueID());
1294 
1295  TString beamSide = word->GetTitle();
1296  DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1297 
1298  if (!beamSide.EqualTo("B")) continue;
1299  triggers |= AliAODForwardMult::kB;
1300 
1301  // break; // No more to do here
1302  }
1303  TIter nextBg(&fBgWords);
1304  UInt_t all = (AliAODForwardMult::kA |
1307  while ((word = nextBg())) {
1308  DMSG(fDebug,10,"Checking if %s (BG) trigger %s is fired (%d)",
1309  word->GetTitle(), word->GetName(), word->GetUniqueID());
1310  if (!esd.IsTriggerClassFired(word->GetName())) continue;
1311  // if (fHWords) fHWords->Fill(word->GetName(), 1.);
1312  if (fHWords) fHWords->Fill(word->GetUniqueID());
1313 
1314  TString beamSide = word->GetTitle();
1315  DMSG(fDebug,10,"Found it - this is a %s trigger", beamSide.Data());
1316 
1317  if (beamSide.Contains("A")) triggers |= AliAODForwardMult::kA;
1318  if (beamSide.Contains("C")) triggers |= AliAODForwardMult::kC;
1319  if (beamSide.Contains("E")) triggers |= AliAODForwardMult::kE;
1320 
1321  if ((triggers & all) == all) break; // No more to do
1322  }
1323  return true;
1324 }
1325 
1326 
1327 //____________________________________________________________________
1328 Bool_t
1330 {
1331  //
1332  // Read the vertex information from the ESD event
1333  //
1334  // Parameters:
1335  // esd ESD event
1336  // vz On return, the vertex Z position
1337  //
1338  // Return:
1339  // @c true on success, @c false otherwise
1340  //
1341  DGUARD(fDebug,2,"Read the vertex in AliFMDEventInspector");
1342  ip.SetXYZ(1024, 1024, 0);
1343 
1344  EVtxStatus s = kNoVtx;
1346  s = kVtxOK;
1347  ip.SetZ(fDisplacedVertex.GetVertexZ());
1348  }
1349  else if (fVtxMethod == kPWGUD)
1350  s = CheckPWGUDVertex(esd, ip);
1351  else if (fVtxMethod == kpA2012)
1352  s = CheckpA2012Vertex(esd,ip);
1353  else if (fVtxMethod == kpA2013)
1354  s = CheckpA2013Vertex(esd,ip);
1355  else
1356  s = CheckVertex(esd, ip);
1357 
1358  fHVtxStatus->Fill(s);
1359 
1360  return s == kVtxOK;
1361 }
1362 
1363 //____________________________________________________________________
1366  TVector3& ip) const
1367 {
1368  // This is the code used by the 1st physics people
1369  const AliESDVertex* vertex = esd.GetPrimaryVertex();
1370  if (!vertex || !vertex->GetStatus()) {
1371  DMSG(fDebug,2,"No primary vertex (%p) or bad status %d",
1372  vertex, (vertex ? vertex->GetStatus() : -1));
1373  return kNoVtx;
1374  }
1375  const AliESDVertex* vertexSPD = esd.GetPrimaryVertexSPD();
1376  if (!vertexSPD || !vertexSPD->GetStatus()) {
1377  DMSG(fDebug,2,"No primary SPD vertex (%p) or bad status %d",
1378  vertexSPD, (vertexSPD ? vertexSPD->GetStatus() : -1));
1379  return kNoSPDVtx;
1380  }
1381 
1382  // if vertex is from SPD vertexZ, require more stringent cuts
1383  if (vertex->IsFromVertexerZ()) {
1384  if (vertex->GetDispersion() > fMaxVzErr ||
1385  vertex->GetZRes() > 1.25 * fMaxVzErr) {
1386  DMSG(fDebug,2,"Dispersion %f > %f or resolution %f > %f",
1387  vertex->GetDispersion(), fMaxVzErr,
1388  vertex->GetZRes(), 1.25 * fMaxVzErr);
1389  return kUncertain;
1390  }
1391  }
1392  ip.SetZ(vertex->GetZ());
1393 
1394  if(!vertex->IsFromVertexerZ()) {
1395  ip.SetX(vertex->GetX());
1396  ip.SetY(vertex->GetY());
1397  }
1398  return kVtxOK;
1399 }
1400 //____________________________________________________________________
1403  TVector3& ip) const
1404 {
1405  const Int_t nMinContrib = 0;
1406  const AliESDVertex *vertex = esd.GetPrimaryVertexSPD();
1407  if (!vertex) return kNoSPDVtx;
1408  if (vertex->GetNContributors() <= nMinContrib) return kFewContrib;
1409 
1410  if (vertex->IsFromVertexerZ()) return kNotVtxZ;
1411 
1412  if (vertex->GetDispersion() >= 0.04 || vertex->GetZRes()>=0.25)
1413  return kUncertain;
1414 
1415  ip.SetX(vertex->GetX());
1416  ip.SetY(vertex->GetY());
1417  ip.SetZ(vertex->GetZ());
1418 
1419  return kVtxOK;
1420 }
1421 //____________________________________________________________________
1424  TVector3& ip) const
1425 {
1426  // This code is adopted from
1427  //
1428  // AliAnalysisUtils::IsVertexSelected2013pA(AliVEvent *event)
1429  //
1430  const Int_t nMinContrib = 0;
1431  const AliESDVertex* primVtx = esd.GetPrimaryVertex();
1432  if (!primVtx) return kNoVtx;
1433  if (primVtx->GetNContributors() <= nMinContrib) return kFewContrib;
1434 
1435  const AliESDVertex* spdVtx = esd.GetPrimaryVertexSPD();
1436  if (!spdVtx) return kNoSPDVtx;
1437  if (spdVtx->GetNContributors() <= nMinContrib) return kFewContrib;
1438 
1439  // TString vtxTyp = spdVtx->GetTitle();
1440  // if (spdVtx->IsFromVertexerZ()) {
1441  if (spdVtx->IsFromVertexerZ() && spdVtx->GetZRes()>0.25) return kUncertain;
1442 
1443  Bool_t correlateVtx = true;
1444  if (correlateVtx) {
1445  if (TMath::Abs(spdVtx->GetZ() - primVtx->GetZ()) > 0.5) return kUncertain;
1446  }
1447 
1448  ip.SetX(primVtx->GetX());
1449  ip.SetY(primVtx->GetY());
1450  ip.SetZ(primVtx->GetZ());
1451 
1452  return kVtxOK;
1453 }
1454 
1455 //____________________________________________________________________
1458  TVector3& ip) const
1459 {
1460  // Use standard SPD vertex (perhaps preferable for Pb+Pb)
1461  // Get the vertex
1462  const AliESDVertex* vertex = esd.GetPrimaryVertexSPD();
1463  if (!vertex) {
1464  if (fDebug > 2) {
1465  AliWarning("No SPD vertex found in ESD"); }
1466  return kNoSPDVtx;
1467  }
1468 
1469  // #if 0 // Check disabled - seem to kill a lot of PbPb events
1470  // Check that enough tracklets contributed
1471  if(vertex->GetNContributors() <= 0) {
1472  DMSG(fDebug,2,"Number of contributors to vertex is %d<=0",
1473  vertex->GetNContributors());
1474  ip.SetZ(0);
1475  return kFewContrib;
1476  }
1477  // #endif
1478 
1479  // Check that the uncertainty isn't too large
1480  if (vertex->GetZRes() > fMaxVzErr) {
1481  DMSG(fDebug,2,"Uncertaintity in Z of vertex is too large %f > %f",
1482  vertex->GetZRes(), fMaxVzErr);
1483  return kUncertain;
1484  }
1485 
1486  // Get the z coordiante
1487  ip.SetZ(vertex->GetZ());
1488  const AliESDVertex* vertexXY = esd.GetPrimaryVertex();
1489 
1490  DMSG(fDebug, 2, "Read vertex: Primary (%f,%f,%f) [%3s] SPD (%f,%f,%f) [%3s]",
1491  vertexXY->GetX(), vertexXY->GetY(), vertexXY->GetZ(),
1492  vertexXY->IsFromVertexerZ() ? "Z" : "XYZ",
1493  vertex ->GetX(), vertex ->GetY(), vertex ->GetZ(),
1494  vertex ->IsFromVertexerZ() ? "Z" : "XYZ");
1495 
1496  // If from vertexer Z, we get the mean of X and Y
1497  ip.SetX(vertexXY->GetX());
1498  ip.SetY(vertexXY->GetY());
1499  // if(!vertexXY->IsFromVertexerZ()) {
1500  // ip.SetX(vertexXY->GetX());
1501  // ip.SetY(vertexXY->GetY());
1502  // }
1503  return kVtxOK;
1504 }
1505 
1506 //____________________________________________________________________
1507 Bool_t
1509 {
1510  //
1511  // Read the collision system, collision energy, and L3 field setting
1512  // from the ESD
1513  //
1514  // Parameters:
1515  // esd ESD to get information from
1516  //
1517  // Return:
1518  // true on success, false
1519  //
1520  // AliInfo(Form("Parameters from 1st ESD event: cms=%s, sNN=%f, field=%f",
1521  // esd->GetBeamType(), 2*esd->GetBeamEnergy(),
1522  // esd->GetMagneticField()));
1523  DGUARD(fDebug,2,"Read the run details in AliFMDEventInspector");
1524  const char* sys = esd->GetBeamType();
1525  Float_t cms = 2 * esd->GetBeamEnergy();
1526  Float_t fld = esd->GetMagneticField();
1529  cms);
1531  fRunNumber = esd->GetRunNumber();
1532 
1533  Int_t a0 = esd->GetBeamParticleA(0);
1534  Int_t a1 = esd->GetBeamParticleA(1);
1535  Int_t z0 = esd->GetBeamParticleZ(0);
1536  Int_t z1 = esd->GetBeamParticleZ(1);
1537  if (fCentMethod.EqualTo("default", TString::kIgnoreCase)) {
1538  if (fCollisionSystem == 1 || fCollisionSystem == 2)
1539  SetCentralityMethod("V0M");
1540  if (fCollisionSystem == 3) {
1541  if (a0 != 0 && a1 != 0)
1542  SetCentralityMethod(a0 > a1 ? "ZNC" : "ZNA");
1543  else
1544  SetCentralityMethod("V0M");
1545  }
1546  }
1547  // store the beam species
1548  fList->Add(AliForwardUtil::MakeParameter("beam0a", a0));
1549  fList->Add(AliForwardUtil::MakeParameter("beam0z", z0));
1550  fList->Add(AliForwardUtil::MakeParameter("beam1a", a1));
1551  fList->Add(AliForwardUtil::MakeParameter("beam1z", z1));
1552 
1553  StoreProduction();
1554  StoreInformation();
1556  AliWarningF("Unknown collision system: %s - please check", sys);
1557  return false;
1558  }
1559  if (fEnergy <= 0) {
1560  AliWarningF("Unknown CMS energy: %f (%d) - please check", cms, fEnergy);
1561  return false;
1562  }
1563  if (TMath::Abs(fField) > 10) {
1564  AliWarningF("Unknown L3 field setting: %f (%d) - please check", fld,fField);
1565  return false;
1566  }
1567 
1568  return true;
1569 }
1570 
1571 
1572 //____________________________________________________________________
1573 const Char_t*
1575 {
1576  static TString s;
1577  s = "";
1578  if (code & kNoEvent) s.Append("NOEVENT ");
1579  if (code & kNoTriggers) s.Append("NOTRIGGERS ");
1580  if (code & kNoSPD) s.Append("NOSPD ");
1581  if (code & kNoFMD) s.Append("NOFMD ");
1582  if (code & kNoVertex) s.Append("NOVERTEX ");
1583  if (code & kBadVertex) s.Append("BADVERTEX ");
1584  return s.Data();
1585 }
1586 #define PF(N,V,...) \
1587  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1588 #define PFB(N,FLAG) \
1589  do { \
1590  AliForwardUtil::PrintName(N); \
1591  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1592  } while(false)
1593 #define PFV(N,VALUE) \
1594  do { \
1595  AliForwardUtil::PrintName(N); \
1596  std::cout << (VALUE) << std::endl; } while(false)
1597 
1598 //____________________________________________________________________
1599 void
1601 {
1602  //
1603  // Print information
1604  //
1605  // option Not used
1606  //
1609  sNN.Strip(TString::kBoth, '0');
1610  sNN.ReplaceAll("GeV", " GeV");
1612  field.ReplaceAll("p", "+");
1613  field.ReplaceAll("m", "-");
1614  field.ReplaceAll("kG", " kG");
1615 
1616  gROOT->IncreaseDirLevel();
1617  PFV("Production", "");
1618  gROOT->IncreaseDirLevel();
1619  PF("Period", "LHC%02d%c", fProdYear, fProdLetter);
1620  PFB("MC anchor", fProdMC);
1621  PFV("AliROOT Revision", fProdSVN);
1622  gROOT->DecreaseDirLevel();
1623  PFV("Vertex bins", fVtxAxis.GetNbins());
1624  PF("Vertex Range", "[%+6.1f,%+6.1f]", fVtxAxis.GetXmin(), fVtxAxis.GetXmax());
1625  PFV("Low flux cut", fLowFluxCut );
1626  PFV("Max(delta v_z)", fMaxVzErr );
1627  PF("Pile-up methods", "spd:%s,tracks:%s,out-of-bunch:%s",
1628  (fPileupFlags & kSPD ? "yes" : "no"),
1629  (fPileupFlags & kTracks ? "yes" : "no"),
1630  (fPileupFlags & kOutOfBunch ? "yes" : "no"));
1631  PFV("Min(nContrib_pileup)", fMinPileupContrib );
1632  PFV("Min(v-pileup)", fMinPileupDistance );
1634  PFV("CMS energy per nucleon", sNN);
1635  PFV("Field", field);
1636  PFB("Satellite events", AllowDisplaced());
1637  // PFB("Use 2012 pA vertex", fUsepA2012Vertex );
1638  // PFB("Use PWG-UD vertex", fUseFirstPhysicsVertex);
1639  TString vtxMethod("normal");
1640  switch(fVtxMethod) {
1641  case kNormal: vtxMethod = "normal" ; break;
1642  case kpA2012: vtxMethod = "pA 2012"; break;
1643  case kpA2013: vtxMethod = "pA 2013"; break;
1644  case kPWGUD: vtxMethod = "PWG-UD"; break;
1645  case kDisplaced: vtxMethod = "Satellite"; break;
1646  }
1647  PFV("IP method", vtxMethod);
1648  PFB("Simulation input", fMC );
1649  PFV("Centrality method", fCentMethod);
1650  PFV("Centrality axis", (!fCentAxis ? "none" : ""));
1651  if (!fCentAxis) {return; }
1652  Int_t nBin = fCentAxis->GetNbins();
1653  for (Int_t i = 0; i < nBin; i++) {
1654  if (i != 0 && (i % 10) == 0) std::cout << '\n';
1655  if ((i % 10) == 0) std::cout << " ";
1656  std::cout << std::setw(5) << fCentAxis->GetBinLowEdge(i+1) << '-';
1657  }
1658  std::cout << std::setw(5) << fCentAxis->GetBinUpEdge(nBin) << std::endl;
1659  gROOT->DecreaseDirLevel();
1660 }
1661 
1662 
1663 //
1664 // EOF
1665 //
1666 
static Short_t ParseMagneticField(Float_t field)
UInt_t kTriggerMask
Definition: QA.C:42
Bool_t Process(const AliESDEvent *esd)
double Double_t
Definition: External.C:58
static void FillTriggerHistogram(UInt_t triggerMask, UInt_t trg, TH1 *hist)
UShort_t fEnergy
Histogram container.
virtual EVtxStatus CheckVertex(const AliESDEvent &esd, TVector3 &ip) const
Definition: External.C:236
static const char * CenterOfMassEnergyString(UShort_t cms)
virtual Bool_t CheckPileup(const AliESDEvent &esd, UInt_t &triggers) const
TH1I * fHEventsTrVtx
Histogram of events w/trigger.
static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms)
Int_t fLowFluxCut
Pile-up status.
char Char_t
Definition: External.C:18
virtual Bool_t CheckCosmics(const TString &trigStri) const
TH1I * fHTriggers
XY vtx with trigger and Z vertex in range.
static TH1I * MakeTriggerHistogram(const char *name="triggers", UInt_t mask=0)
#define DMSG(L, N, F,...)
Int_t fCollisionSystem
void SetMaxCentrality(Double_t maxcent=-1.0)
TH1I * fHPileup
Trigger processing status.
TH2D * fHEventsAcceptedXY
Events w/trigger and vertex in range.
Bool_t ReadRunDetails(const AliESDEvent *esd)
void CacheConfiguredTriggerClasses(TList &cache, const TList *classes, AliOADBPhysicsSelection *o)
Bool_t FetchHistograms(const TList *d, TH1I *&hEventsTr, TH1I *&hEventsTrVtx, TH1I *&hEventsAcc, TH1I *&hTriggers) const
virtual Bool_t CheckFastPartition(bool fastonly) const
TH1I * fHVtxStatus
Event processing status.
UInt_t Process(const AliESDEvent *event, UInt_t &triggers, Bool_t &lowFlux, UShort_t &ivz, TVector3 &ip, Double_t &cent, UShort_t &nClusters)
Per-event per bin.
static const char * MagneticFieldString(Short_t field)
int Int_t
Definition: External.C:63
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
static const char * fgkFolderName
virtual EVtxStatus CheckpA2012Vertex(const AliESDEvent &esd, TVector3 &ip) const
Various utilities used in PWGLF/FORWARD.
void SetCentralityMethod(const TString &m)
TList fBgWords
Configured collision words.
Definition: External.C:228
Bool_t ReadVertex(const AliESDEvent &esd, TVector3 &ip)
virtual EVtxStatus CheckpA2013Vertex(const AliESDEvent &esd, TVector3 &ip) const
virtual Bool_t ReadCentrality(const AliESDEvent &esd, Double_t &cent, UShort_t &qual) const
#define PFB(N, FLAG)
TH1I * fHEventsAccepted
Events w/trigger and vertex.
void CreateOutputObjects(TList *dir)
virtual Bool_t CheckINELGT0(const AliESDEvent &esd, UShort_t &nClusters, UInt_t &triggers) const
#define DGUARD(L, N, F,...)
TH1I * fHStatus
Centrality vs quality.
static void PrintTask(const TObject &o)
static TObject * MakeParameter(const char *name, UShort_t value)
TString fCentMethod
Configured background words.
static UShort_t ParseCollisionSystem(const char *sys)
void SetMinCentrality(Double_t mincent=-1.0)
AliDisplacedVertexSelection fDisplacedVertex
TH1I * fHTrgStatus
Vertex processing status.
static Float_t GetCentrality(const AliVEvent &event, const TString &method, Int_t &qual, Bool_t verbose=false)
TH1I * fHWords
Type (low/high flux) of event.
static ULong_t AliROOTBranch()
virtual Bool_t CheckWords(const AliESDEvent &esd, UInt_t &triggers) const
void Print(Option_t *option="") const
unsigned short UShort_t
Definition: External.C:28
virtual Bool_t CheckMultiVertex(const AliESDEvent &esd, Bool_t checkOtherBC=false) const
const char Option_t
Definition: External.C:48
TH1F * fHCent
Trigger words.
static const char * CollisionSystemString(UShort_t sys)
bool Bool_t
Definition: External.C:53
#define PF(N, V,...)
void SetupForData(TList *l, const char *name=0, Bool_t mc=false)
TH2F * fHCentVsQual
Centrality.
static const char * CodeString(UInt_t mask)
#define PFV(N, VALUE)
static ULong_t AliROOTRevision()
Bool_t AllowDisplaced() const
Bool_t ReadTriggers(const AliESDEvent &esd, UInt_t &triggers, UShort_t &nClusters)
virtual Bool_t CheckEmpty(const TString &trigStr, UInt_t &triggers) const
Bool_t CheckpAExtraV0(const AliESDEvent &esd) const
virtual void SetupForData(const TAxis &vtxAxis)
virtual EVtxStatus CheckPWGUDVertex(const AliESDEvent &esd, TVector3 &ip) const