AliPhysics  a5cd6b6 (a5cd6b6)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakeEvaluateTriggers.C
Go to the documentation of this file.
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 //====================================================================
28 class EvaluateTrigger : public AliAnalysisTaskSE
29 {
30 public:
31  enum {
32  kNone = 0x0,
33  kESD = 0x1,
34  kMC = 0x2
35  };
36  enum {
37  kMaxN = 9
38  };
42  EvaluateTrigger()
43  : AliAnalysisTaskSE(),
44  fInel(),
45  fInelGt0(),
46  fNSD(),
47  fNClusterGt0(),
48  fInspector(),
49  fFirstEvent(true),
50  fData(0),
51  fTriggers(0),
52  fTrackletRequirement(kESD),
53  fVertexRequirement(kESD),
54  fVertexAxis(0, 0, 0),
55  fVertexESD(0),
56  fVertexMC(0),
57  fM(0)
58  {}
62  EvaluateTrigger(const char* /*name*/)
63  : AliAnalysisTaskSE("evaluateTrigger"),
64  fInel(AliAODForwardMult::kInel),
65  fInelGt0(AliAODForwardMult::kInelGt0),
66  fNSD(AliAODForwardMult::kNSD),
67  fNClusterGt0(AliAODForwardMult::kNClusterGt0),
68  fInspector("eventInspector"),
69  fFirstEvent(true),
70  fData(0),
71  fTriggers(0),
72  fTrackletRequirement(kESD),
73  fVertexRequirement(kESD),
74  fVertexAxis(10, -10, 10),
75  fVertexESD(0),
76  fVertexMC(0),
77  fM(0)
78  {
79  DefineOutput(1, TList::Class());
80  DefineOutput(2, TList::Class());
81  }
82  void SetVertexRequirement(UShort_t m) { fVertexRequirement = m; }
83  void SetTrackletRequirement(UShort_t m) { fTrackletRequirement = m; }
84  void SetVertexAxis(Int_t nBins, Double_t low, Double_t high)
85  {
86  fVertexAxis.Set(nBins, low, high);
87  }
91  void Init() {}
95  void UserCreateOutputObjects()
96  {
97  fList = new TList;
98  fList->SetOwner();
99  fList->SetName("triggerSums");
100 
101  // Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 11 };
102  // Int_t nM = 10;
103  TAxis eAxis(200, -4, 6);
104  TAxis pAxis(40, 0, 2*TMath::Pi());
105 
106  fData = new TH2D("data", "Cache",
107  eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
108  pAxis.GetNbins(), pAxis.GetXmin(), pAxis.GetXmax());
109  fData->SetDirectory(0);
110  fData->SetXTitle("#eta");
111  fData->SetYTitle("#varphi [radians]");
112  fData->SetZTitle("N_{ch}(#eta,#varphi)");
113  fData->Sumw2();
114 
115  fM = new TH1D("m", "Distribution of N_{ch}|_{|#eta|<1}", kMaxN+1,0,kMaxN+1);
116  fM->SetXTitle("N_{ch}|_{|#eta|<1}");
117  fM->SetYTitle("Events");
118  fM->SetFillColor(kRed+1);
119  fM->SetFillStyle(3001);
120  fM->SetDirectory(0);
121  fList->Add(fM);
122 
123  for (Int_t i = 0; i <= kMaxN; i++) {
124  TString lbl;
125  if (i == 0) lbl = "all";
126  else if (i == kMaxN) lbl = Form("%d+",i-1);
127  else lbl = Form("<%d",i);
128  fM->GetXaxis()->SetBinLabel(i+1, lbl);
129  }
130 
131  fTriggers = new TH1I("triggers", "Triggers", 8, -.5, 7.5);
132  fTriggers->SetDirectory(0);
133  fTriggers->GetXaxis()->SetBinLabel(1, "INEL (MC)");
134  fTriggers->GetXaxis()->SetBinLabel(2, "INEL (ESD)");
135  fTriggers->GetXaxis()->SetBinLabel(3, "INEL & N_{cluster}>0 (MC)");
136  fTriggers->GetXaxis()->SetBinLabel(4, "INEL & N_{cluster}>0 (ESD)");
137  fTriggers->GetXaxis()->SetBinLabel(5, "INEL>0 (MC)");
138  fTriggers->GetXaxis()->SetBinLabel(6, "INEL>0 (ESD)");
139  fTriggers->GetXaxis()->SetBinLabel(7, "NSD (MC)");
140  fTriggers->GetXaxis()->SetBinLabel(8, "NSD (ESD)");
141  fTriggers->SetFillColor(kYellow+1);
142  fTriggers->SetFillStyle(3001);
143  fList->Add(fTriggers);
144 
145  fVertexESD = new TH1D("vertexESD", "ESD vertex distribution",
146  fVertexAxis.GetNbins(),
147  fVertexAxis.GetXmin(),
148  fVertexAxis.GetXmax());
149  fVertexESD->SetDirectory(0);
150  fVertexESD->SetFillColor(kRed+1);
151  fVertexESD->SetFillStyle(3001);
152  fVertexESD->SetXTitle("v_{z} [cm]");
153  fVertexESD->SetYTitle("P(v_{z})");
154  fList->Add(fVertexESD);
155 
156  fVertexMC = new TH1D("vertexMC", "MC vertex distribution",
157  fVertexAxis.GetNbins(),
158  fVertexAxis.GetXmin(),
159  fVertexAxis.GetXmax());
160  fVertexMC->SetDirectory(0);
161  fVertexMC->SetFillColor(kBlue+1);
162  fVertexMC->SetFillStyle(3001);
163  fVertexMC->SetXTitle("v_{z} [cm]");
164  fVertexMC->SetYTitle("P(v_{z})");
165  fList->Add(fVertexMC);
166 
167  fInel.CreateObjects(fList, fM, fData);
168  fInelGt0.CreateObjects(fList, fM, fData);
169  fNSD.CreateObjects(fList, fM, fData);
170  fNClusterGt0.CreateObjects(fList, fM, fData);
171 
172 
173  fInspector.DefineOutput(fList);
174  fInspector.Init(fVertexAxis);
175 
176  PostData(1, fList);
177  }
181  void UserExec(Option_t*)
182  {
183  // Get the input data - MC event
184  AliMCEvent* mcEvent = MCEvent();
185  if (!mcEvent) {
186  AliWarning("No MC event found");
187  return;
188  }
189 
190  // Get the input data - ESD event
191  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
192  if (!esd) {
193  AliWarning("No ESD event found for input event");
194  return;
195  }
196 
197  if (fFirstEvent && esd->GetESDRun()) {
198  fInspector.ReadRunDetails(esd);
199 
200  AliInfo(Form("Initializing with parameters from the ESD:\n"
201  " AliESDEvent::GetBeamEnergy() ->%f\n"
202  " AliESDEvent::GetBeamType() ->%s\n"
203  " AliESDEvent::GetCurrentL3() ->%f\n"
204  " AliESDEvent::GetMagneticField()->%f\n"
205  " AliESDEvent::GetRunNumber() ->%d\n",
206  esd->GetBeamEnergy(),
207  esd->GetBeamType(),
208  esd->GetCurrentL3(),
209  esd->GetMagneticField(),
210  esd->GetRunNumber()));
211 
212  fFirstEvent = false;
213  }
214 
215  // Get the particle stack
216  AliStack* stack = mcEvent->Stack();
217 
218  // Some variables
219  UInt_t triggers; // Trigger bits
220  Bool_t lowFlux; // Low flux flag
221  UShort_t iVz; // Vertex bin from ESD
222  Double_t vZ; // Z coordinate from ESD
223  Double_t cent; // Centrality
224  UShort_t iVzMc; // Vertex bin from MC
225  Double_t vZMc; // Z coordinate of IP vertex from MC
226  Double_t b; // Impact parameter
227  Int_t nPart; // Number of participants
228  Int_t nBin; // Number of binary collisions
229  Double_t phiR; // Reaction plane from MC
230  UShort_t nClusters;// Number of clisters
231  // Process the data
232  Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent,
233  nClusters);
234  Int_t retMC = fInspector.ProcessMC(mcEvent, triggers, iVzMc,
235  vZMc, b, nPart, nBin, phiR);
236 
237  Bool_t hasESDVtx = retESD == AliFMDEventInspector::kOk;
238  Bool_t hasMCVtx = retMC == AliFMDEventInspector::kOk;
239  if (hasESDVtx) fVertexESD->Fill(vZ);
240  if (hasMCVtx) fVertexMC->Fill(vZMc);
241 
242  Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
243  Bool_t isMcNSD = (triggers & AliAODForwardMult::kMCNSD);
244 
245  Int_t mESD = 0;
246  const AliMultiplicity* spdmult = esd->GetMultiplicity();
247  if (!spdmult) {
248  AliWarning("No SPD multiplicity");
249  }
250  else {
251  // Check if we have one or more tracklets
252  // in the range -1 < eta < 1 to set the INEL>0
253  // trigger flag.
254  Int_t n = spdmult->GetNumberOfTracklets();
255  for (Int_t j = 0; j < n; j++)
256  if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
257  }
258 
259  // Reset cache
260  fData->Reset();
261  Int_t mMC = 0; // Number of particles in |eta|<1
262 
263  // Loop over all tracks
264  Int_t nTracks = mcEvent->GetNumberOfTracks();
265  for (Int_t iTr = 0; iTr < nTracks; iTr++) {
266  AliMCParticle* particle =
267  static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr));
268 
269  // Check the returned particle
270  if (!particle) continue;
271 
272  // Check if this charged and a primary
273  Bool_t isCharged = particle->Charge() != 0;
274  Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
275 
276  if (!isCharged || !isPrimary) continue;
277 
278 
279  // Fill (eta,phi) of the particle into histograsm for b
280  Double_t eta = particle->Eta();
281  Double_t phi = particle->Phi();
282 
283  fData->Fill(eta, phi);
284  if (TMath::Abs(eta) <= 1) mMC++;
285  }
286  Int_t m = mESD;
287  if (fTrackletRequirement == kMC) m = mMC;
288  fM->Fill(m);
289 
290  bool isMcInelGt0 = isMcInel && (mMC > 0);
291 
292  bool hasVertex = true;
293  if (fVertexRequirement & kMC) hasVertex = hasVertex && hasMCVtx;
294  if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;
295 
296  if (isMcInel) {
297  fTriggers->Fill(0);
298  bool triggered = (triggers & AliAODForwardMult::kInel);
299  if (triggered) fTriggers->Fill(1);
300  fInel.AddEvent(triggered, hasVertex, m, fData);
301  }
302  if (isMcInel) { // && nClusters > 0) {
303  fTriggers->Fill(2);
304  bool triggered = (triggers & AliAODForwardMult::kNClusterGt0);
305  if (triggered) fTriggers->Fill(3);
306  fNClusterGt0.AddEvent(triggered, hasVertex, m, fData);
307  }
308  if (isMcInelGt0) {
309  fTriggers->Fill(4);
310  bool triggered = (triggers & AliAODForwardMult::kInelGt0);
311  if (triggered) fTriggers->Fill(5);
312  fInelGt0.AddEvent(triggered, hasVertex, m, fData);
313  }
314  if (isMcNSD) {
315  fTriggers->Fill(6);
316  bool triggered = (triggers & AliAODForwardMult::kNSD);
317  if (triggered) fTriggers->Fill(7);
318  fNSD.AddEvent(triggered, hasVertex, m, fData);
319  }
320  PostData(1, fList);
321  }
325  void Terminate(Option_t*)
326  {
327  fList = dynamic_cast<TList*>(GetOutputData(1));
328  if (!fList) {
329  AliError(Form("No output list defined (%p)", GetOutputData(1)));
330  if (GetOutputData(1)) GetOutputData(1)->Print();
331  return;
332  }
333 
334 
335  TList* output = new TList;
336  output->SetName("triggerResults");
337  output->SetOwner();
338 
339  fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
340  fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
341  fM = static_cast<TH1D*>(fList->FindObject("m"));
342  if (fVertexMC) {
343  TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
344  vtxMC->SetDirectory(0);
345  if (vtxMC->GetEntries() > 0)
346  vtxMC->Scale(1. / vtxMC->GetEntries());
347  else
348  vtxMC->Scale(0);
349  output->Add(vtxMC);
350  }
351  if (fVertexESD) {
352  TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
353  vtxESD->SetDirectory(0);
354  if (vtxESD->GetEntries() > 0)
355  vtxESD->Scale(1. / vtxESD->GetEntries());
356  else
357  vtxESD->Scale(0);
358  output->Add(vtxESD);
359  }
360  if (fM) {
361  TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
362  m->SetDirectory(0);
363  m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
364  if (m->GetBinContent(1) > 0)
365  m->Scale(1. / m->GetBinContent(1));
366  else
367  m->Scale(0);
368  output->Add(m);
369  }
370 
371  TString vtxReq;
372  if (fVertexRequirement & kMC) vtxReq.Append("MC ");
373  if (fVertexRequirement & kESD) vtxReq.Append("ESD ");
374  output->Add(new TNamed("vtxReq", vtxReq.Data()));
375  output->Add(new TNamed("trkReq",
376  fTrackletRequirement == kMC ? "MC" : "ESD"));
377 
378  fInel.Finish(fList, output);
379  fInelGt0.Finish(fList, output);
380  fNSD.Finish(fList, output);
381  fNClusterGt0.Finish(fList, output);
382 
383  PostData(2, output);
384  }
385 
386 protected:
387  //__________________________________________________________________
391  struct TriggerType : public TNamed
392  {
393  //________________________________________________________________
397  struct MBin : public TNamed
398  {
399  TH2D* fTruth;
400  TH2D* fTriggered;
401  TH2D* fAccepted;
402  TH1I* fCounts;
403  UShort_t fLow;
404  UShort_t fHigh;
405  Bool_t IsAll() const { return fLow > fHigh; }
406  Bool_t IsLast() const { return fHigh >= kMaxN; }
410  MBin()
411  : fTruth(0), fTriggered(0), fAccepted(0),
412  fCounts(0), fLow(0), fHigh(1000) {}
422  MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist)
423  : fTruth(0),
424  fTriggered(0),
425  fAccepted(0),
426  fCounts(0),
427  fLow(low),
428  fHigh(high)
429  {
430  if (IsAll()) {
431  SetName("all");
432  SetTitle("All");
433  }
434  else if (IsLast()) {
435  SetName(Form("m%03dplus", fLow));
436  SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1}", fLow));
437  }
438  else {
439  SetName(Form("m%03d_%03d", fLow, fHigh));
440  SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", fLow,fHigh));
441  }
442 
443  TList* l = new TList;
444  l->SetOwner();
445  l->SetName(GetName());
446  p->Add(l);
447 
448  fTruth = static_cast<TH2D*>(dHist->Clone(("truth")));
449  fTruth->SetTitle("MC truth");
450  fTruth->SetDirectory(0);
451  fTruth->SetZTitle("#sum_i^{N_X} N_{ch}(#eta,#varphi)");
452  // fTruth->Sumw2();
453  fTruth->Reset();
454  l->Add(fTruth);
455 
456  fTriggered = static_cast<TH2D*>(fTruth->Clone(("triggered")));
457  fTriggered->SetTitle("Triggered");
458  fTriggered->SetDirectory(0);
459  fTriggered->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
460  // fTriggered->Sumw2();
461  fTriggered->Reset();
462  l->Add(fTriggered);
463 
464  fAccepted = static_cast<TH2D*>(fTruth->Clone(("accepted")));
465  fAccepted->SetTitle("Accepted");
466  fAccepted->SetDirectory(0);
467  fAccepted->SetZTitle("#sum_i^{N_T} N_{ch}(#eta,#varphi)");
468  // fAccepted->Sumw2();
469  fAccepted->Reset();
470  l->Add(fAccepted);
471 
472  fCounts = new TH1I("counts", "Event counts", 3, -.5, 2.5);
473  fCounts->SetDirectory(0);
474  fCounts->GetXaxis()->SetBinLabel(1, "Truth");
475  fCounts->GetXaxis()->SetBinLabel(2, "Triggered");
476  fCounts->GetXaxis()->SetBinLabel(3, "Accepted");
477  fCounts->SetYTitle("# events");
478  l->Add(fCounts);
479  }
486  void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event)
487  {
488  fCounts->Fill(0);
489  fTruth->Add(event);
490  if (triggered) {
491  fCounts->Fill(1);
492  fTriggered->Add(event);
493  if (hasVtx) {
494  fCounts->Fill(2);
495  fAccepted->Add(event);
496  }
497  }
498  }
508  Double_t Finish(const TList* p, TList* o, THStack* stack)
509  {
510  TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
511  if (!l) {
512  Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
513  return 0;
514  }
515  fTruth = static_cast<TH2D*>(l->FindObject("truth"));
516  fTriggered = static_cast<TH2D*>(l->FindObject("triggered"));
517  fAccepted = static_cast<TH2D*>(l->FindObject("accepted"));
518  fCounts = static_cast<TH1I*>(l->FindObject("counts"));
519 
520  Int_t nTruth = fCounts->GetBinContent(1);
521  Int_t nTriggered = fCounts->GetBinContent(2);
522  Int_t nAccepted = fCounts->GetBinContent(3);
523  Double_t eff = 0;
524  if (nTruth > 0) eff = double(nTriggered) / nTruth;
525  else if (nTriggered == nTruth) eff = 1;
526 
527  if (nTruth > 0) fTruth->Scale(1. / nTruth);
528  if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
529  if (nAccepted > 0) fAccepted->Scale(1. / nAccepted);
530 
531  if (IsAll())
532  Info("Finish", "%-12s [all] E_X=N_T/N_X=%9d/%-9d=%f "
533  "E_V=N_A/N_T=%9d/%-9d=%f",
534  p->GetName(), nTriggered, nTruth, eff, nAccepted, nTriggered,
535  (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
536  else if (IsLast())
537  Info("Finish", "%-12s [%2d+] E_X=N_T/N_X=%9d/%-9d=%f "
538  "E_V=N_A/N_T=%9d/%-9d=%f",
539  p->GetName(), fLow, nTriggered, nTruth, eff,
540  nAccepted, nTriggered,
541  (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
542  else
543  Info("Finish", "%-12s [%2d-%2d] E_X=N_T/N_X=%9d/%-9d=%f "
544  "E_V=N_A/N_T=%9d/%-9d=%f",
545  p->GetName(), fLow, fHigh, nTriggered, nTruth, eff,
546  nAccepted, nTriggered,
547  (nTriggered > 0 ? double(nAccepted) / nTriggered: 0));
548 
549  TList* out = new TList;
550  out->SetName(GetName());
551  out->SetOwner();
552  o->Add(out);
553 
554  out->Add(fTruth);
555  out->Add(fTriggered);
556  out->Add(new TParameter<double>("eff", eff));
557 
558  TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
559  bias->Divide(fTruth);
560  bias->SetDirectory(0);
561  bias->SetZTitle("Trigger bias (accepted/truth)");
562  out->Add(bias);
563 
564  TString title = GetTitle();
565  TH1D* truth_px = static_cast<TH1D*>(fTruth->ProjectionX("truth_eta"));
566  truth_px->SetTitle(title);
567  truth_px->Scale(1. / fTruth->GetNbinsY());
568  truth_px->SetDirectory(0);
569  truth_px->SetLineColor(kRed+1);
570  truth_px->SetMarkerColor(kRed+1);
571  truth_px->SetFillColor(kRed+1);
572  truth_px->SetFillStyle(0);
573  out->Add(truth_px);
574 
575  TH1D* triggered_px =
576  static_cast<TH1D*>(fTriggered->ProjectionX("triggered_eta"));
577  triggered_px->SetTitle(title);
578  triggered_px->Scale(1. / fTriggered->GetNbinsY());
579  triggered_px->SetDirectory(0);
580  triggered_px->SetLineColor(kGreen+1);
581  triggered_px->SetMarkerColor(kGreen+1);
582  triggered_px->SetFillColor(kGreen+1);
583  triggered_px->SetFillStyle(0);
584  out->Add(triggered_px);
585 
586  TH1D* accepted_px =
587  static_cast<TH1D*>(fAccepted->ProjectionX("accepted_eta"));
588  accepted_px->SetTitle(title);
589  accepted_px->Scale(1. / fAccepted->GetNbinsY());
590  accepted_px->SetLineColor(kBlue+1);
591  accepted_px->SetMarkerColor(kBlue+1);
592  accepted_px->SetFillColor(kBlue+1);
593  accepted_px->SetDirectory(0);
594  out->Add(accepted_px);
595 
596  THStack* data = new THStack("data", "Data distributions");
597  data->Add(truth_px);
598  data->Add(triggered_px);
599  data->Add(accepted_px);
600  out->Add(data);
601 
602 #if 0
603  TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
604  bias_px->SetTitle(title);
605  bias_px->Scale(1. / bias->GetNbinsY());
606 #else
607  TH1D* bias_px = static_cast<TH1D*>(accepted_px->Clone("bias_px"));
608  bias_px->Divide(truth_px);
609  bias_px->SetYTitle("Trigger bias (triggered/truth)");
610 #endif
611  bias_px->SetDirectory(0);
612  bias_px->SetMarkerStyle(20);
613  bias_px->SetFillStyle(0);
614  bias_px->SetMinimum(0);
615  out->Add(bias_px);
616 
617  stack->Add(bias_px);
618 
619  return eff;
620  }
621  ClassDef(MBin,1);
622  };
623 
627  TriggerType()
628  : TNamed(),
629  fMask(0),
630  fM(0),
631  fBins(0)
632  {}
633  //--- Constructor ------------------------------------------------
639  TriggerType(UShort_t mask)
640  : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
641  fMask(mask),
642  fM(0),
643  fBins(0)
644  {
645  TString n(GetName());
646  n = n.Strip(TString::kBoth);
647  SetName(n);
648  }
657  void CreateObjects(TList* list,
658  const TH1D* mHist,
659  const TH2D* dHist)
660  {
661  TList* ours = new TList;
662  ours->SetName(GetName());
663  ours->SetOwner();
664  list->Add(ours);
665 
666  fM = static_cast<TH1D*>(mHist->Clone("m"));
667  fM->SetDirectory(0);
668  ours->Add(fM);
669 
670  fBins = new TObjArray;
671  fBins->AddAt(new MBin(ours, 1000, 0, dHist), 0);
672  for (Int_t i = 1; i < fM->GetNbinsX(); i++) {
673  Double_t low = fM->GetXaxis()->GetBinLowEdge(i);
674  Double_t high = fM->GetXaxis()->GetBinUpEdge(i);
675  Info("CreatePbjects", "Adding bin [%3d,%3d] @ %d",
676  int(low), int(high), i);
677  fBins->AddAt(new MBin(ours, UShort_t(low), UShort_t(high), dHist), i);
678  }
679  }
687  MBin* FindBin(UShort_t m)
688  {
689 #if 0
690  for (Int_t i = 1; i < fBins->GetEntries(); i++) {
691  MBin* b = static_cast<MBin*>(fBins->At(i));
692  if (m >= b->fLow && m < b->fHigh) return b;
693  }
694  Warning("FindBin", "Found no bin for m=%d", m);
695  return 0;
696 #else
697  Int_t b = fM->GetXaxis()->FindBin(m);
698  if (b <= 0) return 0;
699  if (b >= fM->GetNbinsX()+1) b = fM->GetNbinsX();
700  MBin* bin = static_cast<MBin*>(fBins->At(b));
701  return bin;
702 #endif
703  }
711  void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
712  {
713  fM->AddBinContent(1);
714  fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));
715 
716  MBin* all = static_cast<MBin*>(fBins->At(0));
717  all->AddEvent(triggered, hasVtx, data);
718 
719  MBin* bin = FindBin(m);
720  if (!bin) return;
721  bin->AddEvent(triggered, hasVtx, data);
722  }
729  void Finish(const TList* p, TList* o)
730  {
731  TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
732  if (!l) {
733  Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
734  return;
735  }
736 
737  TList* ours = new TList;
738  ours->SetName(GetName());
739  ours->SetOwner();
740  o->Add(ours);
741 
742  fM = static_cast<TH1D*>(l->FindObject("m"));
743  if (!fM) {
744  Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
745  return;
746  }
747  TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
748  m->SetDirectory(0);
749  m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
750  if (m->GetBinContent(1) > 0)
751  m->Scale(1. / m->GetBinContent(1));
752  ours->Add(m);
753 
754  Int_t nBin = fM->GetNbinsX();
755  TH1D* effs = static_cast<TH1D*>(fM->Clone("effs"));
756  effs->SetYTitle("#epsilon_{X}");
757  effs->SetFillColor(kRed+1);
758  effs->SetDirectory(0);
759  effs->SetMinimum(0);
760 
761  gStyle->SetPalette(1);
762  Int_t nCol = gStyle->GetNumberOfColors();
763  THStack* stack = new THStack("biases", "Trigger biases");
764  for (Int_t i = 0; i < nBin; i++) {
765  MBin* bin = static_cast<MBin*>(fBins->At(i));
766  effs->SetBinContent(i+1, bin->Finish(l, ours, stack));
767  TH1* h = static_cast<TH1*>(stack->GetHists()->At(i));
768  Int_t col = kBlack;
769  if (i != 0) {
770  Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
771  col = gStyle->GetColorPalette(icol);
772  }
773  h->SetMarkerColor(col);
774  h->SetFillColor(col);
775  h->SetLineColor(col);
776  }
777 
778  ours->Add(stack);
779  ours->Add(effs);
780  }
781  UShort_t fMask;
782  TH1D* fM;
783  TObjArray* fBins;
784  ClassDef(TriggerType,1);
785  };
786  TriggerType fInel;
787  TriggerType fInelGt0;
788  TriggerType fNSD;
789  TriggerType fNClusterGt0;
790  AliFMDMCEventInspector fInspector;
791  TList* fList;
792  Bool_t fFirstEvent;
793  TH2D* fData;
794  TH1I* fTriggers;
795  UInt_t fTrackletRequirement;
796  UInt_t fVertexRequirement;
797  TAxis fVertexAxis;
798  TH1D* fVertexESD;
799  TH1D* fVertexMC;
800  TH1D* fM;
801  ClassDef(EvaluateTrigger,1);
802 };
803 #else
804 //====================================================================
805 void MakeEvaluateTriggers(const char* esddir,
806  Int_t nEvents = -1,
807  UInt_t vtx = 0x2,
808  UInt_t trk = 0x1,
809  Int_t vz = 10,
810  Int_t proof = 0)
811 {
812  // --- Libraries to load -------------------------------------------
813  gROOT->Macro("$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
814 
815  // --- Check for proof mode, and possibly upload pars --------------
816  if (proof> 0) {
817  gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/scripts/LoadPars.C");
818  if (!LoadPars(proof)) {
819  Error("MakeAOD", "Failed to load PARs");
820  return;
821  }
822  }
823 
824  // --- Our data chain ----------------------------------------------
825  gROOT->LoadMacro("$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/scripts/MakeChain.C");
826  TChain* chain = MakeChain("ESD", esddir,true);
827  // If 0 or less events is select, choose all
828  if (nEvents <= 0) nEvents = chain->GetEntries();
829 
830  // --- Set the macro path ------------------------------------------
831  gROOT->SetMacroPath(Form("%s:$(ALICE_PHYSICS)/PWGLF/FORWARD/analysis2:"
832  "$ALICE_PHYSICS/ANALYSIS/macros",
833  gROOT->GetMacroPath()));
834 
835  // --- Creating the manager and handlers ---------------------------
836  AliAnalysisManager *mgr = new AliAnalysisManager("Triggers",
837  "Forward multiplicity");
838  AliAnalysisManager::SetCommonFileName("triggers.root");
839 
840  // --- ESD input handler -------------------------------------------
841  AliESDInputHandler *esdHandler = new AliESDInputHandler();
842  mgr->SetInputEventHandler(esdHandler);
843 
844  // --- Monte Carlo handler -----------------------------------------
845  AliMCEventHandler* mcHandler = new AliMCEventHandler();
846  mgr->SetMCtruthEventHandler(mcHandler);
847  mcHandler->SetReadTR(true);
848 
849  // --- Add tasks ---------------------------------------------------
850  // Physics selection
851  gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,0)"));
852 
853 #if 0
854  // --- Fix up physics selection to give proper A,C, and E triggers -
855  AliInputEventHandler* ih =
856  static_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
857  AliPhysicsSelection* ps =
858  static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
859  // Ignore trigger class when selecting events. This mean that we
860  // get offline+(A,C,E) events too
861  ps->SetSkipTriggerClassSelection(true);
862 #endif
863 
864  // --- compile our code --------------------------------------------
865  gSystem->AddIncludePath("-I${ALICE_PHYSICS}/PWGLF/FORWARD/analysis2 "
866  "-I${ALICE_PHYSICS}/include "
867  "-I${ALICE_ROOT}/ANALYSIS "
868  "-I${ALICE_ROOT}/include -DBUILD=1");
869  gROOT->LoadMacro("${ALICE_PHYSICS}/PWGLF/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
870 
871  // --- Make our object ---------------------------------------------
872  EvaluateTrigger* task = new EvaluateTrigger("triggers");
873  mgr->AddTask(task);
874  task->SetVertexRequirement(vtx);
875  task->SetTrackletRequirement(trk);
876  task->SetVertexAxis(10, -vz, vz);
877 
878  // --- create containers for input/output --------------------------
879  AliAnalysisDataContainer *sums =
880  mgr->CreateContainer("triggerSums", TList::Class(),
881  AliAnalysisManager::kOutputContainer,
882  AliAnalysisManager::GetCommonFileName());
883  AliAnalysisDataContainer *output =
884  mgr->CreateContainer("triggerResults", TList::Class(),
885  AliAnalysisManager::kParamContainer,
886  AliAnalysisManager::GetCommonFileName());
887 
888  // --- connect input/output ----------------------------------------
889  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
890  mgr->ConnectOutput(task, 1, sums);
891  mgr->ConnectOutput(task, 2, output);
892 
893  // --- Run the analysis --------------------------------------------
894  TStopwatch t;
895  if (!mgr->InitAnalysis()) {
896  Error("MakeAOD", "Failed to initialize analysis train!");
897  return;
898  }
899  // Skip terminate if we're so requested and not in Proof or full mode
900  mgr->SetSkipTerminate(false);
901  // Some informative output
902  mgr->PrintStatus();
903  if (proof) mgr->SetDebugLevel(3);
904  if (mgr->GetDebugLevel() < 1 && !proof)
905  mgr->SetUseProgressBar(kTRUE,100);
906 
907  // Run the train
908  t.Start();
909  Printf("=== RUNNING ANALYSIS on %9d events ==========================",
910  nEvents);
911  mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
912  t.Stop();
913  t.Print();
914 }
915 //====================================================================
916 void DrawEvaluateTriggers(const char* filename="triggers.root")
917 {
918  TFile* file = TFile::Open(filename, "READ");
919  if (!file) {
920  Error("DrawEvaluteTriggers", "Failed to open %s", filename);
921  return;
922  }
923 
924  TList* list = static_cast<TList*>(file->Get("triggerResults"));
925  if (!list) {
926  Error("DrawEvaluateTriggers", "Faield to get triggerResults from %s",
927  list->GetName());
928  return;
929  }
930 
931  TCanvas* c = new TCanvas("c", "c", 1200, 900);
932  c->SetFillColor(0);
933  c->SetBorderSize(0);
934  c->SetBorderMode(0);
935 
936  TPad* p = new TPad("top", "Top", 0, .9, 1, 1, 0, 0, 0);
937  p->SetFillColor(0);
938  p->SetFillStyle(0);
939  p->SetBorderSize(0);
940  p->SetBorderMode(0);
941  c->cd();
942  p->Draw();
943  p->cd();
944  TObject* vtxReq = list->FindObject("vtxReq");
945  TObject* trkReq = list->FindObject("trkReq");
946  if (!vtxReq) vtxReq = new TNamed("vtxReq", "Unknown");
947  if (!trkReq) trkReq = new TNamed("trkReq", "Unknown");
948  TLatex* ltx = new TLatex(0.5, 0.5,
949  Form("Trigger bias and efficiency. "
950  "Vertex requirement: %s, "
951  "Tracklet requirement %s",
952  vtxReq->GetTitle(), trkReq->GetTitle()));
953  ltx->Draw();
954  ltx->SetTextAlign(22);
955  ltx->SetTextSize(0.4);
956 
957  const char* trigs[] = { "INEL", "NCluster>0", "INEL>0", "NSD", 0 };
958  const char** trig = trigs;
959  Int_t n = 4;
960  Int_t i = 0;
961  while (*trig) {
962  TList* lt = static_cast<TList*>(list->FindObject(*trig));
963  if (!lt) {
964  Warning("DrawEvaluteTriggers", "No list for '%s' in '%s'",
965  *trig, list->GetName());
966  // list->ls();
967  TIter next(triggerResults);
968  TObject* o ;
969  while ((o = next())) Printf("'%s'", o->GetName());
970  trig++;
971  i++;
972  continue;
973  }
974 
975  Double_t y1 = 1-double(i+1)/n;
976  Double_t y2 = 1-double(i)/n;
977  Double_t x1 = double(i)/n;
978  Double_t x2 = double(i+1)/n;
979  p = new TPad(Form("p%d", 2*i), Form("%s biases", *trig),
980  x1, .3, x2, .9, 0, 0);
981  p->SetFillColor(0);
982  p->SetFillStyle(0);
983  p->SetBorderSize(0);
984  p->SetBorderMode(0);
985  p->SetTopMargin(0.01);
986  p->SetBottomMargin(0.05);
987  p->SetRightMargin(0.01);
988  c->cd();
989  p->Draw();
990  p->cd();
991 
992  THStack* biases = static_cast<THStack*>(lt->FindObject("biases"));
993  biases->SetTitle(*trig);
994  TLegend* l = new TLegend(.1, .15, .95, .35,
995  Form("1/N_{T}#sum N_{ch}}/"
996  "1/N_{X}#sum N_{ch} - %s", *trig));
997  l->SetFillColor(0);
998  l->SetFillStyle(0);
999  l->SetBorderSize(0);
1000  l->SetNColumns(2);
1001 
1002  gStyle->SetOptTitle(0);
1003  TIter next(biases->GetHists());
1004  TH1* frame = 0;
1005  TH1* hist = 0;
1006  Int_t j = 1;
1007  while ((hist = static_cast<TH1*>(next()))) {
1008  if (!frame) {
1009 
1010  frame = new TH2D("frame", "",
1011  hist->GetXaxis()->GetNbins(),
1012  hist->GetXaxis()->GetXmin(),
1013  hist->GetXaxis()->GetXmax(),
1014  100, 0, 4);
1015  frame->SetDirectory(0);
1016  frame->SetXTitle(hist->GetXaxis()->GetTitle());
1017  frame->SetYTitle(hist->GetYaxis()->GetTitle());
1018  frame->SetStats(0);
1019  frame->Draw();
1020  }
1021  // hist->Scale(j);
1022  // hist->Draw("same p hist");
1023  // if (j == 1) hist->SetTitle("all");
1024  l->AddEntry(hist, hist->GetTitle(), "p");
1025  j++;
1026  }
1027  // biases->SetMaximum(3.5);
1028  biases->Draw("p hist nostack");
1029  l->Draw();
1030 
1031  p = new TPad(Form("p%d", 2*i+1), Form("%s efficiencies", *trig),
1032  x1, 0, x2, .3, 0, 0, 0);
1033  p->SetFillColor(0);
1034  p->SetFillStyle(0);
1035  p->SetBorderSize(0);
1036  p->SetBorderMode(0);
1037  p->SetTopMargin(0.01);
1038  p->SetRightMargin(0.01);
1039  c->cd();
1040  p->Draw();
1041  p->cd();
1042 
1043  TH1* effs = static_cast<TH1*>(lt->FindObject("effs"));
1044  effs->SetTitle(*trig);
1045  effs->SetStats(0);
1046  effs->SetMinimum(0);
1047  effs->SetMaximum(1.4);
1048  effs->SetMarkerSize(2.5);
1049  effs->Draw("hist text");
1050 
1051 
1052  i++;
1053  trig++;
1054  }
1055  c->cd();
1056 
1057  c->SaveAs("triggers.png");
1058 }
1059 
1060 #endif
filename
const char * filename
Definition: TestFCM.C:1
MakeChain
TChain * MakeChain(const char *what, const char *datadir, bool recursive=false)
Definition: MakeChain.C:151
Double_t
double Double_t
Definition: External.C:58
AliAODForwardMult::kMCNSD
Definition: AliAODForwardMult.h:129
AliAODForwardMult::kInel
Definition: AliAODForwardMult.h:111
title
const char * title
Definition: MakeQAPdf.C:26
MakeEvaluateTriggers
void MakeEvaluateTriggers(const char *esddir, Int_t nEvents=-1, UInt_t vtx=0x2, UInt_t trk=0x1, Int_t vz=10, Int_t proof=0)
Definition: MakeEvaluateTriggers.C:805
gSystem
TSystem * gSystem
list
TList * list
Definition: DrawAnaCaloTrackQA.C:18
c
TCanvas * c
Definition: TestFitELoss.C:172
kBoth
Definition: FitMassSpectra.C:39
AliFMDMCEventInspector.h
kMC
Bool_t kMC
Definition: ana.C:62
AliAODForwardMult.h
Per-event per bin.
AliAODForwardMult::kNSD
Definition: AliAODForwardMult.h:115
Int_t
int Int_t
Definition: External.C:63
TH1I
Definition: External.C:204
UInt_t
unsigned int UInt_t
Definition: External.C:33
TH2D
Definition: External.C:228
TH1D
Definition: External.C:212
DrawEvaluateTriggers
void DrawEvaluateTriggers(const char *filename="triggers.root")
Definition: MakeEvaluateTriggers.C:916
AliAODForwardMult::kNClusterGt0
Definition: AliAODForwardMult.h:133
nEvents
Float_t nEvents[nProd]
Definition: DrawProductionComparison.C:36
file
TFile * file
Definition: DrawAnaCaloTrackQA.C:19
AliAODForwardMult::kInelGt0
Definition: AliAODForwardMult.h:113
UShort_t
unsigned short UShort_t
Definition: External.C:28
Option_t
const char Option_t
Definition: External.C:48
Bool_t
bool Bool_t
Definition: External.C:53
LoadPars
Bool_t LoadPars(Int_t nWorkers=4)
Definition: LoadPars.C:13
TH1
Definition: External.C:196