Skip to content

memory leak fixes for Hydjets #37053

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

Merged
cmsbuild merged 1 commit into from
Mar 1, 2022
Merged

memory leak fixes for Hydjets #37053

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
89 changes: 43 additions & 46 deletions GeneratorInterface/Hydjet2Interface/src/Hydjet2Hadronizer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
/**
* \brief Interface to the HYDJET++ (Hydjet2) generator (since core v. 2.4.3), produces HepMC events
* \version 1.1
* \version 1.2
* \author Andrey Belyaev
*/

Expand Down Expand Up @@ -329,29 +329,30 @@ bool Hydjet2Hadronizer::generatePartonsAndHadronize() {
nsub_++;

// event information
HepMC::GenEvent *evt = new HepMC::GenEvent();
std::unique_ptr<HepMC::GenEvent> evt = std::make_unique<HepMC::GenEvent>();
std::unique_ptr<edm::HepMCProduct> HepMCEvt = std::make_unique<edm::HepMCProduct>();

if (nhard_ > 0 || nsoft_ > 0)
get_particles(evt);
get_particles(evt.get());

evt->set_signal_process_id(pypars.msti[0]); // type of the process
evt->set_event_scale(pypars.pari[16]); // Q^2
add_heavy_ion_rec(evt);
add_heavy_ion_rec(evt.get());

if (fVertex_) {
// generate new vertex & apply the shift
// Copy the HepMC::GenEvent
std::unique_ptr<edm::HepMCProduct> HepMCEvt(new edm::HepMCProduct(evt));
HepMCEvt = std::make_unique<edm::HepMCProduct>(evt.get());
HepMCEvt->applyVtxGen(fVertex_);
evt = new HepMC::GenEvent((*HepMCEvt->GetEvent()));
evt = std::make_unique<HepMC::GenEvent>(*HepMCEvt->GetEvent());
}

HepMC::HEPEVT_Wrapper::check_hepevt_consistency();
LogDebug("HEPEVT_info") << "Ev numb: " << HepMC::HEPEVT_Wrapper::event_number()
<< " Entries number: " << HepMC::HEPEVT_Wrapper::number_entries() << " Max. entries "
<< HepMC::HEPEVT_Wrapper::max_number_entries() << std::endl;

event().reset(evt);
event() = std::move(evt);
return kTRUE;
}

Expand Down Expand Up @@ -385,83 +386,79 @@ void Hydjet2Hadronizer::rotateEvtPlane() {

//_____________________________________________________________________
bool Hydjet2Hadronizer::get_particles(HepMC::GenEvent *evt) {
LogDebug("SubEvent") << " Number of sub events " << nsub_;
LogDebug("Hydjet") << " Number of hard events " << hj2->GetNjet();
LogDebug("Hydjet") << " Number of hard particles " << nhard_;
LogDebug("Hydjet") << " Number of soft particles " << nsoft_;
LogDebug("Hydjet") << " nhard_ + nsoft_ = " << nhard_ + nsoft_ << " Ntot = " << hj2->GetNtot() << endl;
LogDebug("Hydjet2") << " Number of sub events " << nsub_;
LogDebug("Hydjet2") << " Number of hard events " << hj2->GetNjet();
LogDebug("Hydjet2") << " Number of hard particles " << nhard_;
LogDebug("Hydjet2") << " Number of soft particles " << nsoft_;
LogDebug("Hydjet2") << " nhard_ + nsoft_ = " << nhard_ + nsoft_ << " Ntot = " << hj2->GetNtot() << endl;

int ihy = 0;
int isub_l = -1;
int stab = 0;

vector<HepMC::GenParticle *> primary_particle(hj2->GetNtot());
vector<HepMC::GenParticle *> particle(hj2->GetNtot());

HepMC::GenVertex *sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), 0); // just initialization
HepMC::GenVertex *sub_vertices = nullptr;

while (ihy < hj2->GetNtot()) {
if ((hj2->GetiJet().at(ihy)) != isub_l) {
sub_vertices = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, 0), hj2->GetiJet().at(ihy));

evt->add_vertex(sub_vertices);

if (!evt->signal_process_vertex())
evt->set_signal_process_vertex(sub_vertices);
isub_l = hj2->GetiJet().at(ihy);
}

