CMS 3D CMS Logo

HiFJRhoFlowModulationProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer
4 // Class: HiFJRhoFlowModulationProducer
5 
25 
26 #include "TF1.h"
27 #include "TF2.h"
28 #include "TH1.h"
29 #include "TMath.h"
30 #include "Math/ProbFuncMathCore.h"
31 #include "TMinuitMinimizer.h"
32 
33 #include <algorithm>
34 #include <cmath>
35 #include <memory>
36 #include <string>
37 #include <utility>
38 #include <vector>
39 
40 namespace {
41 
42  double flowFunction(double* x, double* par) {
43  unsigned int nFlow = par[0]; // Number of fitted flow components is defined in the first parameter
44  unsigned int firstFittedVn = par[1]; // The first fitted flow component is defined in the second parameter
45 
46  // Add each component separately to the total fit value
47  double flowModulation = par[2];
48  for (unsigned int iFlow = 0; iFlow < nFlow; iFlow++) {
49  flowModulation +=
50  par[2] * 2.0 * par[2 * iFlow + 3] * std::cos((iFlow + firstFittedVn) * (x[0] - par[2 * iFlow + 4]));
51  }
52  return flowModulation;
53  }
54 
55  /*
56  * Weighing function for particle flow candidates in case jetty regions are exluded from the flow fit
57  * x[0] = DeltaPhi between jet axis and particle flow candidate in range [0,2Pi]
58  * x[1] = Absolute value of the jet eta
59  * par[0] = Eta cut applied for the particle flow candidates
60  * par[1] = Exclusion radius around the jet axis
61  *
62  * return: The weight to be used with this particle flow candidate is (1 + number provided by this function)
63  */
64  double weightFunction(double* x, double* par) {
65  // 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
66  if (x[0] > par[1])
67  return 0;
68 
69  // 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
70  double exclusionArea = std::sqrt(par[1] * par[1] - x[0] * x[0]);
71 
72  // Particle flow candidates are only considered until eta = par[0]. Check that the exclusion area does not go beyond that
73 
74  // First, check cases where the jet is found outside of the particle flow candidate acceptance in eta
75  if (x[1] > par[0]) {
76  // 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.
77  if (x[1] - exclusionArea > par[0])
78  return 0;
79 
80  // Now we know that part of the exclusion area will be inside of the acceptance. Check how big this is.
81  exclusionArea += par[0] - x[1];
82 
83  } else {
84  // 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
85  if (x[1] + exclusionArea > par[0]) {
86  exclusionArea += par[0] - x[1];
87  } else {
88  // If not, the the exclusion area is two times half of the nominal exclusion area
89  exclusionArea *= 2;
90  }
91  }
92 
93  // Normalize the exclusion area to the total acceptance. This number should be added to the total weight for the particle flow candidate
94  return (2 * par[0]) / (2 * par[0] - exclusionArea) - 1;
95  }
96 }; // namespace
97 
99 public:
101  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
102 
103 private:
104  void beginStream(edm::StreamID) override;
105  void produce(edm::Event&, const edm::EventSetup&) override;
106  void endStream() override;
107 
109  const bool doEvtPlane_;
110  const bool doFreePlaneFit_;
111  const bool doJettyExclusion_;
112  const double exclusionRadius_;
113  const int evtPlaneLevel_;
114  const double pfCandidateEtaCut_;
115  const double pfCandidateMinPtCut_;
116  const double pfCandidateMaxPtCut_;
122  std::unique_ptr<TF1> flowFit_p_;
123  std::unique_ptr<TF2> pfWeightFunction_;
125 };
127  : minPfCandidatesPerEvent_(iConfig.getParameter<int>("minPfCandidatesPerEvent")),
128  doEvtPlane_(iConfig.getParameter<bool>("doEvtPlane")),
129  doFreePlaneFit_(iConfig.getParameter<bool>("doFreePlaneFit")),
130  doJettyExclusion_(iConfig.getParameter<bool>("doJettyExclusion")),
131  exclusionRadius_(iConfig.getParameter<double>("exclusionRadius")),
132  evtPlaneLevel_(iConfig.getParameter<int>("evtPlaneLevel")),
133  pfCandidateEtaCut_(iConfig.getParameter<double>("pfCandidateEtaCut")),
134  pfCandidateMinPtCut_(iConfig.getParameter<double>("pfCandidateMinPtCut")),
135  pfCandidateMaxPtCut_(iConfig.getParameter<double>("pfCandidateMaxPtCut")),
136  firstFittedVn_(iConfig.getParameter<int>("firstFittedVn")),
137  lastFittedVn_(iConfig.getParameter<int>("lastFittedVn")),
138  jetToken_(doJettyExclusion_ ? consumes<reco::JetView>(iConfig.getParameter<edm::InputTag>("jetTag"))
139  : edm::EDGetTokenT<reco::JetView>()),
140  pfCandidateToken_(consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"))),
141  evtPlaneToken_(consumes<reco::EvtPlaneCollection>(iConfig.getParameter<edm::InputTag>("EvtPlane"))) {
142  produces<std::vector<double>>("rhoFlowFitParams");
143 
144  // Converter to get PF candidate ID from packed candidates
146 
147  // Sanity check for the fitted flow component input
148  if (firstFittedVn_ < 1)
149  firstFittedVn_ = 2;
150  if (lastFittedVn_ < 1)
151  lastFittedVn_ = 2;
153  int flipper = lastFittedVn_;
155  firstFittedVn_ = flipper;
156  }
157 
158  // Calculate the number of flow components
159  const int nFlow = lastFittedVn_ - firstFittedVn_ + 1;
160 
161  // Define all fit and weight functions
162  TMinuitMinimizer::UseStaticMinuit(false);
163  flowFit_p_ = std::make_unique<TF1>("flowFit", flowFunction, -TMath::Pi(), TMath::Pi(), nFlow * 2 + 3);
164  flowFit_p_->FixParameter(0, nFlow); // The first parameter defines the number of fitted flow components
165  flowFit_p_->FixParameter(1, firstFittedVn_); // The second parameter defines the first fitted flow component
166  pfWeightFunction_ = std::make_unique<TF2>("weightFunction", weightFunction, 0, TMath::Pi(), -5, 5, 2);
167  pfWeightFunction_->SetParameter(0, pfCandidateEtaCut_); // Set the allowed eta range for particle flow candidates
168  pfWeightFunction_->SetParameter(1, exclusionRadius_); // Set the exclusion radius around the jets
169 }
170 
171 // ------------ method called to produce the data ------------
173  // Get the particle flow candidate collection
174  auto const& pfCands = iEvent.get(pfCandidateToken_);
175 
176  // If we are reading the event plane information for the forest, find the event planes
177  std::array<float, hi::NumEPNames> hiEvtPlane;
178  if (doEvtPlane_) {
179  auto const& evtPlanes = iEvent.get(evtPlaneToken_);
180  assert(evtPlanes.size() == hi::NumEPNames);
181  std::transform(evtPlanes.begin(), evtPlanes.end(), hiEvtPlane.begin(), [this](auto const& ePlane) -> float {
182  return ePlane.angle(evtPlaneLevel_);
183  });
184  }
185 
186  // If jetty regions are excluded from the flow fits, read the jet collection
188  if (doJettyExclusion_) {
189  iEvent.getByToken(jetToken_, jets);
190  }
191 
192  int nFill = 0;
193 
194  // Initialize arrays for event planes
195  // nFlow: Number of flow components in the fit to particle flow condidate distribution
196  const int nFlow = lastFittedVn_ - firstFittedVn_ + 1;
197  double eventPlane[nFlow];
198  for (int iFlow = 0; iFlow < nFlow; iFlow++)
199  eventPlane[iFlow] = -100;
200 
201  // Initialize the output vector with flow fit components
202  // 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
203  const int nParamVals = nFlow * 2 + 4;
204  auto rhoFlowFitParamsOut = std::make_unique<std::vector<double>>(nParamVals, 1e-6);
205 
206  // Set the parameters related to flow fit to zero
207  for (int iParameter = 0; iParameter < nFlow * 2 + 1; iParameter++) {
208  rhoFlowFitParamsOut->at(iParameter) = 0;
209  }
210  // Set the first fitted flow component to the last index of the output vector
211  rhoFlowFitParamsOut->at(nFlow * 2 + 3) = firstFittedVn_;
212 
213  // Initialize arrays for manual event plane calculation
214  double eventPlaneCos[nFlow];
215  double eventPlaneSin[nFlow];
216  for (int iFlow = 0; iFlow < nFlow; iFlow++) {
217  eventPlaneCos[iFlow] = 0;
218  eventPlaneSin[iFlow] = 0;
219  }
220 
221  // Define helper variables for looping over particle flow candidates
222  std::vector<bool> pfcuts(pfCands.size(), false);
223  int iCand = -1;
224  double thisPfCandidateWeight = 0;
225  double deltaPhiJetPf = 0;
226  std::vector<double> pfCandidateWeight(pfCands.size(), 1);
227 
228  // Loop over the particle flow candidates
229  for (auto const& pfCandidate : pfCands) {
230  iCand++; // Update the particle flow candidate index
231 
232  // Find the PF candidate ID from the packed candidates collection
233  auto particleFlowCandidateID = converter_.translatePdgIdToType(pfCandidate.pdgId());
234 
235  // This cut selects particle flow candidates that are charged hadrons. The ID numbers are:
236  // 0 = undefined
237  // 1 = charged hadron
238  // 2 = electron
239  // 3 = muon
240  // 4 = photon
241  // 5 = neutral hadron
242  // 6 = HF tower identified as a hadron
243  // 7 = HF tower identified as an EM particle
244  if (particleFlowCandidateID != 1)
245  continue;
246 
247  // Kinematic cuts for particle flow candidate pT and eta
248  if (std::abs(pfCandidate.eta()) > pfCandidateEtaCut_)
249  continue;
250  if (pfCandidate.pt() < pfCandidateMinPtCut_)
251  continue;
252  if (pfCandidate.pt() > pfCandidateMaxPtCut_)
253  continue;
254 
255  nFill++; // Use same nFill criterion with and without jetty region subtraction
256 
257  // If the jetty regions are excluded from the fit, find which particle flow candidates are close to jets
258  thisPfCandidateWeight = 1;
259  if (doJettyExclusion_) {
260  bool isGood = true;
261  for (auto const& jet : *jets) {
262  if (deltaR2(jet, pfCandidate) < exclusionRadius_ * exclusionRadius_) {
263  isGood = false;
264  break;
265  } else {
266  // 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.
267  deltaPhiJetPf = std::abs(pfCandidate.phi() - jet.phi());
268  if (deltaPhiJetPf > TMath::Pi())
269  deltaPhiJetPf = TMath::Pi() * 2 - deltaPhiJetPf;
270  // Weight currently does not take into account overlapping jets
271  thisPfCandidateWeight += pfWeightFunction_->Eval(deltaPhiJetPf, std::abs(jet.eta()));
272  }
273  }
274  // Do not use this particle flow candidate in the manual event plane calculation
275  if (!isGood) {
276  continue;
277  }
278  }
279 
280  // Update the weight for this particle flow candidate
281  pfCandidateWeight[iCand] = thisPfCandidateWeight;
282 
283  // This particle flow candidate passes all the cuts
284  pfcuts[iCand] = true;
285 
286  // If the event plane is calculated manually, add this particle flow candidate to the calculation
287  if (!doEvtPlane_) {
288  for (int iFlow = 0; iFlow < nFlow; iFlow++) {
289  eventPlaneCos[iFlow] += std::cos((iFlow + firstFittedVn_) * pfCandidate.phi()) * thisPfCandidateWeight;
290  eventPlaneSin[iFlow] += std::sin((iFlow + firstFittedVn_) * pfCandidate.phi()) * thisPfCandidateWeight;
291  }
292  }
293  }
294 
295  // Determine the event plane angle
296  if (!doEvtPlane_) {
297  // Manual calculation for the event plane angle using the particle flow candidates
298  for (int iFlow = 0; iFlow < nFlow; iFlow++) {
299  eventPlane[iFlow] = std::atan2(eventPlaneSin[iFlow], eventPlaneCos[iFlow]) / (iFlow + firstFittedVn_);
300  }
301  } else {
302  // Read the event plane angle determined from the HF calorimeters from the HiForest
303  // Only v2 and v3 are available in the HiForest. This option should not be used if other flow components are fitted
304  int halfWay = nFlow / 2;
305  for (int iFlow = 0; iFlow < halfWay; iFlow++) {
306  eventPlane[iFlow] = hiEvtPlane[hi::HF2];
307  }
308  for (int iFlow = halfWay; iFlow < nFlow; iFlow++) {
309  eventPlane[iFlow] = hiEvtPlane[hi::HF3];
310  }
311  }
312 
313  // Do the flow fit provided that there are enough particle flow candidates in the event and that the event planes have been determined successfully
314  int pfcuts_count = 0;
315  if (nFill >= minPfCandidatesPerEvent_ && eventPlane[0] > -99) {
316  // Create a particle flow candidate phi-histogram
317  const int nPhiBins = std::max(10, nFill / 30);
318  std::string name = "phiTestIEta4_" + std::to_string(iEvent.id().event()) + "_h";
319  std::unique_ptr<TH1F> phi_h = std::make_unique<TH1F>(name.data(), "", nPhiBins, -TMath::Pi(), TMath::Pi());
320  phi_h->SetDirectory(nullptr);
321  for (auto const& pfCandidate : pfCands) {
322  if (pfcuts.at(pfcuts_count)) { // Only use particle flow candidates that pass all cuts
323  phi_h->Fill(pfCandidate.phi(), pfCandidateWeight.at(pfcuts_count));
324  }
325  pfcuts_count++;
326  }
327 
328  // Set initial values for the fit parameters
329  flowFit_p_->SetParameter(2, 10);
330  for (int iFlow = 0; iFlow < nFlow; iFlow++) {
331  flowFit_p_->SetParameter(3 + 2 * iFlow, 0);
332  flowFit_p_->SetParameter(4 + 2 * iFlow, eventPlane[iFlow]);
333  // If we are not allowing the event plane angles to be free parameters in the fit, fix them
334  if (!doFreePlaneFit_) {
335  flowFit_p_->FixParameter(4 + 2 * iFlow, eventPlane[iFlow]);
336  }
337  }
338 
339  // Do the Fourier fit
340  phi_h->Fit(flowFit_p_.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi());
341 
342  // Put the fit parameters to the output vector
343  for (int iParameter = 0; iParameter < nFlow * 2 + 1; iParameter++) {
344  rhoFlowFitParamsOut->at(iParameter) = flowFit_p_->GetParameter(iParameter + 2);
345  }
346 
347  // Also add chi2 and ndf information to the output vector
348  rhoFlowFitParamsOut->at(nFlow * 2 + 1) = flowFit_p_->GetChisquare();
349  rhoFlowFitParamsOut->at(nFlow * 2 + 2) = flowFit_p_->GetNDF();
350 
351  phi_h.reset();
352  pfcuts.clear();
353  }
354 
355  // Define the produced output
356  iEvent.put(std::move(rhoFlowFitParamsOut), "rhoFlowFitParams");
357 }
358 
359 // ------------ method called once each stream before processing any runs, lumis or events ------------
361 
362 // ------------ method called once each stream after processing all runs, lumis and events ------------
364 
365 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
368  desc.add<int>("minPfCandidatesPerEvent", 100);
369  desc.add<bool>("doEvtPlane", false);
370  desc.add<edm::InputTag>("EvtPlane", edm::InputTag("hiEvtPlane"));
371  desc.add<bool>("doJettyExclusion", false);
372  desc.add<double>("exclusionRadius", 0.4);
373  desc.add<bool>("doFreePlaneFit", false);
374  desc.add<edm::InputTag>("jetTag", edm::InputTag("ak4PFJetsForFlow"));
375  desc.add<edm::InputTag>("pfCandSource", edm::InputTag("packedPFCandidates"));
376  desc.add<int>("evtPlaneLevel", 0);
377  desc.add<double>("pfCandidateEtaCut", 1.0);
378  desc.add<double>("pfCandidateMinPtCut", 0.3);
379  desc.add<double>("pfCandidateMaxPtCut", 3.0);
380  desc.add<int>("firstFittedVn", 2);
381  desc.add<int>("lastFittedVn", 3);
382  descriptions.add("hiFJRhoFlowModulationProducer", desc);
383 }
384 
const double Pi
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const edm::EDGetTokenT< reco::EvtPlaneCollection > evtPlaneToken_
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< pat::PackedCandidate > PackedCandidateCollection
void beginStream(edm::StreamID) override
edm::View< Jet > JetView
edm references
Definition: JetCollection.h:11
assert(be >=bs)
static std::string to_string(const XMLCh *ch)
Definition: HeavyIon.h:7
void produce(edm::Event &, const edm::EventSetup &) override
int iEvent
Definition: GenABIO.cc:224
T sqrt(T t)
Definition: SSEVec.h:23
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::vector< EvtPlane > EvtPlaneCollection
Definition: EvtPlane.h:62
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
const edm::EDGetTokenT< reco::JetView > jetToken_
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
fixed size matrix
ParticleType translatePdgIdToType(int pdgid) const
Definition: PFCandidate.cc:231
HLT enums.
static constexpr int nPhiBins
float x
static const int NumEPNames
HiFJRhoFlowModulationProducer(const edm::ParameterSet &)
def move(src, dest)
Definition: eostools.py:511
MPlex< T, D1, D2, N > atan2(const MPlex< T, D1, D2, N > &y, const MPlex< T, D1, D2, N > &x)
Definition: Matriplex.h:648
const edm::EDGetTokenT< pat::PackedCandidateCollection > pfCandidateToken_
unsigned transform(const HcalDetId &id, unsigned transformCode)