text: #include "FWCore/Framework/interface/MakerMacros.h"
  #include "RecoJets/JetProducers/plugins/CSJetProducer.h"
  #include "FWCore/Utilities/interface/Exception.h"
  #include "RecoJets/JetProducers/interface/JetSpecific.h"
  
  #include "Math/ProbFuncMathCore.h"
  
  #include <iostream>
  
  #include "fastjet/contrib/ConstituentSubtractor.hh"
  
  using namespace std;
  using namespace reco;
  using namespace edm;
  using namespace cms;
  
  CSJetProducer::CSJetProducer(edm::ParameterSet const& conf):
    VirtualJetProducer( conf ),
    csRParam_(-1.0),
    csAlpha_(0.),
    useModulatedRho_(false)
  {
    //get eta range, rho and rhom map
    etaToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>( "etaMap" ));
    rhoToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>( "rho" ));
    rhomToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>( "rhom" ));
    csRParam_ = conf.getParameter<double>("csRParam");
    csAlpha_ = conf.getParameter<double>("csAlpha");
    useModulatedRho_ = conf.getParameter<bool>("useModulatedRho");
    minFlowChi2Prob_ = conf.getParameter<double>("minFlowChi2Prob");
    maxFlowChi2Prob_ = conf.getParameter<double>("maxFlowChi2Prob");
    if(useModulatedRho_) rhoFlowFitParamsToken_ = consumes<std::vector<double>>(conf.getParameter<edm::InputTag>( "rhoFlowFitParams" ));
  }
  
  void CSJetProducer::produce( edm::Event & iEvent, const edm::EventSetup & iSetup )
  {
    // use the default production from one collection
    VirtualJetProducer::produce( iEvent, iSetup );
    //use runAlgorithm of this class
  
    //Delete allocated memory. It is allocated every time runAlgorithm is called
    fjClusterSeq_.reset();
  }
  
  //______________________________________________________________________________
  void CSJetProducer::runAlgorithm( edm::Event & iEvent, edm::EventSetup const& iSetup)
  {
    // run algorithm
    if ( !doAreaFastjet_ && !doRhoFastjet_) {
      fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs_, *fjJetDefinition_ ) );
    } else if (voronoiRfact_ <= 0) {
      fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequenceArea( fjInputs_, *fjJetDefinition_ , *fjAreaDefinition_ ) );
    } else {
      fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequenceVoronoiArea( fjInputs_, *fjJetDefinition_ , fastjet::VoronoiAreaSpec(voronoiRfact_) ) );
    }
  
    fjJets_.clear();
    std::vector<fastjet::PseudoJet> tempJets = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
  
    //Get local rho and rhom map
    edm::Handle<std::vector<double>> etaRanges;
    edm::Handle<std::vector<double>> rhoRanges;
    edm::Handle<std::vector<double>> rhomRanges;
    edm::Handle<std::vector<double>> rhoFlowFitParams;
  
    iEvent.getByToken(etaToken_, etaRanges);
    iEvent.getByToken(rhoToken_, rhoRanges);
    iEvent.getByToken(rhomToken_, rhomRanges);
    if(useModulatedRho_) iEvent.getByToken(rhoFlowFitParamsToken_, rhoFlowFitParams);
  
    //Check if size of eta and background density vectors is the same
    unsigned int bkgVecSize = etaRanges->size();
    if(bkgVecSize<1) { throw cms::Exception("WrongBkgInput") << "Producer needs vector with background estimates\n"; }
    if(bkgVecSize != (rhoRanges->size()+1) || bkgVecSize != (rhomRanges->size()+1) ) {
      throw cms::Exception("WrongBkgInput") << "Size of etaRanges (" << bkgVecSize << ") and rhoRanges (" << rhoRanges->size() << ") and/or rhomRanges (" << rhomRanges->size() << ") vectors inconsistent\n";
    }
  
    
    //Allow the background densities to change within the jet
  
    for(fastjet::PseudoJet& ijet : tempJets ) {
    
      //----------------------------------------------------------------------
      // sift ghosts and particles in the input jet
      std::vector<fastjet::PseudoJet> particles, ghosts;
      fastjet::SelectorIsPureGhost().sift(ijet.constituents(), ghosts, particles);
      unsigned long nParticles=particles.size();
      if(nParticles==0) continue; //don't subtract ghost jets
      
      //assign rho and rhom to ghosts according to local eta-dependent map + modulation as function of phi
      for(fastjet::PseudoJet& ighost : ghosts) {
        double rhoModulationFactor = 1.;
        double ghostPhi = ighost.phi_std();
  
        if(useModulatedRho_){
  	if(rhoFlowFitParams->size() > 0){
  	  double val = ROOT::Math::chisquared_cdf_c(rhoFlowFitParams->at(5), rhoFlowFitParams->at(6));
  	  bool minProb = val > minFlowChi2Prob_;
  	  bool maxProb = val < maxFlowChi2Prob_;
  
  	  if(minProb && maxProb)
  	    rhoModulationFactor = getModulatedRhoFactor(ghostPhi,
  							rhoFlowFitParams->at(2),
  							rhoFlowFitParams->at(4),
  							rhoFlowFitParams->at(1),
  							rhoFlowFitParams->at(3)
  							);
  	}
        }
  
        int ghostPos = -1;
        if(ighost.eta()<=etaRanges->at(0)) ghostPos = 0;
        else if(ighost.eta()>=etaRanges->at(etaRanges->size()-1)) ghostPos = rhoRanges->size()-1;
        else{
          for(unsigned int ie = 0; ie < etaRanges->size()-1; ++ie){
            if(ighost.eta()>=etaRanges->at(ie) && ighost.eta()<etaRanges->at(ie+1)){
              ghostPos = ie;
              break;
            }
          }
        }
  
  
        double pt = rhoRanges->at(ghostPos)*ighost.area()*rhoModulationFactor;
        double mass_squared=std::pow(rhoModulationFactor*rhomRanges->at(ghostPos)*ighost.area()+pt,2)-std::pow(pt,2);
        double mass=0;
        if (mass_squared>0) mass=sqrt(mass_squared);
        ighost.reset_momentum_PtYPhiM(pt,ighost.rap(),ighost.phi(),mass);
      }
  
      //----------------------------------------------------------------------
      //from here use official fastjet contrib package
      
      fastjet::contrib::ConstituentSubtractor subtractor;
      subtractor.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR); // distance in eta-phi plane
      subtractor.set_max_distance(csRParam_); // free parameter for the maximal allowed distance between particle i and ghost k
      subtractor.set_alpha(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
      subtractor.set_do_mass_subtraction(true);
      subtractor.set_remove_all_zero_pt_particles(true);
  
      std::vector<fastjet::PseudoJet> subtracted_particles = subtractor.do_subtraction(particles,ghosts);
  
      //Create subtracted jets
      fastjet::PseudoJet subtracted_jet=join(subtracted_particles);
      if(subtracted_jet.perp()>0.)
        fjJets_.push_back( subtracted_jet );
    }
    fjJets_ = fastjet::sorted_by_pt(fjJets_); 
  }
  
  void CSJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  
    edm::ParameterSetDescription descCSJetProducer;
    ////// From CSJetProducer
    fillDescriptionsFromCSJetProducer(descCSJetProducer);
    ///// From VirtualJetProducer
    descCSJetProducer.add<string>("jetCollInstanceName", ""    );
    VirtualJetProducer::fillDescriptionsFromVirtualJetProducer(descCSJetProducer);
    descCSJetProducer.add<bool> ("sumRecHits", false);
    
    /////////////////////
    descriptions.add("CSJetProducer",descCSJetProducer);
    
  }
  
  void CSJetProducer::fillDescriptionsFromCSJetProducer(edm::ParameterSetDescription& desc) {
  
    desc.add<double>("csRParam",-1.);
    desc.add<double>("csAlpha",2.);
    desc.add<bool> ("useModulatedRho", false);
    desc.add<double> ("minFlowChi2Prob", false);
    desc.add<double> ("maxFlowChi2Prob", false);
  
    desc.add<edm::InputTag>("etaMap",edm::InputTag("hiFJRhoProducer","mapEtaEdges") );
    desc.add<edm::InputTag>("rho",edm::InputTag("hiFJRhoProducer","mapToRho") );
    desc.add<edm::InputTag>("rhom",edm::InputTag("hiFJRhoProducer","mapToRhoM") );
  
    desc.add<edm::InputTag>("rhoFlowFitParams", edm::InputTag("hiFJRhoFlowModulationProducer","rhoFlowFitParams") );
  }
  
  double CSJetProducer::getModulatedRhoFactor(const double phi, const double eventPlane2, const double eventPlane3, const double par1, const double par2) {
    //get the rho modulation as function of phi
    //flow modulation fit is done in HiJetBackground/HiFJRhoFlowModulationProducer
    double mod = 1. + 2.*(par1*cos(2.*(phi - eventPlane2))) + par2*cos(3.*(phi - eventPlane3));
  
    return mod;
  }
  
  ////////////////////////////////////////////////////////////////////////////////
  // define as cmssw plugin
  ////////////////////////////////////////////////////////////////////////////////
  
  DEFINE_FWK_MODULE(CSJetProducer);
  