Belle II Software  release-08-01-10
KLMTimeAlgorithm.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 /* Own header. */
10 #include <klm/calibration/KLMTimeAlgorithm.h>
11 
12 /* KLM headers. */
13 #include <klm/dataobjects/bklm/BKLMElementNumbers.h>
14 
15 /* Basf2 headers. */
16 #include <framework/database/Database.h>
17 #include <framework/database/DBObjPtr.h>
18 #include <framework/database/DBStore.h>
19 #include <framework/dataobjects/EventMetaData.h>
20 #include <framework/gearbox/Const.h>
21 #include <framework/logging/Logger.h>
22 
23 /* ROOT headers. */
24 #include <Math/MinimizerOptions.h>
25 #include <TFile.h>
26 #include <TFitResult.h>
27 #include <TMinuit.h>
28 #include <TROOT.h>
29 #include <TString.h>
30 #include <TTree.h>
31 #include <TVector3.h>
32 
33 using namespace Belle2;
34 using namespace ROOT::Math;
35 
37 const int c_NBinsTime = 100;
38 
40 const int c_NBinsDistance = 100;
41 
43 static double s_BinnedData[c_NBinsTime][c_NBinsDistance];
44 
46 static double s_LowerTimeBoundary = 0;
47 
49 static double s_UpperTimeBoundary = 0;
50 
52 static double s_StripLength = 0;
53 
54 static bool compareEventNumber(const std::pair<KLMChannelNumber, unsigned int>& pair1,
55  const std::pair<KLMChannelNumber, unsigned int>& pair2)
56 {
57  return pair1.second < pair2.second;
58 }
59 
60 static double timeDensity(const double x[2], const double* par)
61 {
62  double polynomial, t0, gauss;
63  polynomial = par[0];
64  t0 = par[2] + par[4] * x[1];
65  gauss = par[1] / (sqrt(2.0 * M_PI) * par[3]) *
66  exp(-0.5 * pow((x[0] - t0) / par[3], 2));
67  return fabs(polynomial + gauss);
68 }
69 
70 static void fcn(int& npar, double* grad, double& fval, double* par, int iflag)
71 {
72  (void)npar;
73  (void)grad;
74  (void)iflag;
75  double x[2];
76  fval = 0;
77  for (int i = 0; i < c_NBinsTime; ++i) {
78  x[0] = s_LowerTimeBoundary +
79  (s_UpperTimeBoundary - s_LowerTimeBoundary) *
80  (double(i) + 0.5) / c_NBinsTime;
81  for (int j = 0; j < c_NBinsDistance; ++j) {
82  x[1] = s_StripLength * (double(j) + 0.5) / c_NBinsDistance;
83  double f = timeDensity(x, par);
84  if (s_BinnedData[i][j] == 0)
85  fval = fval + 2.0 * f;
86  else
87  fval = fval + 2.0 * (f - s_BinnedData[i][j] *
88  (1.0 - log(s_BinnedData[i][j] / f)));
89  }
90  }
91 }
92 
93 KLMTimeAlgorithm::KLMTimeAlgorithm() :
94  CalibrationAlgorithm("KLMTimeCollector")
95 {
96  m_ElementNumbers = &(KLMElementNumbers::Instance());
97  m_minimizerOptions = ROOT::Math::MinimizerOptions();
98 }
99 
100 KLMTimeAlgorithm::~KLMTimeAlgorithm()
101 {
102 }
103 
104 void KLMTimeAlgorithm::setupDatabase()
105 {
106  const std::vector<Calibration::ExpRun>& runs = getRunList();
107  int firstExperiment = runs[0].first;
108  int lastExperiment = runs[runs.size() - 1].first;
109  if (firstExperiment != lastExperiment) {
110  B2FATAL("Runs from different experiments are used "
111  "for KLM time calibration (single algorithm run).");
112  }
113  /* DataStore. */
114  DataStore::Instance().setInitializeActive(true);
115  StoreObjPtr<EventMetaData> eventMetaData;
116  eventMetaData.registerInDataStore();
117  DataStore::Instance().setInitializeActive(false);
118  /* Database. */
119  if (eventMetaData.isValid()) {
120  if (eventMetaData->getExperiment() != firstExperiment) {
121  B2FATAL("Runs from different experiments are used "
122  "for KLM time calibration (consecutive algorithm runs).");
123  }
124  eventMetaData->setExperiment(firstExperiment);
125  eventMetaData->setRun(runs[0].second);
126  } else {
127  eventMetaData.construct(1, runs[0].second, firstExperiment);
128  }
129  DBStore& dbStore = DBStore::Instance();
130  dbStore.update();
131  dbStore.updateEvent();
132  /*
133  * For calibration algorithms, the database is not initialized on class
134  * creation. Do not move the database object to class members.
135  */
136  DBObjPtr<BKLMGeometryPar> bklmGeometry;
137  m_BKLMGeometry = bklm::GeometryPar::instance(*bklmGeometry);
138  m_EKLMGeometry = &(EKLM::GeometryData::Instance());
139 }
140 
141 CalibrationAlgorithm::EResult KLMTimeAlgorithm::readCalibrationData()
142 {
143  B2INFO("Read tree entries and seprate events by module id.");
144  Event event;
145  std::shared_ptr<TTree> timeCalibrationData;
146  timeCalibrationData = getObjectPtr<TTree>("time_calibration_data");
147  timeCalibrationData->SetBranchAddress("t0", &event.t0);
148  timeCalibrationData->SetBranchAddress("flyTime", &event.flyTime);
149  timeCalibrationData->SetBranchAddress("recTime", &event.recTime);
150  timeCalibrationData->SetBranchAddress("dist", &event.dist);
151  timeCalibrationData->SetBranchAddress("diffDistX", &event.diffDistX);
152  timeCalibrationData->SetBranchAddress("diffDistY", &event.diffDistY);
153  timeCalibrationData->SetBranchAddress("diffDistZ", &event.diffDistZ);
154  timeCalibrationData->SetBranchAddress("eDep", &event.eDep);
155  timeCalibrationData->SetBranchAddress("nPE", &event.nPE);
156  timeCalibrationData->SetBranchAddress("channelId", &event.channelId);
157  timeCalibrationData->SetBranchAddress("inRPC", &event.inRPC);
158  timeCalibrationData->SetBranchAddress("isFlipped", &event.isFlipped);
159 
160  B2INFO(LogVar("Total number of digits:", timeCalibrationData->GetEntries()));
161  m_evts.clear();
162 
163  int n = timeCalibrationData->GetEntries();
164  if (n < m_MinimalDigitNumber)
165  return CalibrationAlgorithm::c_NotEnoughData;
166  for (int i = 0; i < n; ++i) {
167  timeCalibrationData->GetEntry(i);
168  m_evts[event.channelId].push_back(event);
169  }
170  B2INFO("Events packing finish.");
171  return CalibrationAlgorithm::c_OK;
172 }
173 
174 void KLMTimeAlgorithm::createHistograms()
175 {
176  if (m_mc) {
177  m_LowerTimeBoundaryRPC = -10.0;
178  m_UpperTimeBoundaryRPC = 10.0;
179  m_LowerTimeBoundaryScintilltorsBKLM = 20.0;
180  m_UpperTimeBoundaryScintilltorsBKLM = 70.0;
181  m_LowerTimeBoundaryScintilltorsEKLM = 20.0;
182  m_UpperTimeBoundaryScintilltorsEKLM = 70.0;
183  } else {
184  m_LowerTimeBoundaryRPC = -800.0;
185  m_UpperTimeBoundaryRPC = -600.0;
186  m_LowerTimeBoundaryScintilltorsBKLM = -4800.0;
187  m_UpperTimeBoundaryScintilltorsBKLM = -4400.0;
188  m_LowerTimeBoundaryScintilltorsEKLM = -4950.0;
189  m_UpperTimeBoundaryScintilltorsEKLM = -4650.0;
190 
191  }
192  int nBin = 200;
193  int nBin_scint = 400;
194 
195  TString iFstring[2] = {"Backward", "Forward"};
196  TString iPstring[2] = {"ZReadout", "PhiReadout"};
197  TString hn, ht;
198 
199  h_diff = new TH1F("h_diff", "Position difference between bklmHit2d and extHit;position difference", 100, 0, 10);
200  h_calibrated = new TH1I("h_calibrated_summary", "h_calibrated_summary;calibrated or not", 3, 0, 3);
201  hc_calibrated = new TH1I("hc_calibrated_summary", "hc_calibrated_summary;calibrated or not", 3, 0, 3);
202 
203  gre_time_channel_scint = new TGraphErrors();
204  gre_time_channel_rpc = new TGraphErrors();
205  gre_time_channel_scint_end = new TGraphErrors();
206 
207  gr_timeShift_channel_scint = new TGraph();
208  gr_timeShift_channel_rpc = new TGraph();
209  gr_timeShift_channel_scint_end = new TGraph();
210 
211  gre_ctime_channel_scint = new TGraphErrors();
212  gre_ctime_channel_rpc = new TGraphErrors();
213  gre_ctime_channel_scint_end = new TGraphErrors();
214 
215  gr_timeRes_channel_scint = new TGraph();
216  gr_timeRes_channel_rpc = new TGraph();
217  gr_timeRes_channel_scint_end = new TGraph();
218 
219  double maximalPhiStripLengthBKLM =
220  m_BKLMGeometry->getMaximalPhiStripLength();
221  double maximalZStripLengthBKLM =
222  m_BKLMGeometry->getMaximalZStripLength();
223  double maximalStripLengthEKLM =
224  m_EKLMGeometry->getMaximalStripLength() / CLHEP::cm * Unit::cm;
225 
226  m_ProfileRpcPhi = new TProfile("hprf_rpc_phi_effC",
227  "Time over propagation length for RPCs (Phi_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
228  400.0);
229  m_ProfileRpcZ = new TProfile("hprf_rpc_z_effC",
230  "Time over propagation length for RPCs (Z_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
231  400.0);
232  m_ProfileBKLMScintillatorPhi = new TProfile("hprf_scint_phi_effC",
233  "Time over propagation length for scintillators (Phi_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
234  200, 0.0, maximalPhiStripLengthBKLM);
235  m_ProfileBKLMScintillatorZ = new TProfile("hprf_scint_z_effC",
236  "Time over propagation length for scintillators (Z_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
237  200, 0.0, maximalZStripLengthBKLM);
238  m_ProfileEKLMScintillatorPlane1 = new TProfile("hprf_scint_plane1_effC_end",
239  "Time over propagation length for scintillators (plane1, Endcap); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
240  200, 0.0, maximalStripLengthEKLM);
241  m_ProfileEKLMScintillatorPlane2 = new TProfile("hprf_scint_plane2_effC_end",
242  "Time over propagation length for scintillators (plane2, Endcap); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
243  200, 0.0, maximalStripLengthEKLM);
244 
245  m_Profile2RpcPhi = new TProfile("hprf2_rpc_phi_effC",
246  "Time over propagation length for RPCs (Phi_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
247  400.0);
248  m_Profile2RpcZ = new TProfile("hprf2_rpc_z_effC",
249  "Time over propagation length for RPCs (Z_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
250  400.0);
251  m_Profile2BKLMScintillatorPhi = new TProfile("hprf2_scint_phi_effC",
252  "Time over propagation length for scintillators (Phi_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
253  200, 0.0, maximalPhiStripLengthBKLM);
254  m_Profile2BKLMScintillatorZ = new TProfile("hprf2_scint_z_effC",
255  "Time over propagation length for scintillators (Z_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
256  200, 0.0, maximalZStripLengthBKLM);
257  m_Profile2EKLMScintillatorPlane1 = new TProfile("hprf2_scint_plane1_effC_end",
258  "Time over propagation length for scintillators (plane1, Endcap); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
259  200, 0.0, maximalStripLengthEKLM);
260  m_Profile2EKLMScintillatorPlane2 = new TProfile("hprf2_scint_plane2_effC_end",
261  "Time over propagation length for scintillators (plane2, Endcap); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]",
262  200, 0.0, maximalStripLengthEKLM);
263 
264  h_time_rpc_tc = new TH1F("h_time_rpc_tc", "time distribtution for RPC", nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
265  h_time_scint_tc = new TH1F("h_time_scint_tc", "time distribtution for Scintillator", nBin_scint,
266  m_LowerTimeBoundaryScintilltorsBKLM, m_UpperTimeBoundaryScintilltorsBKLM);
267  h_time_scint_tc_end = new TH1F("h_time_scint_tc_end", "time distribtution for Scintillator (Endcap)", nBin_scint,
268  m_LowerTimeBoundaryScintilltorsEKLM,
269  m_UpperTimeBoundaryScintilltorsEKLM);
270 
272  h_time_rpc = new TH1F("h_time_rpc", "time distribtution for RPC; T_rec-T_0-T_fly-T_propagation[ns]", nBin, m_LowerTimeBoundaryRPC,
273  m_UpperTimeBoundaryRPC);
274  h_time_scint = new TH1F("h_time_scint", "time distribtution for Scintillator; T_rec-T_0-T_fly-T_propagation[ns]", nBin_scint,
275  m_LowerTimeBoundaryScintilltorsBKLM, m_UpperTimeBoundaryScintilltorsBKLM);
276  h_time_scint_end = new TH1F("h_time_scint_end", "time distribtution for Scintillator (Endcap); T_rec-T_0-T_fly-T_propagation[ns]",
277  nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM, m_UpperTimeBoundaryScintilltorsEKLM);
278 
279  hc_time_rpc = new TH1F("hc_time_rpc", "Calibrated time distribtution for RPC; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
280  nBin, m_LowerTimeBoundaryCalibratedRPC, m_UpperTimeBoundaryCalibratedRPC);
281  hc_time_scint = new TH1F("hc_time_scint",
282  "Calibrated time distribtution for Scintillator; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", nBin_scint,
283  m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
284  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
285  hc_time_scint_end = new TH1F("hc_time_scint_end",
286  "Calibrated time distribtution for Scintillator (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", nBin_scint,
287  m_LowerTimeBoundaryCalibratedScintilltorsEKLM, m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
288 
289  for (int iF = 0; iF < 2; ++iF) {
290  hn = Form("h_timeF%d_rpc", iF);
291  ht = Form("Time distribtution for RPC of %s; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
292  h_timeF_rpc[iF] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
293  hn = Form("h_timeF%d_scint", iF);
294  ht = Form("Time distribtution for Scintillator of %s; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
295  h_timeF_scint[iF] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
296  m_UpperTimeBoundaryScintilltorsBKLM);
297  hn = Form("h_timeF%d_scint_end", iF);
298  ht = Form("Time distribtution for Scintillator of %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
299  h_timeF_scint_end[iF] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
300  m_UpperTimeBoundaryScintilltorsEKLM);
301 
302  hn = Form("h2_timeF%d_rpc", iF);
303  ht = Form("Time distribtution for RPC of %s; Sector Index; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
304  h2_timeF_rpc[iF] = new TH2F(hn.Data(), ht.Data(), 8, 0, 8, nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
305  hn = Form("h2_timeF%d_scint", iF);
306  ht = Form("Time distribtution for Scintillator of %s; Sector Index; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
307  h2_timeF_scint[iF] = new TH2F(hn.Data(), ht.Data(), 8, 0, 8, nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
308  m_UpperTimeBoundaryScintilltorsBKLM);
309  hn = Form("h2_timeF%d_scint_end", iF);
310  ht = Form("Time distribtution for Scintillator of %s (Endcap); Sector Index; T_rec-T_0-T_fly-T_propagation[ns]",
311  iFstring[iF].Data());
312  h2_timeF_scint_end[iF] = new TH2F(hn.Data(), ht.Data(), 4, 0, 4, nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
313  m_UpperTimeBoundaryScintilltorsEKLM);
314 
315  hn = Form("hc_timeF%d_rpc", iF);
316  ht = Form("Calibrated time distribtution for RPC of %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iFstring[iF].Data());
317  hc_timeF_rpc[iF] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryCalibratedRPC, m_UpperTimeBoundaryCalibratedRPC);
318  hn = Form("hc_timeF%d_scint", iF);
319  ht = Form("Calibrated time distribtution for Scintillator of %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
320  iFstring[iF].Data());
321  hc_timeF_scint[iF] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
322  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
323  hn = Form("hc_timeF%d_scint_end", iF);
324  ht = Form("Calibrated time distribtution for Scintillator of %s (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
325  iFstring[iF].Data());
326  hc_timeF_scint_end[iF] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
327  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
328 
329  hn = Form("h2c_timeF%d_rpc", iF);
330  ht = Form("Calibrated time distribtution for RPC of %s; Sector Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
331  iFstring[iF].Data());
332  h2c_timeF_rpc[iF] = new TH2F(hn.Data(), ht.Data(), 8, 0, 8, nBin, m_LowerTimeBoundaryCalibratedRPC,
333  m_UpperTimeBoundaryCalibratedRPC);
334  hn = Form("h2c_timeF%d_scint", iF);
335  ht = Form("Calibrated time distribtution for Scintillator of %s; Sector Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
336  iFstring[iF].Data());
337  h2c_timeF_scint[iF] = new TH2F(hn.Data(), ht.Data(), 8, 0, 8, nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
338  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
339  hn = Form("h2c_timeF%d_scint_end", iF);
340  ht = Form("Calibrated time distribtution for Scintillator of %s (Endcap) ; Sector Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
341  iFstring[iF].Data());
342  h2c_timeF_scint_end[iF] = new TH2F(hn.Data(), ht.Data(), 4, 0, 4, nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
343  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
344 
345  for (int iS = 0; iS < 8; ++iS) {
346  // Barrel parts
347  hn = Form("h_timeF%d_S%d_scint", iF, iS);
348  ht = Form("Time distribtution for Scintillator of Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
349  h_timeFS_scint[iF][iS] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
350  m_UpperTimeBoundaryScintilltorsBKLM);
351  hn = Form("h_timeF%d_S%d_rpc", iF, iS);
352  ht = Form("Time distribtution for RPC of Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
353  h_timeFS_rpc[iF][iS] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
354  hn = Form("h2_timeF%d_S%d", iF, iS);
355  ht = Form("Time distribtution of Sector%d, %s; Layer Index; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
356  h2_timeFS[iF][iS] = new TH2F(hn.Data(), ht.Data(), 15, 0, 15, nBin_scint, m_LowerTimeBoundaryRPC,
357  m_UpperTimeBoundaryScintilltorsBKLM);
358 
359  hn = Form("hc_timeF%d_S%d_scint", iF, iS);
360  ht = Form("Calibrated time distribtution for Scintillator of Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iS,
361  iFstring[iF].Data());
362  hc_timeFS_scint[iF][iS] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
363  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
364  hn = Form("hc_timeF%d_S%d_rpc", iF, iS);
365  ht = Form("Calibrated time distribtution for RPC of Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iS,
366  iFstring[iF].Data());
367  hc_timeFS_rpc[iF][iS] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedRPC,
368  m_UpperTimeBoundaryCalibratedRPC);
369  hn = Form("h2c_timeF%d_S%d", iF, iS);
370  ht = Form("Calibrated time distribtution of Sector%d, %s; Layer Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iS,
371  iFstring[iF].Data());
372  h2c_timeFS[iF][iS] = new TH2F(hn.Data(), ht.Data(), 15, 0, 15, nBin_scint, m_LowerTimeBoundaryCalibratedRPC,
373  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
374 
375  // Inner 2 layers --> Scintillators
376  for (int iL = 0; iL < 2; ++iL) {
377  hn = Form("h_timeF%d_S%d_L%d", iF, iS, iL);
378  ht = Form("Time distribtution for Scintillator of Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iL, iS,
379  iFstring[iF].Data());
380  h_timeFSL[iF][iS][iL] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
381  m_UpperTimeBoundaryScintilltorsBKLM);
382  hn = Form("hc_timeF%d_S%d_L%d", iF, iS, iL);
383  ht = Form("Calibrated time distribtution for Scintillator of Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
384  iL, iS, iFstring[iF].Data());
385  hc_timeFSL[iF][iS][iL] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
386  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
387 
388  for (int iP = 0; iP < 2; ++iP) {
389  hn = Form("h_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
390  ht = Form("Time distribtution for Scintillator of %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]",
391  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
392  h_timeFSLP[iF][iS][iL][iP] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
393  m_UpperTimeBoundaryScintilltorsBKLM);
394  hn = Form("h2_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
395  ht = Form("Time distribtution for Scintillator of %s, Layer%d, Sector%d, %s; Channel Index; T_rec-T_0-T_fly-T_propagation[ns]",
396  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
397  h2_timeFSLP[iF][iS][iL][iP] = new TH2F(hn.Data(), ht.Data(), 54, 0, 54, nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
398  m_UpperTimeBoundaryScintilltorsBKLM);
399 
400  hn = Form("hc_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
401  ht = Form("Calibrated time distribtution for Scintillator of %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
402  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
403  hc_timeFSLP[iF][iS][iL][iP] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
404  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
405  hn = Form("h2c_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
406  ht = Form("Calibrated time distribtution for Scintillator of %s, Layer%d, Sector%d, %s; Channel Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
407  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
408  h2c_timeFSLP[iF][iS][iL][iP] = new TH2F(hn.Data(), ht.Data(), 54, 0, 54, nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
409  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
410 
411  int nchannel_max = BKLMElementNumbers::getNStrips(iF, iS + 1, iL + 1, iP);
412  for (int iC = 0; iC < nchannel_max; ++iC) {
413  hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_tc", iF, iS, iL, iP, iC);
414  ht = Form("time distribtution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iC,
415  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
416  h_timeFSLPC_tc[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
417  m_UpperTimeBoundaryScintilltorsBKLM);
418 
419  hn = Form("h_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
420  ht = Form("time distribtution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iC,
421  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
422  h_timeFSLPC[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsBKLM,
423  m_UpperTimeBoundaryScintilltorsBKLM);
424 
425  hn = Form("hc_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
426  ht = Form("Calibrated time distribtution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
427  iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
428  hc_timeFSLPC[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
429  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
430  hn = Form("time_length_bklm_F%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
431  double stripLength = 200;
432  m_HistTimeLengthBKLM[iF][iS][iL][iP][iC] =
433  new TH2F(hn.Data(),
434  "Time versus propagation length; "
435  "propagation distance[cm]; "
436  "T_rec-T_0-T_fly-'T_calibration'[ns]",
437  200, 0.0, stripLength,
438  400, m_LowerTimeBoundaryCalibratedScintilltorsBKLM,
439  m_UpperTimeBoundaryCalibratedScintilltorsBKLM);
440  }
441  }
442  }
443 
444  for (int iL = 2; iL < 15; ++iL) {
445  hn = Form("h_timeF%d_S%d_L%d", iF, iS, iL);
446  ht = Form("time distribtution for RPC of Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iL, iS, iFstring[iF].Data());
447  h_timeFSL[iF][iS][iL] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
448 
449  hn = Form("hc_timeF%d_S%d_L%d", iF, iS, iL);
450  ht = Form("Calibrated time distribtution for RPC of Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iL, iS,
451  iFstring[iF].Data());
452  hc_timeFSL[iF][iS][iL] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryCalibratedRPC, m_UpperTimeBoundaryCalibratedRPC);
453 
454  for (int iP = 0; iP < 2; ++iP) {
455  hn = Form("h_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
456  ht = Form("time distribtution for RPC of %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iPstring[iP].Data(), iL, iS,
457  iFstring[iF].Data());
458  h_timeFSLP[iF][iS][iL][iP] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
459 
460  hn = Form("h2_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
461  ht = Form("time distribtution for RPC of %s, Layer%d, Sector%d, %s; Channel Index; T_rec-T_0-T_fly-T_propagation[ns]",
462  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
463  h2_timeFSLP[iF][iS][iL][iP] = new TH2F(hn.Data(), ht.Data(), 48, 0, 48, nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
464 
465  hn = Form("hc_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
466  ht = Form("Calibrated time distribtution for RPC of %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
467  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
468  hc_timeFSLP[iF][iS][iL][iP] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryCalibratedRPC,
469  m_UpperTimeBoundaryCalibratedRPC);
470 
471  hn = Form("h2c_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
472  ht = Form("Calibrated time distribtution for RPC of %s, Layer%d, Sector%d, %s; Channel Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
473  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
474  h2c_timeFSLP[iF][iS][iL][iP] = new TH2F(hn.Data(), ht.Data(), 48, 0, 48, nBin, m_LowerTimeBoundaryCalibratedRPC,
475  m_UpperTimeBoundaryCalibratedRPC);
476 
477  int nchannel_max = BKLMElementNumbers::getNStrips(iF, iS + 1, iL + 1, iP);
478  for (int iC = 0; iC < nchannel_max; ++iC) {
479  hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_tc", iF, iS, iL, iP, iC);
480  ht = Form("Time distribtution for RPC of Channel%d, %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iC,
481  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
482  h_timeFSLPC_tc[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
483 
484  hn = Form("h_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
485  ht = Form("Time distribtution for RPC of Channel%d, %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iC,
486  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
487  h_timeFSLPC[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryRPC, m_UpperTimeBoundaryRPC);
488 
489  hn = Form("hc_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
490  ht = Form("Calibrated time distribtution for RPC of Channel%d, %s, Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
491  iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
492  hc_timeFSLPC[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin, m_LowerTimeBoundaryCalibratedRPC,
493  m_UpperTimeBoundaryCalibratedRPC);
494  }
495  }
496  }
497  }
498  // Endcap part
499  int maxLay = 12 + 2 * iF;
500  for (int iS = 0; iS < 4; ++iS) {
501  hn = Form("h_timeF%d_S%d_scint_end", iF, iS);
502  ht = Form("Time distribtution for Scintillator of Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]", iS,
503  iFstring[iF].Data());
504  h_timeFS_scint_end[iF][iS] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
505  m_UpperTimeBoundaryScintilltorsEKLM);
506  hn = Form("h2_timeF%d_S%d_end", iF, iS);
507  ht = Form("Time distribtution of Sector%d, %s (Endcap); Layer Index; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
508  h2_timeFS_end[iF][iS] = new TH2F(hn.Data(), ht.Data(), maxLay, 0, maxLay, nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
509  m_UpperTimeBoundaryScintilltorsEKLM);
510  hn = Form("hc_timeF%d_S%d_scint_end", iF, iS);
511  ht = Form("Calibrated time distribtution for Scintillator of Sector%d (Endcap), %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
512  iS, iFstring[iF].Data());
513  hc_timeFS_scint_end[iF][iS] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
514  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
515  hn = Form("h2c_timeF%d_S%d_end", iF, iS);
516  ht = Form("Calibrated time distribtution of Sector%d, %s (Endcap); Layer Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
517  iS, iFstring[iF].Data());
518  h2c_timeFS_end[iF][iS] = new TH2F(hn.Data(), ht.Data(), maxLay, 0, maxLay, nBin_scint,
519  m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
520  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
521 
522  for (int iL = 0; iL < maxLay; ++iL) {
523  hn = Form("h_timeF%d_S%d_L%d_end", iF, iS, iL);
524  ht = Form("Time distribtution for Scintillator of Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]", iL, iS,
525  iFstring[iF].Data());
526  h_timeFSL_end[iF][iS][iL] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
527  m_UpperTimeBoundaryScintilltorsEKLM);
528  hn = Form("hc_timeF%d_S%d_L%d_end", iF, iS, iL);
529  ht = Form("Calibrated time distribtution for Scintillator of Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
530  iL, iS, iFstring[iF].Data());
531  hc_timeFSL_end[iF][iS][iL] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
532  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
533 
534  for (int iP = 0; iP < 2; ++iP) {
535  hn = Form("h_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
536  ht = Form("Time distribtution for Scintillator of %s, Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]",
537  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
538  h_timeFSLP_end[iF][iS][iL][iP] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
539  m_UpperTimeBoundaryScintilltorsEKLM);
540 
541  hn = Form("h2_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
542  ht = Form("Time distribtution for Scintillator of %s, Layer%d, Sector%d, %s (Endcap); Channel Index; T_rec-T_0-T_fly-T_propagation[ns]",
543  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
544  h2_timeFSLP_end[iF][iS][iL][iP] = new TH2F(hn.Data(), ht.Data(), 75, 0, 75, nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
545  m_UpperTimeBoundaryScintilltorsEKLM);
546 
547  hn = Form("hc_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
548  ht = Form("Calibrated time distribtution for Scintillator of %s, Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
549  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
550  hc_timeFSLP_end[iF][iS][iL][iP] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
551  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
552 
553  hn = Form("h2c_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
554  ht = Form("Calibrated time distribtution for Scintillator of %s, Layer%d, Sector%d, %s (Endcap); Channel Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
555  iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
556  h2c_timeFSLP_end[iF][iS][iL][iP] = new TH2F(hn.Data(), ht.Data(), 75, 0, 75, nBin_scint,
557  m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
558  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
559 
560  for (int iC = 0; iC < 75; ++iC) {
561  hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_tc_end", iF, iS, iL, iP, iC);
562  ht = Form("Time distribtution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]",
563  iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
564  h_timeFSLPC_tc_end[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
565  m_UpperTimeBoundaryScintilltorsEKLM);
566 
567  hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_end", iF, iS, iL, iP, iC);
568  ht = Form("Time distribtution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]",
569  iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
570  h_timeFSLPC_end[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryScintilltorsEKLM,
571  m_UpperTimeBoundaryScintilltorsEKLM);
572  hn = Form("hc_timeF%d_S%d_L%d_P%d_C%d_end", iF, iS, iL, iP, iC);
573  ht = Form("Calibrated time distribtution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
574  iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
575  hc_timeFSLPC_end[iF][iS][iL][iP][iC] = new TH1F(hn.Data(), ht.Data(), nBin_scint, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
576  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
577  hn = Form("time_length_eklm_F%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
578  double stripLength = m_EKLMGeometry->getStripLength(iC + 1) /
579  CLHEP::cm * Unit::cm;
580  m_HistTimeLengthEKLM[iF][iS][iL][iP][iC] =
581  new TH2F(hn.Data(),
582  "Time versus propagation length; "
583  "propagation distance[cm]; "
584  "T_rec-T_0-T_fly-'T_calibration'[ns]",
585  200, 0.0, stripLength,
586  400, m_LowerTimeBoundaryCalibratedScintilltorsEKLM,
587  m_UpperTimeBoundaryCalibratedScintilltorsEKLM);
588  }
589  }
590  }
591  }
592  }
593 }
594 
595 void KLMTimeAlgorithm::fillTimeDistanceProfiles(
596  TProfile* profileRpcPhi, TProfile* profileRpcZ,
597  TProfile* profileBKLMScintillatorPhi, TProfile* profileBKLMScintillatorZ,
598  TProfile* profileEKLMScintillatorPlane1,
599  TProfile* profileEKLMScintillatorPlane2, bool fill2dHistograms)
600 {
601  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
602  KLMChannelNumber channel = klmChannel.getKLMChannelNumber();
603  if (m_cFlag[channel] == ChannelCalibrationStatus::c_NotEnoughData)
604  continue;
605 
606  std::vector<struct Event>::iterator it;
607  std::vector<struct Event> eventsChannel;
608  eventsChannel = m_evts[channel];
609  int iSub = klmChannel.getSubdetector();
610 
611  for (it = eventsChannel.begin(); it != eventsChannel.end(); ++it) {
612  double timeHit = it->time() - m_timeShift[channel];
613  if (m_useEventT0)
614  timeHit = timeHit - it->t0;
615  double distHit = it->dist;
616 
617  if (iSub == KLMElementNumbers::c_BKLM) {
618  int iF = klmChannel.getSection();
619  int iS = klmChannel.getSector() - 1;
620  int iL = klmChannel.getLayer() - 1;
621  int iP = klmChannel.getPlane();
622  int iC = klmChannel.getStrip() - 1;
623  if (iL > 1) {
624  if (iP) {
625  profileRpcPhi->Fill(distHit, timeHit);
626  } else {
627  profileRpcZ->Fill(distHit, timeHit);
628  }
629  } else {
630  if (fill2dHistograms)
631  m_HistTimeLengthBKLM[iF][iS][iL][iP][iC]->Fill(distHit, timeHit);
632  if (iP) {
633  profileBKLMScintillatorPhi->Fill(distHit, timeHit);
634  } else {
635  profileBKLMScintillatorZ->Fill(distHit, timeHit);
636  }
637  }
638  } else {
639  int iF = klmChannel.getSection() - 1;
640  int iS = klmChannel.getSector() - 1;
641  int iL = klmChannel.getLayer() - 1;
642  int iP = klmChannel.getPlane() - 1;
643  int iC = klmChannel.getStrip() - 1;
644  if (fill2dHistograms)
645  m_HistTimeLengthEKLM[iF][iS][iL][iP][iC]->Fill(distHit, timeHit);
646  if (iP) {
647  profileEKLMScintillatorPlane1->Fill(distHit, timeHit);
648  } else {
649  profileEKLMScintillatorPlane2->Fill(distHit, timeHit);
650  }
651  }
652  }
653  }
654 }
655 
656 void KLMTimeAlgorithm::timeDistance2dFit(
657  const std::vector< std::pair<KLMChannelNumber, unsigned int> >& channels,
658  double& delay, double& delayError)
659 {
660  std::vector<struct Event>::iterator it;
661  std::vector<struct Event> eventsChannel;
662  int nFits = 1000;
663  int nConvergedFits = 0;
664  delay = 0;
665  delayError = 0;
666  if (nFits > (int)channels.size())
667  nFits = channels.size();
668  for (int i = 0; i < nFits; ++i) {
669  int subdetector, section, sector, layer, plane, strip;
670  m_ElementNumbers->channelNumberToElementNumbers(
671  channels[i].first, &subdetector, &section, &sector, &layer, &plane,
672  &strip);
673  if (subdetector == KLMElementNumbers::c_BKLM) {
674  s_LowerTimeBoundary = m_LowerTimeBoundaryScintilltorsBKLM;
675  s_UpperTimeBoundary = m_UpperTimeBoundaryScintilltorsBKLM;
676  const bklm::Module* module =
677  m_BKLMGeometry->findModule(section, sector, layer);
678  s_StripLength = module->getStripLength(plane, strip);
679  } else {
680  s_LowerTimeBoundary = m_LowerTimeBoundaryScintilltorsEKLM;
681  s_UpperTimeBoundary = m_UpperTimeBoundaryScintilltorsEKLM;
682  s_StripLength = m_EKLMGeometry->getStripLength(strip) /
683  CLHEP::cm * Unit::cm;
684  }
685  for (int j = 0; j < c_NBinsTime; ++j) {
686  for (int k = 0; k < c_NBinsDistance; ++k)
687  s_BinnedData[j][k] = 0;
688  }
689  eventsChannel = m_evts[channels[i].first];
690  double averageTime = 0;
691  for (it = eventsChannel.begin(); it != eventsChannel.end(); ++it) {
692  double timeHit = it->time();
693  if (m_useEventT0)
694  timeHit = timeHit - it->t0;
695  averageTime = averageTime + timeHit;
696  int timeBin = std::floor((timeHit - s_LowerTimeBoundary) * c_NBinsTime /
697  (s_UpperTimeBoundary - s_LowerTimeBoundary));
698  if (timeBin < 0 || timeBin >= c_NBinsTime)
699  continue;
700  int distanceBin = std::floor(it->dist * c_NBinsDistance / s_StripLength);
701  if (distanceBin < 0 || distanceBin >= c_NBinsDistance) {
702  B2ERROR("The distance to SiPM is greater than the strip length.");
703  continue;
704  }
705  s_BinnedData[timeBin][distanceBin] += 1;
706  }
707  averageTime = averageTime / eventsChannel.size();
708  TMinuit minuit(5);
709  minuit.SetPrintLevel(-1);
710  int minuitResult;
711  minuit.mnparm(0, "P0", 1, 0.001, 0, 0, minuitResult);
712  minuit.mnparm(1, "N", 10, 0.001, 0, 0, minuitResult);
713  minuit.mnparm(2, "T0", averageTime, 0.001, 0, 0, minuitResult);
714  minuit.mnparm(3, "SIGMA", 10, 0.001, 0, 0, minuitResult);
715  minuit.mnparm(4, "DELAY", 0.0, 0.001, 0, 0, minuitResult);
716  minuit.SetFCN(fcn);
717  minuit.mncomd("FIX 2 3 4 5", minuitResult);
718  minuit.mncomd("MIGRAD 10000", minuitResult);
719  minuit.mncomd("RELEASE 2", minuitResult);
720  minuit.mncomd("MIGRAD 10000", minuitResult);
721  minuit.mncomd("RELEASE 3", minuitResult);
722  minuit.mncomd("MIGRAD 10000", minuitResult);
723  minuit.mncomd("RELEASE 4", minuitResult);
724  minuit.mncomd("MIGRAD 10000", minuitResult);
725  minuit.mncomd("RELEASE 5", minuitResult);
726  minuit.mncomd("MIGRAD 10000", minuitResult);
727  /* Require converged fit with accurate error matrix. */
728  if (minuit.fISW[1] != 3)
729  continue;
730  nConvergedFits++;
731  double channelDelay, channelDelayError;
732  minuit.GetParameter(4, channelDelay, channelDelayError);
733  delay = delay + channelDelay;
734  delayError = delayError + channelDelayError * channelDelayError;
735  }
736  delay = delay / nConvergedFits;
737  delayError = sqrt(delayError) / (nConvergedFits - 1);
738 }
739 
740 CalibrationAlgorithm::EResult KLMTimeAlgorithm::calibrate()
741 {
742  int channelId;
743  gROOT->SetBatch(kTRUE);
744  setupDatabase();
745  m_timeCableDelay = new KLMTimeCableDelay();
746  m_timeConstants = new KLMTimeConstants();
747  m_timeResolution = new KLMTimeResolution();
748 
749  fcn_gaus = new TF1("fcn_gaus", "gaus");
750  fcn_land = new TF1("fcn_land", "landau");
751  fcn_pol1 = new TF1("fcn_pol1", "pol1");
752  fcn_const = new TF1("fcn_const", "pol0");
753 
754  CalibrationAlgorithm::EResult result = readCalibrationData();
755  if (result != CalibrationAlgorithm::c_OK)
756  return result;
757 
758  /* Choose non-existing file name. */
759  std::string name = "time_calibration.root";
760  int i = 1;
761  while (1) {
762  struct stat buffer;
763  if (stat(name.c_str(), &buffer) != 0)
764  break;
765  name = "time_calibration_" + std::to_string(i) + ".root";
766  i = i + 1;
767  /* Overflow. */
768  if (i < 0)
769  break;
770  }
771  m_outFile = new TFile(name.c_str(), "recreate");
772  createHistograms();
773 
774  std::vector<struct Event>::iterator it;
775  std::vector<struct Event> eventsChannel;
776 
777  eventsChannel.clear();
778  m_cFlag.clear();
779  m_minimizerOptions.SetDefaultStrategy(2);
780 
781  /* Sort channels by number of events. */
782  std::vector< std::pair<KLMChannelNumber, unsigned int> > channelsBKLM;
783  std::vector< std::pair<KLMChannelNumber, unsigned int> > channelsEKLM;
784  KLMChannelIndex klmChannels;
785  for (KLMChannelIndex& klmChannel : klmChannels) {
786  KLMChannelNumber channel = klmChannel.getKLMChannelNumber();
787  m_cFlag[channel] = ChannelCalibrationStatus::c_NotEnoughData;
788  if (m_evts.find(channel) == m_evts.end())
789  continue;
790  int nEvents = m_evts[channel].size();
791  if (nEvents < m_lower_limit_counts) {
792  B2WARNING("Not enough calibration data collected."
793  << LogVar("channel", channel)
794  << LogVar("number of digit", nEvents));
795  continue;
796  }
797  m_cFlag[channel] = ChannelCalibrationStatus::c_FailedFit;
798  if (klmChannel.getSubdetector() == KLMElementNumbers::c_BKLM &&
799  klmChannel.getLayer() < BKLMElementNumbers::c_FirstRPCLayer) {
800  channelsBKLM.push_back(
801  std::pair<KLMChannelNumber, unsigned int>(channel, nEvents));
802  }
803  if (klmChannel.getSubdetector() == KLMElementNumbers::c_EKLM) {
804  channelsEKLM.push_back(
805  std::pair<KLMChannelNumber, unsigned int>(channel, nEvents));
806  }
807  }
808  std::sort(channelsBKLM.begin(), channelsBKLM.end(), compareEventNumber);
809  std::sort(channelsEKLM.begin(), channelsEKLM.end(), compareEventNumber);
810 
811  /* Two-dimensional fit for the channel with the maximal number of events. */
812  double delayBKLM, delayBKLMError;
813  double delayEKLM, delayEKLMError;
814  timeDistance2dFit(channelsBKLM, delayBKLM, delayBKLMError);
815  timeDistance2dFit(channelsEKLM, delayEKLM, delayEKLMError);
816 
817  /**********************************
818  * First loop
819  * Estimation of effective light speed for Scintillators and RPCs, separately.
820  **********************************/
821  B2INFO("Effective light speed Estimation.");
822  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
823  channelId = klmChannel.getKLMChannelNumber();
824  if (m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
825  continue;
826 
827  eventsChannel = m_evts[channelId];
828  int iSub = klmChannel.getSubdetector();
829 
830  for (it = eventsChannel.begin(); it != eventsChannel.end(); ++it) {
831  TVector3 diffD = TVector3(it->diffDistX, it->diffDistY, it->diffDistZ);
832  h_diff->Fill(diffD.Mag());
833  double timeHit = it->time();
834  if (m_useEventT0)
835  timeHit = timeHit - it->t0;
836  if (iSub == KLMElementNumbers::c_BKLM) {
837  int iF = klmChannel.getSection();
838  int iS = klmChannel.getSector() - 1;
839  int iL = klmChannel.getLayer() - 1;
840  int iP = klmChannel.getPlane();
841  int iC = klmChannel.getStrip() - 1;
842  h_timeFSLPC_tc[iF][iS][iL][iP][iC]->Fill(timeHit);
843  if (iL > 1) {
844  h_time_rpc_tc->Fill(timeHit);
845  } else {
846  h_time_scint_tc->Fill(timeHit);
847  }
848  } else {
849  int iF = klmChannel.getSection() - 1;
850  int iS = klmChannel.getSector() - 1;
851  int iL = klmChannel.getLayer() - 1;
852  int iP = klmChannel.getPlane() - 1;
853  int iC = klmChannel.getStrip() - 1;
854  h_timeFSLPC_tc_end[iF][iS][iL][iP][iC]->Fill(timeHit);
855  h_time_scint_tc_end->Fill(timeHit);
856  }
857  }
858  }
859  B2INFO("Effective light speed Estimation! Hists and Graph filling done.");
860 
861  m_timeShift.clear();
862 
863  double tmpMean_rpc_global = h_time_rpc_tc->GetMean();
864  double tmpMean_scint_global = h_time_scint_tc->GetMean();
865  double tmpMean_scint_global_end = h_time_scint_tc_end->GetMean();
866 
867  B2INFO("Global Mean for Raw." << LogVar("RPC", tmpMean_rpc_global) << LogVar("Scint BKLM",
868  tmpMean_scint_global) << LogVar("Scint EKLM", tmpMean_scint_global_end));
869 
870  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
871  channelId = klmChannel.getKLMChannelNumber();
872  if (m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
873  continue;
874 
875  int iSub = klmChannel.getSubdetector();
876  if (iSub == KLMElementNumbers::c_BKLM) {
877  int iF = klmChannel.getSection();
878  int iS = klmChannel.getSector() - 1;
879  int iL = klmChannel.getLayer() - 1;
880  int iP = klmChannel.getPlane();
881  int iC = klmChannel.getStrip() - 1;
882  h_timeFSLPC_tc[iF][iS][iL][iP][iC]->Fit(fcn_gaus, "LESQ");
883  double tmpMean_channel = fcn_gaus->GetParameter(1);
884  if (iL > 1) {
885  m_timeShift[channelId] = tmpMean_channel - tmpMean_rpc_global;
886  } else {
887  m_timeShift[channelId] = tmpMean_channel - tmpMean_scint_global;
888  }
889  } else {
890  int iF = klmChannel.getSection() - 1;
891  int iS = klmChannel.getSector() - 1;
892  int iL = klmChannel.getLayer() - 1;
893  int iP = klmChannel.getPlane() - 1;
894  int iC = klmChannel.getStrip() - 1;
895  h_timeFSLPC_tc_end[iF][iS][iL][iP][iC]->Fit(fcn_gaus, "LESQ");
896  double tmpMean_channel = fcn_gaus->GetParameter(1);
897  m_timeShift[channelId] = tmpMean_channel - tmpMean_scint_global_end;
898  }
899  }
900 
901  delete h_time_scint_tc;
902  delete h_time_scint_tc_end;
903  delete h_time_rpc_tc;
904  B2INFO("Effective Light m_timeShift obtained. done.");
905 
906  fillTimeDistanceProfiles(
907  m_ProfileRpcPhi, m_ProfileRpcZ,
908  m_ProfileBKLMScintillatorPhi, m_ProfileBKLMScintillatorZ,
909  m_ProfileEKLMScintillatorPlane1, m_ProfileEKLMScintillatorPlane2, false);
910 
911  B2INFO("Effective light speed fitting.");
912  m_ProfileRpcPhi->Fit("fcn_pol1", "EMQ");
913  double delayRPCPhi = fcn_pol1->GetParameter(1);
914  double e_slope_rpc_phi = fcn_pol1->GetParError(1);
915 
916  m_ProfileRpcZ->Fit("fcn_pol1", "EMQ");
917  double delayRPCZ = fcn_pol1->GetParameter(1);
918  double e_slope_rpc_z = fcn_pol1->GetParError(1);
919 
920  m_ProfileBKLMScintillatorPhi->Fit("fcn_pol1", "EMQ");
921  double slope_scint_phi = fcn_pol1->GetParameter(1);
922  double e_slope_scint_phi = fcn_pol1->GetParError(1);
923 
924  m_ProfileBKLMScintillatorZ->Fit("fcn_pol1", "EMQ");
925  double slope_scint_z = fcn_pol1->GetParameter(1);
926  double e_slope_scint_z = fcn_pol1->GetParError(1);
927 
928  m_ProfileEKLMScintillatorPlane1->Fit("fcn_pol1", "EMQ");
929  double slope_scint_plane1_end = fcn_pol1->GetParameter(1);
930  double e_slope_scint_plane1_end = fcn_pol1->GetParError(1);
931 
932  m_ProfileEKLMScintillatorPlane2->Fit("fcn_pol1", "EMQ");
933  double slope_scint_plane2_end = fcn_pol1->GetParameter(1);
934  double e_slope_scint_plane2_end = fcn_pol1->GetParError(1);
935 
936  TString logStr_phi, logStr_z;
937  logStr_phi = Form("%.4f +/- %.4f ns/cm", delayRPCPhi, e_slope_rpc_phi);
938  logStr_z = Form("%.4f +/- %.4f ns/cm", delayRPCZ, e_slope_rpc_z);
939  B2INFO("Delay in RPCs:"
940  << LogVar("Fitted Value (phi readout) ", logStr_phi.Data())
941  << LogVar("Fitted Value (z readout) ", logStr_z.Data()));
942  logStr_phi = Form("%.4f +/- %.4f ns/cm", slope_scint_phi, e_slope_scint_phi);
943  logStr_z = Form("%.4f +/- %.4f ns/cm", slope_scint_z, e_slope_scint_z);
944  B2INFO("Delay in BKLM scintillators:"
945  << LogVar("Fitted Value (phi readout) ", logStr_phi.Data())
946  << LogVar("Fitted Value (z readout) ", logStr_z.Data()));
947  logStr_phi = Form("%.4f +/- %.4f ns/cm", slope_scint_plane1_end,
948  e_slope_scint_plane1_end);
949  logStr_z = Form("%.4f +/- %.4f ns/cm", slope_scint_plane2_end,
950  e_slope_scint_plane2_end);
951  B2INFO("Delay in EKLM scintillators:"
952  << LogVar("Fitted Value (plane1 readout) ", logStr_phi.Data())
953  << LogVar("Fitted Value (plane2 readout) ", logStr_z.Data()));
954 
955  logStr_z = Form("%.4f +/- %.4f ns/cm", delayBKLM, delayBKLMError);
956  B2INFO("Delay in BKLM scintillators:"
957  << LogVar("Fitted Value (2d fit) ", logStr_z.Data()));
958  logStr_z = Form("%.4f +/- %.4f ns/cm", delayEKLM, delayEKLMError);
959  B2INFO("Delay in EKLM scintillators:"
960  << LogVar("Fitted Value (2d fit) ", logStr_z.Data()));
961 
962  // Default Effective light speed in current Database
963  //delayEKLM = 0.5 * (slope_scint_plane1_end + slope_scint_plane2_end);
964  //delayBKLM = 0.5 * (slope_scint_phi + slope_scint_z);
965 
966  m_timeConstants->setDelay(delayEKLM, KLMTimeConstants::c_EKLM);
967  m_timeConstants->setDelay(delayBKLM, KLMTimeConstants::c_BKLM);
968  m_timeConstants->setDelay(delayRPCPhi, KLMTimeConstants::c_RPCPhi);
969  m_timeConstants->setDelay(delayRPCZ, KLMTimeConstants::c_RPCZ);
970 
972  B2INFO("Time distribution filling begins.");
973  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
974  channelId = klmChannel.getKLMChannelNumber();
975  int iSub = klmChannel.getSubdetector();
976 
977  if (m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
978  continue;
979  eventsChannel = m_evts[channelId];
980 
981  for (it = eventsChannel.begin(); it != eventsChannel.end(); ++it) {
982  double timeHit = it->time();
983  if (m_useEventT0)
984  timeHit = timeHit - it->t0;
985  if (iSub == KLMElementNumbers::c_BKLM) {
986  int iF = klmChannel.getSection();
987  int iS = klmChannel.getSector() - 1;
988  int iL = klmChannel.getLayer() - 1;
989  int iP = klmChannel.getPlane();
990  int iC = klmChannel.getStrip() - 1;
991  if (iL > 1) {
992  double propgationT;
993  if (iP == BKLMElementNumbers::c_ZPlane)
994  propgationT = it->dist * delayRPCZ;
995  else
996  propgationT = it->dist * delayRPCPhi;
997  double time = timeHit - propgationT;
998  h_time_rpc->Fill(time);
999  h_timeF_rpc[iF]->Fill(time);
1000  h_timeFS_rpc[iF][iS]->Fill(time);
1001  h_timeFSL[iF][iS][iL]->Fill(time);
1002  h_timeFSLP[iF][iS][iL][iP]->Fill(time);
1003  h_timeFSLPC[iF][iS][iL][iP][iC]->Fill(time);
1004  h2_timeF_rpc[iF]->Fill(iS, time);
1005  h2_timeFS[iF][iS]->Fill(iL, time);
1006  h2_timeFSLP[iF][iS][iL][iP]->Fill(iC, time);
1007  } else {
1008  double propgationT = it->dist * delayBKLM;
1009  double time = timeHit - propgationT;
1010  h_time_scint->Fill(time);
1011  h_timeF_scint[iF]->Fill(time);
1012  h_timeFS_scint[iF][iS]->Fill(time);
1013  h_timeFSL[iF][iS][iL]->Fill(time);
1014  h_timeFSLP[iF][iS][iL][iP]->Fill(time);
1015  h_timeFSLPC[iF][iS][iL][iP][iC]->Fill(time);
1016  h2_timeF_scint[iF]->Fill(iS, time);
1017  h2_timeFS[iF][iS]->Fill(iL, time);
1018  h2_timeFSLP[iF][iS][iL][iP]->Fill(iC, time);
1019  }
1020  } else {
1021  int iF = klmChannel.getSection() - 1;
1022  int iS = klmChannel.getSector() - 1;
1023  int iL = klmChannel.getLayer() - 1;
1024  int iP = klmChannel.getPlane() - 1;
1025  int iC = klmChannel.getStrip() - 1;
1026  double propgationT = it->dist * delayEKLM;
1027  double time = timeHit - propgationT;
1028  h_time_scint_end->Fill(time);
1029  h_timeF_scint_end[iF]->Fill(time);
1030  h_timeFS_scint_end[iF][iS]->Fill(time);
1031  h_timeFSL_end[iF][iS][iL]->Fill(time);
1032  h_timeFSLP_end[iF][iS][iL][iP]->Fill(time);
1033  h_timeFSLPC_end[iF][iS][iL][iP][iC]->Fill(time);
1034  h2_timeF_scint_end[iF]->Fill(iS, time);
1035  h2_timeFS_end[iF][iS]->Fill(iL, time);
1036  h2_timeFSLP_end[iF][iS][iL][iP]->Fill(iC, time);
1037  }
1038  }
1039  }
1040 
1041  B2INFO("Orignal filling done.");
1042 
1043  int iChannel_rpc = 0;
1044  int iChannel = 0;
1045  int iChannel_end = 0;
1046  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1047  channelId = klmChannel.getKLMChannelNumber();
1048  if (m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1049  continue;
1050  int iSub = klmChannel.getSubdetector();
1051 
1052  if (iSub == KLMElementNumbers::c_BKLM) {
1053  int iF = klmChannel.getSection();
1054  int iS = klmChannel.getSector() - 1;
1055  int iL = klmChannel.getLayer() - 1;
1056  int iP = klmChannel.getPlane();
1057  int iC = klmChannel.getStrip() - 1;
1058 
1059  TFitResultPtr r = h_timeFSLPC[iF][iS][iL][iP][iC]->Fit(fcn_gaus, "LESQ");
1060  if (int(r) != 0)
1061  continue;
1062  if (int(r) == 0)
1063  m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1064  m_time_channel[channelId] = fcn_gaus->GetParameter(1);
1065  m_etime_channel[channelId] = fcn_gaus->GetParError(1);
1066  if (iL > 1) {
1067  gre_time_channel_rpc->SetPoint(iChannel_rpc, channelId, m_time_channel[channelId]);
1068  gre_time_channel_rpc->SetPointError(iChannel_rpc, 0., m_etime_channel[channelId]);
1069  iChannel_rpc++;
1070  } else {
1071  gre_time_channel_scint->SetPoint(iChannel, channelId, m_time_channel[channelId]);
1072  gre_time_channel_scint->SetPointError(iChannel, 0., m_etime_channel[channelId]);
1073  iChannel++;
1074  }
1075  } else {
1076  int iF = klmChannel.getSection() - 1;
1077  int iS = klmChannel.getSector() - 1;
1078  int iL = klmChannel.getLayer() - 1;
1079  int iP = klmChannel.getPlane() - 1;
1080  int iC = klmChannel.getStrip() - 1;
1081 
1082  TFitResultPtr r = h_timeFSLPC_end[iF][iS][iL][iP][iC]->Fit(fcn_gaus, "LESQ");
1083  if (int(r) != 0)
1084  continue;
1085  if (int(r) == 0)
1086  m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1087  m_time_channel[channelId] = fcn_gaus->GetParameter(1);
1088  m_etime_channel[channelId] = fcn_gaus->GetParError(1);
1089  gre_time_channel_scint_end->SetPoint(iChannel_end, channelId, m_time_channel[channelId]);
1090  gre_time_channel_scint_end->SetPointError(iChannel_end, 0., m_etime_channel[channelId]);
1091  iChannel_end++;
1092  }
1093  }
1094 
1095  gre_time_channel_scint->Fit("fcn_const", "EMQ");
1096  m_time_channelAvg_scint = fcn_const->GetParameter(0);
1097  m_etime_channelAvg_scint = fcn_const->GetParError(0);
1098 
1099  gre_time_channel_scint_end->Fit("fcn_const", "EMQ");
1100  m_time_channelAvg_scint_end = fcn_const->GetParameter(0);
1101  m_etime_channelAvg_scint_end = fcn_const->GetParError(0);
1102 
1103  gre_time_channel_rpc->Fit("fcn_const", "EMQ");
1104  m_time_channelAvg_rpc = fcn_const->GetParameter(0);
1105  m_etime_channelAvg_rpc = fcn_const->GetParError(0);
1106 
1107  B2INFO("Channel's time distribution fitting done.");
1108  B2DEBUG(20, LogVar("Average time (RPC)", m_time_channelAvg_rpc)
1109  << LogVar("Average time (BKLM scintillators)", m_time_channelAvg_scint)
1110  << LogVar("Average time (EKLM scintillators)", m_time_channelAvg_scint_end));
1111 
1112  B2INFO("Calibrated channel's time distribution filling begins.");
1113 
1114  m_timeShift.clear();
1115  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1116  channelId = klmChannel.getKLMChannelNumber();
1117  h_calibrated->Fill(m_cFlag[channelId]);
1118  if (m_time_channel.find(channelId) == m_time_channel.end())
1119  continue;
1120  double timeShift = m_time_channel[channelId];
1121  int iSub = klmChannel.getSubdetector();
1122  if (iSub == KLMElementNumbers::c_BKLM) {
1123  int iL = klmChannel.getLayer() - 1;
1124  if (iL > 1)
1125  m_timeShift[channelId] = timeShift;
1126  else
1127  m_timeShift[channelId] = timeShift;
1128  } else {
1129  m_timeShift[channelId] = timeShift;
1130  }
1131  m_timeCableDelay->setTimeDelay(channelId, m_timeShift[channelId]);
1132  }
1133 
1134  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1135  channelId = klmChannel.getKLMChannelNumber();
1136  if (m_timeShift.find(channelId) != m_timeShift.end())
1137  continue;
1138  m_timeShift[channelId] = esti_timeShift(klmChannel);
1139  m_timeCableDelay->setTimeDelay(channelId, m_timeShift[channelId]);
1140  B2DEBUG(20, "Uncalibrated Estimation " << LogVar("Channel", channelId) << LogVar("Estimated value", m_timeShift[channelId]));
1141  }
1142 
1143  iChannel_rpc = 0;
1144  iChannel = 0;
1145  iChannel_end = 0;
1146  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1147  channelId = klmChannel.getKLMChannelNumber();
1148  if (m_timeShift.find(channelId) == m_timeShift.end()) {
1149  B2ERROR("!!! Not All Channels Calibration Constant Set. Error Happended on " << LogVar("Channel", channelId));
1150  continue;
1151  }
1152  int iSub = klmChannel.getSubdetector();
1153  if (iSub == KLMElementNumbers::c_BKLM) {
1154  int iL = klmChannel.getLayer();
1155  if (iL > 2) {
1156  gr_timeShift_channel_rpc->SetPoint(iChannel_rpc, channelId, m_timeShift[channelId]);
1157  iChannel_rpc++;
1158  } else {
1159  gr_timeShift_channel_scint->SetPoint(iChannel, channelId, m_timeShift[channelId]);
1160  iChannel++;
1161  }
1162  } else {
1163  gr_timeShift_channel_scint_end->SetPoint(iChannel_end, channelId, m_timeShift[channelId]);
1164  iChannel_end++;
1165  }
1166  }
1167 
1168  fillTimeDistanceProfiles(
1169  m_Profile2RpcPhi, m_Profile2RpcZ,
1170  m_Profile2BKLMScintillatorPhi, m_Profile2BKLMScintillatorZ,
1171  m_Profile2EKLMScintillatorPlane1, m_Profile2EKLMScintillatorPlane2, true);
1172 
1173  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1174  channelId = klmChannel.getKLMChannelNumber();
1175  int iSub = klmChannel.getSubdetector();
1176  eventsChannel = m_evts[channelId];
1177  for (it = eventsChannel.begin(); it != eventsChannel.end(); ++it) {
1178  double timeHit = it->time();
1179  if (m_useEventT0)
1180  timeHit = timeHit - it->t0;
1181  if (iSub == KLMElementNumbers::c_BKLM) {
1182  int iF = klmChannel.getSection();
1183  int iS = klmChannel.getSector() - 1;
1184  int iL = klmChannel.getLayer() - 1;
1185  int iP = klmChannel.getPlane();
1186  int iC = klmChannel.getStrip() - 1;
1187  if (iL > 1) {
1188  double propgationT;
1189  if (iP == BKLMElementNumbers::c_ZPlane)
1190  propgationT = it->dist * delayRPCZ;
1191  else
1192  propgationT = it->dist * delayRPCPhi;
1193  double time = timeHit - propgationT - m_timeShift[channelId];
1194  hc_time_rpc->Fill(time);
1195  hc_timeF_rpc[iF]->Fill(time);
1196  hc_timeFS_rpc[iF][iS]->Fill(time);
1197  hc_timeFSL[iF][iS][iL]->Fill(time);
1198  hc_timeFSLP[iF][iS][iL][iP]->Fill(time);
1199  hc_timeFSLPC[iF][iS][iL][iP][iC]->Fill(time);
1200  h2c_timeF_rpc[iF]->Fill(iS, time);
1201  h2c_timeFS[iF][iS]->Fill(iL, time);
1202  h2c_timeFSLP[iF][iS][iL][iP]->Fill(iC, time);
1203  } else {
1204  double propgationT = it->dist * delayBKLM;
1205  double time = timeHit - propgationT - m_timeShift[channelId];
1206  hc_time_scint->Fill(time);
1207  hc_timeF_scint[iF]->Fill(time);
1208  hc_timeFS_scint[iF][iS]->Fill(time);
1209  hc_timeFSL[iF][iS][iL]->Fill(time);
1210  hc_timeFSLP[iF][iS][iL][iP]->Fill(time);
1211  hc_timeFSLPC[iF][iS][iL][iP][iC]->Fill(time);
1212  h2c_timeF_scint[iF]->Fill(iS, time);
1213  h2c_timeFS[iF][iS]->Fill(iL, time);
1214  h2c_timeFSLP[iF][iS][iL][iP]->Fill(iC, time);
1215  }
1216  } else {
1217  int iF = klmChannel.getSection() - 1;
1218  int iS = klmChannel.getSector() - 1;
1219  int iL = klmChannel.getLayer() - 1;
1220  int iP = klmChannel.getPlane() - 1;
1221  int iC = klmChannel.getStrip() - 1;
1222  double propgationT = it->dist * delayEKLM;
1223  double time = timeHit - propgationT - m_timeShift[channelId];
1224  hc_time_scint_end->Fill(time);
1225  hc_timeF_scint_end[iF]->Fill(time);
1226  hc_timeFS_scint_end[iF][iS]->Fill(time);
1227  hc_timeFSL_end[iF][iS][iL]->Fill(time);
1228  hc_timeFSLP_end[iF][iS][iL][iP]->Fill(time);
1229  hc_timeFSLPC_end[iF][iS][iL][iP][iC]->Fill(time);
1230  h2c_timeF_scint_end[iF]->Fill(iS, time);
1231  h2c_timeFS_end[iF][iS]->Fill(iL, time);
1232  h2c_timeFSLP_end[iF][iS][iL][iP]->Fill(iC, time);
1233  }
1234  }
1235  }
1236 
1237  int icChannel_rpc = 0;
1238  int icChannel = 0;
1239  int icChannel_end = 0;
1240  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1241  channelId = klmChannel.getKLMChannelNumber();
1242  if (m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1243  continue;
1244  int iSub = klmChannel.getSubdetector();
1245 
1246  if (iSub == KLMElementNumbers::c_BKLM) {
1247  int iF = klmChannel.getSection();
1248  int iS = klmChannel.getSector() - 1;
1249  int iL = klmChannel.getLayer() - 1;
1250  int iP = klmChannel.getPlane();
1251  int iC = klmChannel.getStrip() - 1;
1252 
1253  TFitResultPtr rc = hc_timeFSLPC[iF][iS][iL][iP][iC]->Fit(fcn_gaus, "LESQ");
1254  if (int(rc) != 0)
1255  continue;
1256  if (int(rc) == 0)
1257  m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1258  m_ctime_channel[channelId] = fcn_gaus->GetParameter(1);
1259  mc_etime_channel[channelId] = fcn_gaus->GetParError(1);
1260  if (iL > 1) {
1261  gre_ctime_channel_rpc->SetPoint(icChannel_rpc, channelId, m_ctime_channel[channelId]);
1262  gre_ctime_channel_rpc->SetPointError(icChannel_rpc, 0., mc_etime_channel[channelId]);
1263  icChannel_rpc++;
1264  } else {
1265  gre_ctime_channel_scint->SetPoint(icChannel, channelId, m_ctime_channel[channelId]);
1266  gre_ctime_channel_scint->SetPointError(icChannel, 0., mc_etime_channel[channelId]);
1267  icChannel++;
1268  }
1269  } else {
1270  int iF = klmChannel.getSection() - 1;
1271  int iS = klmChannel.getSector() - 1;
1272  int iL = klmChannel.getLayer() - 1;
1273  int iP = klmChannel.getPlane() - 1;
1274  int iC = klmChannel.getStrip() - 1;
1275 
1276  TFitResultPtr rc = hc_timeFSLPC_end[iF][iS][iL][iP][iC]->Fit(fcn_gaus, "LESQ");
1277  if (int(rc) != 0)
1278  continue;
1279  if (int(rc) == 0)
1280  m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1281  m_ctime_channel[channelId] = fcn_gaus->GetParameter(1);
1282  mc_etime_channel[channelId] = fcn_gaus->GetParError(1);
1283  gre_ctime_channel_scint_end->SetPoint(icChannel_end, channelId, m_ctime_channel[channelId]);
1284  gre_ctime_channel_scint_end->SetPointError(icChannel_end, 0., mc_etime_channel[channelId]);
1285  icChannel_end++;
1286  }
1287  }
1288 
1289  gre_ctime_channel_scint->Fit("fcn_const", "EMQ");
1290  m_ctime_channelAvg_scint = fcn_const->GetParameter(0);
1291  mc_etime_channelAvg_scint = fcn_const->GetParError(0);
1292 
1293  gre_ctime_channel_scint_end->Fit("fcn_const", "EMQ");
1294  m_ctime_channelAvg_scint_end = fcn_const->GetParameter(0);
1295  mc_etime_channelAvg_scint_end = fcn_const->GetParError(0);
1296 
1297  gre_ctime_channel_rpc->Fit("fcn_const", "EMQ");
1298  m_ctime_channelAvg_rpc = fcn_const->GetParameter(0);
1299  mc_etime_channelAvg_rpc = fcn_const->GetParError(0);
1300 
1301  B2INFO("Channel's time distribution fitting done.");
1302  B2DEBUG(20, LogVar("Average calibrated time (RPC)", m_ctime_channelAvg_rpc)
1303  << LogVar("Average calibrated time (BKLM scintillators)", m_ctime_channelAvg_scint)
1304  << LogVar("Average calibrated time (EKLM scintillators)", m_ctime_channelAvg_scint_end));
1305 
1306  B2INFO("Calibrated channel's time distribution filling begins.");
1307 
1308  m_timeRes.clear();
1309  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1310  channelId = klmChannel.getKLMChannelNumber();
1311  hc_calibrated->Fill(m_cFlag[channelId]);
1312  if (m_ctime_channel.find(channelId) == m_ctime_channel.end())
1313  continue;
1314  double timeRes = m_ctime_channel[channelId];
1315  int iSub = klmChannel.getSubdetector();
1316  if (iSub == KLMElementNumbers::c_BKLM) {
1317  int iL = klmChannel.getLayer() - 1;
1318  if (iL > 1)
1319  m_timeRes[channelId] = timeRes;
1320  else
1321  m_timeRes[channelId] = timeRes;
1322  } else {
1323  m_timeRes[channelId] = timeRes;
1324  }
1325  m_timeResolution->setTimeResolution(channelId, m_timeRes[channelId]);
1326  }
1327 
1328  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1329  channelId = klmChannel.getKLMChannelNumber();
1330  if (m_timeRes.find(channelId) != m_timeRes.end())
1331  continue;
1332  m_timeRes[channelId] = esti_timeRes(klmChannel);
1333  m_timeResolution->setTimeResolution(channelId, m_timeRes[channelId]);
1334  B2DEBUG(20, "Calibrated Estimation " << LogVar("Channel", channelId) << LogVar("Estimated value", m_timeRes[channelId]));
1335  }
1336 
1337  icChannel_rpc = 0;
1338  icChannel = 0;
1339  icChannel_end = 0;
1340  for (KLMChannelIndex klmChannel = m_klmChannels.begin(); klmChannel != m_klmChannels.end(); ++klmChannel) {
1341  channelId = klmChannel.getKLMChannelNumber();
1342  if (m_timeRes.find(channelId) == m_timeRes.end()) {
1343  B2ERROR("!!! Not All Channels Calibration Constant Set. Error Happended on " << LogVar("Channel", channelId));
1344  continue;
1345  }
1346  int iSub = klmChannel.getSubdetector();
1347  if (iSub == KLMElementNumbers::c_BKLM) {
1348  int iL = klmChannel.getLayer();
1349  if (iL > 2) {
1350  gr_timeRes_channel_rpc->SetPoint(icChannel_rpc, channelId, m_timeRes[channelId]);
1351  icChannel_rpc++;
1352  } else {
1353  gr_timeRes_channel_scint->SetPoint(icChannel, channelId, m_timeRes[channelId]);
1354  icChannel++;
1355  }
1356  } else {
1357  gr_timeRes_channel_scint_end->SetPoint(icChannel_end, channelId, m_timeRes[channelId]);
1358  icChannel_end++;
1359  }
1360  }
1361  saveHist();
1362 
1363  delete fcn_const;
1364  m_evts.clear();
1365  m_timeShift.clear();
1366  m_timeRes.clear();
1367  m_cFlag.clear();
1368 
1369  saveCalibration(m_timeCableDelay, "KLMTimeCableDelay");
1370  saveCalibration(m_timeConstants, "KLMTimeConstants");
1371  saveCalibration(m_timeResolution, "KLMTimeResolution");
1372  return CalibrationAlgorithm::c_OK;
1373 }
1374 
1375 
1376 void KLMTimeAlgorithm::saveHist()
1377 {
1378  m_outFile->cd();
1379  B2INFO("Save Histograms into Files.");
1380  TDirectory* dir_monitor = m_outFile->mkdir("monitor_Hists");
1381  dir_monitor->cd();
1382  h_calibrated->SetDirectory(dir_monitor);
1383  hc_calibrated->SetDirectory(dir_monitor);
1384  h_diff->SetDirectory(dir_monitor);
1385 
1386  m_outFile->cd();
1387  TDirectory* dir_effC = m_outFile->mkdir("effC_Hists");
1388  dir_effC->cd();
1389  m_ProfileRpcPhi->SetDirectory(dir_effC);
1390  m_ProfileRpcZ->SetDirectory(dir_effC);
1391  m_ProfileBKLMScintillatorPhi->SetDirectory(dir_effC);
1392  m_ProfileBKLMScintillatorZ->SetDirectory(dir_effC);
1393  m_ProfileEKLMScintillatorPlane1->SetDirectory(dir_effC);
1394  m_ProfileEKLMScintillatorPlane2->SetDirectory(dir_effC);
1395  m_Profile2RpcPhi->SetDirectory(dir_effC);
1396  m_Profile2RpcZ->SetDirectory(dir_effC);
1397  m_Profile2BKLMScintillatorPhi->SetDirectory(dir_effC);
1398  m_Profile2BKLMScintillatorZ->SetDirectory(dir_effC);
1399  m_Profile2EKLMScintillatorPlane1->SetDirectory(dir_effC);
1400  m_Profile2EKLMScintillatorPlane2->SetDirectory(dir_effC);
1401 
1402  m_outFile->cd();
1403  TDirectory* dir_time = m_outFile->mkdir("time");
1404  dir_time->cd();
1405 
1406  h_time_scint->SetDirectory(dir_time);
1407  hc_time_scint->SetDirectory(dir_time);
1408 
1409  h_time_scint_end->SetDirectory(dir_time);
1410  hc_time_scint_end->SetDirectory(dir_time);
1411 
1412  h_time_rpc->SetDirectory(dir_time);
1413  hc_time_rpc->SetDirectory(dir_time);
1414 
1415  gre_time_channel_rpc->Write("gre_time_channel_rpc");
1416  gre_time_channel_scint->Write("gre_time_channel_scint");
1417  gre_time_channel_scint_end->Write("gre_time_channel_scint_end");
1418  gr_timeShift_channel_rpc->Write("gr_timeShift_channel_rpc");
1419  gr_timeShift_channel_scint->Write("gr_timeShift_channel_scint");
1420  gr_timeShift_channel_scint_end->Write("gr_timeShift_channel_scint_end");
1421  gre_ctime_channel_rpc->Write("gre_ctime_channel_rpc");
1422  gre_ctime_channel_scint->Write("gre_ctime_channel_scint");
1423  gre_ctime_channel_scint_end->Write("gre_ctime_channel_scint_end");
1424  gr_timeRes_channel_rpc->Write("gr_timeRes_channel_rpc");
1425  gr_timeRes_channel_scint->Write("gr_timeRes_channel_scint");
1426  gr_timeRes_channel_scint_end->Write("gr_timeRes_channel_scint_end");
1427 
1428  B2INFO("Top file setup Done.");
1429 
1430  TDirectory* dir_time_F[2];
1431  TDirectory* dir_time_FS[2][8];
1432  TDirectory* dir_time_FSL[2][8][15];
1433  TDirectory* dir_time_FSLP[2][8][15][2];
1434  TDirectory* dir_time_F_end[2];
1435  TDirectory* dir_time_FS_end[2][4];
1436  TDirectory* dir_time_FSL_end[2][4][14];
1437  TDirectory* dir_time_FSLP_end[2][4][14][2];
1438  char dirname[50];
1439  B2INFO("Sub files declare Done.");
1440  for (int iF = 0; iF < 2; ++iF) {
1441  h_timeF_rpc[iF]->SetDirectory(dir_time);
1442  hc_timeF_rpc[iF]->SetDirectory(dir_time);
1443 
1444  h2_timeF_rpc[iF]->SetDirectory(dir_time);
1445  h2c_timeF_rpc[iF]->SetDirectory(dir_time);
1446 
1447  h_timeF_scint[iF]->SetDirectory(dir_time);
1448  hc_timeF_scint[iF]->SetDirectory(dir_time);
1449 
1450  h2_timeF_scint[iF]->SetDirectory(dir_time);
1451  h2c_timeF_scint[iF]->SetDirectory(dir_time);
1452 
1453  h_timeF_scint_end[iF]->SetDirectory(dir_time);
1454  hc_timeF_scint_end[iF]->SetDirectory(dir_time);
1455 
1456  h2_timeF_scint_end[iF]->SetDirectory(dir_time);
1457  h2c_timeF_scint_end[iF]->SetDirectory(dir_time);
1458 
1459  sprintf(dirname, "isForward_%d", iF);
1460  dir_time_F[iF] = dir_time->mkdir(dirname);
1461  dir_time_F[iF]->cd();
1462 
1463  for (int iS = 0; iS < 8; ++iS) {
1464  h_timeFS_rpc[iF][iS]->SetDirectory(dir_time_F[iF]);
1465  hc_timeFS_rpc[iF][iS]->SetDirectory(dir_time_F[iF]);
1466 
1467  h_timeFS_scint[iF][iS]->SetDirectory(dir_time_F[iF]);
1468  hc_timeFS_scint[iF][iS]->SetDirectory(dir_time_F[iF]);
1469 
1470  h2_timeFS[iF][iS]->SetDirectory(dir_time_F[iF]);
1471  h2c_timeFS[iF][iS]->SetDirectory(dir_time_F[iF]);
1472 
1473  sprintf(dirname, "Sector_%d", iS + 1);
1474  dir_time_FS[iF][iS] = dir_time_F[iF]->mkdir(dirname);
1475  dir_time_FS[iF][iS]->cd();
1476 
1477  for (int iL = 0; iL < 15; ++iL) {
1478  h_timeFSL[iF][iS][iL]->SetDirectory(dir_time_FS[iF][iS]);
1479  hc_timeFSL[iF][iS][iL]->SetDirectory(dir_time_FS[iF][iS]);
1480 
1481  sprintf(dirname, "Layer_%d", iL + 1);
1482  dir_time_FSL[iF][iS][iL] = dir_time_FS[iF][iS]->mkdir(dirname);
1483  dir_time_FSL[iF][iS][iL]->cd();
1484  for (int iP = 0; iP < 2; ++iP) {
1485  h_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1486  hc_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1487  h2_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1488  h2c_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1489 
1490  sprintf(dirname, "Plane_%d", iP);
1491  dir_time_FSLP[iF][iS][iL][iP] = dir_time_FSL[iF][iS][iL]->mkdir(dirname);
1492  dir_time_FSLP[iF][iS][iL][iP]->cd();
1493 
1494  int nchannel_max = BKLMElementNumbers::getNStrips(iF, iS + 1, iL + 1, iP);
1495  for (int iC = 0; iC < nchannel_max; ++iC) {
1496  if (iL < 2)
1497  m_HistTimeLengthBKLM[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP[iF][iS][iL][iP]);
1498  h_timeFSLPC[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP[iF][iS][iL][iP]);
1499  hc_timeFSLPC[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP[iF][iS][iL][iP]);
1500  delete h_timeFSLPC_tc[iF][iS][iL][iP][iC];
1501  }
1502  }
1503  }
1504  }
1505 
1506  sprintf(dirname, "isForward_%d_end", iF + 1);
1507  dir_time_F_end[iF] = dir_time->mkdir(dirname);
1508  dir_time_F_end[iF]->cd();
1509  int maxLayer = 12 + 2 * iF;
1510  for (int iS = 0; iS < 4; ++iS) {
1511  h_timeFS_scint_end[iF][iS]->SetDirectory(dir_time_F_end[iF]);
1512  hc_timeFS_scint_end[iF][iS]->SetDirectory(dir_time_F_end[iF]);
1513 
1514  h2_timeFS_end[iF][iS]->SetDirectory(dir_time_F_end[iF]);
1515  h2c_timeFS_end[iF][iS]->SetDirectory(dir_time_F_end[iF]);
1516 
1517  sprintf(dirname, "Sector_%d_end", iS + 1);
1518  dir_time_FS_end[iF][iS] = dir_time_F_end[iF]->mkdir(dirname);
1519  dir_time_FS_end[iF][iS]->cd();
1520  for (int iL = 0; iL < maxLayer; ++iL) {
1521  h_timeFSL_end[iF][iS][iL]->SetDirectory(dir_time_FS_end[iF][iS]);
1522  hc_timeFSL_end[iF][iS][iL]->SetDirectory(dir_time_FS_end[iF][iS]);
1523 
1524  sprintf(dirname, "Layer_%d_end", iL + 1);
1525  dir_time_FSL_end[iF][iS][iL] = dir_time_FS_end[iF][iS]->mkdir(dirname);
1526  dir_time_FSL_end[iF][iS][iL]->cd();
1527  for (int iP = 0; iP < 2; ++iP) {
1528  h_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1529  hc_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1530  h2_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1531  h2c_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1532 
1533  sprintf(dirname, "plane_%d_end", iP);
1534  dir_time_FSLP_end[iF][iS][iL][iP] = dir_time_FSL_end[iF][iS][iL]->mkdir(dirname);
1535  dir_time_FSLP_end[iF][iS][iL][iP]->cd();
1536 
1537  for (int iC = 0; iC < 75; ++iC) {
1538  h_timeFSLPC_end[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP_end[iF][iS][iL][iP]);
1539  hc_timeFSLPC_end[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP_end[iF][iS][iL][iP]);
1540  m_HistTimeLengthEKLM[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP_end[iF][iS][iL][iP]);
1541  delete h_timeFSLPC_tc_end[iF][iS][iL][iP][iC];
1542  }
1543  }
1544  }
1545  }
1546  }
1547  m_outFile->cd();
1548  m_outFile->Write();
1549  m_outFile->Close();
1550  B2INFO("File Write and Close. Done.");
1551 }
1552 
1553 double KLMTimeAlgorithm::esti_timeShift(const KLMChannelIndex& klmChannel)
1554 {
1555  double tS = 0.0;
1556  int iSub = klmChannel.getSubdetector();
1557  int iF = klmChannel.getSection();
1558  int iS = klmChannel.getSector();
1559  int iL = klmChannel.getLayer();
1560  int iP = klmChannel.getPlane();
1561  int iC = klmChannel.getStrip();
1562  int totNStrips = EKLMElementNumbers::getMaximalStripNumber();
1563  if (iSub == KLMElementNumbers::c_BKLM)
1564  totNStrips = BKLMElementNumbers::getNStrips(iF, iS, iL, iP);
1565  if (iC == 1) {
1566  KLMChannelIndex kCIndex_upper(iSub, iF, iS, iL, iP, iC + 1);
1567  tS = tS_upperStrip(kCIndex_upper).second;
1568  } else if (iC == totNStrips) {
1569  KLMChannelIndex kCIndex_lower(iSub, iF, iS, iL, iP, iC - 1);
1570  tS = tS_lowerStrip(kCIndex_lower).second;
1571  } else {
1572  KLMChannelIndex kCIndex_upper(iSub, iF, iS, iL, iP, iC + 1);
1573  KLMChannelIndex kCIndex_lower(iSub, iF, iS, iL, iP, iC - 1);
1574  std::pair<int, double> tS_upper = tS_upperStrip(kCIndex_upper);
1575  std::pair<int, double> tS_lower = tS_lowerStrip(kCIndex_lower);
1576  unsigned int td_upper = tS_upper.first - iC;
1577  unsigned int td_lower = iC - tS_lower.first;
1578  unsigned int td = tS_upper.first - tS_lower.first;
1579  tS = (double(td_upper) * tS_lower.second + double(td_lower) * tS_upper.second) / double(td);
1580  }
1581  return tS;
1582 }
1583 
1584 std::pair<int, double> KLMTimeAlgorithm::tS_upperStrip(const KLMChannelIndex& klmChannel)
1585 {
1586  std::pair<int, double> tS;
1587  int cId = klmChannel.getKLMChannelNumber();
1588  int iSub = klmChannel.getSubdetector();
1589  int iF = klmChannel.getSection();
1590  int iS = klmChannel.getSector();
1591  int iL = klmChannel.getLayer();
1592  int iP = klmChannel.getPlane();
1593  int iC = klmChannel.getStrip();
1594  int totNStrips = EKLMElementNumbers::getMaximalStripNumber();
1595  if (iSub == KLMElementNumbers::c_BKLM)
1596  totNStrips = BKLMElementNumbers::getNStrips(iF, iS, iL, iP);
1597  if (m_timeShift.find(cId) != m_timeShift.end()) {
1598  tS.first = iC;
1599  tS.second = m_timeShift[cId];
1600  } else if (iC == totNStrips) {
1601  tS.first = iC;
1602  tS.second = 0.0;
1603  } else {
1604  KLMChannelIndex kCIndex(iSub, iF, iS, iL, iP, iC + 1);
1605  tS = tS_upperStrip(kCIndex);
1606  }
1607  return tS;
1608 }
1609 
1610 std::pair<int, double> KLMTimeAlgorithm::tS_lowerStrip(const KLMChannelIndex& klmChannel)
1611 {
1612  std::pair<int, double> tS;
1613  int cId = klmChannel.getKLMChannelNumber();
1614  int iSub = klmChannel.getSubdetector();
1615  int iF = klmChannel.getSection();
1616  int iS = klmChannel.getSector();
1617  int iL = klmChannel.getLayer();
1618  int iP = klmChannel.getPlane();
1619  int iC = klmChannel.getStrip();
1620  if (m_timeShift.find(cId) != m_timeShift.end()) {
1621  tS.first = iC;
1622  tS.second = m_timeShift[cId];
1623  } else if (iC == 1) {
1624  tS.first = iC;
1625  tS.second = 0.0;
1626  } else {
1627  KLMChannelIndex kCIndex(iSub, iF, iS, iL, iP, iC - 1);
1628  tS = tS_lowerStrip(kCIndex);
1629  }
1630  return tS;
1631 }
1632 
1633 double KLMTimeAlgorithm::esti_timeRes(const KLMChannelIndex& klmChannel)
1634 {
1635  double tR = 0.0;
1636  int iSub = klmChannel.getSubdetector();
1637  int iF = klmChannel.getSection();
1638  int iS = klmChannel.getSector();
1639  int iL = klmChannel.getLayer();
1640  int iP = klmChannel.getPlane();
1641  int iC = klmChannel.getStrip();
1642  int totNStrips = EKLMElementNumbers::getMaximalStripNumber();
1643  if (iSub == KLMElementNumbers::c_BKLM)
1644  totNStrips = BKLMElementNumbers::getNStrips(iF, iS, iL, iP);
1645  if (iC == 1) {
1646  KLMChannelIndex kCIndex_upper(iSub, iF, iS, iL, iP, iC + 1);
1647  tR = tR_upperStrip(kCIndex_upper).second;
1648  } else if (iC == totNStrips) {
1649  KLMChannelIndex kCIndex_lower(iSub, iF, iS, iL, iP, iC - 1);
1650  tR = tR_lowerStrip(kCIndex_lower).second;
1651  } else {
1652  KLMChannelIndex kCIndex_upper(iSub, iF, iS, iL, iP, iC + 1);
1653  KLMChannelIndex kCIndex_lower(iSub, iF, iS, iL, iP, iC - 1);
1654  std::pair<int, double> tR_upper = tR_upperStrip(kCIndex_upper);
1655  std::pair<int, double> tR_lower = tR_lowerStrip(kCIndex_lower);
1656  unsigned int tr_upper = tR_upper.first - iC;
1657  unsigned int tr_lower = iC - tR_lower.first;
1658  unsigned int tr = tR_upper.first - tR_lower.first;
1659  tR = (double(tr_upper) * tR_lower.second + double(tr_lower) * tR_upper.second) / double(tr);
1660  }
1661  return tR;
1662 }
1663 
1664 std::pair<int, double> KLMTimeAlgorithm::tR_upperStrip(const KLMChannelIndex& klmChannel)
1665 {
1666  std::pair<int, double> tR;
1667  int cId = klmChannel.getKLMChannelNumber();
1668  int iSub = klmChannel.getSubdetector();
1669  int iF = klmChannel.getSection();
1670  int iS = klmChannel.getSector();
1671  int iL = klmChannel.getLayer();
1672  int iP = klmChannel.getPlane();
1673  int iC = klmChannel.getStrip();
1674  int totNStrips = EKLMElementNumbers::getMaximalStripNumber();
1675  if (iSub == KLMElementNumbers::c_BKLM)
1676  totNStrips = BKLMElementNumbers::getNStrips(iF, iS, iL, iP);
1677  if (m_timeRes.find(cId) != m_timeRes.end()) {
1678  tR.first = iC;
1679  tR.second = m_timeRes[cId];
1680  } else if (iC == totNStrips) {
1681  tR.first = iC;
1682  tR.second = 0.0;
1683  } else {
1684  KLMChannelIndex kCIndex(iSub, iF, iS, iL, iP, iC + 1);
1685  tR = tR_upperStrip(kCIndex);
1686  }
1687  return tR;
1688 }
1689 
1690 std::pair<int, double> KLMTimeAlgorithm::tR_lowerStrip(const KLMChannelIndex& klmChannel)
1691 {
1692  std::pair<int, double> tR;
1693  int cId = klmChannel.getKLMChannelNumber();
1694  int iSub = klmChannel.getSubdetector();
1695  int iF = klmChannel.getSection();
1696  int iS = klmChannel.getSector();
1697  int iL = klmChannel.getLayer();
1698  int iP = klmChannel.getPlane();
1699  int iC = klmChannel.getStrip();
1700  if (m_timeRes.find(cId) != m_timeRes.end()) {
1701  tR.first = iC;
1702  tR.second = m_timeRes[cId];
1703  } else if (iC == 1) {
1704  tR.first = iC;
1705  tR.second = 0.0;
1706  } else {
1707  KLMChannelIndex kCIndex(iSub, iF, iS, iL, iP, iC - 1);
1708  tR = tR_lowerStrip(kCIndex);
1709  }
1710  return tR;
1711 }
1712 
Belle2::BKLMElementNumbers::c_FirstRPCLayer
@ c_FirstRPCLayer
First RPC layer.
Definition: BKLMElementNumbers.h:64
Belle2::BKLMElementNumbers::getNStrips
static int getNStrips(int section, int sector, int layer, int plane)
Get number of strips.
Definition: BKLMElementNumbers.cc:107
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:37
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:297
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:40
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:41
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:43
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:266
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:21
Belle2::DBStore
Singleton class to cache database objects.
Definition: DBStore.h:31
Belle2::DataStore::Instance
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
Belle2::DataStore::setInitializeActive
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
Belle2::EKLMElementNumbers::getMaximalStripNumber
static constexpr int getMaximalStripNumber()
Get maximal strip number.
Definition: EKLMElementNumbers.h:319
Belle2::EKLM::GeometryData::Instance
static const GeometryData & Instance(enum DataSource dataSource=c_Database, const GearDir *gearDir=nullptr)
Instantiation.
Definition: GeometryData.cc:33
Belle2::EKLM::GeometryData::getStripLength
double getStripLength(int strip) const
Get strip length.
Definition: GeometryData.h:71
Belle2::EKLM::GeometryData::getMaximalStripLength
double getMaximalStripLength() const
Get maximal strip length.
Definition: GeometryData.h:79
Belle2::KLMChannelIndex
KLM channel index.
Definition: KLMChannelIndex.h:23
Belle2::KLMChannelIndex::getSubdetector
int getSubdetector() const
Get subdetector.
Definition: KLMChannelIndex.h:106
Belle2::KLMChannelIndex::getLayer
int getLayer() const
Get layer.
Definition: KLMChannelIndex.h:130
Belle2::KLMChannelIndex::begin
KLMChannelIndex begin()
First channel.
Definition: KLMChannelIndex.h:199
Belle2::KLMChannelIndex::end
KLMChannelIndex & end()
Last channel.
Definition: KLMChannelIndex.h:207
Belle2::KLMChannelIndex::getSection
int getSection() const
Get section.
Definition: KLMChannelIndex.h:114
Belle2::KLMChannelIndex::getPlane
int getPlane() const
Get plane.
Definition: KLMChannelIndex.h:138
Belle2::KLMChannelIndex::getStrip
int getStrip() const
Get strip.
Definition: KLMChannelIndex.h:146
Belle2::KLMChannelIndex::getSector
int getSector() const
Get sector.
Definition: KLMChannelIndex.h:122
Belle2::KLMChannelIndex::getKLMChannelNumber
KLMChannelNumber getKLMChannelNumber() const
Get KLM channel number.
Definition: KLMChannelIndex.cc:143
Belle2::KLMElementNumbers::Instance
static const KLMElementNumbers & Instance()
Instantiation.
Definition: KLMElementNumbers.cc:29
Belle2::KLMElementNumbers::channelNumberToElementNumbers
void channelNumberToElementNumbers(KLMChannelNumber channel, int *subdetector, int *section, int *sector, int *layer, int *plane, int *strip) const
Get element numbers by channel number.
Definition: KLMElementNumbers.cc:104
Belle2::KLMTimeAlgorithm::m_Profile2EKLMScintillatorPlane2
TProfile * m_Profile2EKLMScintillatorPlane2
For EKLM scintillator plane2.
Definition: KLMTimeAlgorithm.h:518
Belle2::KLMTimeAlgorithm::mc_etime_channelAvg_rpc
double mc_etime_channelAvg_rpc
Calibrated central value error of the global time distribution (BKLM RPC part).
Definition: KLMTimeAlgorithm.h:367
Belle2::KLMTimeAlgorithm::m_HistTimeLengthEKLM
TH2F * m_HistTimeLengthEKLM[2][4][14][2][75]
Two-dimensional distributions of time versus propagation length.
Definition: KLMTimeAlgorithm.h:728
Belle2::KLMTimeAlgorithm::h_timeFSLPC_tc
TH1F * h_timeFSLPC_tc[2][8][15][2][54]
BKLM part, used for effective light speed estimation.
Definition: KLMTimeAlgorithm.h:713
Belle2::KLMTimeAlgorithm::h2c_timeF_scint_end
TH2F * h2c_timeF_scint_end[2]
EKLM part.
Definition: KLMTimeAlgorithm.h:610
Belle2::KLMTimeAlgorithm::h_timeFSLPC_tc_end
TH1F * h_timeFSLPC_tc_end[2][4][14][2][75]
EKLM part, used for effective light speed estimation.
Definition: KLMTimeAlgorithm.h:722
Belle2::KLMTimeAlgorithm::m_timeResolution
KLMTimeResolution * m_timeResolution
DBObject of time resolution.
Definition: KLMTimeAlgorithm.h:416
Belle2::KLMTimeAlgorithm::h_time_scint_tc_end
TH1F * h_time_scint_tc_end
EKLM part.
Definition: KLMTimeAlgorithm.h:532
Belle2::KLMTimeAlgorithm::createHistograms
void createHistograms()
Create histograms.
Definition: KLMTimeAlgorithm.cc:174
Belle2::KLMTimeAlgorithm::gre_time_channel_scint
TGraphErrors * gre_time_channel_scint
BKLM Scintillator.
Definition: KLMTimeAlgorithm.h:442
Belle2::KLMTimeAlgorithm::m_LowerTimeBoundaryCalibratedScintilltorsEKLM
double m_LowerTimeBoundaryCalibratedScintilltorsEKLM
Lower time boundary for EKLM scintillators (calibrated data).
Definition: KLMTimeAlgorithm.h:336
Belle2::KLMTimeAlgorithm::m_LowerTimeBoundaryScintilltorsEKLM
double m_LowerTimeBoundaryScintilltorsEKLM
Lower time boundary for EKLM scintillators.
Definition: KLMTimeAlgorithm.h:318
Belle2::KLMTimeAlgorithm::h_timeFSL
TH1F * h_timeFSL[2][8][15]
BKLM part.
Definition: KLMTimeAlgorithm.h:659
Belle2::KLMTimeAlgorithm::hc_timeFSLPC_end
TH1F * hc_timeFSLPC_end[2][4][14][2][75]
EKLM part.
Definition: KLMTimeAlgorithm.h:736
Belle2::KLMTimeAlgorithm::hc_timeFSL_end
TH1F * hc_timeFSL_end[2][4][14]
EKLM part.
Definition: KLMTimeAlgorithm.h:670
Belle2::KLMTimeAlgorithm::h_timeFSLP_end
TH1F * h_timeFSLP_end[2][4][14][2]
EKLM part.
Definition: KLMTimeAlgorithm.h:678
Belle2::KLMTimeAlgorithm::gr_timeRes_channel_rpc
TGraph * gr_timeRes_channel_rpc
BKLM RPC.
Definition: KLMTimeAlgorithm.h:472
Belle2::KLMTimeAlgorithm::hc_timeFSLP_end
TH1F * hc_timeFSLP_end[2][4][14][2]
EKLM part.
Definition: KLMTimeAlgorithm.h:686
Belle2::KLMTimeAlgorithm::h_timeFSLP
TH1F * h_timeFSLP[2][8][15][2]
BKLM part.
Definition: KLMTimeAlgorithm.h:675
Belle2::KLMTimeAlgorithm::hc_timeF_scint_end
TH1F * hc_timeF_scint_end[2]
EKLM part.
Definition: KLMTimeAlgorithm.h:582
Belle2::KLMTimeAlgorithm::m_timeShift
std::map< KLMChannelNumber, double > m_timeShift
Shift values of each channel.
Definition: KLMTimeAlgorithm.h:288
Belle2::KLMTimeAlgorithm::h_time_scint
TH1F * h_time_scint
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:540
Belle2::KLMTimeAlgorithm::m_time_channelAvg_scint
double m_time_channelAvg_scint
Central value of the global time distribution (BKLM scintillator part).
Definition: KLMTimeAlgorithm.h:349
Belle2::KLMTimeAlgorithm::hc_timeFS_scint_end
TH1F * hc_timeFS_scint_end[2][4]
EKLM part.
Definition: KLMTimeAlgorithm.h:632
Belle2::KLMTimeAlgorithm::esti_timeRes
double esti_timeRes(const KLMChannelIndex &klmChannel)
Estimate value of calibration constant for calibrated channels.
Definition: KLMTimeAlgorithm.cc:1633
Belle2::KLMTimeAlgorithm::m_UpperTimeBoundaryCalibratedRPC
double m_UpperTimeBoundaryCalibratedRPC
Upper time boundary for RPC (calibrated data).
Definition: KLMTimeAlgorithm.h:327
Belle2::KLMTimeAlgorithm::m_ctime_channelAvg_rpc
double m_ctime_channelAvg_rpc
Calibrated central value of the global time distribution (BKLM RPC part).
Definition: KLMTimeAlgorithm.h:364
Belle2::KLMTimeAlgorithm::m_timeConstants
KLMTimeConstants * m_timeConstants
DBObject of time cost on some parts of the detector.
Definition: KLMTimeAlgorithm.h:407
Belle2::KLMTimeAlgorithm::setupDatabase
void setupDatabase()
Setup the database.
Definition: KLMTimeAlgorithm.cc:104
Belle2::KLMTimeAlgorithm::m_LowerTimeBoundaryCalibratedScintilltorsBKLM
double m_LowerTimeBoundaryCalibratedScintilltorsBKLM
Lower time boundary for BKLM scintillators (calibrated data).
Definition: KLMTimeAlgorithm.h:330
Belle2::KLMTimeAlgorithm::m_UpperTimeBoundaryCalibratedScintilltorsBKLM
double m_UpperTimeBoundaryCalibratedScintilltorsBKLM
Upper time boundary for BKLM scintillators (calibrated data).
Definition: KLMTimeAlgorithm.h:333
Belle2::KLMTimeAlgorithm::hc_timeFS_scint
TH1F * hc_timeFS_scint[2][8]
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:629
Belle2::KLMTimeAlgorithm::m_time_channel
std::map< KLMChannelNumber, double > m_time_channel
Time distribution central value of each channel.
Definition: KLMTimeAlgorithm.h:294
Belle2::KLMTimeAlgorithm::m_ctime_channelAvg_scint_end
double m_ctime_channelAvg_scint_end
Calibrated central value of the global time distribution (EKLM scintillator part).
Definition: KLMTimeAlgorithm.h:379
Belle2::KLMTimeAlgorithm::readCalibrationData
CalibrationAlgorithm::EResult readCalibrationData()
Read calibration data.
Definition: KLMTimeAlgorithm.cc:141
Belle2::KLMTimeAlgorithm::gr_timeShift_channel_scint_end
TGraph * gr_timeShift_channel_scint_end
EKLM.
Definition: KLMTimeAlgorithm.h:467
Belle2::KLMTimeAlgorithm::gr_timeRes_channel_scint
TGraph * gr_timeRes_channel_scint
BKLM scintillator.
Definition: KLMTimeAlgorithm.h:475
Belle2::KLMTimeAlgorithm::hc_timeF_scint
TH1F * hc_timeF_scint[2]
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:579
Belle2::KLMTimeAlgorithm::h_timeFS_scint
TH1F * h_timeFS_scint[2][8]
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:618
Belle2::KLMTimeAlgorithm::m_ElementNumbers
const KLMElementNumbers * m_ElementNumbers
Element numbers.
Definition: KLMTimeAlgorithm.h:392
Belle2::KLMTimeAlgorithm::tR_upperStrip
std::pair< int, double > tR_upperStrip(const KLMChannelIndex &klmChannel)
Tracing avaiable channels with increasing strip number.
Definition: KLMTimeAlgorithm.cc:1664
Belle2::KLMTimeAlgorithm::hc_timeF_rpc
TH1F * hc_timeF_rpc[2]
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:576
Belle2::KLMTimeAlgorithm::h2c_timeFS_end
TH2F * h2c_timeFS_end[2][4]
EKLM part.
Definition: KLMTimeAlgorithm.h:654
Belle2::KLMTimeAlgorithm::h_timeFSLPC_end
TH1F * h_timeFSLPC_end[2][4][14][2][75]
EKLM part.
Definition: KLMTimeAlgorithm.h:725
Belle2::KLMTimeAlgorithm::m_EKLMGeometry
const EKLM::GeometryData * m_EKLMGeometry
EKLM geometry data.
Definition: KLMTimeAlgorithm.h:398
Belle2::KLMTimeAlgorithm::gre_ctime_channel_scint_end
TGraphErrors * gre_ctime_channel_scint_end
EKLM.
Definition: KLMTimeAlgorithm.h:456
Belle2::KLMTimeAlgorithm::m_Profile2BKLMScintillatorPhi
TProfile * m_Profile2BKLMScintillatorPhi
For BKLM scintillator phi plane.
Definition: KLMTimeAlgorithm.h:509
Belle2::KLMTimeAlgorithm::hc_time_scint_end
TH1F * hc_time_scint_end
EKLM part.
Definition: KLMTimeAlgorithm.h:554
Belle2::KLMTimeAlgorithm::gre_time_channel_scint_end
TGraphErrors * gre_time_channel_scint_end
EKLM.
Definition: KLMTimeAlgorithm.h:445
Belle2::KLMTimeAlgorithm::h2_timeFSLP
TH2F * h2_timeFSLP[2][8][15][2]
BKLM part.
Definition: KLMTimeAlgorithm.h:694
Belle2::KLMTimeAlgorithm::gr_timeShift_channel_scint
TGraph * gr_timeShift_channel_scint
BKLM scintillator.
Definition: KLMTimeAlgorithm.h:464
Belle2::KLMTimeAlgorithm::m_time_channelAvg_scint_end
double m_time_channelAvg_scint_end
Central value of the global time distribution (EKLM scintillator part).
Definition: KLMTimeAlgorithm.h:357
Belle2::KLMTimeAlgorithm::m_Profile2EKLMScintillatorPlane1
TProfile * m_Profile2EKLMScintillatorPlane1
For EKLM scintillator plane1.
Definition: KLMTimeAlgorithm.h:515
Belle2::KLMTimeAlgorithm::hc_timeFS_rpc
TH1F * hc_timeFS_rpc[2][8]
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:626
Belle2::KLMTimeAlgorithm::fillTimeDistanceProfiles
void fillTimeDistanceProfiles(TProfile *profileRpcPhi, TProfile *profileRpcZ, TProfile *profileBKLMScintillatorPhi, TProfile *profileBKLMScintillatorZ, TProfile *profileEKLMScintillatorPlane1, TProfile *profileEKLMScintillatorPlane2, bool fill2dHistograms)
Fill profiles of time versus distance.
Definition: KLMTimeAlgorithm.cc:595
Belle2::KLMTimeAlgorithm::m_outFile
TFile * m_outFile
Output file.
Definition: KLMTimeAlgorithm.h:753
Belle2::KLMTimeAlgorithm::esti_timeShift
double esti_timeShift(const KLMChannelIndex &klmChannel)
Estimate value of calibration constant for uncalibrated channels.
Definition: KLMTimeAlgorithm.cc:1553
Belle2::KLMTimeAlgorithm::tS_upperStrip
std::pair< int, double > tS_upperStrip(const KLMChannelIndex &klmChannel)
Tracing avaiable channels with increasing strip number.
Definition: KLMTimeAlgorithm.cc:1584
Belle2::KLMTimeAlgorithm::hc_timeFSLP
TH1F * hc_timeFSLP[2][8][15][2]
BKLM part.
Definition: KLMTimeAlgorithm.h:683
Belle2::KLMTimeAlgorithm::gre_ctime_channel_rpc
TGraphErrors * gre_ctime_channel_rpc
BKLM RPC.
Definition: KLMTimeAlgorithm.h:450
Belle2::KLMTimeAlgorithm::saveHist
void saveHist()
Save histograms to file.
Definition: KLMTimeAlgorithm.cc:1376
Belle2::KLMTimeAlgorithm::m_BKLMGeometry
const bklm::GeometryPar * m_BKLMGeometry
BKLM geometry data.
Definition: KLMTimeAlgorithm.h:395
Belle2::KLMTimeAlgorithm::h2c_timeFSLP
TH2F * h2c_timeFSLP[2][8][15][2]
BKLM part.
Definition: KLMTimeAlgorithm.h:705
Belle2::KLMTimeAlgorithm::m_ctime_channelAvg_scint
double m_ctime_channelAvg_scint
Calibrated central value of the global time distribution (BKLM scintillator part).
Definition: KLMTimeAlgorithm.h:371
Belle2::KLMTimeAlgorithm::fcn_const
TF1 * fcn_const
Const function.
Definition: KLMTimeAlgorithm.h:744
Belle2::KLMTimeAlgorithm::m_Profile2RpcZ
TProfile * m_Profile2RpcZ
For BKLM RPC z plane.
Definition: KLMTimeAlgorithm.h:506
Belle2::KLMTimeAlgorithm::h_timeFSLPC
TH1F * h_timeFSLPC[2][8][15][2][54]
BKLM part.
Definition: KLMTimeAlgorithm.h:716
Belle2::KLMTimeAlgorithm::m_HistTimeLengthBKLM
TH2F * m_HistTimeLengthBKLM[2][8][15][2][54]
Two-dimensional distributions of time versus propagation length.
Definition: KLMTimeAlgorithm.h:719
Belle2::KLMTimeAlgorithm::h2_timeFSLP_end
TH2F * h2_timeFSLP_end[2][4][14][2]
EKLM part.
Definition: KLMTimeAlgorithm.h:697
Belle2::KLMTimeAlgorithm::hc_calibrated
TH1I * hc_calibrated
Calibration statistics for each channel.
Definition: KLMTimeAlgorithm.h:431
Belle2::KLMTimeAlgorithm::timeDistance2dFit
void timeDistance2dFit(const std::vector< std::pair< KLMChannelNumber, unsigned int > > &channels, double &delay, double &delayError)
Two-dimensional fit for individual channels.
Definition: KLMTimeAlgorithm.cc:656
Belle2::KLMTimeAlgorithm::gr_timeRes_channel_scint_end
TGraph * gr_timeRes_channel_scint_end
EKLM.
Definition: KLMTimeAlgorithm.h:478
Belle2::KLMTimeAlgorithm::h_calibrated
TH1I * h_calibrated
Calibration statistics for each channel.
Definition: KLMTimeAlgorithm.h:428
Belle2::KLMTimeAlgorithm::m_ProfileBKLMScintillatorZ
TProfile * m_ProfileBKLMScintillatorZ
For BKLM scintillator z plane.
Definition: KLMTimeAlgorithm.h:492
Belle2::KLMTimeAlgorithm::h_time_rpc_tc
TH1F * h_time_rpc_tc
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:526
Belle2::KLMTimeAlgorithm::h_time_scint_end
TH1F * h_time_scint_end
EKLM part.
Definition: KLMTimeAlgorithm.h:543
Belle2::KLMTimeAlgorithm::h2c_timeF_scint
TH2F * h2c_timeF_scint[2]
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:607
Belle2::KLMTimeAlgorithm::fcn_pol1
TF1 * fcn_pol1
Pol1 function.
Definition: KLMTimeAlgorithm.h:741
Belle2::KLMTimeAlgorithm::m_etime_channelAvg_scint_end
double m_etime_channelAvg_scint_end
Central value error of the global time distribution (EKLM scintillator part).
Definition: KLMTimeAlgorithm.h:361
Belle2::KLMTimeAlgorithm::hc_time_rpc
TH1F * hc_time_rpc
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:548
Belle2::KLMTimeAlgorithm::mc_etime_channelAvg_scint
double mc_etime_channelAvg_scint
Calibrated central value error of the global time distribution (BKLM scintillator part).
Definition: KLMTimeAlgorithm.h:375
Belle2::KLMTimeAlgorithm::h2c_timeFSLP_end
TH2F * h2c_timeFSLP_end[2][4][14][2]
EKLM part.
Definition: KLMTimeAlgorithm.h:708
Belle2::KLMTimeAlgorithm::mc_etime_channelAvg_scint_end
double mc_etime_channelAvg_scint_end
Calibrated central value error of the global time distribution (EKLM scintillator part).
Definition: KLMTimeAlgorithm.h:383
Belle2::KLMTimeAlgorithm::h2_timeF_scint_end
TH2F * h2_timeF_scint_end[2]
EKLM part.
Definition: KLMTimeAlgorithm.h:596
Belle2::KLMTimeAlgorithm::m_UpperTimeBoundaryScintilltorsBKLM
double m_UpperTimeBoundaryScintilltorsBKLM
Upper time boundary for BKLM scintillators.
Definition: KLMTimeAlgorithm.h:315
Belle2::KLMTimeAlgorithm::fcn_gaus
TF1 * fcn_gaus
Gaussian function.
Definition: KLMTimeAlgorithm.h:747
Belle2::KLMTimeAlgorithm::h_timeF_rpc
TH1F * h_timeF_rpc[2]
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:562
Belle2::KLMTimeAlgorithm::m_ProfileBKLMScintillatorPhi
TProfile * m_ProfileBKLMScintillatorPhi
For BKLM scintillator phi plane.
Definition: KLMTimeAlgorithm.h:489
Belle2::KLMTimeAlgorithm::hc_time_scint
TH1F * hc_time_scint
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:551
Belle2::KLMTimeAlgorithm::h2_timeFS
TH2F * h2_timeFS[2][8]
BKLM part.
Definition: KLMTimeAlgorithm.h:640
Belle2::KLMTimeAlgorithm::m_etime_channelAvg_scint
double m_etime_channelAvg_scint
Central value error of the global time distribution (BKLM scintillator part).
Definition: KLMTimeAlgorithm.h:353
Belle2::KLMTimeAlgorithm::m_UpperTimeBoundaryCalibratedScintilltorsEKLM
double m_UpperTimeBoundaryCalibratedScintilltorsEKLM
Upper time boundary for BKLM scintillators (calibrated data).
Definition: KLMTimeAlgorithm.h:339
Belle2::KLMTimeAlgorithm::h_timeF_scint_end
TH1F * h_timeF_scint_end[2]
EKLM part.
Definition: KLMTimeAlgorithm.h:568
Belle2::KLMTimeAlgorithm::m_ProfileRpcPhi
TProfile * m_ProfileRpcPhi
For BKLM RPC phi plane.
Definition: KLMTimeAlgorithm.h:483
Belle2::KLMTimeAlgorithm::gre_ctime_channel_scint
TGraphErrors * gre_ctime_channel_scint
BKLM Scintillator.
Definition: KLMTimeAlgorithm.h:453
Belle2::KLMTimeAlgorithm::m_ProfileEKLMScintillatorPlane2
TProfile * m_ProfileEKLMScintillatorPlane2
For EKLM scintillator plane2.
Definition: KLMTimeAlgorithm.h:498
Belle2::KLMTimeAlgorithm::h_time_rpc
TH1F * h_time_rpc
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:537
Belle2::KLMTimeAlgorithm::h_timeFSL_end
TH1F * h_timeFSL_end[2][4][14]
EKLM part.
Definition: KLMTimeAlgorithm.h:662
Belle2::KLMTimeAlgorithm::m_Profile2RpcPhi
TProfile * m_Profile2RpcPhi
For BKLM RPC phi plane.
Definition: KLMTimeAlgorithm.h:503
Belle2::KLMTimeAlgorithm::h_timeF_scint
TH1F * h_timeF_scint[2]
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:565
Belle2::KLMTimeAlgorithm::hc_timeFSL
TH1F * hc_timeFSL[2][8][15]
BKLM part.
Definition: KLMTimeAlgorithm.h:667
Belle2::KLMTimeAlgorithm::m_Profile2BKLMScintillatorZ
TProfile * m_Profile2BKLMScintillatorZ
For BKLM scintillator z plane.
Definition: KLMTimeAlgorithm.h:512
Belle2::KLMTimeAlgorithm::h_timeFS_rpc
TH1F * h_timeFS_rpc[2][8]
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:615
Belle2::KLMTimeAlgorithm::m_klmChannels
KLMChannelIndex m_klmChannels
KLM ChannelIndex object.
Definition: KLMTimeAlgorithm.h:401
Belle2::KLMTimeAlgorithm::gr_timeShift_channel_rpc
TGraph * gr_timeShift_channel_rpc
BKLM RPC.
Definition: KLMTimeAlgorithm.h:461
Belle2::KLMTimeAlgorithm::m_timeRes
std::map< KLMChannelNumber, double > m_timeRes
Resolution values of each channel.
Definition: KLMTimeAlgorithm.h:291
Belle2::KLMTimeAlgorithm::h_diff
TH1F * h_diff
Distance between global and local position.
Definition: KLMTimeAlgorithm.h:434
Belle2::KLMTimeAlgorithm::h2_timeF_scint
TH2F * h2_timeF_scint[2]
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:593
Belle2::KLMTimeAlgorithm::h_time_scint_tc
TH1F * h_time_scint_tc
BKLM scintillator part.
Definition: KLMTimeAlgorithm.h:529
Belle2::KLMTimeAlgorithm::m_LowerTimeBoundaryRPC
double m_LowerTimeBoundaryRPC
Lower time boundary for RPC.
Definition: KLMTimeAlgorithm.h:306
Belle2::KLMTimeAlgorithm::calibrate
virtual EResult calibrate() override
Run algorithm on data.
Definition: KLMTimeAlgorithm.cc:740
Belle2::KLMTimeAlgorithm::m_ctime_channel
std::map< KLMChannelNumber, double > m_ctime_channel
Calibrated time distribution central value of each channel.
Definition: KLMTimeAlgorithm.h:300
Belle2::KLMTimeAlgorithm::m_LowerTimeBoundaryCalibratedRPC
double m_LowerTimeBoundaryCalibratedRPC
Lower time boundary for RPC (calibrated data).
Definition: KLMTimeAlgorithm.h:324
Belle2::KLMTimeAlgorithm::m_useEventT0
bool m_useEventT0
Whether to use event T0 from CDC.
Definition: KLMTimeAlgorithm.h:425
Belle2::KLMTimeAlgorithm::m_MinimalDigitNumber
int m_MinimalDigitNumber
Minimal digit number (total).
Definition: KLMTimeAlgorithm.h:386
Belle2::KLMTimeAlgorithm::m_ProfileEKLMScintillatorPlane1
TProfile * m_ProfileEKLMScintillatorPlane1
For EKLM scintillator plane1.
Definition: KLMTimeAlgorithm.h:495
Belle2::KLMTimeAlgorithm::m_UpperTimeBoundaryRPC
double m_UpperTimeBoundaryRPC
Upper time boundary for RPC.
Definition: KLMTimeAlgorithm.h:309
Belle2::KLMTimeAlgorithm::h2c_timeF_rpc
TH2F * h2c_timeF_rpc[2]
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:604
Belle2::KLMTimeAlgorithm::fcn_land
TF1 * fcn_land
Landau function.
Definition: KLMTimeAlgorithm.h:750
Belle2::KLMTimeAlgorithm::m_timeCableDelay
KLMTimeCableDelay * m_timeCableDelay
DBObject of the calibration constant of each channel due to cable decay.
Definition: KLMTimeAlgorithm.h:413
Belle2::KLMTimeAlgorithm::tS_lowerStrip
std::pair< int, double > tS_lowerStrip(const KLMChannelIndex &klmChannel)
Tracing avaiable channels with decreasing strip number.
Definition: KLMTimeAlgorithm.cc:1610
Belle2::KLMTimeAlgorithm::m_ProfileRpcZ
TProfile * m_ProfileRpcZ
For BKLM RPC z plane.
Definition: KLMTimeAlgorithm.h:486
Belle2::KLMTimeAlgorithm::mc_etime_channel
std::map< KLMChannelNumber, double > mc_etime_channel
Calibrated time distribution central value Error of each channel.
Definition: KLMTimeAlgorithm.h:303
Belle2::KLMTimeAlgorithm::tR_lowerStrip
std::pair< int, double > tR_lowerStrip(const KLMChannelIndex &klmChannel)
Tracing avaiable channels with decreasing strip number.
Definition: KLMTimeAlgorithm.cc:1690
Belle2::KLMTimeAlgorithm::m_LowerTimeBoundaryScintilltorsBKLM
double m_LowerTimeBoundaryScintilltorsBKLM
Lower time boundary for BKLM scintillators.
Definition: KLMTimeAlgorithm.h:312
Belle2::KLMTimeAlgorithm::hc_timeFSLPC
TH1F * hc_timeFSLPC[2][8][15][2][54]
BKLM part.
Definition: KLMTimeAlgorithm.h:733
Belle2::KLMTimeAlgorithm::m_minimizerOptions
ROOT::Math::MinimizerOptions m_minimizerOptions
Minimization options.
Definition: KLMTimeAlgorithm.h:404
Belle2::KLMTimeAlgorithm::m_cFlag
std::map< KLMChannelNumber, int > m_cFlag
Calibration flag if the channel has enough hits collected and fitted OK.
Definition: KLMTimeAlgorithm.h:285
Belle2::KLMTimeAlgorithm::gre_time_channel_rpc
TGraphErrors * gre_time_channel_rpc
BKLM RPC.
Definition: KLMTimeAlgorithm.h:439
Belle2::KLMTimeAlgorithm::h2_timeFS_end
TH2F * h2_timeFS_end[2][4]
EKLM part.
Definition: KLMTimeAlgorithm.h:643
Belle2::KLMTimeAlgorithm::h2_timeF_rpc
TH2F * h2_timeF_rpc[2]
BKLM RPC part.
Definition: KLMTimeAlgorithm.h:590
Belle2::KLMTimeAlgorithm::m_evts
std::map< KLMChannelNumber, std::vector< struct Event > > m_evts
Container of hit information.
Definition: KLMTimeAlgorithm.h:279
Belle2::KLMTimeAlgorithm::h_timeFS_scint_end
TH1F * h_timeFS_scint_end[2][4]
EKLM part.
Definition: KLMTimeAlgorithm.h:621
Belle2::KLMTimeAlgorithm::m_time_channelAvg_rpc
double m_time_channelAvg_rpc
Central value of the global time distribution (BKLM RPC part).
Definition: KLMTimeAlgorithm.h:342
Belle2::KLMTimeAlgorithm::h2c_timeFS
TH2F * h2c_timeFS[2][8]
BKLM part.
Definition: KLMTimeAlgorithm.h:651
Belle2::KLMTimeAlgorithm::m_etime_channelAvg_rpc
double m_etime_channelAvg_rpc
Central value error of the global time distribution (BKLM RPC part).
Definition: KLMTimeAlgorithm.h:345
Belle2::KLMTimeAlgorithm::m_etime_channel
std::map< KLMChannelNumber, double > m_etime_channel
Time distribution central value Error of each channel.
Definition: KLMTimeAlgorithm.h:297
Belle2::KLMTimeAlgorithm::m_UpperTimeBoundaryScintilltorsEKLM
double m_UpperTimeBoundaryScintilltorsEKLM
Upper time boundary for BKLM scintillators.
Definition: KLMTimeAlgorithm.h:321
Belle2::KLMTimeAlgorithm::m_lower_limit_counts
int m_lower_limit_counts
Lower limit of hits collected for on single channel.
Definition: KLMTimeAlgorithm.h:389
Belle2::KLMTimeCableDelay
Class to store BKLM delay time coused by cable in the database.
Definition: KLMTimeCableDelay.h:29
Belle2::KLMTimeCableDelay::setTimeDelay
void setTimeDelay(KLMChannelNumber channel, float timeDelay)
Set the time calibration constant.
Definition: KLMTimeCableDelay.cc:25
Belle2::KLMTimeConstants
Class to store KLM constants related to time.
Definition: KLMTimeConstants.h:23
Belle2::KLMTimeConstants::c_BKLM
@ c_BKLM
BKLM scintillator.
Definition: KLMTimeConstants.h:36
Belle2::KLMTimeConstants::c_EKLM
@ c_EKLM
EKLM scintillator.
Definition: KLMTimeConstants.h:33
Belle2::KLMTimeConstants::setDelay
void setDelay(float delay, int cType)
Set effective light speed of scintillators.
Definition: KLMTimeConstants.cc:32
Belle2::KLMTimeResolution
Class to store KLM time resolution in the database.
Definition: KLMTimeResolution.h:30
Belle2::KLMTimeResolution::setTimeResolution
void setTimeResolution(KLMChannelNumber channel, float resolution)
Set time resolution.
Definition: KLMTimeResolution.cc:26
Belle2::StoreAccessorBase::registerInDataStore
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Definition: StoreAccessorBase.h:52
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:111
Belle2::StoreObjPtr::construct
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
Definition: StoreObjPtr.h:119
Belle2::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:47
Belle2::bklm::GeometryPar::findModule
const Module * findModule(int section, int sector, int layer) const
Get the pointer to the definition of a module.
Definition: GeometryPar.cc:721
Belle2::bklm::GeometryPar::getMaximalZStripLength
double getMaximalZStripLength() const
Get maximal Z strip length (for scintillators).
Definition: GeometryPar.h:289
Belle2::bklm::GeometryPar::instance
static GeometryPar * instance(void)
Static method to get a reference to the singleton GeometryPar instance.
Definition: GeometryPar.cc:27
Belle2::bklm::GeometryPar::getMaximalPhiStripLength
double getMaximalPhiStripLength() const
Get maximal phi strip length (for scintillators).
Definition: GeometryPar.h:283
Belle2::bklm::Module
Define the geometry of a BKLM module Each sector [octant] contains Modules.
Definition: Module.h:76
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
Belle2::DBStore::Instance
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:28
Belle2::DBStore::updateEvent
void updateEvent()
Updates all intra-run dependent objects.
Definition: DBStore.cc:142
Belle2::DBStore::update
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:79
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::KLMChannelNumber
uint16_t KLMChannelNumber
Channel number.
Definition: KLMElementNumberDefinitions.h:36
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17