Run algorithm on data.
1135{
1136 int channelId;
1137 gROOT->SetBatch(kTRUE);
1143
1144 fcn_gaus =
new TF1(
"fcn_gaus",
"gaus");
1145 fcn_land =
new TF1(
"fcn_land",
"landau");
1146 fcn_pol1 =
new TF1(
"fcn_pol1",
"pol1");
1147 fcn_const =
new TF1(
"fcn_const",
"pol0");
1148
1149
1152 return result;
1153
1154
1155 std::string name = "time_calibration.root";
1156 int i = 1;
1157 while (1) {
1158 struct stat buffer;
1159 if (stat(name.c_str(), &buffer) != 0)
1160 break;
1161 name = "time_calibration_" + std::to_string(i) + ".root";
1162 i = i + 1;
1163 if (i < 0)
1164 break;
1165 }
1166 m_outFile =
new TFile(name.c_str(),
"recreate");
1168
1169 std::vector<struct Event> eventsChannel;
1170 eventsChannel.clear();
1173
1174 B2INFO("Counting events per channel...");
1175 std::map<KLMChannelNumber, unsigned int> eventCounts;
1177
1178
1179 std::vector< std::pair<KLMChannelNumber, unsigned int> > channelsBKLM;
1180 std::vector< std::pair<KLMChannelNumber, unsigned int> > channelsEKLM;
1181 KLMChannelIndex klmChannels;
1182
1183 for (const KLMChannelIndex& klmChannel : klmChannels) {
1185 m_cFlag[channel] = ChannelCalibrationStatus::c_NotEnoughData;
1186
1187 if (eventCounts.find(channel) == eventCounts.end())
1188 continue;
1189
1190 int nEvents = eventCounts[channel];
1192 B2WARNING("Not enough calibration data collected."
1193 << LogVar("channel", channel)
1194 << LogVar("number of digit", nEvents));
1195 continue;
1196 }
1197
1198 m_cFlag[channel] = ChannelCalibrationStatus::c_FailedFit;
1199
1202 channelsBKLM.push_back(std::pair<KLMChannelNumber, unsigned int>(channel, nEvents));
1203 }
1205 channelsEKLM.push_back(std::pair<KLMChannelNumber, unsigned int>(channel, nEvents));
1206 }
1207 }
1208
1209 std::sort(channelsBKLM.begin(), channelsBKLM.end(), compareEventNumber);
1210 std::sort(channelsEKLM.begin(), channelsEKLM.end(), compareEventNumber);
1211
1212
1213 double delayBKLM, delayBKLMError;
1214 double delayEKLM, delayEKLMError;
1215
1216
1221
1222 B2INFO("2D fits complete, data cleared.");
1223
1224
1225 auto isRPCBackward = [](const KLMChannelIndex & ch) {
1229 };
1230
1231 auto isRPCForward = [](const KLMChannelIndex & ch) {
1235 };
1236
1237 auto isBKLMScintillatorBackward = [](const KLMChannelIndex & ch) {
1241 };
1242
1243 auto isBKLMScintillatorForward = [](const KLMChannelIndex & ch) {
1247 };
1248
1249 auto isEKLMScintillatorBackward = [](const KLMChannelIndex & ch) {
1252 };
1253
1254 auto isEKLMScintillatorForward = [](const KLMChannelIndex & ch) {
1257 };
1258
1259 std::vector<std::pair<std::string, std::function<bool(const KLMChannelIndex&)>>> batches = {
1260 {"RPC Backward", isRPCBackward},
1261 {"RPC Forward", isRPCForward},
1262 {"BKLM Scintillator Backward", isBKLMScintillatorBackward},
1263 {"BKLM Scintillator Forward", isBKLMScintillatorForward},
1264 {"EKLM Scintillator Backward", isEKLMScintillatorBackward},
1265 {"EKLM Scintillator Forward", isEKLMScintillatorForward}
1266 };
1267
1268
1269
1270
1271
1272 B2INFO("First loop: Computing global statistics (batched processing)...");
1273
1274 TString iFstring[2] = {"Backward", "Forward"};
1275 TString iPstring[2] = {"ZReadout", "PhiReadout"};
1276 int nBin = 80;
1277 int nBin_scint = 80;
1278
1279 for (const auto& batch : batches) {
1280 B2INFO("Processing batch for global stats: " << batch.first);
1282
1284 channelId = klmChannel.getKLMChannelNumber();
1285
1286 if (!batch.second(klmChannel))
1287 continue;
1288
1289 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1290 continue;
1291
1293 continue;
1294
1295 eventsChannel =
m_evts[channelId];
1296 int iSub = klmChannel.getSubdetector();
1298
1299
1300 for (
const Event& event : eventsChannel) {
1301
1303 continue;
1304
1305 XYZVector diffD = XYZVector(event.diffDistX, event.diffDistY, event.diffDistZ);
1307
1308 double timeHit = event.time();
1310 timeHit = timeHit - event.t0;
1311
1312 if (timeHit <= -400e3)
1313 continue;
1314
1316
1319 } else {
1321 }
1322 } else {
1324 }
1325 }
1326 }
1327
1329 B2INFO("Batch processed and cleared: " << batch.first);
1330 }
1331
1332
1337
1338 B2INFO("Global Mean for Raw." << LogVar("RPC", tmpMean_rpc_global)
1339 << LogVar("Scint BKLM", tmpMean_scint_global)
1340 << LogVar("Scint EKLM", tmpMean_scint_global_end));
1341
1342
1343
1344
1345
1346 B2INFO("Second pass: Computing per-channel time shifts (batched processing)...");
1347
1348 for (const auto& batch : batches) {
1349 B2INFO("Processing batch for time shifts: " << batch.first);
1351
1353 channelId = klmChannel.getKLMChannelNumber();
1354
1355 if (!batch.second(klmChannel))
1356 continue;
1357
1358 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1359 continue;
1360
1362 continue;
1363
1364 eventsChannel =
m_evts[channelId];
1365 int iSub = klmChannel.getSubdetector();
1366 int iF, iS, iL, iP, iC;
1367
1369 iF = klmChannel.getSection();
1370 iS = klmChannel.getSector() - 1;
1371 iL = klmChannel.getLayer() - 1;
1372 iP = klmChannel.getPlane();
1373 iC = klmChannel.getStrip() - 1;
1374 } else {
1375 iF = klmChannel.getSection() - 1;
1376 iS = klmChannel.getSector() - 1;
1377 iL = klmChannel.getLayer() - 1;
1378 iP = klmChannel.getPlane() - 1;
1379 iC = klmChannel.getStrip() - 1;
1380 }
1381
1382
1383 TString hn, ht;
1384 TH1F* h_temp_tc = nullptr;
1385
1387
1389 hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_tc", iF, iS, iL, iP, iC);
1390 ht = Form("Time distribution for RPC of Channel%d, %s, Layer%d, Sector%d, %s",
1391 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1393 } else {
1394 hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_tc", iF, iS, iL, iP, iC);
1395 ht = Form("time distribution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s",
1396 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1399 }
1400 } else {
1401 hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_tc_end", iF, iS, iL, iP, iC);
1402 ht = Form("Time distribution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s (Endcap)",
1403 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1406 }
1407
1408 for (
const Event& event : eventsChannel) {
1409
1411 continue;
1412
1413 double timeHit = event.time();
1415 timeHit = timeHit - event.t0;
1416 if (timeHit <= -400e3)
1417 continue;
1418 h_temp_tc->Fill(timeHit);
1419 }
1420
1422 double tmpMean_channel =
fcn_gaus->GetParameter(1);
1423
1425
1427 m_timeShift[channelId] = tmpMean_channel - tmpMean_rpc_global;
1428 } else {
1429 m_timeShift[channelId] = tmpMean_channel - tmpMean_scint_global;
1430 }
1431 } else {
1432 m_timeShift[channelId] = tmpMean_channel - tmpMean_scint_global_end;
1433 }
1434
1435 delete h_temp_tc;
1436 }
1437
1439 B2INFO("Batch processed and cleared: " << batch.first);
1440 }
1441
1445 B2INFO("Effective Light m_timeShift obtained.");
1446
1447
1452
1453 B2INFO("Effective light speed fitting.");
1454
1455
1457 double fittedDelayRPCPhi =
fcn_pol1->GetParameter(1);
1458 double e_slope_rpc_phi =
fcn_pol1->GetParError(1);
1459
1461 double fittedDelayRPCZ =
fcn_pol1->GetParameter(1);
1462 double e_slope_rpc_z =
fcn_pol1->GetParError(1);
1463
1464
1465 double delayRPCPhi, delayRPCZ;
1469 B2INFO(
"Using fixed RPC propagation delay: " <<
m_fixedRPCDelay <<
" ns/cm (c_eff = 0.5c)"
1470 << LogVar("Fitted phi (not used)", fittedDelayRPCPhi)
1471 << LogVar("Fitted Z (not used)", fittedDelayRPCZ));
1472 } else {
1473 delayRPCPhi = fittedDelayRPCPhi;
1474 delayRPCZ = fittedDelayRPCZ;
1475 }
1476
1478 double slope_scint_phi =
fcn_pol1->GetParameter(1);
1479 double e_slope_scint_phi =
fcn_pol1->GetParError(1);
1480
1482 double slope_scint_z =
fcn_pol1->GetParameter(1);
1483 double e_slope_scint_z =
fcn_pol1->GetParError(1);
1484
1486 double slope_scint_plane1_end =
fcn_pol1->GetParameter(1);
1487 double e_slope_scint_plane1_end =
fcn_pol1->GetParError(1);
1488
1490 double slope_scint_plane2_end =
fcn_pol1->GetParameter(1);
1491 double e_slope_scint_plane2_end =
fcn_pol1->GetParError(1);
1492
1493 TString logStr_phi, logStr_z;
1495 logStr_phi = Form("%.4f ns/cm (fixed)", delayRPCPhi);
1496 logStr_z = Form("%.4f ns/cm (fixed)", delayRPCZ);
1497 B2INFO("Delay in RPCs (using fixed value):"
1498 << LogVar("Used Value (phi readout)", logStr_phi.Data())
1499 << LogVar("Used Value (z readout)", logStr_z.Data()));
1500 } else {
1501 logStr_phi = Form("%.4f +/- %.4f ns/cm", delayRPCPhi, e_slope_rpc_phi);
1502 logStr_z = Form("%.4f +/- %.4f ns/cm", delayRPCZ, e_slope_rpc_z);
1503 B2INFO("Delay in RPCs:"
1504 << LogVar("Fitted Value (phi readout)", logStr_phi.Data())
1505 << LogVar("Fitted Value (z readout)", logStr_z.Data()));
1506 }
1507 logStr_phi = Form("%.4f +/- %.4f ns/cm", slope_scint_phi, e_slope_scint_phi);
1508 logStr_z = Form("%.4f +/- %.4f ns/cm", slope_scint_z, e_slope_scint_z);
1509 B2INFO("Delay in BKLM scintillators:"
1510 << LogVar("Fitted Value (phi readout) ", logStr_phi.Data())
1511 << LogVar("Fitted Value (z readout) ", logStr_z.Data()));
1512 logStr_phi = Form("%.4f +/- %.4f ns/cm", slope_scint_plane1_end,
1513 e_slope_scint_plane1_end);
1514 logStr_z = Form("%.4f +/- %.4f ns/cm", slope_scint_plane2_end,
1515 e_slope_scint_plane2_end);
1516 B2INFO("Delay in EKLM scintillators:"
1517 << LogVar("Fitted Value (plane1 readout) ", logStr_phi.Data())
1518 << LogVar("Fitted Value (plane2 readout) ", logStr_z.Data()));
1519
1520 logStr_z = Form("%.4f +/- %.4f ns/cm", delayBKLM, delayBKLMError);
1521 B2INFO("Delay in BKLM scintillators:"
1522 << LogVar("Fitted Value (2d fit) ", logStr_z.Data()));
1523 logStr_z = Form("%.4f +/- %.4f ns/cm", delayEKLM, delayEKLMError);
1524 B2INFO("Delay in EKLM scintillators:"
1525 << LogVar("Fitted Value (2d fit) ", logStr_z.Data()));
1526
1531
1532
1533
1534
1535
1536 B2INFO("Third loop: Time distribution filling (batched processing)...");
1537
1538 for (const auto& batch : batches) {
1539 B2INFO("Processing batch: " << batch.first);
1541
1543 channelId = klmChannel.getKLMChannelNumber();
1544
1545 if (!batch.second(klmChannel))
1546 continue;
1547
1548 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData)
1549 continue;
1550
1552 continue;
1553
1554 eventsChannel =
m_evts[channelId];
1555 int iSub = klmChannel.getSubdetector();
1556 int iF, iS, iL, iP, iC;
1557
1559 iF = klmChannel.getSection();
1560 iS = klmChannel.getSector() - 1;
1561 iL = klmChannel.getLayer() - 1;
1562 iP = klmChannel.getPlane();
1563 iC = klmChannel.getStrip() - 1;
1564 } else {
1565 iF = klmChannel.getSection() - 1;
1566 iS = klmChannel.getSector() - 1;
1567 iL = klmChannel.getLayer() - 1;
1568 iP = klmChannel.getPlane() - 1;
1569 iC = klmChannel.getStrip() - 1;
1570 }
1571
1572
1573 TString hn, ht;
1574 TH1F* h_temp = nullptr;
1575
1577
1579 hn = Form("h_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
1580 ht = Form("Time distribution for RPC of Channel%d, %s, Layer%d, Sector%d, %s",
1581 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1583 } else {
1584 hn = Form("h_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
1585 ht = Form("time distribution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s",
1586 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1589 }
1590 } else {
1591 hn = Form("h_timeF%d_S%d_L%d_P%d_C%d_end", iF, iS, iL, iP, iC);
1592 ht = Form("Time distribution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s (Endcap)",
1593 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1596 }
1597
1598
1599 for (
const Event& event : eventsChannel) {
1600
1602 continue;
1603
1604 double timeHit = event.time();
1606 timeHit = timeHit - event.t0;
1607 if (timeHit <= -400e3)
1608 continue;
1609
1611
1613 double propgationT;
1615 propgationT = event.dist * delayRPCZ;
1616 else
1617 propgationT = event.dist * delayRPCPhi;
1618 double time = timeHit - propgationT;
1619
1621 h_temp->Fill(time);
1622
1631 }
1632 } else {
1633 double propgationT = event.dist * delayBKLM;
1634 double time = timeHit - propgationT;
1635
1637 h_temp->Fill(time);
1638
1647 }
1648 }
1649 } else {
1650 double propgationT = event.dist * delayEKLM;
1651 double time = timeHit - propgationT;
1652
1654 h_temp->Fill(time);
1655
1664 }
1665 }
1666 }
1667
1668 TFitResultPtr r = h_temp->Fit(
fcn_gaus,
"LESQ");
1669 if (int(r) == 0) {
1670 m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1673 }
1674
1675
1680 }
1681
1683 B2INFO("Batch processed and cleared: " << batch.first);
1684 }
1685
1686 B2INFO("Original filling done.");
1687
1688
1689 int iChannel_rpc = 0;
1690 int iChannel = 0;
1691 int iChannel_end = 0;
1693 channelId = klmChannel.getKLMChannelNumber();
1694 if (
m_cFlag[channelId] != ChannelCalibrationStatus::c_SuccessfulCalibration)
1695 continue;
1696
1697 int iSub = klmChannel.getSubdetector();
1699 int iL = klmChannel.getLayer() - 1;
1700
1704 iChannel_rpc++;
1705 } else {
1708 iChannel++;
1709 }
1710 } else {
1713 iChannel_end++;
1714 }
1715 }
1716
1720
1724
1728
1729 B2INFO("Channel's time distribution fitting done.");
1733
1734 B2INFO("Calibrated channel's time distribution filling begins.");
1735
1738 channelId = klmChannel.getKLMChannelNumber();
1741 continue;
1745 }
1746
1748 channelId = klmChannel.getKLMChannelNumber();
1750 continue;
1753 B2DEBUG(20,
"Uncalibrated Estimation " << LogVar(
"Channel", channelId) << LogVar(
"Estimated value",
m_timeShift[channelId]));
1754 }
1755
1756 iChannel_rpc = 0;
1757 iChannel = 0;
1758 iChannel_end = 0;
1760 channelId = klmChannel.getKLMChannelNumber();
1762 B2ERROR("!!! Not All Channels Calibration Constant Set. Error Happened on " << LogVar("Channel", channelId));
1763 continue;
1764 }
1765 int iSub = klmChannel.getSubdetector();
1767
1768 int iL = klmChannel.getLayer() - 1;
1771 iChannel_rpc++;
1772 } else {
1774 iChannel++;
1775 }
1776 } else {
1778 iChannel_end++;
1779 }
1780 }
1781
1782
1787
1788
1789
1790
1791
1792 B2INFO("Fourth loop: Calibrated time distribution filling (batched processing)...");
1793
1794 for (const auto& batch : batches) {
1795 B2INFO("Processing batch: " << batch.first);
1797
1799 channelId = klmChannel.getKLMChannelNumber();
1800
1801 if (!batch.second(klmChannel))
1802 continue;
1803
1805 continue;
1806
1807 eventsChannel =
m_evts[channelId];
1808 int iSub = klmChannel.getSubdetector();
1809 int iF, iS, iL, iP, iC;
1810
1812 iF = klmChannel.getSection();
1813 iS = klmChannel.getSector() - 1;
1814 iL = klmChannel.getLayer() - 1;
1815 iP = klmChannel.getPlane();
1816 iC = klmChannel.getStrip() - 1;
1817 } else {
1818 iF = klmChannel.getSection() - 1;
1819 iS = klmChannel.getSector() - 1;
1820 iL = klmChannel.getLayer() - 1;
1821 iP = klmChannel.getPlane() - 1;
1822 iC = klmChannel.getStrip() - 1;
1823 }
1824
1825 TString hn, ht;
1826 TH1F* hc_temp = nullptr;
1827
1829
1831 hn = Form("hc_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
1832 ht = Form("Calibrated time distribution for RPC of Channel%d, %s, Layer%d, Sector%d, %s",
1833 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1836 } else {
1837 hn = Form("hc_timeF%d_S%d_L%d_P%d_C%d", iF, iS, iL, iP, iC);
1838 ht = Form("Calibrated time distribution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s",
1839 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1842 }
1843 } else {
1844 hn = Form("hc_timeF%d_S%d_L%d_P%d_C%d_end", iF, iS, iL, iP, iC);
1845 ht = Form("Calibrated time distribution for Scintillator of Channel%d, %s, Layer%d, Sector%d, %s (Endcap)",
1846 iC, iPstring[iP].Data(), iL, iS, iFstring[iF].Data());
1849 }
1850
1851 for (
const Event& event : eventsChannel) {
1852
1854 continue;
1855
1856 double timeHit = event.time();
1858 timeHit = timeHit - event.t0;
1859 if (timeHit <= -400e3)
1860 continue;
1861
1863
1865 double propgationT;
1867 propgationT = event.dist * delayRPCZ;
1868 else
1869 propgationT = event.dist * delayRPCPhi;
1870 double time = timeHit - propgationT -
m_timeShift[channelId];
1871
1873 hc_temp->Fill(time);
1874
1883 }
1884 } else {
1885 double propgationT = event.dist * delayBKLM;
1886 double time = timeHit - propgationT -
m_timeShift[channelId];
1887
1889 hc_temp->Fill(time);
1890
1899 }
1900 }
1901 } else {
1902 double propgationT = event.dist * delayEKLM;
1903 double time = timeHit - propgationT -
m_timeShift[channelId];
1904
1906 hc_temp->Fill(time);
1907
1916 }
1917 }
1918 }
1919
1920 if (
m_cFlag[channelId] == ChannelCalibrationStatus::c_NotEnoughData) {
1921 delete hc_temp;
1922 continue;
1923 }
1924
1925 TFitResultPtr rc = hc_temp->Fit(
fcn_gaus,
"LESQ");
1926 if (int(rc) == 0) {
1927 m_cFlag[channelId] = ChannelCalibrationStatus::c_SuccessfulCalibration;
1930 }
1931
1932
1937 }
1938
1940 B2INFO("Batch processed and cleared: " << batch.first);
1941 }
1942
1943
1944 int icChannel_rpc = 0;
1945 int icChannel = 0;
1946 int icChannel_end = 0;
1948 channelId = klmChannel.getKLMChannelNumber();
1949 if (
m_cFlag[channelId] != ChannelCalibrationStatus::c_SuccessfulCalibration)
1950 continue;
1951
1952 int iSub = klmChannel.getSubdetector();
1954 int iL = klmChannel.getLayer() - 1;
1955
1959 icChannel_rpc++;
1960 } else {
1963 icChannel++;
1964 }
1965 } else {
1968 icChannel_end++;
1969 }
1970 }
1971
1975
1979
1983
1984 B2INFO("Channel's time distribution fitting done.");
1988
1989 B2INFO("Calibrated channel's time distribution filling begins.");
1990
1993 channelId = klmChannel.getKLMChannelNumber();
1996 continue;
2000 }
2001
2003 channelId = klmChannel.getKLMChannelNumber();
2005 continue;
2008 B2DEBUG(20,
"Calibrated Estimation " << LogVar(
"Channel", channelId) << LogVar(
"Estimated value",
m_timeRes[channelId]));
2009 }
2010
2011 icChannel_rpc = 0;
2012 icChannel = 0;
2013 icChannel_end = 0;
2015 channelId = klmChannel.getKLMChannelNumber();
2017 B2ERROR("!!! Not All Channels Calibration Constant Set. Error Happened on " << LogVar("Channel", channelId));
2018 continue;
2019 }
2020 int iSub = klmChannel.getSubdetector();
2022
2023 int iL = klmChannel.getLayer() - 1;
2026 icChannel_rpc++;
2027 } else {
2029 icChannel++;
2030 }
2031 } else {
2033 icChannel_end++;
2034 }
2035 }
2036
2037
2038
2039
2040 B2INFO("Fifth pass: Computing di-muon ΔT0 for EventT0 hit resolution calibration...");
2041
2042
2043 struct TrackT0Info {
2044 int charge;
2045 int nHits_BKLM_Scint;
2046 int nHits_BKLM_RPC_Phi;
2047 int nHits_BKLM_RPC_Z;
2048 int nHits_EKLM_Scint;
2049 int nHits_BKLM_RPC;
2050 double sumT0_BKLM_Scint;
2051 double sumT0_BKLM_RPC_Phi;
2052 double sumT0_BKLM_RPC_Z;
2053 double sumT0_EKLM_Scint;
2054 double sumT0_BKLM_RPC;
2055
2056 TrackT0Info() : charge(0),
2057 nHits_BKLM_Scint(0),
2058 nHits_BKLM_RPC_Phi(0),
2059 nHits_BKLM_RPC_Z(0),
2060 nHits_EKLM_Scint(0),
2061 nHits_BKLM_RPC(0),
2062 sumT0_BKLM_Scint(0.0),
2063 sumT0_BKLM_RPC_Phi(0.0),
2064 sumT0_BKLM_RPC_Z(0.0),
2065 sumT0_EKLM_Scint(0.0),
2066 sumT0_BKLM_RPC(0.0) {}
2067 };
2068
2069
2070 std::map<std::pair<int, int>, std::map<int, TrackT0Info>> eventTrackMap;
2071
2072
2073 for (const auto& batch : batches) {
2074 B2INFO("Processing batch for di-muon analysis: " << batch.first);
2076
2077 for (
const auto& channelPair :
m_evts) {
2079 const std::vector<Event>& chEvents = channelPair.second;
2080
2081
2082 int subdetector, section, sector, layer, plane, strip;
2084 chId, &subdetector, §ion, §or, &layer, &plane, &strip);
2085
2086
2087 int iSub = subdetector;
2088 int iL = layer - 1;
2089
2090 for (
const Event& event : chEvents) {
2091
2093 continue;
2094
2095
2096 std::pair<int, int> eventKey(event.Run, event.Events);
2097 int trackIdx = event.nTrack;
2098 int charge = event.Track_Charge;
2099
2100
2101 double timeHit =
event.time() -
m_timeShift[chId];
2102
2103 if (timeHit <= -400e3)
2104 continue;
2105
2106
2107 double propT = 0.0;
2110
2112 propT = event.dist * delayRPCZ;
2113 else
2114 propT = event.dist * delayRPCPhi;
2115 } else {
2116
2117 propT = event.dist * delayBKLM;
2118 }
2119 } else {
2120
2121 propT = event.dist * delayEKLM;
2122 }
2123
2124 double t0_estimate = timeHit - propT;
2125
2126
2127 TrackT0Info& trackInfo = eventTrackMap[eventKey][trackIdx];
2128 trackInfo.charge = charge;
2129
2132
2133 trackInfo.nHits_BKLM_RPC++;
2134 trackInfo.sumT0_BKLM_RPC += t0_estimate;
2135
2137 trackInfo.nHits_BKLM_RPC_Z++;
2138 trackInfo.sumT0_BKLM_RPC_Z += t0_estimate;
2139 } else {
2140 trackInfo.nHits_BKLM_RPC_Phi++;
2141 trackInfo.sumT0_BKLM_RPC_Phi += t0_estimate;
2142 }
2143 } else {
2144
2145 trackInfo.nHits_BKLM_Scint++;
2146 trackInfo.sumT0_BKLM_Scint += t0_estimate;
2147 }
2148 } else {
2149 trackInfo.nHits_EKLM_Scint++;
2150 trackInfo.sumT0_EKLM_Scint += t0_estimate;
2151 }
2152 }
2153 }
2154
2156 }
2157
2158 B2INFO("Event-track map built. Processing events for EventT0 histograms...");
2159
2160
2161
2162 double sum_delta2_over_v_BKLM_Scint = 0.0;
2163 double sum_delta2_over_v_BKLM_RPC_Phi = 0.0;
2164 double sum_delta2_over_v_BKLM_RPC_Z = 0.0;
2165 double sum_delta2_over_v_EKLM_Scint = 0.0;
2166 double sum_delta2_over_v_BKLM_RPC = 0.0;
2167
2168 int nDimuon_BKLM_Scint = 0;
2169 int nDimuon_BKLM_RPC_Phi = 0;
2170 int nDimuon_BKLM_RPC_Z = 0;
2171 int nDimuon_EKLM_Scint = 0;
2172 int nDimuon_BKLM_RPC = 0;
2173
2174 for (const auto& eventPair : eventTrackMap) {
2175 const auto& trackMap = eventPair.second;
2176
2177
2178 if (trackMap.size() != 2)
2179 continue;
2180
2181 auto it1 = trackMap.begin();
2182 auto it2 = trackMap.begin();
2183 ++it2;
2184 const TrackT0Info& track1 = it1->second;
2185 const TrackT0Info& track2 = it2->second;
2186
2187
2188 if (track1.charge * track2.charge >= 0)
2189 continue;
2190
2191
2192 const TrackT0Info& muPlus = (track1.charge > 0) ? track1 : track2;
2193 const TrackT0Info& muMinus = (track1.charge > 0) ? track2 : track1;
2194
2195
2196 if (muPlus.nHits_BKLM_Scint > 0 && muMinus.nHits_BKLM_Scint > 0) {
2197 double t0_plus = muPlus.sumT0_BKLM_Scint / muPlus.nHits_BKLM_Scint;
2198 double t0_minus = muMinus.sumT0_BKLM_Scint / muMinus.nHits_BKLM_Scint;
2199 double deltaT0 = t0_plus - t0_minus;
2200
2201
2202 double v = 1.0 / muPlus.nHits_BKLM_Scint + 1.0 / muMinus.nHits_BKLM_Scint;
2203 int nTotal = muPlus.nHits_BKLM_Scint + muMinus.nHits_BKLM_Scint;
2204
2205 if (v <= 0.0)
2206 continue;
2207
2208
2209 sum_delta2_over_v_BKLM_Scint += (deltaT0 * deltaT0) / v;
2211 nDimuon_BKLM_Scint++;
2212
2213
2217
2218
2224
2225
2226 if (nTotal < 5) {
2228 } else if (nTotal < 15) {
2230 } else {
2232 }
2233 }
2234
2235
2236 if (muPlus.nHits_BKLM_RPC_Phi > 0 && muMinus.nHits_BKLM_RPC_Phi > 0) {
2237 double t0_plus = muPlus.sumT0_BKLM_RPC_Phi / muPlus.nHits_BKLM_RPC_Phi;
2238 double t0_minus = muMinus.sumT0_BKLM_RPC_Phi / muMinus.nHits_BKLM_RPC_Phi;
2239 double deltaT0 = t0_plus - t0_minus;
2240
2241 double v = 1.0 / muPlus.nHits_BKLM_RPC_Phi + 1.0 / muMinus.nHits_BKLM_RPC_Phi;
2242 int nTotal = muPlus.nHits_BKLM_RPC_Phi + muMinus.nHits_BKLM_RPC_Phi;
2243
2244 if (v <= 0.0)
2245 continue;
2246
2247 sum_delta2_over_v_BKLM_RPC_Phi += (deltaT0 * deltaT0) / v;
2249 nDimuon_BKLM_RPC_Phi++;
2250
2254
2255
2261
2262 if (nTotal < 10) {
2264 } else if (nTotal < 30) {
2266 } else {
2268 }
2269 }
2270
2271
2272 if (muPlus.nHits_BKLM_RPC_Z > 0 && muMinus.nHits_BKLM_RPC_Z > 0) {
2273 double t0_plus = muPlus.sumT0_BKLM_RPC_Z / muPlus.nHits_BKLM_RPC_Z;
2274 double t0_minus = muMinus.sumT0_BKLM_RPC_Z / muMinus.nHits_BKLM_RPC_Z;
2275 double deltaT0 = t0_plus - t0_minus;
2276
2277 double v = 1.0 / muPlus.nHits_BKLM_RPC_Z + 1.0 / muMinus.nHits_BKLM_RPC_Z;
2278 int nTotal = muPlus.nHits_BKLM_RPC_Z + muMinus.nHits_BKLM_RPC_Z;
2279
2280 if (v <= 0.0)
2281 continue;
2282
2283 sum_delta2_over_v_BKLM_RPC_Z += (deltaT0 * deltaT0) / v;
2285 nDimuon_BKLM_RPC_Z++;
2286
2290
2291
2297
2298 if (nTotal < 10) {
2300 } else if (nTotal < 30) {
2302 } else {
2304 }
2305 }
2306
2307
2308 if (muPlus.nHits_BKLM_RPC > 0 && muMinus.nHits_BKLM_RPC > 0) {
2309 double deltaT0 = muPlus.sumT0_BKLM_RPC / muPlus.nHits_BKLM_RPC
2310 - muMinus.sumT0_BKLM_RPC / muMinus.nHits_BKLM_RPC;
2311 double v = 1.0 / muPlus.nHits_BKLM_RPC + 1.0 / muMinus.nHits_BKLM_RPC;
2312 sum_delta2_over_v_BKLM_RPC += (deltaT0 * deltaT0) / v;
2314 nDimuon_BKLM_RPC++;
2315 }
2316
2317
2318 if (muPlus.nHits_EKLM_Scint > 0 && muMinus.nHits_EKLM_Scint > 0) {
2319 double t0_plus = muPlus.sumT0_EKLM_Scint / muPlus.nHits_EKLM_Scint;
2320 double t0_minus = muMinus.sumT0_EKLM_Scint / muMinus.nHits_EKLM_Scint;
2321 double deltaT0 = t0_plus - t0_minus;
2322
2323 double v = 1.0 / muPlus.nHits_EKLM_Scint + 1.0 / muMinus.nHits_EKLM_Scint;
2324 int nTotal = muPlus.nHits_EKLM_Scint + muMinus.nHits_EKLM_Scint;
2325
2326 if (v <= 0.0)
2327 continue;
2328
2329 sum_delta2_over_v_EKLM_Scint += (deltaT0 * deltaT0) / v;
2331 nDimuon_EKLM_Scint++;
2332
2336
2337
2343
2344
2345 if (nTotal < 5) {
2347 } else if (nTotal < 15) {
2349 } else {
2351 }
2352 }
2353 }
2354
2355 B2INFO("Di-muon ΔT0 data collected."
2356 << LogVar("BKLM Scint di-muon events", nDimuon_BKLM_Scint)
2357 << LogVar("BKLM RPC Phi di-muon events", nDimuon_BKLM_RPC_Phi)
2358 << LogVar("BKLM RPC (phi and z) di-muon events", nDimuon_BKLM_RPC)
2359 << LogVar("BKLM RPC Z di-muon events", nDimuon_BKLM_RPC_Z)
2360 << LogVar("EKLM Scint di-muon events", nDimuon_EKLM_Scint));
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379 auto coreSigma = [](TH1F * h, double fallback, double & err, int nEntries) -> double {
2380 err = (nEntries > 0) ? fallback / std::sqrt(2.0 * nEntries) : 0.0;
2381
2382 auto useMoment = [&](const char* reason)
2383 {
2384 B2WARNING("sigma_hit: " << reason << ", using the second moment"
2385 << LogVar("histogram", std::string((h != nullptr) ? h->GetName() : "none"))
2386 << LogVar("second moment [ns]", fallback));
2387 return fallback;
2388 };
2389 if (h == nullptr || h->GetEntries() < 100)
2390 return useMoment("fewer than 100 events");
2391 double prob[3] = {0.25, 0.50, 0.75};
2392 double q[3];
2393 h->GetQuantiles(3, q, prob);
2394 const double median = q[1];
2395 const double sigma0 = (q[2] - q[0]) / 1.349;
2396 if (!(sigma0 > 0.0))
2397 return useMoment("zero interquartile range");
2398 TF1 f("fcore", "gaus", median - 2.5 * sigma0, median + 2.5 * sigma0);
2399 f.SetParameters(h->GetMaximum(), median, sigma0);
2400 if (h->Fit(&f, "RQN0") != 0)
2401 return useMoment("core fit did not converge");
2402 const double sg = std::fabs(f.GetParameter(2));
2403
2404 if (sg > 5.0 * fallback || sg < 0.1 * fallback)
2405 return useMoment("core fit width far from the second moment");
2406 err = f.GetParError(2);
2407 return sg;
2408 };
2409
2410 float sigma_BKLM_Scint = 10.0f;
2411 float sigma_BKLM_Scint_err = 1.0f;
2412 if (nDimuon_BKLM_Scint > 0) {
2413 double sigma2 = sum_delta2_over_v_BKLM_Scint / static_cast<double>(nDimuon_BKLM_Scint);
2414 double moment = std::sqrt(sigma2);
2415 double fitErr = 0.0;
2416 double core = coreSigma(
h_z_bklm_scint, moment, fitErr, nDimuon_BKLM_Scint);
2417 sigma_BKLM_Scint = static_cast<float>(core);
2418 B2INFO("sigma_hit BKLM Scint" << LogVar("core fit [ns]", core)
2419 << LogVar("second moment [ns]", moment) << LogVar("ratio", moment / core));
2420 sigma_BKLM_Scint_err = static_cast<float>(fitErr);
2421 }
2422
2423 float sigma_RPC_Phi = 10.0f;
2424 float sigma_RPC_Phi_err = 1.0f;
2425 if (nDimuon_BKLM_RPC_Phi > 0) {
2426 double sigma2 = sum_delta2_over_v_BKLM_RPC_Phi / static_cast<double>(nDimuon_BKLM_RPC_Phi);
2427 double moment = std::sqrt(sigma2);
2428 double fitErr = 0.0;
2429 double core = coreSigma(
h_z_bklm_rpc_phi, moment, fitErr, nDimuon_BKLM_RPC_Phi);
2430 sigma_RPC_Phi = static_cast<float>(core);
2431 B2INFO("sigma_hit RPC phi" << LogVar("core fit [ns]", core)
2432 << LogVar("second moment [ns]", moment) << LogVar("ratio", moment / core));
2433 sigma_RPC_Phi_err = static_cast<float>(fitErr);
2434 }
2435
2436 float sigma_RPC_Z = 10.0f;
2437 float sigma_RPC_Z_err = 1.0f;
2438 if (nDimuon_BKLM_RPC_Z > 0) {
2439 double sigma2 = sum_delta2_over_v_BKLM_RPC_Z / static_cast<double>(nDimuon_BKLM_RPC_Z);
2440 double moment = std::sqrt(sigma2);
2441 double fitErr = 0.0;
2442 double core = coreSigma(
h_z_bklm_rpc_z, moment, fitErr, nDimuon_BKLM_RPC_Z);
2443 sigma_RPC_Z = static_cast<float>(core);
2444 B2INFO("sigma_hit RPC z" << LogVar("core fit [ns]", core)
2445 << LogVar("second moment [ns]", moment) << LogVar("ratio", moment / core));
2446 sigma_RPC_Z_err = static_cast<float>(fitErr);
2447 }
2448
2449
2450
2451 float sigma_RPC = 10.0f;
2452 float sigma_RPC_err = 1.0f;
2453 if (nDimuon_BKLM_RPC > 0) {
2454 double sigma2 = sum_delta2_over_v_BKLM_RPC / static_cast<double>(nDimuon_BKLM_RPC);
2455 double moment = std::sqrt(sigma2);
2456 double fitErr = 0.0;
2457 double core = coreSigma(
h_z_bklm_rpc, moment, fitErr, nDimuon_BKLM_RPC);
2458 sigma_RPC = static_cast<float>(core);
2459 B2INFO("sigma_hit RPC (phi and z together)" << LogVar("core fit [ns]", core)
2460 << LogVar("second moment [ns]", moment) << LogVar("ratio", moment / core));
2461 sigma_RPC_err = static_cast<float>(fitErr);
2462 }
2463
2464 float sigma_EKLM_Scint = 10.0f;
2465 float sigma_EKLM_Scint_err = 1.0f;
2466 if (nDimuon_EKLM_Scint > 0) {
2467 double sigma2 = sum_delta2_over_v_EKLM_Scint / static_cast<double>(nDimuon_EKLM_Scint);
2468 double moment = std::sqrt(sigma2);
2469 double fitErr = 0.0;
2470 double core = coreSigma(
h_z_eklm_scint, moment, fitErr, nDimuon_EKLM_Scint);
2471 sigma_EKLM_Scint = static_cast<float>(core);
2472 B2INFO("sigma_hit EKLM Scint" << LogVar("core fit [ns]", core)
2473 << LogVar("second moment [ns]", moment) << LogVar("ratio", moment / core));
2474 sigma_EKLM_Scint_err = static_cast<float>(fitErr);
2475 }
2476
2477 B2INFO("Extracted per-hit resolutions using event-by-event weighting:"
2478 << LogVar("σ_BKLM_Scint [ns]", sigma_BKLM_Scint) << LogVar("±", sigma_BKLM_Scint_err)
2479 << LogVar("σ_RPC_Phi [ns]", sigma_RPC_Phi) << LogVar("±", sigma_RPC_Phi_err)
2480 << LogVar("σ_RPC_Z [ns]", sigma_RPC_Z) << LogVar("±", sigma_RPC_Z_err)
2481 << LogVar("σ_RPC_combined [ns]", sigma_RPC) << LogVar("±", sigma_RPC_err)
2482 << LogVar("σ_EKLM_Scint [ns]", sigma_EKLM_Scint) << LogVar("±", sigma_EKLM_Scint_err));
2483
2484
2490
2491 B2INFO("EventT0 hit resolution calibration complete and stored in payload.");
2492
2493
2494 eventTrackMap.clear();
2495
2501
2503
2508
2510}
@ c_FirstRPCLayer
First RPC layer.
@ c_ForwardSection
Forward.
@ c_BackwardSection
Backward.
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_ForwardSection
Forward.
@ c_BackwardSection
Backward.
TProfile * m_Profile2EKLMScintillatorPlane2
For EKLM scintillator plane2.
TDirectory * m_channelHistDir_BKLM[2][8][15][2]
Directory structure for per-channel histograms (BKLM).
double mc_etime_channelAvg_rpc
Calibrated central value error of the global time distribution (BKLM RPC part).
TH2F * h2c_timeF_scint_end[2]
EKLM part.
KLMTimeResolution * m_timeResolution
DBObject of time resolution.
TH1F * hc_eventT0_scint_end_highN
Corrected EventT0 for EKLM scintillator with high hit count.
TProfile * prof_deltaT0_rms_vs_v_scint_end
DeltaT0 RMS vs inverse hit count profile for EKLM scintillator.
TH1F * h_time_scint_tc_end
EKLM part.
TH1F * h_z_bklm_rpc_z
Normalised di-muon time difference, BKLM RPC z readout.
TH1F * hc_eventT0_scint_midN
Corrected EventT0 for BKLM scintillator with medium hit count.
void createHistograms()
Create histograms.
TGraphErrors * gre_time_channel_scint
BKLM Scintillator.
TH1F * h_timeFSL[2][8][15]
BKLM part.
TH1F * hc_timeFSL_end[2][4][14]
EKLM part.
TH1F * h_timeFSLP_end[2][4][14][2]
EKLM part.
TGraph * gr_timeRes_channel_rpc
BKLM RPC.
TH1F * hc_timeFSLP_end[2][4][14][2]
EKLM part.
TH1F * h_timeFSLP[2][8][15][2]
BKLM part.
TH1F * hc_timeF_scint_end[2]
EKLM part.
std::map< KLMChannelNumber, double > m_timeShift
Shift values of each channel.
TH1F * h_time_scint
BKLM scintillator part.
double m_time_channelAvg_scint
Central value of the global time distribution (BKLM scintillator part).
TH1F * hc_timeFS_scint_end[2][4]
EKLM part.
double esti_timeRes(const KLMChannelIndex &klmChannel)
Estimate value of calibration constant for calibrated channels.
double m_UpperTimeBoundaryCalibratedRPC
Upper time boundary for RPC (calibrated data).
double m_ctime_channelAvg_rpc
Calibrated central value of the global time distribution (BKLM RPC part).
KLMTimeConstants * m_timeConstants
DBObject of time cost on some parts of the detector.
void setupDatabase()
Setup the database.
TH1F * h_z_bklm_rpc
Normalised di-muon time difference, BKLM RPC, phi and z readouts together.
TProfile * prof_deltaT0_rms_vs_v_rpc
DeltaT0 RMS vs inverse hit count profile for RPC.
TH1F * h_eventT0_scint
EventT0 seen by BKLM scintillator hits.
TH1F * hc_timeFS_scint[2][8]
BKLM scintillator part.
std::map< KLMChannelNumber, double > m_time_channel
Time distribution central value of each channel.
TH1F * hc_eventT0_rpc_lowN
Corrected EventT0 for RPC with low hit count.
double m_ctime_channelAvg_scint_end
Calibrated central value of the global time distribution (EKLM scintillator part).
CalibrationAlgorithm::EResult readCalibrationData()
Read calibration data.
TH1F * hc_eventT0_scint_lowN
Corrected EventT0 for BKLM scintillator with low hit count.
TGraph * gr_timeShift_channel_scint_end
EKLM.
TGraph * gr_timeRes_channel_scint
BKLM scintillator.
TH2F * h2_deltaT0_vs_v_scint
DeltaT0 vs inverse hit count for BKLM scintillator.
TH1F * hc_timeF_scint[2]
BKLM scintillator part.
TH1F * h_timeFS_scint[2][8]
BKLM scintillator part.
bool m_saveChannelHists
Write per-channel temporary histograms (tc/raw/hc) in minimal mode.
double m_UpperTimeBoundaryScintillatorsBKLM
Upper time boundary for BKLM scintillators.
void writeThenDelete_(TH1 *h, bool write, TDirectory *dir=nullptr)
Optionally write a histogram, then delete it to free memory.
TH1F * hc_timeF_rpc[2]
BKLM RPC part.
TH2F * h2c_timeFS_end[2][4]
EKLM part.
TGraphErrors * gre_ctime_channel_scint_end
EKLM.
TH1F * h_nHits_plus_scint_end
Number of EKLM scintillator hits per mu+ track.
TProfile * m_Profile2BKLMScintillatorPhi
For BKLM scintillator phi plane.
TH1F * hc_time_scint_end
EKLM part.
TH1F * hc_eventT0_scint
Corrected EventT0 for BKLM scintillator hits.
TGraphErrors * gre_time_channel_scint_end
EKLM.
TH2F * h2_timeFSLP[2][8][15][2]
BKLM part.
double m_UpperTimeBoundaryCalibratedScintillatorsEKLM
Upper time boundary for BKLM scintillators (calibrated data).
TGraph * gr_timeShift_channel_scint
BKLM scintillator.
double m_time_channelAvg_scint_end
Central value of the global time distribution (EKLM scintillator part).
TProfile * m_Profile2EKLMScintillatorPlane1
For EKLM scintillator plane1.
TH1F * hc_timeFS_rpc[2][8]
BKLM RPC part.
double m_UpperTimeBoundaryScintillatorsEKLM
Upper time boundary for BKLM scintillators.
void fillTimeDistanceProfiles(TProfile *profileRpcPhi, TProfile *profileRpcZ, TProfile *profileBKLMScintillatorPhi, TProfile *profileBKLMScintillatorZ, TProfile *profileEKLMScintillatorPlane1, TProfile *profileEKLMScintillatorPlane2, bool fill2dHistograms)
Fill profiles of time versus distance.
TFile * m_outFile
Output file.
double esti_timeShift(const KLMChannelIndex &klmChannel)
Estimate value of calibration constant for uncalibrated channels.
TH1F * h_z_bklm_scint
Normalised di-muon time difference z = dT0 / sqrt(1/N+ + 1/N-).
double m_LowerTimeBoundaryCalibratedScintillatorsEKLM
Lower time boundary for EKLM scintillators (calibrated data).
TH1F * hc_timeFSLP[2][8][15][2]
BKLM part.
TGraphErrors * gre_ctime_channel_rpc
BKLM RPC.
void saveHist()
Save histograms to file.
bool m_saveAllPlots
Default minimal unless you set true in your header script.
TH2F * h2c_timeFSLP[2][8][15][2]
BKLM part.
double m_ctime_channelAvg_scint
Calibrated central value of the global time distribution (BKLM scintillator part).
TF1 * fcn_const
Const function.
double m_UpperTimeBoundaryCalibratedScintillatorsBKLM
Upper time boundary for BKLM scintillators (calibrated data).
TProfile * m_Profile2RpcZ
For BKLM RPC z plane.
TH1F * h_nHits_minus_rpc
Number of RPC hits per mu- track.
TH2F * h2_timeFSLP_end[2][4][14][2]
EKLM part.
TH1I * hc_calibrated
Calibration statistics for each channel.
TH1F * h_eventT0_rpc
EventT0 seen by RPC hits.
void timeDistance2dFit(const std::vector< std::pair< KLMChannelNumber, unsigned int > > &channels, double &delay, double &delayError)
Two-dimensional fit for individual channels.
TGraph * gr_timeRes_channel_scint_end
EKLM.
TH1I * h_calibrated
Calibration statistics for each channel.
TProfile * m_ProfileBKLMScintillatorZ
For BKLM scintillator z plane.
double m_LowerTimeBoundaryScintillatorsBKLM
Lower time boundary for BKLM scintillators.
TH1F * hc_eventT0_scint_end_midN
Corrected EventT0 for EKLM scintillator with medium hit count.
TH1F * h_time_rpc_tc
BKLM RPC part.
TH1F * h_time_scint_end
EKLM part.
TH2F * h2c_timeF_scint[2]
BKLM scintillator part.
TF1 * fcn_pol1
Pol1 function.
double m_etime_channelAvg_scint_end
Central value error of the global time distribution (EKLM scintillator part).
TH1F * hc_time_rpc
BKLM RPC part.
double mc_etime_channelAvg_scint
Calibrated central value error of the global time distribution (BKLM scintillator part).
TH2F * h2c_timeFSLP_end[2][4][14][2]
EKLM part.
double mc_etime_channelAvg_scint_end
Calibrated central value error of the global time distribution (EKLM scintillator part).
bool passesADCCut(const Event &event, int subdetector, int layer) const
Check if event passes ADC count cuts for quality selection.
TH2F * h2_timeF_scint_end[2]
EKLM part.
KLMEventT0HitResolution * m_eventT0HitResolution
DBObject of per-hit time resolution for EventT0.
TF1 * fcn_gaus
Gaussian function.
double m_LowerTimeBoundaryCalibratedScintillatorsBKLM
Lower time boundary for BKLM scintillators (calibrated data).
TH1F * h_timeF_rpc[2]
BKLM RPC part.
TProfile * m_ProfileBKLMScintillatorPhi
For BKLM scintillator phi plane.
TH1F * hc_time_scint
BKLM scintillator part.
TH2F * h2_timeFS[2][8]
BKLM part.
double m_fixedRPCDelay
Fixed propagation delay for RPCs (ns/cm).
double m_etime_channelAvg_scint
Central value error of the global time distribution (BKLM scintillator part).
TH1F * h_timeF_scint_end[2]
EKLM part.
TH1F * hc_eventT0_scint_end
Corrected EventT0 for EKLM scintillator hits.
TH1F * h_z_bklm_rpc_phi
Normalised di-muon time difference, BKLM RPC phi readout.
TH1F * h_z_eklm_scint
Normalised di-muon time difference, EKLM scintillator.
TProfile * m_ProfileRpcPhi
For BKLM RPC phi plane.
TGraphErrors * gre_ctime_channel_scint
BKLM Scintillator.
TH2F * h2_deltaT0_vs_nhits_scint
DeltaT0 vs total hit count for BKLM scintillator.
TH2F * h2_deltaT0_vs_nhits_rpc
DeltaT0 vs total hit count for RPC.
TProfile * m_ProfileEKLMScintillatorPlane2
For EKLM scintillator plane2.
TH1F * h_time_rpc
BKLM RPC part.
TH1F * h_timeFSL_end[2][4][14]
EKLM part.
TProfile * m_Profile2RpcPhi
For BKLM RPC phi plane.
TH1F * h_timeF_scint[2]
BKLM scintillator part.
TH1F * hc_timeFSL[2][8][15]
BKLM part.
TProfile * m_Profile2BKLMScintillatorZ
For BKLM scintillator z plane.
TH2F * h2_deltaT0_vs_v_scint_end
DeltaT0 vs inverse hit count for EKLM scintillator.
TH1F * h_timeFS_rpc[2][8]
BKLM RPC part.
KLMChannelIndex m_klmChannels
KLM ChannelIndex object.
TGraph * gr_timeShift_channel_rpc
BKLM RPC.
std::map< KLMChannelNumber, double > m_timeRes
Resolution values of each channel.
TH1F * h_eventT0_scint_end
EventT0 seen by EKLM scintillator hits.
TH1F * h_diff
Distance between global and local position.
TH1F * hc_eventT0_rpc_highN
Corrected EventT0 for RPC with high hit count.
TH1F * h_nHits_plus_scint
Number of BKLM scintillator hits per mu+ track.
TH2F * h2_timeF_scint[2]
BKLM scintillator part.
TH1F * h_time_scint_tc
BKLM scintillator part.
double m_LowerTimeBoundaryRPC
Lower time boundary for RPC.
std::map< KLMChannelNumber, double > m_ctime_channel
Calibrated time distribution central value of each channel.
double m_LowerTimeBoundaryCalibratedRPC
Lower time boundary for RPC (calibrated data).
TH2F * h2_deltaT0_vs_nhits_scint_end
DeltaT0 vs total hit count for EKLM scintillator.
bool m_useEventT0
Whether to use event T0 from CDC.
TProfile * m_ProfileEKLMScintillatorPlane1
For EKLM scintillator plane1.
double m_UpperTimeBoundaryRPC
Upper time boundary for RPC.
TH2F * h2c_timeF_rpc[2]
BKLM RPC part.
TH1F * hc_eventT0_scint_highN
Corrected EventT0 for BKLM scintillator with high hit count.
TF1 * fcn_land
Landau function.
KLMTimeCableDelay * m_timeCableDelay
DBObject of the calibration constant of each channel due to cable decay.
TH2F * h2_deltaT0_vs_v_rpc
DeltaT0 vs inverse hit count for RPC.
TProfile * m_ProfileRpcZ
For BKLM RPC z plane.
bool m_useFixedRPCDelay
Whether to use fixed propagation delay for RPCs.
std::map< KLMChannelNumber, double > mc_etime_channel
Calibrated time distribution central value Error of each channel.
void readCalibrationDataBatch(std::function< bool(const KLMChannelIndex &)> channelFilter)
Load calibration data for a specific batch of channels.
TH1F * hc_eventT0_rpc_midN
Corrected EventT0 for RPC with medium hit count.
TH1F * h_nHits_minus_scint
Number of BKLM scintillator hits per mu- track.
TProfile * prof_deltaT0_rms_vs_v_scint
DeltaT0 RMS vs inverse hit count profile for BKLM scintillator.
std::map< KLMChannelNumber, int > m_cFlag
Calibration flag if the channel has enough hits collected and fitted OK.
TGraphErrors * gre_time_channel_rpc
BKLM RPC.
TH2F * h2_timeFS_end[2][4]
EKLM part.
TH2F * h2_timeF_rpc[2]
BKLM RPC part.
TH1F * h_nHits_plus_rpc
Number of RPC hits per mu+ track.
std::map< KLMChannelNumber, std::vector< struct Event > > m_evts
Container of hit information.
TDirectory * m_channelHistDir_EKLM[2][4][14][2]
Directory structure for per-channel histograms (EKLM).
TH1F * hc_eventT0_rpc
Corrected EventT0 for RPC hits.
TH1F * h_nHits_minus_scint_end
Number of EKLM scintillator hits per mu- track.
TH1F * h_timeFS_scint_end[2][4]
EKLM part.
double m_time_channelAvg_rpc
Central value of the global time distribution (BKLM RPC part).
TH2F * h2c_timeFS[2][8]
BKLM part.
double m_etime_channelAvg_rpc
Central value error of the global time distribution (BKLM RPC part).
double m_LowerTimeBoundaryScintillatorsEKLM
Lower time boundary for EKLM scintillators.
void readCalibrationDataCounts(std::map< KLMChannelNumber, unsigned int > &eventCounts)
Count events per channel (lightweight scan without loading full data).
std::map< KLMChannelNumber, double > m_etime_channel
Time distribution central value Error of each channel.
int m_lower_limit_counts
Lower limit of hits collected for on single channel.
TH1F * hc_eventT0_scint_end_lowN
Corrected EventT0 for EKLM scintillator with low hit count.
void readCalibrationDataFor2DFit(const std::vector< std::pair< KLMChannelNumber, unsigned int > > &channelsBKLM, const std::vector< std::pair< KLMChannelNumber, unsigned int > > &channelsEKLM)
Load calibration data only for channels needed for 2D fit.
@ c_BKLM
BKLM scintillator.
@ c_EKLM
EKLM scintillator.
uint16_t KLMChannelNumber
Channel number.
1.13.2