CommonAlignmentProducer/src/AlignmentTwoBoyDecayTrackSelector.cc

0001 //Framework
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 //DataFormats
0009 #include <DataFormats/TrackReco/interface/Track.h>
0010 #include <DataFormats/METReco/interface/CaloMET.h>
0011 #include <DataFormats/Math/interface/deltaPhi.h>
0012
0013 //STL
0014 #include <cmath>
0015 //ROOT
0016 #include "TLorentzVector.h"
0017
0018 #include "Alignment/CommonAlignmentProducer/interface/AlignmentTwoBodyDecayTrackSelector.h"
0019 //TODO put those namespaces into functions?
0020 using namespace std;
0021 using namespace edm;
0022 // constructor ----------------------------------------------------------------
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");  //Number of candidates to keep
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   // Mass range filter
0069   desc.add<bool>("applyMassrangeFilter", false);
0070   desc.add<double>("daughterMass", 0.105);  // GeV
0071
0072   // Charge-related parameters
0073   desc.add<bool>("useUnsignedCharge", true);
0074   desc.add<int>("charge", 0);
0075
0076   // Missing ET source
0077   desc.add<edm::InputTag>("missingETSource", edm::InputTag("met"));
0078   desc.add<double>("maxXMass", 15000.0);  // GeV
0079   desc.add<double>("minXMass", 0.0);      // GeV
0080
0081   // Acoplanarity settings
0082   desc.add<double>("acoplanarDistance", 1.0);  // Radian
0083
0084   // Filters
0085   desc.add<bool>("applyChargeFilter", false);
0086   desc.add<bool>("applyAcoplanarityFilter", false);
0087   desc.add<bool>("applyMissingETFilter", false);
0088
0089   // Candidate selection
0090   desc.add<unsigned int>("numberOfCandidates", 1);
0091   desc.add<bool>("applySecThreshold", false);
0092   desc.add<double>("secondThreshold", 6.0);
0093 }
0094
0095 // destructor -----------------------------------------------------------------
0096
0097 AlignmentTwoBodyDecayTrackSelector::~AlignmentTwoBodyDecayTrackSelector() {}
0098
0099 ///returns if any of the Filters is used.
0100 bool AlignmentTwoBodyDecayTrackSelector::useThisFilter() {
0101   return theMassrangeSwitch || theChargeSwitch || theAcoplanarityFilterSwitch;
0102 }
0103
0104 // do selection ---------------------------------------------------------------
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 ///checks if the mass of the X is in the mass region
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 ///checks if the mass of the X is in the mass region adding missing E_T
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 ///checks if the mother has charge = [theCharge]
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 ///checks if the [cands] are acoplanar (returns empty set if not)
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 ///checks if the [cands] are acoplanar (returns empty set if not)
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 //===================HELPERS===================
0294
0295 ///print Information on Track-Collection
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 }