Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Flow modulated CS subtraction for heavy ions #30057

Merged
merged 24 commits into from
Jul 22, 2020
Merged

Flow modulated CS subtraction for heavy ions #30057

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
ab78b11
flow cs subtraction
stepobr May 22, 2020
72a0a56
flow cs subtraction pt2
stepobr May 22, 2020
793ec3c
fix code and style checks
stepobr May 22, 2020
273940e
drop HiPFCandCleaner for a different PR
stepobr May 23, 2020
f18011b
add new cs jets as default for reco sequence
stepobr May 23, 2020
84041f2
add back
stepobr May 30, 2020
8c0db1b
add missing task
stepobr May 30, 2020
d4b72a6
Merged csflow from repository stepobr with cms-merge-topic
stepobr Jun 15, 2020
a0e1538
remove headers, move to thread-aware producers
stepobr Jun 15, 2020
1839096
remove cfi
stepobr Jun 29, 2020
e6189c3
remove HiPFCandCleaner
stepobr Jul 2, 2020
e865781
apply suggestions
stepobr Jul 3, 2020
f2c6d16
declare const in a common place
stepobr Jul 10, 2020
daf657b
apply suggestions
stepobr Jul 10, 2020
62a1138
apply code check
stepobr Jul 10, 2020
c226d35
change etaedge to constexpr
stepobr Jul 16, 2020
1f6f809
remove unnecessary type conversion
stepobr Jul 16, 2020
027d512
style check fix
stepobr Jul 16, 2020
98d7ea0
move constant to HITowerHelper
stepobr Jul 16, 2020
5c8bc8e
Merge branch 'master' into csflow
stepobr Jul 16, 2020
922e835
Update RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc
stepobr Jul 17, 2020
c2b6c81
fix style check
stepobr Jul 17, 2020
05f2e98
set useHF = True
stepobr Jul 17, 2020
968daf9
minor optimization of .eta() calls
stepobr Jul 21, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
251 changes: 251 additions & 0 deletions RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,251 @@
// -*- C++ -*-
//
// Package: RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer
// Class: HiFJRhoFlowModulationProducer

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/Provenance/interface/EventID.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/src/WorkerMaker.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
#include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
#include "DataFormats/JetReco/interface/JetCollection.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Utilities/interface/StreamID.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

#include "TF1.h"
#include "TH1.h"
#include "TMath.h"

#include <algorithm>
#include <cmath>
#include <memory>
#include <string>
#include <utility>
#include <vector>

class HiFJRhoFlowModulationProducer : public edm::stream::EDProducer<> {
public:
explicit HiFJRhoFlowModulationProducer(const edm::ParameterSet&);

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
void beginStream(edm::StreamID) override;
void produce(edm::Event&, const edm::EventSetup&) override;
void endStream() override;

bool doEvtPlane_;
bool doFreePlaneFit_;
bool doJettyExclusion_;
int evtPlaneLevel_;
edm::EDGetTokenT<reco::JetView> jetTag_;
stepobr marked this conversation as resolved.
Show resolved Hide resolved
edm::EDGetTokenT<reco::PFCandidateCollection> pfCandsToken_;
edm::EDGetTokenT<reco::EvtPlaneCollection> evtPlaneToken_;
stepobr marked this conversation as resolved.
Show resolved Hide resolved

float* hiEvtPlane;
stepobr marked this conversation as resolved.
Show resolved Hide resolved
};

HiFJRhoFlowModulationProducer::HiFJRhoFlowModulationProducer(const edm::ParameterSet& iConfig)
: doEvtPlane_(iConfig.getParameter<bool>("doEvtPlane")),
doFreePlaneFit_(iConfig.getParameter<bool>("doFreePlaneFit")),
doJettyExclusion_(iConfig.getParameter<bool>("doJettyExclusion")),
evtPlaneLevel_(iConfig.getParameter<int>("evtPlaneLevel")),
jetTag_(consumes<reco::JetView>(iConfig.getParameter<edm::InputTag>("jetTag"))),
stepobr marked this conversation as resolved.
Show resolved Hide resolved
pfCandsToken_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"))),
evtPlaneToken_(consumes<reco::EvtPlaneCollection>(iConfig.getParameter<edm::InputTag>("EvtPlane"))) {
produces<std::vector<double>>("rhoFlowFitParams");
stepobr marked this conversation as resolved.
Show resolved Hide resolved
}
stepobr marked this conversation as resolved.
Show resolved Hide resolved

