Belle II Software  release-05-02-19
T0CalibrationAlgorithm.cc
1 #include <cdc/calibration/T0CalibrationAlgorithm.h>
2 #include <calibration/CalibrationAlgorithm.h>
3 #include <cdc/dbobjects/CDCTimeZeros.h>
4 #include <cdc/dataobjects/WireID.h>
5 #include <cdc/geometry/CDCGeometryPar.h>
6 
7 #include <TError.h>
8 #include <TROOT.h>
9 #include <TGraphErrors.h>
10 #include <TF1.h>
11 #include <TFile.h>
12 #include <TTree.h>
13 
14 #include <framework/database/DBObjPtr.h>
15 #include <framework/database/IntervalOfValidity.h>
16 #include <framework/logging/Logger.h>
17 
18 using namespace std;
19 using namespace Belle2;
20 using namespace CDC;
21 
22 T0CalibrationAlgorithm::T0CalibrationAlgorithm(): CalibrationAlgorithm("CDCCalibrationCollector")
23 {
24 
25  setDescription(
26  " -------------------------- T0 Calibration Algorithm -------------------------\n"
27  );
28 }
29 
30 void T0CalibrationAlgorithm::createHisto()
31 {
32 
33  B2INFO("Creating histograms");
34  Float_t x;
35  Float_t t_mea;
36  Float_t t_fit;
37  Float_t ndf;
38  Float_t Pval;
39  UShort_t IWire;
40  UChar_t lay;
41 
42  auto tree = getObjectPtr<TTree>("tree");
43  tree->SetBranchAddress("lay", &lay);
44  tree->SetBranchAddress("IWire", &IWire);
45  tree->SetBranchAddress("x_u", &x);
46  tree->SetBranchAddress("t", &t_mea);
47  tree->SetBranchAddress("t_fit", &t_fit);
48  tree->SetBranchAddress("ndf", &ndf);
49  tree->SetBranchAddress("Pval", &Pval);
50 
51 
52  /* Disable unused branch */
53  std::vector<TString> list_vars = {"lay", "IWire", "x_u", "t", "t_fit", "Pval", "ndf"};
54  tree->SetBranchStatus("*", 0);
55 
56  for (TString brname : list_vars) {
57  tree->SetBranchStatus(brname, 1);
58  }
59 
60 
61  double halfCSize[56];
62 
63  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
64 
65  for (int i = 0; i < 56; ++i) {
66  double R = cdcgeo.senseWireR(i);
67  double nW = cdcgeo.nWiresInLayer(i);
68  halfCSize[i] = M_PI * R / nW;
69  }
70 
71  m_hTotal = new TH1F("hTotal", "hTotal", 30, -15, 15);
72 
73  //for each channel
74  for (int il = 0; il < 56; il++) {
75  const int nW = cdcgeo.nWiresInLayer(il);
76  for (int ic = 0; ic < nW; ++ic) {
77  m_h1[il][ic] = new TH1F(Form("hdT_L%d_W%d", il, ic), Form("L%d_cell%d", il, ic), 30, -15, 15);
78  }
79  }
80 
81  //for each board
82  for (int ib = 0; ib < 300; ++ib) {
83  m_hT0b[ib] = new TH1F(Form("hdT_b%d", ib), Form("boardID_%d", ib), 100, -20, 20);
84  }
85 
86  //read data
87  const Long64_t nEntries = tree->GetEntries();
88  B2INFO("Number of entries: " << nEntries);
89  for (Long64_t i = 0; i < nEntries; ++i) {
90  tree->GetEntry(i);
91  double xmax = halfCSize[lay] - 0.1;
92  if ((fabs(x) < m_xmin) || (fabs(x) > xmax)
93  || (ndf < m_ndfmin)
94  || (Pval < m_Pvalmin)) continue; /*select good region*/
95  //each channel
96  m_hTotal->Fill(t_mea - t_fit);
97  m_h1[lay][IWire]->Fill(t_mea - t_fit);
98  //each board
99  int boardID = cdcgeo.getBoardID(WireID(lay, IWire));
100  m_hT0b[boardID]->Fill(t_mea - t_fit);
101  }
102 
103  B2INFO("Finish making histogram for all channels");
104 }
105 
106 CalibrationAlgorithm::EResult T0CalibrationAlgorithm::calibrate()
107 {
108  B2INFO("Start calibration");
109 
110  gROOT->SetBatch(1);
111  gErrorIgnoreLevel = 3001;
112 
113  // We are potentially using data from several runs at once during execution
114  // (which may have different DBObject values). So in general you would need to
115  // average them, or aply them to the correct collector data.
116 
117  // However since this is the geometry lets assume it is fixed for now.
118  const auto exprun = getRunList()[0];
119  B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
120  updateDBObjPtrs(1, exprun.second, exprun.first);
121 
122  // CDCGeometryPar basically constructs a ton of objects and other DB objects.
123  // Normally we'd call updateDBObjPtrs to set the values of the requested DB objects.
124  // But in CDCGeometryPar the DB objects get used during the constructor so they must
125  // be set before/during the constructor.
126 
127  // Since we are avoiding using the DataStore EventMetaData, we need to pass in
128  // an EventMetaData object to be used when constructing the DB objects.
129 
130  B2INFO("Creating CDCGeometryPar object");
131  CDC::CDCGeometryPar::Instance(&(*m_cdcGeo));
132 
133  auto hEvtT0 = getObjectPtr<TH1F>("hEventT0");
134  double dEventT0 = hEvtT0->GetMean();
135  createHisto();
136  TH1F* hm_All = new TH1F("hm_All", "mean of #DeltaT distribution for all chanels", 500, -10, 10);
137  TH1F* hs_All = new TH1F("hs_All", "#sigma of #DeltaT distribution for all chanels", 100, 0, 10);
138  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
139 
140  TF1* g1 = new TF1("g1", "gaus", -100, 100);
141  g1->SetParLimits(1, -20, 20);
142  vector<double> b, db, Sb, dSb;
143  vector<double> c[56];
144  vector<double> dc[56];
145  vector<double> s[56];
146  vector<double> ds[56];
147 
148  B2INFO("Gaus fitting for whole channel");
149  double par[3];
150  m_hTotal->SetDirectory(0);
151  double mean = m_hTotal->GetMean();
152  m_hTotal->Fit("g1", "Q", "", mean - 15, mean + 15);
153  g1->GetParameters(par);
154 
155  B2INFO("Gaus fitting for each board");
156  for (int ib = 1; ib < 300; ++ib) {
157  // Set Delta_T0=0 again to make sure there is no strange case
158  // in which T0 might be initialed with a strange value
159  dtb[ib] = 0;
160  err_dtb[ib] = 0;
161 
162  if (m_hT0b[ib]->Integral(1, m_hT0b[ib]->GetNbinsX()) < 50) {
163  //set error to large number as a flag of bad value
164  err_dtb[ib] = 50; continue;
165  }
166  mean = m_hT0b[ib]->GetMean();
167  m_hT0b[ib]->SetDirectory(0);
168  m_hT0b[ib]->Fit("g1", "Q", "", mean - 15, mean + 15);
169  g1->GetParameters(par);
170  b.push_back(ib);
171  db.push_back(0.0);
172  Sb.push_back(par[1]);
173  dSb.push_back(g1->GetParError(1));
174  dtb[ib] = par[1] + dEventT0;// add dEvtT0 here
175  err_dtb[ib] = g1->GetParError(1);
176  }
177  B2INFO("Gaus fitting for each cell");
178  for (int ilay = 0; ilay < 56; ++ilay) {
179  const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
180  for (unsigned int iwire = 0; iwire < nW; ++iwire) {
181 
182  const int n = m_h1[ilay][iwire]->Integral(1, m_h1[ilay][iwire]->GetNbinsX()) ;
183  B2DEBUG(21, "layer " << ilay << " wire " << iwire << " entries " << n);
184  // Set Delta_T0=0 again to make sure there is no strange case
185  // in which T0 might be initialed with a strange value
186  dt[ilay][iwire] = 0;
187  err_dt[ilay][iwire] = 0;
188 
189  if (n < 30) {
190  //set error to large number as a flag of bad value
191  err_dt[ilay][iwire] = 50; continue;
192  }
193 
194  mean = m_h1[ilay][iwire]->GetMean();
195  m_h1[ilay][iwire]->SetDirectory(0);
196  m_h1[ilay][iwire]->Fit("g1", "Q", "", mean - 15, mean + 15);
197  g1->GetParameters(par);
198  c[ilay].push_back(iwire);
199  dc[ilay].push_back(0.0);
200  s[ilay].push_back(par[1]);
201  ds[ilay].push_back(g1->GetParError(1));
202 
203  dt[ilay][iwire] = par[1] + dEventT0; // add dEvtT0 here;
204  err_dt[ilay][iwire] = g1->GetParError(1);
205  hm_All->Fill(par[1]);// mean of gauss fitting.
206  hs_All->Fill(par[2]); // sigma of gauss fitting.
207  }
208  }
209 
210  if (m_storeHisto) {
211  B2INFO("Storing histograms");
212  auto hNDF = getObjectPtr<TH1F>("hNDF");
213  auto hPval = getObjectPtr<TH1F>("hPval");
214  TFile* fout = new TFile(m_histName.c_str(), "RECREATE");
215  fout->cd();
216  TGraphErrors* gr[56];
217  TDirectory* top = gDirectory;
218 
219  //store NDF, P-val. EventT0 histogram for monitoring during calibration
220  if (hNDF && hPval && hEvtT0) {
221  hEvtT0->Write();
222  hPval->Write();
223  hNDF->Write();
224  }
225  m_hTotal->Write();
226  hm_All->Write();
227  hs_All->Write();
228  TDirectory* subDir[56];
229  for (int ilay = 0; ilay < 56; ++ilay) {
230  subDir[ilay] = top ->mkdir(Form("lay_%d", ilay));
231  subDir[ilay]->cd();
232  const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
233  for (unsigned int iwire = 0; iwire < nW; ++iwire) {
234  m_h1[ilay][iwire]->Write();
235  }
236  }
237 
238  top->cd();
239  TDirectory* corrT0 = top->mkdir("DeltaT0");
240  corrT0->cd();
241 
242  if (b.size() > 2) {
243  TGraphErrors* grb = new TGraphErrors(b.size(), &b.at(0), &Sb.at(0), &db.at(0), &dSb.at(0));
244  grb->SetMarkerColor(2);
245  grb->SetMarkerSize(1.0);
246  grb->SetTitle("#DeltaT0;BoardID;#DeltaT0[ns]");
247  grb->SetMaximum(10);
248  grb->SetMinimum(-10);
249  grb->SetName("Board");
250  grb->Write();
251  }
252  for (int sl = 0; sl < 56; ++sl) {
253  if (c[sl].size() < 2) continue;
254  gr[sl] = new TGraphErrors(c[sl].size(), &c[sl].at(0), &s[sl].at(0), &dc[sl].at(0), &ds[sl].at(0));
255  gr[sl]->SetMarkerColor(2);
256  gr[sl]->SetMarkerSize(1.0);
257  gr[sl]->SetTitle(Form("Layer_%d;IWire;#LT t_{mea}-t_{fit} #GT [ns]", sl));
258  gr[sl]->SetMaximum(10);
259  gr[sl]->SetMinimum(-10);
260  gr[sl]->SetName(Form("lay%d", sl));
261  gr[sl]->Write();
262  }
263  fout->Close();
264  }
265  B2INFO("Writing constants");
266  write();
267 
268 
269  if (fabs(hm_All->GetMean()) < m_maxMeanDt && fabs(hm_All->GetRMS()) < m_maxRMSDt) {
270  B2INFO("mean " << fabs(hm_All->GetMean()) << " " << m_maxMeanDt);
271  B2INFO("sigma " << fabs(hm_All->GetRMS()) << " " << m_maxRMSDt);
272  return c_OK;
273  } else {
274  B2INFO("mean " << fabs(hm_All->GetMean()) << " " << m_maxMeanDt);
275  B2INFO("sigma " << fabs(hm_All->GetRMS()) << " " << m_maxRMSDt);
276  return c_Iterate;
277  }
278 }
279 
280 void T0CalibrationAlgorithm::write()
281 {
282  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
283  CDCTimeZeros* tz = new CDCTimeZeros();
284 
285  TH1F* hT0B[300];
286  for (int ib = 0; ib < 300; ++ib) {
287  hT0B[ib] = new TH1F(Form("hT0B%d", ib), Form("boardID_%d", ib), 9000, 0, 9000);
288  }
289 
290  TH1F* hT0_all = new TH1F("hT0_all", "T0 distribution", 9000, 0, 9000);
291  double T0;
292  // get old T0 to calculate T0 mean of each board
293  for (int ilay = 0; ilay < 56; ++ilay) {
294  const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
295  for (unsigned int iwire = 0; iwire < nW; ++iwire) {
296  WireID wireid(ilay, iwire);
297  T0 = cdcgeo.getT0(wireid);
298  hT0_all->Fill(T0);
299  hT0B[cdcgeo.getBoardID(wireid)]->Fill(T0);
300  }
301  }
302 
303  //get Nominal T0
304  double T0_average = getMeanT0(hT0_all);
305  if (fabs(T0_average - m_commonT0) > 50) {B2WARNING("Large difference between common T0 (" << m_commonT0 << ") and aveage value" << T0_average);}
306  // get average T0 for each board and also apply T0 correction for T0 of that board
307  double T0B[300];
308  for (int i = 0; i < 300; ++i) {T0B[i] = m_commonT0;}
309 
310 
311  for (int ib = 1; ib < 300; ++ib) {
312  T0B[ib] = getMeanT0(hT0B[ib]);
313  if (fabs(T0B[ib] - T0_average) > 25) {
314  B2WARNING("T0 of Board " << ib << " (= " << T0B[ib] << ") is too different with common T0: " << T0_average <<
315  "\n It will be replaced by common T0");
316  T0B[ib] = T0_average;
317  continue;
318  }
319  //correct T0 board
320  if (abs(err_dtb[ib]) < 2 && abs(dtb[ib]) < 20) {
321  T0B[ib] -= dtb[ib];
322  }
323  }
324 
325  //correct T0 and write
326  double dT;
327  for (int ilay = 0; ilay < 56; ++ilay) {
328  const unsigned int nW = cdcgeo.nWiresInLayer(ilay);
329  for (unsigned int iwire = 0; iwire < nW; ++iwire) {
330  WireID wireid(ilay, iwire);
331  int bID = cdcgeo.getBoardID(wireid);
332 
333  //get old T0, replace with common T0 of board or average over all channel.
334  T0 = cdcgeo.getT0(wireid);
335  if (fabs(T0 - T0B[bID]) > 25 || fabs(T0 - T0_average) > 25) {
336  B2WARNING("T0 of wireID L-W: " << ilay << "--" << iwire << " (= " << T0
337  << ") is too different with common T0 of board: " << bID << " = " << T0B[bID] <<
338  "\n It will be replaced by common T0 of the Board");
339  T0 = T0B[bID];
340  }
341 
342  //select DeltaT
343  dT = 0;
344  if (abs(err_dt[ilay][iwire]) < 2 && abs(dt[ilay][iwire]) < 20) {
345  dT = dt[ilay][iwire];
346  }
347 
348  //export
349  tz->setT0(wireid, T0 - dT);
350  }
351  }
352 
353  if (m_textOutput == true) {
354  tz->outputToFile(m_outputT0FileName);
355  }
356  saveCalibration(tz, "CDCTimeZeros");
357 }
358 double T0CalibrationAlgorithm::getMeanT0(TH1F* h1)
359 {
360  double mean1 = h1->GetMean();
361  h1->GetXaxis()->SetRangeUser(mean1 - 50, 8000);
362  double mean2 = h1->GetMean();
363  while (fabs(mean1 - mean2) > 0.5) {
364  mean1 = mean2;
365  h1->GetXaxis()->SetRangeUser(mean1 - 50, 8000);
366  mean2 = h1->GetMean();
367  }
368  return mean2;
369 
370 }
Belle2::CDC::T0CalibrationAlgorithm::m_Pvalmin
double m_Pvalmin
minimum pvalue required
Definition: T0CalibrationAlgorithm.h:78
Belle2::CDC::T0CalibrationAlgorithm::m_ndfmin
double m_ndfmin
minimum ndf required
Definition: T0CalibrationAlgorithm.h:77
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:44
Belle2::CDC::T0CalibrationAlgorithm::m_maxRMSDt
double m_maxRMSDt
RMS of dT distribution of all channels.
Definition: T0CalibrationAlgorithm.h:81
Belle2::CDC::T0CalibrationAlgorithm::err_dt
double err_dt[56][385]
error of dt of each channel
Definition: T0CalibrationAlgorithm.h:83
Belle2::CDC::T0CalibrationAlgorithm::m_commonT0
double m_commonT0
A common T0 of all channels.
Definition: T0CalibrationAlgorithm.h:86
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::CDC::T0CalibrationAlgorithm::m_storeHisto
bool m_storeHisto
store histo or not
Definition: T0CalibrationAlgorithm.h:88
Belle2::CDC::T0CalibrationAlgorithm::m_hTotal
TH1F * m_hTotal
1D histogram of delta T whole channel
Definition: T0CalibrationAlgorithm.h:73
Belle2::CDCTimeZeros::setT0
void setT0(unsigned short iCLayer, unsigned short iWire, float t0)
Set t0 in the list.
Definition: CDCTimeZeros.h:50
Belle2::CalibrationAlgorithm::c_Iterate
@ c_Iterate
Needs iteration =1 in Python.
Definition: CalibrationAlgorithm.h:52
Belle2::CDC::T0CalibrationAlgorithm::m_maxMeanDt
double m_maxMeanDt
Mean of dT distribution of all channels;.
Definition: T0CalibrationAlgorithm.h:80
Belle2::CDCTimeZeros::outputToFile
void outputToFile(std::string fileName) const
Output the contents in text file format.
Definition: CDCTimeZeros.h:133
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2::CDC::T0CalibrationAlgorithm::m_outputT0FileName
std::string m_outputT0FileName
output t0 file name for text file
Definition: T0CalibrationAlgorithm.h:90
Belle2::CalibrationAlgorithm::updateDBObjPtrs
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition: CalibrationAlgorithm.cc:398
Belle2::CDC::T0CalibrationAlgorithm::m_xmin
double m_xmin
minimum drift length
Definition: T0CalibrationAlgorithm.h:76
Belle2::CDC::T0CalibrationAlgorithm::m_h1
TH1F * m_h1[56][385]
1D histogram for each channel
Definition: T0CalibrationAlgorithm.h:74
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::CDC::T0CalibrationAlgorithm::m_histName
std::string m_histName
root file name
Definition: T0CalibrationAlgorithm.h:91
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:41
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDC::T0CalibrationAlgorithm::createHisto
void createHisto()
create histo for each channel
Definition: T0CalibrationAlgorithm.cc:30
Belle2::CDCTimeZeros
Database object for timing offset (t0).
Definition: CDCTimeZeros.h:36
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::T0CalibrationAlgorithm::dt
double dt[56][385]
dt of each channel
Definition: T0CalibrationAlgorithm.h:82
Belle2::CDC::T0CalibrationAlgorithm::getMeanT0
double getMeanT0(TH1F *h1)
calculate mean of the T0 distribution
Definition: T0CalibrationAlgorithm.cc:358
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:276
Belle2::CDC::CDCGeometryPar::getT0
float getT0(const WireID &wireID) const
Returns t0 parameter of the specified sense wire.
Definition: CDCGeometryPar.h:565
Belle2::CDC::T0CalibrationAlgorithm::m_cdcGeo
DBObjPtr< CDCGeometry > m_cdcGeo
Geometry of CDC.
Definition: T0CalibrationAlgorithm.h:92
Belle2::CDC::T0CalibrationAlgorithm::m_textOutput
bool m_textOutput
output text file if true
Definition: T0CalibrationAlgorithm.h:89
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1211
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::CDC::T0CalibrationAlgorithm::dtb
double dtb[300]
dt of each board
Definition: T0CalibrationAlgorithm.h:84
Belle2::CDC::T0CalibrationAlgorithm::err_dtb
double err_dtb[300]
error of dt of board
Definition: T0CalibrationAlgorithm.h:85
Belle2::CDC::T0CalibrationAlgorithm::calibrate
EResult calibrate() override
Run algo on data.
Definition: T0CalibrationAlgorithm.cc:106
Belle2::CDC::T0CalibrationAlgorithm::write
void write()
write outut or store db
Definition: T0CalibrationAlgorithm.cc:280
Belle2::CDC::T0CalibrationAlgorithm::m_hT0b
TH1F * m_hT0b[300]
1D histogram for each board
Definition: T0CalibrationAlgorithm.h:75