TagAndProbe/src/TagProbePairMaker.cc

0001 #include "PhysicsTools/TagAndProbe/interface/TagProbePairMaker.h"
0002 #include "DataFormats/Candidate/interface/Candidate.h"
0003
0004 tnp::TagProbePairMaker::TagProbePairMaker(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)
0005     : srcToken_(iC.consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("tagProbePairs"))),
0006       randGen_(nullptr) {
0007   std::string arbitration = iConfig.getParameter<std::string>("arbitration");
0008   if (arbitration == "None") {
0009     arbitration_ = None;
0010   } else if (arbitration == "OneProbe") {
0011     arbitration_ = OneProbe;
0012   } else if (arbitration == "OnePair") {
0013     arbitration_ = OnePair;
0014   } else if (arbitration == "NonDuplicate") {
0015     arbitration_ = NonDuplicate;
0016   } else if (arbitration == "BestMass") {
0017     arbitration_ = BestMass;
0018     arbitrationMass_ = iConfig.getParameter<double>("massForArbitration");
0019   } else if (arbitration == "Random2") {
0020     arbitration_ = Random2;
0021     randGen_ = new TRandom2(7777);
0022   } else if (arbitration == "HighestPt") {
0023     arbitration_ = HighestPt;
0024   } else
0025     throw cms::Exception("Configuration") << "TagProbePairMakerOnTheFly: the only currently "
0026                                           << "allowed values for 'arbitration' are "
0027                                           << "'None', 'OneProbe', 'BestMass', 'Random2'\n";
0028
0029   if (iConfig.existsAs<bool>("phiCutForTwoLeg")) {
0030     phiCutForTwoLeg_ = iConfig.getParameter<bool>("phiCutForTwoLeg");
0031     //std::cout << "Set phiCutForTwoLeg_ to " << phiCutForTwoLeg_ << std::endl;
0032   } else {
0033     phiCutForTwoLeg_ = false;
0034     //std::cout << "Set phiCutForTwoLeg_ to default " << phiCutForTwoLeg_ << std::endl;
0035   }
0036 }
0037
0038 tnp::TagProbePairs tnp::TagProbePairMaker::run(const edm::Event &iEvent) const {
0039   // declare output
0040   tnp::TagProbePairs pairs;
0041
0042   // read from event
0043   edm::Handle<reco::CandidateView> src;
0044   iEvent.getByToken(srcToken_, src);
0045
0046   // convert
0047   for (reco::CandidateView::const_iterator it = src->begin(), ed = src->end(); it != ed; ++it) {
0048     const reco::Candidate &mother = *it;
0049     if (mother.numberOfDaughters() != 2)
0050       throw cms::Exception("CorruptData") << "Tag&Probe pair with " << mother.numberOfDaughters() << " daughters\n";
0051     pairs.push_back(tnp::TagProbePair(mother.daughter(0)->masterClone(),
0052                                       mother.daughter(1)->masterClone(),
0053                                       src->refAt(it - src->begin()),
0054                                       mother.mass()));
0055   }
0056
0057   if ((arbitration_ != None) && (pairs.size() > 1)) {
0058     // might need to clean up
0059     arbitrate(pairs);
0060   }
0061
0062   if (phiCutForTwoLeg_ && !pairs.empty()) {
0063     int eventNum = iEvent.id().event();
0064     std::cout << "Calling phiCutByEventNumber on eventNum=" << eventNum << std::endl;
0065     phiCutByEventNumber(pairs, eventNum);
0066   }
0067
0068   // return
0069   return pairs;
0070 }
0071
0072 void tnp::TagProbePairMaker::phiCutByEventNumber(TagProbePairs &pairs, int eventNumber) const {
0073   unsigned int currentNum = 0;
0074
0075   size_t nclean = pairs.size();
0076   for (TagProbePairs::iterator it = pairs.begin(), ed = pairs.end(); it != ed; ++it) {
0077     if (it->tag.isNull())
0078       continue;  // skip already invalidated pairs
0079     if (eventNumber % 2) {
0080       std::cout << "Odd event number " << eventNumber << ", require 0 < phi(tag) < pi... ";
0081       if (!(it->tag->phi() > 0. && it->tag->phi() < 3.141592654)) {
0082         std::cout << "Rejecting pair number " << currentNum++ << " with tag phi " << it->tag->phi();
0083         nclean--;
0084         it->tag = reco::CandidateBaseRef();
0085         --nclean;
0086       } else {
0087         std::cout << "Keeping pair number " << currentNum++ << " with tag phi " << it->tag->phi();
0088       }
0089     } else {
0090       std::cout << "Even event number " << eventNumber << ", require -pi < phi(tag) < 0... ";
0091       //      if (!(it->tag->phi() > 3.141592654 && it->tag->phi() < 2*3.141592654)) {
0092       if (!(it->tag->phi() > -3.141592654 && it->tag->phi() < 0)) {
0093         std::cout << "Rejecting pair number " << currentNum++ << " with tag phi " << it->tag->phi();
0094         nclean--;
0095         it->tag = reco::CandidateBaseRef();
0096         --nclean;
0097       } else {
0098         std::cout << "Keeping pair number " << currentNum++ << " with tag phi " << it->tag->phi();
0099       }
0100     }
0101     std::cout << std::endl;
0102   }
0103
0104   if (nclean == 0) {
0105     pairs.clear();
0106   } else if (nclean < pairs.size()) {
0107     TagProbePairs cleaned;
0108     cleaned.reserve(nclean);
0109     for (TagProbePairs::iterator it = pairs.begin(), ed = pairs.end(); it != ed; ++it) {
0110       if (it->tag.isNonnull())
0111         cleaned.push_back(*it);
0112     }
0113     pairs.swap(cleaned);
0114   }
0115 }
0116
0117 void tnp::TagProbePairMaker::arbitrate(TagProbePairs &pairs) const {
0118   size_t nclean = pairs.size();
0119   for (TagProbePairs::iterator it = pairs.begin(), ed = pairs.end(); it != ed; ++it) {
0120     if (it->tag.isNull())
0121       continue;  // skip already invalidated pairs
0122
0123     bool TTpair = false;
0124     for (TagProbePairs::iterator it2 = pairs.begin(); it2 != ed; ++it2) {  // first check for Tag-Tag pairs
0125       if (it2->tag.isNull())
0126         continue;  // skip already invalidated pairs
0127       if (it != it2 && it->probe == it2->tag && it->tag == it2->probe) {
0128         //std::cout << "----------> probe is tag!!!!" << std::endl;
0129         TTpair = true;
0130       }
0131     }
0132     //if(!TTpair) std::cout << "this is not TT!" << std::endl;
0133     bool invalidateThis = false;
0134     int numberOfProbes = 0;
0135     for (TagProbePairs::iterator it2 = it + 1; it2 != ed; ++it2) {  // it+1 <= ed, otherwise we don't arrive here
0136       // arbitrate the case where multiple probes are matched to the same tag
0137       if ((arbitration_ != NonDuplicate) && (it->tag == it2->tag)) {
0138         if (TTpair) {  // we already have found a TT pair, no need to guess
0139           it2->tag = reco::CandidateBaseRef();
0140           --nclean;
0141           //std::cout << "remove unnecessary pair! -----------" << std::endl;
0142           continue;
0143         }
0144
0145         if (arbitration_ == OneProbe) {
0146           // invalidate this one
0147           it2->tag = reco::CandidateBaseRef();
0148           --nclean;
0149           // remember to invalidate also the other one (but after finishing to check for duplicates)
0150           invalidateThis = true;
0151         } else if (arbitration_ == BestMass) {
0152           // but the best one in the first  iterator
0153           if (fabs(it2->mass - arbitrationMass_) < fabs(it->mass - arbitrationMass_)) {
0154             std::swap(*it, *it2);
0155           }
0156           // and invalidate it2
0157           it2->tag = reco::CandidateBaseRef();
0158           --nclean;
0159         } else if (arbitration_ == HighestPt) {
0160           // but the best one in the first  iterator
0161           if (it2->probe->pt() > it->probe->pt()) {
0162             std::swap(*it, *it2);
0163           }
0164           // and invalidate it2
0165           it2->tag = reco::CandidateBaseRef();
0166           --nclean;
0167         } else if (arbitration_ == Random2) {
0168           numberOfProbes++;
0169           if (numberOfProbes > 1) {
0170             //std::cout << "more than 2 probes!" << std::endl;
0171             invalidateThis = true;
0172             it2->tag = reco::CandidateBaseRef();
0173             --nclean;
0174           } else {
0175             // do a coin toss to decide if we want to swap them
0176             if (randGen_->Rndm() > 0.5) {
0177               std::swap(*it, *it2);
0178             }
0179             // and invalidate it2
0180             it2->tag = reco::CandidateBaseRef();
0181             --nclean;
0182           }
0183         }
0184       }
0185       // arbitrate the case where the same probe is associated to more then one tag
0186       if ((arbitration_ == NonDuplicate) && (it->probe == it2->probe)) {
0187         // invalidate the pair in which the tag has lower pT
0188         if (it2->tag->pt() > it->tag->pt())
0189           std::swap(*it, *it2);
0190         it2->tag = reco::CandidateBaseRef();
0191         --nclean;
0192       }
0193       // arbitrate the OnePair case: disallow the same pair to enter the lineshape twice
0194       // this can't be done by reference, unfortunately, so we resort to a simple matching
0195       if ((arbitration_ == OnePair) && std::abs(it->mass - it2->mass) < 1e-4 &&
0196           std::abs(it->probe->phi() - it2->tag->phi()) < 1e-5 && std::abs(it->probe->eta() - it2->tag->eta()) < 1e-5 &&
0197           std::abs(it->probe->pt() - it2->tag->pt()) < 1e-5 && std::abs(it2->probe->phi() - it->tag->phi()) < 1e-5 &&
0198           std::abs(it2->probe->eta() - it->tag->eta()) < 1e-5 && std::abs(it2->probe->pt() - it->tag->pt()) < 1e-5) {
0199         it2->tag = reco::CandidateBaseRef();
0200         --nclean;
0201       }
0202     }
0203     if (invalidateThis) {
0204       it->tag = reco::CandidateBaseRef();
0205       --nclean;
0206     }
0207   }
0208
0209   if (nclean == 0) {
0210     pairs.clear();
0211   } else if (nclean < pairs.size()) {
0212     TagProbePairs cleaned;
0213     cleaned.reserve(nclean);
0214     for (TagProbePairs::iterator it = pairs.begin(), ed = pairs.end(); it != ed; ++it) {
0215       if (it->tag.isNonnull())
0216         cleaned.push_back(*it);
0217     }
0218     pairs.swap(cleaned);
0219   }
0220 }