JetProducers/plugins/CSJetProducer.cc

0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "RecoJets/JetProducers/plugins/CSJetProducer.h"
0003
0004 #include <memory>
0005
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0008
0009 #include "Math/ProbFuncMathCore.h"
0010
0011 #include "fastjet/contrib/ConstituentSubtractor.hh"
0012
0013 using namespace std;
0014 using namespace reco;
0015 using namespace edm;
0016 using namespace cms;
0017
0018 CSJetProducer::CSJetProducer(edm::ParameterSet const& conf)
0019     : VirtualJetProducer(conf),
0020       csRParam_(-1.0),
0021       csAlpha_(0.),
0022       useModulatedRho_(conf.getParameter<bool>("useModulatedRho")),
0023       minFlowChi2Prob_(conf.getParameter<double>("minFlowChi2Prob")),
0024       maxFlowChi2Prob_(conf.getParameter<double>("maxFlowChi2Prob")) {
0025   //get eta range, rho and rhom map
0026   etaToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>("etaMap"));
0027   rhoToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>("rho"));
0028   rhomToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>("rhom"));
0029   csRParam_ = conf.getParameter<double>("csRParam");
0030   csAlpha_ = conf.getParameter<double>("csAlpha");
0031   if (useModulatedRho_)
0032     rhoFlowFitParamsToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>("rhoFlowFitParams"));
0033 }
0034
0035 void CSJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0036   // use the default production from one collection
0037   VirtualJetProducer::produce(iEvent, iSetup);
0038   //use runAlgorithm of this class
0039
0040   //Delete allocated memory. It is allocated every time runAlgorithm is called
0041   fjClusterSeq_.reset();
0042 }
0043
0044 //______________________________________________________________________________
0045 void CSJetProducer::runAlgorithm(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0046   // run algorithm
0047   if (!doAreaFastjet_ && !doRhoFastjet_) {
0048     fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
0049   } else if (voronoiRfact_ <= 0) {
0050     fjClusterSeq_ =
0051         ClusterSequencePtr(new fastjet::ClusterSequenceArea(fjInputs_, *fjJetDefinition_, *fjAreaDefinition_));
0052   } else {
0053     fjClusterSeq_ = ClusterSequencePtr(
0054         new fastjet::ClusterSequenceVoronoiArea(fjInputs_, *fjJetDefinition_, fastjet::VoronoiAreaSpec(voronoiRfact_)));
0055   }
0056
0057   fjJets_.clear();
0058   std::vector<fastjet::PseudoJet> tempJets = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
0059
0060   //Get local rho and rhom map
0061   edm::Handle<std::vector<double>> etaRanges;
0062   edm::Handle<std::vector<double>> rhoRanges;
0063   edm::Handle<std::vector<double>> rhomRanges;
0064   edm::Handle<std::vector<double>> rhoFlowFitParams;
0065   iEvent.getByToken(etaToken_, etaRanges);
0066   iEvent.getByToken(rhoToken_, rhoRanges);
0067   iEvent.getByToken(rhomToken_, rhomRanges);
0068   if (useModulatedRho_)
0069     iEvent.getByToken(rhoFlowFitParamsToken_, rhoFlowFitParams);
0070
0071   //Check if size of eta and background density vectors is the same
0072   unsigned int bkgVecSize = etaRanges->size();
0073   if (bkgVecSize < 1) {
0074     throw cms::Exception("WrongBkgInput") << "Producer needs vector with background estimates\n";
0075   }
0076   if (bkgVecSize != (rhoRanges->size() + 1) || bkgVecSize != (rhomRanges->size() + 1)) {
0077     throw cms::Exception("WrongBkgInput")
0078         << "Size of etaRanges (" << bkgVecSize << ") and rhoRanges (" << rhoRanges->size() << ") and/or rhomRanges ("
0079         << rhomRanges->size() << ") vectors inconsistent\n";
0080   }
0081   bool minProb = false;
0082   bool maxProb = false;
0083   if (useModulatedRho_) {
0084     if (!rhoFlowFitParams->empty()) {
0085       int chi2index = rhoFlowFitParams->size() - 3;
0086       double val = ROOT::Math::chisquared_cdf_c(rhoFlowFitParams->at(chi2index), rhoFlowFitParams->at(chi2index + 1));
0087       minProb = val > minFlowChi2Prob_;
0088       maxProb = val < maxFlowChi2Prob_;
0089     }
0090   }
0091
0092   //Allow the background densities to change within the jet
0093
0094   for (fastjet::PseudoJet& ijet : tempJets) {
0095     //----------------------------------------------------------------------
0096     // sift ghosts and particles in the input jet
0097     std::vector<fastjet::PseudoJet> particles, ghosts;
0098     fastjet::SelectorIsPureGhost().sift(ijet.constituents(), ghosts, particles);
0099     unsigned long nParticles = particles.size();
0100     if (nParticles == 0)
0101       continue;  //don't subtract ghost jets
0102
0103     //assign rho and rhom to ghosts according to local eta-dependent map + modulation as function of phi
0104     for (fastjet::PseudoJet& ighost : ghosts) {
0105       double rhoModulationFactor = 1.;
0106       double ghostPhi = ighost.phi_std();
0107
0108       if (useModulatedRho_) {
0109         if (!rhoFlowFitParams->empty()) {
0110           if (minProb && maxProb) {
0111             rhoModulationFactor = getModulatedRhoFactor(ghostPhi, rhoFlowFitParams);
0112           }
0113         }
0114       }
0115
0116       int32_t ghostPos = -1;
0117       auto const ighostEta = ighost.eta();
0118       if (ighostEta <= etaRanges->at(0) || bkgVecSize == 1) {
0119         ghostPos = 0;
0120       } else if (ighostEta >= etaRanges->at(bkgVecSize - 1)) {
0121         ghostPos = rhoRanges->size() - 1;
0122       } else {
0123         for (unsigned int ie = 0; ie < (bkgVecSize - 1); ie++) {
0124           if (ighostEta >= etaRanges->at(ie) && ighostEta < etaRanges->at(ie + 1)) {
0125             ghostPos = ie;
0126             break;
0127           }
0128         }
0129       }
0130       double pt = rhoRanges->at(ghostPos) * ighost.area() * rhoModulationFactor;
0131       double mass_squared =
0132           std::pow(rhoModulationFactor * rhomRanges->at(ghostPos) * ighost.area() + pt, 2) - std::pow(pt, 2);
0133       double mass = (mass_squared > 0) ? sqrt(mass_squared) : 0;
0134       ighost.reset_momentum_PtYPhiM(pt, ighost.rap(), ighost.phi(), mass);
0135     }
0136
0137     //----------------------------------------------------------------------
0138     //from here use official fastjet contrib package
0139
0140     fastjet::contrib::ConstituentSubtractor subtractor;
0141     subtractor.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR);  // distance in eta-phi plane
0142     subtractor.set_max_distance(
0143         csRParam_);  // free parameter for the maximal allowed distance between particle i and ghost k
0144     subtractor.set_alpha(
0145         csAlpha_);  // free parameter for the distance measure (the exponent of particle pt). Note that in older versions of the package alpha was multiplied by two but in newer versions this is not the case anymore
0146     subtractor.set_do_mass_subtraction();
0147     subtractor.set_remove_all_zero_pt_particles(true);
0148
0149     std::vector<fastjet::PseudoJet> subtracted_particles = subtractor.do_subtraction(particles, ghosts);
0150
0151     //Create subtracted jets
0152     fastjet::PseudoJet subtracted_jet = join(subtracted_particles);
0153     if (subtracted_jet.perp() > jetPtMin_)
0154       fjJets_.push_back(subtracted_jet);
0155   }
0156   fjJets_ = fastjet::sorted_by_pt(fjJets_);
0157 }
0158
0159 void CSJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0160   edm::ParameterSetDescription descCSJetProducer;
0161   ////// From CSJetProducer
0162   fillDescriptionsFromCSJetProducer(descCSJetProducer);
0163   ///// From VirtualJetProducer
0164   descCSJetProducer.add<string>("jetCollInstanceName", "");
0165   VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(descCSJetProducer);
0166   descCSJetProducer.add<bool>("sumRecHits", false);
0167
0168   /////////////////////
0169   descriptions.add("CSJetProducer", descCSJetProducer);
0170 }
0171
0172 void CSJetProducer::fillDescriptionsFromCSJetProducer(edm::ParameterSetDescription& desc) {
0173   desc.add<double>("csRParam", -1.);
0174   desc.add<double>("csAlpha", 2.);
0175   desc.add<bool>("useModulatedRho", false);
0176   desc.add<double>("minFlowChi2Prob", 0.05);
0177   desc.add<double>("maxFlowChi2Prob", 0.95);
0178   desc.add<edm::InputTag>("etaMap", {"hiFJRhoProducer", "mapEtaEdges"});
0179   desc.add<edm::InputTag>("rho", {"hiFJRhoProducer", "mapToRho"});
0180   desc.add<edm::InputTag>("rhom", {"hiFJRhoProducer", "mapToRhoM"});
0181   desc.add<edm::InputTag>("rhoFlowFitParams", {"hiFJRhoFlowModulationProducer", "rhoFlowFitParams"});
0182 }
0183
0184 double CSJetProducer::getModulatedRhoFactor(const double phi, const edm::Handle<std::vector<double>>& flowComponents) {
0185   // get the rho modulation as function of phi
0186   // flow modulation fit is done in RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer
0187   int nFlow = (flowComponents->size() - 4) / 2;
0188   double modulationValue = 1;
0189   double firstFlowComponent = flowComponents->at(nFlow * 2 + 3);
0190   for (int iFlow = 0; iFlow < nFlow; iFlow++) {
0191     modulationValue += 2.0 * flowComponents->at(2 * iFlow + 1) *
0192                        cos((iFlow + firstFlowComponent) * (phi - flowComponents->at(2 * iFlow + 2)));
0193   }
0194   return modulationValue;
0195 }
0196
0197 ////////////////////////////////////////////////////////////////////////////////
0198 // define as cmssw plugin
0199 ////////////////////////////////////////////////////////////////////////////////
0200
0201 DEFINE_FWK_MODULE(CSJetProducer);