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       double val = ROOT::Math::chisquared_cdf_c(rhoFlowFitParams->at(5), rhoFlowFitParams->at(6));
0086       minProb = val > minFlowChi2Prob_;
0087       maxProb = val < maxFlowChi2Prob_;
0088     }
0089   }
0090
0091   //Allow the background densities to change within the jet
0092
0093   for (fastjet::PseudoJet& ijet : tempJets) {
0094     //----------------------------------------------------------------------
0095     // sift ghosts and particles in the input jet
0096     std::vector<fastjet::PseudoJet> particles, ghosts;
0097     fastjet::SelectorIsPureGhost().sift(ijet.constituents(), ghosts, particles);
0098     unsigned long nParticles = particles.size();
0099     if (nParticles == 0)
0100       continue;  //don't subtract ghost jets
0101
0102     //assign rho and rhom to ghosts according to local eta-dependent map + modulation as function of phi
0103     for (fastjet::PseudoJet& ighost : ghosts) {
0104       double rhoModulationFactor = 1.;
0105       double ghostPhi = ighost.phi_std();
0106
0107       if (useModulatedRho_) {
0108         if (!rhoFlowFitParams->empty()) {
0109           if (minProb && maxProb) {
0110             rhoModulationFactor = getModulatedRhoFactor(ghostPhi,
0111                                                         rhoFlowFitParams->at(2),
0112                                                         rhoFlowFitParams->at(4),
0113                                                         rhoFlowFitParams->at(1),
0114                                                         rhoFlowFitParams->at(3));
0115           }
0116         }
0117       }
0118
0119       int32_t ghostPos = -1;
0120       auto const ighostEta = ighost.eta();
0121       if (ighostEta <= etaRanges->at(0) || bkgVecSize == 1) {
0122         ghostPos = 0;
0123       } else if (ighostEta >= etaRanges->at(bkgVecSize - 1)) {
0124         ghostPos = rhoRanges->size() - 1;
0125       } else {
0126         for (unsigned int ie = 0; ie < (bkgVecSize - 1); ie++) {
0127           if (ighostEta >= etaRanges->at(ie) && ighostEta < etaRanges->at(ie + 1)) {
0128             ghostPos = ie;
0129             break;
0130           }
0131         }
0132       }
0133       double pt = rhoRanges->at(ghostPos) * ighost.area() * rhoModulationFactor;
0134       double mass_squared =
0135           std::pow(rhoModulationFactor * rhomRanges->at(ghostPos) * ighost.area() + pt, 2) - std::pow(pt, 2);
0136       double mass = (mass_squared > 0) ? sqrt(mass_squared) : 0;
0137       ighost.reset_momentum_PtYPhiM(pt, ighost.rap(), ighost.phi(), mass);
0138     }
0139
0140     //----------------------------------------------------------------------
0141     //from here use official fastjet contrib package
0142
0143     fastjet::contrib::ConstituentSubtractor subtractor;
0144     subtractor.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR);  // distance in eta-phi plane
0145     subtractor.set_max_distance(
0146         csRParam_);  // free parameter for the maximal allowed distance between particle i and ghost k
0147     subtractor.set_alpha(
0148         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
0149     subtractor.set_do_mass_subtraction();
0150     subtractor.set_remove_all_zero_pt_particles(true);
0151
0152     std::vector<fastjet::PseudoJet> subtracted_particles = subtractor.do_subtraction(particles, ghosts);
0153
0154     //Create subtracted jets
0155     fastjet::PseudoJet subtracted_jet = join(subtracted_particles);
0156     if (subtracted_jet.perp() > 0.)
0157       fjJets_.push_back(subtracted_jet);
0158   }
0159   fjJets_ = fastjet::sorted_by_pt(fjJets_);
0160 }
0161
0162 void CSJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0163   edm::ParameterSetDescription descCSJetProducer;
0164   ////// From CSJetProducer
0165   fillDescriptionsFromCSJetProducer(descCSJetProducer);
0166   ///// From VirtualJetProducer
0167   descCSJetProducer.add<string>("jetCollInstanceName", "");
0168   VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(descCSJetProducer);
0169   descCSJetProducer.add<bool>("sumRecHits", false);
0170
0171   /////////////////////
0172   descriptions.add("CSJetProducer", descCSJetProducer);
0173 }
0174
0175 void CSJetProducer::fillDescriptionsFromCSJetProducer(edm::ParameterSetDescription& desc) {
0176   desc.add<double>("csRParam", -1.);
0177   desc.add<double>("csAlpha", 2.);
0178   desc.add<bool>("useModulatedRho", false);
0179   desc.add<double>("minFlowChi2Prob", 0.05);
0180   desc.add<double>("maxFlowChi2Prob", 0.95);
0181   desc.add<edm::InputTag>("etaMap", {"hiFJRhoProducer", "mapEtaEdges"});
0182   desc.add<edm::InputTag>("rho", {"hiFJRhoProducer", "mapToRho"});
0183   desc.add<edm::InputTag>("rhom", {"hiFJRhoProducer", "mapToRhoM"});
0184   desc.add<edm::InputTag>("rhoFlowFitParams", {"hiFJRhoFlowModulationProducer", "rhoFlowFitParams"});
0185 }
0186
0187 double CSJetProducer::getModulatedRhoFactor(
0188     const double phi, const double eventPlane2, const double eventPlane3, const double par1, const double par2) {
0189   // get the rho modulation as function of phi
0190   // flow modulation fit is done in RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer
0191   return 1. + 2. * par1 * cos(2. * (phi - eventPlane2)) + par2 * cos(3. * (phi - eventPlane3));
0192 }
0193
0194 ////////////////////////////////////////////////////////////////////////////////
0195 // define as cmssw plugin
0196 ////////////////////////////////////////////////////////////////////////////////
0197
0198 DEFINE_FWK_MODULE(CSJetProducer);