Skip to content

Fix TAGH/TAGM tagger-overlap efficiency double-counting (closes #850) #1153

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sdobbs merged 3 commits into from
Jul 29, 2026
Merged

Fix TAGH/TAGM tagger-overlap efficiency double-counting (closes #850) #1153

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
987e324
Fix TAGGEDMCGEN denominator collapse in TAGH/TAGM overlap regions
lihaoahil Jul 15, 2026
de1a751
Fill_ThrownTree: write one row per TAGGEDMCGEN entry
lihaoahil Jul 15, 2026
57136b3
Fill_ThrownTree: add ThrownBeam__IsTAGH branch and old/new toggle
lihaoahil Jul 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 35 additions & 10 deletions src/libraries/ANALYSIS/DEventWriterROOT.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ static bool STORE_MC_TRAJECTORIES = false;

static bool STORE_SC_VETO_INFO = false;
static bool STORE_THROWN_DECAYING_PARTICLES = true;
static bool STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS = true; //true = current (post-c3deea4b6) behavior: 1 row per TAGGEDMCGEN entry; false = pre-c3deea4b6: collapse to TAGGEDMCGEN[0]/MCGEN[0]

void DEventWriterROOT::Initialize(const std::shared_ptr<const JEvent>& locEvent)
{
Expand Down Expand Up @@ -129,6 +130,7 @@ void DEventWriterROOT::Create_ThrownTree(const std::shared_ptr<const JEvent>& lo
// set parameters specifically for thrown trees
// if(japp->Exists("ANALYSIS:STORE_THROWN_DECAYING_PARTICLES"))
japp->GetParameter("ANALYSIS:STORE_THROWN_DECAYING_PARTICLES",STORE_THROWN_DECAYING_PARTICLES); cout << "ANALYSIS:STORE_THROWN_DECAYING_PARTICLES set to " << STORE_THROWN_DECAYING_PARTICLES << ", IGNORE the \"<-- NO DEFAULT! (TYPO?)\" message " << endl;
japp->GetParameter("ANALYSIS:STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS",STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS); cout << "ANALYSIS:STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS set to " << STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS << ", IGNORE the \"<-- NO DEFAULT! (TYPO?)\" message " << endl;

//TTREE BRANCHES
DTreeBranchRegister locBranchRegister;
Expand Down Expand Up @@ -637,6 +639,7 @@ void DEventWriterROOT::Create_Branches_Thrown(DTreeBranchRegister& locBranchRegi
locBranchRegister.Register_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "X4")); //reported at target center
locBranchRegister.Register_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "P4"));
locBranchRegister.Register_Single<Float_t>(Build_BranchName("ThrownBeam", "GeneratedEnergy"));
locBranchRegister.Register_Single<Bool_t>(Build_BranchName("ThrownBeam", "IsTAGH")); //true: beam photon tagged in TAGH; false: TAGM or SYS_NULL (untagged/MCGEN fallback)

//EVENT-WIDE INFO
locBranchRegister.Register_Single<ULong64_t>("NumPIDThrown_FinalState"); //19 digits
Expand Down Expand Up @@ -1188,22 +1191,43 @@ void DEventWriterROOT::Fill_ThrownTree(const std::shared_ptr<const JEvent>& locE
vector<const DMCTrajectoryPoint*> locDMCTrajectoryPoints;
locEvent->Get(locDMCTrajectoryPoints);

const DBeamPhoton* locTaggedMCGenBeam = locTaggedMCGenBeams.empty() ? locMCGenBeams[0] : locTaggedMCGenBeams[0]; //if empty: will have to do.
//if empty: will have to do.
//If STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS: write one row per TAGGEDMCGEN entry instead of collapsing
//to locTaggedMCGenBeams[0]: in TAGH/TAGM overlap regions TAGGEDMCGEN correctly holds one entry per
//truth-tagged system (mirroring DBeamPhoton_factory_TRUTH), and the reconstructed side is likewise
//left un-merged, so collapsing the thrown side to a single row would undercount the "thrown"
//denominator relative to the (correctly doubled) reconstructed numerator in any analysis that
//divides by this tree's row count instead of querying DBeamPhoton:TAGGEDMCGEN object counts
//directly. Otherwise (flag off), reproduce the pre-c3deea4b6 single-row-per-event behavior for
//consumers that expect exactly one Thrown_Tree row per event.
vector<const DBeamPhoton*> locThrownTreeBeams;
if(STORE_TAGGEDMCGEN_PER_SYSTEM_ROWS)
{
locThrownTreeBeams = locTaggedMCGenBeams.empty() ?
vector<const DBeamPhoton*>(1, locMCGenBeams[0]) : locTaggedMCGenBeams;
}
else
{
locThrownTreeBeams = vector<const DBeamPhoton*>(1, locTaggedMCGenBeams.empty() ? locMCGenBeams[0] : locTaggedMCGenBeams[0]);
}

DEvent::GetLockService(locEvent)->RootWriteLock();

//primary event info
dThrownTreeFillData.Fill_Single<UInt_t>("RunNumber", locEvent->GetRunNumber());
dThrownTreeFillData.Fill_Single<ULong64_t>("EventNumber", locEvent->GetEventNumber());
for(const DBeamPhoton* locTaggedMCGenBeam : locThrownTreeBeams)
{
//primary event info
dThrownTreeFillData.Fill_Single<UInt_t>("RunNumber", locEvent->GetRunNumber());
dThrownTreeFillData.Fill_Single<ULong64_t>("EventNumber", locEvent->GetEventNumber());

//throwns
Fill_ThrownInfo(&dThrownTreeFillData, locMCReaction, locTaggedMCGenBeam, locMCThrownsToSave, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying, locDMCTrajectoryPoints);
//throwns
Fill_ThrownInfo(&dThrownTreeFillData, locMCReaction, locTaggedMCGenBeam, locMCThrownsToSave, locThrownIndexMap, locNumPIDThrown_FinalState, locPIDThrown_Decaying, locDMCTrajectoryPoints);

//Custom Branches
Fill_CustomBranches_ThrownTree(&dThrownTreeFillData, locEvent, locMCReaction, locMCThrownsToSave);
//Custom Branches
Fill_CustomBranches_ThrownTree(&dThrownTreeFillData, locEvent, locMCReaction, locMCThrownsToSave);

//FILL TTREE
dThrownTreeInterface->Fill(dThrownTreeFillData);
//FILL TTREE
dThrownTreeInterface->Fill(dThrownTreeFillData);
}


DEvent::GetLockService(locEvent)->RootUnLock();
Expand Down Expand Up @@ -1701,6 +1725,7 @@ void DEventWriterROOT::Fill_ThrownInfo(DTreeFillData* locTreeFillData, const DMC
DLorentzVector locThrownBeamP4 = locTaggedMCGenBeam->lorentzMomentum();
TLorentzVector locThrownBeamTP4(locThrownBeamP4.Px(), locThrownBeamP4.Py(), locThrownBeamP4.Pz(), locThrownBeamP4.E());
locTreeFillData->Fill_Single<TLorentzVector>(Build_BranchName("ThrownBeam", "P4"), locThrownBeamTP4);
locTreeFillData->Fill_Single<Bool_t>(Build_BranchName("ThrownBeam", "IsTAGH"), locTaggedMCGenBeam->dSystem == SYS_TAGH);

//THROWN PRODUCTS
locTreeFillData->Fill_Single<UInt_t>("NumThrown", locMCThrowns.size());
Expand Down
73 changes: 50 additions & 23 deletions src/libraries/PID/DBeamPhoton_factory_TAGGEDMCGEN.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,39 +27,66 @@ void DBeamPhoton_factory_TAGGEDMCGEN::Process(const std::shared_ptr<const JEvent
if(locMCReactions.empty())
return; //Not a thrown event

//Get the MCGEN beam
const DBeamPhoton* locMCGenBeam;
event->GetSingle(locMCGenBeam, "MCGEN");

//See if it was tagged
auto locSystem = locMCGenBeam->dSystem;
if(locSystem == SYS_NULL)
return; //Nope, no objects to create
//Get every primary (bg == 0) truth-tagged system for the thrown photon.
//In TAGH/TAGM overlap regions the photon can legitimately register truth
//hits in BOTH systems (see DBeamPhoton_factory_TRUTH), just as a real
//photon in that energy range produces coincident TAGH+TAGM hits and
//un-merged reconstructed DBeamPhoton candidates (see DBeamPhoton_factory).
//"MCGEN" only reports a single canonical system (preferring TAGM) since
//it exists to represent the one physical photon's kinematics, not to
//enumerate which tagger systems it was tagged in - so it can't be used
//here without silently dropping the second, equally real, truth match.
vector<const DBeamPhoton*> locTruthPhotons;
event->Get(locTruthPhotons, "TRUTH");

//Get reconstructed beam photons
vector<const DBeamPhoton*> locBeamPhotons;
event->Get(locBeamPhotons);

//Loop over beam photons
double locBestDeltaT = 9.9E9;
const DBeamPhoton* locBestPhoton = nullptr;
for(auto& locBeamPhoton : locBeamPhotons)
for(auto& locTruthPhoton : locTruthPhotons)
{
if(locBeamPhoton->dSystem != locSystem)
continue;
if(locBeamPhoton->dCounter != locMCGenBeam->dCounter)
//only consider the primary photon's own hits, not accidentals/pileup
bool locIsPrimary = false;
vector<const DTAGMHit*> locTAGMHits;
locTruthPhoton->Get(locTAGMHits);
for(auto locTAGMHit : locTAGMHits)
{
if(locTAGMHit->bg == 0)
locIsPrimary = true;
}
vector<const DTAGHHit*> locTAGHHits;
locTruthPhoton->Get(locTAGHHits);
for(auto locTAGHHit : locTAGHHits)
{
if(locTAGHHit->bg == 0)
locIsPrimary = true;
}
if(!locIsPrimary)
continue;

auto locDeltaT = fabs(locMCGenBeam->time() - locBeamPhoton->time());
if(locDeltaT >= locBestDeltaT)
auto locSystem = locTruthPhoton->dSystem;
if(locSystem == SYS_NULL)
continue;
locBestDeltaT = locDeltaT;
locBestPhoton = locBeamPhoton;
}

if(locBestPhoton == nullptr)
return; //Uh oh. Shouldn't be possible.
//Loop over beam photons
double locBestDeltaT = 9.9E9;
const DBeamPhoton* locBestPhoton = nullptr;
for(auto& locBeamPhoton : locBeamPhotons)
{
if(locBeamPhoton->dSystem != locSystem)
continue;
if(locBeamPhoton->dCounter != locTruthPhoton->dCounter)
continue;

mData.push_back(new DBeamPhoton(*locBestPhoton));
auto locDeltaT = fabs(locTruthPhoton->time() - locBeamPhoton->time());
if(locDeltaT >= locBestDeltaT)
continue;
locBestDeltaT = locDeltaT;
locBestPhoton = locBeamPhoton;
}

if(locBestPhoton != nullptr)
mData.push_back(new DBeamPhoton(*locBestPhoton));
}
}

2 changes: 2 additions & 0 deletions src/libraries/PID/DBeamPhoton_factory_TAGGEDMCGEN.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@
#include <PID/DBeamPhoton.h>
#include "PID/DMCReaction.h"
#include "DANA/DStatusBits.h"
#include "TAGGER/DTAGHHit.h"
#include "TAGGER/DTAGMHit.h"

class DBeamPhoton_factory_TAGGEDMCGEN:public JFactoryT<DBeamPhoton>{
public:
Expand Down