2 #include "AliAnalysisTaskSE.h"
3 #include "AliAnalysisManager.h"
4 #include "AliAnalysisDataContainer.h"
5 #include "AliMCEvent.h"
6 #include "AliESDEvent.h"
8 #include "AliMultiplicity.h"
9 #include "AliFMDMCEventInspector.h"
10 #include "AliAODForwardMult.h"
16 #include <TObjArray.h>
17 #include <TParameter.h>
18 #include <TStopwatch.h>
23 //====================================================================
25 * Task to evaluate trigger bias in pp
28 class EvaluateTrigger : public AliAnalysisTaskSE
43 : AliAnalysisTaskSE(),
52 fTrackletRequirement(kESD),
53 fVertexRequirement(kESD),
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"),
72 fTrackletRequirement(kESD),
73 fVertexRequirement(kESD),
74 fVertexAxis(10, -10, 10),
79 DefineOutput(1, TList::Class());
80 DefineOutput(2, TList::Class());
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)
86 fVertexAxis.Set(nBins, low, high);
93 * Create output objects
95 void UserCreateOutputObjects()
99 fList->SetName("triggerSums");
101 // Double_t mb[] = { 0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 11 };
103 TAxis eAxis(200, -4, 6);
104 TAxis pAxis(40, 0, 2*TMath::Pi());
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)");
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);
123 for (Int_t i = 0; i <= kMaxN; i++) {
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);
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);
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);
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);
167 fInel.CreateObjects(fList, fM, fData);
168 fInelGt0.CreateObjects(fList, fM, fData);
169 fNSD.CreateObjects(fList, fM, fData);
170 fNClusterGt0.CreateObjects(fList, fM, fData);
173 fInspector.DefineOutput(fList);
174 fInspector.Init(fVertexAxis);
181 void UserExec(Option_t*)
183 // Get the input data - MC event
184 AliMCEvent* mcEvent = MCEvent();
186 AliWarning("No MC event found");
190 // Get the input data - ESD event
191 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
193 AliWarning("No ESD event found for input event");
197 if (fFirstEvent && esd->GetESDRun()) {
198 fInspector.ReadRunDetails(esd);
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(),
209 esd->GetMagneticField(),
210 esd->GetRunNumber()));
215 // Get the particle stack
216 AliStack* stack = mcEvent->Stack();
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
232 Int_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, cent,
234 Int_t retMC = fInspector.ProcessMC(mcEvent, triggers, iVzMc,
235 vZMc, b, nPart, nBin, phiR);
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);
242 Bool_t isMcInel = true; // (triggers & AliAODForwardMult::kB);
243 Bool_t isMcNSD = (triggers & AliAODForwardMult::kMCNSD);
246 const AliMultiplicity* spdmult = esd->GetMultiplicity();
248 AliWarning("No SPD multiplicity");
251 // Check if we have one or more tracklets
252 // in the range -1 < eta < 1 to set the INEL>0
254 Int_t n = spdmult->GetNumberOfTracklets();
255 for (Int_t j = 0; j < n; j++)
256 if(TMath::Abs(spdmult->GetEta(j)) < 1) mESD++;
261 Int_t mMC = 0; // Number of particles in |eta|<1
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));
269 // Check the returned particle
270 if (!particle) continue;
272 // Check if this charged and a primary
273 Bool_t isCharged = particle->Charge() != 0;
274 Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
276 if (!isCharged || !isPrimary) continue;
279 // Fill (eta,phi) of the particle into histograsm for b
280 Double_t eta = particle->Eta();
281 Double_t phi = particle->Phi();
283 fData->Fill(eta, phi);
284 if (TMath::Abs(eta) <= 1) mMC++;
287 if (fTrackletRequirement == kMC) m = mMC;
290 bool isMcInelGt0 = isMcInel && (mMC > 0);
292 bool hasVertex = true;
293 if (fVertexRequirement & kMC) hasVertex = hasVertex && hasMCVtx;
294 if (fVertexRequirement & kESD) hasVertex = hasVertex && hasESDVtx;
298 bool triggered = (triggers & AliAODForwardMult::kInel);
299 if (triggered) fTriggers->Fill(1);
300 fInel.AddEvent(triggered, hasVertex, m, fData);
302 if (isMcInel) { // && nClusters > 0) {
304 bool triggered = (triggers & AliAODForwardMult::kNClusterGt0);
305 if (triggered) fTriggers->Fill(3);
306 fNClusterGt0.AddEvent(triggered, hasVertex, m, fData);
310 bool triggered = (triggers & AliAODForwardMult::kInelGt0);
311 if (triggered) fTriggers->Fill(5);
312 fInelGt0.AddEvent(triggered, hasVertex, m, fData);
316 bool triggered = (triggers & AliAODForwardMult::kNSD);
317 if (triggered) fTriggers->Fill(7);
318 fNSD.AddEvent(triggered, hasVertex, m, fData);
323 * End of job processing
325 void Terminate(Option_t*)
327 fList = dynamic_cast<TList*>(GetOutputData(1));
329 AliError(Form("No output list defined (%p)", GetOutputData(1)));
330 if (GetOutputData(1)) GetOutputData(1)->Print();
335 TList* output = new TList;
336 output->SetName("triggerResults");
339 fVertexMC = static_cast<TH1D*>(fList->FindObject("vertexMC"));
340 fVertexESD = static_cast<TH1D*>(fList->FindObject("vertexESD"));
341 fM = static_cast<TH1D*>(fList->FindObject("m"));
343 TH1D* vtxMC = static_cast<TH1D*>(fVertexMC->Clone("vertexMC"));
344 vtxMC->SetDirectory(0);
345 if (vtxMC->GetEntries() > 0)
346 vtxMC->Scale(1. / vtxMC->GetEntries());
352 TH1D* vtxESD = static_cast<TH1D*>(fVertexESD->Clone("vertexESD"));
353 vtxESD->SetDirectory(0);
354 if (vtxESD->GetEntries() > 0)
355 vtxESD->Scale(1. / vtxESD->GetEntries());
361 TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
363 m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
364 if (m->GetBinContent(1) > 0)
365 m->Scale(1. / m->GetBinContent(1));
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"));
378 fInel.Finish(fList, output);
379 fInelGt0.Finish(fList, output);
380 fNSD.Finish(fList, output);
381 fNClusterGt0.Finish(fList, output);
387 //__________________________________________________________________
389 * Structure to hold per trigger type information
391 struct TriggerType : public TNamed
393 //________________________________________________________________
395 * Structure to hold per multiplicity bin information
397 struct MBin : public TNamed
405 Bool_t IsAll() const { return fLow > fHigh; }
406 Bool_t IsLast() const { return fHigh >= kMaxN; }
411 : fTruth(0), fTriggered(0), fAccepted(0),
412 fCounts(0), fLow(0), fHigh(1000) {}
416 * @param p Parent list
418 * @param high High cut
419 * @param eAxis Eta axis
420 * @param pAxis Phi axis
422 MBin(TList* p, UShort_t low, UShort_t high, const TH2D* dHist)
435 SetName(Form("m%03dplus", fLow));
436 SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1}", fLow));
439 SetName(Form("m%03d_%03d", fLow, fHigh));
440 SetTitle(Form("%d #leq N_{tracklets}|_{|#eta|<1} < %d", fLow,fHigh));
443 TList* l = new TList;
445 l->SetName(GetName());
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)");
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();
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();
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");
481 * Add event observation
483 * @param triggered Whether the event was triggered
484 * @param event Data for this event
486 void AddEvent(Bool_t triggered, Bool_t hasVtx, const TH2D* event)
492 fTriggered->Add(event);
495 fAccepted->Add(event);
500 * End of job processing
502 * @param p Parent list
503 * @param o Output parent list
504 * @param stack Stack of histograms
506 * @return Trigger efficiency
508 Double_t Finish(const TList* p, TList* o, THStack* stack)
510 TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
512 Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
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"));
520 Int_t nTruth = fCounts->GetBinContent(1);
521 Int_t nTriggered = fCounts->GetBinContent(2);
522 Int_t nAccepted = fCounts->GetBinContent(3);
524 if (nTruth > 0) eff = double(nTriggered) / nTruth;
525 else if (nTriggered == nTruth) eff = 1;
527 if (nTruth > 0) fTruth->Scale(1. / nTruth);
528 if (nTriggered > 0) fTriggered->Scale(1. / nTriggered);
529 if (nAccepted > 0) fAccepted->Scale(1. / nAccepted);
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));
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));
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));
549 TList* out = new TList;
550 out->SetName(GetName());
555 out->Add(fTriggered);
556 out->Add(new TParameter<double>("eff", eff));
558 TH2D* bias = static_cast<TH2D*>(fAccepted->Clone("bias"));
559 bias->Divide(fTruth);
560 bias->SetDirectory(0);
561 bias->SetZTitle("Trigger bias (accepted/truth)");
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);
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);
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);
596 THStack* data = new THStack("data", "Data distributions");
598 data->Add(triggered_px);
599 data->Add(accepted_px);
603 TH1D* bias_px = static_cast<TH1D*>(bias->ProjectionX("bias_eta"));
604 bias_px->SetTitle(title);
605 bias_px->Scale(1. / bias->GetNbinsY());
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)");
611 bias_px->SetDirectory(0);
612 bias_px->SetMarkerStyle(20);
613 bias_px->SetFillStyle(0);
614 bias_px->SetMinimum(0);
633 //--- Constructor ------------------------------------------------
637 * @param mask Trigger mask
639 TriggerType(UShort_t mask)
640 : TNamed(AliAODForwardMult::GetTriggerString(mask), ""),
645 TString n(GetName());
646 n = n.Strip(TString::kBoth);
652 * @param list PArent list
653 * @param mAxis Multiplicity axis
654 * @param eAxis Eta axis
655 * @param pAxis Phi axis
657 void CreateObjects(TList* list,
661 TList* ours = new TList;
662 ours->SetName(GetName());
666 fM = static_cast<TH1D*>(mHist->Clone("m"));
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);
681 * Find bin corresponding to m
683 * @param m Multiplicity
687 MBin* FindBin(UShort_t m)
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;
694 Warning("FindBin", "Found no bin for m=%d", m);
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));
705 * Add event observation
707 * @param triggered IF this is triggered
708 * @param m Multiplicity
709 * @param data Observation
711 void AddEvent(Bool_t triggered, Bool_t hasVtx, UShort_t m, const TH2D* data)
713 fM->AddBinContent(1);
714 fM->AddBinContent(TMath::Min(fM->GetNbinsX(), m+2));
716 MBin* all = static_cast<MBin*>(fBins->At(0));
717 all->AddEvent(triggered, hasVtx, data);
719 MBin* bin = FindBin(m);
721 bin->AddEvent(triggered, hasVtx, data);
724 * End of job processing
726 * @param p Parent list
727 * @param o Parent output list
729 void Finish(const TList* p, TList* o)
731 TList* l = dynamic_cast<TList*>(p->FindObject(GetName()));
733 Warning("Finish", "Cannot find %s in %s", GetName(), p->GetName());
737 TList* ours = new TList;
738 ours->SetName(GetName());
742 fM = static_cast<TH1D*>(l->FindObject("m"));
744 Warning("Finish", "Didn't find object 'm' in %s", l->GetName());
747 TH1D* m = static_cast<TH1D*>(fM->Clone("m"));
749 m->SetYTitle("P(N_{ch}|_{|#eta|<1} < X)");
750 if (m->GetBinContent(1) > 0)
751 m->Scale(1. / m->GetBinContent(1));
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);
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));
770 Int_t icol = TMath::Min(nCol-1,int(double(i)/nBin * nCol + .5));
771 col = gStyle->GetColorPalette(icol);
773 h->SetMarkerColor(col);
774 h->SetFillColor(col);
775 h->SetLineColor(col);
784 ClassDef(TriggerType,1);
787 TriggerType fInelGt0;
789 TriggerType fNClusterGt0;
790 AliFMDMCEventInspector fInspector;
795 UInt_t fTrackletRequirement;
796 UInt_t fVertexRequirement;
801 ClassDef(EvaluateTrigger,1);
804 //====================================================================
805 void MakeEvaluateTriggers(const char* esddir,
812 // --- Libraries to load -------------------------------------------
813 gROOT->Macro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C");
815 // --- Check for proof mode, and possibly upload pars --------------
817 gROOT->LoadMacro("$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadPars.C");
818 if (!LoadPars(proof)) {
819 Error("MakeAOD", "Failed to load PARs");
824 // --- Our data chain ----------------------------------------------
825 gROOT->LoadMacro("$ALICE_ROOT/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();
830 // --- Set the macro path ------------------------------------------
831 gROOT->SetMacroPath(Form("%s:$(ALICE_ROOT)/PWGLF/FORWARD/analysis2:"
832 "$ALICE_ROOT/ANALYSIS/macros",
833 gROOT->GetMacroPath()));
835 // --- Creating the manager and handlers ---------------------------
836 AliAnalysisManager *mgr = new AliAnalysisManager("Triggers",
837 "Forward multiplicity");
838 AliAnalysisManager::SetCommonFileName("triggers.root");
840 // --- ESD input handler -------------------------------------------
841 AliESDInputHandler *esdHandler = new AliESDInputHandler();
842 mgr->SetInputEventHandler(esdHandler);
844 // --- Monte Carlo handler -----------------------------------------
845 AliMCEventHandler* mcHandler = new AliMCEventHandler();
846 mgr->SetMCtruthEventHandler(mcHandler);
847 mcHandler->SetReadTR(true);
849 // --- Add tasks ---------------------------------------------------
851 gROOT->Macro(Form("AddTaskPhysicsSelection.C(1,1,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);
864 // --- compile our code --------------------------------------------
865 gSystem->AddIncludePath("-I${ALICE_ROOT}/PWGLF/FORWARD/analysis2 "
866 "-I${ALICE_ROOT}/ANALYSIS "
867 "-I${ALICE_ROOT}/include -DBUILD=1");
868 gROOT->LoadMacro("${ALICE_ROOT}/PWGLF/FORWARD/analysis2/MakeEvaluateTriggers.C++g");
870 // --- Make our object ---------------------------------------------
871 EvaluateTrigger* task = new EvaluateTrigger("triggers");
873 task->SetVertexRequirement(vtx);
874 task->SetTrackletRequirement(trk);
875 task->SetVertexAxis(10, -vz, vz);
877 // --- create containers for input/output --------------------------
878 AliAnalysisDataContainer *sums =
879 mgr->CreateContainer("triggerSums", TList::Class(),
880 AliAnalysisManager::kOutputContainer,
881 AliAnalysisManager::GetCommonFileName());
882 AliAnalysisDataContainer *output =
883 mgr->CreateContainer("triggerResults", TList::Class(),
884 AliAnalysisManager::kParamContainer,
885 AliAnalysisManager::GetCommonFileName());
887 // --- connect input/output ----------------------------------------
888 mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
889 mgr->ConnectOutput(task, 1, sums);
890 mgr->ConnectOutput(task, 2, output);
892 // --- Run the analysis --------------------------------------------
894 if (!mgr->InitAnalysis()) {
895 Error("MakeAOD", "Failed to initialize analysis train!");
898 // Skip terminate if we're so requested and not in Proof or full mode
899 mgr->SetSkipTerminate(false);
900 // Some informative output
902 if (proof) mgr->SetDebugLevel(3);
903 if (mgr->GetDebugLevel() < 1 && !proof)
904 mgr->SetUseProgressBar(kTRUE,100);
908 Printf("=== RUNNING ANALYSIS on %9d events ==========================",
910 mgr->StartAnalysis(proof ? "proof" : "local", chain, nEvents);
914 //====================================================================
915 void DrawEvaluateTriggers(const char* filename="triggers.root")
917 TFile* file = TFile::Open(filename, "READ");
919 Error("DrawEvaluteTriggers", "Failed to open %s", filename);
923 TList* list = static_cast<TList*>(file->Get("triggerResults"));
925 Error("DrawEvaluateTriggers", "Faield to get triggerResults from %s",
930 TCanvas* c = new TCanvas("c", "c", 1200, 900);
935 TPad* p = new TPad("top", "Top", 0, .9, 1, 1, 0, 0, 0);
943 TObject* vtxReq = list->FindObject("vtxReq");
944 TObject* trkReq = list->FindObject("trkReq");
945 if (!vtxReq) vtxReq = new TNamed("vtxReq", "Unknown");
946 if (!trkReq) trkReq = new TNamed("trkReq", "Unknown");
947 TLatex* ltx = new TLatex(0.5, 0.5,
948 Form("Trigger bias and efficiency. "
949 "Vertex requirement: %s, "
950 "Tracklet requirement %s",
951 vtxReq->GetTitle(), trkReq->GetTitle()));
953 ltx->SetTextAlign(22);
954 ltx->SetTextSize(0.4);
956 const char* trigs[] = { "INEL", "NCluster>0", "INEL>0", "NSD", 0 };
957 const char** trig = trigs;
961 TList* lt = static_cast<TList*>(list->FindObject(*trig));
963 Warning("DrawEvaluteTriggers", "No list for '%s' in '%s'",
964 *trig, list->GetName());
966 TIter next(triggerResults);
968 while ((o = next())) Printf("'%s'", o->GetName());
974 Double_t y1 = 1-double(i+1)/n;
975 Double_t y2 = 1-double(i)/n;
976 Double_t x1 = double(i)/n;
977 Double_t x2 = double(i+1)/n;
978 p = new TPad(Form("p%d", 2*i), Form("%s biases", *trig),
979 x1, .3, x2, .9, 0, 0);
984 p->SetTopMargin(0.01);
985 p->SetBottomMargin(0.05);
986 p->SetRightMargin(0.01);
991 THStack* biases = static_cast<THStack*>(lt->FindObject("biases"));
992 biases->SetTitle(*trig);
993 TLegend* l = new TLegend(.1, .15, .95, .35,
994 Form("1/N_{T}#sum N_{ch}}/"
995 "1/N_{X}#sum N_{ch} - %s", *trig));
1001 gStyle->SetOptTitle(0);
1002 TIter next(biases->GetHists());
1006 while ((hist = static_cast<TH1*>(next()))) {
1009 frame = new TH2D("frame", "",
1010 hist->GetXaxis()->GetNbins(),
1011 hist->GetXaxis()->GetXmin(),
1012 hist->GetXaxis()->GetXmax(),
1014 frame->SetDirectory(0);
1015 frame->SetXTitle(hist->GetXaxis()->GetTitle());
1016 frame->SetYTitle(hist->GetYaxis()->GetTitle());
1021 // hist->Draw("same p hist");
1022 // if (j == 1) hist->SetTitle("all");
1023 l->AddEntry(hist, hist->GetTitle(), "p");
1026 // biases->SetMaximum(3.5);
1027 biases->Draw("p hist nostack");
1030 p = new TPad(Form("p%d", 2*i+1), Form("%s efficiencies", *trig),
1031 x1, 0, x2, .3, 0, 0, 0);
1034 p->SetBorderSize(0);
1035 p->SetBorderMode(0);
1036 p->SetTopMargin(0.01);
1037 p->SetRightMargin(0.01);
1042 TH1* effs = static_cast<TH1*>(lt->FindObject("effs"));
1043 effs->SetTitle(*trig);
1045 effs->SetMinimum(0);
1046 effs->SetMaximum(1.4);
1047 effs->SetMarkerSize(2.5);
1048 effs->Draw("hist text");
1056 c->SaveAs("triggers.png");