diff --git a/EventFiltering/macros/checkBCrangesSkimming.C b/EventFiltering/macros/checkBCrangesSkimming.C index 3dea4c56318..bd6d39b41fb 100644 --- a/EventFiltering/macros/checkBCrangesSkimming.C +++ b/EventFiltering/macros/checkBCrangesSkimming.C @@ -10,13 +10,14 @@ // or submit itself to any jurisdiction. // O2 includes -#include -#include -#include #include #include -#include #include +#include +#include +#include +#include +#include #include "CommonDataFormat/InteractionRecord.h" #include "CommonDataFormat/IRFrame.h" @@ -24,9 +25,7 @@ using o2::InteractionRecord; using o2::dataformats::IRFrame; // Set the bit of trigger which need to be checked -const ULong64_t Trigger0BIT = BIT(61); -const ULong64_t Trigger1BIT = 0; -const ULong64_t bcDiffTolerance = 100; +const ULong64_t bcDiffTolerance = 0; const char outputFileName[15] = "output.root"; struct bcTuple { @@ -40,6 +39,7 @@ struct bcTuple { }; struct selectedFrames : public IRFrame { + selectedFrames(ULong64_t bcAO2D, ULong64_t bcEvSel, const IRFrame& frame) : IRFrame(frame), bcAO2D(bcAO2D), bcEvSel(bcEvSel), triMask{0, 0}, selMask{0, 0} {} selectedFrames(ULong64_t bcAO2D, ULong64_t bcEvSel, ULong64_t triMask[2], ULong64_t selMask[2], const IRFrame& frame) : IRFrame(frame), bcAO2D(bcAO2D), bcEvSel(bcEvSel), triMask{triMask[0], triMask[1]}, selMask{selMask[0], selMask[1]} {} ULong64_t triMask[2]{0ull}, selMask[2]{0ull}, bcAO2D, bcEvSel; int numSameTriggerInNearbyBCs = 0; // related to bcDiffTolerance @@ -58,20 +58,32 @@ int DoBCSubraction(ULong64_t bc1, ULong64_t bc2) } } +int DoBCSubraction(selectedFrames bc1, selectedFrames bc2) +{ + if (bc1.getMin() > bc2.getMax()) { + return DoBCSubraction(bc1.getMin().toLong(), bc2.getMax().toLong()); + } else if (bc1.getMax() < bc2.getMin()) { + return DoBCSubraction(bc1.getMax().toLong(), bc2.getMin().toLong()); + } else { + return 0; + } +} + bool isClose(selectedFrames a, selectedFrames b, ULong64_t bcDiffTolerance) { - if (a.getMin() > b.getMax() + bcDiffTolerance || a.getMax() < b.getMin() - bcDiffTolerance) + if (a.getMin() > b.getMax() + bcDiffTolerance || a.getMax() + bcDiffTolerance < b.getMin()) return false; else return true; } -std::vector getSelectedFrames(TFile& file, ULong64_t trigger0Bit, ULong64_t trigger1Bit) +std::vector> getFrames(std::unique_ptr& file, int trgIDStart, int N) { - std::vector selectedFrames; ULong64_t bcAO2D{0ull}, bcEvSel{0ull}, triMask[2]{0ull}, selMask[2]{0ull}; - for (auto key : *file.GetListOfKeys()) { - auto dir = dynamic_cast(file.Get(key->GetName())); + std::vector> frames; + frames.resize(N); + for (auto key : *file->GetListOfKeys()) { + auto dir = dynamic_cast(file->Get(key->GetName())); if (!dir) { continue; } @@ -96,18 +108,36 @@ std::vector getSelectedFrames(TFile& file, ULong64_t trigger0Bit if (!selMask[0] && !selMask[1]) { continue; } - if (selMask[0] & trigger0Bit || selMask[1] & trigger1Bit) { - InteractionRecord irstart, irend; - irstart.setFromLong(std::min(bcAO2D, bcEvSel)); - irend.setFromLong(std::max(bcAO2D, bcEvSel)); - IRFrame frame(irstart, irend); - selectedFrames.push_back({bcAO2D, bcEvSel, triMask, selMask, frame}); + for (int trgID = trgIDStart; trgID < trgIDStart + N; trgID++) { + ULong64_t trigger0Bit = 0, trigger1Bit = 0; + if (trgID < 64) { + trigger0Bit = BIT(trgID); + } else { + trigger1Bit = BIT(trgID - 64); + } + if (selMask[0] & trigger0Bit || selMask[1] & trigger1Bit) { + InteractionRecord irstart, irend; + irstart.setFromLong(std::min(bcAO2D, bcEvSel)); + irend.setFromLong(std::max(bcAO2D, bcEvSel)); + IRFrame frame(irstart, irend); + int index = trgID - trgIDStart; + frames[index].push_back({bcAO2D, bcEvSel, triMask, selMask, frame}); + } } } } - return selectedFrames; + + return frames; +} + +std::vector getSelectedFrames(std::unique_ptr& file, int trgID) +{ + auto frames = getFrames(file, trgID, 1); + return frames[0]; } +// Check how many other triggers are in a compatible BC window with the current one +// Ideally, most of triggers are singles (num = 1), which means none of others is in the same window void checkNearbyBCs(std::vector& frames, ULong64_t bcDiffTolerance) { std::sort(frames.begin(), frames.end(), [](const selectedFrames& a, const selectedFrames& b) { @@ -138,27 +168,163 @@ void checkNearbyBCs(std::vector& frames, ULong64_t bcDiffToleran } } -// Calulate the ratio of duplicate triggers -void checkDuplicateTriggerAndBCs(std::string AnaFileName = "AnalysisResults.root", std::string originalFileName = "bcRanges_fullrun.root", std::string skimmedFileName = "bcRanges_fullrun_skimmed.root") +// Get RunNumber +std::string getRunNumber(std::string fileName) { - - // Get RunNumber std::string runNumber = ""; std::regex re("/5[0-9]*"); std::smatch match; - if (std::regex_search(originalFileName, match, re)) { + if (std::regex_search(fileName, match, re)) { // Remove the leading '/' runNumber = match.str().substr(1); } + return runNumber; +} + +// Detailed checks for specific trigger +void checkBCForSelectedTrg(std::vector& originalFrames, std::vector& skimmedFrames, string runNumber, string triggerLabel) +{ + + TH1D hTriggerCounter("hTriggerCounter", (runNumber + " " + triggerLabel + ";;Total number of trigger").data(), 2, -0.5, 1.5); + hTriggerCounter.GetXaxis()->SetBinLabel(1, "Original"); + hTriggerCounter.GetXaxis()->SetBinLabel(2, "Skimmed"); + TH1D hBCDiffAO2D("hBCDiffAO2D", (runNumber + " " + triggerLabel + ";;#DeltaBC_{AO2D} between paired singles").data(), 201, -100.5, 100.5); + TH1D hBCDiffEvSel("hBCDiffEvSel", (runNumber + " " + triggerLabel + ";;#DeltaBC_{EvSel} between paired singles").data(), 201, -100.5, 100.5); + + TH1D hBCOriginal("hBCOriginal", (runNumber + " " + triggerLabel + " Original;;Trigger counts").data(), 4, -0.5, 3.5); + hBCOriginal.GetXaxis()->SetBinLabel(1, "Total"); + hBCOriginal.GetXaxis()->SetBinLabel(2, "Same AO2D BC"); + hBCOriginal.GetXaxis()->SetBinLabel(3, "Same EvSel BC"); + hBCOriginal.GetXaxis()->SetBinLabel(4, "Same Both BC"); + TH1D hBCSkimmed("hBCSkimmed", (runNumber + " " + triggerLabel + " Skimmed;;Trigger counts").data(), 4, -0.5, 3.5); + hBCSkimmed.GetXaxis()->SetBinLabel(1, "Total"); + hBCSkimmed.GetXaxis()->SetBinLabel(2, "Same AO2D BC"); + hBCSkimmed.GetXaxis()->SetBinLabel(3, "Same EvSel BC"); + hBCSkimmed.GetXaxis()->SetBinLabel(4, "Same Both BC"); - // Checks for BC difference between original and skimming data, and the ratio of triggers which have BCdiff==0 - TH1D hPairedNumCounterTotal("hPairedNumCounterTotal", "hPairedNumCounterTotal", 10, -0.5, 9.5); - TH1D hBCDiffAO2DTotal("hBCDiffAO2DTotal", "hBCDiffAO2DTotal", 201, -100.5, 100.5); - TH1D hBCDiffEvSelTotal("hBCDiffEvSelTotal", "hBCDiffEvSelTotal", 201, -100.5, 100.5); + TH1D hPairedNumCounter("hPairedNumCounter", (runNumber + " " + triggerLabel + ";;Number of matched triggers in skimmed data").data(), 10, -0.5, 9.5); + + checkNearbyBCs(originalFrames, bcDiffTolerance); + checkNearbyBCs(skimmedFrames, bcDiffTolerance); + + std::vector bcSet; + int firstID = 0; + for (auto frame : originalFrames) { + hTriggerCounter.Fill(0); + hBCOriginal.Fill(0); + //------------------------------ Check if there are triggers which have same BC, time-consuming! ------------------------------------------------------- + auto p1 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcAO2D == frame.bcAO2D; }); + if (p1 != bcSet.end()) { + hBCOriginal.Fill(1); + } + auto p2 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcEvSel == frame.bcEvSel; }); + if (p2 != bcSet.end()) { + hBCOriginal.Fill(2); + } + bcTuple currentBC(frame.bcAO2D, frame.bcEvSel); + auto p3 = std::find(bcSet.begin(), bcSet.end(), currentBC); + if (p3 == bcSet.end()) { + bcSet.push_back(currentBC); + } else { + hBCOriginal.Fill(3); + } + //------------------------------------------------------------------------------------- + + if (frame.GetNum() != 1) { + continue; // Only check singles + } + std::vector skimmedbcs; + int n = 0; + bool isFirst = true; + for (int i = firstID; i < skimmedFrames.size(); i++) { + auto& skimmedFrame = skimmedFrames[i]; + if (skimmedFrame.getMin() > frame.getMax()) { + break; + } + if (skimmedFrame.GetNum() != 1) { + continue; // Only check singles + } + if (isClose(frame, skimmedFrame, bcDiffTolerance)) { + bool found = frame.selMask[0] & skimmedFrame.selMask[0] || frame.selMask[1] & skimmedFrame.selMask[1]; + if (found) { + skimmedbcs.push_back({skimmedFrame.bcAO2D, skimmedFrame.bcEvSel}); + n++; + } + } else { + if (isFirst) { + firstID = i; + } + } + } + if (n == 0) { + // std::cout << "Trigger not found!!!" << std::endl; + } else if (n == 1) { + hBCDiffAO2D.Fill(DoBCSubraction(frame.bcAO2D, skimmedbcs[0].bcAO2D)); + hBCDiffEvSel.Fill(DoBCSubraction(frame.bcEvSel, skimmedbcs[0].bcEvSel)); + } + hPairedNumCounter.Fill(n); + } + + //------------------------------ Check if there are triggers which have same BC, time-consuming! ------------------------------------------------------- + bcSet.clear(); + for (auto& skimmedFrame : skimmedFrames) { + hTriggerCounter.Fill(1); + hBCSkimmed.Fill(0); + auto p1 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcAO2D == skimmedFrame.bcAO2D; }); + if (p1 != bcSet.end()) { + hBCSkimmed.Fill(1); + } + auto p2 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcEvSel == skimmedFrame.bcEvSel; }); + if (p2 != bcSet.end()) { + hBCSkimmed.Fill(2); + } + bcTuple currentBC(skimmedFrame.bcAO2D, skimmedFrame.bcEvSel); + auto p3 = std::find(bcSet.begin(), bcSet.end(), currentBC); + if (p3 == bcSet.end()) { + bcSet.push_back(currentBC); + } else { + hBCSkimmed.Fill(3); + } + } + //------------------------------------------------------------------------------------- + + TFile fout(outputFileName, "UPDATE"); + fout.cd(); + TDirectory* dir1 = fout.GetDirectory(runNumber.data()); + if (!dir1) { + dir1 = fout.mkdir(runNumber.data()); + } + dir1->cd(); + TDirectory* dir2 = dir1->GetDirectory(triggerLabel.data()); + if (!dir2) { + dir2 = dir1->mkdir(triggerLabel.data()); + } + dir2->cd(); + + hTriggerCounter.Write(); + hBCOriginal.Write(); + hBCSkimmed.Write(); + hBCDiffAO2D.Write(); + hBCDiffEvSel.Write(); + hPairedNumCounter.Write(); + fout.Close(); +} + +// Detailed checks for specific trigger +void checkBCForSelectedTrg(std::string AnaFileName = "AnalysisResults.root", std::string originalFileName = "bcRanges_fullrun.root", std::string skimmedFileName = "bcRanges_fullrun_skimmed.root", int triggerID = 1, bool useAlien = true) +{ + + string runNumber = getRunNumber(originalFileName); + if (useAlien) { + TGrid::Connect("alien://"); + AnaFileName = "alien://" + AnaFileName; + originalFileName = "alien://" + originalFileName; + skimmedFileName = "alien://" + skimmedFileName; + } // Readin labels - TFile AnaFile(AnaFileName.c_str(), "READ"); - TH1* hist0 = dynamic_cast(AnaFile.Get("central-event-filter-task/scalers/mFiltered;1")); + std::unique_ptr AnaFile{TFile::Open(AnaFileName.c_str(), "READ")}; + TH1* hist0 = dynamic_cast(AnaFile->Get("central-event-filter-task/scalers/mFiltered;1")); std::vector labels; std::vector binNum; for (int i = 1; i <= hist0->GetNbinsX(); i++) { @@ -168,29 +334,69 @@ void checkDuplicateTriggerAndBCs(std::string AnaFileName = "AnalysisResults.root binNum.push_back(i); } } - AnaFile.Close(); + AnaFile->Close(); + std::string triggerLabel = labels[triggerID]; - TFile originalFile(originalFileName.c_str(), "READ"); - TFile skimmedFile(skimmedFileName.c_str(), "READ"); - std::vector sel_labels; - std::vector numOriginal, numSkimmed, numOriginalSingle, numSkimmedSingle, numOriginalDouble, numSkimmedDouble, numOriginalMultiple, numSkimmedMultiple; - std::vector numpair, numpairedBCAO2D, numpairedBCEvSel, maxDeltaBCAO2D, maxDeltaBCEvSel; - for (int i = 0; i < labels.size(); i++) { - // std::cout << "i:" << i << std::endl; - ULong64_t trigger0Bit = 0, trigger1Bit = 0; - int triggerBit = binNum[i] - 2; - if (triggerBit < 64) { - trigger0Bit = BIT(triggerBit); - } else { - trigger1Bit = BIT(triggerBit - 64); + std::unique_ptr originalFile{TFile::Open(originalFileName.c_str(), "READ")}; + std::unique_ptr skimmedFile{TFile::Open(skimmedFileName.c_str(), "READ")}; + auto originalFrames = getSelectedFrames(originalFile, triggerID); + auto skimmedFrames = getSelectedFrames(skimmedFile, triggerID); + originalFile->Close(); + skimmedFile->Close(); + + checkBCForSelectedTrg(originalFrames, skimmedFrames, runNumber, triggerLabel); +} + +// Calulate the ratio of duplicate triggers +void checkBCrangesSkimming(std::string AnaFileName = "AnalysisResults.root", std::string originalFileName = "bcRanges_fullrun.root", std::string skimmedFileName = "bcRanges_fullrun_skimmed.root", bool useAlien = true) +{ + + string runNumber = getRunNumber(originalFileName); + if (useAlien) { + TGrid::Connect("alien://"); + AnaFileName = "alien://" + AnaFileName; + originalFileName = "alien://" + originalFileName; + skimmedFileName = "alien://" + skimmedFileName; + } + + // Readin labels + std::unique_ptr AnaFile{TFile::Open(AnaFileName.c_str(), "READ")}; + TH1* hist0 = dynamic_cast(AnaFile->Get("central-event-filter-task/scalers/mFiltered;1")); + std::vector labels; + std::vector binNum; + for (int i = 1; i <= hist0->GetNbinsX(); i++) { + std::string label = hist0->GetXaxis()->GetBinLabel(i); + if (label != "Total number of events" && label != "Filtered events") { + labels.push_back(label); + binNum.push_back(i); + // std::cout << i - 2 << ": " << label << std::endl; } + } + AnaFile->Close(); + + // Due to potential selection on triggers, histograms should be created later + // for example: skip triggers which have no enrties + std::vector sel_labels; + std::vector numOriginal, numSkimmed, numOriginalSingle, numSkimmedSingle, numOriginalDouble, numSkimmedDouble, numOriginalMultiple, numSkimmedMultiple, numCloseSkimmed; + std::vector numpair, numpairedBCAO2D, numpairedBCEvSel; + std::vector avgDeltaBCAO2D, avgDeltaBCEvSel, avgDeltaBC, rmsDeltaBCAO2D, rmsDeltaBCEvSel, rmsDeltaBC; + std::vector avgNumPairedTrigger, rmsNumPairedTrigger; + + std::unique_ptr originalFile{TFile::Open(originalFileName.c_str(), "READ")}; + std::unique_ptr skimmedFile{TFile::Open(skimmedFileName.c_str(), "READ")}; + std::vector> originalAllFrames = getFrames(originalFile, 0, labels.size()); + std::vector> skimmedAllFrames = getFrames(skimmedFile, 0, labels.size()); + for (int trgID = 0; trgID < labels.size(); trgID++) { // Caculate singles, doubles, and multiples + int noriginal{0}, nskimmed{0}, noriginalsingle{0}, nskimmedsingle{0}, noriginaldouble{0}, nskimmeddouble{0}, noriginalmultiple{0}, nskimmedmultiple{0}; + // Caculate mean and rms of diff BC + TH1D hDiffBCAO2DCount("hDiffBCAO2DCount", "hDiffBCAO2DCount", 21, -10.5, 10.5); + TH1D hDiffBCEvSelCount("hDiffBCEvSelCount", "hDiffBCEvSelCount", 21, -10.5, 10.5); + TH1D hDiffBCCount("hDiffBCCount", "hDiffBCCount", 21, -10.5, 10.5); + TH1D hNumPairedTriggerCount("hNumPairedTriggerCount", "hNumPairedTriggerCount", 10, -0.5, 9.5); // For Original dataset - std::vector bcSet; - std::vector bcFullSet; - int noriginal{0}, nskimmed{0}, noriginalsingle{0}, nskimmedsingle{0}, noriginaldouble{0}, nskimmeddouble{0}, noriginalmultiple{0}, nskimmedmultiple{0}, maxdiffBCAO2D{0}, maxdiffBCEvSel{0}; - auto originalFrames = getSelectedFrames(originalFile, trigger0Bit, trigger1Bit); - checkNearbyBCs(originalFrames, bcDiffTolerance); + auto& originalFrames = originalAllFrames[trgID]; + checkNearbyBCs(originalFrames, bcDiffTolerance); // include sorting noriginal = originalFrames.size(); for (auto originalFrame : originalFrames) { if (originalFrame.GetNum() == 0) { @@ -204,8 +410,8 @@ void checkDuplicateTriggerAndBCs(std::string AnaFileName = "AnalysisResults.root } } // For skimmed dataset - auto skimmedFrames = getSelectedFrames(skimmedFile, trigger0Bit, trigger1Bit); - checkNearbyBCs(skimmedFrames, bcDiffTolerance); + auto& skimmedFrames = skimmedAllFrames[trgID]; + checkNearbyBCs(skimmedFrames, bcDiffTolerance); // include sorting nskimmed = skimmedFrames.size(); for (auto& skimmedFrame : skimmedFrames) { if (skimmedFrame.GetNum() == 0) { @@ -219,344 +425,233 @@ void checkDuplicateTriggerAndBCs(std::string AnaFileName = "AnalysisResults.root } } - sel_labels.push_back(labels[i]); - numOriginal.push_back(noriginal); - numOriginalSingle.push_back(noriginalsingle); - numOriginalDouble.push_back(noriginaldouble); - numOriginalMultiple.push_back(noriginalmultiple); - numSkimmed.push_back(nskimmed); - numSkimmedSingle.push_back(nskimmedsingle); - numSkimmedDouble.push_back(nskimmeddouble); - numSkimmedMultiple.push_back(nskimmedmultiple); - // Check BC differences - int npair{0}, npairedBCAO2D{0}, npairedBCEvSel{0}, maxdeltaBCAO2D{0}, maxdeltaBCEvSel{0}; + int npair{0}, npairedBCAO2D{0}, npairedBCEvSel{0}, ncloseskimmed{0}, maxdeltaBCAO2D{0}, maxdeltaBCEvSel{0}; int firstID = 0; for (auto frame : originalFrames) { - if (frame.selMask[0] & trigger0Bit || frame.selMask[1] & trigger1Bit) { - // std::cout << "------------------------------------------------" << std::endl; - if (frame.GetNum() != 1) { + if (frame.GetNum() != 1) { + continue; // Only check singles + } + std::vector skimmedbcs; + int n = 0; + bool isFirst = true; + for (int i = firstID; i < skimmedFrames.size(); i++) { + auto& skimmedFrame = skimmedFrames[i]; + if (skimmedFrame.getMin() > frame.getMax()) { + break; + } + if (skimmedFrame.GetNum() != 1) { continue; // Only check singles } - std::vector skimmedbcs; - int n = 0; - bool isFirst = true; - for (int i = firstID; i < skimmedFrames.size(); i++) { - auto& skimmedFrame = skimmedFrames[i]; - if (skimmedFrame.getMin() > frame.getMax()) { - break; - } - if (skimmedFrame.GetNum() != 1) { - continue; // Only check singles + if (isClose(frame, skimmedFrame, bcDiffTolerance)) { + isFirst = false; + bool found = frame.selMask[0] & skimmedFrame.selMask[0] || frame.selMask[1] & skimmedFrame.selMask[1]; + if (found) { + InteractionRecord irstart, irend; + irstart.setFromLong(std::min(skimmedFrame.bcAO2D, skimmedFrame.bcEvSel)); + irend.setFromLong(std::max(skimmedFrame.bcAO2D, skimmedFrame.bcEvSel)); + IRFrame frame(irstart, irend); + skimmedbcs.push_back({skimmedFrame.bcAO2D, skimmedFrame.bcEvSel, frame}); + n++; } - if (isClose(frame, skimmedFrame, bcDiffTolerance)) { - isFirst = false; - bool found = frame.selMask[0] & skimmedFrame.selMask[0] || frame.selMask[1] & skimmedFrame.selMask[1]; - // found = found && (frame.bcAO2D == skimmedFrame.bcAO2D || frame.bcEvSel == skimmedFrame.bcEvSel); - if (found) { - skimmedbcs.push_back({skimmedFrame.bcAO2D, skimmedFrame.bcEvSel}); - n++; - } - } else { - if (isFirst) { - firstID = i; - } + } else { + if (isFirst) { + firstID = i; } } - if (n == 1) { - npair++; - int bcdiffAO2D = DoBCSubraction(frame.bcAO2D, skimmedbcs[0].bcAO2D); - int bcdiffEvSel = DoBCSubraction(frame.bcEvSel, skimmedbcs[0].bcEvSel); - hBCDiffAO2DTotal.Fill(bcdiffAO2D); - hBCDiffEvSelTotal.Fill(bcdiffEvSel); - maxdiffBCAO2D = std::max(std::abs(maxdiffBCAO2D), bcdiffAO2D); - maxdiffBCEvSel = std::max(std::abs(maxdiffBCEvSel), bcdiffEvSel); - if (frame.bcAO2D == skimmedbcs[0].bcAO2D) { - npairedBCAO2D++; - } - if (frame.bcEvSel == skimmedbcs[0].bcEvSel) { - npairedBCEvSel++; - } + } + if (n == 1) { + npair++; + int bcdiffAO2D = DoBCSubraction(frame.bcAO2D, skimmedbcs[0].bcAO2D); + int bcdiffEvSel = DoBCSubraction(frame.bcEvSel, skimmedbcs[0].bcEvSel); + hDiffBCAO2DCount.Fill(std::abs(bcdiffAO2D)); + hDiffBCEvSelCount.Fill(std::abs(bcdiffEvSel)); + hDiffBCCount.Fill(std::abs(DoBCSubraction(frame, skimmedbcs[0]))); + if (frame.bcAO2D == skimmedbcs[0].bcAO2D) { + npairedBCAO2D++; + } + if (frame.bcEvSel == skimmedbcs[0].bcEvSel) { + npairedBCEvSel++; } - hPairedNumCounterTotal.Fill(n); } + ncloseskimmed += n; + hNumPairedTriggerCount.Fill(n); + } + + // if (static_cast(ncloseskimmed) / noriginal > 0.95 || noriginal == 0) + if (noriginal == 0) { + // continue; } + sel_labels.push_back(labels[trgID]); + numOriginal.push_back(noriginal); + numOriginalSingle.push_back(noriginalsingle); + numOriginalDouble.push_back(noriginaldouble); + numOriginalMultiple.push_back(noriginalmultiple); + numSkimmed.push_back(nskimmed); + numSkimmedSingle.push_back(nskimmedsingle); + numSkimmedDouble.push_back(nskimmeddouble); + numSkimmedMultiple.push_back(nskimmedmultiple); + numpair.push_back(npair); numpairedBCAO2D.push_back(npairedBCAO2D); numpairedBCEvSel.push_back(npairedBCEvSel); - maxDeltaBCAO2D.push_back(maxdiffBCAO2D); - maxDeltaBCEvSel.push_back(maxdiffBCEvSel); + numCloseSkimmed.push_back(ncloseskimmed); + avgDeltaBCAO2D.push_back(hDiffBCAO2DCount.GetMean()); + avgDeltaBCEvSel.push_back(hDiffBCEvSelCount.GetMean()); + avgDeltaBC.push_back(hDiffBCCount.GetMean()); + rmsDeltaBCAO2D.push_back(hDiffBCAO2DCount.GetRMS()); + rmsDeltaBCEvSel.push_back(hDiffBCEvSelCount.GetRMS()); + rmsDeltaBC.push_back(hDiffBCCount.GetRMS()); + avgNumPairedTrigger.push_back(hNumPairedTriggerCount.GetMean()); + rmsNumPairedTrigger.push_back(hNumPairedTriggerCount.GetRMS()); } - originalFile.Close(); - skimmedFile.Close(); - - TH1D hOriginalTotal("hOriginalTotal", (runNumber + " AO2D Original;;Number of events").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hOriginalSingles("hOriginalSingles", (runNumber + " Original;;Number of Singles").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hOriginalSinglesRatio("hOriginalSinglesRatio", (runNumber + " Original;;Singles / Total").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hOriginalDoubles("hOriginalDoubles", (runNumber + " Original;;Number of Doubles").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hOriginalDoublesRatio("hOriginalDoublesRatio", (runNumber + " Original;;Doubles / Total").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hOriginalMultiples("hOriginalMultiples", (runNumber + " Original;;Number of Multiples").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hOriginalMultiplesRatio("hOriginalMultiplesRatio", (runNumber + " Original;;Multiples / Total").data(), sel_labels.size(), 0, sel_labels.size()); - - TH1D hSkimmedTotal("hSkimmedTotal", (runNumber + " AO2D Skimmed;;Number of events").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hSkimmedSingles("hSkimmedSingles", (runNumber + " Skimmed;;Number of Singles").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hSkimmedSinglesRatio("hSkimmedSinglesRatio", (runNumber + " Skimmed;;Singles / Total").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hSkimmedDoubles("hSkimmedDoubles", (runNumber + " Skimmed;;Number of Doubles").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hSkimmedDoublesRatio("hSkimmedDoublesRatio", (runNumber + " Skimmed;;Doubles / Total").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hSkimmedMultiples("hSkimmedMultiples", (runNumber + " Skimmed;;Number of Multiples").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hSkimmedMultiplesRatio("hSkimmedMultiplesRatio", (runNumber + " Skimmed;;Multiples / Total").data(), sel_labels.size(), 0, sel_labels.size()); - - TH1D hPairedBCAO2DRatio("hPairedBCAO2DRatio", (runNumber + " One-to-One Pairs;; Pairs with same BCAO2D / Total").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hPairedBCEvSelRatio("hPairedBCEvSelRatio", (runNumber + " One-to-One Pairs;; Pairs with same BCEvSel / Total").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hMaxDiffBCAO2D("hMaxDiffBCAO2D", (runNumber + " One-to-One Pairs;;|#DeltaBCAO2D|_{max}").data(), sel_labels.size(), 0, sel_labels.size()); - TH1D hMaxDiffBCEvSel("hMaxDiffBCEvSel", (runNumber + " One-to-One Pairs;;|#DeltaBCEvSel|_{max}").data(), sel_labels.size(), 0, sel_labels.size()); + originalFile->Close(); + skimmedFile->Close(); + + TH1D hOriginalTotal("hOriginalTotal", (runNumber + " Original;;Number of events").data(), sel_labels.size(), 0, sel_labels.size()); + TH1D hOriginalSingles("hOriginalSingles", (runNumber + " Original;;Number of singles").data(), sel_labels.size(), 0, sel_labels.size()); + TH1D hOriginalDoubles("hOriginalDoubles", (runNumber + " Original;;Number of doubles").data(), sel_labels.size(), 0, sel_labels.size()); + TH1D hOriginalMultiples("hOriginalMultiples", (runNumber + " Original;;Number of multiples").data(), sel_labels.size(), 0, sel_labels.size()); + + TH1D hSkimmedTotal("hSkimmedTotal", (runNumber + " Skimmed;;Number of events").data(), sel_labels.size(), 0, sel_labels.size()); + TH1D hSkimmedSingles("hSkimmedSingles", (runNumber + " Skimmed;;Number of singles").data(), sel_labels.size(), 0, sel_labels.size()); + TH1D hSkimmedDoubles("hSkimmedDoubles", (runNumber + " Skimmed;;Number of doubles").data(), sel_labels.size(), 0, sel_labels.size()); + TH1D hSkimmedMultiples("hSkimmedMultiples", (runNumber + " Skimmed;;Number of multiples").data(), sel_labels.size(), 0, sel_labels.size()); + + TH1D hTriggerPairsRatio("hTriggerPairsRatio", (runNumber + " Skimmed Efficiency;; Matched skimmed triggers / Original singles").data(), sel_labels.size(), 0, sel_labels.size()); // the ratio of triggers in skimmed dataset whose BC is compatible with original triggers to the number of original triggers, might be duplicate since we check it based on every trigger in unskimmed data + TH1D hTriggerSinglePairsRatio("hTriggerSinglePairsRatio", (runNumber + " Skimmed Efficiency;; One-to-one pairs / Original singles").data(), sel_labels.size(), 0, sel_labels.size()); // the ratio of 1-1 paired triggers to the number of original triggers + TH1D hPairsSameBCAO2DRatio("hPairsSameBCAO2DRatio", (runNumber + " One-to-one pairs;; Pairs with same BC_{AO2D} / Total").data(), sel_labels.size(), 0, sel_labels.size()); // In 1-1 pairs, the ratio of pairs who have same BCAO2D + TH1D hPairsSameBCEvSelRatio("hPairsSameBCEvSelRatio", (runNumber + " One-to-one pairs;; Pairs with same BC_{EvSel} / Total").data(), sel_labels.size(), 0, sel_labels.size()); // In 1-1 pairs, the ratio of pairs who have same BCEvSel + TH1D hDiffBCAO2D("hDiffBCAO2D", (runNumber + " One-to-one pairs;;#DeltaBC_{AO2D}").data(), sel_labels.size(), 0, sel_labels.size()); // difference in BCAO2D of 1-1 pairs + TH1D hDiffBCEvSel("hDiffBCEvSel", (runNumber + " One-to-one pairs;;#DeltaBC_{EvSel}").data(), sel_labels.size(), 0, sel_labels.size()); // difference in BCEvSel of 1-1 pairs + TH1D hDiffBC("hDiffBC", (runNumber + " One-to-one pairs;;#DeltaBC").data(), sel_labels.size(), 0, sel_labels.size()); // difference between the BC tuple, expected to be 0 if bcDiffTolerance = 0 + TH1D hNumPairsInSkimmed("hNumPairsInSkimmed", (runNumber + " number of matched triggers in skimmed data;;Matched trigger count").data(), sel_labels.size(), 0, sel_labels.size()); // number of triggers in skimmed data which are compatible in the BC ranges of singles in original selection for (int i = 0; i < sel_labels.size(); i++) { + // Original data hOriginalTotal.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); hOriginalSingles.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hOriginalDoubles.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); hOriginalMultiples.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hOriginalSinglesRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hOriginalDoublesRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hOriginalMultiplesRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); hOriginalTotal.SetBinContent(i + 1, numOriginal[i]); + hOriginalTotal.SetBinError(i + 1, std::sqrt(numOriginal[i])); hOriginalSingles.SetBinContent(i + 1, numOriginalSingle[i]); + hOriginalSingles.SetBinError(i + 1, std::sqrt(numOriginalSingle[i])); hOriginalDoubles.SetBinContent(i + 1, numOriginalDouble[i]); + hOriginalDoubles.SetBinError(i + 1, std::sqrt(numOriginalDouble[i])); hOriginalMultiples.SetBinContent(i + 1, numOriginalMultiple[i]); - if (numOriginal[i] > 0) { - hOriginalSinglesRatio.SetBinContent(i + 1, static_cast(numOriginalSingle[i]) / numOriginal[i]); - hOriginalDoublesRatio.SetBinContent(i + 1, static_cast(numOriginalSingle[i]) / numOriginal[i]); - hOriginalMultiplesRatio.SetBinContent(i + 1, static_cast(numOriginalMultiple[i]) / numOriginal[i]); - } else { - hOriginalSinglesRatio.SetBinContent(i + 1, 0); - hOriginalDoublesRatio.SetBinContent(i + 1, 0); - hOriginalMultiplesRatio.SetBinContent(i + 1, 0); - } + hOriginalMultiples.SetBinError(i + 1, std::sqrt(numOriginalMultiple[i])); + // Skimmed data hSkimmedTotal.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); hSkimmedSingles.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hSkimmedDoubles.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); hSkimmedMultiples.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hSkimmedSinglesRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hSkimmedDoublesRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hSkimmedMultiplesRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); hSkimmedTotal.SetBinContent(i + 1, numSkimmed[i]); + hSkimmedTotal.SetBinError(i + 1, std::sqrt(numSkimmed[i])); hSkimmedSingles.SetBinContent(i + 1, numSkimmedSingle[i]); + hSkimmedSingles.SetBinError(i + 1, std::sqrt(numSkimmedSingle[i])); hSkimmedDoubles.SetBinContent(i + 1, numSkimmedDouble[i]); + hSkimmedDoubles.SetBinError(i + 1, std::sqrt(numSkimmedDouble[i])); hSkimmedMultiples.SetBinContent(i + 1, numSkimmedMultiple[i]); - if (numSkimmed[i] > 0) { - hSkimmedSinglesRatio.SetBinContent(i + 1, static_cast(numSkimmedSingle[i]) / numSkimmed[i]); - hSkimmedDoublesRatio.SetBinContent(i + 1, static_cast(numSkimmedDouble[i]) / numSkimmed[i]); - hSkimmedMultiplesRatio.SetBinContent(i + 1, static_cast(numSkimmedMultiple[i]) / numSkimmed[i]); - } else { - hSkimmedSinglesRatio.SetBinContent(i + 1, 0); - hSkimmedDoublesRatio.SetBinContent(i + 1, 0); - hSkimmedMultiplesRatio.SetBinContent(i + 1, 0); - } + hSkimmedMultiples.SetBinError(i + 1, std::sqrt(numSkimmedMultiple[i])); + + // Pairs QA + hTriggerPairsRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hTriggerSinglePairsRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hPairsSameBCAO2DRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hPairsSameBCEvSelRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hDiffBCAO2D.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hDiffBCEvSel.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hDiffBC.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); + hNumPairsInSkimmed.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hPairedBCAO2DRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hPairedBCEvSelRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); if (numpair[i] > 0) { - hPairedBCAO2DRatio.SetBinContent(i + 1, static_cast(numpairedBCAO2D[i]) / numpair[i]); - hPairedBCEvSelRatio.SetBinContent(i + 1, static_cast(numpairedBCEvSel[i]) / numpair[i]); + hPairsSameBCAO2DRatio.SetBinContent(i + 1, static_cast(numpairedBCAO2D[i]) / numpair[i]); + hPairsSameBCEvSelRatio.SetBinContent(i + 1, static_cast(numpairedBCEvSel[i]) / numpair[i]); } - - hMaxDiffBCAO2D.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hMaxDiffBCEvSel.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str()); - hMaxDiffBCAO2D.SetBinContent(i + 1, maxDeltaBCAO2D[i]); - hMaxDiffBCEvSel.SetBinContent(i + 1, maxDeltaBCEvSel[i]); + hTriggerPairsRatio.SetBinContent(i + 1, numCloseSkimmed[i]); + hTriggerPairsRatio.SetBinError(i + 1, std::sqrt(numCloseSkimmed[i])); + hTriggerSinglePairsRatio.SetBinContent(i + 1, numpair[i]); + hTriggerSinglePairsRatio.SetBinError(i + 1, std::sqrt(numpair[i])); + hDiffBCAO2D.SetBinContent(i + 1, avgDeltaBCAO2D[i]); + hDiffBCAO2D.SetBinError(i + 1, rmsDeltaBCAO2D[i]); + hDiffBCEvSel.SetBinContent(i + 1, avgDeltaBCEvSel[i]); + hDiffBCEvSel.SetBinError(i + 1, rmsDeltaBCEvSel[i]); + hDiffBC.SetBinContent(i + 1, avgDeltaBC[i]); + hDiffBC.SetBinError(i + 1, rmsDeltaBC[i]); + hNumPairsInSkimmed.SetBinContent(i + 1, avgNumPairedTrigger[i]); + hNumPairsInSkimmed.SetBinError(i + 1, rmsNumPairedTrigger[i]); } + TH1D* hTriggerEff; // Ratio of the total number of triggers in skimmed data to that in original data + TH1D *hOriginalSinglesRatio, *hOriginalDoublesRatio, *hOriginalMultiplesRatio; + TH1D *hSkimmedSinglesRatio, *hSkimmedDoublesRatio, *hSkimmedMultiplesRatio; + + hTriggerEff = (TH1D*)hSkimmedTotal.Clone("hTriggerEff"); + hTriggerEff->SetTitle((runNumber + " skimmed efficiency;; Skimmed / Original").data()); + hTriggerEff->Divide(&hOriginalTotal); + hTriggerPairsRatio.Divide(&hOriginalSingles); + hTriggerSinglePairsRatio.Divide(&hOriginalSingles); + hOriginalSinglesRatio = (TH1D*)hOriginalSingles.Clone("hOriginalSinglesRatio"); + hOriginalSinglesRatio->SetTitle((runNumber + " Original;;Singles / Total").data()); + hOriginalSinglesRatio->Divide(&hOriginalTotal); + hOriginalDoublesRatio = (TH1D*)hOriginalDoubles.Clone("hOriginalDoublesRatio"); + hOriginalDoublesRatio->SetTitle((runNumber + " Original;;Doubles / Total").data()); + hOriginalDoublesRatio->Divide(&hOriginalTotal); + hOriginalMultiplesRatio = (TH1D*)hOriginalMultiples.Clone("hOriginalMultiplesRatio"); + hOriginalMultiplesRatio->SetTitle((runNumber + " Original;;Multiples / Total").data()); + hOriginalMultiplesRatio->Divide(&hOriginalTotal); + + hSkimmedSinglesRatio = (TH1D*)hSkimmedSingles.Clone("hSkimmedSinglesRatio"); + hSkimmedSinglesRatio->SetTitle((runNumber + " Skimmed;;Singles / Total").data()); + hSkimmedSinglesRatio->Divide(&hSkimmedTotal); + hSkimmedDoublesRatio = (TH1D*)hSkimmedDoubles.Clone("hSkimmedDoublesRatio"); + hSkimmedDoublesRatio->SetTitle((runNumber + " Skimmed;;Doubles / Total").data()); + hSkimmedDoublesRatio->Divide(&hSkimmedTotal); + hSkimmedMultiplesRatio = (TH1D*)hSkimmedMultiples.Clone("hSkimmedMultiplesRatio"); + hSkimmedMultiplesRatio->SetTitle((runNumber + " Skimmed;;Multiples / Total").data()); + hSkimmedMultiplesRatio->Divide(&hSkimmedTotal); + TFile fout(outputFileName, "UPDATE"); fout.cd(); + TDirectory* dir = fout.mkdir(runNumber.data()); + dir->cd(); + hTriggerEff->Write(); + hTriggerPairsRatio.Write(); + hTriggerSinglePairsRatio.Write(); + hDiffBCAO2D.Write(); + hDiffBCEvSel.Write(); + hNumPairsInSkimmed.Write(); + if (bcDiffTolerance > 0) { + hDiffBC.Write(); + } + TDirectory* dirextra = dir->mkdir("ExtraQA"); + dirextra->cd(); hOriginalTotal.Write(); hOriginalSingles.Write(); hOriginalDoubles.Write(); hOriginalMultiples.Write(); - hOriginalSinglesRatio.Write(); - hOriginalDoublesRatio.Write(); - hOriginalMultiplesRatio.Write(); + hOriginalSinglesRatio->Write(); + hOriginalDoublesRatio->Write(); + hOriginalMultiplesRatio->Write(); hSkimmedTotal.Write(); hSkimmedSingles.Write(); hSkimmedDoubles.Write(); hSkimmedMultiples.Write(); - hSkimmedSinglesRatio.Write(); - hSkimmedDoublesRatio.Write(); - hSkimmedMultiplesRatio.Write(); - hPairedNumCounterTotal.Write(); - hBCDiffAO2DTotal.Write(); - hBCDiffEvSelTotal.Write(); - hPairedBCAO2DRatio.Write(); - hPairedBCEvSelRatio.Write(); - hMaxDiffBCAO2D.Write(); - hMaxDiffBCEvSel.Write(); + hSkimmedSinglesRatio->Write(); + hSkimmedDoublesRatio->Write(); + hSkimmedMultiplesRatio->Write(); + hPairsSameBCAO2DRatio.Write(); + hPairsSameBCEvSelRatio.Write(); fout.Close(); -} - -// Detailed checks for specific trigger -void checkBCrangesSkimming(std::string originalFileName = "bcRanges_fullrun.root", std::string skimmedFileName = "bcRanges_fullrun_skimmed.root") -{ - //---------------------------------For specific trigger---------------------------------- - TH1D hTriggerCounter("hTriggerCounter", "hTriggerCounter", 3, 0.5, 3.5); - hTriggerCounter.GetXaxis()->SetBinLabel(1, "Original"); - hTriggerCounter.GetXaxis()->SetBinLabel(2, "Skimmed"); - TH1D hNumCounter("hNumCounter", "hNumCounter", 10, -0.5, 9.5); - TH1D hSinglePairCheck("hSinglePairCheck", "hSinglePairCheck", 4, 0.5, 4.5); - hSinglePairCheck.GetXaxis()->SetBinLabel(1, "Total"); - hSinglePairCheck.GetXaxis()->SetBinLabel(2, "Same AO2D BC"); - hSinglePairCheck.GetXaxis()->SetBinLabel(3, "Same EvSel BC"); - hSinglePairCheck.GetXaxis()->SetBinLabel(4, "Same Both BC"); - TH1D hMultiPairCheck("hMultiPairCheck", "hMultiPairCheck", 5, 0.5, 5.5); - hMultiPairCheck.GetXaxis()->SetBinLabel(1, "Total"); - hMultiPairCheck.GetXaxis()->SetBinLabel(2, "Total Pair"); - hMultiPairCheck.GetXaxis()->SetBinLabel(3, "Same AO2D BC"); - hMultiPairCheck.GetXaxis()->SetBinLabel(4, "Same EvSel BC"); - hMultiPairCheck.GetXaxis()->SetBinLabel(5, "Same Both BC"); - TH1D hBCDiffAO2D("hBCDiffAO2D", "hBCDiffAO2D", 2001, -1000.5, 1000.5); - TH1D hBCDiffEvSel("hBCDiffEvSel", "hBCDiffEvSel", 2001, -1000.5, 1000.5); - - TH1D hBCOriginal("hBCOriginal", "hBCOriginal", 4, 0.5, 4.5); - hBCOriginal.GetXaxis()->SetBinLabel(1, "Total"); - hBCOriginal.GetXaxis()->SetBinLabel(2, "Same AO2D BC"); - hBCOriginal.GetXaxis()->SetBinLabel(3, "Same EvSel BC"); - hBCOriginal.GetXaxis()->SetBinLabel(4, "Same Both BC"); - TH1D hBCSkimmed("hBCSkimmed", "hBCSkimmed", 4, 0.5, 4.5); - hBCSkimmed.GetXaxis()->SetBinLabel(1, "Total"); - hBCSkimmed.GetXaxis()->SetBinLabel(2, "Same AO2D BC"); - hBCSkimmed.GetXaxis()->SetBinLabel(3, "Same EvSel BC"); - hBCSkimmed.GetXaxis()->SetBinLabel(4, "Same Both BC"); - - auto t1 = std::chrono::steady_clock::now(); - TFile originalFile(originalFileName.c_str(), "READ"); - TFile skimmedFile(skimmedFileName.c_str(), "READ"); - auto originalFrames = getSelectedFrames(originalFile, Trigger0BIT, Trigger1BIT); - auto skimmedFrames = getSelectedFrames(skimmedFile, Trigger0BIT, Trigger1BIT); - originalFile.Close(); - skimmedFile.Close(); - auto t2 = std::chrono::steady_clock::now(); - int d1 = std::chrono::duration_cast(t2 - t1).count(); - std::cout << "Readin Time: " << d1 << std::endl; - - auto t3 = std::chrono::steady_clock::now(); - checkNearbyBCs(originalFrames, bcDiffTolerance); - checkNearbyBCs(skimmedFrames, bcDiffTolerance); - auto t4 = std::chrono::steady_clock::now(); - int d2 = std::chrono::duration_cast(t4 - t3).count(); - std::cout << "Sort Time: " << d2 << std::endl; - auto t5 = std::chrono::steady_clock::now(); - std::vector bcSet; - for (auto frame : originalFrames) { - if (frame.selMask[0] & Trigger0BIT || frame.selMask[1] & Trigger1BIT) { - hTriggerCounter.Fill(1); - hBCOriginal.Fill(1); - auto p1 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcAO2D == frame.bcAO2D; }); - if (p1 != bcSet.end()) { - hBCOriginal.Fill(2); - } - auto p2 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcEvSel == frame.bcEvSel; }); - if (p2 != bcSet.end()) { - hBCOriginal.Fill(3); - } - bcTuple currentBC(frame.bcAO2D, frame.bcEvSel); - auto p3 = std::find(bcSet.begin(), bcSet.end(), currentBC); - if (p3 == bcSet.end()) { - bcSet.push_back(currentBC); - } else { - hBCOriginal.Fill(4); - } - // std::cout << "------------------------------------------------" << std::endl; - if (frame.GetNum() != 1) { - continue; // Only check singles - } - std::vector skimmedbcs; - int n = 0; - for (auto& skimmedFrame : skimmedFrames) { - if (skimmedFrame.getMin() > frame.getMax()) { - break; - } - if (skimmedFrame.GetNum() != 1) { - continue; // Only check singles - } - if (isClose(frame, skimmedFrame, bcDiffTolerance)) { - bool found = frame.selMask[0] & skimmedFrame.selMask[0] || frame.selMask[1] & skimmedFrame.selMask[1]; - // found = found && (frame.bcAO2D == skimmedFrame.bcAO2D || frame.bcEvSel == skimmedFrame.bcEvSel); - if (found) { - skimmedbcs.push_back({skimmedFrame.bcAO2D, skimmedFrame.bcEvSel}); - n++; - } - } - } - if (n == 0) { - // std::cout << "Trigger not found!!!" << std::endl; - } else if (n == 1) { - hSinglePairCheck.Fill(1); - hBCDiffAO2D.Fill(DoBCSubraction(frame.bcAO2D, skimmedbcs[0].bcAO2D)); - hBCDiffEvSel.Fill(DoBCSubraction(frame.bcEvSel, skimmedbcs[0].bcEvSel)); - if (frame.bcAO2D == skimmedbcs[0].bcAO2D) { - hSinglePairCheck.Fill(2); - } - if (frame.bcEvSel == skimmedbcs[0].bcEvSel) { - hSinglePairCheck.Fill(3); - if (frame.bcAO2D == skimmedbcs[0].bcAO2D) { - hSinglePairCheck.Fill(4); - } - } - } else { - // std::cout << "Unexpected trigger!!! n=" << n << std::endl; - hMultiPairCheck.Fill(1); - for (auto skimmedbc : skimmedbcs) { - hMultiPairCheck.Fill(2); - if (frame.bcAO2D == skimmedbc.bcAO2D) { - hMultiPairCheck.Fill(3); - } - if (frame.bcEvSel == skimmedbc.bcEvSel) { - hMultiPairCheck.Fill(4); - if (frame.bcAO2D == skimmedbc.bcAO2D) { - hMultiPairCheck.Fill(5); - } - } - } - } - hNumCounter.Fill(n); - } - } - auto t6 = std::chrono::steady_clock::now(); - int d3 = std::chrono::duration_cast(t6 - t5).count(); - std::cout << "Search Time: " << d3 << std::endl; - - bcSet.clear(); - for (auto& skimmedFrame : skimmedFrames) { - if (skimmedFrame.selMask[0] & Trigger0BIT || skimmedFrame.selMask[1] & Trigger1BIT) { - hTriggerCounter.Fill(2); - hBCSkimmed.Fill(1); - auto p1 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcAO2D == skimmedFrame.bcAO2D; }); - if (p1 != bcSet.end()) { - hBCSkimmed.Fill(2); - } - auto p2 = std::find_if(bcSet.begin(), bcSet.end(), [&](const auto& val) { return val.bcEvSel == skimmedFrame.bcEvSel; }); - if (p2 != bcSet.end()) { - hBCSkimmed.Fill(3); - } - bcTuple currentBC(skimmedFrame.bcAO2D, skimmedFrame.bcEvSel); - auto p3 = std::find(bcSet.begin(), bcSet.end(), currentBC); - if (p3 == bcSet.end()) { - bcSet.push_back(currentBC); - } else { - hBCSkimmed.Fill(4); - } + // Do checks for trigger + for (int trgID = 0; trgID < labels.size(); trgID++) { + if (trgID == 77 || trgID == 78 || trgID == 79) { + // checkBCForSelectedTrg(originalAllFrames[trgID], skimmedAllFrames[trgID], runNumber, labels[trgID]); } } - - TFile fout(outputFileName, "RECREATE"); - fout.cd(); - hTriggerCounter.Write(); - hBCOriginal.Write(); - hBCSkimmed.Write(); - hNumCounter.Write(); - hSinglePairCheck.Write(); - hBCDiffAO2D.Write(); - hBCDiffEvSel.Write(); - hMultiPairCheck.Write(); - fout.Close(); }