30 #include "Math/ProbFuncMathCore.h" 31 #include "TMinuitMinimizer.h" 42 double flowFunction(
double*
x,
double* par) {
43 unsigned int nFlow = par[0];
44 unsigned int firstFittedVn = par[1];
47 double flowModulation = par[2];
48 for (
unsigned int iFlow = 0; iFlow < nFlow; iFlow++) {
50 par[2] * 2.0 * par[2 * iFlow + 3] *
std::cos((iFlow + firstFittedVn) * (
x[0] - par[2 * iFlow + 4]));
52 return flowModulation;
64 double weightFunction(
double*
x,
double* par) {
70 double exclusionArea =
std::sqrt(par[1] * par[1] -
x[0] *
x[0]);
77 if (
x[1] - exclusionArea > par[0])
81 exclusionArea += par[0] -
x[1];
85 if (
x[1] + exclusionArea > par[0]) {
86 exclusionArea += par[0] -
x[1];
94 return (2 * par[0]) / (2 * par[0] - exclusionArea) - 1;
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")),
142 produces<std::vector<double>>(
"rhoFlowFitParams");
162 TMinuitMinimizer::UseStaticMinuit(
false);
197 double eventPlane[nFlow];
198 for (
int iFlow = 0; iFlow < nFlow; iFlow++)
199 eventPlane[iFlow] = -100;
203 const int nParamVals = nFlow * 2 + 4;
204 auto rhoFlowFitParamsOut = std::make_unique<std::vector<double>>(nParamVals, 1
e-6);
207 for (
int iParameter = 0; iParameter < nFlow * 2 + 1; iParameter++) {
208 rhoFlowFitParamsOut->at(iParameter) = 0;
214 double eventPlaneCos[nFlow];
215 double eventPlaneSin[nFlow];
216 for (
int iFlow = 0; iFlow < nFlow; iFlow++) {
217 eventPlaneCos[iFlow] = 0;
218 eventPlaneSin[iFlow] = 0;
222 std::vector<bool> pfcuts(pfCands.size(),
false);
224 double thisPfCandidateWeight = 0;
225 double deltaPhiJetPf = 0;
226 std::vector<double> pfCandidateWeight(pfCands.size(), 1);
229 for (
auto const& pfCandidate : pfCands) {
244 if (particleFlowCandidateID != 1)
258 thisPfCandidateWeight = 1;
261 for (
auto const&
jet : *
jets) {
267 deltaPhiJetPf =
std::abs(pfCandidate.phi() -
jet.phi());
269 deltaPhiJetPf =
TMath::Pi() * 2 - deltaPhiJetPf;
281 pfCandidateWeight[iCand] = thisPfCandidateWeight;
284 pfcuts[iCand] =
true;
288 for (
int iFlow = 0; iFlow < nFlow; iFlow++) {
298 for (
int iFlow = 0; iFlow < nFlow; iFlow++) {
304 int halfWay = nFlow / 2;
305 for (
int iFlow = 0; iFlow < halfWay; iFlow++) {
308 for (
int iFlow = halfWay; iFlow < nFlow; iFlow++) {
314 int pfcuts_count = 0;
320 phi_h->SetDirectory(
nullptr);
321 for (
auto const& pfCandidate : pfCands) {
322 if (pfcuts.at(pfcuts_count)) {
323 phi_h->Fill(pfCandidate.phi(), pfCandidateWeight.at(pfcuts_count));
330 for (
int iFlow = 0; iFlow < nFlow; iFlow++) {
332 flowFit_p_->SetParameter(4 + 2 * iFlow, eventPlane[iFlow]);
335 flowFit_p_->FixParameter(4 + 2 * iFlow, eventPlane[iFlow]);
343 for (
int iParameter = 0; iParameter < nFlow * 2 + 1; iParameter++) {
344 rhoFlowFitParamsOut->at(iParameter) =
flowFit_p_->GetParameter(iParameter + 2);
348 rhoFlowFitParamsOut->at(nFlow * 2 + 1) =
flowFit_p_->GetChisquare();
349 rhoFlowFitParamsOut->at(nFlow * 2 + 2) =
flowFit_p_->GetNDF();
368 desc.add<
int>(
"minPfCandidatesPerEvent", 100);
369 desc.add<
bool>(
"doEvtPlane",
false);
371 desc.add<
bool>(
"doJettyExclusion",
false);
372 desc.add<
double>(
"exclusionRadius", 0.4);
373 desc.add<
bool>(
"doFreePlaneFit",
false);
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);
reco::PFCandidate converter_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const edm::EDGetTokenT< reco::EvtPlaneCollection > evtPlaneToken_
Sin< T >::type sin(const T &t)
std::vector< pat::PackedCandidate > PackedCandidateCollection
void beginStream(edm::StreamID) override
void endStream() override
edm::View< Jet > JetView
edm references
std::unique_ptr< TF2 > pfWeightFunction_
static std::string to_string(const XMLCh *ch)
void produce(edm::Event &, const edm::EventSetup &) override
Cos< T >::type cos(const T &t)
std::vector< EvtPlane > EvtPlaneCollection
Abs< T >::type abs(const T &t)
#define DEFINE_FWK_MODULE(type)
const bool doFreePlaneFit_
const edm::EDGetTokenT< reco::JetView > jetToken_
void add(std::string const &label, ParameterSetDescription const &psetDescription)
const double exclusionRadius_
const bool doJettyExclusion_
Particle reconstructed by the particle flow algorithm.
ParticleType translatePdgIdToType(int pdgid) const
const double pfCandidateEtaCut_
static constexpr int nPhiBins
const double pfCandidateMinPtCut_
static const int NumEPNames
const double pfCandidateMaxPtCut_
HiFJRhoFlowModulationProducer(const edm::ParameterSet &)
const int minPfCandidatesPerEvent_
std::unique_ptr< TF1 > flowFit_p_
MPlex< T, D1, D2, N > atan2(const MPlex< T, D1, D2, N > &y, const MPlex< T, D1, D2, N > &x)
const edm::EDGetTokenT< pat::PackedCandidateCollection > pfCandidateToken_