193 int nBin_scint = 400;
195 TString iFstring[2] = {
"Backward",
"Forward"};
196 TString iPstring[2] = {
"ZReadout",
"PhiReadout"};
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);
219 double maximalPhiStripLengthBKLM =
221 double maximalZStripLengthBKLM =
223 double maximalStripLengthEKLM =
227 "Time over propagation length for RPCs (Phi_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
230 "Time over propagation length for RPCs (Z_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
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);
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);
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);
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);
246 "Time over propagation length for RPCs (Phi_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
249 "Time over propagation length for RPCs (Z_Readout); propagation distance[cm]; T_rec-T_0-T_fly-'T_calibration'[ns]", 400, 0.0,
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);
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);
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);
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);
265 h_time_scint_tc =
new TH1F(
"h_time_scint_tc",
"time distribution for Scintillator", nBin_scint,
267 h_time_scint_tc_end =
new TH1F(
"h_time_scint_tc_end",
"time distribution for Scintillator (Endcap)", nBin_scint,
274 h_time_scint =
new TH1F(
"h_time_scint",
"time distribution for Scintillator; T_rec-T_0-T_fly-T_propagation[ns]", nBin_scint,
276 h_time_scint_end =
new TH1F(
"h_time_scint_end",
"time distribution for Scintillator (Endcap); T_rec-T_0-T_fly-T_propagation[ns]",
279 hc_time_rpc =
new TH1F(
"hc_time_rpc",
"Calibrated time distribution for RPC; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
282 "Calibrated time distribution for Scintillator; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", nBin_scint,
286 "Calibrated time distribution for Scintillator (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", nBin_scint,
289 for (
int iF = 0; iF < 2; ++iF) {
290 hn = Form(
"h_timeF%d_rpc", iF);
291 ht = Form(
"Time distribution for RPC of %s; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
293 hn = Form(
"h_timeF%d_scint", iF);
294 ht = Form(
"Time distribution for Scintillator of %s; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
297 hn = Form(
"h_timeF%d_scint_end", iF);
298 ht = Form(
"Time distribution for Scintillator of %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
302 hn = Form(
"h2_timeF%d_rpc", iF);
303 ht = Form(
"Time distribution for RPC of %s; Sector Index; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
305 hn = Form(
"h2_timeF%d_scint", iF);
306 ht = Form(
"Time distribution for Scintillator of %s; Sector Index; T_rec-T_0-T_fly-T_propagation[ns]", iFstring[iF].Data());
309 hn = Form(
"h2_timeF%d_scint_end", iF);
310 ht = Form(
"Time distribution for Scintillator of %s (Endcap); Sector Index; T_rec-T_0-T_fly-T_propagation[ns]",
311 iFstring[iF].Data());
315 hn = Form(
"hc_timeF%d_rpc", iF);
316 ht = Form(
"Calibrated time distribution for RPC of %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iFstring[iF].Data());
318 hn = Form(
"hc_timeF%d_scint", iF);
319 ht = Form(
"Calibrated time distribution for Scintillator of %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
320 iFstring[iF].Data());
323 hn = Form(
"hc_timeF%d_scint_end", iF);
324 ht = Form(
"Calibrated time distribution for Scintillator of %s (Endcap); T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
325 iFstring[iF].Data());
329 hn = Form(
"h2c_timeF%d_rpc", iF);
330 ht = Form(
"Calibrated time distribution for RPC of %s; Sector Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
331 iFstring[iF].Data());
334 hn = Form(
"h2c_timeF%d_scint", iF);
335 ht = Form(
"Calibrated time distribution for Scintillator of %s; Sector Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
336 iFstring[iF].Data());
339 hn = Form(
"h2c_timeF%d_scint_end", iF);
340 ht = Form(
"Calibrated time distribution for Scintillator of %s (Endcap) ; Sector Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
341 iFstring[iF].Data());
345 for (
int iS = 0; iS < 8; ++iS) {
347 hn = Form(
"h_timeF%d_S%d_scint", iF, iS);
348 ht = Form(
"Time distribution for Scintillator of Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
351 hn = Form(
"h_timeF%d_S%d_rpc", iF, iS);
352 ht = Form(
"Time distribution for RPC of Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
354 hn = Form(
"h2_timeF%d_S%d", iF, iS);
355 ht = Form(
"Time distribution of Sector%d, %s; Layer Index; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
359 hn = Form(
"hc_timeF%d_S%d_scint", iF, iS);
360 ht = Form(
"Calibrated time distribution for Scintillator of Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iS,
361 iFstring[iF].Data());
364 hn = Form(
"hc_timeF%d_S%d_rpc", iF, iS);
365 ht = Form(
"Calibrated time distribution for RPC of Sector%d, %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iS,
366 iFstring[iF].Data());
369 hn = Form(
"h2c_timeF%d_S%d", iF, iS);
370 ht = Form(
"Calibrated time distribution of Sector%d, %s; Layer Index; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]", iS,
371 iFstring[iF].Data());
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 distribution for Scintillator of Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iL, iS,
379 iFstring[iF].Data());
382 hn = Form(
"hc_timeF%d_S%d_L%d", iF, iS, iL);
383 ht = Form(
"Calibrated time distribution 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());
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 distribution 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());
394 hn = Form(
"h2_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
395 ht = Form(
"Time distribution 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());
400 hn = Form(
"hc_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
401 ht = Form(
"Calibrated time distribution 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());
405 hn = Form(
"h2c_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
406 ht = Form(
"Calibrated time distribution 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());
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 distribution 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());
419 hn = Form(
"h_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
420 ht = Form(
"time distribution 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());
425 hn = Form(
"hc_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
426 ht = Form(
"Calibrated time distribution 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());
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;
434 "Time versus propagation length; "
435 "propagation distance[cm]; "
436 "T_rec-T_0-T_fly-'T_calibration'[ns]",
437 200, 0.0, stripLength,
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 distribution for RPC of Layer%d, Sector%d, %s; T_rec-T_0-T_fly-T_propagation[ns]", iL, iS, iFstring[iF].Data());
449 hn = Form(
"hc_timeF%d_S%d_L%d", iF, iS, iL);
450 ht = Form(
"Calibrated time distribution 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());
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 distribution 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());
460 hn = Form(
"h2_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
461 ht = Form(
"time distribution 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());
465 hn = Form(
"hc_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
466 ht = Form(
"Calibrated time distribution 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());
471 hn = Form(
"h2c_timeF%d_S%d_L%d_P%d", iF, iS, iL, iP);
472 ht = Form(
"Calibrated time distribution 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());
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 distribution 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());
484 hn = Form(
"h_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
485 ht = Form(
"Time distribution 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());
489 hn = Form(
"hc_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
490 ht = Form(
"Calibrated time distribution 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());
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 distribution for Scintillator of Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]", iS,
503 iFstring[iF].Data());
506 hn = Form(
"h2_timeF%d_S%d_end", iF, iS);
507 ht = Form(
"Time distribution of Sector%d, %s (Endcap); Layer Index; T_rec-T_0-T_fly-T_propagation[ns]", iS, iFstring[iF].Data());
510 hn = Form(
"hc_timeF%d_S%d_scint_end", iF, iS);
511 ht = Form(
"Calibrated time distribution for Scintillator of Sector%d (Endcap), %s; T_rec-T_0-T_fly-T_propagation-T_calibration[ns]",
512 iS, iFstring[iF].Data());
515 hn = Form(
"h2c_timeF%d_S%d_end", iF, iS);
516 ht = Form(
"Calibrated time distribution 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,
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 distribution for Scintillator of Layer%d, Sector%d, %s (Endcap); T_rec-T_0-T_fly-T_propagation[ns]", iL, iS,
525 iFstring[iF].Data());
528 hn = Form(
"hc_timeF%d_S%d_L%d_end", iF, iS, iL);
529 ht = Form(
"Calibrated time distribution 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());
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 distribution 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());
541 hn = Form(
"h2_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
542 ht = Form(
"Time distribution 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());
547 hn = Form(
"hc_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
548 ht = Form(
"Calibrated time distribution 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());
553 hn = Form(
"h2c_timeF%d_S%d_L%d_P%d_end", iF, iS, iL, iP);
554 ht = Form(
"Calibrated time distribution 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,
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 distribution 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());
567 hn = Form(
"h_timeF%d_S%d_L%d_P%d_C%d_end", iF, iS, iL, iP, iC);
568 ht = Form(
"Time distribution 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());
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 distribution 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());
577 hn = Form(
"time_length_eklm_F%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
582 "Time versus propagation length; "
583 "propagation distance[cm]; "
584 "T_rec-T_0-T_fly-'T_calibration'[ns]",
585 200, 0.0, stripLength,
656 const std::vector< std::pair<KLMChannelNumber, unsigned int> >& channels,
657 double& delay,
double& delayError)
659 std::vector<struct Event> eventsChannel;
661 int nConvergedFits = 0;
664 if (nFits > (
int)channels.size())
665 nFits = channels.size();
666 for (
int i = 0; i < nFits; ++i) {
667 int subdetector, section, sector, layer, plane, strip;
669 channels[i].first, &subdetector, §ion, §or, &layer, &plane,
676 s_StripLength =
module->getStripLength(plane, strip);
683 for (
int j = 0; j < c_NBinsTime; ++j) {
684 for (
int k = 0; k < c_NBinsDistance; ++k)
685 s_BinnedData[j][k] = 0;
687 eventsChannel =
m_evts[channels[i].first];
688 double averageTime = 0;
689 for (
const Event& event : eventsChannel) {
690 double timeHit =
event.time();
692 timeHit = timeHit -
event.t0;
693 averageTime = averageTime + timeHit;
694 int timeBin = std::floor((timeHit - s_LowerTimeBoundary) * c_NBinsTime /
695 (s_UpperTimeBoundary - s_LowerTimeBoundary));
696 if (timeBin < 0 || timeBin >= c_NBinsTime)
698 int distanceBin = std::floor(event.dist * c_NBinsDistance / s_StripLength);
699 if (distanceBin < 0 || distanceBin >= c_NBinsDistance) {
700 B2ERROR(
"The distance to SiPM is greater than the strip length.");
703 s_BinnedData[timeBin][distanceBin] += 1;
705 averageTime = averageTime / eventsChannel.size();
707 minuit.SetPrintLevel(-1);
709 minuit.mnparm(0,
"P0", 1, 0.001, 0, 0, minuitResult);
710 minuit.mnparm(1,
"N", 10, 0.001, 0, 0, minuitResult);
711 minuit.mnparm(2,
"T0", averageTime, 0.001, 0, 0, minuitResult);
712 minuit.mnparm(3,
"SIGMA", 10, 0.001, 0, 0, minuitResult);
713 minuit.mnparm(4,
"DELAY", 0.0, 0.001, 0, 0, minuitResult);
715 minuit.mncomd(
"FIX 2 3 4 5", minuitResult);
716 minuit.mncomd(
"MIGRAD 10000", minuitResult);
717 minuit.mncomd(
"RELEASE 2", minuitResult);
718 minuit.mncomd(
"MIGRAD 10000", minuitResult);
719 minuit.mncomd(
"RELEASE 3", minuitResult);
720 minuit.mncomd(
"MIGRAD 10000", minuitResult);
721 minuit.mncomd(
"RELEASE 4", minuitResult);
722 minuit.mncomd(
"MIGRAD 10000", minuitResult);
723 minuit.mncomd(
"RELEASE 5", minuitResult);
724 minuit.mncomd(
"MIGRAD 10000", minuitResult);
726 if (minuit.fISW[1] != 3)
729 double channelDelay, channelDelayError;
730 minuit.GetParameter(4, channelDelay, channelDelayError);
731 delay = delay + channelDelay;
732 delayError = delayError + channelDelayError * channelDelayError;
734 delay = delay / nConvergedFits;
735 delayError =
sqrt(delayError) / (nConvergedFits - 1);
741 gROOT->SetBatch(kTRUE);
747 fcn_gaus =
new TF1(
"fcn_gaus",
"gaus");
748 fcn_land =
new TF1(
"fcn_land",
"landau");
749 fcn_pol1 =
new TF1(
"fcn_pol1",
"pol1");
750 fcn_const =
new TF1(
"fcn_const",
"pol0");
757 std::string name =
"time_calibration.root";
761 if (stat(name.c_str(), &buffer) != 0)
763 name =
"time_calibration_" + std::to_string(i) +
".root";
769 m_outFile =
new TFile(name.c_str(),
"recreate");
772 std::vector<struct Event> eventsChannel;
774 eventsChannel.clear();
779 std::vector< std::pair<KLMChannelNumber, unsigned int> > channelsBKLM;
780 std::vector< std::pair<KLMChannelNumber, unsigned int> > channelsEKLM;
784 m_cFlag[channel] = ChannelCalibrationStatus::c_NotEnoughData;
787 int nEvents =
m_evts[channel].size();
789 B2WARNING(
"Not enough calibration data collected."
790 <<
LogVar(
"channel", channel)
791 <<
LogVar(
"number of digit", nEvents));
794 m_cFlag[channel] = ChannelCalibrationStatus::c_FailedFit;
797 channelsBKLM.push_back(
798 std::pair<KLMChannelNumber, unsigned int>(channel, nEvents));
801 channelsEKLM.push_back(
802 std::pair<KLMChannelNumber, unsigned int>(channel, nEvents));
805 std::sort(channelsBKLM.begin(), channelsBKLM.end(), compareEventNumber);
806 std::sort(channelsEKLM.begin(), channelsEKLM.end(), compareEventNumber);
809 double delayBKLM, delayBKLMError;
810 double delayEKLM, delayEKLMError;
818 B2INFO(
"Effective light speed Estimation.");
820 channelId = klmChannel.getKLMChannelNumber();
821 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
824 eventsChannel =
m_evts[channelId];
825 int iSub = klmChannel.getSubdetector();
827 for (
const Event& event : eventsChannel) {
828 XYZVector diffD = XYZVector(event.diffDistX, event.diffDistY, event.diffDistZ);
830 double timeHit =
event.time();
832 timeHit = timeHit -
event.t0;
834 int iF = klmChannel.getSection();
835 int iS = klmChannel.getSector() - 1;
836 int iL = klmChannel.getLayer() - 1;
837 int iP = klmChannel.getPlane();
838 int iC = klmChannel.getStrip() - 1;
846 int iF = klmChannel.getSection() - 1;
847 int iS = klmChannel.getSector() - 1;
848 int iL = klmChannel.getLayer() - 1;
849 int iP = klmChannel.getPlane() - 1;
850 int iC = klmChannel.getStrip() - 1;
856 B2INFO(
"Effective light speed Estimation! Hists and Graph filling done.");
864 B2INFO(
"Global Mean for Raw." <<
LogVar(
"RPC", tmpMean_rpc_global) <<
LogVar(
"Scint BKLM",
865 tmpMean_scint_global) <<
LogVar(
"Scint EKLM", tmpMean_scint_global_end));
868 channelId = klmChannel.getKLMChannelNumber();
869 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
872 int iSub = klmChannel.getSubdetector();
874 int iF = klmChannel.getSection();
875 int iS = klmChannel.getSector() - 1;
876 int iL = klmChannel.getLayer() - 1;
877 int iP = klmChannel.getPlane();
878 int iC = klmChannel.getStrip() - 1;
880 double tmpMean_channel =
fcn_gaus->GetParameter(1);
882 m_timeShift[channelId] = tmpMean_channel - tmpMean_rpc_global;
884 m_timeShift[channelId] = tmpMean_channel - tmpMean_scint_global;
887 int iF = klmChannel.getSection() - 1;
888 int iS = klmChannel.getSector() - 1;
889 int iL = klmChannel.getLayer() - 1;
890 int iP = klmChannel.getPlane() - 1;
891 int iC = klmChannel.getStrip() - 1;
893 double tmpMean_channel =
fcn_gaus->GetParameter(1);
894 m_timeShift[channelId] = tmpMean_channel - tmpMean_scint_global_end;
901 B2INFO(
"Effective Light m_timeShift obtained. done.");
908 B2INFO(
"Effective light speed fitting.");
910 double delayRPCPhi =
fcn_pol1->GetParameter(1);
911 double e_slope_rpc_phi =
fcn_pol1->GetParError(1);
914 double delayRPCZ =
fcn_pol1->GetParameter(1);
915 double e_slope_rpc_z =
fcn_pol1->GetParError(1);
918 double slope_scint_phi =
fcn_pol1->GetParameter(1);
919 double e_slope_scint_phi =
fcn_pol1->GetParError(1);
922 double slope_scint_z =
fcn_pol1->GetParameter(1);
923 double e_slope_scint_z =
fcn_pol1->GetParError(1);
926 double slope_scint_plane1_end =
fcn_pol1->GetParameter(1);
927 double e_slope_scint_plane1_end =
fcn_pol1->GetParError(1);
930 double slope_scint_plane2_end =
fcn_pol1->GetParameter(1);
931 double e_slope_scint_plane2_end =
fcn_pol1->GetParError(1);
933 TString logStr_phi, logStr_z;
934 logStr_phi = Form(
"%.4f +/- %.4f ns/cm", delayRPCPhi, e_slope_rpc_phi);
935 logStr_z = Form(
"%.4f +/- %.4f ns/cm", delayRPCZ, e_slope_rpc_z);
936 B2INFO(
"Delay in RPCs:"
937 <<
LogVar(
"Fitted Value (phi readout) ", logStr_phi.Data())
938 <<
LogVar(
"Fitted Value (z readout) ", logStr_z.Data()));
939 logStr_phi = Form(
"%.4f +/- %.4f ns/cm", slope_scint_phi, e_slope_scint_phi);
940 logStr_z = Form(
"%.4f +/- %.4f ns/cm", slope_scint_z, e_slope_scint_z);
941 B2INFO(
"Delay in BKLM scintillators:"
942 <<
LogVar(
"Fitted Value (phi readout) ", logStr_phi.Data())
943 <<
LogVar(
"Fitted Value (z readout) ", logStr_z.Data()));
944 logStr_phi = Form(
"%.4f +/- %.4f ns/cm", slope_scint_plane1_end,
945 e_slope_scint_plane1_end);
946 logStr_z = Form(
"%.4f +/- %.4f ns/cm", slope_scint_plane2_end,
947 e_slope_scint_plane2_end);
948 B2INFO(
"Delay in EKLM scintillators:"
949 <<
LogVar(
"Fitted Value (plane1 readout) ", logStr_phi.Data())
950 <<
LogVar(
"Fitted Value (plane2 readout) ", logStr_z.Data()));
952 logStr_z = Form(
"%.4f +/- %.4f ns/cm", delayBKLM, delayBKLMError);
953 B2INFO(
"Delay in BKLM scintillators:"
954 <<
LogVar(
"Fitted Value (2d fit) ", logStr_z.Data()));
955 logStr_z = Form(
"%.4f +/- %.4f ns/cm", delayEKLM, delayEKLMError);
956 B2INFO(
"Delay in EKLM scintillators:"
957 <<
LogVar(
"Fitted Value (2d fit) ", logStr_z.Data()));
969 B2INFO(
"Time distribution filling begins.");
971 channelId = klmChannel.getKLMChannelNumber();
972 int iSub = klmChannel.getSubdetector();
974 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
976 eventsChannel =
m_evts[channelId];
978 for (
const Event& event : eventsChannel) {
979 double timeHit =
event.time();
981 timeHit = timeHit -
event.t0;
983 int iF = klmChannel.getSection();
984 int iS = klmChannel.getSector() - 1;
985 int iL = klmChannel.getLayer() - 1;
986 int iP = klmChannel.getPlane();
987 int iC = klmChannel.getStrip() - 1;
991 propgationT =
event.dist * delayRPCZ;
993 propgationT =
event.dist * delayRPCPhi;
994 double time = timeHit - propgationT;
1005 double propgationT =
event.dist * delayBKLM;
1006 double time = timeHit - propgationT;
1018 int iF = klmChannel.getSection() - 1;
1019 int iS = klmChannel.getSector() - 1;
1020 int iL = klmChannel.getLayer() - 1;
1021 int iP = klmChannel.getPlane() - 1;
1022 int iC = klmChannel.getStrip() - 1;
1023 double propgationT =
event.dist * delayEKLM;
1024 double time = timeHit - propgationT;
1038 B2INFO(
"Original filling done.");
1040 int iChannel_rpc = 0;
1042 int iChannel_end = 0;
1044 channelId = klmChannel.getKLMChannelNumber();
1045 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1047 int iSub = klmChannel.getSubdetector();
1050 int iF = klmChannel.getSection();
1051 int iS = klmChannel.getSector() - 1;
1052 int iL = klmChannel.getLayer() - 1;
1053 int iP = klmChannel.getPlane();
1054 int iC = klmChannel.getStrip() - 1;
1060 m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1073 int iF = klmChannel.getSection() - 1;
1074 int iS = klmChannel.getSector() - 1;
1075 int iL = klmChannel.getLayer() - 1;
1076 int iP = klmChannel.getPlane() - 1;
1077 int iC = klmChannel.getStrip() - 1;
1083 m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1104 B2INFO(
"Channel's time distribution fitting done.");
1109 B2INFO(
"Calibrated channel's time distribution filling begins.");
1113 channelId = klmChannel.getKLMChannelNumber();
1118 int iSub = klmChannel.getSubdetector();
1120 int iL = klmChannel.getLayer() - 1;
1132 channelId = klmChannel.getKLMChannelNumber();
1137 B2DEBUG(20,
"Uncalibrated Estimation " <<
LogVar(
"Channel", channelId) <<
LogVar(
"Estimated value",
m_timeShift[channelId]));
1144 channelId = klmChannel.getKLMChannelNumber();
1146 B2ERROR(
"!!! Not All Channels Calibration Constant Set. Error Happened on " <<
LogVar(
"Channel", channelId));
1149 int iSub = klmChannel.getSubdetector();
1151 int iL = klmChannel.getLayer();
1171 channelId = klmChannel.getKLMChannelNumber();
1172 int iSub = klmChannel.getSubdetector();
1173 eventsChannel =
m_evts[channelId];
1174 for (
const Event& event : eventsChannel) {
1175 double timeHit =
event.time();
1177 timeHit = timeHit -
event.t0;
1179 int iF = klmChannel.getSection();
1180 int iS = klmChannel.getSector() - 1;
1181 int iL = klmChannel.getLayer() - 1;
1182 int iP = klmChannel.getPlane();
1183 int iC = klmChannel.getStrip() - 1;
1187 propgationT =
event.dist * delayRPCZ;
1189 propgationT =
event.dist * delayRPCPhi;
1190 double time = timeHit - propgationT -
m_timeShift[channelId];
1201 double propgationT =
event.dist * delayBKLM;
1202 double time = timeHit - propgationT -
m_timeShift[channelId];
1214 int iF = klmChannel.getSection() - 1;
1215 int iS = klmChannel.getSector() - 1;
1216 int iL = klmChannel.getLayer() - 1;
1217 int iP = klmChannel.getPlane() - 1;
1218 int iC = klmChannel.getStrip() - 1;
1219 double propgationT =
event.dist * delayEKLM;
1220 double time = timeHit - propgationT -
m_timeShift[channelId];
1234 int icChannel_rpc = 0;
1236 int icChannel_end = 0;
1238 channelId = klmChannel.getKLMChannelNumber();
1239 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1241 int iSub = klmChannel.getSubdetector();
1244 int iF = klmChannel.getSection();
1245 int iS = klmChannel.getSector() - 1;
1246 int iL = klmChannel.getLayer() - 1;
1247 int iP = klmChannel.getPlane();
1248 int iC = klmChannel.getStrip() - 1;
1254 m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1267 int iF = klmChannel.getSection() - 1;
1268 int iS = klmChannel.getSector() - 1;
1269 int iL = klmChannel.getLayer() - 1;
1270 int iP = klmChannel.getPlane() - 1;
1271 int iC = klmChannel.getStrip() - 1;
1277 m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1298 B2INFO(
"Channel's time distribution fitting done.");
1303 B2INFO(
"Calibrated channel's time distribution filling begins.");
1307 channelId = klmChannel.getKLMChannelNumber();
1312 int iSub = klmChannel.getSubdetector();
1314 int iL = klmChannel.getLayer() - 1;
1326 channelId = klmChannel.getKLMChannelNumber();
1331 B2DEBUG(20,
"Calibrated Estimation " <<
LogVar(
"Channel", channelId) <<
LogVar(
"Estimated value",
m_timeRes[channelId]));
1338 channelId = klmChannel.getKLMChannelNumber();
1340 B2ERROR(
"!!! Not All Channels Calibration Constant Set. Error Happened on " <<
LogVar(
"Channel", channelId));
1343 int iSub = klmChannel.getSubdetector();
1345 int iL = klmChannel.getLayer();
1376 B2INFO(
"Save Histograms into Files.");
1377 TDirectory* dir_monitor =
m_outFile->mkdir(
"monitor_Hists",
"",
true);
1381 h_diff->SetDirectory(dir_monitor);
1384 TDirectory* dir_effC =
m_outFile->mkdir(
"effC_Hists",
"",
true);
1400 TDirectory* dir_time =
m_outFile->mkdir(
"time",
"",
true);
1425 B2INFO(
"Top file setup Done.");
1427 TDirectory* dir_time_F[2];
1428 TDirectory* dir_time_FS[2][8];
1429 TDirectory* dir_time_FSL[2][8][15];
1430 TDirectory* dir_time_FSLP[2][8][15][2];
1431 TDirectory* dir_time_F_end[2];
1432 TDirectory* dir_time_FS_end[2][4];
1433 TDirectory* dir_time_FSL_end[2][4][14];
1434 TDirectory* dir_time_FSLP_end[2][4][14][2];
1436 B2INFO(
"Sub files declare Done.");
1437 for (
int iF = 0; iF < 2; ++iF) {
1456 sprintf(dirname,
"isForward_%d", iF);
1457 dir_time_F[iF] = dir_time->mkdir(dirname,
"",
true);
1458 dir_time_F[iF]->cd();
1460 for (
int iS = 0; iS < 8; ++iS) {
1467 h2_timeFS[iF][iS]->SetDirectory(dir_time_F[iF]);
1468 h2c_timeFS[iF][iS]->SetDirectory(dir_time_F[iF]);
1470 sprintf(dirname,
"Sector_%d", iS + 1);
1471 dir_time_FS[iF][iS] = dir_time_F[iF]->mkdir(dirname,
"",
true);
1472 dir_time_FS[iF][iS]->cd();
1474 for (
int iL = 0; iL < 15; ++iL) {
1475 h_timeFSL[iF][iS][iL]->SetDirectory(dir_time_FS[iF][iS]);
1476 hc_timeFSL[iF][iS][iL]->SetDirectory(dir_time_FS[iF][iS]);
1478 sprintf(dirname,
"Layer_%d", iL + 1);
1479 dir_time_FSL[iF][iS][iL] = dir_time_FS[iF][iS]->mkdir(dirname,
"",
true);
1480 dir_time_FSL[iF][iS][iL]->cd();
1481 for (
int iP = 0; iP < 2; ++iP) {
1482 h_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1483 hc_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1484 h2_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1485 h2c_timeFSLP[iF][iS][iL][iP]->SetDirectory(dir_time_FSL[iF][iS][iL]);
1487 sprintf(dirname,
"Plane_%d", iP);
1488 dir_time_FSLP[iF][iS][iL][iP] = dir_time_FSL[iF][iS][iL]->mkdir(dirname,
"",
true);
1489 dir_time_FSLP[iF][iS][iL][iP]->cd();
1492 for (
int iC = 0; iC < nchannel_max; ++iC) {
1495 h_timeFSLPC[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP[iF][iS][iL][iP]);
1496 hc_timeFSLPC[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP[iF][iS][iL][iP]);
1503 sprintf(dirname,
"isForward_%d_end", iF + 1);
1504 dir_time_F_end[iF] = dir_time->mkdir(dirname,
"",
true);
1505 dir_time_F_end[iF]->cd();
1506 int maxLayer = 12 + 2 * iF;
1507 for (
int iS = 0; iS < 4; ++iS) {
1514 sprintf(dirname,
"Sector_%d_end", iS + 1);
1515 dir_time_FS_end[iF][iS] = dir_time_F_end[iF]->mkdir(dirname,
"",
true);
1516 dir_time_FS_end[iF][iS]->cd();
1517 for (
int iL = 0; iL < maxLayer; ++iL) {
1518 h_timeFSL_end[iF][iS][iL]->SetDirectory(dir_time_FS_end[iF][iS]);
1519 hc_timeFSL_end[iF][iS][iL]->SetDirectory(dir_time_FS_end[iF][iS]);
1521 sprintf(dirname,
"Layer_%d_end", iL + 1);
1522 dir_time_FSL_end[iF][iS][iL] = dir_time_FS_end[iF][iS]->mkdir(dirname,
"",
true);
1523 dir_time_FSL_end[iF][iS][iL]->cd();
1524 for (
int iP = 0; iP < 2; ++iP) {
1525 h_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1526 hc_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1527 h2_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1528 h2c_timeFSLP_end[iF][iS][iL][iP]->SetDirectory(dir_time_FSL_end[iF][iS][iL]);
1530 sprintf(dirname,
"plane_%d_end", iP);
1531 dir_time_FSLP_end[iF][iS][iL][iP] = dir_time_FSL_end[iF][iS][iL]->mkdir(dirname,
"",
true);
1532 dir_time_FSLP_end[iF][iS][iL][iP]->cd();
1534 for (
int iC = 0; iC < 75; ++iC) {
1535 h_timeFSLPC_end[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP_end[iF][iS][iL][iP]);
1536 hc_timeFSLPC_end[iF][iS][iL][iP][iC]->SetDirectory(dir_time_FSLP_end[iF][iS][iL][iP]);
1547 B2INFO(
"File Write and Close. Done.");