Skip to content

Commit

Permalink
Update trigger BC checks (AliceO2Group#7108)
Browse files Browse the repository at this point in the history 
* Update trigger BC checks

* Delete useless output

* Small fixes

* Small fixes

* Small fixes

* Small fixes
  • Loading branch information
@wang-yuanzhe
wang-yuanzhe authored Aug 1, 2024
1 parent 36a905f commit cbbf994
Showing 1 changed file with 148 additions and 65 deletions.
213 changes: 148 additions & 65 deletions EventFiltering/macros/checkBCrangesSkimming.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,11 +24,10 @@ using o2::InteractionRecord;
using o2::dataformats::IRFrame;

// Set the bit of trigger which need to be checked
const ULong64_t trigger0Bit = BIT(54);
const ULong64_t trigger1Bit = 0;
const int bcDiffTolerance = 0;
const ULong64_t Trigger0BIT = BIT(54);
const ULong64_t Trigger1BIT = 0;
const ULong64_t bcDiffTolerance = 100;
const char outputFileName[15] = "output.root";
bool skipNoneDuplicate = false;

struct bcTuple {
bcTuple(ULong64_t bcAO2D, ULong64_t bcEvSel) : bcAO2D(bcAO2D), bcEvSel(bcEvSel) {}
Expand All @@ -43,8 +42,20 @@ struct bcTuple {
struct selectedFrames : public IRFrame {
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
bool isSingle() { return numSameTriggerInNearbyBCs == 0; }
void SetNum(int n) { numSameTriggerInNearbyBCs = n; }
int GetNum() { return numSameTriggerInNearbyBCs; }
};

bool isClose(selectedFrames a, selectedFrames b, ULong64_t bcDiffTolerance)
{
if (a.getMin() > b.getMax() + bcDiffTolerance || a.getMax() < b.getMin() - bcDiffTolerance)
return false;
else
return true;
}

std::vector<selectedFrames> getSelectedFrames(TFile& file, ULong64_t trigger0Bit, ULong64_t trigger1Bit)
{
std::vector<selectedFrames> selectedFrames;
Expand Down Expand Up @@ -87,10 +98,39 @@ std::vector<selectedFrames> getSelectedFrames(TFile& file, ULong64_t trigger0Bit
return selectedFrames;
}

void checkNearbyBCs(std::vector<selectedFrames>& frames, ULong64_t bcDiffTolerance)
{
std::sort(frames.begin(), frames.end(), [](const selectedFrames& a, const selectedFrames& b) {
return a.getMin() < b.getMin();
});
int firstID = 0;
for (auto& currentFrame : frames) {
int num = 0;
bool isFirst = true;
for (int i = firstID; i < frames.size(); i++) {
auto& frame = frames[i];
if (frame.getMin() > currentFrame.getMax() + bcDiffTolerance) {
break;
}
if (isClose(currentFrame, frame, bcDiffTolerance)) {
isFirst = false;
bool found = currentFrame.selMask[0] & frame.selMask[0] || currentFrame.selMask[1] & frame.selMask[1];
if (found) {
num++;
}
} else {
if (isFirst) {
firstID = i;
}
}
}
currentFrame.SetNum(num);
}
}

// Calulate the ratio of duplicate triggers
void checkDuplicateTrigger(std::string AnaFileName = "AnalysisResults.root", std::string originalFileName = "bcRanges_fullrun.root", std::string skimmedFileName = "bcRanges_fullrun_skimmed.root")
{

std::string runNumber = "";
std::regex re("/5[0-9]*");
std::smatch match;
Expand All @@ -116,7 +156,7 @@ void checkDuplicateTrigger(std::string AnaFileName = "AnalysisResults.root", std
TFile originalFile(originalFileName.c_str(), "READ");
TFile skimmedFile(skimmedFileName.c_str(), "READ");
std::vector<std::string> sel_labels;
std::vector<double> numOriginal, numSkimmed, numOriginalDuplicate, numSkimmedDuplicate;
std::vector<int> numOriginal, numSkimmed, numOriginalSingle, numSkimmedSingle, numOriginalDouble, numSkimmedDouble, numOriginalMultiple, numSkimmedMultiple;
for (int i = 0; i < labels.size(); i++) {
ULong64_t trigger0Bit = 0, trigger1Bit = 0;
int triggerBit = binNum[i] - 2;
Expand All @@ -127,81 +167,127 @@ void checkDuplicateTrigger(std::string AnaFileName = "AnalysisResults.root", std
}
// For Original dataset
std::vector<bcTuple> bcSet;
double noriginal{0}, nskimmed{0}, noriginalduplicate{0}, nskimmedduplicate{0};
auto Frames = getSelectedFrames(originalFile, trigger0Bit, trigger1Bit);
for (auto frame : Frames) {
noriginal++;
bcTuple currentBC(frame.bcAO2D, frame.bcEvSel);
auto p = std::find(bcSet.begin(), bcSet.end(), currentBC);
if (p == bcSet.end()) {
bcSet.push_back(currentBC);
std::vector<ULong64_t> bcFullSet;
int noriginal{0}, nskimmed{0}, noriginalsingle{0}, nskimmedsingle{0}, noriginaldouble{0}, nskimmeddouble{0}, noriginalmultiple{0}, nskimmedmultiple{0};
auto originalFrames = getSelectedFrames(originalFile, trigger0Bit, trigger1Bit);
checkNearbyBCs(originalFrames, bcDiffTolerance);
noriginal = originalFrames.size();
for (auto originalFrame : originalFrames) {
if (originalFrame.GetNum() == 0) {
std::cerr << "Unexpected trigger!!! " << std::endl;
} else if (originalFrame.GetNum() == 1) {
noriginalsingle++;
} else if (originalFrame.GetNum() == 2) {
noriginaldouble++;
} else {
noriginalduplicate++;
noriginalmultiple++;
}
}
// For skimmed dataset
bcSet.clear();
auto skimmedFrames = getSelectedFrames(skimmedFile, trigger0Bit, trigger1Bit);
checkNearbyBCs(skimmedFrames, bcDiffTolerance);
nskimmed = skimmedFrames.size();
for (auto& skimmedFrame : skimmedFrames) {
nskimmed++;
bcTuple currentBC(skimmedFrame.bcAO2D, skimmedFrame.bcEvSel);
auto p = std::find(bcSet.begin(), bcSet.end(), currentBC);
if (p == bcSet.end()) {
bcSet.push_back(currentBC);
if (skimmedFrame.GetNum() == 0) {
std::cerr << "Unexpected trigger!!! " << std::endl;
} else if (skimmedFrame.GetNum() == 1) {
nskimmedsingle++;
} else if (skimmedFrame.GetNum() == 2) {
nskimmeddouble++;
} else {
nskimmedduplicate++;
nskimmedmultiple++;
}
}
if (!skipNoneDuplicate || noriginalduplicate != 0 || nskimmedduplicate != 0) {
sel_labels.push_back(labels[i]);
numOriginal.push_back(noriginal);
numOriginalDuplicate.push_back(noriginalduplicate);
numSkimmed.push_back(nskimmed);
numSkimmedDuplicate.push_back(nskimmedduplicate);
}

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);
}
originalFile.Close();
skimmedFile.Close();

TH1D hOriginalTotal("hOriginalTotal", "AO2D Original;;Number of events", sel_labels.size(), 0, sel_labels.size());
TH1D hOriginalDuplicate("hOriginalDuplicate", "Duplicate Trigger Original;;Number of events", sel_labels.size(), 0, sel_labels.size());
TH1D hOriginalRatio("hOriginalRatio", (runNumber + " Original;;Duplicate / Total").data(), sel_labels.size(), 0, sel_labels.size());
TH1D hSkimmedTotal("hSkimmedTotal", "AO2D Skimmed;;Number of events", sel_labels.size(), 0, sel_labels.size());
TH1D hSkimmedDuplicate("hSkimmedDuplicate", "Duplicate Trigger Skimmed;;Number of events", sel_labels.size(), 0, sel_labels.size());
TH1D hSkimmedRatio("hSkimmedRatio", (runNumber + " Skimmed;;Duplicate / Total").data(), sel_labels.size(), 0, sel_labels.size());
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());

for (int i = 0; i < sel_labels.size(); i++) {
hOriginalTotal.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str());
hOriginalDuplicate.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str());
hOriginalRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str());
hOriginalSingles.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]);
hOriginalDuplicate.SetBinContent(i + 1, numOriginalDuplicate[i]);
if (hOriginalTotal.GetBinContent(i + 1) > 0) {
hOriginalRatio.SetBinContent(i + 1, hOriginalDuplicate.GetBinContent(i + 1) / hOriginalTotal.GetBinContent(i + 1));
hOriginalSingles.SetBinContent(i + 1, numOriginalSingle[i]);
hOriginalDoubles.SetBinContent(i + 1, numOriginalDouble[i]);
hOriginalMultiples.SetBinContent(i + 1, numOriginalMultiple[i]);
if (numOriginal[i] > 0) {
hOriginalSinglesRatio.SetBinContent(i + 1, static_cast<double>(numOriginalSingle[i]) / numOriginal[i]);
hOriginalDoublesRatio.SetBinContent(i + 1, static_cast<double>(numOriginalSingle[i]) / numOriginal[i]);
hOriginalMultiplesRatio.SetBinContent(i + 1, static_cast<double>(numOriginalMultiple[i]) / numOriginal[i]);
} else {
hOriginalRatio.SetBinContent(i + 1, 0);
hOriginalSinglesRatio.SetBinContent(i + 1, 0);
hOriginalDoublesRatio.SetBinContent(i + 1, 0);
hOriginalMultiplesRatio.SetBinContent(i + 1, 0);
}

hSkimmedTotal.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str());
hSkimmedDuplicate.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str());
hSkimmedRatio.GetXaxis()->SetBinLabel(i + 1, sel_labels[i].c_str());
hSkimmedSingles.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]);
hSkimmedDuplicate.SetBinContent(i + 1, numSkimmedDuplicate[i]);
if (hSkimmedTotal.GetBinContent(i + 1) > 0) {
hSkimmedRatio.SetBinContent(i + 1, hSkimmedDuplicate.GetBinContent(i + 1) / hSkimmedTotal.GetBinContent(i + 1));
hSkimmedSingles.SetBinContent(i + 1, numSkimmedSingle[i]);
hSkimmedDoubles.SetBinContent(i + 1, numSkimmedDouble[i]);
hSkimmedMultiples.SetBinContent(i + 1, numSkimmedMultiple[i]);
if (numSkimmed[i] > 0) {
hSkimmedSinglesRatio.SetBinContent(i + 1, static_cast<double>(numSkimmedSingle[i]) / numSkimmed[i]);
hSkimmedDoublesRatio.SetBinContent(i + 1, static_cast<double>(numSkimmedDouble[i]) / numSkimmed[i]);
hSkimmedMultiplesRatio.SetBinContent(i + 1, static_cast<double>(numSkimmedMultiple[i]) / numSkimmed[i]);
} else {
hSkimmedRatio.SetBinContent(i + 1, 0);
hSkimmedSinglesRatio.SetBinContent(i + 1, 0);
hSkimmedDoublesRatio.SetBinContent(i + 1, 0);
hSkimmedMultiplesRatio.SetBinContent(i + 1, 0);
}
}