// ------------ method called to produce the data ------------
void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
edm::Handle<reco::PFCandidateCollection> pfCands;
iEvent.getByToken(pfCandsToken_, pfCands);
stepobr marked this conversation as resolved.
Show resolved Hide resolved

edm::Handle<reco::EvtPlaneCollection> evtPlanes;
if (doEvtPlane_) {
iEvent.getByToken(evtPlaneToken_, evtPlanes);
if (evtPlanes.isValid()) {
for (unsigned int i = 0; i < evtPlanes->size(); ++i) {
hiEvtPlane[i] = (*evtPlanes)[i].angle(evtPlaneLevel_);
stepobr marked this conversation as resolved.
Show resolved Hide resolved
}
}
}

edm::Handle<reco::JetView> jets;
if (doJettyExclusion_)
iEvent.getByToken(jetTag_, jets);

int nFill = 0;

double eventPlane2 = -100;
double eventPlane3 = -100;

constexpr int nParamVals = 9;
auto rhoFlowFitParamsOut = std::make_unique<std::vector<double>>(nParamVals, 1e-6);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

use std::array<double, 9> (or, nParamVals) instead to be more explicit about the size

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry can't figure that out, auto rhoFlowFitParamsOut = std::make_unique<std::array<double, nParamVals>> not seem to work

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what was the error?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably just doing something wrong
error: pointer to a function used in arithmetic [-Werror=pointer-arith] rhoFlowFitParamsOut[0] = 0; error: assignment of read-only location '**rhoFlowFitParamsOut' rhoFlowFitParamsOut[0] = 0;

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

but this looks more like about (*rhoFlowFitParamsOut)[0] vs rhoFlowFitParamsOut[0]

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Adding the array to StdDictionaries gives an error during compilation:
Error: It is not necessary to explicitly select class array<double,9>. I/O is supported for it transparently.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The symptoms (ROOT saying it has a dictionary, framework that it has not) smell similar to #27861. I'll take a look.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The cause is different from #27861, and may take some time to resolve. I opened an issue for it #30493.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In case the resolution will take longer than ~today, I'm willing to accept the current vector solution for this product.
With an assumption/note that it can/should be updated once the framework has support for array.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Almost certainly will take longer than today.


rhoFlowFitParamsOut->at(0) = 0;
stepobr marked this conversation as resolved.
Show resolved Hide resolved
rhoFlowFitParamsOut->at(1) = 0;
rhoFlowFitParamsOut->at(2) = 0;
rhoFlowFitParamsOut->at(3) = 0;
rhoFlowFitParamsOut->at(4) = 0;

double eventPlane2Cos = 0;
double eventPlane2Sin = 0;

double eventPlane3Cos = 0;
double eventPlane3Sin = 0;

for (auto const& pfCandidate : *pfCands) {
if (pfCandidate.eta() < -1.0)
continue;
if (pfCandidate.eta() > 1.0)
continue;
if (pfCandidate.pt() < .3)
continue;
if (pfCandidate.pt() > 3.)
continue;
if (pfCandidate.particleId() != 1)
continue;
stepobr marked this conversation as resolved.
Show resolved Hide resolved

if (doJettyExclusion_) {
bool isGood = true;
for (auto const& jet : *jets) {
if (deltaR2(jet, pfCandidate) < .16) {
isGood = false;
break;
}
}
if (!isGood)
continue;
}

nFill++;

if (!doEvtPlane_) {
eventPlane2Cos += std::cos(2 * pfCandidate.phi());
eventPlane2Sin += std::sin(2 * pfCandidate.phi());

eventPlane3Cos += std::cos(3 * pfCandidate.phi());
eventPlane3Sin += std::sin(3 * pfCandidate.phi());
}
}

