diff --git a/macro/beamtime/mcbm2024/bmon_Q_Factor.C b/macro/beamtime/mcbm2024/bmon_Q_Factor.C
new file mode 100644
index 0000000000000000000000000000000000000000..13dec765574951cdebc25ce463be31dc08b7ea03
--- /dev/null
+++ b/macro/beamtime/mcbm2024/bmon_Q_Factor.C
@@ -0,0 +1,309 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer] */
+
+
+double_t ExtractQFactor(TH1* pHistoIn)
+{
+ // Q-Factor = Max Bin Content / Mean Content of all bin in range
+ // => Tend toward 1 if bins are more identical
+ if (pHistoIn->Integral()) {
+ return (pHistoIn->GetBinContent(pHistoIn->GetMaximumBin())) / (pHistoIn->Integral() / pHistoIn->GetNbinsX());
+ }
+ else {
+ return 0.0;
+ }
+}
+double_t ExtractMean(TH1* pHistoIn)
+{
+ // Q-Factor = Max Bin Content / Mean Content of all bin in range
+ // => Tend toward 1 if bins are more identical
+ if (pHistoIn->Integral()) {
+ return (pHistoIn->Integral() / pHistoIn->GetNbinsX());
+ }
+ else {
+ return 0.0;
+ }
+}
+
+/** @brief Macro for check of Hades-like Q-Factor based on bmon digis, against variations of bin size and integration
+ * length
+ ** @param input Name of input file
+ ** @param nTimeSlices Number of time-slices to process
+ ** @param dTsSizeNs Size of one Ts in ns = Total length of one time disribution histogram
+ **/
+void bmon_Q_Factor(TString inputFileName, uint32_t uRunId, size_t numTimeslices = 0, double_t dTsSizeNs = 128000000.0)
+{
+
+ gROOT->cd();
+ TFile* file = new TFile(inputFileName, "READ");
+ TTree* tree = (TTree*) (file->Get("cbmsim"));
+
+ std::vector<CbmBmonDigi>* vDigisBmon = new std::vector<CbmBmonDigi>();
+ tree->SetBranchAddress("BmonDigi", &vDigisBmon);
+
+ uint32_t nentries = tree->GetEntries();
+ cout << "Entries: " << nentries << endl;
+ nentries = (numTimeslices && numTimeslices < nentries ? numTimeslices : nentries);
+
+ gROOT->cd();
+
+ uint32_t uMinNbBins = 10;
+ uint32_t uMaxNbBins = 100000;
+ /*
+ /// Hint: keep fractions of TS size and under 100 us
+ std::vector<double_t> vdBinSizesNs = {10, 20, 40, 80, 100, 200, 400, 800, 1.28e3, 2.56e3, 5.12e3, 12.8e3, 25.6e3,
+ 51.2e3}; // 14 values
+ /// Hint: keep fractions of TS size + multiples of bin size and above 10 us
+ std::vector<double_t> vdIntegrationNs = {12.8e3, 25.6e3, 51.2e3, 102.4e3, 204.8e3, 512e3, 1.28e6, 2.56e6, 5.12e6,
+ 12.8e6, 25.6e6, 32e6, 64e6}; // 13 values
+ */
+ /// Hint: keep fractions of TS size and under 100 us
+ std::vector<double_t> vdBinSizesNs = {10, 100, 1.28e3, 25.6e3}; // 4 values, biggest ~ HADES style
+ /// Hint: keep fractions of TS size + multiples of bin size and above 10 us
+ std::vector<double_t> vdIntegrationNs = {102.4e3, 512e3, 1.28e6, 32e6, 64e6}; // 4 values, prev to last ~ HADES style
+ /// Dimension: same as BinSizes vector!!
+ std::vector<double_t> vdQfactorHistMax = {2000., 400., 40., 20.};
+
+ std::vector<std::vector<uint32_t>> vuNbBinsHisto(vdIntegrationNs.size(),
+ std::vector<uint32_t>(vdBinSizesNs.size(), 0));
+ std::vector<uint32_t> vuNbHistoCyclesPerTS(vdIntegrationNs.size(), 0);
+ std::vector<std::vector<TH1*>> vHistoOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vHistoNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vQvalOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vQvalNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vMeanOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vMeanNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+
+ uint16_t uNbPlots = 0;
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ /// Pre-check values before in spreadsheet to make sure integer !!!!
+ vuNbHistoCyclesPerTS[uHistSz] = dTsSizeNs / vdIntegrationNs[uHistSz];
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ /// Pre-check values before in spreadsheet to make sure integer !!!!
+ vuNbBinsHisto[uHistSz][uBinSz] = vdIntegrationNs[uHistSz] / vdBinSizesNs[uBinSz];
+ if (uMinNbBins <= vuNbBinsHisto[uHistSz][uBinSz] /*&& vuNbBinsHisto[uHistSz][uBinSz] <= uMaxNbBins*/) {
+ vHistoOld[uHistSz][uBinSz] =
+ new TH1D(Form("binHist_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin in cycle of range %9.0f ns, Old Bmon; Time in Cycle [ns]; Digis []",
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbBinsHisto[uHistSz][uBinSz], 0.0, vdIntegrationNs[uHistSz]);
+ vHistoNew[uHistSz][uBinSz] =
+ new TH1D(Form("binHist_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin in cycle of range %9.0f ns, New Bmon; Time in Cycle [ns]; Digis []",
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbBinsHisto[uHistSz][uBinSz], 0.0, vdIntegrationNs[uHistSz]);
+
+ double_t dBinOffset = 1.0 / (2.0 * vuNbHistoCyclesPerTS[uHistSz]);
+ vQvalOld[uHistSz][uBinSz] =
+ new TH1D(Form("evoQFactor_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, Old Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+ vQvalNew[uHistSz][uBinSz] =
+ new TH1D(Form("evoQFactor_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, New Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+
+ vMeanOld[uHistSz][uBinSz] =
+ new TH1D(Form("evoMean_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, Old Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+ vMeanNew[uHistSz][uBinSz] =
+ new TH1D(Form("evoMean_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, New Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+
+ uNbPlots++;
+ }
+ }
+ }
+
+ std::vector<uint32_t> vuIdxHistoCycleinTS(vdIntegrationNs.size(), 0);
+
+ /// Loop on timeslices
+ for (Int_t iEntry = 0; iEntry < nentries; iEntry++) {
+ tree->GetEntry(iEntry);
+ uint32_t nDigisBmon = vDigisBmon->size();
+
+ if (nDigisBmon < 300) {
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => spill break, skipping this TS"
+ << std::endl;
+ continue;
+ }
+
+ std::vector<CbmBmonDigi> vDigisOld;
+ std::vector<CbmBmonDigi> vDigisNew;
+ vDigisOld.reserve(nDigisBmon);
+ vDigisNew.reserve(nDigisBmon);
+ for (auto& digi : *vDigisBmon) {
+ if (1 == CbmTofAddress::GetChannelSide(digi.GetAddress())) {
+ if (CbmTofAddress::GetChannelId(digi.GetAddress()) < 4) {
+ vDigisNew.push_back(digi);
+ }
+ else {
+ LOG(fatal) << "Bad sCVD channel: " << CbmTofAddress::GetChannelId(digi.GetAddress());
+ }
+ }
+ else {
+ vDigisOld.push_back(digi);
+ }
+ }
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => " << vDigisOld.size() << " Old + "
+ << vDigisNew.size() << " sCVD" << std::endl;
+
+ vuIdxHistoCycleinTS.assign(vdIntegrationNs.size(), 0);
+ for (auto itOld = vDigisOld.begin(); itOld != vDigisOld.end(); ++itOld) {
+ double_t dTime = (*itOld).GetTime();
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ uint32_t uCurrentCycle = std::floor(dTime / vdIntegrationNs[uHistSz]);
+ if (vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle) {
+ for (; vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle; ++vuIdxHistoCycleinTS[uHistSz]) {
+ double_t dTsFractional = (vdIntegrationNs[uHistSz] * vuIdxHistoCycleinTS[uHistSz]) / dTsSizeNs + iEntry;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ double_t dQFactor = ExtractQFactor(vHistoOld[uHistSz][uBinSz]);
+ vQvalOld[uHistSz][uBinSz]->Fill(dTsFractional, dQFactor);
+ vMeanOld[uHistSz][uBinSz]->Fill(dTsFractional, ExtractMean(vHistoOld[uHistSz][uBinSz]));
+ /*
+ if (1280 == vdBinSizesNs[uBinSz]) {
+ std::cout << Form("%9f %9.0f %9.0f %9d %9.0f",
+ dTsFractional, vdIntegrationNs[uHistSz], vHistoOld[uHistSz][uBinSz]->Integral(),
+ vHistoOld[uHistSz][uBinSz]->GetNbinsX(),
+ vHistoOld[uHistSz][uBinSz]->GetBinContent(vHistoOld[uHistSz][uBinSz]->GetMaximumBin()))
+ << std::endl;
+ }
+ */
+ if (0.0 < dQFactor) {
+ vHistoOld[uHistSz][uBinSz]->Reset();
+ }
+ }
+ }
+ }
+ }
+
+ double_t dTimeInCycle = std::fmod(dTime, vdIntegrationNs[uHistSz]);
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ vHistoOld[uHistSz][uBinSz]->Fill(dTimeInCycle);
+ }
+ }
+ }
+ }
+
+ vuIdxHistoCycleinTS.assign(vdIntegrationNs.size(), 0);
+ uint32_t uDigiIdx = 0;
+ for (auto itNew = vDigisNew.begin(); itNew != vDigisNew.end(); ++itNew) {
+ double_t dTime = (*itNew).GetTime();
+ if (dTime < 0) {
+ std::cout << Form("TS %5d Digi %6u %8.0f!!!!!", iEntry, uDigiIdx, dTime) << std::endl;
+ continue;
+ }
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ uint32_t uCurrentCycle = std::floor(dTime / vdIntegrationNs[uHistSz]);
+ if (vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle) {
+ for (; vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle; ++vuIdxHistoCycleinTS[uHistSz]) {
+ double_t dTsFractional = (vdIntegrationNs[uHistSz] * vuIdxHistoCycleinTS[uHistSz]) / dTsSizeNs + iEntry;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalNew[uHistSz][uBinSz]) {
+ double_t dQFactor = ExtractQFactor(vHistoNew[uHistSz][uBinSz]);
+ vQvalNew[uHistSz][uBinSz]->Fill(dTsFractional, dQFactor);
+ vMeanNew[uHistSz][uBinSz]->Fill(dTsFractional, ExtractMean(vHistoNew[uHistSz][uBinSz]));
+ if (0.0 < dQFactor) {
+ vHistoNew[uHistSz][uBinSz]->Reset();
+ }
+ }
+ }
+ }
+ }
+
+ double_t dTimeInCycle = std::fmod(dTime, vdIntegrationNs[uHistSz]);
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalNew[uHistSz][uBinSz]) {
+ vHistoNew[uHistSz][uBinSz]->Fill(dTimeInCycle);
+ }
+ }
+ }
+ uDigiIdx++;
+ }
+ }
+
+ uint16_t nPadX = std::ceil(std::sqrt(uNbPlots));
+ uint16_t nPadY = std::ceil(1.0 * uNbPlots / nPadX);
+ std::cout << Form("Nb plots %3d Nb pads X %2d Y %2d", uNbPlots, nPadX, nPadY) << std::endl;
+
+ TCanvas* test = new TCanvas("test", "test");
+ test->Divide(nPadX, nPadY);
+
+ TCanvas* testMean = new TCanvas("testMean", "testMean");
+ testMean->Divide(nPadX, nPadY);
+
+ uint32_t uPadIdx = 1;
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ THStack* pStack = new THStack(Form("pStack_%2d_%2d", uHistSz, uBinSz),
+ Form("Q-Factor, Old and New Bmon, run %4d, %5.0f ns bin, %9.0f ns range", uRunId,
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]));
+ vQvalOld[uHistSz][uBinSz]->SetLineColor(kBlue);
+ vQvalOld[uHistSz][uBinSz]->SetLineWidth(2);
+ vQvalOld[uHistSz][uBinSz]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[uBinSz]);
+ pStack->Add(vQvalOld[uHistSz][uBinSz]);
+
+ vQvalNew[uHistSz][uBinSz]->SetLineColor(kRed);
+ vQvalNew[uHistSz][uBinSz]->SetLineWidth(2);
+ vQvalNew[uHistSz][uBinSz]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[uBinSz]);
+ pStack->Add(vQvalNew[uHistSz][uBinSz]);
+
+ test->cd(uPadIdx);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStack->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ THStack* pStackMean =
+ new THStack(Form("pStackMean_%2d_%2d", uHistSz, uBinSz),
+ Form("Mean bin content, Old and New Bmon, run %4d, %5.0f ns bin, %9.0f ns range", uRunId,
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]));
+ vMeanOld[uHistSz][uBinSz]->SetLineColor(kBlue);
+ vMeanOld[uHistSz][uBinSz]->SetLineWidth(2);
+ pStackMean->Add(vMeanOld[uHistSz][uBinSz]);
+
+ vMeanNew[uHistSz][uBinSz]->SetLineColor(kRed);
+ vMeanNew[uHistSz][uBinSz]->SetLineWidth(2);
+ pStackMean->Add(vMeanNew[uHistSz][uBinSz]);
+
+ testMean->cd(uPadIdx);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackMean->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ uPadIdx++;
+ }
+ }
+ }
+
+ TFile* outFile = new TFile(Form("data/bmon_q_factor_%04u_%03luTs.root", uRunId, numTimeslices), "RECREATE");
+ outFile->cd();
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ vQvalOld[uHistSz][uBinSz]->Write();
+ vQvalNew[uHistSz][uBinSz]->Write();
+ vMeanOld[uHistSz][uBinSz]->Write();
+ vMeanNew[uHistSz][uBinSz]->Write();
+ }
+ }
+ }
+ test->Write();
+ testMean->Write();
+
+ gROOT->cd();
+
+ outFile->Close();
+}
diff --git a/macro/beamtime/mcbm2024/bmon_Q_Factor_plots_2391_2984_3006.C b/macro/beamtime/mcbm2024/bmon_Q_Factor_plots_2391_2984_3006.C
new file mode 100644
index 0000000000000000000000000000000000000000..4c31fd6045b04fc1a50b539ce70f2f91e8cacbea
--- /dev/null
+++ b/macro/beamtime/mcbm2024/bmon_Q_Factor_plots_2391_2984_3006.C
@@ -0,0 +1,301 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer] */
+
+
+void bmon_Q_Factor_plots_2391_2984_3006()
+{
+ std::vector<uint32_t> vRunId = {2391, 2894, 3006};
+ // std::vector<uint32_t> vRunId = {3006, 3107, 3109};
+
+ /// Hint: keep fractions of TS size and under 100 us
+ std::vector<double_t> vdBinSizesNs = {10, 100, 1.28e3, 25.6e3}; // 4 values, biggest ~ HADES style
+ /// Hint: keep fractions of TS size + multiples of bin size and above 10 us
+ std::vector<double_t> vdIntegrationNs = {102.4e3, 512e3, 1.28e6, 32e6, 64e6}; // 4 values, prev to last ~ HADES style
+ /// Dimension: same as BinSizes vector!!
+ std::vector<double_t> vdQfactorHistMax = {2000., 400., 40., 20.};
+
+ std::map<uint32_t, std::vector<std::vector<TH1*>>> vQvalOld;
+ std::map<uint32_t, std::vector<std::vector<TH1*>>> vQvalNew;
+ std::map<uint32_t, std::vector<std::vector<TH1*>>> vMeanOld;
+ std::map<uint32_t, std::vector<std::vector<TH1*>>> vMeanNew;
+
+ for (uint32_t uRun = 0; uRun < vRunId.size(); ++uRun) {
+ TFile* inFile = new TFile(Form("data/bmon_q_factor_%04u_080Ts.root", vRunId[uRun]), "READ");
+ gROOT->cd();
+ std::cout << vRunId[uRun] << std::endl;
+
+ vQvalOld[uRun].resize(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ vQvalNew[uRun].resize(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ vMeanOld[uRun].resize(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ vMeanNew[uRun].resize(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ TH1* pTemp = dynamic_cast<TH1*>(
+ inFile->FindObjectAny(Form("evoQFactor_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(
+ Form("EvoQFactor_Old_%04d_%9.0f_%5.0f", vRunId[uRun], vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ }
+ vQvalOld[uRun][uHistSz][uBinSz] = pTemp;
+
+ pTemp = dynamic_cast<TH1*>(
+ inFile->FindObjectAny(Form("evoQFactor_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(
+ Form("EvoQFactor_New_%04d_%9.0f_%5.0f", vRunId[uRun], vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ }
+ vQvalNew[uRun][uHistSz][uBinSz] = pTemp;
+
+ pTemp = dynamic_cast<TH1*>(
+ inFile->FindObjectAny(Form("evoMean_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(
+ Form("EvoMean_Old_%04d_%9.0f_%5.0f", vRunId[uRun], vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ }
+ vMeanOld[uRun][uHistSz][uBinSz] = pTemp;
+
+ pTemp = dynamic_cast<TH1*>(
+ inFile->FindObjectAny(Form("evoMean_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(
+ Form("EvoMean_New_%04d_%9.0f_%5.0f", vRunId[uRun], vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz])));
+ }
+ vMeanNew[uRun][uHistSz][uBinSz] = pTemp;
+ }
+ }
+ inFile->Close();
+ }
+
+ /// HADES ---------------------------------------------------------------------------------------------------------///
+ uint32_t uHadesHistSz = 3;
+ uint32_t uHadesBinSz = 3;
+ TCanvas* cHadesQfactor = new TCanvas("cHadesQfactor", Form("HADES-like Q factors: %5.0f ns bin, %9.0f ns range",
+ vdBinSizesNs[uHadesBinSz], vdIntegrationNs[uHadesHistSz]));
+ cHadesQfactor->Divide(vRunId.size(), 2);
+
+ std::vector<THStack*> vStacksHadesQfactor(2 * vRunId.size(), nullptr);
+ for (uint32_t uRun = 0; uRun < vRunId.size(); ++uRun) {
+
+ vStacksHadesQfactor[uRun] =
+ new THStack(Form("stackHades_Qfact_%04d", vRunId[uRun]), Form("Q-Factor, run %4d", vRunId[uRun]));
+ if (nullptr != vQvalOld[uRun][uHadesHistSz][uHadesBinSz]) {
+ vQvalOld[uRun][uHadesHistSz][uHadesBinSz]->SetLineColor(kBlue);
+ vQvalOld[uRun][uHadesHistSz][uHadesBinSz]->SetLineWidth(2);
+ vQvalOld[uRun][uHadesHistSz][uHadesBinSz]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[uHadesBinSz]);
+ vStacksHadesQfactor[uRun]->Add(vQvalOld[uRun][uHadesHistSz][uHadesBinSz]);
+ }
+
+ if (nullptr != vQvalNew[uRun][uHadesHistSz][uHadesBinSz]) {
+ vQvalNew[uRun][uHadesHistSz][uHadesBinSz]->SetLineColor(kRed);
+ vQvalNew[uRun][uHadesHistSz][uHadesBinSz]->SetLineWidth(2);
+ vQvalNew[uRun][uHadesHistSz][uHadesBinSz]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[uHadesBinSz]);
+ vStacksHadesQfactor[uRun]->Add(vQvalNew[uRun][uHadesHistSz][uHadesBinSz]);
+ }
+
+ cHadesQfactor->cd(1 + uRun);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksHadesQfactor[uRun]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ vStacksHadesQfactor[vRunId.size() + uRun] =
+ new THStack(Form("stackHades_Mean_%04d", vRunId[uRun]), Form("Mean bin content, run %4d", vRunId[uRun]));
+ if (nullptr != vMeanOld[uRun][uHadesHistSz][uHadesBinSz]) {
+ vMeanOld[uRun][uHadesHistSz][uHadesBinSz]->SetLineColor(kBlue);
+ vMeanOld[uRun][uHadesHistSz][uHadesBinSz]->SetLineWidth(2);
+ vStacksHadesQfactor[vRunId.size() + uRun]->Add(vMeanOld[uRun][uHadesHistSz][uHadesBinSz]);
+ }
+
+ if (nullptr != vMeanNew[uRun][uHadesHistSz][uHadesBinSz]) {
+ vMeanNew[uRun][uHadesHistSz][uHadesBinSz]->SetLineColor(kRed);
+ vMeanNew[uRun][uHadesHistSz][uHadesBinSz]->SetLineWidth(2);
+ vStacksHadesQfactor[vRunId.size() + uRun]->Add(vMeanNew[uRun][uHadesHistSz][uHadesBinSz]);
+ }
+
+ cHadesQfactor->cd(1 + vRunId.size() + uRun);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksHadesQfactor[vRunId.size() + uRun]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+ }
+ ///----------------------------------------------------------------------------------------------------------------///
+
+ /// Timescale -----------------------------------------------------------------------------------------------------///
+ uint32_t uTimescaleHistSz = 4;
+ TCanvas* cTimescaleQfactor =
+ new TCanvas("cTimescaleQfactor", Form("Q factors vs timescale: %9.0f ns range", vdIntegrationNs[uTimescaleHistSz]));
+ cTimescaleQfactor->Divide(vdBinSizesNs.size(), vRunId.size());
+ TCanvas* cTimescaleMean = new TCanvas(
+ "cTimescaleMean", Form("Mean bin content vs timescale: %9.0f ns range", vdIntegrationNs[uTimescaleHistSz]));
+ cTimescaleMean->Divide(vdBinSizesNs.size(), vRunId.size());
+
+ std::vector<std::vector<THStack*>> vStacksTimescaleQfactor(vRunId.size(),
+ std::vector<THStack*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<THStack*>> vStacksTimescaleMean(vRunId.size(),
+ std::vector<THStack*>(vdBinSizesNs.size(), nullptr));
+ for (uint32_t uRun = 0; uRun < vRunId.size(); ++uRun) {
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ vStacksTimescaleQfactor[uRun][uBinSz] =
+ new THStack(Form("stackBinSz_Qfact_%04d_%02d", vRunId[uRun], uBinSz),
+ Form("Q-Factor, run %4d, %5.0f ns bin", vRunId[uRun], vdBinSizesNs[uBinSz]));
+ if (nullptr != vQvalOld[uRun][uTimescaleHistSz][uBinSz]) {
+ vQvalOld[uRun][uTimescaleHistSz][uBinSz]->SetLineColor(kBlue);
+ vQvalOld[uRun][uTimescaleHistSz][uBinSz]->SetLineWidth(2);
+ vQvalOld[uRun][uTimescaleHistSz][uBinSz]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[uBinSz]);
+ vStacksTimescaleQfactor[uRun][uBinSz]->Add(vQvalOld[uRun][uTimescaleHistSz][uBinSz]);
+ }
+
+ if (nullptr != vQvalNew[uRun][uTimescaleHistSz][uBinSz]) {
+ vQvalNew[uRun][uTimescaleHistSz][uBinSz]->SetLineColor(kRed);
+ vQvalNew[uRun][uTimescaleHistSz][uBinSz]->SetLineWidth(2);
+ vQvalNew[uRun][uTimescaleHistSz][uBinSz]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[uBinSz]);
+ vStacksTimescaleQfactor[uRun][uBinSz]->Add(vQvalNew[uRun][uTimescaleHistSz][uBinSz]);
+ }
+
+ cTimescaleQfactor->cd(1 + uRun * vdBinSizesNs.size() + uBinSz);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksTimescaleQfactor[uRun][uBinSz]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ vStacksTimescaleMean[uRun][uBinSz] =
+ new THStack(Form("stackBinSz_Mean_%04d_%02d", vRunId[uRun], uBinSz),
+ Form("Mean bin content, run %4d, %5.0f ns bin", vRunId[uRun], vdBinSizesNs[uBinSz]));
+ if (nullptr != vMeanOld[uRun][uTimescaleHistSz][uBinSz]) {
+ vMeanOld[uRun][uTimescaleHistSz][uBinSz]->SetLineColor(kBlue);
+ vMeanOld[uRun][uTimescaleHistSz][uBinSz]->SetLineWidth(2);
+ vStacksTimescaleMean[uRun][uBinSz]->Add(vMeanOld[uRun][uTimescaleHistSz][uBinSz]);
+ }
+
+ if (nullptr != vMeanNew[uRun][uTimescaleHistSz][uBinSz]) {
+ vMeanNew[uRun][uTimescaleHistSz][uBinSz]->SetLineColor(kRed);
+ vMeanNew[uRun][uTimescaleHistSz][uBinSz]->SetLineWidth(2);
+ vStacksTimescaleMean[uRun][uBinSz]->Add(vMeanNew[uRun][uTimescaleHistSz][uBinSz]);
+ }
+
+ cTimescaleMean->cd(1 + uRun * vdBinSizesNs.size() + uBinSz);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksTimescaleMean[uRun][uBinSz]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+ }
+ }
+ ///----------------------------------------------------------------------------------------------------------------///
+
+ /// Integration ---------------------------------------------------------------------------------------------------///
+ uint32_t uIntegrationHistSzSmall = 1;
+ uint32_t uIntegrationHistSzBig = 3;
+ TCanvas* cIntegrationQfactor =
+ new TCanvas("cIntegrationQfactor", "Q factors vs Integration length, min/max timescales");
+ cIntegrationQfactor->Divide(4, vRunId.size());
+
+ std::vector<std::vector<THStack*>> vStacksIntegrationQfactor(vRunId.size(), std::vector<THStack*>(4, nullptr));
+ for (uint32_t uRun = 0; uRun < vRunId.size(); ++uRun) {
+ /// Small scale, small integration
+ vStacksIntegrationQfactor[uRun][0] =
+ new THStack(Form("stackIntegr_Qfact_%04d_%02d", vRunId[uRun], 0),
+ Form("Q-Factor, run %4d, %5.0f ns bin, %9.0f ns range", vRunId[uRun], vdBinSizesNs[0],
+ vdIntegrationNs[uIntegrationHistSzSmall]));
+ if (nullptr != vQvalOld[uRun][uIntegrationHistSzSmall][0]) {
+ vQvalOld[uRun][uIntegrationHistSzSmall][0]->SetLineColor(kBlue);
+ vQvalOld[uRun][uIntegrationHistSzSmall][0]->SetLineWidth(2);
+ vQvalOld[uRun][uIntegrationHistSzSmall][0]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[0]);
+ vStacksIntegrationQfactor[uRun][0]->Add(vQvalOld[uRun][uIntegrationHistSzSmall][0]);
+ }
+
+ if (nullptr != vQvalNew[uRun][uIntegrationHistSzSmall][0]) {
+ vQvalNew[uRun][uIntegrationHistSzSmall][0]->SetLineColor(kRed);
+ vQvalNew[uRun][uIntegrationHistSzSmall][0]->SetLineWidth(2);
+ vQvalNew[uRun][uIntegrationHistSzSmall][0]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[0]);
+ vStacksIntegrationQfactor[uRun][0]->Add(vQvalNew[uRun][uIntegrationHistSzSmall][0]);
+ }
+
+ cIntegrationQfactor->cd(1 + uRun * 4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksIntegrationQfactor[uRun][0]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ /// Large scale, small integration
+ vStacksIntegrationQfactor[uRun][1] =
+ new THStack(Form("stackIntegr_Qfact_%04d_%02d", vRunId[uRun], 1),
+ Form("Q-Factor, run %4d, %5.0f ns bin, %9.0f ns range", vRunId[uRun],
+ vdBinSizesNs[vdBinSizesNs.size() - 1], vdIntegrationNs[uIntegrationHistSzSmall]));
+ if (nullptr != vQvalOld[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]) {
+ vQvalOld[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]->SetLineColor(kBlue);
+ vQvalOld[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]->SetLineWidth(2);
+ vQvalOld[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]->GetYaxis()->SetRangeUser(
+ 0., vdQfactorHistMax[vdBinSizesNs.size() - 1]);
+ vStacksIntegrationQfactor[uRun][1]->Add(vQvalOld[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]);
+ }
+
+ if (nullptr != vQvalNew[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]) {
+ vQvalNew[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]->SetLineColor(kRed);
+ vQvalNew[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]->SetLineWidth(2);
+ vQvalNew[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]->GetYaxis()->SetRangeUser(
+ 0., vdQfactorHistMax[vdBinSizesNs.size() - 1]);
+ vStacksIntegrationQfactor[uRun][1]->Add(vQvalNew[uRun][uIntegrationHistSzSmall][vdBinSizesNs.size() - 1]);
+ }
+
+ cIntegrationQfactor->cd(2 + uRun * 4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksIntegrationQfactor[uRun][1]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ /// Small scale, large integration
+ vStacksIntegrationQfactor[uRun][2] =
+ new THStack(Form("stackIntegr_Qfact_%04d_%02d", vRunId[uRun], 2),
+ Form("Q-Factor, run %4d, %5.0f ns bin, %9.0f ns range", vRunId[uRun], vdBinSizesNs[0],
+ vdIntegrationNs[uIntegrationHistSzBig]));
+ if (nullptr != vQvalOld[uRun][uIntegrationHistSzBig][0]) {
+ vQvalOld[uRun][uIntegrationHistSzBig][0]->SetLineColor(kBlue);
+ vQvalOld[uRun][uIntegrationHistSzBig][0]->SetLineWidth(2);
+ vQvalOld[uRun][uIntegrationHistSzBig][0]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[0]);
+ vStacksIntegrationQfactor[uRun][2]->Add(vQvalOld[uRun][uIntegrationHistSzBig][0]);
+ }
+
+ if (nullptr != vQvalNew[uRun][uIntegrationHistSzBig][0]) {
+ vQvalNew[uRun][uIntegrationHistSzBig][0]->SetLineColor(kRed);
+ vQvalNew[uRun][uIntegrationHistSzBig][0]->SetLineWidth(2);
+ vQvalNew[uRun][uIntegrationHistSzBig][0]->GetYaxis()->SetRangeUser(0., vdQfactorHistMax[0]);
+ vStacksIntegrationQfactor[uRun][2]->Add(vQvalNew[uRun][uIntegrationHistSzBig][0]);
+ }
+
+ cIntegrationQfactor->cd(3 + uRun * 4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksIntegrationQfactor[uRun][2]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ /// Large scale, large integration
+ vStacksIntegrationQfactor[uRun][3] =
+ new THStack(Form("stackIntegr_Qfact_%04d_%02d", vRunId[uRun], 3),
+ Form("Q-Factor, run %4d, %5.0f ns bin, %9.0f ns range", vRunId[uRun],
+ vdBinSizesNs[vdBinSizesNs.size() - 1], vdIntegrationNs[uIntegrationHistSzBig]));
+ if (nullptr != vQvalOld[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]) {
+ vQvalOld[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]->SetLineColor(kBlue);
+ vQvalOld[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]->SetLineWidth(2);
+ vQvalOld[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]->GetYaxis()->SetRangeUser(
+ 0., vdQfactorHistMax[vdBinSizesNs.size() - 1]);
+ vStacksIntegrationQfactor[uRun][3]->Add(vQvalOld[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]);
+ }
+
+ if (nullptr != vQvalNew[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]) {
+ vQvalNew[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]->SetLineColor(kRed);
+ vQvalNew[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]->SetLineWidth(2);
+ vQvalNew[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]->GetYaxis()->SetRangeUser(
+ 0., vdQfactorHistMax[vdBinSizesNs.size() - 1]);
+ vStacksIntegrationQfactor[uRun][3]->Add(vQvalNew[uRun][uIntegrationHistSzBig][vdBinSizesNs.size() - 1]);
+ }
+
+ cIntegrationQfactor->cd(4 + uRun * 4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vStacksIntegrationQfactor[uRun][3]->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+ }
+ ///----------------------------------------------------------------------------------------------------------------///
+}
diff --git a/macro/beamtime/mcbm2024/bmon_Q_Factor_timescale.C b/macro/beamtime/mcbm2024/bmon_Q_Factor_timescale.C
new file mode 100644
index 0000000000000000000000000000000000000000..215639addea0410ed6af2547c98e413f15455430
--- /dev/null
+++ b/macro/beamtime/mcbm2024/bmon_Q_Factor_timescale.C
@@ -0,0 +1,471 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer] */
+
+
+std::vector<double> GenerateLogBinArray(uint32_t uNbDecadesLog, uint32_t uNbStepsDecade, uint32_t uNbSubStepsInStep,
+ uint32_t& uNbBinsLog, int32_t iStartExp = 0, bool bAddZeroStart = false)
+{
+ /// Logarithmic bining for self time comparison
+ /// Number of log bins =
+ /// 9 for the sub-unit decade
+ /// + 9 for each unit of each decade * 10 for the subdecade range
+ /// + 1 for the closing bin top edge
+ uNbBinsLog = uNbStepsDecade + uNbStepsDecade * uNbSubStepsInStep * uNbDecadesLog;
+
+ /// Need uNbBinsLog + 1 values as we need to provide the end of last bin
+ uint32_t uArrayLength = uNbBinsLog + 1;
+ double dBinsLog[uArrayLength];
+ /// First fill sub-unit decade
+ for (uint32_t uSubU = 0; uSubU < uNbStepsDecade; uSubU++) {
+ dBinsLog[uSubU] = std::pow(10, iStartExp - 1) * (1 + uSubU);
+ }
+
+ /// Then fill the main decades
+ double dSubstepSize = 1.0 / uNbSubStepsInStep;
+ for (uint32_t uDecade = 0; uDecade < uNbDecadesLog; uDecade++) {
+ double dBase = std::pow(10, iStartExp + static_cast<int32_t>(uDecade));
+ uint32_t uDecadeIdx = uNbStepsDecade + uDecade * uNbStepsDecade * uNbSubStepsInStep;
+ for (uint32_t uStep = 0; uStep < uNbStepsDecade; uStep++) {
+ uint32_t uStepIdx = uDecadeIdx + uStep * uNbSubStepsInStep;
+ for (uint32_t uSubStep = 0; uSubStep < uNbSubStepsInStep; uSubStep++) {
+ dBinsLog[uStepIdx + uSubStep] = dBase * (1 + uStep) + dBase * dSubstepSize * uSubStep;
+ } // for( uint32_t uSubStep = 0; uSubStep < uNbSubStepsInStep; uSubStep++ )
+ } // for( uint32_t uStep = 0; uStep < uNbStepsDecade; uStep++ )
+ } // for( uint32_t uDecade = 0; uDecade < uNbDecadesLog; uDecade ++)
+ dBinsLog[uNbBinsLog] = std::pow(10, iStartExp + uNbDecadesLog);
+
+ /// use vector instead
+ std::vector<double> dBinsLogVect;
+
+ /// + 1 optional if bin [ 0; Min [ should be added
+ if (bAddZeroStart) {
+ uNbBinsLog++;
+ dBinsLogVect.push_back(0);
+ }
+
+ for (uint32_t i = 0; i < uArrayLength; ++i) {
+ dBinsLogVect.push_back(dBinsLog[i]);
+ } // for( uint32_t i = 0; i < uArrayLength; ++i )
+
+ return dBinsLogVect;
+}
+
+double_t ExtractQFactor(TH1* pHistoIn)
+{
+ // Q-Factor = Max Bin Content / Mean Content of all bin in range
+ // => Tend toward 1 if bins are more identical
+ if (pHistoIn->Integral()) {
+ return (pHistoIn->GetBinContent(pHistoIn->GetMaximumBin())) / (pHistoIn->Integral() / pHistoIn->GetNbinsX());
+ }
+ else {
+ return 0.0;
+ }
+}
+double_t ExtractMean(TH1* pHistoIn)
+{
+ // Q-Factor = Max Bin Content / Mean Content of all bin in range
+ // => Tend toward 1 if bins are more identical
+ if (pHistoIn->Integral()) {
+ return (pHistoIn->Integral() / pHistoIn->GetNbinsX());
+ }
+ else {
+ return 0.0;
+ }
+}
+
+/** @brief Macro for check of Hades-like Q-Factor based on bmon digis, against variations of bin size and integration
+ * length
+ ** @param input Name of input file
+ ** @param nTimeSlices Number of time-slices to process
+ ** @param dTsSizeNs Size of one Ts in ns = Total length of one time disribution histogram
+ **/
+void bmon_Q_Factor_timescale(TString inputFileName, uint32_t uRunId, size_t numTimeslices = 0, double_t dStartTs = 0,
+ double_t dStopTs = 0, double_t dTsSizeNs = 128000000.0)
+{
+
+ gROOT->cd();
+ TFile* file = new TFile(inputFileName, "READ");
+ TTree* tree = (TTree*) (file->Get("cbmsim"));
+
+ std::vector<CbmBmonDigi>* vDigisBmon = new std::vector<CbmBmonDigi>();
+ tree->SetBranchAddress("BmonDigi", &vDigisBmon);
+
+ uint32_t nentries = tree->GetEntries();
+ cout << "Entries: " << nentries << endl;
+ nentries = (numTimeslices && numTimeslices < nentries ? numTimeslices : nentries);
+
+ if (0 < dStopTs && dStopTs < nentries) {
+ nentries = dStopTs + 1;
+ }
+ numTimeslices = nentries - dStartTs;
+
+ gROOT->cd();
+
+ uint32_t uMinNbBins = 10;
+ uint32_t uMaxNbBins = 300000;
+ /// Hint: keep fractions of TS size and under 100 us
+ std::vector<double_t> vdBinSizesNs = {10, 20, 40, 80, 100, 200, 400, 800, 1.28e3, 2.56e3,
+ 5.12e3, 6.4e3, 10.24e3, 25.6e3, 51.2e3, 64e3, 102.4e3, 204.8e3}; // 14 values
+ /// Hint: keep fractions of TS size + multiples of bin size and above 10 us
+ std::vector<double_t> vdIntegrationNs = {2.56e6}; // 13 values
+ /// Dimension: same as BinSizes vector!!
+ std::vector<double_t> vdQfactorHistMax = {2000., 400., 40., 20.};
+
+ std::vector<std::vector<uint32_t>> vuNbBinsHisto(vdIntegrationNs.size(),
+ std::vector<uint32_t>(vdBinSizesNs.size(), 0));
+ std::vector<uint32_t> vuNbHistoCyclesPerTS(vdIntegrationNs.size(), 0);
+ std::vector<std::vector<TH1*>> vHistoOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vHistoNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vQvalOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vQvalNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vMeanOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vMeanNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+
+ std::vector<TH1*> vBinCountDistributionOld(vdBinSizesNs.size(), nullptr);
+ std::vector<TH1*> vBinCountDistributionNew(vdBinSizesNs.size(), nullptr);
+
+ uint16_t uNbPlots = 0;
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ /// Pre-check values before in spreadsheet to make sure integer !!!!
+ vuNbHistoCyclesPerTS[uHistSz] = dTsSizeNs / vdIntegrationNs[uHistSz];
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ /// Pre-check values before in spreadsheet to make sure integer !!!!
+ vuNbBinsHisto[uHistSz][uBinSz] = vdIntegrationNs[uHistSz] / vdBinSizesNs[uBinSz];
+ if (uMinNbBins <= vuNbBinsHisto[uHistSz][uBinSz] /*&& vuNbBinsHisto[uHistSz][uBinSz] <= uMaxNbBins*/) {
+ vHistoOld[uHistSz][uBinSz] =
+ new TH1D(Form("binHist_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin in cycle of range %9.0f ns, Old Bmon; Time in Cycle [ns]; Digis []",
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbBinsHisto[uHistSz][uBinSz], 0.0, vdIntegrationNs[uHistSz]);
+ vHistoNew[uHistSz][uBinSz] =
+ new TH1D(Form("binHist_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin in cycle of range %9.0f ns, New Bmon; Time in Cycle [ns]; Digis []",
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbBinsHisto[uHistSz][uBinSz], 0.0, vdIntegrationNs[uHistSz]);
+
+ double_t dBinOffset = 1.0 / (2.0 * vuNbHistoCyclesPerTS[uHistSz]);
+ vQvalOld[uHistSz][uBinSz] =
+ new TH1D(Form("evoQFactor_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, Old Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+ vQvalNew[uHistSz][uBinSz] =
+ new TH1D(Form("evoQFactor_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, New Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+
+ vMeanOld[uHistSz][uBinSz] =
+ new TH1D(Form("evoMean_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, Old Bmon; Time in Run [TS]; Mean []", uRunId,
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+ vMeanNew[uHistSz][uBinSz] =
+ new TH1D(Form("evoMean_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, New Bmon; Time in Run [TS]; Mean []", uRunId,
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+
+ uNbPlots++;
+ }
+ }
+ }
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ vBinCountDistributionOld[uBinSz] =
+ new TH1D(Form("binCntDist_Old_%5.0f", vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin, Old Bmon; Digis []; Bins []", vdBinSizesNs[uBinSz]), //
+ 10000, -0.5, 9999.5);
+ vBinCountDistributionNew[uBinSz] =
+ new TH1D(Form("binCntDist_New_%5.0f", vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin, New Bmon; Digis []; Bins []", vdBinSizesNs[uBinSz]), //
+ 10000, -0.5, 9999.5);
+ }
+
+ std::vector<uint32_t> vuIdxHistoCycleinTS(vdIntegrationNs.size(), 0);
+ uint32_t uTotalCountOld = 0;
+ uint32_t uTotalCountNew = 0;
+
+ /// Loop on timeslices
+ for (Int_t iEntry = dStartTs; iEntry < nentries; iEntry++) {
+ tree->GetEntry(iEntry);
+ uint32_t nDigisBmon = vDigisBmon->size();
+
+ if (nDigisBmon < 500) {
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => spill break, skipping this TS"
+ << std::endl;
+ continue;
+ }
+
+ std::vector<CbmBmonDigi> vDigisOld;
+ std::vector<CbmBmonDigi> vDigisNew;
+ vDigisOld.reserve(nDigisBmon);
+ vDigisNew.reserve(nDigisBmon);
+ for (auto& digi : *vDigisBmon) {
+ if (1 == CbmTofAddress::GetChannelSide(digi.GetAddress())) {
+ if (CbmTofAddress::GetChannelId(digi.GetAddress()) < 4) {
+ vDigisNew.push_back(digi);
+ }
+ else {
+ LOG(fatal) << "Bad sCVD channel: " << CbmTofAddress::GetChannelId(digi.GetAddress());
+ }
+ }
+ else {
+ vDigisOld.push_back(digi);
+ }
+ }
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => " << vDigisOld.size() << " Old + "
+ << vDigisNew.size() << " sCVD" << std::endl;
+ uTotalCountOld += vDigisOld.size();
+ uTotalCountNew += vDigisNew.size();
+
+ vuIdxHistoCycleinTS.assign(vdIntegrationNs.size(), 0);
+ uint32_t uDigiIdx = 0;
+ for (auto itOld = vDigisOld.begin(); itOld != vDigisOld.end(); ++itOld) {
+ double_t dTime = (*itOld).GetTime();
+ if (dTime < 0) {
+ std::cout << Form("TS %5d Old Digi %6u %8.0f!!!!!", iEntry, uDigiIdx, dTime) << std::endl;
+ continue;
+ }
+ if (dTsSizeNs * 1.01 < dTime) {
+ std::cout << Form("TS %5d Old Digi %6u %8.0f out of TS up [%5.3f]!!!!!", iEntry, uDigiIdx, dTime,
+ dTime / dTsSizeNs)
+ << std::endl;
+ uDigiIdx++;
+ continue;
+ }
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ uint32_t uCurrentCycle = std::floor(dTime / vdIntegrationNs[uHistSz]);
+ if (vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle) {
+ for (; vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle; ++vuIdxHistoCycleinTS[uHistSz]) {
+ double_t dTsFractional =
+ (vdIntegrationNs[uHistSz] * vuIdxHistoCycleinTS[uHistSz]) / dTsSizeNs + iEntry - dStartTs;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ double_t dQFactor = ExtractQFactor(vHistoOld[uHistSz][uBinSz]);
+ vQvalOld[uHistSz][uBinSz]->Fill(dTsFractional, dQFactor);
+ vMeanOld[uHistSz][uBinSz]->Fill(dTsFractional, ExtractMean(vHistoOld[uHistSz][uBinSz]));
+ /*
+ if (1280 == vdBinSizesNs[uBinSz]) {
+ std::cout << Form("%9f %9.0f %9.0f %9d %9.0f",
+ dTsFractional, vdIntegrationNs[uHistSz], vHistoOld[uHistSz][uBinSz]->Integral(),
+ vHistoOld[uHistSz][uBinSz]->GetNbinsX(),
+ vHistoOld[uHistSz][uBinSz]->GetBinContent(vHistoOld[uHistSz][uBinSz]->GetMaximumBin()))
+ << std::endl;
+ }
+ */
+ for (uint32_t uBin = 1; uBin <= vHistoOld[uHistSz][uBinSz]->GetNbinsX(); ++uBin) {
+ vBinCountDistributionOld[uBinSz]->Fill(vHistoOld[uHistSz][uBinSz]->GetBinContent(uBin));
+ }
+
+ if (0.0 < dQFactor) {
+ vHistoOld[uHistSz][uBinSz]->Reset();
+ }
+ }
+ }
+ }
+ }
+
+ double_t dTimeInCycle = std::fmod(dTime, vdIntegrationNs[uHistSz]);
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ vHistoOld[uHistSz][uBinSz]->Fill(dTimeInCycle);
+ }
+ }
+ }
+ uDigiIdx++;
+ }
+ /// FIXME: last "slice" is lost or properly take?
+
+ vuIdxHistoCycleinTS.assign(vdIntegrationNs.size(), 0);
+ uDigiIdx = 0;
+ for (auto itNew = vDigisNew.begin(); itNew != vDigisNew.end(); ++itNew) {
+ double_t dTime = (*itNew).GetTime();
+ if (dTime < 0) {
+ std::cout << Form("TS %5d New Digi %6u %8.0f!!!!!", iEntry, uDigiIdx, dTime) << std::endl;
+ continue;
+ }
+ if (dTsSizeNs * 1.01 < dTime) {
+ std::cout << Form("TS %5d New Digi %6u %8.0f out of TS up [%5.3f]!!!!!", iEntry, uDigiIdx, dTime,
+ dTime / dTsSizeNs)
+ << std::endl;
+ uDigiIdx++;
+ continue;
+ }
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ uint32_t uCurrentCycle = std::floor(dTime / vdIntegrationNs[uHistSz]);
+ if (vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle) {
+ for (; vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle; ++vuIdxHistoCycleinTS[uHistSz]) {
+ double_t dTsFractional =
+ (vdIntegrationNs[uHistSz] * vuIdxHistoCycleinTS[uHistSz]) / dTsSizeNs + iEntry - dStartTs;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalNew[uHistSz][uBinSz]) {
+ double_t dQFactor = ExtractQFactor(vHistoNew[uHistSz][uBinSz]);
+ vQvalNew[uHistSz][uBinSz]->Fill(dTsFractional, dQFactor);
+ vMeanNew[uHistSz][uBinSz]->Fill(dTsFractional, ExtractMean(vHistoNew[uHistSz][uBinSz]));
+
+ for (uint32_t uBin = 1; uBin <= vHistoNew[uHistSz][uBinSz]->GetNbinsX(); ++uBin) {
+ vBinCountDistributionNew[uBinSz]->Fill(vHistoNew[uHistSz][uBinSz]->GetBinContent(uBin));
+ }
+
+ if (0.0 < dQFactor) {
+ vHistoNew[uHistSz][uBinSz]->Reset();
+ }
+ }
+ }
+ }
+ }
+
+ double_t dTimeInCycle = std::fmod(dTime, vdIntegrationNs[uHistSz]);
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalNew[uHistSz][uBinSz]) {
+ vHistoNew[uHistSz][uBinSz]->Fill(dTimeInCycle);
+ }
+ }
+ }
+ uDigiIdx++;
+ }
+ /// FIXME: last "slice" is lost or properly take?
+ }
+
+ uint16_t nPadX = std::ceil(std::sqrt(uNbPlots));
+ uint16_t nPadY = std::ceil(1.0 * uNbPlots / nPadX);
+ std::cout << Form("Nb plots %3d Nb pads X %2d Y %2d", uNbPlots, nPadX, nPadY) << std::endl;
+
+ uint32_t uNbBinsLog = 0;
+ std::vector<double> dBinsLogVector = GenerateLogBinArray(4, 9, 1, uNbBinsLog, 2);
+ double* dBinsLog = dBinsLogVector.data();
+
+ TH1* hMeanVsBinSzOld = new TH1D(
+ Form("MeanVsBinSzOld_%04u", uRunId),
+ Form("Fit of the Mean vs Bin Size, Old, run %04u; Bin Size [ns]; Mean [Digis/bin]", uRunId), uNbBinsLog, dBinsLog);
+ TH1* hMeanVsBinSzNew = new TH1D(
+ Form("MeanVsBinSzNew_%04u", uRunId),
+ Form("Fit of the Mean vs Bin Size, New, run %04u; Bin Size [ns]; Mean [Digis/bin]", uRunId), uNbBinsLog, dBinsLog);
+ TH1* hQfactorVsBinSzOld =
+ new TH1D(Form("hQfactorVsBinSzOld_%04u", uRunId),
+ Form("Fit of Q-Factor vs Bin Size, run %04u; Bin Size [ns]; Q-Factor []", uRunId), uNbBinsLog, dBinsLog);
+ TH1* hQfactorVsBinSzNew =
+ new TH1D(Form("hQfactorVsBinSzNew_%04u", uRunId),
+ Form("Fit of Q-Factor vs Bin Size, run %04u; Bin Size [ns]; Q-Factor []", uRunId), uNbBinsLog, dBinsLog);
+
+ TCanvas* tempOld = new TCanvas("tempOld", "temp Old");
+ TCanvas* tempNew = new TCanvas("tempNew", "temp New");
+ tempOld->Divide(5, 4);
+ tempNew->Divide(5, 4);
+ uint32_t uPadIdx = 1;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ /// Fit Mean and Q-Factor
+ if (0 < vMeanOld[0][uBinSz]->GetEntries()) {
+ tempOld->cd(uPadIdx);
+ auto fCstMeanOld = new TF1("fCstMeanOld", "[0]", 0, numTimeslices);
+ vMeanOld[0][uBinSz]->Fit(fCstMeanOld, "", "", 10, 50);
+ hMeanVsBinSzOld->Fill(vdBinSizesNs[uBinSz], fCstMeanOld->GetParameter(0));
+
+ TCanvas temp("tempcanv", "temp");
+ temp.cd();
+ auto fCstQfactorOld = new TF1("fCstQfactorOld", "[0]", 0, numTimeslices);
+ vQvalOld[0][uBinSz]->Fit(fCstQfactorOld, "", "", 10, 50);
+ hQfactorVsBinSzOld->Fill(vdBinSizesNs[uBinSz], fCstQfactorOld->GetParameter(0));
+ delete fCstMeanOld;
+ delete fCstQfactorOld;
+ }
+
+ if (0 < vMeanNew[0][uBinSz]->GetEntries()) {
+ tempNew->cd(uPadIdx);
+ auto fCstMeanNew = new TF1("fCstMeanNew", "[0]", 0, numTimeslices);
+ vMeanNew[0][uBinSz]->Fit(fCstMeanNew, "", "", 10, 50);
+ hMeanVsBinSzNew->Fill(vdBinSizesNs[uBinSz], fCstMeanNew->GetParameter(0));
+
+ TCanvas temp("tempcanv", "temp");
+ temp.cd();
+ auto fCstQfactorNew = new TF1("fCstQfactorNew", "[0]", 0, numTimeslices);
+ vQvalNew[0][uBinSz]->Fit(fCstQfactorNew, "", "", 10, 50);
+ hQfactorVsBinSzNew->Fill(vdBinSizesNs[uBinSz], fCstQfactorNew->GetParameter(0));
+ delete fCstMeanNew;
+ delete fCstQfactorNew;
+ }
+ uPadIdx++;
+ }
+
+ TCanvas* test = new TCanvas("test", "test");
+ test->Divide(2, 2);
+
+ test->cd(1);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hMeanVsBinSzOld->Draw("hist");
+
+ test->cd(2);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hMeanVsBinSzNew->Draw("hist");
+
+ test->cd(3);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hQfactorVsBinSzOld->Draw("hist");
+
+ test->cd(4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hQfactorVsBinSzNew->Draw("hist");
+
+ TCanvas* cBinCntDistOld = new TCanvas("cBinCntDistOld", "BinCntDist Old");
+ TCanvas* cBinCntDistNew = new TCanvas("cBinCntDistNew", "BinCntDist New");
+ cBinCntDistOld->Divide(5, 4);
+ cBinCntDistNew->Divide(5, 4);
+ uPadIdx = 1;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ cBinCntDistOld->cd(uPadIdx);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogy();
+ vBinCountDistributionOld[uBinSz]->Draw("hist");
+
+ cBinCntDistNew->cd(uPadIdx);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogy();
+ vBinCountDistributionNew[uBinSz]->Draw("hist");
+
+ uPadIdx++;
+ }
+
+ TFile* outFile = new TFile(Form("data/bmon_q_factor_timescale_%04u_%03luTs.root", uRunId, numTimeslices), "RECREATE");
+ outFile->cd();
+ test->Write();
+ hMeanVsBinSzOld->Write();
+ hMeanVsBinSzNew->Write();
+ hQfactorVsBinSzOld->Write();
+ hQfactorVsBinSzNew->Write();
+ gROOT->cd();
+
+ outFile->Close();
+
+ outFile = new TFile(
+ Form("data/bmon_bin_counts_vs_bin_size_%04u_%03luTs_%03.0f_%03.0f.root", uRunId, numTimeslices, dStartTs, dStopTs),
+ "RECREATE");
+ outFile->cd();
+ test->Write();
+ cBinCntDistOld->Write();
+ cBinCntDistNew->Write();
+ gROOT->cd();
+
+ outFile->Close();
+
+ std::cout << "Beam particles in spill: " << uTotalCountOld << " (Old) VS " << uTotalCountNew << " (New)" << std::endl;
+ double_t dAvgCountSecOld = (1.0 * uTotalCountOld) / (numTimeslices * 0.128);
+ double_t dAvgCountSecNew = (1.0 * uTotalCountNew) / (numTimeslices * 0.128);
+ std::cout << "Avg Beam particles/s: " << dAvgCountSecOld << " (Old) VS " << dAvgCountSecNew << " (New)" << std::endl;
+}
diff --git a/macro/beamtime/mcbm2024/bmon_Q_Factor_timescale_plots.C b/macro/beamtime/mcbm2024/bmon_Q_Factor_timescale_plots.C
new file mode 100644
index 0000000000000000000000000000000000000000..96acbeda3f7c14a57f0180d8dd30b0c43709a516
--- /dev/null
+++ b/macro/beamtime/mcbm2024/bmon_Q_Factor_timescale_plots.C
@@ -0,0 +1,251 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer] */
+
+
+void bmon_Q_Factor_timescale_plots()
+{
+ std::vector<uint32_t> vRunId = {2391, 2984, 3006, 3107, 3109};
+ std::map<uint32_t, EColor> mColors = {{2391, EColor(kBlue)},
+ {2984, EColor(kRed)},
+ {3006, EColor(kGreen + 2)},
+ {3107, EColor(kViolet)},
+ {3109, EColor(kBlack)}};
+ uint32_t uRefRun = 2391;
+
+ std::map<uint32_t, TH1*> vMeanOld;
+ std::map<uint32_t, TH1*> vMeanNew;
+ std::map<uint32_t, TH1*> vQvalOld;
+ std::map<uint32_t, TH1*> vQvalNew;
+ std::map<uint32_t, TH1*> vMeanOldVsRef;
+ std::map<uint32_t, TH1*> vQvalOldVsRef;
+
+ std::map<uint32_t, TF1*> vMeanFitOld;
+ std::map<uint32_t, TF1*> vMeanFitNew;
+
+ for (auto runIdx : vRunId) {
+ TFile* inFile = new TFile(Form("data/bmon_q_factor_timescale_%04u_Spill.root", runIdx), "READ");
+ gROOT->cd();
+ std::cout << runIdx << std::endl;
+
+ TH1* pTemp = dynamic_cast<TH1*>(inFile->FindObjectAny(Form("MeanVsBinSzOld_%04u", runIdx)));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(Form("MeanVsBinSzOld_%04d_", runIdx)));
+ }
+ vMeanOld[runIdx] = pTemp;
+
+ pTemp = dynamic_cast<TH1*>(inFile->FindObjectAny(Form("MeanVsBinSzNew_%04u", runIdx)));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(Form("MeanVsBinSzNew_%04d_", runIdx)));
+ }
+ vMeanNew[runIdx] = pTemp;
+
+ pTemp = dynamic_cast<TH1*>(inFile->FindObjectAny(Form("hQfactorVsBinSzOld_%04u", runIdx)));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(Form("QfactorVsBinSzOld_%04d", runIdx)));
+ }
+ vQvalOld[runIdx] = pTemp;
+
+ pTemp = dynamic_cast<TH1*>(inFile->FindObjectAny(Form("hQfactorVsBinSzNew_%04u", runIdx)));
+ if (nullptr != pTemp) {
+ pTemp = dynamic_cast<TH1*>(pTemp->Clone(Form("QfactorVsBinSzOld_%04d", runIdx)));
+ }
+ vQvalNew[runIdx] = pTemp;
+
+ inFile->Close();
+ }
+
+ if (vMeanOld[uRefRun] && vQvalOld[uRefRun]) {
+ for (auto runIdx : vRunId) {
+ if (nullptr != vMeanOld[runIdx]) {
+ vMeanOldVsRef[runIdx] = dynamic_cast<TH1*>(vMeanOld[runIdx]->Clone(Form("MeanOldVsRef_%04d_", runIdx)));
+ vMeanOldVsRef[runIdx]->Divide(vMeanOld[uRefRun]);
+
+ vMeanFitOld[runIdx] = new TF1(Form("fMeanFitOld_%04u", runIdx), "pol1", 0, 1e9);
+ TCanvas temp("tempcanv", "temp");
+ temp.cd();
+ vMeanOld[runIdx]->Draw();
+ vMeanOld[runIdx]->Fit(vMeanFitOld[runIdx]);
+ }
+ if (nullptr != vQvalOld[runIdx]) {
+ vQvalOldVsRef[runIdx] = dynamic_cast<TH1*>(vQvalOld[runIdx]->Clone(Form("QvalOldVsRef_%04d_", runIdx)));
+ vQvalOldVsRef[runIdx]->Divide(vQvalOld[uRefRun]);
+ }
+ if (nullptr != vMeanNew[runIdx]) {
+ vMeanFitNew[runIdx] = new TF1(Form("fMeanFitNew_%04u", runIdx), "pol1", 0, 1e9);
+ TCanvas temp("tempcanv", "temp");
+ temp.cd();
+ vMeanNew[runIdx]->Draw();
+ vMeanNew[runIdx]->Fit(vMeanFitNew[runIdx]);
+ }
+ }
+ }
+
+ /// Timescale -----------------------------------------------------------------------------------------------------///
+ uint32_t uTimescaleHistSz = 4;
+ TCanvas* cTimescaleQfactor = new TCanvas("cTimescaleQfactor", "Mean & Q-factors vs timescale: 2.56 ms range");
+ cTimescaleQfactor->Divide(2, 2);
+
+ THStack* stackTimescaleMeanOld =
+ new THStack("stackBinSz_Mean_Old", "Mean bin content, Old BMon; Bin size [ns]; Mean bin content [digis]");
+ THStack* stackTimescaleMeanNew =
+ new THStack("stackBinSz_Mean_New", "Mean bin content, New BMon; Bin size [ns]; Mean bin content [digis]");
+ THStack* stackTimescaleQfactorOld =
+ new THStack("stackBinSz_Qfact_Old", "Q-Factor, Old BMon; Bin size [ns]; Q-Factor []");
+ THStack* stackTimescaleQfactorNew =
+ new THStack("stackBinSz_Qfact_New", "Q-Factor, New BMon; Bin size [ns]; Q-Factor []");
+
+ for (auto runIdx : vRunId) {
+ if (nullptr != vMeanOld[runIdx]) {
+ vMeanOld[runIdx]->SetLineColor(mColors[runIdx]);
+ vMeanOld[runIdx]->SetLineWidth(2);
+ stackTimescaleMeanOld->Add(vMeanOld[runIdx]);
+ }
+
+ if (nullptr != vMeanNew[runIdx]) {
+ vMeanNew[runIdx]->SetLineColor(mColors[runIdx]);
+ vMeanNew[runIdx]->SetLineWidth(2);
+ stackTimescaleMeanNew->Add(vMeanNew[runIdx]);
+ }
+
+ if (nullptr != vQvalOld[runIdx]) {
+ vQvalOld[runIdx]->SetLineColor(mColors[runIdx]);
+ vQvalOld[runIdx]->SetLineWidth(2);
+ stackTimescaleQfactorOld->Add(vQvalOld[runIdx]);
+ }
+
+ if (nullptr != vQvalNew[runIdx]) {
+ vQvalNew[runIdx]->SetLineColor(mColors[runIdx]);
+ vQvalNew[runIdx]->SetLineWidth(2);
+ stackTimescaleQfactorNew->Add(vQvalNew[runIdx]);
+ }
+ }
+
+ cTimescaleQfactor->cd(1);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ stackTimescaleMeanOld->Draw("hist nostack");
+ gPad->BuildLegend(0.15, 0.70, 0.45, 0.90, "");
+
+ cTimescaleQfactor->cd(2);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ stackTimescaleMeanNew->Draw("hist nostack");
+ gPad->BuildLegend(0.15, 0.70, 0.45, 0.90, "");
+
+ cTimescaleQfactor->cd(3);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ stackTimescaleQfactorOld->Draw("hist nostack");
+ gPad->BuildLegend(0.69, 0.70, 0.99, 0.90, "");
+
+ cTimescaleQfactor->cd(4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ stackTimescaleQfactorNew->Draw("hist nostack");
+ gPad->BuildLegend(0.69, 0.70, 0.99, 0.90, "");
+ ///----------------------------------------------------------------------------------------------------------------///
+
+ /// Reference -----------------------------------------------------------------------------------------------------///
+ TCanvas* cTimescaleRef = new TCanvas("cTimescaleRef", "Mean & Q-factors vs timescale: 2.56 ms range");
+ cTimescaleRef->Divide(2);
+
+ THStack* stackTimescaleRefMeanOld = new THStack(
+ "stackBinSzRef_Mean_Old",
+ Form("Mean bin content as fraction of run %04u, Old BMon; Bin size [ns]; Mean / Mean_%04u []", uRefRun, uRefRun));
+ THStack* stackTimescaleRefQfactorOld = new THStack(
+ "stackBinSzRef_Qfact_Old",
+ Form("Q-Factor as fraction of run %04u, Old BMon; Bin size [ns]; Q-Factor / Q-Factor)%04u []", uRefRun, uRefRun));
+
+ for (auto runIdx : vRunId) {
+ if (nullptr != vMeanOldVsRef[runIdx]) {
+ vMeanOldVsRef[runIdx]->SetLineColor(mColors[runIdx]);
+ vMeanOldVsRef[runIdx]->SetLineWidth(2);
+ stackTimescaleRefMeanOld->Add(vMeanOldVsRef[runIdx]);
+ }
+
+ if (nullptr != vQvalOldVsRef[runIdx]) {
+ vQvalOldVsRef[runIdx]->SetLineColor(mColors[runIdx]);
+ vQvalOldVsRef[runIdx]->SetLineWidth(2);
+ stackTimescaleRefQfactorOld->Add(vQvalOldVsRef[runIdx]);
+ }
+ }
+
+ cTimescaleRef->cd(1);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ stackTimescaleRefMeanOld->Draw("hist nostack");
+ gPad->BuildLegend(0.25, 0.20, 0.75, 0.40, "");
+
+ cTimescaleRef->cd(2);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ stackTimescaleRefQfactorOld->Draw("hist nostack");
+ gPad->BuildLegend(0.49, 0.60, 0.99, 0.90, "");
+ ///----------------------------------------------------------------------------------------------------------------///
+
+ /// Check of beam particles estimated counts ----------------------------------------------------------------------///
+ std::cout << "Fitted average flux per seconds: " << std::endl;
+ for (auto runIdx : vRunId) {
+ std::cout << Form("%04d => %9e Digi/s (Old) %9e Digi/s (New)", runIdx, vMeanFitOld[runIdx]->Eval(1e9),
+ vMeanFitNew[runIdx]->Eval(1e9))
+ << std::endl;
+ }
+ std::cout << "Fitted average count per spill: " << std::endl;
+ for (auto runIdx : vRunId) {
+ double_t dSpillDurationNs = 9.6e9; // ns
+ if (2391 == runIdx) {
+ dSpillDurationNs = 7.8e9; // ns
+ }
+ std::cout << Form("%04d => %9e Digis (Old) %9e Digis (New)", runIdx, vMeanFitOld[runIdx]->Eval(dSpillDurationNs),
+ vMeanFitNew[runIdx]->Eval(dSpillDurationNs))
+ << std::endl;
+ }
+ /*
+ ===============================================================
+ Run 2391 [38, 98] (61 TS)
+ Beam particles in spill: 35459706 (Old) VS 0 (New)
+ Avg Beam particles/s: 4.54146e+06 (Old) VS 0 (New)
+ ===============================================================
+ Run 2984 [66, 140] (75 TS, 141 is mostly spill ramp down)
+ Beam particles in spill: 27967704 (Old) VS 21882062 (New)
+ Avg Beam particles/s: 2.9133e+06 (Old) VS 2.27938e+06 (New)
+ ===============================================================
+ Run 3006 [76, 148] (73 TS, 149 is half empty, maybe sharp spill edge)
+ Beam particles in spill: 18387484 (Old) VS 29189274 (New)
+ Avg Beam particles/s: 1.96784e+06 (Old) VS 3.12385e+06 (New)
+ ===============================================================
+ Run 3107 [53, 128] (76 TS)
+ Beam particles in spill: 44813487 (Old) VS 50719284 (New)
+ Avg Beam particles/s: 4.60665e+06 (Old) VS 5.21374e+06 (New)
+ ===============================================================
+ Run 3109 [14, 88] (75 TS)
+ Beam particles in spill: 20582262 (Old) VS 24777405 (New)
+ Avg Beam particles/s: 2.14399e+06 (Old) VS 2.58098e+06 (New)
+ ===============================================================
+ Fitted average flux per seconds:
+ 2391 => 2.102520e+06 Digi/s (Old) 0.000000e+00 Digi/s (New)
+ 2984 => 1.212869e+06 Digi/s (Old) 1.066250e+06 Digi/s (New)
+ 3006 => 1.188376e+06 Digi/s (Old) 2.202738e+06 Digi/s (New)
+ 3107 => 1.761680e+06 Digi/s (Old) 2.133611e+06 Digi/s (New)
+ 3109 => 1.113414e+06 Digi/s (Old) 1.419566e+06 Digi/s (New)
+ Fitted average count per spill:
+ 2391 => 1.639965e+07 Digis (Old) 0.000000e+00 Digis (New)
+ 2984 => 1.164355e+07 Digis (Old) 1.023600e+07 Digis (New)
+ 3006 => 1.140841e+07 Digis (Old) 2.114629e+07 Digis (New)
+ 3107 => 1.691213e+07 Digis (Old) 2.048267e+07 Digis (New)
+ 3109 => 1.068877e+07 Digis (Old) 1.362784e+07 Digis (New)
+ */
+
+ ///----------------------------------------------------------------------------------------------------------------///
+}
diff --git a/macro/beamtime/mcbm2024/first_spill_bmon_bins.C b/macro/beamtime/mcbm2024/first_spill_bmon_bins.C
new file mode 100644
index 0000000000000000000000000000000000000000..ed0cb45cea76c3a40a2329402a40e71eff2ba011
--- /dev/null
+++ b/macro/beamtime/mcbm2024/first_spill_bmon_bins.C
@@ -0,0 +1,469 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer] */
+
+/** @brief Macro for check of bmon digis time spread depending on bin size in direct vector from legacy or online
+ * unpacking
+ ** @param input Name of input file
+ ** @param nTimeSlices Number of time-slices to process
+ ** @param dTsSizeNs Size of one Ts in ns = Total length of one time disribution histogram
+ **/
+void first_spill_bmon_bins(TString inputFileName, uint32_t uRunId, size_t numTimeslices = 0,
+ double_t dTsSizeNs = 128000000.0)
+{
+ gROOT->cd();
+ TFile* file = new TFile(inputFileName, "READ");
+ TTree* tree = (TTree*) (file->Get("cbmsim"));
+
+ std::vector<CbmBmonDigi>* vDigisBmon = new std::vector<CbmBmonDigi>();
+ tree->SetBranchAddress("BmonDigi", &vDigisBmon);
+
+ uint32_t nentries = tree->GetEntries();
+ cout << "Entries: " << nentries << endl;
+ nentries = (numTimeslices && numTimeslices < nentries ? numTimeslices : nentries);
+
+ gROOT->cd();
+
+
+ std::vector<double_t> vdBinSizes = {10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0,
+ 200, 400, 600, 800, 1e3, 2e3, 4e3, 6e3, 8e3, 1e4};
+ std::vector<int32_t> viNbBinsRange(vdBinSizes.size(), 1);
+
+ std::vector<std::vector<TH1*>> vDigiTimeDistBmonOld;
+ std::vector<std::vector<TH1*>> vDigiTimeDistBmonNew;
+ vDigiTimeDistBmonOld.resize(vdBinSizes.size());
+ vDigiTimeDistBmonNew.resize(vdBinSizes.size());
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ vDigiTimeDistBmonOld[uBinSz].resize(nentries);
+ vDigiTimeDistBmonNew[uBinSz].resize(nentries);
+ viNbBinsRange[uBinSz] = std::round(dTsSizeNs / vdBinSizes[uBinSz]);
+ }
+ TH2* vDigiTimeDistBmonOld_tenUs =
+ new TH2I("vDigiTimeDistBmonOld_tenUs", "Time distribution of Old Bmon digis; TS; time in TS [ns]; Nb Digis []",
+ nentries, 0, nentries, 12800, 0.0, dTsSizeNs);
+ TH2* vDigiTimeDistBmonNew_tenUs =
+ new TH2I("vDigiTimeDistBmonNew_tenUs", "Time distribution of Old Bmon digis; TS; time in TS [ns]; Nb Digis []",
+ nentries, 0, nentries, 12800, 0.0, dTsSizeNs);
+
+
+ /// Loop on timeslices
+ for (Int_t iEntry = 0; iEntry < nentries; iEntry++) {
+ //if (iEntry % 10 == 0 ) std::cout << "Entry " << iEntry << " / " << nentries << std::endl;
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ vDigiTimeDistBmonOld[uBinSz][iEntry] = new TH1I(
+ Form("hDistBmonOld_%03d_%03d", uBinSz, iEntry),
+ Form("Time distribution of Old Bmon digis in TS %3d, bin %.0f ns; Time in TS [ns]; Nb Digis []", iEntry,
+ vdBinSizes[uBinSz]), //
+ viNbBinsRange[uBinSz], 0.0, dTsSizeNs);
+ vDigiTimeDistBmonNew[uBinSz][iEntry] = new TH1I(
+ Form("hDistBmonNew_%03d_%03d", uBinSz, iEntry),
+ Form("Time distribution of New Bmon digis in TS %3d, bin %.0f ns; Time in TS [ns]; Nb Digis []", iEntry,
+ vdBinSizes[uBinSz]), //
+ viNbBinsRange[uBinSz], 0.0, dTsSizeNs);
+ }
+
+ tree->GetEntry(iEntry);
+ uint32_t nDigisBmon = vDigisBmon->size();
+
+ std::vector<CbmBmonDigi> vDigisOld;
+ std::vector<CbmBmonDigi> vDigisNew;
+ vDigisOld.reserve(nDigisBmon);
+ vDigisNew.reserve(nDigisBmon);
+ for (auto& digi : *vDigisBmon) {
+ if (1 == CbmTofAddress::GetChannelSide(digi.GetAddress())) {
+ if (CbmTofAddress::GetChannelId(digi.GetAddress()) < 4) {
+ vDigisNew.push_back(digi);
+ }
+ else {
+ LOG(fatal) << "Bad sCVD channel: " << CbmTofAddress::GetChannelId(digi.GetAddress());
+ }
+ }
+ else {
+ vDigisOld.push_back(digi);
+ }
+ }
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => " << vDigisOld.size() << " Old + "
+ << vDigisNew.size() << " sCVD" << std::endl;
+
+ for (auto itOld = vDigisOld.begin(); itOld != vDigisOld.end(); ++itOld) {
+ double_t dTime = (*itOld).GetTime();
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ vDigiTimeDistBmonOld[uBinSz][iEntry]->Fill(dTime);
+ }
+ vDigiTimeDistBmonOld_tenUs->Fill(iEntry, dTime);
+ }
+ for (auto itNew = vDigisNew.begin(); itNew != vDigisNew.end(); ++itNew) {
+ double_t dTime = (*itNew).GetTime();
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ vDigiTimeDistBmonNew[uBinSz][iEntry]->Fill(dTime);
+ }
+ vDigiTimeDistBmonNew_tenUs->Fill(iEntry, dTime);
+ }
+ }
+
+ /// Loop on bins in each timeslice
+ std::vector<TH1*> vhDigiMulOld(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhDigiMulNew(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhGapDurationNsOld(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhGapDurationNsNew(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhGapDurationUsOld(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhGapDurationUsNew(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhGapDurationMsOld(vdBinSizes.size(), nullptr);
+ std::vector<TH1*> vhGapDurationMsNew(vdBinSizes.size(), nullptr);
+
+ constexpr double epsilon = std::numeric_limits<double>::epsilon();
+ int32_t iNbBinsDurationBin0 =
+ std::floor(50000.0 / vdBinSizes[0]) + (epsilon < std::fabs(std::fmod(50000.0, vdBinSizes[0])) ? 1 : 0);
+
+ TH1* vhGapDurationNsOld_bin0 =
+ new TH1I("hGapDurationNsOld_bin0",
+ Form("Gaps duration, N bins * %3.0f ns bin, Old Bmon, %04d; Gap duration [ns]", vdBinSizes[0], uRunId), //
+ iNbBinsDurationBin0, 0, 50000.0 + (std::fmod(50000.0, vdBinSizes[0])));
+ TH1* vhGapDurationNsNew_bin0 =
+ new TH1I("hGapDurationNsNew_bin0",
+ Form("Gaps duration, N bins * %3.0f ns bin, New Bmon, %04d; Gap duration [ns]", vdBinSizes[0], uRunId), //
+ iNbBinsDurationBin0, 0, 50000.0 + (std::fmod(50000.0, vdBinSizes[0])));
+
+ TH2* vhGapDurationVsPrevNsOld_bin0 =
+ new TH2I("hGapDurationVsPrevNsOld_bin0",
+ Form("Gaps duration, N bins * %3.0f ns bin, Old Bmon, %04d; Prev Gap duration [ns]; Gap duration [ns]",
+ vdBinSizes[0], uRunId), //
+ iNbBinsDurationBin0 / 5, 0, 10000.0 + (std::fmod(10000.0, vdBinSizes[0])), //
+ iNbBinsDurationBin0 / 5, 0, 10000.0 + (std::fmod(10000.0, vdBinSizes[0])));
+ TH2* vhGapDurationVsPrevNsNew_bin0 =
+ new TH2I("hGapDurationVsPrevNsNew_bin0",
+ Form("Gaps duration, N bins * %3.0f ns bin, New Bmon, %04d; Prev Gap duration [ns];; Gap duration [ns]",
+ vdBinSizes[0], uRunId), //
+ iNbBinsDurationBin0 / 5, 0, 10000.0 + (std::fmod(10000.0, vdBinSizes[0])), //
+ iNbBinsDurationBin0 / 5, 0, 10000.0 + (std::fmod(10000.0, vdBinSizes[0])));
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ vhDigiMulOld[uBinSz] =
+ new TH1I(Form("hDigiMulOld_%03d", uBinSz),
+ Form("Digi Multiplicity per %3.0f ns bin, Old Bmon, %04d; Old [Digis]", vdBinSizes[uBinSz], uRunId), //
+ 30, -0.5, 29.5);
+ vhDigiMulNew[uBinSz] =
+ new TH1I(Form("hDigiMulNew_%03d", uBinSz),
+ Form("Digi Multiplicity per %3.0f ns bin, New Bmon, %04d; New [Digis]", vdBinSizes[uBinSz], uRunId), //
+ 30, -0.5, 29.5);
+
+ if (vdBinSizes[uBinSz] < 1000.0) {
+ int32_t iNbBinsDuration =
+ std::floor(1000.0 / vdBinSizes[uBinSz]) + (epsilon < std::fabs(std::fmod(1000.0, vdBinSizes[uBinSz])) ? 1 : 0);
+
+ vhGapDurationNsOld[uBinSz] = new TH1I(
+ Form("hGapDurationNsOld_%03d", uBinSz),
+ Form("Gaps duration, N bins * %3.0f ns bin, Old Bmon, %04d; Gap duration [ns]", vdBinSizes[uBinSz], uRunId), //
+ iNbBinsDuration, 0, 1000.0 + (std::fmod(1000.0, vdBinSizes[uBinSz])));
+ vhGapDurationNsNew[uBinSz] = new TH1I(
+ Form("hGapDurationNsNew_%03d", uBinSz),
+ Form("Gaps duration, N bins * %3.0f ns bin, New Bmon, %04d; Gap duration [ns]", vdBinSizes[uBinSz], uRunId), //
+ iNbBinsDuration, 0, 1000.0 + (std::fmod(1000.0, vdBinSizes[uBinSz])));
+ }
+ int32_t iNbBinsDuration = 2000;
+ double_t dDurationMax = 1e6;
+ if (1e3 <= vdBinSizes[uBinSz] && vdBinSizes[uBinSz] < 1e6) {
+ iNbBinsDuration =
+ std::floor(1e6 / vdBinSizes[uBinSz]) + (epsilon < std::fabs(std::fmod(1e6, vdBinSizes[uBinSz])) ? 1 : 0);
+ dDurationMax = 1e6 + (std::fmod(1e6, vdBinSizes[uBinSz]));
+ }
+
+ vhGapDurationUsOld[uBinSz] = new TH1I(
+ Form("hGapDurationUsOld_%03d", uBinSz),
+ Form("Gaps duration, N bins * %3.0f ns bin, Old Bmon, %04d; Gap duration [ns]", vdBinSizes[uBinSz], uRunId), //
+ iNbBinsDuration, 0, dDurationMax);
+ vhGapDurationUsNew[uBinSz] = new TH1I(
+ Form("hGapDurationUsNew_%03d", uBinSz),
+ Form("Gaps duration, N bins * %3.0f ns bin, New Bmon, %04d; Gap duration [ns]", vdBinSizes[uBinSz], uRunId), //
+ iNbBinsDuration, 0, dDurationMax);
+
+ vhGapDurationMsOld[uBinSz] = new TH1I(
+ Form("hGapDurationMsOld_%03d", uBinSz),
+ Form("Gaps duration, N bins * %3.0f ns bin, Old Bmon, %04d; Gap duration [ms]", vdBinSizes[uBinSz], uRunId), //
+ 200, 0, 200);
+ vhGapDurationMsNew[uBinSz] = new TH1I(
+ Form("hGapDuratioMsNew_%03d", uBinSz),
+ Form("Gaps duration, N bins * %3.0f ns bin, New Bmon, %04d; Gap duration [ms]", vdBinSizes[uBinSz], uRunId), //
+ 200, 0, 200);
+ }
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ for (Int_t iEntry = 0; iEntry < nentries; iEntry++) {
+ bool bLastEmptyOld = false;
+ bool bLastEmptyNew = false;
+ int32_t iNbEmptiesOld = 0;
+ int32_t iNbEmptiesNew = 0;
+ double_t dLastGapDurationOld = -1.0;
+ double_t dLastGapDurationNew = -1.0;
+ for (Int_t iBin = 0; iBin < viNbBinsRange[uBinSz]; iBin++) {
+ int32_t iCountsOld = vDigiTimeDistBmonOld[uBinSz][iEntry]->GetBinContent(iBin);
+ int32_t iCountsNew = vDigiTimeDistBmonNew[uBinSz][iEntry]->GetBinContent(iBin);
+
+ if (0 == iCountsOld) {
+ bLastEmptyOld = true;
+ iNbEmptiesOld++;
+ }
+ else {
+ vhDigiMulOld[uBinSz]->Fill(iCountsOld);
+ if (bLastEmptyOld) {
+
+ double_t dGapDuration = iNbEmptiesOld * vdBinSizes[uBinSz];
+ if (dGapDuration < 1e3 && vdBinSizes[uBinSz] < 1000.0) {
+ vhGapDurationNsOld[uBinSz]->Fill(dGapDuration);
+ }
+ if (dGapDuration < 1e6) {
+ vhGapDurationUsOld[uBinSz]->Fill(dGapDuration);
+ }
+ vhGapDurationMsOld[uBinSz]->Fill(dGapDuration / 1e6);
+
+ if (0 == uBinSz) {
+ vhGapDurationNsOld_bin0->Fill(dGapDuration);
+ if (0.0 < dLastGapDurationOld) {
+ vhGapDurationVsPrevNsOld_bin0->Fill(dLastGapDurationOld, dGapDuration);
+ }
+ dLastGapDurationOld = dGapDuration;
+ }
+
+ bLastEmptyOld = false;
+ iNbEmptiesOld = 0;
+ }
+ }
+
+ if (0 == iCountsNew) {
+ bLastEmptyNew = true;
+ iNbEmptiesNew++;
+ }
+ else {
+ vhDigiMulNew[uBinSz]->Fill(iCountsNew);
+ if (bLastEmptyNew) {
+ double_t dGapDuration = iNbEmptiesNew * vdBinSizes[uBinSz];
+ if (dGapDuration < 1e3) {
+ vhGapDurationNsNew[uBinSz]->Fill(dGapDuration);
+ }
+ if (dGapDuration < 1e6) {
+ vhGapDurationUsNew[uBinSz]->Fill(dGapDuration);
+ }
+ vhGapDurationMsNew[uBinSz]->Fill(dGapDuration / 1e6);
+
+ if (0 == uBinSz) {
+ vhGapDurationNsNew_bin0->Fill(dGapDuration);
+ if (0.0 < dLastGapDurationNew) {
+ vhGapDurationVsPrevNsNew_bin0->Fill(dLastGapDurationNew, dGapDuration);
+ }
+ dLastGapDurationNew = dGapDuration;
+ }
+
+ bLastEmptyNew = false;
+ iNbEmptiesNew = 0;
+ }
+ }
+ }
+ }
+ }
+ TCanvas* fCanvBinSizes =
+ new TCanvas("canvBinSizes", "Digi Multiplicity per bin, as function of bin size, Old and New Bmon");
+ fCanvBinSizes->Divide(2);
+
+ const TArrayI& palette = TColor::GetPalette();
+ int32_t iNbColors = palette.GetSize();
+
+ THStack* pStackOld =
+ new THStack("pStackOld", Form("Digi Multiplicity per bin, as function of bin size, Old Bmon, run %4d", uRunId));
+ THStack* pStackNew =
+ new THStack("pStackNew", Form("Digi Multiplicity per bin, as function of bin size, New Bmon, run %4d", uRunId));
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ int32_t iColor = iNbColors * (static_cast<double_t>(uBinSz) / vdBinSizes.size());
+ vhDigiMulOld[uBinSz]->SetLineColor(palette[iColor]);
+ vhDigiMulOld[uBinSz]->SetLineWidth(2);
+ pStackOld->Add(vhDigiMulOld[uBinSz]);
+
+ vhDigiMulNew[uBinSz]->SetLineColor(palette[iColor]);
+ vhDigiMulNew[uBinSz]->SetLineWidth(2);
+ pStackNew->Add(vhDigiMulNew[uBinSz]);
+
+ std::cout << Form("Bin %5.0f => Mean Old %5.2f New %5.2f", vdBinSizes[uBinSz], vhDigiMulOld[uBinSz]->GetMean(),
+ vhDigiMulNew[uBinSz]->GetMean())
+ << std::endl;
+ }
+
+ fCanvBinSizes->cd(1);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackOld->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ fCanvBinSizes->cd(2);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackNew->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ TCanvas* fCanvGapsDuration =
+ new TCanvas("canvGapsDuration", "Gaps durations between digis, as function of bin size, Old and New Bmon");
+ fCanvGapsDuration->Divide(3, 2);
+
+ THStack* pStackGapsNsOld =
+ new THStack("pStackGapsNsOld",
+ Form("Gaps durations between digis, as function of bin size, ns scale, Old Bmon, run %4d", uRunId));
+ THStack* pStackGapsNsNew =
+ new THStack("pStackGapsNsNew",
+ Form("Gaps durations between digis, as function of bin size, ns scale, New Bmon, run %4d", uRunId));
+
+ THStack* pStackGapsUsOld =
+ new THStack("pStackGapsUsOld",
+ Form("Gaps durations between digis, as function of bin size, us scale, Old Bmon, run %4d", uRunId));
+ THStack* pStackGapsUsNew =
+ new THStack("pStackGapsUsNew",
+ Form("Gaps durations between digis, as function of bin size, us scale, New Bmon, run %4d", uRunId));
+
+ THStack* pStackGapsMsOld =
+ new THStack("pStackGapMsOld",
+ Form("Gaps durations between digis, as function of bin size, ms scale, Old Bmon, run %4d", uRunId));
+ THStack* pStackGapsMsNew =
+ new THStack("pStackGapsMsNew",
+ Form("Gaps durations between digis, as function of bin size, ms scale, New Bmon, run %4d", uRunId));
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizes.size(); ++uBinSz) {
+ int32_t iColor = iNbColors * (static_cast<double_t>(uBinSz) / vdBinSizes.size());
+
+ if (vdBinSizes[uBinSz] < 1000.0) {
+ vhGapDurationNsOld[uBinSz]->SetLineColor(palette[iColor]);
+ vhGapDurationNsOld[uBinSz]->SetLineWidth(2);
+ pStackGapsNsOld->Add(vhGapDurationNsOld[uBinSz]);
+
+ vhGapDurationNsNew[uBinSz]->SetLineColor(palette[iColor]);
+ vhGapDurationNsNew[uBinSz]->SetLineWidth(2);
+ pStackGapsNsNew->Add(vhGapDurationNsNew[uBinSz]);
+ }
+
+ vhGapDurationUsOld[uBinSz]->SetLineColor(palette[iColor]);
+ vhGapDurationUsOld[uBinSz]->SetLineWidth(2);
+ pStackGapsUsOld->Add(vhGapDurationUsOld[uBinSz]);
+
+ vhGapDurationUsNew[uBinSz]->SetLineColor(palette[iColor]);
+ vhGapDurationUsNew[uBinSz]->SetLineWidth(2);
+ pStackGapsUsNew->Add(vhGapDurationUsNew[uBinSz]);
+
+ vhGapDurationMsOld[uBinSz]->SetLineColor(palette[iColor]);
+ vhGapDurationMsOld[uBinSz]->SetLineWidth(2);
+ pStackGapsMsOld->Add(vhGapDurationMsOld[uBinSz]);
+
+ vhGapDurationMsNew[uBinSz]->SetLineColor(palette[iColor]);
+ vhGapDurationMsNew[uBinSz]->SetLineWidth(2);
+ pStackGapsMsNew->Add(vhGapDurationMsNew[uBinSz]);
+ }
+
+ fCanvGapsDuration->cd(1);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackGapsNsOld->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ fCanvGapsDuration->cd(4);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackGapsNsNew->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ fCanvGapsDuration->cd(2);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackGapsUsOld->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ fCanvGapsDuration->cd(5);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackGapsUsNew->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ fCanvGapsDuration->cd(3);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackGapsMsOld->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+ fCanvGapsDuration->cd(6);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ pStackGapsMsNew->Draw("hist nostack");
+ gPad->BuildLegend(0.79, 0.79, 0.99, 0.94, "");
+
+
+ TCanvas* fCanvGapsDurationBin0 =
+ new TCanvas("canvGapsDurationBin0", "Gaps durations between digis, Old and New Bmon");
+ fCanvGapsDurationBin0->cd(1);
+ gPad->SetLogy();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vhGapDurationNsOld_bin0->SetLineColor(kBlue);
+ vhGapDurationNsNew_bin0->SetLineColor(kRed);
+ vhGapDurationNsOld_bin0->Draw("hist");
+ vhGapDurationNsNew_bin0->Draw("hist same");
+
+ TCanvas* fCanvGapsDurationVsPrevBin0 =
+ new TCanvas("canvGapsDurationVsPrevBin0", "Gaps durations between digis, Old and New Bmon");
+ fCanvGapsDurationVsPrevBin0->Divide(2, 2);
+
+ fCanvGapsDurationVsPrevBin0->cd(1);
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vhGapDurationVsPrevNsOld_bin0->Draw("colz");
+
+ fCanvGapsDurationVsPrevBin0->cd(2);
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vhGapDurationVsPrevNsNew_bin0->Draw("colz");
+
+ fCanvGapsDurationVsPrevBin0->cd(3);
+ gPad->SetLogx();
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vhGapDurationVsPrevNsOld_bin0->Draw("colz");
+
+ fCanvGapsDurationVsPrevBin0->cd(4);
+ gPad->SetLogx();
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vhGapDurationVsPrevNsNew_bin0->Draw("colz");
+
+ TCanvas* fCanvDigiTimeDist_tenUs =
+ new TCanvas("fCanvDigiTimeDist_tenUs", "Digi time dist, 10 us bining, per TS, Old and New Bmon");
+ fCanvDigiTimeDist_tenUs->Divide(2);
+
+ fCanvDigiTimeDist_tenUs->cd(1);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vDigiTimeDistBmonOld_tenUs->Draw("colz");
+
+ fCanvDigiTimeDist_tenUs->cd(2);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ vDigiTimeDistBmonNew_tenUs->Draw("colz");
+
+ TFile* outFile = new TFile(Form("data/first_spill_bmon_bins_%04u_%03luTs.root", uRunId, numTimeslices), "RECREATE");
+ outFile->cd();
+ fCanvBinSizes->Write();
+ fCanvGapsDuration->Write();
+ fCanvGapsDurationBin0->Write();
+ fCanvGapsDurationVsPrevBin0->Write();
+
+ gROOT->cd();
+
+ outFile->Close();
+}
diff --git a/macro/beamtime/mcbm2024/first_spill_bunches.C b/macro/beamtime/mcbm2024/first_spill_bunches.C
index 3069b2a4774e431ce45a69fa0ae8a00cf43db982..24be04d143fd33b1aab902f3824328df6382b17f 100644
--- a/macro/beamtime/mcbm2024/first_spill_bunches.C
+++ b/macro/beamtime/mcbm2024/first_spill_bunches.C
@@ -8,9 +8,8 @@
** @param dBinSizeNs Size of one bin in the time distribution in ns
** @param dTsSizeNs Size of one Ts in ns = Total length of one time disribution histogram
**/
- /// TODO: Add plots for individual STS and TOF "stations"
-void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimeslices = 0,
- double_t dBinSizeNs = 100.0,
+/// TODO: Add plots for individual STS and TOF "stations"
+void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimeslices = 0, double_t dBinSizeNs = 100.0,
double_t dZoomOffset = 60000000.0, double_t dZoomLength = 10000.0,
double_t dTsSizeNs = 128000000.0)
{
@@ -27,6 +26,12 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
std::vector<CbmTofDigi>* vDigisTof = new std::vector<CbmTofDigi>();
tree->SetBranchAddress("TofDigi", &vDigisTof);
+ std::vector<CbmTrdDigi>* vDigisTrd = new std::vector<CbmTrdDigi>();
+ tree->SetBranchAddress("TrdDigi", &vDigisTrd);
+
+ std::vector<CbmRichDigi>* vDigisRich = new std::vector<CbmRichDigi>();
+ tree->SetBranchAddress("RichDigi", &vDigisRich);
+
uint32_t nentries = tree->GetEntries();
cout << "Entries: " << nentries << endl;
nentries = (numTimeslices && numTimeslices < nentries ? numTimeslices : nentries);
@@ -34,37 +39,37 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
gROOT->cd();
int32_t iNbBinsRange = std::round(dTsSizeNs / dBinSizeNs);
- int32_t iNbBinsZoom = std::round(dZoomLength / dBinSizeNs);
+ int32_t iNbBinsZoom = std::round(dZoomLength / dBinSizeNs);
TH1* hDigiTimeDistZoomBmonOld = new TH1I(
Form("hDistZoomBmonOld_%04d", uRunId),
Form("Time distribution of Old Bmon digis in TS 0 of %4d, offset %.0f us bin %.0f ns; Time in TS [ns]; Nb Digis []",
- uRunId, dZoomOffset/1e3, dBinSizeNs),
- iNbBinsZoom, 0.0, dZoomLength);
+ uRunId, dZoomOffset / 1e3, dBinSizeNs), //
+ iNbBinsZoom, 0.0, dZoomLength);
TH1* hDigiTimeDistZoomBmonNew = new TH1I(
Form("hDistZoomBmonNew_%04d", uRunId),
Form("Time distribution of New Bmon digis in TS 0 of %4d, offset %.0f us bin %.0f ns; Time in TS [ns]; Nb Digis []",
- uRunId, dZoomOffset/1e3, dBinSizeNs),
+ uRunId, dZoomOffset / 1e3, dBinSizeNs), //
iNbBinsZoom, 0.0, dZoomLength);
TH1* hDigiTimeDistZoomStsSta0 = new TH1I(
Form("hDistZoomStsSta0_%04d", uRunId),
Form("Time distribution of Sts St 0 digis in TS 0 of %4d, offset %.0f us bin %.0f ns; Time in TS [ns]; Nb Digis []",
- uRunId, dZoomOffset/1e3, dBinSizeNs),
+ uRunId, dZoomOffset / 1e3, dBinSizeNs), //
iNbBinsZoom, 0.0, dZoomLength);
TH1* hDigiTimeDistZoomStsSta1 = new TH1I(
Form("hDistZoomStsSta1_%04d", uRunId),
Form("Time distribution of Sts St 1 digis in TS 0 of %4d, offset %.0f us bin %.0f ns; Time in TS [ns]; Nb Digis []",
- uRunId, dZoomOffset/1e3, dBinSizeNs),
+ uRunId, dZoomOffset / 1e3, dBinSizeNs), //
iNbBinsZoom, 0.0, dZoomLength);
TH1* hDigiTimeDistZoomStsSta2 = new TH1I(
Form("hDistZoomStsSta2_%04d", uRunId),
Form("Time distribution of Sts St 2 digis in TS 0 of %4d, offset %.0f us bin %.0f ns; Time in TS [ns]; Nb Digis []",
- uRunId, dZoomOffset/1e3, dBinSizeNs),
+ uRunId, dZoomOffset / 1e3, dBinSizeNs), //
iNbBinsZoom, 0.0, dZoomLength);
TH1* hDigiTimeDistZoomTofSta0 = new TH1I(
Form("hDistZoomTofSta0_%04d", uRunId),
Form("Time distribution of Tof St 0 digis in TS 0 of %4d, offset %.0f us bin %.0f ns; Time in TS [ns]; Nb Digis []",
- uRunId, dZoomOffset/1e3, dBinSizeNs),
+ uRunId, dZoomOffset / 1e3, dBinSizeNs), //
iNbBinsZoom, 0.0, dZoomLength);
std::vector<TH1*> vDigiTimeDistBmonOld(nentries, nullptr);
@@ -72,26 +77,58 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
std::vector<TH1*> vDigiTimeDistSts(nentries, nullptr);
std::vector<TH1*> vDigiTimeDistTof(nentries, nullptr);
+ std::vector<TH1*> vDigiTimeDistStsSta0(nentries, nullptr);
+ std::vector<TH1*> vDigiTimeDistStsSta1(nentries, nullptr);
+ std::vector<TH1*> vDigiTimeDistStsSta2(nentries, nullptr);
+
+ std::vector<TH1*> vDigiTimeDistTrd(nentries, nullptr);
+ std::vector<TH1*> vDigiTimeDistRich(nentries, nullptr);
+
/// Loop on timeslices
for (Int_t iEntry = 0; iEntry < nentries; iEntry++) {
//if (iEntry % 10 == 0 ) std::cout << "Entry " << iEntry << " / " << nentries << std::endl;
- vDigiTimeDistBmonOld[iEntry] = new TH1I(
- Form("hDistBmonOld_%03d", iEntry),
- Form("Time distribution of Old Bmon digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry),
- iNbBinsRange, 0.0, dTsSizeNs);
- vDigiTimeDistBmonNew[iEntry] = new TH1I(
- Form("hDistBmonNew_%03d", iEntry),
- Form("Time distribution of New Bmon digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry),
- iNbBinsRange, 0.0, dTsSizeNs);
- vDigiTimeDistSts[iEntry] = new TH1I(
- Form("hDistSts_%03d", iEntry),
- Form("Time distribution of Sts digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry),
- iNbBinsRange, 0.0, dTsSizeNs);
- vDigiTimeDistTof[iEntry] = new TH1I(
- Form("hDistTof_%03d", iEntry),
- Form("Time distribution of Tof digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry),
- iNbBinsRange, 0.0, dTsSizeNs);
+ vDigiTimeDistBmonOld[iEntry] =
+ new TH1I(Form("hDistBmonOld_%03d", iEntry),
+ Form("Time distribution of Old Bmon digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+ vDigiTimeDistBmonNew[iEntry] =
+ new TH1I(Form("hDistBmonNew_%03d", iEntry),
+ Form("Time distribution of New Bmon digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+ vDigiTimeDistSts[iEntry] =
+ new TH1I(Form("hDistSts_%03d", iEntry),
+ Form("Time distribution of Sts digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+ vDigiTimeDistTof[iEntry] =
+ new TH1I(Form("hDistTof_%03d", iEntry),
+ Form("Time distribution of Tof digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+
+ vDigiTimeDistStsSta0[iEntry] =
+ new TH1I(Form("hDistStsSta0_%03d", iEntry),
+ Form("Time distribution of Sts digis in TS %3d, Station 0; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+
+ vDigiTimeDistStsSta1[iEntry] =
+ new TH1I(Form("hDistStsSta1_%03d", iEntry),
+ Form("Time distribution of Sts digis in TS %3d, Station 1; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+
+ vDigiTimeDistStsSta2[iEntry] =
+ new TH1I(Form("hDistStsSta2_%03d", iEntry),
+ Form("Time distribution of Sts digis in TS %3d, Station 2; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+
+ vDigiTimeDistTrd[iEntry] =
+ new TH1I(Form("hDistTrd_%03d", iEntry),
+ Form("Time distribution of Trd digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
+
+ vDigiTimeDistRich[iEntry] =
+ new TH1I(Form("hDistRich_%03d", iEntry),
+ Form("Time distribution of Rich digis in TS %3d; Time in TS [ns]; Nb Digis []", iEntry), //
+ iNbBinsRange, 0.0, dTsSizeNs);
tree->GetEntry(iEntry);
@@ -120,16 +157,39 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
<< vDigisNew.size() << " sCVD, STS Digis: " << nDigisSts << ", TOF Digis: " << nDigisTof << std::endl;
for (auto itOld = vDigisOld.begin(); itOld != vDigisOld.end(); ++itOld) {
- vDigiTimeDistBmonOld[iEntry]->Fill( (*itOld).GetTime() );
+ vDigiTimeDistBmonOld[iEntry]->Fill((*itOld).GetTime());
}
for (auto itNew = vDigisNew.begin(); itNew != vDigisNew.end(); ++itNew) {
- vDigiTimeDistBmonNew[iEntry]->Fill( (*itNew).GetTime() );
+ vDigiTimeDistBmonNew[iEntry]->Fill((*itNew).GetTime());
}
for (auto itSts = vDigisSts->begin(); itSts != vDigisSts->end(); ++itSts) {
- vDigiTimeDistSts[iEntry]->Fill( (*itSts).GetTime() );
+ auto timeInTs = (*itSts).GetTime();
+ vDigiTimeDistSts[iEntry]->Fill(timeInTs);
+
+ auto addr = (*itSts).GetAddress();
+ auto station = CbmStsAddress::GetElementId(addr, EStsElementLevel::kStsUnit);
+
+ if (0 == station) {
+ vDigiTimeDistStsSta0[iEntry]->Fill(timeInTs);
+ }
+ else if (1 == station) {
+ vDigiTimeDistStsSta1[iEntry]->Fill(timeInTs);
+ }
+ else if (2 == station) {
+ vDigiTimeDistStsSta2[iEntry]->Fill(timeInTs);
+ }
}
for (auto itTof = vDigisTof->begin(); itTof != vDigisTof->end(); ++itTof) {
- vDigiTimeDistTof[iEntry]->Fill( (*itTof).GetTime() );
+ vDigiTimeDistTof[iEntry]->Fill((*itTof).GetTime());
+ }
+
+ uint32_t nDigisTrd = vDigisTrd->size();
+ uint32_t nDigisRich = vDigisRich->size();
+ for (auto itTrd = vDigisTrd->begin(); itTrd != vDigisTrd->end(); ++itTrd) {
+ vDigiTimeDistTrd[iEntry]->Fill((*itTrd).GetTime());
+ }
+ for (auto itRich = vDigisRich->begin(); itRich != vDigisRich->end(); ++itRich) {
+ vDigiTimeDistRich[iEntry]->Fill((*itRich).GetTime());
}
if (0 == iEntry) {
@@ -139,7 +199,7 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
if (timeInTs < dZoomOffset) continue;
if (dZoomOffset + dZoomLength < timeInTs) break;
- hDigiTimeDistZoomBmonOld->Fill( timeInTs - dZoomOffset );
+ hDigiTimeDistZoomBmonOld->Fill(timeInTs - dZoomOffset);
}
for (auto itNew = vDigisNew.begin(); itNew != vDigisNew.end(); ++itNew) {
auto timeInTs = (*itNew).GetTime();
@@ -147,7 +207,7 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
if (timeInTs < dZoomOffset) continue;
if (dZoomOffset + dZoomLength < timeInTs) break;
- hDigiTimeDistZoomBmonNew->Fill( timeInTs - dZoomOffset );
+ hDigiTimeDistZoomBmonNew->Fill(timeInTs - dZoomOffset);
}
for (auto itSts = vDigisSts->begin(); itSts != vDigisSts->end(); ++itSts) {
auto timeInTs = (*itSts).GetTime();
@@ -155,17 +215,17 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
if (timeInTs < dZoomOffset) continue;
if (dZoomOffset + dZoomLength < timeInTs) break;
- auto addr = (*itSts).GetAddress();
+ auto addr = (*itSts).GetAddress();
auto station = CbmStsAddress::GetElementId(addr, EStsElementLevel::kStsUnit);
- if( 0 == station) {
- hDigiTimeDistZoomStsSta0->Fill( timeInTs - dZoomOffset );
+ if (0 == station) {
+ hDigiTimeDistZoomStsSta0->Fill(timeInTs - dZoomOffset);
}
- else if( 1 == station) {
- hDigiTimeDistZoomStsSta1->Fill( timeInTs - dZoomOffset );
+ else if (1 == station) {
+ hDigiTimeDistZoomStsSta1->Fill(timeInTs - dZoomOffset);
}
- else if( 2 == station) {
- hDigiTimeDistZoomStsSta2->Fill( timeInTs - dZoomOffset );
+ else if (2 == station) {
+ hDigiTimeDistZoomStsSta2->Fill(timeInTs - dZoomOffset);
}
}
for (auto itTof = vDigisTof->begin(); itTof != vDigisTof->end(); ++itTof) {
@@ -174,13 +234,13 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
if (timeInTs < dZoomOffset) continue;
if (dZoomOffset + dZoomLength < timeInTs) break;
- hDigiTimeDistZoomTofSta0->Fill( timeInTs - dZoomOffset );
+ hDigiTimeDistZoomTofSta0->Fill(timeInTs - dZoomOffset);
}
}
}
TCanvas* fCanvFirstTs = new TCanvas("canvFirstTs", "Digis time distribution first TS, Bmon Old/New + STS + TOF");
-/*
+ /*
THStack* pStackTimeDist = new THStack(
"pStackTimeDist",
Form("Digis time distribution first TS, zoomed %1.0f us at %1.0f us, run %4d", dZoomLength/1e3, dZoomOffset/1e3,
@@ -254,36 +314,87 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
/// Loop on bins in each timeslice
TH2* hDigiMulNewOld = new TH2I(
- "hDigiMulNewOld",
- Form("Digi Multiplicity per %3.0f ns bin, Old VS New Bmon, %04d; New [Digis]; Old [Digis]", dBinSizeNs, uRunId),
- 30, -0.5, 29.5, 30, -0.5, 29.5);
+ "hDigiMulNewOld",
+ Form("Digi Multiplicity per %3.0f ns bin, Old VS New Bmon, %04d; New [Digis]; Old [Digis]", dBinSizeNs, uRunId), //
+ 30, -0.5, 29.5, 30, -0.5, 29.5);
TH2* hDigiMulOldSts = new TH2I(
- "hDigiMulOldSts",
- Form("Digi Multiplicity per %3.0f ns bin, STS VS Old Bmon, %04d; Old [Digis]; STS [Digis]", dBinSizeNs, uRunId),
- 30, -0.5, 29.5, 1000, -0.5, 999.5);
+ "hDigiMulOldSts",
+ Form("Digi Multiplicity per %3.0f ns bin, STS VS Old Bmon, %04d; Old [Digis]; STS [Digis]", dBinSizeNs, uRunId), //
+ 30, -0.5, 29.5, 1000, -0.5, 999.5);
TH2* hDigiMulOldTof = new TH2I(
- "hDigiMulOldTof",
- Form("Digi Multiplicity per %3.0f ns bin, TOF VS Old Bmon, %04d; Old [Digis]; TOF [Digis]", dBinSizeNs, uRunId),
- 30, -0.5, 29.5, 500, -0.5, 499.5);
+ "hDigiMulOldTof",
+ Form("Digi Multiplicity per %3.0f ns bin, TOF VS Old Bmon, %04d; Old [Digis]; TOF [Digis]", dBinSizeNs, uRunId), //
+ 30, -0.5, 29.5, 500, -0.5, 499.5);
TH2* hDigiMulNewSts = new TH2I(
- "hDigiMulNewSts",
- Form("Digi Multiplicity per %3.0f ns bin, STS VS New Bmon, %04d; New [Digis]; STS [Digis]", dBinSizeNs, uRunId),
- 30, -0.5, 29.5, 1000, -0.5, 999.5);
+ "hDigiMulNewSts",
+ Form("Digi Multiplicity per %3.0f ns bin, STS VS New Bmon, %04d; New [Digis]; STS [Digis]", dBinSizeNs, uRunId), //
+ 30, -0.5, 29.5, 1000, -0.5, 999.5);
TH2* hDigiMulNewTof = new TH2I(
- "hDigiMulNewTof",
- Form("Digi Multiplicity per %3.0f ns bin, TOF VS New Bmon, %04d; New [Digis]; TOF [Digis]", dBinSizeNs, uRunId),
- 30, -0.5, 29.5, 500, -0.5, 499.5);
+ "hDigiMulNewTof",
+ Form("Digi Multiplicity per %3.0f ns bin, TOF VS New Bmon, %04d; New [Digis]; TOF [Digis]", dBinSizeNs, uRunId), //
+ 30, -0.5, 29.5, 500, -0.5, 499.5);
TH2* hDigiMulStsTof = new TH2I(
- "hDigiMulStsTof",
- Form("Digi Multiplicity per %3.0f ns bin, TOF VS STS, %04d; Sts [Digis]; TOF [Digis]", dBinSizeNs, uRunId),
- 100, -0.5, 999.5, 50, -0.5, 499.5);
+ "hDigiMulStsTof",
+ Form("Digi Multiplicity per %3.0f ns bin, TOF VS STS, %04d; Sts [Digis]; TOF [Digis]", dBinSizeNs, uRunId), //
+ 100, -0.5, 999.5, 50, -0.5, 499.5);
+ TH2* hDigiMulTofTrd = new TH2I(
+ "hDigiMulTofTrd",
+ Form("Digi Multiplicity per %3.0f ns bin, TRD VS TOF, %04d; TOF [Digis]; TRD [Digis]", dBinSizeNs, uRunId), //
+ 500, -0.5, 499.5, 500, -0.5, 499.5);
+ TH2* hDigiMulTofRich = new TH2I(
+ "hDigiMulTofRich",
+ Form("Digi Multiplicity per %3.0f ns bin, RICH VS TOF, %04d; TOF [Digis]; RICH [Digis]", dBinSizeNs, uRunId), //
+ 500, -0.5, 499.5, 300, -0.5, 299.5);
+
+ TH2* hDigiMulEvoBmonOld =
+ new TH2I(Form("hDigiMulEvoBmonOld_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Old Bmon digis in %4d; TS [ns]; Nb Digis per bin []",
+ dBinSizeNs, uRunId), //
+ nentries * 10, 0.0, nentries, 30, -0.5, 29.5);
+ TH2* hDigiMulEvoBmonNew =
+ new TH2I(Form("hDigiMulEvoBmonNew_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of New Bmon digis in %4d; TS [ns]; Nb Digis per bin []",
+ dBinSizeNs, uRunId), //
+ nentries * 10, 0.0, nentries, 30, -0.5, 29.5);
+ TH2* hDigiMulEvoStsSta0 =
+ new TH2I(Form("hDigiMulEvoStsSta0_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Sts St 0 digis in %4d; TS [ns]; Nb Digis per bin []",
+ dBinSizeNs, uRunId), //
+ nentries * 10, 0.0, nentries, 1000, -0.5, 999.5);
+ TH2* hDigiMulEvoStsSta1 =
+ new TH2I(Form("hDigiMulEvoStsSta1_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Sts St 1 digis in %4d; TS [ns]; Nb Digis per bin []",
+ dBinSizeNs, uRunId), //
+ nentries * 10, 0.0, nentries, 1000, -0.5, 999.5);
+ TH2* hDigiMulEvoStsSta2 =
+ new TH2I(Form("hDigiMulEvoStsSta2_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Sts St 2 digis in %4d; TS [ns]; Nb Digis per bin []",
+ dBinSizeNs, uRunId), //
+ nentries * 10, 0.0, nentries, 1000, -0.5, 999.5);
+ TH2* hDigiMulEvoTofSta0 =
+ new TH2I(Form("hDigiMulEvoTofSta0_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Tof St 0 digis in %4d; TS [ns]; Nb Digis per bin []",
+ dBinSizeNs, uRunId), //
+ nentries * 10, 0.0, nentries, 500, -0.5, 499.5);
+ TH2* hDigiMulEvoTrd =
+ new TH2I(Form("hDigiMulEvoTrd_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Trd digis in %4d; TS [ns]; Nb Digis per bin []", dBinSizeNs,
+ uRunId), //
+ nentries * 10, 0.0, nentries, 500, -0.5, 499.5);
+ TH2* hDigiMulEvoRich =
+ new TH2I(Form("hDigiMulEvoRich_%04d", uRunId),
+ Form("Digi Multiplicity per %.0f ns bin of Rich digis in %4d; TS [ns]; Nb Digis per bin []", dBinSizeNs,
+ uRunId), //
+ nentries * 10, 0.0, nentries, 300, -0.5, 299.5);
for (Int_t iEntry = 0; iEntry < nentries; iEntry++) {
for (Int_t iBin = 0; iBin < iNbBinsRange; iBin++) {
- Int_t iNbDigisNew = vDigiTimeDistBmonNew[iEntry]->GetBinContent(iBin);
- Int_t iNbDigisOld = vDigiTimeDistBmonOld[iEntry]->GetBinContent(iBin);
- Int_t iNbDigisSts = vDigiTimeDistSts[iEntry]->GetBinContent(iBin);
- Int_t iNbDigisTof = vDigiTimeDistTof[iEntry]->GetBinContent(iBin);
+ Int_t iNbDigisNew = vDigiTimeDistBmonNew[iEntry]->GetBinContent(iBin);
+ Int_t iNbDigisOld = vDigiTimeDistBmonOld[iEntry]->GetBinContent(iBin);
+ Int_t iNbDigisSts = vDigiTimeDistSts[iEntry]->GetBinContent(iBin);
+ Int_t iNbDigisTof = vDigiTimeDistTof[iEntry]->GetBinContent(iBin);
+ Int_t iNbDigisTrd = vDigiTimeDistTrd[iEntry]->GetBinContent(iBin);
+ Int_t iNbDigisRich = vDigiTimeDistRich[iEntry]->GetBinContent(iBin);
if (0 < iNbDigisNew || 0 < iNbDigisOld) {
hDigiMulNewOld->Fill(iNbDigisNew, iNbDigisOld);
}
@@ -302,19 +413,34 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
if (0 < iNbDigisSts && 0 < iNbDigisTof) {
hDigiMulStsTof->Fill(iNbDigisSts, iNbDigisTof);
}
+ if (0 < iNbDigisTrd && 0 < iNbDigisTof) {
+ hDigiMulTofTrd->Fill(iNbDigisTof, iNbDigisTrd);
+ }
+ if (0 < iNbDigisRich && 0 < iNbDigisTof) {
+ hDigiMulTofRich->Fill(iNbDigisTof, iNbDigisRich);
+ }
+
+ double_t dTsFractional = (dBinSizeNs * iBin) / dTsSizeNs + iEntry;
+ hDigiMulEvoBmonOld->Fill(dTsFractional, iNbDigisOld);
+ hDigiMulEvoBmonNew->Fill(dTsFractional, iNbDigisNew);
+ hDigiMulEvoStsSta0->Fill(dTsFractional, vDigiTimeDistStsSta0[iEntry]->GetBinContent(iBin));
+ hDigiMulEvoStsSta1->Fill(dTsFractional, vDigiTimeDistStsSta1[iEntry]->GetBinContent(iBin));
+ hDigiMulEvoStsSta2->Fill(dTsFractional, vDigiTimeDistStsSta2[iEntry]->GetBinContent(iBin));
+ hDigiMulEvoTofSta0->Fill(dTsFractional, iNbDigisTof);
+ hDigiMulEvoTrd->Fill(dTsFractional, iNbDigisTrd);
+ hDigiMulEvoRich->Fill(dTsFractional, iNbDigisRich);
}
}
- TCanvas* fCanvCorrels = new TCanvas(
- "canvCorrels",
- Form("Digi Multiplicity per %3.0f ns bin, New Bmon vs Old Bmon/STS/TOF", dBinSizeNs));
- fCanvCorrels->Divide(3, 2);
+ TCanvas* fCanvCorrels =
+ new TCanvas("canvCorrels", Form("Digi Multiplicity per %3.0f ns bin, New Bmon vs Old Bmon/STS/TOF", dBinSizeNs));
+ fCanvCorrels->Divide(4, 2);
fCanvCorrels->cd(1);
gPad->SetLogz();
gPad->SetGridx();
gPad->SetGridy();
hDigiMulNewOld->Draw("colz");
- TProfile * hDigiMulNewOldCorr = hDigiMulNewOld->ProfileX("hDigiMulNewOldCorr");
+ TProfile* hDigiMulNewOldCorr = hDigiMulNewOld->ProfileX("hDigiMulNewOldCorr");
hDigiMulNewOldCorr->Draw("same");
fCanvCorrels->cd(2);
@@ -322,7 +448,7 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
gPad->SetGridx();
gPad->SetGridy();
hDigiMulOldSts->Draw("colz");
- TProfile * hDigiMulOldStsCorr = hDigiMulOldSts->ProfileX("hDigiMulOldStsCorr");
+ TProfile* hDigiMulOldStsCorr = hDigiMulOldSts->ProfileX("hDigiMulOldStsCorr");
hDigiMulOldStsCorr->Draw("same");
fCanvCorrels->cd(3);
@@ -330,30 +456,130 @@ void first_spill_bunches(TString inputFileName, uint32_t uRunId, size_t numTimes
gPad->SetGridx();
gPad->SetGridy();
hDigiMulOldTof->Draw("colz");
- TProfile * hDigiMulOldTofCorr = hDigiMulOldTof->ProfileX("hDigiMulOldTofCorr");
+ TProfile* hDigiMulOldTofCorr = hDigiMulOldTof->ProfileX("hDigiMulOldTofCorr");
hDigiMulOldTofCorr->Draw("same");
- fCanvCorrels->cd(4);
+ fCanvCorrels->cd(5);
gPad->SetLogz();
gPad->SetGridx();
gPad->SetGridy();
hDigiMulStsTof->Draw("colz");
- TProfile * hDigiMulStsTofCorr = hDigiMulStsTof->ProfileX("hDigiMulStsTofCorr");
+ TProfile* hDigiMulStsTofCorr = hDigiMulStsTof->ProfileX("hDigiMulStsTofCorr");
hDigiMulStsTofCorr->Draw("same");
- fCanvCorrels->cd(5);
+ fCanvCorrels->cd(6);
gPad->SetLogz();
gPad->SetGridx();
gPad->SetGridy();
hDigiMulNewSts->Draw("colz");
- TProfile * hDigiMulNewStsCorr = hDigiMulNewSts->ProfileX("hDigiMulNewStsCorr");
+ TProfile* hDigiMulNewStsCorr = hDigiMulNewSts->ProfileX("hDigiMulNewStsCorr");
hDigiMulNewStsCorr->Draw("same");
- fCanvCorrels->cd(6);
+ fCanvCorrels->cd(7);
gPad->SetLogz();
gPad->SetGridx();
gPad->SetGridy();
hDigiMulNewTof->Draw("colz");
- TProfile * hDigiMulNewTofCorr = hDigiMulNewTof->ProfileX("hDigiMulNewTofCorr");
+ TProfile* hDigiMulNewTofCorr = hDigiMulNewTof->ProfileX("hDigiMulNewTofCorr");
hDigiMulNewTofCorr->Draw("same");
+
+ fCanvCorrels->cd(4);
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulTofTrd->Draw("colz");
+ TProfile* hDigiMulTofTrdCorr = hDigiMulTofTrd->ProfileX("hDigiMulTofTrdCorr");
+ hDigiMulTofTrdCorr->Draw("same");
+
+ fCanvCorrels->cd(8);
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulTofRich->Draw("colz");
+ TProfile* hDigiMulTofRichCorr = hDigiMulTofRich->ProfileX("hDigiMulTofRichCorr");
+ hDigiMulTofRichCorr->Draw("same");
+
+ TCanvas* fCanvEvo = new TCanvas("canvEvo", Form("Digi Multiplicity per %3.0f ns bin vs timeslice", dBinSizeNs));
+ fCanvEvo->Divide(4, 2);
+
+ fCanvEvo->cd(1);
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoBmonOld->Draw("colz");
+ TProfile* hDigiMulEvoBmonOldMean = hDigiMulEvoBmonOld->ProfileX("hDigiMulEvoBmonOldMean");
+ hDigiMulEvoBmonOldMean->Draw("same");
+
+ fCanvEvo->cd(2);
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoBmonNew->Draw("colz");
+ TProfile* hDigiMulEvoBmonNewMean = hDigiMulEvoBmonNew->ProfileX("hDigiMulEvoBmonNewMean");
+ hDigiMulEvoBmonNewMean->Draw("same");
+
+ fCanvEvo->cd(3);
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoStsSta0->Draw("colz");
+ TProfile* hDigiMulEvoStsSta0Mean = hDigiMulEvoStsSta0->ProfileX("hDigiMulEvoStsSta0Mean");
+ hDigiMulEvoStsSta0Mean->Draw("same");
+
+ fCanvEvo->cd(4);
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoStsSta1->Draw("colz");
+ TProfile* hDigiMulEvoStsSta1Mean = hDigiMulEvoStsSta1->ProfileX("hDigiMulEvoStsSta1Mean");
+ hDigiMulEvoStsSta1Mean->Draw("same");
+
+ fCanvEvo->cd(5);
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoStsSta2->Draw("colz");
+ TProfile* hDigiMulEvoStsSta2Mean = hDigiMulEvoStsSta2->ProfileX("hDigiMulEvoStsSta2Mean");
+ hDigiMulEvoStsSta2Mean->Draw("same");
+
+ fCanvEvo->cd(6);
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoTofSta0->Draw("colz");
+ TProfile* hDigiMulEvoTofSta0Mean = hDigiMulEvoTofSta0->ProfileX("hDigiMulEvoTofSta0Mean");
+ hDigiMulEvoTofSta0Mean->Draw("same");
+
+ fCanvEvo->cd(7);
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoTrd->Draw("colz");
+ TProfile* hDigiMulEvoTrdMean = hDigiMulEvoTrd->ProfileX("hDigiMulEvoTrdMean");
+ hDigiMulEvoTrdMean->Draw("same");
+
+ fCanvEvo->cd(8);
+ gPad->SetLogy();
+ gPad->SetLogz();
+ gPad->SetGridx();
+ gPad->SetGridy();
+ hDigiMulEvoRich->Draw("colz");
+ TProfile* hDigiMulEvoRichMean = hDigiMulEvoRich->ProfileX("hDigiMulEvoRichMean");
+ hDigiMulEvoRichMean->Draw("same");
+
+ TFile* outFile =
+ new TFile(Form("data/first_spill_bunches_%04u_%03luTs_%.0fNs.root", uRunId, numTimeslices, dBinSizeNs), "RECREATE");
+ outFile->cd();
+ fCanvFirstTs->Write();
+ fCanvCorrels->Write();
+ fCanvEvo->Write();
+
+ gROOT->cd();
+
+ outFile->Close();
}
diff --git a/macro/beamtime/mcbm2024/sts_tof_Q_Factor_timescale.C b/macro/beamtime/mcbm2024/sts_tof_Q_Factor_timescale.C
new file mode 100644
index 0000000000000000000000000000000000000000..0ae73885132bfc751f4bbd35dc6d7ba0b1165031
--- /dev/null
+++ b/macro/beamtime/mcbm2024/sts_tof_Q_Factor_timescale.C
@@ -0,0 +1,508 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau [committer] */
+
+
+std::vector<double> GenerateLogBinArray(uint32_t uNbDecadesLog, uint32_t uNbStepsDecade, uint32_t uNbSubStepsInStep,
+ uint32_t& uNbBinsLog, int32_t iStartExp = 0, bool bAddZeroStart = false)
+{
+ /// Logarithmic bining for self time comparison
+ /// Number of log bins =
+ /// 9 for the sub-unit decade
+ /// + 9 for each unit of each decade * 10 for the subdecade range
+ /// + 1 for the closing bin top edge
+ uNbBinsLog = uNbStepsDecade + uNbStepsDecade * uNbSubStepsInStep * uNbDecadesLog;
+
+ /// Need uNbBinsLog + 1 values as we need to provide the end of last bin
+ uint32_t uArrayLength = uNbBinsLog + 1;
+ double dBinsLog[uArrayLength];
+ /// First fill sub-unit decade
+ for (uint32_t uSubU = 0; uSubU < uNbStepsDecade; uSubU++) {
+ dBinsLog[uSubU] = std::pow(10, iStartExp - 1) * (1 + uSubU);
+ }
+
+ /// Then fill the main decades
+ double dSubstepSize = 1.0 / uNbSubStepsInStep;
+ for (uint32_t uDecade = 0; uDecade < uNbDecadesLog; uDecade++) {
+ double dBase = std::pow(10, iStartExp + static_cast<int32_t>(uDecade));
+ uint32_t uDecadeIdx = uNbStepsDecade + uDecade * uNbStepsDecade * uNbSubStepsInStep;
+ for (uint32_t uStep = 0; uStep < uNbStepsDecade; uStep++) {
+ uint32_t uStepIdx = uDecadeIdx + uStep * uNbSubStepsInStep;
+ for (uint32_t uSubStep = 0; uSubStep < uNbSubStepsInStep; uSubStep++) {
+ dBinsLog[uStepIdx + uSubStep] = dBase * (1 + uStep) + dBase * dSubstepSize * uSubStep;
+ } // for( uint32_t uSubStep = 0; uSubStep < uNbSubStepsInStep; uSubStep++ )
+ } // for( uint32_t uStep = 0; uStep < uNbStepsDecade; uStep++ )
+ } // for( uint32_t uDecade = 0; uDecade < uNbDecadesLog; uDecade ++)
+ dBinsLog[uNbBinsLog] = std::pow(10, iStartExp + uNbDecadesLog);
+
+ /// use vector instead
+ std::vector<double> dBinsLogVect;
+
+ /// + 1 optional if bin [ 0; Min [ should be added
+ if (bAddZeroStart) {
+ uNbBinsLog++;
+ dBinsLogVect.push_back(0);
+ }
+
+ for (uint32_t i = 0; i < uArrayLength; ++i) {
+ dBinsLogVect.push_back(dBinsLog[i]);
+ } // for( uint32_t i = 0; i < uArrayLength; ++i )
+
+ return dBinsLogVect;
+}
+
+double_t ExtractMean(TH1* pHistoIn, double_t dThreshold = 1.0)
+{
+ // Q-Factor = Max Bin Content / Mean Content of all bin in range
+ // => Tend toward 1 if bins are more identical
+ if (pHistoIn->Integral()) {
+ uint32_t uNbBins = pHistoIn->GetNbinsX();
+ uint32_t uNbNonEmpty = 0;
+ double_t dIntegralNonEmpty = 0.0;
+ for (uint32_t uBin = 1; uBin < uNbBins; ++uBin) {
+ double_t dBinCount = pHistoIn->GetBinContent(uBin);
+ if (dThreshold <= dBinCount) {
+ uNbNonEmpty++;
+ dIntegralNonEmpty += dBinCount;
+ }
+ }
+ return (dIntegralNonEmpty / uNbNonEmpty);
+ }
+ else {
+ return 0.0;
+ }
+}
+double_t ExtractQFactor(TH1* pHistoIn, double_t dThreshold = 1.0)
+{
+ // Q-Factor = Max Bin Content / Mean Content of all bin in range
+ // => Tend toward 1 if bins are more identical
+ if (pHistoIn->Integral()) {
+ return (pHistoIn->GetBinContent(pHistoIn->GetMaximumBin())) / ExtractMean(pHistoIn, dThreshold);
+ }
+ else {
+ return 0.0;
+ }
+}
+
+/** @brief Macro for check of Hades-like Q-Factor based on bmon digis, against variations of bin size and integration
+ * length
+ ** @param input Name of input file
+ ** @param nTimeSlices Number of time-slices to process
+ ** @param dTsSizeNs Size of one Ts in ns = Total length of one time disribution histogram
+ **/
+void sts_tof_Q_Factor_timescale(TString inputFileName, uint32_t uRunId, size_t numTimeslices = 0, double_t dStartTs = 0,
+ double_t dStopTs = 0, double_t dTsSizeNs = 128000000.0)
+{
+
+ gROOT->cd();
+ TFile* file = new TFile(inputFileName, "READ");
+ TTree* tree = (TTree*) (file->Get("cbmsim"));
+
+ std::vector<CbmBmonDigi>* vDigisBmon = new std::vector<CbmBmonDigi>();
+ tree->SetBranchAddress("BmonDigi", &vDigisBmon);
+
+ std::vector<CbmStsDigi>* vDigisSts = new std::vector<CbmStsDigi>();
+ tree->SetBranchAddress("StsDigi", &vDigisSts);
+
+ std::vector<CbmTofDigi>* vDigisTof = new std::vector<CbmTofDigi>();
+ tree->SetBranchAddress("TofDigi", &vDigisTof);
+
+ uint32_t nentries = tree->GetEntries();
+ cout << "Entries: " << nentries << endl;
+ nentries = (numTimeslices && numTimeslices < nentries ? numTimeslices : nentries);
+
+ if (0 < dStopTs && dStopTs < nentries) {
+ nentries = dStopTs + 1;
+ }
+ numTimeslices = nentries - dStartTs;
+
+ gROOT->cd();
+
+ uint32_t uMinNbBins = 10;
+ uint32_t uMaxNbBins = 300000;
+ /// Hint: keep fractions of TS size and under 100 us
+ std::vector<double_t> vdBinSizesNs = {40, 100, 1.28e3, 10.24e3}; // 14 values
+ /// Hint: keep fractions of TS size + multiples of bin size and above 10 us
+ std::vector<double_t> vdIntegrationNs = {2.56e6}; // 13 values
+ /// Dimension: same as BinSizes vector!!
+ std::vector<double_t> vdQfactorHistMax = {2000., 400., 40., 20.};
+
+ std::vector<std::vector<uint32_t>> vuNbBinsHisto(vdIntegrationNs.size(),
+ std::vector<uint32_t>(vdBinSizesNs.size(), 0));
+ std::vector<uint32_t> vuNbHistoCyclesPerTS(vdIntegrationNs.size(), 0);
+ std::vector<std::vector<TH1*>> vHistoOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vHistoNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vQvalOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vQvalNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vMeanOld(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+ std::vector<std::vector<TH1*>> vMeanNew(vdIntegrationNs.size(), std::vector<TH1*>(vdBinSizesNs.size(), nullptr));
+
+ std::vector<TH1*> vBinCountDistributionOld(vdBinSizesNs.size(), nullptr);
+ std::vector<TH1*> vBinCountDistributionNew(vdBinSizesNs.size(), nullptr);
+
+ uint16_t uNbPlots = 0;
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ /// Pre-check values before in spreadsheet to make sure integer !!!!
+ vuNbHistoCyclesPerTS[uHistSz] = dTsSizeNs / vdIntegrationNs[uHistSz];
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ /// Pre-check values before in spreadsheet to make sure integer !!!!
+ vuNbBinsHisto[uHistSz][uBinSz] = vdIntegrationNs[uHistSz] / vdBinSizesNs[uBinSz];
+ if (uMinNbBins <= vuNbBinsHisto[uHistSz][uBinSz] /*&& vuNbBinsHisto[uHistSz][uBinSz] <= uMaxNbBins*/) {
+ vHistoOld[uHistSz][uBinSz] =
+ new TH1D(Form("binHist_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin in cycle of range %9.0f ns, Old Bmon; Time in Cycle [ns]; Digis []",
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbBinsHisto[uHistSz][uBinSz], 0.0, vdIntegrationNs[uHistSz]);
+ vHistoNew[uHistSz][uBinSz] =
+ new TH1D(Form("binHist_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin in cycle of range %9.0f ns, New Bmon; Time in Cycle [ns]; Digis []",
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbBinsHisto[uHistSz][uBinSz], 0.0, vdIntegrationNs[uHistSz]);
+
+ double_t dBinOffset = 1.0 / (2.0 * vuNbHistoCyclesPerTS[uHistSz]);
+ vQvalOld[uHistSz][uBinSz] =
+ new TH1D(Form("evoQFactor_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, Old Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+ vQvalNew[uHistSz][uBinSz] =
+ new TH1D(Form("evoQFactor_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, New Bmon; Time in Run [TS]; Q Factor []",
+ uRunId, vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+
+ vMeanOld[uHistSz][uBinSz] =
+ new TH1D(Form("evoMean_Old_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, Old Bmon; Time in Run [TS]; Mean []", uRunId,
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+ vMeanNew[uHistSz][uBinSz] =
+ new TH1D(Form("evoMean_New_%9.0f_%5.0f", vdIntegrationNs[uHistSz], vdBinSizesNs[uBinSz]),
+ Form("Q Factor, run %4d, %5.0f ns bin, %9.0f ns range, New Bmon; Time in Run [TS]; Mean []", uRunId,
+ vdBinSizesNs[uBinSz], vdIntegrationNs[uHistSz]), //
+ vuNbHistoCyclesPerTS[uHistSz] * numTimeslices, 0.0 - dBinOffset, numTimeslices - dBinOffset);
+
+ uNbPlots++;
+ }
+ }
+ }
+
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ vBinCountDistributionOld[uBinSz] =
+ new TH1D(Form("binCntDist_Old_%5.0f", vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin, Old Bmon; Digis []; Bins []", vdBinSizesNs[uBinSz]), //
+ 10000, -0.5, 9999.5);
+ vBinCountDistributionNew[uBinSz] =
+ new TH1D(Form("binCntDist_New_%5.0f", vdBinSizesNs[uBinSz]),
+ Form("Counts per %5.0f ns bin, New Bmon; Digis []; Bins []", vdBinSizesNs[uBinSz]), //
+ 10000, -0.5, 9999.5);
+ }
+
+ std::vector<uint32_t> vuIdxHistoCycleinTS(vdIntegrationNs.size(), 0);
+ uint32_t uTotalCountOld = 0;
+ uint32_t uTotalCountNew = 0;
+
+ /// Loop on timeslices
+ for (Int_t iEntry = dStartTs; iEntry < nentries; iEntry++) {
+ tree->GetEntry(iEntry);
+ uint32_t nDigisBmon = vDigisBmon->size();
+ uint32_t nDigisSts = vDigisSts->size();
+ uint32_t nDigisTof = vDigisTof->size();
+
+
+ if (nDigisBmon < 500) {
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => spill break, skipping this TS"
+ << std::endl;
+ continue;
+ }
+
+ /*
+ std::vector<CbmBmonDigi> vDigisOld;
+ std::vector<CbmBmonDigi> vDigisNew;
+ vDigisOld.reserve(nDigisBmon);
+ vDigisNew.reserve(nDigisBmon);
+ for (auto& digi : *vDigisBmon) {
+ if (1 == CbmTofAddress::GetChannelSide(digi.GetAddress())) {
+ if (CbmTofAddress::GetChannelId(digi.GetAddress()) < 4) {
+ vDigisNew.push_back(digi);
+ }
+ else {
+ LOG(fatal) << "Bad sCVD channel: " << CbmTofAddress::GetChannelId(digi.GetAddress());
+ }
+ }
+ else {
+ vDigisOld.push_back(digi);
+ }
+ }
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " => " << vDigisOld.size() << " Old + "
+ << vDigisNew.size() << " sCVD" << std::endl;
+ uTotalCountOld += vDigisOld.size();
+ uTotalCountNew += vDigisNew.size();
+ */
+ std::cout << "TS " << iEntry << ", BMON Digis: " << nDigisBmon << " & " << nDigisSts << " Sts"
+ << " & " << nDigisTof << " TOF" << std::endl;
+
+ vuIdxHistoCycleinTS.assign(vdIntegrationNs.size(), 0);
+ uint32_t uDigiIdx = 0;
+ for (auto itOld = vDigisTof->begin(); itOld != vDigisTof->end(); ++itOld) {
+ double_t dTime = (*itOld).GetTime();
+ if (dTime < 0) {
+ std::cout << Form("TS %5d TOF Digi %6u %8.0f out of TS down!!!!!", iEntry, uDigiIdx, dTime) << std::endl;
+ uDigiIdx++;
+ continue;
+ }
+ if (dTsSizeNs * 1.01 < dTime) {
+ std::cout << Form("TS %5d TOF Digi %6u %8.0f out of TS up [%5.3f]!!!!!", iEntry, uDigiIdx, dTime,
+ dTime / dTsSizeNs)
+ << std::endl;
+ uDigiIdx++;
+ continue;
+ }
+ if (0 == uDigiIdx % 1000000) {
+ std::cout << Form("TS %5d TOF Digi %8u / %8u", iEntry, uDigiIdx, nDigisTof) << std::endl;
+ }
+
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ uint32_t uCurrentCycle = std::floor(dTime / vdIntegrationNs[uHistSz]);
+ if (vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle) {
+ for (; vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle; ++vuIdxHistoCycleinTS[uHistSz]) {
+ double_t dTsFractional =
+ (vdIntegrationNs[uHistSz] * vuIdxHistoCycleinTS[uHistSz]) / dTsSizeNs + iEntry - dStartTs;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ /*
+ std::cout << Form("Extracting TS %5d Digi %8u / %8u hist sz %2d bin sz %2u cycle %u / %u time %9.0f ns",
+ iEntry, uDigiIdx, nDigisTof, uHistSz, uBinSz, vuIdxHistoCycleinTS[uHistSz], uCurrentCycle, dTime)
+ << std::endl;
+ */
+
+ double_t dQFactor = ExtractQFactor(vHistoOld[uHistSz][uBinSz], 8);
+ vQvalOld[uHistSz][uBinSz]->Fill(dTsFractional, dQFactor);
+ vMeanOld[uHistSz][uBinSz]->Fill(dTsFractional, ExtractMean(vHistoOld[uHistSz][uBinSz], 8));
+ /*
+ if (1280 == vdBinSizesNs[uBinSz]) {
+ std::cout << Form("%9f %9.0f %9.0f %9d %9.0f",
+ dTsFractional, vdIntegrationNs[uHistSz], vHistoOld[uHistSz][uBinSz]->Integral(),
+ vHistoOld[uHistSz][uBinSz]->GetNbinsX(),
+ vHistoOld[uHistSz][uBinSz]->GetBinContent(vHistoOld[uHistSz][uBinSz]->GetMaximumBin()))
+ << std::endl;
+ }
+ */
+ for (uint32_t uBin = 1; uBin <= vHistoOld[uHistSz][uBinSz]->GetNbinsX(); ++uBin) {
+ vBinCountDistributionOld[uBinSz]->Fill(vHistoOld[uHistSz][uBinSz]->GetBinContent(uBin));
+ }
+
+ // if (0.0 < dQFactor) {
+ vHistoOld[uHistSz][uBinSz]->Reset();
+ // }
+ }
+ }
+ }
+ }
+
+ double_t dTimeInCycle = std::fmod(dTime, vdIntegrationNs[uHistSz]);
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalOld[uHistSz][uBinSz]) {
+ vHistoOld[uHistSz][uBinSz]->Fill(dTimeInCycle);
+ }
+ }
+ }
+ uDigiIdx++;
+ }
+ /// FIXME: last "slice" is lost or properly take?
+ vuIdxHistoCycleinTS.assign(vdIntegrationNs.size(), 0);
+ uDigiIdx = 0;
+ for (auto itNew = vDigisSts->begin(); itNew != vDigisSts->end(); ++itNew) {
+ double_t dTime = (*itNew).GetTime();
+ if (dTime < 0) {
+ std::cout << Form("TS %5d Sts Digi %6u %8.0f!!!!!", iEntry, uDigiIdx, dTime) << std::endl;
+ continue;
+ }
+ if (dTsSizeNs * 1.01 < dTime) {
+ std::cout << Form("TS %5d Sts Digi %6u %8.0f out of TS up [%5.3f]!!!!!", iEntry, uDigiIdx, dTime,
+ dTime / dTsSizeNs)
+ << std::endl;
+ uDigiIdx++;
+ continue;
+ }
+ if (0 == uDigiIdx % 1000000) {
+ std::cout << Form("TS %5d Sts Digi %8u / %8u", iEntry, uDigiIdx, nDigisTof) << std::endl;
+ }
+ for (uint32_t uHistSz = 0; uHistSz < vdIntegrationNs.size(); ++uHistSz) {
+ uint32_t uCurrentCycle = std::floor(dTime / vdIntegrationNs[uHistSz]);
+ if (vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle) {
+ for (; vuIdxHistoCycleinTS[uHistSz] < uCurrentCycle; ++vuIdxHistoCycleinTS[uHistSz]) {
+ double_t dTsFractional =
+ (vdIntegrationNs[uHistSz] * vuIdxHistoCycleinTS[uHistSz]) / dTsSizeNs + iEntry - dStartTs;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalNew[uHistSz][uBinSz]) {
+ double_t dQFactor = ExtractQFactor(vHistoNew[uHistSz][uBinSz]);
+ vQvalNew[uHistSz][uBinSz]->Fill(dTsFractional, dQFactor);
+ vMeanNew[uHistSz][uBinSz]->Fill(dTsFractional, ExtractMean(vHistoNew[uHistSz][uBinSz]));
+
+ for (uint32_t uBin = 1; uBin <= vHistoNew[uHistSz][uBinSz]->GetNbinsX(); ++uBin) {
+ vBinCountDistributionNew[uBinSz]->Fill(vHistoNew[uHistSz][uBinSz]->GetBinContent(uBin));
+ }
+
+ if (0.0 < dQFactor) {
+ vHistoNew[uHistSz][uBinSz]->Reset();
+ }
+ }
+ }
+ }
+ }
+
+ double_t dTimeInCycle = std::fmod(dTime, vdIntegrationNs[uHistSz]);
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ if (nullptr != vQvalNew[uHistSz][uBinSz]) {
+ vHistoNew[uHistSz][uBinSz]->Fill(dTimeInCycle);
+ }
+ }
+ }
+ uDigiIdx++;
+ }
+ /// FIXME: last "slice" is lost or properly take?
+ std::cout << Form("Donw with TS %5d Digi %8u / %8u", iEntry, uDigiIdx, nDigisTof) << std::endl;
+ }
+
+ uint16_t nPadX = std::ceil(std::sqrt(uNbPlots));
+ uint16_t nPadY = std::ceil(1.0 * uNbPlots / nPadX);
+ std::cout << Form("Nb plots %3d Nb pads X %2d Y %2d", uNbPlots, nPadX, nPadY) << std::endl;
+
+ uint32_t uNbBinsLog = 0;
+ std::vector<double> dBinsLogVector = GenerateLogBinArray(4, 9, 1, uNbBinsLog, 2);
+ double* dBinsLog = dBinsLogVector.data();
+
+ TH1* hMeanVsBinSzOld = new TH1D(
+ Form("MeanVsBinSzOld_%04u", uRunId),
+ Form("Fit of the Mean vs Bin Size, Old, run %04u; Bin Size [ns]; Mean [Digis/bin]", uRunId), uNbBinsLog, dBinsLog);
+ TH1* hMeanVsBinSzNew = new TH1D(
+ Form("MeanVsBinSzNew_%04u", uRunId),
+ Form("Fit of the Mean vs Bin Size, New, run %04u; Bin Size [ns]; Mean [Digis/bin]", uRunId), uNbBinsLog, dBinsLog);
+ TH1* hQfactorVsBinSzOld =
+ new TH1D(Form("hQfactorVsBinSzOld_%04u", uRunId),
+ Form("Fit of Q-Factor vs Bin Size, run %04u; Bin Size [ns]; Q-Factor []", uRunId), uNbBinsLog, dBinsLog);
+ TH1* hQfactorVsBinSzNew =
+ new TH1D(Form("hQfactorVsBinSzNew_%04u", uRunId),
+ Form("Fit of Q-Factor vs Bin Size, run %04u; Bin Size [ns]; Q-Factor []", uRunId), uNbBinsLog, dBinsLog);
+
+ TCanvas* tempOld = new TCanvas("tempOld", "temp Old");
+ TCanvas* tempNew = new TCanvas("tempNew", "temp New");
+ tempOld->Divide(5, 4);
+ tempNew->Divide(5, 4);
+ uint32_t uPadIdx = 1;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ /// Fit Mean and Q-Factor
+ if (0 < vMeanOld[0][uBinSz]->GetEntries()) {
+ tempOld->cd(uPadIdx);
+ auto fCstMeanOld = new TF1("fCstMeanOld", "[0]", 0, numTimeslices);
+ vMeanOld[0][uBinSz]->Fit(fCstMeanOld, "", "", 10, 50);
+ hMeanVsBinSzOld->Fill(vdBinSizesNs[uBinSz], fCstMeanOld->GetParameter(0));
+
+ TCanvas temp("tempcanv", "temp");
+ temp.cd();
+ auto fCstQfactorOld = new TF1("fCstQfactorOld", "[0]", 0, numTimeslices);
+ vQvalOld[0][uBinSz]->Fit(fCstQfactorOld, "", "", 10, 50);
+ hQfactorVsBinSzOld->Fill(vdBinSizesNs[uBinSz], fCstQfactorOld->GetParameter(0));
+ delete fCstMeanOld;
+ delete fCstQfactorOld;
+ }
+
+ if (0 < vMeanNew[0][uBinSz]->GetEntries()) {
+ tempNew->cd(uPadIdx);
+ auto fCstMeanNew = new TF1("fCstMeanNew", "[0]", 0, numTimeslices);
+ vMeanNew[0][uBinSz]->Fit(fCstMeanNew, "", "", 10, 50);
+ hMeanVsBinSzNew->Fill(vdBinSizesNs[uBinSz], fCstMeanNew->GetParameter(0));
+
+ TCanvas temp("tempcanv", "temp");
+ temp.cd();
+ auto fCstQfactorNew = new TF1("fCstQfactorNew", "[0]", 0, numTimeslices);
+ vQvalNew[0][uBinSz]->Fit(fCstQfactorNew, "", "", 10, 50);
+ hQfactorVsBinSzNew->Fill(vdBinSizesNs[uBinSz], fCstQfactorNew->GetParameter(0));
+ delete fCstMeanNew;
+ delete fCstQfactorNew;
+ }
+ uPadIdx++;
+ }
+
+ TCanvas* test = new TCanvas("test", "test");
+ test->Divide(2, 2);
+
+ test->cd(1);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hMeanVsBinSzOld->Draw("hist");
+
+ test->cd(2);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hMeanVsBinSzNew->Draw("hist");
+
+ test->cd(3);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hQfactorVsBinSzOld->Draw("hist");
+
+ test->cd(4);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogx();
+ gPad->SetLogy();
+ hQfactorVsBinSzNew->Draw("hist");
+
+ TCanvas* cBinCntDistOld = new TCanvas("cBinCntDistOld", "BinCntDist Old");
+ TCanvas* cBinCntDistNew = new TCanvas("cBinCntDistNew", "BinCntDist New");
+ cBinCntDistOld->Divide(5, 4);
+ cBinCntDistNew->Divide(5, 4);
+ uPadIdx = 1;
+ for (uint32_t uBinSz = 0; uBinSz < vdBinSizesNs.size(); ++uBinSz) {
+ cBinCntDistOld->cd(uPadIdx);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogy();
+ vBinCountDistributionOld[uBinSz]->Draw("hist");
+
+ cBinCntDistNew->cd(uPadIdx);
+ gPad->SetGridx();
+ gPad->SetGridy();
+ gPad->SetLogy();
+ vBinCountDistributionNew[uBinSz]->Draw("hist");
+
+ uPadIdx++;
+ }
+
+ TFile* outFile =
+ new TFile(Form("data/sts_tof_q_factor_timescale_%04u_%03luTs.root", uRunId, numTimeslices), "RECREATE");
+ outFile->cd();
+ test->Write();
+ hMeanVsBinSzOld->Write();
+ hMeanVsBinSzNew->Write();
+ hQfactorVsBinSzOld->Write();
+ hQfactorVsBinSzNew->Write();
+ gROOT->cd();
+
+ outFile->Close();
+
+ outFile = new TFile(Form("data/sts_tof_bin_counts_vs_bin_size_%04u_%03luTs_%03.0f_%03.0f.root", uRunId, numTimeslices,
+ dStartTs, dStopTs),
+ "RECREATE");
+ outFile->cd();
+ test->Write();
+ cBinCntDistOld->Write();
+ cBinCntDistNew->Write();
+ gROOT->cd();
+
+ outFile->Close();
+
+ std::cout << "Beam particles in spill: " << uTotalCountOld << " (Old) VS " << uTotalCountNew << " (New)" << std::endl;
+ double_t dAvgCountSecOld = (1.0 * uTotalCountOld) / (numTimeslices * 0.128);
+ double_t dAvgCountSecNew = (1.0 * uTotalCountNew) / (numTimeslices * 0.128);
+ std::cout << "Avg Beam particles/s: " << dAvgCountSecOld << " (Old) VS " << dAvgCountSecNew << " (New)" << std::endl;
+}