MuonSeedGenerator/src/MuonDTSeedFromRecHits.cc

0001 /**
0002  *  See header file for a description of this class.
0003  *
0004  *
0005  *  \author A. Vitelli - INFN Torino, V.Palichik
0006  *  \author porting  R. Bellan
0007  *
0008  */
0009 #include "RecoMuon/MuonSeedGenerator/src/MuonDTSeedFromRecHits.h"
0010 #include "RecoMuon/MuonSeedGenerator/src/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   // these weights are supposed to be 1/sigma^2(pt), but they're so small.
0030   // Instead of the 0.2-0.5 GeV here, we'll add something extra in quadrature later
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   /// compute pts with vertex constraint
0036   computePtWithVtx(pt, spt);
0037
0038   /// now w/o vertex constrain
0039   computePtWithoutVtx(pt, spt);
0040
0041   // some dump...
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   /// now combine all pt estimate
0053   float ptmean = 0.;
0054   float sptmean = 0.;
0055   //@@ Use Shih-Chuan's
0056   computeBestPt(pt, spt, ptmean, sptmean);
0057
0058   // add an extra term to the error in quadrature, 30% of pT per point
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   // take the best candidate
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   }  //  +v.
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     //+vvp !:
0160
0161     float eta1 = (*iter)->globalPosition().eta();
0162     if (fabs(eta1 - eta0) > .2)
0163       continue;  //   !!! +vvp
0164
0165     // assign Pt from MB1 & vtx
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     // assign Pt from MB2 & vtx
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     //+vvp !:
0192     float eta1 = (*iter)->globalPosition().eta();
0193     if (fabs(eta1 - eta0) > .2)
0194       continue;  //   !!! +vvp
0195
0196     for (MuonRecHitContainer::const_iterator iter2 = theRhits.begin(); iter2 != iter; iter2++) {
0197       //+vvp !:
0198       float eta2 = (*iter2)->globalPosition().eta();
0199       if (fabs(eta2 - eta0) > .2)
0200         continue;  //   !!! +vvp
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: {  //MB1-MB2
0214           pt[2] = thispt;
0215           break;
0216         }
0217         case 13: {  //MB1-MB3
0218           pt[3] = thispt;
0219           break;
0220         }
0221         case 14: {  //MB1-MB4
0222           pt[5] = thispt;
0223           break;
0224         }
0225         case 23: {  //MB2-MB3
0226           pt[4] = thispt;
0227           break;
0228         }
0229         case 24: {  //MB2-MB4
0230           pt[6] = thispt;
0231           break;
0232         }
0233         case 34: {  //MB3-MB4
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     /// just one pt estimate:  use it.
0261     ptmean = pt[igood];
0262     sptmean = 1 / sqrt(spt[igood]);
0263   } else if (nTotal == 0) {
0264     // No estimate (e.g. only one rechit!)
0265     ptmean = 50;
0266     sptmean = 30;
0267   } else {
0268     /// more than a pt estimate, do all the job.
0269     // calculate pt with vertex
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     // FIXME: temp hack
0276     if (ptvtx != 0.) {
0277       ptmean = ptvtx;
0278       sptmean = sptvtx;
0279       return;
0280       //
0281     }
0282     // calculate pt w/o vertex
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     // decide wheter the muon comes or not from vertex
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     // now choose the "right" pt => weighted mean
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       // Recompute mean with a cut at 3 sigma
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 }