BmnRoot
Loading...
Searching...
No Matches
BmnTriggerEfficiencyRun7.cxx
Go to the documentation of this file.
2
4{
5
6 if (!fSpectraFile.IsNull()) {
7
8 // Efficiency for trigger that contains either FD @and@ BD conditions (like, BT + BD1 + FD2) ...
9 if (fSpectraFile.Contains("FD") && fSpectraFile.Contains("BD")) {
10 // Triggers with three or more conditions are not supported now ...
11 if (fEffFiles.size() > 2)
12 return;
13
15
16 // Loop over add. files passed ...
17 for (auto effFile : fEffFiles)
18 lAnal->SetTriggerEffData(effFile.ReplaceAll("eff_", "").ReplaceAll(".root", ""),
19 effFile); // FIXME !!!
20
21 lAnal->ReadData();
22
23 TClonesArray* filledData = lAnal->GetFilledArray();
24
25 TClonesArray* out = new TClonesArray("TriggerEfficiency");
26
27 for (Int_t iData1 = 0; iData1 < filledData->GetEntriesFast(); iData1++) {
28 TriggerEfficiency* eff1 = (TriggerEfficiency*)filledData->UncheckedAt(iData1);
29
30 TString trig1 = eff1->trigger();
31 Int_t minMul1 = eff1->multilplicity().first;
32 Int_t maxMul1 = eff1->multilplicity().second;
33
34 // Starting from iData1 ...
35 for (Int_t iData2 = iData1; iData2 < filledData->GetEntriesFast(); iData2++) {
36 TriggerEfficiency* eff2 = (TriggerEfficiency*)filledData->UncheckedAt(iData2);
37
38 TString trig2 = eff2->trigger();
39 Int_t minMul2 = eff2->multilplicity().first;
40 Int_t maxMul2 = eff2->multilplicity().second;
41
42 if (minMul1 != minMul2 || maxMul1 != maxMul2 || trig1 == trig2)
43 continue;
44
45 Double_t eff = eff1->efficiency() * eff2->efficiency();
46
47 Double_t relErr1 = eff1->efficiencyError() / eff1->efficiency();
48 Double_t relErr2 = eff2->efficiencyError() / eff2->efficiency();
49
50 Double_t relErr = TMath::Sqrt(relErr1 * relErr1 + relErr2 * relErr2);
51 Double_t absErr = eff * relErr;
52
53 new ((*out)[out->GetEntriesFast()])
54 TriggerEfficiency(fTrigger, make_pair(minMul1, maxMul1), make_pair(eff, absErr));
55 }
56 }
57
58 // Clear multiplicity map from set values before ...
59 fMultMap.clear();
60
61 // Setting number of multiplicity bins (possibly, could change before by user)
62 nMultBins = out->GetEntriesFast();
63
64 vector<Double_t> triggEffs;
65 vector<Double_t> triggEffErrs;
66
67 for (Int_t iEntry = 0; iEntry < out->GetEntriesFast(); iEntry++) {
68 TriggerEfficiency* eff = (TriggerEfficiency*)out->UncheckedAt(iEntry);
69
70 fMultMap[iEntry] = make_pair(eff->multilplicity().first, eff->multilplicity().second);
71 triggEffs.push_back(eff->efficiency());
72 triggEffErrs.push_back(eff->efficiencyError());
73 }
74
75 delete lAnal;
76
77 TFile* eff = new TFile(Form("eff_%s.root", fTrigger.Data()), "recreate");
78 TH1F* effGr =
79 new TH1F("Eff", Form("Efficiency of %s", fTrigger.Data()), triggEffs.size(), 0., triggEffs.size());
80
81 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++) {
82 effGr->SetBinContent(iMultBin + 1, (triggEffs.at(iMultBin) > 1.) ? 1. : triggEffs.at(iMultBin));
83 effGr->SetBinError(iMultBin + 1, triggEffErrs.at(iMultBin));
84
85 Int_t min = fMultMap.find(iMultBin)->second.first;
86 Int_t max = fMultMap.find(iMultBin)->second.second;
87
88 effGr->GetXaxis()->SetBinLabel(iMultBin + 1, Form("[%d, %d)", min, max));
89 }
90
91 effGr->SetMarkerStyle(20);
92 effGr->SetMarkerColor(kSpring - 6);
93 effGr->SetLineColor(kSpring - 6);
94 effGr->SetLineWidth(1);
95 effGr->GetXaxis()->SetTitle("nTracks in V_{p}");
96 effGr->GetYaxis()->SetTitle("");
97 effGr->GetXaxis()->SetNdivisions(5);
98 effGr->GetYaxis()->SetNdivisions(5);
99 effGr->GetYaxis()->SetLabelOffset(0.01);
100 effGr->GetXaxis()->SetLabelSize(0.06);
101 effGr->GetYaxis()->SetLabelSize(0.04);
102 effGr->GetXaxis()->SetTitleSize(0.05);
103 effGr->GetYaxis()->SetTitleSize(0.05);
104 effGr->GetXaxis()->SetTitleOffset(0.85);
105 effGr->GetYaxis()->SetTitleOffset(1.5);
106 effGr->GetYaxis()->CenterTitle();
107 effGr->Draw("PE1X0");
108
109 eff->Write();
110 delete eff;
111
112 return;
113 }
114
115 // Efficiency for trigger that contains either FD @or@ BD conditions
116
117 TFile* f = new TFile(fSpectraFile.Data());
118 if (!f->IsOpen())
119 return;
120
121 TString targ = "";
122 for (auto target : fTarget)
123 targ += target;
124
125 auto it = fMultMap.begin();
126
127 const Int_t nTrigCond = 2;
128
129 map<Int_t, pair<TH1F*, TH1F*>> histoMap; // mult. bin --> pair of (denom., nom.)
130
131 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++) {
132 auto itShifted = next(it, iMultBin);
133
134 vector<TH1F*> vSet;
135 for (Int_t iCond = 0; iCond < nTrigCond; iCond++) {
136 TString hName = "Target: " + targ
137 + TString::Format(", Trigger: %s (%s)", fTrigger.Data(),
138 (iCond == TriggConditions::severe) ? "nominator" : "denominator")
139 + TString::Format(", minMult (>=) %d, maxMult (<) %d", itShifted->second.first,
140 itShifted->second.second);
141
142 TH1F* h = (TH1F*)f->Get(hName.Data());
143 vSet.push_back(h);
144 }
145
146 histoMap[iMultBin] = make_pair(vSet.at(0), vSet.at(1)); // denom., nom. (!!!)
147 }
148
149 Double_t mean, sigma = -1.;
150 pair<Double_t, Double_t> T = make_pair(-1., -1.);
151 pair<Double_t, Double_t> B = make_pair(-1., -1.);
152
153 TFile* out = new TFile("outH.root", "recreate");
154
155 vector<TH1F*> noms;
156 vector<TH1F*> denoms;
157
158 vector<Double_t> nomsSignal;
159 vector<Double_t> denomsSignal;
160
161 vector<Double_t> errsRelTot;
162
163 for (auto histo : histoMap) {
164 TH1F* hNom = histo.second.second;
165 TH1F* hDenom = histo.second.first;
166
167 TF1* f1 = nullptr;
168 Double_t dS = 0., dB = 0., errTotRelNom = 0., errTotRelDenom = 0.;
169
170 // Fitting spectra for certain mult. bin ...
171 fitSpectrum(hNom, mean, sigma, T, B);
172 nomsSignal.push_back(T.second - B.second);
173
174 delete hNom->GetListOfFunctions()->FindObject("f");
175 delete hNom->GetListOfFunctions()->FindObject("f2");
176 f1 = (TF1*)hNom->GetListOfFunctions()->FindObject("f1");
177 f1->SetLineColor(kRed);
178 dS = TMath::Sqrt(T.second - B.second);
179 dB = deltaB(f1, fitSpectrum((TH1F*)hNom->Clone()));
180
181 errTotRelNom = TMath::Sqrt((dS / (T.second - B.second)) * (dS / (T.second - B.second))
182 + (dB / B.second) * (dB / B.second));
183 hNom->Write();
184
185 //
186 fitSpectrum(hDenom, mean, sigma, T, B);
187 denomsSignal.push_back(T.second - B.second);
188
189 delete hDenom->GetListOfFunctions()->FindObject("f");
190 delete hDenom->GetListOfFunctions()->FindObject("f2");
191 f1 = (TF1*)hDenom->GetListOfFunctions()->FindObject("f1");
192 f1->SetLineColor(kRed);
193 dS = TMath::Sqrt(T.second - B.second);
194 dB = deltaB(f1, fitSpectrum((TH1F*)hDenom->Clone()));
195
196 errTotRelDenom = TMath::Sqrt((dS / (T.second - B.second)) * (dS / (T.second - B.second))
197 + (dB / B.second) * (dB / B.second));
198 hDenom->Write();
199
200 noms.push_back(hNom);
201 denoms.push_back(hDenom);
202
203 errsRelTot.push_back(TMath::Sqrt(errTotRelNom * errTotRelNom + errTotRelDenom * errTotRelDenom));
204 }
205
206 TCanvas* c = new TCanvas("c", "c", 1200, 600);
207 c->Divide(nMultBins, 2);
208
209 gStyle->SetOptStat(0);
210 gStyle->SetTitleFontSize(0.1);
211 gStyle->SetTitleY(1.0);
212
213 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++) {
214 TH1F* n = noms.at(iMultBin);
215 c->cd(iMultBin + 1);
216 n->Draw("PE1X0");
217 n->SetTitle("");
218
219 TLatex* t1 = new TLatex();
220 t1->SetNDC(); // we want to draw this object in the Normalized system [0,1]
221 t1->SetTextSize(.045);
222 TString name = n->GetName();
223 name.ReplaceAll("Target: PbSnCuAl,", "")
224 .ReplaceAll("Trigger: ", "")
225 .ReplaceAll("nominator", "nom.")
226 .ReplaceAll("minMult", "nVpTr:")
227 .ReplaceAll("maxMult", "");
228 t1->DrawLatex(.15, .15, name.Data());
229 t1->DrawLatex(.1, .85, Form("S (bin cont. in sig. reg. - bckgrnd) = %d", (Int_t)nomsSignal.at(iMultBin)));
230
231 TH1F* d = denoms.at(iMultBin);
232 c->cd(iMultBin + 1 + nMultBins);
233 d->Draw("PE1X0");
234 d->SetTitle("");
235
236 t1 = new TLatex();
237 t1->SetNDC(); // we want to draw this object in the Normalized system [0,1]
238 t1->SetTextSize(.045);
239 name = d->GetName();
240 name.ReplaceAll("Target: PbSnCuAl,", "")
241 .ReplaceAll("Trigger: ", "")
242 .ReplaceAll("denominator", "den.")
243 .ReplaceAll("minMult", "nVpTr:")
244 .ReplaceAll("maxMult", "");
245 t1->DrawLatex(.15, .15, name.Data());
246 t1->SetTextSize(.04);
247 t1->DrawLatex(.15, .85,
248 Form("S (bin cont. in sig. reg. - bckgrnd) = %d", (Int_t)denomsSignal.at(iMultBin)));
249 }
250
251 c->SaveAs(Form("triggEff_%s_multBins.pdf", fTrigger.Data()));
252
253 vector<Double_t> triggEffs;
254 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++)
255 triggEffs.push_back(nomsSignal.at(iMultBin) / denomsSignal.at(iMultBin));
256
257 TFile* eff = new TFile(Form("eff_%s.root", fTrigger.Data()), "recreate");
258 TH1F* effGr =
259 new TH1F("Eff", Form("Efficiency of %s", fTrigger.Data()), triggEffs.size(), 0., triggEffs.size());
260
261 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++) {
262 effGr->SetBinContent(iMultBin + 1, (triggEffs.at(iMultBin) > 1.) ? 1. : triggEffs.at(iMultBin));
263 effGr->SetBinError(iMultBin + 1, triggEffs.at(iMultBin) * errsRelTot.at(iMultBin));
264
265 Int_t min = fMultMap.find(iMultBin)->second.first;
266 Int_t max = fMultMap.find(iMultBin)->second.second;
267
268 effGr->GetXaxis()->SetBinLabel(iMultBin + 1, Form("[%d, %d)", min, max));
269 }
270
271 effGr->SetMarkerStyle(20);
272 effGr->SetMarkerColor(kSpring - 6);
273 effGr->SetLineColor(kSpring - 6);
274 effGr->SetLineWidth(1);
275 effGr->GetXaxis()->SetTitle("nTracks in V_{p}");
276 effGr->GetYaxis()->SetTitle("");
277 effGr->GetXaxis()->SetNdivisions(5);
278 effGr->GetYaxis()->SetNdivisions(5);
279 effGr->GetYaxis()->SetLabelOffset(0.01);
280 effGr->GetXaxis()->SetLabelSize(0.06);
281 effGr->GetYaxis()->SetLabelSize(0.04);
282 effGr->GetXaxis()->SetTitleSize(0.05);
283 effGr->GetYaxis()->SetTitleSize(0.05);
284 effGr->GetXaxis()->SetTitleOffset(0.85);
285 effGr->GetYaxis()->SetTitleOffset(1.5);
286 effGr->GetYaxis()->CenterTitle();
287 effGr->Draw("PE1X0");
288
289 eff->Write();
290 delete eff;
291
292 delete f;
293 delete out;
294
295 return;
296 }
297
298 // Getting geom. cut map (for all targets) ...
299 map<TString, BmnParticlePairCut*> tCutMap = GetTargetCutsMap();
300
301 BmnParticlePairCut* cuts = nullptr;
302
303 TString targ = "";
304 for (auto target : fTarget) {
305 targ += target;
306 cuts = tCutMap.find(target)->second;
307 }
308
309 // Preparing histos ...
310 const Int_t nTrigCond = 2; // Nominator and denominator (always two) ...
311 hSpectra = new TH1F**[nTrigCond];
312 hMult = new TH1F*[nTrigCond];
313
314 auto it = fMultMap.begin();
315
316 for (Int_t iCond = 0; iCond < nTrigCond; iCond++) {
317 hSpectra[iCond] = new TH1F*[nMultBins];
318
319 TString hName = "Target: " + targ
320 + TString::Format(", Trigger: %s (%s)", fTrigger.Data(),
321 (iCond == TriggConditions::severe) ? "nominator" : "denominator");
322
323 hMult[iCond] = new TH1F(hName.Data(), hName.Data(), 100, 0., 100.);
324
325 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++) {
326 auto itShifted = next(it, iMultBin);
327 hName = "Target: " + targ
328 + TString::Format(", Trigger: %s (%s)", fTrigger.Data(),
329 (iCond == TriggConditions::severe) ? "nominator" : "denominator")
330 + TString::Format(", minMult (>=) %d, maxMult (<) %d", itShifted->second.first,
331 itShifted->second.second);
332 hSpectra[iCond][iMultBin] = new TH1F(hName.Data(), hName.Data(), 75, xLow, xUp);
333 }
334 }
335
336 // Processing dst's ...
337 if (cuts) {
338 vector<TString> triggsSet;
339
340 // Setting datasets to be processed ...
341 if (fTrigger == "BT+FD3" || fTrigger == "BT+FD2") {
342 triggsSet.push_back("BT+BD2");
343 triggsSet.push_back("BT+BD3");
344 } else if (fTrigger == "BT+BD3" || fTrigger == "BT+BD2" || fTrigger == "BT+BD1") {
345 triggsSet.push_back("BT+FD2");
346 triggsSet.push_back("BT+FD3");
347 } else if (fTrigger == "BT+BD1+FD2")
348 return;
349
350 else
351 Fatal("BmnTriggerEfficiencyRun7::triggerEfficiency", "Trigger condition not supported!!!");
352
353 for (auto dst : createFilelist(triggsSet)) {
354 // Getting denom. histos (base) ...
355 isAddTriggerCondition = kFALSE;
356 ReadFile(dst, cuts);
357
358 // Getting nomin. histos (severe) ...
359 isAddTriggerCondition = kTRUE;
360 ReadFile(dst, cuts);
361 }
362 }
363
364 // Writing filled histos to out file ...
365 TFile* f = new TFile(Form("out_%s_start%d_finish%d.root", fTrigger.Data(), fStartRun, fFinishRun), "recreate");
366
367 for (Int_t iCond = 0; iCond < nTrigCond; iCond++)
368 for (Int_t iMultBin = 0; iMultBin < nMultBins; iMultBin++)
369 hSpectra[iCond][iMultBin]->Write();
370
371 delete f;
372}
373
374vector<TString> BmnTriggerEfficiencyRun7::createFilelist(vector<TString> triggsSet)
375{
376 vector<TString> list;
377
378 // Access to trigger information ...
380
381 for (Int_t iFile = fStartRun; iFile < fFinishRun; iFile++) {
382 UniRun* pCurrentRun = UniRun::GetRun(fPeriod, iFile);
383
384 // Check presense of the current run in DB ...
385 if (!pCurrentRun)
386 continue;
387
388 if (pCurrentRun->GetBeamParticle() != fBeam)
389 continue;
390
391 // Trigger check ...
392 Bool_t isNeededTrigger = kFALSE;
393
394 for (auto it : triggsSet) {
395 if (trInfo->GetTrigger(iFile) != it)
396 continue;
397 else {
398 isNeededTrigger = kTRUE;
399 break;
400 }
401 }
402
403 if (!isNeededTrigger)
404 continue;
405
406 TString targ = pCurrentRun->GetTargetParticle();
407
408 // Skipping runs w/o target ...
409 if (targ == "")
410 continue;
411
412 Bool_t isNeededTarget = kFALSE;
413
414 for (auto it : fTarget) {
415 if (it != targ)
416 continue;
417 else {
418 isNeededTarget = kTRUE;
419 break;
420 }
421 }
422
423 if (!isNeededTarget)
424 continue;
425
426 for (TString tmp : fInFiles) {
427 if (!tmp.Contains(Form("%d", iFile)))
428 continue;
429
430 list.push_back(tmp);
431 break;
432 }
433 }
434
435 delete trInfo;
436
437 return list;
438}
439
440void BmnTriggerEfficiencyRun7::ReadFile(TString f, BmnParticlePairCut* cut0)
441{
442 TChain* out = new TChain("bmndata");
443 out->Add(f.Data());
444
445 if (!out->GetFile()) {
446 delete out;
447 return;
448 }
449
450 TClonesArray* particlePair = nullptr;
451 out->SetBranchAddress("ParticlePair", &particlePair);
452 if (!particlePair)
453 return;
454
455 DstEventHeaderExtended* eHeaderExt = nullptr;
456 out->SetBranchAddress("DstEventHeaderExtended.", &eHeaderExt);
457 if (!eHeaderExt)
458 return;
459
460 for (Int_t iEv = 0; iEv < out->GetEntries(); iEv++) {
461 out->GetEntry(iEv);
462
463 Int_t nVpTracks = eHeaderExt->nTracks();
464 Double_t z = eHeaderExt->Vp().Z();
465
466 // -6. < VpZ < +6. [cm]
467 const Double_t zCor = 6.;
468
469 if (TMath::Abs(z) > zCor || nVpTracks < 2)
470 continue;
471
472 // Getting trigger conditions and track multiplicity ...
473 if ((fTrigger.Contains("BT+FD3") || fTrigger.Contains("BT+FD2")) && isAddTriggerCondition) {
474 // Nominator ...
475 nFdCounts = (fTrigger.Contains("BT+FD3")) ? 3 : 2;
476 if (eHeaderExt->fd() <= nFdCounts)
477 continue;
478 } else if ((fTrigger.Contains("BT+BD3") || fTrigger.Contains("BT+BD2") || fTrigger.Contains("BT+BD1"))
479 && isAddTriggerCondition)
480 {
481 // Nominator ...
482 nBdCounts = (fTrigger.Contains("BT+BD3")) ? 3 : (fTrigger.Contains("BT+BD2")) ? 2 : 1;
483 if (eHeaderExt->bd() <= nBdCounts)
484 continue;
485 }
486
487 for (Int_t iPair = 0; iPair < particlePair->GetEntriesFast(); iPair++) {
488 BmnParticlePair* pair = (BmnParticlePair*)particlePair->UncheckedAt(iPair);
489
490 // V0 ...
491 Double_t v0z = pair->GetV0Z();
492
493 // Skipping so far reconstructed V0's ...
494 if (v0z < 0. || v0z > 50.)
495 continue;
496
497 // Geometry topology ...
498 Double_t DCA12 = pair->GetDCA12();
499 Double_t DCA1 = pair->GetDCA1();
500 Double_t DCA2 = pair->GetDCA2();
501 Double_t DCA0 = pair->GetDCA0();
502
503 Double_t PATH = pair->GetPath();
504
505 // Track length (GEM part of glob. track) ...
506 Int_t protonNGemHits = pair->GetNHitsPart1("GEM");
507 Int_t pionNGemHits = pair->GetNHitsPart2("GEM");
508
509 BmnParticlePairCut* cut = cut0;
510
511 Double_t dca0 = cut->dca0();
512 Double_t dca1 = cut->dca1();
513 Double_t dca2 = cut->dca2();
514 Double_t dca12 = cut->dca12();
515
516 Double_t path = cut->path();
517
518 Int_t nPos = cut->nHitsGemPos();
519 Int_t nNeg = cut->nHitsGemNeg();
520
521 if (DCA0 > dca0 || DCA12 > dca12) // <=
522 continue;
523
524 if (DCA1 < dca1 || DCA2 < dca2 || PATH < path) // >=
525 continue;
526
527 if (protonNGemHits <= nPos || pionNGemHits <= nNeg) // >
528 continue;
529
530 if (FinMultBin(nVpTracks) != -1)
531 hSpectra[isAddTriggerCondition ? (TriggConditions::severe) : (TriggConditions::base)]
532 [FinMultBin(nVpTracks)]
533 ->Fill(pair->GetInvMass());
534 }
535 }
536
537 cout << "File " << f << " processed# " << endl;
538 delete out;
539}
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
float f
Definition P4_F32vec4.h:21
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:31
map< TString, BmnParticlePairCut * > GetTargetCutsMap()
void SetTriggerEffData(TString trigger, TString data)
TClonesArray * GetFilledArray()
vector< TString > createFilelist()
Double_t deltaB(TF1 *f, TFitResultPtr fitRes)
vector< TString > fInFiles
TFitResultPtr fitSpectrum(TH1F *)
vector< TString > fTarget
Int_t GetNHitsPart1(TString det="")
Int_t GetNHitsPart2(TString det="")
Double_t GetInvMass()
pair< Int_t, Int_t > multilplicity()
static UniRun * GetRun(int period_number, int run_number)
get run from the database
Definition UniRun.cxx:177
TString GetBeamParticle()
get beam particle of the current run
Definition UniRun.h:113
TString GetTargetParticle()
get target particle of the current run
Definition UniRun.h:115