CommonAlignmentProducer/src/AlignmentMuonSelector.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "Alignment/CommonAlignmentProducer/interface/AlignmentMuonSelector.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005
0006 #include "TLorentzVector.h"
0007
0008 // constructor ----------------------------------------------------------------
0009
0010 AlignmentMuonSelector::AlignmentMuonSelector(const edm::ParameterSet& cfg)
0011     : applyBasicCuts(cfg.getParameter<bool>("applyBasicCuts")),
0012       applyNHighestPt(cfg.getParameter<bool>("applyNHighestPt")),
0013       applyMultiplicityFilter(cfg.getParameter<bool>("applyMultiplicityFilter")),
0014       applyMassPairFilter(cfg.getParameter<bool>("applyMassPairFilter")),
0015       nHighestPt(cfg.getParameter<int>("nHighestPt")),
0016       minMultiplicity(cfg.getParameter<int>("minMultiplicity")),
0017       pMin(cfg.getParameter<double>("pMin")),
0018       pMax(cfg.getParameter<double>("pMax")),
0019       ptMin(cfg.getParameter<double>("ptMin")),
0020       ptMax(cfg.getParameter<double>("ptMax")),
0021       etaMin(cfg.getParameter<double>("etaMin")),
0022       etaMax(cfg.getParameter<double>("etaMax")),
0023       phiMin(cfg.getParameter<double>("phiMin")),
0024       phiMax(cfg.getParameter<double>("phiMax")),
0025       nHitMinSA(cfg.getParameter<double>("nHitMinSA")),
0026       nHitMaxSA(cfg.getParameter<double>("nHitMaxSA")),
0027       chi2nMaxSA(cfg.getParameter<double>("chi2nMaxSA")),
0028       nHitMinGB(cfg.getParameter<double>("nHitMinGB")),
0029       nHitMaxGB(cfg.getParameter<double>("nHitMaxGB")),
0030       chi2nMaxGB(cfg.getParameter<double>("chi2nMaxGB")),
0031       nHitMinTO(cfg.getParameter<double>("nHitMinTO")),
0032       nHitMaxTO(cfg.getParameter<double>("nHitMaxTO")),
0033       chi2nMaxTO(cfg.getParameter<double>("chi2nMaxTO")),
0034       minMassPair(cfg.getParameter<double>("minMassPair")),
0035       maxMassPair(cfg.getParameter<double>("maxMassPair")) {
0036   if (applyBasicCuts)
0037     edm::LogInfo("AlignmentMuonSelector")
0038         << "applying basic muon cuts ..."
0039         << "\npmin,pmax:           " << pMin << "," << pMax << "\nptmin,ptmax:         " << ptMin << "," << ptMax
0040         << "\netamin,etamax:       " << etaMin << "," << etaMax << "\nphimin,phimax:       " << phiMin << "," << phiMax
0041         << "\nnhitminSA,nhitmaxSA: " << nHitMinSA << "," << nHitMaxSA << "\nchi2nmaxSA:          " << chi2nMaxSA << ","
0042         << "\nnhitminGB,nhitmaxGB: " << nHitMinGB << "," << nHitMaxGB << "\nchi2nmaxGB:          " << chi2nMaxGB << ","
0043         << "\nnhitminTO,nhitmaxTO: " << nHitMinTO << "," << nHitMaxTO << "\nchi2nmaxTO:          " << chi2nMaxTO;
0044
0045   if (applyNHighestPt)
0046     edm::LogInfo("AlignmentMuonSelector") << "filter N muons with highest Pt N=" << nHighestPt;
0047
0048   if (applyMultiplicityFilter)
0049     edm::LogInfo("AlignmentMuonSelector") << "apply multiplicity filter N>=" << minMultiplicity;
0050
0051   if (applyMassPairFilter)
0052     edm::LogInfo("AlignmentMuonSelector")
0053         << "apply Mass Pair filter minMassPair=" << minMassPair << " maxMassPair=" << maxMassPair;
0054 }
0055
0056 // destructor -----------------------------------------------------------------
0057
0058 AlignmentMuonSelector::~AlignmentMuonSelector() {}
0059
0060 // do selection ---------------------------------------------------------------
0061
0062 AlignmentMuonSelector::Muons AlignmentMuonSelector::select(const Muons& muons, const edm::Event& evt) const {
0063   Muons result = muons;
0064
0065   // apply basic muon cuts (if selected)
0066   if (applyBasicCuts)
0067     result = this->basicCuts(result);
0068
0069   // filter N muons with highest Pt (if selected)
0070   if (applyNHighestPt)
0071     result = this->theNHighestPtMuons(result);
0072
0073   // apply minimum multiplicity requirement (if selected)
0074   if (applyMultiplicityFilter) {
0075     if (result.size() < (unsigned int)minMultiplicity)
0076       result.clear();
0077   }
0078
0079   // apply mass pair requirement (if selected)
0080   if (applyMassPairFilter) {
0081     if (result.size() < 2)
0082       result.clear();  // at least 2 muons are require for a mass pair...
0083     else
0084       result = this->theBestMassPairCombinationMuons(result);
0085   }
0086
0087   edm::LogInfo("AlignmentMuonSelector") << "muons all,kept: " << muons.size() << "," << result.size();
0088
0089   return result;
0090 }
0091
0092 // make basic cuts ------------------------------------------------------------
0093
0094 AlignmentMuonSelector::Muons AlignmentMuonSelector::basicCuts(const Muons& muons) const {
0095   Muons result;
0096
0097   for (Muons::const_iterator it = muons.begin(); it != muons.end(); ++it) {
0098     const reco::Muon* muonp = *it;
0099     float p = muonp->p();
0100     float pt = muonp->pt();
0101     float eta = muonp->eta();
0102     float phi = muonp->phi();
0103
0104     int nhitSA = 0;
0105     float chi2nSA = 9999.;
0106     if (muonp->isStandAloneMuon()) {
0107       nhitSA = muonp->standAloneMuon()->numberOfValidHits();  // standAlone Muon
0108       chi2nSA = muonp->standAloneMuon()->normalizedChi2();    // standAlone Muon
0109     }
0110     int nhitGB = 0;
0111     float chi2nGB = 9999.;
0112     if (muonp->isGlobalMuon()) {
0113       nhitGB = muonp->combinedMuon()->numberOfValidHits();  // global Muon
0114       chi2nGB = muonp->combinedMuon()->normalizedChi2();    // global Muon
0115     }
0116     int nhitTO = 0;
0117     float chi2nTO = 9999.;
0118     if (muonp->isTrackerMuon()) {
0119       nhitTO = muonp->track()->numberOfValidHits();  // Tracker Only
0120       chi2nTO = muonp->track()->normalizedChi2();    // Tracker Only
0121     }
0122     edm::LogInfo("AlignmentMuonSelector")
0123         << " pt,eta,phi,nhitSA,chi2nSA,nhitGB,chi2nGB,nhitTO,chi2nTO: " << pt << "," << eta << "," << phi << ","
0124         << nhitSA << "," << chi2nSA << "," << nhitGB << "," << chi2nGB << "," << nhitTO << "," << chi2nTO;
0125
0126     if (p > pMin && p < pMax && pt > ptMin && pt < ptMax && eta > etaMin && eta < etaMax && phi > phiMin &&
0127         phi < phiMax && nhitSA >= nHitMinSA && nhitSA <= nHitMaxSA && chi2nSA < chi2nMaxSA && nhitGB >= nHitMinGB &&
0128         nhitGB <= nHitMaxGB && chi2nGB < chi2nMaxGB && nhitTO >= nHitMinTO && nhitTO <= nHitMaxTO &&
0129         chi2nTO < chi2nMaxTO) {
0130       result.push_back(muonp);
0131     }
0132   }
0133
0134   return result;
0135 }
0136
0137 //-----------------------------------------------------------------------------
0138
0139 AlignmentMuonSelector::Muons AlignmentMuonSelector::theNHighestPtMuons(const Muons& muons) const {
0140   Muons sortedMuons = muons;
0141   Muons result;
0142
0143   // sort in pt
0144   std::sort(sortedMuons.begin(), sortedMuons.end(), ptComparator);
0145
0146   // copy theMuonMult highest pt muons to result vector
0147   int n = 0;
0148   for (Muons::const_iterator it = sortedMuons.begin(); it != sortedMuons.end(); ++it) {
0149     if (n < nHighestPt) {
0150       result.push_back(*it);
0151       n++;
0152     }
0153   }
0154
0155   return result;
0156 }
0157
0158 //-----------------------------------------------------------------------------
0159
0160 AlignmentMuonSelector::Muons AlignmentMuonSelector::theBestMassPairCombinationMuons(const Muons& muons) const {
0161   Muons sortedMuons = muons;
0162   Muons result;
0163   TLorentzVector mu1, mu2, pair;
0164   double mass = 0, minDiff = 999999.;
0165
0166   // sort in pt
0167   std::sort(sortedMuons.begin(), sortedMuons.end(), ptComparator);
0168
0169   // copy best mass pair combination muons to result vector
0170   // Criteria:
0171   // a) maxMassPair !=    minMassPair: the two highest pt muons with mass pair inside the given mass window
0172   // b) maxMassPair ==    minMassPair: the muon pair with massPair closest to given mass value
0173   for (Muons::const_iterator it1 = sortedMuons.begin(); it1 != sortedMuons.end(); ++it1) {
0174     for (Muons::const_iterator it2 = it1 + 1; it2 != sortedMuons.end(); ++it2) {
0175       mu1 = TLorentzVector((*it1)->momentum().x(), (*it1)->momentum().y(), (*it1)->momentum().z(), (*it1)->p());
0176       mu2 = TLorentzVector((*it2)->momentum().x(), (*it2)->momentum().y(), (*it2)->momentum().z(), (*it2)->p());
0177       pair = mu1 + mu2;
0178       mass = pair.M();
0179
0180       if (maxMassPair != minMassPair) {
0181         if (mass < maxMassPair && mass > minMassPair) {
0182           result.push_back(*it1);
0183           result.push_back(*it2);
0184           break;
0185         }
0186       } else {
0187         if (fabs(mass - maxMassPair) < minDiff) {
0188           minDiff = fabs(mass - maxMassPair);
0189           result.clear();
0190           result.push_back(*it1);
0191           result.push_back(*it2);
0192         }
0193       }
0194     }
0195   }
0196
0197   return result;
0198 }