8#include <cdc/calibration/XTCalibrationAlgorithm.h>
9#include <cdc/calibration/XTFunction.h>
10#include <cdc/geometry/CDCGeometryPar.h>
11#include <cdc/dbobjects/CDCXtRelations.h>
14#include <TStopwatch.h>
23#include <framework/database/DBObjPtr.h>
29typedef std::array<float, 3> array3;
33 " -------------------------- XT Calibration Algorithm -------------------------\n"
40 B2INFO(
"create and fill histo");
42 for (
int i = 0; i < 56; ++i) {
43 for (
int lr = 0; lr < 2; ++lr) {
46 m_hProf[i][lr][al][th] =
new TProfile(Form(
"m_hProf%d_%d_%d_%d", i, lr, al, th),
47 Form(
"(L=%d)-(lr=%d)-(#alpha=%3.0f)-(#theta=%3.0f); Drift time (ns);Drift Length (cm)",
49 m_hist2d[i][lr][al][th] =
new TH2F(Form(
"h%d_%d_%d_%d", i, lr, al, th),
50 Form(
"(L=%d)-(lr=%d)-(#alpha=%3.0f)-(#theta=%3.0f); Drift time (ns);Drift Length (cm)",
53 m_hist2dDraw[i][al][th] =
new TH2F(Form(
"h_draw%d_%d_%d", i, al, th),
54 Form(
"(L=%d)-(#alpha=%3.0f)-(#theta=%3.0f); Drift time (ns);Drift Length (cm)",
68 Float_t Pval, alpha, theta;
71 tree->SetBranchAddress(
"lay", &lay);
72 tree->SetBranchAddress(
"t", &dt);
73 tree->SetBranchAddress(
"x_u", &dx);
74 tree->SetBranchAddress(
"alpha", &alpha);
75 tree->SetBranchAddress(
"theta", &theta);
76 tree->SetBranchAddress(
"Pval", &Pval);
77 tree->SetBranchAddress(
"ndf", &ndf);
80 std::vector<TString> list_vars = {
"lay",
"t",
"x_u",
"alpha",
"theta",
"Pval",
"ndf"};
81 tree->SetBranchStatus(
"*", 0);
83 for (TString brname : list_vars) {
84 tree->SetBranchStatus(brname, 1);
92 const Long64_t nEntries = tree->GetEntries();
93 B2INFO(
"Number of entries " << nEntries);
94 for (Long64_t i = 0; i < nEntries; ++i) {
111 int lr = dx > 0 ? c_Right : c_Left;
113 m_hProf[lay][lr][al][th]->Fill(dt, abs(dx));
114 m_hist2d[lay][lr][al][th]->Fill(dt, abs(dx));
116 m_hProf[lay][0][al][th]->Fill(dt, abs(dx));
117 m_hist2d[lay][0][al][th]->Fill(dt, abs(dx));
118 m_hProf[lay][1][al][th]->Fill(dt, abs(dx));
119 m_hist2d[lay][1][al][th]->Fill(dt, abs(dx));
124 B2INFO(
"Time to fill histograms: " << time.RealTime() <<
"s");
131 gPrintViaErrorHandler =
true;
132 gErrorIgnoreLevel = 3001;
133 B2INFO(
"Start calibration");
137 B2INFO(
"ExpRun used for DB Geometry : " << exprun.first <<
" " << exprun.second);
139 B2INFO(
"Creating CDCGeometryPar object");
145 B2INFO(
"Start Fitting");
146 std::unique_ptr<XTFunction> xt;
147 for (
int l = 0; l < 56; ++l) {
148 for (
int lr = 0; lr < 2; ++lr) {
158 m_hist2d[l][lr][al][th]->FitSlicesY(0, 0, -1, 5);
159 m_hist2d_1[l][lr][al][th] = (TH1F*)gDirectory->Get(Form(
"h%d_%d_%d_%d_1", l, lr, al, th));
162 B2WARNING(
"Error, not found results of slices fit");
165 m_hist2d_1[l][lr][al][th]->Fit(
"pol1",
"Q",
"", 30, 60);
166 fpol1 = (TF1*)
m_hProf[l][lr][al][th]->GetFunction(
"pol1");
169 for (
int n = 0; n <
m_hProf[l][lr][al][th]->GetNbinsX(); ++n) {
170 if (
m_hProf[l][lr][al][th]->GetBinEntries(n) < 5 &&
m_hProf[l][lr][al][th]->GetBinEntries(n) > 1) {
171 m_hProf[l][lr][al][th]->SetBinError(n, 0.3 /
m_hProf[l][lr][al][th]->GetBinEntries(n));
174 m_hProf[l][lr][al][th]->Fit(
"pol1",
"Q",
"", 30, 60);
175 fpol1 = (TF1*)
m_hProf[l][lr][al][th]->GetFunction(
"pol1");
180 p0 = fpol1->GetParameter(0);
181 p1 = fpol1->GetParameter(1);
182 tmin = -1 * p0 / p1 + 15;
212 double p6 =
m_xtPrior[l][lr][ial_old][ith_old][6];
218 xt->setXTParams(
m_xtPrior[l][lr][ial_old][ith_old]);
221 xt->setXTParams(p0, p1, 0., 0., 0., 0., p6,
m_xtPrior[l][lr][ial_old][ith_old][7]);
223 xt->setFitRange(tmin, p6 + 100);
225 xt->setXTParams(p0, p1, 0., 0., 0., 0.,
m_par6[l], 0.0001);
226 xt->setFitRange(tmin,
m_par6[l] + 100);
231 if (xt->isValid() ==
false) {
232 B2WARNING(
"Empty xt");
236 if (xt->getFitStatus() != 1) {
237 B2WARNING(
"Fit failed");
241 if (xt->validate() ==
true) {
243 m_xtFunc[l][lr][al][th] = (TF1*)xt->getXTFunction();
246 m_hist2d_1[l][lr][al][th] = (TH1F*)xt->getFittedHisto();
248 m_hProf[l][lr][al][th] = (TProfile*)xt->getFittedHisto();
265 const double tMax = 500;
266 for (
int l = 0; l < 56; ++l) {
267 for (
int lr = 0; lr < 2; ++lr) {
270 if (
m_fitStatus[l][lr][al][th] == FitStatus::c_OK) {
272 double y = fun->Eval(tMax);
274 B2INFO(
"Strange XT function l " << l <<
" lr " << lr <<
" alpha " << al <<
" theta " << th
275 <<
", replaced by initial one");
276 fun->SetParameters(
m_xtPrior[l][lr][al][th]);
288 int nFitCompleted = 0;
289 for (
int l = 0; l < 56; ++l) {
290 for (
int lr = 0; lr < 2; ++lr) {
293 if (
m_fitStatus[l][lr][al][th] == FitStatus::c_OK) {
301 if (
static_cast<double>(nFitCompleted) / nTotal <
m_threshold) {
302 B2WARNING(
"Less than " <<
m_threshold * 100 <<
" % of XTs were fitted.");
310 B2INFO(
"Prepare calibration of XT");
311 const double rad2deg = 180 / M_PI;
317 array3 alpha = dbXT->getAlphaBin(i);
324 array3 theta = dbXT->getThetaBin(i);
332 B2FATAL(
"Function type before calibration is wrong " <<
m_xtModePrior);
337 B2INFO(
"Function type " <<
m_xtMode);
339 for (
unsigned short iCL = 0; iCL < 56; ++iCL) {
340 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
343 const std::vector<float> params = dbXT->getXtParams(iCL, iLR, iA, iT);
344 unsigned short np = params.size();
345 for (
unsigned short i = 0; i < np; ++i) {
346 m_xtPrior[iCL][iLR][iA][iT][i] = params[i];
356 B2INFO(
"write calibrated XT");
368 const float deg2rad =
static_cast<float>(
Unit::deg);
371 std::array<float, 3> alpha3 = {
m_lowerAlpha[i]* deg2rad,
379 std::array<float, 3> theta3 = {
m_lowerTheta[i]* deg2rad,
390 for (
int l = 0; l < 56; ++l) {
391 for (
int lr = 0; lr < 2; ++lr) {
392 if (
m_fitStatus[l][lr][al][th] != FitStatus::c_OK) {
394 B2DEBUG(21,
"fit failure status = " <<
m_fitStatus[l][lr][al][th]);
395 B2DEBUG(21,
"layer " << l <<
", r " << lr <<
", alpha " <<
m_iAlpha[al] <<
", theta " <<
m_iTheta[th]);
396 B2DEBUG(21,
"number of event: " <<
m_hProf[l][lr][al][th]->GetEntries());
397 if (
m_fitStatus[l][lr][al][th] != FitStatus::c_lowStat) {
399 B2DEBUG(21,
"Probability of fit: " <<
m_xtFunc[l][lr][al][th]->GetProb());
405 for (
int i = 0; i < 8; ++i) {
409 B2FATAL(
"XT mode before/after calibration is different!");
414 for (
int i = 0; i < 8; ++i) {
418 m_xtFunc[l][lr][al][th]->GetParameters(par);
422 std::vector<float> xtbuff;
423 for (
int i = 0; i < 8; ++i) {
424 xtbuff.push_back(par[i]);
439 B2RESULT(
"Successfully Fitted: " << nfitted);
440 B2RESULT(
"Failure Fit: " << nfailure);
450 B2INFO(
"saving histograms");
451 TFile* fout =
new TFile(
m_histName.c_str(),
"RECREATE");
452 TDirectory* top = gDirectory;
454 if (hNDF && hPval && hEvtT0) {
461 TDirectory* Direct[56];
463 for (
int l = 0; l < 56; ++l) {
465 Direct[l] = gDirectory->mkdir(Form(
"lay_%d", l));
470 for (
int lr = 0; lr < 2; ++lr) {
479 m_hProf[l][lr][al][th]->Write();
489 B2RESULT(
" " << nhisto <<
" histograms was stored.");
Database object for xt-relations.
void setThetaBin(const array3 &theta)
Set theta-angle bin (rad)
void outputToFile(std::string fileName) const
Output the contents in test file format.
void setXtParamMode(unsigned short mode)
Set xt parameterization mode.
void setXtParams(const XtID xtID, const std::vector< float > ¶ms)
Set xt parameters for the specified id.
void setAlphaBin(const array3 &alpha)
Set alpha-angle bin (rad)
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
void storeHisto()
Store histogram to file.
void prepare()
Prepare the calibration of XT.
void sanitaryCheck()
Check if there are any wrong xt functions.
double m_par6[56]
boundary parameter for fitting, semi-experiment number
float m_lowerTheta[7]
Lower boundaries of theta bins.
double m_xtPrior[56][2][18][7][8]
parameters of XT before calibration
std::string m_histName
root file name
TF1 * m_xtFunc[56][2][20][10]
XTFunction.
double m_threshold
minimal requirement for the fraction of fitted results
double m_minNdf
minimum ndf required
bool m_LRseparate
Separate LR in calibration or mix.
void createHisto()
Create histogram for calibration.
bool m_debug
run in debug or silent
int m_minEntriesRequired
minimum number of hit per hitosgram.
int m_fitStatus[56][2][20][10]
Fit flag.
TH2F * m_hist2d[56][2][20][10]
2D histo of xt
int m_nAlphaBins
number of alpha bins
int m_nThetaBins
number of theta bins
double m_minPval
minimum pvalue required
TProfile * m_hProf[56][2][20][10]
Profile xt histo.
float m_iAlpha[18]
Represented alpha in alpha bins.
float m_lowerAlpha[18]
Lower boundaries of alpha bins.
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
XTCalibrationAlgorithm()
Constructor.
bool m_useSliceFit
Use slice fit or profile.
bool m_bField
with b field or none
void write()
Store calibrated constant.
float m_upperAlpha[18]
Upper boundaries of alpha bins.
bool m_textOutput
output text file if true
EResult calibrate() override
Run algo on data.
EResult checkConvergence()
Check the convergence of XT fit.
int m_xtModePrior
Mode of xt before calibration; 0 is polynomial;1 is Chebyshev.
TH2F * m_hist2dDraw[56][20][10]
2d histo for draw
float m_upperTheta[7]
Upper boundaries of theta bins.
int m_xtMode
Mode of xt; 0 is polynomial;1 is Chebyshev.
TH1F * m_hist2d_1[56][2][20][10]
1D xt histo, results of slice fit
std::string m_outputFileName
Output xt filename.
float m_iTheta[7]
Represented theta in theta bins.
Class to perform fitting for each xt function.
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_NotEnoughData
Needs more data =2 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.
static const double deg
degree to radians
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.