Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:36 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackAssociator
0004 // Class:      CachedTrajectory
0005 //
0006 //
0007 //
0008 
0009 #include "TrackingTools/TrackAssociator/interface/CachedTrajectory.h"
0010 // #include "Utilities/Timing/interface/TimerStack.h"
0011 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0012 #include "DataFormats/GeometrySurface/interface/Plane.h"
0013 #include "Geometry/DTGeometry/interface/DTChamber.h"
0014 #include "Geometry/CSCGeometry/interface/CSCChamber.h"
0015 #include <deque>
0016 #include <algorithm>
0017 
0018 std::vector<SteppingHelixStateInfo> propagateThoughFromIP(const SteppingHelixStateInfo& state,
0019                                                           const Propagator* prop,
0020                                                           const FiducialVolume& volume,
0021                                                           int nsteps,
0022                                                           float step,
0023                                                           float minR,
0024                                                           float minZ,
0025                                                           float maxR,
0026                                                           float maxZ) {
0027   CachedTrajectory neckLace;
0028   neckLace.setStateAtIP(state);
0029   neckLace.reset_trajectory();
0030   neckLace.setPropagator(prop);
0031   neckLace.setPropagationStep(0.1);
0032   neckLace.setMinDetectorRadius(minR);
0033   neckLace.setMinDetectorLength(minZ * 2.);
0034   neckLace.setMaxDetectorRadius(maxR);
0035   neckLace.setMaxDetectorLength(maxZ * 2.);
0036 
0037   // Propagate track
0038   bool isPropagationSuccessful = neckLace.propagateAll(state);
0039 
0040   if (!isPropagationSuccessful)
0041     return std::vector<SteppingHelixStateInfo>();
0042 
0043   std::vector<SteppingHelixStateInfo> complicatePoints;
0044   neckLace.getTrajectory(complicatePoints, volume, nsteps);
0045 
0046   return complicatePoints;
0047 }
0048 
0049 CachedTrajectory::CachedTrajectory() : propagator_(nullptr) {
0050   reset_trajectory();
0051   setMaxDetectorRadius();
0052   setMaxDetectorLength();
0053   setMinDetectorRadius();
0054   setMinDetectorLength();
0055   setPropagationStep();
0056 }
0057 
0058 void CachedTrajectory::propagateForward(SteppingHelixStateInfo& state, float distance) {
0059   // defined a normal plane wrt the particle trajectory direction
0060   // let's hope that I computed the rotation matrix correctly.
0061   GlobalVector vector(state.momentum().unit());
0062   float r21 = 0;
0063   float r22 = vector.z() / sqrt(1 - pow(vector.x(), 2));
0064   float r23 = -vector.y() / sqrt(1 - pow(vector.x(), 2));
0065   float r31 = vector.x();
0066   float r32 = vector.y();
0067   float r33 = vector.z();
0068   float r11 = r22 * r33 - r23 * r32;
0069   float r12 = r23 * r31;
0070   float r13 = -r22 * r31;
0071 
0072   Surface::RotationType rotation(r11, r12, r13, r21, r22, r23, r31, r32, r33);
0073   Plane::PlanePointer target = Plane::build(state.position() + vector * distance, rotation);
0074   propagate(state, *target);
0075 }
0076 
0077 void CachedTrajectory::propagate(SteppingHelixStateInfo& state, const Plane& plane) {
0078   if (const SteppingHelixPropagator* shp = dynamic_cast<const SteppingHelixPropagator*>(propagator_)) {
0079     try {
0080       shp->propagate(state, plane, state);
0081     } catch (cms::Exception& ex) {
0082       edm::LogWarning("TrackAssociator") << "Caught exception " << ex.category() << ": " << ex.explainSelf();
0083       edm::LogWarning("TrackAssociator") << "An exception is caught during the track propagation\n"
0084                                          << state.momentum().x() << ", " << state.momentum().y() << ", "
0085                                          << state.momentum().z();
0086       state = SteppingHelixStateInfo();
0087     }
0088   } else {
0089     FreeTrajectoryState fts;
0090     state.getFreeState(fts);
0091     TrajectoryStateOnSurface stateOnSurface = propagator_->propagate(fts, plane);
0092     state = SteppingHelixStateInfo(*(stateOnSurface.freeState()));
0093   }
0094 }
0095 
0096 void CachedTrajectory::propagate(SteppingHelixStateInfo& state, const Cylinder& cylinder) {
0097   if (const SteppingHelixPropagator* shp = dynamic_cast<const SteppingHelixPropagator*>(propagator_)) {
0098     try {
0099       shp->propagate(state, cylinder, state);
0100     } catch (cms::Exception& ex) {
0101       edm::LogWarning("TrackAssociator") << "Caught exception " << ex.category() << ": " << ex.explainSelf();
0102       edm::LogWarning("TrackAssociator") << "An exception is caught during the track propagation\n"
0103                                          << state.momentum().x() << ", " << state.momentum().y() << ", "
0104                                          << state.momentum().z();
0105       state = SteppingHelixStateInfo();
0106     }
0107   } else {
0108     FreeTrajectoryState fts;
0109     state.getFreeState(fts);
0110     TrajectoryStateOnSurface stateOnSurface = propagator_->propagate(fts, cylinder);
0111     state = SteppingHelixStateInfo(*(stateOnSurface.freeState()));
0112   }
0113 }
0114 
0115 bool CachedTrajectory::propagateAll(const SteppingHelixStateInfo& initialState) {
0116   if (fullTrajectoryFilled_) {
0117     edm::LogWarning("TrackAssociator")
0118         << "Reseting all trajectories. Please call reset_trajectory() explicitely to avoid this message";
0119     reset_trajectory();
0120   }
0121 
0122   //   TimerStack timers(TimerStack::Disableable);
0123 
0124   reset_trajectory();
0125   if (propagator_ == nullptr)
0126     throw cms::Exception("FatalError") << "Track propagator is not defined\n";
0127   SteppingHelixStateInfo currentState(initialState);
0128   fullTrajectory_.push_back(currentState);
0129 
0130   while (currentState.position().perp() < maxRho_ && fabs(currentState.position().z()) < maxZ_) {
0131     LogTrace("TrackAssociator") << "[propagateAll] Propagate outward from (rho, r, z, phi) ("
0132                                 << currentState.position().perp() << ", " << currentState.position().mag() << ", "
0133                                 << currentState.position().z() << ", " << currentState.position().phi() << ")";
0134     propagateForward(currentState, step_);
0135     if (!currentState.isValid()) {
0136       LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n";
0137       break;
0138     }
0139     LogTrace("TrackAssociator") << "\treached (rho, r, z, phi) (" << currentState.position().perp() << ", "
0140                                 << currentState.position().mag() << ", " << currentState.position().z() << ", "
0141                                 << currentState.position().phi() << ")";
0142     fullTrajectory_.push_back(currentState);
0143   }
0144 
0145   SteppingHelixStateInfo currentState2(initialState);
0146   SteppingHelixStateInfo previousState;
0147   while (currentState2.position().perp() > minRho_ || fabs(currentState2.position().z()) > minZ_) {
0148     previousState = currentState2;
0149     propagateForward(currentState2, -step_);
0150     if (!currentState2.isValid()) {
0151       LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n";
0152       break;
0153     }
0154     if (previousState.position().perp() - currentState2.position().perp() < 0) {
0155       LogTrace("TrackAssociator")
0156           << "Error: TrackAssociator has propogated the particle past the point of closest approach to IP" << std::endl;
0157       break;
0158     }
0159     LogTrace("TrackAssociator") << "[propagateAll] Propagated inward from (rho, r, z, phi) ("
0160                                 << previousState.position().perp() << ", " << previousState.position().mag() << ", "
0161                                 << previousState.position().z() << "," << previousState.position().phi() << ") to ("
0162                                 << currentState2.position().perp() << ", " << currentState2.position().mag() << ", "
0163                                 << currentState2.position().z() << ", " << currentState2.position().phi() << ")";
0164     fullTrajectory_.push_front(currentState2);
0165   }
0166 
0167   // LogTrace("TrackAssociator") << "fullTrajectory_ has " << fullTrajectory_.size() << " states with (R, z):\n";
0168   // for(unsigned int i=0; i<fullTrajectory_.size(); i++) {
0169   //  LogTrace("TrackAssociator") << "state " << i << ": (" << fullTrajectory_[i].position().perp() << ", "
0170   //    << fullTrajectory_[i].position().z() << ")\n";
0171   // }
0172 
0173   LogTrace("TrackAssociator") << "Done with the track propagation in the detector. Number of steps: "
0174                               << fullTrajectory_.size();
0175   fullTrajectoryFilled_ = true;
0176   return !fullTrajectory_.empty();
0177 }
0178 
0179 TrajectoryStateOnSurface CachedTrajectory::propagate(const Plane* plane) {
0180   // TimerStack timers(TimerStack::Disableable);
0181   // timers.benchmark("CachedTrajectory::propagate::benchmark");
0182   // timers.push("CachedTrajectory::propagate",TimerStack::FastMonitoring);
0183   // timers.push("CachedTrajectory::propagate::findClosestPoint",TimerStack::FastMonitoring);
0184 
0185   // Assume that all points along the trajectory are equally spread out.
0186   // For simplication assume that the trajectory is just a straight
0187   // line and find a point closest to the target plane. Propagate to
0188   // the plane from the point.
0189 
0190   const float matchingDistance = 1;
0191   // find the closest point to the plane
0192   int leftIndex = 0;
0193   int rightIndex = fullTrajectory_.size() - 1;
0194   int closestPointOnLeft = 0;
0195 
0196   // check whether the trajectory crossed the plane (signs should be different)
0197   if (sign(distance(plane, leftIndex)) * sign(distance(plane, rightIndex)) != -1) {
0198     LogTrace("TrackAssociator") << "Track didn't cross the plane:\n\tleft distance: " << distance(plane, leftIndex)
0199                                 << "\n\tright distance: " << distance(plane, rightIndex);
0200     return TrajectoryStateOnSurface();
0201   }
0202 
0203   while (leftIndex + 1 < rightIndex) {
0204     closestPointOnLeft = int((leftIndex + rightIndex) / 2);
0205     float dist = distance(plane, closestPointOnLeft);
0206     // LogTrace("TrackAssociator") << "Closest point on left: " << closestPointOnLeft << "\n"
0207     //    << "Distance to the plane: " << dist;
0208     if (fabs(dist) < matchingDistance) {
0209       // found close match, verify that we are on the left side
0210       if (closestPointOnLeft > 0 && sign(distance(plane, closestPointOnLeft - 1)) * dist == -1)
0211         closestPointOnLeft--;
0212       break;
0213     }
0214 
0215     // check where is the plane
0216     if (sign(distance(plane, leftIndex) * dist) == -1)
0217       rightIndex = closestPointOnLeft;
0218     else
0219       leftIndex = closestPointOnLeft;
0220 
0221     // LogTrace("TrackAssociator") << "Distance on left: " << distance(plane, leftIndex) << "\n"
0222     //  << "Distance to closest point: " <<  distance(plane, closestPointOnLeft) << "\n"
0223     //  << "Left index: " << leftIndex << "\n"
0224     //  << "Right index: " << rightIndex;
0225   }
0226   LogTrace("TrackAssociator") << "closestPointOnLeft: " << closestPointOnLeft << "\n\ttrajectory point (z,R,eta,phi): "
0227                               << fullTrajectory_[closestPointOnLeft].position().z() << ", "
0228                               << fullTrajectory_[closestPointOnLeft].position().perp() << " , "
0229                               << fullTrajectory_[closestPointOnLeft].position().eta() << " , "
0230                               << fullTrajectory_[closestPointOnLeft].position().phi()
0231                               << "\n\tplane center (z,R,eta,phi): " << plane->position().z() << ", "
0232                               << plane->position().perp() << " , " << plane->position().eta() << " , "
0233                               << plane->position().phi();
0234 
0235   // propagate to the plane
0236   // timers.pop_and_push("CachedTrajectory::propagate::localPropagation",TimerStack::FastMonitoring);
0237   if (const SteppingHelixPropagator* shp = dynamic_cast<const SteppingHelixPropagator*>(propagator_)) {
0238     SteppingHelixStateInfo state;
0239     try {
0240       shp->propagate(fullTrajectory_[closestPointOnLeft], *plane, state);
0241     } catch (cms::Exception& ex) {
0242       edm::LogWarning("TrackAssociator") << "Caught exception " << ex.category() << ": " << ex.explainSelf();
0243       edm::LogWarning("TrackAssociator") << "An exception is caught during the track propagation\n"
0244                                          << state.momentum().x() << ", " << state.momentum().y() << ", "
0245                                          << state.momentum().z();
0246       return TrajectoryStateOnSurface();
0247     }
0248     return state.getStateOnSurface(*plane);
0249   } else {
0250     FreeTrajectoryState fts;
0251     fullTrajectory_[closestPointOnLeft].getFreeState(fts);
0252     return propagator_->propagate(fts, *plane);
0253   }
0254 }
0255 
0256 std::pair<float, float> CachedTrajectory::trajectoryDelta(TrajectorType trajectoryType) {
0257   // MEaning of trajectory change depends on its usage. In most cases we measure
0258   // change in a trajectory as difference between final track position and initial
0259   // direction. In some cases such as change of trajectory in the muon detector we
0260   // might want to compare theta-phi of two points or even find local maximum and
0261   // mimimum. In general it's not essential what defenition of the trajectory change
0262   // is used since we use these numbers only as a rough estimate on how much wider
0263   // we should make the preselection region.
0264   std::pair<float, float> result(0, 0);
0265   if (!stateAtIP_.isValid()) {
0266     edm::LogWarning("TrackAssociator") << "State at IP is not known set. Cannot estimate trajectory change. "
0267                                        << "Trajectory change is not taken into account in matching";
0268     return result;
0269   }
0270   switch (trajectoryType) {
0271     case IpToEcal:
0272       if (ecalTrajectory_.empty())
0273         edm::LogWarning("TrackAssociator") << "ECAL trajector is empty. Cannot estimate trajectory change. "
0274                                            << "Trajectory change is not taken into account in matching";
0275       else
0276         return delta(stateAtIP_.momentum().theta(),
0277                      ecalTrajectory_.front().position().theta(),
0278                      stateAtIP_.momentum().phi(),
0279                      ecalTrajectory_.front().position().phi());
0280       break;
0281     case IpToHcal:
0282       if (hcalTrajectory_.empty())
0283         edm::LogWarning("TrackAssociator") << "HCAL trajector is empty. Cannot estimate trajectory change. "
0284                                            << "Trajectory change is not taken into account in matching";
0285       else
0286         return delta(stateAtIP_.momentum().theta(),
0287                      hcalTrajectory_.front().position().theta(),
0288                      stateAtIP_.momentum().phi(),
0289                      hcalTrajectory_.front().position().phi());
0290       break;
0291     case IpToHO:
0292       if (hoTrajectory_.empty())
0293         edm::LogWarning("TrackAssociator") << "HO trajector is empty. Cannot estimate trajectory change. "
0294                                            << "Trajectory change is not taken into account in matching";
0295       else
0296         return delta(stateAtIP_.momentum().theta(),
0297                      hoTrajectory_.front().position().theta(),
0298                      stateAtIP_.momentum().phi(),
0299                      hoTrajectory_.front().position().phi());
0300       break;
0301     case FullTrajectory:
0302       if (fullTrajectory_.empty())
0303         edm::LogWarning("TrackAssociator") << "Full trajector is empty. Cannot estimate trajectory change. "
0304                                            << "Trajectory change is not taken into account in matching";
0305       else
0306         return delta(stateAtIP_.momentum().theta(),
0307                      fullTrajectory_.back().position().theta(),
0308                      stateAtIP_.momentum().phi(),
0309                      fullTrajectory_.back().position().phi());
0310       break;
0311     default:
0312       edm::LogWarning("TrackAssociator")
0313           << "Unkown or not supported trajector type. Cannot estimate trajectory change. "
0314           << "Trajectory change is not taken into account in matching";
0315   }
0316   return result;
0317 }
0318 
0319 std::pair<float, float> CachedTrajectory::delta(const double& theta1,
0320                                                 const double& theta2,
0321                                                 const double& phi1,
0322                                                 const double& phi2) {
0323   std::pair<float, float> result(theta2 - theta1, phi2 - phi1);
0324   // this won't work for loopers, since deltaPhi cannot be larger than Pi.
0325   if (fabs(result.second) > 2 * M_PI - fabs(result.second)) {
0326     if (result.second > 0)
0327       result.second -= 2 * M_PI;
0328     else
0329       result.second += 2 * M_PI;
0330   }
0331   return result;
0332 }
0333 
0334 void CachedTrajectory::getTrajectory(std::vector<SteppingHelixStateInfo>& trajectory,
0335                                      const FiducialVolume& volume,
0336                                      int steps) {
0337   if (!fullTrajectoryFilled_)
0338     throw cms::Exception("FatalError") << "trajectory is not defined yet. Please use propagateAll first.";
0339   if (fullTrajectory_.empty()) {
0340     LogTrace("TrackAssociator") << "Trajectory is empty. Move on";
0341     return;
0342   }
0343   if (!volume.isValid()) {
0344     LogTrace("TrackAssociator") << "no trajectory is expected to be found since the fiducial volume is not valid";
0345     return;
0346   }
0347   double step = std::max(volume.maxR() - volume.minR(), volume.maxZ() - volume.minZ()) / steps;
0348 
0349   int closestPointOnLeft = -1;
0350 
0351   // check whether the trajectory crossed the region
0352   if (!((fullTrajectory_.front().position().perp() < volume.maxR() &&
0353          fabs(fullTrajectory_.front().position().z()) < volume.maxZ()) &&
0354         (fullTrajectory_.back().position().perp() > volume.minR() ||
0355          fabs(fullTrajectory_.back().position().z()) > volume.minZ()))) {
0356     LogTrace("TrackAssociator") << "Track didn't cross the region (R1,R2,L1,L2): " << volume.minR() << ", "
0357                                 << volume.maxR() << ", " << volume.minZ() << ", " << volume.maxZ();
0358     return;
0359   }
0360 
0361   // get distance along momentum to the surface.
0362 
0363   // the following code can be made faster, but it'll hardly be a significant improvement
0364   // simplifications:
0365   //   1) direct loop over stored trajectory points instead of some sort
0366   //      of fast root search (Newton method)
0367   //   2) propagate from the closest point outside the region with the
0368   //      requested step ignoring stored trajectory points.
0369   double dZ(-1.);
0370   double dR(-1.);
0371   int firstPointInside(-1);
0372   for (unsigned int i = 0; i < fullTrajectory_.size(); i++) {
0373     // LogTrace("TrackAssociator") << "Trajectory info (i,perp,r1,r2,z,z1,z2): " << i << ", " << fullTrajectory_[i].position().perp() <<
0374     //  ", " << volume.minR() << ", " << volume.maxR() << ", " << fullTrajectory_[i].position().z() << ", " << volume.minZ() << ", " <<
0375     //  volume.maxZ() << ", " << closestPointOnLeft;
0376     dR = fullTrajectory_[i].position().perp() - volume.minR();
0377     dZ = fabs(fullTrajectory_[i].position().z()) - volume.minZ();
0378     if (dR > 0 || dZ > 0) {
0379       if (i > 0) {
0380         firstPointInside = i;
0381         closestPointOnLeft = i - 1;
0382       } else {
0383         firstPointInside = 0;
0384         closestPointOnLeft = 0;
0385       }
0386       break;
0387     }
0388   }
0389   if (closestPointOnLeft == -1)
0390     throw cms::Exception("FatalError") << "This shouls never happen - internal logic error";
0391 
0392   SteppingHelixStateInfo currentState(fullTrajectory_[closestPointOnLeft]);
0393   if (currentState.position().x() * currentState.momentum().x() +
0394           currentState.position().y() * currentState.momentum().y() +
0395           currentState.position().z() * currentState.momentum().z() <
0396       0)
0397     step = -step;
0398 
0399   // propagate to the inner surface of the active volume
0400 
0401   if (firstPointInside != closestPointOnLeft) {
0402     if (dR > 0) {
0403       Cylinder::CylinderPointer barrel =
0404           Cylinder::build(volume.minR(), Cylinder::PositionType(0, 0, 0), Cylinder::RotationType());
0405       propagate(currentState, *barrel);
0406     } else {
0407       Plane::PlanePointer endcap =
0408           Plane::build(Plane::PositionType(0, 0, currentState.position().z() > 0 ? volume.minZ() : -volume.minZ()),
0409                        Plane::RotationType());
0410       propagate(currentState, *endcap);
0411     }
0412     if (currentState.isValid())
0413       trajectory.push_back(currentState);
0414   } else
0415     LogTrace("TrackAssociator") << "Weird message\n";
0416 
0417   while (currentState.isValid() && currentState.position().perp() < volume.maxR() &&
0418          fabs(currentState.position().z()) < volume.maxZ()) {
0419     propagateForward(currentState, step);
0420     if (!currentState.isValid()) {
0421       LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n";
0422       break;
0423     }
0424     // LogTrace("TrackAssociator") << "New state (perp, z): " << currentState.position().perp() << ", " << currentState.position().z();
0425     //if ( ( currentState.position().perp() < volume.maxR() && fabs(currentState.position().z()) < volume.maxZ() ) &&
0426     //     ( currentState.position().perp()-volume.minR() > 0  || fabs(currentState.position().z()) - volume.minZ() >0 ) )
0427     trajectory.push_back(currentState);
0428   }
0429 }
0430 
0431 void CachedTrajectory::reset_trajectory() {
0432   fullTrajectory_.clear();
0433   ecalTrajectory_.clear();
0434   hcalTrajectory_.clear();
0435   hoTrajectory_.clear();
0436   preshowerTrajectory_.clear();
0437   wideEcalTrajectory_.clear();
0438   wideHcalTrajectory_.clear();
0439   wideHOTrajectory_.clear();
0440   fullTrajectoryFilled_ = false;
0441 }
0442 
0443 void CachedTrajectory::findEcalTrajectory(const FiducialVolume& volume) {
0444   LogTrace("TrackAssociator") << "getting trajectory in ECAL";
0445   getTrajectory(ecalTrajectory_, volume, 4);
0446   LogTrace("TrackAssociator") << "# of points in ECAL trajectory:" << ecalTrajectory_.size();
0447 }
0448 
0449 void CachedTrajectory::findPreshowerTrajectory(const FiducialVolume& volume) {
0450   LogTrace("TrackAssociator") << "getting trajectory in Preshower";
0451   getTrajectory(preshowerTrajectory_, volume, 2);
0452   LogTrace("TrackAssociator") << "# of points in Preshower trajectory:" << preshowerTrajectory_.size();
0453 }
0454 
0455 const std::vector<SteppingHelixStateInfo>& CachedTrajectory::getEcalTrajectory() const { return ecalTrajectory_; }
0456 
0457 const std::vector<SteppingHelixStateInfo>& CachedTrajectory::getPreshowerTrajectory() const {
0458   return preshowerTrajectory_;
0459 }
0460 
0461 void CachedTrajectory::findHcalTrajectory(const FiducialVolume& volume) {
0462   LogTrace("TrackAssociator") << "getting trajectory in HCAL";
0463   getTrajectory(hcalTrajectory_, volume, 4);  // more steps to account for different depth
0464   LogTrace("TrackAssociator") << "# of points in HCAL trajectory:" << hcalTrajectory_.size();
0465 }
0466 
0467 const std::vector<SteppingHelixStateInfo>& CachedTrajectory::getHcalTrajectory() const { return hcalTrajectory_; }
0468 
0469 void CachedTrajectory::findHOTrajectory(const FiducialVolume& volume) {
0470   LogTrace("TrackAssociator") << "getting trajectory in HO";
0471   getTrajectory(hoTrajectory_, volume, 2);
0472   LogTrace("TrackAssociator") << "# of points in HO trajectory:" << hoTrajectory_.size();
0473 }
0474 
0475 const std::vector<SteppingHelixStateInfo>& CachedTrajectory::getHOTrajectory() const { return hoTrajectory_; }
0476 
0477 std::vector<GlobalPoint>* CachedTrajectory::getWideTrajectory(const std::vector<SteppingHelixStateInfo>& states,
0478                                                               WideTrajectoryType wideTrajectoryType) {
0479   std::vector<GlobalPoint>* wideTrajectory = nullptr;
0480   switch (wideTrajectoryType) {
0481     case Ecal:
0482       LogTrace("TrackAssociator") << "Filling ellipses in Ecal trajectory";
0483       wideTrajectory = &wideEcalTrajectory_;
0484       break;
0485     case Hcal:
0486       LogTrace("TrackAssociator") << "Filling ellipses in Hcal trajectory";
0487       wideTrajectory = &wideHcalTrajectory_;
0488       break;
0489     case HO:
0490       LogTrace("TrackAssociator") << "Filling ellipses in HO trajectory";
0491       wideTrajectory = &wideHOTrajectory_;
0492       break;
0493   }
0494   if (!wideTrajectory)
0495     return nullptr;
0496 
0497   for (std::vector<SteppingHelixStateInfo>::const_iterator state = states.begin(); state != states.end(); state++) {
0498     // defined a normal plane wrt the particle trajectory direction
0499     // let's hope that I computed the rotation matrix correctly.
0500     GlobalVector vector(state->momentum().unit());
0501     float r21 = 0;
0502     float r22 = vector.z() / sqrt(1 - pow(vector.x(), 2));
0503     float r23 = -vector.y() / sqrt(1 - pow(vector.x(), 2));
0504     float r31 = vector.x();
0505     float r32 = vector.y();
0506     float r33 = vector.z();
0507     float r11 = r22 * r33 - r23 * r32;
0508     float r12 = r23 * r31;
0509     float r13 = -r22 * r31;
0510 
0511     Plane::RotationType rotation(r11, r12, r13, r21, r22, r23, r31, r32, r33);
0512     Plane* target = new Plane(state->position(), rotation);
0513 
0514     TrajectoryStateOnSurface tsos = state->getStateOnSurface(*target);
0515 
0516     if (!tsos.isValid()) {
0517       LogTrace("TrackAssociator") << "[getWideTrajectory] TSOS not valid";
0518       continue;
0519     }
0520     if (!tsos.hasError()) {
0521       LogTrace("TrackAssociator") << "[getWideTrajectory] TSOS does not have Errors";
0522       continue;
0523     }
0524     LocalError localErr = tsos.localError().positionError();
0525     localErr.scale(2);  // get the 2 sigma ellipse
0526     float xx = localErr.xx();
0527     float xy = localErr.xy();
0528     float yy = localErr.yy();
0529 
0530     float denom = yy - xx;
0531     float phi = 0.;
0532     if (xy == 0 && denom == 0)
0533       phi = M_PI_4;
0534     else
0535       phi = 0.5 * atan2(2. * xy, denom);  // angle of MAJOR axis
0536     // Unrotate the error ellipse to get the semimajor and minor axes. Then place points on
0537     // the endpoints of semiminor an seminajor axes on original(rotated) error ellipse.
0538     LocalError rotErr = localErr.rotate(-phi);  // xy covariance of rotErr should be zero
0539     float semi1 = sqrt(rotErr.xx());
0540     float semi2 = sqrt(rotErr.yy());
0541 
0542     // Just use one point if the ellipse is small
0543     // if(semi1 < 0.1 && semi2 < 0.1) {
0544     //   LogTrace("TrackAssociator") << "[getWideTrajectory] Error ellipse is small, using one trajectory point";
0545     //   wideTrajectory->push_back(state->position());
0546     //   continue;
0547     // }
0548 
0549     Local2DPoint bounds[4];
0550     bounds[0] = Local2DPoint(semi1 * cos(phi), semi1 * sin(phi));
0551     bounds[1] = Local2DPoint(semi1 * cos(phi + M_PI), semi1 * sin(phi + M_PI));
0552     phi += M_PI_2;  // add pi/2 for the semi2 axis
0553     bounds[2] = Local2DPoint(semi2 * cos(phi), semi2 * sin(phi));
0554     bounds[3] = Local2DPoint(semi2 * cos(phi + M_PI), semi2 * sin(phi + M_PI));
0555 
0556     // LogTrace("TrackAssociator") << "Axes " << semi1 <<","<< semi2 <<"   phi "<< phi;
0557     // LogTrace("TrackAssociator") << "Local error ellipse: " << bounds[0] << bounds[1] << bounds[2] << bounds[3];
0558 
0559     wideTrajectory->push_back(state->position());
0560     for (int index = 0; index < 4; ++index)
0561       wideTrajectory->push_back(target->toGlobal(bounds[index]));
0562 
0563     //LogTrace("TrackAssociator") <<"Global error ellipse: (" << target->toGlobal(bounds[0]) <<","<< target->toGlobal(bounds[1])
0564     //         <<","<< target->toGlobal(bounds[2]) <<","<< target->toGlobal(bounds[3]) <<","<<state->position() <<")";
0565   }
0566 
0567   return wideTrajectory;
0568 }
0569 
0570 SteppingHelixStateInfo CachedTrajectory::getStateAtEcal() {
0571   if (ecalTrajectory_.empty())
0572     return SteppingHelixStateInfo();
0573   else
0574     return ecalTrajectory_.front();
0575 }
0576 
0577 SteppingHelixStateInfo CachedTrajectory::getStateAtPreshower() {
0578   if (preshowerTrajectory_.empty())
0579     return SteppingHelixStateInfo();
0580   else
0581     return preshowerTrajectory_.front();
0582 }
0583 
0584 SteppingHelixStateInfo CachedTrajectory::getStateAtHcal() {
0585   if (hcalTrajectory_.empty())
0586     return SteppingHelixStateInfo();
0587   else
0588     return hcalTrajectory_.front();
0589 }
0590 
0591 SteppingHelixStateInfo CachedTrajectory::getStateAtHO() {
0592   if (hoTrajectory_.empty())
0593     return SteppingHelixStateInfo();
0594   else
0595     return hoTrajectory_.front();
0596 }
0597 
0598 SteppingHelixStateInfo CachedTrajectory::getInnerState() {
0599   if (fullTrajectory_.empty())
0600     return SteppingHelixStateInfo();
0601   else
0602     return fullTrajectory_.front();
0603 }
0604 
0605 SteppingHelixStateInfo CachedTrajectory::getOuterState() {
0606   if (fullTrajectory_.empty())
0607     return SteppingHelixStateInfo();
0608   else
0609     return fullTrajectory_.back();
0610 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!