CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
SeedFromNuclearInteraction.cc
Go to the documentation of this file.
1 #include <memory>
2 
4 
7 
11 
13 
15  const TrackerGeometry* geom,
16  double iPtMin)
17  : ptMin(iPtMin), thePropagator(prop), theTrackerGeom(geom) {
18  isValid_ = true;
19  initialTSOS_ = std::make_shared<TrajectoryStateOnSurface>();
20  updatedTSOS_ = std::make_shared<TrajectoryStateOnSurface>();
21  freeTS_ = std::make_shared<FreeTrajectoryState>();
22 }
23 
24 //----------------------------------------------------------------------
26  ConstRecHitPointer ihit,
27  ConstRecHitPointer ohit) {
28  // delete pointer to TrackingRecHits
29  theHits.clear();
30 
31  // get the inner and outer transient TrackingRecHits
32  innerHit_ = ihit;
33  outerHit_ = ohit;
34 
35  //theHits.push_back( inner_TM.recHit() ); // put temporarily - TODO: remove this line
36  theHits.push_back(outerHit_);
37 
38  initialTSOS_.reset(new TrajectoryStateOnSurface(inner_TSOS));
39 
40  // calculate the initial FreeTrajectoryState.
41  freeTS_.reset(stateWithError());
42 
43  // check transverse momentum
44  if (freeTS_->momentum().perp() < ptMin) {
45  isValid_ = false;
46  } else {
47  // convert freeTS_ into a persistent TSOS on the outer surface
48  isValid_ = construct();
49  }
50 }
51 //----------------------------------------------------------------------
53  const TSOS& inner_TSOS,
54  ConstRecHitPointer ihit,
55  ConstRecHitPointer ohit) {
56  // delete pointer to TrackingRecHits
57  theHits.clear();
58 
59  // get the inner and outer transient TrackingRecHits
60  innerHit_ = ihit;
61  outerHit_ = ohit;
62 
63  GlobalPoint innerPos =
64  theTrackerGeom->idToDet(innerHit_->geographicalId())->surface().toGlobal(innerHit_->localPosition());
65  GlobalPoint outerPos =
66  theTrackerGeom->idToDet(outerHit_->geographicalId())->surface().toGlobal(outerHit_->localPosition());
67 
68  TangentHelix helix(thePrimaryHelix, outerPos, innerPos);
69 
70  theHits.push_back(innerHit_);
71  theHits.push_back(outerHit_);
72 
73  initialTSOS_.reset(new TrajectoryStateOnSurface(inner_TSOS));
74 
75  // calculate the initial FreeTrajectoryState from the inner and outer TM assuming that the helix equation is already known.
76  freeTS_.reset(stateWithError(helix));
77 
78  if (freeTS_->momentum().perp() < ptMin) {
79  isValid_ = false;
80  } else {
81  // convert freeTS_ into a persistent TSOS on the outer surface
82  isValid_ = construct();
83  }
84 }
85 //----------------------------------------------------------------------
87  // Calculation of the helix assuming that the secondary track has the same direction
88  // than the primary track and pass through the inner and outer hits.
89  GlobalVector direction = initialTSOS_->globalDirection();
90  GlobalPoint inner = initialTSOS_->globalPosition();
91  TangentHelix helix(direction, inner, outerHitPosition());
92 
93  return stateWithError(helix);
94 }
95 //----------------------------------------------------------------------
97  // typedef TkRotation<float> Rotation;
98 
99  GlobalVector dirAtVtx = helix.directionAtVertex();
100  const MagneticField& mag = initialTSOS_->globalParameters().magneticField();
101 
102  // Get the global parameters of the trajectory
103  // we assume that the magnetic field at the vertex is equal to the magnetic field at the inner TM.
105  helix.vertexPoint(), dirAtVtx, helix.charge(mag.inTesla(helix.vertexPoint()).z()) / helix.rho(), 0, &mag);
106 
107  // Error matrix in a frame where z is in the direction of the track at the vertex
108  AlgebraicSymMatrix66 primaryError(initialTSOS_->cartesianError().matrix());
109  double p_max = initialTSOS_->globalParameters().momentum().mag();
110  AlgebraicMatrix33 rot = this->rotationMatrix(dirAtVtx);
111 
112  AlgebraicMatrix66 globalRotation;
113  globalRotation.Place_at(rot, 0, 0);
114  globalRotation.Place_at(rot, 3, 3);
115  AlgebraicSymMatrix66 primaryErrorInNewFrame = ROOT::Math::Similarity(globalRotation, primaryError);
116 
117  AlgebraicSymMatrix66 secondaryErrorInNewFrame = AlgebraicMatrixID();
118  double p_perp_max = 2; // energy max of a secondary track emited perpendicularly to the
119  // primary track is +/- 2 GeV
120  secondaryErrorInNewFrame(0, 0) = primaryErrorInNewFrame(0, 0) + helix.vertexError() * p_perp_max / p_max;
121  secondaryErrorInNewFrame(1, 1) = primaryErrorInNewFrame(1, 1) + helix.vertexError() * p_perp_max / p_max;
122  secondaryErrorInNewFrame(2, 2) = helix.vertexError() * helix.vertexError();
123  secondaryErrorInNewFrame(3, 3) = p_perp_max * p_perp_max;
124  secondaryErrorInNewFrame(4, 4) = p_perp_max * p_perp_max;
125  secondaryErrorInNewFrame(5, 5) = p_max * p_max;
126 
127  AlgebraicSymMatrix66 secondaryError = ROOT::Math::SimilarityT(globalRotation, secondaryErrorInNewFrame);
128 
129  return new FreeTrajectoryState(gtp, CartesianTrajectoryError(secondaryError));
130 }
131 
132 //----------------------------------------------------------------------
134  // loop on all hits in theHits
135  KFUpdator theUpdator;
136 
137  const TrackingRecHit* hit = nullptr;
138 
139  LogDebug("NuclearSeedGenerator") << "Seed ** initial state " << freeTS_->cartesianError().matrix();
140 
141  for (unsigned int iHit = 0; iHit < theHits.size(); iHit++) {
142  hit = theHits[iHit]->hit();
144  (iHit == 0)
147 
148  if (!state.isValid())
149  return false;
150 
152  updatedTSOS_.reset(new TrajectoryStateOnSurface(theUpdator.update(state, *tth)));
153  }
154 
155  LogDebug("NuclearSeedGenerator") << "Seed ** updated state " << updatedTSOS_->cartesianError().matrix();
156 
158  return true;
159 }
160 
161 //----------------------------------------------------------------------
163  recHitContainer _hits;
164  for (ConstRecHitContainer::const_iterator it = theHits.begin(); it != theHits.end(); it++) {
165  _hits.push_back(it->get()->hit()->clone());
166  }
167  return _hits;
168 }
169 //----------------------------------------------------------------------
172 
173  // z axis coincides with perp
174  GlobalVector zAxis = perp.unit();
175 
176  // x axis has no global Z component
177  GlobalVector xAxis;
178  if (zAxis.x() != 0 || zAxis.y() != 0) {
179  // precision is not an issue here, just protect against divizion by zero
180  xAxis = GlobalVector(-zAxis.y(), zAxis.x(), 0).unit();
181  } else { // perp coincides with global Z
182  xAxis = GlobalVector(1, 0, 0);
183  }
184 
185  // y axis obtained by cross product
186  GlobalVector yAxis(zAxis.cross(xAxis));
187 
188  result(0, 0) = xAxis.x();
189  result(0, 1) = xAxis.y();
190  result(0, 2) = xAxis.z();
191  result(1, 0) = yAxis.x();
192  result(1, 1) = yAxis.y();
193  result(1, 2) = yAxis.z();
194  result(2, 0) = zAxis.x();
195  result(2, 1) = zAxis.y();
196  result(2, 2) = zAxis.z();
197  return result;
198 }
const TrackerGeometry * theTrackerGeom
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:110
PropagationDirection direction() const
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepStd< double, 3, 3 > > AlgebraicMatrix33
GlobalVector directionAtVertex()
Definition: TangentHelix.cc:17
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
constexpr float ptMin
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
T y() const
Definition: PV3DBase.h:60
SeedFromNuclearInteraction(const Propagator *prop, const TrackerGeometry *geom, double ptMin)
void setMeasurements(const TSOS &tsosAtInteractionPoint, ConstRecHitPointer ihit, ConstRecHitPointer ohit)
Fill all data members from 2 TM&#39;s where the first one is supposed to be at the interaction point...
tuple result
Definition: mps_fire.py:311
std::shared_ptr< TSOS > initialTSOS_
void push_back(D *&d)
Definition: OwnVector.h:326
TrajectoryStateOnSurface update(const TrajectoryStateOnSurface &, const TrackingRecHit &) const override
Definition: KFUpdator.cc:177
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer
T z() const
Definition: PV3DBase.h:61
double vertexError()
Definition: TangentHelix.h:47
const TrackerGeomDet * idToDet(DetId) const override
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepStd< double, 6, 6 > > AlgebraicMatrix66
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
AlgebraicMatrix33 rotationMatrix(const GlobalVector &perp) const
Vector3DBase unit() const
Definition: Vector3DBase.h:54
std::shared_ptr< TSOS > updatedTSOS_
std::shared_ptr< FreeTrajectoryState > freeTS_
int charge(float magz)
Definition: TangentHelix.h:41
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:50
T perp() const
Magnitude of transverse component.
tuple zAxis
Definition: MetAnalyzer.py:57
GlobalPoint vertexPoint() const
Definition: TangentHelix.h:35
DetId geographicalId() const
double rho() const
Definition: TangentHelix.h:43
T x() const
Definition: PV3DBase.h:59
Global3DVector GlobalVector
Definition: GlobalVector.h:10
Basic3DVector unit() const
#define LogDebug(id)
FreeTrajectoryState * stateWithError() const