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