Skip to content

ClusterSearcher and RecoCluster#365

Open
flemmons wants to merge 6 commits intoANNIEsoft:Applicationfrom
flemmons:Application
Open

ClusterSearcher and RecoCluster#365
flemmons wants to merge 6 commits intoANNIEsoft:Applicationfrom
flemmons:Application

Conversation

@flemmons
Copy link
Contributor

-Updated DigitBuilder to be more consistent with and independent of ClusterFinder's corresponding functions
--Added some extra features, such as hit LAPPD count and strip-by-strip LAPPD hits. (some of these are considered experimental/preliminary)
-Added ClusterSearcher to replace ClusterFinder using RecoDigit and RecoCluster classes
-Updated RecoDigit and RecoCluster classes for corresponding use and new features, such as various cluster parameters
-Added NeutronCheck tool as output for RecoCluster information
-Added sample toolchain configfolder for using the new tools

…ClusterFinder's corresponding functions

--Added some extra features, such as hit LAPPD count and strip-by-strip LAPPD hits.
-Added ClusterSearcher to replace ClusterFinder using RecoDigit and RecoCluster classes
-Updated RecoDigit and RecoCluster classes for corresponding use and new features, such as various cluster parameters
-Added NeutronCheck tool as output for RecoCluster information
-Added sample toolchain configfolder for using the new tools
@jminock jminock self-requested a review October 30, 2025 13:03
@jminock jminock self-assigned this Oct 30, 2025
@S81D S81D self-assigned this Oct 30, 2025
}


