Belle II Software  release-05-01-25
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  //for each channel
73  for (int il = 0; il < 56; il++) {
74  const int nW = cdcgeo.nWiresInLayer(il);
75  for (int ic = 0; ic < nW; ++ic) {
76  m_h1[il][ic] = new TH1F(Form("h1%d@%d", il, ic), Form("L%d_cell%d", il, ic), 30, -15, 15);
77  }
78  }
79  //for each board
80  for (int ib = 0; ib < 300; ++ib) {
81  m_hT0b[ib] = new TH1F(Form("hT0b%d", ib), Form("boardID_%d", ib), 100, -20, 20);
82  }
83  //read data
84  const Long64_t nEntries = tree->GetEntries();
85  B2INFO("Number of entries: " << nEntries);
86  for (Long64_t i = 0; i < nEntries; ++i) {
87  tree->GetEntry(i);
88  double xmax = halfCSize[lay] - 0.1;
89  if ((fabs(x) < m_xmin) || (fabs(x) > xmax)
90  || (ndf < m_ndfmin)
91  || (Pval < m_Pvalmin)) continue; /*select good region*/
92  //each channel
93  m_hTotal->Fill(t_mea - t_fit);
94  m_h1[lay][IWire]->Fill(t_mea - t_fit);
95  //each board
96  int boardID = cdcgeo.getBoardID(WireID(lay, IWire));
97  m_hT0b[boardID]->Fill(t_mea - t_fit);
98  }
99  B2INFO("Finish making histogram for all channels");
100 }
101 
102 CalibrationAlgorithm::EResult T0CalibrationAlgorithm::calibrate()
103 {
104  B2INFO("Start calibration");
105 
106  gROOT->SetBatch(1);
107  gErrorIgnoreLevel = 3001;
108 
109  // We are potentially using data from several runs at once during execution
110  // (which may have different DBObject values). So in general you would need to
111  // average them, or aply them to the correct collector data.
112 
113  // However since this is the geometry lets assume it is fixed for now.
114  const auto exprun = getRunList()[0];
115  B2INFO("ExpRun used for DB Geometry : " << exprun.first << " " << exprun.second);
116  updateDBObjPtrs(1, exprun.second, exprun.first);
117 
118  // CDCGeometryPar basically constructs a ton of objects and other DB objects.
119  // Normally we'd call updateDBObjPtrs to set the values of the requested DB objects.
120  // But in CDCGeometryPar the DB objects get used during the constructor so they must
121  // be set before/during the constructor.
122 
123  // Since we are avoiding using the DataStore EventMetaData, we need to pass in
124  // an EventMetaData object to be used when constructing the DB objects.
125 
126  B2INFO("Creating CDCGeometryPar object");
127  CDC::CDCGeometryPar::Instance(&(*m_cdcGeo));
128 
129  auto hEvtT0 = getObjectPtr<TH1F>("hEventT0");
130  double dEventT0 = hEvtT0->GetMean();
131  createHisto();
132  TH1F* hm_All = new TH1F("hm_All", "mean of #DeltaT distribution for all chanels", 500, -10, 10);
133  TH1F* hs_All = new TH1F("hs_All", "#sigma of #DeltaT distribution for all chanels", 100, 0, 10);
134  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
135 
136  TF1* g1 = new TF1("g1", "gaus", -100, 100);
137  vector<double> b, db, Sb, dSb;
138  vector<double> c[56];
139  vector<double> dc[56];
140  vector<double> s[56];
141  vector<double> ds[56];
142 
143  B2INFO("Gaus fitting for whole channel");
144  double par[3];
145  m_hTotal->SetDirectory(0);
146  double mean = m_hTotal->GetMean();
147  m_hTotal->Fit("g1", "Q", "", mean - 15, mean + 15);
148  g1->GetParameters(par);
149 
150  B2INFO("Gaus fitting for each board");
151  for (int ib = 1; ib < 300; ++ib) {
152  if (m_hT0b[ib]->GetEntries() < 10) continue;
153  mean = m_hT0b[ib]->GetMean();
154  m_hT0b[ib]->SetDirectory(0);
155  m_hT0b[ib]->Fit("g1", "Q", "", mean - 15, mean + 15);
156  g1->GetParameters(par);
157  b.push_back(ib);
158  db.push_back(0.0);
159  Sb.push_back(par[1]);
160  dSb.push_back(g1->GetParError(1));
161  dtb[ib] = par[1];
162  err_dtb[ib] = g1->GetParError(1);
163  }
164  B2INFO("Gaus fitting for each cell");
165  for (int ilay = 0; ilay < 56; ++ilay) {
166  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
167  const int n = m_h1[ilay][iwire]->GetEntries();
168  B2DEBUG(21, "layer " << ilay << " wire " << iwire << " entries " << n);
169  if (n < 10) continue;
170  mean = m_h1[ilay][iwire]->GetMean();
171  m_h1[ilay][iwire]->SetDirectory(0);
172  m_h1[ilay][iwire]->Fit("g1", "Q", "", mean - 15, mean + 15);
173  g1->GetParameters(par);
174  c[ilay].push_back(iwire);
175  dc[ilay].push_back(0.0);
176  s[ilay].push_back(par[1]);
177  ds[ilay].push_back(g1->GetParError(1));
178 
179  dt[ilay][iwire] = par[1];
180  err_dt[ilay][iwire] = g1->GetParError(1);
181  hm_All->Fill(par[1]);// mean of gauss fitting.
182  hs_All->Fill(par[2]); // sigma of gauss fitting.
183  }
184  }
185 
186  // mean shift
187  // const double dt0Mean = hm_All->GetMean();
188  for (int ilay = 0; ilay < 56; ++ilay) {
189  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
190  dt[ilay][iwire] += dEventT0;
191  }
192  }
193 
194 
195 
196  if (m_storeHisto) {
197  B2INFO("Storing histograms");
198  auto hNDF = getObjectPtr<TH1F>("hNDF");
199  auto hPval = getObjectPtr<TH1F>("hPval");
200  TFile* fout = new TFile(m_histName.c_str(), "RECREATE");
201  fout->cd();
202  TGraphErrors* gr[56];
203  TDirectory* top = gDirectory;
204 
205  //store NDF, P-val. EventT0 histogram for monitoring during calibration
206  if (hNDF && hPval && hEvtT0) {
207  hEvtT0->Write();
208  hPval->Write();
209  hNDF->Write();
210  }
211  m_hTotal->Write();
212  hm_All->Write();
213  hs_All->Write();
214  TDirectory* subDir[56];
215  for (int ilay = 0; ilay < 56; ++ilay) {
216  subDir[ilay] = top ->mkdir(Form("lay_%d", ilay));
217  subDir[ilay]->cd();
218  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
219  m_h1[ilay][iwire]->Write();
220  }
221  }
222  top->cd();
223  TDirectory* corrT0 = top->mkdir("DeltaT0");
224  corrT0->cd();
225 
226  if (b.size() > 2) {
227  TGraphErrors* grb = new TGraphErrors(b.size(), &b.at(0), &Sb.at(0), &db.at(0), &dSb.at(0));
228  grb->SetMarkerColor(2);
229  grb->SetMarkerSize(1.0);
230  grb->SetTitle("#DeltaT0;BoardID;#DeltaT0[ns]");
231  grb->SetMaximum(10);
232  grb->SetMinimum(-10);
233  grb->SetName("Board");
234  grb->Write();
235  }
236  for (int sl = 0; sl < 56; ++sl) {
237  if (c[sl].size() < 2) continue;
238  gr[sl] = new TGraphErrors(c[sl].size(), &c[sl].at(0), &s[sl].at(0), &dc[sl].at(0), &ds[sl].at(0));
239  gr[sl]->SetMarkerColor(2);
240  gr[sl]->SetMarkerSize(1.0);
241  gr[sl]->SetTitle(Form("Layer_%d;IWire;#LT t_{mea}-t_{fit} #GT [ns]", sl));
242  gr[sl]->SetMaximum(10);
243  gr[sl]->SetMinimum(-10);
244  gr[sl]->SetName(Form("lay%d", sl));
245  gr[sl]->Write();
246  }
247  fout->Close();
248  }
249  B2INFO("Writing constants");
250  write();
251 
252 
253  if (fabs(hm_All->GetMean()) < m_maxMeanDt && fabs(hm_All->GetRMS()) < m_maxRMSDt) {
254  B2INFO("mean " << fabs(hm_All->GetMean()) << " " << m_maxMeanDt);
255  B2INFO("sigma " << fabs(hm_All->GetRMS()) << " " << m_maxRMSDt);
256  return c_OK;
257  } else {
258  B2INFO("mean " << fabs(hm_All->GetMean()) << " " << m_maxMeanDt);
259  B2INFO("sigma " << fabs(hm_All->GetRMS()) << " " << m_maxRMSDt);
260  return c_Iterate;
261  }
262 }
263 
264 void T0CalibrationAlgorithm::write()
265 {
266  CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
267  CDCTimeZeros* tz = new CDCTimeZeros();
268  double T0;
269  TH1F* T0B[300];
270  for (int ib = 0; ib < 300; ++ib) {
271  T0B[ib] = new TH1F(Form("T0B%d", ib), Form("boardID_%d", ib), 9000, 0, 9000);
272  }
273  //for calculate T0 mean of each board
274  for (int ilay = 0; ilay < 56; ++ilay) {
275  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
276  WireID wireid(ilay, iwire);
277  T0 = cdcgeo.getT0(wireid);
278  T0B[cdcgeo.getBoardID(wireid)]->Fill(T0);
279  }
280  }
281 
282  //correct T0 and write
283  for (int ilay = 0; ilay < 56; ++ilay) {
284  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
285  WireID wireid(ilay, iwire);
286  int bID = cdcgeo.getBoardID(wireid);
287  if (abs(err_dt[ilay][iwire]) > 2 || abs(dt[ilay][iwire]) > 1.e3) {
288  T0 = T0B[bID]->GetMean();
289  dt[ilay][iwire] = dtb[bID];
290  } else {
291  T0 = cdcgeo.getT0(wireid);
292  }
293  tz->setT0(wireid, T0 - dt[ilay][iwire]);
294  }
295  }
296 
297  if (m_textOutput == true) {
298  tz->outputToFile(m_outputT0FileName);
299  }
300  saveCalibration(tz, "CDCTimeZeros");
301 }
Belle2::CDC::T0CalibrationAlgorithm::m_Pvalmin
double m_Pvalmin
minimum pvalue required
Definition: T0CalibrationAlgorithm.h:75
Belle2::CDC::T0CalibrationAlgorithm::m_ndfmin
double m_ndfmin
minimum ndf required
Definition: T0CalibrationAlgorithm.h:74
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:78
Belle2::CDC::T0CalibrationAlgorithm::err_dt
double err_dt[56][385]
error of dt of each channel
Definition: T0CalibrationAlgorithm.h:80
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:85
Belle2::CDC::T0CalibrationAlgorithm::m_hTotal
TH1F * m_hTotal
1D histogram of delta T whole channel
Definition: T0CalibrationAlgorithm.h:70
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:77
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:87
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:73
Belle2::CDC::T0CalibrationAlgorithm::m_h1
TH1F * m_h1[56][385]
1D histogram for each channel
Definition: T0CalibrationAlgorithm.h:71
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:88
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:79
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:89
Belle2::CDC::T0CalibrationAlgorithm::m_textOutput
bool m_textOutput
output text file if true
Definition: T0CalibrationAlgorithm.h:86
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:81
Belle2::CDC::T0CalibrationAlgorithm::err_dtb
double err_dtb[300]
error of dt of board
Definition: T0CalibrationAlgorithm.h:82
Belle2::CDC::T0CalibrationAlgorithm::calibrate
EResult calibrate() override
Run algo on data.
Definition: T0CalibrationAlgorithm.cc:102
Belle2::CDC::T0CalibrationAlgorithm::write
void write()
write outut or store db
Definition: T0CalibrationAlgorithm.cc:264
Belle2::CDC::T0CalibrationAlgorithm::m_hT0b
TH1F * m_hT0b[300]
1D histogram for each board
Definition: T0CalibrationAlgorithm.h:72