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