HiJetAlgos/plugins/HiFJRhoFlowModulationProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer
0004 // Class:      HiFJRhoFlowModulationProducer
0005
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/Common/interface/View.h"
0008 #include "DataFormats/JetReco/interface/Jet.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0011 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0012 #include "DataFormats/Provenance/interface/EventID.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
0019 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
0020 #include "DataFormats/JetReco/interface/JetCollection.h"
0021 #include "FWCore/Framework/interface/stream/EDProducer.h"
0022 #include "FWCore/Utilities/interface/StreamID.h"
0023 #include "FWCore/Utilities/interface/EDGetToken.h"
0024 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h"
0025
0026 #include "TF1.h"
0027 #include "TF2.h"
0028 #include "TH1.h"
0029 #include "TMath.h"
0030 #include "Math/ProbFuncMathCore.h"
0031 #include "TMinuitMinimizer.h"
0032
0033 #include <algorithm>
0034 #include <cmath>
0035 #include <memory>
0036 #include <string>
0037 #include <utility>
0038 #include <vector>
0039
0040 namespace {
0041
0042   double flowFunction(double* x, double* par) {
0043     unsigned int nFlow = par[0];          // Number of fitted flow components is defined in the first parameter
0044     unsigned int firstFittedVn = par[1];  // The first fitted flow component is defined in the second parameter
0045
0046     // Add each component separately to the total fit value
0047     double flowModulation = par[2];
0048     for (unsigned int iFlow = 0; iFlow < nFlow; iFlow++) {
0049       flowModulation +=
0050           par[2] * 2.0 * par[2 * iFlow + 3] * std::cos((iFlow + firstFittedVn) * (x[0] - par[2 * iFlow + 4]));
0051     }
0052     return flowModulation;
0053   }
0054
0055   /*
0056    * Weighing function for particle flow candidates in case jetty regions are exluded from the flow fit
0057    *   x[0] = DeltaPhi between jet axis and particle flow candidate in range [0,2Pi]
0058    *   x[1] = Absolute value of the jet eta
0059    *   par[0] = Eta cut applied for the particle flow candidates
0060    *   par[1] = Exclusion radius around the jet axis
0061    *
0062    *   return: The weight to be used with this particle flow candidate is (1 + number provided by this function)
0063    */
0064   double weightFunction(double* x, double* par) {
0065     // If the particle flow candidate is farther than the exclusion radius from the jet axis, no need for weighting due to area where no particle flow candidates can be found
0066     if (x[0] > par[1])
0067       return 0;
0068
0069     // If the particle flow candidate is closer than 0.4 in phi from the jet axis, then there is part of phi acceptance from which no particle flow candidates can be found. Calculate the half of the size of that strip in eta
0070     double exclusionArea = std::sqrt(par[1] * par[1] - x[0] * x[0]);
0071
0072     // Particle flow candidates are only considered until eta = par[0]. Check that the exclusion area does not go beyond that
0073
0074     // First, check cases where the jet is found outside of the particle flow candidate acceptance in eta
0075     if (x[1] > par[0]) {
0076       // Check if the whole exclusion area is outside of the acceptance. If this is the case, we should not add anything to the weight for this particle flow candidate.
0077       if (x[1] - exclusionArea > par[0])
0078         return 0;
0079
0080       // Now we know that part of the exclusion area will be inside of the acceptance. Check how big this is.
0081       exclusionArea += par[0] - x[1];
0082
0083     } else {
0084       // In the next case, the jet is found inside of the particle flow candidate acceptance. In this case, we need to check if some of the exclusion area goes outside of the acceptance. If is does, we need to exclude that from the exclusion area
0085       if (x[1] + exclusionArea > par[0]) {
0086         exclusionArea += par[0] - x[1];
0087       } else {
0088         // If not, the the exclusion area is two times half of the nominal exclusion area
0089         exclusionArea *= 2;
0090       }
0091     }
0092
0093     // Normalize the exclusion area to the total acceptance. This number should be added to the total weight for the particle flow candidate
0094     return (2 * par[0]) / (2 * par[0] - exclusionArea) - 1;
0095   }
0096 };  // namespace
0097
0098 class HiFJRhoFlowModulationProducer : public edm::stream::EDProducer<> {
0099 public:
0100   explicit HiFJRhoFlowModulationProducer(const edm::ParameterSet&);
0101   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0102
0103 private:
0104   void beginStream(edm::StreamID) override;
0105   void produce(edm::Event&, const edm::EventSetup&) override;
0106   void endStream() override;
0107
0108   const int minPfCandidatesPerEvent_;
0109   const bool doEvtPlane_;
0110   const bool doFreePlaneFit_;
0111   const bool doJettyExclusion_;
0112   const double exclusionRadius_;
0113   const int evtPlaneLevel_;
0114   const double pfCandidateEtaCut_;
0115   const double pfCandidateMinPtCut_;
0116   const double pfCandidateMaxPtCut_;
0117   int firstFittedVn_;
0118   int lastFittedVn_;
0119   const edm::EDGetTokenT<reco::JetView> jetToken_;
0120   const edm::EDGetTokenT<pat::PackedCandidateCollection> pfCandidateToken_;
0121   const edm::EDGetTokenT<reco::EvtPlaneCollection> evtPlaneToken_;
0122   std::unique_ptr<TF1> flowFit_p_;
0123   std::unique_ptr<TF2> pfWeightFunction_;
0124   reco::PFCandidate converter_;
0125 };
0126 HiFJRhoFlowModulationProducer::HiFJRhoFlowModulationProducer(const edm::ParameterSet& iConfig)
0127     : minPfCandidatesPerEvent_(iConfig.getParameter<int>("minPfCandidatesPerEvent")),
0128       doEvtPlane_(iConfig.getParameter<bool>("doEvtPlane")),
0129       doFreePlaneFit_(iConfig.getParameter<bool>("doFreePlaneFit")),
0130       doJettyExclusion_(iConfig.getParameter<bool>("doJettyExclusion")),
0131       exclusionRadius_(iConfig.getParameter<double>("exclusionRadius")),
0132       evtPlaneLevel_(iConfig.getParameter<int>("evtPlaneLevel")),
0133       pfCandidateEtaCut_(iConfig.getParameter<double>("pfCandidateEtaCut")),
0134       pfCandidateMinPtCut_(iConfig.getParameter<double>("pfCandidateMinPtCut")),
0135       pfCandidateMaxPtCut_(iConfig.getParameter<double>("pfCandidateMaxPtCut")),
0136       firstFittedVn_(iConfig.getParameter<int>("firstFittedVn")),
0137       lastFittedVn_(iConfig.getParameter<int>("lastFittedVn")),
0138       jetToken_(doJettyExclusion_ ? consumes<reco::JetView>(iConfig.getParameter<edm::InputTag>("jetTag"))
0139                                   : edm::EDGetTokenT<reco::JetView>()),
0140       pfCandidateToken_(consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"))),
0141       evtPlaneToken_(consumes<reco::EvtPlaneCollection>(iConfig.getParameter<edm::InputTag>("EvtPlane"))) {
0142   produces<std::vector<double>>("rhoFlowFitParams");
0143
0144   // Converter to get PF candidate ID from packed candidates
0145   converter_ = reco::PFCandidate();
0146
0147   // Sanity check for the fitted flow component input
0148   if (firstFittedVn_ < 1)
0149     firstFittedVn_ = 2;
0150   if (lastFittedVn_ < 1)
0151     lastFittedVn_ = 2;
0152   if (firstFittedVn_ > lastFittedVn_) {
0153     int flipper = lastFittedVn_;
0154     lastFittedVn_ = firstFittedVn_;
0155     firstFittedVn_ = flipper;
0156   }
0157
0158   // Calculate the number of flow components
0159   const int nFlow = lastFittedVn_ - firstFittedVn_ + 1;
0160
0161   // Define all fit and weight functions
0162   TMinuitMinimizer::UseStaticMinuit(false);
0163   flowFit_p_ = std::make_unique<TF1>("flowFit", flowFunction, -TMath::Pi(), TMath::Pi(), nFlow * 2 + 3);
0164   flowFit_p_->FixParameter(0, nFlow);           // The first parameter defines the number of fitted flow components
0165   flowFit_p_->FixParameter(1, firstFittedVn_);  // The second parameter defines the first fitted flow component
0166   pfWeightFunction_ = std::make_unique<TF2>("weightFunction", weightFunction, 0, TMath::Pi(), -5, 5, 2);
0167   pfWeightFunction_->SetParameter(0, pfCandidateEtaCut_);  // Set the allowed eta range for particle flow candidates
0168   pfWeightFunction_->SetParameter(1, exclusionRadius_);    // Set the exclusion radius around the jets
0169 }
0170
0171 // ------------ method called to produce the data  ------------
0172 void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0173   // Get the particle flow candidate collection
0174   auto const& pfCands = iEvent.get(pfCandidateToken_);
0175
0176   // If we are reading the event plane information for the forest, find the event planes
0177   std::array<float, hi::NumEPNames> hiEvtPlane;
0178   if (doEvtPlane_) {
0179     auto const& evtPlanes = iEvent.get(evtPlaneToken_);
0180     assert(evtPlanes.size() == hi::NumEPNames);
0181     std::transform(evtPlanes.begin(), evtPlanes.end(), hiEvtPlane.begin(), [this](auto const& ePlane) -> float {
0182       return ePlane.angle(evtPlaneLevel_);
0183     });
0184   }
0185
0186   // If jetty regions are excluded from the flow fits, read the jet collection
0187   edm::Handle<reco::JetView> jets;
0188   if (doJettyExclusion_) {
0189     iEvent.getByToken(jetToken_, jets);
0190   }
0191
0192   int nFill = 0;
0193
0194   // Initialize arrays for event planes
0195   // nFlow: Number of flow components in the fit to particle flow condidate distribution
0196   const int nFlow = lastFittedVn_ - firstFittedVn_ + 1;
0197   double eventPlane[nFlow];
0198   for (int iFlow = 0; iFlow < nFlow; iFlow++)
0199     eventPlane[iFlow] = -100;
0200
0201   // Initialize the output vector with flow fit components
0202   // v_n and Psi_n for each fitted flow component, plus overall normalization factor, chi^2 and ndf for flow fit, and the first fitted flow component
0203   const int nParamVals = nFlow * 2 + 4;
0204   auto rhoFlowFitParamsOut = std::make_unique<std::vector<double>>(nParamVals, 1e-6);
0205
0206   // Set the parameters related to flow fit to zero
0207   for (int iParameter = 0; iParameter < nFlow * 2 + 1; iParameter++) {
0208     rhoFlowFitParamsOut->at(iParameter) = 0;
0209   }
0210   // Set the first fitted flow component to the last index of the output vector
0211   rhoFlowFitParamsOut->at(nFlow * 2 + 3) = firstFittedVn_;
0212
0213   // Initialize arrays for manual event plane calculation
0214   double eventPlaneCos[nFlow];
0215   double eventPlaneSin[nFlow];
0216   for (int iFlow = 0; iFlow < nFlow; iFlow++) {
0217     eventPlaneCos[iFlow] = 0;
0218     eventPlaneSin[iFlow] = 0;
0219   }
0220
0221   // Define helper variables for looping over particle flow candidates
0222   std::vector<bool> pfcuts(pfCands.size(), false);
0223   int iCand = -1;
0224   double thisPfCandidateWeight = 0;
0225   double deltaPhiJetPf = 0;
0226   std::vector<double> pfCandidateWeight(pfCands.size(), 1);
0227
0228   // Loop over the particle flow candidates
0229   for (auto const& pfCandidate : pfCands) {
0230     iCand++;  // Update the particle flow candidate index
0231
0232     // Find the PF candidate ID from the packed candidates collection
0233     auto particleFlowCandidateID = converter_.translatePdgIdToType(pfCandidate.pdgId());
0234
0235     // This cut selects particle flow candidates that are charged hadrons. The ID numbers are:
0236     // 0 = undefined
0237     // 1 = charged hadron
0238     // 2 = electron
0239     // 3 = muon
0240     // 4 = photon
0241     // 5 = neutral hadron
0242     // 6 = HF tower identified as a hadron
0243     // 7 = HF tower identified as an EM particle
0244     if (particleFlowCandidateID != 1)
0245       continue;
0246
0247     // Kinematic cuts for particle flow candidate pT and eta
0248     if (std::abs(pfCandidate.eta()) > pfCandidateEtaCut_)
0249       continue;
0250     if (pfCandidate.pt() < pfCandidateMinPtCut_)
0251       continue;
0252     if (pfCandidate.pt() > pfCandidateMaxPtCut_)
0253       continue;
0254
0255     nFill++;  // Use same nFill criterion with and without jetty region subtraction
0256
0257     // If the jetty regions are excluded from the fit, find which particle flow candidates are close to jets
0258     thisPfCandidateWeight = 1;
0259     if (doJettyExclusion_) {
0260       bool isGood = true;
0261       for (auto const& jet : *jets) {
0262         if (deltaR2(jet, pfCandidate) < exclusionRadius_ * exclusionRadius_) {
0263           isGood = false;
0264           break;
0265         } else {
0266           // If the particle flow candidates are not excluded from the fit, check if there are any jets in the same phi-slice as the particle flow candidate. If this is the case, add a weight for the particle flow candidate taking into account the lost acceptance for underlaying event particle flow candidates.
0267           deltaPhiJetPf = std::abs(pfCandidate.phi() - jet.phi());
0268           if (deltaPhiJetPf > TMath::Pi())
0269             deltaPhiJetPf = TMath::Pi() * 2 - deltaPhiJetPf;
0270           // Weight currently does not take into account overlapping jets
0271           thisPfCandidateWeight += pfWeightFunction_->Eval(deltaPhiJetPf, std::abs(jet.eta()));
0272         }
0273       }
0274       // Do not use this particle flow candidate in the manual event plane calculation
0275       if (!isGood) {
0276         continue;
0277       }
0278     }
0279
0280     // Update the weight for this particle flow candidate
0281     pfCandidateWeight[iCand] = thisPfCandidateWeight;
0282
0283     // This particle flow candidate passes all the cuts
0284     pfcuts[iCand] = true;
0285
0286     // If the event plane is calculated manually, add this particle flow candidate to the calculation
0287     if (!doEvtPlane_) {
0288       for (int iFlow = 0; iFlow < nFlow; iFlow++) {
0289         eventPlaneCos[iFlow] += std::cos((iFlow + firstFittedVn_) * pfCandidate.phi()) * thisPfCandidateWeight;
0290         eventPlaneSin[iFlow] += std::sin((iFlow + firstFittedVn_) * pfCandidate.phi()) * thisPfCandidateWeight;
0291       }
0292     }
0293   }
0294
0295   // Determine the event plane angle
0296   if (!doEvtPlane_) {
0297     // Manual calculation for the event plane angle using the particle flow candidates
0298     for (int iFlow = 0; iFlow < nFlow; iFlow++) {
0299       eventPlane[iFlow] = std::atan2(eventPlaneSin[iFlow], eventPlaneCos[iFlow]) / (iFlow + firstFittedVn_);
0300     }
0301   } else {
0302     // Read the event plane angle determined from the HF calorimeters from the HiForest
0303     // Only v2 and v3 are available in the HiForest. This option should not be used if other flow components are fitted
0304     int halfWay = nFlow / 2;
0305     for (int iFlow = 0; iFlow < halfWay; iFlow++) {
0306       eventPlane[iFlow] = hiEvtPlane[hi::HF2];
0307     }
0308     for (int iFlow = halfWay; iFlow < nFlow; iFlow++) {
0309       eventPlane[iFlow] = hiEvtPlane[hi::HF3];
0310     }
0311   }
0312
0313   // Do the flow fit provided that there are enough particle flow candidates in the event and that the event planes have been determined successfully
0314   int pfcuts_count = 0;
0315   if (nFill >= minPfCandidatesPerEvent_ && eventPlane[0] > -99) {
0316     // Create a particle flow candidate phi-histogram
0317     const int nPhiBins = std::max(10, nFill / 30);
0318     std::string name = "phiTestIEta4_" + std::to_string(iEvent.id().event()) + "_h";
0319     std::unique_ptr<TH1F> phi_h = std::make_unique<TH1F>(name.data(), "", nPhiBins, -TMath::Pi(), TMath::Pi());
0320     phi_h->SetDirectory(nullptr);
0321     for (auto const& pfCandidate : pfCands) {
0322       if (pfcuts.at(pfcuts_count)) {  // Only use particle flow candidates that pass all cuts
0323         phi_h->Fill(pfCandidate.phi(), pfCandidateWeight.at(pfcuts_count));
0324       }
0325       pfcuts_count++;
0326     }
0327
0328     // Set initial values for the fit parameters
0329     flowFit_p_->SetParameter(2, 10);
0330     for (int iFlow = 0; iFlow < nFlow; iFlow++) {
0331       flowFit_p_->SetParameter(3 + 2 * iFlow, 0);
0332       flowFit_p_->SetParameter(4 + 2 * iFlow, eventPlane[iFlow]);
0333       // If we are not allowing the event plane angles to be free parameters in the fit, fix them
0334       if (!doFreePlaneFit_) {
0335         flowFit_p_->FixParameter(4 + 2 * iFlow, eventPlane[iFlow]);
0336       }
0337     }
0338
0339     // Do the Fourier fit
0340     phi_h->Fit(flowFit_p_.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi());
0341
0342     // Put the fit parameters to the output vector
0343     for (int iParameter = 0; iParameter < nFlow * 2 + 1; iParameter++) {
0344       rhoFlowFitParamsOut->at(iParameter) = flowFit_p_->GetParameter(iParameter + 2);
0345     }
0346
0347     // Also add chi2 and ndf information to the output vector
0348     rhoFlowFitParamsOut->at(nFlow * 2 + 1) = flowFit_p_->GetChisquare();
0349     rhoFlowFitParamsOut->at(nFlow * 2 + 2) = flowFit_p_->GetNDF();
0350
0351     phi_h.reset();
0352     pfcuts.clear();
0353   }
0354
0355   // Define the produced output
0356   iEvent.put(std::move(rhoFlowFitParamsOut), "rhoFlowFitParams");
0357 }
0358
0359 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0360 void HiFJRhoFlowModulationProducer::beginStream(edm::StreamID) {}
0361
0362 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0363 void HiFJRhoFlowModulationProducer::endStream() {}
0364
0365 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0366 void HiFJRhoFlowModulationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0367   edm::ParameterSetDescription desc;
0368   desc.add<int>("minPfCandidatesPerEvent", 100);
0369   desc.add<bool>("doEvtPlane", false);
0370   desc.add<edm::InputTag>("EvtPlane", edm::InputTag("hiEvtPlane"));
0371   desc.add<bool>("doJettyExclusion", false);
0372   desc.add<double>("exclusionRadius", 0.4);
0373   desc.add<bool>("doFreePlaneFit", false);
0374   desc.add<edm::InputTag>("jetTag", edm::InputTag("ak4PFJetsForFlow"));
0375   desc.add<edm::InputTag>("pfCandSource", edm::InputTag("packedPFCandidates"));
0376   desc.add<int>("evtPlaneLevel", 0);
0377   desc.add<double>("pfCandidateEtaCut", 1.0);
0378   desc.add<double>("pfCandidateMinPtCut", 0.3);
0379   desc.add<double>("pfCandidateMaxPtCut", 3.0);
0380   desc.add<int>("firstFittedVn", 2);
0381   desc.add<int>("lastFittedVn", 3);
0382   descriptions.add("hiFJRhoFlowModulationProducer", desc);
0383 }
0384
0385 DEFINE_FWK_MODULE(HiFJRhoFlowModulationProducer);