Skip to content

Commit 7e038b4

Browse files
f3schf3sch
committed
ITS: modify staggering macro
Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>
1 parent 6cecb46 commit 7e038b4

File tree

1 file changed

+65
-226
lines changed

1 file changed

+65
-226
lines changed

‎Detectors/ITSMFT/ITS/macros/test/CheckStaggering.C‎

Lines changed: 65 additions & 226 deletions
Original file line numberDiff line numberDiff line change
@@ -13,9 +13,11 @@
1313
#include <cmath>
1414
#include <vector>
1515
#include <filesystem>
16-
#include <vector>
16+
#include <algorithm>
17+
#include <random>
1718

1819
#include <TGeoManager.h>
20+
#include <TRandom.h>
1921
#include <TFile.h>
2022
#include <TTree.h>
2123
#include <TF1.h>
@@ -30,16 +32,6 @@
3032
#include <TLine.h>
3133
#include <TLorentzVector.h>
3234
#include <TPaveText.h>
33-
#include <RooRealVar.h>
34-
#include <RooDataHist.h>
35-
#include <RooGaussian.h>
36-
#include <RooAddPdf.h>
37-
#include "RooVoigtian.h"
38-
#include <RooChebychev.h>
39-
#include <RooPlot.h>
40-
#include <RooFitResult.h>
41-
#include <RooHist.h>
42-
#include <RooFit.h>
4335

4436
#include "DetectorsBase/Propagator.h"
4537
#include "DataFormatsITSMFT/CompCluster.h"
@@ -56,7 +48,6 @@
5648
#include "DetectorsVertexing/SVertexHypothesis.h"
5749
#include "DCAFitter/DCAFitterN.h"
5850
#endif
59-
using namespace RooFit;
6051

6152
constexpr const char* tracFile = "o2trac_its.root";
6253
constexpr const char* clsFile = "o2clus_its.root";
@@ -69,8 +60,6 @@ struct PairCandidate {
6960
};
7061

7162
std::vector<std::filesystem::path> findDirs(const std::string&);
72-
void fitK0(TH1D*, int);
73-
void fitPhiMeson(TH1D*, int);
7463

7564
void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
7665
{
@@ -82,6 +71,9 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
8271
return;
8372
}
8473
if (max > 0 && (int)dirs.size() > max) {
74+
std::random_device rd;
75+
std::mt19937 g(rd());
76+
std::shuffle(dirs.begin(), dirs.end(), g);
8577
dirs.resize(max);
8678
LOGP(info, "restricting to {} dirs", max);
8779
}
@@ -111,11 +103,12 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
111103
auto hTrkTSE = new TH1D("hTrkTSE", "assigned time errors; tE (BC)", 198, -0.5, 198 - 0.5);
112104

113105
// K0 && Phi-Meson
114-
const float mMinITSPt{0.3};
106+
const float mMinITSPt{0.15};
115107
//
116108
const int phiMinNCls{7};
117-
const float phiMaxDCAR2MVTX{0.1}; // max distance to mean vtx
118-
auto hPhiMeson = new TH1D("hPhiMeson", "#phi meson;mass (GeV/c^{2})", 100, 0.96, 1.1);
109+
const float phiMaxDCAR2MVTX{0.05}; // max distance to mean vtx
110+
auto hPhiMeson = new TH1D("hPhiMeson", "#phi meson;mass (GeV/c^{2})", 200, 0.96, 1.1);
111+
auto hPhiMesonBkg = new TH1D("hPhiMesonBkg", "#phi meson background;mass (GeV/c^{2})", 200, 0.96, 1.1);
119112

