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 // destructor -----------------------------------------------------------------
0068
0069 AlignmentTwoBodyDecayTrackSelector::~AlignmentTwoBodyDecayTrackSelector() {}
0070
0071 ///returns if any of the Filters is used.
0072 bool AlignmentTwoBodyDecayTrackSelector::useThisFilter() {
0073   return theMassrangeSwitch || theChargeSwitch || theAcoplanarityFilterSwitch;
0074 }
0075
0076 // do selection ---------------------------------------------------------------
0077
0078 AlignmentTwoBodyDecayTrackSelector::Tracks AlignmentTwoBodyDecayTrackSelector::select(const Tracks& tracks,
0079                                                                                       const edm::Event& iEvent,
0080                                                                                       const edm::EventSetup& iSetup) {
0081   Tracks result = tracks;
0082
0083   if (theMassrangeSwitch) {
0084     if (theMissingETSwitch)
0085       result = checkMETMass(result, iEvent);
0086     else
0087       result = checkMass(result);
0088   }
0089
0090   LogDebug("Alignment") << ">  TwoBodyDecay tracks all,kept: " << tracks.size() << "," << result.size();
0091   return result;
0092 }
0093
0094 ///checks if the mass of the X is in the mass region
0095 AlignmentTwoBodyDecayTrackSelector::Tracks AlignmentTwoBodyDecayTrackSelector::checkMass(const Tracks& cands) const {
0096   Tracks result;
0097
0098   LogDebug("Alignment") << ">  cands size : " << cands.size();
0099
0100   if (cands.size() < 2)
0101     return result;
0102
0103   TLorentzVector track0;
0104   TLorentzVector track1;
0105   TLorentzVector mother;
0106   typedef pair<const reco::Track*, const reco::Track*> constTrackPair;
0107   typedef pair<double, constTrackPair> candCollectionItem;
0108   vector<candCollectionItem> candCollection;
0109
0110   for (unsigned int iCand = 0; iCand < cands.size(); iCand++) {
0111     track0.SetXYZT(cands.at(iCand)->px(),
0112                    cands.at(iCand)->py(),
0113                    cands.at(iCand)->pz(),
0114                    sqrt(cands.at(iCand)->p() * cands.at(iCand)->p() + theDaughterMass * theDaughterMass));
0115
0116     for (unsigned int jCand = iCand + 1; jCand < cands.size(); jCand++) {
0117       track1.SetXYZT(cands.at(jCand)->px(),
0118                      cands.at(jCand)->py(),
0119                      cands.at(jCand)->pz(),
0120                      sqrt(cands.at(jCand)->p() * cands.at(jCand)->p() + theDaughterMass * theDaughterMass));
0121       if (secThrBool == true && track1.Pt() < thesecThr && track0.Pt() < thesecThr)
0122         continue;
0123       mother = track0 + track1;
0124
0125       const reco::Track* trk1 = cands.at(iCand);
0126       const reco::Track* trk2 = cands.at(jCand);
0127
0128       bool correctCharge = true;
0129       if (theChargeSwitch)
0130         correctCharge = this->checkCharge(trk1, trk2);
0131
0132       bool acoplanarTracks = true;
0133       if (theAcoplanarityFilterSwitch)
0134         acoplanarTracks = this->checkAcoplanarity(trk1, trk2);
0135
0136       if (mother.M() > theMinMass && mother.M() < theMaxMass && correctCharge && acoplanarTracks) {
0137         candCollection.push_back(candCollectionItem(mother.Pt(), constTrackPair(trk1, trk2)));
0138       }
0139     }
0140   }
0141
0142   if (candCollection.empty())
0143     return result;
0144
0145   sort(candCollection.begin(), candCollection.end(), [](auto& a, auto& b) { return a.first > b.first; });
0146
0147   std::map<const reco::Track*, unsigned int> uniqueTrackIndex;
0148   std::map<const reco::Track*, unsigned int>::iterator it;
0149   for (unsigned int i = 0; i < candCollection.size() && i < theCandNumber; i++) {
0150     constTrackPair& trackPair = candCollection[i].second;
0151
0152     it = uniqueTrackIndex.find(trackPair.first);
0153     if (it == uniqueTrackIndex.end()) {
0154       result.push_back(trackPair.first);
0155       uniqueTrackIndex[trackPair.first] = i;
0156     }
0157
0158     it = uniqueTrackIndex.find(trackPair.second);
0159     if (it == uniqueTrackIndex.end()) {
0160       result.push_back(trackPair.second);
0161       uniqueTrackIndex[trackPair.second] = i;
0162     }
0163   }
0164
0165   return result;
0166 }
0167
0168 ///checks if the mass of the X is in the mass region adding missing E_T
0169 AlignmentTwoBodyDecayTrackSelector::Tracks AlignmentTwoBodyDecayTrackSelector::checkMETMass(
0170     const Tracks& cands, const edm::Event& iEvent) const {
0171   Tracks result;
0172
0173   LogDebug("Alignment") << ">  cands size : " << cands.size();
0174
0175   if (cands.empty())
0176     return result;
0177
0178   TLorentzVector track;
0179   TLorentzVector met4;
0180   TLorentzVector mother;
0181
0182   Handle<reco::CaloMETCollection> missingET;
0183   iEvent.getByToken(theMissingETToken, missingET);
0184   if (!missingET.isValid()) {
0185     LogError("Alignment") << "@SUB=AlignmentTwoBodyDecayTrackSelector::checkMETMass"
0186                           << ">  could not optain missingET Collection!";
0187     return result;
0188   }
0189
0190   typedef pair<double, const reco::Track*> candCollectionItem;
0191   vector<candCollectionItem> candCollection;
0192
0193   for (reco::CaloMETCollection::const_iterator itMET = missingET->begin(); itMET != missingET->end(); ++itMET) {
0194     met4.SetXYZT((*itMET).px(), (*itMET).py(), (*itMET).pz(), (*itMET).p());
0195
0196     for (unsigned int iCand = 0; iCand < cands.size(); iCand++) {
0197       track.SetXYZT(cands.at(iCand)->px(),
0198                     cands.at(iCand)->py(),
0199                     cands.at(iCand)->pz(),
0200                     sqrt(cands.at(iCand)->p() * cands.at(iCand)->p() + theDaughterMass * theDaughterMass));
0201
0202       mother = track + met4;
0203
0204       const reco::Track* trk = cands.at(iCand);
0205       const reco::CaloMET* met = &(*itMET);
0206
0207       bool correctCharge = true;
0208       if (theChargeSwitch)
0209         correctCharge = this->checkCharge(trk);
0210
0211       bool acoplanarTracks = true;
0212       if (theAcoplanarityFilterSwitch)
0213         acoplanarTracks = this->checkMETAcoplanarity(trk, met);
0214
0215       if (mother.M() > theMinMass && mother.M() < theMaxMass && correctCharge && acoplanarTracks) {
0216         candCollection.push_back(candCollectionItem(mother.Pt(), trk));
0217       }
0218     }
0219   }
0220
0221   if (candCollection.empty())
0222     return result;
0223
0224   sort(candCollection.begin(), candCollection.end(), [](auto& a, auto& b) { return a.first > b.first; });
0225
0226   std::map<const reco::Track*, unsigned int> uniqueTrackIndex;
0227   std::map<const reco::Track*, unsigned int>::iterator it;
0228   for (unsigned int i = 0; i < candCollection.size() && i < theCandNumber; i++) {
0229     it = uniqueTrackIndex.find(candCollection[i].second);
0230     if (it == uniqueTrackIndex.end()) {
0231       result.push_back(candCollection[i].second);
0232       uniqueTrackIndex[candCollection[i].second] = i;
0233     }
0234   }
0235
0236   return result;
0237 }
0238
0239 ///checks if the mother has charge = [theCharge]
0240 bool AlignmentTwoBodyDecayTrackSelector::checkCharge(const reco::Track* trk1, const reco::Track* trk2) const {
0241   int sumCharge = trk1->charge();
0242   if (trk2)
0243     sumCharge += trk2->charge();
0244   if (theUnsignedSwitch)
0245     sumCharge = std::abs(sumCharge);
0246   if (sumCharge == theCharge)
0247     return true;
0248   return false;
0249 }
0250
0251 ///checks if the [cands] are acoplanar (returns empty set if not)
0252 bool AlignmentTwoBodyDecayTrackSelector::checkAcoplanarity(const reco::Track* trk1, const reco::Track* trk2) const {
0253   if (fabs(deltaPhi(trk1->phi(), trk2->phi() - M_PI)) < theAcoplanarDistance)
0254     return true;
0255   return false;
0256 }
0257
0258 ///checks if the [cands] are acoplanar (returns empty set if not)
0259 bool AlignmentTwoBodyDecayTrackSelector::checkMETAcoplanarity(const reco::Track* trk1, const reco::CaloMET* met) const {
0260   if (fabs(deltaPhi(trk1->phi(), met->phi() - M_PI)) < theAcoplanarDistance)
0261     return true;
0262   return false;
0263 }
0264
0265 //===================HELPERS===================
0266
0267 ///print Information on Track-Collection
0268 void AlignmentTwoBodyDecayTrackSelector::printTracks(const Tracks& col) const {
0269   int count = 0;
0270   LogDebug("Alignment") << ">......................................";
0271   for (Tracks::const_iterator it = col.begin(); it < col.end(); ++it, ++count) {
0272     LogDebug("Alignment") << ">  Track No. " << count << ": p = (" << (*it)->px() << "," << (*it)->py() << ","
0273                           << (*it)->pz() << ")\n"
0274                           << ">                        pT = " << (*it)->pt() << " eta = " << (*it)->eta()
0275                           << " charge = " << (*it)->charge();
0276   }
0277   LogDebug("Alignment") << ">......................................";
0278 }