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 void AlignmentMuonSelector::fillPSetDescription(edm::ParameterSetDescription& desc) {
0057   desc.add<bool>("applyBasicCuts", true);
0058   desc.add<bool>("applyNHighestPt", false);
0059   desc.add<bool>("applyMultiplicityFilter", false);
0060   desc.add<bool>("applyMassPairFilter", false);
0061   desc.add<int>("nHighestPt", 2);
0062   desc.add<int>("minMultiplicity", 1);
0063   desc.add<double>("pMin", 0.0);
0064   desc.add<double>("pMax", 999999.0);
0065   desc.add<double>("ptMin", 10.);
0066   desc.add<double>("ptMax", 999999.0);
0067   desc.add<double>("etaMin", -2.4);
0068   desc.add<double>("etaMax", 2.4);
0069   desc.add<double>("phiMin", -3.1416);
0070   desc.add<double>("phiMax", 3.1416);
0071   desc.add<double>("nHitMinSA", 0.0);
0072   desc.add<double>("nHitMaxSA", 999999.0);
0073   desc.add<double>("chi2nMaxSA", 999999.0);
0074   desc.add<double>("nHitMinGB", 0.0);
0075   desc.add<double>("nHitMaxGB", 999999.0);
0076   desc.add<double>("chi2nMaxGB", 999999.0);
0077   desc.add<double>("nHitMinTO", 0.0);
0078   desc.add<double>("nHitMaxTO", 999999.0);
0079   desc.add<double>("chi2nMaxTO", 999999.0);
0080   desc.add<double>("minMassPair", 89.0)
0081       ->setComment(
0082           "copy best mass pair combination muons to result vector \n Criteria: \n a) maxMassPair != minMassPair: the "
0083           "two highest pt muons with mass pair inside the given mass window \n b) maxMassPair == minMassPair: the muon "
0084           "pair with mass pair closest to given mass value");
0085   desc.add<double>("maxMassPair", 90.0);
0086 }
0087
0088 // destructor -----------------------------------------------------------------
0089
0090 AlignmentMuonSelector::~AlignmentMuonSelector() {}
0091
0092 // do selection ---------------------------------------------------------------
0093
0094 AlignmentMuonSelector::Muons AlignmentMuonSelector::select(const Muons& muons, const edm::Event& evt) const {
0095   Muons result = muons;
0096
0097   // apply basic muon cuts (if selected)
0098   if (applyBasicCuts)
0099     result = this->basicCuts(result);
0100
0101   // filter N muons with highest Pt (if selected)
0102   if (applyNHighestPt)
0103     result = this->theNHighestPtMuons(result);
0104
0105   // apply minimum multiplicity requirement (if selected)
0106   if (applyMultiplicityFilter) {
0107     if (result.size() < (unsigned int)minMultiplicity)
0108       result.clear();
0109   }
0110
0111   // apply mass pair requirement (if selected)
0112   if (applyMassPairFilter) {
0113     if (result.size() < 2)
0114       result.clear();  // at least 2 muons are require for a mass pair...
0115     else
0116       result = this->theBestMassPairCombinationMuons(result);
0117   }
0118
0119   edm::LogInfo("AlignmentMuonSelector") << "muons all,kept: " << muons.size() << "," << result.size();
0120
0121   return result;
0122 }
0123
0124 // make basic cuts ------------------------------------------------------------
0125
0126 AlignmentMuonSelector::Muons AlignmentMuonSelector::basicCuts(const Muons& muons) const {
0127   Muons result;
0128
0129   for (Muons::const_iterator it = muons.begin(); it != muons.end(); ++it) {
0130     const reco::Muon* muonp = *it;
0131     float p = muonp->p();
0132     float pt = muonp->pt();
0133     float eta = muonp->eta();
0134     float phi = muonp->phi();
0135
0136     int nhitSA = 0;
0137     float chi2nSA = 9999.;
0138     if (muonp->isStandAloneMuon()) {
0139       nhitSA = muonp->standAloneMuon()->numberOfValidHits();  // standAlone Muon
0140       chi2nSA = muonp->standAloneMuon()->normalizedChi2();    // standAlone Muon
0141     }
0142     int nhitGB = 0;
0143     float chi2nGB = 9999.;
0144     if (muonp->isGlobalMuon()) {
0145       nhitGB = muonp->combinedMuon()->numberOfValidHits();  // global Muon
0146       chi2nGB = muonp->combinedMuon()->normalizedChi2();    // global Muon
0147     }
0148     int nhitTO = 0;
0149     float chi2nTO = 9999.;
0150     if (muonp->isTrackerMuon()) {
0151       nhitTO = muonp->track()->numberOfValidHits();  // Tracker Only
0152       chi2nTO = muonp->track()->normalizedChi2();    // Tracker Only
0153     }
0154     edm::LogInfo("AlignmentMuonSelector")
0155         << " pt,eta,phi,nhitSA,chi2nSA,nhitGB,chi2nGB,nhitTO,chi2nTO: " << pt << "," << eta << "," << phi << ","
0156         << nhitSA << "," << chi2nSA << "," << nhitGB << "," << chi2nGB << "," << nhitTO << "," << chi2nTO;
0157
0158     if (p > pMin && p < pMax && pt > ptMin && pt < ptMax && eta > etaMin && eta < etaMax && phi > phiMin &&
0159         phi < phiMax && nhitSA >= nHitMinSA && nhitSA <= nHitMaxSA && chi2nSA < chi2nMaxSA && nhitGB >= nHitMinGB &&
0160         nhitGB <= nHitMaxGB && chi2nGB < chi2nMaxGB && nhitTO >= nHitMinTO && nhitTO <= nHitMaxTO &&
0161         chi2nTO < chi2nMaxTO) {
0162       result.push_back(muonp);
0163     }
0164   }
0165
0166   return result;
0167 }
0168
0169 //-----------------------------------------------------------------------------
0170
0171 AlignmentMuonSelector::Muons AlignmentMuonSelector::theNHighestPtMuons(const Muons& muons) const {
0172   Muons sortedMuons = muons;
0173   Muons result;
0174
0175   // sort in pt
0176   std::sort(sortedMuons.begin(), sortedMuons.end(), ptComparator);
0177
0178   // copy theMuonMult highest pt muons to result vector
0179   int n = 0;
0180   for (Muons::const_iterator it = sortedMuons.begin(); it != sortedMuons.end(); ++it) {
0181     if (n < nHighestPt) {
0182       result.push_back(*it);
0183       n++;
0184     }
0185   }
0186
0187   return result;
0188 }
0189
0190 //-----------------------------------------------------------------------------
0191
0192 AlignmentMuonSelector::Muons AlignmentMuonSelector::theBestMassPairCombinationMuons(const Muons& muons) const {
0193   Muons sortedMuons = muons;
0194   Muons result;
0195   TLorentzVector mu1, mu2, pair;
0196   double mass = 0, minDiff = 999999.;
0197
0198   // sort in pt
0199   std::sort(sortedMuons.begin(), sortedMuons.end(), ptComparator);
0200
0201   // copy best mass pair combination muons to result vector
0202   // Criteria:
0203   // a) maxMassPair !=    minMassPair: the two highest pt muons with mass pair inside the given mass window
0204   // b) maxMassPair ==    minMassPair: the muon pair with massPair closest to given mass value
0205   for (Muons::const_iterator it1 = sortedMuons.begin(); it1 != sortedMuons.end(); ++it1) {
0206     for (Muons::const_iterator it2 = it1 + 1; it2 != sortedMuons.end(); ++it2) {
0207       mu1 = TLorentzVector((*it1)->momentum().x(), (*it1)->momentum().y(), (*it1)->momentum().z(), (*it1)->p());
0208       mu2 = TLorentzVector((*it2)->momentum().x(), (*it2)->momentum().y(), (*it2)->momentum().z(), (*it2)->p());
0209       pair = mu1 + mu2;
0210       mass = pair.M();
0211
0212       if (maxMassPair != minMassPair) {
0213         if (mass < maxMassPair && mass > minMassPair) {
0214           result.push_back(*it1);
0215           result.push_back(*it2);
0216           break;
0217         }
0218       } else {
0219         if (fabs(mass - maxMassPair) < minDiff) {
0220           minDiff = fabs(mass - maxMassPair);
0221           result.clear();
0222           result.push_back(*it1);
0223           result.push_back(*it2);
0224         }
0225       }
0226     }
0227   }
0228
0229   return result;
0230 }