#include "SHERPA/Main/Filter.H" #include "SHERPA/Tools/Output_Base.H" #include "ATOOLS/Org/Scoped_Settings.H" using namespace SHERPA; using namespace ATOOLS; Filter::Filter() : m_on(false) { } Filter::~Filter() { for (std::map::iterator fit=m_filters.begin(); fit!=m_filters.end();fit++) delete fit->second; m_filters.clear(); } bool Filter::Init() { for (auto s : Settings::GetMainSettings()["FILTERS"].GetItems()) { Add(s); } return (m_on = m_filters.size()>0); } void Filter::Add(Scoped_Settings& s) { if (!s["Kf"].IsSetExplicitly()) THROW(fatal_error, "Every filter needs to specify \"Kf: \"."); FilterCriterion * crit = new FilterCriterion; crit->m_flav = Flavour(s["Kf"].SetDefault(0).Get()); crit->m_etamin = s["EtaMin"].SetDefault(-1e20).Get(); crit->m_etamax = s["EtaMax"].SetDefault( 1e20).Get(); crit->m_pTmin = s["PTMin"] .SetDefault( 0).Get(); crit->m_pTmax = s["PTMax"] .SetDefault( 1e20).Get(); crit->m_Nmin = s["NMin"] .SetDefault( 0).Get(); crit->m_Nmax = s["NMax"] .SetDefault( 1e20).Get(); m_filters[crit->m_flav] = crit; } void Filter::ShowSyntax(int mode) { if (!msg_LevelIsInfo() || mode == 0) return; msg_Out() << METHOD << "(): {\n\n"; msg_Out() << " FILTERS:\n"; msg_Out() << " - Kf: \n"; msg_Out() << " # optional:\n"; msg_Out() << " EtaMin: \n"; msg_Out() << " EtaMax: \n"; msg_Out() << " PTMin: \n"; msg_Out() << " PTMax: \n"; msg_Out() << " NMin: \n"; msg_Out() << " NMax: \n"; msg_Out() << " - ...\n"; msg_Out() << "\n}" << std::endl; } bool Filter::operator()(Blob_List * blobs) { if (!m_on) return true; Reset(); HarvestActiveParticles(blobs); FilterAccepted(); return Check(); } bool Filter::Check() { for (std::map::iterator flit=m_filters.begin(); flit!=m_filters.end();flit++) { //msg_Out()<first; std::map::iterator accit = m_accepted.find(flit->first); if (accit==m_accepted.end()) { //msg_Out()<<" --> not found at all.\n"; return false; } if (accit->second < flit->second->m_Nmin || accit->second > flit->second->m_Nmax) { //msg_Out()<<" --> wrong multiplicity.\n"; return false; } //msg_Out()<<" --> ok.\n"; } return true; } void Filter::FilterAccepted() { for (std::list::iterator pit=m_particles.begin(); pit!=m_particles.end();pit++) { Flavour flav = (*pit)->Flav(); std::map::iterator flit = m_filters.find(flav); if (flit==m_filters.end()) continue; Vec4D mom = (*pit)->Momentum(); double eta = mom.Eta(), pT = mom.PPerp(); FilterCriterion * crit = flit->second; if (eta>=crit->m_etamin && eta<=crit->m_etamax && pT>=crit->m_pTmin && pT<=crit->m_pTmax) { std::map::iterator accit = m_accepted.find(flav); if (accit==m_accepted.end()) m_accepted[flav]=1; else accit->second++; } } //for (std::map::iterator accit=m_accepted.begin(); // accit!=m_accepted.end();accit++) { // msg_Out()<<" --> found "<first<<": "<second<<"\n"; //} } void Filter::HarvestActiveParticles(Blob_List * blobs) { for (Blob_List::iterator bit=blobs->begin();bit!=blobs->end();bit++) { const Particle_Vector & particles = (*bit)->GetOutParticles(); for (Particle_Vector::const_iterator pit=particles.begin(); pit!=particles.end();pit++) { if (!(*pit)->DecayBlob()) m_particles.push_back((*pit)); } } } void Filter::Reset() { m_particles.clear(); m_accepted.clear(); }