TFile fout(outputFileName, "UPDATE");
fout.cd();
hOriginalTotal.Write();
hOriginalDuplicate.Write();
hOriginalRatio.Write();
hOriginalSingles.Write();
hOriginalDoubles.Write();
hOriginalMultiples.Write();
hOriginalSinglesRatio.Write();
hOriginalDoublesRatio.Write();
hOriginalMultiplesRatio.Write();
hSkimmedTotal.Write();
hSkimmedDuplicate.Write();
hSkimmedRatio.Write();
hSkimmedSingles.Write();
hSkimmedDoubles.Write();
hSkimmedMultiples.Write();
hSkimmedSinglesRatio.Write();
hSkimmedDoublesRatio.Write();
hSkimmedMultiplesRatio.Write();
fout.Close();
}

Expand Down Expand Up @@ -244,29 +330,25 @@ void checkBCrangesSkimming(std::string originalFileName = "bcRanges_fullrun.root
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);
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<std::chrono::milliseconds>(t2 - t1).count();
std::cout << "Readin Time: " << d1 << std::endl;

auto t3 = std::chrono::steady_clock::now();
std::sort(originalFrames.begin(), originalFrames.end(), [](const selectedFrames& a, const selectedFrames& b) {
return a.getMin() < b.getMin();
});
std::sort(skimmedFrames.begin(), skimmedFrames.end(), [](const selectedFrames& a, const selectedFrames& b) {
return a.getMin() < b.getMin();
});
checkNearbyBCs(originalFrames, bcDiffTolerance);
checkNearbyBCs(skimmedFrames, bcDiffTolerance);
auto t4 = std::chrono::steady_clock::now();
int d2 = std::chrono::duration_cast<std::chrono::milliseconds>(t4 - t3).count();
std::cout << "Sort Time: " << d2 << std::endl;

auto t5 = std::chrono::steady_clock::now();
std::vector<bcTuple> bcSet;
for (auto frame : originalFrames) {
if (frame.selMask[0] & trigger0Bit) {
if (frame.selMask[0] & Trigger0BIT || frame.selMask[1] & Trigger1BIT) {
bool found = false;
hTriggerCounter.Fill(1);
hBCOriginal.Fill(1);
Expand All @@ -286,17 +368,18 @@ void checkBCrangesSkimming(std::string originalFileName = "bcRanges_fullrun.root
hBCOriginal.Fill(4);
}
// std::cout << "------------------------------------------------" << std::endl;
frame.getMin() -= bcDiffTolerance;
frame.getMax() += bcDiffTolerance;
if (frame.GetNum() != 1) {
continue; // Only check singles
}
std::vector<bcTuple> skimmedbcs;
int n = 0;
for (auto& skimmedFrame : skimmedFrames) {
if (skimmedFrame.getMin() > frame.getMax()) {
break;
}
if (!frame.isOutside(skimmedFrame)) {
if (isClose(frame, skimmedFrame, bcDiffTolerance)) {
found = frame.selMask[0] & skimmedFrame.selMask[0] || frame.selMask[1] & skimmedFrame.selMask[1];
found = found && (frame.bcAO2D == skimmedFrame.bcAO2D || frame.bcEvSel == skimmedFrame.bcEvSel);
// found = found && (frame.bcAO2D == skimmedFrame.bcAO2D || frame.bcEvSel == skimmedFrame.bcEvSel);
if (found) {
hPairedTriggerCounter.Fill(1);
if (frame.bcAO2D == skimmedFrame.bcAO2D) {
Expand Down Expand Up @@ -353,7 +436,7 @@ void checkBCrangesSkimming(std::string originalFileName = "bcRanges_fullrun.root

bcSet.clear();
for (auto& skimmedFrame : skimmedFrames) {
if (skimmedFrame.selMask[0] & trigger0Bit || skimmedFrame.selMask[1] & trigger1Bit) {
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; });
Expand Down

0 comments on commit cbbf994

Please sign in to comment.