8#include <cdc/calibration/TimeWalkCalibrationAlgorithm.h>
9#include <cdc/dbobjects/CDCTimeWalks.h>
10#include <cdc/geometry/CDCGeometryPar.h>
11#include <cdc/dataobjects/WireID.h>
15#include <TDirectory.h>
18#include <TStopwatch.h>
20#include <framework/database/DBObjPtr.h>
21#include <framework/logging/Logger.h>
30 " -------------------------- Time Walk Calibration Algorithm -------------------------\n"
36 B2INFO(
"Creating and filling histograms");
41 for (
int i = 0; i < 56; ++i) {
44 halfCSize[i] = M_PI *
R / nW;
48 for (
int i = 0; i < 300; ++i) {
49 m_h2[i] =
new TH2D(Form(
"board_%d", i), Form(
"board_%d", i), 50, 0., 500, 60, -30, 30);
65 tree->SetBranchAddress(
"lay", &lay);
66 tree->SetBranchAddress(
"IWire", &IWire);
67 tree->SetBranchAddress(
"x_u", &x);
68 tree->SetBranchAddress(
"t", &t_mea);
69 tree->SetBranchAddress(
"t_fit", &t_fit);
70 tree->SetBranchAddress(
"ndf", &ndf);
71 tree->SetBranchAddress(
"Pval", &Pval);
72 tree->SetBranchAddress(
"adc", &adc);
75 std::vector<TString> list_vars = {
"lay",
"IWire",
"x_u",
"t",
"t_fit",
"Pval",
"ndf",
"adc"};
76 tree->SetBranchStatus(
"*", 0);
78 for (TString brname : list_vars) {
79 tree->SetBranchStatus(brname, 1);
82 const Long64_t nEntries = tree->GetEntries();
83 B2INFO(
"Number of entries: " << nEntries);
86 for (Long64_t i = 0; i < nEntries; ++i) {
88 const double xmax = halfCSize[lay] - 0.12;
89 if ((fabs(x) <
m_xmin) || (fabs(x) > xmax)
96 B2INFO(
"Time to fill histograms: " << time.RealTime() <<
"s");
102 B2INFO(
"Start calibration");
106 B2INFO(
"ExpRun used for DB Geometry : " << exprun.first <<
" " << exprun.second);
108 B2INFO(
"Creating CDCGeometryPar object");
116 fold =
new TF1(
"fold",
"[0]/sqrt(x)", 0, 500);
118 fold =
new TF1(
"fold",
"[0]*exp(-1*[1]*x)", 0, 500);
120 if (fold ==
nullptr) {
121 B2FATAL(
"Old fitting function is not defined.");
125 B2INFO(
"time walk formular: ");
126 [](TF1 * f) {
auto title = f->GetTitle(); B2INFO(
"Expression : " << title);}(fold);
129 for (
int ib = 1; ib < 300; ++ib) {
131 B2DEBUG(21,
"Board ID:" << ib);
132 m_h2[ib]->SetDirectory(0);
135 if (
m_h2[ib]->GetEntries() < 500) {
139 m_h2[ib]->FitSlicesY(0, 1, -1, 10);
140 TString name =
m_h2[ib]->GetName();
141 TString hm_name = name +
"_1";
142 m_h1[ib] = (TH1D*)gDirectory->Get(hm_name);
144 m_h1[ib]->SetDirectory(0);
145 if (
m_h1[ib]->GetEntries() < 5) {
147 B2WARNING(
"Low statistic, number of points after slice fit: " <<
m_h1[ib]->GetEntries());
153 fold->SetParameter(p,
m_tw_old[ib][p]);
158 TF1*
func =
new TF1(
"func",
"[0]+[1]/sqrt(x)", 0, 500);
159 func->SetParameters(-4, 28);
160 m_h1[ib]->Fit(
"func",
"MQ",
"", 20, 150);
165 B2FATAL(
"Mode " <<
m_twParamMode <<
" is not available, please check again");
169 TF1* f1 =
m_h1[ib]->GetFunction(
"func");
170 if (!f1) {
m_flag[ib] = 0;
continue;}
173 m_tw_new[ib][i - 1] = f1->GetParameter(i);
177 B2DEBUG(21,
"Prob of fitting:" << f1->GetProb());
178 B2DEBUG(21,
"Fitting Param 0-1:" << f1->GetParameter(0) <<
" - " << f1->GetParameter(1));
195 B2INFO(
"Storing histogram");
196 B2DEBUG(21,
"Store 1D histogram");
197 TFile* fhist =
new TFile(
m_histName.c_str(),
"RECREATE");
202 if (hNDF && hPval && hEvtT0) {
208 TDirectory* old = gDirectory;
209 TDirectory* h1D = old->mkdir(
"h1D");
210 TDirectory* h2D = old->mkdir(
"h2D");
212 for (
int ib = 1; ib < 300; ++ib) {
214 if (
m_h1[ib]->GetEntries() < 5)
continue;
215 m_h1[ib]->SetMinimum(-5);
216 m_h1[ib]->SetMaximum(15);
220 B2DEBUG(21,
"Store 2D histogram");
222 for (
int ib = 1; ib < 300; ++ib) {
223 if (
m_h2[ib] ==
nullptr)
continue;
224 if (
m_h2[ib]->GetEntries() < 5)
continue;
230 B2INFO(
"Hitograms were stored");
235 TH1F* hDtw =
new TH1F(
"hDtw",
"", 100, -1, 1);
236 for (
int ib = 0; ib < 300; ++ib) {
238 if (std::isnan(dtw) == 0) {
243 B2INFO(hDtw->GetRMS());
245 if (hDtw->GetRMS() < 0.02) {
255 B2INFO(
"Save to the local DB");
260 for (
int ib = 0; ib < 300; ++ib) {
265 const int num =
static_cast<int>(
m_tw_old[ib].size());
266 for (
int i = 0; i < num; ++i) {
276 B2RESULT(
"Time-walk coeff. table has been written to " <<
m_outputTWFileName.c_str());
280 B2RESULT(
"Failure to calibrate time-walk for " << nfailure <<
" board");
286 B2INFO(
"Prepare calibration");
292 B2FATAL(
"Mode " <<
m_twParamMode <<
" haven't implemented yet.");
296 const int nEntries = dbTw->getEntries();
297 for (
int ib = 0; ib < nEntries; ++ib) {
298 m_tw_old[ib] = dbTw->getTimeWalkParams(ib);
305 B2INFO(
"Add constant term into T0 database");
309 for (
int ilay = 0; ilay < 56; ++ilay) {
311 for (
unsigned int iwire = 0; iwire < nW; ++iwire) {
312 WireID wireid(ilay, iwire);
314 T0 = cdcgeo.
getT0(wireid);
318 tz->
setT0(wireid, T0);
331 int max = h1->GetMaximumBin();
332 double maxX = h1->GetBinCenter(max);
333 double maxY = h1->GetBinContent(max);
334 B2DEBUG(21,
"Max: id - x - y : " << max <<
" " << maxX <<
" " << maxY);
338 h1->Fit(
"pol0",
"MQE",
"", maxX + 125, 500);
339 if (h1->GetFunction(
"pol0")) {
340 p0 = h1->GetFunction(
"pol0")->GetParameter(0);
344 TH1D* hshift =
new TH1D(
"hshift",
"shift", h1->GetNbinsX(), 0, 500);
345 hshift->SetDirectory(0);
346 for (
int i = 0; i <= h1->GetNbinsX(); ++i) {
347 hshift->SetBinContent(i, h1->GetBinContent(i) - p0);
349 hshift->Fit(
"expo",
"MQE",
"", 0, 400);
350 double p1 = maxY + p0;
352 if (hshift->GetFunction(
"expo")) {
353 p1 = exp(hshift->GetFunction(
"expo")->GetParameter(0));
354 p2 = hshift->GetFunction(
"expo")->GetParameter(1);
358 TF1*
func =
new TF1(
"func",
"[0]+ [1]*exp(-1*[2]*x)", 0, 500);
359 func->SetParameters(p0, p1, -1 * p2);
360 func->SetParLimits(0, -5, 5);
361 func->SetParLimits(1, -5, 500);
362 func->SetParLimits(2, -0.01, 0.1);
363 h1->Fit(
"func",
"MQE",
"", 10, 400);
Database object for time-walk.
void outputToFile(std::string fileName) const
Output the contents in text file format.
void setTwParamMode(unsigned short mode)
Set tw parameterization mode mode=0: tw (in ns) = p0/sqrt(FADCsum); mode=1: tw( in ns) = p0*exp(-p1*F...
void setTimeWalkParams(unsigned short boardID, const std::vector< float > ¶ms)
Set the time-walk coefficients in the list.
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 value cut of drift length.
void prepare()
prepare calibration.
std::string m_outputT0FileName
t0 file name after calibration.
TH2D * m_h2[300]
2D histogram of residual vs ADC for each board.
unsigned short m_twParamMode
=0 for P0/Sqrt(ADC); =1 for P0*Exp(-P1*ADC).
void storeHist()
Store histograms.
std::string m_histName
root file name
double m_constTerm[300]
const term in fitting, it will be added to T0 instead tw
bool m_storeHisto
Store all Histogram or not.
TimeWalkCalibrationAlgorithm()
Constructor.
double m_minNdf
minimum number of degree of freedom required for track.
unsigned short m_flag[300]
flag for fit status
void createHisto()
Create histo for calibrate.
double m_minPval
minimum number of Prob(chi2) of fitted track.
void updateT0()
update constant term to t0 database.
std::vector< float > m_tw_old[300]
tw list before calibration.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
std::string m_outputTWFileName
tw file name after calibration.
std::vector< float > m_tw_new[300]
tw list after calibration.
void write()
save calibration.
unsigned short m_nTwParams
No.
void fitToExponentialFunc(TH1D *h1)
fit tw histogram
bool m_textOutput
output text file if true
EResult calibrate() override
Run algo on data.
EResult checkConvergence()
check convergence.
TH1D * m_h1[300]
Mean of residual as function of ADC of each board.
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.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
Class for accessing objects in the database.
Class to identify a wire inside the CDC.
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Abstract base class for different kinds of events.