if ((hj2->GetFinal().at(ihy)) == 1) //convertStatus(hj2->GetPythiaStatus().at(ihy)) == 1)
stab++;
LogDebug("Hydjet_array") << ihy << " MULTin ev.:" << hj2->GetNtot() << " SubEv.#" << hj2->GetiJet().at(ihy)
<< " Part #" << ihy + 1 << ", PDG: " << hj2->GetPdg().at(ihy) << " (st. "
<< convertStatus(hj2->GetPythiaStatus().at(ihy))
<< ") mother=" << hj2->GetMotherIndex().at(ihy) + 1 << ", childs ("
<< hj2->GetFirstDaughterIndex().at(ihy) + 1 << "-"
<< hj2->GetLastDaughterIndex().at(ihy) + 1 << "), vtx (" << hj2->GetX().at(ihy) << ","
<< hj2->GetY().at(ihy) << "," << hj2->GetZ().at(ihy) << ") " << std::endl;
LogDebug("Hydjet2_array") << ihy << " MULTin ev.:" << hj2->GetNtot() << " SubEv.#" << hj2->GetiJet().at(ihy)
<< " Part #" << ihy + 1 << ", PDG: " << hj2->GetPdg().at(ihy) << " (st. "
<< convertStatus(hj2->GetPythiaStatus().at(ihy))
<< ") mother=" << hj2->GetMotherIndex().at(ihy) + 1 << ", childs ("
<< hj2->GetFirstDaughterIndex().at(ihy) + 1 << "-"
<< hj2->GetLastDaughterIndex().at(ihy) + 1 << "), vtx (" << hj2->GetX().at(ihy) << ","
<< hj2->GetY().at(ihy) << "," << hj2->GetZ().at(ihy) << ") " << std::endl;

if ((hj2->GetMotherIndex().at(ihy)) <= 0) {
primary_particle.at(ihy) = build_hyjet2(ihy, ihy + 1);
sub_vertices->add_particle_out(primary_particle.at(ihy));
LogDebug("Hydjet_array") << " ---> " << ihy + 1 << std::endl;
particle.at(ihy) = build_hyjet2(ihy, ihy + 1);
if (!sub_vertices)
LogError("Hydjet2_array") << "##### HYDJET2: Vertex not initialized!";
else
sub_vertices->add_particle_out(particle.at(ihy));
LogDebug("Hydjet2_array") << " ---> " << ihy + 1 << std::endl;
} else {
particle.at(ihy) = build_hyjet2(ihy, ihy + 1);
int mid = hj2->GetMotherIndex().at(ihy);
int mid_t = mid;

while ((mid < ihy) && ((hj2->GetPdg().at(ihy)) < 100) && ((hj2->GetFirstDaughterIndex().at(ihy)) == ihy))
while (((mid + 1) < ihy) && (std::abs(hj2->GetPdg().at(mid)) < 100) &&
((hj2->GetFirstDaughterIndex().at(mid + 1)) <= ihy)) {
mid++;
LogDebug("Hydjet2_array") << "======== MID changed to " << mid
<< " ======== PDG(mid) = " << hj2->GetPdg().at(mid) << std::endl;
}

if ((hj2->GetPdg().at(ihy)) < 100)
mid = mid_t;

HepMC::GenParticle *mother = primary_particle.at(mid);
HepMC::GenVertex *prods = build_hyjet2_vertex(ihy, (hj2->GetiJet().at(ihy)));

if (!mother) {
mother = particle.at(mid);
primary_particle.at(mid) = mother;
if (std::abs(hj2->GetPdg().at(mid)) < 100) {
mid = hj2->GetMotherIndex().at(ihy);
LogDebug("Hydjet2_array") << "======== MID changed BACK to " << mid
<< " ======== PDG(mid) = " << hj2->GetPdg().at(mid) << std::endl;
}

HepMC::GenParticle *mother = particle.at(mid);
HepMC::GenVertex *prod_vertex = mother->end_vertex();

if (!prod_vertex) {
prod_vertex = prods;
prod_vertex = build_hyjet2_vertex(ihy, (hj2->GetiJet().at(ihy)));
prod_vertex->add_particle_in(mother);
LogDebug("Hydjet_array") << " <--- " << mid + 1 << std::endl;
LogDebug("Hydjet2_array") << " <--- " << mid + 1 << std::endl;
evt->add_vertex(prod_vertex);
prods = nullptr;
}
prod_vertex->add_particle_out(particle.at(ihy));
LogDebug("Hydjet_array") << " ---" << mid + 1 << "---> " << ihy + 1 << std::endl;
delete prods;
LogDebug("Hydjet2_array") << " ---" << mid + 1 << "---> " << ihy + 1 << std::endl;
}
ihy++;
}
delete sub_vertices;

