Skip to content

Commit b807420

Browse files
asanchezcastilloasanchezcastillo
committed
Cleanup
1 parent 3414827 commit b807420

1 file changed

Lines changed: 10 additions & 41 deletions

File tree

‎sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc‎

Lines changed: 10 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -240,7 +240,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer {
240240
double fZError;
241241
double fAngleError;
242242
std::vector<double> fFlashVetoWindow;
243-
int fDebugEvent;
244243
// Event-level data members
245244
int fRun; ///< Number of the run being processed
246245
int fEvent; ///< Number of the event being processed
@@ -270,7 +269,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer {
270269
double fDeltaY_Trigger; ///< | Triggering flash Y center - charge Y center | (cm)
271270
double fDeltaZ_Trigger; ///< | Triggering flash Z center - charge Z center | (cm)
272271
double fRadius_Trigger; ///< Hypotenuse of DeltaY_Trigger and DeltaZ_Trigger (cm)
273-
double fNuScore; ///< NuScore of the slice, if available, otherwise -9999
274272
TTree* fMatchTree; ///< Tree to store all match information
275273
// Geometry service
276274
geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle<geo::WireReadout>()->Get();
@@ -310,7 +308,6 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet
310308
fZError(p.get<double>("ZError")), // cm
311309
fAngleError(p.get<double>("AngleError")), // deg
312310
fFlashVetoWindow(p.get<std::vector<double>>("FlashVetoWindow")), // us
313-
fDebugEvent(p.get<int>("DebugEvent", -1)) // Event to debug, -1 means no debug event
314311
{
315312
// Call appropriate produces<>() functions here.
316313

@@ -343,7 +340,6 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet
343340
fMatchTree->Branch("chargeWidthX", &fChargeWidthX, "chargeWidthX/d" );
344341
fMatchTree->Branch("chargeWidthY", &fChargeWidthY, "chargeWidthY/d" );
345342
fMatchTree->Branch("chargeWidthZ", &fChargeWidthZ, "chargeWidthZ/d" );
346-
fMatchTree->Branch("nuScore", &fNuScore, "nuScore/d" );
347343

348344
//Matched Flash Info
349345
fMatchTree->Branch("flashTime", &fFlashTime, "flashTime/d" );
@@ -442,7 +438,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
442438
unsigned nSlices = (*sliceHandle).size();
443439
::art::Handle<std::vector<recob::PFParticle>> pfpHandle;
444440
e.getByLabel(fPandoraLabel, pfpHandle);
445-
//unsigned nPFPs = (*pfpHandle).size();
446441
art::FindManyP<larpandoraobj::PFParticleMetadata> pfp_to_metadata(pfpHandle, e, fPandoraLabel);
447442
art::FindManyP<recob::PFParticle> slice_pfp_assns (sliceHandle, e, fPandoraLabel);
448443

@@ -465,27 +460,17 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
465460
const art::Ptr<recob::Slice> slicePtr { sliceHandle, j };
466461
// -----------------> This is to get the nuscore of the slice, only for debugging purposes <--------------------- //
467462
//Slice to PFParticles association
468-
fNuScore = -9999.; //Default value
469463
//Vector for recob PFParticles
470-
471464
std::vector<art::Ptr<recob::Hit>> tpcHitsVec;
472465
std::vector<art::Ptr<recob::PFParticle>> pfpVect = slice_pfp_assns.at(j);
473466

474-
// PFP Metadata
467+
// Get the hits associated to the PFParticles in this slice
475468
for(const art::Ptr<recob::PFParticle> &pfp : pfpVect){
476469
std::vector<art::Ptr<recob::Cluster>> cluster_v = fmClusterPfp.at(pfp.key());
477470
for(size_t i=0; i<cluster_v.size(); i++){
478471
std::vector<art::Ptr<recob::Hit>> hitVect = cluster_hit_assns.at(cluster_v[i].key());
479472
tpcHitsVec.insert(tpcHitsVec.end(), hitVect.begin(), hitVect.end());
480473
}
481-
if(pfp->IsPrimary() && ( std::abs(pfp->PdgCode())==12 || std::abs(pfp->PdgCode())==14 ) )
482-
{
483-
const std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> pfpMetaVec = pfp_to_metadata.at(pfp.key());
484-
for (auto const pfpMeta : pfpMetaVec) {
485-
larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap();
486-
fNuScore = propertiesMap.at("NuScore");
487-
}
488-
}
489474
}
490475

491476
for ( size_t tpc=0; tpc<2 ; tpc ++) {
@@ -591,16 +576,14 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
591576
e.getByLabel(fOpFlashesModuleLabel[tpc], flashHandle);
592577
std::vector< art::Ptr<recob::OpFlash> > flashVect;
593578
art::fill_ptr_vector(flashVect, flashHandle);
579+
594580
//OpHit OpFlash assns
595581
art::FindManyP<recob::OpHit> flash_ophit_assns(flashHandle, e, fOpFlashesModuleLabel[tpc]);
596582

597583
//Vector to store the idxs of the candidate flashes
598584
std::vector<int> candidateFlashIdxs;
599-
600585
int nFlashes = flashVect.size();
601586

602-
//Vector to store the distance between the barycenter of the charge and the barycenter of the light
603-
std::vector<double> chargeLightDistance;
604587
//For flash...
605588
for ( int m = 0; m < nFlashes; m++ ) {
606589
auto & flash = flashVect[m];
@@ -612,7 +595,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
612595
if(fDo3DMatching)
613596
thisDistance = std::hypot( (thisFlashCenterX - fChargeCenterX), (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) );
614597
else thisDistance = std::hypot( (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) );
615-
chargeLightDistance.push_back(thisDistance);
598+
616599
if ( thisDistance < minDistance ) {
617600
minDistance = thisDistance;
618601
matchIndex = m;
@@ -621,10 +604,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
621604
if(thisDistance<fDistanceCandidateFlashes){
622605
candidateFlashIdxs.push_back(m);
623606
}
624-
if(fEvent == fDebugEvent )
625-
{
626-
std::cout << " Flash " << m << " at time " << flash->Time() <<" has barycenter " << thisFlashCenterX << ", " << thisFlashCenterY << ", " << thisFlashCenterZ << " and distance " << thisDistance << " and PE " << flash->TotalPE() << std::endl;
627-
}
628607
} //End for flash
629608

630609
//No valid match found...
@@ -664,16 +643,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
664643
else
665644
chi2 = std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2);
666645

667-
if(fNuScore>0.55)
668-
{
669-
std::cout << " flash " << i << " at time " << flash->Time() <<" has barycenter " << _currentFlashX << ", " << _currentFlashY << ", " << _currentFlashZ << " and distance " << chargeLightDistance[idx] << " with angle " << angle << " and charge light ratio " << lightChargeRatio << " and PE " << flash->TotalPE() <<std::endl;
670-
std::cout << " chi2 " << chi2 << std::endl;
671-
}
672-
673-
/*if(fNuScore>0.55)
674-
{
675-
std::cout << " --------> Candidate flash " << idx << " at time " << flash->Time() << " has LightChargeRatio " << lightChargeRatio <<" distance " << chargeLightDistance[idx] << " angle " << angle << " and chi2 " << chi2 << std::endl;
676-
}*/
677646
if(lightChargeRatio < fLightChargeRatioBounds[0] || lightChargeRatio > fLightChargeRatioBounds[1]) continue;
678647
if(chi2 < minChi2)
679648
{
@@ -682,9 +651,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
682651
}
683652
}
684653

685-
//std::cout << " Chosen flash is " << matchIndex << std::endl;
686-
//Now we have a pool of flashes so we can have a look at the PCAs and the calorimetry
687-
688654
//Best match flash pointer
689655
unsigned unsignedMatchIndex = matchIndex;
690656
const art::Ptr<recob::OpFlash> flashPtr { flashHandle, unsignedMatchIndex };
@@ -717,7 +683,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e)
717683
flashAssns->addSingle(infoPtrVector[maxFlashIdx], flashPtrVector[maxFlashIdx]);
718684
matchInfoVector->push_back(std::move(sliceMatchInfoVector[maxFlashIdx]));
719685
}
720-
if(fNuScore>0.55) std::cout << " -----------------------> Selecting flash at time " << flashPtrVector[maxFlashIdx]->Time() << std::endl;
721686
} //End for slice
722687

