CMS 3D CMS Logo

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