]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Use new sub-algos to do hits from track-refs
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / MakeEvaluateTriggers.C
CommitLineData
5bb5d1f6 1#ifdef BUILD
2#include "AliAnalysisTaskSE.h"
3#include "AliAnalysisManager.h"
4#include "AliAnalysisDataContainer.h"
5#include "AliMCEvent.h"
6#include "AliESDEvent.h"
7#include "AliStack.h"
8#include "AliMultiplicity.h"
9#include "AliFMDMCEventInspector.h"
10#include "AliAODForwardMult.h"
11#include "AliLog.h"
12#include <TH1I.h>
13#include <TH2D.h>
14#include <TAxis.h>
15#include <TList.h>
16#include <TObjArray.h>
17#include <TParameter.h>
18#include <TStopwatch.h>
19#include <TROOT.h>
20#include <THStack.h>
21#include <TStyle.h>
22
23//====================================================================
24/**
25 * Task to evaluate trigger bias in pp
26 *
27 */
28class EvaluateTrigger : public AliAnalysisTaskSE
29{
30public:
31 enum {
32 kNone = 0x0,
33 kESD = 0x1,
34 kMC = 0x2
35 };
36 /**
37 * Constructor
38 */
39 EvaluateTrigger()
40 : AliAnalysisTaskSE(),
41 fInel(),
42 fInelGt0(),
43 fNSD(),
44 fInspector(),
45 fFirstEvent(true),
46 fData(0),
47 fTriggers(0),
48 fTrackletRequirement(kESD),
49 fVertexRequirement(kESD),
50 fVertexAxis(0, 0, 0),
51 fVertexESD(0),
52 fVertexMC(0),
53 fM(0)
54 {}
55 /**
56 * Constructor
57 */
58 EvaluateTrigger(const char* /*name*/)
59 : AliAnalysisTaskSE("evaluateTrigger"),
60 fInel(AliAODForwardMult::kInel),
61 fInelGt0(AliAODForwardMult::kInelGt0),
62 fNSD(AliAODForwardMult::kNSD),
63 fInspector("eventInspector"),
64 fFirstEvent(true),
65 fData(0),
66 fTriggers(0),
67 fTrackletRequirement(kESD),
68 fVertexRequirement(kESD),
69 fVertexAxis(10, -10, 10),
70 fVertexESD(0),
71 fVertexMC(0),
72 fM(0)
73 {
74 DefineOutput(1, TList::Class());
75 DefineOutput(2, TList::Class());
76 }
77 void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; }
78 void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; }
79 void SetVertexAxis(Int_t nBins, Double_t low, Double_t high)
80 {
81 fVertexAxis.Set(nBins, low, high);
82 }
83 /**
84 * Intialize
85 */
86 void Init() {}
87 /**
88 * Create output objects
89 */
90 void UserCreateOutputObjects()
91 {
92 fList = new TList;
93 fList->SetOwner();
94 fList->SetName(GetName());
95
96 Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 1000 };
97 Int_t nM = 10;
98 TAxis mAxis(nM, mb);
99 TAxis eAxis(200, -4, 6);
100 TAxis pAxis(40, 0, 2*TMath::Pi());
101
102 fData = new TH2D("data", "Cache",
103 eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
104 pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax());
105 fData->SetDirectory(0);
106 fData->SetXTitle("#eta");
107 fData->SetYTitle("#varphi [radians]");
108 fData->SetZTitle("N_{ch}(#eta,#varphi)");
109 fData->Sumw2();
110
111 fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}",nM+1, -0.5, nM+.5);
112 fM->SetXTitle("N_{ch}|_{|#eta|<1}");
113 fM->SetYTitle("Events");
114 fM->SetFillColor(kRed+1);
115 fM->SetFillStyle(3001);
116 fM->SetDirectory(0);
117 fList->Add(fM);
118
119 for (Int_t i = 0; i <= nM; i++) {
120 TString lbl;
121 if (i == 0) lbl = "all";
122 else if (i == nM) lbl = Form("%d+",int(mAxis.GetBinLowEdge(i)+.5));
123 else lbl = Form("<%d",int(mAxis.GetBinUpEdge(i)+.5));
124 fM->GetXaxis()->SetBinLabel(i+1, lbl);
125 }
126
127 fTriggers = new TH1I("triggers", "Triggers", 6, -.5, 5.5);
128 fTriggers->SetDirectory(0);
129 fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)");
130 fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)");
131 fTriggers->GetXaxis()->SetBinLabel(3, "INEL>0 (MC)");
132 fTriggers->GetXaxis()->SetBinLabel(4, "INEL>0 (ESD)");
133 fTriggers->GetXaxis()->SetBinLabel(5, "NSD (MC)");
134 fTriggers->GetXaxis()->SetBinLabel(6, "NSD (ESD)");
135 fTriggers->SetFillColor(kYellow+1);
136 fTriggers->SetFillStyle(3001);
137 fList->Add(fTriggers);
138
139 fVertexESD = new TH1D("vertexESD", "ESD vertex distribution",
140 fVertexAxis.GetNbins(),
141 fVertexAxis.GetXmin(),
142 fVertexAxis.GetXmax());
143 fVertexESD->SetDirectory(0);
144 fVertexESD->SetFillColor(kRed+1);
145 fVertexESD->SetFillStyle(3001);
146 fVertexESD->SetXTitle("v_{z} [cm]");
147 fVertexESD->SetYTitle("P(v_{z})");
148 fList->Add(fVertexESD);
149
150 fVertexMC = new TH1D("vertexMC", "MC vertex distribution",
151 fVertexAxis.GetNbins(),
152 fVertexAxis.GetXmin(),
153 fVertexAxis.GetXmax());
154 fVertexMC->SetDirectory(0);
155 fVertexMC->SetFillColor(kBlue+1);
156 fVertexMC->SetFillStyle(3001);
157 fVertexMC->SetXTitle("v_{z} [cm]");
158 fVertexMC->SetYTitle("P(v_{z})");
159 fList->Add(fVertexMC);
160
161 fInel.CreateObjects(fList, fM, fData);
162 fInelGt0.CreateObjects(fList, fM, fData);
163 fNSD.CreateObjects(fList, fM, fData);
164
165
166 fInspector.DefineOutput(fList);
167 fInspector.Init(fVertexAxis);
168
169 PostData(1, fList);
170 }
171 /**
172 * Event processing
173 */
174 void UserExec(Option_t*)
175 {
176 // Get the input data - MC event
177 AliMCEvent* mcEvent = MCEvent();
178 if (!mcEvent) {
179 AliWarning("No MC event found");
180 return;
181 }
182
183 // Get the input data - ESD event
184 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
185 if (!esd) {
186 AliWarning("No ESD event found for input event");
187 return;
188 }
189
190 if (fFirstEvent && esd->GetESDRun()) {
191 fInspector.ReadRunDetails(esd);
192
193 AliInfo(Form("Initializing with parameters from the ESD:\n"
194 " AliESDEvent::GetBeamEnergy() ->%f\n"
195 " AliESDEvent::GetBeamType() ->%s\n"
196 " AliESDEvent::GetCurrentL3() ->%f\n"
197 " AliESDEvent::GetMagneticField()->%f\n"
198 " AliESDEvent::GetRunNumber() ->%d\n",
199 esd->GetBeamEnergy(),
200 esd->GetBeamType(),
201 esd->GetCurrentL3(),
202 esd->GetMagneticField(),
203 esd->GetRunNumber()));
204
205 fFirstEvent = false;
206 }
207
208 // Get the particle stack
209 AliStack* stack = mcEvent->Stack();
210
211 // Some variables
212 UInt_t triggers; // Trigger bits
213 Bool_t lowFlux; // Low flux flag
214 UShort_t iVz; // Vertex bin from ESD
215 Double_t vZ; // Z coordinate from ESD
216 Double_t cent; // Centrality
217 UShort_t iVzMc; // Vertex bin from MC
218 Double_t vZMc; // Z coordinate of IP vertex from MC
219 Double_t b; // Impact parameter
220 Int_t nPart; // Number of participants
221 Int_t nBin; // Number of binary collisions
222 Double_t phiR; // Reaction plane from MC
223
224 // Process the data
225 Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent);
226 Int_t retMC = fInspector.ProcessMC(mcEvent, triggers, iVzMc,
227 vZMc, b, nPart, nBin, phiR);
228
229 Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk;
230 Bool_t hasMCVtx = retMC == AliFMDEventInspector::kOk;
231 if (hasESDVtx) fVertexESD->Fill(vZ);
232 if (hasMCVtx) fVertexMC->Fill(vZMc);
233
234 Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
235 Bool_t isMcNSD = (triggers & AliAODForwardMult::kMCNSD);
236
237 Int_t mESD = 0;
238 const AliMultiplicity* spdmult = esd->GetMultiplicity();
239 if (!spdmult) {
240 AliWarning("No SPD multiplicity");
241 }
242 else {
243 // Check if we have one or more tracklets
244 // in the range -1 < eta < 1 to set the INEL>0
245 // trigger flag.
246 Int_t n = spdmult->GetNumberOfTracklets();
247 for (Int_t j = 0; j < n; j++)
248 if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
249 }
250
251 // Reset cache
252 fData->Reset();
253 Int_t mMC = 0; // Number of particles in |eta|<1
254
255 // Loop over all tracks
256 Int_t nTracks = mcEvent->GetNumberOfTracks();
257 for (Int_t iTr = 0; iTr < nTracks; iTr++) {
258 AliMCParticle* particle =
259 static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
260
261 // Check the returned particle
262 if (!particle) continue;
263
264 // Check if this charged and a primary
265 Bool_t isCharged = particle->Charge() != 0;
266 Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
267
268 if (!isCharged || !isPrimary) continue;
269
270
271 // Fill (eta,phi) of the particle into histograsm for b
272 Double_t eta = particle->Eta();
273 Double_t phi = particle->Phi();
274
275 fData->Fill(eta, phi);
276 if (TMath::Abs(eta) <= 1) mMC++;
277 }
278 Int_t m = mESD;
279 if (fTrackletRequirement == kMC) m = mMC;
280 fM->Fill(m);
281
282 bool isMcInelGt0 = isMcInel && (mMC > 0);
283
284 bool hasVertex = true;
285 if (fVertexRequirement & kMC) hasVertex = hasVertex && hasMCVtx;
286 if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;
287
288 if (isMcInel) {
289 fTriggers->Fill(0);
290 bool triggered = (triggers & AliAODForwardMult::kInel);
291 if (triggered) fTriggers->Fill(1);
292 fInel.AddEvent(triggered, hasVertex, m, fData);
293 }
294 if (isMcInelGt0) {
295 fTriggers->Fill(2);
296 bool triggered = (triggers & AliAODForwardMult::kInelGt0);
297 if (triggered) fTriggers->Fill(3);
298 fInelGt0.AddEvent(triggered, hasVertex, m, fData);
299 }
300 if (isMcNSD) {
301 fTriggers->Fill(4);
302 bool triggered = (triggers & AliAODForwardMult::kNSD);
303 if (triggered) fTriggers->Fill(5);
304 fNSD.AddEvent(triggered, hasVertex, m, fData);
305 }
306 PostData(1, fList);
307 }
308 /**
309 * End of job processing
310 */
311 void Terminate(Option_t*)
312 {
313 fList = dynamic_cast<TList*>(GetOutputData(1));
314 if (!fList) {
315 AliError(Form("No output list defined (%p)", GetOutputData(1)));
316 if (GetOutputData(1)) GetOutputData(1)->Print();
317 return;
318 }
319
320
321 TList* output = new TList;
322 output->SetName(GetName());
323 output->SetOwner();
324
325 fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
326 fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
327 fM = static_cast<TH1D*>(fList->FindObject("m"));
328 if (fVertexMC) {
329 TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
330 vtxMC->SetDirectory(0);
331 vtxMC->Scale(1. / vtxMC->GetEntries());
332 output->Add(vtxMC);
333 }
334 if (fVertexESD) {
335 TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
336 vtxESD->SetDirectory(0);
337 vtxESD->Scale(1. / vtxESD->GetEntries());
338 output->Add(vtxESD);
339 }
340 if (fM) {
341 TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
342 m->SetDirectory(0);
343 m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
344 m->Scale(1. / m->GetBinContent(1));
345 output->Add(m);
346 }
347
348 fInel.Finish(fList, output);
349 fInelGt0.Finish(fList, output);
350 fNSD.Finish(fList, output);
351
352 PostData(2, output);
353 }
354
355protected:
356 //__________________________________________________________________
357 /**
358 * Structure to hold per trigger type information
359 */
360 struct TriggerType : public TNamed
361 {
362 //________________________________________________________________
363 /**
364 * Structure to hold per multiplicity bin information
365 */
366 struct MBin : public TNamed
367 {
368 TH2D* fTruth;
369 TH2D* fTriggered;
370 TH2D* fAccepted;
371 TH1I* fCounts;
372 UShort_t fLow;
373 UShort_t fHigh;
374 /**
375 * Constructor
376 */
377 MBin()
378 : fTruth(0), fTriggered(0), fAccepted(0),
379 fCounts(0), fLow(0), fHigh(1000) {}
380 /**
381 * Constructor
382 *
383 * @param p Parent list
384 * @param low Low cut
385 * @param high High cut
386 * @param eAxis Eta axis
387 * @param pAxis Phi axis
388 */
389 MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist)
390 : fTruth(0),
391 fTriggered(0),
392 fAccepted(0),
393 fCounts(0),
394 fLow(low),
395 fHigh(high)
396 {
397 if (low >= high) SetName("all");
398 else SetName(Form("m%03d_%03d", fLow, fHigh));
399 TList* l = new TList;
400 l->SetOwner();
401 l->SetName(GetName());
402 p->Add(l);
403
404 fTruth = static_cast<TH2D*>(dHist->Clone(("truth")));
405 fTruth->SetTitle("MC truth");
406 fTruth->SetDirectory(0);
407 fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)");
408 fTruth->Sumw2();
409 fTruth->Reset();
410 l->Add(fTruth);
411
412 fTriggered = static_cast<TH2D*>(fTruth->Clone(("triggered")));
413 fTriggered->SetTitle("Triggered");
414 fTriggered->SetDirectory(0);
415 fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
416 fTriggered->Sumw2();
417 fTriggered->Reset();
418 l->Add(fTriggered);
419
420 fAccepted = static_cast<TH2D*>(fTruth->Clone(("accepted")));
421 fAccepted->SetTitle("Accepted");
422 fAccepted->SetDirectory(0);
423 fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
424 fAccepted->Sumw2();
425 fAccepted->Reset();
426 l->Add(fAccepted);
427
428 fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5);
429 fCounts->SetDirectory(0);
430 fCounts->GetXaxis()->SetBinLabel(1, "Truth");
431 fCounts->GetXaxis()->SetBinLabel(2, "Triggered");
432 fCounts->GetXaxis()->SetBinLabel(3, "Accepted");
433 fCounts->SetYTitle("# events");
434 l->Add(fCounts);
435 }
436 /**
437 * Add event observation
438 *
439 * @param triggered Whether the event was triggered
440 * @param event Data for this event
441 */
442 void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event)
443 {
444 fCounts->Fill(0);
445 fTruth->Add(event);
446 if (triggered) {
447 fCounts->Fill(1);
448 fTriggered->Add(event);
449 if (hasVtx) {
450 fCounts->Fill(2);
451 fAccepted->Add(event);
452 }
453 }
454 }
455 /**
456 * End of job processing
457 *
458 * @param p Parent list
459 * @param o Output parent list
460 * @param stack Stack of histograms
461 *
462 * @return Trigger efficiency
463 */
464 Double_t Finish(const TList* p, TList* o, THStack* stack)
465 {
466 TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
467 if (!l) {
468 Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
469 return 0;
470 }
471 fTruth = static_cast<TH2D*>(l->FindObject("truth"));
472 fTriggered = static_cast<TH2D*>(l->FindObject("triggered"));
473 fAccepted = static_cast<TH2D*>(l->FindObject("accepted"));
474 fCounts = static_cast<TH1I*>(l->FindObject("counts"));
475
476 Int_t nTruth = fCounts->GetBinContent(1);
477 Int_t nTriggered = fCounts->GetBinContent(2);
478 Int_t nAccepted = fCounts->GetBinContent(3);
479 Double_t eff = 0;
480 if (nTruth > 0) eff = double(nTriggered) / nTruth;
481 else if (nTriggered == nTruth) eff = 1;
482
483 if (nTruth > 0) fTruth->Scale(1. / nTruth);
484 if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
485 if (nAccepted > 0) fAccepted->Scale(1. / nAccepted);
486
487 if (fLow >= fHigh)
488 Info("Finish", "%-6s [all] E_X=N_T/N_X=%9d/%-9d=%f "
489 "E_V=N_A/N_T=%9d/%-9d=%f",
490 p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered,
491 (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
492 else
493 Info("Finish", "%-6s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f "
494 "E_V=N_A/N_T=%9d/%-9d=%f",
495 p->GetName(), fLow, fHigh, nTriggered, nTruth, eff,
496 nAccepted, nTriggered,
497 (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
498
499 TList* out = new TList;
500 out->SetName(GetName());
501 out->SetOwner();
502 o->Add(out);
503
504 out->Add(fTruth);
505 out->Add(fTriggered);
506 out->Add(new TParameter<double>("eff", eff));
507
508 TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
509 bias->Divide(fTruth);
510 bias->SetDirectory(0);
511 bias->SetZTitle("Trigger bias (accepted/truth)");
512 out->Add(bias);
513
514 TH1D* truth_px = static_cast<TH1D*>(fTruth->ProjectionX("truth_eta"));
515 truth_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
516 fLow, fHigh));
517 truth_px->Scale(1. / fTruth->GetNbinsY());
518 truth_px->SetDirectory(0);
519 truth_px->SetLineColor(kRed+1);
520 truth_px->SetMarkerColor(kRed+1);
521 truth_px->SetFillColor(kRed+1);
522 truth_px->SetFillStyle(0);
523 out->Add(truth_px);
524
525 TH1D* triggered_px =
526 static_cast<TH1D*>(fTriggered->ProjectionX("triggered_eta"));
527 triggered_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
528 fLow, fHigh));
529 triggered_px->Scale(1. / fTriggered->GetNbinsY());
530 triggered_px->SetDirectory(0);
531 triggered_px->SetLineColor(kGreen+1);
532 triggered_px->SetMarkerColor(kGreen+1);
533 triggered_px->SetFillColor(kGreen+1);
534 triggered_px->SetFillStyle(0);
535 out->Add(triggered_px);
536
537 TH1D* accepted_px =
538 static_cast<TH1D*>(fAccepted->ProjectionX("accepted_eta"));
539 accepted_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
540 fLow, fHigh));
541 accepted_px->Scale(1. / fAccepted->GetNbinsY());
542 accepted_px->SetLineColor(kBlue+1);
543 accepted_px->SetMarkerColor(kBlue+1);
544 accepted_px->SetFillColor(kBlue+1);
545 accepted_px->SetDirectory(0);
546 out->Add(accepted_px);
547
548 THStack* data = new THStack("data", "Data distributions");
549 data->Add(truth_px);
550 data->Add(triggered_px);
551 data->Add(accepted_px);
552 out->Add(data);
553
554#if 0
555 TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
556 bias_px->SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d",
557 fLow, fHigh));
558 bias_px->Scale(1. / bias->GetNbinsY());
559#else
560 TH1D* bias_px = static_cast<TH1D*>(accepted_px->Clone("bias_px"));
561 bias_px->Divide(truth_px);
562 bias_px->SetYTitle("Trigger bias (triggered/truth)");
563#endif
564 bias_px->SetDirectory(0);
565 bias_px->SetMarkerStyle(20);
566 bias_px->SetFillStyle(0);
567 bias_px->SetMinimum(0);
568 out->Add(bias_px);
569
570 stack->Add(bias_px);
571
572 return eff;
573 }
574 ClassDef(MBin,1);
575 };
576
577 /**
578 * Constructor
579 */
580 TriggerType()
581 : TNamed(),
582 fMask(0),
583 fM(0),
584 fBins(0)
585 {}
586 //--- Constructor ------------------------------------------------
587 /**
588 * Constructor
589 *
590 * @param mask Trigger mask
591 */
592 TriggerType(UShort_t mask)
593 : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
594 fMask(mask),
595 fM(0),
596 fBins(0)
597 {
598 }
599 /**
600 * Create objects
601 *
602 * @param list PArent list
603 * @param mAxis Multiplicity axis
604 * @param eAxis Eta axis
605 * @param pAxis Phi axis
606 */
607 void CreateObjects(TList* list,
608 const TH1D* mHist,
609 const TH2D* dHist)
610 {
611 TList* ours = new TList;
612 ours->SetName(GetName());
613 ours->SetOwner();
614 list->Add(ours);
615
616 fM = static_cast<TH1D*>(mHist->Clone("m"));
617 fM->SetDirectory(0);
618 ours->Add(fM);
619
620 fBins = new TObjArray;
621 fBins->AddAt(new MBin(ours, 0, 0, dHist), 0);
622 for (Int_t i = 1; i <= fM->GetNbinsX(); i++) {
623 Double_t low = fM->GetXaxis()->GetBinLowEdge(i);
624 Double_t high = fM->GetXaxis()->GetBinUpEdge(i);
625
626 fBins->AddAt(new MBin(ours, low, high, dHist), i);
627 }
628 }
629 /**
630 * Find bin corresponding to m
631 *
632 * @param m Multiplicity
633 *
634 * @return Bin.
635 */
636 MBin* FindBin(UShort_t m)
637 {
638 Int_t b = fM->GetXaxis()->FindBin(m);
639 if (b <= 0) return 0;
640 if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX();
641 return static_cast<MBin*>(fBins->At(b));
642 }
643 /**
644 * Add event observation
645 *
646 * @param triggered IF this is triggered
647 * @param m Multiplicity
648 * @param data Observation
649 */
650 void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
651 {
652 fM->AddBinContent(1);
653 fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));
654
655 MBin* all = static_cast<MBin*>(fBins->At(0));
656 all->AddEvent(triggered, hasVtx, data);
657
658 MBin* bin = FindBin(m);
659 bin->AddEvent(triggered, hasVtx, data);
660 }
661 /**
662 * End of job processing
663 *
664 * @param p Parent list
665 * @param o Parent output list
666 */
667 void Finish(const TList* p, TList* o)
668 {
669 TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
670 if (!l) {
671 Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
672 return;
673 }
674
675 TList* ours = new TList;
676 ours->SetName(GetName());
677 ours->SetOwner();
678 o->Add(ours);
679
680 fM = static_cast<TH1D*>(l->FindObject("m"));
681 if (!fM) {
682 Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
683 return;
684 }
685 TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
686 m->SetDirectory(0);
687 m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
688 if (m->GetBinContent(1) > 0)
689 m->Scale(1. / m->GetBinContent(1));
690 ours->Add(m);
691
692 Int_t nBin = fM->GetNbinsX();
693 TH1D* effs = static_cast<TH1D*>(fM->Clone("effs"));
694 effs->SetYTitle("#epsilon_{X}");
695 effs->SetFillColor(kRed+1);
696 effs->SetDirectory(0);
697 effs->SetMinimum(0);
698
699 gStyle->SetPalette(1);
700 Int_t nCol = gStyle->GetNumberOfColors();
701 THStack* stack = new THStack("biases", "Trigger biases");
702 for (Int_t i = 0; i <= nBin; i++) {
703 MBin* bin = static_cast<MBin*>(fBins->At(i));
704 effs->SetBinContent(i+1, bin->Finish(l, ours, stack));
705 TH1* h = static_cast<TH1*>(stack->GetHists()->At(i));
706 Int_t col = kBlack;
707 if (i != 0) {
708 Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
709 col = gStyle->GetColorPalette(icol);
710 }
711 h->SetMarkerColor(col);
712 h->SetFillColor(col);
713 h->SetLineColor(col);
714 }
715
716 ours->Add(stack);
717 ours->Add(effs);
718 }
719 UShort_t fMask;
720 TH1D* fM;
721 TObjArray* fBins;
722 ClassDef(TriggerType,1);
723 };
724 TriggerType fInel;
725 TriggerType fInelGt0;
726 TriggerType fNSD;
727 AliFMDMCEventInspector fInspector;
728 TList* fList;
729 Bool_t fFirstEvent;
730 TH2D* fData;
731 TH1I* fTriggers;
732 UInt_t fTrackletRequirement;
733 UInt_t fVertexRequirement;
734 TAxis fVertexAxis;
735 TH1D* fVertexESD;
736 TH1D* fVertexMC;
737 TH1D* fM;
738 ClassDef(EvaluateTrigger,1);
739};
740#else
741//====================================================================
742void MakeEvaluateTriggers(const char* esddir,
743 Int_t nEvents = -1,
744 UInt_t vtx = 0x1,
745 UInt_t trk = 0x1,
746 UInt_t vz = 10,
747 Int_t proof = 0)
748{
749 // --- Libraries to load -------------------------------------------
750 gROOT->Macro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C");
751
752 // --- Check for proof mode, and possibly upload pars --------------
753 if (proof> 0) {
754 gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadPars.C");
755 if (!LoadPars(proof)) {
756 Error("MakeAOD", "Failed to load PARs");
757 return;
758 }
759 }
760
761 // --- Our data chain ----------------------------------------------
762 gROOT->LoadMacro("$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/MakeChain.C");
763 TChain* chain = MakeChain("ESD", esddir,true);
764 // If 0 or less events is select, choose all
765 if (nEvents <= 0) nEvents = chain->GetEntries();
766
767 // --- Set the macro path ------------------------------------------
768 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWG2/FORWARD/analysis2:"
769 "$ALICE_ROOT/ANALYSIS/macros",
770 gROOT->GetMacroPath()));
771
772 // --- Creating the manager and handlers ---------------------------
773 AliAnalysisManager *mgr = new AliAnalysisManager("Triggers",
774 "Forward multiplicity");
775 AliAnalysisManager::SetCommonFileName("triggers.root");
776
777 // --- ESD input handler -------------------------------------------
778 AliESDInputHandler *esdHandler = new AliESDInputHandler();
779 mgr->SetInputEventHandler(esdHandler);
780
781 // --- Monte Carlo handler -----------------------------------------
782 AliMCEventHandler* mcHandler = new AliMCEventHandler();
783 mgr->SetMCtruthEventHandler(mcHandler);
784 mcHandler->SetReadTR(true);
785
786 // --- Add tasks ---------------------------------------------------
787 // Physics selection
788 gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)"));
789
790#if 0
791 // --- Fix up physics selection to give proper A,C, and E triggers -
792 AliInputEventHandler* ih =
793 static_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
794 AliPhysicsSelection* ps =
795 static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
796 // Ignore trigger class when selecting events. This mean that we
797 // get offline+(A,C,E) events too
798 ps->SetSkipTriggerClassSelection(true);
799#endif
800
801 // --- compile our code --------------------------------------------
802 gSystem->AddIncludePath("-I${ALICE_ROOT}/PWG2/FORWARD/analysis2 "
803 "-I${ALICE_ROOT}/ANALYSIS "
804 "-I${ALICE_ROOT}/include -DBUILD=1");
805 gROOT->LoadMacro("${ALICE_ROOT}/PWG2/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
806
807 // --- Make our object ---------------------------------------------
808 EvaluateTrigger* task = new EvaluateTrigger("triggers");
809 mgr->AddTask(task);
810 task->SetVertexRequirement(vtx);
811 task->SetTrackletRequirement(trk);
812 task->SetVertexAxis(10, -vz, vz);
813
814 // --- create containers for input/output --------------------------
815 AliAnalysisDataContainer *sums =
816 mgr->CreateContainer("triggerSums", TList::Class(),
817 AliAnalysisManager::kOutputContainer,
818 AliAnalysisManager::GetCommonFileName());
819 AliAnalysisDataContainer *output =
820 mgr->CreateContainer("triggerResults", TList::Class(),
821 AliAnalysisManager::kParamContainer,
822 AliAnalysisManager::GetCommonFileName());
823
824 // --- connect input/output ----------------------------------------
825 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
826 mgr->ConnectOutput(task, 1, sums);
827 mgr->ConnectOutput(task, 2, output);
828
829 // --- Run the analysis --------------------------------------------
830 TStopwatch t;
831 if (!mgr->InitAnalysis()) {
832 Error("MakeAOD", "Failed to initialize analysis train!");
833 return;
834 }
835 // Skip terminate if we're so requested and not in Proof or full mode
836 mgr->SetSkipTerminate(false);
837 // Some informative output
838 mgr->PrintStatus();
839 if (proof) mgr->SetDebugLevel(3);
840 if (mgr->GetDebugLevel() < 1 && !proof)
841 mgr->SetUseProgressBar(kTRUE,100);
842
843 // Run the train
844 t.Start();
845 Printf("=== RUNNING ANALYSIS on %9d events ==========================",
846 nEvents);
847 mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
848 t.Stop();
849 t.Print();
850}
851#endif