From e7aadd7d75952d7e93b794789c06bfa1ac770d99 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 16 Jan 2024 10:45:17 -0500 Subject: [PATCH 01/13] Add matching to FMWPC clusters to DDetectorMatches --- src/libraries/PID/DChargedTrackHypothesis.h | 6 +- .../PID/DChargedTrackHypothesis_factory.cc | 15 ++++ src/libraries/PID/DDetectorMatches.h | 53 ++++++++++-- src/libraries/PID/DDetectorMatches_factory.cc | 80 ++++++++++++++++++- src/libraries/PID/DDetectorMatches_factory.h | 2 + 5 files changed, 147 insertions(+), 9 deletions(-) diff --git a/src/libraries/PID/DChargedTrackHypothesis.h b/src/libraries/PID/DChargedTrackHypothesis.h index 4fba5ecef..74e33c72a 100644 --- a/src/libraries/PID/DChargedTrackHypothesis.h +++ b/src/libraries/PID/DChargedTrackHypothesis.h @@ -69,6 +69,7 @@ class DChargedTrackHypothesis : public DKinematicData shared_ptr Get_FCALSingleHitMatchParams(void) const{return dTrackingInfo->dFCALSingleHitMatchParams;} shared_ptr Get_DIRCMatchParams(void) const{return dTrackingInfo->dDIRCMatchParams;} shared_ptr Get_CTOFHitMatchParams(void) const{return dTrackingInfo->dCTOFHitMatchParams;} + shared_ptr Get_FMWPCMatchParams(void) const{return dTrackingInfo->dFMWPCMatchParams;} //SETTERS @@ -93,6 +94,7 @@ class DChargedTrackHypothesis : public DKinematicData void Set_FCALSingleHitMatchParams(shared_ptr locMatchParams){dTrackingInfo->dFCALSingleHitMatchParams = locMatchParams;} void Set_DIRCMatchParams(shared_ptr locMatchParams){dTrackingInfo->dDIRCMatchParams = locMatchParams;} void Set_CTOFHitMatchParams(shared_ptr locMatchParams){dTrackingInfo->dCTOFHitMatchParams = locMatchParams;} + void Set_FMWPCMatchParams(shared_ptr locMatchParams){dTrackingInfo->dFMWPCMatchParams = locMatchParams;} void toStrings(vector > &items) const { @@ -171,6 +173,7 @@ class DChargedTrackHypothesis : public DKinematicData shared_ptr dFCALSingleHitMatchParams = nullptr; shared_ptr dDIRCMatchParams = nullptr; shared_ptr dCTOFHitMatchParams = nullptr; + shared_ptr dFMWPCMatchParams=nullptr; }; private: @@ -405,11 +408,12 @@ inline void DChargedTrackHypothesis::DTrackingInfo::Reset(void) dTrackTimeBased = nullptr; dSCHitMatchParams = nullptr; dTOFHitMatchParams = nullptr; - dCTOFHitMatchParams = nullptr; + dCTOFHitMatchParams = nullptr; dBCALShowerMatchParams = nullptr; dFCALShowerMatchParams = nullptr; dFCALSingleHitMatchParams = nullptr; dDIRCMatchParams = nullptr; + dFMWPCMatchParams = nullptr; } inline void DChargedTrackHypothesis::DEOverPInfo::Reset(void) diff --git a/src/libraries/PID/DChargedTrackHypothesis_factory.cc b/src/libraries/PID/DChargedTrackHypothesis_factory.cc index 5cdeb8ae2..f226ba2d2 100644 --- a/src/libraries/PID/DChargedTrackHypothesis_factory.cc +++ b/src/libraries/PID/DChargedTrackHypothesis_factory.cc @@ -283,6 +283,21 @@ DChargedTrackHypothesis* DChargedTrackHypothesis_factory::Create_ChargedTrackHyp locChargedTrackHypothesis->Set_DIRCMatchParams(locDIRCMatchParams); if(dPIDAlgorithm->Get_BestCTOFMatchParams(locTrackTimeBased, locDetectorMatches, locCTOFHitMatchParams)) locChargedTrackHypothesis->Set_CTOFHitMatchParams(locCTOFHitMatchParams); + // Matching to CPP wire chambers + vector > locFMWPCMatchParamsVec; + if (locDetectorMatches->Get_FMWPCMatchParams(locTrackTimeBased, + locFMWPCMatchParamsVec)){ + locChargedTrackHypothesis->Set_FMWPCMatchParams(locFMWPCMatchParamsVec[0]); + shared_ptr myshared=locChargedTrackHypothesis->Get_FMWPCMatchParams(); + if (myshared){ + for (int k=0;kdLayers.size();k++){ + cout << myshared->dLayers[k]<<": " <dNhits[k] + << " " << myshared->dDists[k] + <<" " << myshared->dClosestWires[k]<t1_detector() == SYS_BCAL) diff --git a/src/libraries/PID/DDetectorMatches.h b/src/libraries/PID/DDetectorMatches.h index 68607d610..63f38b666 100644 --- a/src/libraries/PID/DDetectorMatches.h +++ b/src/libraries/PID/DDetectorMatches.h @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -155,6 +156,15 @@ class DDIRCMatchParams double dExtrapolatedTime; }; +class DFMWPCMatchParams +{ + public: + vectordLayers; + vectordNhits; + vectordDists; + vectordClosestWires; +}; + class DDetectorMatches : public JObject { public: @@ -166,9 +176,10 @@ class DDetectorMatches : public JObject inline bool Get_FCALMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const; inline bool Get_TOFMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const; inline bool Get_CTOFMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const; + inline bool Get_FMWPCMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const; inline bool Get_SCMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const; inline bool Get_DIRCMatchParams(const DTrackingData* locTrack, shared_ptr& locMatchParams) const; - + inline bool Get_IsMatchedToTrack(const DBCALShower* locBCALShower) const; inline bool Get_IsMatchedToTrack(const DFCALShower* locFCALShower) const; inline bool Get_IsMatchedToHit(const DTrackingData* locTrack) const; @@ -179,7 +190,7 @@ class DDetectorMatches : public JObject inline bool Get_TrackMatchParams(const DTOFPoint* locTOFPoint, vector >& locMatchParams) const; inline bool Get_TrackMatchParams(const DSCHit* locSCHit, vector >& locMatchParams) const; inline bool Get_DIRCTrackMatchParamsMap(map, vector >& locDIRCTrackMatchParamsMap); - inline bool Get_TrackMatchParams(const DCTOFPoint* locCTOFPoint, vector >& locMatchParams) const; + inline bool Get_TrackMatchParams(const DCTOFPoint* locCTOFPoint, vector >& locMatchParams) const; inline bool Get_DistanceToNearestTrack(const DBCALShower* locBCALShower, double& locDistance) const; inline bool Get_DistanceToNearestTrack(const DBCALShower* locBCALShower, double& locDeltaPhi, double& locDeltaZ) const; @@ -194,6 +205,8 @@ class DDetectorMatches : public JObject inline size_t Get_NumTrackSCMatches(void) const; inline size_t Get_NumTrackDIRCMatches(void) const; inline size_t Get_NumTrackCTOFMatches(void) const; + inline size_t Get_NumTrackFMWPCMatches(void) const; + //SETTERS: inline void Add_Match(const DTrackingData* locTrack, const DBCALShower* locBCALShower, const shared_ptr& locShowerMatchParams); @@ -203,6 +216,9 @@ class DDetectorMatches : public JObject inline void Add_Match(const DTrackingData* locTrack, const DCTOFPoint* locCTOFPoint, const shared_ptr& locHitMatchParams); inline void Add_Match(const DTrackingData* locTrack, const DSCHit* locSCHit, const shared_ptr& locHitMatchParams); inline void Add_Match(const DTrackingData* locTrack, const shared_ptr& locDIRCMatchParams); + inline void Add_Match(const DTrackingData *locTrack, + const shared_ptr& locFMWPCMatchParams); + inline void Set_DistanceToNearestTrack(const DBCALShower* locBCALShower, double locDeltaPhi, double locDeltaZ); inline void Set_DistanceToNearestTrack(const DFCALShower* locFCALShower, double locDistanceToNearestTrack); inline void Set_FlightTimePCorrelation(const DTrackingData* locTrack, DetectorSystem_t locDetectorSystem, double locCorrelation); @@ -214,7 +230,8 @@ class DDetectorMatches : public JObject AddString(items, "#_Track_TOF_Matches", "%d", Get_NumTrackTOFMatches()); AddString(items, "#_Track_SC_Matches", "%d", Get_NumTrackSCMatches()); AddString(items, "#_Track_DIRC_Matches", "%d", Get_NumTrackDIRCMatches()); - AddString(items, "#_Track_CTOF_Matches", "%d", Get_NumTrackCTOFMatches()); + AddString(items, "#_Track_CTOF_Matches", "%d", Get_NumTrackCTOFMatches()); + AddString(items, "#_Track_FMWPC_Matches", "%d", Get_NumTrackFMWPCMatches()); } private: @@ -225,6 +242,7 @@ class DDetectorMatches : public JObject map > > dTrackFCALMatchParams; map > > dTrackTOFMatchParams; map > > dTrackCTOFMatchParams; + map > > dTrackFMWPCMatchParams; map > > dTrackSCMatchParams; map > dTrackDIRCMatchParams; @@ -285,6 +303,16 @@ inline bool DDetectorMatches::Get_CTOFMatchParams(const DTrackingData* locTrack, return true; } +inline bool DDetectorMatches::Get_FMWPCMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const +{ + locMatchParams.clear(); + auto locIterator = dTrackFMWPCMatchParams.find(locTrack); + if(locIterator == dTrackFMWPCMatchParams.end()) + return false; + locMatchParams = locIterator->second; + return true; +} + inline bool DDetectorMatches::Get_TOFMatchParams(const DTrackingData* locTrack, vector >& locMatchParams) const { locMatchParams.clear(); @@ -334,8 +362,6 @@ inline bool DDetectorMatches::Get_IsMatchedToHit(const DTrackingData* locTrack) return true; if(dTrackSCMatchParams.find(locTrack) != dTrackSCMatchParams.end()) return true; - if(dTrackCTOFMatchParams.find(locTrack) != dTrackCTOFMatchParams.end()) - return true; return false; } @@ -351,6 +377,8 @@ inline bool DDetectorMatches::Get_IsMatchedToDetector(const DTrackingData* locTr return (dTrackSCMatchParams.find(locTrack) != dTrackSCMatchParams.end()); else if(locDetectorSystem == SYS_CTOF) return (dTrackCTOFMatchParams.find(locTrack) != dTrackCTOFMatchParams.end()); + else if(locDetectorSystem == SYS_FMWPC) + return (dTrackFMWPCMatchParams.find(locTrack) != dTrackFMWPCMatchParams.end()); else return false; } @@ -385,8 +413,6 @@ inline bool DDetectorMatches::Get_TrackMatchParams(const DCTOFPoint* locCTOFPoin return true; } - - inline bool DDetectorMatches::Get_TrackMatchParams(const DTOFPoint* locTOFPoint, vector >& locMatchParams) const { locMatchParams.clear(); @@ -492,6 +518,15 @@ inline size_t DDetectorMatches::Get_NumTrackCTOFMatches(void) const return locNumTrackMatches; } +inline size_t DDetectorMatches::Get_NumTrackFMWPCMatches(void) const +{ + auto locIterator = dTrackFMWPCMatchParams.begin(); + unsigned int locNumTrackMatches = 0; + for(; locIterator != dTrackFMWPCMatchParams.end(); ++locIterator) + locNumTrackMatches += locIterator->second.size(); + return locNumTrackMatches; +} + inline size_t DDetectorMatches::Get_NumTrackSCMatches(void) const { auto locIterator = dTrackSCMatchParams.begin(); @@ -530,6 +565,10 @@ inline void DDetectorMatches::Add_Match(const DTrackingData* locTrack, const DCT dTrackCTOFMatchParams[locTrack].push_back(locHitMatchParams); dCTOFTrackMatchParams[locCTOFPoint].push_back(locHitMatchParams); } +inline void DDetectorMatches::Add_Match(const DTrackingData* locTrack, const shared_ptr& locHitMatchParams) +{ + dTrackFMWPCMatchParams[locTrack].push_back(locHitMatchParams); +} inline void DDetectorMatches::Add_Match(const DTrackingData* locTrack, const DSCHit* locSCHit, const shared_ptr& locHitMatchParams) { dTrackSCMatchParams[locTrack].push_back(locHitMatchParams); diff --git a/src/libraries/PID/DDetectorMatches_factory.cc b/src/libraries/PID/DDetectorMatches_factory.cc index d1f2462a0..50f9e778c 100644 --- a/src/libraries/PID/DDetectorMatches_factory.cc +++ b/src/libraries/PID/DDetectorMatches_factory.cc @@ -69,6 +69,9 @@ DDetectorMatches* DDetectorMatches_factory::Create_DDetectorMatches(jana::JEvent vector locCTOFPoints; locEventLoop->Get(locCTOFPoints); + vector locFMWPCClusters; + locEventLoop->Get(locFMWPCClusters); + DDetectorMatches* locDetectorMatches = new DDetectorMatches(); //Match tracks to showers/hits @@ -79,7 +82,10 @@ DDetectorMatches* DDetectorMatches_factory::Create_DDetectorMatches(jana::JEvent MatchToFCAL(locParticleID, locTrackTimeBasedVector[loc_i], locFCALShowers, locDetectorMatches); MatchToSC(locParticleID, locTrackTimeBasedVector[loc_i], locSCHits, locDetectorMatches); MatchToDIRC(locParticleID, locTrackTimeBasedVector[loc_i], locDIRCHits, locDetectorMatches, locDIRCBarHits); - MatchToCTOF(locParticleID, locTrackTimeBasedVector[loc_i], locCTOFPoints, locDetectorMatches); + if (locTrackTimeBasedVector[loc_i]->PID()<10){ // GEANT ids; ignore proton=14 and kaons=11+12 + MatchToCTOF(locParticleID, locTrackTimeBasedVector[loc_i], locCTOFPoints, locDetectorMatches); + MatchToFMWPC(locTrackTimeBasedVector[loc_i], locFMWPCClusters, locDetectorMatches); + } } //Find nearest tracks to showers @@ -141,6 +147,78 @@ void DDetectorMatches_factory::MatchToBCAL(const DParticleID* locParticleID, con } } +void DDetectorMatches_factory::MatchToFMWPC(const DTrackTimeBased* locTrackTimeBased, const vector& locFMWPCClusters, DDetectorMatches* locDetectorMatches) const{ + auto fmwpc_projections=locTrackTimeBased->extrapolations.at(SYS_FMWPC); + if (fmwpc_projections.size()==0) return; + + // Loop over projections filling in arrays of matching info + vectorlocLayers; + vectorlocNhits; + vectorlocDists; + vectorlocClosestWires; + const double FMWPC_WIRE_SPACING=1.016; //cm + bool got_match=false; + for( int layer=1; layer<=(int)fmwpc_projections.size(); layer++){ + auto proj = fmwpc_projections[layer-1]; + // x and y projections from track + double xpos=proj.position.x(); + double ypos=proj.position.y(); + + // Loop over DFMWPCClusters and find closest match to this projection + int min_dist = 1000000; + int wire_trk_proj=0; + const DFMWPCCluster* closest_fmwpc_cluster= nullptr; + for(auto fmwpccluster : locFMWPCClusters){ + if( fmwpccluster->layer != layer ) continue; + + // Convert into local coordinates so we can work with wire numbers + double s=fmwpccluster->orientation==DGeometry::kFMWPC_WIRE_ORIENTATION_VERTICAL ? xpos+fmwpccluster->xoffset : ypos+fmwpccluster->yoffset; + wire_trk_proj = round(71.5 + s/FMWPC_WIRE_SPACING) + 1; // 1-144 + + // If the projection is outside of the wire range then bail now + if( (wire_trk_proj<1) || (wire_trk_proj>144) ) continue; + + int dist=1000000; + if( wire_trk_proj >= fmwpccluster->first_wire ){ + dist = wire_trk_proj - fmwpccluster->last_wire; // distance beyond last wire (will be negative if inside cluster) + if( dist < 0 ) dist = 0; // force dist to 0 if inside cluster + }else{ + dist = fmwpccluster->first_wire - wire_trk_proj; // distance before first wire + } + + if( dist < min_dist ){ + min_dist = dist; + closest_fmwpc_cluster = fmwpccluster; + } + } + + // If a DFMWPCCluster was found, add the match info to the track + if( closest_fmwpc_cluster ){ + int closest_wire=wire_trk_proj; + if (wire_trk_proj < closest_fmwpc_cluster->first_wire ) { + closest_wire = closest_fmwpc_cluster->first_wire; + } + else if (wire_trk_proj > closest_fmwpc_cluster->last_wire){ + closest_wire = closest_fmwpc_cluster->last_wire; + } + locLayers.push_back(layer); + locNhits.push_back(closest_fmwpc_cluster->Nhits); + locDists.push_back(min_dist); + locClosestWires.push_back(closest_wire); + + got_match=true; + } + } + if (got_match){ + shared_ptrlocFMWPCMatchParams=std::make_shared(); + locFMWPCMatchParams->dLayers=locLayers; + locFMWPCMatchParams->dNhits=locNhits; + locFMWPCMatchParams->dDists=locDists; + locFMWPCMatchParams->dClosestWires=locClosestWires; + locDetectorMatches->Add_Match(locTrackTimeBased,locFMWPCMatchParams); + } +} + void DDetectorMatches_factory::MatchToCTOF(const DParticleID* locParticleID, const DTrackTimeBased* locTrackTimeBased, const vector& locCTOFPoints, DDetectorMatches* locDetectorMatches) const { vector extrapolations=locTrackTimeBased->extrapolations.at(SYS_CTOF); diff --git a/src/libraries/PID/DDetectorMatches_factory.h b/src/libraries/PID/DDetectorMatches_factory.h index da0fb1d5b..38324eac7 100644 --- a/src/libraries/PID/DDetectorMatches_factory.h +++ b/src/libraries/PID/DDetectorMatches_factory.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -48,6 +49,7 @@ class DDetectorMatches_factory : public jana::JFactory void MatchToSC(const DParticleID* locParticleID, const DTrackTimeBased* locTrackTimeBased, const vector& locSCHits, DDetectorMatches* locDetectorMatches) const; void MatchToDIRC(const DParticleID* locParticleID, const DTrackTimeBased* locTrackTimeBased, const vector& locDIRCHits, DDetectorMatches* locDetectorMatches, const vector& locDIRCBarHits) const; void MatchToCTOF(const DParticleID* locParticleID, const DTrackTimeBased* locTrackTimeBased, const vector& locCTOFPoints, DDetectorMatches* locDetectorMatches) const; + void MatchToFMWPC(const DTrackTimeBased* locTrackTimeBased, const vector& locFMWPCClusters, DDetectorMatches* locDetectorMatches) const; void MatchToFCAL(const DParticleID* locParticleID, const DTrackTimeBased *locTrackTimeBased, From fc133b1f5644b1c7e85e57674b40657b5de6738d Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 16 Jan 2024 10:48:45 -0500 Subject: [PATCH 02/13] add orientation and offset data to FMWPC clusters --- src/libraries/FMWPC/DFMWPCCluster.h | 5 ++++- src/libraries/FMWPC/DFMWPCCluster_factory.cc | 3 +++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/libraries/FMWPC/DFMWPCCluster.h b/src/libraries/FMWPC/DFMWPCCluster.h index 17ca14b99..9f59fe592 100644 --- a/src/libraries/FMWPC/DFMWPCCluster.h +++ b/src/libraries/FMWPC/DFMWPCCluster.h @@ -20,18 +20,21 @@ class DFMWPCCluster:public jana::JObject{ vector members; ///< DFMWPCHits that make up this cluster int layer; ///< #1-6 FMWPC layer + int orientation; ///< Vertical/Horizontal + float xoffset,yoffset; ///< x and y offsets of wires float q; ///< total charge in the cluster float u; ///< center of gravity of cluster in wire coordinates int first_wire; ///< first wire in cluster 1-144 int last_wire; ///< last wire in cluster 1-144 int Nhits; ///< Number of wire hits in cluster - DVector3 pos; ///< position (x,y,z) in global coodinates + DVector3 pos; ///< position (x,y,z) in global coordinates float t; // time in ns // This method is used primarily for pretty printing // the second argument to AddString is printf style format void toStrings(vector > &items)const{ AddString(items, "layer", "%d", layer); + AddString(items, "orientation", "%d", orientation); AddString(items, "q", "%10.2f", q); AddString(items, "u", "%3.4f", u); AddString(items, "first_wire", "%3d", first_wire); diff --git a/src/libraries/FMWPC/DFMWPCCluster_factory.cc b/src/libraries/FMWPC/DFMWPCCluster_factory.cc index bf57b2a10..c99e87a44 100644 --- a/src/libraries/FMWPC/DFMWPCCluster_factory.cc +++ b/src/libraries/FMWPC/DFMWPCCluster_factory.cc @@ -184,6 +184,9 @@ void DFMWPCCluster_factory::pique(vector& H) // global coordinate system // set to -777 for not measured coordinate auto orientation = fmwpc_wire_orientation[newCluster->layer-1]; + newCluster->orientation=orientation; + newCluster->xoffset=xvec[newCluster->layer-1]; + newCluster->yoffset=yvec[newCluster->layer-1]; double x = (orientation==DGeometry::kFMWPC_WIRE_ORIENTATION_VERTICAL ) ? xvec[newCluster->layer-1]+(newCluster->u-72.5)*FMWPC_WIRE_SPACING : -777 ; double y = (orientation==DGeometry::kFMWPC_WIRE_ORIENTATION_HORIZONTAL) ? yvec[newCluster->layer-1]+(newCluster->u-72.5)*FMWPC_WIRE_SPACING : -777 ; double z = zvec[newCluster->layer-1]; From 203419e555fa06b45a789c59ffd6d53a5e04fae7 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 16 Jan 2024 10:57:40 -0500 Subject: [PATCH 03/13] add fmwpc track match info to rest xml --- src/libraries/HDDM/DEventWriterREST.cc | 61 +++++++++++++++++--------- src/libraries/HDDM/rest.xml | 6 +++ 2 files changed, 46 insertions(+), 21 deletions(-) diff --git a/src/libraries/HDDM/DEventWriterREST.cc b/src/libraries/HDDM/DEventWriterREST.cc index efd0273fa..669abc1e5 100644 --- a/src/libraries/HDDM/DEventWriterREST.cc +++ b/src/libraries/HDDM/DEventWriterREST.cc @@ -606,27 +606,6 @@ bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutpu fcalList().setTflightvar(locFCALShowerMatchParamsVector[loc_k]->dFlightTimeVariance); } - vector> locCTOFHitMatchParamsVector; - locDetectorMatches[loc_i]->Get_CTOFMatchParams(tracks[loc_j], locCTOFHitMatchParamsVector); - for(size_t loc_k = 0; loc_k < locCTOFHitMatchParamsVector.size(); ++loc_k) - { - hddm_r::CtofMatchParamsList ctofList = matches().addCtofMatchParamses(1); - ctofList().setTrack(loc_j); - - size_t locCTOFindex = 0; - for(; locCTOFindex < ctofpoints.size(); ++locCTOFindex) - { - if(ctofpoints[locCTOFindex] == locCTOFHitMatchParamsVector[loc_k]->dCTOFPoint) - break; - } - ctofList().setHit(locCTOFindex); - ctofList().setDEdx(locCTOFHitMatchParamsVector[loc_k]->dEdx); - ctofList().setTflight(locCTOFHitMatchParamsVector[loc_k]->dFlightTime); - - ctofList().setDeltax(locCTOFHitMatchParamsVector[loc_k]->dDeltaXToHit); - ctofList().setDeltay(locCTOFHitMatchParamsVector[loc_k]->dDeltaYToHit); - } - vector> locFCALSingleHitMatchParamsVector; locDetectorMatches[loc_i]->Get_FCALSingleHitMatchParams(tracks[loc_j], locFCALSingleHitMatchParamsVector); for (size_t loc_k = 0; loc_k < locFCALSingleHitMatchParamsVector.size(); ++loc_k) @@ -772,6 +751,46 @@ bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutpu correlationList().setSystem(SYS_START); correlationList().setCorrelation(locFlightTimePCorrelation); } + + //--------- The following are for CPP -------------// + vector> locFMWPCMatchParamsVector; + locDetectorMatches[loc_i]->Get_FMWPCMatchParams(tracks[loc_j], locFMWPCMatchParamsVector); + for(size_t loc_k = 0; loc_k < locFMWPCMatchParamsVector.size(); ++loc_k){ + hddm_r::FmwpcMatchParamsList fmwpcList = matches().addFmwpcMatchParamses(1); + fmwpcList().setTrack(loc_j); + vectorlocLayers=locFMWPCMatchParamsVector[loc_k]->dLayers; + vectorlocNhits=locFMWPCMatchParamsVector[loc_k]->dNhits; + vectorlocDists=locFMWPCMatchParamsVector[loc_k]->dDists; + vectorlocClosestWires=locFMWPCMatchParamsVector[loc_k]->dClosestWires; + for (size_t loc_m=0;loc_m> locCTOFHitMatchParamsVector; + locDetectorMatches[loc_i]->Get_CTOFMatchParams(tracks[loc_j], locCTOFHitMatchParamsVector); + for(size_t loc_k = 0; loc_k < locCTOFHitMatchParamsVector.size(); ++loc_k) + { + hddm_r::CtofMatchParamsList ctofList = matches().addCtofMatchParamses(1); + ctofList().setTrack(loc_j); + + size_t locCTOFindex = 0; + for(; locCTOFindex < ctofpoints.size(); ++locCTOFindex) + { + if(ctofpoints[locCTOFindex] == locCTOFHitMatchParamsVector[loc_k]->dCTOFPoint) + break; + } + ctofList().setHit(locCTOFindex); + ctofList().setDEdx(locCTOFHitMatchParamsVector[loc_k]->dEdx); + ctofList().setTflight(locCTOFHitMatchParamsVector[loc_k]->dFlightTime); + + ctofList().setDeltax(locCTOFHitMatchParamsVector[loc_k]->dDeltaXToHit); + ctofList().setDeltay(locCTOFHitMatchParamsVector[loc_k]->dDeltaYToHit); + } } for(size_t loc_j = 0; loc_j < bcalshowers.size(); ++loc_j) diff --git a/src/libraries/HDDM/rest.xml b/src/libraries/HDDM/rest.xml index 506f4e8d2..b9a2e5e06 100644 --- a/src/libraries/HDDM/rest.xml +++ b/src/libraries/HDDM/rest.xml @@ -140,6 +140,12 @@ track="int" hit="int" dEdx="float" deltax="float" deltay="float" tflight="float" tunit="ns" lunit="cm" dEdx_unit="GeV/cm"/> + + + Date: Tue, 16 Jan 2024 11:35:08 -0500 Subject: [PATCH 04/13] retrieve FMWPC track matching data from REST file --- src/libraries/HDDM/DEventSourceREST.cc | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/libraries/HDDM/DEventSourceREST.cc b/src/libraries/HDDM/DEventSourceREST.cc index 7151c2ccb..b981bbcf7 100644 --- a/src/libraries/HDDM/DEventSourceREST.cc +++ b/src/libraries/HDDM/DEventSourceREST.cc @@ -1794,7 +1794,26 @@ jerror_t DEventSourceREST::Extract_DDetectorMatches(JEventLoop* locEventLoop, hd } } } - + + // Extract track matching data for FMPWCs + const hddm_r::FmwpcMatchParamsList &fmwpcList = iter->getFmwpcMatchParamses(); + hddm_r::FmwpcMatchParamsList::iterator fmwpcIter = fmwpcList.begin(); + for(; fmwpcIter != fmwpcList.end(); ++fmwpcIter) + { + size_t locTrackIndex = fmwpcIter->getTrack(); + const hddm_r::FmwpcDataList &fmwpcDataList = fmwpcIter->getFmwpcDatas(); + hddm_r::FmwpcDataList::iterator fmwpcDataIter = fmwpcDataList.begin(); + + auto locFMWPCMatchParams = std::make_shared(); + for(; fmwpcDataIter != fmwpcDataList.end(); ++fmwpcDataIter){ + locFMWPCMatchParams->dLayers.push_back(fmwpcDataIter->getLayer()); + locFMWPCMatchParams->dNhits.push_back(fmwpcDataIter->getNhits()); + locFMWPCMatchParams->dDists.push_back(fmwpcDataIter->getDist()); + locFMWPCMatchParams->dClosestWires.push_back(fmwpcDataIter->getClosestwire()); + } + locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], std::const_pointer_cast(locFMWPCMatchParams)); + } + // Extract track matching data for CTOF const hddm_r::CtofMatchParamsList &ctofList = iter->getCtofMatchParamses(); hddm_r::CtofMatchParamsList::iterator ctofIter = ctofList.begin(); From d8e2af758a8840b316bffa4d20d94ab9469d0813 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Tue, 16 Jan 2024 14:47:37 -0500 Subject: [PATCH 05/13] remove debug printout --- src/libraries/PID/DChargedTrackHypothesis_factory.cc | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/libraries/PID/DChargedTrackHypothesis_factory.cc b/src/libraries/PID/DChargedTrackHypothesis_factory.cc index f226ba2d2..de1f163fb 100644 --- a/src/libraries/PID/DChargedTrackHypothesis_factory.cc +++ b/src/libraries/PID/DChargedTrackHypothesis_factory.cc @@ -288,15 +288,6 @@ DChargedTrackHypothesis* DChargedTrackHypothesis_factory::Create_ChargedTrackHyp if (locDetectorMatches->Get_FMWPCMatchParams(locTrackTimeBased, locFMWPCMatchParamsVec)){ locChargedTrackHypothesis->Set_FMWPCMatchParams(locFMWPCMatchParamsVec[0]); - shared_ptr myshared=locChargedTrackHypothesis->Get_FMWPCMatchParams(); - if (myshared){ - for (int k=0;kdLayers.size();k++){ - cout << myshared->dLayers[k]<<": " <dNhits[k] - << " " << myshared->dDists[k] - <<" " << myshared->dClosestWires[k]< Date: Tue, 16 Jan 2024 16:40:06 -0500 Subject: [PATCH 06/13] add additional fcal energy info to match class, needed for CPP --- src/libraries/HDDM/DEventSourceREST.cc | 11 ++++++++++- src/libraries/HDDM/DEventWriterREST.cc | 5 +++++ src/libraries/HDDM/rest.xml | 5 ++++- 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/libraries/HDDM/DEventSourceREST.cc b/src/libraries/HDDM/DEventSourceREST.cc index b981bbcf7..9b44aa0f9 100644 --- a/src/libraries/HDDM/DEventSourceREST.cc +++ b/src/libraries/HDDM/DEventSourceREST.cc @@ -1699,7 +1699,16 @@ jerror_t DEventSourceREST::Extract_DDetectorMatches(JEventLoop* locEventLoop, hd locShowerMatchParams->dFlightTimeVariance = fcalIter->getTflightvar(); locShowerMatchParams->dPathLength = fcalIter->getPathlength(); locShowerMatchParams->dDOCAToShower = fcalIter->getDoca(); - + locShowerMatchParams->dEcenter=0.; + locShowerMatchParams->dE3x3=0.; + locShowerMatchParams->dE5x5=0.; + const hddm_r::FcalEnergyParamsList &fcalEnergyList = fcalIter->getFcalEnergyParamses(); + hddm_r::FcalEnergyParamsList::iterator fcalEnergyIter = fcalEnergyList.begin(); + for(; fcalEnergyIter != fcalEnergyList.end(); ++fcalEnergyIter){ + locShowerMatchParams->dEcenter=fcalEnergyIter->getEcenter(); + locShowerMatchParams->dE3x3=fcalEnergyIter->getE3x3(); + locShowerMatchParams->dE5x5=fcalEnergyIter->getE5x5(); + } locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], locFCALShowers[locShowerIndex], std::const_pointer_cast(locShowerMatchParams)); } diff --git a/src/libraries/HDDM/DEventWriterREST.cc b/src/libraries/HDDM/DEventWriterREST.cc index 669abc1e5..d28e2f135 100644 --- a/src/libraries/HDDM/DEventWriterREST.cc +++ b/src/libraries/HDDM/DEventWriterREST.cc @@ -604,6 +604,11 @@ bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutpu fcalList().setPathlength(locFCALShowerMatchParamsVector[loc_k]->dPathLength); fcalList().setTflight(locFCALShowerMatchParamsVector[loc_k]->dFlightTime); fcalList().setTflightvar(locFCALShowerMatchParamsVector[loc_k]->dFlightTimeVariance); + // Additional energy information + hddm_r::FcalEnergyParamsList fcalEnergyParamsList = fcalList().addFcalEnergyParamses(1); + fcalEnergyParamsList().setEcenter(locFCALShowerMatchParamsVector[loc_k]->dEcenter); + fcalEnergyParamsList().setE3x3(locFCALShowerMatchParamsVector[loc_k]->dE3x3); + fcalEnergyParamsList().setE5x5(locFCALShowerMatchParamsVector[loc_k]->dE5x5); } vector> locFCALSingleHitMatchParamsVector; diff --git a/src/libraries/HDDM/rest.xml b/src/libraries/HDDM/rest.xml index b9a2e5e06..ae8ae554b 100644 --- a/src/libraries/HDDM/rest.xml +++ b/src/libraries/HDDM/rest.xml @@ -131,7 +131,10 @@ + tunit="ns" lunit="cm"> + + Date: Wed, 17 Jan 2024 14:37:02 -0500 Subject: [PATCH 07/13] add additional FCAL data needed for CPP analysis to the FCAL/track match class. These additional data are not by default added to the rest files; this can be turned on by switching the ADD_FCAL_DATA_FOR_CPP flag to true --- src/libraries/HDDM/DEventWriterREST.cc | 10 +++++--- src/libraries/HDDM/DEventWriterREST.h | 1 + src/libraries/PID/DParticleID.cc | 34 ++++++++++++++++++++++++++ src/libraries/PID/DParticleID.h | 1 + 4 files changed, 42 insertions(+), 4 deletions(-) diff --git a/src/libraries/HDDM/DEventWriterREST.cc b/src/libraries/HDDM/DEventWriterREST.cc index d28e2f135..626ade589 100644 --- a/src/libraries/HDDM/DEventWriterREST.cc +++ b/src/libraries/HDDM/DEventWriterREST.cc @@ -605,10 +605,12 @@ bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutpu fcalList().setTflight(locFCALShowerMatchParamsVector[loc_k]->dFlightTime); fcalList().setTflightvar(locFCALShowerMatchParamsVector[loc_k]->dFlightTimeVariance); // Additional energy information - hddm_r::FcalEnergyParamsList fcalEnergyParamsList = fcalList().addFcalEnergyParamses(1); - fcalEnergyParamsList().setEcenter(locFCALShowerMatchParamsVector[loc_k]->dEcenter); - fcalEnergyParamsList().setE3x3(locFCALShowerMatchParamsVector[loc_k]->dE3x3); - fcalEnergyParamsList().setE5x5(locFCALShowerMatchParamsVector[loc_k]->dE5x5); + if (ADD_FCAL_DATA_FOR_CPP){ + hddm_r::FcalEnergyParamsList fcalEnergyParamsList = fcalList().addFcalEnergyParamses(1); + fcalEnergyParamsList().setEcenter(locFCALShowerMatchParamsVector[loc_k]->dEcenter); + fcalEnergyParamsList().setE3x3(locFCALShowerMatchParamsVector[loc_k]->dE3x3); + fcalEnergyParamsList().setE5x5(locFCALShowerMatchParamsVector[loc_k]->dE5x5); + } } vector> locFCALSingleHitMatchParamsVector; diff --git a/src/libraries/HDDM/DEventWriterREST.h b/src/libraries/HDDM/DEventWriterREST.h index 3af022d59..98b360165 100644 --- a/src/libraries/HDDM/DEventWriterREST.h +++ b/src/libraries/HDDM/DEventWriterREST.h @@ -59,6 +59,7 @@ class DEventWriterREST : public JObject bool REST_WRITE_DIRC_HITS; bool REST_WRITE_CCAL_SHOWERS; bool REST_WRITE_TRACK_EXIT_PARAMS; + bool ADD_FCAL_DATA_FOR_CPP; // metadata to save in the REST file // these should be consistent during program execution diff --git a/src/libraries/PID/DParticleID.cc b/src/libraries/PID/DParticleID.cc index 1e6781b32..3bec876a0 100644 --- a/src/libraries/PID/DParticleID.cc +++ b/src/libraries/PID/DParticleID.cc @@ -46,6 +46,8 @@ DParticleID::DParticleID(JEventLoop *loop) CDC_TRUNCATE_DEDX = true; gPARMS->SetDefaultParameter("PID:CDC_TRUNCATE_DEDX",CDC_TRUNCATE_DEDX); + ADD_FCAL_DATA_FOR_CPP=false; + gPARMS->SetDefaultParameter("PID:ADD_FCAL_DATA_FOR_CPP",ADD_FCAL_DATA_FOR_CPP); DApplication* dapp = dynamic_cast(loop->GetJApplication()); if(!dapp){ @@ -1302,6 +1304,38 @@ bool DParticleID::Distance_ToTrack(const vector & locShowerMatchParams->dPathLength = locPathLength; locShowerMatchParams->dDOCAToShower = d; + // The following additional information is needed for CPP ML mu/pi separation + if (ADD_FCAL_DATA_FOR_CPP){ + // Row/column corresponding to projected position + auto row = dFCALGeometry->row( (float)locProjPos.y() ); + auto col = dFCALGeometry->column( (float)locProjPos.x() ); + if (row>=0 && col>=0){ + locShowerMatchParams->dE5x5=0.; + locShowerMatchParams->dE3x3=0.; + locShowerMatchParams->dEcenter=0.; + // Get the list of hits from the cluster associated with the shower + const DFCALCluster*cluster=NULL; + locFCALShower->GetSingle(cluster); + if (cluster){ + vectorhits=cluster->GetHits(); + for( auto hit : hits ){ + auto delta_row = abs( dFCALGeometry->row(hit.ch) - row ); + auto delta_col = abs( dFCALGeometry->column(hit.ch) - col ); + cout << "Deltas " << delta_row <<" " << delta_col << endl; + if( (delta_row<=2) && (delta_col<=2) ){ + locShowerMatchParams->dE5x5 += hit.E; + if( (delta_row<=1) && (delta_col<=1) ){ + locShowerMatchParams->dE3x3 += hit.E; + if( (delta_row==0) && (delta_col==0) ){ + locShowerMatchParams->dEcenter = hit.E; + } + } + } + } + } + } + } // for CPP + return true; } bool DParticleID::Distance_ToTrack(const vector&extrapolations, const DTOFPoint* locTOFPoint, double locInputStartTime,shared_ptr& locTOFHitMatchParams, DVector3* locOutputProjPos, DVector3* locOutputProjMom) const diff --git a/src/libraries/PID/DParticleID.h b/src/libraries/PID/DParticleID.h index 56ca03eaa..5602c296b 100644 --- a/src/libraries/PID/DParticleID.h +++ b/src/libraries/PID/DParticleID.h @@ -333,6 +333,7 @@ class DParticleID:public jana::JObject double CDC_TIME_CUT_FOR_DEDX; bool CDC_TRUNCATE_DEDX; // dE/dx truncation: ignore hits with highest dE + bool ADD_FCAL_DATA_FOR_CPP; double dTargetZCenter; From 9db4bb05cf4f7bc8e0fca726fc16d23bea4677e5 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Wed, 17 Jan 2024 14:39:33 -0500 Subject: [PATCH 08/13] add flag to optionally add additional information to FCAL track matches for CPP --- src/libraries/HDDM/DEventWriterREST.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/libraries/HDDM/DEventWriterREST.cc b/src/libraries/HDDM/DEventWriterREST.cc index 626ade589..44f4352e4 100644 --- a/src/libraries/HDDM/DEventWriterREST.cc +++ b/src/libraries/HDDM/DEventWriterREST.cc @@ -51,6 +51,10 @@ DEventWriterREST::DEventWriterREST(JEventLoop* locEventLoop, string locOutputFil REST_WRITE_CCAL_SHOWERS = true; gPARMS->SetDefaultParameter("REST:WRITE_CCAL_SHOWERS", REST_WRITE_CCAL_SHOWERS); + ADD_FCAL_DATA_FOR_CPP=false; + gPARMS->SetDefaultParameter("PID:ADD_FCAL_DATA_FOR_CPP",ADD_FCAL_DATA_FOR_CPP); + + CCDB_CONTEXT_STRING = ""; // if we can get the calibration context from the DANA interface, then save this as well DApplication *dapp = dynamic_cast(locEventLoop->GetJApplication()); From 1e9ebde9675c6911cf77c7601a19f4bf35739b64 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Wed, 17 Jan 2024 14:40:45 -0500 Subject: [PATCH 09/13] add additional FCAL data to track match params for CPP --- src/libraries/PID/DDetectorMatches.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/libraries/PID/DDetectorMatches.h b/src/libraries/PID/DDetectorMatches.h index 63f38b666..d957b6dd6 100644 --- a/src/libraries/PID/DDetectorMatches.h +++ b/src/libraries/PID/DDetectorMatches.h @@ -65,7 +65,7 @@ class DFCALShowerMatchParams { public: DFCALShowerMatchParams(void) : dFCALShower(NULL), - dx(0.0), dFlightTime(0.0), dFlightTimeVariance(0.0), dPathLength(0.0), dDOCAToShower(0.0){} + dx(0.0), dFlightTime(0.0), dFlightTimeVariance(0.0), dPathLength(0.0), dDOCAToShower(0.0), dE5x5(0.0), dE3x3(0.0), dEcenter(0.0){} const DFCALShower* dFCALShower; @@ -74,6 +74,9 @@ class DFCALShowerMatchParams double dFlightTimeVariance; double dPathLength; //path length from DKinematicData::position() to the shower double dDOCAToShower; //DOCA of track to shower + double dE5x5; // energy sum over 5x5 array centered on track projection + double dE3x3; // energy sum over 3x3 array centered on track projection + double dEcenter; // energy in "center" block }; class DCTOFHitMatchParams From 7e9dc4ae8811b29ad98a226cc35b20f58eb84e19 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Thu, 18 Jan 2024 10:03:40 -0500 Subject: [PATCH 10/13] extract features for mu/pi ML from detector match data for pi+ and pi- tracks. --- src/libraries/FMWPC/DCPPEpEm_factory.cc | 120 +++++++++++++----------- src/libraries/FMWPC/DCPPEpEm_factory.h | 3 +- 2 files changed, 68 insertions(+), 55 deletions(-) diff --git a/src/libraries/FMWPC/DCPPEpEm_factory.cc b/src/libraries/FMWPC/DCPPEpEm_factory.cc index ee9bd1546..527a3283d 100644 --- a/src/libraries/FMWPC/DCPPEpEm_factory.cc +++ b/src/libraries/FMWPC/DCPPEpEm_factory.cc @@ -123,13 +123,13 @@ jerror_t DCPPEpEm_factory::evnt(JEventLoop *loop, uint64_t eventnumber) unsigned int ip=(q1>q2)?0:1; unsigned int in=(q1>q2)?1:0; - hyp=tracks[ip]->Get_Hypothesis(PiPlus); - if (hyp==NULL) return RESOURCE_UNAVAILABLE; - const DTrackTimeBased *piplus=hyp->Get_TrackTimeBased(); + const DChargedTrackHypothesis *PiPhyp=tracks[ip]->Get_Hypothesis(PiPlus); + if (PiPhyp==NULL) return RESOURCE_UNAVAILABLE; + const DTrackTimeBased *piplus=PiPhyp->Get_TrackTimeBased(); - hyp=tracks[in]->Get_Hypothesis(PiMinus); - if (hyp==NULL) return RESOURCE_UNAVAILABLE; - const DTrackTimeBased *piminus=hyp->Get_TrackTimeBased(); + const DChargedTrackHypothesis *PiMhyp=tracks[in]->Get_Hypothesis(PiMinus); + if (PiMhyp==NULL) return RESOURCE_UNAVAILABLE; + const DTrackTimeBased *piminus=PiMhyp->Get_TrackTimeBased(); hyp=tracks[ip]->Get_Hypothesis(KPlus); if (hyp==NULL) return RESOURCE_UNAVAILABLE; @@ -285,7 +285,7 @@ jerror_t DCPPEpEm_factory::evnt(JEventLoop *loop, uint64_t eventnumber) // n.b. if we need to use a mutex then we should pass a local // array for "input" and the lock the mutex just for the copy // to the tflite tensor. - if( PiMuFillFeatures(loop, piplus, piminus, pimu_input) ){ + if( PiMuFillFeatures(loop, tracks.size(), PiPhyp, PiMhyp, pimu_input) ){ // Run inference if( pimu_interpreter->Invoke() == kTfLiteOk){ @@ -418,26 +418,9 @@ bool DCPPEpEm_factory::VetoNeutrals(double t0_rf,const DVector3 &vect, // Return true if values are valid, false otherwise. // e.g. return false if there is not at least 1 pi+ // and 1 pi- candidate. -bool DCPPEpEm_factory::PiMuFillFeatures(jana::JEventLoop *loop, const DTrackTimeBased *piplus, const DTrackTimeBased *piminus, float *features){ - - vector matchedtracks; - loop->Get( matchedtracks ); - - // Find the DFMWPCMatchedTrack objects corresponding to - // the piplus, piminus tracks used in for the kinematic fit - // (i.e. the ones passed into this method as arguments) - const DFMWPCMatchedTrack *piplus_mt = nullptr; - const DFMWPCMatchedTrack *piminus_mt = nullptr; - for( auto mt : matchedtracks ){ - const DTrackTimeBased *trk; - mt->GetSingleT(trk); - if( trk == piplus ) piplus_mt = mt; - if( trk == piminus ) piminus_mt = mt; - } - - // Must have a DFMWPCMatchedTrack for both a pi+ and pi- - if( (piplus_mt==nullptr) || (piminus_mt==nullptr) ) return false; - +bool DCPPEpEm_factory::PiMuFillFeatures(jana::JEventLoop *loop, unsigned int nChargedTracks,const DChargedTrackHypothesis *PiPhyp, const DChargedTrackHypothesis *PiMhyp, float *features){ + memset(features,0,47*sizeof(float)); + // Features list is the following: // 0 nChargedTracks // 1 nFCALShowers @@ -487,53 +470,82 @@ bool DCPPEpEm_factory::PiMuFillFeatures(jana::JEventLoop *loop, const DTrackTime // 45 FMWPC_dist_closest_wire6_9 // 46 FMWPC_Nhits_cluster6_9 - vector chargedtracks; vector fcalshowers; vector fcalhits; vector fmwpchits; - loop->Get( chargedtracks ); + vectorstats; + loop->Get( stats ); loop->Get( fcalshowers ); loop->Get( fcalhits ); loop->Get( fmwpchits ); - features[ 0] = chargedtracks.size(); + features[ 0] = nChargedTracks; features[ 1] = fcalshowers.size(); - features[ 2] = fcalhits.size(); + features[ 2] = (stats.size()>0) ? stats[0]->fcal_blocks : fcalhits.size(); features[ 3] = fmwpchits.size(); - features[ 4] = matchedtracks.size(); - features[ 5] = piplus_mt->FCAL_E_center; - features[ 6] = piplus_mt->FCAL_E_3x3; - features[ 7] = piplus_mt->FCAL_E_5x5; - for(int ilayer=0; ilayer<6; ilayer++){ - features[ 8+3*ilayer] = piplus_mt->FMWPC_closest_wire[ilayer]; - features[ 9+3*ilayer] = piplus_mt->FMWPC_dist_closest_wire[ilayer]; - features[10+3*ilayer] = piplus_mt->FMWPC_Nhits_cluster[ilayer]; - } - features[26] = piminus_mt->FCAL_E_center; - features[27] = piminus_mt->FCAL_E_3x3; - features[28] = piminus_mt->FCAL_E_5x5; - for(int ilayer=0; ilayer<6; ilayer++){ - features[29+3*ilayer] = piminus_mt->FMWPC_closest_wire[ilayer]; - features[30+3*ilayer] = piminus_mt->FMWPC_dist_closest_wire[ilayer]; - features[31+3*ilayer] = piminus_mt->FMWPC_Nhits_cluster[ilayer]; + features[ 4] = 2; // piPlus + piMinus tracks found ?? + + // Match to FCAL for pi+ hypothesis + shared_ptrfcalparms=PiPhyp->Get_FCALShowerMatchParams(); + if (fcalparms!=nullptr){ + features[ 5] = fcalparms->dEcenter; + features[ 6] = fcalparms->dE3x3; + features[ 7] = fcalparms->dE5x5; } - + // Match to FMWPCs for pi+ hypothesis + shared_ptrfmwpcparms=PiPhyp->Get_FMWPCMatchParams(); // Before training the model, Nikhil's code replaced feature values // where the distance to the closest wire was >30 with values used // to indicate no wire hit. - for(int ilayer=0; ilayer<6; ilayer++){ - if( piminus_mt->FMWPC_dist_closest_wire[ilayer] >30.0 ){ - features[29+3*ilayer] = -1000.0; - features[30+3*ilayer] = 1000000; - features[31+3*ilayer] = 0; + for (int ilayer=0; ilayer<6; ilayer++){ + features[ 8+3*ilayer] = -1000.0; + features[ 9+3*ilayer] = 1000000; + features[10+3*ilayer] = 0; + } + if (fmwpcparms!=nullptr){ + for (unsigned int i=0;idLayers.size();i++){ + if (fmwpcparms->dDists[i]<30){ + int ilayer=fmwpcparms->dLayers[i]-1; + features[ 8+3*ilayer] = fmwpcparms->dClosestWires[i]; + features[ 9+3*ilayer] = fmwpcparms->dDists[i]; + features[10+3*ilayer] = fmwpcparms->dNhits[i]; + } } } + // Match to FCAL for pi- hypothesis + fcalparms=PiMhyp->Get_FCALShowerMatchParams(); + if (fcalparms!=nullptr){ + features[26] = fcalparms->dEcenter; + features[27] = fcalparms->dE3x3; + features[28] = fcalparms->dE5x5; + } + // Match to FMWPCs for pi- hypothesis + fmwpcparms=PiMhyp->Get_FMWPCMatchParams(); + // Before training the model, Nikhil's code replaced feature values + // where the distance to the closest wire was >30 with values used + // to indicate no wire hit. + for (int ilayer=0; ilayer<6; ilayer++){ + features[29+3*ilayer] = -1000.0; + features[30+3*ilayer] = 1000000; + features[31+3*ilayer] = 0; + } + if (fmwpcparms!=nullptr){ + for (unsigned int i=0;idLayers.size();i++){ + if (fmwpcparms->dDists[i]<30){ + int ilayer=fmwpcparms->dLayers[i]-1; + features[29+3*ilayer] = fmwpcparms->dClosestWires[i]; + features[30+3*ilayer] = fmwpcparms->dDists[i]; + features[31+3*ilayer] = fmwpcparms->dNhits[i]; + } + } + } + // These are values Nikhil sent that were used for normalizing the // features before training the model. static const float feature_min[] = {2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,0.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,-1000.0,0.0,0.0,0.0}; static const float feature_max[] = {6.0,10.0,20.0,94.0,8.0,3.924656391143799,5.177245497703552,5.349521217867732,144.0,1000000.0,39.0,144.0,1000000.0,17.0,144.0,1000000.0,12.0,144.0,1000000.0,11.0,144.0,1000000.0,8.0,144.0,1000000.0,7.0,4.154212951660156,5.578885164111853,5.9553504548966885,144.0,1000000.0,39.0,144.0,1000000.0,32.0,144.0,1000000.0,14.0,144.0,1000000.0,35.0,144.0,1000000.0,7.0,144.0,1000000.0,11.0,1.0}; - for(int i=0; i<48; i++){ + for(int i=0; i<47; i++){ features[i] = (features[i] - feature_min[i])/(feature_max[i]-feature_min[i]); } diff --git a/src/libraries/FMWPC/DCPPEpEm_factory.h b/src/libraries/FMWPC/DCPPEpEm_factory.h index 4712d6c54..d957f01c2 100644 --- a/src/libraries/FMWPC/DCPPEpEm_factory.h +++ b/src/libraries/FMWPC/DCPPEpEm_factory.h @@ -28,6 +28,7 @@ #include #include #include +#include class ReadMLPMinus; @@ -56,7 +57,7 @@ void DoKinematicFit(const DBeamPhoton *beamphoton, bool VetoNeutrals(double t0_rf,const DVector3 &vect, vector&neutrals) const; -bool PiMuFillFeatures(jana::JEventLoop *loop, const DTrackTimeBased *piplus, const DTrackTimeBased *piminus, float *features); + bool PiMuFillFeatures(jana::JEventLoop *loop, unsigned int nChargedTracks,const DChargedTrackHypothesis *piplus, const DChargedTrackHypothesis *piminus, float *features); double SPLIT_CUT,FCAL_THRESHOLD,BCAL_THRESHOLD,GAMMA_DT_CUT; string PIMU_MODEL_FILE; From b38561f1e353a2fd2084c734f7ac572f58555dbc Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Fri, 2 Feb 2024 13:11:42 -0500 Subject: [PATCH 11/13] remove debug printout --- src/libraries/PID/DParticleID.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/libraries/PID/DParticleID.cc b/src/libraries/PID/DParticleID.cc index 3bec876a0..017f172e4 100644 --- a/src/libraries/PID/DParticleID.cc +++ b/src/libraries/PID/DParticleID.cc @@ -1321,7 +1321,6 @@ bool DParticleID::Distance_ToTrack(const vector & for( auto hit : hits ){ auto delta_row = abs( dFCALGeometry->row(hit.ch) - row ); auto delta_col = abs( dFCALGeometry->column(hit.ch) - col ); - cout << "Deltas " << delta_row <<" " << delta_col << endl; if( (delta_row<=2) && (delta_col<=2) ){ locShowerMatchParams->dE5x5 += hit.E; if( (delta_row<=1) && (delta_col<=1) ){ From 5b6822c4940a3b01a427aae0e3a0f220e3d74013 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Fri, 2 Feb 2024 13:12:40 -0500 Subject: [PATCH 12/13] for CPP, check that the additional FCAL info added to the output makes sense --- src/libraries/HDDM/DEventWriterREST.cc | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/libraries/HDDM/DEventWriterREST.cc b/src/libraries/HDDM/DEventWriterREST.cc index 2d0d9c592..9ca81e1fa 100644 --- a/src/libraries/HDDM/DEventWriterREST.cc +++ b/src/libraries/HDDM/DEventWriterREST.cc @@ -633,10 +633,16 @@ bool DEventWriterREST::Write_RESTEvent(JEventLoop* locEventLoop, string locOutpu fcalList().setTflightvar(locFCALShowerMatchParamsVector[loc_k]->dFlightTimeVariance); // Additional energy information if (ADD_FCAL_DATA_FOR_CPP){ - hddm_r::FcalEnergyParamsList fcalEnergyParamsList = fcalList().addFcalEnergyParamses(1); - fcalEnergyParamsList().setEcenter(locFCALShowerMatchParamsVector[loc_k]->dEcenter); - fcalEnergyParamsList().setE3x3(locFCALShowerMatchParamsVector[loc_k]->dE3x3); - fcalEnergyParamsList().setE5x5(locFCALShowerMatchParamsVector[loc_k]->dE5x5); + // Sanity check for this additional info + double myE5x5=locFCALShowerMatchParamsVector[loc_k]->dE5x5; + double myE3x3=locFCALShowerMatchParamsVector[loc_k]->dE3x3; + double myEcenter=locFCALShowerMatchParamsVector[loc_k]->dEcenter; + if (myEcenter>0. || myE3x3>0. || myE5x5>0.){ + hddm_r::FcalEnergyParamsList fcalEnergyParamsList = fcalList().addFcalEnergyParamses(1); + fcalEnergyParamsList().setEcenter(myEcenter); + fcalEnergyParamsList().setE3x3(myE3x3); + fcalEnergyParamsList().setE5x5(myE5x5); + } } } From acc669c2641e9e7f7a5db2e28c36c1a0e7ec8726 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Wed, 14 Feb 2024 08:36:40 -0500 Subject: [PATCH 13/13] make features[4]=4 --- src/libraries/FMWPC/DCPPEpEm_factory.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libraries/FMWPC/DCPPEpEm_factory.cc b/src/libraries/FMWPC/DCPPEpEm_factory.cc index 527a3283d..f89f3f09c 100644 --- a/src/libraries/FMWPC/DCPPEpEm_factory.cc +++ b/src/libraries/FMWPC/DCPPEpEm_factory.cc @@ -483,7 +483,7 @@ bool DCPPEpEm_factory::PiMuFillFeatures(jana::JEventLoop *loop, unsigned int nCh features[ 1] = fcalshowers.size(); features[ 2] = (stats.size()>0) ? stats[0]->fcal_blocks : fcalhits.size(); features[ 3] = fmwpchits.size(); - features[ 4] = 2; // piPlus + piMinus tracks found ?? + features[ 4] = 4; // lepton and pion tracks found // Match to FCAL for pi+ hypothesis shared_ptrfcalparms=PiPhyp->Get_FCALShowerMatchParams();