8#include <cdc/calibration/T0CalibrationAlgorithm.h>
9#include <calibration/CalibrationAlgorithm.h>
10#include <cdc/dbobjects/CDCTimeZeros.h>
11#include <cdc/dataobjects/WireID.h>
12#include <cdc/geometry/CDCGeometryPar.h>
16#include <TGraphErrors.h>
20#include <TStopwatch.h>
22#include <framework/database/DBObjPtr.h>
23#include <framework/database/IntervalOfValidity.h>
24#include <framework/logging/Logger.h>
34 " -------------------------- T0 Calibration Algorithm -------------------------\n"
41 B2INFO(
"Creating histograms");
50 auto tree = getObjectPtr<TTree>(
"tree");
51 tree->SetBranchAddress(
"lay", &lay);
52 tree->SetBranchAddress(
"IWire", &IWire);
53 tree->SetBranchAddress(
"x_u", &x);
54 tree->SetBranchAddress(
"t", &t_mea);
55 tree->SetBranchAddress(
"t_fit", &t_fit);
56 tree->SetBranchAddress(
"ndf", &ndf);
57 tree->SetBranchAddress(
"Pval", &Pval);
61 std::vector<TString> list_vars = {
"lay",
"IWire",
"x_u",
"t",
"t_fit",
"Pval",
"ndf"};
62 tree->SetBranchStatus(
"*", 0);
64 for (TString brname : list_vars) {
65 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);
82 for (
int il = 0; il < 56; il++) {
84 for (
int ic = 0; ic < nW; ++ic) {
85 m_h1[il][ic] =
new TH1F(Form(
"hdT_L%d_W%d", il, ic), Form(
"L%d_cell%d", il, ic), 30, -15, 15);
90 for (
int ib = 0; ib < 300; ++ib) {
91 m_hT0b[ib] =
new TH1F(Form(
"hdT_b%d", ib), Form(
"boardID_%d", ib), 100, -20, 20);
95 const Long64_t nEntries = tree->GetEntries();
96 B2INFO(
"Number of entries: " << nEntries);
99 for (Long64_t i = 0; i < nEntries; ++i) {
101 double xmax = halfCSize[lay] - 0.1;
102 if ((fabs(x) <
m_xmin) || (fabs(x) > xmax)
107 m_h1[lay][IWire]->Fill(t_mea - t_fit);
110 m_hT0b[boardID]->Fill(t_mea - t_fit);
113 B2INFO(
"Finish making histogram for all channels");
114 B2INFO(
"Time to fill histograms: " << timer.RealTime() <<
"s");
120 B2INFO(
"Start calibration");
123 gErrorIgnoreLevel = 3001;
131 B2INFO(
"ExpRun used for DB Geometry : " << exprun.first <<
" " << exprun.second);
142 B2INFO(
"Creating CDCGeometryPar object");
145 auto hEvtT0 = getObjectPtr<TH1F>(
"hEventT0");
146 double dEventT0 = hEvtT0->GetMean();
148 TH1F* hm_All =
new TH1F(
"hm_All",
"mean of #DeltaT distribution for all chanels", 500, -10, 10);
149 TH1F* hs_All =
new TH1F(
"hs_All",
"#sigma of #DeltaT distribution for all chanels", 100, 0, 10);
152 TF1* g1 =
new TF1(
"g1",
"gaus", -100, 100);
153 g1->SetParLimits(1, -20, 20);
154 vector<double> b, db, Sb, dSb;
155 vector<double> c[56];
156 vector<double> dc[56];
157 vector<double> s[56];
158 vector<double> ds[56];
160 B2INFO(
"Gaus fitting for whole channel");
164 m_hTotal->Fit(
"g1",
"Q",
"", mean - 15, mean + 15);
165 g1->GetParameters(par);
167 B2INFO(
"Gaus fitting for each board");
168 for (
int ib = 1; ib < 300; ++ib) {
174 if (
m_hT0b[ib]->Integral(1,
m_hT0b[ib]->GetNbinsX()) < 50) {
178 mean =
m_hT0b[ib]->GetMean();
179 m_hT0b[ib]->SetDirectory(0);
180 m_hT0b[ib]->Fit(
"g1",
"Q",
"", mean - 15, mean + 15);
181 g1->GetParameters(par);
184 Sb.push_back(par[1]);
185 dSb.push_back(g1->GetParError(1));
186 dtb[ib] = par[1] + dEventT0;
187 err_dtb[ib] = g1->GetParError(1);
189 B2INFO(
"Gaus fitting for each cell");
190 for (
int ilay = 0; ilay < 56; ++ilay) {
192 for (
unsigned int iwire = 0; iwire < nW; ++iwire) {
194 const int n =
m_h1[ilay][iwire]->Integral(1,
m_h1[ilay][iwire]->GetNbinsX()) ;
195 B2DEBUG(21,
"layer " << ilay <<
" wire " << iwire <<
" entries " << n);
203 err_dt[ilay][iwire] = 50;
continue;
206 mean =
m_h1[ilay][iwire]->GetMean();
207 m_h1[ilay][iwire]->SetDirectory(0);
208 m_h1[ilay][iwire]->Fit(
"g1",
"Q",
"", mean - 15, mean + 15);
209 g1->GetParameters(par);
210 c[ilay].push_back(iwire);
211 dc[ilay].push_back(0.0);
212 s[ilay].push_back(par[1]);
213 ds[ilay].push_back(g1->GetParError(1));
215 dt[ilay][iwire] = par[1] + dEventT0;
216 err_dt[ilay][iwire] = g1->GetParError(1);
217 hm_All->Fill(par[1]);
218 hs_All->Fill(par[2]);
223 B2INFO(
"Storing histograms");
224 auto hNDF = getObjectPtr<TH1F>(
"hNDF");
225 auto hPval = getObjectPtr<TH1F>(
"hPval");
226 TFile* fout =
new TFile(
m_histName.c_str(),
"RECREATE");
228 TGraphErrors* gr[56];
229 TDirectory* top = gDirectory;
232 if (hNDF && hPval && hEvtT0) {
240 TDirectory* subDir[56];
241 for (
int ilay = 0; ilay < 56; ++ilay) {
242 subDir[ilay] = top ->mkdir(Form(
"lay_%d", ilay));
245 for (
unsigned int iwire = 0; iwire < nW; ++iwire) {
246 m_h1[ilay][iwire]->Write();
251 TDirectory* corrT0 = top->mkdir(
"DeltaT0");
255 TGraphErrors* grb =
new TGraphErrors(b.size(), &b.at(0), &Sb.at(0), &db.at(0), &dSb.at(0));
256 grb->SetMarkerColor(2);
257 grb->SetMarkerSize(1.0);
258 grb->SetTitle(
"#DeltaT0;BoardID;#DeltaT0[ns]");
260 grb->SetMinimum(-10);
261 grb->SetName(
"Board");
264 for (
int sl = 0; sl < 56; ++sl) {
265 if (c[sl].size() < 2)
continue;
266 gr[sl] =
new TGraphErrors(c[sl].size(), &c[sl].at(0), &s[sl].at(0), &dc[sl].at(0), &ds[sl].at(0));
267 gr[sl]->SetMarkerColor(2);
268 gr[sl]->SetMarkerSize(1.0);
269 gr[sl]->SetTitle(Form(
"Layer_%d;IWire;#LT t_{mea}-t_{fit} #GT [ns]", sl));
270 gr[sl]->SetMaximum(10);
271 gr[sl]->SetMinimum(-10);
272 gr[sl]->SetName(Form(
"lay%d", sl));
277 B2INFO(
"Writing constants...");
278 int Ngood_T0 =
write();
279 B2INFO(
"Checking conversion conditons...");
280 double rms_dT = hm_All->GetRMS();
281 double n_below = hm_All->Integral(0, hm_All->GetXaxis()->FindBin(
m_offsetMeanDt - 5 * rms_dT));
282 double n_upper = hm_All->Integral(hm_All->GetXaxis()->FindBin(
m_offsetMeanDt + 5 * rms_dT), hm_All->GetXaxis()->GetNbins() - 1);
283 int N_remaining = n_below + n_upper - (14112 - Ngood_T0);
284 if (N_remaining < 0) N_remaining = 0;
285 B2INFO(
"+ Number of channel which are properly calibrated : " << Ngood_T0);
286 B2INFO(
"+ Number of channel which still need to be calibrated are: " << N_remaining <<
"(requirement <" <<
m_maxBadChannel <<
")");
287 B2INFO(
"+ Median of Delta_T - offset:" << hm_All->GetMean() -
m_offsetMeanDt <<
"(requirement: <" <<
m_maxMeanDt <<
")");
288 B2INFO(
"+ RMS of Delta_T dist. :" << rms_dT <<
" (Requirement: <" <<
m_maxRMSDt <<
")");
293 B2INFO(
"T0 Calibration Finished:");
296 B2INFO(
"Need more iteration ...");
307 for (
int ib = 0; ib < 300; ++ib) {
308 hT0B[ib] =
new TH1F(Form(
"hT0B%d", ib), Form(
"boardID_%d", ib), 9000, 0, 9000);
311 TH1F* hT0_all =
new TH1F(
"hT0_all",
"T0 distribution", 9000, 0, 9000);
314 for (
int ilay = 0; ilay < 56; ++ilay) {
316 for (
unsigned int iwire = 0; iwire < nW; ++iwire) {
317 WireID wireid(ilay, iwire);
318 T0 = cdcgeo.
getT0(wireid);
326 if (fabs(T0_average -
m_commonT0) > 50) {B2WARNING(
"Large difference between common T0 (" <<
m_commonT0 <<
") and aveage value" << T0_average);}
329 for (
int i = 0; i < 300; ++i) {T0B[i] =
m_commonT0;}
332 for (
int ib = 1; ib < 300; ++ib) {
334 if (fabs(T0B[ib] - T0_average) > 25) {
335 B2WARNING(
"T0 of Board " << ib <<
" (= " << T0B[ib] <<
") is too different with common T0: " << T0_average <<
336 "\n It will be replaced by common T0");
337 T0B[ib] = T0_average;
341 if (abs(
err_dtb[ib]) < 2 && abs(
dtb[ib]) < 20) {
349 for (
int ilay = 0; ilay < 56; ++ilay) {
351 for (
unsigned int iwire = 0; iwire < nW; ++iwire) {
352 WireID wireid(ilay, iwire);
356 T0 = cdcgeo.
getT0(wireid);
357 if (fabs(T0 - T0B[bID]) > 25 || fabs(T0 - T0_average) > 25) {
358 B2WARNING(
"T0 of wireID L-W: " << ilay <<
"--" << iwire <<
" (= " << T0
359 <<
") is too different with common T0 of board: " << bID <<
" = " << T0B[bID] <<
360 "\n It will be replaced by common T0 of the Board");
366 if (abs(
err_dt[ilay][iwire]) < 2 && abs(
dt[ilay][iwire]) < 20) {
367 dT =
dt[ilay][iwire];
372 tz->
setT0(wireid, T0 - dT);
384 double mean1 = h1->GetMean();
385 h1->GetXaxis()->SetRangeUser(mean1 - 50, 8000);
386 double mean2 = h1->GetMean();
387 while (fabs(mean1 - mean2) > 0.5) {
389 h1->GetXaxis()->SetRangeUser(mean1 - 50, 8000);
390 mean2 = h1->GetMean();
Database object for timing offset (t0).
void outputToFile(std::string fileName) const
Output the contents in text file format.
void setT0(unsigned short iCLayer, unsigned short iWire, float t0)
Set t0 in the list.
The Class for CDC Geometry Parameters.
unsigned short getBoardID(const WireID &wID) const
Returns frontend board id. corresponding to the wire id.
float getT0(const WireID &wireID) const
Returns t0 parameter of the specified sense wire.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
double m_xmin
minimum drift length
std::string m_outputT0FileName
output t0 file name for text file
double m_offsetMeanDt
offset dT distribution caused by event timing extraction;
double dtb[300]
dt of each board
std::string m_histName
root file name
TH1F * m_hTotal
1D histogram of delta T whole channel
double dt[56][385]
dt of each channel
bool m_storeHisto
store histo or not
double m_ndfmin
minimum ndf required
void createHisto()
create histo for each channel
double getMeanT0(TH1F *h1)
calculate mean of the T0 distribution
T0CalibrationAlgorithm()
Constructor.
double m_maxRMSDt
RMS of dT distribution of all channels.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
double m_commonT0
A common T0 of all channels.
double m_maxMeanDt
Mean of dT distribution of all channels;.
TH1F * m_h1[56][385]
1D histogram for each channel
double err_dt[56][385]
error of dt of each channel
bool m_textOutput
output text file if true
EResult calibrate() override
Run algo on data.
TH1F * m_hT0b[300]
1D histogram for each board
double m_Pvalmin
minimum pvalue required
int write()
write outut or store db
double err_dtb[300]
error of dt of board
double m_maxBadChannel
Number of channels which has DeltaT0 larger then 0.5.
Base class for calibration algorithms.
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_Iterate
Needs iteration =1 in Python.
Class to identify a wire inside the CDC.
Abstract base class for different kinds of events.