Fix phi flip and abs() bugs

iLCSoft · Aug 27, 2021 · 254db39 · 254db39
1 parent 0a9ff65
commit 254db39
Show file tree

Hide file tree

Showing 2 changed files with 63 additions and 53 deletions.
diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc
@@ -97,8 +97,8 @@ void TOFEstimators::init() {
         const auto& detector = dd4hep::Detector::getInstance();
         detector.field().magneticField({0., 0., 0.}, _bField);
     }
-    else { 
-        throw EVENT::Exception(std::string("Invalid ProcessorVersion parameter passed: \'") + _procVersion + std::string("\'\n Viable options are idr (default), dev")); 
+    else {
+        throw EVENT::Exception(std::string("Invalid ProcessorVersion parameter passed: \'") + _procVersion + std::string("\'\n Viable options are idr (default), dev"));
     }
 }
 
@@ -497,13 +497,18 @@ dd4hep::rec::Vector3D TOFEstimators::getMomAtCalo(const Track* track){
 
 double TOFEstimators::getFlightLength(const Track* track){
     const TrackState* ts = track->getTrackState(TrackState::AtIP);
-    double phiIP = ts->getPhi();
+    double phiIp = ts->getPhi();
     ts = track->getTrackState(TrackState::AtCalorimeter);
     double phiCalo = ts->getPhi();
     double omegaCalo = ts->getOmega();
     double tanLCalo = ts->getTanLambda();
 
-    return abs((phiIP - phiCalo)/omegaCalo)*sqrt(1. + tanLCalo*tanLCalo);
+    double dPhi = std::abs(phiIp - phiCalo);
+
+    // if dPhi > PI assume phi is flipped
+    if (dPhi > M_PI) dPhi = 2*M_PI - dPhi;
+
+    return dPhi/std::abs(omegaCalo)*std::sqrt(1. + tanLCalo*tanLCalo);
 }
 
 

diff --git a/TimeOfFlight/src/TOFUtils.cc b/TimeOfFlight/src/TOFUtils.cc
@@ -15,50 +15,55 @@
 
 namespace TOFUtils{
 
-  using namespace lcio ;
-
-
-  float computeFlightLength( EVENT::Track* trk){
-
-    const TrackState* tsIP = trk->getTrackState( TrackState::AtIP ) ; 	
-    const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; 	
-
-    float Phicalo = tscalo->getPhi() ;
-    float PhiIP = tsIP->getPhi() ;
-
-    float Omega = tsIP->getOmega()  ;
-    float tanL = tsIP->getTanLambda() ;
-
-    float length = (PhiIP-Phicalo)*(1/Omega) * sqrt( 1 + tanL * tanL ) ;
-
-    return length ;
-  }
+using namespace lcio ;
+
+
+float computeFlightLength( EVENT::Track* trk){
+
+    const TrackState* tsIp = trk->getTrackState( TrackState::AtIP ) ;
+    const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ;
+
+    float phiIp = tsIp->getPhi() ;
+    float phiCalo = tscalo->getPhi() ;
+
+    float omega = tsIp->getOmega()  ;
+    float tanL = tsIp->getTanLambda() ;
+
+    float dPhi = std::abs(phiIp - phiCalo);
 
-  float computeFlightLength(const TrackState* ts0,  const TrackState* ts1 ){
-
-    float Phicalo = ts1->getPhi() ;
-    float PhiIP = ts0->getPhi() ;
-
-    float Omega = ts0->getOmega()  ;
+    // if dPhi > PI assume phi is flipped
+    if (dPhi > M_PI) dPhi = 2*M_PI - dPhi;
+
+    return dPhi/std::abs(omega) * std::sqrt( 1. + tanL * tanL ) ;
+}
+
+float computeFlightLength(const TrackState* ts0,  const TrackState* ts1 ){
+    float phi1 = ts1->getPhi() ;
+    float phi0 = ts0->getPhi() ;
+
+    float omega = ts0->getOmega()  ;
     float tanL = ts0->getTanLambda() ;
-
-    float length = (PhiIP-Phicalo)*(1/Omega) * sqrt( 1 + tanL * tanL ) ;
 
-    return length ;
-  }
+    float dPhi = std::abs(phi0 - phi1);
+
+    // if dPhi > PI assume phi is flipped
+    if (dPhi > M_PI) dPhi = 2*M_PI - dPhi;
+
+    return dPhi/std::abs(omega) * std::sqrt( 1. + tanL * tanL ) ;
+}
 
 
   int layer( EVENT::CalorimeterHit* h ) {
-    
+
     return CHT( h->getType() ).layer() ;
   }
-  
-  
-  /// helper function to check if this is an Ecal hit 
+
+
+  /// helper function to check if this is an Ecal hit
   bool isEcal( EVENT::CalorimeterHit* h ) {
-    
-    return CHT( h->getType() ).caloID() == CHT::ecal ; 
-  } 
+
+    return CHT( h->getType() ).caloID() == CHT::ecal ;
+  }
 
 
   std::string caloTypeStr(  EVENT::CalorimeterHit* h ) {
@@ -72,7 +77,7 @@ namespace TOFUtils{
   std::string CaloHitData::toString(){
 
     std::stringstream s ;
-    s << "  l= " << layer ;  
+    s << "  l= " << layer ;
     s << ", st= " << smearedTime ;
     s << ", tr= " << timeResolution  ;
     s << ", dIP=" << distanceFromIP  ;
@@ -87,30 +92,30 @@ namespace TOFUtils{
   float computeDistanceFromLine( EVENT::CalorimeterHit* h, const dd4hep::rec::Vector3D& point,
 				 const dd4hep::rec::Vector3D& unitDir) {
 
-    dd4hep::rec::Vector3D pos( h->getPosition()[0], 
-			       h->getPosition()[1], 
+    dd4hep::rec::Vector3D pos( h->getPosition()[0],
+			       h->getPosition()[1],
 			       h->getPosition()[2] ) ;
-    
+
     dd4hep::rec::Vector3D diff = pos - point ;
 
-    
+
     return diff.cross( unitDir ).r() ;
-    
+
   }
- 
+
 
   CaloHitDataVec findHitsClosestToLine( const CaloHitLayerMap& layerMap ){
-    
+
     CaloHitDataVec hitVec ;
-    
+
     for( auto m : layerMap ){
-      
+
       const CaloHitDataVec& chv = m.second ;
-      
+
       CaloHitData* closestHit =
 	*std::min_element( chv.begin() , chv.end () ,
 			   [](CaloHitData* c0, CaloHitData* c1 ){ return c0->distancefromStraightline < c1->distancefromStraightline  ; }
-	  ) ; 
+	  ) ;
 
       hitVec.push_back( closestHit ) ;
     }
@@ -128,14 +133,14 @@ namespace TOFUtils{
     float meansq = 0. ;
 
     for( auto ch : chv ){
-      float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; 
+      float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ;
       mean   += t   ;
       ++nHit ;
     }
     mean /= nHit ;
 
     for( auto ch : chv ){
-      float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; 
+      float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ;
       meansq += ( t - mean ) * ( t - mean )   ;
     }