CMS 3D CMS Logo
RPCSeedPattern.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  *
5  * \author Haiyun.Teng - Peking University
6  *
7  */
8 
9 #include "RecoMuon/MuonSeedGenerator/src/RPCSeedPattern.h"
10 #include <RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h>
11 #include <MagneticField/Engine/interface/MagneticField.h>
12 #include <MagneticField/Records/interface/IdealMagneticFieldRecord.h>
13 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
14 #include <TrackingTools/DetLayers/interface/DetLayer.h>
15 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
16 #include <DataFormats/Common/interface/OwnVector.h>
17 #include <DataFormats/MuonDetId/interface/RPCDetId.h>
18 #include <FWCore/Framework/interface/ESHandle.h>
19 #include <FWCore/MessageLogger/interface/MessageLogger.h>
20 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
21 #include <Geometry/RPCGeometry/interface/RPCChamber.h>
22 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
23 #include <Geometry/Records/interface/MuonGeometryRecord.h>
24 
25 #include "gsl/gsl_statistics.h"
26 #include "TH1F.h"
27 #include <cmath>
28 
29 using namespace std;
30 using namespace edm;
31 
32 RPCSeedPattern::RPCSeedPattern() {
33  isPatternChecked = false;
34  isConfigured = false;
35  MagnecticFieldThreshold = 0.5;
36 }
37 
38 RPCSeedPattern::~RPCSeedPattern() {}
39 
40 void RPCSeedPattern::configure(const edm::ParameterSet& iConfig) {
41  MaxRSD = iConfig.getParameter<double>("MaxRSD");
42  deltaRThreshold = iConfig.getParameter<double>("deltaRThreshold");
43  AlgorithmType = iConfig.getParameter<unsigned int>("AlgorithmType");
44  autoAlgorithmChoose = iConfig.getParameter<bool>("autoAlgorithmChoose");
45  ZError = iConfig.getParameter<double>("ZError");
46  MinDeltaPhi = iConfig.getParameter<double>("MinDeltaPhi");
47  MagnecticFieldThreshold = iConfig.getParameter<double>("MagnecticFieldThreshold");
48  stepLength = iConfig.getParameter<double>("stepLength");
49  sampleCount = iConfig.getParameter<unsigned int>("sampleCount");
50  isConfigured = true;
51 }
52 
53 RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::seed(const edm::EventSetup& eSetup, int& isGoodSeed) {
54  if (isConfigured == false) {
55  cout << "Configuration not set yet" << endl;
56  return createFakeSeed(isGoodSeed, eSetup);
57  }
58 
59  // Check recHit number, if fail we return a fake seed and set pattern to "wrong"
60  unsigned int NumberofHitsinSeed = nrhit();
61  if (NumberofHitsinSeed < 3)
62  return createFakeSeed(isGoodSeed, eSetup);
63  // If only three recHits, we don't have other choice
64  if (NumberofHitsinSeed == 3)
65  ThreePointsAlgorithm();
66 
67  if (NumberofHitsinSeed > 3) {
68  if (autoAlgorithmChoose == false) {
69  cout << "computePtWithmorerecHits" << endl;
70  if (AlgorithmType == 0)
71  ThreePointsAlgorithm();
72  if (AlgorithmType == 1)
73  MiddlePointsAlgorithm();
74  if (AlgorithmType == 2)
75  SegmentAlgorithm();
76  if (AlgorithmType == 3) {
77  if (checkSegment())
78  SegmentAlgorithmSpecial(eSetup);
79  else {
80  cout << "Not enough recHits for Special Segment Algorithm" << endl;
81  return createFakeSeed(isGoodSeed, eSetup);
82  }
83  }
84  } else {
85  if (checkSegment()) {
86  AlgorithmType = 3;
87  SegmentAlgorithmSpecial(eSetup);
88  } else {
89  AlgorithmType = 2;
90  SegmentAlgorithm();
91  }
92  }
93  }
94 
95  // Check the pattern
96  if (isPatternChecked == false) {
97  if (AlgorithmType != 3) {
98  checkSimplePattern(eSetup);
99  } else {
100  checkSegmentAlgorithmSpecial(eSetup);
101  }
102  }
103 
104  return createSeed(isGoodSeed, eSetup);
105 }
106 
107 void RPCSeedPattern::ThreePointsAlgorithm() {
108  cout << "computePtWith3recHits" << endl;
109  unsigned int NumberofHitsinSeed = nrhit();
110  // Check recHit number, if fail we set the pattern to "wrong"
111  if (NumberofHitsinSeed < 3) {
112  isPatternChecked = true;
113  isGoodPattern = -1;
114  return;
115  }
116  // Choose every 3 recHits to form a part
117  unsigned int NumberofPart = NumberofHitsinSeed * (NumberofHitsinSeed - 1) * (NumberofHitsinSeed - 2) / (3 * 2);
118  ;
119  double* pt = new double[NumberofPart];
120  double* pt_err = new double[NumberofPart];
121  // Loop for each three-recHits part
122  ConstMuonRecHitPointer precHit[3];
123  unsigned int n = 0;
124  unsigned int NumberofStraight = 0;
125  for (unsigned int i = 0; i < (NumberofHitsinSeed - 2); i++)
126  for (unsigned int j = (i + 1); j < (NumberofHitsinSeed - 1); j++)
127  for (unsigned int k = (j + 1); k < NumberofHitsinSeed; k++) {
128  precHit[0] = theRecHits[i];
129  precHit[1] = theRecHits[j];
130  precHit[2] = theRecHits[k];
131  bool checkStraight = checkStraightwithThreerecHits(precHit, MinDeltaPhi);
132  if (!checkStraight) {
133  GlobalVector Center_temp = computePtwithThreerecHits(pt[n], pt_err[n], precHit);
134  // For simple pattern
135  Center += Center_temp;
136  } else {
137  // For simple pattern
138  NumberofStraight++;
139  pt[n] = upper_limit_pt;
140  pt_err[n] = 0;
141  }
142  n++;
143  }
144  // For simple pattern, only one general parameter for pattern
145  if (NumberofStraight == NumberofPart) {
146  isStraight = true;
147  meanRadius = -1;
148  } else {
149  isStraight = false;
150  Center /= (NumberofPart - NumberofStraight);
151  double meanR = 0.;
152  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
153  meanR += getDistance(*iter, Center);
154  meanR /= NumberofHitsinSeed;
155  meanRadius = meanR;
156  }
157 
158  // Unset the pattern estimation signa
159  isPatternChecked = false;
160 
161  //double ptmean0 = 0;
162  //double sptmean0 = 0;
163  //computeBestPt(pt, pt_err, ptmean0, sptmean0, (NumberofPart - NumberofStraight));
164 
165  delete[] pt;
166  delete[] pt_err;
167 }
168 
169 void RPCSeedPattern::MiddlePointsAlgorithm() {
170  cout << "Using middle points algorithm" << endl;
171  unsigned int NumberofHitsinSeed = nrhit();
172  // Check recHit number, if fail we set the pattern to "wrong"
173  if (NumberofHitsinSeed < 4) {
174  isPatternChecked = true;
175  isGoodPattern = -1;
176  return;
177  }
178  double* X = new double[NumberofHitsinSeed];
179  double* Y = new double[NumberofHitsinSeed];
180  unsigned int n = 0;
181  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
182  X[n] = (*iter)->globalPosition().x();
183  Y[n] = (*iter)->globalPosition().y();
184  cout << "X[" << n << "] = " << X[n] << ", Y[" << n << "]= " << Y[n] << endl;
185  n++;
186  }
187  unsigned int NumberofPoints = NumberofHitsinSeed;
188  while (NumberofPoints > 3) {
189  for (unsigned int i = 0; i <= (NumberofPoints - 2); i++) {
190  X[i] = (X[i] + X[i + 1]) / 2;
191  Y[i] = (Y[i] + Y[i + 1]) / 2;
192  }
193  NumberofPoints--;
194  }
195  double x[3], y[3];
196  for (unsigned int i = 0; i < 3; i++) {
197  x[i] = X[i];
198  y[i] = Y[i];
199  }
200  double pt = 0;
201  double pt_err = 0;
202  bool checkStraight = checkStraightwithThreerecHits(x, y, MinDeltaPhi);
203  if (!checkStraight) {
204  GlobalVector Center_temp = computePtWithThreerecHits(pt, pt_err, x, y);
205  double meanR = 0.;
206  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
207  meanR += getDistance(*iter, Center_temp);
208  meanR /= NumberofHitsinSeed;
209  // For simple pattern
210  isStraight = false;
211  Center = Center_temp;
212  meanRadius = meanR;
213  } else {
214  // For simple pattern
215  isStraight = true;
216  meanRadius = -1;
217  }
218 
219  // Unset the pattern estimation signa
220  isPatternChecked = false;
221 
222  delete[] X;
223  delete[] Y;
224 }
225 
226 void RPCSeedPattern::SegmentAlgorithm() {
227  cout << "Using segments algorithm" << endl;
228  unsigned int NumberofHitsinSeed = nrhit();
229  // Check recHit number, if fail we set the pattern to "wrong"
230  if (NumberofHitsinSeed < 4) {
231  isPatternChecked = true;
232  isGoodPattern = -1;
233  return;
234  }
235 
236  RPCSegment* Segment;
237  unsigned int NumberofSegment = NumberofHitsinSeed - 2;
238  Segment = new RPCSegment[NumberofSegment];
239  unsigned int n = 0;
240  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 2); iter++) {
241  Segment[n].first = (*iter);
242  Segment[n].second = (*(iter + 2));
243  n++;
244  }
245  unsigned int NumberofStraight = 0;
246  for (unsigned int i = 0; i < NumberofSegment - 1; i++) {
247  bool checkStraight = checkStraightwithSegment(Segment[i], Segment[i + 1], MinDeltaPhi);
248  if (checkStraight == true) {
249  // For simple patterm
250  NumberofStraight++;
251  } else {
252  GlobalVector Center_temp = computePtwithSegment(Segment[i], Segment[i + 1]);
253  // For simple patterm
254  Center += Center_temp;
255  }
256  }
257  // For simple pattern, only one general parameter for pattern
258  if ((NumberofSegment - 1 - NumberofStraight) > 0) {
259  isStraight = false;
260  Center /= (NumberofSegment - 1 - NumberofStraight);
261  double meanR = 0.;
262  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
263  meanR += getDistance(*iter, Center);
264  meanR /= NumberofHitsinSeed;
265  meanRadius = meanR;
266  } else {
267  isStraight = true;
268  meanRadius = -1;
269  }
270 
271  // Unset the pattern estimation signal
272  isPatternChecked = false;
273 
274  delete[] Segment;
275 }
276 
277 void RPCSeedPattern::SegmentAlgorithmSpecial(const edm::EventSetup& eSetup) {
278  // Get magnetic field
279  edm::ESHandle<MagneticField> Field;
280  eSetup.get<IdealMagneticFieldRecord>().get(Field);
281 
282  //unsigned int NumberofHitsinSeed = nrhit();
283  if (!checkSegment()) {
284  isPatternChecked = true;
285  isGoodPattern = -1;
286  return;
287  }
288 
289  // Get magnetice field sampling information, recHit's position is not the border of Chamber and Iron
290  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
291  GlobalPoint gpFirst = (*iter)->globalPosition();
292  GlobalPoint gpLast = (*(iter + 1))->globalPosition();
293  GlobalPoint* gp = new GlobalPoint[sampleCount];
294  double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
295  double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
296  double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
297  for (unsigned int index = 0; index < sampleCount; index++) {
298  gp[index] = GlobalPoint(
299  (gpFirst.x() + dx * (index + 1)), (gpFirst.y() + dy * (index + 1)), (gpFirst.z() + dz * (index + 1)));
300  GlobalVector MagneticVec_temp = Field->inTesla(gp[index]);
301  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
302  //BValue.push_back(MagneticVec_temp);
303  }
304  delete[] gp;
305  }
306 
307  // form two segments
308  ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin();
309  for (unsigned int n = 0; n <= 1; n++) {
310  SegmentRB[n].first = (*iter);
311  cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
312  iter++;
313  SegmentRB[n].second = (*iter);
314  cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
315  iter++;
316  }
317  GlobalVector segvec1 = (SegmentRB[0].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
318  GlobalVector segvec2 = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
319 
320  // extrapolate the segment to find the Iron border which magnetic field is at large value
321  entryPosition = (SegmentRB[0].second)->globalPosition();
322  leavePosition = (SegmentRB[1].first)->globalPosition();
323  while (fabs(Field->inTesla(entryPosition).z()) < MagnecticFieldThreshold) {
324  cout << "Entry position is : " << entryPosition << ", and stepping into next point" << endl;
325  entryPosition += segvec1.unit() * stepLength;
326  }
327  // Loop back for more accurate by stepLength/10
328  while (fabs(Field->inTesla(entryPosition).z()) >= MagnecticFieldThreshold) {
329  cout << "Entry position is : " << entryPosition << ", and stepping back into next point" << endl;
330  entryPosition -= segvec1.unit() * stepLength / 10;
331  }
332  entryPosition += 0.5 * segvec1.unit() * stepLength / 10;
333  cout << "Final entry position is : " << entryPosition << endl;
334 
335  while (fabs(Field->inTesla(leavePosition).z()) < MagnecticFieldThreshold) {
336  cout << "Leave position is : " << leavePosition << ", and stepping into next point" << endl;
337  leavePosition -= segvec2.unit() * stepLength;
338  }
339  // Loop back for more accurate by stepLength/10
340  while (fabs(Field->inTesla(leavePosition).z()) >= MagnecticFieldThreshold) {
341  cout << "Leave position is : " << leavePosition << ", and stepping back into next point" << endl;
342  leavePosition += segvec2.unit() * stepLength / 10;
343  }
344  leavePosition -= 0.5 * segvec2.unit() * stepLength / 10;
345  cout << "Final leave position is : " << leavePosition << endl;
346 
347  // Sampling magnetic field in Iron region
348  GlobalPoint* gp = new GlobalPoint[sampleCount];
349  double dx = (leavePosition.x() - entryPosition.x()) / (sampleCount + 1);
350  double dy = (leavePosition.y() - entryPosition.y()) / (sampleCount + 1);
351  double dz = (leavePosition.z() - entryPosition.z()) / (sampleCount + 1);
352  std::vector<GlobalVector> BValue;
353  BValue.clear();
354  for (unsigned int index = 0; index < sampleCount; index++) {
355  gp[index] = GlobalPoint((entryPosition.x() + dx * (index + 1)),
356  (entryPosition.y() + dy * (index + 1)),
357  (entryPosition.z() + dz * (index + 1)));
358  GlobalVector MagneticVec_temp = Field->inTesla(gp[index]);
359  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
360  BValue.push_back(MagneticVec_temp);
361  }
362  delete[] gp;
363  GlobalVector meanB2(0, 0, 0);
364  for (std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
365  meanB2 += (*BIter);
366  meanB2 /= BValue.size();
367  cout << "Mean B field is " << meanB2 << endl;
368  meanMagneticField2 = meanB2;
369 
370  double meanBz2 = meanB2.z();
371  double deltaBz2 = 0.;
372  for (std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
373  deltaBz2 += (BIter->z() - meanBz2) * (BIter->z() - meanBz2);
374  deltaBz2 /= BValue.size();
375  deltaBz2 = sqrt(deltaBz2);
376  cout << "delta Bz is " << deltaBz2 << endl;
377 
378  // Distance of the initial 3 segment
379  S = 0;
380  bool checkStraight = checkStraightwithSegment(SegmentRB[0], SegmentRB[1], MinDeltaPhi);
381  if (checkStraight == true) {
382  // Just for complex pattern
383  isStraight2 = checkStraight;
384  Center2 = GlobalVector(0, 0, 0);
385  meanRadius2 = -1;
386  GlobalVector MomentumVec = (SegmentRB[1].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
387  S += MomentumVec.perp();
388  lastPhi = MomentumVec.phi().value();
389  } else {
390  GlobalVector seg1 = entryPosition - (SegmentRB[0].first)->globalPosition();
391  S += seg1.perp();
392  GlobalVector seg2 = (SegmentRB[1].second)->globalPosition() - leavePosition;
393  S += seg2.perp();
394  GlobalVector vecZ(0, 0, 1);
395  GlobalVector gvec1 = seg1.cross(vecZ);
396  GlobalVector gvec2 = seg2.cross(vecZ);
397  double A1 = gvec1.x();
398  double B1 = gvec1.y();
399  double A2 = gvec2.x();
400  double B2 = gvec2.y();
401  double X1 = entryPosition.x();
402  double Y1 = entryPosition.y();
403  double X2 = leavePosition.x();
404  double Y2 = leavePosition.y();
405  double XO = (A1 * A2 * (Y2 - Y1) + A2 * B1 * X1 - A1 * B2 * X2) / (A2 * B1 - A1 * B2);
406  double YO = (B1 * B2 * (X2 - X1) + B2 * A1 * Y1 - B1 * A2 * Y2) / (B2 * A1 - B1 * A2);
407  GlobalVector Center_temp(XO, YO, 0);
408  // Just for complex pattern
409  isStraight2 = checkStraight;
410  Center2 = Center_temp;
411 
412  cout << "entryPosition: " << entryPosition << endl;
413  cout << "leavePosition: " << leavePosition << endl;
414  cout << "Center2 is : " << Center_temp << endl;
415 
416  double R1 = GlobalVector((entryPosition.x() - Center_temp.x()),
417  (entryPosition.y() - Center_temp.y()),
418  (entryPosition.z() - Center_temp.z()))
419  .perp();
420  double R2 = GlobalVector((leavePosition.x() - Center_temp.x()),
421  (leavePosition.y() - Center_temp.y()),
422  (leavePosition.z() - Center_temp.z()))
423  .perp();
424  double meanR = (R1 + R2) / 2;
425  double deltaR = sqrt(((R1 - meanR) * (R1 - meanR) + (R2 - meanR) * (R2 - meanR)) / 2);
426  meanRadius2 = meanR;
427  cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
428  cout << "Mean radius is " << meanR << endl;
429  cout << "Delta R is " << deltaR << endl;
430  double deltaPhi =
431  fabs(((leavePosition - GlobalPoint(XO, YO, 0)).phi() - (entryPosition - GlobalPoint(XO, YO, 0)).phi()).value());
432  S += meanR * deltaPhi;
433  lastPhi = seg2.phi().value();
434  }
435 
436  // Unset the pattern estimation signa
437  isPatternChecked = false;
438 }
439 
440 bool RPCSeedPattern::checkSegment() const {
441  bool isFit = true;
442  unsigned int count = 0;
443  // first 4 recHits should be located in RB1 and RB2
444  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
445  count++;
446  const GeomDet* Detector = (*iter)->det();
447  if (dynamic_cast<const RPCChamber*>(Detector) != nullptr) {
448  const RPCChamber* RPCCh = dynamic_cast<const RPCChamber*>(Detector);
449  RPCDetId RPCId = RPCCh->id();
450  int Region = RPCId.region();
451  unsigned int Station = RPCId.station();
452  //int Layer = RPCId.layer();
453  if (count <= 4) {
454  if (Region != 0)
455  isFit = false;
456  if (Station > 2)
457  isFit = false;
458  }
459  }
460  }
461  // more than 4 recHits for pattern building
462  if (count <= 4)
463  isFit = false;
464  cout << "Check for segment fit: " << isFit << endl;
465  return isFit;
466 }
467 
468 MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::FirstRecHit() const { return theRecHits.front(); }
469 
470 MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::BestRefRecHit() const {
471  ConstMuonRecHitPointer best;
472  int index = 0;
473  // Use the last one for recHit on last layer has minmum delta Z for barrel or delta R for endcap while calculating the momentum
474  // But for Algorithm 3 we use the 4th recHit on the 2nd segment for more accurate
475  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
476  if (AlgorithmType != 3)
477  best = (*iter);
478  else if (index < 4)
479  best = (*iter);
480  index++;
481  }
482  return best;
483 }
484 
485 double RPCSeedPattern::getDistance(const ConstMuonRecHitPointer& precHit, const GlobalVector& Center) const {
486  return sqrt((precHit->globalPosition().x() - Center.x()) * (precHit->globalPosition().x() - Center.x()) +
487  (precHit->globalPosition().y() - Center.y()) * (precHit->globalPosition().y() - Center.y()));
488 }
489 
490 bool RPCSeedPattern::checkStraightwithThreerecHits(ConstMuonRecHitPointer (&precHit)[3], double MinDeltaPhi) const {
491  GlobalVector segvec1 = precHit[1]->globalPosition() - precHit[0]->globalPosition();
492  GlobalVector segvec2 = precHit[2]->globalPosition() - precHit[1]->globalPosition();
493  double dPhi = (segvec2.phi() - segvec1.phi()).value();
494  if (fabs(dPhi) > MinDeltaPhi) {
495  cout << "Part is estimate to be not straight" << endl;
496  return false;
497  } else {
498  cout << "Part is estimate to be straight" << endl;
499  return true;
500  }
501 }
502 
503 GlobalVector RPCSeedPattern::computePtwithThreerecHits(double& pt,
504  double& pt_err,
505  ConstMuonRecHitPointer (&precHit)[3]) const {
506  double x[3], y[3];
507  x[0] = precHit[0]->globalPosition().x();
508  y[0] = precHit[0]->globalPosition().y();
509  x[1] = precHit[1]->globalPosition().x();
510  y[1] = precHit[1]->globalPosition().y();
511  x[2] = precHit[2]->globalPosition().x();
512  y[2] = precHit[2]->globalPosition().y();
513  double A = (y[2] - y[1]) / (x[2] - x[1]) - (y[1] - y[0]) / (x[1] - x[0]);
514  double TYO = (x[2] - x[0]) / A + (y[2] * y[2] - y[1] * y[1]) / ((x[2] - x[1]) * A) -
515  (y[1] * y[1] - y[0] * y[0]) / ((x[1] - x[0]) * A);
516  double TXO = (x[2] + x[1]) + (y[2] * y[2] - y[1] * y[1]) / (x[2] - x[1]) - TYO * (y[2] - y[1]) / (x[2] - x[1]);
517  double XO = 0.5 * TXO;
518  double YO = 0.5 * TYO;
519  double R2 = (x[0] - XO) * (x[0] - XO) + (y[0] - YO) * (y[0] - YO);
520  cout << "R2 is " << R2 << endl;
521  // How this algorithm get the pt without magnetic field??
522  pt = 0.01 * sqrt(R2) * 2 * 0.3;
523  cout << "pt is " << pt << endl;
524  GlobalVector Center(XO, YO, 0);
525  return Center;
526 }
527 
528 bool RPCSeedPattern::checkStraightwithSegment(const RPCSegment& Segment1,
529  const RPCSegment& Segment2,
530  double MinDeltaPhi) const {
531  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
532  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
533  GlobalVector segvec3 = (Segment2.first)->globalPosition() - (Segment1.first)->globalPosition();
534  // compare segvec 1&2 for paralle, 1&3 for straight
535  double dPhi1 = (segvec2.phi() - segvec1.phi()).value();
536  double dPhi2 = (segvec3.phi() - segvec1.phi()).value();
537  cout << "Checking straight with 2 segments. dPhi1: " << dPhi1 << ", dPhi2: " << dPhi2 << endl;
538  cout << "Checking straight with 2 segments. dPhi1 in degree: " << dPhi1 * 180 / 3.1415926
539  << ", dPhi2 in degree: " << dPhi2 * 180 / 3.1415926 << endl;
540  if (fabs(dPhi1) > MinDeltaPhi || fabs(dPhi2) > MinDeltaPhi) {
541  cout << "Segment is estimate to be not straight" << endl;
542  return false;
543  } else {
544  cout << "Segment is estimate to be straight" << endl;
545  return true;
546  }
547 }
548 
549 GlobalVector RPCSeedPattern::computePtwithSegment(const RPCSegment& Segment1, const RPCSegment& Segment2) const {
550  GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
551  GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
552  GlobalPoint Point1(((Segment1.second)->globalPosition().x() + (Segment1.first)->globalPosition().x()) / 2,
553  ((Segment1.second)->globalPosition().y() + (Segment1.first)->globalPosition().y()) / 2,
554  ((Segment1.second)->globalPosition().z() + (Segment1.first)->globalPosition().z()) / 2);
555  GlobalPoint Point2(((Segment2.second)->globalPosition().x() + (Segment2.first)->globalPosition().x()) / 2,
556  ((Segment2.second)->globalPosition().y() + (Segment2.first)->globalPosition().y()) / 2,
557  ((Segment2.second)->globalPosition().z() + (Segment2.first)->globalPosition().z()) / 2);
558  GlobalVector vecZ(0, 0, 1);
559  GlobalVector gvec1 = segvec1.cross(vecZ);
560  GlobalVector gvec2 = segvec2.cross(vecZ);
561  double A1 = gvec1.x();
562  double B1 = gvec1.y();
563  double A2 = gvec2.x();
564  double B2 = gvec2.y();
565  double X1 = Point1.x();
566  double Y1 = Point1.y();
567  double X2 = Point2.x();
568  double Y2 = Point2.y();
569  double XO = (A1 * A2 * (Y2 - Y1) + A2 * B1 * X1 - A1 * B2 * X2) / (A2 * B1 - A1 * B2);
570  double YO = (B1 * B2 * (X2 - X1) + B2 * A1 * Y1 - B1 * A2 * Y2) / (B2 * A1 - B1 * A2);
571  GlobalVector Center(XO, YO, 0);
572  return Center;
573 }
574 
575 bool RPCSeedPattern::checkStraightwithThreerecHits(double (&x)[3], double (&y)[3], double MinDeltaPhi) const {
576  GlobalVector segvec1((x[1] - x[0]), (y[1] - y[0]), 0);
577  GlobalVector segvec2((x[2] - x[1]), (y[2] - y[1]), 0);
578  double dPhi = (segvec2.phi() - segvec1.phi()).value();
579  if (fabs(dPhi) > MinDeltaPhi) {
580  cout << "Part is estimate to be not straight" << endl;
581  return false;
582  } else {
583  cout << "Part is estimate to be straight" << endl;
584  return true;
585  }
586 }
587 
588 GlobalVector RPCSeedPattern::computePtWithThreerecHits(double& pt,
589  double& pt_err,
590  double (&x)[3],
591  double (&y)[3]) const {
592  double A = (y[2] - y[1]) / (x[2] - x[1]) - (y[1] - y[0]) / (x[1] - x[0]);
593  double TYO = (x[2] - x[0]) / A + (y[2] * y[2] - y[1] * y[1]) / ((x[2] - x[1]) * A) -
594  (y[1] * y[1] - y[0] * y[0]) / ((x[1] - x[0]) * A);
595  double TXO = (x[2] + x[1]) + (y[2] * y[2] - y[1] * y[1]) / (x[2] - x[1]) - TYO * (y[2] - y[1]) / (x[2] - x[1]);
596  double XO = 0.5 * TXO;
597  double YO = 0.5 * TYO;
598  double R2 = (x[0] - XO) * (x[0] - XO) + (y[0] - YO) * (y[0] - YO);
599  cout << "R2 is " << R2 << endl;
600  // How this algorithm get the pt without magnetic field??
601  pt = 0.01 * sqrt(R2) * 2 * 0.3;
602  cout << "pt is " << pt << endl;
603  GlobalVector Center(XO, YO, 0);
604  return Center;
605 }
606 
607 void RPCSeedPattern::checkSimplePattern(const edm::EventSetup& eSetup) {
608  if (isPatternChecked == true)
609  return;
610 
611  // Get magnetic field
612  edm::ESHandle<MagneticField> Field;
613  eSetup.get<IdealMagneticFieldRecord>().get(Field);
614 
615  unsigned int NumberofHitsinSeed = nrhit();
616 
617  // Print the recHit's position
618  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
619  cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
620 
621  // Get magnetice field information
622  std::vector<double> BzValue;
623  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
624  GlobalPoint gpFirst = (*iter)->globalPosition();
625  GlobalPoint gpLast = (*(iter + 1))->globalPosition();
626  GlobalPoint* gp = new GlobalPoint[sampleCount];
627  double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
628  double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
629  double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
630  for (unsigned int index = 0; index < sampleCount; index++) {
631  gp[index] = GlobalPoint(
632  (gpFirst.x() + dx * (index + 1)), (gpFirst.y() + dy * (index + 1)), (gpFirst.z() + dz * (index + 1)));
633  GlobalVector MagneticVec_temp = Field->inTesla(gp[index]);
634  cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
635  BzValue.push_back(MagneticVec_temp.z());
636  }
637  delete[] gp;
638  }
639  meanBz = 0.;
640  for (unsigned int index = 0; index < BzValue.size(); index++)
641  meanBz += BzValue[index];
642  meanBz /= BzValue.size();
643  cout << "Mean Bz is " << meanBz << endl;
644  deltaBz = 0.;
645  for (unsigned int index = 0; index < BzValue.size(); index++)
646  deltaBz += (BzValue[index] - meanBz) * (BzValue[index] - meanBz);
647  deltaBz /= BzValue.size();
648  deltaBz = sqrt(deltaBz);
649  cout << "delata Bz is " << deltaBz << endl;
650 
651  // Set isGoodPattern to default true and check the failure
652  isGoodPattern = 1;
653 
654  // Check the Z direction
655  if (fabs((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError) {
656  if (((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
657  isParralZ = 1;
658  else
659  isParralZ = -1;
660  } else
661  isParralZ = 0;
662 
663  cout << " Check isParralZ is :" << isParralZ << endl;
664  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
665  if (isParralZ == 0) {
666  if (fabs((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) > ZError) {
667  cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
668  isGoodPattern = 0;
669  }
670  } else {
671  if ((int)(((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) / ZError) * isParralZ < 0) {
672  cout << "Pattern find error in Z direction: wrong Z direction" << endl;
673  isGoodPattern = 0;
674  }
675  }
676  }
677 
678  // Check pattern
679  if (isStraight == false) {
680  // Check clockwise direction
681  GlobalVector* vec = new GlobalVector[NumberofHitsinSeed];
682  unsigned int index = 0;
683  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
684  GlobalVector vec_temp(((*iter)->globalPosition().x() - Center.x()),
685  ((*iter)->globalPosition().y() - Center.y()),
686  ((*iter)->globalPosition().z() - Center.z()));
687  vec[index] = vec_temp;
688  index++;
689  }
690  isClockwise = 0;
691  for (unsigned int index = 0; index < (NumberofHitsinSeed - 1); index++) {
692  // Check phi direction, all sub-dphi direction should be the same
693  if ((vec[index + 1].phi() - vec[index].phi()) > 0)
694  isClockwise--;
695  else
696  isClockwise++;
697  cout << "Current isClockwise is : " << isClockwise << endl;
698  }
699  cout << "Check isClockwise is : " << isClockwise << endl;
700  if ((unsigned int)abs(isClockwise) != (NumberofHitsinSeed - 1)) {
701  cout << "Pattern find error in Phi direction" << endl;
702  isGoodPattern = 0;
703  isClockwise = 0;
704  } else
705  isClockwise /= abs(isClockwise);
706  delete[] vec;
707 
708  // Get meanPt and meanSpt
709  double deltaRwithBz = fabs(deltaBz * meanRadius / meanBz);
710  cout << "deltaR with Bz is " << deltaRwithBz << endl;
711 
712  if (isClockwise == 0)
713  meanPt = upper_limit_pt;
714  else
715  meanPt = 0.01 * meanRadius * meanBz * 0.3 * isClockwise;
716  if (fabs(meanPt) > upper_limit_pt)
717  meanPt = upper_limit_pt * meanPt / fabs(meanPt);
718  cout << " meanRadius is " << meanRadius << endl;
719  cout << " meanPt is " << meanPt << endl;
720 
721  double deltaR = 0.;
722  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
723  deltaR += (getDistance(*iter, Center) - meanRadius) * (getDistance(*iter, Center) - meanRadius);
724  }
725  deltaR = deltaR / NumberofHitsinSeed;
726  deltaR = sqrt(deltaR);
727  //meanSpt = 0.01 * deltaR * meanBz * 0.3;
728  meanSpt = deltaR;
729  cout << "DeltaR is " << deltaR << endl;
730  if (deltaR > deltaRThreshold) {
731  cout << "Pattern find error: deltaR over threshold" << endl;
732  isGoodPattern = 0;
733  }
734  } else {
735  // Just set pattern to be straight
736  isClockwise = 0;
737  meanPt = upper_limit_pt;
738  // Set the straight pattern with lowest priority among good pattern
739  meanSpt = deltaRThreshold;
740  }
741  cout << "III--> Seed Pt : " << meanPt << endl;
742  cout << "III--> Pattern is: " << isGoodPattern << endl;
743 
744  // Set the pattern estimation signal
745  isPatternChecked = true;
746 }
747 
748 void RPCSeedPattern::checkSegmentAlgorithmSpecial(const edm::EventSetup& eSetup) {
749  if (isPatternChecked == true)
750  return;
751 
752  if (!checkSegment()) {
753  isPatternChecked = true;
754  isGoodPattern = -1;
755  return;
756  }
757 
758  // Set isGoodPattern to default true and check the failure
759  isGoodPattern = 1;
760 
761  // Print the recHit's position
762  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
763  cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
764 
765  // Check the Z direction
766  if (fabs((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError) {
767  if (((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
768  isParralZ = 1;
769  else
770  isParralZ = -1;
771  } else
772  isParralZ = 0;
773 
774  cout << " Check isParralZ is :" << isParralZ << endl;
775  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
776  if (isParralZ == 0) {
777  if (fabs((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) > ZError) {
778  cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
779  isGoodPattern = 0;
780  }
781  } else {
782  if ((int)(((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) / ZError) * isParralZ < 0) {
783  cout << "Pattern find error in Z direction: wrong Z direction" << endl;
784  isGoodPattern = 0;
785  }
786  }
787  }
788 
789  // Check the pattern
790  if (isStraight2 == true) {
791  // Set pattern to be straight
792  isClockwise = 0;
793  meanPt = upper_limit_pt;
794  // Set the straight pattern with lowest priority among good pattern
795  meanSpt = deltaRThreshold;
796 
797  // Extrapolate to other recHits and check deltaR
798  GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
799  GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
800  GlobalVector startMomentum = startSegment * (meanPt / startSegment.perp());
801  unsigned int index = 0;
802  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
803  if (index < 4) {
804  index++;
805  continue;
806  }
807  double tracklength = 0;
808  cout << "Now checking recHit " << index << endl;
809  double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, eSetup);
810  cout << "Final distance is " << Distance << endl;
811  if (Distance > MaxRSD) {
812  cout << "Pattern find error in distance for other recHits: " << Distance << endl;
813  isGoodPattern = 0;
814  }
815  index++;
816  }
817  } else {
818  // Get clockwise direction
819  GlobalVector vec[2];
820  vec[0] = GlobalVector(
821  (entryPosition.x() - Center2.x()), (entryPosition.y() - Center2.y()), (entryPosition.z() - Center2.z()));
822  vec[1] = GlobalVector(
823  (leavePosition.x() - Center2.x()), (leavePosition.y() - Center2.y()), (leavePosition.z() - Center2.z()));
824  isClockwise = 0;
825  if ((vec[1].phi() - vec[0].phi()).value() > 0)
826  isClockwise = -1;
827  else
828  isClockwise = 1;
829 
830  cout << "Check isClockwise is : " << isClockwise << endl;
831 
832  // Get meanPt
833  meanPt = 0.01 * meanRadius2 * meanMagneticField2.z() * 0.3 * isClockwise;
834  //meanPt = 0.01 * meanRadius2[0] * (-3.8) * 0.3 * isClockwise;
835  cout << " meanRadius is " << meanRadius2 << ", with meanBz " << meanMagneticField2.z() << endl;
836  cout << " meanPt is " << meanPt << endl;
837  if (fabs(meanPt) > upper_limit_pt)
838  meanPt = upper_limit_pt * meanPt / fabs(meanPt);
839 
840  // Check the initial 3 segments
841  cout << "entryPosition: " << entryPosition << endl;
842  cout << "leavePosition: " << leavePosition << endl;
843  cout << "Center2 is : " << Center2 << endl;
844  double R1 = vec[0].perp();
845  double R2 = vec[1].perp();
846  double deltaR = sqrt(((R1 - meanRadius2) * (R1 - meanRadius2) + (R2 - meanRadius2) * (R2 - meanRadius2)) / 2);
847  meanSpt = deltaR;
848  cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
849  cout << "Delta R for the initial 3 segments is " << deltaR << endl;
850  if (deltaR > deltaRThreshold) {
851  cout << "Pattern find error in delta R for the initial 3 segments" << endl;
852  isGoodPattern = 0;
853  }
854 
855  // Extrapolate to other recHits and check deltaR
856  GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
857  GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
858  GlobalVector startMomentum = startSegment * (fabs(meanPt) / startSegment.perp());
859  unsigned int index = 0;
860  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
861  if (index < 4) {
862  index++;
863  continue;
864  }
865  double tracklength = 0;
866  cout << "Now checking recHit " << index << endl;
867  double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, eSetup);
868  cout << "Final distance is " << Distance << endl;
869  if (Distance > MaxRSD) {
870  cout << "Pattern find error in distance for other recHits: " << Distance << endl;
871  isGoodPattern = 0;
872  }
873  index++;
874  }
875  }
876 
877  cout << "Checking finish, isGoodPattern now is " << isGoodPattern << endl;
878  // Set the pattern estimation signal
879  isPatternChecked = true;
880 }
881 
882 double RPCSeedPattern::extropolateStep(const GlobalPoint& startPosition,
883  const GlobalVector& startMomentum,
884  ConstMuonRecHitContainer::const_iterator iter,
885  const int ClockwiseDirection,
886  double& tracklength,
887  const edm::EventSetup& eSetup) {
888  // Get magnetic field
889  edm::ESHandle<MagneticField> Field;
890  eSetup.get<IdealMagneticFieldRecord>().get(Field);
891 
892  cout << "Extrapolating the track to check the pattern" << endl;
893  tracklength = 0;
894  // Get the iter recHit's detector geometry
895  DetId hitDet = (*iter)->hit()->geographicalId();
896  RPCDetId RPCId = RPCDetId(hitDet.rawId());
897  //const RPCChamber* hitRPC = dynamic_cast<const RPCChamber*>(hitDet);
898  edm::ESHandle<RPCGeometry> pRPCGeom;
899  eSetup.get<MuonGeometryRecord>().get(pRPCGeom);
900  const RPCGeometry* rpcGeometry = (const RPCGeometry*)&*pRPCGeom;
901 
902  const BoundPlane RPCSurface = rpcGeometry->chamber(RPCId)->surface();
903  double startSide = RPCSurface.localZ(startPosition);
904  cout << "Start side : " << startSide;
905 
906  GlobalPoint currentPosition = startPosition;
907  GlobalVector currentMomentum = startMomentum;
908  GlobalVector ZDirection(0, 0, 1);
909 
910  // Use the perp other than mag, since initial segment might have small value while final recHit have large difference value at Z direction
911  double currentDistance = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
912  cout << "Start current position is : " << currentPosition << endl;
913  cout << "Start current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
914  cout << "Start current distance is " << currentDistance << endl;
915  cout << "Start current radius is " << currentPosition.perp() << endl;
916  cout << "Destination radius is " << (*iter)->globalPosition().perp() << endl;
917 
918  // Judge roughly if the stepping cross the Det surface of the recHit
919  //while((currentPosition.perp() < ((*iter)->globalPosition().perp())))
920  double currentDistance_next = currentDistance;
921  do {
922  currentDistance = currentDistance_next;
923  if (ClockwiseDirection == 0) {
924  currentPosition += currentMomentum.unit() * stepLength;
925  } else {
926  double Bz = Field->inTesla(currentPosition).z();
927  double Radius = currentMomentum.perp() / fabs(Bz * 0.01 * 0.3);
928  double deltaPhi = (stepLength * currentMomentum.perp() / currentMomentum.mag()) / Radius;
929 
930  // Get the center for current step
931  GlobalVector currentPositiontoCenter = currentMomentum.unit().cross(ZDirection);
932  currentPositiontoCenter *= Radius;
933  // correction of ClockwiseDirection correction
934  currentPositiontoCenter *= ClockwiseDirection;
935  // continue to get the center for current step
936  GlobalPoint currentCenter = currentPosition;
937  currentCenter += currentPositiontoCenter;
938 
939  // Get the next step position
940  GlobalVector CentertocurrentPosition = (GlobalVector)(currentPosition - currentCenter);
941  double Phi = CentertocurrentPosition.phi().value();
942  Phi += deltaPhi * (-1) * ClockwiseDirection;
943  double deltaZ = stepLength * currentMomentum.z() / currentMomentum.mag();
944  GlobalVector CentertonewPosition(GlobalVector::Cylindrical(CentertocurrentPosition.perp(), Phi, deltaZ));
945  double PtPhi = currentMomentum.phi().value();
946  PtPhi += deltaPhi * (-1) * ClockwiseDirection;
947  currentMomentum = GlobalVector(GlobalVector::Cylindrical(currentMomentum.perp(), PtPhi, currentMomentum.z()));
948  currentPosition = currentCenter + CentertonewPosition;
949  }
950 
951  // count the total step length
952  tracklength += stepLength * currentMomentum.perp() / currentMomentum.mag();
953 
954  // Get the next step distance
955  double currentSide = RPCSurface.localZ(currentPosition);
956  cout << "Stepping current side : " << currentSide << endl;
957  cout << "Stepping current position is: " << currentPosition << endl;
958  cout << "Stepping current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
959  currentDistance_next = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
960  cout << "Stepping current distance is " << currentDistance << endl;
961  cout << "Stepping current radius is " << currentPosition.perp() << endl;
962  } while (currentDistance_next < currentDistance);
963 
964  return currentDistance;
965 }
966 
967 RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createFakeSeed(int& isGoodSeed, const edm::EventSetup& eSetup) {
968  // Create a fake seed and return
969  cout << "Now create a fake seed" << endl;
970  isPatternChecked = true;
971  isGoodPattern = -1;
972  isStraight = true;
973  meanPt = upper_limit_pt;
974  meanSpt = 0;
975  Charge = 0;
976  isClockwise = 0;
977  isParralZ = 0;
978  meanRadius = -1;
979  //return createSeed(isGoodSeed, eSetup);
980 
981  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
982  const ConstMuonRecHitPointer best = BestRefRecHit();
983 
984  Momentum = GlobalVector(0, 0, 0);
985  LocalPoint segPos = best->localPosition();
986  LocalVector segDirFromPos = best->det()->toLocal(Momentum);
987  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
988 
989  //AlgebraicVector t(4);
990  AlgebraicSymMatrix mat(5, 0);
991  mat = best->parametersError().similarityT(best->projectionMatrix());
992  mat[0][0] = meanSpt;
993  LocalTrajectoryError error(asSMatrix<5>(mat));
994 
995  edm::ESHandle<MagneticField> Field;
996  eSetup.get<IdealMagneticFieldRecord>().get(Field);
997 
998  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &*Field);
999 
1000  DetId id = best->geographicalId();
1001 
1002  PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
1003 
1004  edm::OwnVector<TrackingRecHit> container;
1005  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
1006  container.push_back((*iter)->hit()->clone());
1007 
1008  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
1009  weightedTrajectorySeed theweightedSeed;
1010  theweightedSeed.first = theSeed;
1011  theweightedSeed.second = meanSpt;
1012  isGoodSeed = isGoodPattern;
1013 
1014  return theweightedSeed;
1015 }
1016 
1017 RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createSeed(int& isGoodSeed, const edm::EventSetup& eSetup) {
1018  if (isPatternChecked == false || isGoodPattern == -1) {
1019  cout << "Pattern is not yet checked! Create a fake seed instead!" << endl;
1020  return createFakeSeed(isGoodSeed, eSetup);
1021  }
1022 
1023  edm::ESHandle<MagneticField> Field;
1024  eSetup.get<IdealMagneticFieldRecord>().get(Field);
1025 
1026  MuonPatternRecoDumper debug;
1027 
1028  //double theMinMomentum = 3.0;
1029  //if(fabs(meanPt) < lower_limit_pt)
1030  //meanPt = lower_limit_pt * meanPt / fabs(meanPt);
1031 
1032  // For pattern we use is Clockwise other than isStraight to estimate charge
1033  if (isClockwise == 0)
1034  Charge = 0;
1035  else
1036  Charge = (int)(meanPt / fabs(meanPt));
1037 
1038  // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
1039  const ConstMuonRecHitPointer best = BestRefRecHit();
1040  const ConstMuonRecHitPointer first = FirstRecHit();
1041 
1042  if (isClockwise != 0) {
1043  if (AlgorithmType != 3) {
1044  // Get the momentum on reference recHit
1045  GlobalVector vecRef1((first->globalPosition().x() - Center.x()),
1046  (first->globalPosition().y() - Center.y()),
1047  (first->globalPosition().z() - Center.z()));
1048  GlobalVector vecRef2((best->globalPosition().x() - Center.x()),
1049  (best->globalPosition().y() - Center.y()),
1050  (best->globalPosition().z() - Center.z()));
1051 
1052  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1053  double deltaS = meanRadius * fabs(deltaPhi);
1054  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1055 
1056  GlobalVector vecZ(0, 0, 1);
1057  GlobalVector vecPt = (vecRef2.unit()).cross(vecZ);
1058  if (isClockwise == -1)
1059  vecPt *= -1;
1060  vecPt *= deltaS;
1061  Momentum = GlobalVector(0, 0, deltaZ);
1062  Momentum += vecPt;
1063  Momentum *= fabs(meanPt / deltaS);
1064  } else {
1065  double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1066  Momentum = GlobalVector(GlobalVector::Cylindrical(S, lastPhi, deltaZ));
1067  Momentum *= fabs(meanPt / S);
1068  }
1069  } else {
1070  Momentum = best->globalPosition() - first->globalPosition();
1071  double deltaS = Momentum.perp();
1072  Momentum *= fabs(meanPt / deltaS);
1073  }
1074  LocalPoint segPos = best->localPosition();
1075  LocalVector segDirFromPos = best->det()->toLocal(Momentum);
1076  LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
1077 
1078  LocalTrajectoryError error = getSpecialAlgorithmErrorMatrix(first, best);
1079 
1080  TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &*Field);
1081  cout << "Trajectory State on Surface before the extrapolation" << endl;
1082  cout << debug.dumpTSOS(tsos);
1083  DetId id = best->geographicalId();
1084  cout << "The RecSegment relies on: " << endl;
1085  cout << debug.dumpMuonId(id);
1086  cout << debug.dumpTSOS(tsos);
1087 
1088  PTrajectoryStateOnDet const& seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
1089 
1090  edm::OwnVector<TrackingRecHit> container;
1091  for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
1092  // This casting withou clone will cause memory overflow when used in push_back
1093  // Since container's deconstructor functiion free the pointer menber!
1094  //TrackingRecHit* pt = dynamic_cast<TrackingRecHit*>(&*(*iter));
1095  //cout << "Push recHit type " << pt->getType() << endl;
1096  container.push_back((*iter)->hit()->clone());
1097  }
1098 
1099  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
1100  weightedTrajectorySeed theweightedSeed;
1101  theweightedSeed.first = theSeed;
1102  theweightedSeed.second = meanSpt;
1103  isGoodSeed = isGoodPattern;
1104 
1105  return theweightedSeed;
1106 }
1107 
1108 LocalTrajectoryError RPCSeedPattern::getSpecialAlgorithmErrorMatrix(const ConstMuonRecHitPointer& first,
1109  const ConstMuonRecHitPointer& best) {
1110  LocalTrajectoryError Error;
1111  double dXdZ = 0;
1112  double dYdZ = 0;
1113  double dP = 0;
1114  AlgebraicSymMatrix mat(5, 0);
1115  mat = best->parametersError().similarityT(best->projectionMatrix());
1116  if (AlgorithmType != 3) {
1117  GlobalVector vecRef1((first->globalPosition().x() - Center.x()),
1118  (first->globalPosition().y() - Center.y()),
1119  (first->globalPosition().z() - Center.z()));
1120  GlobalVector vecRef2((best->globalPosition().x() - Center.x()),
1121  (best->globalPosition().y() - Center.y()),
1122  (best->globalPosition().z() - Center.z()));
1123  double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1124  double L = meanRadius * fabs(deltaPhi);
1125  double N = nrhit();
1126  double A_N = 180 * N * N * N / ((N - 1) * (N + 1) * (N + 2) * (N + 3));
1127  double sigma_x = sqrt(mat[3][3]);
1128  double betaovergame = Momentum.mag() / 0.1066;
1129  double beta = sqrt((betaovergame * betaovergame) / (1 + betaovergame * betaovergame));
1130  double dPt = meanPt * (0.0136 * sqrt(1 / 100) * sqrt(4 * A_N / N) / (beta * 0.3 * meanBz * L) +
1131  sigma_x * meanPt * sqrt(4 * A_N) / (0.3 * meanBz * L * L));
1132  double dP = dPt * Momentum.mag() / meanPt;
1133  mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1134  mat[1][1] = dXdZ * dXdZ;
1135  mat[2][2] = dYdZ * dYdZ;
1136  Error = LocalTrajectoryError(asSMatrix<5>(mat));
1137  } else {
1138  AlgebraicSymMatrix mat0(5, 0);
1139  mat0 = (SegmentRB[0].first)->parametersError().similarityT((SegmentRB[0].first)->projectionMatrix());
1140  double dX0 = sqrt(mat0[3][3]);
1141  double dY0 = sqrt(mat0[4][4]);
1142  AlgebraicSymMatrix mat1(5, 0);
1143  mat1 = (SegmentRB[0].second)->parametersError().similarityT((SegmentRB[0].second)->projectionMatrix());
1144  double dX1 = sqrt(mat1[3][3]);
1145  double dY1 = sqrt(mat1[4][4]);
1146  AlgebraicSymMatrix mat2(5, 0);
1147  mat2 = (SegmentRB[1].first)->parametersError().similarityT((SegmentRB[1].first)->projectionMatrix());
1148  double dX2 = sqrt(mat2[3][3]);
1149  double dY2 = sqrt(mat2[4][4]);
1150  AlgebraicSymMatrix mat3(5, 0);
1151  mat3 = (SegmentRB[1].second)->parametersError().similarityT((SegmentRB[1].second)->projectionMatrix());
1152  double dX3 = sqrt(mat3[3][3]);
1153  double dY3 = sqrt(mat3[4][4]);
1154  cout << "Local error for 4 recHits are: " << dX0 << ", " << dY0 << ", " << dX1 << ", " << dY1 << ", " << dX2 << ", "
1155  << dY2 << ", " << dX3 << ", " << dY3 << endl;
1156  const GeomDetUnit* refRPC1 = (SegmentRB[0].second)->detUnit();
1157  LocalPoint recHit0 = refRPC1->toLocal((SegmentRB[0].first)->globalPosition());
1158  LocalPoint recHit1 = refRPC1->toLocal((SegmentRB[0].second)->globalPosition());
1159  LocalVector localSegment00 = (LocalVector)(recHit1 - recHit0);
1160  LocalVector localSegment01 = LocalVector(localSegment00.x() + dX0 + dX1, localSegment00.y(), localSegment00.z());
1161  LocalVector localSegment02 = LocalVector(localSegment00.x() - dX0 - dX1, localSegment00.y(), localSegment00.z());
1162  GlobalVector globalSegment00 = refRPC1->toGlobal(localSegment00);
1163  GlobalVector globalSegment01 = refRPC1->toGlobal(localSegment01);
1164  GlobalVector globalSegment02 = refRPC1->toGlobal(localSegment02);
1165 
1166  const GeomDetUnit* refRPC2 = (SegmentRB[1].first)->detUnit();
1167  LocalPoint recHit2 = refRPC2->toLocal((SegmentRB[1].first)->globalPosition());
1168  LocalPoint recHit3 = refRPC2->toLocal((SegmentRB[1].second)->globalPosition());
1169  LocalVector localSegment10 = (LocalVector)(recHit3 - recHit2);
1170  LocalVector localSegment11 = LocalVector(localSegment10.x() + dX2 + dX3, localSegment10.y(), localSegment10.z());
1171  LocalVector localSegment12 = LocalVector(localSegment10.x() - dX2 - dX3, localSegment10.y(), localSegment10.z());
1172  GlobalVector globalSegment10 = refRPC2->toGlobal(localSegment10);
1173  GlobalVector globalSegment11 = refRPC2->toGlobal(localSegment11);
1174  GlobalVector globalSegment12 = refRPC2->toGlobal(localSegment12);
1175 
1176  if (isClockwise != 0) {
1177  GlobalVector vec[2];
1178  vec[0] = GlobalVector(
1179  (entryPosition.x() - Center2.x()), (entryPosition.y() - Center2.y()), (entryPosition.z() - Center2.z()));
1180  vec[1] = GlobalVector(
1181  (leavePosition.x() - Center2.x()), (leavePosition.y() - Center2.y()), (leavePosition.z() - Center2.z()));
1182  double halfPhiCenter = fabs((vec[1].phi() - vec[0].phi()).value()) / 2;
1183  // dPhi0 shoule be the same clockwise direction, while dPhi1 should be opposite clockwise direction, w.r.t to the track clockwise
1184  double dPhi0 = (((globalSegment00.phi() - globalSegment01.phi()).value() * isClockwise) > 0)
1185  ? fabs((globalSegment00.phi() - globalSegment01.phi()).value())
1186  : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1187  double dPhi1 = (((globalSegment10.phi() - globalSegment11.phi()).value() * isClockwise) < 0)
1188  ? fabs((globalSegment10.phi() - globalSegment11.phi()).value())
1189  : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1190  // For the deltaR should be kept small, we assume the delta Phi0/Phi1 should be in a same limit value
1191  double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1192  cout << "DPhi for new Segment is about " << dPhi << endl;
1193  // Check the variance of halfPhiCenter
1194  double newhalfPhiCenter = ((halfPhiCenter - dPhi) > 0 ? (halfPhiCenter - dPhi) : 0);
1195  if (newhalfPhiCenter != 0) {
1196  double newmeanPt = meanPt * halfPhiCenter / newhalfPhiCenter;
1197  if (fabs(newmeanPt) > upper_limit_pt)
1198  newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1199  cout << "The error is inside range. Max meanPt could be " << newmeanPt << endl;
1200  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1201  } else {
1202  double newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1203  cout << "The error is outside range. Max meanPt could be " << newmeanPt << endl;
1204  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1205  }
1206  } else {
1207  double dPhi0 = (fabs((globalSegment00.phi() - globalSegment01.phi()).value()) <=
1208  fabs((globalSegment00.phi() - globalSegment02.phi()).value()))
1209  ? fabs((globalSegment00.phi() - globalSegment01.phi()).value())
1210  : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1211  double dPhi1 = (fabs((globalSegment10.phi() - globalSegment11.phi()).value()) <=
1212  fabs((globalSegment10.phi() - globalSegment12.phi()).value()))
1213  ? fabs((globalSegment10.phi() - globalSegment11.phi()).value())
1214  : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1215  double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1216  GlobalVector middleSegment = leavePosition - entryPosition;
1217  double halfDistance = middleSegment.perp() / 2;
1218  double newmeanPt = halfDistance / dPhi;
1219  cout << "The error is for straight. Max meanPt could be " << newmeanPt << endl;
1220  dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1221  }
1222 
1223  double dXdZ1 = globalSegment11.x() / globalSegment11.z() - globalSegment10.x() / globalSegment10.z();
1224  double dXdZ2 = globalSegment12.x() / globalSegment12.z() - globalSegment10.x() / globalSegment10.z();
1225  dXdZ = (fabs(dXdZ1) >= fabs(dXdZ2)) ? dXdZ1 : dXdZ2;
1226 
1227  LocalVector localSegment13 = LocalVector(localSegment10.x(), localSegment10.y() + dY2 + dY3, localSegment10.z());
1228  LocalVector localSegment14 = LocalVector(localSegment10.x(), localSegment10.y() - dY2 - dY3, localSegment10.z());
1229  GlobalVector globalSegment13 = refRPC2->toGlobal(localSegment13);
1230  GlobalVector globalSegment14 = refRPC2->toGlobal(localSegment14);
1231  double dYdZ1 = globalSegment13.y() / globalSegment13.z() - globalSegment10.y() / globalSegment10.z();
1232  double dYdZ2 = globalSegment14.y() / globalSegment14.z() - globalSegment10.y() / globalSegment10.z();
1233  dYdZ = (fabs(dYdZ1) >= fabs(dYdZ2)) ? dYdZ1 : dYdZ2;
1234 
1235  mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1236  mat[1][1] = dXdZ * dXdZ;
1237  mat[2][2] = dYdZ * dYdZ;
1238  Error = LocalTrajectoryError(asSMatrix<5>(mat));
1239  }
1240  return Error;
1241 }
Vector3DBase
Definition: Vector3DBase.h:8
RPCSeedPattern::getSpecialAlgorithmErrorMatrix
LocalTrajectoryError getSpecialAlgorithmErrorMatrix(const ConstMuonRecHitPointer &first, const ConstMuonRecHitPointer &best)
Definition: RPCSeedPattern.cc:1108
dttmaxenums::L
Definition: DTTMax.h:29
MagneticField::inTesla
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
mps_fire.i
i
Definition: mps_fire.py:428
RPCSeedPattern::ThreePointsAlgorithm
void ThreePointsAlgorithm()
Definition: RPCSeedPattern.cc:107
MuonTransientTrackingRecHit::ConstMuonRecHitPointer
std::shared_ptr< MuonTransientTrackingRecHit const > ConstMuonRecHitPointer
Definition: MuonTransientTrackingRecHit.h:22
MessageLogger.h
MuonPatternRecoDumper.h
RPCDetId::station
int station() const
Definition: RPCDetId.h:78
RPCDetId::region
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:53
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
GeomDet
Definition: GeomDet.h:27
ESHandle.h
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
X
#define X(str)
Definition: MuonsGrabber.cc:38
RPCSeedPattern::createSeed
weightedTrajectorySeed createSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
Definition: RPCSeedPattern.cc:1017
RPCSeedPattern::extropolateStep
double extropolateStep(const GlobalPoint &startPosition, const GlobalVector &startMomentum, ConstMuonRecHitContainer::const_iterator iter, const int ClockwiseDirection, double &tracklength, const edm::EventSetup &eSetup)
Definition: RPCSeedPattern.cc:882
edm
HLT enums.
Definition: AlignableModifier.h:19
zMuMuMuonUserData.beta
beta
Definition: zMuMuMuonUserData.py:10
trajectoryStateTransform::persistentState
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
Definition: TrajectoryStateTransform.cc:14
RPCDetId
Definition: RPCDetId.h:16
RPCChamber::id
RPCDetId id() const
Return the RPCChamberId of this chamber.
Definition: RPCChamber.cc:26
gather_cfg.cout
cout
Definition: gather_cfg.py:144
RPCChamber.h
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:222
JetPartonCorrections_cff.Radius
Radius
Definition: JetPartonCorrections_cff.py:11
RPCSeedPattern::weightedTrajectorySeed
std::pair< TrajectorySeed, double > weightedTrajectorySeed
Definition: RPCSeedPattern.h:36
GlobalVector
Global3DVector GlobalVector
Definition: GlobalVector.h:10
HLT_FULL_cff.dPhi
dPhi
Definition: HLT_FULL_cff.py:13765
RPCSeedPattern::checkStraightwithSegment
bool checkStraightwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2, double MinDeltaPhi) const
Definition: RPCSeedPattern.cc:528
relativeConstraints.error
error
Definition: relativeConstraints.py:53
dqmdumpme.first
first
Definition: dqmdumpme.py:55
RPCChamber
Definition: RPCChamber.h:19
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
LocalTrajectoryParameters
Definition: LocalTrajectoryParameters.h:25
IdealMagneticFieldRecord
Definition: IdealMagneticFieldRecord.h:11
DetId
Definition: DetId.h:17
GeomDet::surface
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
TrajectoryStateOnSurface
Definition: TrajectoryStateOnSurface.h:16
MuonME0Digis_cfi.sigma_x
sigma_x
Definition: MuonME0Digis_cfi.py:9
debug
#define debug
Definition: HDRShower.cc:19
Vector3DBase::unit
Vector3DBase unit() const
Definition: Vector3DBase.h:54
edm::EventSetup::get
T get() const
Definition: EventSetup.h:80
upper_limit_pt
#define upper_limit_pt
Definition: RPCSeedPattern.h:21
PTrajectoryStateOnDet.h
RPCSeedPattern::createFakeSeed
weightedTrajectorySeed createFakeSeed(int &isGoodSeed, const edm::EventSetup &eSetup)
Definition: RPCSeedPattern.cc:967
MuonPatternRecoDumper
Definition: MuonPatternRecoDumper.h:18
SiPixelRawToDigiRegional_cfi.deltaPhi
deltaPhi
Definition: SiPixelRawToDigiRegional_cfi.py:9
PythiaFilterGammaJetWithBg_cfi.MinDeltaPhi
MinDeltaPhi
Definition: PythiaFilterGammaJetWithBg_cfi.py:8
RPCSeedPattern::computePtwithSegment
GlobalVector computePtwithSegment(const RPCSegment &Segment1, const RPCSegment &Segment2) const
Definition: RPCSeedPattern.cc:549
RPCSeedPattern::FirstRecHit
ConstMuonRecHitPointer FirstRecHit() const
Definition: RPCSeedPattern.cc:468
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
Geom::Cylindrical2Cartesian
Definition: CoordinateSets.h:34
RPCSeedPattern::RPCSeedPattern
RPCSeedPattern()
Definition: RPCSeedPattern.cc:32
IdealMagneticFieldRecord.h
edm::ESHandle< MagneticField >
RPCSeedPattern::computePtwithThreerecHits
GlobalVector computePtwithThreerecHits(double &pt, double &pt_err, ConstMuonRecHitPointer(&precHit)[3]) const
Definition: RPCSeedPattern.cc:503
N
#define N
Definition: blowfish.cc:9
submitPVResolutionJobs.count
count
Definition: submitPVResolutionJobs.py:352
VtxSmearedParameters_cfi.Phi
Phi
Definition: VtxSmearedParameters_cfi.py:112
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
dqmdumpme.k
k
Definition: dqmdumpme.py:60
Point3DBase< float, GlobalTag >
GeomDet::toLocal
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:58
PbPb_ZMuSkimMuonDPG_cff.deltaR
deltaR
Definition: PbPb_ZMuSkimMuonDPG_cff.py:63
RPCDetId.h
LocalTrajectoryError
Definition: LocalTrajectoryError.h:20
runTauDisplay.gp
gp
Definition: runTauDisplay.py:431
cross
Basic3DVector cross(const Basic3DVector &v) const
Vector product, or "cross" product, with a vector of same type.
Definition: Basic3DVectorLD.h:225
RPCSeedPattern::BestRefRecHit
ConstMuonRecHitPointer BestRefRecHit() const
Definition: RPCSeedPattern.cc:470
GeomDet::toGlobal
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
edm::ParameterSet
Definition: ParameterSet.h:47
RPCSeedPattern::~RPCSeedPattern
~RPCSeedPattern()
Definition: RPCSeedPattern.cc:38
RPCSeedPattern::RPCSegment
std::pair< ConstMuonRecHitPointer, ConstMuonRecHitPointer > RPCSegment
Definition: RPCSeedPattern.h:35
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
leef::Error
edm::ErrorSummaryEntry Error
Definition: LogErrorEventFilter.cc:29
LocalVector
Local3DVector LocalVector
Definition: LocalVector.h:12
createfilelist.int
int
Definition: createfilelist.py:10
RPCSeedPattern::seed
weightedTrajectorySeed seed(const edm::EventSetup &eSetup, int &isGoodSeed)
Definition: RPCSeedPattern.cc:53
RPCSeedPattern::checkSegmentAlgorithmSpecial
void checkSegmentAlgorithmSpecial(const edm::EventSetup &eSetup)
Definition: RPCSeedPattern.cc:748
Vector3DBase::cross
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:110
g4SimHits_cfi.Region
Region
Definition: g4SimHits_cfi.py:152
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:15
PVValHelper::dy
Definition: PVValidationHelpers.h:49
MagneticField.h
edm::EventSetup
Definition: EventSetup.h:57
get
#define get
RPCSeedPattern::ConstMuonRecHitPointer
MuonTransientTrackingRecHit::ConstMuonRecHitPointer ConstMuonRecHitPointer
Definition: RPCSeedPattern.h:30
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
align::Detector
Definition: StructureType.h:86
RPCGeometry::chamber
const RPCChamber * chamber(RPCDetId id) const
Definition: RPCGeometry.cc:46
GeomDet.h
RPCSeedPattern::SegmentAlgorithmSpecial
void SegmentAlgorithmSpecial(const edm::EventSetup &eSetup)
Definition: RPCSeedPattern.cc:277
MaterialEffects_cfi.A
A
Definition: MaterialEffects_cfi.py:11
std
Definition: JetResolutionObject.h:76
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
PVValHelper::dz
Definition: PVValidationHelpers.h:50
relativeConstraints.value
value
Definition: relativeConstraints.py:53
RPCSeedPattern::checkSegment
bool checkSegment() const
Definition: RPCSeedPattern.cc:440
BoundPlane
muon::AlgorithmType
AlgorithmType
Definition: MuonSelectors.h:154
RPCSeedPattern::checkSimplePattern
void checkSimplePattern(const edm::EventSetup &eSetup)
Definition: RPCSeedPattern.cc:607
PixelTestBeamValidation_cfi.Charge
Charge
Definition: PixelTestBeamValidation_cfi.py:91
DetLayer.h
TrajectorySeed
Definition: TrajectorySeed.h:18
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
TrajectoryStateTransform.h
BeamSpotPI::Y
Definition: BeamSpotPayloadInspectorHelper.h:31
edm::OwnVector::push_back
void push_back(D *&d)
Definition: OwnVector.h:326
S
Definition: CSCDBL1TPParametersExtended.h:16
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
RPCSeedPattern::getDistance
double getDistance(const ConstMuonRecHitPointer &precHit, const GlobalVector &Center) const
Definition: RPCSeedPattern.cc:485
RPCSeedPattern::SegmentAlgorithm
void SegmentAlgorithm()
Definition: RPCSeedPattern.cc:226
RPCGeometry
Definition: RPCGeometry.h:20
RPCSeedPattern::computePtWithThreerecHits
GlobalVector computePtWithThreerecHits(double &pt, double &pt_err, double(&x)[3], double(&y)[3]) const
Definition: RPCSeedPattern.cc:588
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
PTrajectoryStateOnDet
Definition: PTrajectoryStateOnDet.h:10
OwnVector.h
MuonGeometryRecord.h
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
RPCSeedPattern.h
volumeBasedMagneticField_160812_cfi.BValue
BValue
Definition: volumeBasedMagneticField_160812_cfi.py:28
A
MuonGeometryRecord
Definition: MuonGeometryRecord.h:34
RPCGeometry.h
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
PVValHelper::dx
Definition: PVValidationHelpers.h:48
RPCSeedPattern::MiddlePointsAlgorithm
void MiddlePointsAlgorithm()
Definition: RPCSeedPattern.cc:169
alongMomentum
Definition: PropagationDirection.h:4
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
RPCSeedPattern::checkStraightwithThreerecHits
bool checkStraightwithThreerecHits(ConstMuonRecHitPointer(&precHit)[3], double MinDeltaPhi) const
Definition: RPCSeedPattern.cc:490
edm::OwnVector< TrackingRecHit >
RPCSeedPattern::configure
void configure(const edm::ParameterSet &iConfig)
Definition: RPCSeedPattern.cc:40
Generated for CMSSW Reference Manual by doxygen 1.8.16