diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index cb66fbd1..03c84e8d 100644 --- 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 index 95168ed9..0ef9d235 100644 --- 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 ) ; }