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 "TH1.h"
28 #include "TMath.h"
29 #include "TMinuitMinimizer.h"
30 
31 #include <algorithm>
32 #include <cmath>
33 #include <memory>
34 #include <string>
35 #include <utility>
36 #include <vector>
37 
38 namespace {
39  double lineFunction(double* x, double* par) { return par[0]; }
40  double flowFunction(double* x, double* par) {
41  return par[0] * (1. + 2. * (par[1] * std::cos(2. * (x[0] - par[2])) + par[3] * std::cos(3. * (x[0] - par[4]))));
42  }
43 }; // namespace
44 
46 public:
48  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
49 
50 private:
51  void produce(edm::Event&, const edm::EventSetup&) override;
52 
53  const bool doEvtPlane_;
54  const bool doFreePlaneFit_;
55  const bool doJettyExclusion_;
56  const int evtPlaneLevel_;
60  std::unique_ptr<TF1> lineFit_p_;
61  std::unique_ptr<TF1> flowFit_p_;
62 };
64  : doEvtPlane_(iConfig.getParameter<bool>("doEvtPlane")),
65  doFreePlaneFit_(iConfig.getParameter<bool>("doFreePlaneFit")),
66  doJettyExclusion_(iConfig.getParameter<bool>("doJettyExclusion")),
67  evtPlaneLevel_(iConfig.getParameter<int>("evtPlaneLevel")),
68  jetToken_(doJettyExclusion_ ? consumes<reco::JetView>(iConfig.getParameter<edm::InputTag>("jetTag"))
69  : edm::EDGetTokenT<reco::JetView>()),
70  pfCandsToken_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"))),
71  evtPlaneToken_(consumes<reco::EvtPlaneCollection>(iConfig.getParameter<edm::InputTag>("EvtPlane"))) {
72  produces<std::vector<double>>("rhoFlowFitParams");
73  TMinuitMinimizer::UseStaticMinuit(false);
74  lineFit_p_ = std::make_unique<TF1>("lineFit", lineFunction, -TMath::Pi(), TMath::Pi());
75  flowFit_p_ = std::make_unique<TF1>("flowFit", flowFunction, -TMath::Pi(), TMath::Pi());
76 }
77 
78 // ------------ method called to produce the data ------------
80  int nPhiBins = 10;
81  auto const& pfCands = iEvent.get(pfCandsToken_);
82  static constexpr int kMaxEvtPlane = 29;
83  std::array<float, kMaxEvtPlane> hiEvtPlane;
84 
85  if (doEvtPlane_) {
86  auto const& evtPlanes = iEvent.get(evtPlaneToken_);
87  assert(evtPlanes.size() == hi::NumEPNames);
88  std::transform(evtPlanes.begin(), evtPlanes.end(), hiEvtPlane.begin(), [this](auto const& ePlane) -> float {
89  return ePlane.angle(evtPlaneLevel_);
90  });
91  }
92 
95  iEvent.getByToken(jetToken_, jets);
96 
97  int nFill = 0;
98 
99  double eventPlane2 = -100;
100  double eventPlane3 = -100;
101 
102  constexpr int nParamVals = 9;
103  auto rhoFlowFitParamsOut = std::make_unique<std::vector<double>>(nParamVals, 1e-6);
104 
105  rhoFlowFitParamsOut->at(0) = 0;
106  rhoFlowFitParamsOut->at(1) = 0;
107  rhoFlowFitParamsOut->at(2) = 0;
108  rhoFlowFitParamsOut->at(3) = 0;
109  rhoFlowFitParamsOut->at(4) = 0;
110 
111  double eventPlane2Cos = 0;
112  double eventPlane2Sin = 0;
113 
114  double eventPlane3Cos = 0;
115  double eventPlane3Sin = 0;
116 
117  std::vector<bool> pfcuts(pfCands.size(), false);
118  int iCand = -1;
119  for (auto const& pfCandidate : pfCands) {
120  iCand++;
121  if (pfCandidate.particleId() != 1 || std::abs(pfCandidate.eta()) > 1.0 || pfCandidate.pt() < .3 ||
122  pfCandidate.pt() > 3.) {
123  continue;
124  }
125 
126  if (doJettyExclusion_) {
127  bool isGood = true;
128  for (auto const& jet : *jets) {
129  if (deltaR2(jet, pfCandidate) < .16) {
130  isGood = false;
131  break;
132  }
133  }
134  if (!isGood) {
135  continue;
136  }
137  }
138 
139  nFill++;
140  pfcuts[iCand] = true;
141 
142  if (!doEvtPlane_) {
143  eventPlane2Cos += std::cos(2 * pfCandidate.phi());
144  eventPlane2Sin += std::sin(2 * pfCandidate.phi());
145 
146  eventPlane3Cos += std::cos(3 * pfCandidate.phi());
147  eventPlane3Sin += std::sin(3 * pfCandidate.phi());
148  }
149  }
150 
151  if (!doEvtPlane_) {
152  eventPlane2 = std::atan2(eventPlane2Sin, eventPlane2Cos) / 2.;
153  eventPlane3 = std::atan2(eventPlane3Sin, eventPlane3Cos) / 3.;
154  } else {
155  eventPlane2 = hiEvtPlane[hi::HF2];
156  eventPlane3 = hiEvtPlane[hi::HF3];
157  }
158  int pfcuts_count = 0;
159  if (nFill >= 100 && eventPlane2 > -99) {
160  nPhiBins = std::max(10, nFill / 30);
161 
162  std::string name = "phiTestIEta4_" + std::to_string(iEvent.id().event()) + "_h";
163  std::string nameFlat = "phiTestIEta4_Flat_" + std::to_string(iEvent.id().event()) + "_h";
164  std::unique_ptr<TH1F> phi_h = std::make_unique<TH1F>(name.data(), "", nPhiBins, -TMath::Pi(), TMath::Pi());
165  phi_h->SetDirectory(nullptr);
166  for (auto const& pfCandidate : pfCands) {
167  if (pfcuts.at(pfcuts_count))
168  phi_h->Fill(pfCandidate.phi());
169  pfcuts_count++;
170  }
171  flowFit_p_->SetParameter(0, 10);
172  flowFit_p_->SetParameter(1, 0);
173  flowFit_p_->SetParameter(2, eventPlane2);
174  flowFit_p_->SetParameter(3, 0);
175  flowFit_p_->SetParameter(4, eventPlane3);
176  if (!doFreePlaneFit_) {
177  flowFit_p_->FixParameter(2, eventPlane2);
178  flowFit_p_->FixParameter(4, eventPlane3);
179  }
180 
181  lineFit_p_->SetParameter(0, 10);
182 
183  phi_h->Fit(flowFit_p_.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi());
184  phi_h->Fit(lineFit_p_.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi());
185  rhoFlowFitParamsOut->at(0) = flowFit_p_->GetParameter(0);
186  rhoFlowFitParamsOut->at(1) = flowFit_p_->GetParameter(1);
187  rhoFlowFitParamsOut->at(2) = flowFit_p_->GetParameter(2);
188  rhoFlowFitParamsOut->at(3) = flowFit_p_->GetParameter(3);
189  rhoFlowFitParamsOut->at(4) = flowFit_p_->GetParameter(4);
190 
191  rhoFlowFitParamsOut->at(5) = flowFit_p_->GetChisquare();
192  rhoFlowFitParamsOut->at(6) = flowFit_p_->GetNDF();
193 
194  rhoFlowFitParamsOut->at(7) = lineFit_p_->GetChisquare();
195  rhoFlowFitParamsOut->at(8) = lineFit_p_->GetNDF();
196 
197  phi_h.reset();
198  pfcuts.clear();
199  }
200 
201  iEvent.put(std::move(rhoFlowFitParamsOut), "rhoFlowFitParams");
202 }
203 
204 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
207  desc.add<bool>("doEvtPlane", false);
208  desc.add<edm::InputTag>("EvtPlane", edm::InputTag("hiEvtPlane"));
209  desc.add<bool>("doJettyExclusion", false);
210  desc.add<bool>("doFreePlaneFit", false);
211  desc.add<bool>("doFlatTest", false);
212  desc.add<edm::InputTag>("jetTag", edm::InputTag("ak4PFJets"));
213  desc.add<edm::InputTag>("pfCandSource", edm::InputTag("particleFlow"));
214  desc.add<int>("evtPlaneLevel", 0);
215  descriptions.add("hiFJRhoFlowModulationProducer", desc);
216 }
217 
ConfigurationDescriptions.h
EvtPlane.h
Handle.h
electrons_cff.bool
bool
Definition: electrons_cff.py:366
HiFJRhoFlowModulationProducer::HiFJRhoFlowModulationProducer
HiFJRhoFlowModulationProducer(const edm::ParameterSet &)
Definition: HiFJRhoFlowModulationProducer.cc:63
PFCandidate.h
edm::EDGetTokenT
Definition: EDGetToken.h:33
edm
HLT enums.
Definition: AlignableModifier.h:19
HLT_FULL_cff.InputTag
InputTag
Definition: HLT_FULL_cff.py:89281
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
HiFJRhoFlowModulationProducer
Definition: HiFJRhoFlowModulationProducer.cc:45
cms::cuda::assert
assert(be >=bs)
HiFJRhoFlowModulationProducer::doEvtPlane_
const bool doEvtPlane_
Definition: HiFJRhoFlowModulationProducer.cc:53
EDProducer.h
Jet.h
HiFJRhoFlowModulationProducer::produce
void produce(edm::Event &, const edm::EventSetup &) override
Definition: HiFJRhoFlowModulationProducer.cc:79
reco::JetView
edm::View< Jet > JetView
edm references
Definition: JetCollection.h:11
singleTopDQM_cfi.jets
jets
Definition: singleTopDQM_cfi.py:42
HiFJRhoFlowModulationProducer::doJettyExclusion_
const bool doJettyExclusion_
Definition: HiFJRhoFlowModulationProducer.cc:55
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
HiFJRhoFlowModulationProducer::jetToken_
const edm::EDGetTokenT< reco::JetView > jetToken_
Definition: HiFJRhoFlowModulationProducer.cc:57
edm::Handle
Definition: AssociativeIterator.h:50
HiFJRhoFlowModulationProducer::pfCandsToken_
const edm::EDGetTokenT< reco::PFCandidateCollection > pfCandsToken_
Definition: HiFJRhoFlowModulationProducer.cc:58
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MakerMacros.h
HiFJRhoFlowModulationProducer::doFreePlaneFit_
const bool doFreePlaneFit_
Definition: HiFJRhoFlowModulationProducer.cc:54
HiFJRhoFlowModulationProducer::lineFit_p_
std::unique_ptr< TF1 > lineFit_p_
Definition: HiFJRhoFlowModulationProducer.cc:60
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
EventID.h
hi::HF2
Definition: HiEvtPlaneList.h:23
hi::HF3
Definition: HiEvtPlaneList.h:24
HcalDetIdTransform::transform
unsigned transform(const HcalDetId &id, unsigned transformCode)
Definition: HcalDetIdTransform.cc:7
JetCollection.h
ParameterSetDescription.h
EDGetToken.h
reco::EvtPlaneCollection
std::vector< EvtPlane > EvtPlaneCollection
Definition: EvtPlane.h:62
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::ParameterSet
Definition: ParameterSet.h:47
HiEvtPlane_cfi.hiEvtPlane
hiEvtPlane
Definition: HiEvtPlane_cfi.py:3
Event.h
deltaR.h
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
HiFJRhoFlowModulationProducer::flowFit_p_
std::unique_ptr< TF1 > flowFit_p_
Definition: HiFJRhoFlowModulationProducer.cc:61
createfilelist.int
int
Definition: createfilelist.py:10
iEvent
int iEvent
Definition: GenABIO.cc:224
edm::stream::EDProducer
Definition: EDProducer.h:38
edm::EventSetup
Definition: EventSetup.h:58
HiEvtPlaneList.h
l1t::PFCandidateCollection
std::vector< l1t::PFCandidate > PFCandidateCollection
Definition: PFCandidate.h:57
HiFJRhoFlowModulationProducer::evtPlaneLevel_
const int evtPlaneLevel_
Definition: HiFJRhoFlowModulationProducer.cc:56
submitPVResolutionJobs.desc
string desc
Definition: submitPVResolutionJobs.py:251
eostools.move
def move(src, dest)
Definition: eostools.py:511
HLTMuonOfflineAnalyzer_cfi.deltaR2
deltaR2
Definition: HLTMuonOfflineAnalyzer_cfi.py:105
metsig::jet
Definition: SignAlgoResolutions.h:47
genVertex_cff.x
x
Definition: genVertex_cff.py:12
HiFJRhoFlowModulationProducer::evtPlaneToken_
const edm::EDGetTokenT< reco::EvtPlaneCollection > evtPlaneToken_
Definition: HiFJRhoFlowModulationProducer.cc:59
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
hi::NumEPNames
static const int NumEPNames
Definition: HiEvtPlaneList.h:28
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
View.h
ParameterSet.h
ParameterSetDescriptionFiller.h
ecaldqm::binning::nPhiBins
Definition: MESetBinningUtils.h:70
edm::Event
Definition: Event.h:73
StreamID.h
edm::InputTag
Definition: InputTag.h:15
HiFJRhoFlowModulationProducer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: HiFJRhoFlowModulationProducer.cc:205
PFCandidateFwd.h
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37