Belle II Software  release-05-01-25
T0Correction.cc
1 #include <cdc/calibration/T0Correction.h>
2 #include <cdc/dbobjects/CDCTimeZeros.h>
3 #include <cdc/geometry/CDCGeometryPar.h>
4 #include <cdc/dataobjects/WireID.h>
5 
6 #include <TError.h>
7 #include <TROOT.h>
8 #include <TGraphErrors.h>
9 #include <TF1.h>
10 #include <TFile.h>
11 #include <TChain.h>
12 
13 #include <framework/database/IntervalOfValidity.h>
14 #include <framework/database/DBImportObjPtr.h>
15 #include <framework/logging/Logger.h>
16 
17 using namespace std;
18 using namespace Belle2;
19 using namespace CDC;
20 
21 T0Correction::T0Correction():
22  m_firstExperiment(0), m_firstRun(0),
23  m_lastExperiment(-1), m_lastRun(-1)
24 {
25  /*
26  setDescription(
27  " -------------------------- Test Calibration Algoritm -------------------------\n"
28  " \n"
29  " Testing algorithm which just gets mean of a test histogram collected by \n"
30  " CaTest module and provides a DB object with another histogram with one \n"
31  " entry at calibrated value. \n"
32  " ------------------------------------------------------------------------------\n"
33  );
34  */
35 }
36 
38 {
39 
40  B2INFO("CreateHisto");
41  double x;
42  double t_mea;
43  double w;
44  double t_fit;
45  double ndf;
46  double Pval;
47  int IWire;
48  int lay;
49 
50  // auto& tree = getObject<TTree>("cdcCalib");
51  TChain* tree = new TChain("tree");
52  tree->Add(m_inputRootFileName.c_str());
53  tree->SetBranchAddress("lay", &lay);
54  tree->SetBranchAddress("IWire", &IWire);
55  tree->SetBranchAddress("x_u", &x);
56  tree->SetBranchAddress("t", &t_mea);
57  tree->SetBranchAddress("t_fit", &t_fit);
58  tree->SetBranchAddress("weight", &w);
59  tree->SetBranchAddress("ndf", &ndf);
60  tree->SetBranchAddress("Pval", &Pval);
61 
62  /* Disable unused branch */
63  std::vector<TString> list_vars = {"lay", "IWire", "x_u", "t", "t_fit", "weight", "Pval", "ndf"};
64  tree->SetBranchStatus("*", 0);
65 
66  for (TString brname : list_vars) {
67  tree->SetBranchStatus(brname, 1);
68  }
69 
70 
71  double halfCSize[56];
72  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
73  for (int i = 0; i < 56; ++i) {
74  double R = cdcgeo.senseWireR(i);
75  double nW = cdcgeo.nWiresInLayer(i);
76  halfCSize[i] = M_PI * R / nW;
77  }
78 
79  m_hTotal = new TH1F("hTotal", "hTotal", 30, -15, 15);
80  //for each channel
81  for (int il = 0; il < 56; il++) {
82  const int nW = cdcgeo.nWiresInLayer(il);
83  for (int ic = 0; ic < nW; ++ic) {
84  m_h1[il][ic] = new TH1F(Form("h1%d@%d", il, ic), Form("L%d_cell%d", il, ic), 30, -15, 15);
85  }
86  }
87  //for each board
88  for (int ib = 0; ib < 300; ++ib) {
89  m_hT0b[ib] = new TH1F(Form("hT0b%d", ib), Form("boardID_%d", ib), 100, -20, 20);
90  }
91  //read data
92  const int nEntries = tree->GetEntries();
93  B2INFO("Number of entries: " << nEntries);
94  for (int i = 0; i < nEntries; ++i) {
95  tree->GetEntry(i);
96  double xmax = halfCSize[lay] - 0.1;
97  if ((fabs(x) < m_xmin) || (fabs(x) > xmax)
98  || (ndf < m_ndfmin)
99  || (Pval < m_Pvalmin)) continue; /*select good region*/
100  //each channel
101  m_hTotal->Fill(t_mea - t_fit);
102  m_h1[lay][IWire]->Fill(t_mea - t_fit);
103  //each board
104  m_hT0b[cdcgeo.getBoardID(WireID(lay, IWire))]->Fill(t_mea - t_fit);
105  }
106  B2INFO("Finish making histogram for all channels");
107 }
108 
110 {
111  B2INFO("Start calibration");
112 
113  gROOT->SetBatch(1);
114  gErrorIgnoreLevel = 3001;
115 
116  CreateHisto();
117  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
118 
119  TF1* g1 = new TF1("g1", "gaus", -100, 100);
120  vector<double> b, db, Sb, dSb;
121  vector<double> c[56];
122  vector<double> dc[56];
123  vector<double> s[56];
124  vector<double> ds[56];
125 
126  B2INFO("Gaus fitting for whole channel");
127  double par[3];
128  m_hTotal->SetDirectory(0);
129  double mean = m_hTotal->GetMean();
130  m_hTotal->Fit("g1", "Q", "", mean - 15, mean + 15);
131  g1->GetParameters(par);
132 
133  B2INFO("Gaus fitting for each board");
134  for (int ib = 1; ib < 300; ++ib) {
135  if (m_hT0b[ib]->GetEntries() < 10) continue;
136  mean = m_hT0b[ib]->GetMean();
137  m_hT0b[ib]->SetDirectory(0);
138  m_hT0b[ib]->Fit("g1", "Q", "", mean - 15, mean + 15);
139  g1->GetParameters(par);
140  b.push_back(ib);
141  db.push_back(0.0);
142  Sb.push_back(par[1]);
143  dSb.push_back(g1->GetParError(1));
144  dtb[ib] = par[1];
145  err_dtb[ib] = g1->GetParError(1);
146  }
147  B2INFO("Gaus fitting for each cell");
148  for (int ilay = 0; ilay < 56; ++ilay) {
149  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
150  const int n = m_h1[ilay][iwire]->GetEntries();
151  B2DEBUG(99, "layer " << ilay << " wire " << iwire << " entries " << n);
152  if (n < 10) continue;
153  mean = m_h1[ilay][iwire]->GetMean();
154  m_h1[ilay][iwire]->SetDirectory(0);
155  m_h1[ilay][iwire]->Fit("g1", "Q", "", mean - 15, mean + 15);
156  g1->GetParameters(par);
157  c[ilay].push_back(iwire);
158  dc[ilay].push_back(0.0);
159  s[ilay].push_back(par[1]);
160  ds[ilay].push_back(g1->GetParError(1));
161 
162  dt[ilay][iwire] = par[1];
163  err_dt[ilay][iwire] = g1->GetParError(1);
164  }
165  }
166 
167 
168  if (m_storeHisto) {
169  B2INFO("Store histo");
170  TFile* fout = new TFile("Correct_T0.root", "RECREATE");
171  fout->cd();
172  TGraphErrors* gr[56];
173  TDirectory* top = gDirectory;
174  m_hTotal->Write();
175  TDirectory* subDir[56];
176  for (int ilay = 0; ilay < 56; ++ilay) {
177  subDir[ilay] = top ->mkdir(Form("lay_%d", ilay));
178  subDir[ilay]->cd();
179  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
180  m_h1[ilay][iwire]->Write();
181  }
182  }
183  top->cd();
184  TDirectory* corrT0 = top->mkdir("DeltaT0");
185  corrT0->cd();
186 
187  if (b.size() > 2) {
188  TGraphErrors* grb = new TGraphErrors(b.size(), &b.at(0), &Sb.at(0), &db.at(0), &dSb.at(0));
189  grb->SetMarkerColor(2);
190  grb->SetMarkerSize(1.0);
191  grb->SetTitle("#DeltaT0;BoardID;#DeltaT0[ns]");
192  grb->SetMaximum(10);
193  grb->SetMinimum(-10);
194  grb->SetName("Board");
195  grb->Write();
196  }
197  for (int sl = 0; sl < 56; ++sl) {
198  if (c[sl].size() < 2) continue;
199  gr[sl] = new TGraphErrors(c[sl].size(), &c[sl].at(0), &s[sl].at(0), &dc[sl].at(0), &ds[sl].at(0));
200  gr[sl]->SetMarkerColor(2);
201  gr[sl]->SetMarkerSize(1.0);
202  gr[sl]->SetTitle(Form("Layer_%d;IWire;#LT t_{mea}-t_{fit} #GT [ns]", sl));
203  gr[sl]->SetMaximum(10);
204  gr[sl]->SetMinimum(-10);
205  gr[sl]->SetName(Form("lay%d", sl));
206  gr[sl]->Write();
207  }
208  fout->Close();
209  }
210  B2INFO("Write constants");
211  Write();
212  return true;
213 }
214 void T0Correction::Write()
215 {
216  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
217  ofstream ofs(m_outputT0FileName.c_str());
219  tz.construct();
220 
221  double T0;
222  TH1F* T0B[300];
223  for (int ib = 0; ib < 300; ++ib) {
224  T0B[ib] = new TH1F(Form("T0B%d", ib), Form("boardID_%d", ib), 9000, 0, 9000);
225  }
226  //for calculate T0 mean of each board
227  for (int ilay = 0; ilay < 56; ++ilay) {
228  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
229  WireID wireid(ilay, iwire);
230  T0 = cdcgeo.getT0(wireid);
231  T0B[cdcgeo.getBoardID(wireid)]->Fill(T0);
232  }
233  }
234 
235  //correct T0 and write
236  for (int ilay = 0; ilay < 56; ++ilay) {
237  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
238  WireID wireid(ilay, iwire);
239  int bID = cdcgeo.getBoardID(wireid);
240  if (abs(err_dt[ilay][iwire]) > 2 || abs(dt[ilay][iwire]) > 1.e3) {
241  // if (abs(err_dt[ilay][iwire]) > 2) {
242  T0 = T0B[bID]->GetMean();
243  dt[ilay][iwire] = dtb[bID];
244  } else {
245  T0 = cdcgeo.getT0(wireid);
246  }
247  ofs << ilay << "\t" << iwire << "\t" << T0 - dt[ilay][iwire] << std::endl;
248  if (m_useDB) {
249  tz->setT0(wireid, T0 - dt[ilay][iwire]);
250  }
251  }
252  }
253  ofs.close();
254  if (m_useDB) {
257  tz.import(iov);
258  B2RESULT("T0 was calibrated and imported to database.");
259  }
260 
261 }
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::CDC::T0Correction::m_xmin
double m_xmin
minimum drift length
Definition: T0Correction.h:68
Belle2::CDC::T0Correction::calibrate
virtual bool calibrate()
Run algo on data.
Belle2::CDC::T0Correction::dtb
double dtb[300]
dt of each board
Definition: T0Correction.h:74
Belle2::CDC::T0Correction::m_firstRun
int m_firstRun
First run.
Definition: T0Correction.h:83
Belle2::CDCTimeZeros::setT0
void setT0(unsigned short iCLayer, unsigned short iWire, float t0)
Set t0 in the list.
Definition: CDCTimeZeros.h:50
Belle2::CDC::T0Correction::m_firstExperiment
int m_firstExperiment
First experiment.
Definition: T0Correction.h:82
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::CDC::T0Correction::m_hTotal
TH1F * m_hTotal
1D histogram of delta T whole channel
Definition: T0Correction.h:65
Belle2::CDC::T0Correction::m_Pvalmin
double m_Pvalmin
minimum pvalue required
Definition: T0Correction.h:70
Belle2::CDC::T0Correction::dt
double dt[56][385]
dt of each channel
Definition: T0Correction.h:72
Belle2::CDC::T0Correction::m_inputRootFileName
std::string m_inputRootFileName
input file names
Definition: T0Correction.h:81
Belle2::CDC::T0Correction::err_dt
double err_dt[56][385]
error of dt of each channel
Definition: T0Correction.h:73
Belle2::CDC::T0Correction::Write
virtual void Write()
write outut or store db
Belle2::CDC::T0Correction::m_ndfmin
double m_ndfmin
minimum ndf required
Definition: T0Correction.h:69
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:75
Belle2::CDC::T0Correction::m_outputT0FileName
std::string m_outputT0FileName
output t0 file name for text file
Definition: T0Correction.h:80
Belle2::CDC::T0Correction::m_useDB
bool m_useDB
use DB or text mode
Definition: T0Correction.h:79
Belle2::CDC::CDCGeometryPar::senseWireR
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
Definition: CDCGeometryPar.h:1231
Belle2::CDC::T0Correction::err_dtb
double err_dtb[300]
error of dt of board
Definition: T0Correction.h:75
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::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::CDC::T0Correction::m_h1
TH1F * m_h1[56][385]
1D histogram for each channel
Definition: T0Correction.h:66
Belle2::CDC::T0Correction::m_lastRun
int m_lastRun
Last run.
Definition: T0Correction.h:85
Belle2::CDC::T0Correction::m_storeHisto
bool m_storeHisto
store histo or not
Definition: T0Correction.h:78
Belle2::CDC::T0Correction::m_lastExperiment
int m_lastExperiment
Last experiment.
Definition: T0Correction.h:84
Belle2::CDC::T0Correction::CreateHisto
virtual void CreateHisto()
create histo for each channel
Belle2::CDC::T0Correction::m_hT0b
TH1F * m_hT0b[300]
1D histogram for each board
Definition: T0Correction.h:67