CMS 3D CMS Logo

L1MuonPixelTrackFitter.cc
Go to the documentation of this file.
2 
4 
12 
13 #include <math.h>
14 
15 template <class T> T sqr( T t) {return t*t;}
16 
18  : theConfig(cfg),
19  invPtErrorScale{theConfig.getParameter<double>("invPtErrorScale")},
20  phiErrorScale{theConfig.getParameter<double>("phiErrorScale")},
21  cotThetaErrorScale{theConfig.getParameter<double>("cotThetaErrorScale")},
22  tipErrorScale{theConfig.getParameter<double>("tipErrorScale")},
23  zipErrorScale{theConfig.getParameter<double>("zipErrorScale")}
24 { }
25 
27 {
28  thePhiL1 = muon.phiValue()+0.021817;
29  theEtaL1 = muon.etaValue();
30  theChargeL1 = muon.charge();
31 }
32 
34 {
35  theHit1 = hits[0]->globalPosition();
36  theHit1 = hits[1]->globalPosition();
37 }
38 
40  const std::vector<const TrackingRecHit *>& hits, const TrackingRegion& region) const
41 {
43  es.get<IdealMagneticFieldRecord>().get(fieldESH);
44 
45  double phi_vtx_fromHits = (theHit2-theHit1).phi();
46 
47  double invPt = valInversePt(phi_vtx_fromHits, thePhiL1, theEtaL1);
48  double invPtErr = errInversePt(invPt, theEtaL1);
49 
50  int charge = (invPt > 0.) ? 1 : -1;
51  double curvature = PixelRecoUtilities::curvature(invPt, es);
52 
53  Circle circle(theHit1, theHit2, curvature);
54 
55  double valPhi = this->valPhi(circle, charge);
56  double valTip = this->valTip(circle, curvature);
57  double valZip = this->valZip(curvature, theHit1,theHit2);
59 
60 // if ( (fabs(invPt)-0.1)/invPtErr > 3.) return 0;
61 
62  PixelTrackBuilder builder;
63  double valPt = (fabs(invPt) > 1.e-4) ? 1./fabs(invPt) : 1.e4;
64  double errPt = invPtErrorScale*invPtErr * valPt*valPt;
65  double eta = asinh(valCotTheta);
66  double errPhi = this->errPhi(invPt, eta) * phiErrorScale;
67  double errCotTheta = this->errCotTheta(invPt, eta) * cotThetaErrorScale;
68  double errTip = this->errTip(invPt, eta) * tipErrorScale;
69  double errZip = this->errZip(invPt, eta) * zipErrorScale;
70 
71  Measurement1D pt(valPt, errPt);
72  Measurement1D phi(valPhi, errPhi);
73  Measurement1D cotTheta(valCotTheta, errCotTheta);
74  Measurement1D tip(valTip, errTip);
75  Measurement1D zip(valZip, errZip);
76 
77  float chi2 = 0.;
78 
79 /*
80  std::ostringstream str;
81  str <<"\t ipt: " << invPt <<"+/-"<<invPtErr
82  <<"\t pt: " << pt.value() <<"+/-"<<pt.error()
83  <<"\t phi: " << phi.value() <<"+/-"<<phi.error()
84  <<"\t cot: " << cotTheta.value() <<"+/-"<<cotTheta.error()
85  <<"\t tip: " << tip.value() <<"+/-"<<tip.error()
86  <<"\t zip: " << zip.value() <<"+/-"<<zip.error()
87  <<"\t charge: " << charge;
88  std::cout <<str.str()<< std::endl;
89 */
90 
91 
92  return builder.build(pt, phi, cotTheta, tip, zip, chi2, charge, hits, fieldESH.product());
93 }
94 
95 
96 double L1MuonPixelTrackFitter::valPhi(const Circle &circle, int charge) const
97 {
98  Circle::Point zero(0.,0.,0.);
99  Circle::Vector center = circle.center()-zero;
100  long double radius = center.perp();
101  center /= radius;
102  Circle::Vector dir = center.cross(-charge*Circle::Vector(0.,0.,1.));
103  return dir.phi();
104 }
105 
106 double L1MuonPixelTrackFitter::errPhi(double invPt, double eta) const {
107  //
108  // sigma = p0+p1/|pt|;
109  //
110  double p0,p1;
111  int ieta = int (10*fabs(eta));
112  switch (ieta) {
113  case 0: { p0 = 0.000597506; p1 = 0.00221057; break; }
114  case 1: { p0 = 0.000591867; p1 = 0.00278744; break; }
115  case 2: { p0 = 0.000635666; p1 = 0.00207433; break; }
116  case 3: { p0 = 0.000619086; p1 = 0.00255121; break; }
117  case 4: { p0 = 0.000572067; p1 = 0.00310618; break; }
118  case 5: { p0 = 0.000596239; p1 = 0.00288442; break; }
119  case 6: { p0 = 0.000607608; p1 = 0.00282996; break; }
120  case 7: { p0 = 0.000606446; p1 = 0.00281118; break; }
121  case 8: { p0 = 0.000594076; p1 = 0.00280546; break; }
122  case 9: { p0 = 0.000579615; p1 = 0.00335534; break; }
123  case 10: { p0 = 0.000659546; p1 = 0.00340443; break; }
124  case 11: { p0 = 0.000659031; p1 = 0.00343151; break; }
125  case 12: { p0 = 0.000738391; p1 = 0.00337297; break; }
126  case 13: { p0 = 0.000798966; p1 = 0.00330008; break; }
127  case 14: { p0 = 0.000702997; p1 = 0.00562643; break; }
128  case 15: { p0 = 0.000973417; p1 = 0.00312666; break; }
129  case 16: { p0 = 0.000995213; p1 = 0.00564278; break; }
130  case 17: { p0 = 0.00121436; p1 = 0.00572704; break; }
131  case 18: { p0 = 0.00119216; p1 = 0.00760204; break; }
132  case 19: { p0 = 0.00141204; p1 = 0.0093777; break; }
133  default: { p0 = 0.00153161; p1 = 0.00940265; break; }
134  }
135  return p0+p1*fabs(invPt);
136 }
137 
138 
140 {
141  return line.cotLine();
142 }
143 
144 double L1MuonPixelTrackFitter::errCotTheta(double invPt, double eta) const
145 {
146  //
147  // sigma = p0+p1/|pt|;
148  //
149  double p0,p1;
150  int ieta = int (10*fabs(eta));
151  switch (ieta) {
152  case 0: { p0 = 0.00166115; p1 = 5.75533e-05; break; }
153  case 1: { p0 = 0.00157525; p1 = 0.000707437; break; }
154  case 2: { p0 = 0.00122246; p1 = 0.000325456; break; }
155  case 3: { p0 = 0.000852422; p1 = 0.000429216; break; }
156  case 4: { p0 = 0.000637561; p1 = 0.00122298; break; }
157  case 5: { p0 = 0.000555766; p1 = 0.00158096; break; }
158  case 6: { p0 = 0.000641202; p1 = 0.00143339; break; }
159  case 7: { p0 = 0.000803207; p1 = 0.000648816; break; }
160  case 8: { p0 = 0.000741394; p1 = 0.0015289; break; }
161  case 9: { p0 = 0.000652019; p1 = 0.00168873; break; }
162  case 10: { p0 = 0.000716902; p1 = 0.00257556; break; }
163  case 11: { p0 = 0.000800409; p1 = 0.00190563; break; }
164  case 12: { p0 = 0.000808778; p1 = 0.00264139; break; }
165  case 13: { p0 = 0.000775757; p1 = 0.00318478; break; }
166  case 14: { p0 = 0.000705781; p1 = 0.00460576; break; }
167  case 15: { p0 = 0.000580679; p1 = 0.00748248; break; }
168  case 16: { p0 = 0.000561667; p1 = 0.00767487; break; }
169  case 17: { p0 = 0.000521626; p1 = 0.0100178; break; }
170  case 18: { p0 = 0.00064253; p1 = 0.0106062; break; }
171  case 19: { p0 = 0.000636868; p1 = 0.0140047; break; }
172  default: { p0 = 0.000682478; p1 = 0.0163569; break; }
173  }
174  return p0+p1*fabs(invPt);
175 }
176 
177 
178 double L1MuonPixelTrackFitter::valTip(const Circle &circle, double curvature) const
179 {
180  Circle::Point zero(0.,0.,0.);
181  Circle::Vector center = circle.center()-zero;
182  long double radius = center.perp();
183  return radius-1./fabs(curvature);
184 }
185 
186 double L1MuonPixelTrackFitter::errTip(double invPt, double eta) const {
187  //
188  // sigma = p0+p1/|pt|;
189  //
190  double p0,p1;
191  int ieta = int (10*fabs(eta));
192  switch (ieta) {
193  case 0: { p0 = 0.00392416; p1 = 0.00551809; break; }
194  case 1: { p0 = 0.00390391; p1 = 0.00543244; break; }
195  case 2: { p0 = 0.0040651; p1 = 0.00406496; break; }
196  case 3: { p0 = 0.00387782; p1 = 0.00797637; break; }
197  case 4: { p0 = 0.00376798; p1 = 0.00866894; break; }
198  case 5: { p0 = 0.0042131; p1 = 0.00462184; break; }
199  case 6: { p0 = 0.00392579; p1 = 0.00784685; break; }
200  case 7: { p0 = 0.00370472; p1 = 0.00790174; break; }
201  case 8: { p0 = 0.00364433; p1 = 0.00928368; break; }
202  case 9: { p0 = 0.00387578; p1 = 0.00640431; break; }
203  case 10: { p0 = 0.00382464; p1 = 0.00960763; break; }
204  case 11: { p0 = 0.0038907; p1 = 0.0104562; break; }
205  case 12: { p0 = 0.00392525; p1 = 0.0106442; break; }
206  case 13: { p0 = 0.00400634; p1 = 0.011218; break; }
207  case 14: { p0 = 0.0036229; p1 = 0.0156403; break; }
208  case 15: { p0 = 0.00444317; p1 = 0.00832987; break; }
209  case 16: { p0 = 0.00465492; p1 = 0.0179908; break; }
210  case 17: { p0 = 0.0049652; p1 = 0.0216647; break; }
211  case 18: { p0 = 0.0051395; p1 = 0.0233692; break; }
212  case 19: { p0 = 0.0062917; p1 = 0.0262175; break; }
213  default: { p0 = 0.00714444; p1 = 0.0253856; break; }
214  }
215  return p0+p1*fabs(invPt);
216 }
217 
218 
220  const GlobalPoint& pinner, const GlobalPoint& pouter) const
221 {
222  //
223  // phi = asin(r*rho/2) with asin(x) ~= x+x**3/(2*3)
224  //
225  double rho3 = curv*curv*curv;
226  double r1 = pinner.perp();
227  double phi1 = r1*curv/2 + pinner.perp2()*r1*rho3/48.;
228  double r2 = pouter.perp();
229  double phi2 = r2*curv/2 + pouter.perp2()*r2*rho3/48.;
230  double z1 = pinner.z();
231  double z2 = pouter.z();
232 
233  return z1 - phi1/(phi1-phi2)*(z1-z2);
234 }
235 
236 
237 double L1MuonPixelTrackFitter::errZip(double invPt, double eta) const {
238  //
239  // sigma = p0+p1/pt;
240  //
241  double p0,p1;
242  int ieta = int (10*fabs(eta));
243  switch (ieta) {
244  case 0: { p0 = 0.0120743; p1 = 0; break; }
245  case 1: { p0 = 0.0110343; p1 = 0.0051199; break; }
246  case 2: { p0 = 0.00846487; p1 = 0.00570084; break; }
247  case 3: { p0 = 0.00668726; p1 = 0.00331165; break; }
248  case 4: { p0 = 0.00467126; p1 = 0.00578239; break; }
249  case 5: { p0 = 0.0043042; p1 = 0.00598517; break; }
250  case 6: { p0 = 0.00515392; p1 = 0.00495422; break; }
251  case 7: { p0 = 0.0060843; p1 = 0.00320512; break; }
252  case 8: { p0 = 0.00564942; p1 = 0.00478876; break; }
253  case 9: { p0 = 0.00532111; p1 = 0.0073239; break; }
254  case 10: { p0 = 0.00579429; p1 = 0.00952782; break; }
255  case 11: { p0 = 0.00614229; p1 = 0.00977795; break; }
256  case 12: { p0 = 0.00714661; p1 = 0.00550482; break; }
257  case 13: { p0 = 0.0066593; p1 = 0.00999362; break; }
258  case 14: { p0 = 0.00634922; p1 = 0.0148156; break; }
259  case 15: { p0 = 0.00600586; p1 = 0.0318022; break; }
260  case 16: { p0 = 0.00676919; p1 = 0.027456; break; }
261  case 17: { p0 = 0.00670066; p1 = 0.0317005; break; }
262  case 18: { p0 = 0.00752392; p1 = 0.0347714; break; }
263  case 19: { p0 = 0.00791425; p1 = 0.0566665; break; }
264  default: { p0 = 0.00882372; p1 = 0.0596858; break; }
265  }
266  return p0+p1*fabs(invPt);
267 }
268 
269 
270 double L1MuonPixelTrackFitter::valInversePt( double phi0, double phiL1, double eta) const
271 {
272  double result = 0.;
273 
274  // solve equtaion p3*result^3 + p1*result + dphi = 0;
275  // where: result = 1/pt
276  // dphi = phiL1 - phi0
277  // Cardan way is used
278 
279  double p1,p2,p3; //parameters p1,p3 are negative by parameter fix in fit!
280  param(eta, p1, p2, p3);
281 
282  double dphi = deltaPhi(phiL1,phi0); // phi0-phiL1
283  if ( fabs(dphi) < 0.01) {
284  result = -dphi/p1;
285  } else {
286  double q = dphi/2./p3;
287  double p = p1/3./p3; // positive
288  double D = q*q + p*p*p; // positive
289  double u = pow(-q+sqrt(D), 1./3.);
290  double v = -pow(q+sqrt(D), 1./3.);
291  result = u+v;
292  }
293  return result;
294 }
295 
296 
297 double L1MuonPixelTrackFitter::errInversePt(double invPt, double eta) const {
298  //
299  // pt*sigma(1/pt) = p0+p1*pt;
300  //
301  double p0,p1;
302  int ieta = int (10*fabs(eta));
303  switch (ieta) {
304  case 0: { p0 = 0.0196835; p1 = 0.00517533; break; }
305  case 1: { p0 = 0.0266583; p1 = 0.00478101; break; }
306  case 2: { p0 = 0.0217164; p1 = 0.00545425; break; }
307  case 3: { p0 = 0.0197547; p1 = 0.00552263; break; }
308  case 4: { p0 = 0.0208778; p1 = 0.00536009; break; }
309  case 5: { p0 = 0.024192; p1 = 0.00521709; break; }
310  case 6: { p0 = 0.0265315; p1 = 0.0051897; break; }
311  case 7: { p0 = 0.0198071; p1 = 0.00566822; break; }
312  case 8: { p0 = 0.0361955; p1 = 0.00486352; break; }
313  case 9: { p0 = 0.037864; p1 = 0.00509094; break; }
314  case 10: { p0 = 0.0382968; p1 = 0.00612354; break; }
315  case 11: { p0 = 0.0308326; p1 = 0.0074234; break; }
316  case 12: { p0 = 0.0248577; p1 = 0.00883049; break; }
317  case 13: { p0 = 0.0279965; p1 = 0.00888293; break; }
318  case 14: { p0 = 0.0372582; p1 = 0.00950252; break; }
319  case 15: { p0 = 0.0281366; p1 = 0.0111501; break; }
320  case 16: { p0 = 0.0421483; p1 = 0.0109413; break; }
321  case 17: { p0 = 0.0461798; p1 = 0.0125824; break; }
322  case 18: { p0 = 0.0530603; p1 = 0.0132638; break; }
323  case 19: { p0 = 0.0601148; p1 = 0.0147911; break; }
324  default: { p0 = 0.0552377; p1 = 0.0155574; break; }
325  }
326  return p1+p0*fabs(invPt);
327 }
328 
329 double L1MuonPixelTrackFitter::findPt(double phi0, double phiL1, double eta, int charge) const {
330 
331  double dphi_min = fabs(deltaPhi(phi0,phiL1));
332  double pt_best = 1.;
333  double pt_cur = 1;
334  while ( pt_cur < 10000.) {
335  double phi_exp = phi0+getBending(1./pt_cur, eta, charge);
336  double dphi = fabs(deltaPhi(phi_exp,phiL1));
337  if ( dphi < dphi_min) {
338  pt_best = pt_cur;
339  dphi_min = dphi;
340  }
341  if (pt_cur < 10.) pt_cur += 0.01;
342  else if (pt_cur < 20.) pt_cur+= 0.025;
343  else if( pt_cur < 100.) pt_cur+= 0.1;
344  else pt_cur +=1;
345  };
346  return pt_best;
347 }
348 
349 
350 double L1MuonPixelTrackFitter::getBending(double invPt, double eta, int charge)
351 {
352  double p1, p2, p3;
353  param(eta,p1,p2,p3);
354  return charge*p1*invPt + charge*p2*invPt*invPt + charge*p3*invPt*invPt*invPt;
355 }
356 
357 double L1MuonPixelTrackFitter::getBendingError(double invPt, double eta)
358 {
359  int ieta = int (10*fabs(eta));
360  double p0,p1;
361  switch (ieta) {
362  case 0: { p0 = 0.0196835; p1 = 0.00517533; break; }
363  case 1: { p0 = 0.0266583; p1 = 0.00478101; break; }
364  case 2: { p0 = 0.0217164; p1 = 0.00545425; break; }
365  case 3: { p0 = 0.0197547; p1 = 0.00552263; break; }
366  case 4: { p0 = 0.0208778; p1 = 0.00536009; break; }
367  case 5: { p0 = 0.024192; p1 = 0.00521709; break; }
368  case 6: { p0 = 0.0265315; p1 = 0.0051897; break; }
369  case 7: { p0 = 0.0198071; p1 = 0.00566822; break; }
370  case 8: { p0 = 0.0361955; p1 = 0.00486352; break; }
371  case 9: { p0 = 0.037864; p1 = 0.00509094; break; }
372  case 10: { p0 = 0.0382968; p1 = 0.00612354; break; }
373  case 11: { p0 = 0.0308326; p1 = 0.0074234; break; }
374  case 12: { p0 = 0.0248577; p1 = 0.00883049; break; }
375  case 13: { p0 = 0.0279965; p1 = 0.00888293; break; }
376  case 14: { p0 = 0.0372582; p1 = 0.00950252; break; }
377  case 15: { p0 = 0.0281366; p1 = 0.0111501; break; }
378  case 16: { p0 = 0.0421483; p1 = 0.0109413; break; }
379  case 17: { p0 = 0.0461798; p1 = 0.0125824; break; }
380  case 18: { p0 = 0.0530603; p1 = 0.0132638; break; }
381  case 19: { p0 = 0.0601148; p1 = 0.0147911; break; }
382  default: { p0 = 0.0552377; p1 = 0.0155574; break; }
383  }
384  return p0+p1*sqr(invPt);
385 }
386 
387 
388 void L1MuonPixelTrackFitter::param(double eta, double &p1, double& p2, double& p3)
389 {
390  int ieta = int (10*fabs(eta));
391  switch (ieta) {
392  case 0: { p1 = -2.68016; p2 = 0; p3 = -12.9653; break; }
393  case 1: { p1 = -2.67864; p2 = 0; p3 = -12.0036; break; }
394  case 2: { p1 = -2.72997; p2 = 0; p3 = -10.3468; break; }
395  case 3: { p1 = -2.68836; p2 = 0; p3 = -12.3369; break; }
396  case 4: { p1 = -2.66885; p2 = 0; p3 = -11.589; break; }
397  case 5: { p1 = -2.64932; p2 = 0; p3 = -12.7176; break; }
398  case 6: { p1 = -2.64513; p2 = 0; p3 = -11.3912; break; }
399  case 7: { p1 = -2.64487; p2 = 0; p3 = -9.83239; break; }
400  case 8: { p1 = -2.58875; p2 = 0; p3 = -10.0541; break; }
401  case 9: { p1 = -2.5551; p2 = 0; p3 = -9.75757; break; }
402  case 10: { p1 = -2.55294; p2 = 0; p3 = -7.25238; break; }
403  case 11: { p1 = -2.45333; p2 = 0; p3 = -7.966; break; }
404  case 12: { p1 = -2.29283; p2 = 0; p3 = -7.03231; break; }
405  case 13: { p1 = -2.13167; p2 = 0; p3 = -4.29182; break; }
406  case 14: { p1 = -1.9102; p2 = 0; p3 = -3.21295; break; }
407  case 15: { p1 = -1.74552; p2 = 0; p3 = -0.531827; break; }
408  case 16: { p1 = -1.53642; p2 = 0; p3 = -1.74057; break; }
409  case 17: { p1 = -1.39446; p2 = 0; p3 = -1.56819; break; }
410  case 18: { p1 = -1.26176; p2 = 0; p3 = -0.598631; break; }
411  case 19: { p1 = -1.14133; p2 = 0; p3 = -0.0941055; break; }
412  default: { p1 = -1.02086; p2 = 0; p3 = -0.100491; break; }
413  }
414 }
415 
416 double L1MuonPixelTrackFitter::deltaPhi(double phi1, double phi2) const
417 {
418  while ( phi1 >= 2*M_PI) phi1 -= 2*M_PI;
419  while ( phi2 >= 2*M_PI) phi2 -= 2*M_PI;
420  while ( phi1 < 0) phi1 += 2*M_PI;
421  while ( phi2 < 0) phi2 += 2*M_PI;
422  double dPhi = phi2-phi1;
423 
424  if ( dPhi > M_PI ) dPhi -= 2*M_PI;
425  if ( dPhi < -M_PI ) dPhi += 2*M_PI;
426 
427  return dPhi;
428 }
429 
float etaValue() const
Definition: L1MuGMTCand.cc:114
double errZip(double invPt, double eta) const
T getParameter(std::string const &) const
double errTip(double invPt, double eta) const
float phiValue() const
Definition: L1MuGMTCand.cc:100
T perp() const
Definition: PV3DBase.h:72
void setPxConstraint(const SeedingHitSet &hits)
double findPt(double phi0, double phiL1, double eta, int charge) const
double valTip(const Circle &c, double curvature) const
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
static double getBendingError(double invPt, double eta)
T perp2() const
Definition: PV3DBase.h:71
double valPhi(const Circle &c, int charge) const
T sqr(T t)
float cotLine() const
T curvature(T InversePt, const edm::EventSetup &iSetup)
static void param(double eta, double &p1, double &p2, double &p3)
OutputIterator zip(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2, OutputIterator result, Compare comp)
T sqrt(T t)
Definition: SSEVec.h:18
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
virtual reco::Track * run(const edm::EventSetup &es, const std::vector< const TrackingRecHit * > &hits, const TrackingRegion &region) const
T z() const
Definition: PV3DBase.h:64
reco::Track * build(const Measurement1D &pt, const Measurement1D &phi, const Measurement1D &cotTheta, const Measurement1D &tip, const Measurement1D &zip, float chi2, int charge, const std::vector< const TrackingRecHit * > &hits, const MagneticField *mf, const GlobalPoint &reference=GlobalPoint(0, 0, 0)) const
double valZip(double curvature, const GlobalPoint &p0, const GlobalPoint &p1) const
static double getBending(double invPt, double eta, int charge)
double p2[4]
Definition: TauolaWrapper.h:90
#define M_PI
double deltaPhi(double phi1, double phi2) const
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:151
double valInversePt(double phi0, double phiL1, double eta) const
const T & get() const
Definition: EventSetup.h:55
double errInversePt(double invPt, double eta) const
double p1[4]
Definition: TauolaWrapper.h:89
void setL1Constraint(const L1MuGMTCand &muon)
double valCotTheta(const PixelRecoLineRZ &line) const
dbl *** dir
Definition: mlp_gen.cc:35
long double T
double errPhi(double invPt, double eta) const
T const * product() const
Definition: ESHandle.h:86
L1MuonPixelTrackFitter(const edm::ParameterSet &cfg)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
double errCotTheta(double invPt, double eta) const
double p3[4]
Definition: TauolaWrapper.h:91
int charge() const
get charge (+1 -1)
Definition: L1MuGMTCand.h:135