@@ -52,7 +52,8 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo){
5252 // Pathlength (in AV) and start/end point
5353 pinfo.pathLength = PathLength ( part, pinfo.startPoint , pinfo.endPoint );
5454 // Central position of trajectory
55- pinfo.position = 0.5 *(pinfo.startPoint +pinfo.endPoint );
55+ pinfo.Position = geo::vect::middlePoint ({ pinfo.startPoint , pinfo.endPoint });
56+
5657 // Energy/charge deposited by this particle, found using SimEnergyDeposits
5758 pinfo.depEnergy = 0 ;
5859 pinfo.depElectrons = 0 ;
@@ -157,7 +158,7 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
157158
158159 // If this is a new blip, initialize
159160 if ( !tblip.G4ChargeMap .size () ) {
160- tblip.Position = pinfo.position ;
161+ tblip.Position = pinfo.Position ;
161162 tblip.Time = pinfo.time ;
162163
163164 // .. otherwise, check that the new particle
@@ -167,7 +168,9 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
167168 float totE = tblip.Energy + pinfo.depEnergy ;
168169 float w1 = tblip.Energy /totE;
169170 float w2 = pinfo.depEnergy /totE;
170- tblip.Position = w1*tblip.Position + w2*pinfo.position ;
171+ tblip.Position .SetXYZ ( w1*tblip.Position .X () + w2*pinfo.Position .X (),
172+ w1*tblip.Position .Y () + w2*pinfo.Position .Y (),
173+ w1*tblip.Position .Z () + w2*pinfo.Position .Z ());
171174 tblip.Time = w1*tblip.Time + w2*pinfo.time ;
172175 tblip.LeadCharge = pinfo.depElectrons ;
173176 // ... if the particle isn't a match, show's over
@@ -208,15 +211,17 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
208211 // check that the times are similar (we don't want to merge
209212 // together a blip that happened much later but in the same spot)
210213 if ( fabs (blip_i.Time - blip_j.Time ) > 5 ) continue ;
211- float d = (blip_i.Position -blip_j.Position ).Mag ();
214+ float d = (blip_i.Position -blip_j.Position ).R (); // Size of vector spanning two blips
212215 if ( d < dmin ) {
213216 isGrouped.at (j) = true ;
214217 // float totE = blip_i.Energy + blip_j.Energy;
215218 float totQ = blip_i.DepElectrons + blip_j.DepElectrons ;
216219 float w1 = blip_i.DepElectrons /totQ;
217220 float w2 = blip_j.DepElectrons /totQ;
218221 blip_i.Energy += blip_j.Energy ;
219- blip_i.Position = w1*blip_i.Position + w2*blip_j.Position ;
222+ blip_i.Position .SetXYZ ( w1*blip_i.Position .X () + w2*blip_j.Position .X (),
223+ w1*blip_i.Position .Y () + w2*blip_j.Position .Y (),
224+ w1*blip_i.Position .Z () + w2*blip_j.Position .Z ());
220225 blip_i.DriftTime = w1*blip_i.DriftTime + w2*blip_j.DriftTime ;
221226 blip_i.Time = w1*blip_i.Time + w2*blip_j.Time ;
222227 blip_i.DepElectrons += blip_j.DepElectrons ;
@@ -387,7 +392,7 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
387392 // ------------------------------------------------
388393 // / Look for valid wire intersections between
389394 // central-most hits in each cluster
390- std::vector<TVector3 > wirex;
395+ std::vector<geo::Point_t > wirex;
391396 for (size_t i=0 ; i<hcs.size (); i++) {
392397 int pli = hcs[i].Plane ;
393398 auto const & planegeo = wireReadoutGeom->Get ().Plane (geo::PlaneID{(unsigned int )hcs[i].Cryostat , (unsigned int )hcs[i].TPC , (unsigned int )hcs[i].Plane });
@@ -415,7 +420,7 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
415420 }
416421
417422 if ( match3d ) {
418- TVector3 a ( 0 ., intsec_p.Y (), intsec_p.Z ()) ;
423+ geo::Point_t a{ 0 ., intsec_p.Y (), intsec_p.Z ()} ;
419424 wirex.push_back (a);
420425 newblip.clusters [pli] = hcs[i];
421426 newblip.clusters [plj] = hcs[j];
@@ -430,12 +435,16 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
430435 // YZ-plane, as well as the mean difference between intersection points.
431436 newblip.Position .SetXYZ (0 ,0 ,0 );
432437 if ( wirex.size () == 1 ) {
433- newblip.Position = wirex[0 ];
438+ newblip.Position = wirex[0 ];
434439 } else {
435440 newblip.SigmaYZ = 0 ;
436441 double fact = 1 ./wirex.size ();
437- for (auto & v : wirex ) newblip.Position += v * fact;
438- for (auto & v : wirex ) newblip.SigmaYZ += (v-newblip.Position ).Mag () * fact;
442+ for (auto & v : wirex ) newblip.Position .SetXYZ ( newblip.Position .X () + v.X () * fact,
443+ newblip.Position .Y () + v.Y () * fact,
444+ newblip.Position .Z () + v.Z () * fact);
445+ for (auto & v : wirex ) newblip.SigmaYZ += TMath::Sqrt ( pow (v.X ()-newblip.Position .X (), 2 ) +
446+ pow (v.Y ()-newblip.Position .Y (), 2 ) +
447+ pow (v.Z ()-newblip.Position .Z (), 2 )) * fact;
439448 // Ensure that difference between intersection points is
440449 // consistent with the maximal wire extent
441450 if ( newblip.SigmaYZ > std::max (1 .,0.5 *newblip.dYZ ) ) return newblip;
@@ -625,7 +634,7 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
625634
626635 // =============================================================================
627636 // Length of particle trajectory
628- double PathLength (const simb::MCParticle& part, TVector3 & start, TVector3 & end)
637+ double PathLength (const simb::MCParticle& part, geo::Point_t & start, geo::Point_t & end)
629638 {
630639 int n = part.NumberTrajectoryPoints ();
631640 if ( n <= 1 ) return 0 .;
@@ -645,7 +654,7 @@ void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo,
645654 return L;
646655 }
647656 double PathLength (const simb::MCParticle& part){
648- TVector3 a,b;
657+ geo::Point_t a,b;
649658 return PathLength (part,a,b);
650659 }
651660
0 commit comments