File indexing completed on 2025-02-05 23:50:58
0001
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Utilities/interface/EDMException.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007
0008
0009 #include <DataFormats/TrackReco/interface/Track.h>
0010 #include <DataFormats/METReco/interface/CaloMET.h>
0011 #include <DataFormats/Math/interface/deltaPhi.h>
0012
0013
0014 #include <cmath>
0015
0016 #include "TLorentzVector.h"
0017
0018 #include "Alignment/CommonAlignmentProducer/interface/AlignmentTwoBodyDecayTrackSelector.h"
0019
0020 using namespace std;
0021 using namespace edm;
0022
0023
0024 AlignmentTwoBodyDecayTrackSelector::AlignmentTwoBodyDecayTrackSelector(const edm::ParameterSet& cfg,
0025 edm::ConsumesCollector& iC) {
0026 LogDebug("Alignment") << "> applying two body decay Trackfilter ...";
0027 theMassrangeSwitch = cfg.getParameter<bool>("applyMassrangeFilter");
0028 if (theMassrangeSwitch) {
0029 theMinMass = cfg.getParameter<double>("minXMass");
0030 theMaxMass = cfg.getParameter<double>("maxXMass");
0031 theDaughterMass = cfg.getParameter<double>("daughterMass");
0032 theCandNumber = cfg.getParameter<unsigned int>("numberOfCandidates");
0033 secThrBool = cfg.getParameter<bool>("applySecThreshold");
0034 thesecThr = cfg.getParameter<double>("secondThreshold");
0035 LogDebug("Alignment") << "> Massrange min,max : " << theMinMass << "," << theMaxMass
0036 << "\n> Mass of daughter Particle : " << theDaughterMass;
0037
0038 } else {
0039 theMinMass = 0;
0040 theMaxMass = 0;
0041 theDaughterMass = 0;
0042 }
0043 theChargeSwitch = cfg.getParameter<bool>("applyChargeFilter");
0044 if (theChargeSwitch) {
0045 theCharge = cfg.getParameter<int>("charge");
0046 theUnsignedSwitch = cfg.getParameter<bool>("useUnsignedCharge");
0047 if (theUnsignedSwitch)
0048 theCharge = std::abs(theCharge);
0049 LogDebug("Alignment") << "> Desired Charge, unsigned: " << theCharge << " , " << theUnsignedSwitch;
0050 } else {
0051 theCharge = 0;
0052 theUnsignedSwitch = true;
0053 }
0054 theMissingETSwitch = cfg.getParameter<bool>("applyMissingETFilter");
0055 if (theMissingETSwitch) {
0056 edm::InputTag theMissingETSource = cfg.getParameter<InputTag>("missingETSource");
0057 theMissingETToken = iC.consumes<reco::CaloMETCollection>(theMissingETSource);
0058 LogDebug("Alignment") << "> missing Et Source: " << theMissingETSource;
0059 }
0060 theAcoplanarityFilterSwitch = cfg.getParameter<bool>("applyAcoplanarityFilter");
0061 if (theAcoplanarityFilterSwitch) {
0062 theAcoplanarDistance = cfg.getParameter<double>("acoplanarDistance");
0063 LogDebug("Alignment") << "> Acoplanar Distance: " << theAcoplanarDistance;
0064 }
0065 }
0066
0067 void AlignmentTwoBodyDecayTrackSelector::fillPSetDescription(edm::ParameterSetDescription& desc) {
0068
0069 desc.add<bool>("applyMassrangeFilter", false);
0070 desc.add<double>("daughterMass", 0.105);
0071
0072
0073 desc.add<bool>("useUnsignedCharge", true);
0074 desc.add<int>("charge", 0);
0075
0076
0077 desc.add<edm::InputTag>("missingETSource", edm::InputTag("met"));
0078 desc.add<double>("maxXMass", 15000.0);
0079 desc.add<double>("minXMass", 0.0);
0080
0081
0082 desc.add<double>("acoplanarDistance", 1.0);
0083
0084
0085 desc.add<bool>("applyChargeFilter", false);
0086 desc.add<bool>("applyAcoplanarityFilter", false);
0087 desc.add<bool>("applyMissingETFilter", false);
0088
0089
0090 desc.add<unsigned int>("numberOfCandidates", 1);
0091 desc.add<bool>("applySecThreshold", false);
0092 desc.add<double>("secondThreshold", 6.0);
0093 }
0094
0095
0096
0097 AlignmentTwoBodyDecayTrackSelector::~AlignmentTwoBodyDecayTrackSelector() {}
0098
0099
0100 bool AlignmentTwoBodyDecayTrackSelector::useThisFilter() {
0101 return theMassrangeSwitch || theChargeSwitch || theAcoplanarityFilterSwitch;
0102 }
0103
0104
0105
0106 AlignmentTwoBodyDecayTrackSelector::Tracks AlignmentTwoBodyDecayTrackSelector::select(const Tracks& tracks,
0107 const edm::Event& iEvent,
0108 const edm::EventSetup& iSetup) {
0109 Tracks result = tracks;
0110
0111 if (theMassrangeSwitch) {
0112 if (theMissingETSwitch)
0113 result = checkMETMass(result, iEvent);
0114 else
0115 result = checkMass(result);
0116 }
0117
0118 LogDebug("Alignment") << "> TwoBodyDecay tracks all,kept: " << tracks.size() << "," << result.size();
0119 return result;
0120 }
0121
0122
0123 AlignmentTwoBodyDecayTrackSelector::Tracks AlignmentTwoBodyDecayTrackSelector::checkMass(const Tracks& cands) const {
0124 Tracks result;
0125
0126 LogDebug("Alignment") << "> cands size : " << cands.size();
0127
0128 if (cands.size() < 2)
0129 return result;
0130
0131 TLorentzVector track0;
0132 TLorentzVector track1;
0133 TLorentzVector mother;
0134 typedef pair<const reco::Track*, const reco::Track*> constTrackPair;
0135 typedef pair<double, constTrackPair> candCollectionItem;
0136 vector<candCollectionItem> candCollection;
0137
0138 for (unsigned int iCand = 0; iCand < cands.size(); iCand++) {
0139 track0.SetXYZT(cands.at(iCand)->px(),
0140 cands.at(iCand)->py(),
0141 cands.at(iCand)->pz(),
0142 sqrt(cands.at(iCand)->p() * cands.at(iCand)->p() + theDaughterMass * theDaughterMass));
0143
0144 for (unsigned int jCand = iCand + 1; jCand < cands.size(); jCand++) {
0145 track1.SetXYZT(cands.at(jCand)->px(),
0146 cands.at(jCand)->py(),
0147 cands.at(jCand)->pz(),
0148 sqrt(cands.at(jCand)->p() * cands.at(jCand)->p() + theDaughterMass * theDaughterMass));
0149 if (secThrBool == true && track1.Pt() < thesecThr && track0.Pt() < thesecThr)
0150 continue;
0151 mother = track0 + track1;
0152
0153 const reco::Track* trk1 = cands.at(iCand);
0154 const reco::Track* trk2 = cands.at(jCand);
0155
0156 bool correctCharge = true;
0157 if (theChargeSwitch)
0158 correctCharge = this->checkCharge(trk1, trk2);
0159
0160 bool acoplanarTracks = true;
0161 if (theAcoplanarityFilterSwitch)
0162 acoplanarTracks = this->checkAcoplanarity(trk1, trk2);
0163
0164 if (mother.M() > theMinMass && mother.M() < theMaxMass && correctCharge && acoplanarTracks) {
0165 candCollection.push_back(candCollectionItem(mother.Pt(), constTrackPair(trk1, trk2)));
0166 }
0167 }
0168 }
0169
0170 if (candCollection.empty())
0171 return result;
0172
0173 sort(candCollection.begin(), candCollection.end(), [](auto& a, auto& b) { return a.first > b.first; });
0174
0175 std::map<const reco::Track*, unsigned int> uniqueTrackIndex;
0176 std::map<const reco::Track*, unsigned int>::iterator it;
0177 for (unsigned int i = 0; i < candCollection.size() && i < theCandNumber; i++) {
0178 constTrackPair& trackPair = candCollection[i].second;
0179
0180 it = uniqueTrackIndex.find(trackPair.first);
0181 if (it == uniqueTrackIndex.end()) {
0182 result.push_back(trackPair.first);
0183 uniqueTrackIndex[trackPair.first] = i;
0184 }
0185
0186 it = uniqueTrackIndex.find(trackPair.second);
0187 if (it == uniqueTrackIndex.end()) {
0188 result.push_back(trackPair.second);
0189 uniqueTrackIndex[trackPair.second] = i;
0190 }
0191 }
0192
0193 return result;
0194 }
0195
0196
0197 AlignmentTwoBodyDecayTrackSelector::Tracks AlignmentTwoBodyDecayTrackSelector::checkMETMass(
0198 const Tracks& cands, const edm::Event& iEvent) const {
0199 Tracks result;
0200
0201 LogDebug("Alignment") << "> cands size : " << cands.size();
0202
0203 if (cands.empty())
0204 return result;
0205
0206 TLorentzVector track;
0207 TLorentzVector met4;
0208 TLorentzVector mother;
0209
0210 Handle<reco::CaloMETCollection> missingET;
0211 iEvent.getByToken(theMissingETToken, missingET);
0212 if (!missingET.isValid()) {
0213 LogError("Alignment") << "@SUB=AlignmentTwoBodyDecayTrackSelector::checkMETMass"
0214 << "> could not optain missingET Collection!";
0215 return result;
0216 }
0217
0218 typedef pair<double, const reco::Track*> candCollectionItem;
0219 vector<candCollectionItem> candCollection;
0220
0221 for (reco::CaloMETCollection::const_iterator itMET = missingET->begin(); itMET != missingET->end(); ++itMET) {
0222 met4.SetXYZT((*itMET).px(), (*itMET).py(), (*itMET).pz(), (*itMET).p());
0223
0224 for (unsigned int iCand = 0; iCand < cands.size(); iCand++) {
0225 track.SetXYZT(cands.at(iCand)->px(),
0226 cands.at(iCand)->py(),
0227 cands.at(iCand)->pz(),
0228 sqrt(cands.at(iCand)->p() * cands.at(iCand)->p() + theDaughterMass * theDaughterMass));
0229
0230 mother = track + met4;
0231
0232 const reco::Track* trk = cands.at(iCand);
0233 const reco::CaloMET* met = &(*itMET);
0234
0235 bool correctCharge = true;
0236 if (theChargeSwitch)
0237 correctCharge = this->checkCharge(trk);
0238
0239 bool acoplanarTracks = true;
0240 if (theAcoplanarityFilterSwitch)
0241 acoplanarTracks = this->checkMETAcoplanarity(trk, met);
0242
0243 if (mother.M() > theMinMass && mother.M() < theMaxMass && correctCharge && acoplanarTracks) {
0244 candCollection.push_back(candCollectionItem(mother.Pt(), trk));
0245 }
0246 }
0247 }
0248
0249 if (candCollection.empty())
0250 return result;
0251
0252 sort(candCollection.begin(), candCollection.end(), [](auto& a, auto& b) { return a.first > b.first; });
0253
0254 std::map<const reco::Track*, unsigned int> uniqueTrackIndex;
0255 std::map<const reco::Track*, unsigned int>::iterator it;
0256 for (unsigned int i = 0; i < candCollection.size() && i < theCandNumber; i++) {
0257 it = uniqueTrackIndex.find(candCollection[i].second);
0258 if (it == uniqueTrackIndex.end()) {
0259 result.push_back(candCollection[i].second);
0260 uniqueTrackIndex[candCollection[i].second] = i;
0261 }
0262 }
0263
0264 return result;
0265 }
0266
0267
0268 bool AlignmentTwoBodyDecayTrackSelector::checkCharge(const reco::Track* trk1, const reco::Track* trk2) const {
0269 int sumCharge = trk1->charge();
0270 if (trk2)
0271 sumCharge += trk2->charge();
0272 if (theUnsignedSwitch)
0273 sumCharge = std::abs(sumCharge);
0274 if (sumCharge == theCharge)
0275 return true;
0276 return false;
0277 }
0278
0279
0280 bool AlignmentTwoBodyDecayTrackSelector::checkAcoplanarity(const reco::Track* trk1, const reco::Track* trk2) const {
0281 if (fabs(deltaPhi(trk1->phi(), trk2->phi() - M_PI)) < theAcoplanarDistance)
0282 return true;
0283 return false;
0284 }
0285
0286
0287 bool AlignmentTwoBodyDecayTrackSelector::checkMETAcoplanarity(const reco::Track* trk1, const reco::CaloMET* met) const {
0288 if (fabs(deltaPhi(trk1->phi(), met->phi() - M_PI)) < theAcoplanarDistance)
0289 return true;
0290 return false;
0291 }
0292
0293
0294
0295
0296 void AlignmentTwoBodyDecayTrackSelector::printTracks(const Tracks& col) const {
0297 int count = 0;
0298 LogDebug("Alignment") << ">......................................";
0299 for (Tracks::const_iterator it = col.begin(); it < col.end(); ++it, ++count) {
0300 LogDebug("Alignment") << "> Track No. " << count << ": p = (" << (*it)->px() << "," << (*it)->py() << ","
0301 << (*it)->pz() << ")\n"
0302 << "> pT = " << (*it)->pt() << " eta = " << (*it)->eta()
0303 << " charge = " << (*it)->charge();
0304 }
0305 LogDebug("Alignment") << ">......................................";
0306 }