723688
//Store new products at the end of the event
@@ -816,28 +781,32 @@ void TPCPMTBarycenterMatchProducer::GetPCA(std::vector<double> const&x, std::vec
816781
points(i, 0) = x[i];
817782
points(i, 1) = y[i];
818783
}
784+
819785
//Get the weighted centroid
820786
Eigen::Vector2d weighted_sum(0.0, 0.0);
821787
for (size_t i = 0; i < n; ++i) {
822788
weighted_sum(0) += weight[i] * x[i];
823789
weighted_sum(1) += weight[i] * y[i];
824790
}
825-
Eigen::Vector2d centroid = weighted_sum / weight_sum;
791+
826792
//Center the points wrt centroid
793+
Eigen::Vector2d centroid = weighted_sum / weight_sum;
827794
Eigen::MatrixXd centered_points = points.rowwise() - centroid.transpose();
795+
828796
//Get the weighted covariance matrix
829797
Eigen::Matrix2d cov = Eigen::Matrix2d::Zero();
830798
for (size_t i = 0; i < n; ++i) {
831799
Eigen::Vector2d pt = centered_points.row(i);
832800
cov += std::abs(weight[i]) * (pt * pt.transpose());
833801
}
834802
cov /= weight_sum;
835-
// Obtener autovalores y autovectores
803+
804+
// Get eigenvalues / eigenvectors
836805
Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> solver(cov);
837806
Eigen::VectorXd eigvals = solver.eigenvalues();
838807
Eigen::Matrix2d eigvecs = solver.eigenvectors();
839808

840-
// Dirección principal (mayor autovalor)
809+
// Get Principal component (largest eigenvalue)
841810
int maxIndex;
842811
eigvals.maxCoeff(&maxIndex);
843812
Eigen::Vector2d principal_direction = eigvecs.col(maxIndex);

0 commit comments

Comments
 (0)