Commit 38f79cdd authored by Ewelina Maria Lobodzinska's avatar Ewelina Maria Lobodzinska
Browse files

Merge branch '21.6-direct-photon-filter-remove-overlap' into '21.6'

Add option to avoid overlap in DirectPhotonFilter slices

See merge request !27585
parents 9e2f3334 54c05cef
......@@ -13,15 +13,17 @@ class DirectPhotonFilter : public GenFilter {
public:
DirectPhotonFilter(const std::string& name, ISvcLocator* pSvcLocator);
virtual StatusCode filterInitialize();
virtual StatusCode filterEvent();
private:
double m_Ptmin;
double m_Ptmax;
std::vector<double> m_Ptmin;
std::vector<double> m_Ptmax;
double m_EtaRange;
int m_NPhotons;
size_t m_NPhotons;
bool m_AllowSUSYDecay;
bool m_OrderPhotons;
};
......
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/DirectPhotonFilter.h"
#include <limits>
#include <algorithm>
DirectPhotonFilter::DirectPhotonFilter(const std::string& name, ISvcLocator* pSvcLocator)
: GenFilter(name, pSvcLocator)
{
declareProperty("Ptmin",m_Ptmin = 10000.);
declareProperty("Ptmax",m_Ptmax = 100000000.);
declareProperty("Etacut", m_EtaRange = 2.50);
declareProperty("NPhotons", m_NPhotons = 1);
declareProperty("OrderPhotons",m_OrderPhotons = true);
declareProperty("Ptmin",m_Ptmin = std::vector<double>(m_NPhotons, 10000.));
declareProperty("Ptmax",m_Ptmax = std::vector<double>(m_NPhotons, std::numeric_limits<double>::max()));
declareProperty("Etacut", m_EtaRange = 2.50);
declareProperty("AllowSUSYDecay",m_AllowSUSYDecay = false);
// Backward compatibility aliases
declareProperty("Ptcut", m_Ptmin = 10000.);
declareProperty("Ptcut", m_Ptmin = std::vector<double>(m_NPhotons, 10000.));
}
StatusCode DirectPhotonFilter::filterInitialize() {
StatusCode DirectPhotonFilter::filterEvent() {
int NPhotons = 0;
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent* genEvt = (*itr);
ATH_MSG_DEBUG("----->>> Process : " << genEvt->signal_process_id());
for (HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr) {
if (((*pitr)->pdg_id() == 22)) {
if( ((*pitr)->momentum().perp() >= m_Ptmin) && ((*pitr)->momentum().perp() <= m_Ptmax) &&
fabs((*pitr)->momentum().pseudoRapidity()) <= m_EtaRange){
ATH_MSG_DEBUG("Generic photon found with status = " << (*pitr)->status() << " " << (*pitr)->barcode());
ATH_MSG_INFO("Initialising DirectPhoton filter with OrderPhotons="<<m_OrderPhotons);
HepMC::GenVertex* CandProdVertex = (*pitr)->production_vertex();
ATH_MSG_DEBUG("Candidate production vertex ID = " << CandProdVertex->id());
ATH_MSG_DEBUG("Candidate production vertex barcode = " << CandProdVertex->barcode());
if (m_Ptmin.size()>m_NPhotons || m_Ptmax.size()>m_NPhotons) {
ATH_MSG_ERROR("Too many Ptmin/max for given NPhotons");
return StatusCode::FAILURE;
}
HepMC::GenVertex::particle_iterator firstChild = (*pitr)->production_vertex()->particles_begin(HepMC::children);
HepMC::GenVertex::particle_iterator endChild = (*pitr)->production_vertex()->particles_end(HepMC::children);
HepMC::GenVertex::particle_iterator thisChild = firstChild;
for (; thisChild != endChild; ++thisChild) {
ATH_MSG_DEBUG("Looping on Production (children) vertex : " << (*thisChild)->pdg_id() << " " << (*thisChild)->barcode());
}
if (!m_OrderPhotons && (m_Ptmin.size()>1 || m_Ptmax.size()>1)) {
ATH_MSG_ERROR("Too many Ptmin/max requirements for OrderPhotons=false case.");
return StatusCode::FAILURE;
}
HepMC::GenVertex::particle_iterator firstChild1 = (*pitr)->production_vertex()->particles_begin(HepMC::parents);
HepMC::GenVertex::particle_iterator endChild1 = (*pitr)->production_vertex()->particles_end(HepMC::parents);
HepMC::GenVertex::particle_iterator thisChild1 = firstChild1;
for (; thisChild1 != endChild1; ++thisChild1) {
ATH_MSG_DEBUG("Looping on Production (parents) vertex : " << (*thisChild1)->pdg_id() << " " << (*thisChild1)->barcode());
}
}
}
// allow specifying only one pTmin/max to be applied to all (further) photons
// for backward compatibility
if (m_Ptmin.size()<m_NPhotons) {
size_t origsize = m_Ptmin.size();
double lastPt = m_Ptmin.back();
m_Ptmin.resize(m_NPhotons);
for (size_t i=origsize; i<m_NPhotons; ++i) m_Ptmin[i]=lastPt;
}
if (m_Ptmax.size()<m_NPhotons) {
size_t origsize = m_Ptmax.size();
double lastPt = m_Ptmax.back();
m_Ptmax.resize(m_NPhotons);
for (size_t i=origsize; i<m_NPhotons; ++i) m_Ptmax[i]=lastPt;
}
return StatusCode::SUCCESS;
}
// Check for a photon with desired kinematics
if ( ((*pitr)->pdg_id() == 22) ){
if( (*pitr)->status() == 1 &&
((*pitr)->momentum().perp() >= m_Ptmin) && ((*pitr)->momentum().perp() <= m_Ptmax) &&
fabs((*pitr)->momentum().pseudoRapidity()) <= m_EtaRange){
bool DirectPhotonFilterCmpByPt(HepMC::GenParticle* p1, HepMC::GenParticle* p2) {
return (p1->momentum().perp()<p2->momentum().perp());
}
// The following lines are for cross checking purpose when using different generators
HepMC::GenVertex* CandProdVertex = (*pitr)->production_vertex();
ATH_MSG_DEBUG("Candidate production vertex ID = " << CandProdVertex->id());
ATH_MSG_DEBUG("Candidate production vertex barcode = " << CandProdVertex->barcode());
StatusCode DirectPhotonFilter::filterEvent() {
std::vector<HepMC::GenParticle*> promptPhotonsInEta;
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent* genEvt = (*itr);
ATH_MSG_DEBUG("----->>> Process : " << genEvt->signal_process_id());
HepMC::GenVertex::particle_iterator firstChild = (*pitr)->production_vertex()->particles_begin(HepMC::children);
HepMC::GenVertex::particle_iterator endChild = (*pitr)->production_vertex()->particles_end(HepMC::children);
HepMC::GenVertex::particle_iterator thisChild = firstChild;
for(; thisChild != endChild; ++thisChild) {
ATH_MSG_DEBUG("Looping on Production (children) vertex : " << (*thisChild)->pdg_id() << " " << (*thisChild)->barcode());
// Find all prompt photons with within given eta range
for (HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr) {
if ((*pitr)->pdg_id() == 22 &&
(*pitr)->status() == 1 &&
fabs((*pitr)->momentum().pseudoRapidity()) <= m_EtaRange) {
// iterate over parent particles to exclude photons from hadron decays
HepMC::GenVertex* prodVtx = (*pitr)->production_vertex();
bool fromHadron(false);
for (auto parent = prodVtx->particles_begin(HepMC::parents);
parent != prodVtx->particles_end(HepMC::parents); ++parent) {
int pdgindex = abs((*parent)->pdg_id());
ATH_MSG_DEBUG("Looping on Production (parents) vertex : " << (*parent)->pdg_id() << " " << (*parent)->barcode());
if (pdgindex > 100) {
fromHadron = true;
if (m_AllowSUSYDecay && ( (pdgindex > 1000000 && pdgindex < 1000040) || (pdgindex > 2000000 && pdgindex < 2000016) ) ) fromHadron = false;
}
}
///////////////////////////////////////////////////////////////////////////////////////////////
//
// 1) once a status = 1 photon is found check where it comes from : loop on incoming particle
//
// Pythia : - a photon from HP comes out with status=3 and then turned into a status = 1 photon
// - coming from the status-3 photon.
// : - a photon from brem comes out directly with a status=1
// Herwig : - NOT CHECKED on HP !
// - a photon from brem comes out with different status but always turned into a
// status=1 photon at the end coming from a photon with a different status.
//
// So requiring a status=1 photon coming from q/CLHEP::g saves brem photons in Pythia. Requiring
// a status=1 photon coming from a photon should save HP photons in Pythia and Brem in
// Herwig
//
// 2) the second option is to ask for a status = 1 photon which doesn't come from a PDG>100
// particle. In this way we should veto 'background photon'
//
//////////////////////////////////////////////////////////////////////////////////////////////
HepMC::GenVertex::particle_iterator firstChild1 = (*pitr)->production_vertex()->particles_begin(HepMC::parents);
HepMC::GenVertex::particle_iterator endChild1 = (*pitr)->production_vertex()->particles_end(HepMC::parents);
HepMC::GenVertex::particle_iterator thisChild1 = firstChild1;
if (!fromHadron) promptPhotonsInEta.push_back((*pitr));
else ATH_MSG_DEBUG("non-prompt photon ignored");
}
}
}
bool fromHadron(false);
for (; thisChild1 != endChild1; ++thisChild1) {
int pdgindex = abs((*thisChild1)->pdg_id());
ATH_MSG_DEBUG("Looping on Production (parents) vertex : " << (*thisChild1)->pdg_id() << " " << (*thisChild1)->barcode());
if (pdgindex > 100) {
fromHadron = true;
if (m_AllowSUSYDecay && ( (pdgindex > 1000000 && pdgindex < 1000040) || (pdgindex > 2000000 && pdgindex < 2000016) ) ) fromHadron = false;
ATH_MSG_DEBUG("event kept");
}
}
if (!fromHadron) NPhotons++;
} // if( (*pitr)->status()==1 && ...
} // if( ((*pitr)->pdg_id() == 22) )
if (promptPhotonsInEta.size()<m_NPhotons) {
setFilterPassed(false);
}
else {
for (auto photon: promptPhotonsInEta) {
ATH_MSG_DEBUG("Found prompt photon with pt="<<photon->momentum().perp());
}
if (m_OrderPhotons) { // apply cuts to leading/subleading/... photon as specified
std::sort(promptPhotonsInEta.begin(), promptPhotonsInEta.end(), DirectPhotonFilterCmpByPt);
bool pass = true;
for (size_t i = 0; i < m_NPhotons; ++i) {
double pt = promptPhotonsInEta[i]->momentum().perp();
if (pt < m_Ptmin[i] || pt > m_Ptmax[i]) {
ATH_MSG_DEBUG(" rejected pt="<<pt);
pass = false;
}
}
if (pass) ATH_MSG_DEBUG("Passed!");
setFilterPassed(pass);
}
else { // just require NPhotons to pass m_Ptmin/max[0]
size_t NPhotons=0;
for (size_t i = 0; i < promptPhotonsInEta.size(); ++i) {
double pt = promptPhotonsInEta[i]->momentum().perp();
if (pt > m_Ptmin[0] && pt < m_Ptmax[0]) ++NPhotons;
}
if (NPhotons>=m_NPhotons) ATH_MSG_DEBUG("Passed!");
setFilterPassed(NPhotons>=m_NPhotons);
}
}
if (NPhotons >= m_NPhotons) return StatusCode::SUCCESS;
setFilterPassed(false);
return StatusCode::SUCCESS;
}
Supports Markdown
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