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