CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonSeedCreator.cc
Go to the documentation of this file.
1 
10 
12 
14 
17 
20 
26 
30 
31 #include "gsl/gsl_statistics.h"
32 
33 
34 /*
35  * Constructor
36  */
38 
39  theMinMomentum = pset.getParameter<double>("minimumSeedPt");
40  theMaxMomentum = pset.getParameter<double>("maximumSeedPt");
41  defaultMomentum = pset.getParameter<double>("defaultSeedPt");
42  debug = pset.getParameter<bool>("DebugMuonSeed");
43  sysError = pset.getParameter<double>("SeedPtSystematics");
44  // load seed PT parameters
45  DT12 = pset.getParameter<std::vector<double> >("DT_12");
46  DT13 = pset.getParameter<std::vector<double> >("DT_13");
47  DT14 = pset.getParameter<std::vector<double> >("DT_14");
48  DT23 = pset.getParameter<std::vector<double> >("DT_23");
49  DT24 = pset.getParameter<std::vector<double> >("DT_24");
50  DT34 = pset.getParameter<std::vector<double> >("DT_34");
51 
52  CSC01 = pset.getParameter<std::vector<double> >("CSC_01");
53  CSC12 = pset.getParameter<std::vector<double> >("CSC_12");
54  CSC02 = pset.getParameter<std::vector<double> >("CSC_02");
55  CSC13 = pset.getParameter<std::vector<double> >("CSC_13");
56  CSC03 = pset.getParameter<std::vector<double> >("CSC_03");
57  CSC14 = pset.getParameter<std::vector<double> >("CSC_14");
58  CSC23 = pset.getParameter<std::vector<double> >("CSC_23");
59  CSC24 = pset.getParameter<std::vector<double> >("CSC_24");
60  CSC34 = pset.getParameter<std::vector<double> >("CSC_34");
61 
62  OL1213 = pset.getParameter<std::vector<double> >("OL_1213");
63  OL1222 = pset.getParameter<std::vector<double> >("OL_1222");
64  OL1232 = pset.getParameter<std::vector<double> >("OL_1232");
65  OL1213 = pset.getParameter<std::vector<double> >("OL_1213");
66  OL2222 = pset.getParameter<std::vector<double> >("OL_1222");
67 
68  SME11 = pset.getParameter<std::vector<double> >("SME_11");
69  SME12 = pset.getParameter<std::vector<double> >("SME_12");
70  SME13 = pset.getParameter<std::vector<double> >("SME_13");
71  SME21 = pset.getParameter<std::vector<double> >("SME_21");
72  SME22 = pset.getParameter<std::vector<double> >("SME_22");
73  SME31 = pset.getParameter<std::vector<double> >("SME_31");
74  SME32 = pset.getParameter<std::vector<double> >("SME_32");
75  SME41 = pset.getParameter<std::vector<double> >("SME_41");
76 
77  SMB10 = pset.getParameter<std::vector<double> >("SMB_10");
78  SMB11 = pset.getParameter<std::vector<double> >("SMB_11");
79  SMB12 = pset.getParameter<std::vector<double> >("SMB_12");
80  SMB20 = pset.getParameter<std::vector<double> >("SMB_20");
81  SMB21 = pset.getParameter<std::vector<double> >("SMB_21");
82  SMB22 = pset.getParameter<std::vector<double> >("SMB_22");
83  SMB30 = pset.getParameter<std::vector<double> >("SMB_30");
84  SMB31 = pset.getParameter<std::vector<double> >("SMB_31");
85  SMB32 = pset.getParameter<std::vector<double> >("SMB_32");
86 
87  // Load dphi scale parameters
88  CSC01_1 = pset.getParameter<std::vector<double> >("CSC_01_1_scale");
89  CSC12_1 = pset.getParameter<std::vector<double> >("CSC_12_1_scale");
90  CSC12_2 = pset.getParameter<std::vector<double> >("CSC_12_2_scale");
91  CSC12_3 = pset.getParameter<std::vector<double> >("CSC_12_3_scale");
92  CSC13_2 = pset.getParameter<std::vector<double> >("CSC_13_2_scale");
93  CSC13_3 = pset.getParameter<std::vector<double> >("CSC_13_3_scale");
94  CSC14_3 = pset.getParameter<std::vector<double> >("CSC_14_3_scale");
95  CSC23_1 = pset.getParameter<std::vector<double> >("CSC_23_1_scale");
96  CSC23_2 = pset.getParameter<std::vector<double> >("CSC_23_2_scale");
97  CSC24_1 = pset.getParameter<std::vector<double> >("CSC_24_1_scale");
98  CSC34_1 = pset.getParameter<std::vector<double> >("CSC_34_1_scale");
99 
100  DT12_1 = pset.getParameter<std::vector<double> >("DT_12_1_scale");
101  DT12_2 = pset.getParameter<std::vector<double> >("DT_12_2_scale");
102  DT13_1 = pset.getParameter<std::vector<double> >("DT_13_1_scale");
103  DT13_2 = pset.getParameter<std::vector<double> >("DT_13_2_scale");
104  DT14_1 = pset.getParameter<std::vector<double> >("DT_14_1_scale");
105  DT14_2 = pset.getParameter<std::vector<double> >("DT_14_2_scale");
106  DT23_1 = pset.getParameter<std::vector<double> >("DT_23_1_scale");
107  DT23_2 = pset.getParameter<std::vector<double> >("DT_23_2_scale");
108  DT24_1 = pset.getParameter<std::vector<double> >("DT_24_1_scale");
109  DT24_2 = pset.getParameter<std::vector<double> >("DT_24_2_scale");
110  DT34_1 = pset.getParameter<std::vector<double> >("DT_34_1_scale");
111  DT34_2 = pset.getParameter<std::vector<double> >("DT_34_2_scale");
112 
113  OL_1213 = pset.getParameter<std::vector<double> >("OL_1213_0_scale");
114  OL_1222 = pset.getParameter<std::vector<double> >("OL_1222_0_scale");
115  OL_1232 = pset.getParameter<std::vector<double> >("OL_1232_0_scale");
116  OL_2213 = pset.getParameter<std::vector<double> >("OL_2213_0_scale");
117  OL_2222 = pset.getParameter<std::vector<double> >("OL_2222_0_scale");
118 
119  SMB_10S = pset.getParameter<std::vector<double> >("SMB_10_0_scale");
120  SMB_11S = pset.getParameter<std::vector<double> >("SMB_11_0_scale");
121  SMB_12S = pset.getParameter<std::vector<double> >("SMB_12_0_scale");
122  SMB_20S = pset.getParameter<std::vector<double> >("SMB_20_0_scale");
123  SMB_21S = pset.getParameter<std::vector<double> >("SMB_21_0_scale");
124  SMB_22S = pset.getParameter<std::vector<double> >("SMB_22_0_scale");
125  SMB_30S = pset.getParameter<std::vector<double> >("SMB_30_0_scale");
126  SMB_31S = pset.getParameter<std::vector<double> >("SMB_31_0_scale");
127  SMB_32S = pset.getParameter<std::vector<double> >("SMB_32_0_scale");
128 
129  SME_11S = pset.getParameter<std::vector<double> >("SME_11_0_scale");
130  SME_12S = pset.getParameter<std::vector<double> >("SME_12_0_scale");
131  SME_13S = pset.getParameter<std::vector<double> >("SME_13_0_scale");
132  SME_21S = pset.getParameter<std::vector<double> >("SME_21_0_scale");
133  SME_22S = pset.getParameter<std::vector<double> >("SME_22_0_scale");
134 
135 }
136 
137 
138 /*
139  * Destructor
140  */
142 
143 }
144 
145 
146 /*
147  * createSeed
148  *
149  * Note type = 1 --> CSC
150  * = 2 --> Overlap
151  * = 3 --> DT
152  */
153 TrajectorySeed MuonSeedCreator::createSeed(int type, SegmentContainer seg, std::vector<int> layers, int NShowers, int NShowerSegments ) {
154 
155  // The index of the station closest to the IP
156  int last = 0;
157 
158  double ptmean = theMinMomentum;
159  double sptmean = theMinMomentum;
160 
161  AlgebraicVector t(4);
162  AlgebraicSymMatrix mat(5,0) ;
163  LocalPoint segPos;
164 
165  // Compute the pt according to station types used;
166  if (type == 1 ) estimatePtCSC(seg, layers, ptmean, sptmean);
167  if (type == 2 ) estimatePtOverlap(seg, layers, ptmean, sptmean);
168  if (type == 3 ) estimatePtDT(seg, layers, ptmean, sptmean);
169  if (type == 4 ) estimatePtSingle(seg, layers, ptmean, sptmean);
170  // type 5 are the seeding for ME1/4
171  if (type == 5 ) estimatePtCSC(seg, layers, ptmean, sptmean);
172 
173  // for certain clear showering case, set-up the minimum value
174  if ( NShowers > 0 ) estimatePtShowering( NShowers, NShowerSegments, ptmean, sptmean );
175  //if ( NShowers > 0 ) std::cout<<" Showering happened "<<NShowers<<" times w/ "<< NShowerSegments<<std::endl; ;
176 
177 
178  // Minimal pt
179  double charge = 1.0;
180  if (ptmean < 0.) charge = -1.0;
181  if ( (charge * ptmean) < theMinMomentum ) {
182  ptmean = theMinMomentum * charge;
183  sptmean = theMinMomentum ;
184  }
185  else if ( (charge * ptmean) > theMaxMomentum ) {
186  ptmean = theMaxMomentum * charge;
187  sptmean = theMaxMomentum * 0.25 ;
188  }
189 
191 
192  double p_err =0.0;
193  // determine the seed layer
194  int best_seg= 0;
195  double chi2_dof = 9999.0;
196  unsigned int ini_seg = 0;
197  // avoid generating seed from 1st layer(ME1/1)
198  if ( type == 5 ) ini_seg = 1;
199  for (size_t i = ini_seg ; i < seg.size(); i++) {
200  double dof = static_cast<double>(seg[i]->degreesOfFreedom());
201  if ( chi2_dof < ( seg[i]->chi2()/dof ) ) continue;
202  chi2_dof = seg[i]->chi2() / dof ;
203  best_seg = static_cast<int>(i);
204  }
205 
206 
207  if ( type==1 || type==5 || type== 4) {
208  // Fill the LocalTrajectoryParameters
210  last = best_seg;
211  // last = 0;
212  GlobalVector mom = seg[last]->globalPosition()-GlobalPoint();
213  segPos = seg[last]->localPosition();
215 
216  GlobalVector polar(GlobalVector::Spherical(mom.theta(),seg[last]->globalDirection().phi(),1.));
217 
219  polar *= fabs(ptmean)/polar.perp();
221  LocalVector segDirFromPos = seg[last]->det()->toLocal(polar);
222  int chargeI = static_cast<int>(charge);
223  LocalTrajectoryParameters param1(segPos, segDirFromPos, chargeI);
224  param = param1;
225  p_err = (sptmean*sptmean)/(polar.mag()*polar.mag()*ptmean*ptmean) ;
226  mat = seg[last]->parametersError().similarityT( seg[last]->projectionMatrix() );
227  mat[0][0]= p_err;
228  if ( type == 5 ) {
229  mat[0][0] = mat[0][0]/fabs( tan(mom.theta()) );
230  mat[1][1] = mat[1][1]/fabs( tan(mom.theta()) );
231  mat[3][3] = 2.25*mat[3][3];
232  mat[4][4] = 2.25*mat[4][4];
233  }
234  if ( type == 4 ) {
235  mat[0][0] = mat[0][0]/fabs( tan(mom.theta()) );
236  mat[1][1] = mat[1][1]/fabs( tan(mom.theta()) );
237  mat[2][2] = 2.25*mat[2][2];
238  mat[3][3] = 2.25*mat[3][3];
239  mat[4][4] = 2.25*mat[4][4];
240  }
241  double dh = fabs( seg[last]->globalPosition().eta() ) - 1.6 ;
242  if ( fabs(dh) < 0.1 && type == 1 ) {
243  mat[1][1] = 4.*mat[1][1];
244  mat[2][2] = 4.*mat[2][2];
245  mat[3][3] = 9.*mat[3][3];
246  mat[4][4] = 9.*mat[4][4];
247  }
248 
249  //if ( !highPt && type != 1 ) mat[1][1]= 2.25*mat[1][1];
250  //if ( highPt && type != 1 ) mat[3][3]= 2.25*mat[1][1];
251  //mat[2][2]= 3.*mat[2][2];
252  //mat[3][3]= 2.*mat[3][3];
253  //mat[4][4]= 2.*mat[4][4];
254  }
255  else {
256  // Fill the LocalTrajectoryParameters
258  last = 0;
259  segPos = seg[last]->localPosition();
260  GlobalVector mom = seg[last]->globalPosition()-GlobalPoint();
262  GlobalVector polar(GlobalVector::Spherical(mom.theta(),seg[last]->globalDirection().phi(),1.));
263  //GlobalVector polar(GlobalVector::Spherical(seg[last]->globalDirection().theta(),seg[last]->globalDirection().phi(),1.));
264 
266  //double ptRatio = 1.0 - (2.808/(fabs(ptmean) -1)) + (4.546/( (fabs(ptmean)-1)*(fabs(ptmean)-1)) );
267  //ptmean = ptmean*ptRatio ;
268 
270  polar *= fabs(ptmean)/polar.perp();
272  LocalVector segDirFromPos = seg[last]->det()->toLocal(polar);
273  int chargeI = static_cast<int>(charge);
274  LocalTrajectoryParameters param1(segPos, segDirFromPos, chargeI);
275  param = param1;
276  p_err = (sptmean*sptmean)/(polar.mag()*polar.mag()*ptmean*ptmean) ;
277  mat = seg[last]->parametersError().similarityT( seg[last]->projectionMatrix() );
278  //mat[0][0]= 1.44 * p_err;
279  mat[0][0]= p_err;
280  }
281 
282  if ( debug ) {
283  GlobalPoint gp = seg[last]->globalPosition();
284  float Geta = gp.eta();
285 
286  std::cout << "Type " << type << " Nsegments " << layers.size() << " ";
287  std::cout << "pt " << ptmean << " errpt " << sptmean
288  << " eta " << Geta << " charge " << charge
289  << std::endl;
290  }
291 
292  // this perform H.T() * parErr * H, which is the projection of the
293  // the measurement error (rechit rf) to the state error (TSOS rf)
294  // Legend:
295  // H => is the 4x5 projection matrix
296  // parError the 4x4 parameter error matrix of the RecHit
297 
298 
300 
301  // Create the TrajectoryStateOnSurface
302  TrajectoryStateOnSurface tsos(param, error, seg[last]->det()->surface(),&*BField);
303 
304  // Take the DetLayer on which relies the segment
305  DetId id = seg[last]->geographicalId();
306 
307  // Transform it in a TrajectoryStateOnSurface
308  TrajectoryStateTransform tsTransform;
309 
310  //PTrajectoryStateOnDet *seedTSOS = tsTransform.persistentState( tsos, id.rawId());
311  std::auto_ptr<PTrajectoryStateOnDet> seedTSOS(tsTransform.persistentState( tsos, id.rawId()));
312 
314  for (unsigned l=0; l<seg.size(); l++) {
315  container.push_back( seg[l]->hit()->clone() );
316  //container.push_back(seg[l]->hit());
317  }
318 
319  TrajectorySeed theSeed(*seedTSOS,container,alongMomentum);
320  return theSeed;
321 }
322 
323 
324 /*
325  * estimatePtCSC
326  *
327  * Look at delta phi to determine pt as:
328  * pt = (c_1 * eta + c_2) / dphi
329  *
330  * Note that only segment pairs with one segment in the ME1 layers give sensitive results
331  *
332  */
333 void MuonSeedCreator::estimatePtCSC(SegmentContainer seg, std::vector<int> layers, double& thePt, double& theSpt ) {
334 
335  unsigned size = seg.size();
336  if (size < 2) return;
337 
338  // reverse the segment and layer container first for pure CSC case
339  //if ( layers[0] > layers[ layers.size()-1 ] ) {
340  // reverse( layers.begin(), layers.end() );
341  // reverse( seg.begin(), seg.end() );
342  //}
343 
344  std::vector<double> ptEstimate;
345  std::vector<double> sptEstimate;
346 
347  thePt = defaultMomentum;
348  theSpt = defaultMomentum;
349 
350  double pt = 0.;
351  double spt = 0.;
352  GlobalPoint segPos[2];
353 
354  int layer0 = layers[0];
355  segPos[0] = seg[0]->globalPosition();
356  float eta = fabs( segPos[0].eta() );
357  //float corr = fabs( tan(segPos[0].theta()) );
358  // use pt from vertex information
359  /*
360  if ( layer0 == 0 ) {
361  SegmentContainer seg0;
362  seg0.push_back(seg[0]);
363  std::vector<int> lyr0(1,0);
364  estimatePtSingle( seg0, lyr0, thePt, theSpt);
365  ptEstimate.push_back( thePt );
366  sptEstimate.push_back( theSpt );
367  }
368  */
369 
370  //std::cout<<" estimate CSC "<<std::endl;
371 
372  unsigned idx1 = size;
373  if (size > 1) {
374  while ( idx1 > 1 ) {
375  idx1--;
376  int layer1 = layers[idx1];
377  if (layer0 == layer1) continue;
378  segPos[1] = seg[idx1]->globalPosition();
379 
380  double dphi = segPos[0].phi() - segPos[1].phi();
381  //double temp_dphi = dphi/corr;
382  double temp_dphi = dphi;
383 
384  double sign = 1.0;
385  if (temp_dphi < 0.) {
386  temp_dphi = -1.0*temp_dphi;
387  sign = -1.0;
388  }
389 
390  // Ensure that delta phi is not too small to prevent pt from blowing up
391  if (temp_dphi < 0.0001 ) {
392  temp_dphi = 0.0001 ;
393  pt = theMaxMomentum ;
394  spt = theMaxMomentum*0.25 ;
395  ptEstimate.push_back( pt*sign );
396  sptEstimate.push_back( spt );
397  }
398  // ME1 is inner-most
399  if ( layer0 == 0 && temp_dphi >= 0.0001 ) {
400 
401  // ME1/2 is outer-most
402  if ( layer1 == 1 ) {
403  //temp_dphi = scaledPhi(temp_dphi, CSC01_1[3] );
404  pt = getPt( CSC01, eta , temp_dphi )[0];
405  spt = getPt( CSC01, eta , temp_dphi )[1];
406  }
407  // ME2 is outer-most
408  else if ( layer1 == 2 ) {
409  //temp_dphi = scaledPhi(temp_dphi, CSC12_3[3] );
410  pt = getPt( CSC02, eta , temp_dphi )[0];
411  spt = getPt( CSC02, eta , temp_dphi )[1];
412  }
413  // ME3 is outer-most
414  else if ( layer1 == 3 ) {
415  //temp_dphi = scaledPhi(temp_dphi, CSC13_3[3] );
416  pt = getPt( CSC03, eta , temp_dphi )[0];
417  spt = getPt( CSC03, eta , temp_dphi )[1];
418  }
419  // ME4 is outer-most
420  else {
421  //temp_dphi = scaledPhi(temp_dphi, CSC14_3[3]);
422  pt = getPt( CSC14, eta , temp_dphi )[0];
423  spt = getPt( CSC14, eta , temp_dphi )[1];
424  }
425  ptEstimate.push_back( pt*sign );
426  sptEstimate.push_back( spt );
427  }
428 
429  // ME1/2,ME1/3 is inner-most
430  if ( layer0 == 1 ) {
431  // ME2 is outer-most
432  if ( layer1 == 2 ) {
433 
434  //if ( eta <= 1.2 ) { temp_dphi = scaledPhi(temp_dphi, CSC12_1[3]); }
435  //if ( eta > 1.2 ) { temp_dphi = scaledPhi(temp_dphi, CSC12_2[3]); }
436  pt = getPt( CSC12, eta , temp_dphi )[0];
437  spt = getPt( CSC12, eta , temp_dphi )[1];
438  }
439  // ME3 is outer-most
440  else if ( layer1 == 3 ) {
441  temp_dphi = scaledPhi(temp_dphi, CSC13_2[3]);
442  pt = getPt( CSC13, eta , temp_dphi )[0];
443  spt = getPt( CSC13, eta , temp_dphi )[1];
444  }
445  // ME4 is outer-most
446  else {
447  temp_dphi = scaledPhi(temp_dphi, CSC14_3[3]);
448  pt = getPt( CSC14, eta , temp_dphi )[0];
449  spt = getPt( CSC14, eta , temp_dphi )[1];
450  }
451  ptEstimate.push_back( pt*sign );
452  sptEstimate.push_back( spt );
453  }
454 
455  // ME2 is inner-most
456  if ( layer0 == 2 && temp_dphi > 0.0001 ) {
457 
458  // ME4 is outer-most
459  bool ME4av =false;
460  if ( layer1 == 4 ) {
461  temp_dphi = scaledPhi(temp_dphi, CSC24_1[3]);
462  pt = getPt( CSC24, eta , temp_dphi )[0];
463  spt = getPt( CSC24, eta , temp_dphi )[1];
464  ME4av = true;
465  }
466  // ME3 is outer-most
467  else {
468  // if ME2-4 is availabe , discard ME2-3
469  if ( !ME4av ) {
470  if ( eta <= 1.7 ) { temp_dphi = scaledPhi(temp_dphi, CSC23_1[3]); }
471  if ( eta > 1.7 ) { temp_dphi = scaledPhi(temp_dphi, CSC23_2[3]); }
472  pt = getPt( CSC23, eta , temp_dphi )[0];
473  spt = getPt( CSC23, eta , temp_dphi )[1];
474  }
475  }
476  ptEstimate.push_back( pt*sign );
477  sptEstimate.push_back( spt );
478  }
479 
480  // ME3 is inner-most
481  if ( layer0 == 3 && temp_dphi > 0.0001 ) {
482 
483  temp_dphi = scaledPhi(temp_dphi, CSC34_1[3]);
484  pt = getPt( CSC34, eta , temp_dphi )[0];
485  spt = getPt( CSC34, eta , temp_dphi )[1];
486  ptEstimate.push_back( pt*sign );
487  sptEstimate.push_back( spt );
488  }
489 
490  }
491  }
492 
493  // Compute weighted average if have more than one estimator
494  if ( ptEstimate.size() > 0 ) weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
495 
496 }
497 
498 
499 /*
500  * estimatePtDT
501  *
502  * Look at delta phi between segments to determine pt as:
503  * pt = (c_1 * eta + c_2) / dphi
504  */
505 void MuonSeedCreator::estimatePtDT(SegmentContainer seg, std::vector<int> layers, double& thePt, double& theSpt) {
506 
507  unsigned size = seg.size();
508  if (size < 2) return;
509 
510  std::vector<double> ptEstimate;
511  std::vector<double> sptEstimate;
512 
513  thePt = defaultMomentum;
514  theSpt = defaultMomentum;
515 
516  double pt = 0.;
517  double spt = 0.;
518  GlobalPoint segPos[2];
519 
520  int layer0 = layers[0];
521  segPos[0] = seg[0]->globalPosition();
522  float eta = fabs(segPos[0].eta());
523 
524  //std::cout<<" estimate DT "<<std::endl;
525  // inner-most layer
526  //for ( unsigned idx0 = 0; idx0 < size-1; ++idx0 ) {
527  // layer0 = layers[idx0];
528  // segPos[0] = seg[idx0]->globalPosition();
529  // outer-most layer
530  // for ( unsigned idx1 = idx0+1; idx1 < size; ++idx1 ) {
531  for ( unsigned idx1 = 1; idx1 <size ; ++idx1 ) {
532 
533  int layer1 = layers[idx1];
534  segPos[1] = seg[idx1]->globalPosition();
535 
536  //eta = fabs(segPos[1].eta());
537  //if (layer1 == -4) eta = fabs(segPos[0].eta());
538 
539  double dphi = segPos[0].phi() - segPos[1].phi();
540  double temp_dphi = dphi;
541 
542  // Ensure that delta phi is not too small to prevent pt from blowing up
543 
544  double sign = 1.0;
545  if (temp_dphi < 0.) {
546  temp_dphi = -temp_dphi;
547  sign = -1.0;
548  }
549 
550  if (temp_dphi < 0.0001 ) {
551  temp_dphi = 0.0001 ;
552  pt = theMaxMomentum ;
553  spt = theMaxMomentum*0.25 ;
554  ptEstimate.push_back( pt*sign );
555  sptEstimate.push_back( spt );
556  }
557 
558  // MB1 is inner-most
559  bool MB23av = false;
560  if (layer0 == -1 && temp_dphi > 0.0001 ) {
561  // MB2 is outer-most
562  if (layer1 == -2) {
563 
564  if ( eta <= 0.7 ) { temp_dphi = scaledPhi(temp_dphi, DT12_1[3]); }
565  if ( eta > 0.7 ) { temp_dphi = scaledPhi(temp_dphi, DT12_2[3]); }
566  pt = getPt( DT12, eta , temp_dphi )[0];
567  spt = getPt( DT12, eta , temp_dphi )[1];
568  MB23av = true;
569  }
570  // MB3 is outer-most
571  else if (layer1 == -3) {
572 
573  if ( eta <= 0.6 ) { temp_dphi = scaledPhi(temp_dphi, DT13_1[3]); }
574  if ( eta > 0.6 ) { temp_dphi = scaledPhi(temp_dphi, DT13_2[3]); }
575  pt = getPt( DT13, eta , temp_dphi )[0];
576  spt = getPt( DT13, eta , temp_dphi )[1];
577  MB23av = true;
578  }
579  // MB4 is outer-most
580  else {
581  if ( !MB23av ) {
582  if ( eta <= 0.52 ) { temp_dphi = scaledPhi(temp_dphi, DT14_1[3]); }
583  if ( eta > 0.52 ) { temp_dphi = scaledPhi(temp_dphi, DT14_2[3]); }
584  pt = getPt( DT14, eta , temp_dphi )[0];
585  spt = getPt( DT14, eta , temp_dphi )[1];
586  }
587  }
588  ptEstimate.push_back( pt*sign );
589  sptEstimate.push_back( spt );
590  }
591 
592  // MB2 is inner-most
593  if (layer0 == -2 && temp_dphi > 0.0001 ) {
594  // MB3 is outer-most
595  if ( layer1 == -3) {
596 
597  if ( eta <= 0.6 ) { temp_dphi = scaledPhi(temp_dphi, DT23_1[3]); }
598  if ( eta > 0.6 ) { temp_dphi = scaledPhi(temp_dphi, DT23_2[3]); }
599  pt = getPt( DT23, eta , temp_dphi )[0];
600  spt = getPt( DT23, eta , temp_dphi )[1];
601  }
602  // MB4 is outer-most
603  else {
604 
605  if ( eta <= 0.52 ) { temp_dphi = scaledPhi(temp_dphi, DT24_1[3]); }
606  if ( eta > 0.52 ) { temp_dphi = scaledPhi(temp_dphi, DT24_2[3]); }
607  pt = getPt( DT24, eta , temp_dphi )[0];
608  spt = getPt( DT24, eta , temp_dphi )[1];
609  }
610  ptEstimate.push_back( pt*sign );
611  sptEstimate.push_back( spt );
612  }
613 
614  // MB3 is inner-most -> only marginally useful to pick up the charge
615  if (layer0 == -3 && temp_dphi > 0.0001 ) {
616  // MB4 is outer-most
617 
618  if ( eta <= 0.51 ) { temp_dphi = scaledPhi(temp_dphi, DT34_1[3]); }
619  if ( eta > 0.51 ) { temp_dphi = scaledPhi(temp_dphi, DT34_2[3]); }
620  pt = getPt( DT34, eta , temp_dphi )[0];
621  spt = getPt( DT34, eta , temp_dphi )[1];
622  ptEstimate.push_back( pt*sign );
623  sptEstimate.push_back( spt );
624  }
625  }
626  //}
627 
628 
629  // Compute weighted average if have more than one estimator
630  if (ptEstimate.size() > 0 ) weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
631 
632 }
633 
634 
635 /*
636  * estimatePtOverlap
637  *
638  */
639 void MuonSeedCreator::estimatePtOverlap(SegmentContainer seg, std::vector<int> layers, double& thePt, double& theSpt) {
640 
641  int size = layers.size();
642 
643  thePt = defaultMomentum;
644  theSpt = defaultMomentum;
645 
646  SegmentContainer segCSC;
647  std::vector<int> layersCSC;
648  SegmentContainer segDT;
649  std::vector<int> layersDT;
650 
651  // DT layers are numbered as -4 to -1, whereas CSC layers go from 0 to 4:
652  for ( unsigned j = 0; j < layers.size(); ++j ) {
653  if ( layers[j] > -1 ) {
654  segCSC.push_back(seg[j]);
655  layersCSC.push_back(layers[j]);
656  }
657  else {
658  segDT.push_back(seg[j]);
659  layersDT.push_back(layers[j]);
660  }
661  }
662 
663  std::vector<double> ptEstimate;
664  std::vector<double> sptEstimate;
665 
666  GlobalPoint segPos[2];
667  int layer0 = layers[0];
668  segPos[0] = seg[0]->globalPosition();
669  float eta = fabs(segPos[0].eta());
670  //std::cout<<" estimate OL "<<std::endl;
671 
672  if ( segDT.size() > 0 && segCSC.size() > 0 ) {
673  int layer1 = layers[size-1];
674  segPos[1] = seg[size-1]->globalPosition();
675 
676  double dphi = segPos[0].phi() - segPos[1].phi();
677  double temp_dphi = dphi;
678 
679  // Ensure that delta phi is not too small to prevent pt from blowing up
680 
681  double sign = 1.0;
682  if (temp_dphi < 0.) {
683  temp_dphi = -temp_dphi;
684  sign = -1.0;
685  }
686 
687  if (temp_dphi < 0.0001 ) {
688  temp_dphi = 0.0001 ;
689  thePt = theMaxMomentum ;
690  theSpt = theMaxMomentum*0.25 ;
691  ptEstimate.push_back( thePt*sign );
692  sptEstimate.push_back( theSpt );
693  }
694 
695  // MB1 is inner-most
696  if ( layer0 == -1 && temp_dphi > 0.0001 ) {
697  // ME1/3 is outer-most
698  if ( layer1 == 1 ) {
699  temp_dphi = scaledPhi(temp_dphi, OL_1213[3]);
700  thePt = getPt( OL1213, eta , temp_dphi )[0];
701  theSpt = getPt( OL1213, eta , temp_dphi )[1];
702  }
703  // ME2 is outer-most
704  else if ( layer1 == 2) {
705  temp_dphi = scaledPhi(temp_dphi, OL_1222[3]);
706  thePt = getPt( OL1222, eta , temp_dphi )[0];
707  theSpt = getPt( OL1222, eta , temp_dphi )[1];
708  }
709  // ME3 is outer-most
710  else {
711  temp_dphi = scaledPhi(temp_dphi, OL_1232[3]);
712  thePt = getPt( OL1232, eta , temp_dphi )[0];
713  theSpt = getPt( OL1232, eta , temp_dphi )[1];
714  }
715  ptEstimate.push_back(thePt*sign);
716  sptEstimate.push_back(theSpt);
717  }
718  // MB2 is inner-most
719  if ( layer0 == -2 && temp_dphi > 0.0001 ) {
720  // ME1/3 is outer-most
721  if ( layer1 == 1 ) {
722  temp_dphi = scaledPhi(temp_dphi, OL_2213[3]);
723  thePt = getPt( OL2213, eta , temp_dphi )[0];
724  theSpt = getPt( OL2213, eta , temp_dphi )[1];
725  ptEstimate.push_back(thePt*sign);
726  sptEstimate.push_back(theSpt);
727  }
728  // ME2 is outer-most
729  if ( layer1 == 2) {
730  temp_dphi = scaledPhi(temp_dphi, OL_2222[3]);
731  thePt = getPt( OL2222, eta , temp_dphi )[0];
732  theSpt = getPt( OL2222, eta , temp_dphi )[1];
733  }
734  }
735  }
736 
737  if ( segDT.size() > 1 ) {
738  estimatePtDT(segDT, layersDT, thePt, theSpt);
739  ptEstimate.push_back(thePt);
740  sptEstimate.push_back(theSpt);
741  }
742 
743  /*
744  // not useful ....and pt estimation is bad
745  if ( segCSC.size() > 1 ) {
746  // don't estimate pt without ME1 information
747  bool CSCLayer1=false;
748  for (unsigned i=0; i< layersCSC.size(); i++) {
749  if ( layersCSC[i]==0 || layersCSC[i]==1 ) CSCLayer1 = true;
750  }
751  if (CSCLayer1) {
752  estimatePtCSC(segCSC, layersCSC, thePt, theSpt);
753  ptEstimate.push_back(thePt);
754  sptEstimate.push_back(theSpt);
755  }
756  }
757  */
758 
759  // Compute weighted average if have more than one estimator
760  if (ptEstimate.size() > 0 ) weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
761 
762 }
763 /*
764  *
765  * estimate Pt for single segment events
766  *
767  */
768 void MuonSeedCreator::estimatePtSingle(SegmentContainer seg, std::vector<int> layers, double& thePt, double& theSpt) {
769 
770  thePt = defaultMomentum;
771  theSpt = defaultMomentum;
772 
773  GlobalPoint segPos = seg[0]->globalPosition();
774  double eta = segPos.eta();
775  GlobalVector gv = seg[0]->globalDirection();
776 
777  // Psi is angle between the segment origin and segment direction
778  // Use dot product between two vectors to get Psi in global x-y plane
779  double cosDpsi = (gv.x()*segPos.x() + gv.y()*segPos.y());
780  cosDpsi /= sqrt(segPos.x()*segPos.x() + segPos.y()*segPos.y());
781  cosDpsi /= sqrt(gv.x()*gv.x() + gv.y()*gv.y());
782 
783  double axb = ( segPos.x()*gv.y() ) - ( segPos.y()*gv.x() ) ;
784  double sign = (axb < 0.) ? 1.0 : -1.0;
785 
786  double dpsi = fabs(acos(cosDpsi)) ;
787  if ( dpsi > 1.570796 ) {
788  dpsi = 3.141592 - dpsi;
789  sign = -1.*sign ;
790  }
791  if (fabs(dpsi) < 0.00005) {
792  dpsi = 0.00005;
793  }
794 
795  // the 1st layer
796  if ( layers[0] == -1 ) {
797  // MB10
798  if ( fabs(eta) < 0.3 ) {
799  dpsi = scaledPhi(dpsi, SMB_10S[3] );
800  thePt = getPt( SMB10, eta , dpsi )[0];
801  theSpt = getPt( SMB10, eta , dpsi )[1];
802  }
803  // MB11
804  if ( fabs(eta) >= 0.3 && fabs(eta) < 0.82 ) {
805  dpsi = scaledPhi(dpsi, SMB_11S[3] );
806  thePt = getPt( SMB11, eta , dpsi )[0];
807  theSpt = getPt( SMB11, eta , dpsi )[1];
808  }
809  // MB12
810  if ( fabs(eta) >= 0.82 && fabs(eta) < 1.2 ) {
811  dpsi = scaledPhi(dpsi, SMB_12S[3] );
812  thePt = getPt( SMB12, eta , dpsi )[0];
813  theSpt = getPt( SMB12, eta , dpsi )[1];
814  }
815  }
816  if ( layers[0] == 1 ) {
817  // ME13
818  if ( fabs(eta) > 0.92 && fabs(eta) < 1.16 ) {
819  dpsi = scaledPhi(dpsi, SME_13S[3] );
820  thePt = getPt( SME13, eta , dpsi )[0];
821  theSpt = getPt( SME13, eta , dpsi )[1];
822  }
823  // ME12
824  if ( fabs(eta) >= 1.16 && fabs(eta) <= 1.6 ) {
825  dpsi = scaledPhi(dpsi, SME_12S[3] );
826  thePt = getPt( SME12, eta , dpsi )[0];
827  theSpt = getPt( SME12, eta , dpsi )[1];
828  }
829  }
830  if ( layers[0] == 0 ) {
831  // ME11
832  if ( fabs(eta) > 1.6 ) {
833  dpsi = scaledPhi(dpsi, SMB_11S[3] );
834  thePt = getPt( SME11, eta , dpsi )[0];
835  theSpt = getPt( SME11, eta , dpsi )[1];
836  }
837  }
838  // the 2nd layer
839  if ( layers[0] == -2 ) {
840  // MB20
841  if ( fabs(eta) < 0.25 ) {
842  dpsi = scaledPhi(dpsi, SMB_20S[3] );
843  thePt = getPt( SMB20, eta , dpsi )[0];
844  theSpt = getPt( SMB20, eta , dpsi )[1];
845  }
846  // MB21
847  if ( fabs(eta) >= 0.25 && fabs(eta) < 0.72 ) {
848  dpsi = scaledPhi(dpsi, SMB_21S[3] );
849  thePt = getPt( SMB21, eta , dpsi )[0];
850  theSpt = getPt( SMB21, eta , dpsi )[1];
851  }
852  // MB22
853  if ( fabs(eta) >= 0.72 && fabs(eta) < 1.04 ) {
854  dpsi = scaledPhi(dpsi, SMB_22S[3] );
855  thePt = getPt( SMB22, eta , dpsi )[0];
856  theSpt = getPt( SMB22, eta , dpsi )[1];
857  }
858  }
859  if ( layers[0] == 2 ) {
860  // ME22
861  if ( fabs(eta) > 0.95 && fabs(eta) <= 1.6 ) {
862  dpsi = scaledPhi(dpsi, SME_22S[3] );
863  thePt = getPt( SME22, eta , dpsi )[0];
864  theSpt = getPt( SME22, eta , dpsi )[1];
865  }
866  // ME21
867  if ( fabs(eta) > 1.6 && fabs(eta) < 2.45 ) {
868  dpsi = scaledPhi(dpsi, SME_21S[3] );
869  thePt = getPt( SME21, eta , dpsi )[0];
870  theSpt = getPt( SME21, eta , dpsi )[1];
871  }
872  }
873 
874  // the 3rd layer
875  if ( layers[0] == -3 ) {
876  // MB30
877  if ( fabs(eta) <= 0.22 ) {
878  dpsi = scaledPhi(dpsi, SMB_30S[3] );
879  thePt = getPt( SMB30, eta , dpsi )[0];
880  theSpt = getPt( SMB30, eta , dpsi )[1];
881  }
882  // MB31
883  if ( fabs(eta) > 0.22 && fabs(eta) <= 0.6 ) {
884  dpsi = scaledPhi(dpsi, SMB_31S[3] );
885  thePt = getPt( SMB31, eta , dpsi )[0];
886  theSpt = getPt( SMB31, eta , dpsi )[1];
887  }
888  // MB32
889  if ( fabs(eta) > 0.6 && fabs(eta) < 0.95 ) {
890  dpsi = scaledPhi(dpsi, SMB_32S[3] );
891  thePt = getPt( SMB32, eta , dpsi )[0];
892  theSpt = getPt( SMB32, eta , dpsi )[1];
893  }
894  }
895  thePt = fabs(thePt)*sign;
896  theSpt = fabs(theSpt);
897 
898  return;
899 }
900 
901 // setup the minimum value for obvious showering cases
902 void MuonSeedCreator::estimatePtShowering(int& NShowers, int& NShowerSegments, double& thePt, double& theSpt) {
903 
904  if ( NShowers > 2 && thePt < 300. ) {
905  thePt = 800. ;
906  theSpt = 200. ;
907  }
908  if ( NShowers == 2 && NShowerSegments > 11 && thePt < 150. ) {
909  thePt = 280. ;
910  theSpt = 70. ;
911  }
912  if ( NShowers == 2 && NShowerSegments <= 11 && thePt < 50. ) {
913  thePt = 80.;
914  theSpt = 40. ;
915  }
916  if ( NShowers == 1 && NShowerSegments <= 5 && thePt < 10. ) {
917  thePt = 16. ;
918  theSpt = 8. ;
919  }
920 
921 }
922 
923 /*
924  * weightedPt
925  *
926  * Look at delta phi between segments to determine pt as:
927  * pt = (c_1 * eta + c_2) / dphi
928  */
929 void MuonSeedCreator::weightedPt(std::vector<double> ptEstimate, std::vector<double> sptEstimate, double& thePt, double& theSpt) {
930 
931 
932  int size = ptEstimate.size();
933 
934  // If only one element, by-pass computation below
935  if (size < 2) {
936  thePt = ptEstimate[0];
937  theSpt = sptEstimate[0];
938  return;
939  }
940 
941  double charge = 0.;
942  // If have more than one pt estimator, first look if any estimator is biased
943  // by comparing the charge of each estimator
944 
945  for ( unsigned j = 0; j < ptEstimate.size(); j++ ) {
946  //std::cout<<" weighting pt: "<< ptEstimate[j] <<std::endl;
947  if ( ptEstimate[j] < 0. ) {
948  // To prevent from blowing up, add 0.1
949  charge -= 1. * (ptEstimate[j]*ptEstimate[j]) / (sptEstimate[j]*sptEstimate[j] ); // weight by relative error on pt
950  } else {
951  charge += 1. * (ptEstimate[j]*ptEstimate[j]) / (sptEstimate[j]*sptEstimate[j] ); // weight by relative error on pt
952  }
953  }
954 
955  // No need to normalize as we want to know only sign ( + or - )
956  if (charge < 0.) {
957  charge = -1.;
958  } else {
959  charge = 1.;
960  }
961 
962  //int n = 0;
963  double weightPtSum = 0.;
964  double sigmaSqr_sum = 0.;
965 
966  // Now, we want to compute average Pt using estimators with "correct" charge
967  // This is to remove biases
968  for ( unsigned j = 0; j < ptEstimate.size(); ++j ) {
969  //if ( (minpt_ratio < 0.5) && (fabs(ptEstimate[j]) < 5.0) ) continue;
970  //if ( ptEstimate[j] * charge > 0. ) {
971  //n++;
972  sigmaSqr_sum += 1.0 / (sptEstimate[j]*sptEstimate[j]);
973  weightPtSum += fabs(ptEstimate[j])/(sptEstimate[j]*sptEstimate[j]);
974  //}
975  }
976  /*
977  if (n < 1) {
978  thePt = defaultMomentum*charge;
979  theSpt = defaultMomentum;
980  return;
981  }
982  */
983  // Compute weighted mean and error
984 
985  thePt = (charge*weightPtSum) / sigmaSqr_sum;
986  theSpt = sqrt( 1.0 / sigmaSqr_sum ) ;
987 
988  //std::cout<<" final weighting : "<< thePt <<" ~ "<< fabs( theSpt/thePt ) <<std::endl;
989 
990  return;
991 }
992 
993 std::vector<double> MuonSeedCreator::getPt(std::vector<double> vPara, double eta, double dPhi ) {
994 
995  double h = fabs(eta);
996  double estPt = ( vPara[0] + vPara[1]*h + vPara[2]*h*h ) / dPhi;
997  double estSPt = ( vPara[3] + vPara[4]*h + vPara[5]*h*h ) * estPt;
998  std::vector<double> paraPt ;
999  paraPt.push_back( estPt );
1000  paraPt.push_back( estSPt ) ;
1001 
1002  //std::cout<<" pt:"<<estPt<<" +/-"<< estSPt<<" h:"<<eta<<" df:"<<dPhi<<std::endl;
1003  return paraPt ;
1004 }
1005 
1006 double MuonSeedCreator::scaledPhi( double dphi, double t1) {
1007 
1008  if (dphi != 0. ) {
1009 
1010  double oPhi = 1./dphi ;
1011  dphi = dphi /( 1. + t1/( oPhi + 10. ) ) ;
1012  return dphi ;
1013 
1014  } else {
1015  return dphi ;
1016  }
1017 
1018 }
1019 
std::vector< double > OL2213
type
Definition: HCALResponse.h:22
T getParameter(std::string const &) const
std::vector< double > DT13_2
int i
Definition: DBlmapReader.cc:9
TrajectorySeed createSeed(int type, SegmentContainer seg, std::vector< int > layers, int NShower, int NShowerSeg)
Create a seed from set of segments.
std::vector< double > SMB32
void estimatePtCSC(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from CSC measurements.
std::vector< double > CSC23_1
std::vector< double > SMB_30S
std::vector< double > SMB_21S
std::vector< double > DT14_2
std::vector< double > DT23_2
std::vector< double > DT34_1
T perp() const
Definition: PV3DBase.h:66
std::vector< double > DT14_1
std::vector< double > SME_22S
std::vector< double > SMB22
std::vector< double > OL1232
std::vector< double > DT34_2
std::vector< double > CSC24
std::vector< double > CSC23_2
std::vector< double > SMB20
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
std::vector< double > getPt(std::vector< double > vParameters, double eta, double dPhi)
Compute pt from parameters.
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< double > CSC14
void estimatePtSingle(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from single segment.
void estimatePtOverlap(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from CSC + DT measurements.
T y() const
Definition: PV3DBase.h:57
std::vector< double > DT34
std::vector< double > CSC12_2
std::vector< double > CSC13_2
std::vector< double > CSC12_1
std::vector< double > SME12
std::vector< double > SME22
void weightedPt(std::vector< double > ptEstimate, std::vector< double > sptEstimate, double &ptAvg, double &sptAvg)
Compute weighted mean pt from different pt estimators.
MuonTransientTrackingRecHit::MuonRecHitContainer SegmentContainer
std::vector< double > SME21
double charge(const std::vector< uint8_t > &Ampls)
T eta() const
std::vector< double > CSC13_3
MuonSeedCreator(const edm::ParameterSet &pset)
Constructor.
const MagneticField * BField
std::vector< double > CSC34
std::vector< double > SME13
Geom::Theta< T > theta() const
Definition: PV3DBase.h:69
std::vector< double > DT13
std::vector< double > CSC12_3
void push_back(D *&d)
Definition: OwnVector.h:288
std::vector< double > SMB10
std::vector< double > CSC12
std::vector< double > DT12_2
std::vector< double > OL_1222
std::vector< double > OL1222
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
void estimatePtShowering(int &NShowers, int &NShowerSeg, double &pt, double &spt)
Estimate transverse momentum of track from showering segment.
T sqrt(T t)
Definition: SSEVec.h:28
std::vector< double > CSC23
std::vector< double > CSC13
std::vector< double > DT23_1
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
std::vector< double > SMB30
int j
Definition: DBlmapReader.cc:9
std::vector< double > SMB_10S
tuple pset
Definition: CrabTask.py:85
PTrajectoryStateOnDet * persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid) const
std::vector< double > DT24_2
std::vector< double > SMB_31S
std::vector< double > SMB_32S
std::vector< double > DT12_1
std::vector< double > CSC14_3
std::vector< double > DT12
void estimatePtDT(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from DT measurements.
std::vector< double > SMB_20S
std::vector< double > SMB_11S
std::vector< double > OL1213
Definition: DetId.h:20
CLHEP::HepVector AlgebraicVector
std::vector< double > SMB12
std::vector< double > SME_21S
T * clone(const T *tp)
Definition: Ptr.h:42
std::vector< double > SME_12S
std::vector< double > OL2222
std::vector< double > SME_11S
std::vector< double > DT14
T eta() const
Definition: PV3DBase.h:70
std::vector< double > DT13_1
std::vector< double > SME31
std::vector< double > DT24
std::vector< double > OL_2213
std::vector< double > DT23
std::vector< double > SME32
std::vector< double > CSC24_1
std::vector< double > SMB21
CLHEP::HepSymMatrix AlgebraicSymMatrix
std::vector< double > CSC01_1
tuple cout
Definition: gather_cfg.py:41
double scaledPhi(double dphi, double t1)
Scale the dPhi from segment position.
~MuonSeedCreator()
Destructor.
std::vector< double > SME_13S
std::vector< double > OL_2222
std::vector< double > SMB31
std::vector< double > CSC03
std::vector< double > CSC02
std::vector< double > DT24_1
std::vector< double > CSC34_1
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
std::vector< double > SME11
T x() const
Definition: PV3DBase.h:56
std::vector< double > SMB11
std::vector< double > OL_1232
std::vector< double > SMB_12S
std::vector< double > SME41
tuple size
Write out results.
Definition: QTest.h:564
std::vector< double > OL_1213
std::vector< double > SMB_22S