CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SeedFromGenericPairOrTriplet.cc
Go to the documentation of this file.
10 
12  const TrackerGeometry* geom,
13  const TransientTrackingRecHitBuilder* builder,
14  const Propagator* propagatorAlong,
15  const Propagator* propagatorOpposite,
16  const std::vector<int>& charges,
17  bool momFromPSet,
18  double errorRescaling):theMagfield(mf), theTracker(geom), theBuilder(builder), thePropagatorAlong(propagatorAlong), thePropagatorOpposite(propagatorOpposite), theSetMomentum(momFromPSet), theCharges(charges), theErrorRescaling(errorRescaling){}
19 
20 
21 std::vector<TrajectorySeed*> SeedFromGenericPairOrTriplet::seed(const SeedingHitSet& hits,
23  const NavigationDirection& seedDir,
24  const edm::EventSetup& iSetup){
25  std::vector<TrajectorySeed*> seeds;
26  if (hits.size() == 3) {
27  if(!theSetMomentum){
28  TrajectorySeed* seed = seedFromTriplet(hits, dir, seedDir, iSetup);
29  if (seed) seeds.push_back(seed);
30  } else {
31  for (std::vector<int>::const_iterator iCh = theCharges.begin(); iCh != theCharges.end(); iCh++){
32  TrajectorySeed* seed = seedFromTriplet(hits, dir, seedDir, iSetup, *iCh);
33  if(seed) seeds.push_back(seed);
34  }
35  }
36 
37  } else if (hits.size() == 2){
38  if(!theSetMomentum){
39  TrajectorySeed* seed = seedFromPair(hits, dir, seedDir);
40  if (seed) seeds.push_back(seed);
41  } else {
42  for (std::vector<int>::const_iterator iCh = theCharges.begin(); iCh != theCharges.end(); iCh++){
43  TrajectorySeed* seed = seedFromPair(hits, dir, seedDir, *iCh);
44  if (seed) seeds.push_back(seed);
45  }
46  }
47  } else {
48  throw cms::Exception("CombinatorialSeedGeneratorForCosmics") << " Wrong number of hits in Set: "
49  << hits.size() << ", should be 2 or 3 ";
50  }
51  return seeds;
52 }
53 
56  const NavigationDirection& seedDir,
57  const edm::EventSetup& iSetup, int charge) const{
58  if (hits.size() != 3) {
59  throw cms::Exception("CombinatorialSeedGeneratorForCosmics") <<
60  "call to SeedFromGenericPairOrTriplet::seedFromTriplet with " << hits.size() << " hits ";
61  }
62 
63  const TrackingRecHit* innerHit = hits[0]->hit();
64  const TrackingRecHit* middleHit = hits[1]->hit();
65  const TrackingRecHit* outerHit = hits[2]->hit();
66  GlobalPoint inner = theTracker->idToDet(innerHit->geographicalId() )->surface().toGlobal(innerHit->localPosition() );
67  GlobalPoint middle = theTracker->idToDet(middleHit->geographicalId())->surface().toGlobal(middleHit->localPosition());
68  GlobalPoint outer = theTracker->idToDet(outerHit->geographicalId() )->surface().toGlobal(outerHit->localPosition() );
69  if (theSetMomentum){
70  LogDebug("SeedFromGenericPairOrTriplet") <<
71  "Using the following hits: outer(r, phi, theta) " << outerHit->geographicalId().subdetId()<< " (" << outer.perp()
72  << ", " << outer.phi()
73  << "," << outer.theta()
74  << ") " << middleHit->geographicalId().subdetId() << " middle ("
75  << middle.perp() << ", "
76  << middle.phi() << ","
77  << middle.theta() << ") "
78  <<innerHit->geographicalId().subdetId() <<" inner ("
79  << inner.perp()
80  << ", " << inner.phi()
81  << "," << inner.theta() <<") "
82  << " (x, y, z) outer ("
83  << inner.x() << ", "
84  << inner.y() << ", "
85  << inner.z() << ") middle ("
86  << middle.x() << ", "
87  << middle.y() << ", "
88  << middle.z() << ")";
89  SeedingHitSet newSet;
90  if (seedDir == outsideIn){
91  newSet.add(hits[1]);
92  newSet.add(hits[2]);
93  } else {
94  newSet.add(hits[0]);
95  newSet.add(hits[1]);
96  }
97  if (!qualityFilter(hits)) return 0;
98  TrajectorySeed* seed = seedFromPair(newSet, dir, seedDir, charge);
99  return seed;
100 
101  }
102  GlobalPoint* firstPoint = 0;
103  GlobalPoint* secondPoint = 0;
104  GlobalPoint* thirdPoint = 0;
105  int momentumSign = 1;
106  //const TrackingRecHit* firstHit = 0;
107  //const TrackingRecHit* secondHit = 0;
108  //choose the prop dir and hit order accordingly to where the seed is made
109  std::vector<const TrackingRecHit*> trHits;
110  if (seedDir == outsideIn){
111  LogDebug("SeedFromGenericPairOrTriplet")
112  << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
113  firstPoint = &outer;
114  secondPoint = &middle;
115  thirdPoint = &inner;
116  trHits.push_back(outerHit);
117  trHits.push_back(middleHit);
118  //firstHit = outerHit;
119  //secondHit = middleHit;
120  } else {
121  LogDebug("SeedFromGenericPairOrTriplet")
122  << "Seed from outsideIn oppositeToMomentum OR insideOut alongMomentum";
123  firstPoint = &inner;
124  secondPoint = &middle;
125  thirdPoint = &outer;
126  trHits.push_back(innerHit);
127  trHits.push_back(middleHit);
128  //firstHit = innerHit;
129  //secondHit = middleHit;
130  }
131  if (dir == oppositeToMomentum) momentumSign = -1;
132  FastHelix helix(*thirdPoint, *secondPoint, *firstPoint, iSetup);
133  FreeTrajectoryState originalStartingState = helix.stateAtVertex();
134  LogDebug("SeedFromGenericPairOrTriplet") << "originalStartingState " << originalStartingState;
135  /*GlobalTrajectoryParameters originalPar = originalStartingState.parameters();
136  GlobalTrajectoryParameters newPar = GlobalTrajectoryParameters(originalPar.position(),
137  momentumSign*originalPar.momentum(),
138  originalPar.charge(),
139  &originalPar.magneticField());
140  */
141 
142  /*FastCircle helix(*thirdPoint, *secondPoint, *firstPoint);
143  GlobalTrajectoryParameters newPar = GlobalTrajectoryParameters(*secondPoint,
144  momentumSign*originalPar.momentum(),
145  originalPar.charge(),
146  &originalPar.magneticField());*/
147  //FreeTrajectoryState* startingState = new FreeTrajectoryState(newPar, initialError(trHits[0]));//trHits[1]));
148  if (!qualityFilter(momentumSign*originalStartingState.momentum())) return 0;
149  //TrajectorySeed* seed = buildSeed(startingState, trHits, dir);
150  TrajectorySeed* seed = buildSeed(momentumSign*originalStartingState.momentum(),originalStartingState.charge(), trHits, dir);
151  return seed;
152 }
153 
155  const PropagationDirection& dir,
156  const NavigationDirection& seedDir, int charge) const{
157  if (hits.size() != 2) {
158  throw cms::Exception("CombinatorialSeedGeneratorForCosmics") <<
159  "call to SeedFromGenericPairOrTriplet::seedFromPair with " << hits.size() << " hits ";
160  }
161  const TrackingRecHit* innerHit = hits[0]->hit();
162  const TrackingRecHit* outerHit = hits[1]->hit();
163  GlobalPoint inner = theTracker->idToDet(innerHit->geographicalId() )->surface().toGlobal(innerHit->localPosition() );
164  GlobalPoint outer = theTracker->idToDet(outerHit->geographicalId() )->surface().toGlobal(outerHit->localPosition() );
165  LogDebug("SeedFromGenericPairOrTriplet") <<
166  "Using the following hits: outer(r, phi, theta) (" << outer.perp()
167  << ", " << outer.phi()
168  << "," << outer.theta()
169  <<") inner ("
170  << inner.perp()
171  << ", " << inner.phi()
172  << "," << inner.theta() <<")";
173  GlobalPoint* firstPoint = 0;
174  GlobalPoint* secondPoint = 0;
175  int momentumSign = 1;
176  //const TrackingRecHit* firstHit = 0;
177  //const TrackingRecHit* secondHit = 0;
178  std::vector<const TrackingRecHit*> trHits;
179  //choose the prop dir and hit order accordingly to where the seed is made
180  if (seedDir == outsideIn){
181  LogDebug("SeedFromGenericPairOrTriplet")
182  << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
183  firstPoint = &outer;
184  secondPoint = &inner;
185  //firstHit = outerHit;
186  //secondHit = innerHit;
187  trHits.push_back(outerHit);
188  trHits.push_back(innerHit);
189  } else {
190  LogDebug("SeedFromGenericPairOrTriplet")
191  << "Seed from outsideIn oppositeToMomentum OR insideOut alongMomentum";
192  firstPoint = &inner;
193  secondPoint = &outer;
194  momentumSign = -1;
195  //firstHit = innerHit;
196  //secondHit = outerHit;
197  trHits.push_back(innerHit);
198  trHits.push_back(outerHit);
199  }
200  if (dir == oppositeToMomentum) momentumSign = -1;
201  GlobalVector momentum = momentumSign*theP*(*secondPoint-*firstPoint).unit();
202  /*GlobalTrajectoryParameters gtp(*firstPoint,
203  momentum,
204  charge,
205  &(*theMagfield));*/
206  //FreeTrajectoryState* startingState = new FreeTrajectoryState(gtp,initialError(trHits[0]));//trHits[1]));
207  //if (!qualityFilter(startingState, hits)) return 0;
208  //TrajectorySeed* seed = buildSeed(startingState, firstHit, dir);
209  //TrajectorySeed* seed = buildSeed(startingState, trHits, dir);
210  TrajectorySeed* seed = buildSeed(momentum, charge, trHits, dir);
211  return seed;
212 }
213 
214 
216  int charge,
217  //const TrackingRecHit* firsthit,
218  std::vector<const TrackingRecHit*>& trHits,
219  const PropagationDirection& dir) const {
221  if (dir == oppositeToMomentum) propagator = thePropagatorOpposite;
222 
223  //debug
224  GlobalPoint first = theTracker->idToDet(trHits[0]->geographicalId() )->surface().toGlobal(trHits[0]->localPosition() );
225  GlobalPoint second = theTracker->idToDet(trHits[1]->geographicalId() )->surface().toGlobal(trHits[1]->localPosition() );
226  LogDebug("SeedFromGenericPairOrTriplet") <<
227  "Using the following hits: first(r, phi, theta) (" << first.perp()
228  << ", " << first.phi()
229  << "," << first.theta()
230  <<") second ("
231  << second.perp()
232  << ", " << second.phi()
233  << "," << second.theta() <<")";
234  //build an initial trajectory state with big errors,
235  //then update it with the first hit
237  LocalVector lmom = transHit->surface()->toLocal(momentum);
238  LocalTrajectoryParameters ltp(charge/momentum.mag(),
239  lmom.x()/lmom.z(),lmom.y()/lmom.z(),
240  transHit->localPosition().x(),
241  transHit->localPosition().y(),
242  lmom.z()>0.?1.:-1.);
243 
244  LocalTrajectoryError lte(100.,100.,
245  (1.+(lmom.x()/lmom.z())*(lmom.x()/lmom.z())),
246  (1.+(lmom.y()/lmom.z())*(lmom.y()/lmom.z())),
247  1./momentum.mag());
248 
249  TrajectoryStateOnSurface initialState(ltp,lte,*transHit->surface(),&(*theMagfield));
250  KFUpdator updator;
251  TrajectoryStateOnSurface startingState = updator.update(initialState,*transHit);
252 
254 
255  //TrajectoryStateOnSurface seedTSOS(*startingState, *plane);
256  const TrajectoryStateOnSurface propTSOS = propagator->propagate(startingState,*(transHit2->surface()));
257  if (!propTSOS.isValid()){
258  LogDebug("SeedFromGenericPairOrTriplet") << "first propagation failed";
259  return 0;
260  }
261  TrajectoryStateOnSurface seedTSOS = updator.update(propTSOS, *transHit2);
262  if (!seedTSOS.isValid()){
263  LogDebug("SeedFromGenericPairOrTriplet") << "first update failed";
264  return 0;
265  }
266  LogDebug("SeedFromGenericPairOrTriplet") << "starting TSOS " << seedTSOS ;
268  PTrajectoryStateOnDet *PTraj=
269  theTransformer.persistentState(seedTSOS, trHits[1]->geographicalId().rawId());
271  //build the transientTrackingRecHit for the starting hit, then call the method clone to rematch if needed
272  std::vector<const TrackingRecHit*>::const_iterator ihits;
273  for (ihits = trHits.begin(); ihits != trHits.end(); ihits++){
274  seed_hits.push_back((*ihits)->clone());
275  }
276  TrajectorySeed* seed = new TrajectorySeed(*PTraj,seed_hits,dir);
277  delete PTraj;
278  return seed;
279 }
280 /*
281 CurvilinearTrajectoryError SeedFromGenericPairOrTriplet::initialError(const TrackingRecHit* rechit) {
282  TransientTrackingRecHit::RecHitPointer transHit = theBuilder->build(rechit);
283  //AlgebraicSymMatrix C(5,1);
284  AlgebraicSymMatrix55 C = AlgebraicMatrixID();
285  //C*=0.1;
286  if (theSetMomentum){
287  C(0,0)=1/theP;
288  C(3,3)=transHit->globalPositionError().cxx();
289  C(4,4)=transHit->globalPositionError().czz();
290 
291  //C[0][0]=1/theP;
292  //C[3][3]=transHit->globalPositionError().cxx();
293  //C[4][4]=transHit->globalPositionError().czz();
294  } else {
295  float zErr = transHit->globalPositionError().czz();
296  float transverseErr = transHit->globalPositionError().cxx(); // assume equal cxx cyy
297  C(3,3) = transverseErr;
298  C(4,4) = zErr;
299  //C[3][3] = transverseErr;
300  //C[4][4] = zErr;
301  }
302 
303  return CurvilinearTrajectoryError(C);
304 }
305 */
307  //if we are setting the momentum with the PSet we look for aligned hits
308  if (theSetMomentum){
309  if (hits.size()==3){
310  std::vector<GlobalPoint> gPoints;
311  unsigned int nHits = hits.size();
312  for (unsigned int iHit = 0; iHit < nHits; ++iHit) gPoints.push_back(hits[iHit]->globalPosition());
313  unsigned int subid=(*hits[0]).geographicalId().subdetId();
314  if(subid == StripSubdetector::TEC || subid == StripSubdetector::TID){
315  LogDebug("SeedFromGenericPairOrTriplet")
316  << "In the endcaps we cannot decide if hits are aligned with only phi and z";
317  return true;
318  }
319  FastCircle circle(gPoints[0],
320  gPoints[1],
321  gPoints[2]);
322  if (circle.rho() < 500 && circle.rho() != 0) {
323  LogDebug("SeedFromGenericPairOrTriplet") <<
324  "Seed qualityFilter rejected because rho = " << circle.rho();
325  return false;
326  }
327 
328  }
329  }
330  return true;
331 }
332 
334  if (!theSetMomentum){
335  if (momentum.perp() < theP) return false;
336  }
337  return true;
338 }
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:78
FTS stateAtVertex() const
Definition: FastHelix.cc:7
T perp() const
Definition: PV3DBase.h:66
void update(const LocalTrajectoryParameters &p, const Surface &aSurface, const MagneticField *field, const SurfaceSide side=SurfaceSideDefinition::atCenterOfSurface)
TrajectorySeed * seedFromPair(const SeedingHitSet &hits, const PropagationDirection &dir, const NavigationDirection &seedDir, int charge=-1) const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
T y() const
Definition: PV3DBase.h:57
PropagationDirection
TrackCharge charge() const
double charge(const std::vector< uint8_t > &Ampls)
double rho() const
Definition: FastCircle.h:54
Geom::Theta< T > theta() const
Definition: PV3DBase.h:69
U second(std::pair< T, U > const &p)
void push_back(D *&d)
Definition: OwnVector.h:290
T mag() const
Definition: PV3DBase.h:61
const TransientTrackingRecHitBuilder * theBuilder
string unit
Definition: csvLumiCalc.py:46
virtual RecHitPointer build(const TrackingRecHit *p) const =0
build a tracking rechit from an existing rechit
T z() const
Definition: PV3DBase.h:58
TrajectorySeed * seedFromTriplet(const SeedingHitSet &hits, const PropagationDirection &dir, const NavigationDirection &seedDir, const edm::EventSetup &iSetup, int charge=-1) const
PTrajectoryStateOnDet * persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid) const
bool first
Definition: L1TdeRCT.cc:79
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
virtual const GeomDet * idToDet(DetId) const
GlobalVector momentum() const
TrajectorySeed * buildSeed(const GlobalVector &momentum, int charge, std::vector< const TrackingRecHit * > &trHits, const PropagationDirection &dir) const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:9
std::vector< TrajectorySeed * > seed(const SeedingHitSet &hits, const PropagationDirection &dir, const NavigationDirection &seedDir, const edm::EventSetup &iSetup)
SeedFromGenericPairOrTriplet(const MagneticField *mf, const TrackerGeometry *geom, const TransientTrackingRecHitBuilder *builder, const Propagator *propagatorAlong, const Propagator *propagatorOpposite, const std::vector< int > &charges, bool momFromPSet, double errorRescaling)
bool qualityFilter(const SeedingHitSet &hits) const
unsigned int size() const
Definition: SeedingHitSet.h:14
void add(TransientTrackingRecHit::ConstRecHitPointer pHit)
Definition: SeedingHitSet.h:15
DetId geographicalId() const
dbl *** dir
Definition: mlp_gen.cc:35
T x() const
Definition: PV3DBase.h:56
virtual LocalPoint localPosition() const =0
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37