Belle II Software  release-05-01-25
TimeWalkCalibration.cc
1 #include <cdc/calibration/TimeWalkCalibration.h>
2 #include <cdc/dbobjects/CDCTimeWalks.h>
3 #include <cdc/geometry/CDCGeometryPar.h>
4 #include <cdc/dataobjects/WireID.h>
5 
6 #include <TF1.h>
7 #include <TFile.h>
8 #include <TChain.h>
9 #include <TDirectory.h>
10 #include <TROOT.h>
11 #include <framework/utilities/FileSystem.h>
12 #include <framework/database/DBImportObjPtr.h>
13 #include <framework/database/DBObjPtr.h>
14 #include <framework/database/IntervalOfValidity.h>
15 #include <framework/logging/Logger.h>
16 
17 using namespace std;
18 using namespace Belle2;
19 using namespace CDC;
20 
21 TimeWalkCalibration::TimeWalkCalibration()
22 {
23  // setDescription("CDC Time Walk Calibration");
24 }
25 void TimeWalkCalibration::CreateHisto()
26 {
27  double x;
28  double t_mea;
29  double w;
30  double t_fit;
31  double ndf;
32  double Pval;
33  unsigned short adc;
34  // int board;
35  int IWire;
36  int lay;
37 
38  // auto& tree = getObject<TTree>("cdcCalib");
39  TChain* tree = new TChain("tree");
40  tree->Add(m_InputRootFileName.c_str());
41 
42  tree->SetBranchAddress("lay", &lay);
43  tree->SetBranchAddress("IWire", &IWire);
44  tree->SetBranchAddress("x_u", &x);
45  tree->SetBranchAddress("t", &t_mea);
46  tree->SetBranchAddress("t_fit", &t_fit);
47  tree->SetBranchAddress("weight", &w);
48  tree->SetBranchAddress("ndf", &ndf);
49  tree->SetBranchAddress("Pval", &Pval);
50  tree->SetBranchAddress("adc", &adc);
51 
52  /* Disable unused branch */
53  std::vector<TString> list_vars = {"lay", "IWire", "x_u", "t", "t_fit", "weight", "Pval", "ndf", "adc"};
54  tree->SetBranchStatus("*", 0);
55 
56  for (TString brname : list_vars) {
57  tree->SetBranchStatus(brname, 1);
58  }
59 
60 
61  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
62  double halfCSize[56];
63  for (int i = 0; i < 56; ++i) {
64  double R = cdcgeo.senseWireR(i);
65  double nW = cdcgeo.nWiresInLayer(i);
66  halfCSize[i] = M_PI * R / nW;
67  }
68  int nEntries = tree->GetEntries();
69  B2INFO("Number of entry: " << nEntries);
70  //Histogram of each board
71  for (int i = 0; i < 300; ++i) {
72  m_h2[i] = new TH2D(Form("board_%d", i), Form("board_%d", i), 50, 0., 500, 60, -30, 30);
73  }
74 
75  //read data
76  for (int i = 0; i < nEntries; ++i) {
77  tree->GetEntry(i);
78  double xmax = halfCSize[lay] - 0.12;
79  if ((fabs(x) < m_xmin) || (fabs(x) > xmax)
80  || (ndf < m_ndfmin)
81  || (Pval < m_Pvalmin)) continue; /*select good region*/
82 
83  m_h2[cdcgeo.getBoardID(WireID(lay, IWire))]->Fill(adc, fabs(t_mea) - fabs(t_fit));
84  }
85  B2INFO("Reading data and filling histogram are done");
86 }
87 
88 bool TimeWalkCalibration::calibrate()
89 {
90  gROOT->SetBatch(1);
91  readTW();
92  CreateHisto();
93  TF1* fold = nullptr;
94  if (m_twParamMode_old == 0)
95  fold = new TF1("fold", "[0]/sqrt(x)", 0, 500);
96  else if (m_twParamMode_old == 1)
97  fold = new TF1("fold", "[0]*exp(-1*[1]*x)", 0, 500);
98  else
99  B2FATAL("Mode " << m_twParamMode_new << " haven't implemented yet.");
100 
101  B2INFO("Old time walk formular: ");
102  fold->Print();
103  B2INFO("Time walk mode new: " << m_twParamMode_new << " with " << m_nTwParams_new << " parameters");
104  // double p3, p1;
105  m_tw_new[0].resize(m_nTwParams_new, 0.0); //for board 0, no available
106  for (int ib = 1; ib < 300; ++ib) {
107  m_flag[ib] = 1;
108  B2DEBUG(199, "Board ID:" << ib);
109  m_h2[ib]->SetDirectory(0);
110 
111  // ignore if histogram is low statistic
112  if (m_h2[ib]->GetEntries() < 500) { m_flag[ib] = 0; continue;}
113  // do slide fit
114  m_h2[ib]->FitSlicesY(0, 1, -1, 10);
115  TString name = m_h2[ib]->GetName();
116  TString hm_name = name + "_1";
117  m_h1[ib] = (TH1D*)gDirectory->Get(hm_name);
118  if (!m_h1[ib]) {m_flag[ib] = 0; continue;}
119  m_h1[ib]->SetDirectory(0);
120 
121  //if number of poist after slide fit quite low, ignore
122  if (m_h1[ib]->GetEntries() < 5) {
123  m_flag[ib] = 0;
124  B2WARNING("Low statistic, number of points after slice fit: " << m_h1[ib]->GetEntries());
125  continue;
126  }
127 
128  // Add previous correction to this
129  for (int p = 0; p < m_nTwParams_old; ++p) {
130  fold->SetParameter(p, m_tw_old[ib][p]);
131  }
132  m_h1[ib]->Add(fold);
133  if (m_twParamMode_new == 0) {
134  TF1* func = new TF1("func", "[0]+[1]/sqrt(x)", 0, 500);
135  func->SetParameters(-4, 28);
136  m_h1[ib]->Fit("func", "MQ", "", 20, 150);
137  }//fit for mode 0
138  else if (m_twParamMode_new == 1) {
139  fitToExponentialFunc(m_h1[ib]);//do fitting for mode 1
140  } else {
141  B2FATAL("Mode " << m_twParamMode_new << " is not available, please check again");
142  }
143 
144  // read fitted parameters
145  TF1* f1 = m_h1[ib]->GetFunction("func");
146  if (!f1) { m_flag[ib] = 0; continue;}
147  m_constTerm[ib] = f1->GetParameter(0);
148  m_tw_new[ib].resize(m_nTwParams_new, 0.0);
149  for (int i = 1; i <= m_nTwParams_new; ++i) {
150  m_tw_new[ib][i - 1] = f1->GetParameter(i);
151  }
152  B2DEBUG(199, "Prob of fitting:" << f1->GetProb());
153  B2DEBUG(199, "Fitting Param 0-1:" << f1->GetParameter(0) << " - " << f1->GetParameter(1));
154  }
155  //Write histogram to file
156  if (m_storeHisto) {
157  B2INFO("Storing histogram");
158  TFile* fhist = new TFile("tw_histo.root", "recreate");
159  TDirectory* old = gDirectory;
160  TDirectory* h1D = old->mkdir("h1D");
161  TDirectory* h2D = old->mkdir("h2D");
162  h1D->cd();
163  for (int ib = 1; ib < 300; ++ib) {
164  if (!m_h1[ib] || m_flag[ib] != 1) continue;
165  m_h1[ib]->SetDirectory(0);
166  if (m_h1[ib]->GetEntries() < 5) continue;
167 
168  m_h1[ib]->SetMinimum(-5);
169  m_h1[ib]->SetMaximum(15);
170  m_h1[ib]->Write();
171  }
172  B2DEBUG(199, "Store 2D histogram");
173  h2D->cd();
174  for (int ib = 1; ib < 300; ++ib) {
175  if (m_h2[ib]) {
176  m_h2[ib]->SetDirectory(0);
177  m_h2[ib]->Write();
178  }
179  }
180  fhist->Close();
181  B2INFO("All hitograms were stored");
182  }
183  Write();
184  updateT0();
185  return true;
186 }
187 void TimeWalkCalibration::Write()
188 {
189  B2INFO("Store calibrated time walk parameters");
190  int nfailure(-1);//for no exist board; board 0;
191  // CDCTimeWalks* dbTw = new CDCTimeWalks();
192  DBImportObjPtr<CDCTimeWalks> dbTw;
193  dbTw.construct();
194 
195  for (int ib = 0; ib < 300; ++ib) {
196  // If fitting fail and new mode is same as previous input mode. use old param
197  if (m_flag[ib] != 1) {
198  nfailure += 1;
199  if (m_twParamMode_old == m_twParamMode_new) {
200  dbTw->setTimeWalkParams(ib, m_tw_old[ib]);
201  } else {
202  //and even calibrated fail but mode is different from previous.
203  //in this case, param is zero
204  m_tw_new[ib].resize(m_nTwParams_new, 0.0);
205  dbTw->setTimeWalkParams(ib, m_tw_new[ib]);
206  }
207  } else {
208  //Use new param if board is successfuly calibrated
209  dbTw->setTimeWalkParams(ib, m_tw_new[ib]);
210  }
211  }
212  dbTw->setTwParamMode(m_twParamMode_new);
213  dbTw->outputToFile(m_outputTWFileName);
214  if (m_useDB) {
215  IntervalOfValidity iov(m_firstExperiment, m_firstRun,
216  m_lastExperiment, m_lastRun);
217  dbTw.import(iov);
218  }
219 
220  B2RESULT("Failure to calibrate time-walk for " << nfailure << " board");
221  B2RESULT("Time-walk coeff. table wrote to " << m_outputTWFileName.c_str());
222 }
223 void TimeWalkCalibration::updateT0()
224 {
225  B2INFO("Add constant term into T0 database");
226  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
227  ofstream ofs(m_outputT0FileName.c_str());
228  DBImportObjPtr<CDCTimeZeros> tz;// = new CDCTimeZeros();
229  tz.construct();
230  double T0;
231  //correct T0 and write
232  for (int ilay = 0; ilay < 56; ++ilay) {
233  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
234  WireID wireid(ilay, iwire);
235  int bID = cdcgeo.getBoardID(wireid);
236  T0 = cdcgeo.getT0(wireid);
237  if (m_flag[bID] == 1) {
238  ofs << ilay << "\t" << iwire << "\t" << T0 - m_constTerm[bID] << std::endl;
239  if (m_useDB)
240  tz->setT0(wireid, T0 - m_constTerm[bID]);
241  } else {
242  ofs << ilay << "\t" << iwire << "\t" << T0 << std::endl;
243  if (m_useDB)
244  tz->setT0(wireid, T0 - m_constTerm[bID]);
245  }
246  }
247  }
248  ofs.close();
249  if (m_useDB) {
250  IntervalOfValidity iov(m_firstExperiment, m_firstRun,
251  m_lastExperiment, m_lastRun);
252  tz.import(iov);
253  }
254  B2RESULT("updated constant term of tw correction to T0 constant.");
255 }
256 void TimeWalkCalibration::fitToExponentialFunc(TH1D* h1)
257 {
258  h1->SetDirectory(0);
259  int max = h1->GetMaximumBin();
260  double maxX = h1->GetBinCenter(max);
261  double maxY = h1->GetBinContent(max);
262  B2DEBUG(199, "Max: id - x - y : " << max << " " << maxX << " " << maxY);
263 
264  //search for p0
265  double p0 = -1;
266  h1->Fit("pol0", "MQE", "", maxX + 125, 500);
267  if (h1->GetFunction("pol0")) {
268  p0 = h1->GetFunction("pol0")->GetParameter(0);
269  }
270  //create histo = old-p0;
271  // fit with expo function to find intial parameters
272  TH1D* hshift = new TH1D("hshift", "shift", h1->GetNbinsX(), 0, 500);
273  hshift->SetDirectory(0);
274  for (int i = 0; i <= h1->GetNbinsX(); ++i) {
275  hshift->SetBinContent(i, h1->GetBinContent(i) - p0);
276  }
277  hshift->Fit("expo", "MQE", "", 0, 400);
278  double p1 = maxY + p0;
279  double p2 = -0.04;
280  if (hshift->GetFunction("expo")) {
281  p1 = exp(hshift->GetFunction("expo")->GetParameter(0));
282  p2 = hshift->GetFunction("expo")->GetParameter(1);
283  }
284 
285  // fit with function
286  TF1* func = new TF1("func", "[0]+ [1]*exp(-1*[2]*x)", 0, 500);
287  func->SetParameters(p0, p1, -1 * p2);
288  func->SetParLimits(0, -5, 5);
289  func->SetParLimits(1, -5, 500);
290  func->SetParLimits(2, -0.01, 0.1);
291  h1->Fit("func", "MQE", "", 10, 400);
292 }
293 void TimeWalkCalibration::readTW()
294 {
295  if (m_useDB) {
296  DBObjPtr<CDCTimeWalks> oldDB;
297  m_twParamMode_old = oldDB->getTwParamMode();
298  for (int ib = 0; ib < 300; ++ib) {
299  m_tw_old[ib] = oldDB->getTimeWalkParams(ib);
300  }
301  } else {
302  //For text mode
303  std::string fileName1 = "/data/cdc" + m_InputTWFileName;
304  std::string fileName = FileSystem::findFile(fileName1);
305  ifstream ifs;
306  if (fileName == "") {
307  fileName = FileSystem::findFile(m_InputTWFileName);
308  }
309  if (fileName == "") {
310  B2FATAL("CDCGeometryPar: " << fileName1 << " not exist!");
311  } else {
312  B2INFO("Time Walk Calibration: open " << fileName);
313  ifs.open(fileName.c_str());
314  if (!ifs) B2FATAL("CDCGeometryPar: cannot open " << fileName1 << " !");
315  }
316 
317  unsigned short iBoard(0);
318  int nRead(0);
319  double dumy(0.);// dumy2(0.), dumy3(0.);
320  ifs >> m_twParamMode_old >> m_nTwParams_old;
321  while (ifs >> iBoard) {
322  m_tw_old[iBoard].resize(m_nTwParams_old);
323  for (int i = 0; i < m_nTwParams_old; ++i) {
324  ifs >> dumy;
325  m_tw_old[iBoard][i] = dumy;
326  }
327  ++nRead;
328  }
329  if (nRead - 1 != 299) B2FATAL("#lines read-in (=" << nRead << ") is not equal #boards 299 !");
330  ifs.close();
331  B2INFO("TW from database: mode = " << m_twParamMode_old << " with " << m_nTwParams_old << " parameters");
332  }
333 }
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:35
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:44
Belle2::CDCTimeZeros::setT0
void setT0(unsigned short iCLayer, unsigned short iWire, float t0)
Set t0 in the list.
Definition: CDCTimeZeros.h:50
Belle2::DBImportObjPtr::construct
void construct(Args &&... params)
Construct an object of type T in this DBImportObjPtr using the provided constructor arguments.
Definition: DBImportObjPtr.h:57
Belle2::DBImportBase::import
bool import(const IntervalOfValidity &iov)
Import the object to database.
Definition: DBImportBase.cc:38
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:75
Belle2::CDC::CDCGeometryPar::senseWireR
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
Definition: CDCGeometryPar.h:1231
Belle2::CDCTimeWalks::setTimeWalkParams
void setTimeWalkParams(unsigned short boardID, const std::vector< float > &params)
Set the time-walk coefficients in the list.
Definition: CDCTimeWalks.h:58
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::DBImportObjPtr
Class for importing a single object to the database.
Definition: DBImportObjPtr.h:33
Belle2::CDC::CDCGeometryPar::getBoardID
unsigned short getBoardID(const WireID &wID) const
Returns frontend board id. corresponding to the wire id.
Definition: CDCGeometryPar.h:576
Belle2::CDC::CDCGeometryPar::getT0
float getT0(const WireID &wireID) const
Returns t0 parameter of the specified sense wire.
Definition: CDCGeometryPar.h:565
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1211
Belle2::CDCTimeWalks::setTwParamMode
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...
Definition: CDCTimeWalks.h:48
Belle2::CDCTimeWalks::outputToFile
void outputToFile(std::string fileName) const
Output the contents in text file format.
Definition: CDCTimeWalks.h:144