/*while (!file_singlepe.eof()) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please delete the old commented out code. It unnecessarily clutters the file.

//}
}
}
//FIXME: Need a method to have the 123 be equal to the number of operating detectors
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be fixed now?

double max_angle = 0;
double angle;
Position i_position, j_position;
for (int i = 0; i < fDigitList.size(); i++) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are many separate loops of the same size with the same index going on. Could we combine them? If this is a DataModel function, it could be used across multiple Tools so it would be best to improve efficiency within reason.

LappdNumStrips 60 ## num channels to construct from each LAPPD
LappdStripLength 100 ## relative x position of each LAPPD strip, for dual-sided readout [mm]
LappdStripSeparation 10 ## stripline separation, for calculating relative y position of each LAPPD strip [mm]
PMTMask configfiles/BeamClusterAnalysisMC/DeadPMTIDs_p2v7.txt ## Which PMTs should be masked out? / are dead?
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can point to the most up to date path configfiles/LoadWCSim/DeadPMTIDs_p2v7.txt

if (!fisMC){
int pmtid = recoDigit->GetDetectorID();
unsigned long chankey = pmt_tubeid_to_channelkey[pmtid];
if (pmt_gains[chankey]>0) qep/=pmt_gains[chankey];
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In line 77 (m_variables.Get("SinglePEGains",singlePEgains);) you grab the SPE gains from the gains file then use it to convert from nQ --> pe for data. It doesn't look like you ever specify the gains file in the Config file. Could you add it to ClusterSearcherConfig? Alternatively you can grab that directly from the Store since the LoadGeometry tool populates it for downstream tools (like this one) to use. See PhaseIITreeMaker for an example.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The current configuration in the provided toolchain is for use on MC simulations, which do not require the conversion. So this variable is not necessary within the configuration. I will add it to the readme file.
Though as I start looking at Data in the next couple of weeks, I'll likely move into taking the gains from the store in the next update to this tool. Thank you for the suggestion.

if(fMCParticles->at(i).GetPdgCode()==2112 && fMCParticles->at(i).GetParentPdg()==0) {
fTrueNeutronMult++;

if(fMCParticles->at(i).GetStopTime()>10000) fTrueNeutronDelayed++;
Copy link
Collaborator

@S81D S81D Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might be wise to make this a configurable variable in case someone needs to look for neutrons in a different region of interest. Especially considering 10us is used extensively throughout the code.

//true_Emu*=1000; //GeV->MeV to match other energies(unneeded, possibly)

double theta = truevtx->GetDirection().GetTheta();
double p = sqrt(pow(true_Emu,2)-pow(105.7,2));
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jminock is there a way to grab the momentum and Q^2 from the Store? I thought the LoadGenieEvent tool will store those values for use.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There absolutely is via: m_data->Stores["GenieInfo"]->Get("EventQ2",TrueQ2). Magnitude of muon momentum is not saved to the Store so this is the intended method to get the true muon momentum

Copy link
Collaborator

@jminock jminock left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please make the changes listed. I worry with the large number of pointers if memory is being handled properly. I don't know how necessary all of the pointers are, and I am not enough of an expert on it to make a certain statement. I recommend double checking all are fitting general best practice before this gets sent off to Level 0 review

double maxCharge=0;
int tempPDG = -5;
for (RecoDigit* i_digit : fDigitList) {
cout<<"Digit parent list size: "<<i_digit->GetParents().size()<<endl;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add verbosity conditions to the cout's throughout this file and others. Preferably Log functions if you would like Marcus to give approval down the road


// Default Clustering parameters
fConfig = ClusterSearcher::kPulseHeightAndClusters;
fPmtMinPulseHeight = 5.0; // minimum pulse height (PEs) //Ioana... initial 1.0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The configuration file should take care of the initialization. Leave the default parameters in the configuration file. If any of these NEED to be initialized outside of the configuration file, please do so in the header file upon declaration. It cuts down on clutter

double temp_gain;
while (!file_singlepe.eof()){
file_singlepe >> temp_chankey >> temp_gain;
if (file_singlepe.eof()) break;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line seems unnecessary. While loop already includes this condition

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think it would also prevent the emplace call for the last line of the file.


carryon = 1;
while( carryon ){
carryon = 0;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this for loop inside of a while loop?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seemingly because only one of the two loops had a well-defined range over which to run. This code was inherited from the HitCleaner tool, so I can't answer the why of that decision with complete certainty, but the code flows to my reading, and works as needed.
This seems like a pet peeve, more than a problem, so for the time being, I'll stand by it as is.


/////////////////// Usefull header ///////////////////////
if(verbosity) cout<<"Initializing Tool DigitBuilder"<<endl;
cout<<"Initializing Tool DigitBuilder"<<endl;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this has been mentioned before. Please include verbosity conditions, preferably as Log functions

double sumY=0;
Position pos;
//std::vector<RecoDigits> hullDigits;
std::sort(fDigitList.begin(), fDigitList.end());
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is sorting a vector of pointers with the default comparator going to do what you expect? I see you have a operator< for RecoDigit, but think you still need to provide a comparitor lambda to perform the pointer de-reference:

std::sort(fDigitList.begin(), fDigitList.end(), [](const RecoDigit* a, const RecoDigit* b){ return a < b; });

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm... just going to remove the RecoCluster::ConvexHull function, at least for now. It's not actually useful right now.

}

if( !fgClusterSearcher ){
assert(fgClusterSearcher);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this seems redundant. If the previous allocation didn't succeed it would have thrown bad_alloc.

}

if( fgClusterSearcher ){

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

}

if (!fisMC){
ifstream file_singlepe(singlePEgains.c_str());
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

recommend adding a check file_singlepe.is_open() to catch typos in filename.


// run clustering algorithm
// ========================
std::vector<RecoCluster*>* ClusterList = (std::vector<RecoCluster*>*)(this->RecoClusters(DigitList));
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

again shouldn't need the c-style cast - in fact the local variable ClusterList is just an alias for the member variable fRecoClusters, so seems redundant. The aliasing in this Tool makes it hard to track objects.

// ===================
for(int idigit=0; idigit<int(fSelectByNeighbours->size()); idigit++ ){
RecoDigit* recoDigit = (RecoDigit*)(fSelectByNeighbours->at(idigit));
RecoClusterDigit* clusterDigit = new RecoClusterDigit(recoDigit);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why do the digits need to be on the heap?

//Loop over lines, collect all detector data (should only be one line here)
while (getline(file_singlepe, line)) {
if (verbosity > 3) std::cout << line << std::endl; //has our stuff;
if (line.find("#") != std::string::npos) continue;
Copy link
Collaborator

@marc1uk marc1uk Nov 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

be wary that this will skip lines with trailing comments, as well as commented-out lines.
perhaps if(line.empty() || line[0]=='#') continue is a safer/clearer check

if (verbosity > 3) std::cout << line << std::endl; //has our stuff;
if (line.find("#") != std::string::npos) continue;
std::vector<std::string> DataEntries;
boost::split(DataEntries, line, boost::is_any_of(","), boost::token_compress_on);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i think token_compress_on merges repeated tokens? Is this desirable? If there are repeated tokens, doesn't this suggest an empty column, and by compressing that, your later columns will be shifted?

Log("This HIT'S TIME AND CHARGE: " + to_string(ahit.GetTime()) + ", " + to_string(ahit.GetCharge()),v_debug,verbosity);
double hitTime = ahit.GetTime()*1.0;
if(hitTime>-10 && hitTime<40) {
if(hitTime>-10 && hitTime<70) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

better for magic numbers to be configuration variables, especially if they are subject to change

-Changed RecoCluster and RecoDigit lists in several tools to be vectors of objects, rather than vectors of pointers, to avoid memory complications
-removed unused convex hull function from RecoCluster class
-simplified CalcAS function in RecoCluster class by consolidating for loops
-removed several instances of commented-out code from older versions
-removed several debug outputs
-altered hard-coded time window values in several instances to rely on configuration input
-added use of the true Q2 value from the GenieInfo store to NeutronCheck's output.
-removed several uncontrolled cout lines, and replaced useful ones with Log.
-Tidied the Instance() function of ClusterSearcher
-Added vertex information to NeutronCheck
@jminock
Copy link
Collaborator

jminock commented Jan 13, 2026

Thank you @flemmons for the updates! Unfortunately, this branch has conflicts with the main branch that must be resolved in order to be merged. I will wait for the conflicts to be resolved and for the workflow check to be performed such that ToolAnalysis is confirmed to compile before I review the actual changes to the files

@flemmons
Copy link
Contributor Author

Conflicts were just Factory and Unity listing different added tools. I've cleared them.

@jminock
Copy link
Collaborator

jminock commented Jan 14, 2026

Okay, but could you still resolve those conflicts? Also, it's great that ToolAnalysis can compile in your workspace, but it would be better to get confirmation that ToolAnalysis can compile in a general workspace (on GitHub via workflow). Given that this PR is self contained, there is no reason to not resolve these conflicts. I will look at the other files in the meantime.

@jminock jminock added bug and removed Conflicts labels Jan 16, 2026
@jminock jminock removed the bug label Jan 26, 2026
MRDRecoCut 0
RecoPMTVolCut 0
RecoFVCut 0
ArgonFV 0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not a hold up for merger, but what is ArgonFV? That isn't expected for current the current EventSelector

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ArgonFV is a holdover from a previous experimental study. It...shouldn't be there. Good catch.

@@ -0,0 +1,12 @@
verbosity 0
FluxVersion 0 # use 0 to load genie files based on bnb_annie_0000.root etc files
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The flux files corresponding to version 0 are very outdated. None of the GENIE files I produced use those form of flux files. I don't even know where those flux files are. Pointing this out more as a courtesy for your consideration

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's good to know... And likely explains why I've seen strange behavior in the extracted true q2 values I've been working with since adding them from LoadGenieEvent. I appreciate it.

MCParticleProperties MCParticleProperties ./configfiles/ClusterSearcher/MCParticlePropertiesConfig
DigitBuilder DigitBuilder ./configfiles/ClusterSearcher/DigitBuilderConfig
ClusterSearcher ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcherConfig
ClusterSearcher2 ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcher2Config
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You have to use ClusterSearcher twice? Would a second run of ClusterSearcher overwrite the first use? Sorry, just confused why this is here. Also, how does the information get written out to a file? Does NeutronCheck double as an output file maker? Tools are meant to be singular in purpose, and there already exists output file makers (PhaseIITreeMaker)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ClusterSearcher is designed to be used multiple times, yes. It will add the new clusters to the prior list. This enables different sorts of clusters to be made for muli-level analysis. The current setup is running ClusterSearcher to make concentrated Cherenkov ring-like clusters to analyze the initial muon emission for vertexing and CC verification, and then ClusterSearcher2 collects other clusters without any spatial density requirements to be tested for neutron identification. The different types of clusters are separated by the ClusterMode tag.

As to NeutronCheck, it was meant to be more of a debug tool, akin to VertexGeometryCheck but for the neutron analysis. Its scope grew as I needed more relevant values, but it still doesn't represent a final-state of analysis. In the current ToolAnalysis, it's the only 'end of toolchain' tool that uses RecoClusters, and so adding it was the best option for this toolchain. However, I'm also not sure that PhaseIITreeMaker should be the be-all, end-all output tool for all of ANNIE. Having all output centered in a single tool will make it very cumbersome for focused analyses that only need certain variables. It may be worth considering having separate output tools for separate analyses. As I finalize the neutron multiplicity analysis, there will likely be a more concrete output, whether that's adding the important variables to PhaseIITreeMaker, advancing NeutronCheck to a more immediately useful form, or something else.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gotcha, please document this and how to use ClusterSearcher in a central, accessible place.

As far as the opinion on separate output Tools: having multiple Tools with similar functionality that could be altered via a configuration file is against the intended direction of ANNIE analysis and software. Different sets of Tools for different analyses will have duplicated work, divergent branches, and different standards for different analyses. Yes, PhaseIITreeMaker is cumbersome to work with. But I can guarentee standardizing different Tools with the sane functionality is much more cumbersome. There already is ANNIETreeMaker as another output file maker. I'm not opposed to that Tool being the new standard, but there should be a singular standard across analyses. Making a new separate output Tool will result in: it continuing to advance and develop after you, diverging from other analysis ToolChains, silo'ing off neutron analysis from everything else; or it getting scrapped in favor of shared output in a Tool, losing important work done now, making unnecessary work in the future. Do what you need for your present analysis, but please keep in mind that existing resources are to here to help each of us and the work we do serves ourselves, each other, and those who come after.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is detailed in the ClusterSearcher tool's readme file.

You're right. I'll plan to fold the final analysis outputs into PhaseIITreeMaker or the like as I complete them, and hold NeutronCheck to the level of a debugging check of mid-level values. But, in that role, and for the time being, it is useful, and should not hold back this pull request. Take it as a demonstration of what ClusterSearcher is doing, more than any analysis output.

//}
}
//FIXME: Need a method to have the 123 be equal to the number of operating detectors
double ucharge_balance = sqrt((total_QSquared) / (total_Q * total_Q) - (1. / 123.));
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

123 is a magic number. While it doesn't impact functionality for now, it can be troublesome for documentation and updating it. What does 123 represent? What is an "operating detector"?

Copy link
Collaborator

@jminock jminock left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please address all the comments. I left some that are more general questions that would be nice to have accessible documentation addressing, and shouldn't be an issue to create. There are some comments that ask about memory. Please address them as they do impact functionality

for(int idigit=0; idigit<myCluster->GetNDigits(); idigit++ ){
RecoDigit* myDigit = (RecoDigit*)(myCluster->GetDigit(idigit));
for(int idigit=0; idigit<myCluster.GetNDigits(); idigit++ ){
RecoDigit* myDigit = new RecoDigit;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are creating pointers in a nested loop. How are you freeing the memory once they're done? This looks like it can and will cause memory issues.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does this Tool exist when there already exists output file maker Tools?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As to NeutronCheck, it was meant to be more of a debug tool, akin to VertexGeometryCheck but for the neutron analysis. Its scope grew as I needed more relevant values, but it still doesn't represent a final-state of analysis. In the current ToolAnalysis, it's the only 'end of toolchain' tool that uses RecoClusters, and so adding it was the best option for this toolchain. However, I'm also not sure that PhaseIITreeMaker should be the be-all, end-all output tool for all of ANNIE. Having all output centered in a single tool will make it very cumbersome for focused analyses that only need certain variables. It may be worth considering having separate output tools for separate analyses. As I finalize the neutron multiplicity analysis, there will likely be a more concrete output, whether that's adding the important variables to PhaseIITreeMaker, advancing NeutronCheck to a more immediately useful form, or something else.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants