File indexing completed on 2025-01-09 23:33:55
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "RecoMuon/MuonSeedGenerator/interface/MuonDTSeedFromRecHits.h"
0010 #include "RecoMuon/MuonSeedGenerator/interface/MuonSeedPtExtractor.h"
0011 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0012
0013 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0014
0015 #include "DataFormats/Common/interface/OwnVector.h"
0016 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0017 #include "DataFormats/Math/interface/deltaPhi.h"
0018
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020
0021 #include "gsl/gsl_statistics.h"
0022
0023 using namespace std;
0024
0025 MuonDTSeedFromRecHits::MuonDTSeedFromRecHits() : MuonSeedFromRecHits() {}
0026
0027 TrajectorySeed MuonDTSeedFromRecHits::seed() const {
0028 double pt[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0029
0030
0031 double spt[8] = {1 / 0.048, 1 / 0.075, 1 / 0.226, 1 / 0.132, 1 / 0.106, 1 / 0.175, 1 / 0.125, 1 / 0.185};
0032
0033 const std::string metname = "Muon|RecoMuon|MuonDTSeedFromRecHits";
0034
0035
0036 computePtWithVtx(pt, spt);
0037
0038
0039 computePtWithoutVtx(pt, spt);
0040
0041
0042 LogTrace(metname) << " Pt MB1 @vtx: " << pt[0] << " w: " << spt[0] << "\n"
0043 << " Pt MB2 @vtx: " << pt[1] << " w: " << spt[1] << endl;
0044
0045 LogTrace(metname) << " Pt MB2-MB1 " << pt[2] << " w: " << spt[2] << "\n"
0046 << " Pt MB3-MB1 " << pt[3] << " w: " << spt[3] << "\n"
0047 << " Pt MB3-MB2 " << pt[4] << " w: " << spt[4] << "\n"
0048 << " Pt MB4-MB1 " << pt[5] << " w: " << spt[5] << "\n"
0049 << " Pt MB4-MB2 " << pt[6] << " w: " << spt[6] << "\n"
0050 << " Pt MB4-MB3 " << pt[7] << " w: " << spt[7] << endl;
0051
0052
0053 float ptmean = 0.;
0054 float sptmean = 0.;
0055
0056 computeBestPt(pt, spt, ptmean, sptmean);
0057
0058
0059 int npoints = 0;
0060 for (int i = 0; i < 8; ++i)
0061 if (pt[i] != 0)
0062 ++npoints;
0063 if (npoints != 0) {
0064 sptmean = sqrt(sptmean * sptmean + 0.09 * ptmean * ptmean / npoints);
0065 }
0066
0067 LogTrace(metname) << " Seed Pt: " << ptmean << " +/- " << sptmean << endl;
0068
0069 ConstMuonRecHitPointer last = bestBarrelHit(theRhits);
0070 return createSeed(ptmean, sptmean, last);
0071 }
0072
0073 MuonDTSeedFromRecHits::ConstMuonRecHitPointer MuonDTSeedFromRecHits::bestBarrelHit(
0074 const MuonRecHitContainer& barrelHits) const {
0075 int alt_npt = 0;
0076 int best_npt = 0;
0077 int cur_npt = 0;
0078 MuonRecHitPointer best = nullptr;
0079 MuonRecHitPointer alter = nullptr;
0080
0081 for (MuonRecHitContainer::const_iterator iter = barrelHits.begin(); iter != barrelHits.end(); iter++) {
0082 bool hasZed = ((*iter)->projectionMatrix()[1][1] == 1);
0083
0084 cur_npt = 0;
0085 vector<TrackingRecHit*> slrhs = (*iter)->recHits();
0086 for (vector<TrackingRecHit*>::const_iterator slrh = slrhs.begin(); slrh != slrhs.end(); ++slrh) {
0087 cur_npt += (*slrh)->recHits().size();
0088 }
0089 float radius1 = (*iter)->det()->position().perp();
0090
0091 if (hasZed) {
0092 if (cur_npt > best_npt) {
0093 best = (*iter);
0094 best_npt = cur_npt;
0095 } else if (best && cur_npt == best_npt) {
0096 float radius2 = best->det()->position().perp();
0097 if (radius1 < radius2) {
0098 best = (*iter);
0099 best_npt = cur_npt;
0100 }
0101 }
0102 }
0103
0104 if (cur_npt > alt_npt) {
0105 alter = (*iter);
0106 alt_npt = cur_npt;
0107 } else if (alter && cur_npt == alt_npt) {
0108 float radius2 = alter->det()->position().perp();
0109 if (radius1 < radius2) {
0110 alter = (*iter);
0111 alt_npt = cur_npt;
0112 }
0113 }
0114 }
0115
0116 if (!best)
0117 best = alter;
0118
0119 return best;
0120 }
0121
0122 float MuonDTSeedFromRecHits::bestEta() const {
0123 int Maxseg = 0;
0124 float Msdeta = 100.;
0125 float result = (*theRhits.begin())->globalPosition().eta();
0126
0127 for (MuonRecHitContainer::const_iterator iter = theRhits.begin(); iter != theRhits.end(); iter++) {
0128 float eta1 = (*iter)->globalPosition().eta();
0129
0130 int Nseg = 0;
0131 float sdeta = .0;
0132
0133 for (MuonRecHitContainer::const_iterator iter2 = theRhits.begin(); iter2 != theRhits.end(); iter2++) {
0134 if (iter2 == iter)
0135 continue;
0136
0137 float eta2 = (*iter2)->globalPosition().eta();
0138
0139 if (fabs(eta1 - eta2) > .2)
0140 continue;
0141
0142 Nseg++;
0143 sdeta += fabs(eta1 - eta2);
0144
0145 if (Nseg > Maxseg || (Nseg == Maxseg && sdeta < Msdeta)) {
0146 Maxseg = Nseg;
0147 Msdeta = sdeta;
0148 result = eta1;
0149 }
0150 }
0151 }
0152 return result;
0153 }
0154
0155 void MuonDTSeedFromRecHits::computePtWithVtx(double* pt, double* spt) const {
0156 float eta0 = bestEta();
0157
0158 for (MuonRecHitContainer::const_iterator iter = theRhits.begin(); iter != theRhits.end(); iter++) {
0159
0160
0161 float eta1 = (*iter)->globalPosition().eta();
0162 if (fabs(eta1 - eta0) > .2)
0163 continue;
0164
0165
0166 unsigned int stat = DTChamberId((**iter).geographicalId()).station();
0167 if (stat > 2)
0168 continue;
0169
0170 double thispt = thePtExtractor->pT_extract(*iter, *iter)[0];
0171 float ptmax = 2000.;
0172 if (thispt > ptmax)
0173 thispt = ptmax;
0174 if (thispt < -ptmax)
0175 thispt = -ptmax;
0176
0177 if (stat == 1) {
0178 pt[0] = thispt;
0179 }
0180
0181 if (stat == 2) {
0182 pt[1] = thispt;
0183 }
0184 }
0185 }
0186
0187 void MuonDTSeedFromRecHits::computePtWithoutVtx(double* pt, double* spt) const {
0188 float eta0 = bestEta();
0189
0190 for (MuonRecHitContainer::const_iterator iter = theRhits.begin(); iter != theRhits.end(); iter++) {
0191
0192 float eta1 = (*iter)->globalPosition().eta();
0193 if (fabs(eta1 - eta0) > .2)
0194 continue;
0195
0196 for (MuonRecHitContainer::const_iterator iter2 = theRhits.begin(); iter2 != iter; iter2++) {
0197
0198 float eta2 = (*iter2)->globalPosition().eta();
0199 if (fabs(eta2 - eta0) > .2)
0200 continue;
0201
0202 unsigned int stat1 = DTChamberId((*iter)->geographicalId()).station();
0203 unsigned int stat2 = DTChamberId((*iter2)->geographicalId()).station();
0204
0205 if (stat1 > stat2) {
0206 int tmp = stat1;
0207 stat1 = stat2;
0208 stat2 = tmp;
0209 }
0210 unsigned int st = stat1 * 10 + stat2;
0211 double thispt = thePtExtractor->pT_extract(*iter, *iter2)[0];
0212 switch (st) {
0213 case 12: {
0214 pt[2] = thispt;
0215 break;
0216 }
0217 case 13: {
0218 pt[3] = thispt;
0219 break;
0220 }
0221 case 14: {
0222 pt[5] = thispt;
0223 break;
0224 }
0225 case 23: {
0226 pt[4] = thispt;
0227 break;
0228 }
0229 case 24: {
0230 pt[6] = thispt;
0231 break;
0232 }
0233 case 34: {
0234 pt[7] = thispt;
0235 break;
0236 }
0237 default: {
0238 break;
0239 }
0240 }
0241 }
0242 }
0243 }
0244
0245 void MuonDTSeedFromRecHits::computeBestPt(double* pt, double* spt, float& ptmean, float& sptmean) const {
0246 const std::string metname = "Muon|RecoMuon|MuonDTSeedFromRecHits";
0247
0248 int nTotal = 8;
0249 int igood = -1;
0250 for (int i = 0; i <= 7; i++) {
0251 if (pt[i] == 0) {
0252 spt[i] = 0.;
0253 nTotal--;
0254 } else {
0255 igood = i;
0256 }
0257 }
0258
0259 if (nTotal == 1) {
0260
0261 ptmean = pt[igood];
0262 sptmean = 1 / sqrt(spt[igood]);
0263 } else if (nTotal == 0) {
0264
0265 ptmean = 50;
0266 sptmean = 30;
0267 } else {
0268
0269
0270 float ptvtx = 0.0;
0271 float sptvtx = 0.0;
0272 computeMean(pt, spt, 2, false, ptvtx, sptvtx);
0273 LogTrace(metname) << " GSL: Pt w vtx :" << ptvtx << "+/-" << sptvtx << endl;
0274
0275
0276 if (ptvtx != 0.) {
0277 ptmean = ptvtx;
0278 sptmean = sptvtx;
0279 return;
0280
0281 }
0282
0283 float ptMB = 0.0;
0284 float sptMB = 0.0;
0285 computeMean(pt + 2, spt + 2, 6, false, ptMB, sptMB);
0286 LogTrace(metname) << " GSL: Pt w/o vtx :" << ptMB << "+/-" << sptMB << endl;
0287
0288
0289 float ptpool = 0.0;
0290 if ((ptvtx + ptMB) != 0.0)
0291 ptpool = (ptvtx - ptMB) / (ptvtx + ptMB);
0292 bool fromvtx = true;
0293 if (fabs(ptpool) > 0.2)
0294 fromvtx = false;
0295 LogTrace(metname) << "From vtx? " << fromvtx << " ptpool " << ptpool << endl;
0296
0297
0298 computeMean(pt, spt, 8, true, ptmean, sptmean);
0299 LogTrace(metname) << " GSL Pt :" << ptmean << "+/-" << sptmean << endl;
0300 }
0301 }
0302
0303 void MuonDTSeedFromRecHits::computeMean(
0304 const double* pt, const double* weights, int sz, bool tossOutlyers, float& ptmean, float& sptmean) const {
0305 int n = 0;
0306 double ptTmp[8];
0307 double wtTmp[8];
0308 assert(sz <= 8);
0309
0310 for (int i = 0; i < 8; ++i) {
0311 ptTmp[i] = 0.;
0312 wtTmp[i] = 0;
0313 if (i < sz && pt[i] != 0) {
0314 ptTmp[n] = pt[i];
0315 wtTmp[n] = weights[i];
0316 ++n;
0317 }
0318 }
0319 if (n != 0) {
0320 if (n == 1) {
0321 ptmean = ptTmp[0];
0322 sptmean = 1 / sqrt(wtTmp[0]);
0323 } else {
0324 ptmean = gsl_stats_wmean(wtTmp, 1, ptTmp, 1, n);
0325 sptmean = sqrt(gsl_stats_wvariance_m(wtTmp, 1, ptTmp, 1, n, ptmean));
0326 }
0327
0328 if (tossOutlyers) {
0329
0330 for (int nm = 0; nm < 8; nm++) {
0331 if (ptTmp[nm] != 0 && fabs(ptTmp[nm] - ptmean) > 3 * (sptmean)) {
0332 wtTmp[nm] = 0.;
0333 }
0334 }
0335 ptmean = gsl_stats_wmean(wtTmp, 1, ptTmp, 1, n);
0336 sptmean = sqrt(gsl_stats_wvariance_m(wtTmp, 1, ptTmp, 1, n, ptmean));
0337 }
0338 }
0339 }