Belle II Software  release-05-01-25
CrudeT0CalibrationAlgorithm.cc
1 #include <cdc/calibration/CrudeT0CalibrationAlgorithm.h>
2 #include <calibration/CalibrationAlgorithm.h>
3 #include <cdc/geometry/CDCGeometryPar.h>
4 #include <cdc/dbobjects/CDCTimeZeros.h>
5 #include <framework/gearbox/Const.h>
6 #include <TError.h>
7 #include <TROOT.h>
8 #include <TF1.h>
9 #include <TFile.h>
10 #include <TTree.h>
11 
12 #include <framework/datastore/StoreObjPtr.h>
13 #include <framework/database/Database.h>
14 #include <framework/database/DBObjPtr.h>
15 #include <framework/database/IntervalOfValidity.h>
16 #include <framework/database/DBImportObjPtr.h>
17 #include <framework/logging/Logger.h>
18 using namespace std;
19 using namespace Belle2;
20 using namespace CDC;
21 
22 CrudeT0CalibrationAlgorithm::CrudeT0CalibrationAlgorithm(): CalibrationAlgorithm("CDCCrudeT0Collector")
23 {
24 
25  setDescription(
26  " -------------------------- T0 Calibration Algorithm -------------------------\n"
27  );
28 }
29 
30 void CrudeT0CalibrationAlgorithm::createHisto(StoreObjPtr<EventMetaData>& evtPtr)
31 
32 {
33 
34  B2INFO("CreateHisto");
35 
36  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
37  const auto exprun = getRunList();
38  B2INFO("Changed ExpRun to: " << exprun[0].first << " " << exprun[0].second);
39  evtPtr->setExperiment(exprun[0].first);
40  evtPtr->setRun(exprun[0].second);
41  DBStore::Instance().update();
42  B2INFO("create TDChist");
43 
44  for (int il = 0; il < 56; ++il) {
45  for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
46  m_hTDC[il][w] = new TH1D(Form("hLay%d_ch%d", il, w), "tdc", m_tdcMax - m_tdcMin, m_tdcMin, m_tdcMax);
47  }
48  }
49 
50  B2INFO("create TDChist(board)");
51  for (int ib = 0; ib < 300; ++ib) {
52  m_hTDCBoard[ib] = new TH1D(Form("hTDCBoard%d", ib), "", m_tdcMax - m_tdcMin, m_tdcMin, m_tdcMax);
53  }
54  m_hT0All = new TH1D("hT0All", "", 8500, 0, 8500);
55 
56  unsigned short lay;
57  unsigned short wire;
58  unsigned short tdc;
59  auto tree = getObjectPtr<TTree>("tree");
60  tree->SetBranchAddress("lay", &lay);
61  tree->SetBranchAddress("wire", &wire);
62  tree->SetBranchAddress("tdc", &tdc);
63 
64  const int nEntries = tree->GetEntries();
65  B2INFO("fill hist");
66  for (int i = 0; i < nEntries ; ++i) {
67  tree->GetEntry(i);
68  m_hTDC[lay][wire]->Fill(tdc);
69  m_hTDCBoard[cdcgeo.getBoardID(WireID(lay, wire))]->Fill(tdc);
70  m_hT0All->Fill(tdc);
71  }
72  B2INFO("end of filling hist");
73 }
74 
75 CalibrationAlgorithm::EResult CrudeT0CalibrationAlgorithm::calibrate()
76 {
77  B2INFO("Start calibration");
78 
79  gROOT->SetBatch(1);
80  gErrorIgnoreLevel = 3001;
81 
82  // We create an EventMetaData object. But since it's possible we're re-running this algorithm inside a process
83  // that has already created a DataStore, we need to check if it's already valid, or if it needs registering.
84  StoreObjPtr<EventMetaData> evtPtr;
85  if (!evtPtr.isValid()) {
86  // Construct an EventMetaData object in the Datastore so that the DB objects in CDCGeometryPar can work
87  DataStore::Instance().setInitializeActive(true);
88  B2INFO("Registering EventMetaData object in DataStore");
89  evtPtr.registerInDataStore();
90  DataStore::Instance().setInitializeActive(false);
91  B2INFO("Creating EventMetaData object");
92  evtPtr.create();
93  } else {
94  B2INFO("A valid EventMetaData object already exists.");
95  }
96  // Construct a CDCGeometryPar object which will update to the correct DB values when we change the EventMetaData and update
97  // the Database instance
98  DBObjPtr<CDCGeometry> cdcGeometry;
99  CDC::CDCGeometryPar::Instance(&(*cdcGeometry));
100  B2INFO("ExpRun at init : " << evtPtr->getExperiment() << " " << evtPtr->getRun());
101 
102  createHisto(evtPtr);
103 
104  TH1D* hs = new TH1D("hs", "sigma", 100, 0, 20);
105  std::vector<double> sb;
106  std::vector<double> dsb;
107  std::vector<double> t0b;
108  std::vector<double> dt0b;
109  std::vector<double> b;
110  std::vector<double> db;
111 
112  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
113  TF1* f1 = new TF1("f1", "[0]+[1]*(exp([2]*(x-[3]))/(1+exp(-([4]-x)/[5])))", m_tdcMin, m_tdcMax);
114  f1->SetParLimits(0, 0., 1000.);
115  f1->SetLineColor(kRed);
116  const double tdcBinWidth = cdcgeo.getTdcBinWidth();
117  bool bflag[300] = {false};
118 
119  for (int ib = 1; ib < 300; ++ib) {
120  if (m_hTDCBoard[ib]->GetEntries() < m_minEntries) {
121  B2DEBUG(199, "Warning: this board low statistic: " << m_hTDCBoard[ib]->GetEntries());
122  bflag[ib] = false;
123  m_t0b[ib] = m_initT0;
124  continue;
125  }
126  double p3 = m_hTDCBoard[ib]->GetXaxis()->GetBinCenter(m_hTDCBoard[ib]->GetMaximumBin());
127  f1->SetParameters(0, m_hTDCBoard[ib]->GetMaximum(), -0.001, p3, m_initT0, 2.5);
128  m_hTDCBoard[ib]->Fit("f1", "QM", "", m_initT0 - 60, m_initT0 + 60);
129 
130  if ((fabs(f1->GetParameter(4) - m_initT0) > 100)
131  || (fabs(f1->GetParameter(5)) < 0.01)
132  || (fabs(f1->GetParameter(5)) > 16)) {
133 
134  bflag[ib] = false;
135  m_t0b[ib] = m_initT0;
136  continue;
137  }
138 
139  bflag[ib] = true;
140  m_t0b[ib] = f1->GetParameter(4) * tdcBinWidth;
141 
142  sb.push_back(f1->GetParameter(5));
143  dsb.push_back(f1->GetParError(5));
144  t0b.push_back(f1->GetParameter(4));
145  dt0b.push_back(f1->GetParError(4));
146  b.push_back(ib);
147  db.push_back(0);
148  }
149 
150  for (int il = 0; il < 56; ++il) {
151  for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
152  B2DEBUG(99, "fitting for channel: " << il << " - " << w);
153  B2DEBUG(99, "number of entries" << m_hTDC[il][w]->GetEntries());
154  m_t0[il][w] = m_initT0 * tdcBinWidth;
155  int bid = cdcgeo.getBoardID(WireID(il, w));
156  if (m_hTDC[il][w]->GetEntries() < m_minEntries) {
157  B2DEBUG(99, "Warning: low statistic channel: " << m_hTDC[il][w]->GetEntries());
158  if (bflag[bid] != false) {
159  m_t0[il][w] = m_t0b[bid];
160  m_flag[il][w] = true;
161  } else {
162  m_flag[il][w] = false;
163  }
164  } else {
165  double p3 = m_hTDC[il][w]->GetXaxis()->GetBinCenter(m_hTDC[il][w]->GetMaximumBin());
166  f1->SetParameters(0, m_hTDC[il][w]->GetMaximum(), -0.001, p3, m_initT0, 2.5);
167  m_hTDC[il][w]->Fit("f1", "QM", "", m_initT0 - 60, m_initT0 + 60);
168  B2DEBUG(99, "prob of fit : " << f1->GetProb());
169  if ((f1->GetProb() < 1E-150) || (fabs(f1->GetParameter(4) - m_initT0) > 100) || (f1->GetParameter(5) < 0.1)
170  || (f1->GetParameter(5) > 20)) {
171  if (bflag[bid] != 0) {
172  m_t0[il][w] = m_t0b[bid];
173  m_flag[il][w] = true;
174  } else {
175  m_flag[il][w] = false;
176  }
177  } else {
178  m_t0[il][w] = f1->GetParameter(4) * tdcBinWidth;
179  hs->Fill(f1->GetParameter(5));
180  m_flag[il][w] = true;
181  }
182  }
183 
184  B2DEBUG(99, "P4 = " << m_t0[il][w]);
185  if (m_cosmic && cdcgeo.wireBackwardPosition(il, w).Y() > 0) {
186  m_t0[il][w] -= cdcgeo.senseWireR(il) / Const::speedOfLight;
187  } else {
188  m_t0[il][w] += cdcgeo.senseWireR(il) / Const::speedOfLight;
189  }
190  m_t0[il][w] += (m_zOffset - cdcgeo.wireBackwardPosition(il, w).Z()) / 27.25;
191  m_t0[il][w] += 6.122;
192  m_hT0All->Fill(m_t0[il][w]);
193  }
194  }
195 
196  B2INFO("Write constants");
197  write(evtPtr);
198  saveHisto();
199  return c_OK;
200 }
201 
202 void CrudeT0CalibrationAlgorithm::write(StoreObjPtr<EventMetaData>& evtPtr)
203 {
204  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
205  const auto exprun = getRunList();
206  B2INFO("Changed ExpRun to: " << exprun[0].first << " " << exprun[0].second);
207  evtPtr->setExperiment(exprun[0].first);
208  evtPtr->setRun(exprun[0].second);
209  DBStore::Instance().update();
210 
211  CDCTimeZeros* tz = new CDCTimeZeros();
212  for (int ilay = 0; ilay < 56; ++ilay) {
213  for (unsigned int iwire = 0; iwire < cdcgeo.nWiresInLayer(ilay); ++iwire) {
214  WireID wireid(ilay, iwire);
215  tz->setT0(wireid, m_t0[ilay][iwire]);
216  }
217  }
218  saveCalibration(tz, "CDCTimeZeros");
219 }
220 
221 void CrudeT0CalibrationAlgorithm::saveHisto()
222 {
223  static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
224  TFile* fhist = new TFile("histCrudeT0.root", "recreate");
225  fhist->cd();
226  TDirectory* top = gDirectory;
227  TDirectory* Direct[56];
228  for (int il = 0; il < 56; ++il) {
229  top->cd();
230  Direct[il] = gDirectory->mkdir(Form("lay_%d", il));
231  Direct[il]->cd();
232  for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
233  if (m_flag[il][w] == 1) {
234  m_hTDC[il][w]->Write();
235  }
236  }
237  }
238  top->cd();
239  TDirectory* board = gDirectory->mkdir("board");
240  board->cd();
241  for (int ib = 0; ib < 300; ++ib) {
242  if (m_hTDCBoard[ib]) {
243  m_hTDCBoard[ib]->Write();
244  }
245  }
246  top->cd();
247  m_hT0All->Write();
248 
249  fhist->Close();
250 }
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:44
Belle2::CDC::CrudeT0CalibrationAlgorithm::createHisto
virtual void createHisto(StoreObjPtr< EventMetaData > &evtPtr)
create histo for each channel
Definition: CrudeT0CalibrationAlgorithm.cc:30
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_cosmic
bool m_cosmic
for cosmic case, tof of upper sector will be negative
Definition: CrudeT0CalibrationAlgorithm.h:106
Belle2::DataStore::Instance
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_hTDC
TH1D * m_hTDC[56][400]
TDC distribution histo.
Definition: CrudeT0CalibrationAlgorithm.h:98
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::DataStore::setInitializeActive
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_minEntries
unsigned short m_minEntries
minimum entries required by histo.
Definition: CrudeT0CalibrationAlgorithm.h:104
Belle2::CDCTimeZeros::setT0
void setT0(unsigned short iCLayer, unsigned short iWire, float t0)
Set t0 in the list.
Definition: CDCTimeZeros.h:50
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_initT0
float m_initT0
Common initial T0 for fitting.
Definition: CrudeT0CalibrationAlgorithm.h:103
Belle2::CDC::CrudeT0CalibrationAlgorithm::calibrate
EResult calibrate() override
Run algo on data.
Definition: CrudeT0CalibrationAlgorithm.cc:75
Belle2::CDC::CDCGeometryPar::getTdcBinWidth
double getTdcBinWidth() const
Return TDC bin width (nsec).
Definition: CDCGeometryPar.h:732
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::CrudeT0CalibrationAlgorithm::m_zOffset
float m_zOffset
z offset for calculate prop time, it is position of trigger counter,
Definition: CrudeT0CalibrationAlgorithm.h:105
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
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::CDCGeometryPar::wireBackwardPosition
const TVector3 wireBackwardPosition(int layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
Definition: CDCGeometryPar.cc:1662
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_tdcMax
unsigned short m_tdcMax
maximum of TDC hist for fitting
Definition: CrudeT0CalibrationAlgorithm.h:102
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::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
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::DBStore::Instance
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:36
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_t0
double m_t0[56][400]
T0 of each channel.
Definition: CrudeT0CalibrationAlgorithm.h:95
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_t0b
double m_t0b[300]
T0 for each board.
Definition: CrudeT0CalibrationAlgorithm.h:94
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_tdcMin
unsigned short m_tdcMin
minimum of TDC hist for fitting
Definition: CrudeT0CalibrationAlgorithm.h:101
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::CrudeT0CalibrationAlgorithm::m_flag
bool m_flag[56][400]
flag =1 for good, =0 for low statistic or bad fit
Definition: CrudeT0CalibrationAlgorithm.h:96
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_hT0All
TH1D * m_hT0All
T0 distribution of all channel.
Definition: CrudeT0CalibrationAlgorithm.h:100
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1211
Belle2::CDC::CrudeT0CalibrationAlgorithm::m_hTDCBoard
TH1D * m_hTDCBoard[300]
T0 distribution of each board.
Definition: CrudeT0CalibrationAlgorithm.h:99
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::CDC::CrudeT0CalibrationAlgorithm::write
virtual void write(StoreObjPtr< EventMetaData > &evtPtr)
write outut or store db
Definition: CrudeT0CalibrationAlgorithm.cc:202
Belle2::DBStore::update
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:87
Belle2::CDC::CrudeT0CalibrationAlgorithm::saveHisto
void saveHisto()
Save hitograms of the calibration results.
Definition: CrudeT0CalibrationAlgorithm.cc:221
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:120