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