Commit a3b92d1c authored by Petr Jacka's avatar Petr Jacka Committed by Nils Erik Krumnack
Browse files

Updating rc jets substructure

parent 61347f34
......@@ -751,12 +751,6 @@ namespace top {
//re-clustered jets
// -> need unordered map for systematics
bool m_makeRCJets; // making re-clustered jets
bool m_makeVarRCJets; // making VarRC jets
bool m_useRCJSS; // write RCJSS variables
bool m_useRCAdditionalJSS; // write RCJSS additional variables
bool m_useVarRCJSS; // write Variable-R RCJSS variables
bool m_useVarRCAdditionalJSS; // write Variable-R RCJSS additional variables
bool m_useElectronChargeIDSelection; // write ECID tool output variables
std::string m_RCJetContainer; // name for RC jets container in TStore
std::vector<std::string> m_VarRCJetRho;
......@@ -781,42 +775,9 @@ namespace top {
std::vector<std::vector<float> > m_rcjetsub_phi;
std::vector<std::vector<float> > m_rcjetsub_e;
std::vector<std::vector<float> > m_rcjetsub_mv2c10;
std::vector<float> m_rrcjet_pt;
std::vector<float> m_rrcjet_eta;
std::vector<float> m_rrcjet_phi;
std::vector<float> m_rrcjet_e;
std::vector<float> m_rcjet_tau32_clstr;
std::vector<float> m_rcjet_tau21_clstr;
std::vector<float> m_rcjet_tau3_clstr;
std::vector<float> m_rcjet_tau2_clstr;
std::vector<float> m_rcjet_tau1_clstr;
std::vector<float> m_rcjet_D2_clstr;
std::vector<float> m_rcjet_ECF1_clstr;
std::vector<float> m_rcjet_ECF2_clstr;
std::vector<float> m_rcjet_ECF3_clstr;
std::vector<float> m_rcjet_d12_clstr;
std::vector<float> m_rcjet_d23_clstr;
std::vector<float> m_rcjet_Qw_clstr;
std::vector<float> m_rcjet_nconstituent_clstr;
std::vector<float> m_rcjet_gECF332_clstr;
std::vector<float> m_rcjet_gECF461_clstr;
std::vector<float> m_rcjet_gECF322_clstr;
std::vector<float> m_rcjet_gECF331_clstr;
std::vector<float> m_rcjet_gECF422_clstr;
std::vector<float> m_rcjet_gECF441_clstr;
std::vector<float> m_rcjet_gECF212_clstr;
std::vector<float> m_rcjet_gECF321_clstr;
std::vector<float> m_rcjet_gECF311_clstr;
std::vector<float> m_rcjet_L1_clstr;
std::vector<float> m_rcjet_L2_clstr;
std::vector<float> m_rcjet_L3_clstr;
std::vector<float> m_rcjet_L4_clstr;
std::vector<float> m_rcjet_L5_clstr;
// maps containing rc jet substructure variables
std::map<std::string,std::vector<float>> m_rcjetJSSVariables;
std::map<std::string,std::map<std::string,std::vector<float>>> m_VarRCjetJSSVariables;
//met
......@@ -1426,8 +1387,6 @@ namespace top {
//re-clustered jets
// -> need unordered map for systematics
const bool& makeRCJets() const {return m_makeRCJets;} // making re-clustered jets
const bool& makeVarRCJets() const {return m_makeVarRCJets;} // making VarRC jets
const std::string& RCJetContainer() const {return m_RCJetContainer;} // name for RC jets container in TStore
const std::vector<std::string>& VarRCJetRho() const {return m_VarRCJetRho;}
const std::vector<std::string>& VarRCJetMassScale() const {return m_VarRCJetMassScale;}
......@@ -1449,34 +1408,10 @@ namespace top {
const std::vector<std::vector<float> >& rcjetsub_phi() const {return m_rcjetsub_phi;}
const std::vector<std::vector<float> >& rcjetsub_e() const {return m_rcjetsub_e;}
const std::vector<std::vector<float> >& rcjetsub_mv2c10() const {return m_rcjetsub_mv2c10;}
const std::vector<float>& rcjet_tau32_clstr() const {return m_rcjet_tau32_clstr;}
const std::vector<float>& rcjet_tau21_clstr() const {return m_rcjet_tau21_clstr;}
const std::vector<float>& rcjet_tau3_clstr() const {return m_rcjet_tau3_clstr;}
const std::vector<float>& rcjet_tau2_clstr() const {return m_rcjet_tau2_clstr;}
const std::vector<float>& rcjet_tau1_clstr() const {return m_rcjet_tau1_clstr;}
const std::vector<float>& rcjet_D2_clstr() const {return m_rcjet_D2_clstr;}
const std::vector<float>& rcjet_ECF1_clstr() const {return m_rcjet_ECF1_clstr;}
const std::vector<float>& rcjet_ECF2_clstr() const {return m_rcjet_ECF2_clstr;}
const std::vector<float>& rcjet_ECF3_clstr() const {return m_rcjet_ECF3_clstr;}
const std::vector<float>& rcjet_d12_clstr() const {return m_rcjet_d12_clstr;}
const std::vector<float>& rcjet_d23_clstr() const {return m_rcjet_d23_clstr;}
const std::vector<float>& rcjet_Qw_clstr() const {return m_rcjet_Qw_clstr;}
const std::vector<float>& rcjet_nconstituent_clstr() const {return m_rcjet_nconstituent_clstr;}
const std::vector<float>& rcjet_gECF332_clstr() const {return m_rcjet_gECF332_clstr;}
const std::vector<float>& rcjet_gECF461_clstr() const {return m_rcjet_gECF461_clstr;}
const std::vector<float>& rcjet_gECF322_clstr() const {return m_rcjet_gECF322_clstr;}
const std::vector<float>& rcjet_gECF331_clstr() const {return m_rcjet_gECF331_clstr;}
const std::vector<float>& rcjet_gECF422_clstr() const {return m_rcjet_gECF422_clstr;}
const std::vector<float>& rcjet_gECF441_clstr() const {return m_rcjet_gECF441_clstr;}
const std::vector<float>& rcjet_gECF212_clstr() const {return m_rcjet_gECF212_clstr;}
const std::vector<float>& rcjet_gECF321_clstr() const {return m_rcjet_gECF321_clstr;}
const std::vector<float>& rcjet_gECF311_clstr() const {return m_rcjet_gECF311_clstr;}
const std::vector<float>& rcjet_L1_clstr() const {return m_rcjet_L1_clstr;}
const std::vector<float>& rcjet_L2_clstr() const {return m_rcjet_L2_clstr;}
const std::vector<float>& rcjet_L3_clstr() const {return m_rcjet_L3_clstr;}
const std::vector<float>& rcjet_L4_clstr() const {return m_rcjet_L4_clstr;}
const std::vector<float>& rcjet_L5_clstr() const {return m_rcjet_L5_clstr;}
const std::map<std::string,std::vector<float>>& rcjetJSSVariables() const {return m_rcjetJSSVariables;}
const std::map<std::string,std::map<std::string,std::vector<float>>>& VarRCjetJSSVariables() const {return m_VarRCjetJSSVariables;}
//met
const float& met_met() const {return m_met_met;}
......
......@@ -247,8 +247,10 @@ namespace top {
registerParameter("RCJetRadius", "Reclustered Jet radius for object selection. Default 1.0", "1.0");
registerParameter("UseRCJetSubstructure", "Calculate Reclustered Jet Substructure Variables. Default False",
"False");
registerParameter("UseRCJetAdditionalSubstructure",
"Calculate Additional Reclustered Jet Substructure Variables. Default False", "False");
registerParameter("RCJetSubstructureVariables",
R"("List of reclustered jet substructure variables stored in the output separated by commas.
Available variables: Tau32_clstr,Tau21_clstr,Tau3_clstr,Tau2_clstr,Tau1_clstr,d12_clstr,d23_clstr,Qw_clstr,nconstituent_clstr,D2_clstr,ECF1_clstr,ECF2_clstr,ECF3_clstr,gECF332_clstr,gECF461_clstr,gECF322_clstr,gECF331_clstr,gECF422_clstr,gECF441_clstr,gECF212_clstr,gECF321_clstr,gECF311_clstr,L1_clstr,L2_clstr,L3_clstr,L4_clstr,L5_clstr,rrcjet_pt,rrcjet_eta,rrcjet_phi,rrcjet_e)",
" ");
registerParameter("UseRCJets", "Use Reclustered Jets. Default False.", "False");
......@@ -267,8 +269,11 @@ namespace top {
registerParameter("UseVarRCJets", "Use Reclustered Jets (Variable-R Jets). Default False.", "False");
registerParameter("UseVarRCJetSubstructure",
"Calculate Variable-R Reclustered Jet Substructure Variables. Default False", "False");
registerParameter("UseVarRCJetAdditionalSubstructure",
"Calculate Additional Variable-R Reclustered Jet Substructure Variables. Default False", "False");
registerParameter("VarRCJetSubstructureVariables",
R"(List of variable-R reclustered jet substructure variables stored in the output separated by commas.
Available variables: Tau32_clstr,Tau21_clstr,Tau3_clstr,Tau2_clstr,Tau1_clstr,d12_clstr,d23_clstr,Qw_clstr,nconstituent_clstr,D2_clstr,ECF1_clstr,ECF2_clstr,ECF3_clstr,gECF332_clstr,gECF461_clstr,gECF322_clstr,gECF331_clstr,gECF422_clstr,gECF441_clstr,gECF212_clstr,gECF321_clstr,gECF311_clstr,L1_clstr,L2_clstr,L3_clstr,L4_clstr,L5_clstr,rrcjet_pt,rrcjet_eta,rrcjet_phi,rrcjet_e)",
" ");
registerParameter("TauPt",
"Pt cut applied to both tight and loose taus (in MeV)."
......
......@@ -288,7 +288,6 @@ namespace top {
m_RCJetTrimcut(0.05),
m_RCJetRadius(1.0),
m_useRCJetSubstructure(false),
m_useRCJetAdditionalSubstructure(false),
m_VarRCJetPtcut(100000.),
m_VarRCJetEtacut(2.0),
......@@ -297,7 +296,6 @@ namespace top {
m_VarRCJetRho("2"),
m_VarRCJetMassScale("m_w,m_z,m_h,m_t"),
m_useVarRCJetSubstructure(false),
m_useVarRCJetAdditionalSubstructure(false),
m_trackQuality("SetMe"),
......@@ -1342,6 +1340,8 @@ namespace top {
this->largeRJetEtacut(std::stof(settings->value("LargeRJetEta")));
// now get all substructure variables from the config file.
std::string strSubstructure = settings->value("LargeRJetSubstructureVariables");
// Making vector of strings with "," used as separator
......@@ -1393,10 +1393,10 @@ namespace top {
if (settings->value("UseRCJetSubstructure") == "True" ||
settings->value("UseRCJetSubstructure") == "true") this->m_useRCJetSubstructure = true;
else this->m_useRCJetSubstructure = false;
if (settings->value("UseRCJetAdditionalSubstructure") == "True" ||
settings->value("UseRCJetAdditionalSubstructure") == "true") this->m_useRCJetAdditionalSubstructure = true;
else this->m_useRCJetAdditionalSubstructure = false;
if(this->m_useRCJetSubstructure) {
m_rcJetSubstructureVariables = readListOfSubstructureVariables(settings, "RCJetSubstructureVariables");
}
this->VarRCJetPtcut(std::stof(settings->value("VarRCJetPt")));
this->VarRCJetEtacut(std::stof(settings->value("VarRCJetEta")));
......@@ -1409,10 +1409,10 @@ namespace top {
if (settings->value("UseVarRCJetSubstructure") == "True" ||
settings->value("UseVarRCJetSubstructure") == "true") this->m_useVarRCJetSubstructure = true;
else this->m_useVarRCJetSubstructure = false;
if (settings->value("UseVarRCJetAdditionalSubstructure") == "True" ||
settings->value("UseVarRCJetAdditionalSubstructure") ==
"true") this->m_useVarRCJetAdditionalSubstructure = true;
else this->m_useVarRCJetAdditionalSubstructure = false;
if(this->m_useVarRCJetSubstructure) {
m_VarRCJetSubstructureVariables = readListOfSubstructureVariables(settings, "VarRCJetSubstructureVariables");
}
if (settings->value("StoreJetTruthLabels") == "False") {
this->jetStoreTruthLabels(false);
......@@ -3958,4 +3958,34 @@ namespace top {
os << "\n";
return os;
}
std::vector<std::string> TopConfig::readListOfSubstructureVariables(top::ConfigurationSettings* const& settings, const std::string& in) const {
const std::string& message = settings->stringData().at(in).m_human_explanation;
const std::string& value = settings->value(in);
const std::string substr = "Available variables:";
size_t found = message.find(substr);
if(found==std::string::npos) throw std::runtime_error("TopConfig: Expected 'Available variables:' string in " + in + " option description!");
const std::string strDefinedVariables = message.substr(found+substr.length());
std::vector<std::string> vecDefinedVariables;
tokenize(strDefinedVariables,vecDefinedVariables,",",true); // splitting string
for(std::string& defVar : vecDefinedVariables) {// removing spaces
defVar.erase(std::remove_if(defVar.begin(), defVar.end(), isspace), defVar.end());
}
std::vector<std::string> result;
tokenize(value,result,",",true); // splitting string
for(std::string& var : result) {
var.erase(std::remove_if(var.begin(), var.end(), isspace), var.end()); // removing empty spaces
if(std::find(std::begin(vecDefinedVariables),std::end(vecDefinedVariables),var) == std::end(vecDefinedVariables)) // Compare with list of defined variables
throw std::runtime_error("TopConfig: Unknown variable " + var + " found in " + in);
}
return result;
}
}
......@@ -1191,7 +1191,7 @@ namespace top {
inline virtual float trackEtacut() const {return m_trackEtacut;}
virtual std::vector<std::string> readListOfSubstructureVariables(top::ConfigurationSettings* const& settings, const std::string& in) const;
inline virtual float RCJetPtcut() const {return m_RCJetPtcut;}
inline virtual float RCJetEtacut() const {return m_RCJetEtacut;}
......@@ -1200,8 +1200,11 @@ namespace top {
inline virtual float RCJetTrimcut() const {return m_RCJetTrimcut;}
inline virtual float RCJetRadius() const {return m_RCJetRadius;}
inline virtual bool useRCJetSubstructure() const {return m_useRCJetSubstructure;}
inline virtual bool useRCJetAdditionalSubstructure() const {return m_useRCJetAdditionalSubstructure;}
inline virtual const std::vector<std::string>& rcJetSubstructureVariables() const {return m_rcJetSubstructureVariables;}
inline virtual void RCJetPtcut(const float pt) {
if (!m_configFixed) {
m_RCJetPtcut = pt;
......@@ -1244,12 +1247,13 @@ namespace top {
}
}
inline virtual void useRCJetAdditionalSubstructure(const bool use) {
inline virtual void rcJetSubstructureVariables(const std::vector<std::string>& use) {
if (!m_configFixed) {
m_useRCJetAdditionalSubstructure = use;
m_rcJetSubstructureVariables = use;
}
}
inline virtual float VarRCJetPtcut() const {return m_VarRCJetPtcut;}
inline virtual float VarRCJetEtacut() const {return m_VarRCJetEtacut;}
inline virtual float VarRCJetTrimcut() const {return m_VarRCJetTrimcut;}
......@@ -1257,7 +1261,7 @@ namespace top {
inline virtual const std::string& VarRCJetRho() const {return m_VarRCJetRho;}
inline virtual const std::string& VarRCJetMassScale() const {return m_VarRCJetMassScale;}
inline virtual bool useVarRCJetSubstructure() const {return m_useVarRCJetSubstructure;}
inline virtual bool useVarRCJetAdditionalSubstructure() const {return m_useVarRCJetAdditionalSubstructure;}
inline virtual const std::vector<std::string>& VarRCJetSubstructureVariables() const {return m_VarRCJetSubstructureVariables;}
inline virtual void VarRCJetPtcut(const float pt) {
if (!m_configFixed) {
......@@ -1301,9 +1305,9 @@ namespace top {
}
}
inline virtual void useVarRCJetAdditionalSubstructure(const bool use) {
inline virtual void VarRCJetSubstructureVariables(const std::vector<std::string>& use) {
if (!m_configFixed) {
m_useVarRCJetAdditionalSubstructure = use;
m_VarRCJetSubstructureVariables = use;
}
}
......@@ -2324,7 +2328,7 @@ namespace top {
float m_RCJetTrimcut;
float m_RCJetRadius;
bool m_useRCJetSubstructure;
bool m_useRCJetAdditionalSubstructure;
std::vector<std::string> m_rcJetSubstructureVariables;
// Jet configuration for variable large-R jets
float m_VarRCJetPtcut;
......@@ -2334,7 +2338,7 @@ namespace top {
std::string m_VarRCJetRho;
std::string m_VarRCJetMassScale;
bool m_useVarRCJetSubstructure;
bool m_useVarRCJetAdditionalSubstructure;
std::vector<std::string> m_VarRCJetSubstructureVariables;
std::string m_trackQuality; // track quality to be used in track selection
......
......@@ -48,7 +48,6 @@ RCJetMC15::RCJetMC15(const std::string& name) :
m_minradius(0.),
m_massscale(0.),
m_useJSS(false),
m_useAdditionalJSS(false),
m_egamma("EG_"),
m_jetsyst("JET_"),
m_muonsyst("MUON_"),
......@@ -108,7 +107,7 @@ StatusCode RCJetMC15::initialize() {
m_massscale = rho * m_scale * 1e-3; // e.g., 2*m_top; in [GeV]!
m_useJSS = m_config->useVarRCJetSubstructure();
m_useAdditionalJSS = m_config->useVarRCJetAdditionalSubstructure();
m_substructureVariables = m_config->VarRCJetSubstructureVariables();
} else {
m_ptcut = std::stof(configSettings->value("RCJetPt")); // for initialize [GeV] & passSelection
m_etamax = std::stof(configSettings->value("RCJetEta")); // for passSelection
......@@ -117,63 +116,58 @@ StatusCode RCJetMC15::initialize() {
m_minradius = -1.0;
m_massscale = -1.0;
m_useJSS = m_config->useRCJetSubstructure();
m_useAdditionalJSS = m_config->useRCJetAdditionalSubstructure();
m_substructureVariables = m_config->rcJetSubstructureVariables();
}
m_inputJetPtMin = std::stof(configSettings->value("RCInputJetPtMin"));
m_inputJetEtaMax = std::stof(configSettings->value("RCInputJetEtaMax"));
if (m_useJSS || m_useAdditionalJSS) {
if (m_useJSS) {
ATH_MSG_INFO("Calculating RCJet Substructure");
// Setup a bunch of FastJet stuff
//define the type of jets you will build (http://fastjet.fr/repo/doxygen-3.0.3/classfastjet_1_1JetDefinition.html)
m_jet_def_rebuild = std::make_shared<fastjet::JetDefinition>(fastjet::antikt_algorithm, 1.0, fastjet::E_scheme,
m_jet_def_rebuild = std::make_unique<fastjet::JetDefinition>(fastjet::antikt_algorithm, 1.0, fastjet::E_scheme,
fastjet::Best);
}
if (m_useJSS) {
//Substructure tool definitions
m_nSub1_beta1 = std::make_shared<fastjet::contrib::Nsubjettiness>(1,
m_nSub1_beta1 = std::make_unique<fastjet::contrib::Nsubjettiness>(1,
fastjet::contrib::OnePass_WTA_KT_Axes(),
fastjet::contrib::UnnormalizedMeasure(1.0));
m_nSub2_beta1 = std::make_shared<fastjet::contrib::Nsubjettiness>(2,
m_nSub2_beta1 = std::make_unique<fastjet::contrib::Nsubjettiness>(2,
fastjet::contrib::OnePass_WTA_KT_Axes(),
fastjet::contrib::UnnormalizedMeasure(1.0));
m_nSub3_beta1 = std::make_shared<fastjet::contrib::Nsubjettiness>(3,
m_nSub3_beta1 = std::make_unique<fastjet::contrib::Nsubjettiness>(3,
fastjet::contrib::OnePass_WTA_KT_Axes(),
fastjet::contrib::UnnormalizedMeasure(1.0));
m_split12 = std::make_shared<JetSubStructureUtils::KtSplittingScale>(1);
m_split23 = std::make_shared<JetSubStructureUtils::KtSplittingScale>(2);
m_split12 = std::make_unique<JetSubStructureUtils::KtSplittingScale>(1);
m_split23 = std::make_unique<JetSubStructureUtils::KtSplittingScale>(2);
m_qw = std::make_shared<JetSubStructureUtils::Qw>();
m_qw = std::make_unique<JetSubStructureUtils::Qw>();
m_ECF1 = std::make_shared<fastjet::contrib::EnergyCorrelator>(1, 1.0, fastjet::contrib::EnergyCorrelator::pt_R);
m_ECF2 = std::make_shared<fastjet::contrib::EnergyCorrelator>(2, 1.0, fastjet::contrib::EnergyCorrelator::pt_R);
m_ECF3 = std::make_shared<fastjet::contrib::EnergyCorrelator>(3, 1.0, fastjet::contrib::EnergyCorrelator::pt_R);
m_ECF1 = std::make_unique<fastjet::contrib::EnergyCorrelator>(1, 1.0, fastjet::contrib::EnergyCorrelator::pt_R);
m_ECF2 = std::make_unique<fastjet::contrib::EnergyCorrelator>(2, 1.0, fastjet::contrib::EnergyCorrelator::pt_R);
m_ECF3 = std::make_unique<fastjet::contrib::EnergyCorrelator>(3, 1.0, fastjet::contrib::EnergyCorrelator::pt_R);
}
if (m_useAdditionalJSS) {
m_gECF332 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(3, 3, 2,
m_gECF332 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(3, 3, 2,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF461 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(6, 4, 1,
m_gECF461 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(6, 4, 1,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF322 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(2, 3, 2,
m_gECF322 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(2, 3, 2,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF331 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(3, 3, 1,
m_gECF331 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(3, 3, 1,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF422 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(2, 4, 2,
m_gECF422 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(2, 4, 2,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF441 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(4, 4, 1,
m_gECF441 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(4, 4, 1,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF212 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(1, 2, 2,
m_gECF212 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(1, 2, 2,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF321 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(2, 3, 1,
m_gECF321 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(2, 3, 1,
JetSubStructureUtils::EnergyCorrelator::pt_R);
m_gECF311 = std::make_shared<JetSubStructureUtils::EnergyCorrelatorGeneralized>(1, 3, 1,
m_gECF311 = std::make_unique<JetSubStructureUtils::EnergyCorrelatorGeneralized>(1, 3, 1,
JetSubStructureUtils::EnergyCorrelator::pt_R);
}
......@@ -309,7 +303,7 @@ StatusCode RCJetMC15::execute(const top::Event& event) {
rcjet->auxdecor<bool>("PassedSelection") = passSelection(*rcjet);
}
if (m_useJSS || m_useAdditionalJSS) {
if (m_useJSS) {
static const SG::AuxElement::ConstAccessor<bool> passedSelection("PassedSelection");
for (auto rcjet : *myJets) {
......@@ -339,127 +333,110 @@ StatusCode RCJetMC15::execute(const top::Event& event) {
}
if (clusters.size() != 0) {
// Now rebuild the large jet from the small jet constituents aka the original clusters
fastjet::ClusterSequence clust_seq_rebuild = fastjet::ClusterSequence(clusters, *m_jet_def_rebuild);
std::vector<fastjet::PseudoJet> my_pjets = fastjet::sorted_by_pt(clust_seq_rebuild.inclusive_jets(0.0));
fastjet::PseudoJet correctedJet;
correctedJet = my_pjets[0];
//Sometimes fastjet splits the jet into two, so need to correct for that!!
if (my_pjets.size() > 1) correctedJet += my_pjets[1];
if (m_useJSS) {
// Now finally we can calculate some substructure!
double tau32 = -1, tau21 = -1;
double tau1 = m_nSub1_beta1->result(correctedJet);
double tau2 = m_nSub2_beta1->result(correctedJet);
double tau3 = m_nSub3_beta1->result(correctedJet);
if (std::abs(tau1) > 1e-8) tau21 = tau2 / tau1;
else tau21 = -999.0;
if (std::abs(tau2) > 1e-8) tau32 = tau3 / tau2;
else tau32 = -999.0;
double split12 = m_split12->result(correctedJet);
double split23 = m_split23->result(correctedJet);
double qw = m_qw->result(correctedJet);
double D2 = -1;
double vECF1 = m_ECF1->result(correctedJet);
double vECF2 = m_ECF2->result(correctedJet);
double vECF3 = m_ECF3->result(correctedJet);
if (std::abs(vECF2) > 1e-8) D2 = vECF3 * vECF1* vECF1* vECF1 / (vECF2 * vECF2 * vECF2);
else D2 = -999.0;
// now attach the results to the original jet
rcjet->auxdecor<float>("Tau32_clstr") = tau32;
rcjet->auxdecor<float>("Tau21_clstr") = tau21;
// lets also write out the components so we can play with them later
rcjet->auxdecor<float>("Tau3_clstr") = tau3;
rcjet->auxdecor<float>("Tau2_clstr") = tau2;
rcjet->auxdecor<float>("Tau1_clstr") = tau1;
rcjet->auxdecor<float>("d12_clstr") = split12;
rcjet->auxdecor<float>("d23_clstr") = split23;
rcjet->auxdecor<float>("Qw_clstr") = qw;
rcjet->auxdecor<float>("nconstituent_clstr") = clusters.size();
rcjet->auxdecor<float>("ECF1_clstr") = vECF1;
rcjet->auxdecor<float>("ECF2_clstr") = vECF2;
rcjet->auxdecor<float>("ECF3_clstr") = vECF3;
rcjet->auxdecor<float>("D2_clstr") = D2;
} // end of if useJSS
if (m_useAdditionalJSS) {
// MlB's t/H discriminators
// E = (a*n) / (b*m)
// for an ECFG_X_Y_Z, a=Y, n=Z -> dimenionless variable
double gECF332 = m_gECF332->result(correctedJet);
double gECF461 = m_gECF461->result(correctedJet);
double gECF322 = m_gECF322->result(correctedJet);
double gECF331 = m_gECF331->result(correctedJet);
double gECF422 = m_gECF422->result(correctedJet);
double gECF441 = m_gECF441->result(correctedJet);
double gECF212 = m_gECF212->result(correctedJet);
double gECF321 = m_gECF321->result(correctedJet);
double gECF311 = m_gECF311->result(correctedJet);
double L1 = -999.0, L2 = -999.0, L3 = -999.0, L4 = -999.0, L5 = -999.0;
if (std::abs(gECF212) > 1e-12) {
L1 = gECF321 / gECF212;
L2 = gECF331 / sqrt(gECF212*gECF212*gECF212);
}
if (std::abs(gECF331) > 1e-12) {
L3 = gECF311 / pow(gECF331,1./3.);
L4 = gECF322 / pow(gECF331,4./3.);
}
if (std::abs(gECF441) > 1e-12) {
L5 = gECF422/gECF441;
}
rcjet->auxdecor<float>("gECF332_clstr") = gECF332;
rcjet->auxdecor<float>("gECF461_clstr") = gECF461;
rcjet->auxdecor<float>("gECF322_clstr") = gECF322;
rcjet->auxdecor<float>("gECF331_clstr") = gECF331;
rcjet->auxdecor<float>("gECF422_clstr") = gECF422;
rcjet->auxdecor<float>("gECF441_clstr") = gECF441;
rcjet->auxdecor<float>("gECF212_clstr") = gECF212;
rcjet->auxdecor<float>("gECF321_clstr") = gECF321;
rcjet->auxdecor<float>("gECF311_clstr") = gECF311;
rcjet->auxdecor<float>("L1_clstr") = L1;
rcjet->auxdecor<float>("L2_clstr") = L2;
rcjet->auxdecor<float>("L3_clstr") = L3;
rcjet->auxdecor<float>("L4_clstr") = L4;
rcjet->auxdecor<float>("L5_clstr") = L5;
// lets also store the rebuilt jet incase we need it later
rcjet->auxdecor<float>("RRCJet_pt") = correctedJet.pt();
rcjet->auxdecor<float>("RRCJet_eta") = correctedJet.eta();
rcjet->auxdecor<float>("RRCJet_phi") = correctedJet.phi();
rcjet->auxdecor<float>("RRCJet_e") = correctedJet.e();
}// end of if useAdditional JSS
// Now rebuild the large jet from the small jet constituents aka the original clusters
fastjet::ClusterSequence clust_seq_rebuild = fastjet::ClusterSequence(clusters, *m_jet_def_rebuild);
std::vector<fastjet::PseudoJet> my_pjets = fastjet::sorted_by_pt(clust_seq_rebuild.inclusive_jets(0.0));
fastjet::PseudoJet correctedJet;
correctedJet = my_pjets[0];
//Sometimes fastjet splits the jet into two, so need to correct for that!!
if (my_pjets.size() > 1) correctedJet += my_pjets[1];
// Now finally we can calculate some substructure!
auto it1 = std::begin(m_substructureVariables);
auto it2 = std::end(m_substructureVariables);
float tau1(-999.),tau2(-999.),tau3(-999.);
if(std::find(it1,it2,"Tau1_clstr")!=it2) rcjet->auxdecor<float>("Tau1_clstr") = tau1 = m_nSub1_beta1->result(correctedJet);
if(std::find(it1,it2,"Tau2_clstr")!=it2) rcjet->auxdecor<float>("Tau2_clstr") = tau2 = m_nSub2_beta1->result(correctedJet);
if(std::find(it1,it2,"Tau3_clstr")!=it2) rcjet->auxdecor<float>("Tau3_clstr") = tau3 = m_nSub3_beta1->result(correctedJet);
if(std::find(it1,it2,"Tau21_clstr")!=it2) {
if(tau1<0.) tau1 = m_nSub1_beta1->result(correctedJet);
if(tau2<0.) tau2 = m_nSub2_beta1->result(correctedJet);
rcjet->auxdecor<float>("Tau21_clstr") = std::abs(tau1) > 1e-8 ? tau2/tau1 : -999.;
}
if(std::find(it1,it2,"Tau32_clstr")!=it2) {
if(tau3<0.) tau3 = m_nSub3_beta1->result(correctedJet);
if(tau2<0.) tau2 = m_nSub2_beta1->result(correctedJet);