if (!doEvtPlane_) {
eventPlane2 = std::atan2(eventPlane2Sin, eventPlane2Cos) / 2.;
eventPlane3 = std::atan2(eventPlane3Sin, eventPlane3Cos) / 3.;
} else {
eventPlane2 = hiEvtPlane[8];
eventPlane3 = hiEvtPlane[15];
}

if (nFill >= 100 && eventPlane2 > -99) {
const int nPhiBins = std::max(10, nFill / 30);

std::string name = "phiTestIEta4_" + std::to_string(iEvent.id().event()) + "_h";
std::string nameFlat = "phiTestIEta4_Flat_" + std::to_string(iEvent.id().event()) + "_h";

TH1F* phi_h = new TH1F(name.data(), "", nPhiBins, -TMath::Pi(), TMath::Pi());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that the histogram re-creation per event can be rather heavy.
Moving this to the module creation would be more practical, even considering that we need variants with 3 to 10 bins.
Please add the following to the instances

phi_h->SetDirectory(nullptr);

also, use unique_ptr instead of plain pointers.

Copy link
Contributor Author

@stepobr stepobr Jul 3, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

variants with 3 to 10 bins

it's not like that actually, it's std::max(10, nFill / 30); so not sure that we can avoid histogram re-creation

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's try to be more specific/practical then.
what's the current cost per event?
please get a timing measurement per event (TimeReport is good enough) and a profiler report for CPU https://twiki.cern.ch/twiki/bin/viewauth/CMS/RecoIntegration#Run_profiler_igprof
this should help to understand if something other than TH1 should be used or if TH1 is acceptable.
Thank you.

Copy link
Contributor Author

@stepobr stepobr Jul 3, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's time in seconds averaged over 10 events

hiPuRhoProducer 0.3736367 hiFJRhoFlowModulationProducer 0.001974513 CSJetProducer 0.630398

For profiler report, hiFJRhoFlowModulationProducer is mentioned a couple of times like:
0.0 ......... 0.00 / 0.00 HiFJRhoFlowModulationProducer::produce(edm::Event&, edm::EventSetup const&) [11411]

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thank you.
please share the full report.
Also, please check/clarify what input data is used; it is important to have most events to pass if (nFill >= 100 && eventPlane2 > -99) { for this check to make sense

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

profiler report link
time report link

10 out of 10 events have passed this check, but I'll take a look at a larger sample

stepobr marked this conversation as resolved.
Show resolved Hide resolved

for (auto const& pfCandidate : *pfCands) {
if (pfCandidate.eta() < -1.0)
continue;
if (pfCandidate.eta() > 1.0)
continue;
if (pfCandidate.pt() < .3)
continue;
if (pfCandidate.pt() > 3.)
continue;
if (pfCandidate.particleId() != 1)
continue;

if (doJettyExclusion_) {
bool isGood = true;
for (auto const& jet : *jets) {
if (deltaR2(jet, pfCandidate) < .16) {
isGood = false;
break;
}
}
if (!isGood)
continue;
}
stepobr marked this conversation as resolved.
Show resolved Hide resolved

phi_h->Fill(pfCandidate.phi());
}

std::string flowFitForm = "[0]*(1.+2.*([1]*TMath::Cos(2.*(x-[2]))+[3]*TMath::Cos(3.*(x-[4]))))";

TF1* flowFit_p = new TF1("flowFit", flowFitForm.c_str(), -TMath::Pi(), TMath::Pi());
stepobr marked this conversation as resolved.
Show resolved Hide resolved
flowFit_p->SetParameter(0, 10);
flowFit_p->SetParameter(1, 0);
flowFit_p->SetParameter(2, eventPlane2);
flowFit_p->SetParameter(3, 0);
flowFit_p->SetParameter(4, eventPlane3);

if (!doFreePlaneFit_) {
flowFit_p->FixParameter(2, eventPlane2);
flowFit_p->FixParameter(4, eventPlane3);
}

TF1* lineFit_p = new TF1("lineFit", "[0]", -TMath::Pi(), TMath::Pi());
lineFit_p->SetParameter(0, 10);

phi_h->Fit(flowFit_p, "Q", "", -TMath::Pi(), TMath::Pi());
phi_h->Fit(lineFit_p, "Q", "", -TMath::Pi(), TMath::Pi());
stepobr marked this conversation as resolved.
Show resolved Hide resolved

rhoFlowFitParamsOut->at(0) = flowFit_p->GetParameter(0);
rhoFlowFitParamsOut->at(1) = flowFit_p->GetParameter(1);
rhoFlowFitParamsOut->at(2) = flowFit_p->GetParameter(2);
rhoFlowFitParamsOut->at(3) = flowFit_p->GetParameter(3);
rhoFlowFitParamsOut->at(4) = flowFit_p->GetParameter(4);

rhoFlowFitParamsOut->at(5) = flowFit_p->GetChisquare();
rhoFlowFitParamsOut->at(6) = flowFit_p->GetNDF();

rhoFlowFitParamsOut->at(7) = lineFit_p->GetChisquare();
rhoFlowFitParamsOut->at(8) = lineFit_p->GetNDF();

delete phi_h;

delete flowFit_p;
delete lineFit_p;
}

iEvent.put(std::move(rhoFlowFitParamsOut), "rhoFlowFitParams");
stepobr marked this conversation as resolved.
Show resolved Hide resolved
}