LogDebug("Hydjet_array") << " MULTin ev.:" << hj2->GetNtot() << ", last index: " << ihy - 1
<< ", stable particles: " << stab << std::endl;
LogDebug("Hydjet2_array") << " MULTin ev.:" << hj2->GetNtot() << ", last index: " << ihy - 1
<< ", stable particles: " << stab << std::endl;
return kTRUE;
}

Expand Down
50 changes: 28 additions & 22 deletions GeneratorInterface/Hydjet2Interface/test/Hydjet2Analyzer.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/**
\class Hydjet2Analyzer
\brief HepMC events analyzer
\version 2.0
\version 2.1
\authors Yetkin Yilmaz, Andrey Belyaev
*/

Expand All @@ -13,7 +13,7 @@
#include <string>

// user include files
#include "FWCore/Framework/interface/EDAnalyzer.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"

#include "FWCore/Framework/interface/Event.h"
Expand Down Expand Up @@ -42,16 +42,18 @@
// root include file
#include "TFile.h"
#include "TNtuple.h"
#include "TH1.h"

#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

using namespace edm;
using namespace std;
static const int MAXPARTICLES = 5000000;
static const int ETABINS = 3; // Fix also in branch string

namespace {
static const int MAXPARTICLES = 5000000;
static const int ETABINS = 3; // Fix also in branch string
} // namespace
struct Hydjet2Event {
int event;
float b;
Expand Down Expand Up @@ -79,16 +81,15 @@ struct Hydjet2Event {
float vz;
float vr;
};
class Hydjet2Analyzer : public edm::EDAnalyzer {
class Hydjet2Analyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
public:
explicit Hydjet2Analyzer(const edm::ParameterSet &);
~Hydjet2Analyzer();

private:
virtual void beginRun(const edm::Run &, const edm::EventSetup &);
virtual void beginJob();
virtual void analyze(const edm::Event &, const edm::EventSetup &);
virtual void endJob();
void beginJob() final;
void analyze(const edm::Event &, const edm::EventSetup &) final;
void endJob() final;
// ----------member data ---------------------------
std::ofstream out_b;
std::string fBFileName;
Expand All @@ -100,6 +101,7 @@ class Hydjet2Analyzer : public edm::EDAnalyzer {
Hydjet2Event hev_;
TNtuple *nt;
std::string output; // Output filename
bool doTestEvent_;
bool doAnalysis_;
bool printLists_;
bool doCF_;
Expand Down Expand Up @@ -155,9 +157,7 @@ class Hydjet2Analyzer : public edm::EDAnalyzer {
edm::InputTag genParticleSrc_;
edm::InputTag genHIsrc_;
edm::InputTag simVerticesTag_;
edm::ESGetToken<HepPDT::ParticleDataTable, PDTRecord> pdt_;

edm::Service<TFileService> f;
edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdtToken_;

//common
TH1D *dhpt;
Expand Down Expand Up @@ -272,7 +272,6 @@ class Hydjet2Analyzer : public edm::EDAnalyzer {
};

Hydjet2Analyzer::Hydjet2Analyzer(const edm::ParameterSet &iConfig) {
pdt_ = esConsumes<HepPDT::ParticleDataTable, PDTRecord>();
fBFileName = iConfig.getUntrackedParameter<std::string>("output_b", "b_values.txt");
fNFileName = iConfig.getUntrackedParameter<std::string>("output_n", "n_values.txt");
fMFileName = iConfig.getUntrackedParameter<std::string>("output_m", "m_values.txt");
Expand All @@ -293,10 +292,12 @@ Hydjet2Analyzer::Hydjet2Analyzer(const edm::ParameterSet &iConfig) {

genParticleSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("hiGenParticles"));
genHIsrc_ = iConfig.getUntrackedParameter<edm::InputTag>("src", edm::InputTag("heavyIon"));
doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", false);

if (useHepMCProduct_)
pdtToken_ = esConsumes();
doTestEvent_ = iConfig.getUntrackedParameter<bool>("doTestEvent", false);
doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", false);
doHistos_ = iConfig.getUntrackedParameter<bool>("doHistos", false);

if (doHistos_) {
userHistos_ = iConfig.getUntrackedParameter<bool>("userHistos", false);
if (userHistos_) {
Expand Down Expand Up @@ -408,7 +409,6 @@ Hydjet2Analyzer::~Hydjet2Analyzer() {}

// ------------ method called to for each event ------------
void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
auto const &pdt = iSetup.getData(pdt_);
using namespace edm;
using namespace HepMC;
hev_.event = iEvent.id().event();
Expand All @@ -435,7 +435,8 @@ void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &i
int sig = -1;
int src = -1;
if (useHepMCProduct_) {
//_______________________________________________________________________________________
edm::ESHandle<ParticleDataTable> pdt = iSetup.getHandle(pdtToken_);

if (doCF_) {
Handle<CrossingFrame<HepMCProduct>> cf;
iEvent.getByToken(srcTmix_, cf);
Expand All @@ -457,7 +458,7 @@ void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &i
float eta = (*it)->momentum().eta();
float phi = (*it)->momentum().phi();
float pt = (*it)->momentum().perp();
const ParticleData *part = pdt.particle(pdg_id);
const ParticleData *part = pdt->particle(pdg_id);
int charge = static_cast<int>(part->charge());
hev_.pt[hev_.mult] = pt;
hev_.eta[hev_.mult] = eta;
Expand Down Expand Up @@ -502,6 +503,10 @@ void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &i
}

src = evt->particles_size();
if (doTestEvent_) {
std::cout << "Event size: " << src << std::endl;
mc->GetEvent()->print();
}

HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
Expand All @@ -519,8 +524,8 @@ void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &i
float e = (*it)->momentum().e();
float pseudoRapidity = (*it)->momentum().pseudoRapidity();
int charge = -1;
if ((pdt.particle(pdg_id))) {
part = pdt.particle(pdg_id);
if ((pdt->particle(pdg_id))) {
part = pdt->particle(pdg_id);
charge = static_cast<int>(part->charge());
}

Expand Down Expand Up @@ -795,7 +800,6 @@ void Hydjet2Analyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &i
}
}
// ------------ method called once each job just before starting event loop ------------
void Hydjet2Analyzer::beginRun(const edm::Run &, const edm::EventSetup &iSetup) {}
void Hydjet2Analyzer::beginJob() {
if (printLists_) {
out_b.open(fBFileName.c_str());
Expand Down Expand Up @@ -944,6 +948,8 @@ void Hydjet2Analyzer::beginJob() {
}

if (doAnalysis_) {
usesResource(TFileService::kSharedResource);
edm::Service<TFileService> f;
nt = f->make<TNtuple>("nt", "Mixing Analysis", "mix:np:src:sig");
hydjetTree_ = f->make<TTree>("hi", "Tree of Hydjet2 Events");
hydjetTree_->Branch("event", &hev_.event, "event/I");
Expand Down
37 changes: 35 additions & 2 deletions GeneratorInterface/Hydjet2Interface/test/testHydjet.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,14 @@
process.load("Configuration.StandardSequences.Services_cff")
process.load("GeneratorInterface.Hydjet2Interface.hydjet2Default_cfi")

process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(100))
process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(10))

process.ana = cms.EDAnalyzer('Hydjet2Analyzer',

doHistos = cms.untracked.bool(True),
userHistos = cms.untracked.bool(False),
#doAnalysis = cms.untracked.bool(True),
doAnalysis = cms.untracked.bool(False),
doTestEvent = cms.untracked.bool(False), # for debuging event output information

###Settings for USER histos

Expand Down Expand Up @@ -58,6 +59,38 @@

#to separate hydro and jet parts of hydjet2
process.generator.separateHydjetComponents = cms.untracked.bool(True)
Debug = None

if Debug:
process.load("FWCore.MessageLogger.MessageLogger_cfi")

process.MessageLogger = cms.Service("MessageLogger",

destinations = cms.untracked.vstring('LogDebug_Hydjet2'),
categories = cms.untracked.vstring(
'Hydjet2',
'Hydjet2_array'
),
LogDebug_Hydjet2 = cms.untracked.PSet(
threshold = cms.untracked.string('DEBUG'), #Priority: DEBUG < INFO < WARNING < ERROR
DEBUG = cms.untracked.PSet(limit = cms.untracked.int32(-1)),
INFO = cms.untracked.PSet(limit = cms.untracked.int32(0)),
WARNING = cms.untracked.PSet(limit = cms.untracked.int32(0)),
ERROR = cms.untracked.PSet(limit = cms.untracked.int32(0)),
#Categores
Hydjet2 = cms.untracked.PSet(
limit = cms.untracked.int32(-1), # number of masseges
timespan = cms.untracked.int32(0) #time to resete limit counter in seconds
),

Hydjet2_array = cms.untracked.PSet(
limit = cms.untracked.int32(-1), # number of masseges
timespan = cms.untracked.int32(0) #time to resete limit counter in seconds
)

),
debugModules = cms.untracked.vstring('*')
)

process.TFileService = cms.Service('TFileService',
fileName = cms.string('Hydjet2_MB_5020GeV.root')
Expand Down
Loading