1 #include <cdc/calibration/T0Correction.h>
2 #include <cdc/dbobjects/CDCTimeZeros.h>
3 #include <cdc/geometry/CDCGeometryPar.h>
4 #include <cdc/dataobjects/WireID.h>
8 #include <TGraphErrors.h>
13 #include <framework/database/IntervalOfValidity.h>
14 #include <framework/database/DBImportObjPtr.h>
15 #include <framework/logging/Logger.h>
21 T0Correction::T0Correction():
22 m_firstExperiment(0), m_firstRun(0),
23 m_lastExperiment(-1), m_lastRun(-1)
40 B2INFO(
"CreateHisto");
51 TChain* tree =
new TChain(
"tree");
53 tree->SetBranchAddress(
"lay", &lay);
54 tree->SetBranchAddress(
"IWire", &IWire);
55 tree->SetBranchAddress(
"x_u", &x);
56 tree->SetBranchAddress(
"t", &t_mea);
57 tree->SetBranchAddress(
"t_fit", &t_fit);
58 tree->SetBranchAddress(
"weight", &w);
59 tree->SetBranchAddress(
"ndf", &ndf);
60 tree->SetBranchAddress(
"Pval", &Pval);
63 std::vector<TString> list_vars = {
"lay",
"IWire",
"x_u",
"t",
"t_fit",
"weight",
"Pval",
"ndf"};
64 tree->SetBranchStatus(
"*", 0);
66 for (TString brname : list_vars) {
67 tree->SetBranchStatus(brname, 1);
73 for (
int i = 0; i < 56; ++i) {
76 halfCSize[i] = M_PI * R / nW;
79 m_hTotal =
new TH1F(
"hTotal",
"hTotal", 30, -15, 15);
81 for (
int il = 0; il < 56; il++) {
83 for (
int ic = 0; ic < nW; ++ic) {
84 m_h1[il][ic] =
new TH1F(Form(
"h1%d@%d", il, ic), Form(
"L%d_cell%d", il, ic), 30, -15, 15);
88 for (
int ib = 0; ib < 300; ++ib) {
89 m_hT0b[ib] =
new TH1F(Form(
"hT0b%d", ib), Form(
"boardID_%d", ib), 100, -20, 20);
92 const int nEntries = tree->GetEntries();
93 B2INFO(
"Number of entries: " << nEntries);
94 for (
int i = 0; i < nEntries; ++i) {
96 double xmax = halfCSize[lay] - 0.1;
97 if ((fabs(x) <
m_xmin) || (fabs(x) > xmax)
102 m_h1[lay][IWire]->Fill(t_mea - t_fit);
106 B2INFO(
"Finish making histogram for all channels");
111 B2INFO(
"Start calibration");
114 gErrorIgnoreLevel = 3001;
119 TF1* g1 =
new TF1(
"g1",
"gaus", -100, 100);
120 vector<double> b, db, Sb, dSb;
121 vector<double> c[56];
122 vector<double> dc[56];
123 vector<double> s[56];
124 vector<double> ds[56];
126 B2INFO(
"Gaus fitting for whole channel");
130 m_hTotal->Fit(
"g1",
"Q",
"", mean - 15, mean + 15);
131 g1->GetParameters(par);
133 B2INFO(
"Gaus fitting for each board");
134 for (
int ib = 1; ib < 300; ++ib) {
135 if (
m_hT0b[ib]->GetEntries() < 10)
continue;
136 mean =
m_hT0b[ib]->GetMean();
137 m_hT0b[ib]->SetDirectory(0);
138 m_hT0b[ib]->Fit(
"g1",
"Q",
"", mean - 15, mean + 15);
139 g1->GetParameters(par);
142 Sb.push_back(par[1]);
143 dSb.push_back(g1->GetParError(1));
145 err_dtb[ib] = g1->GetParError(1);
147 B2INFO(
"Gaus fitting for each cell");
148 for (
int ilay = 0; ilay < 56; ++ilay) {
149 for (
unsigned int iwire = 0; iwire < cdcgeo.
nWiresInLayer(ilay); ++iwire) {
150 const int n =
m_h1[ilay][iwire]->GetEntries();
151 B2DEBUG(99,
"layer " << ilay <<
" wire " << iwire <<
" entries " << n);
152 if (n < 10)
continue;
153 mean =
m_h1[ilay][iwire]->GetMean();
154 m_h1[ilay][iwire]->SetDirectory(0);
155 m_h1[ilay][iwire]->Fit(
"g1",
"Q",
"", mean - 15, mean + 15);
156 g1->GetParameters(par);
157 c[ilay].push_back(iwire);
158 dc[ilay].push_back(0.0);
159 s[ilay].push_back(par[1]);
160 ds[ilay].push_back(g1->GetParError(1));
162 dt[ilay][iwire] = par[1];
163 err_dt[ilay][iwire] = g1->GetParError(1);
169 B2INFO(
"Store histo");
170 TFile* fout =
new TFile(
"Correct_T0.root",
"RECREATE");
172 TGraphErrors* gr[56];
173 TDirectory* top = gDirectory;
175 TDirectory* subDir[56];
176 for (
int ilay = 0; ilay < 56; ++ilay) {
177 subDir[ilay] = top ->mkdir(Form(
"lay_%d", ilay));
179 for (
unsigned int iwire = 0; iwire < cdcgeo.
nWiresInLayer(ilay); ++iwire) {
180 m_h1[ilay][iwire]->Write();
184 TDirectory* corrT0 = top->mkdir(
"DeltaT0");
188 TGraphErrors* grb =
new TGraphErrors(b.size(), &b.at(0), &Sb.at(0), &db.at(0), &dSb.at(0));
189 grb->SetMarkerColor(2);
190 grb->SetMarkerSize(1.0);
191 grb->SetTitle(
"#DeltaT0;BoardID;#DeltaT0[ns]");
193 grb->SetMinimum(-10);
194 grb->SetName(
"Board");
197 for (
int sl = 0; sl < 56; ++sl) {
198 if (c[sl].size() < 2)
continue;
199 gr[sl] =
new TGraphErrors(c[sl].size(), &c[sl].at(0), &s[sl].at(0), &dc[sl].at(0), &ds[sl].at(0));
200 gr[sl]->SetMarkerColor(2);
201 gr[sl]->SetMarkerSize(1.0);
202 gr[sl]->SetTitle(Form(
"Layer_%d;IWire;#LT t_{mea}-t_{fit} #GT [ns]", sl));
203 gr[sl]->SetMaximum(10);
204 gr[sl]->SetMinimum(-10);
205 gr[sl]->SetName(Form(
"lay%d", sl));
210 B2INFO(
"Write constants");
223 for (
int ib = 0; ib < 300; ++ib) {
224 T0B[ib] =
new TH1F(Form(
"T0B%d", ib), Form(
"boardID_%d", ib), 9000, 0, 9000);
227 for (
int ilay = 0; ilay < 56; ++ilay) {
228 for (
unsigned int iwire = 0; iwire < cdcgeo.
nWiresInLayer(ilay); ++iwire) {
229 WireID wireid(ilay, iwire);
230 T0 = cdcgeo.
getT0(wireid);
236 for (
int ilay = 0; ilay < 56; ++ilay) {
237 for (
unsigned int iwire = 0; iwire < cdcgeo.
nWiresInLayer(ilay); ++iwire) {
238 WireID wireid(ilay, iwire);
240 if (abs(
err_dt[ilay][iwire]) > 2 || abs(
dt[ilay][iwire]) > 1.e3) {
242 T0 = T0B[bID]->GetMean();
243 dt[ilay][iwire] =
dtb[bID];
245 T0 = cdcgeo.
getT0(wireid);
247 ofs << ilay <<
"\t" << iwire <<
"\t" << T0 -
dt[ilay][iwire] << std::endl;
249 tz->
setT0(wireid, T0 -
dt[ilay][iwire]);
258 B2RESULT(
"T0 was calibrated and imported to database.");