// ------------ method called once each stream just before starting event loop ------------
void HiFJRhoFlowModulationProducer::beginStream(edm::StreamID) {
if (doEvtPlane_) {
constexpr int kMaxEvtPlanes = 1000;
hiEvtPlane = new float[kMaxEvtPlanes];
}
}

// ------------ method called once each stream just after ending the event loop ------------
void HiFJRhoFlowModulationProducer::endStream() {
if (doEvtPlane_)
delete hiEvtPlane;
}

// ------------ method fills 'descriptions' with the allowed parameters for the module ------------
void HiFJRhoFlowModulationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
// The following says we do not know what parameters are allowed so do no validation
// Please change this to state exactly what you do use, even if it is no parameters
edm::ParameterSetDescription desc;
desc.setUnknown();
descriptions.addDefault(desc);
stepobr marked this conversation as resolved.
Show resolved Hide resolved
}

DEFINE_FWK_MODULE(HiFJRhoFlowModulationProducer);
63 changes: 63 additions & 0 deletions RecoHI/HiJetAlgos/plugins/HiPFCandCleaner.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/src/WorkerMaker.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Utilities/interface/EDGetToken.h"
#include <cstdlib>
#include <memory>
#include <string>
#include <utility>
#include <vector>

class HiPFCandCleaner : public edm::global::EDProducer<> {
public:
explicit HiPFCandCleaner(const edm::ParameterSet&);

// class methods

private:
void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;

// ----------member data ---------------------------
edm::EDGetTokenT<reco::PFCandidateCollection> candidatesToken_;

double ptMin_;
double absEtaMax_;
};
//
// constructors and destructor
//
HiPFCandCleaner::HiPFCandCleaner(const edm::ParameterSet& iConfig)
stepobr marked this conversation as resolved.
Show resolved Hide resolved
: candidatesToken_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("candidatesSrc"))),
ptMin_(iConfig.getParameter<double>("ptMin")),
absEtaMax_(iConfig.getParameter<double>("absEtaMax")) {
produces<reco::PFCandidateCollection>("particleFlowCleaned");
}

// ------------ method called to for each event ------------
void HiPFCandCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
edm::Handle<reco::PFCandidateCollection> candidates;
iEvent.getByToken(candidatesToken_, candidates);

auto prod = std::make_unique<reco::PFCandidateCollection>();

for (auto const& cand : *candidates) {
if (cand.pt() < ptMin_)
continue;
if (std::abs(cand.eta()) > absEtaMax_)
continue;
if (cand.particleId() != 1)
continue;

prod->push_back(cand);
}

iEvent.put(std::move(prod), "particleFlowCleaned");
}

DEFINE_FWK_MODULE(HiPFCandCleaner);
Loading