From 9db8364450085910917fde93efb2e9defd6dd492 Mon Sep 17 00:00:00 2001 From: root Date: Mon, 27 Oct 2025 15:49:33 +0000 Subject: [PATCH 1/5] -Updated DigitBuilder to be more consistent with, but independent of 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 --- DataModel/RecoCluster.cpp | 437 ++++++++- DataModel/RecoCluster.h | 57 ++ DataModel/RecoDigit.h | 7 + UserTools/ClusterSearcher/ClusterSearcher.cpp | 827 ++++++++++++++++++ UserTools/ClusterSearcher/ClusterSearcher.h | 157 ++++ UserTools/ClusterSearcher/README.md | 50 ++ UserTools/DigitBuilder/DigitBuilder.cpp | 536 +++++++++--- UserTools/DigitBuilder/DigitBuilder.h | 17 +- UserTools/Factory/Factory.cpp | 2 + UserTools/NeutronCheck/NeutronCheck.cpp | 491 +++++++++++ UserTools/NeutronCheck/NeutronCheck.h | 177 ++++ UserTools/NeutronCheck/README.md | 71 ++ UserTools/Unity.h | 2 + .../ClusterSearcher/ClusterSearcher2Config | 21 + .../ClusterSearcher/ClusterSearcherConfig | 21 + .../ClusterSearcher/DigitBuilderConfig | 12 + .../ClusterSearcher/LoadGeometryConfig | 10 + configfiles/ClusterSearcher/LoadWCSimConfig | 17 + .../ClusterSearcher/LoadWCSimLAPPDConfig | 11 + .../MCParticlePropertiesConfig | 3 + .../ClusterSearcher/NeutronCheckConfig | 5 + configfiles/ClusterSearcher/ToolChainConfig | 23 + configfiles/ClusterSearcher/ToolsConfig | 9 + 23 files changed, 2860 insertions(+), 103 deletions(-) create mode 100644 UserTools/ClusterSearcher/ClusterSearcher.cpp create mode 100644 UserTools/ClusterSearcher/ClusterSearcher.h create mode 100644 UserTools/ClusterSearcher/README.md create mode 100644 UserTools/NeutronCheck/NeutronCheck.cpp create mode 100644 UserTools/NeutronCheck/NeutronCheck.h create mode 100644 UserTools/NeutronCheck/README.md create mode 100755 configfiles/ClusterSearcher/ClusterSearcher2Config create mode 100755 configfiles/ClusterSearcher/ClusterSearcherConfig create mode 100755 configfiles/ClusterSearcher/DigitBuilderConfig create mode 100644 configfiles/ClusterSearcher/LoadGeometryConfig create mode 100644 configfiles/ClusterSearcher/LoadWCSimConfig create mode 100644 configfiles/ClusterSearcher/LoadWCSimLAPPDConfig create mode 100644 configfiles/ClusterSearcher/MCParticlePropertiesConfig create mode 100644 configfiles/ClusterSearcher/NeutronCheckConfig create mode 100644 configfiles/ClusterSearcher/ToolChainConfig create mode 100755 configfiles/ClusterSearcher/ToolsConfig diff --git a/DataModel/RecoCluster.cpp b/DataModel/RecoCluster.cpp index c63855600..bea56c09f 100644 --- a/DataModel/RecoCluster.cpp +++ b/DataModel/RecoCluster.cpp @@ -1,4 +1,4 @@ -#include "RecoCluster.h" +#include "RecoCluster.h" #include "ANNIERecoObjectTable.h" #include @@ -48,5 +48,440 @@ int RecoCluster::GetNDigits() return fDigitList.size(); } +void RecoCluster::SetClusterMode(int cmode) +{ + fClusterMode = cmode; +} + +int RecoCluster::GetClusterMode() +{ + return fClusterMode; +} + +void RecoCluster::CalcTime() { + double sum = 0; + for (RecoDigit* i_digit : fDigitList) { + sum += i_digit->GetCalTime(); + } + clusterTime = sum / GetNDigits(); + +} + +void RecoCluster::CalcCharge() { + double sum=0; + for (RecoDigit* i_digit : fDigitList) { + sum += i_digit->GetCalCharge(); + } + clusterCharge = sum; +} + +void RecoCluster::CalcCB() { + //calculate unfiltered CB + double total_Q = 0; + double total_QSquared = 0; + for (int i = 0; i < fDigitList.size(); i++) { + //if(unfilteredDigits->at(i)->GetDigitType()==RecoDigit::PMT8inch){ + if (fDigitList.at(i)->GetFilterStatus()) { + double tube_charge = fDigitList.at(i)->GetCalCharge(); + total_Q += tube_charge; + total_QSquared += (tube_charge * tube_charge); + //} + } + } + //FIXME: Need a method to have the 123 be equal to the number of operating detectors + double charge_balance = sqrt((total_QSquared) / (total_Q * total_Q) - (1. / 123.)); + + clusterCB=charge_balance; +} + +//Calculate Angular span +// - AS0 is the maximum angle between any two hit PMTs +// - AS1 is the largest angle between the greatest hit by charge and any of the other hits +// - ASC is a charge-weighted version of AS0 +void RecoCluster::CalcAS() { + + double max_angle = 0; + double angle; + Position i_position, j_position; + for (int i = 0; i < fDigitList.size(); i++) { + i_position = fDigitList.at(i)->GetPosition(); + for (int j = 0; j < fDigitList.size(); j++) { + if (j == i)continue; + j_position = fDigitList.at(j)->GetPosition(); + + angle = j_position.Angle(i_position); + if (angle > max_angle)max_angle = angle; + + } + } + AS0 = max_angle; + + + + double max_charge = 0; + double max_angle2 = 0; + int max_index = 0; + Position max_position; + for (int i = 0; i < fDigitList.size(); i++) { + if (fDigitList.at(i)->GetCalCharge() > max_charge) { + max_charge = fDigitList.at(i)->GetCalCharge(); + max_index = i; + max_position = fDigitList.at(i)->GetPosition(); + } + } + + + for (int i = 0; i < fDigitList.size(); i++) { + if (i == max_index)continue; + i_position = fDigitList.at(i)->GetPosition(); + angle = i_position.Angle(max_position); + if (angle > max_angle2)max_angle2 = angle; + } + AS1 = max_angle2; + + + double max_chargeangle = 0; + double chargeangle; + double iq,jq; + for (int i = 0; i < fDigitList.size(); i++) { + i_position = fDigitList.at(i)->GetPosition(); + iq=fDigitList.at(i)->GetCalCharge(); + for(int j = 0; j < fDigitList.size(); j++) { + if(j==i)continue; + j_position = fDigitList.at(j)->GetPosition(); + jq = fDigitList.at(j)->GetCalCharge(); + chargeangle=iq*jq*j_position.Angle(i_position); + if(chargeangle>max_chargeangle)max_chargeangle=chargeangle; + } + } + + ASC=max_chargeangle; + +} + +//Retrieve Angular Span value +double RecoCluster::GetAS(int mode) { + if(mode==0)return AS0; + else if(mode==1)return AS1; + else return 0; +} + +void RecoCluster::CalcParameters() { + this->CalcTime(); + this->CalcCharge(); + this->CalcCB(); + this->CalcAS(); + this->CalcSA(); + this->CalcCA(); + this->CalcCV(); + this->CalcAMD(); + + //Other needed parameters here +} + +void RecoCluster::SetParticle(int pID, int pPDG, double eff, double pur) { + bestParticleID=pID; + bestParticlePDG=pPDG; + efficiency=eff; + purity=pur; +} + +int RecoCluster::calcBestParent() { + std::map ParticletoClusterCharge; + std::vector DigiParents; + int parentCandidate; + int maxIndex=0; + double digitCharge; + double maxCharge=0; + int tempPDG = -5; + for (RecoDigit* i_digit : fDigitList) { + cout<<"Digit parent list size: "<GetParents().size()<GetParents().size() == 0) { + cout<<"Digit found with no parents!!!\n"; + } + for(int i=0; iGetParents().size(); i++){ + parentCandidate=i_digit->GetParents().at(i); + cout<<"ClusterParent candidate ID: "<GetCalCharge(); + if (ParticletoClusterCharge.count(parentCandidate) == 0) { + ParticletoClusterCharge.emplace(parentCandidate,digitCharge); + } + else { + ParticletoClusterCharge.at(parentCandidate)+=digitCharge; + } + } + } + for (std::pairapair : ParticletoClusterCharge) { + if (apair.second > maxCharge) { + maxCharge=apair.second; + maxIndex=apair.first; + + } + } + + /*for (std::pair apair : *fMCParticleIndexMap) { + if (apair.second == maxIndex) { + tempPDG=apair.first; + } + }*/ + + //if(tempPDG==-5) return false; + + cout << "Particle " << maxIndex << " wins with charge " << maxCharge << endl; + + bestParticleID=maxIndex; + //bestParticlePDG=tempPDG; + purity=maxCharge/clusterCharge; + + return maxIndex; +} + +std::vector RecoCluster::convexHull() { + + double tempX, tempY, tempPhi, tempTheta; + std::vector effX, effY; + double sumX=0; + double sumY=0; + Position pos; + //std::vector hullDigits; + std::sort(fDigitList.begin(), fDigitList.end()); + for (RecoDigit* adigit : fDigitList) { + pos = adigit->GetPosition(); + tempX = sin(pos.GetTheta()) * cos(pos.GetPhi()); + tempY = sin(pos.GetTheta()) * sin(pos.GetPhi()); + effX.push_back(tempX); + effY.push_back(tempY); + sumX+=tempX; + sumY+=tempY; + } + + TwoDCenter.SetX(sumX/fDigitList.size()); + TwoDCenter.SetY(sumY/fDigitList.size()); + + TwoDCenter.SetZ(0); //2D projection center; used for circularity test. + + std::vector lower, upper; + double diffX, diffY; + double diffprevX, diffprevY; + double cross; + + + // Lower hull + for (int i = 0; i < fDigitList.size(); i++) { + if (lower.size() >= 2) { + diffprevX = effX[lower[lower.size() - 1]] - effX[lower[lower.size() - 2]]; + diffprevY = effY[lower[lower.size() - 1]] - effY[lower[lower.size() - 2]]; + diffX = effX[i] - effX[lower[lower.size() - 1]]; + diffY = effY[i] - effY[lower[lower.size() - 1]]; + cross = diffprevX * diffY - diffprevY * diffX; + } + while (lower.size() >= 2 && cross <= 0) { + //(lower.size() >= 2 && (lower[lower.size() - 1] - lower[lower.size() - 2]).cross(p - lower[lower.size() - 1]) <= 0) { + lower.pop_back(); + if (lower.size() >= 2) { + diffprevX = effX[lower[lower.size() - 1]] - effX[lower[lower.size() - 2]]; + diffprevY = effY[lower[lower.size() - 1]] - effY[lower[lower.size() - 2]]; + diffX = effX[i] - effX[lower[lower.size() - 1]]; + diffY = effY[i] - effY[lower[lower.size() - 1]]; + cross = diffprevX * diffY - diffprevY * diffX; + } + + } + lower.push_back(i); + std::cout<<"lowersize,effX,effY,finalcross: "<= 0; --i) { + if (upper.size() >= 2) { + diffprevX = effX[upper[upper.size() - 1]] - effX[upper[upper.size() - 2]]; + diffprevY = effY[upper[upper.size() - 1]] - effY[upper[upper.size() - 2]]; + diffX = effX[i] - effX[upper[upper.size() - 1]]; + diffY = effY[i] - effY[upper[upper.size() - 1]]; + cross = diffprevX * diffY - diffprevY * diffX; + } + while (upper.size() >= 2 && cross <= 0) { + //(upper[upper.size() - 1] - upper[upper.size() - 2]).cross(points[i] - upper[upper.size() - 1]) <= 0) { + upper.pop_back(); + if (upper.size() >= 2) { + diffprevX = effX[upper[upper.size() - 1]] - effX[upper[upper.size() - 2]]; + diffprevY = effY[upper[upper.size() - 1]] - effY[upper[upper.size() - 2]]; + diffX = effX[i] - effX[upper[upper.size() - 1]]; + diffY = effY[i] - effY[upper[upper.size() - 1]]; + cross = diffprevX * diffY - diffprevY * diffX; + } + } + upper.push_back(i); + std::cout << "uppersize,effX,effY,finalcross: " << upper.size() << "," << effX[i] << "," << effY[i] << "," << cross << endl; + } + + // Remove the last point of each half because it is repeated at the beginning of the other half + lower.pop_back(); + upper.pop_back(); + + // Concatenate lower and upper hull to get the convex hull + lower.insert(lower.end(), upper.begin(), upper.end()); + + for (int i = 0; i < fDigitList.size(); i++) { + fDigitList.at(i)->ResetFilter(); + } + + for (int i = 0; i < lower.size(); i++) { + hullDigits.push_back(fDigitList[lower[i]]); + fDigitList.at(lower[i])->PassFilter(); + std::cout<<"Adding digit index "<GetPosition() - fDigitList.at(j)->GetPosition()).Unit(); + ChargeVector+=(fDigitList.at(i)->GetCalCharge()*fDigitList.at(j)->GetCalCharge())*ud; + } + } +} + +void RecoCluster::CalcAMD() { + double min_distance = 9999; + double ave_min_distance=0; + double distance; + for (int i = 0; i < fDigitList.size(); i++) { + Position ipos=fDigitList.at(i)->GetPosition(); + for (int j = 0; j < fDigitList.size(); j++) { + if(j == i) continue; + Position jpos=fDigitList.at(j)->GetPosition(); + distance = (jpos-ipos).Mag(); + if(distanceGetPosition(); + double iq=fDigitList.at(i)->GetCalCharge(); + aveX += ipos.X() * iq; + aveY += ipos.Y() * iq; + aveZ += ipos.Z() * iq; + } + aveX/=clusterCharge; + aveY/=clusterCharge; + aveZ/=clusterCharge; + SpatialAverage=Position(aveX,aveY,aveZ); +} + +//Containing Angle: angle from SpatialAverage unit vector which contains 68% of the total charge of the cluster +void RecoCluster::CalcCA() { //AW: Angular Width? 1-sigma angular width? + std::map angle_charge_map; + double angle,charge; + double contained=0; + for (RecoDigit* idigit: fDigitList) { + angle=SpatialAverage.Angle(idigit->GetPosition()); + charge = idigit->GetCalCharge(); + if(angle_charge_map.count(angle)==0) + angle_charge_map.emplace(angle,charge); + else + angle_charge_map.at(angle)+=charge; + } + for (std::pair apair: angle_charge_map) { + ContainingAngle=apair.first; + contained+=apair.second/clusterCharge; + if(contained>0.68)break; + } + return; +} + +//Planarity and Sphericity: measure of how planar or spherical cluster is; currently in-process of translating from reference +/*void ComputePlanaritySphericity_ROOT(const std::vector& digits, + double& P, double& S) +{ + P = 0.0; S = 0.0; + + // 1) Charge-weighted centroid + TVector3 mean(0, 0, 0); + double wsum = 0.0; + for (int i = 0; i lambda = { evals[0], evals[1], evals[2] }; + std::sort(lambda.begin(), lambda.end(), std::greater()); + + const double lambda1 = lambda[0]; // largest + const double lambda2 = lambda[1]; + const double lambda3 = lambda[2]; // smallest + const double sum = lambda1 + lambda2 + lambda3; + + if (sum <= 0.0) return; + + // 4) Planarity and Sphericity + // P = (lambda2 - lambda3) / (lambda1 + lambda2 + lambda3) + // S = (3/2) * lambda3 / (lambda1 + lambda2 + lambda3) + P = (lambda2 - lambda3) / sum; + S = 1.5 * (lambda3 / sum); +} + +// ----------------------- +// Example usage +// ----------------------- +int main() { + std::vector digits; + digits.push_back({ TVector3(1.0, 0.0, 0.0), 2.0 }); + digits.push_back({ TVector3(0.0, 1.0, 0.1), 1.0 }); + digits.push_back({ TVector3(0.0,-1.0,-0.1), 1.5 }); + digits.push_back({ TVector3(-1.0, 0.0, 0.0), 1.2 }); + double P = 0, S = 0; + ComputePlanaritySphericity_ROOT(digits, P, S); + std::cout << "Planarity P = " << P << "\n"; + std::cout << "Sphericity S = " << S << "\n"; + return 0; +}*/ \ No newline at end of file diff --git a/DataModel/RecoCluster.h b/DataModel/RecoCluster.h index 4afb14268..5986a6669 100644 --- a/DataModel/RecoCluster.h +++ b/DataModel/RecoCluster.h @@ -3,6 +3,7 @@ #ifndef RECOCLUSTER_H #define RECOCLUSTER_H #include "RecoDigit.h" +#include "Position.h" #include class RecoCluster : public SerialisableObject { @@ -19,16 +20,72 @@ class RecoCluster : public SerialisableObject { void AddDigit(RecoDigit* digit); RecoDigit* GetDigit(int n); + + void SetClusterMode(int cmode); + + int GetClusterMode(); + + std::vector GetDigitList() {return this->fDigitList;} + int GetNDigits(); bool Print() { cout<<"Number of digits in this cluster : "< convexHull(); + + void SetParticle(int pID,int pPDG, double eff, double pur); + inline int GetBestParent(){return bestParticleID;} + inline void SetPDG(int pPDG) {bestParticlePDG=pPDG;} + inline int GetPDG(){return bestParticlePDG;} + inline double Efficiency(){return efficiency;} + inline double Purity(){return purity;} + private: + double clusterTime; + double clusterCharge; + double clusterCB; + double AS0, AS1, ASC; + double AMD; + Position ChargeVector; + Position SpatialAverage; + double ContainingAngle; + int bestParticleID; + int bestParticlePDG; + double efficiency; + double purity; + + int fClusterMode = -999; std::vector fDigitList; + std::vector hullDigits; + Position TwoDCenter; + protected: template void serialize(Archive & ar, const unsigned int version){ diff --git a/DataModel/RecoDigit.h b/DataModel/RecoDigit.h index 13d12c4d0..af64ee7cf 100644 --- a/DataModel/RecoDigit.h +++ b/DataModel/RecoDigit.h @@ -35,6 +35,10 @@ class RecoDigit : public SerialisableObject{ } ~RecoDigit() {/*ANNIERecoObjectTable::Instance()->DeleteDigit();*/} + bool operator<(const RecoDigit other) const { + return (fPosition.X() == other.GetPosition().X()) ? fPosition.Y() < other.GetPosition().Y() : fPosition.X() < other.GetPosition().X(); + } + inline int GetRegion() const {return fRegion;} inline Position GetPosition() const {return fPosition;} inline double GetCalTime() const {return fCalTime;} @@ -42,6 +46,7 @@ class RecoDigit : public SerialisableObject{ inline bool GetFilterStatus() const {return fIsFiltered;} inline int GetDigitType() const {return fDigitType;} inline int GetDetectorID() const {return fDetectorID;} + inline std::vector GetParents() const {return Parents;} inline void SetRegion(int reg) {fRegion = reg;} inline void SetPosition(Position pos){fPosition = pos;} @@ -52,6 +57,7 @@ class RecoDigit : public SerialisableObject{ inline void SetFilter(bool pass = 1) { fIsFiltered = pass;} inline void ResetFilter() {SetFilter(0);} inline void PassFilter() {SetFilter(1);} + inline void SetParents(std::vector parentsin) { Parents = parentsin; } bool Print() { cout<<"Region : "< Parents; template void serialize(Archive & ar, const unsigned int version){ if(serialise){ diff --git a/UserTools/ClusterSearcher/ClusterSearcher.cpp b/UserTools/ClusterSearcher/ClusterSearcher.cpp new file mode 100644 index 000000000..dda70f2ee --- /dev/null +++ b/UserTools/ClusterSearcher/ClusterSearcher.cpp @@ -0,0 +1,827 @@ +#include "ClusterSearcher.h" + +static ClusterSearcher* fgClusterSearcher = 0; + +ClusterSearcher* ClusterSearcher::Instance() +{ + if( !fgClusterSearcher ){ + fgClusterSearcher = new ClusterSearcher(); + } + + if( !fgClusterSearcher ){ + assert(fgClusterSearcher); + } + + if( fgClusterSearcher ){ + + } + + return fgClusterSearcher; +} + +ClusterSearcher::ClusterSearcher():Tool(){} + +ClusterSearcher::~ClusterSearcher() { + +} + +bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){ + + /////////////////// Usefull header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); //loading config file + //m_variables.Print(); + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + + // Default Clustering parameters + fConfig = ClusterSearcher::kPulseHeightAndClusters; + fPmtMinPulseHeight = 5.0; // minimum pulse height (PEs) //Ioana... initial 1.0 + fPmtNeighbourRadius = 50.0; // clustering window (cm) //Ioana... intial 300.0 + fPmtMinNeighbourDigits = 2; // minimum neighbouring digits //Ioana.... initial 2 + fPmtClusterRadius = 50.0; // clustering window (cm) //Ioana... initia 300.0 + fMinClusterDigits = 2; // minimum clustered digits + fPmtTimeWindowN = 10; // timing window for neighbours (ns) + fPmtTimeWindowC = 10; // timing window for clusters (ns) + fPmtMinHitsPerCluster = -1; //min # of hits per cluster //Ioana + fisMC = 1; //default: MC + + fLappdMinPulseHeight = -1.0; // minimum pulse height (PEs) //Ioana... initial 1.0 + fLappdNeighbourRadius = 25.0; // clustering window (cm) //Ioana... intial 300.0 + fLappdMinNeighbourDigits = 5; // minimum neighbouring digits //Ioana.... initial 2 + fLappdClusterRadius = 25.0; // clustering window (cm) //Ioana... initia 300.0 + fLappdTimeWindowN = 10; // timing window for neighbours (ns) + fLappdTimeWindowC = 10; // timing window for clusters (ns) + fLappdMinHitsPerCluster = 5; //min # of hits per cluster //Ioana + + /// Get the Tool configuration variables + m_variables.Get("verbosity",verbosity); + m_variables.Get("IsMC",fisMC); + m_variables.Get("ClusterMode",fClusterMode); + m_variables.Get("Config",fConfig); + m_variables.Get("PmtMinPulseHeight", fPmtMinPulseHeight); + m_variables.Get("PmtNeighbourRadius", fPmtNeighbourRadius); + m_variables.Get("PmtMinNeighbourDigits", fPmtMinNeighbourDigits); + m_variables.Get("PmtClusterRadius", fPmtClusterRadius); + m_variables.Get("PmtTimeWindowN", fPmtTimeWindowN); + m_variables.Get("PmtTimeWindowC", fPmtTimeWindowC); + m_variables.Get("PmtMinHitsPerCluster", fPmtMinHitsPerCluster); + m_variables.Get("LappdMinPulseHeight", fLappdMinPulseHeight); + m_variables.Get("LappdNeighbourRadius", fLappdNeighbourRadius); + m_variables.Get("LappdMinNeighbourDigits", fLappdMinNeighbourDigits); + m_variables.Get("LappdClusterRadius", fLappdClusterRadius ); + m_variables.Get("LappdTimeWindowN", fLappdTimeWindowN ); + m_variables.Get("LappdTimeWindowC", fLappdTimeWindowC ); + m_variables.Get("LappdMinHitsPerCluster", fLappdMinHitsPerCluster ); + m_variables.Get("MinClusterDigits", fMinClusterDigits ); + m_variables.Get("SinglePEGains",singlePEgains); + m_variables.Get("FirstRun",fFirstRun); + + /// Fill map with settings of ClusterSearcher + fClusteringParam = new std::map; + fClusteringParam->emplace("ClusterMode",fClusterMode); + fClusteringParam->emplace("Config",fConfig); + fClusteringParam->emplace("PmtMinPulseHeight",fPmtMinPulseHeight); + fClusteringParam->emplace("PmtNeighbourRadius",fPmtNeighbourRadius); + fClusteringParam->emplace("PmtMinNeighbourDigits",fPmtMinNeighbourDigits); + fClusteringParam->emplace("PmtClusterRadius",fPmtClusterRadius); + fClusteringParam->emplace("PmtTimeWindowN",fPmtTimeWindowN); + fClusteringParam->emplace("PmtTimeWindowC",fPmtTimeWindowC); + fClusteringParam->emplace("LappdMinPulseHeight",fLappdMinPulseHeight); + fClusteringParam->emplace("LappdNeighbourRadius",fLappdNeighbourRadius); + fClusteringParam->emplace("LappdMinNeighbourDigits",fLappdMinNeighbourDigits); + fClusteringParam->emplace("LappdClusterRadius",fLappdClusterRadius); + fClusteringParam->emplace("LappdTimeWindowN",fLappdTimeWindowN); + fClusteringParam->emplace("LappdTimeWindowC",fLappdTimeWindowC); + fClusteringParam->emplace("MinClusterDigits",fMinClusterDigits); + + if (fConfig!=0 && fConfig !=1 && fConfig !=2 && fConfig !=3 && fConfig !=4){ + Log("ClusterSearcher tool: Configuration <"+std::to_string(fConfig)+"> not recognized. Setting Config 3 (kPulseHeightAndClusters)",v_error,verbosity); + fConfig = ClusterSearcher::kPulseHeightAndClusters; + } + + if (!fisMC){ + ifstream file_singlepe(singlePEgains.c_str()); + unsigned long temp_chankey; + double temp_gain; + while (!file_singlepe.eof()){ + file_singlepe >> temp_chankey >> temp_gain; + if (file_singlepe.eof()) break; + pmt_gains.emplace(temp_chankey,temp_gain); + } + file_singlepe.close(); + m_data->CStore.Get("pmt_tubeid_to_channelkey",pmt_tubeid_to_channelkey); + } + + // vector of selected digits + fSelectAll = new std::vector; + fSelectByPulseHeight = new std::vector; + fSelectByNeighbours = new std::vector; + fSelectByClusters = new std::vector; + // only for test + fSelectByTruthInfo = new std::vector; + // vector of clusters + fClusterList = new std::vector; + fRecoClusters = new std::vector; + + //Set clustering parameters in the RecoEvent store + std::map* pre_ClusteringParam = nullptr; // read existing container for parameters + bool cluster_parameter_status = m_data->Stores.at("RecoEvent")->Get("ClusteringParameters", pre_ClusteringParam); + if(!cluster_parameter_status) m_data->Stores.at("RecoEvent")->Set("ClusteringParameters", fClusteringParam); + else { + pre_ClusteringParam->insert(fClusteringParam->begin(), fClusteringParam->end()); + //m_data->Stores.at("RecoEvent")->Set("ClusteringParameters", pre_ClusteringParam); + } + + return true; +} + + +bool ClusterSearcher::Execute(){ + if (fRecoClusters) Log("Prelim check1: fRecoClusters size: " + to_string(fRecoClusters->size()), v_debug, verbosity); + if (fRecoClusters->size()!=0) fRecoClusters->clear(); + + if(fFirstRun && pre_RecoClusters) pre_RecoClusters->clear(); + + std::string name = "ClusterSearcher::Execute()"; + Log(name + ": Executing",v_error,verbosity); + + // print Selecting parameters + if(verbosity>v_message) this->PrintParameters(); + + // see if "ANNIEEvent" exists + auto get_annieevent = m_data->Stores.count("ANNIEEvent"); + if(!get_annieevent){ + Log(name + ": No ANNIEEvent store!",v_error,verbosity); + return false; + }; + + /// see if "RecoEvent" exists + auto get_recoevent = m_data->Stores.count("RecoEvent"); + if(!get_recoevent){ + Log(name + ": No RecoEvent store!",v_error,verbosity); + return false; + }; + + if (fisMC) { + // get true vertex + auto get_truevtx = m_data->Stores.at("RecoEvent")->Get("TrueVertex", fTrueVertex); + if(!get_truevtx){ + Log(name + ": Error retrieving TrueVertex from RecoEvent!",v_error,verbosity); + return false; + } + } + + // get digit list + auto get_recodigit = m_data->Stores.at("RecoEvent")->Get("RecoDigit",fDigitList); ///> Get digits from "RecoEvent" + if(!get_recodigit){ + Log("VtxSeedGenerator Tool: Error retrieving RecoDigits,no digit from the RecoEvent store!",v_error,verbosity); + return false; + } + // copy to a new digit list + std::vector* digits = new std::vector; + for(int i=0;isize());i++) { + RecoDigit* recodigitptr = &(fDigitList->at(i)); + digits->push_back((RecoDigit*)recodigitptr); + } + + // Reset Digit Filter + // ============ + for(int n=0; nsize()); n++ ) { + digits->at(n)->ResetFilter(); + } + + // Run Clustering + // ================ + //std::vector* SelectDigitList = Run(digits); // Cluster digits for muon candidates + + // Set Digit Filter + // ========== + //for(int n=0; nsize()); n++ ) { + // RecoDigit* SelectDigit = (RecoDigit*)(SelectDigitList->at(n)); + // SelectDigit->PassFilter(); + //} + SelectDigits(digits); + + fRecoClusters = this->RecoClusters(digits); + // Digit clustering done! + // ===== + + /*RecoCluster* dummyCluster1 = new RecoCluster(); + Position pos(10,10,10); + RecoDigit* dummyDigit1 = new RecoDigit(0,pos,1,10,10,5); + dummyCluster1->AddDigit(dummyDigit1); + pos.SetX(20); + RecoDigit* dummyDigit2 = new RecoDigit(0,pos,1,11,10,7); + dummyCluster1->AddDigit(dummyDigit2); + dummyCluster1->CalcParameters(); + + fRecoClusters->push_back(dummyCluster1); + Log("dummyCluster1 loaded in. Size of fRecoCLusters: "+to_string(fRecoClusters->size()),v_debug,verbosity); + + RecoCluster* dummyCluster2 = new RecoCluster(); + Position pos2(20, 20, 10); + RecoDigit* dummyDigit3 = new RecoDigit(0, pos2, 1, 10, 10, 9); + dummyCluster2->AddDigit(dummyDigit3); + pos2.SetX(30); + RecoDigit* dummyDigit4 = new RecoDigit(0, pos2, 1, 11, 10, 11); + dummyCluster2->AddDigit(dummyDigit4); + dummyCluster2->CalcParameters(); + + fRecoClusters->push_back(dummyCluster2); + Log("dummyCluster2 loaded in. Size of fRecoClusters: "+to_string(fRecoClusters->size()),v_debug,verbosity); + */ + + pre_RecoClusters = nullptr; // to read existing clusters + bool cluster_status = m_data->Stores.at("RecoEvent")->Get("RecoClusters", pre_RecoClusters); + + Log("ClusterSearcher Tool: Final count of clusters: "+to_string(fRecoClusters->size()),v_debug,verbosity); + if(!cluster_status) { + pre_RecoClusters=new std::vector; + for (int i = 0; i < fRecoClusters->size(); i++) { + pre_RecoClusters->push_back(fRecoClusters->at(i)); + } + fIsHitClusteringDone = true; + m_data->Stores.at("RecoEvent")->Set("HitClusteringDone", fIsHitClusteringDone); + m_data->Stores.at("RecoEvent")->Set("RecoClusters", pre_RecoClusters); + } + else { + pre_RecoClusters->insert(pre_RecoClusters->end(), fRecoClusters->cbegin(), fRecoClusters->cend()); // append new clusters to the existing cluster list + //m_data->Stores.at("RecoEvent")->Set("RecoClusters", pre_RecoClusters); // save updated clusters + //fRecoClusters->clear(); + } + + if (!fFirstRun && pre_RecoClusters->size() == 0) { + Log("ClusterSearcher: NO CLUSTERS FOUND IN EVENT!",v_debug,verbosity); + Log("Digit list size: "+to_string(fDigitList->size()),v_debug,verbosity); + } + + delete digits; digits = 0; + return true; +} + +bool ClusterSearcher::Finalise(){ + //delete fClusteringParam; fClusteringParam = 0; //Will be deleted by the store, don't manually delete + delete fSelectAll; fSelectAll = 0; + delete fSelectByPulseHeight; fSelectByPulseHeight = 0; + delete fSelectByNeighbours; fSelectByNeighbours = 0; + delete fSelectByClusters; fSelectByClusters = 0; + //delete fRecoClusters; fRecoClusters = 0; //Will be deleted by the store, don't manually delete + delete fClusterList; fClusterList = 0; + // for test + delete fSelectByTruthInfo; fSelectByTruthInfo = 0; + return true; +} + +void ClusterSearcher::Config(int config) +{ + ClusterSearcher::Instance()->SetConfig(config); +} + +void ClusterSearcher::PmtMinPulseHeight(double min) +{ + ClusterSearcher::Instance()->SetPmtMinPulseHeight(min); +} + +void ClusterSearcher::PmtNeighbourRadius(double radius) +{ + ClusterSearcher::Instance()->SetPmtNeighbourRadius(radius); +} + +void ClusterSearcher::PmtNeighbourDigits(int digits) +{ + ClusterSearcher::Instance()->SetPmtNeighbourDigits(digits); +} + +void ClusterSearcher::PmtClusterRadius(double radius) +{ + ClusterSearcher::Instance()->SetPmtClusterRadius(radius); +} + +void ClusterSearcher::MinClusterDigits(int digits) +{ + ClusterSearcher::Instance()->SetMinClusterDigits(digits); +} + +void ClusterSearcher::PmtTimeWindowN(double windowN) +{ + ClusterSearcher::Instance()->SetPmtTimeWindowNeighbours(windowN); +} + +void ClusterSearcher::PmtTimeWindowC(double windowC) +{ + ClusterSearcher::Instance()->SetPmtTimeWindowClusters(windowC); +} + +void ClusterSearcher::PrintParameters() +{ + std::cout << " *** ClusterSearcher::PrintParameters() *** " << std::endl; + + std::cout << " Clustering Parameters: " << std::endl + << " ClusterMode = " << fClusterMode<< std::endl + << " Config = " << fConfig<< std::endl + << " PmtMinPulseHeight = " << fPmtMinPulseHeight << std::endl + << " PmtNeighbourRadius = " << fPmtNeighbourRadius << std::endl + << " PmtMinNeighbourDigits = " << fPmtMinNeighbourDigits << std::endl + << " PmtClusterRadius = " << fPmtClusterRadius << std::endl + << " PmtTimeWindowN = " << fPmtTimeWindowN << std::endl + << " PmtTimeWindowC = " << fPmtTimeWindowC << std::endl + << " LappdMinPulseHeight = " << fLappdMinPulseHeight << std::endl + << " LappdNeighbourRadius = " << fLappdNeighbourRadius << std::endl + << " LappdMinNeighbourDigits = " << fLappdMinNeighbourDigits << std::endl + << " LappdClusterRadius = " << fLappdClusterRadius << std::endl + << " LappdTimeWindowN = " << fLappdTimeWindowN << std::endl + << " LappdTimeWindowC = " << fLappdTimeWindowC << std::endl + << " MinClusterDigits = " << fMinClusterDigits << std::endl; +} + +void ClusterSearcher::Reset() +{ + + return; +} + +std::vector* ClusterSearcher::Run(std::vector* DigitList) +{ + if(verbosity>v_debug) std::cout << " *** ClusterSearcher::Run(...) *** " << std::endl; + + // input digit list + // ================ + std::vector* InputList = DigitList; + std::vector* OutputList = DigitList; + + // Select all digits + // ================= + InputList = ResetDigits(OutputList); + OutputList = (std::vector*)(this->SelectAll(InputList)); + OutputList = SelectDigits(OutputList); + if( fConfig==ClusterSearcher::kNone ) return OutputList; + + // Select by pulse height + // ====================== + InputList = ResetDigits(OutputList); + OutputList = (std::vector*)(this->SelectByPulseHeight(InputList)); + OutputList = SelectDigits(OutputList); + if( fConfig==ClusterSearcher::kPulseHeight ) return OutputList; + + // Select using neighbouring digits + // ================================ + InputList = ResetDigits(OutputList); + OutputList = (std::vector*)(this->SelectByNeighbours(InputList)); + OutputList = SelectDigits(OutputList); + if( fConfig==ClusterSearcher::kPulseHeightAndNeighbours ) return OutputList; + + // Select using clustered digits + // ============================= + InputList = ResetDigits(OutputList); + OutputList = (std::vector*)(this->SelectByClusters(InputList)); + OutputList = SelectDigits(OutputList); + if( fConfig==ClusterSearcher::kPulseHeightAndClusters ) return OutputList; + + if (fisMC){ + // Select using truth information (for simulation test only) + // ============================= + InputList = ResetDigits(OutputList); + OutputList = (std::vector*)(this->SelectByTruthInfo(InputList)); + OutputList = SelectDigits(OutputList); + if( fConfig==ClusterSearcher::kPulseHeightAndTruthInfo ) return OutputList; + } + + // return vector of selected digits + // ================================ + return OutputList; +} + +// Set all filter status to 0 (default is 1) +std::vector* ClusterSearcher::ResetDigits(std::vector* DigitList) +{ + for(int idigit=0; idigitsize()); idigit++ ){ + RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit)); + recoDigit->ResetFilter(); + } + + return DigitList; +} + +std::vector* ClusterSearcher::SelectDigits(std::vector* DigitList) +{ + for(int idigit=0; idigitsize()); idigit++ ){ + RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit)); + recoDigit->PassFilter(); + } + + return DigitList; +} + +void ClusterSearcher::SelectDigits(std::vector* DigitList) +{ + for (int idigit = 0; idigitsize()); idigit++) { + RecoDigit recoDigit = (RecoDigit)(DigitList->at(idigit)); + recoDigit.PassFilter(); + } + + return; +} + +std::vector* ClusterSearcher::SelectAll(std::vector* DigitList) +{ + std::string name = "ClusterSearcher::SelectAll() "; + // clear vector of selected digits + // ============================== + fSelectAll->clear(); + // select all digits + // ================= + for(int idigit=0; idigitsize()); idigit++ ){ + RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit)); + fSelectAll->push_back(recoDigit); + } + + // return vector of selected digits + // ================================ + if(verbosity>v_message) std::cout << name << " select all: " << fSelectAll->size() << std::endl; + + return fSelectAll; +} + +std::vector* ClusterSearcher::SelectByPulseHeight(std::vector* DigitList) +{ + std::string name = "ClusterSearcher::SelectByPulseHeight() "; + // clear vector of selected digits + // =============================== + fSelectByPulseHeight->clear(); + + // Select by pulse height + // ====================== + for(int idigit=0; idigitsize()); idigit++ ){ + RecoDigit* recoDigit = (RecoDigit*)(DigitList->at(idigit)); + double qep = recoDigit->GetCalCharge(); + if (!fisMC){ + int pmtid = recoDigit->GetDetectorID(); + unsigned long chankey = pmt_tubeid_to_channelkey[pmtid]; + if (pmt_gains[chankey]>0) qep/=pmt_gains[chankey]; + } + int detType = recoDigit->GetDigitType(); + if(detType == RecoDigit::lappd_v0) { + if( qep>fLappdMinPulseHeight ){ + fSelectByPulseHeight->push_back(recoDigit); + } + } + else if(detType == RecoDigit::PMT8inch) { + if( qep>fPmtMinPulseHeight ){ + fSelectByPulseHeight->push_back(recoDigit); + } + } + else Log(name + " detector type doesn't exist!",v_error,verbosity); + } + + + // return vector of selected digits + // ================================ + if(verbosity>v_message) std::cout <size() << std::endl; + + return fSelectByPulseHeight; +} + +std::vector* ClusterSearcher::SelectByNeighbours(std::vector* DigitList) +{ + std::string name = "ClusterSearcher::SelectByNeighbours() "; + // clear vector of Selected digits + // =============================== + fSelectByNeighbours->clear(); + + // create array of neighbours + // ========================== + int Ndigits = DigitList->size(); + + if( Ndigits<=0 ){ + return fSelectByNeighbours; + } + + int* numNeighbours = new int[Ndigits]; + + for( int idigit=0; idigitsize()); idigit1++ ){ + RecoDigit* fdigit1 = (RecoDigit*)(DigitList->at(idigit1)); + TString digit1Type = fdigit1->GetDigitType(); + for(int idigit2=idigit1+1; idigit2size()); idigit2++ ){ + RecoDigit* fdigit2 = (RecoDigit*)(DigitList->at(idigit2)); + TString digit2Type = fdigit2->GetDigitType(); + + double dx = fdigit1->GetPosition().X() - fdigit2->GetPosition().X(); + double dy = fdigit1->GetPosition().Y() - fdigit2->GetPosition().Y(); + double dz = fdigit1->GetPosition().Z() - fdigit2->GetPosition().Z(); + double dt = fdigit1->GetCalTime() - fdigit2->GetCalTime(); + double drsq = dx*dx + dy*dy + dz*dz; + + if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::PMT8inch) { + if( drsq>0.0 + && drsq0.0 + && drsq0.0 + && drsq0.0 + && drsq0.0 + && drsq0.0 + && drsqsize()); idigit++){ + RecoDigit* fdigit = (RecoDigit*)(DigitList->at(idigit)); + //std::cout << "numNeighbours[" << idigit << "] = " << numNeighbours[idigit] << std::endl; + if( numNeighbours[idigit]>=fPmtMinNeighbourDigits ){ + fSelectByNeighbours->push_back(fdigit); + } + } + + // delete array of neighbours + // ========================== + delete [] numNeighbours; + + // return vector of Selected digits + // ================================ + if(verbosity>v_message) std::cout << name << " Select by neighbours: " << fSelectByNeighbours->size() << std::endl; + + + return fSelectByNeighbours; +} + +std::vector* ClusterSearcher::SelectByClusters(std::vector* DigitList) +{ + std::string name = "ClusterSearcher::SelectByClusters() "; + // clear vector of Selected digits + // =============================== + fSelectByClusters->clear(); + //fRecoClusters->clear(); + + // run clustering algorithm + // ======================== + std::vector* ClusterList = (std::vector*)(this->RecoClusters(DigitList)); + + for(int icluster=0; iclustersize()); icluster++ ){ + RecoCluster* Cluster = (RecoCluster*)(ClusterList->at(icluster)); + fRecoClusters->push_back(Cluster); + + for(int idigit=0; idigitGetNDigits(); idigit++ ){ + RecoDigit* Digit = (RecoDigit*)(Cluster->GetDigit(idigit)); + fSelectByClusters->push_back(Digit); + } + } + + // return vector of Selected digits + // ================================ + if(verbosity>v_message) std::cout << name <<" Select by clusters: " << fSelectByClusters->size() << std::endl; + + + return fSelectByClusters; +} + +std::vector* ClusterSearcher::RecoClusters(std::vector* DigitList) +{ + + // delete cluster digits + // ===================== + for(int i=0; iclear(); + + //Digit Selection + SelectByPulseHeight(DigitList); + SelectByNeighbours(fSelectByPulseHeight); + + // make cluster digits + // =================== + for(int idigit=0; idigitsize()); idigit++ ){ + RecoDigit* recoDigit = (RecoDigit*)(fSelectByNeighbours->at(idigit)); + RecoClusterDigit* clusterDigit = new RecoClusterDigit(recoDigit); + vClusterDigitList.push_back(clusterDigit); + } + + // run clustering algorithm + // ======================== + for(int idigit1=0; idigit1GetDigitType(); + for(int idigit2=idigit1+1; idigit2GetDigitType(); + + double dx = fdigit1->GetX() - fdigit2->GetX(); + double dy = fdigit1->GetY() - fdigit2->GetY(); + double dz = fdigit1->GetZ() - fdigit2->GetZ(); + double dt = fdigit1->GetTime() - fdigit2->GetTime(); + double drsq = dx*dx + dy*dy + dz*dz; + if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::PMT8inch) { + if( drsq>0.0 + && drsqAddClusterDigit(fdigit2); + fdigit2->AddClusterDigit(fdigit1); + } + if (drsq == 0 && fabs(dt) < fPmtTimeWindowC) { + fdigit2->SetClustered(); + } + } + + if(digit1Type == RecoDigit::lappd_v0 && digit2Type == RecoDigit::lappd_v0) { + if( drsq>0.0 + && drsqAddClusterDigit(fdigit2); + fdigit2->AddClusterDigit(fdigit1); + } + } + + if(digit1Type == RecoDigit::PMT8inch && digit2Type == RecoDigit::lappd_v0) { + if( drsq>0.0 + && drsqAddClusterDigit(fdigit2); + } + if( drsq>0.0 + && drsqAddClusterDigit(fdigit1); + } + } + + if(digit1Type == RecoDigit::lappd_v0 && digit2Type == RecoDigit::PMT8inch) { + if( drsq>0.0 + && drsqAddClusterDigit(fdigit2); + } + if( drsq>0.0 + && drsqAddClusterDigit(fdigit1); + } + } + + } + } + + // collect up clusters + // =================== + Bool_t carryon = 0; + for(int idigit=0; idigitIsClustered()==0 && fdigit->GetNClusterDigits()>0 ){ + Log("ClusterSearcher: valid fdigit index: "+to_string(idigit),v_debug,verbosity); + Log("with time: "+to_string(fdigit->GetTime()),v_debug,verbosity); + Log("And location: "+to_string(fdigit->GetX())+", "+to_string(fdigit->GetY())+", "+to_string(fdigit->GetZ()),v_debug,verbosity); + vClusterDigitCollection.clear(); + vClusterDigitCollection.push_back(fdigit); + fdigit->SetClustered(); + + carryon = 1; + while( carryon ){ + carryon = 0; + for(int jdigit=0; jdigit=fMinClusterDigits ){ + RecoCluster* cluster = new RecoCluster(); + cluster->SetClusterMode(fClusterMode); + fClusterList->push_back(cluster); + + for(int jdigit=0; jdigitGetRecoDigit()); + cluster->AddDigit(recodigit); + } + cluster->CalcParameters(); + Log("ClusterSearcher: Clusters made: "+to_string(fClusterList->size()),v_debug,verbosity); + } + } + } + + std::cout <<"Number of clusters = "<size()<* ClusterSearcher::SelectByTruthInfo(std::vector* DigitList) +{ + std::string name = " ClusterSearcher::SelectByTruthInfo(() "; + // clear vector of Selected digits + // =============================== + fSelectByTruthInfo->clear(); + + //===================== true info. + double x0 = fTrueVertex->GetPosition().X(); + double y0 = fTrueVertex->GetPosition().Y(); + double z0 = fTrueVertex->GetPosition().Z(); + double dirx = fTrueVertex->GetDirection().X(); + double diry = fTrueVertex->GetDirection().Y(); + double dirz = fTrueVertex->GetDirection().Z(); + double t0 = 0.0; + + // Select by truth + // ====================== + for(int idigit=0; idigitsize()); idigit++ ){ + RecoDigit* recoDigit = (RecoDigit*)DigitList->at(idigit); + double x = recoDigit->GetPosition().X(); + double y = recoDigit->GetPosition().Y(); + double z = recoDigit->GetPosition().Z(); + double dx = x-x0; + double dy = y-y0; + double dz = z-z0; + double ds = sqrt(dx*dx + dy*dy + dz*dz); + double px = dx/ds; + double py = dy/ds; + double pz = dz/ds; + double cosphi = px*dirx + py*diry + pz*dirz; + double phideg = acos(cosphi)/3.14*180; + if(phideg>38) fSelectByTruthInfo->push_back(recoDigit); + } + + // return vector of Selected digits + // ================================ + if(verbosity>v_message) std::cout << name << " Select by opening angle: " << fSelectByTruthInfo->size() << std::endl; + + return fSelectByTruthInfo; +} + diff --git a/UserTools/ClusterSearcher/ClusterSearcher.h b/UserTools/ClusterSearcher/ClusterSearcher.h new file mode 100644 index 000000000..1fcce0cec --- /dev/null +++ b/UserTools/ClusterSearcher/ClusterSearcher.h @@ -0,0 +1,157 @@ +#ifndef CLUSTERSEARCHER_H +#define CLUSTERSEARCHER_H + +#include +#include +#include + +#include "Tool.h" +#include "RecoCluster.h" +#include "RecoClusterDigit.h" +#include "TString.h" + +class ClusterSearcher: public Tool { + + public: + + ClusterSearcher(); + ~ClusterSearcher(); + bool Initialise(std::string configfile,DataModel &data); + bool Execute(); + bool Finalise(); + + typedef enum EFilterConfig { + kNone = 0, + kPulseHeight = 1, + kPulseHeightAndNeighbours = 2, + kPulseHeightAndClusters = 3, + kPulseHeightAndTruthInfo = 4 + } FilterConfig_t; + + static ClusterSearcher* Instance(); + + static void Config(int config); + static void ClusterType(int ctype); + static void PmtMinPulseHeight(double min); + static void PmtNeighbourRadius(double radius); + static void PmtNeighbourDigits(int digits); + static void PmtClusterRadius(double radius); + static void MinClusterDigits(int digits); + static void PmtTimeWindowN(double windowN); + static void PmtTimeWindowC(double windowC); + + void PrintParameters(); + + void SetConfig(int config) { fConfig = config; } + void SetClusterMode(int cmode) {fClusterMode = cmode;} + void SetPmtMinPulseHeight(double min) { fPmtMinPulseHeight = min; } + void SetPmtNeighbourRadius(double radius) { fPmtNeighbourRadius = radius; } + void SetPmtNeighbourDigits(int digits) { fPmtMinNeighbourDigits = digits; } + void SetPmtClusterRadius(double radius) { fPmtClusterRadius = radius; } + void SetPmtTimeWindowNeighbours(double windowN) { fPmtTimeWindowN = windowN; } + void SetPmtTimeWindowClusters(double windowC) { fPmtTimeWindowC = windowC; } + + void SetLappdMinPulseHeight(double min) { fLappdMinPulseHeight = min; } + void SetLappdNeighbourRadius(double radius) { fLappdNeighbourRadius = radius; } + void SetLappdNeighbourDigits(int digits) { fLappdMinNeighbourDigits = digits; } + void SetLappdClusterRadius(double radius) { fLappdClusterRadius = radius; } + void SetLappdTimeWindowNeighbours(double windowN) { fLappdTimeWindowN = windowN; } + void SetLappdTimeWindowClusters(double windowC) { fLappdTimeWindowC = windowC; } + void LoadConfigFile(string configfilename); + void SetMinClusterDigits(int digits) { fMinClusterDigits = digits; } + + + std::vector* Run(std::vector* digitlist); + std::vector* ResetDigits(std::vector* digitlist); + std::vector* SelectDigits(std::vector* digitlist); + void SelectDigits(std::vector* digitlist); + std::vector* SelectAll(std::vector* digitlist); + std::vector* SelectByPulseHeight(std::vector* digitlist); + std::vector* SelectByNeighbours(std::vector* digitlist); + std::vector* SelectByClusters(std::vector* digitlist); + std::vector* SelectByTruthInfo(std::vector* digitlist); //use truth information. Only for testing the code + std::vector* RecoClusters(std::vector* digitlist); + + + + private: + void Reset(); + + // running mode + int fConfig; + + // clustering mode + int fClusterMode; + + // cleaning parameters + double fPmtMinPulseHeight; + double fPmtNeighbourRadius; + int fPmtMinNeighbourDigits; + double fPmtClusterRadius; + double fPmtTimeWindowN; + double fPmtTimeWindowC; + double fPmtMinHitsPerCluster; + + double fLappdMinPulseHeight; + double fLappdNeighbourRadius; + int fLappdMinNeighbourDigits; + double fLappdClusterRadius; + + double fLappdTimeWindowN; + double fLappdTimeWindowC; + double fLappdMinHitsPerCluster; + + int fMinClusterDigits; + bool fisMC; + bool fFirstRun; //Likely temp; would be nice to have it automatically set. + + // p.e. conversion parameters + std::map pmt_tubeid_to_channelkey; + std::map pmt_gains; + std::string singlePEgains; + + // Container for parameters + std::map* fClusteringParam = nullptr; + + // internal containers + std::vector vNdigitsCluster; + std::vector vClusterDigitList; + std::vector vClusterDigitCollection; + + // vectors of selected digitss + std::vector* fSelectAll; + std::vector* fSelectByPulseHeight; + std::vector* fSelectByNeighbours; + std::vector* fSelectByClusters; + + // for test only + std::vector* fSelectByTruthInfo; + + // vectors of clusters + std::vector* fClusterList; + + // vector of clusters (accessible to the CStore) + std::vector* fRecoClusters = nullptr; + std::vector* pre_RecoClusters = nullptr; + + // true vertex + RecoVertex* fTrueVertex = 0; + + // digit list + std::vector* fDigitList = 0; + + // hit clustering status; + bool fIsHitClusteringDone = false; + + /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged. + int verbosity=1; + int v_error=0; + int v_warning=1; + int v_message=2; + int v_debug=3; + std::string logmessage; + +}; + + +#endif diff --git a/UserTools/ClusterSearcher/README.md b/UserTools/ClusterSearcher/README.md new file mode 100644 index 000000000..1fd30a030 --- /dev/null +++ b/UserTools/ClusterSearcher/README.md @@ -0,0 +1,50 @@ +# ClusterSearcher + +ClusterSearcher +Creates list of RecoClusters according to provided time- and neighboring- parameters. Clusters are tagged with an int ClusterMode variable for use to differentiate different forms of clusters, created from multiple iterations of the tool within the ToolChain, using separate configfiles (e.g. ClusterSearcherConfig, ClusterSearcherConfig2, etc., using Clustermode 1, 2, etc. respectively) + +## Data + +Describe any data formats ClusterSearch creates, destroys, changes, or analyzes. E.G. + +**pmt_tubeid_to_chanelkey** from CStore +* Identifies individual PMTs in order to find effective charge in p.e. + +**ClusteringParameters** from RecoEvent +* Compares execution parameters to other iterations of this tool within the toolchain + +**TrueVertex** from RecoEvent +* True vertex position, used for + +**RecoClusters** Into and from RecoEvent +* List of RecoClusters found based on the configuration parameters provided +* Any RecoClusters created by previous iteration(s) of ClusterSearcher in the toolchain, so that new clusters can be added to the singular list + + + + +## Configuration + + +verbosity 5 +ClusterMode 1 #searching track-like clusters +Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +PmtMinPulseHeight 20 #minimum pulse height +PmtNeighbourRadius 60 #digit neighbouring distance [cm] +PmtMinNeighbourDigits 2 #minimum neighbour digits +PmtClusterRadius 60 #digit clustering distance [cm] +PmtTimeWindowN 10 #neighbouring time window [ns] +PmtTimeWindowC 100 #clustering time window [ns] +PmtMinHitsPerCluster 4 #number of hits per cluster +LappdMinPulseHeight 0 #minimum pulse height +LappdNeighbourRadius 25 #digit neighbouring distance [cm] +LappdMinNeighbourDigits 20 #minimum neighbour digits +LappdClusterRadius 25 #digit clustering distance [cm] +LappdTimeWindowN 1 #neighbouring time window [ns] +LappdTimeWindowC 1 #clustering time window [ns] +LappdMinHitsPerCluster 5 #number of digits per cluster +MinClusterDigits 50 #minimum clustered digits (LAPPD+PMT) +MinClusterDigits 4 #minimum clustered digits (PMT-only) +FirstRun 1 #Debug input for differentiating first iteration from later iterations + +``` diff --git a/UserTools/DigitBuilder/DigitBuilder.cpp b/UserTools/DigitBuilder/DigitBuilder.cpp index b3a5a3e4e..b74f6c2ec 100644 --- a/UserTools/DigitBuilder/DigitBuilder.cpp +++ b/UserTools/DigitBuilder/DigitBuilder.cpp @@ -18,7 +18,7 @@ DigitBuilder::~DigitBuilder() { bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ /////////////////// Usefull header /////////////////////// - if(verbosity) cout<<"Initializing Tool DigitBuilder"<; + fHitLAPPDs = new std::vector; // Make the RecoDigit Store if it doesn't exist int recoeventexists = m_data->Stores.count("RecoEvent"); @@ -57,6 +63,7 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ Log("DigitBuilder Tool: Error retrieving Geometry from ANNIEEvent!",v_error,verbosity); return false; } + Log("Strip Hit Mode:" + to_string(striphit),v_debug,verbosity); // Some hard-coded values of old WCSim LAPPDIDs are in this Tool // I would recommend moving away from the use of WCSim IDs if possible as they are liable to change @@ -66,6 +73,10 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ m_data->CStore.Get("channelkey_to_pmtid",channelkey_to_pmtid); } else { ifstream file_pmtid(path_chankeymap.c_str()); + if (!file_pmtid) { + Log("DigitBuilder Tool: Did not find chankeymap file",v_error,verbosity); + return false; + } while (!file_pmtid.eof()){ unsigned long chankey; int pmtid; @@ -73,19 +84,44 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ channelkey_to_pmtid.emplace(chankey,pmtid); pmtid_to_channelkey.emplace(pmtid,chankey); if (file_pmtid.eof()) break; + Log("DigitBuilder Tool: still gathering chankeys",v_debug,verbosity); } - + Log("DigitBuilder Tool: chankeys done",v_debug,verbosity); file_pmtid.close(); m_data->CStore.Set("pmt_tubeid_to_channelkey",pmtid_to_channelkey); ifstream file_singlepe(singlePEgains.c_str()); + if (!file_singlepe) { + Log("DigitBuilder Tool: Did not find SinglePEgains file",v_error,verbosity); + return false; + } unsigned long temp_chankey; double temp_gain; - while (!file_singlepe.eof()){ + + std::string line; + if (file_singlepe.is_open()) { + //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; + std::vector DataEntries; + boost::split(DataEntries, line, boost::is_any_of(","), boost::token_compress_on); + int channelkey = -9999; + double SPECharge = -9999.; + channelkey = std::stoi(DataEntries.at(0)); + SPECharge = std::stod(DataEntries.at(1)); + pmt_gains.emplace(channelkey, SPECharge); + } + } + + + /*while (!file_singlepe.eof()) { file_singlepe >> temp_chankey >> temp_gain; if (file_singlepe.eof()) break; pmt_gains.emplace(temp_chankey,temp_gain); - } + Log("DigitBuilder Tool: still collecting SPE gains: "+to_string(temp_gain), v_debug, verbosity); + }*/ + Log("DigitBuilder Tool: SPE gains done",v_debug,verbosity); file_singlepe.close(); } @@ -93,10 +129,10 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ //Read the LAPPDID file, if given if(fLAPPDIDFile!="none"){ - if(verbosity>2) std::cout << "Loading digits from LAPPD IDs in file " << fLAPPDIDFile << std::endl; + Log("Loading digits from LAPPD IDs in file " + fLAPPDIDFile,v_debug,verbosity); this->ReadLAPPDIDFile(); } else { - if(verbosity>2) std::cout << "Loading digits from all LAPPDs" << std::endl; + Log("Loading digits from all LAPPDs",v_debug,verbosity); } return true; } @@ -135,7 +171,12 @@ bool DigitBuilder::Execute(){ return false; } } else { - auto get_clusters = m_data->CStore.Get("ClusterMap",m_all_clusters); + auto get_dhits = m_data->Stores.at("ANNIEEvent")->Get("Hits",Hits); + if (!get_dhits) { + Log("DigitBuilder Tool: ERROR retrieving hits in Data mode!",v_error,verbosity); + return false; + } + /*auto get_clusters = m_data->CStore.Get("ClusterMap", m_all_clusters); if (!get_clusters){ Log("DigitBuilder Tool: ERROR retrieving clustered hits (ClusterMap) in Data mode!",v_error,verbosity); return false; @@ -144,7 +185,7 @@ bool DigitBuilder::Execute(){ if (!get_clusters_chankey){ Log("DigitBuilder Tool: ERROR retrieving clustered chankeys (ClusterMapDetkey) in Data mode!",v_error,verbosity); return false; - } + }*/ } /// Build RecoDigit @@ -161,7 +202,7 @@ bool DigitBuilder::Execute(){ bool DigitBuilder::Finalise(){ //delete fDigitList; fDigitList = 0; //Don't delete pointer to fDigitList, will be deleted by the BoostStore! - if(verbosity>0) cout<<"DigitBuilder exitting"<> + + if(fMCPMTHits){ + if (fCollectHits) { + Log("DigitBuilder Tool: Collecting Hits to make Digits",v_message,verbosity); + CollectMCPMTHits(); + Log("DigitBuilder Tool: # of digits in this event: "+to_string(fDigitList->size()),v_debug,verbosity); + } + else { Log("DigitBuilder Tool: Num PMT Digits = "+to_string(fMCPMTHits->size()),v_message, verbosity); /// iterate over the map of sensors with a measurement for(std::pair>&& apair : *fMCPMTHits){ @@ -243,69 +293,88 @@ bool DigitBuilder::BuildMCPMTRecoDigit() { if(det->GetDetectorElement()=="Tank"){ std::vector& hits = apair.second; if(fParametricModel){ - if(verbosity>2) std::cout << "Using parametric model to build PMT hits" << std::endl; + Log("Using parametric model to build PMT hits",v_debug,verbosity); //We'll get all hit info and then define a time/charge for each digit std::vector hitTimes; std::vector hitCharges; + //std::vector hitIDs; for(MCHit& ahit : hits){ - if(verbosity>3){ - std::cout << "This HIT'S TIME AND CHARGE: " << ahit.GetTime() << - "," << ahit.GetCharge() << std::endl; - } + 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) { hitTimes.push_back(ahit.GetTime()*1.0); hitCharges.push_back(ahit.GetCharge()); + //hitIDs.push_back(ahit.GetHitID()); } } // Do median and sum std::sort(hitTimes.begin(), hitTimes.end()); size_t timesize = hitTimes.size(); if (timesize == 0) continue; - if (timesize % 2 == 0){ - calT = (hitTimes.at(timesize/2 - 1) + hitTimes.at(timesize/2))/2; - } else { - calT = hitTimes.at(timesize/2); + if (fParametricModel == 1) { //Mean hit time + if (timesize % 2 == 0) { + calT = (hitTimes.at(timesize / 2 - 1) + hitTimes.at(timesize / 2)) / 2; + } + else { + calT = hitTimes.at(timesize / 2); + } + } + else if (fParametricModel == 2) { //first hit time + calT = hitTimes.at(0); } + else if (fParametricModel == 3) { //Average of first 20% of hit times + for (int i = 0; i < timesize / 5; i++) { + calT += hitTimes.at(i); + } + if ((int)(timesize / 5) > 0)calT = calT / ((int)(timesize / 5)); + } + else if (fParametricModel == 4) { //Average hit time + for (int i = 0; i < timesize; i++) { + calT += hitTimes.at(i); + } + calT = calT / timesize; + } + else if (fParametricModel == 5) { + for (int i = 0; i < timesize; i++) { + if (hitCharges.at(i) > maxT) maxT = hitTimes.at(i); + } + calT = maxT; + } + if (MCPMTResolution > 0) calT = frand.Gaus(calT, MCPMTResolution); calQ = 0.; for(std::vector::iterator it = hitCharges.begin(); it != hitCharges.end(); ++it){ calQ += *it; } - if (verbosity>4) { - std::cout << "PMT position (XfDigitChargeThr) { //changed to 0 for cross-checks with other tools, change back later! digitType = RecoDigit::PMT8inch; RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType, PMTId); + //recoDigit.SetHitIDs(hitIDs); fDigitList->push_back(recoDigit); } - } else { + }else{ for(MCHit& ahit : hits){ //if(v_message4) { - std::cout << "PMT position (Xpush_back(recoDigit); } - } + } } } // end loop over MCHits - } else { - cout<<"No MCHits"<size() > 0) { + fHitLAPPDs->clear(); + } + // repeat for LAPPD hits // MCLAPPDHits is a std::map> - if(fMCLAPPDHits){ + if(fMCLAPPDHits && fMCLAPPDHits->size() > 0){ Log("DigitBuilder Tool: Num LAPPD Digits = "+to_string(fMCLAPPDHits->size()),v_message,verbosity); // iterate over the map of sensors with a measurement for(std::pair>&& apair : *fMCLAPPDHits){ @@ -341,46 +415,122 @@ bool DigitBuilder::BuildMCLAPPDRecoDigit() { } if(!isSelectedLAPPD && fLAPPDId.size()>0) continue; if(verbosity>2){ - std::cout << "Loading in digits for LAPPDID " << LAPPDId << std::endl; + Log("Loading in digits for LAPPDID " + to_string(LAPPDId), v_message,verbosity); + Log("located: ",v_message,verbosity); + det->GetPositionInTank().Print(); + Log("directed: ",v_message,verbosity); + det->GetDetectorDirection().Print(); } - if(det->GetDetectorElement()=="LAPPD"){ // redundant, MCLAPPDHits are LAPPD hitss - std::vector& hits = apair.second; - for(MCLAPPDHit& ahit : hits){ - //if(v_message4) { - std::cout << "LAPPD position (XGetDetectorElement() == "LAPPD") { // redundant, MCLAPPDHits are LAPPD hitss + std::vector& hits = apair.second; + + std::vector sumposX; + std::vector sumposY; + std::vector sumposZ; + std::vector sumT; + double posX; + std::vector nHitsOnStrip; + int a; + for (int i = 0; i < 28; i++) { + sumposX.push_back(0); + sumposY.push_back(0); + sumposZ.push_back(0); + sumT.push_back(0); + nHitsOnStrip.push_back(0); } - // I found the charge is 0 for all the hits. In order to test the code, - // here I just set the charge to 1. We should come back to this later. (Jingbo Wang) - calQ = 1.; - digitType = RecoDigit::lappd_v0; - RecoDigit recoDigit(region, pos_reco, calT, calQ, digitType,LAPPDId); - //if(v_message40 || calT<-10) continue; // cut off delayed hits - fDigitList->push_back(recoDigit); + for(MCLAPPDHit& ahit : hits){ + if (LAPPDId != currentLAPPD && hits.size() > 3) { + Log("THERE ARE HITS ON LAPPD " + to_string(LAPPDId),v_message,verbosity); + currentLAPPD = LAPPDId; + fHitLAPPDs->push_back(currentLAPPD); + Log("fHitLAPPDs size: " + to_string(fHitLAPPDs->size()),v_message,verbosity); + } + if (striphit == 0) { + //if(v_message 40 || calT < -10) continue; // cut off delayed hits + //recoDigit.SetHitIDs(hitIDs); + fDigitList->push_back(recoDigit); + + } + else { + posX = (ahit.GetPosition().at(0) * 100 + xshift) - (det->GetPositionInTank().X() * 100 + xshift); + double dirX = det->GetDetectorDirection().X(); + double AngX = std::asin(dirX); + int strip = (int)((posX / (20.* std::sin(AngX) / 28.)) + 14 - 1); + if (strip > 27 || strip < 0) { + Log("warning: LAPPD hit found off LAPPD: ", v_warning, verbosity); + Log("hit found on strip " + to_string(strip),v_warning,verbosity); + Log("LAPPD truehit X,Y,Z: " + to_string(ahit.GetPosition().at(0)) + ", " + to_string(ahit.GetPosition().at(1)) + ", " + to_string(ahit.GetPosition().at(2)),v_warning,verbosity); + Log("continuing loop without this hit.",v_warning,verbosity); + continue; + } + Log("hit found on strip " + to_string(strip),v_debug,verbosity); + sumposX.at(strip) += ahit.GetPosition().at(0) * 100 + xshift; + sumposY.at(strip) += ahit.GetPosition().at(1) * 100 + yshift; + sumposZ.at(strip) += ahit.GetPosition().at(2) * 100 + zshift; + Log("LAPPD truehit X,Y,Z: " + to_string(ahit.GetPosition().at(0)) + ", " + to_string(ahit.GetPosition().at(1)) + ", " + to_string(ahit.GetPosition().at(2)),v_debug,verbosity); + + if ((striphit == 3 && nHitsOnStrip.at(strip) < 3) || (striphit==2 && nHitsOnStrip.at(strip) == 0) || striphit == 1) + sumT.at(strip) += frand.Gaus(calT, 0.1); // time is smeared with 100 ps time resolution. Harded-coded for now. + nHitsOnStrip.at(strip) += 1; + + } } + if (striphit != 0) { + for (int strip = 0; strip < 28; strip++) { + if (nHitsOnStrip.at(strip) != 0) { + pos_reco.SetX(sumposX.at(strip) / nHitsOnStrip.at(strip)); + pos_reco.SetY(sumposY.at(strip) / nHitsOnStrip.at(strip)); + pos_reco.SetZ(sumposZ.at(strip) / nHitsOnStrip.at(strip)); + if(striphit == 1) calT = sumT.at(strip) / nHitsOnStrip.at(strip); + if (striphit == 3) calT = sumT.at(strip) / 3; + calQ = nHitsOnStrip.at(strip); //set to the number of hits for now. + digitType = RecoDigit::lappd_v0; + + Log("LAPPD ID: " + to_string(LAPPDId),v_message,verbosity); + Log("LAPPD strip-hit position (Xpush_back(recoDigit); + } + } + } + sumposX.clear(); + sumposY.clear(); + sumposZ.clear(); + sumT.clear(); + nHitsOnStrip.clear(); } } // end loop over MCLAPPDHits } else { - cout<<"No MCLAPPDHits"<>&& apair : *Hits) { + unsigned long chankey = apair.first; + Detector* thistube = fGeometry->ChannelToDetector(chankey); + det = fGeometry->ChannelToDetector(chankey); + int detectorkey = thistube->GetDetectorID(); + int PMTId = channelkey_to_pmtid.at(chankey); + if (thistube->GetDetectorElement() == "Tank") { + std::vector& ThisPMTHits = apair.second; + for (Hit& ahit : ThisPMTHits) { + pos_reco=det->GetDetectorPosition(); + calT = ahit.GetTime(); + + calQ = ahit.GetCharge(); + + if (pmt_gains.find(chankey) != pmt_gains.end() && pmt_gains.at(chankey) > 0.0) { + calQ_temp = calQ / pmt_gains.at(chankey); + } + if (calQ_temp > fDigitChargeThr) { + digitType = RecoDigit::PMT8inch; + RecoDigit recoDigit(region, pos_reco, calT, calQ_temp, digitType, PMTId); + + //recoDigit.SetHitIDs(hitIDs); + fDigitList->push_back(recoDigit); + } + + } + + } + } + + + + + + /* + /// m_all_clusters is a std::map> if (m_all_clusters && m_all_clusters_detkey){ int clustersize = m_all_clusters->size(); - std::cout <<"Clustersize of m_all_clusters: "<::iterator it3 = hitcharges.begin(); it3 != hitcharges.end(); ++it3){ calQ += *it3; } - if (verbosity>4) { - std::cout << "PMT position (X 0.0){ calQ_temp = calQ / pmt_gains.at(chankey); @@ -500,6 +685,8 @@ bool DigitBuilder::BuildDataPMTRecoDigit(){ if(calQ_temp>fDigitChargeThr) { digitType = RecoDigit::PMT8inch; RecoDigit recoDigit(region, pos_reco, calT, calQ_temp, digitType, PMTId); + + //recoDigit.SetHitIDs(hitIDs); fDigitList->push_back(recoDigit); } } @@ -524,15 +711,13 @@ bool DigitBuilder::BuildDataPMTRecoDigit(){ for (int i=0; i< int(hitcharges.size()); i++){ calT = hittimes.at(i); calQ = hitcharges.at(i); - if (verbosity>4) { - std::cout << "PMT position (Xpush_back(recoDigit); } } @@ -540,9 +725,9 @@ bool DigitBuilder::BuildDataPMTRecoDigit(){ } } } else { - cout<<"No Clustered Hits found."<Stores.at("RecoEvent")->Set("RecoDigit", fDigitList, savetodisk); ///> Add digits to RecoEvent + m_data->Stores.at("RecoEvent")->Set("HitLAPPDs", fHitLAPPDs, savetodisk); } void DigitBuilder::Reset() { // Reset fDigitList->clear(); + fHitLAPPDs->clear(); } void DigitBuilder::ReadLAPPDIDFile() { @@ -564,7 +751,7 @@ void DigitBuilder::ReadLAPPDIDFile() { if (myfile.is_open()){ while(getline(myfile,line)){ if(verbosity>0){ - std::cout << "DigitBuilder tool: Loading hits from LAPPD ID " << line << std::endl; + Log("DigitBuilder tool: Loading hits from LAPPD ID " + line,v_message,verbosity); } int thisID = std::atoi(line.c_str()); fLAPPDId.push_back(thisID); @@ -573,3 +760,156 @@ void DigitBuilder::ReadLAPPDIDFile() { Log("Unable to open given LAPPD ID File. Using all LAPPDs",v_error,verbosity); } } + +void DigitBuilder::CollectMCPMTHits() { + double calT = 0; + double calQ = 0; + Position pos_reco, pos_sim; + std::map PMT_ishit; + + + for (std::pair>&& apair : *fMCPMTHits) { + unsigned long chankey = apair.first; + Detector* thistube = fGeometry->ChannelToDetector(chankey); + int detectorkey = thistube->GetDetectorID(); + int PMTId = channelkey_to_pmtid.at(chankey); + Log("DigitBuilder Tool: collecting hits for PMT "+to_string(detectorkey)+", which corresponds to PMTId "+to_string(PMTId),v_debug,verbosity); + if (thistube->GetDetectorElement() == "Tank") { + std::vector& hits = apair.second; + PMT_ishit[detectorkey] = 1; + std::vector hit_times; + std::vector temp_times; + double first_time=0; + double mid_time=0; + std::vector datalike_hits; + std::vector datalike_hits_charge; + std::vector Parents_by_hit; + std::vector temp_parents; + double max_charge=-999; + pos_sim = thistube->GetDetectorPosition(); + pos_sim.UnitToCentimeter(); + pos_reco.SetX(pos_sim.X() + xshift); + pos_reco.SetY(pos_sim.Y() + yshift); + pos_reco.SetZ(pos_sim.Z() + zshift); + + + for (MCHit& ahit : hits) { + //if (ahit.GetTime() < end_of_window_time_cut * AcqTimeWindow) { + hit_times.push_back(ahit.GetTime()); + //} + } + /*for (MCHit& ahit : hits) { + //std::cout <<"Key: "< 2000.) std::cout <<"Found hit later than 2us! Hit time : "< combine multiple photons if they are within a 10ns range + //hit times can only be recorded with 2ns precision --> possible times are 0ns, 2ns, 4ns, ... + hits_2ns_res.push_back(2 * (int(ahit.GetTime()) / 2.) + (int(ahit.GetTime()) % 2)); + hits_2ns_res_charge.push_back(ahit.GetCharge()); + } + }*/ + + if (hit_times.size() == 0) { + Log("DigitBuilder tool: no hits in window.",v_message,verbosity); + return; + } + + //Combine multiple MC hits to one pulse + std::sort(hit_times.begin(), hit_times.end()); + for (int i_hit = 0; i_hit < (int)hit_times.size(); i_hit++) { + Log("check 0: hits, hit_times, datalike_hits, parents size: "+to_string(hits.size())+" "+to_string(hit_times.size())+" "+to_string(datalike_hits.size())+" " + to_string(hits.at(i_hit).GetParents()->size()), v_debug, verbosity); + double hit1 = hit_times.at(i_hit); + if(hit1>end_of_window_time_cut * AcqTimeWindow) continue; + int j_hit=0; + if (datalike_hits.size() == 0) { + Log("check 0-b", v_debug, verbosity); + + datalike_hits.push_back(hit1); + first_time=hit1; + + datalike_hits_charge.push_back(hits.at(i_hit).GetCharge()); + Log("check 0-c", v_debug, verbosity); + if (hits.at(i_hit).GetParents()->size() == 0) { + Parents_by_hit.push_back(-5); + Log("found hit with no parent at time "+to_string(hits.at(i_hit).GetTime()),v_debug,verbosity); + } + else { + Parents_by_hit.push_back(hits.at(i_hit).GetParents()->at(0)); + Log("check 0-d: hit parent " + to_string(hits.at(i_hit).GetParents()->at(0)) + " and time: " + to_string(hits.at(i_hit).GetTime()), v_debug, verbosity); + } + + temp_times.push_back(hit1); + } + else { + bool new_pulse = false; + Log("check 0-a", v_debug, verbosity); + for (int j_hit = 0; j_hit < (int)datalike_hits.size(); j_hit++) { + + if (fabs(first_time - hit1) < 10.) { + new_pulse = false; + datalike_hits_charge.at(j_hit) += hits.at(i_hit).GetCharge(); + temp_times.push_back(hit1); + + } + else { + new_pulse = true; + temp_times.push_back(hit1); + + } + Log("check 1-"+to_string(j_hit)+" new_pulse: "+to_string(new_pulse), v_debug, verbosity); + } + + if (new_pulse) { + first_time=hit1; + + // following the DigitBuilder tool --> take median photon hit time as the hit time of the "pulse" + sort(temp_times.begin(),temp_times.end()); + if (temp_times.size() % 2 == 0) { + mid_time = (temp_times.at(temp_times.size() / 2 - 1) + temp_times.at(temp_times.size() / 2)) / 2; + } + else { + mid_time = temp_times.at(temp_times.size() / 2); + } + Log("check 1-a", v_debug, verbosity); + datalike_hits.at(j_hit)=mid_time; //datalike_hits.push_back(mid_time); //Only count as a new pulse if it was 10ns away from every other pulse + datalike_hits_charge.push_back(hits.at(i_hit).GetCharge()); + j_hit++; + + datalike_hits.push_back(hit1); + + if (hits.at(i_hit).GetParents()->size() == 0) Parents_by_hit.push_back(-5); + else { + Parents_by_hit.push_back(hits.at(i_hit).GetParents()->at(0)); + Log("check 0-e: hit parent " + to_string(hits.at(i_hit).GetParents()->at(0)) + " and time: " + to_string(hits.at(i_hit).GetTime()), v_debug, verbosity); + } + temp_times.clear(); + + } + } + } + + Log("check 2", v_debug, verbosity); + //hits.clear(); + + for (int i_hit = 0; i_hit < (int)datalike_hits.size(); i_hit++) { + Log("check 1", v_debug, verbosity); + calT = datalike_hits.at(i_hit); + Log("check 2", v_debug, verbosity); + if (MCPMTResolution > 0) calT= frand.Gaus(calT, MCPMTResolution); + calQ = datalike_hits_charge.at(i_hit); + Log("check 3",v_debug,verbosity); + RecoDigit recoDigit(-999 /*region*/, pos_reco, calT, calQ, RecoDigit::PMT8inch, PMTId); + temp_parents.push_back(Parents_by_hit.at(i_hit)); + Log("check 4",v_debug,verbosity); + recoDigit.SetParents(temp_parents); + fDigitList->push_back(recoDigit); + temp_parents.clear(); + + Log("PMT position (X class DigitBuilder: public Tool { @@ -43,6 +44,7 @@ class DigitBuilder: public Tool { /// It adds PMT hits to the RecoDigit list bool BuildMCPMTRecoDigit(); bool BuildDataPMTRecoDigit(); ///> Same as BuildMCPMTRecoDigit, but applicable for data. + void CollectMCPMTHits(); /// \brief Build LAPPD digits /// @@ -72,15 +74,21 @@ class DigitBuilder: public Tool { int verbosity=1; std::string fInputfile; unsigned long fNumEvents; - + std::vector* fHitLAPPDs; std::vector fLAPPDId; ///< selected LAPPDs std::string fPhotodetectorConfiguration; ///< "PMTs_Only", "LAPPDs_Only", "All_Detectors" - bool fParametricModel; ///< configures if PMTs hits for each event are accumulated into one hit per PMT + int fParametricModel; ///< configures how PMTs hits for each event are accumulated into one hit per PMT 0: they are not, 1: median hit time, 2: first hit time, 3: average first 20% of hits, 4: average all hits, 5: highest charge hit time bool fIsMC; ///< Configure whether to load from MCHits or Hits in boost store std::string fLAPPDIDFile="none"; double fDigitChargeThr; std::string path_chankeymap; std::string singlePEgains; + int striphit; //0 all LAPPD Hits as true; 1 average all times per strip; 2 use first time for each strip + double MCPMTResolution = 0; + + bool fCollectHits; + int AcqTimeWindow; + double end_of_window_time_cut; Geometry* fGeometry=nullptr; ///< ANNIE Geometry TRandom3 frand; ///< Random number generator @@ -104,10 +112,11 @@ class DigitBuilder: public Tool { /// Reconstructed information std::vector* fDigitList; ///< Reconstructed Hits including both LAPPD hits and PMT hits std::map>* fMCPMTHits=nullptr; ///< PMT hits + std::map>* Hits = nullptr; std::map>* fMCLAPPDHits=nullptr; ///< LAPPD hits std::map>* fTDCData=nullptr; ///< MRD & veto hits - std::map>* m_all_clusters=nullptr; ///< Clusters, from ClusterFinder tool - std::map>* m_all_clusters_detkey=nullptr; ///< Chankeys corresponding to clusters, from ClusterFinder tool + std::map>* m_all_clusters=nullptr; ///< Clusters, from ClusterFinder tool - deprecated and to be removed + std::map>* m_all_clusters_detkey=nullptr; ///< Chankeys corresponding to clusters, from ClusterFinder tool - deprecated and to be removed std::map pmt_gains; diff --git a/UserTools/Factory/Factory.cpp b/UserTools/Factory/Factory.cpp index 1c8e2d29c..facb35018 100644 --- a/UserTools/Factory/Factory.cpp +++ b/UserTools/Factory/Factory.cpp @@ -178,5 +178,7 @@ if (tool=="EBPMT") ret=new EBPMT; if (tool=="LAPPDLoadStore") ret=new LAPPDLoadStore; if (tool=="PMTWaveformSim") ret=new PMTWaveformSim; if (tool=="LAPPDWaveformDisplay") ret=new LAPPDWaveformDisplay; +if (tool=="ClusterSearcher") ret=new ClusterSearcher; +if (tool=="NeutronCheck") ret=new NeutronCheck; return ret; } diff --git a/UserTools/NeutronCheck/NeutronCheck.cpp b/UserTools/NeutronCheck/NeutronCheck.cpp new file mode 100644 index 000000000..50391f1af --- /dev/null +++ b/UserTools/NeutronCheck/NeutronCheck.cpp @@ -0,0 +1,491 @@ +#include "NeutronCheck.h" + +NeutronCheck::NeutronCheck():Tool(){} + + +bool NeutronCheck::Initialise(std::string configfile, DataModel &data){ + + /////////////////// Useful header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); // loading config file + //m_variables.Print(); + + m_variables.Get("verbosity",verbosity); + m_variables.Get("outfile",outfile); + m_variables.Get("UseClean",useClean); + m_variables.Get("FinderCompare",fFinderCompare); + m_variables.Get("ParticleInfo",fParticleInfo); + + outfile+=".root"; + + fOutput_tfile = new TFile(outfile.c_str(), "recreate"); + NeutCheckTree = new TTree("NeutCheckTree", "Neutron Check Tree"); + + NeutCheckTree->Branch("eventNumber", &fMCEventNum, "eventNumber/I"); + + NeutCheckTree->Branch("eventCutStatus", &fEventStatusFlagged,"eventStatusFlagged/I"); + NeutCheckTree->Branch("ClusterNumber",&fClusterNum); + NeutCheckTree->Branch("TrueEnu",&true_Enu,"true_Enu/D"); + NeutCheckTree->Branch("TrueEmu",&true_Emu,"true_Emu/D"); + NeutCheckTree->Branch("TrueQ2",&TrueQ2,"true_Q2/D"); + NeutCheckTree->Branch("ClusterCount", &fClusterCount,"clusterCount/I"); + NeutCheckTree->Branch("RecoClusters",&fRecoClusters,"RecoClusters"); + + NeutCheckTree->Branch("TrueNeutronMult",&fTrueNeutronMult,"trueNeutronMult/I"); + NeutCheckTree->Branch("TrueNeutronDelayed",&fTrueNeutronDelayed,"trueNeutronDel/I"); + NeutCheckTree->Branch("TrueNeutCapT",&fMCNeutCapTimes); + + + NeutCheckTree->Branch("ClusterNeutronMult",&fNeutronMult,"neutronMult/I"); + + NeutCheckTree->Branch("ClusterMode",&fClusterMode); + NeutCheckTree->Branch("ClusterPDG",&fClusterPDG); + NeutCheckTree->Branch("ClusterParentPDG",&fClusterParentPDG); + NeutCheckTree->Branch("ClusterParticleEnergy",&fClusterParticleEnergy); + NeutCheckTree->Branch("ClusterHits", &fClusterHits); + NeutCheckTree->Branch("ClusterTime",&fClusterTime); + NeutCheckTree->Branch("ClusterCharge",&fClusterCharge); + NeutCheckTree->Branch("ClusterPurity",&fClusterPurity); + NeutCheckTree->Branch("ClusterCB",&fClusterCB); + NeutCheckTree->Branch("ClusterAS0",&fClusterAS0); + NeutCheckTree->Branch("ClusterAS1",&fClusterAS1); + NeutCheckTree->Branch("ClusterAS2",&fClusterAS2); + //NeutCheckTree->Branch("ClusterSA",&fClusterSA); + NeutCheckTree->Branch("ClusterCVX",&fClusterCVX); + NeutCheckTree->Branch("ClusterCVY", &fClusterCVY); + NeutCheckTree->Branch("ClusterCVZ", &fClusterCVZ); + NeutCheckTree->Branch("ClusterCVR", &fClusterCVR); + NeutCheckTree->Branch("ClusterAMD", &fClusterAMD); + NeutCheckTree->Branch("ClusterCA", &fClusterCA); + + + if(fFinderCompare){ + NeutCheckTree->Branch("FinderClusterNumber",&fFinderClusterNum); + NeutCheckTree->Branch("FinderClusterCount",&fFinderClusterCount,"FinderCount/I"); + NeutCheckTree->Branch("FinderPDG",&fFinderPDG); + NeutCheckTree->Branch("FinderParentPDG",&fFinderParentPDG); + NeutCheckTree->Branch("FinderCharge",&fFinderCharge); + NeutCheckTree->Branch("FinderPurity",&fFinderPurity); + NeutCheckTree->Branch("FinderTime",&fFinderTime); + } + + if (fParticleInfo) { + //Log("Particle Information to be added.",v_debug,verbosity); + NeutCheckTree->Branch("ParticlePDG",&fParticlePDG); + NeutCheckTree->Branch("ParticleParent",&fParticleParent); + NeutCheckTree->Branch("ParticleStartEnergy",&fParticleStartEnergy); + NeutCheckTree->Branch("ParticleStartTime", &fParticleStartTime); + NeutCheckTree->Branch("TrueNeutCapX", &fMCNeutCapX); + NeutCheckTree->Branch("TrueNeutCapY", &fMCNeutCapY); + NeutCheckTree->Branch("TrueNeutCapZ", &fMCNeutCapZ); + } + + + + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + + return true; +} + + +bool NeutronCheck::Execute(){ + + Log("NeutronCheck Tool Executing:",v_message,verbosity); + + ResetVariables(); + + // MC entry number + m_data->Stores.at("ANNIEEvent")->Get("MCEventNum", fMCEventNum); + + // MC trigger number + m_data->Stores.at("ANNIEEvent")->Get("MCTriggernum", fMCTriggerNum); + + // ANNIE Event number + m_data->Stores.at("ANNIEEvent")->Get("EventNumber", fEventNumber); + + auto get_flags = m_data->Stores.at("RecoEvent")->Get("EventFlagged", fEventStatusFlagged); + + if ((fEventStatusFlagged) != 0 && useClean) { + // if (!fEventCutStatus){ + Log("NeutronCheck Tool: Event was flagged with one of the active cuts.", v_debug, verbosity); + return true; + } + + + + + + + //if(EventStore=="ANNIEEvent"){ + + GetClusterInformation(); + + m_data->Stores.at("ANNIEEvent")->Get("ClusterToBestParticleID", fClusterToBestParticleID); + m_data->Stores.at("ANNIEEvent")->Get("ClusterToBestParticlePDG", fClusterToBestParticlePDG); + m_data->Stores.at("ANNIEEvent")->Get("ClusterEfficiency", fClusterEfficiency); + m_data->Stores.at("ANNIEEvent")->Get("ClusterPurity", fClusterPurityMap); + m_data->Stores.at("ANNIEEvent")->Get("ClusterTotalCharge", fClusterTotalCharge); + //bool get_q2 = m_data->Stores["GenieInfo"]->Get("EventQ2", TrueQ2); + + //if(get_q2) Log("True Q2 from GENIE: "+to_string(TrueQ2),v_debug,verbosity); + //else Log("Error finding Q2 info from GENIE",v_debug,verbosity); + //m_data->CStore.Get("ClausterMapMC", m_all_clusters); + + int bestPartId; + for (int i=0;iat(cluster_times.at(i)); + fFinderParentPDG.push_back(fMCParticles->at(bestPartId).GetParentPdg()); + fFinderPDG.push_back(fClusterToBestParticlePDG->at(cluster_times.at(i))); + fFinderPurity.push_back(fClusterPurityMap->at(cluster_times.at(i))); + + } + fFinderClusterCount=fFinderClusterNum.size()+1; + + fTrueNeutronMult=0; + fTrueNeutronDelayed=0; + Log("NeutCheck Tool: Scanning "+to_string(fMCParticles->size())+" MCParticles",v_debug,verbosity); + for (int i = 0; i < fMCParticles->size(); i++) { + if(fMCParticles->at(i).GetPdgCode()==2112 && fMCParticles->at(i).GetParentPdg()==0) { + fTrueNeutronMult++; + + if(fMCParticles->at(i).GetStopTime()>10000) fTrueNeutronDelayed++; + } + else if(fMCParticles->at(i).GetPdgCode()==2112 && fMCParticles->at(i).GetParentPdg()!=0){ + Log("Found non-primary neutron at particle "+to_string(i)+" with parent PDG "+ to_string(fMCParticles->at(i).GetParentPdg()), v_debug, verbosity); + } + if(fParticleInfo){ + fParticlePDG.push_back(fMCParticles->at(i).GetPdgCode()); + fParticleParent.push_back(fMCParticles->at(i).GetParentPdg()); + fParticleStartEnergy.push_back(fMCParticles->at(i).GetStartEnergy()); + fParticleStartTime.push_back(fMCParticles->at(i).GetStartTime()); + if (fMCParticles->at(i).GetPdgCode() == 2112) { + fMCNeutCapTimes.push_back(fMCParticles->at(i).GetStopTime()); + fMCNeutCapX.push_back(fMCParticles->at(i).GetStopVertex().X()); + fMCNeutCapY.push_back(fMCParticles->at(i).GetStopVertex().Y()); + fMCNeutCapZ.push_back(fMCParticles->at(i).GetStopVertex().Z()); + } + } + } + + auto get_muonMCEnergy = m_data->Stores.at("RecoEvent")->Get("TrueMuonEnergy", true_Emu); + bool get_neutrino_energy = m_data->Stores["GenieInfo"]->Get("NeutrinoEnergy", true_Enu); + + RecoVertex* truevtx = 0; + auto get_muonMC = m_data->Stores.at("RecoEvent")->Get("TrueVertex", truevtx); + + //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)); + TrueQ2=2*true_Enu*(true_Emu-p*cos(theta))-pow(105.7,2); + Log("NeutCheck TrueQ2: "+to_string(TrueQ2), v_debug, verbosity); + + CSCheck(); + + //MapRecoCheck->Fill(fClusterToBestParticleID->size(), fRecoClusters->size()); + //} + + NeutCheckTree->Fill(); + return true; +} + + +bool NeutronCheck::Finalise(){ + + Log("NeutronCheck Tool: Finalizing",v_message,verbosity); + fOutput_tfile->cd(); + + + NeutCheckTree->Write(); + + + + fOutput_tfile->Close(); + + + return true; +} + + +bool NeutronCheck::GetClusterInformation() { + + bool return_value = true; + + bool get_cluster_idxN = m_data->Stores["RecoEvent"]->Get("ClusterIndicesNeutron", cluster_neutron); + if (!get_cluster_idxN) Log("NeutronCheck tool: No ClusterIndicesNeutron In RecoEvent!", v_error, verbosity); + bool get_cluster_tN = m_data->Stores["RecoEvent"]->Get("ClusterTimesNeutron", cluster_times_neutron); + if (!get_cluster_tN) Log("NeutronCheck tool: No ClusterTimesNeutron In RecoEvent!", v_error, verbosity); + bool get_cluster_qN = m_data->Stores["RecoEvent"]->Get("ClusterChargesNeutron", cluster_charges_neutron); + if (!get_cluster_qN) Log("NeutronCheck tool: No ClusterChargesNeutron In RecoEvent!", v_error, verbosity); + bool get_cluster_cbN = m_data->Stores["RecoEvent"]->Get("ClusterCBNeutron", cluster_cb_neutron); + if (!get_cluster_cbN) Log("NeutronCheck tool: No ClusterCBNeutron In RecoEvent!", v_error, verbosity); + bool get_cluster_t = m_data->Stores["RecoEvent"]->Get("ClusterTimes", cluster_times); + if (!get_cluster_t) Log("NeutronCheck tool: No ClusterTimes In RecoEvent!", v_error, verbosity); + bool get_cluster_q = m_data->Stores["RecoEvent"]->Get("ClusterCharges", cluster_charges); + if (!get_cluster_q) Log("NeutronCheck tool: No ClusterCharges In RecoEvent!", v_error, verbosity); + bool get_cluster_cb = m_data->Stores["RecoEvent"]->Get("ClusterCB", cluster_cb); + if (!get_cluster_cb) Log("NeutronCheck tool: No ClusterCB In RecoEvent!", v_error, verbosity); + + bool goodMCParticles = m_data->Stores.at("ANNIEEvent")->Get("MCParticles", fMCParticles); + if (!goodMCParticles) { + std::cerr << "BackTracker: no MCParticles in the ANNIEEvent!" << endl; + return false; + } + + return_value = (get_cluster_idxN && get_cluster_tN && get_cluster_qN && get_cluster_cbN && get_cluster_t && get_cluster_q && get_cluster_cb); + + + + reco_Clusters = cluster_cb.size(); + + + + return return_value; + +} + +double NeutronCheck::ASCheck(vector*& digits, int mode) { + //mode config. 0 is maximal distance; 1 is max distance from highest charge + + double max_angle = 0; + double max_angle2 = 0; + if (mode == 0||mode==2) { + + double angle; + Position i_position,j_position; + for (int i = 0; i < digits->size(); i++) { + if (useClean && !digits->at(i).GetFilterStatus())continue; + i_position = digits->at(i).GetPosition(); + for (int j = 0; j < digits->size(); j++) { + if (j == i)continue; + if (useClean && !digits->at(j).GetFilterStatus())continue; + j_position = digits->at(j).GetPosition(); + + angle = j_position.Angle(i_position); + if (angle > max_angle)max_angle = angle; + + } + } + cout<<"Neutcheck max_angle: "<size(); i++) { + if (digits->at(i).GetCalCharge() > max_charge) { + if(useClean && !digits->at(i).GetFilterStatus())continue; + max_charge = digits->at(i).GetCalCharge(); + max_index = i; + max_position = digits->at(i).GetPosition(); + } + } + + + for (int i = 0; i < digits->size(); i++) { + if (i == max_index)continue; + if (useClean && !digits->at(i).GetFilterStatus())continue; + i_position = digits->at(i).GetPosition(); + angle = i_position.Angle(max_position); + if (angle > max_angle2)max_angle2 = angle; + } + cout<<"Neutcheck max_angle2: "<0) ratio = max_angle2/max_angle; + if(mode==2)return ratio; + + return 0; +} + +void NeutronCheck::CSCheck() { + + std::map* fMCParticleIndexMap; + + Log("Cluster Check!",v_debug,verbosity); + bool goodMCParticleIndexMap = m_data->Stores.at("ANNIEEvent")->Get("TrackId_to_MCParticleIndex", fMCParticleIndexMap); + if (!goodMCParticleIndexMap) { + std::cerr << "NeutCheck: no TrackId_to_MCParticleIndex in the ANNIEEvent!" << endl; + return; + } + + + + bool cluster_status = m_data->Stores.at("RecoEvent")->Get("RecoClusters", fRecoClusters); + if (!cluster_status) { + Log("Neutcheck tool found no recoclusters.",v_debug,verbosity); + return; + } + + int cluster_size = fRecoClusters->size(); + if (cluster_size == 0) { + Log("Neutcheck tool found empty recoclusterlist.",v_debug,verbosity); + return; + } + Log("Neutcheck tool found this many clusters: "+to_string(cluster_size),v_debug,verbosity); + + int bestParent; + int bestParticleID=-5; + int bestPDG=-5; + double CVX,CVY,CVZ,CVR; + + for (int i = 0; i < fRecoClusters->size(); i++) { + fClusterNum.push_back(i); + fClusterCount++; + bestParent = fRecoClusters->at(i)->calcBestParent(); + //bestParent = fRecoClusters->at(i)->GetBestParent(); + Log("Check 1: Particle ID "+to_string(bestParent),v_debug,verbosity); + if (bestParent >= fMCParticles->size()) { + Log("Invalid particle ID " + to_string(bestParent) + " vs number of particles: " + to_string(fMCParticles->size()) + "; skipping.",v_error,verbosity); + continue; + } + + //for (std::pair apair : *fMCParticleIndexMap) { + //Log("Testing match for particle: "+to_string(apair.first) + " corresponding to ID " + to_string(apair.second), v_debug, verbosity); + //if(apair.second==bestParent){ + //bestParticleID=apair.first; + bestPDG=fMCParticles->at(bestParent).GetPdgCode(); + //if (bestPDG == 22) { + + if (fMCParticles->at(bestParent).GetFlag() != 0) { + Log("Flagged particle! PDG "+to_string(bestPDG)+" excluding.",v_message,verbosity); + bestPDG=-5; + } + if (fMCParticles->at(bestParent).GetParentPdg() != 0) { + + Log("NeutCheck: PDG "+to_string(bestPDG)+" Finding parent.",v_message,verbosity); + fClusterParentPDG.push_back(fMCParticles->at(bestParent).GetParentPdg()); + Log("Parent PDG "+to_string(fClusterParentPDG.at(fClusterParentPDG.size()-1)),v_message,verbosity); + } + else { + Log("NeutCheck tool: no parent found. Treating as primary. PDG "+to_string(bestPDG), v_message, verbosity); + fClusterParentPDG.push_back(0); + if (bestPDG == 2212 && fRecoClusters->at(i)->GetTime()>10000) { + + Log("Primary proton found ID " + to_string(bestParent) + " at time " +to_string(fRecoClusters->at(i)->GetTime())+" but particle start, stop "+to_string(fMCParticles->at(bestParent).GetStartTime())+", "+to_string(fMCParticles->at(bestParent).GetStopTime()),v_debug,verbosity); + } + if (bestPDG == 13 && fRecoClusters->at(i)->GetTime() > 10000) { + + Log("Primary muon found ID " + to_string(bestParent) + " at time " + to_string(fRecoClusters->at(i)->GetTime()) + " but particle start, stop " + to_string(fMCParticles->at(bestParent).GetStartTime()) + ", " + to_string(fMCParticles->at(bestParent).GetStopTime()), v_debug, verbosity); + } + if (bestPDG == 2112 && fRecoClusters->at(i)->GetTime() > 10000) { + + Log("Primary neutron found ID " + to_string(bestParent) + " at time " + to_string(fRecoClusters->at(i)->GetTime()) + " but particle start, stop " + to_string(fMCParticles->at(bestParent).GetStartTime()) + ", " + to_string(fMCParticles->at(bestParent).GetStopTime()), v_debug, verbosity); + } + } + + //} + fRecoClusters->at(i)->SetPDG(bestPDG); + fClusterParticleEnergy.push_back(fMCParticles->at(bestParent).GetStartEnergy()); + + //} + + //} + /*if (matchCheck == false) { + Log("NeutCheck: Failed to match Particle. IndexMap size: "+to_string(fMCParticleIndexMap->size()), v_error, verbosity); + //continue; + }*/ + + + //ClusterSearchCheck->Fill(fRecoClusters->at(i)->GetClusterMode()); + //recoPDGHist->Fill(fRecoClusters->at(i)->GetPDG()); + fClusterMode.push_back(fRecoClusters->at(i)->GetClusterMode()); + fClusterPDG.push_back(fRecoClusters->at(i)->GetPDG()); + if (fClusterPDG.at(fClusterPDG.size() - 1) == 2112) { + fNeutronMult++; + + } + + fClusterCharge.push_back(fRecoClusters->at(i)->GetCharge()); + fClusterPurity.push_back(fRecoClusters->at(i)->Purity()); + fClusterTime.push_back(fRecoClusters->at(i)->GetTime()); + fClusterCB.push_back(fRecoClusters->at(i)->GetCB()); + fClusterHits.push_back(fRecoClusters->at(i)->GetNDigits()); + //if(fRecoClusters->at(i)->GetTime()>10000) recoPDGHistDelayed->Fill(fRecoClusters->at(i)->GetPDG()); + + //AS0RecoCheck->Fill(fRecoClusters->at(i)->GetAS(0)); + //AS1RecoCheck->Fill(fRecoClusters->at(i)->GetAS(1)); + + fClusterAS0.push_back(fRecoClusters->at(i)->GetAS(0)); + fClusterAS1.push_back(fRecoClusters->at(i)->GetAS(1)); + fClusterAMD.push_back(fRecoClusters->at(i)->GetAMD()); + fClusterCA.push_back(fRecoClusters->at(i)->GetCA()); + if(fRecoClusters->at(i)->GetAS(0)>0){ + //AS2RecoCheck->Fill(fRecoClusters->at(i)->GetAS(1)/fRecoClusters->at(i)->GetAS(0)); + fClusterAS2.push_back(fClusterAS1.at(fClusterAS1.size()-1)/fClusterAS0.at(fClusterAS0.size()-1)); + } + fClusterASC.push_back(fRecoClusters->at(i)->GetASC()); + CVX= fRecoClusters->at(i)->GetCV().X(); + CVY = fRecoClusters->at(i)->GetCV().Y(); + CVZ = fRecoClusters->at(i)->GetCV().Z(); + CVR = pow(pow(CVX,2)+pow(CVY,2)+pow(CVZ,2),0.5); + fClusterCVX.push_back(CVX); + fClusterCVY.push_back(CVY); + fClusterCVZ.push_back(CVZ); + fClusterCVR.push_back(CVR); + //fClusterSA.push_back(fRecoClusters->at(i)->GetSA()); + + + } + + return; +} + +void NeutronCheck::ResetVariables() { + fMCEventNum=-9999; + fClusterNum.clear(); + fClusterCount=0; + fClusterMode.clear(); + fClusterPDG.clear(); + fClusterParentPDG.clear(); + fClusterParticleEnergy.clear(); + fClusterCharge.clear(); + fClusterPurity.clear(); + fClusterCB.clear(); + fClusterTime.clear(); + fClusterAS0.clear(); + fClusterAS1.clear(); + fClusterAS2.clear(); + fClusterASC.clear(); + //fClusterSA.clear(); + fClusterCVX.clear(); + fClusterCVY.clear(); + fClusterCVZ.clear(); + fClusterCVR.clear(); + fClusterAMD.clear(); + fClusterCA.clear(); + fNeutronMult=0; + fMCNeutCapTimes.clear(); + fMCNeutCapX.clear(); + fMCNeutCapY.clear(); + fMCNeutCapZ.clear(); + + fParticlePDG.clear(); + fParticleParent.clear(); + fParticleStartEnergy.clear(); + fParticleStartTime.clear(); + + + fFinderClusterNum.clear(); + fFinderClusterCount=0; + fFinderPDG.clear(); + fFinderParentPDG.clear(); + fFinderCharge.clear(); + fFinderPurity.clear(); + fFinderCB.clear(); + fFinderTime.clear(); +} \ No newline at end of file diff --git a/UserTools/NeutronCheck/NeutronCheck.h b/UserTools/NeutronCheck/NeutronCheck.h new file mode 100644 index 000000000..c446ee6de --- /dev/null +++ b/UserTools/NeutronCheck/NeutronCheck.h @@ -0,0 +1,177 @@ +#ifndef NeutronCheck_H +#define NeutronCheck_H + +#include +#include +#include +#include + +#include "Tool.h" +#include "Particle.h" +#include "Position.h" +#include "RecoVertex.h" +#include "RecoCluster.h" + +#include "TTree.h" +#include "TFile.h" +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TMath.h" + + + +/** + * \class NeutronCheck + * + * This is a blank template for a Tool used by the script to generate a new custom tool. Please fill out the description and author information. +* +* $Author: F. A. Lemmons $ +* $Date: 2025/02/28 10:44:00 $ +* Contact: franklin.lemmons@mines.sdsmt.edu +*/ +class NeutronCheck: public Tool { + + + public: + + NeutronCheck(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise function used to clean up resources. + + + private: + + bool GetClusterInformation(); + double ASCheck(std::vector*& digits,int mode); + void CSCheck(); + void ResetVariables(); + + int verbosity; + std::string outfile; + bool useClean=0; + std::string EventStore; + bool fFinderCompare; + bool fParticleInfo; + + TTree* NeutCheckTree = nullptr; + + /* + -EventNumber + -Clusternumber + -number of clusters + -neutron multiplicity + + + -Cluster PDG + -Cluster Charge + -Cluster Time + -CB + -AS 0,1,2 + + */ + + std::vector cluster_neutron; + std::vector cluster_times_neutron; + std::vector cluster_charges_neutron; + std::vector cluster_cb_neutron; + std::vector cluster_times; + std::vector cluster_charges; + std::vector cluster_cb; + + std::map* fClusterToBestParticleID = nullptr; + std::map* fClusterToBestParticlePDG = nullptr; + std::map* fClusterEfficiency = nullptr; + std::map* fClusterPurityMap = nullptr; + std::map* fClusterTotalCharge = nullptr; + std::vector* fDigitList=nullptr; + vector* fRecoClusters; + + std::vector* fMCParticles; + std::vector fMCNeutCapTimes; + std::vector fMCNeutCapX; + std::vector fMCNeutCapY; + std::vector fMCNeutCapZ; + + + std::map>* m_all_clusters; + + Geometry* fGeometry = nullptr; ///< ANNIE Geometry + + /// \brief MC entry number + uint32_t fMCEventNum; + vector fClusterNum; + int fClusterCount; + + vector fClusterPDG; + vector fClusterParentPDG; + vector fClusterParticleEnergy; + vector fClusterCharge; + vector fClusterPurity; + vector fClusterCB; + vector fClusterTime; + vector fClusterAS0; + vector fClusterAS1; + vector fClusterAS2; + vector fClusterASC; + vector fClusterAMD; + vector fClusterCA; + //vector fClusterSA; + //vector fClusterCV; + vector fClusterCVX,fClusterCVY,fClusterCVZ,fClusterCVR; + + vector fClusterMode; + vector fClusterHits; + + vector fParticlePDG; + vector fParticleParent; + vector fParticleStartEnergy; + vector fParticleStartTime; + + vector fFinderClusterNum; + int fFinderClusterCount; + vector fFinderPDG; + vector fFinderParentPDG; + vector fFinderCharge; + vector fFinderPurity; + vector fFinderCB; + vector fFinderTime; + + + int fTrueNeutronMult; + int fTrueNeutronDelayed; + int fNeutronMult; + double true_Emu; + double true_Enu; + double TrueQ2; + + /// \brief trigger number + uint16_t fMCTriggerNum; + + /// \brief ANNIE event number + uint32_t fEventNumber; + + int fEventStatusApplied; + int fEventStatusFlagged; + + + + int reco_Clusters; + + TFile* fOutput_tfile; + + //verbosity variables + int v_error = 0; + int v_warning = 1; + int v_message = 2; + int v_debug = 3; + int vv_debug = 4; + + +}; + + +#endif diff --git a/UserTools/NeutronCheck/README.md b/UserTools/NeutronCheck/README.md new file mode 100644 index 000000000..d2be28214 --- /dev/null +++ b/UserTools/NeutronCheck/README.md @@ -0,0 +1,71 @@ +# NeutronCheck + +NeutronCheck +Output tool for RecoCluster information and parameters. Creates NeutCheckTree in the outfile.root, for use making analysis and comparison plots. + +## Data + + +From ANNIEEvent store +**MCEventNum** +*simple event number + +**ClusterToBestParticleID** +**ClusterToBestParticlePDG** +**ClusterEfficiency** +**ClusterPurity** +**ClusterTotalCharge** +**ClusterMapMC** +*ClusterFinder-made cluster maps from BackTracker tool for comparison + +**MCParticles** +*vector of all MC Particles in the event + +**TrueMuonEnergy** +*True energy of the muon for inclusion in the output + +**TrackId_to_MCParticleIndex** DEPRECATED +*map of track ID to MCParticle Index. No longer used, as MCParticle Index is used directly in all relevant cases. + + + +From RecoEvent store +**EventFlagged** +*Which event selection cuts were triggered by the current event. + +**ClusterIndiccesNeutron** +**ClusterTimesNeutron** +**ClusterChargesNeutron** +**ClusterCBNeutron** +*ClusterFinder-made cluster maps of neutrons specifically for comparison + +**ClusterTimes** +**ClusterCharges** +**ClusterCB** +*ClusterFinder-made cluster maps of all particles for comparison + +**RecoClusters** +*vector list of all RecoClusters generated from all iterations of ClusterSearcher + + + +From GenieInfo Store +**NeutrinoEnergy +*True neutrino energy for inclusion in the output + +Output .root file +**NeutCheckTree** +*Tree written into output file + + +## Configuration + +Describe any configuration variables for NeutronCheck. + +``` +verbosity 1 +outfile /path/to/output/file #'.root' is added to the end. +UseClean 0 #whether to exclude events that do not pass event selection +FinderCompare 0 #whether to add similar plots from ClusterFinder map-clusters to output for comparison purposes (not all parameters are written into Finder clusters) +ParticleInfo 1 #whether to add true particle information directly from MCParticles list to output +``` diff --git a/UserTools/Unity.h b/UserTools/Unity.h index b25c2ab57..26b165117 100644 --- a/UserTools/Unity.h +++ b/UserTools/Unity.h @@ -187,3 +187,5 @@ #include "LAPPDLoadStore.h" #include "PMTWaveformSim.h" #include "LAPPDWaveformDisplay.h" +#include "ClusterSearcher.h" +#include "NeutronCheck.h" diff --git a/configfiles/ClusterSearcher/ClusterSearcher2Config b/configfiles/ClusterSearcher/ClusterSearcher2Config new file mode 100755 index 000000000..20672cf1c --- /dev/null +++ b/configfiles/ClusterSearcher/ClusterSearcher2Config @@ -0,0 +1,21 @@ +# ClusterSearcher2 config file + +verbosity 5 +ClusterMode 2 #searching track-like clusters +Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +PmtMinPulseHeight 20 #minimum pulse height +PmtNeighbourRadius 600 #digit neighbouring distance [cm] +PmtMinNeighbourDigits 2 #minimum neighbour digits +PmtClusterRadius 600 #digit clustering distance [cm] +PmtTimeWindowN 100 #neighbouring time window [ns] +PmtTimeWindowC 100 #clustering time window [ns] +PmtMinHitsPerCluster 1 #number of hits per cluster +LappdMinPulseHeight 0 #minimum pulse height +LappdNeighbourRadius 250 #digit neighbouring distance [cm] +LappdMinNeighbourDigits 20 #minimum neighbour digits +LappdClusterRadius 250 #digit clustering distance [cm] +LappdTimeWindowN 1 #neighbouring time window [ns] +LappdTimeWindowC 1 #clustering time window [ns] +LappdMinHitsPerCluster 5 #number of digits per cluster +#MinClusterDigits 50 #minimum clustered digits (LAPPD+PMT) +MinClusterDigits 1 #minimum clustered digits (PMT-only) diff --git a/configfiles/ClusterSearcher/ClusterSearcherConfig b/configfiles/ClusterSearcher/ClusterSearcherConfig new file mode 100755 index 000000000..c135e4b6f --- /dev/null +++ b/configfiles/ClusterSearcher/ClusterSearcherConfig @@ -0,0 +1,21 @@ +# ClusterSearcher config file + +verbosity 5 +ClusterMode 1 #searching track-like clusters +Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +PmtMinPulseHeight 20 #minimum pulse height +PmtNeighbourRadius 60 #digit neighbouring distance [cm] +PmtMinNeighbourDigits 2 #minimum neighbour digits +PmtClusterRadius 60 #digit clustering distance [cm] +PmtTimeWindowN 10 #neighbouring time window [ns] +PmtTimeWindowC 100 #clustering time window [ns] +PmtMinHitsPerCluster 4 #number of hits per cluster +LappdMinPulseHeight 0 #minimum pulse height +LappdNeighbourRadius 25 #digit neighbouring distance [cm] +LappdMinNeighbourDigits 20 #minimum neighbour digits +LappdClusterRadius 25 #digit clustering distance [cm] +LappdTimeWindowN 1 #neighbouring time window [ns] +LappdTimeWindowC 1 #clustering time window [ns] +LappdMinHitsPerCluster 5 #number of digits per cluster +#MinClusterDigits 50 #minimum clustered digits (LAPPD+PMT) +MinClusterDigits 4 #minimum clustered digits (PMT-only) diff --git a/configfiles/ClusterSearcher/DigitBuilderConfig b/configfiles/ClusterSearcher/DigitBuilderConfig new file mode 100755 index 000000000..a2aac4e86 --- /dev/null +++ b/configfiles/ClusterSearcher/DigitBuilderConfig @@ -0,0 +1,12 @@ +# DigitBuilder config file + +verbosity 5 +ParametricModel 1 +#Reading in MC files +IsMC 1 +DigitChargeThr 20 #need to change the name +# There are three configurations: "PMT_only", "LAPPD_only", "All" +PhotoDetectorConfiguration PMT_only +#File must be in /pnfs/ space when loading on Fermilab cluster +LAPPDIDFile ./configfiles/VertexReco/PhaseIIRecoFirinne/LAPPDIDs.txt +StripHit 0 diff --git a/configfiles/ClusterSearcher/LoadGeometryConfig b/configfiles/ClusterSearcher/LoadGeometryConfig new file mode 100644 index 000000000..79efc3811 --- /dev/null +++ b/configfiles/ClusterSearcher/LoadGeometryConfig @@ -0,0 +1,10 @@ +verbosity 0 +LAPPDChannelCount 60 +FACCMRDGeoFile ./configfiles/LoadGeometry/FullMRDGeometry_09_29_20.csv +#FACCMRDGeoFile ./configfiles/LoadGeometry/FullMRDGeometry.csv +DetectorGeoFile ./configfiles/LoadGeometry/DetectorGeometrySpecs.csv +LAPPDGeoFile ./configfiles/LoadGeometry/LAPPDGeometry.csv +TankPMTGeoFile ./configfiles/LoadGeometry/FullTankPMTGeometry.csv +TankPMTGainFile ./configfiles/LoadGeometry/ChannelSPEGains2023.csv +TankPMTTimingOffsetFile ./configfiles/LoadGeometry/TankPMTTimingOffsets.csv +AuxiliaryChannelFile ./configfiles/LoadGeometry/AuxChannels.csv diff --git a/configfiles/ClusterSearcher/LoadWCSimConfig b/configfiles/ClusterSearcher/LoadWCSimConfig new file mode 100644 index 000000000..9bd2001a0 --- /dev/null +++ b/configfiles/ClusterSearcher/LoadWCSimConfig @@ -0,0 +1,17 @@ +#LoadWCSim Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 0 + +#InputFile /Data/Ambesim/pmt/*.root +#InputFile /Data/HighEnergy/wcsim-thousand-high_0.root +InputFile /Data/0824Beamsim/pmt/*.root + +WCSimVersion 3 ## should reflect the WCSim version of the files being loaded +HistoricTriggeroffset 0 ## time offset of digits relative to the trigger +UseDigitSmearedTime 1 ## whether to use smeared digit time (T), or true time of first photon (F) +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? +ChankeyToPMTIDMap ./configfiles/BeamClusterAnalysisMC/Chankey_WCSimID_v7.txt +SplitSubTriggers 0 # should subtriggers be loaded in separate Execute steps? diff --git a/configfiles/ClusterSearcher/LoadWCSimLAPPDConfig b/configfiles/ClusterSearcher/LoadWCSimLAPPDConfig new file mode 100644 index 000000000..f39db46d8 --- /dev/null +++ b/configfiles/ClusterSearcher/LoadWCSimLAPPDConfig @@ -0,0 +1,11 @@ +#LoadWCSimLAPPD Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 1 + +#InputFile /Data/Ambesim/lappd/*.root +#InputFile /Data/HighEnergy/wcsim-thousand-high_lappd_0.root +InputFile /Data/0824Beamsim/lappd/*.root + +WCSimVersion 3 ## should reflect the WCSim version of the files being loaded +InnerStructureRadius 1.3545 ## octagonal inner structure radius in m (from drawings 106.64") +DrawDebugGraphs 0 ## whether to draw TPolyMarker3D's of hits diff --git a/configfiles/ClusterSearcher/MCParticlePropertiesConfig b/configfiles/ClusterSearcher/MCParticlePropertiesConfig new file mode 100644 index 000000000..714c40f32 --- /dev/null +++ b/configfiles/ClusterSearcher/MCParticlePropertiesConfig @@ -0,0 +1,3 @@ +# MCParticleProperties configuration file + +verbosity 5 diff --git a/configfiles/ClusterSearcher/NeutronCheckConfig b/configfiles/ClusterSearcher/NeutronCheckConfig new file mode 100644 index 000000000..ccb5d89a3 --- /dev/null +++ b/configfiles/ClusterSearcher/NeutronCheckConfig @@ -0,0 +1,5 @@ +verbosity 5 +outfile /Data/0824Beamsim/NeutronCheck_CS +UseClean 1 +FinderCompare 0 +ParticleInfo 0 \ No newline at end of file diff --git a/configfiles/ClusterSearcher/ToolChainConfig b/configfiles/ClusterSearcher/ToolChainConfig new file mode 100644 index 000000000..df77bb067 --- /dev/null +++ b/configfiles/ClusterSearcher/ToolChainConfig @@ -0,0 +1,23 @@ +#ToolChain dynamic setup file + +##### Runtime Paramiters ##### +verbose 1 +error_level 0 # 0= do not exit, 1= exit on unhandeled errors only, 2= exit on unhandeled errors and handeled errors +attempt_recover 1 + +###### Logging ##### +log_mode Interactive # Interactive=cout , Remote= remote logging system "serservice_name Remote_Logging" , Local = local file log; +log_local_path ./log +log_service LogStore + +###### Service discovery ##### +service_publish_sec -1 +service_kick_sec -1 + +##### Tools To Add ##### +Tools_File ./configfiles/ClusterSearcher/ToolsConfig + +##### Run Type ##### +Inline -1 +Interactive 0 + diff --git a/configfiles/ClusterSearcher/ToolsConfig b/configfiles/ClusterSearcher/ToolsConfig new file mode 100755 index 000000000..adff3b930 --- /dev/null +++ b/configfiles/ClusterSearcher/ToolsConfig @@ -0,0 +1,9 @@ +LoadGeometry LoadGeometry ./configfiles/ClusterSearcher/LoadGeometryConfig +LoadWCSim LoadWCSim ./configfiles/ClusterSearcher/LoadWCSimConfig +LoadWCSimLAPPD LoadWCSimLAPPD ./configfiles/ClusterSearcher/LoadWCSimLAPPDConfig +MCParticleProperties MCParticleProperties ./configfiles/ClusterSearcher/MCParticlePropertiesConfig +DigitBuilder DigitBuilder ./configfiles/ClusterSearcher/DigitBuilderConfig +ClusterSearcher ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcherConfig +ClusterSearcher2 ClusterSearcher2 ./configfiles/ClusterSearcher/ClusterSearcher2Config +EventSelector EventSelector ./configfiles/ClusterSearcher/EventSelectorConfig +NeutronCheck NeutronCheck ./configfiles/ClusterSearcher/NeutronCheckConfig From 36a62bfdefa5c4becbd2f3cab63439fdb7102b52 Mon Sep 17 00:00:00 2001 From: root Date: Tue, 13 Jan 2026 02:49:23 +0000 Subject: [PATCH 2/5] Update to previous commit in accordance with feedback: -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 --- DataModel/RecoCluster.cpp | 393 ++++++++---------- DataModel/RecoCluster.h | 46 +- DataModel/RecoDigit.h | 15 + UserTools/ClusterSearcher/ClusterSearcher.cpp | 176 ++------ UserTools/ClusterSearcher/ClusterSearcher.h | 18 +- UserTools/ClusterSearcher/README.md | 3 +- UserTools/DigitBuilder/DigitBuilder.cpp | 217 +--------- .../ClusterSearcher/ClusterSearcher2Config | 3 +- .../ClusterSearcher/ClusterSearcherConfig | 3 +- 9 files changed, 294 insertions(+), 580 deletions(-) diff --git a/DataModel/RecoCluster.cpp b/DataModel/RecoCluster.cpp index bea56c09f..db6693efd 100644 --- a/DataModel/RecoCluster.cpp +++ b/DataModel/RecoCluster.cpp @@ -16,36 +16,37 @@ RecoCluster::~RecoCluster() void RecoCluster::Reset() { - int j=fDigitList.size(); - for(int i=0;isize(); + /*for (int i = 0; iclear(); } -static bool CompareTimes(RecoDigit *rd1, RecoDigit *rd2) +static bool CompareTimes(RecoDigit rd1, RecoDigit rd2) { - return ( rd1->GetCalTime() > rd2->GetCalTime() ); + return ( rd1.GetCalTime() > rd2.GetCalTime() ); } void RecoCluster::SortCluster() { - sort(fDigitList.begin(), fDigitList.end(), CompareTimes); + sort(fDigitList->begin(), fDigitList->end(), CompareTimes); } -void RecoCluster::AddDigit(RecoDigit* digit) +void RecoCluster::AddDigit(RecoDigit digit) { - fDigitList.push_back(digit); + //digit.AddCluster(fClusterMode); + fDigitList->push_back(digit); } -RecoDigit* RecoCluster::GetDigit(int n) +RecoDigit RecoCluster::GetDigit(int n) { - return (RecoDigit*)(fDigitList.at(n)); + return (RecoDigit)(fDigitList->at(n)); } int RecoCluster::GetNDigits() { - return fDigitList.size(); + return fDigitList->size(); } void RecoCluster::SetClusterMode(int cmode) @@ -58,10 +59,34 @@ int RecoCluster::GetClusterMode() return fClusterMode; } +bool RecoCluster::CheckFilter() { + map ModeDigitMap; + vector DigitModes; + double baseline; + + fIsFiltered=true; + for (int i=0;isize();i++) { + RecoDigit adigit=fDigitList->at(i); + DigitModes=adigit.GetClusteredModes(); + for (int j = 0; j < DigitModes.size(); j++) { + if(!ModeDigitMap.count(DigitModes.at(j))) ModeDigitMap.emplace(DigitModes.at(j),adigit.GetCalCharge()); + else ModeDigitMap.at(DigitModes.at(j))+=adigit.GetCalCharge(); + } + + } + if(ModeDigitMap.size()==0) return fIsFiltered; + baseline=ModeDigitMap.at(fClusterMode); + for (pair apair : ModeDigitMap) { + if (apair.first>=fClusterMode) continue; + if(apair.second/baseline>0.6) fIsFiltered=false; + } + return fIsFiltered; +} + void RecoCluster::CalcTime() { double sum = 0; - for (RecoDigit* i_digit : fDigitList) { - sum += i_digit->GetCalTime(); + for (int i=0; isize(); i++) { + sum += fDigitList->at(i).GetCalTime(); } clusterTime = sum / GetNDigits(); @@ -69,8 +94,8 @@ void RecoCluster::CalcTime() { void RecoCluster::CalcCharge() { double sum=0; - for (RecoDigit* i_digit : fDigitList) { - sum += i_digit->GetCalCharge(); + for (int i = 0; i < fDigitList->size(); i++) { + sum += fDigitList->at(i).GetCalCharge(); } clusterCharge = sum; } @@ -79,10 +104,10 @@ void RecoCluster::CalcCB() { //calculate unfiltered CB double total_Q = 0; double total_QSquared = 0; - for (int i = 0; i < fDigitList.size(); i++) { + for (int i = 0; i < fDigitList->size(); i++) { //if(unfilteredDigits->at(i)->GetDigitType()==RecoDigit::PMT8inch){ - if (fDigitList.at(i)->GetFilterStatus()) { - double tube_charge = fDigitList.at(i)->GetCalCharge(); + if (fDigitList->at(i).GetFilterStatus()) { + double tube_charge = fDigitList->at(i).GetCalCharge(); total_Q += tube_charge; total_QSquared += (tube_charge * tube_charge); //} @@ -103,58 +128,50 @@ void RecoCluster::CalcAS() { double max_angle = 0; double angle; Position i_position, j_position; - for (int i = 0; i < fDigitList.size(); i++) { - i_position = fDigitList.at(i)->GetPosition(); - for (int j = 0; j < fDigitList.size(); j++) { + double max_charge = 0; + double max_angle2 = 0; + int max_index = 0; + Position max_position; + double max_chargeangle = 0; + double chargeangle; + double iq, jq; + + for (int i = 0; i < fDigitList->size(); i++) { + i_position = fDigitList->at(i).GetPosition(); + for (int j = 0; j < fDigitList->size(); j++) { if (j == i)continue; - j_position = fDigitList.at(j)->GetPosition(); + j_position = fDigitList->at(j).GetPosition(); angle = j_position.Angle(i_position); if (angle > max_angle)max_angle = angle; } - } - AS0 = max_angle; - - - - double max_charge = 0; - double max_angle2 = 0; - int max_index = 0; - Position max_position; - for (int i = 0; i < fDigitList.size(); i++) { - if (fDigitList.at(i)->GetCalCharge() > max_charge) { - max_charge = fDigitList.at(i)->GetCalCharge(); + + if (fDigitList->at(i).GetCalCharge() > max_charge) { + max_charge = fDigitList->at(i).GetCalCharge(); max_index = i; - max_position = fDigitList.at(i)->GetPosition(); + max_position = fDigitList->at(i).GetPosition(); } - } + iq = fDigitList->at(i).GetCalCharge(); + for (int j = 0; j < fDigitList->size(); j++) { + if (j == i)continue; + j_position = fDigitList->at(j).GetPosition(); + jq = fDigitList->at(j).GetCalCharge(); + chargeangle = iq * jq * j_position.Angle(i_position); + if (chargeangle > max_chargeangle)max_chargeangle = chargeangle; + } + } + AS0 = max_angle; - for (int i = 0; i < fDigitList.size(); i++) { + for (int i = 0; i < fDigitList->size(); i++) { if (i == max_index)continue; - i_position = fDigitList.at(i)->GetPosition(); + i_position = fDigitList->at(i).GetPosition(); angle = i_position.Angle(max_position); if (angle > max_angle2)max_angle2 = angle; } AS1 = max_angle2; - - double max_chargeangle = 0; - double chargeangle; - double iq,jq; - for (int i = 0; i < fDigitList.size(); i++) { - i_position = fDigitList.at(i)->GetPosition(); - iq=fDigitList.at(i)->GetCalCharge(); - for(int j = 0; j < fDigitList.size(); j++) { - if(j==i)continue; - j_position = fDigitList.at(j)->GetPosition(); - jq = fDigitList.at(j)->GetCalCharge(); - chargeangle=iq*jq*j_position.Angle(i_position); - if(chargeangle>max_chargeangle)max_chargeangle=chargeangle; - } - } - ASC=max_chargeangle; } @@ -172,9 +189,11 @@ void RecoCluster::CalcParameters() { this->CalcCB(); this->CalcAS(); this->CalcSA(); - this->CalcCA(); + this->CalcAW(); this->CalcCV(); this->CalcAMD(); + this->CalcPlanaritySphericity(); + this->CalcTR(); //Other needed parameters here } @@ -188,21 +207,16 @@ void RecoCluster::SetParticle(int pID, int pPDG, double eff, double pur) { int RecoCluster::calcBestParent() { std::map ParticletoClusterCharge; - std::vector DigiParents; int parentCandidate; int maxIndex=0; double digitCharge; double maxCharge=0; int tempPDG = -5; - for (RecoDigit* i_digit : fDigitList) { - cout<<"Digit parent list size: "<GetParents().size()<GetParents().size() == 0) { - cout<<"Digit found with no parents!!!\n"; - } - for(int i=0; iGetParents().size(); i++){ - parentCandidate=i_digit->GetParents().at(i); - cout<<"ClusterParent candidate ID: "<GetCalCharge(); + for (int i = 0; i < fDigitList->size(); i++) { + RecoDigit i_digit=fDigitList->at(i); + for(int j=0; j apair : *fMCParticleIndexMap) { - if (apair.second == maxIndex) { - tempPDG=apair.first; - } - }*/ - - //if(tempPDG==-5) return false; - - cout << "Particle " << maxIndex << " wins with charge " << maxCharge << endl; + //cout << "Particle " << maxIndex << " wins with charge " << maxCharge << endl; bestParticleID=maxIndex; - //bestParticlePDG=tempPDG; purity=maxCharge/clusterCharge; return maxIndex; } -std::vector RecoCluster::convexHull() { - - double tempX, tempY, tempPhi, tempTheta; - std::vector effX, effY; - double sumX=0; - double sumY=0; - Position pos; - //std::vector hullDigits; - std::sort(fDigitList.begin(), fDigitList.end()); - for (RecoDigit* adigit : fDigitList) { - pos = adigit->GetPosition(); - tempX = sin(pos.GetTheta()) * cos(pos.GetPhi()); - tempY = sin(pos.GetTheta()) * sin(pos.GetPhi()); - effX.push_back(tempX); - effY.push_back(tempY); - sumX+=tempX; - sumY+=tempY; - } - - TwoDCenter.SetX(sumX/fDigitList.size()); - TwoDCenter.SetY(sumY/fDigitList.size()); - - TwoDCenter.SetZ(0); //2D projection center; used for circularity test. - - std::vector lower, upper; - double diffX, diffY; - double diffprevX, diffprevY; - double cross; - - - // Lower hull - for (int i = 0; i < fDigitList.size(); i++) { - if (lower.size() >= 2) { - diffprevX = effX[lower[lower.size() - 1]] - effX[lower[lower.size() - 2]]; - diffprevY = effY[lower[lower.size() - 1]] - effY[lower[lower.size() - 2]]; - diffX = effX[i] - effX[lower[lower.size() - 1]]; - diffY = effY[i] - effY[lower[lower.size() - 1]]; - cross = diffprevX * diffY - diffprevY * diffX; - } - while (lower.size() >= 2 && cross <= 0) { - //(lower.size() >= 2 && (lower[lower.size() - 1] - lower[lower.size() - 2]).cross(p - lower[lower.size() - 1]) <= 0) { - lower.pop_back(); - if (lower.size() >= 2) { - diffprevX = effX[lower[lower.size() - 1]] - effX[lower[lower.size() - 2]]; - diffprevY = effY[lower[lower.size() - 1]] - effY[lower[lower.size() - 2]]; - diffX = effX[i] - effX[lower[lower.size() - 1]]; - diffY = effY[i] - effY[lower[lower.size() - 1]]; - cross = diffprevX * diffY - diffprevY * diffX; - } - - } - lower.push_back(i); - std::cout<<"lowersize,effX,effY,finalcross: "<= 0; --i) { - if (upper.size() >= 2) { - diffprevX = effX[upper[upper.size() - 1]] - effX[upper[upper.size() - 2]]; - diffprevY = effY[upper[upper.size() - 1]] - effY[upper[upper.size() - 2]]; - diffX = effX[i] - effX[upper[upper.size() - 1]]; - diffY = effY[i] - effY[upper[upper.size() - 1]]; - cross = diffprevX * diffY - diffprevY * diffX; - } - while (upper.size() >= 2 && cross <= 0) { - //(upper[upper.size() - 1] - upper[upper.size() - 2]).cross(points[i] - upper[upper.size() - 1]) <= 0) { - upper.pop_back(); - if (upper.size() >= 2) { - diffprevX = effX[upper[upper.size() - 1]] - effX[upper[upper.size() - 2]]; - diffprevY = effY[upper[upper.size() - 1]] - effY[upper[upper.size() - 2]]; - diffX = effX[i] - effX[upper[upper.size() - 1]]; - diffY = effY[i] - effY[upper[upper.size() - 1]]; - cross = diffprevX * diffY - diffprevY * diffX; - } - } - upper.push_back(i); - std::cout << "uppersize,effX,effY,finalcross: " << upper.size() << "," << effX[i] << "," << effY[i] << "," << cross << endl; - } - - // Remove the last point of each half because it is repeated at the beginning of the other half - lower.pop_back(); - upper.pop_back(); - - // Concatenate lower and upper hull to get the convex hull - lower.insert(lower.end(), upper.begin(), upper.end()); - - for (int i = 0; i < fDigitList.size(); i++) { - fDigitList.at(i)->ResetFilter(); - } - - for (int i = 0; i < lower.size(); i++) { - hullDigits.push_back(fDigitList[lower[i]]); - fDigitList.at(lower[i])->PassFilter(); - std::cout<<"Adding digit index "<GetPosition() - fDigitList.at(j)->GetPosition()).Unit(); - ChargeVector+=(fDigitList.at(i)->GetCalCharge()*fDigitList.at(j)->GetCalCharge())*ud; + for (int i = 0; i < fDigitList->size(); i++) { + for (int j = i + 1; j < fDigitList->size(); j++) { + ud=(fDigitList->at(i).GetPosition() - fDigitList->at(j).GetPosition()).Unit(); + ChargeVector+=(fDigitList->at(i).GetCalCharge()*fDigitList->at(j).GetCalCharge())*ud; } } } @@ -354,11 +258,11 @@ void RecoCluster::CalcAMD() { double min_distance = 9999; double ave_min_distance=0; double distance; - for (int i = 0; i < fDigitList.size(); i++) { - Position ipos=fDigitList.at(i)->GetPosition(); - for (int j = 0; j < fDigitList.size(); j++) { + for (int i = 0; i < fDigitList->size(); i++) { + Position ipos=fDigitList->at(i).GetPosition(); + for (int j = 0; j < fDigitList->size(); j++) { if(j == i) continue; - Position jpos=fDigitList.at(j)->GetPosition(); + Position jpos=fDigitList->at(j).GetPosition(); distance = (jpos-ipos).Mag(); if(distancesize(); AMD=ave_min_distance; } @@ -375,9 +279,9 @@ void RecoCluster::CalcSA() { double aveX=0; double aveY=0; double aveZ=0; - for (int i = 0; i < fDigitList.size(); i++) { - Position ipos=fDigitList.at(i)->GetPosition(); - double iq=fDigitList.at(i)->GetCalCharge(); + for (int i = 0; i < fDigitList->size(); i++) { + Position ipos=fDigitList->at(i).GetPosition(); + double iq=fDigitList->at(i).GetCalCharge(); aveX += ipos.X() * iq; aveY += ipos.Y() * iq; aveZ += ipos.Z() * iq; @@ -389,20 +293,21 @@ void RecoCluster::CalcSA() { } //Containing Angle: angle from SpatialAverage unit vector which contains 68% of the total charge of the cluster -void RecoCluster::CalcCA() { //AW: Angular Width? 1-sigma angular width? +void RecoCluster::CalcAW() { //AW: Angular Width? 1-sigma angular width? std::map angle_charge_map; double angle,charge; double contained=0; - for (RecoDigit* idigit: fDigitList) { - angle=SpatialAverage.Angle(idigit->GetPosition()); - charge = idigit->GetCalCharge(); + for (int i = 0; i < fDigitList->size(); i++) { + RecoDigit idigit=fDigitList->at(i); + angle=SpatialAverage.Angle(idigit.GetPosition()); + charge = idigit.GetCalCharge(); if(angle_charge_map.count(angle)==0) angle_charge_map.emplace(angle,charge); else angle_charge_map.at(angle)+=charge; } for (std::pair apair: angle_charge_map) { - ContainingAngle=apair.first; + AngularWidth=apair.first; contained+=apair.second/clusterCharge; if(contained>0.68)break; } @@ -410,17 +315,16 @@ void RecoCluster::CalcCA() { //AW: Angular Width? 1-sigma angular width? } //Planarity and Sphericity: measure of how planar or spherical cluster is; currently in-process of translating from reference -/*void ComputePlanaritySphericity_ROOT(const std::vector& digits, - double& P, double& S) +void RecoCluster::CalcPlanaritySphericity() { - P = 0.0; S = 0.0; + double P = 0.0; double S = 0.0; // 1) Charge-weighted centroid - TVector3 mean(0, 0, 0); + Position mean(0, 0, 0); double wsum = 0.0; - for (int i = 0; isize(); i++) { + double w = fDigitList->at(i).GetCalCharge(); + mean += w * fDigitList->at(i).GetPosition(); wsum += w; } if (wsum <= 0.0) return; @@ -430,9 +334,10 @@ void RecoCluster::CalcCA() { //AW: Angular Width? 1-sigma angular width? // M = (1/wsum) * sum_i w_i * (r_i)(r_i)^T, with r_i = pos_i - mean TMatrixDSym cov(3); cov.Zero(); - for (const auto& d : digits) { - const double w = std::max(1e-12, d.charge); - const TVector3 r = d.pos - mean; + for (int i = 0; i < fDigitList->size(); i++) { + RecoDigit d = fDigitList->at(i); + double w = d.GetCalCharge(); + Position r = d.GetPosition() - mean; // Outer product r*r^T scaled by w cov(0, 0) += w * r.X() * r.X(); cov(0, 1) += w * r.X() * r.Y(); @@ -466,22 +371,70 @@ void RecoCluster::CalcCA() { //AW: Angular Width? 1-sigma angular width? // S = (3/2) * lambda3 / (lambda1 + lambda2 + lambda3) P = (lambda2 - lambda3) / sum; S = 1.5 * (lambda3 / sum); + + Planarity=P; + Sphericity=S; +} + +double RecoCluster::CalcBeta(int order) { + double Beta=0; + Position ipos,jpos; + double angle; + int NDigits=fDigitList->size(); + TF1 LegendreFunction("legendrefunction","ROOT::Math::legendre([0],x)",-1,1); + LegendreFunction.SetParameter(0,order); + for (int i = 0; i < fDigitList->size(); i++) { + for (int j = i + 1; j < fDigitList->size(); j++) { + ipos=fDigitList->at(i).GetPosition(); + jpos=fDigitList->at(j).GetPosition(); + angle=ipos.Angle(jpos); + Beta+=2 * LegendreFunction.Eval(cos(angle)) / (NDigits*(NDigits-1)); + } + } + + if(BetaParameters.count(order)) BetaParameters.at(order)=Beta; + else BetaParameters.emplace(order, Beta); + + return Beta; } -// ----------------------- -// Example usage -// ----------------------- -int main() { - std::vector digits; - digits.push_back({ TVector3(1.0, 0.0, 0.0), 2.0 }); - digits.push_back({ TVector3(0.0, 1.0, 0.1), 1.0 }); - digits.push_back({ TVector3(0.0,-1.0,-0.1), 1.5 }); - digits.push_back({ TVector3(-1.0, 0.0, 0.0), 1.2 }); - - double P = 0, S = 0; - ComputePlanaritySphericity_ROOT(digits, P, S); - - std::cout << "Planarity P = " << P << "\n"; - std::cout << "Sphericity S = " << S << "\n"; - return 0; +double RecoCluster::GetBeta(int order) { + if(BetaParameters.count(order))return BetaParameters.at(order); + else return -999; +} + +void RecoCluster::CalcTR() { + double firsttime=fDigitList->at(0).GetCalTime(); + double lasttime=fDigitList->at(fDigitList->size()-1).GetCalTime(); + double maxtime=fDigitList->at(0).GetCalTime(); + double maxcharge=fDigitList->at(0).GetCalCharge(); + + for (int i = 0; i < fDigitList->size(); i++) { + RecoDigit adigit=fDigitList->at(i); + if(adigit.GetCalTime()lasttime)lasttime=adigit.GetCalTime(); + if (adigit.GetCalCharge() > maxcharge) { + maxtime=adigit.GetCalTime(); + maxcharge=adigit.GetCalCharge(); + } + } + + TRTotal=lasttime-firsttime; + if(clusterTime-firsttime>lasttime-clusterTime)TRCluster=clusterTime-firsttime; + else TRCluster=lasttime-clusterTime; + if(maxtime-firsttime>lasttime-maxtime)TRQ=maxtime-firsttime; + else TRQ=lasttime-maxtime; + + TRRatioTQ=TRTotal/TRQ; + TRRatioTC=TRTotal/TRCluster; + TRRatioQC=TRQ/TRCluster; +} + +/*void RecoCluster::CleanDigits() { + for (int i = 0; i < fDigitList->size(); i++) { + if (fDigitList->at(i) != nullptr) { + delete fDigitList->at(i); + fDigitList->at(i)=nullptr; + } + } }*/ \ No newline at end of file diff --git a/DataModel/RecoCluster.h b/DataModel/RecoCluster.h index 5986a6669..d338c392a 100644 --- a/DataModel/RecoCluster.h +++ b/DataModel/RecoCluster.h @@ -5,6 +5,13 @@ #include "RecoDigit.h" #include "Position.h" #include +#include "TVector3.h" +#include "TMatrixDSym.h" +#include "TMatrixDSymEigen.h" +#include "TVectorD.h" +#include "TMath.h" +#include "TF1.h" +#include "Math/IFunction.h" class RecoCluster : public SerialisableObject { @@ -17,15 +24,16 @@ class RecoCluster : public SerialisableObject { void Reset(); void SortCluster(); - void AddDigit(RecoDigit* digit); + void AddDigit(RecoDigit digit); + inline void SetDigits(vector* indigits){fDigitList=indigits;} - RecoDigit* GetDigit(int n); + RecoDigit GetDigit(int n); void SetClusterMode(int cmode); int GetClusterMode(); - std::vector GetDigitList() {return this->fDigitList;} + std::vector* GetDigitList() {return this->fDigitList;} int GetNDigits(); @@ -49,15 +57,26 @@ class RecoCluster : public SerialisableObject { void CalcCV(); void CalcAMD(); void CalcSA(); - void CalcCA(); + void CalcAW(); + void CalcPlanaritySphericity(); + double CalcBeta(int order); + void CalcTR(); void CalcParameters(); int calcBestParent(); double GetAS(int mode); inline double GetASC(){return ASC;} inline double GetAMD(){return AMD;} inline Position GetSA(){return SpatialAverage;} - inline double GetCA(){return ContainingAngle;} - std::vector convexHull(); + inline double GetAW(){return AngularWidth;} + inline double GetPlanarity(){return Planarity;} + inline double GetSphericity(){return Sphericity;} + double GetBeta(int order); + inline double GetTimeRangeT() {return TRTotal;} + inline double GetTimeRangeQ() {return TRQ;} + inline double GetTimeRangeC() {return TRCluster;} + inline double GetTRRTQ(){return TRRatioTQ;} + inline double GetTRRTC() { return TRRatioTC; } + inline double GetTRRQC() { return TRRatioQC; } void SetParticle(int pID,int pPDG, double eff, double pur); inline int GetBestParent(){return bestParticleID;} @@ -66,6 +85,10 @@ class RecoCluster : public SerialisableObject { inline double Efficiency(){return efficiency;} inline double Purity(){return purity;} + bool CheckFilter(); + inline bool GetFilterStatus(){return fIsFiltered;} + //void CleanDigits(); + private: double clusterTime; double clusterCharge; @@ -74,16 +97,21 @@ class RecoCluster : public SerialisableObject { double AMD; Position ChargeVector; Position SpatialAverage; - double ContainingAngle; + double AngularWidth; + double Planarity=-999; + double Sphericity=-999; + std::map BetaParameters; + double TRTotal, TRQ, TRCluster; + double TRRatioTQ, TRRatioTC,TRRatioQC; int bestParticleID; int bestParticlePDG; double efficiency; double purity; + bool fIsFiltered=0; int fClusterMode = -999; - std::vector fDigitList; - std::vector hullDigits; + std::vector* fDigitList; Position TwoDCenter; diff --git a/DataModel/RecoDigit.h b/DataModel/RecoDigit.h index af64ee7cf..f766d0e45 100644 --- a/DataModel/RecoDigit.h +++ b/DataModel/RecoDigit.h @@ -33,6 +33,17 @@ class RecoDigit : public SerialisableObject{ fIsFiltered = 1; ANNIERecoObjectTable::Instance()->NewDigit(); } + RecoDigit(RecoDigit* origin) { + serialise=true; + fRegion=origin->GetRegion(); + fPosition = origin->GetPosition(); + fCalTime = origin->GetCalTime(); + fCalCharge=origin->GetCalCharge(); + fDigitType=origin->GetDigitType(); + fDetectorID=origin->GetDetectorID(); + fIsFiltered=origin->GetFilterStatus(); + + } ~RecoDigit() {/*ANNIERecoObjectTable::Instance()->DeleteDigit();*/} bool operator<(const RecoDigit other) const { @@ -44,6 +55,7 @@ class RecoDigit : public SerialisableObject{ inline double GetCalTime() const {return fCalTime;} inline double GetCalCharge() const {return fCalCharge;} inline bool GetFilterStatus() const {return fIsFiltered;} + inline vector GetClusteredModes() { return fClusterMode; } inline int GetDigitType() const {return fDigitType;} inline int GetDetectorID() const {return fDetectorID;} inline std::vector GetParents() const {return Parents;} @@ -57,6 +69,8 @@ class RecoDigit : public SerialisableObject{ inline void SetFilter(bool pass = 1) { fIsFiltered = pass;} inline void ResetFilter() {SetFilter(0);} inline void PassFilter() {SetFilter(1);} + inline void AddCluster(int InClusterMode) {fClusterMode.push_back(InClusterMode); } + inline void UnCluster(){fClusterMode.pop_back(); } inline void SetParents(std::vector parentsin) { Parents = parentsin; } bool Print() { @@ -76,6 +90,7 @@ class RecoDigit : public SerialisableObject{ double fCalTime; double fCalCharge; bool fIsFiltered; + vector fClusterMode; int fDigitType; int fDetectorID; std::vector Parents; diff --git a/UserTools/ClusterSearcher/ClusterSearcher.cpp b/UserTools/ClusterSearcher/ClusterSearcher.cpp index dda70f2ee..a882513fa 100644 --- a/UserTools/ClusterSearcher/ClusterSearcher.cpp +++ b/UserTools/ClusterSearcher/ClusterSearcher.cpp @@ -8,14 +8,6 @@ ClusterSearcher* ClusterSearcher::Instance() fgClusterSearcher = new ClusterSearcher(); } - if( !fgClusterSearcher ){ - assert(fgClusterSearcher); - } - - if( fgClusterSearcher ){ - - } - return fgClusterSearcher; } @@ -26,7 +18,7 @@ ClusterSearcher::~ClusterSearcher() { } bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){ - + if(verbosity)cout<<"Initializing ClusterSearcher"<; fClusteringParam->emplace("ClusterMode",fClusterMode); - fClusteringParam->emplace("Config",fConfig); fClusteringParam->emplace("PmtMinPulseHeight",fPmtMinPulseHeight); fClusteringParam->emplace("PmtNeighbourRadius",fPmtNeighbourRadius); fClusteringParam->emplace("PmtMinNeighbourDigits",fPmtMinNeighbourDigits); @@ -95,11 +67,7 @@ bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){ fClusteringParam->emplace("LappdTimeWindowC",fLappdTimeWindowC); fClusteringParam->emplace("MinClusterDigits",fMinClusterDigits); - if (fConfig!=0 && fConfig !=1 && fConfig !=2 && fConfig !=3 && fConfig !=4){ - Log("ClusterSearcher tool: Configuration <"+std::to_string(fConfig)+"> not recognized. Setting Config 3 (kPulseHeightAndClusters)",v_error,verbosity); - fConfig = ClusterSearcher::kPulseHeightAndClusters; - } - + Log("ClusterSearcher " + to_string(fClusterMode) + " configs Initialized", v_debug, verbosity); if (!fisMC){ ifstream file_singlepe(singlePEgains.c_str()); unsigned long temp_chankey; @@ -112,6 +80,7 @@ bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){ file_singlepe.close(); m_data->CStore.Get("pmt_tubeid_to_channelkey",pmt_tubeid_to_channelkey); } + Log("ClusterSearcher " + to_string(fClusterMode) + " !MC Initialized", v_debug, verbosity); // vector of selected digits fSelectAll = new std::vector; @@ -121,8 +90,8 @@ bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){ // only for test fSelectByTruthInfo = new std::vector; // vector of clusters - fClusterList = new std::vector; - fRecoClusters = new std::vector; + fClusterList = new std::vector; + fRecoClusters = new std::vector; //Set clustering parameters in the RecoEvent store std::map* pre_ClusteringParam = nullptr; // read existing container for parameters @@ -132,7 +101,7 @@ bool ClusterSearcher::Initialise(std::string configfile, DataModel &data){ pre_ClusteringParam->insert(fClusteringParam->begin(), fClusteringParam->end()); //m_data->Stores.at("RecoEvent")->Set("ClusteringParameters", pre_ClusteringParam); } - + Log("ClusterSearcher "+to_string(fClusterMode)+" Initialized",v_debug,verbosity); return true; } @@ -191,57 +160,25 @@ bool ClusterSearcher::Execute(){ digits->at(n)->ResetFilter(); } - // Run Clustering - // ================ - //std::vector* SelectDigitList = Run(digits); // Cluster digits for muon candidates - // Set Digit Filter // ========== - //for(int n=0; nsize()); n++ ) { - // RecoDigit* SelectDigit = (RecoDigit*)(SelectDigitList->at(n)); - // SelectDigit->PassFilter(); - //} SelectDigits(digits); fRecoClusters = this->RecoClusters(digits); // Digit clustering done! // ===== - - /*RecoCluster* dummyCluster1 = new RecoCluster(); - Position pos(10,10,10); - RecoDigit* dummyDigit1 = new RecoDigit(0,pos,1,10,10,5); - dummyCluster1->AddDigit(dummyDigit1); - pos.SetX(20); - RecoDigit* dummyDigit2 = new RecoDigit(0,pos,1,11,10,7); - dummyCluster1->AddDigit(dummyDigit2); - dummyCluster1->CalcParameters(); - - fRecoClusters->push_back(dummyCluster1); - Log("dummyCluster1 loaded in. Size of fRecoCLusters: "+to_string(fRecoClusters->size()),v_debug,verbosity); - - RecoCluster* dummyCluster2 = new RecoCluster(); - Position pos2(20, 20, 10); - RecoDigit* dummyDigit3 = new RecoDigit(0, pos2, 1, 10, 10, 9); - dummyCluster2->AddDigit(dummyDigit3); - pos2.SetX(30); - RecoDigit* dummyDigit4 = new RecoDigit(0, pos2, 1, 11, 10, 11); - dummyCluster2->AddDigit(dummyDigit4); - dummyCluster2->CalcParameters(); - - fRecoClusters->push_back(dummyCluster2); - Log("dummyCluster2 loaded in. Size of fRecoClusters: "+to_string(fRecoClusters->size()),v_debug,verbosity); - */ pre_RecoClusters = nullptr; // to read existing clusters bool cluster_status = m_data->Stores.at("RecoEvent")->Get("RecoClusters", pre_RecoClusters); Log("ClusterSearcher Tool: Final count of clusters: "+to_string(fRecoClusters->size()),v_debug,verbosity); if(!cluster_status) { - pre_RecoClusters=new std::vector; + pre_RecoClusters=new std::vector; for (int i = 0; i < fRecoClusters->size(); i++) { pre_RecoClusters->push_back(fRecoClusters->at(i)); } fIsHitClusteringDone = true; + Log("Hit Clustering done",v_debug,verbosity); m_data->Stores.at("RecoEvent")->Set("HitClusteringDone", fIsHitClusteringDone); m_data->Stores.at("RecoEvent")->Set("RecoClusters", pre_RecoClusters); } @@ -273,10 +210,10 @@ bool ClusterSearcher::Finalise(){ return true; } -void ClusterSearcher::Config(int config) +/*void ClusterSearcher::Config(int config) { ClusterSearcher::Instance()->SetConfig(config); -} +}*/ void ClusterSearcher::PmtMinPulseHeight(double min) { @@ -319,7 +256,6 @@ void ClusterSearcher::PrintParameters() std::cout << " Clustering Parameters: " << std::endl << " ClusterMode = " << fClusterMode<< std::endl - << " Config = " << fConfig<< std::endl << " PmtMinPulseHeight = " << fPmtMinPulseHeight << std::endl << " PmtNeighbourRadius = " << fPmtNeighbourRadius << std::endl << " PmtMinNeighbourDigits = " << fPmtMinNeighbourDigits << std::endl @@ -341,57 +277,6 @@ void ClusterSearcher::Reset() return; } -std::vector* ClusterSearcher::Run(std::vector* DigitList) -{ - if(verbosity>v_debug) std::cout << " *** ClusterSearcher::Run(...) *** " << std::endl; - - // input digit list - // ================ - std::vector* InputList = DigitList; - std::vector* OutputList = DigitList; - - // Select all digits - // ================= - InputList = ResetDigits(OutputList); - OutputList = (std::vector*)(this->SelectAll(InputList)); - OutputList = SelectDigits(OutputList); - if( fConfig==ClusterSearcher::kNone ) return OutputList; - - // Select by pulse height - // ====================== - InputList = ResetDigits(OutputList); - OutputList = (std::vector*)(this->SelectByPulseHeight(InputList)); - OutputList = SelectDigits(OutputList); - if( fConfig==ClusterSearcher::kPulseHeight ) return OutputList; - - // Select using neighbouring digits - // ================================ - InputList = ResetDigits(OutputList); - OutputList = (std::vector*)(this->SelectByNeighbours(InputList)); - OutputList = SelectDigits(OutputList); - if( fConfig==ClusterSearcher::kPulseHeightAndNeighbours ) return OutputList; - - // Select using clustered digits - // ============================= - InputList = ResetDigits(OutputList); - OutputList = (std::vector*)(this->SelectByClusters(InputList)); - OutputList = SelectDigits(OutputList); - if( fConfig==ClusterSearcher::kPulseHeightAndClusters ) return OutputList; - - if (fisMC){ - // Select using truth information (for simulation test only) - // ============================= - InputList = ResetDigits(OutputList); - OutputList = (std::vector*)(this->SelectByTruthInfo(InputList)); - OutputList = SelectDigits(OutputList); - if( fConfig==ClusterSearcher::kPulseHeightAndTruthInfo ) return OutputList; - } - - // return vector of selected digits - // ================================ - return OutputList; -} - // Set all filter status to 0 (default is 1) std::vector* ClusterSearcher::ResetDigits(std::vector* DigitList) { @@ -596,14 +481,15 @@ std::vector* ClusterSearcher::SelectByClusters(std::vector* ClusterList = (std::vector*)(this->RecoClusters(DigitList)); + std::vector* ClusterList = RecoClusters(DigitList); for(int icluster=0; iclustersize()); icluster++ ){ - RecoCluster* Cluster = (RecoCluster*)(ClusterList->at(icluster)); + RecoCluster Cluster = (ClusterList->at(icluster)); fRecoClusters->push_back(Cluster); - for(int idigit=0; idigitGetNDigits(); idigit++ ){ - RecoDigit* Digit = (RecoDigit*)(Cluster->GetDigit(idigit)); + for(int idigit=0; idigitpush_back(Digit); } } @@ -616,7 +502,7 @@ std::vector* ClusterSearcher::SelectByClusters(std::vector* ClusterSearcher::RecoClusters(std::vector* DigitList) +std::vector* ClusterSearcher::RecoClusters(std::vector* DigitList) { // delete cluster digits @@ -760,22 +646,38 @@ std::vector* ClusterSearcher::RecoClusters(std::vector } } + //std::cout <<"vClusterDigitCollection.size() == "<=fMinClusterDigits ){ - RecoCluster* cluster = new RecoCluster(); - cluster->SetClusterMode(fClusterMode); - fClusterList->push_back(cluster); + RecoCluster cluster; + cluster.SetClusterMode(fClusterMode); + + Log("Adding Digits",v_debug,verbosity); + vector* ClusteredDigits=new vector; for(int jdigit=0; jdigitGetRecoDigit()); - cluster->AddDigit(recodigit); + cout<<"check 2\n"; + RecoDigit* arecodigit=cdigit->GetRecoDigit(); + arecodigit->AddCluster(fClusterMode); + RecoDigit brecodigit(arecodigit); + ClusteredDigits->push_back(brecodigit); + cout<<" check 3\n"; + //cluster.AddDigit(*arecodigit); } - cluster->CalcParameters(); + cluster.SetDigits(ClusteredDigits); + ClusteredDigits = nullptr; + Log("Cluster has " + to_string(cluster.GetNDigits()) + " digits", v_debug, verbosity); + Log("ready to caluclate parameters.",v_debug,verbosity); + cluster.CalcParameters(); Log("ClusterSearcher: Clusters made: "+to_string(fClusterList->size()),v_debug,verbosity); + fClusterList->push_back(cluster); + } } } + std::cout <<"Number of clusters = "<size()<* Run(std::vector* digitlist); std::vector* ResetDigits(std::vector* digitlist); std::vector* SelectDigits(std::vector* digitlist); void SelectDigits(std::vector* digitlist); @@ -70,15 +68,13 @@ class ClusterSearcher: public Tool { std::vector* SelectByNeighbours(std::vector* digitlist); std::vector* SelectByClusters(std::vector* digitlist); std::vector* SelectByTruthInfo(std::vector* digitlist); //use truth information. Only for testing the code - std::vector* RecoClusters(std::vector* digitlist); + std::vector* RecoClusters(std::vector* digitlist); private: void Reset(); - // running mode - int fConfig; // clustering mode int fClusterMode; @@ -128,11 +124,11 @@ class ClusterSearcher: public Tool { std::vector* fSelectByTruthInfo; // vectors of clusters - std::vector* fClusterList; + std::vector* fClusterList; // vector of clusters (accessible to the CStore) - std::vector* fRecoClusters = nullptr; - std::vector* pre_RecoClusters = nullptr; + std::vector* fRecoClusters = nullptr; + std::vector* pre_RecoClusters = nullptr; // true vertex RecoVertex* fTrueVertex = 0; diff --git a/UserTools/ClusterSearcher/README.md b/UserTools/ClusterSearcher/README.md index 1fd30a030..7444b12d7 100644 --- a/UserTools/ClusterSearcher/README.md +++ b/UserTools/ClusterSearcher/README.md @@ -44,7 +44,8 @@ LappdTimeWindowN 1 #neighbouring time window [ns] LappdTimeWindowC 1 #clustering time window [ns] LappdMinHitsPerCluster 5 #number of digits per cluster MinClusterDigits 50 #minimum clustered digits (LAPPD+PMT) -MinClusterDigits 4 #minimum clustered digits (PMT-only) +MinClusterDigits 4 #minimum clustered digits (PMT-only) +SinglePEGains /path/to/gains/file #file in which SPE information can be found for PMT charge conversions; only needed on Data FirstRun 1 #Debug input for differentiating first iteration from later iterations ``` diff --git a/UserTools/DigitBuilder/DigitBuilder.cpp b/UserTools/DigitBuilder/DigitBuilder.cpp index b74f6c2ec..3bf94b7f8 100644 --- a/UserTools/DigitBuilder/DigitBuilder.cpp +++ b/UserTools/DigitBuilder/DigitBuilder.cpp @@ -18,7 +18,6 @@ DigitBuilder::~DigitBuilder() { bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ /////////////////// Usefull header /////////////////////// - cout<<"Initializing Tool DigitBuilder"< 3) std::cout << line << std::endl; //has our stuff; - if (line.find("#") != std::string::npos) continue; + if (line.empty() || line[0] == '#') continue; std::vector DataEntries; boost::split(DataEntries, line, boost::is_any_of(","), boost::token_compress_on); int channelkey = -9999; @@ -115,12 +114,7 @@ bool DigitBuilder::Initialise(std::string configfile, DataModel &data){ } - /*while (!file_singlepe.eof()) { - file_singlepe >> temp_chankey >> temp_gain; - if (file_singlepe.eof()) break; - pmt_gains.emplace(temp_chankey,temp_gain); - Log("DigitBuilder Tool: still collecting SPE gains: "+to_string(temp_gain), v_debug, verbosity); - }*/ + Log("DigitBuilder Tool: SPE gains done",v_debug,verbosity); file_singlepe.close(); @@ -176,16 +170,6 @@ bool DigitBuilder::Execute(){ Log("DigitBuilder Tool: ERROR retrieving hits in Data mode!",v_error,verbosity); return false; } - /*auto get_clusters = m_data->CStore.Get("ClusterMap", m_all_clusters); - if (!get_clusters){ - Log("DigitBuilder Tool: ERROR retrieving clustered hits (ClusterMap) in Data mode!",v_error,verbosity); - return false; - } - auto get_clusters_chankey = m_data->CStore.Get("ClusterMapDetkey",m_all_clusters_detkey); - if (!get_clusters_chankey){ - Log("DigitBuilder Tool: ERROR retrieving clustered chankeys (ClusterMapDetkey) in Data mode!",v_error,verbosity); - return false; - }*/ } /// Build RecoDigit @@ -201,7 +185,6 @@ bool DigitBuilder::Execute(){ } bool DigitBuilder::Finalise(){ - //delete fDigitList; fDigitList = 0; //Don't delete pointer to fDigitList, will be deleted by the BoostStore! Log("DigitBuilder exitting",v_message,verbosity); return true; } @@ -301,7 +284,7 @@ bool DigitBuilder::BuildMCPMTRecoDigit() { for(MCHit& ahit : hits){ 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<70) { + if(hitTime>-10 && hitTimepush_back(recoDigit); } @@ -575,159 +556,6 @@ bool DigitBuilder::BuildDataPMTRecoDigit(){ } } - - - - - - /* - - /// m_all_clusters is a std::map> - - if (m_all_clusters && m_all_clusters_detkey){ - int clustersize = m_all_clusters->size(); - Log("Clustersize of m_all_clusters: " + to_string(clustersize),v_debug,verbosity); - bool clusters_available = false; - bool muon_available = false; - if (clustersize != 0) clusters_available = true; - if (clusters_available){ - //determine the main cluster (max charge and in [0 ... 2000ns] time window) - double max_cluster = 0; - double max_charge = 0; - for(std::pair>&& apair : *m_all_clusters){ - std::vector&Hits = apair.second; - double time = 0; - int hits=0; - double charge = 0; - for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){ - hits++; - time+=Hits.at(i_hit).GetTime(); - charge+=Hits.at(i_hit).GetCharge(); - } - if (hits>0) { - time/=hits; - } - if (time > 2000.) continue; //not a beam muon if not in primary window - if (charge > max_charge) { - muon_available = true; - max_charge = charge; - max_cluster = apair.first; - } - } - if (muon_available){ - std::vector& Hits = m_all_clusters->at(max_cluster); - std::vector detkeys = m_all_clusters_detkey->at(max_cluster); - int hits_pmt = 0; - - std::map> hitTimes; - std::map> hitCharges; - - Log("DigitBuilder Tool: Num PMT Clustered Digits = "+to_string(Hits.size()),v_message, verbosity); - for (unsigned int i_hit = 0; i_hit < Hits.size(); i_hit++){ - Hit ahit = Hits.at(i_hit); - unsigned long chankey = detkeys.at(i_hit); - - - if (hitTimes.find(chankey)!=hitTimes.end()){ - hitTimes.at(chankey).push_back(ahit.GetTime()); - hitCharges.at(chankey).push_back(ahit.GetCharge()); - } - else { - std::vector temp_hittimes{ahit.GetTime()}; - std::vector temp_hitcharges{ahit.GetCharge()}; - hitTimes.emplace(chankey,temp_hittimes); - hitCharges.emplace(chankey,temp_hitcharges); - } - } - - if(fParametricModel){ - Log("DigitBuilder tool: Use Parametric Model to create digits",v_message,verbosity); - // Do median and sum - std::map>::iterator it, it2; - for (it=hitTimes.begin(),it2 = hitCharges.begin(); it != hitTimes.end(), it2 != hitCharges.end(); it++, it2++){ - unsigned long chankey = it->first; - std::vector hittimes = it->second; - std::vector hitcharges = it2->second; - det = fGeometry->ChannelToDetector(chankey); - int PMTId = channelkey_to_pmtid.at(chankey); //PMTID In WCSim - if(det==nullptr){ - Log("DigitBuilder Tool: Detector not found! ",v_message,verbosity); - continue; - } - // convert the WCSim coordinates to the ANNIEreco coordinates - // convert the unit from m to cm - pos_sim = det->GetDetectorPosition(); - pos_sim.UnitToCentimeter(); - pos_reco.SetX(pos_sim.X()+xshift); - pos_reco.SetY(pos_sim.Y()+yshift); - pos_reco.SetZ(pos_sim.Z()+zshift); - - std::sort(hittimes.begin(), hittimes.end()); - size_t timesize = hittimes.size(); - if (timesize == 0) continue; - if (timesize % 2 == 0){ - calT = (hittimes.at(timesize/2 - 1) + hittimes.at(timesize/2))/2; - } else { - calT = hittimes.at(timesize/2); - } - calQ = 0.; - for(std::vector::iterator it3 = hitcharges.begin(); it3 != hitcharges.end(); ++it3){ - calQ += *it3; - } - - Log("PMT position (X 0.0){ - calQ_temp = calQ / pmt_gains.at(chankey); - } - if(calQ_temp>fDigitChargeThr) { - digitType = RecoDigit::PMT8inch; - RecoDigit recoDigit(region, pos_reco, calT, calQ_temp, digitType, PMTId); - - //recoDigit.SetHitIDs(hitIDs); - fDigitList->push_back(recoDigit); - } - } - } else { - std::map>::iterator it, it2; - for (it=hitTimes.begin(), it2 = hitCharges.begin(); it != hitTimes.end(), it2 != hitCharges.end(); it++, it2++){ - unsigned long chankey = it->first; - std::vector hittimes = it->second; - std::vector hitcharges = it2->second; - det = fGeometry->ChannelToDetector(chankey); - int PMTId = channelkey_to_pmtid.at(chankey); //PMTID In WCSim - if(det==nullptr){ - Log("DigitBuilder Tool: Detector not found! ",v_message,verbosity); - continue; - } - pos_sim = det->GetDetectorPosition(); - pos_sim.UnitToCentimeter(); - pos_reco.SetX(pos_sim.X()+xshift); - pos_reco.SetY(pos_sim.Y()+yshift); - pos_reco.SetZ(pos_sim.Z()+zshift); - - for (int i=0; i< int(hitcharges.size()); i++){ - calT = hittimes.at(i); - calQ = hitcharges.at(i); - - Log("PMT position (Xpush_back(recoDigit); - } - } - } - } - } - } else { - Log("No Clustered Hits found.",v_warning,verbosity); - return false; - }*/ return true; @@ -794,20 +622,11 @@ void DigitBuilder::CollectMCPMTHits() { for (MCHit& ahit : hits) { - //if (ahit.GetTime() < end_of_window_time_cut * AcqTimeWindow) { + hit_times.push_back(ahit.GetTime()); - //} + } - /*for (MCHit& ahit : hits) { - //std::cout <<"Key: "< 2000.) std::cout <<"Found hit later than 2us! Hit time : "< combine multiple photons if they are within a 10ns range - //hit times can only be recorded with 2ns precision --> possible times are 0ns, 2ns, 4ns, ... - hits_2ns_res.push_back(2 * (int(ahit.GetTime()) / 2.) + (int(ahit.GetTime()) % 2)); - hits_2ns_res_charge.push_back(ahit.GetCharge()); - } - }*/ + if (hit_times.size() == 0) { Log("DigitBuilder tool: no hits in window.",v_message,verbosity); @@ -817,32 +636,30 @@ void DigitBuilder::CollectMCPMTHits() { //Combine multiple MC hits to one pulse std::sort(hit_times.begin(), hit_times.end()); for (int i_hit = 0; i_hit < (int)hit_times.size(); i_hit++) { - Log("check 0: hits, hit_times, datalike_hits, parents size: "+to_string(hits.size())+" "+to_string(hit_times.size())+" "+to_string(datalike_hits.size())+" " + to_string(hits.at(i_hit).GetParents()->size()), v_debug, verbosity); + Log("Hits, hit_times, datalike_hits, parents size: "+to_string(hits.size())+" "+to_string(hit_times.size())+" "+to_string(datalike_hits.size())+" " + to_string(hits.at(i_hit).GetParents()->size()), v_debug, verbosity); double hit1 = hit_times.at(i_hit); if(hit1>end_of_window_time_cut * AcqTimeWindow) continue; int j_hit=0; if (datalike_hits.size() == 0) { - Log("check 0-b", v_debug, verbosity); + datalike_hits.push_back(hit1); first_time=hit1; datalike_hits_charge.push_back(hits.at(i_hit).GetCharge()); - Log("check 0-c", v_debug, verbosity); if (hits.at(i_hit).GetParents()->size() == 0) { Parents_by_hit.push_back(-5); Log("found hit with no parent at time "+to_string(hits.at(i_hit).GetTime()),v_debug,verbosity); } else { Parents_by_hit.push_back(hits.at(i_hit).GetParents()->at(0)); - Log("check 0-d: hit parent " + to_string(hits.at(i_hit).GetParents()->at(0)) + " and time: " + to_string(hits.at(i_hit).GetTime()), v_debug, verbosity); } temp_times.push_back(hit1); } else { bool new_pulse = false; - Log("check 0-a", v_debug, verbosity); + for (int j_hit = 0; j_hit < (int)datalike_hits.size(); j_hit++) { if (fabs(first_time - hit1) < 10.) { @@ -856,7 +673,7 @@ void DigitBuilder::CollectMCPMTHits() { temp_times.push_back(hit1); } - Log("check 1-"+to_string(j_hit)+" new_pulse: "+to_string(new_pulse), v_debug, verbosity); + } if (new_pulse) { @@ -870,7 +687,7 @@ void DigitBuilder::CollectMCPMTHits() { else { mid_time = temp_times.at(temp_times.size() / 2); } - Log("check 1-a", v_debug, verbosity); + datalike_hits.at(j_hit)=mid_time; //datalike_hits.push_back(mid_time); //Only count as a new pulse if it was 10ns away from every other pulse datalike_hits_charge.push_back(hits.at(i_hit).GetCharge()); j_hit++; @@ -880,7 +697,7 @@ void DigitBuilder::CollectMCPMTHits() { if (hits.at(i_hit).GetParents()->size() == 0) Parents_by_hit.push_back(-5); else { Parents_by_hit.push_back(hits.at(i_hit).GetParents()->at(0)); - Log("check 0-e: hit parent " + to_string(hits.at(i_hit).GetParents()->at(0)) + " and time: " + to_string(hits.at(i_hit).GetTime()), v_debug, verbosity); + Log("Hit parent " + to_string(hits.at(i_hit).GetParents()->at(0)) + " and time: " + to_string(hits.at(i_hit).GetTime()), v_debug, verbosity); } temp_times.clear(); @@ -888,19 +705,19 @@ void DigitBuilder::CollectMCPMTHits() { } } - Log("check 2", v_debug, verbosity); + //hits.clear(); for (int i_hit = 0; i_hit < (int)datalike_hits.size(); i_hit++) { - Log("check 1", v_debug, verbosity); + calT = datalike_hits.at(i_hit); - Log("check 2", v_debug, verbosity); + if (MCPMTResolution > 0) calT= frand.Gaus(calT, MCPMTResolution); calQ = datalike_hits_charge.at(i_hit); - Log("check 3",v_debug,verbosity); + RecoDigit recoDigit(-999 /*region*/, pos_reco, calT, calQ, RecoDigit::PMT8inch, PMTId); temp_parents.push_back(Parents_by_hit.at(i_hit)); - Log("check 4",v_debug,verbosity); + recoDigit.SetParents(temp_parents); fDigitList->push_back(recoDigit); temp_parents.clear(); diff --git a/configfiles/ClusterSearcher/ClusterSearcher2Config b/configfiles/ClusterSearcher/ClusterSearcher2Config index 20672cf1c..d49a09cc9 100755 --- a/configfiles/ClusterSearcher/ClusterSearcher2Config +++ b/configfiles/ClusterSearcher/ClusterSearcher2Config @@ -2,7 +2,8 @@ verbosity 5 ClusterMode 2 #searching track-like clusters -Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +IsMC 1 #whether the input data is simulation or data +SinglePEGains /Path/to/Gains/File #path to the s.p.e. gains file; only needed for running over data PmtMinPulseHeight 20 #minimum pulse height PmtNeighbourRadius 600 #digit neighbouring distance [cm] PmtMinNeighbourDigits 2 #minimum neighbour digits diff --git a/configfiles/ClusterSearcher/ClusterSearcherConfig b/configfiles/ClusterSearcher/ClusterSearcherConfig index c135e4b6f..10b2eb0bf 100755 --- a/configfiles/ClusterSearcher/ClusterSearcherConfig +++ b/configfiles/ClusterSearcher/ClusterSearcherConfig @@ -2,7 +2,8 @@ verbosity 5 ClusterMode 1 #searching track-like clusters -Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +IsMC 1 #whether the input data is simulation or data +SinglePEGains /Path/to/Gains/File #path to the s.p.e. gains file; only needed for running over data PmtMinPulseHeight 20 #minimum pulse height PmtNeighbourRadius 60 #digit neighbouring distance [cm] PmtMinNeighbourDigits 2 #minimum neighbour digits From 8a6499f8a4e3747baed76593ddfd0d67b0fb32b5 Mon Sep 17 00:00:00 2001 From: root Date: Fri, 23 Jan 2026 18:51:24 +0000 Subject: [PATCH 3/5] Removed ClusterCA from NeutronCheck to fix bug. --- UserTools/NeutronCheck/NeutronCheck.cpp | 4 +--- UserTools/NeutronCheck/NeutronCheck.h | 1 - 2 files changed, 1 insertion(+), 4 deletions(-) mode change 100644 => 100755 UserTools/NeutronCheck/NeutronCheck.cpp mode change 100644 => 100755 UserTools/NeutronCheck/NeutronCheck.h diff --git a/UserTools/NeutronCheck/NeutronCheck.cpp b/UserTools/NeutronCheck/NeutronCheck.cpp old mode 100644 new mode 100755 index 50391f1af..443ad40cf --- a/UserTools/NeutronCheck/NeutronCheck.cpp +++ b/UserTools/NeutronCheck/NeutronCheck.cpp @@ -55,7 +55,6 @@ bool NeutronCheck::Initialise(std::string configfile, DataModel &data){ NeutCheckTree->Branch("ClusterCVZ", &fClusterCVZ); NeutCheckTree->Branch("ClusterCVR", &fClusterCVR); NeutCheckTree->Branch("ClusterAMD", &fClusterAMD); - NeutCheckTree->Branch("ClusterCA", &fClusterCA); if(fFinderCompare){ @@ -423,7 +422,7 @@ void NeutronCheck::CSCheck() { fClusterAS0.push_back(fRecoClusters->at(i)->GetAS(0)); fClusterAS1.push_back(fRecoClusters->at(i)->GetAS(1)); fClusterAMD.push_back(fRecoClusters->at(i)->GetAMD()); - fClusterCA.push_back(fRecoClusters->at(i)->GetCA()); + if(fRecoClusters->at(i)->GetAS(0)>0){ //AS2RecoCheck->Fill(fRecoClusters->at(i)->GetAS(1)/fRecoClusters->at(i)->GetAS(0)); fClusterAS2.push_back(fClusterAS1.at(fClusterAS1.size()-1)/fClusterAS0.at(fClusterAS0.size()-1)); @@ -467,7 +466,6 @@ void NeutronCheck::ResetVariables() { fClusterCVZ.clear(); fClusterCVR.clear(); fClusterAMD.clear(); - fClusterCA.clear(); fNeutronMult=0; fMCNeutCapTimes.clear(); fMCNeutCapX.clear(); diff --git a/UserTools/NeutronCheck/NeutronCheck.h b/UserTools/NeutronCheck/NeutronCheck.h old mode 100644 new mode 100755 index c446ee6de..16a01af72 --- a/UserTools/NeutronCheck/NeutronCheck.h +++ b/UserTools/NeutronCheck/NeutronCheck.h @@ -118,7 +118,6 @@ class NeutronCheck: public Tool { vector fClusterAS2; vector fClusterASC; vector fClusterAMD; - vector fClusterCA; //vector fClusterSA; //vector fClusterCV; vector fClusterCVX,fClusterCVY,fClusterCVZ,fClusterCVR; From 0236d079b9cf5bf864fc3b448c5be756b64c85a0 Mon Sep 17 00:00:00 2001 From: root Date: Fri, 23 Jan 2026 19:34:55 +0000 Subject: [PATCH 4/5] Modified HitCleaner for non-pointer Digits list. --- UserTools/HitCleaner/HitCleaner.cpp | 49 +++++++++++++++++++++++------ UserTools/HitCleaner/HitCleaner.h | 7 +++-- 2 files changed, 44 insertions(+), 12 deletions(-) diff --git a/UserTools/HitCleaner/HitCleaner.cpp b/UserTools/HitCleaner/HitCleaner.cpp index 3475b4ce2..4adc803d9 100644 --- a/UserTools/HitCleaner/HitCleaner.cpp +++ b/UserTools/HitCleaner/HitCleaner.cpp @@ -118,8 +118,8 @@ bool HitCleaner::Initialise(std::string configfile, DataModel &data){ // only for test fFilterByTruthInfo = new std::vector; // vector of clusters - fClusterList = new std::vector; - fHitCleaningClusters = new std::vector; + fClusterList = new std::vector; + fHitCleaningClusters = new std::vector; //Set hit cleaner parameters in the RecoEvent store m_data->Stores.at("RecoEvent")->Set("HitCleaningParameters", fHitCleaningParam); @@ -194,6 +194,7 @@ bool HitCleaner::Execute(){ fIsHitCleaningDone = true; m_data->Stores.at("RecoEvent")->Set("HitCleaningDone", fIsHitCleaningDone); m_data->Stores.at("RecoEvent")->Set("HitCleaningClusters", fHitCleaningClusters); + CBCheck(digits,FilterDigitList); delete digits; digits = 0; return true; @@ -523,14 +524,15 @@ std::vector* HitCleaner::FilterByClusters(std::vector* m // run clustering algorithm // ======================== - std::vector* myClusterList = (std::vector*)(this->RecoClusters(myDigitList)); + std::vector* myClusterList = (std::vector*)(this->RecoClusters(myDigitList)); for(int icluster=0; iclustersize()); icluster++ ){ - RecoCluster* myCluster = (RecoCluster*)(myClusterList->at(icluster)); + RecoCluster myCluster = (myClusterList->at(icluster)); fHitCleaningClusters->push_back(myCluster); - for(int idigit=0; idigitGetNDigits(); idigit++ ){ - RecoDigit* myDigit = (RecoDigit*)(myCluster->GetDigit(idigit)); + for(int idigit=0; idigitpush_back(myDigit); } } @@ -543,7 +545,7 @@ std::vector* HitCleaner::FilterByClusters(std::vector* m return fFilterByClusters; } -std::vector* HitCleaner::RecoClusters(std::vector* myDigitList) +std::vector* HitCleaner::RecoClusters(std::vector* myDigitList) { // delete cluster digits @@ -679,13 +681,13 @@ std::vector* HitCleaner::RecoClusters(std::vector* myD } //std::cout <<"vClusterDigitCollection.size() == "<=fMinClusterDigits ){ - RecoCluster* cluster = new RecoCluster(); + RecoCluster cluster; fClusterList->push_back(cluster); for(int jdigit=0; jdigitGetRecoDigit()); - cluster->AddDigit(recodigit); + cluster.AddDigit(recodigit); } } } @@ -738,3 +740,32 @@ std::vector* HitCleaner::FilterByTruthInfo(std::vector* return fFilterByTruthInfo; } + +void HitCleaner::CBCheck(std::vector* unfilteredDigits, std::vector* filteredDigits) { + //calculate unfiltered CB + double total_Q = 0; + double total_QSquared = 0; + for (int i=0;isize();i++) { + //if(unfilteredDigits->at(i)->GetDigitType()==RecoDigit::PMT8inch){ + double tube_charge = unfilteredDigits->at(i)->GetCalCharge(); + total_Q += tube_charge; + total_QSquared += (tube_charge * tube_charge); + //} + } + //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.)); + if (verbosity > 4) std::cout << "HitCleaner Tool: Unfiltered CB: " << ucharge_balance << std::endl; + + total_Q = 0; + total_QSquared = 0; + for (int i = 0; i < filteredDigits->size(); i++) { + //if (filteredDigits->at(i)->GetDigitType() == RecoDigit::PMT8inch) { + double tube_charge = filteredDigits->at(i)->GetCalCharge(); + total_Q += tube_charge; + total_QSquared += (tube_charge * tube_charge); + //} + } + double fcharge_balance = sqrt((total_QSquared) / (total_Q * total_Q) - (1. / 123.)); + if (verbosity > 4) std::cout << "HitCleaner Tool: filtered CB: " << fcharge_balance << std::endl; + +} diff --git a/UserTools/HitCleaner/HitCleaner.h b/UserTools/HitCleaner/HitCleaner.h index e20cef7a0..019459d76 100644 --- a/UserTools/HitCleaner/HitCleaner.h +++ b/UserTools/HitCleaner/HitCleaner.h @@ -67,11 +67,12 @@ class HitCleaner: public Tool { std::vector* FilterByNeighbours(std::vector* digitlist); std::vector* FilterByClusters(std::vector* digitlist); std::vector* FilterByTruthInfo(std::vector* digitlist); //use truth information. Only for testing the code - std::vector* RecoClusters(std::vector* digitlist); + std::vector* RecoClusters(std::vector* digitlist); private: void Reset(); + void CBCheck(std::vector* unfilteredDigits, std::vector* filteredDigits); // running mode int fConfig; @@ -120,10 +121,10 @@ class HitCleaner: public Tool { std::vector* fFilterByTruthInfo; // vectors of clusters - std::vector* fClusterList; + std::vector* fClusterList; // vector of clusters (accessible to the CStore) - std::vector* fHitCleaningClusters = nullptr; + std::vector* fHitCleaningClusters = nullptr; // true vertex RecoVertex* fTrueVertex = 0; From efe448ab077af5a89e1a4f9bdd28302da025ff71 Mon Sep 17 00:00:00 2001 From: root Date: Mon, 26 Jan 2026 18:14:44 +0000 Subject: [PATCH 5/5] Update/fix to ClusterSearcher sample toolchain. --- .../ClusterSearcher/ClusterSearcher2Config | 4 ++-- .../ClusterSearcher/ClusterSearcherConfig | 4 ++-- .../ClusterSearcher/EventSelectorConfig | 18 ++++++++++++++++++ .../ClusterSearcher/LoadGenieEventConfig | 12 ++++++++++++ configfiles/ClusterSearcher/LoadGeometryConfig | 3 +-- .../ClusterSearcher/MCRecoEventLoaderConfig | 11 +++++++++++ configfiles/ClusterSearcher/NeutronCheckConfig | 3 ++- configfiles/ClusterSearcher/ToolsConfig | 4 +++- 8 files changed, 51 insertions(+), 8 deletions(-) create mode 100755 configfiles/ClusterSearcher/EventSelectorConfig create mode 100644 configfiles/ClusterSearcher/LoadGenieEventConfig create mode 100644 configfiles/ClusterSearcher/MCRecoEventLoaderConfig diff --git a/configfiles/ClusterSearcher/ClusterSearcher2Config b/configfiles/ClusterSearcher/ClusterSearcher2Config index d49a09cc9..20af28b1b 100755 --- a/configfiles/ClusterSearcher/ClusterSearcher2Config +++ b/configfiles/ClusterSearcher/ClusterSearcher2Config @@ -2,8 +2,8 @@ verbosity 5 ClusterMode 2 #searching track-like clusters -IsMC 1 #whether the input data is simulation or data -SinglePEGains /Path/to/Gains/File #path to the s.p.e. gains file; only needed for running over data +Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +IsMC 1 PmtMinPulseHeight 20 #minimum pulse height PmtNeighbourRadius 600 #digit neighbouring distance [cm] PmtMinNeighbourDigits 2 #minimum neighbour digits diff --git a/configfiles/ClusterSearcher/ClusterSearcherConfig b/configfiles/ClusterSearcher/ClusterSearcherConfig index 10b2eb0bf..22020dbbf 100755 --- a/configfiles/ClusterSearcher/ClusterSearcherConfig +++ b/configfiles/ClusterSearcher/ClusterSearcherConfig @@ -2,8 +2,8 @@ verbosity 5 ClusterMode 1 #searching track-like clusters -IsMC 1 #whether the input data is simulation or data -SinglePEGains /Path/to/Gains/File #path to the s.p.e. gains file; only needed for running over data +Config 3 #config type: 1 = pulse height cut, 2 = +neighbour cut, 3 = +cluster cut +IsMC 1 PmtMinPulseHeight 20 #minimum pulse height PmtNeighbourRadius 60 #digit neighbouring distance [cm] PmtMinNeighbourDigits 2 #minimum neighbour digits diff --git a/configfiles/ClusterSearcher/EventSelectorConfig b/configfiles/ClusterSearcher/EventSelectorConfig new file mode 100755 index 000000000..6fd425abd --- /dev/null +++ b/configfiles/ClusterSearcher/EventSelectorConfig @@ -0,0 +1,18 @@ +# EventSelector config file + +verbosity 0 +MCPMTVolCut 0 +MCFVCut 1 +MCMRDCut 1 +MCPiKCut 0 +MCIsMuonCut 1 +MRDRecoCut 0 +RecoPMTVolCut 0 +RecoFVCut 0 +ArgonFV 0 +MCIsSingleRingCut 0 +NHitCut 1 +LAPPDMult 0 +TriggerWord 0 +PromptTrigOnly 0 +SaveStatusToStore 1 diff --git a/configfiles/ClusterSearcher/LoadGenieEventConfig b/configfiles/ClusterSearcher/LoadGenieEventConfig new file mode 100644 index 000000000..26c9051d8 --- /dev/null +++ b/configfiles/ClusterSearcher/LoadGenieEventConfig @@ -0,0 +1,12 @@ +verbosity 0 +FluxVersion 0 # use 0 to load genie files based on bnb_annie_0000.root etc files + # use 1 to load files based on beammc_annie_0000.root etc files +#FileDir NA # specify "NA" for newer files: full path is saved in WCSim +#FileDir /pnfs/annie/persistent/users/vfischer/genie_files/BNB_Water_10k_22-05-17 +#FileDir /pnfs/annie/persistent/users/moflaher/genie/BNB_World_10k_11-03-18_gsimpleflux +FileDir /Data/Neutronsim +FilePattern gntp.1078.ghep.root +#FilePattern LoadWCSimTool ## use this pattern to load corresponding genie info with the LoadWCSimTool + ## N.B: FileDir must still be specified for now! (WCSim files do not record their directory) +ManualFileMatching 0 +EventOffset 0 #9900 diff --git a/configfiles/ClusterSearcher/LoadGeometryConfig b/configfiles/ClusterSearcher/LoadGeometryConfig index 79efc3811..3d8ce96bc 100644 --- a/configfiles/ClusterSearcher/LoadGeometryConfig +++ b/configfiles/ClusterSearcher/LoadGeometryConfig @@ -5,6 +5,5 @@ FACCMRDGeoFile ./configfiles/LoadGeometry/FullMRDGeometry_09_29_20.csv DetectorGeoFile ./configfiles/LoadGeometry/DetectorGeometrySpecs.csv LAPPDGeoFile ./configfiles/LoadGeometry/LAPPDGeometry.csv TankPMTGeoFile ./configfiles/LoadGeometry/FullTankPMTGeometry.csv -TankPMTGainFile ./configfiles/LoadGeometry/ChannelSPEGains2023.csv -TankPMTTimingOffsetFile ./configfiles/LoadGeometry/TankPMTTimingOffsets.csv +TankPMTGainFile ./configfiles/LoadGeometry/ChannelSPEGains_BeamRun20192020.csv AuxiliaryChannelFile ./configfiles/LoadGeometry/AuxChannels.csv diff --git a/configfiles/ClusterSearcher/MCRecoEventLoaderConfig b/configfiles/ClusterSearcher/MCRecoEventLoaderConfig new file mode 100644 index 000000000..cba5524e7 --- /dev/null +++ b/configfiles/ClusterSearcher/MCRecoEventLoaderConfig @@ -0,0 +1,11 @@ +# MCRecoEventLoader config file + +verbosity 5 +GetPionKaonInfo 1 +GetNRings 1 +ParticleID 13 +DoParticleSelection 0 +xshift 0.0 +yshift 14.46469 +zshift -168.1 +IsMC 1 diff --git a/configfiles/ClusterSearcher/NeutronCheckConfig b/configfiles/ClusterSearcher/NeutronCheckConfig index ccb5d89a3..0658b1a35 100644 --- a/configfiles/ClusterSearcher/NeutronCheckConfig +++ b/configfiles/ClusterSearcher/NeutronCheckConfig @@ -2,4 +2,5 @@ verbosity 5 outfile /Data/0824Beamsim/NeutronCheck_CS UseClean 1 FinderCompare 0 -ParticleInfo 0 \ No newline at end of file +ParticleInfo 0 +VertexInfo 0 diff --git a/configfiles/ClusterSearcher/ToolsConfig b/configfiles/ClusterSearcher/ToolsConfig index adff3b930..7da2cdb68 100755 --- a/configfiles/ClusterSearcher/ToolsConfig +++ b/configfiles/ClusterSearcher/ToolsConfig @@ -1,9 +1,11 @@ LoadGeometry LoadGeometry ./configfiles/ClusterSearcher/LoadGeometryConfig LoadWCSim LoadWCSim ./configfiles/ClusterSearcher/LoadWCSimConfig LoadWCSimLAPPD LoadWCSimLAPPD ./configfiles/ClusterSearcher/LoadWCSimLAPPDConfig +LoadGenieEvent LoadGenieEvent ./configfiles/ClusterSearcher/LoadGenieEventConfig +MCRecoEventLoader MCRecoEventLoader ./configfiles/ClusterSearcher/MCRecoEventLoaderConfig MCParticleProperties MCParticleProperties ./configfiles/ClusterSearcher/MCParticlePropertiesConfig DigitBuilder DigitBuilder ./configfiles/ClusterSearcher/DigitBuilderConfig ClusterSearcher ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcherConfig -ClusterSearcher2 ClusterSearcher2 ./configfiles/ClusterSearcher/ClusterSearcher2Config +ClusterSearcher2 ClusterSearcher ./configfiles/ClusterSearcher/ClusterSearcher2Config EventSelector EventSelector ./configfiles/ClusterSearcher/EventSelectorConfig NeutronCheck NeutronCheck ./configfiles/ClusterSearcher/NeutronCheckConfig