120113
const int mK0MinNCls{7};
121114
const float mK0minCosPAXYMeanVertex = 0.98;
@@ -133,7 +126,7 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
133126
k0Ft.setMinRelChi2Change(0.9);
134127
k0Ft.setMaxChi2(5);
135128
k0Ft.setUseAbsDCA(true);
136-
auto hK0 = new TH1D("hK0Meson", "K0;mass (GeV/c^{2})", 100, 0.4, 0.6);
129+
auto hK0 = new TH1D("hK0", "K0;mass (GeV/c^{2})", 100, 0.4, 0.6);
137130
o2::vertexing::SVertexHypothesis k0Hyp;
138131
const float k0Par[] = {0., 20, 0., 5.0, 0.0, 1.09004e-03, 2.62291e-04, 8.93179e-03, 2.83121};
139132
k0Hyp.set(o2::track::PID::K0, o2::track::PID::Pion, o2::track::PID::Pion, k0Par, bz);
@@ -160,18 +153,34 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
160153
std::vector<o2::its::TrackITS> trkArr, *trkArrPtr{&trkArr};
161154
std::vector<o2::its::Vertex> vtxArr, *vtxArrPtr{&vtxArr};
162155
std::array<std::vector<o2::itsmft::CompClusterExt>*, 7> clsArr{nullptr};
163-
for (const auto& d : dirs) {
164-
LOGP(info, "Entering {}", d.c_str());
156+
for (size_t iDir{0}; iDir < dirs.size(); ++iDir) {
157+
int progress = static_cast<int>((iDir + 1) * 100 / dirs.size());
158+
printf("\rProgress: [");
159+
int barWidth = 50;
160+
int pos = barWidth * progress / 100;
161+
for (int j = 0; j < barWidth; ++j) {
162+
if (j < pos)
163+
printf("=");
164+
else if (j == pos)
165+
printf(">");
166+
else
167+
printf(" ");
168+
}
169+
printf("] %d%%", progress);
170+
fflush(stdout);
171+
172+
const auto& d = dirs[iDir];
173+
LOGP(debug, "\nEntering {}", d.c_str());
165174
auto fTrks = TFile::Open((d / tracFile).c_str());
166175
auto fCls = TFile::Open((d / clsFile).c_str());
167176
if (!fTrks || !fCls || fTrks->IsZombie() || fCls->IsZombie()) {
168-
LOGP(error, "Could not open files");
177+
LOGP(error, "\nCould not open files");
169178
continue;
170179
}
171180
auto tTrks = fTrks->Get<TTree>("o2sim");
172181
auto tCls = fCls->Get<TTree>("o2sim");
173182
if (!tTrks || !tCls) {
174-
LOGP(error, "Could not retrieve trees");
183+
LOGP(error, "\nCould not retrieve trees");
175184
continue;
176185
}
177186

@@ -189,7 +198,7 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
189198
for (int i{0}; i < 7; ++i) {
190199
ncls += clsArr[i]->size();
191200
}
192-
LOGP(info, "\tTF:{:03} {} trks {} vtx {} cls", iTF, trkArr.size(), vtxArr.size(), ncls);
201+
LOGP(debug, "\n\tTF:{:03} {} trks {} vtx {} cls", iTF, trkArr.size(), vtxArr.size(), ncls);
193202

194203
// for each TF built pool of positive and negaitve tracks
195204
std::vector<const o2::its::TrackITS*> posPool, negPool;
@@ -248,8 +257,19 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
248257
p1.SetXYZM(pP[0], pP[1], pP[2], posPhi.getPID().getMass());
249258
p2.SetXYZM(pN[0], pN[1], pN[2], negPhi.getPID().getMass());
250259
TLorentzVector mother = p1 + p2;
251-
double mass = mother.M();
252-
hPhiMeson->Fill(mass);
260+
hPhiMeson->Fill(mother.M());
261+
// rotate on daughter track to estimate background
262+
for (int i{0}; i < 10; ++i) {
263+
double theta = gRandom->Uniform(165.f, 195.f) * TMath::DegToRad();
264+
double pxN = pN[0] * cos(theta) - pN[1] * sin(theta);
265+
double pyN = pN[0] * sin(theta) + pN[1] * cos(theta);
266+
double pxP = pP[0] * cos(-theta) - pP[1] * sin(-theta);
267+
double pyP = pP[0] * sin(-theta) + pP[1] * cos(-theta);
268+
p1.SetXYZM(pxP, pyP, pP[2], posPhi.getPID().getMass());
269+
p2.SetXYZM(pxN, pyN, pN[2], negPhi.getPID().getMass());
270+
mother = p1 + p2;
271+
hPhiMesonBkg->Fill(mother.M());
272+
}
253273
}
254274
}
255275
}
@@ -359,8 +379,26 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
359379
return lastLine;
360380
};
361381

362-
fitK0(hK0, runNumber);
363-
fitPhiMeson(hPhiMeson, runNumber);
382+
{
383+
TFile out(Form("check_%d.root", runNumber), "RECREATE");
384+
out.WriteTObject(hNTrkCls);
385+
for (int i{0}; i < 5; ++i) {
386+
out.WriteTObject(hTrkTS[i]);
387+
}
388+
out.WriteTObject(hTrkTSE);
389+
out.WriteTObject(hPhiMeson);
390+
out.WriteTObject(hPhiMesonBkg);
391+
out.WriteTObject(hK0);
392+
out.WriteTObject(hVtxXY);
393+
out.WriteTObject(hVtxZ);
394+
out.WriteTObject(hVtxNCont);
395+
out.WriteTObject(hVtxZNCont);
396+
out.WriteTObject(hVtxTS);
397+
out.WriteTObject(hVtxCls);
398+
}
399+
400+
// fitK0(hK0, runNumber);
401+
// fitPhiMeson(hPhiMeson, runNumber);
364402
{
365403
auto c = new TCanvas();
366404
hNTrkCls->Draw();
@@ -485,202 +523,3 @@ std::vector<std::filesystem::path> findDirs(const std::string& roots)
485523
}
486524
return result;
487525
}
488-
489-
void fitK0(TH1D* h, int runNumber)
490-
{
491-
RooRealVar mass("mass", "K^{0} mass (GeV/c^{2})", 0.46, 0.54);
492-
RooDataHist data("data", "data", RooArgList(mass), Import(*h));
493-
RooRealVar mean("mean", "mean", 0.4976, 0.490, 0.503);
494-
RooRealVar sigma1("sigma1", "core sigma", 0.003, 0.0005, 0.02);
495-
RooRealVar sigma2("sigma2", "tail sigma", 0.006, 0.001, 0.03);
496-
RooGaussian gauss1("gauss1", "gauss1", mass, mean, sigma1);
497-
RooGaussian gauss2("gauss2", "gauss2", mass, mean, sigma2);
498-
RooRealVar frac2("frac2", "fraction in second gaussian", 0.15, 0.0, 1.0);
499-
RooAddPdf signal("signal", "signal PDF", RooArgList(gauss1, gauss2), RooArgList(frac2));
500-
RooRealVar c1("c1", "cheby coeff1", 0.0, -5.0, 5.0);
501-
RooChebychev bkg("bkg", "background", mass, RooArgList(c1));
502-
double totalIntegral = h->Integral();
503-
RooRealVar nsig("nsig", "signal yield", totalIntegral * 0.4, 0., totalIntegral * 10.0);
504-
RooRealVar nbkg("nbkg", "background yield", totalIntegral * 0.6, 0., totalIntegral * 10.0);
505-
RooAddPdf model("model", "sig + bkg", RooArgList(signal, bkg), RooArgList(nsig, nbkg));
506-
RooFitResult* fitRes = model.fitTo(data,
507-
Extended(true),
508-
Save(true),
509-
Strategy(1),
510-
Hesse(true),
511-
Minos(false),
512-
PrintLevel(-1));
513-
514-
TCanvas* c = new TCanvas("c", "K0 fit", 800, 800);
515-
c->cd();
516-
TPad* pad1 = new TPad("pad1", "pad1", 0, 0.30, 1, 1.0);
517-
TPad* pad2 = new TPad("pad2", "pad2", 0, 0.0, 1, 0.30);
518-
pad1->SetBottomMargin(0.0);
519-
pad2->SetTopMargin(0.0);
520-
pad2->SetBottomMargin(0.30);
521-
pad1->Draw();
522-
pad2->Draw();
523-
524-
pad1->cd();
525-
RooPlot* frame = mass.frame(Title("K^{0}_{S} mass fit"));
526-
data.plotOn(frame, Name("data"));
527-
model.plotOn(frame, Name("model"));
528-
// draw background component
529-
model.plotOn(frame, Components("bkg"), LineStyle(kDashed), LineColor(kRed), Name("bkg"));
530-
// draw total signal component
531-
model.plotOn(frame, Components("signal"), LineStyle(kDashed), LineColor(kBlue), Name("signal"));
532-
frame->GetYaxis()->SetTitle("Events / bin");
533-
frame->GetXaxis()->SetLabelSize(0.05);
534-
frame->GetXaxis()->SetTitleSize(0.06);
535-
frame->GetYaxis()->SetTitleSize(0.05);
536-
frame->Draw();
537-
538-
TLegend leg(0.60, 0.65, 0.88, 0.88);
539-
leg.SetBorderSize(0);
540-
leg.AddEntry((TObject*)0, Form("mean = %.3f #pm %.3f GeV/c^{2}", mean.getVal(), mean.getError()), "");
541-
leg.AddEntry((TObject*)0, Form("sigma_{1} = %.3f #pm %.3f GeV/c^{2}", sigma1.getVal(), sigma1.getError()), "");
542-
leg.AddEntry((TObject*)0, Form("sigma_{2} = %.3f #pm %.3f GeV/c^{2}", sigma2.getVal(), sigma2.getError()), "");
543-
leg.AddEntry((TObject*)0, Form("signal yield = %.0f #pm %.0f", nsig.getVal(), nsig.getError()), "");
544-
leg.Draw();
545-
546-
int npar = fitRes->floatParsFinal().getSize();
547-
double chi2ndf = frame->chiSquare(npar);
548-
TPaveText* pt = new TPaveText(0.15, 0.75, 0.45, 0.88, "NDC");
549-
pt->SetFillStyle(0);
550-
pt->SetBorderSize(0);
551-
pt->AddText(Form("#chi^{2}/ndf = %.3g", chi2ndf));
552-
pt->Draw();
553-
554-
pad2->cd();
555-
RooPlot* frame_pull = mass.frame(Title("Pull (data - fit) / uncertainty"));
556-
RooHist* hpull = frame->pullHist();
557-
frame_pull->addPlotable(hpull, "P");
558-
frame_pull->GetXaxis()->SetTitle("mass (GeV/c^{2})");
559-
frame_pull->GetXaxis()->SetTitleSize(0.12);
560-
frame_pull->GetXaxis()->SetLabelSize(0.10);
561-
frame_pull->GetYaxis()->SetLabelSize(0.10);
562-
frame_pull->GetYaxis()->SetTitleSize(0.12);
563-
frame_pull->Draw();
564-
frame_pull->GetYaxis()->SetRangeUser(-5, 5);
565-
TLine l(frame_pull->GetXaxis()->GetXmin(), 0, frame_pull->GetXaxis()->GetXmax(), 0);
566-
l.SetLineStyle(kDashed);
567-
l.Draw();
568-
c->SaveAs(Form("k0_%d.pdf", runNumber));
569-
}
570-
571-
void fitPhiMeson(TH1D* h, int runNumber)
572-
{
573-
double hxmin = h->GetXaxis()->GetXmin();
574-
double hxmax = h->GetXaxis()->GetXmax();
575-
int nbins = h->GetNbinsX();
576-
double sig_low = 1.015;
577-
double sig_high = 1.025;
578-
TH1D* h_sb = (TH1D*)h->Clone(Form("%s_sidebands", h->GetName()));
579-
for (int ib = 1; ib <= nbins; ++ib) {
580-
double x = h_sb->GetBinCenter(ib);
581-
if (x > sig_low && x < sig_high) {
582-
h_sb->SetBinContent(ib, 0.0);
583-
h_sb->SetBinError(ib, 0.0);
584-
}
585-
}
586-
RooRealVar mass("mass", "m_{K^{+}K^{-}} (GeV/c^{2})", hxmin, hxmax);
587-
RooDataHist data("data", "data", RooArgList(mass), Import(*h));
588-
RooDataHist data_sb("data_sb", "sideband data", RooArgList(mass), Import(*h_sb));
589-
RooRealVar mean("mean", "phi mass", 1.019, 1.015, 1.025);
590-
RooRealVar gamma("gamma", "natural width (GeV)", 0.00426);
591-
gamma.setConstant(true);
592-
RooRealVar sigma("sigma", "detector resolution (GeV)", 0.001, 0.0002, 0.01);
593-
RooVoigtian signal("signal", "Voigtian(signal)", mass, mean, gamma, sigma);
594-
RooRealVar c1("c1", "cheby coeff1", 0.0, -5.0, 5.0);
595-
RooRealVar c2("c2", "cheby coeff2", 0.0, -5.0, 5.0);
596-
RooChebychev bkg("bkg", "background", mass, RooArgList(c1, c2));
597-
RooFitResult* fbkg = bkg.fitTo(data_sb,
598-
Range(hmin, hxmax),
599-
Save(true),
600-
Verbose(false),
601-
PrintLevel(-1));
602-
603-
if (!fbkg) {
604-
LOGP(error, "Warning: background sideband fit failed to return RooFitResult.");
605-
}
606-
double totalIntegral = h->Integral();
607-
RooRealVar nsig("nsig", "signal yield", totalIntegral * 0.05, 0., totalIntegral * 10.0); // small initial guess
608-
RooRealVar nbkg("nbkg", "background yield", totalIntegral * 0.95, 0., totalIntegral * 10.0);
609-
RooAddPdf model("model", "sig + bkg", RooArgList(signal, bkg), RooArgList(nsig, nbkg));
610-
RooFitResult* fitRes = model.fitTo(data,
611-
Extended(true),
612-
Save(true),
613-
Strategy(1),
614-
Hesse(true),
615-
Minos(false),
616-
Verbose(true),
617-
PrintLevel(-1));
618-
619-
if (!fitRes) {
620-
LOGP(error, "fitPhiMeson: fitTo returned null RooFitResult.");
621-
}
622-
TCanvas* c = new TCanvas("c_phi", "Phi mass fit", 800, 800);
623-
c->cd();
624-
TPad* pad1 = new TPad("pad1", "pad1", 0, 0.30, 1, 1.0);
625-
TPad* pad2 = new TPad("pad2", "pad2", 0, 0.0, 1, 0.30);
626-
pad1->SetBottomMargin(0.0);
627-
pad2->SetTopMargin(0.0);
628-
pad2->SetBottomMargin(0.30);
629-
pad1->Draw();
630-
pad2->Draw();
631-
pad1->cd();
632-
633-
RooPlot* frame = mass.frame(Title("Phi (#phi) meson fit"));
634-
frame->GetXaxis()->SetTitle("m_{K^{+}K^{-}} (GeV/c^{2})");
635-
frame->GetYaxis()->SetTitle("Events / bin");
636-
data.plotOn(frame, Binning(nbins), Name("data"));
637-
model.plotOn(frame, Name("model"));
638-
model.plotOn(frame, Components("bkg"), LineStyle(kDashed), LineColor(kRed), Name("bkg"));
639-
model.plotOn(frame, Components("signal"), LineStyle(kSolid), LineColor(kBlue), Name("signal"));
640-
frame->Draw();
641-
642-
TLegend leg(0.55, 0.60, 0.88, 0.88);
643-
leg.SetBorderSize(0);
644-
leg.SetFillStyle(0);
645-
leg.AddEntry("data", "Data", "lep");
646-
leg.AddEntry("model", "Total fit", "l");
647-
leg.AddEntry("signal", "Signal (Voigtian)", "l");
648-
leg.AddEntry("bkg", "Background (Chebychev)", "l");
649-
leg.Draw();
650-
651-
TPaveText* pv = new TPaveText(0.15, 0.60, 0.45, 0.88, "NDC");
652-
pv->SetFillStyle(0);
653-
pv->SetBorderSize(0);
654-
655-
pv->AddText(Form("m = %.6g #pm %.6g GeV", mean.getVal(), mean.getError()));
656-
pv->AddText(Form("#sigma_{res} = %.6g #pm %.6g GeV", sigma.getVal(), sigma.getError()));
657-
pv->AddText(Form("width (fixed) = %.6g GeV", gamma.getVal()));
658-
pv->AddText(Form("Signal yield = %.1f #pm %.1f", nsig.getVal(), nsig.getError()));
659-
pv->AddText(Form("Background yield = %.1f #pm %.1f", nbkg.getVal(), nbkg.getError()));
660-
pv->Draw();
661-
662-
// chi^2 / ndf (approx via RooPlot::chiSquare)
663-
int nfloat = fitRes ? fitRes->floatParsFinal().getSize() : 0;
664-
double chi2ndf = frame->chiSquare(nfloat);
665-
TPaveText* pv2 = new TPaveText(0.15, 0.52, 0.45, 0.58, "NDC");
666-
pv2->SetFillStyle(0);
667-
pv2->SetBorderSize(0);
668-
pv2->AddText(Form("#chi^{2}/ndf = %.3g", chi2ndf));
669-
pv2->Draw();
670-
pad2->cd();
671-
RooPlot* frame_pull = mass.frame(Title("Pull (data-fit)/#sigma"));
672-
RooHist* hpull = frame->pullHist(); // pulls computed from frame
673-
frame_pull->addPlotable(hpull, "P");
674-
frame_pull->GetXaxis()->SetTitle("m_{K^{+}K^{-}} (GeV/c^{2})");
675-
frame_pull->GetXaxis()->SetTitleSize(0.12);
676-
frame_pull->GetXaxis()->SetLabelSize(0.10);
677-
frame_pull->GetYaxis()->SetLabelSize(0.10);
678-
frame_pull->GetYaxis()->SetTitleSize(0.12);
679-
frame_pull->Draw();
680-
double xmin = frame_pull->GetXaxis()->GetXmin();
681-
double xmax = frame_pull->GetXaxis()->GetXmax();
682-
TLine l(xmin, 0.0, xmax, 0.0);
683-
l.SetLineStyle(kDashed);
684-
l.Draw();
685-
c->SaveAs(Form("phi_%d.pdf", runNumber));
686-
}

0 commit comments

Comments
 (0)