Skip to content
Snippets Groups Projects
Commit 8d622b5a authored by Frank Winklmeier's avatar Frank Winklmeier
Browse files

Merge branch 'hepmc3_EvgenProdTools' into 'master' 

Migration of  EvgenProdTools, first part

See merge request !34341
parents fd0ac69a fdf140f0
No related branches found
No related tags found
6 merge requests!58791DataQualityConfigurations: Modify L1Calo config for web display,!46784MuonCondInterface: Enable thread-safety checking.,!46776Updated LArMonitoring config file for WD to match new files produced using MT,!45405updated ART test cron job,!42417Draft: DIRE and VINCIA Base Fragments for Pythia 8.3,!34341Migration of EvgenProdTools, first part
Hide whitespace changes
Inline Side-by-side
Generators/EvgenProdTools/src/FixHepMC.cxx
+ 9
− 0
View file @ 8d622b5a
......@@ -151,6 +151,14 @@ bool FixHepMC::isLoop(const HepMC::GenParticlePtr p) {
if (p->production_vertex() == p->end_vertex() && p->end_vertex() != NULL) return true;
if (m_loopByBC && p->production_vertex()) {
/// @todo Use new particle MC::parents(...) tool
#ifdef HEPMC3
for (auto itrParent: p->production_vertex()->particles_out()) {
if ( HepMC::barcode(itrParent) > HepMC::barcode(p) ) {
ATH_MSG_VERBOSE("Found a loop (a la Sherpa sample) via barcode.");
return true; // Cannot vectorize, but this is a pretty short loop
} // Check on barcodes
} // Loop over parent particles
#else
for (HepMC::GenVertex::particle_iterator itrParent = p->production_vertex()->particles_begin(HepMC::parents);
itrParent != p->production_vertex()->particles_end(HepMC::parents); ++itrParent) {
if ( (*itrParent)->barcode() > p->barcode() ) {
......@@ -158,6 +166,7 @@ bool FixHepMC::isLoop(const HepMC::GenParticlePtr p) {
return true; // Cannot vectorize, but this is a pretty short loop
} // Check on barcodes
} // Loop over parent particles
#endif
} // Has a production vertex
return false;
}
......
Generators/EvgenProdTools/src/TestHepMC.cxx
+ 89
− 63
View file @ 8d622b5a
......@@ -278,7 +278,14 @@ StatusCode TestHepMC::execute() {
vector<int> undisplaceds;
// Check beams and work out per-event beam energy
const pair<HepMC::GenParticlePtr, HepMC::GenParticlePtr> beams = evt->beam_particles();
#ifdef HEPMC3
auto beams_t = evt->beams();
std::pair<std::shared_ptr<const HepMC3::GenParticle>,std::shared_ptr<const HepMC3::GenParticle>> beams;
beams.first=beams_t.at(0);
beams.second=beams_t.at(1);
#else
auto beams = evt->beam_particles();
#endif
double cmenergy = m_cm_energy;
if (!HepMC::valid_beam_particles(evt)) {
ATH_MSG_WARNING("Invalid beam particle pointers -- this generator interface should be fixed");
......@@ -291,7 +298,7 @@ StatusCode TestHepMC::execute() {
}
const double sumE = beams.first->momentum().e() + beams.second->momentum().e();
const double sumP = beams.first->momentum().pz() + beams.second->momentum().pz();
cmenergy = sqrt(sumE*sumE - sumP*sumP);
cmenergy = std::sqrt(sumE*sumE - sumP*sumP);
if (m_cm_energy > 0 && fabs(cmenergy - m_cm_energy) > m_cme_diff) {
ATH_MSG_FATAL("Beam particles have incorrect energy: " << m_cm_energy/1000. << " GeV expected, vs. " << cmenergy/1000. << " GeV found");
......@@ -310,8 +317,12 @@ StatusCode TestHepMC::execute() {
// Check vertices
int vtxDisplacedstatuscode12CheckRateCnt=0;
int vtxDisplacedstatuscodenot12CheckRateCnt=0;
for (HepMC::GenEvent::vertex_const_iterator vitr = evt->vertices_begin(); vitr != evt->vertices_end(); ++vitr ) {
#ifdef HEPMC3
for (auto vtx: evt->vertices()) {
#else
for (auto vitr = evt->vertices_begin(); vitr != evt->vertices_end(); ++vitr ) {
const HepMC::GenVertex* vtx = *vitr;
#endif
const HepMC::FourVector pos = vtx->position();
// Check for NaNs and infs in vertex position components
......@@ -331,8 +342,8 @@ StatusCode TestHepMC::execute() {
// Anything which propagates macroscopically should be set stable in evgen for G4 to handle
const double dist_trans2 = pos.x()*pos.x() + pos.y()*pos.y(); // in mm2
const double dist2 = dist_trans2 + pos.z()*pos.z(); // in mm2
const double dist_trans = sqrt(dist_trans2); // in mm
const double dist = sqrt(dist2); // in mm
const double dist_trans = std::sqrt(dist_trans2); // in mm
const double dist = std::sqrt(dist2); // in mm
if (dist2 > m_max_dist*m_max_dist) {
ATH_MSG_WARNING("Found vertex position displaced by more than " << m_max_dist << "mm: " << dist << "mm");
++m_vtxDisplacedMoreThan_1m_CheckRateCnt;
......@@ -344,50 +355,59 @@ StatusCode TestHepMC::execute() {
if (dist_trans2 > m_max_dist_trans*m_max_dist_trans) {
ATH_MSG_WARNING("Found vertex position displaced by more than " << m_max_dist_trans << "mm in transverse distance: " << dist_trans << "mm");
for (auto par = vtx->particles_in_const_begin(); par != vtx->particles_in_const_end(); ++par) {
#ifdef HEPMC3
for (auto part: vtx->particles_in()) {
#else
for (auto part_it = vtx->particles_in_const_begin(); part_it != vtx->particles_in_const_end(); ++part_it) {
auto part=(*part_it);
#endif
ATH_MSG_WARNING("Outgoing particle : ");
if (m_dumpEvent) HepMC::Print::line(std::cout,*par);
ATH_MSG_WARNING("production vertex = " << (*par)->production_vertex()->position().x() << ", " << (*par)->production_vertex()->position().y() << ", " << (*par)->production_vertex()->position().z());
ATH_MSG_WARNING("end vertex = " << (*par)->end_vertex()->position().x() << ", " << (*par)->end_vertex()->position().y() << ", " << (*par)->end_vertex()->position().z());
if (m_dumpEvent) HepMC::Print::line(std::cout,part);
ATH_MSG_WARNING("production vertex = " << part->production_vertex()->position().x() << ", " << part->production_vertex()->position().y() << ", " << part->production_vertex()->position().z());
ATH_MSG_WARNING("end vertex = " << part->end_vertex()->position().x() << ", " << part->end_vertex()->position().y() << ", " << part->end_vertex()->position().z());
ATH_MSG_WARNING("parents info: ");
if ((*par)->production_vertex()) {
for(auto p_parents = (*par)->production_vertex()->particles_in_const_begin(); p_parents != (*par)->production_vertex()->particles_in_const_end(); ++p_parents) {
// cout << "\t";
if (m_dumpEvent) HepMC::Print::line(std::cout,*p_parents);
if (part->production_vertex()) {
#ifdef HEPMC3
for(auto p_parents: part->production_vertex()->particles_in()) {
#else
for(auto p_parents_it = part->production_vertex()->particles_in_const_begin(); p_parents_it != part->production_vertex()->particles_in_const_end(); ++p_parents_it) {
auto p_parents=(*p_parents_it);
#endif
if (m_dumpEvent) HepMC::Print::line(std::cout,p_parents);
ATH_MSG_WARNING("\t");
}
} // Done with fancy print
if ((*par)->status()==1 || (*par)->status()==2){
if (part->status()==1 || part->status()==2){
vtxDisplacedstatuscode12CheckRateCnt += 1;
} else {
vtxDisplacedstatuscodenot12CheckRateCnt += 1;
}
if (m_doHist){
m_h_energy_dispVtxCheck->Fill((*par)->momentum().e()*1e-3);
if ((*par)->momentum().e()*1e-3 < 10.) {
m_h_energy_dispVtxCheck_lt10->Fill((*par)->momentum().e()*1e-3);
m_h_energy_dispVtxCheck->Fill(part->momentum().e()*1e-3);
if (part->momentum().e()*1e-3 < 10.) {
m_h_energy_dispVtxCheck_lt10->Fill(part->momentum().e()*1e-3);
}
m_h_pdgid_dispVtxCheck->Fill((*par)->pdg_id());
m_h_status_dispVtxCheck->Fill((*par)->status());
m_h_px_dispVtxCheck->Fill((*par)->momentum().px()*1e-3);
m_h_py_dispVtxCheck->Fill((*par)->momentum().py()*1e-3);
m_h_pz_dispVtxCheck->Fill((*par)->momentum().pz()*1e-3);
m_h_vx_dispVtxCheck->Fill((*par)->end_vertex()->position().x());
m_h_vy_dispVtxCheck->Fill((*par)->end_vertex()->position().y());
m_h_vz_dispVtxCheck->Fill((*par)->end_vertex()->position().z());
m_h_vxprod_dispVtxCheck->Fill((*par)->production_vertex()->position().x());
m_h_vyprod_dispVtxCheck->Fill((*par)->production_vertex()->position().y());
m_h_vzprod_dispVtxCheck->Fill((*par)->production_vertex()->position().z());
double endvx = (*par)->end_vertex()->position().x();
double endvy = (*par)->end_vertex()->position().y();
double endvz = (*par)->end_vertex()->position().z();
double prodvx = (*par)->production_vertex()->position().x();
double prodvy = (*par)->production_vertex()->position().y();
double prodvz = (*par)->production_vertex()->position().z();
double enddis = sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
double proddis = sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
m_h_pdgid_dispVtxCheck->Fill(part->pdg_id());
m_h_status_dispVtxCheck->Fill(part->status());
m_h_px_dispVtxCheck->Fill(part->momentum().px()*1e-3);
m_h_py_dispVtxCheck->Fill(part->momentum().py()*1e-3);
m_h_pz_dispVtxCheck->Fill(part->momentum().pz()*1e-3);
m_h_vx_dispVtxCheck->Fill(part->end_vertex()->position().x());
m_h_vy_dispVtxCheck->Fill(part->end_vertex()->position().y());
m_h_vz_dispVtxCheck->Fill(part->end_vertex()->position().z());
m_h_vxprod_dispVtxCheck->Fill(part->production_vertex()->position().x());
m_h_vyprod_dispVtxCheck->Fill(part->production_vertex()->position().y());
m_h_vzprod_dispVtxCheck->Fill(part->production_vertex()->position().z());
double endvx = part->end_vertex()->position().x();
double endvy = part->end_vertex()->position().y();
double endvz = part->end_vertex()->position().z();
double prodvx = part->production_vertex()->position().x();
double prodvy = part->production_vertex()->position().y();
double prodvz = part->production_vertex()->position().z();
double enddis = std::sqrt(endvx*endvx + endvy*endvy + endvz*endvz);
double proddis = std::sqrt(prodvx*prodvx + prodvy*prodvy + prodvz*prodvz);
m_h_vtxend_dispVtxCheck->Fill(enddis);
m_h_vtxprod_dispVtxCheck->Fill(proddis);
} // End of the filling of histograms for bad vertices
......@@ -398,13 +418,13 @@ StatusCode TestHepMC::execute() {
if (vtxDisplacedstatuscodenot12CheckRateCnt>0) ++m_vtxDisplacedstatuscodenot12CheckRate;
// Check particles
for (auto pitr = evt->particles_begin(); pitr != evt->particles_end(); ++pitr ) {
for (auto pitr: *evt) {
// Local loop variables to clean up the check code
const HepMC::FourVector pmom = (*pitr)->momentum();
const int pstatus = (*pitr)->status();
const int ppdgid = (*pitr)->pdg_id();
const int pbarcode = HepMC::barcode(*pitr);
const HepMC::FourVector pmom = pitr->momentum();
const int pstatus = pitr->status();
const int ppdgid = pitr->pdg_id();
const int pbarcode = HepMC::barcode(pitr);
// Check for NaNs and infs in momentum components
if ( std::isnan(pmom.px()) || std::isinf(pmom.px()) ||
......@@ -414,7 +434,7 @@ StatusCode TestHepMC::execute() {
ATH_MSG_WARNING("NaN (Not A Number) or inf found in the event record momenta");
++m_partMomentumNANandINFCheckRate;
if (m_dumpEvent) HepMC::Print::line(std::cout,*pitr);
if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
if (m_momNaNTest) {
filter_pass = false;
}
......@@ -422,7 +442,7 @@ StatusCode TestHepMC::execute() {
// Check for undecayed pi0s
if (pstatus == 1 || pstatus == 2) {
if (ppdgid == 111 && !(*pitr)->end_vertex() ) {
if (ppdgid == 111 && !pitr->end_vertex() ) {
unDecPi0.push_back( pbarcode);
++m_undecayedPi0CheckRate;
}
......@@ -435,7 +455,7 @@ StatusCode TestHepMC::execute() {
double plifetime = pd->lifetime()*1e+12; // why lifetime doesn't come in common units???
if (plifetime != 0 && plifetime < m_min_tau) { // particles with infinite lifetime get a 0 in the PDT
ATH_MSG_WARNING("Stable particle found with lifetime = " << plifetime << "~ns!!");
if (m_dumpEvent) HepMC::Print::line(std::cout,*pitr);
if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
++m_Status1ShortLifetime;
......@@ -449,14 +469,14 @@ StatusCode TestHepMC::execute() {
vector<int>::size_type count = 0;
while (susyPart==0 && (count < m_SusyPdgID_tab.size() )){
// no warning for SUSY particles from the list susyParticlePdgid.txt
if (m_SusyPdgID_tab[count] == abs(ppdgid)) {
if (m_SusyPdgID_tab[count] == std::abs(ppdgid)) {
susyPart=1;
}
count++;
} // Look through the SUSY table to see if this one should be counted
if (susyPart==0){
ATH_MSG_WARNING("Stable particle not found in PDT, no lifetime check done");
if (m_dumpEvent) HepMC::Print::line(std::cout,*pitr);
if (m_dumpEvent) HepMC::Print::line(std::cout,pitr);
} // It's a SUSY particle -- skip the lifetime check
} // The particle has no data table
} // Test if the particle is stable
......@@ -466,7 +486,7 @@ StatusCode TestHepMC::execute() {
int first_dig = ppdgid;
while (first_dig > 9) first_dig /= 10;
if ((pstatus == 1 ) && (!(*pitr)->end_vertex()) && (!m_nonint.operator()(*pitr)) && (!pid.isNucleus()) && (first_dig != 9) ) {
if ((pstatus == 1 ) && (!pitr->end_vertex()) && (!m_nonint.operator()(pitr)) && (!pid.isNucleus()) && (first_dig != 9) ) {
int known_byG4 = 0;
vector<int>::size_type count =0;
......@@ -490,14 +510,14 @@ StatusCode TestHepMC::execute() {
} // End of check for invalid PDG IDs
// Check for unstables with no end vertex,
if (!(*pitr)->end_vertex() && pstatus == 2) {
if (!pitr->end_vertex() && pstatus == 2) {
unstNoEnd.push_back(pbarcode);
++m_unstableNoEndVtxCheckRate;
} // End of check for unstable with no end vertex
// Sum final state mom/energy, and note negative energy / tachyonic particles
// std::cout << "status " << pstatus << " e " << pmom.e() << " pz " << pmom.pz()<< std::endl;
if ( pstatus == 1 && !(*pitr)->end_vertex() ) {
if ( pstatus == 1 && !pitr->end_vertex() ) {
totalPx += pmom.px();
totalPy += pmom.py();
totalPz += pmom.pz();
......@@ -507,7 +527,7 @@ StatusCode TestHepMC::execute() {
++m_negativeEnergyTachyonicCheckRate;
}
const double aener = fabs(pmom.e());
if ( aener+m_accur_margin < fabs(pmom.px()) || aener+m_accur_margin < fabs(pmom.py()) || aener+m_accur_margin < fabs(pmom.pz()) ) {
if ( aener+m_accur_margin < std::abs(pmom.px()) || aener+m_accur_margin < std::abs(pmom.py()) || aener+m_accur_margin < std::abs(pmom.pz()) ) {
tachyons.push_back(pbarcode);
++m_negativeEnergyTachyonicCheckRate;
}
......@@ -516,16 +536,17 @@ StatusCode TestHepMC::execute() {
// Decay checks (uses PdgToSearch attr value, for tau by default)
/// @todo Clean up / improve / apply to *all* decaying species
int tau_child = 0;
if (abs(ppdgid) == m_pdg && (pstatus == 1 || pstatus == 2)) {
if (std::abs(ppdgid) == m_pdg && (pstatus == 1 || pstatus == 2)) {
++m_TotalTaus;
auto vtx = (*pitr)->end_vertex();
auto vtx = pitr->end_vertex();
if (vtx) {
double p_energy = 0;
for (auto desc = vtx->particles_begin(HepMC::descendants); desc != vtx->particles_end(HepMC::descendants); ++desc) {
if (abs((*desc)->pdg_id()) == m_pdg) tau_child = 1;
if ((*desc)->status() == 1) p_energy += (*desc)->momentum().e();
for (auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
auto desc=(*desc_it);
if (std::abs(desc->pdg_id()) == m_pdg) tau_child = 1;
if (desc->status() == 1) p_energy += desc->momentum().e();
}
if (fabs( p_energy - pmom.e()) > m_energy_diff && !tau_child) {
if (std::abs( p_energy - pmom.e()) > m_energy_diff && !tau_child) {
ATH_MSG_WARNING("Energy sum (decay products): "
<< "Energy (original particle) > " << m_energy_diff << " MeV, "
<< "Event #" << evt->event_number() << ", "
......@@ -546,19 +567,24 @@ StatusCode TestHepMC::execute() {
} // End of checks for specific particle (tau by default)
// Check for undisplaced decay daughters from long-lived hadrons
if ((*pitr)->end_vertex()) {
auto decayvtx = (*pitr)->end_vertex();
if (pitr->end_vertex()) {
auto decayvtx = pitr->end_vertex();
const HepMC::FourVector decaypos = decayvtx->position();
const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
if (displacement > 1e-6) {
#ifdef HEPMC3
for (auto ip:decayvtx->particles_out()) {
#else
HepMC::GenVertex::particle_iterator ipbegin = decayvtx->particles_begin(HepMC::children);
HepMC::GenVertex::particle_iterator ipend = decayvtx->particles_end(HepMC::children);
for (HepMC::GenVertex::particle_iterator ip = ipbegin; ip != ipend; ++ip) {
const HepMC::FourVector pos2 = (*ip)->production_vertex()->position();
for (HepMC::GenVertex::particle_iterator ip_it = ipbegin; ip_it != ipend; ++ip_it) {
auto ip=(*ip_it);
#endif
const HepMC::FourVector pos2 = ip->production_vertex()->position();
const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
if (displacement2 < 1e-6) {
const int pbarcode2 = HepMC::barcode(*ip);
const int vbarcode2 = HepMC::barcode((*ip)->production_vertex());
const int pbarcode2 = HepMC::barcode(ip);
const int vbarcode2 = HepMC::barcode(ip->production_vertex());
ATH_MSG_WARNING("Decay child " << pbarcode2 << " from " << pbarcode
<< " has undisplaced vertex (" << vbarcode2
<< " @ " << displacement2 << "mm) "
......@@ -574,7 +600,7 @@ StatusCode TestHepMC::execute() {
// Check for photons with non-zero masses
/// @todo Persuade generator authors to set proper generated masses in HepMC, then *really* require mass = 0
if (ppdgid == 22 && pstatus == 1) {
const double mass = (*pitr)->generated_mass();
const double mass = pitr->generated_mass();
if (fabs(mass) > 1.0) { // in MeV
ATH_MSG_WARNING("Photon with non-zero mass found! Mass: " << mass << " MeV, BARCODE=" << pbarcode);
++m_nonZeroPhotonMassCheckRate;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment