CMS 3D CMS Logo
MuonSeedPtExtractor.cc
Go to the documentation of this file.
1 #include "RecoMuon/MuonSeedGenerator/src/MuonSeedPtExtractor.h"
2 #include "FWCore/MessageLogger/interface/MessageLogger.h"
3 #include "FWCore/ParameterSet/interface/ParameterSet.h"
4 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
5 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
6 
7 #include "TMath.h"
8 #include <sstream>
9 
10 MuonSeedPtExtractor::MuonSeedPtExtractor(const edm::ParameterSet& par)
11  : theBeamSpot(0., 0., 0.), scaleDT_(par.getParameter<bool>("scaleDT")) {
12  init(par);
13 }
14 
15 void MuonSeedPtExtractor::init(const edm::ParameterSet& par) {
16  // load pT seed parameters
17  // DT combinations
18  fillParametersForCombo("DT_12", par);
19  fillParametersForCombo("DT_13", par);
20  fillParametersForCombo("DT_14", par);
21  fillParametersForCombo("DT_23", par);
22  fillParametersForCombo("DT_24", par);
23  fillParametersForCombo("DT_34", par);
24  // CSC combinations
25  fillParametersForCombo("CSC_01", par);
26  fillParametersForCombo("CSC_12", par);
27  fillParametersForCombo("CSC_02", par);
28  fillParametersForCombo("CSC_13", par);
29  fillParametersForCombo("CSC_03", par);
30  fillParametersForCombo("CSC_14", par);
31  fillParametersForCombo("CSC_23", par);
32  fillParametersForCombo("CSC_24", par);
33  fillParametersForCombo("CSC_34", par);
34 
35  // Overlap combinations
36  fillParametersForCombo("OL_1213", par);
37  fillParametersForCombo("OL_1222", par);
38  fillParametersForCombo("OL_1232", par);
39  fillParametersForCombo("OL_2213", par);
40  fillParametersForCombo("OL_2222", par);
41 
42  // Single segments (CSC)
43  fillParametersForCombo("SME_11", par);
44  fillParametersForCombo("SME_12", par);
45  fillParametersForCombo("SME_13", par);
46  fillParametersForCombo("SME_21", par);
47  fillParametersForCombo("SME_22", par);
48  fillParametersForCombo("SME_31", par);
49  fillParametersForCombo("SME_32", par);
50  fillParametersForCombo("SME_41", par);
51  fillParametersForCombo("SME_42", par);
52 
53  // Single segments (DT)
54  fillParametersForCombo("SMB_10", par);
55  fillParametersForCombo("SMB_11", par);
56  fillParametersForCombo("SMB_12", par);
57  fillParametersForCombo("SMB_20", par);
58  fillParametersForCombo("SMB_21", par);
59  fillParametersForCombo("SMB_22", par);
60  fillParametersForCombo("SMB_30", par);
61  fillParametersForCombo("SMB_31", par);
62  fillParametersForCombo("SMB_32", par);
63 
64  fillScalesForCombo("CSC_01_1_scale", par);
65  fillScalesForCombo("CSC_12_1_scale", par);
66  fillScalesForCombo("CSC_12_2_scale", par);
67  fillScalesForCombo("CSC_12_3_scale", par);
68  fillScalesForCombo("CSC_13_2_scale", par);
69  fillScalesForCombo("CSC_13_3_scale", par);
70  fillScalesForCombo("CSC_14_3_scale", par);
71  fillScalesForCombo("CSC_23_1_scale", par);
72  fillScalesForCombo("CSC_23_2_scale", par);
73  fillScalesForCombo("CSC_24_1_scale", par);
74  fillScalesForCombo("CSC_34_1_scale", par);
75  fillScalesForCombo("DT_12_1_scale", par);
76  fillScalesForCombo("DT_12_2_scale", par);
77  fillScalesForCombo("DT_13_1_scale", par);
78  fillScalesForCombo("DT_13_2_scale", par);
79  fillScalesForCombo("DT_14_1_scale", par);
80  fillScalesForCombo("DT_14_2_scale", par);
81  fillScalesForCombo("DT_23_1_scale", par);
82  fillScalesForCombo("DT_23_2_scale", par);
83  fillScalesForCombo("DT_24_1_scale", par);
84  fillScalesForCombo("DT_24_2_scale", par);
85  fillScalesForCombo("DT_34_1_scale", par);
86  fillScalesForCombo("DT_34_2_scale", par);
87  fillScalesForCombo("OL_1213_0_scale", par);
88  fillScalesForCombo("OL_1222_0_scale", par);
89  fillScalesForCombo("OL_1232_0_scale", par);
90  fillScalesForCombo("OL_2213_0_scale", par);
91  fillScalesForCombo("OL_2222_0_scale", par);
92  fillScalesForCombo("SMB_10_0_scale", par);
93  fillScalesForCombo("SMB_11_0_scale", par);
94  fillScalesForCombo("SMB_12_0_scale", par);
95  fillScalesForCombo("SMB_20_0_scale", par);
96  fillScalesForCombo("SMB_21_0_scale", par);
97  fillScalesForCombo("SMB_22_0_scale", par);
98  fillScalesForCombo("SMB_30_0_scale", par);
99  fillScalesForCombo("SMB_31_0_scale", par);
100  fillScalesForCombo("SMB_32_0_scale", par);
101  fillScalesForCombo("SME_11_0_scale", par);
102  fillScalesForCombo("SME_12_0_scale", par);
103  fillScalesForCombo("SME_13_0_scale", par);
104  fillScalesForCombo("SME_21_0_scale", par);
105  fillScalesForCombo("SME_22_0_scale", par);
106 }
107 
108 MuonSeedPtExtractor::~MuonSeedPtExtractor() {}
109 
110 void MuonSeedPtExtractor::fillParametersForCombo(const std::string& name, const edm::ParameterSet& pset) {
111  theParametersForCombo[name] = pset.getParameter<std::vector<double> >(name);
112 }
113 
114 void MuonSeedPtExtractor::fillScalesForCombo(const std::string& name, const edm::ParameterSet& pset) {
115  theScalesForCombo[name] = pset.getParameter<std::vector<double> >(name);
116 }
117 
118 std::vector<double> MuonSeedPtExtractor::pT_extract(
119  MuonTransientTrackingRecHit::ConstMuonRecHitPointer firstHit,
120  MuonTransientTrackingRecHit::ConstMuonRecHitPointer secondHit) const {
121  GlobalPoint innerPoint = firstHit->globalPosition() - theBeamSpot;
122  GlobalPoint outerPoint = secondHit->globalPosition() - theBeamSpot;
123  MuonTransientTrackingRecHit::ConstMuonRecHitPointer innerHit = firstHit;
124  MuonTransientTrackingRecHit::ConstMuonRecHitPointer outerHit = secondHit;
125 
126  // ways in which the hits could be in the wrong order
127  if ((outerHit->isDT() && innerHit->isCSC()) ||
128  (outerHit->isDT() && innerHit->isDT() && (innerPoint.perp() > outerPoint.perp())) ||
129  (outerHit->isCSC() && innerHit->isCSC() && (fabs(innerPoint.z()) > fabs(outerPoint.z())))) {
130  innerHit = secondHit;
131  outerHit = firstHit;
132  innerPoint = innerHit->globalPosition() - theBeamSpot;
133  outerPoint = outerHit->globalPosition() - theBeamSpot;
134  }
135 
136  double phiInner = innerPoint.phi();
137  double phiOuter = outerPoint.phi();
138 
139  double etaInner = innerPoint.eta();
140  double etaOuter = outerPoint.eta();
141  //std::cout<<" inner pos = "<< innerPoint << " phi eta " << phiInner << " " << etaInner << std::endl;
142  //std::cout<<" outer pos = "<< outerPoint << " phi eta " << phiOuter << " " << etaOuter << std::endl;
143  //double thetaInner = firstHit->globalPosition().theta();
144  // if some of the segments is missing r-phi measurement then we should
145  // use only the 4D phi estimate (also use 4D eta estimate only)
146  // the direction is not so important (it will be corrected)
147  /*
148  bool firstOK = (4==allValidSets[iSet][firstMeasurement]->hit->dimension());
149  bool lastOK = (4==allValidSets[iSet][lastMeasurement]->hit->dimension());
150  if(!(firstOK * lastOK)){
151  if(!firstOK){
152  }
153  if(!firstOK){
154  }
155  }
156  */
157  double dPhi = phiInner - phiOuter;
158  if (dPhi < -TMath::Pi()) {
159  dPhi += 2 * TMath::Pi();
160  } else if (dPhi > TMath::Pi()) {
161  dPhi -= 2 * TMath::Pi();
162  }
163  int sign = 1;
164  if (dPhi < 0.) {
165  dPhi = -dPhi;
166  sign = -1;
167  }
168 
169  if (dPhi < 1.0e-6) {
170  dPhi = 1.0e-6;
171  }
172  double eta = fabs(etaOuter); // what if it is 2D segment? use inner?
173 
174  std::vector<int> stationCoded(2);
175  stationCoded[0] = stationCoded[1] = 999;
176 
177  DetId detId_inner = innerHit->hit()->geographicalId();
178  DetId detId_outer = outerHit->hit()->geographicalId();
179 
180  stationCoded[0] = stationCode(innerHit);
181  stationCoded[1] = stationCode(outerHit);
182 
183  std::ostringstream os0;
184  std::ostringstream os1;
185  os0 << abs(stationCoded[0]);
186  os1 << abs(stationCoded[1]);
187 
188  //std::cout<<" st1 = "<<stationCoded[0]<<" st2 = "<<stationCoded[1]<<std::endl;
189  //std::cout<<" detId_inner = "<<detId_inner.rawId()<<" detId_outer = "<< detId_outer.rawId()<<std::endl;
190  std::string combination = "0";
191  std::string init_combination = combination;
192  bool singleSegment = false;
193  //if(detId_first == detId_second){
194  if (stationCoded[0] == stationCoded[1]) {
195  // single segment - DT or CSC
196  singleSegment = true;
197  //eta = innerPoint.eta();
198  GlobalVector gv = innerHit->globalDirection();
199  double gvPerpPos = gv.x() * gv.x() + gv.y() * gv.y();
200  if (gvPerpPos < 1e-32)
201  gvPerpPos = 1e-32;
202  gvPerpPos = sqrt(gvPerpPos);
203  // Psi is angle between the segment origin and segment direction
204  // Use dot product between two vectors to get Psi in global x-y plane
205  double cosDpsi = (gv.x() * innerPoint.x() + gv.y() * innerPoint.y());
206  if (cosDpsi != 0) {
207  cosDpsi /= sqrt(innerPoint.x() * innerPoint.x() + innerPoint.y() * innerPoint.y());
208  cosDpsi /= gvPerpPos;
209  cosDpsi = cosDpsi > 1 ? 1 : cosDpsi;
210  cosDpsi = cosDpsi < -1 ? -1 : cosDpsi;
211  }
212 
213  double axb = (innerPoint.x() * gv.y()) - (innerPoint.y() * gv.x());
214  sign = (axb < 0.) ? 1 : -1;
215 
216  double dpsi = fabs(acos(cosDpsi));
217  if (dpsi > TMath::Pi() / 2.) {
218  dpsi = TMath::Pi() - dpsi;
219  }
220  if (fabs(dpsi) < 0.00005) {
221  dpsi = 0.00005;
222  }
223  dPhi = dpsi;
224 
225  if (innerHit->isDT()) {
226  DTChamberId dtCh(detId_inner);
227  std::ostringstream os;
228  os << dtCh.station() << abs(dtCh.wheel());
229  combination = "SMB_" + os.str();
230  } else if (innerHit->isCSC()) {
231  CSCDetId cscId(detId_inner);
232  std::ostringstream os;
233  int ring = cscId.ring();
234  if (ring == 4)
235  ring = 1;
236  os << cscId.station() << ring;
237  combination = "SME_" + os.str();
238  } else if (innerHit->isME0()) {
239  ME0DetId me0Id(detId_inner);
240  std::ostringstream os;
241  int ring = 1; //me0Id.ring(); me0 only in ring 1
242  os << me0Id.station() << ring;
243  combination = "SME_" + os.str();
244  } else {
245  throw cms::Exception("MuonSeedPtExtractor") << "Bad hit DetId";
246  }
247  } else {
248  if (stationCoded[0] < 0) {
249  if (stationCoded[1] < 0) {
250  // DT-DT
251  combination = "DT_" + os0.str() + os1.str();
252  } else {
253  // DT-CSC
254  eta = fabs(etaInner);
255  if (-1 == stationCoded[0]) {
256  switch (stationCoded[1]) {
257  case 1:
258  combination = "OL_1213";
259  break;
260  case 2:
261  combination = "OL_1222";
262  break;
263  case 3:
264  combination = "OL_1232";
265  break;
266  default:
267  // can not be
268  break;
269  }
270  } else if (-2 == stationCoded[0]) {
271  if (1 == stationCoded[1]) {
272  combination = "OL_2213";
273  } else {
274  // can not be (not coded?)
275  combination = "OL_2222"; // in case
276  }
277  } else {
278  // can not be
279  }
280  }
281  } else {
282  if (stationCoded[1] < 0) {
283  // CSC-DT
284  // can not be
285  } else {
286  // CSC-CSC
287  combination = "CSC_" + os0.str() + os1.str();
288  if ("CSC_04" == combination) {
289  combination = "CSC_14";
290  }
291  }
292  }
293  }
294 
295  std::vector<double> pTestimate(2); //
296  //std::cout<<" combination = "<<combination<<std::endl;
297  if (init_combination != combination) {
298  //std::cout<<" combination = "<<combination<<" eta = "<<eta<<" dPhi = "<<dPhi<<std::endl;
299  ParametersMap::const_iterator parametersItr = theParametersForCombo.find(combination);
300  if (parametersItr == theParametersForCombo.end()) {
301  // edm::LogWarning("RecoMuon|MuonSeedGenerator|MuonSeedPtExtractor") << "Cannot find parameters for combo " << combination;
302  edm::LogWarning("BadSegmentCombination") << "Cannot find parameters for combo " << combination;
303  pTestimate[0] = pTestimate[1] = 100;
304  // throw cms::Exception("MuonSeedPtEstimator") << "Cannot find parameters for combo " << combination;
305  } else {
306  if (scaleDT_ && outerHit->isDT()) {
307  pTestimate = getPt(parametersItr->second, eta, dPhi, combination, detId_outer);
308  } else {
309  pTestimate = getPt(parametersItr->second, eta, dPhi);
310  }
311 
312  if (singleSegment) {
313  pTestimate[0] = fabs(pTestimate[0]);
314  pTestimate[1] = fabs(pTestimate[1]);
315  }
316  pTestimate[0] *= double(sign);
317  }
318  } else {
319  // often a MB3 - ME1/3 seed
320  pTestimate[0] = pTestimate[1] = 100;
321  // hmm
322  }
323  // std::cout<<" pTestimate[0] = "<<pTestimate[0]<<" pTestimate[1] = "<<pTestimate[1]<<std::endl;
324  /*
325  MuonRecHitContainer_clusters[cluster][iHit+1]->isDT());
326  if(specialCase){
327  DTChamberId dtCh(detId);
328  DTChamberId dtCh_2(detId_2);
329  specialCase = (dtCh.station() == dtCh_2.station());
330  }
331  */
332  //return vPara;
333  return pTestimate;
334 }
335 
336 int MuonSeedPtExtractor::stationCode(MuonTransientTrackingRecHit::ConstMuonRecHitPointer hit) const {
337  DetId detId(hit->hit()->geographicalId());
338  int result = -999;
339  if (hit->isDT()) {
340  DTChamberId dtCh(detId);
341  //std::cout<<"first (DT) St/W/S = "<<dtCh.station()<<"/"<<dtCh.wheel()<<"/"<<dtCh.sector()<<"/"<<std::endl;
342  result = -1 * dtCh.station();
343  } else if (hit->isCSC()) {
344  CSCDetId cscID(detId);
345  //std::cout<<"first (CSC) E/S/R/C = "<<cscID.endcap()<<"/"<<cscID.station()<<"/"<<cscID.ring()<<"/"<<cscID.chamber()<<std::endl;
346  result = cscID.station();
347  if (result == 1 && (1 == cscID.ring() || 4 == cscID.ring()))
348  result = 0;
349  } else if (hit->isME0()) {
350  result = 0;
351  } else if (hit->isRPC()) {
352  }
353  return result;
354 }
355 
356 std::vector<double> MuonSeedPtExtractor::getPt(const std::vector<double>& vPara, double eta, double dPhi) const {
357  //std::cout<<" eta = "<<eta<<" dPhi = "<<dPhi<<" vPara[0] = "<<vPara[0]<<" vPara[1] = "<<vPara[1]<<" vPara[2] = "<<vPara[2]<<std::endl;
358  double h = fabs(eta);
359  double estPt = (vPara[0] + vPara[1] * h + vPara[2] * h * h) / dPhi;
360  double estSPt = (vPara[3] + vPara[4] * h + vPara[5] * h * h) / dPhi;
361  // std::cout<<"estPt = "<<estPt<<std::endl;
362  std::vector<double> paraPt;
363  paraPt.push_back(estPt);
364  paraPt.push_back(estSPt);
365  return paraPt;
366 }
367 
368 // combined pt parameterization and pt scale
369 std::vector<double> MuonSeedPtExtractor::getPt(const std::vector<double>& vPara,
370  double eta,
371  double dPhi,
372  const std::string& combination,
373  const DTChamberId& outerDetId) const {
374  double h = fabs(eta);
375  double estPt = (vPara[0] + vPara[1] * h + vPara[2] * h * h) / dPhi;
376  // changed by S.C.
377  double estSPt = (vPara[3] + vPara[4] * h + vPara[5] * h * h) / dPhi;
378 
379  // scale the pt and spt , changed by S.C.
380  int wheel = 0;
381  if (combination[0] == 'D') {
382  wheel = abs(outerDetId.wheel());
383  }
384 
385  std::ostringstream os;
386  os << combination << "_" << wheel << "_scale";
387 
388  ScalesMap::const_iterator scalesItr = theScalesForCombo.find(os.str());
389  if (scalesItr != theScalesForCombo.end()) {
390  double t1 = scalesItr->second[0];
391  double scaleFactor = 1. / (1. + t1 / (estPt + 10.));
392 
393  estSPt = estSPt * scaleFactor;
394  estPt = estPt * scaleFactor;
395  }
396 
397  // std::cout<<"estPt = "<<estPt<<std::endl;
398  std::vector<double> paraPt;
399  paraPt.push_back(estPt);
400  paraPt.push_back(estSPt);
401  return paraPt;
402 }
Vector3DBase
Definition: Vector3DBase.h:8
electrons_cff.bool
bool
Definition: electrons_cff.py:393
MuonTransientTrackingRecHit::ConstMuonRecHitPointer
std::shared_ptr< MuonTransientTrackingRecHit const > ConstMuonRecHitPointer
Definition: MuonTransientTrackingRecHit.h:22
MessageLogger.h
MuonSeedPtExtractor::~MuonSeedPtExtractor
virtual ~MuonSeedPtExtractor()
Destructor.
Definition: MuonSeedPtExtractor.cc:108
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
CSCDetId::ring
int ring() const
Definition: CSCDetId.h:68
h
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
Definition: L1TUtmAlgorithmRcd.h:4
ME0DetId.h
Validation_hcalonly_cfi.sign
sign
Definition: Validation_hcalonly_cfi.py:32
MuonSeedPtExtractor::scaleDT_
bool scaleDT_
Definition: MuonSeedPtExtractor.h:51
MuonSeedPtExtractor::theScalesForCombo
ScalesMap theScalesForCombo
Definition: MuonSeedPtExtractor.h:49
HLT_FULL_cff.dPhi
dPhi
Definition: HLT_FULL_cff.py:13702
MuonSeedPtExtractor::stationCode
int stationCode(MuonTransientTrackingRecHit::ConstMuonRecHitPointer hit) const
Definition: MuonSeedPtExtractor.cc:336
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
CSCDetId.h
MuonSeedPtExtractor::init
void init(const edm::ParameterSet &par)
Definition: MuonSeedPtExtractor.cc:15
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
DetId
Definition: DetId.h:17
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
MuonSeedPtExtractor::theParametersForCombo
ParametersMap theParametersForCombo
Definition: MuonSeedPtExtractor.h:48
PVValHelper::eta
Definition: PVValidationHelpers.h:69
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
MuonSeedPtExtractor::theBeamSpot
GlobalVector theBeamSpot
Definition: MuonSeedPtExtractor.h:50
MuonSeedPtExtractor::getPt
std::vector< double > getPt(const std::vector< double > &vPara, double eta, double dPhi) const
Definition: MuonSeedPtExtractor.cc:356
h
Point3DBase< float, GlobalTag >
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
MuonSeedPtExtractor::MuonSeedPtExtractor
MuonSeedPtExtractor(const edm::ParameterSet &)
Constructor with Parameter set and MuonServiceProxy.
Definition: MuonSeedPtExtractor.cc:10
edm::ParameterSet
Definition: ParameterSet.h:47
CSCDetId
Definition: CSCDetId.h:26
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:73
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
makeMuonMisalignmentScenario.wheel
wheel
Definition: makeMuonMisalignmentScenario.py:319
ME0DetId::station
int station() const
Definition: ME0DetId.h:58
ME0DetId
Definition: ME0DetId.h:16
MuonSeedPtExtractor::pT_extract
virtual std::vector< double > pT_extract(MuonTransientTrackingRecHit::ConstMuonRecHitPointer firstHit, MuonTransientTrackingRecHit::ConstMuonRecHitPointer secondHit) const
Definition: MuonSeedPtExtractor.cc:118
MuonSeedPtExtractor.h
MuonSeedPtExtractor::fillScalesForCombo
void fillScalesForCombo(const std::string &name, const edm::ParameterSet &pset)
Definition: MuonSeedPtExtractor.cc:114
Exception
Definition: hltDiff.cc:246
multiplicitycorr_cfi.scaleFactor
scaleFactor
Definition: multiplicitycorr_cfi.py:7
relativeConstraints.ring
ring
Definition: relativeConstraints.py:68
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
MuonSeedPtExtractor::fillParametersForCombo
void fillParametersForCombo(const std::string &name, const edm::ParameterSet &pset)
Definition: MuonSeedPtExtractor.cc:110
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
mps_fire.result
result
Definition: mps_fire.py:311
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
DTChamberId
Definition: DTChamberId.h:14
ParameterSet.h
CSCDetId::station
int station() const
Definition: CSCDetId.h:79
DTChamberId::wheel
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:39
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
DTChamberId::station
int station() const
Return the station number.
Definition: DTChamberId.h:42
hit
Definition: SiStripHitEffFromCalibTree.cc:88
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16