diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/Root/EventSaverFlatNtuple.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/Root/EventSaverFlatNtuple.cxx
index c790fa5fdb02db63f954ed29ced2fa369e9f33f5..bc63a435057a17ed0886fb48714dc72d7c8a7085 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/Root/EventSaverFlatNtuple.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/Root/EventSaverFlatNtuple.cxx
@@ -17,7 +17,6 @@
#include <boost/algorithm/string.hpp>
#include "TopParticleLevel/ParticleLevelEvent.h"
-#include "TopParticleLevel/ParticleLevelRCJetObjectLoader.h"
#include "TopFakes/TopFakesMMWeightCalculator.h"
@@ -201,38 +200,24 @@ namespace top {
// fixed-R re-clustering (RC)
- if (config->useRCJets() == true){
+ if (config->useRCJets()){
m_makeRCJets = true;
-
- if (config->useRCJetSubstructure() == true)
- m_useRCJSS = true;
+ m_useRCJSS = config->useRCJetSubstructure();
- m_rc = std::unique_ptr<RCJetMC15> ( new RCJetMC15( "RCJetMC15" ) );
- top::check(m_rc->setProperty( "config" , config ) , "Failed to set config property of RCJetMC15");
- top::check(m_rc->initialize(),"Failed to initialize RCJetMC15");
}
// variable-R re-clustering (VarRC)
- if (config->useVarRCJets() == true){
+ if (config->useVarRCJets()){
m_makeVarRCJets = true;
m_VarRCjetBranches.clear(); // clear map of branches just in case
m_VarRCjetsubBranches.clear();
+ m_VarRCjetBranchesParticle.clear();
+ m_VarRCjetsubBranchesParticle.clear();
+
boost::split(m_VarRCJetRho, config->VarRCJetRho(), boost::is_any_of(","));
boost::split(m_VarRCJetMassScale, config->VarRCJetMassScale(), boost::is_any_of(","));
- for (auto& rho : m_VarRCJetRho){
- for (auto& mass_scale : m_VarRCJetMassScale){
- std::replace( rho.begin(), rho.end(), '.', '_');
- std::string name = rho+mass_scale;
- m_VarRC[name] = std::unique_ptr<RCJetMC15> ( new RCJetMC15( "VarRCJetMC15_"+name ) );
- top::check(m_VarRC[name]->setProperty( "config" , config ) , "Failed to set config property of VarRCJetMC15");
- top::check(m_VarRC[name]->setProperty( "VarRCjets", true ) , "Failed to set VarRCjets property of VarRCJetMC15");
- top::check(m_VarRC[name]->setProperty( "VarRCjets_rho", rho ) , "Failed to set VarRCjets rho property of VarRCJetMC15");
- top::check(m_VarRC[name]->setProperty( "VarRCjets_mass_scale", mass_scale ) , "Failed to set VarRCjets mass scale property of VarRCJetMC15");
- top::check(m_VarRC[name]->initialize(),"Failed to initialize VarRCJetMC15");
- } // end loop over mass scale parameters (e.g., top mass, w mass, etc.)
- } // end loop over mass scale multiplies (e.g., 1.,2.,etc.)
} // end make VarRC jets
@@ -1158,10 +1143,7 @@ namespace top {
// RC branches
if (m_makeRCJets){
- // first initialise the ParticleLevelRCJetObjectLoader
- m_rcjet_particle = std::unique_ptr<ParticleLevelRCJetObjectLoader> ( new ParticleLevelRCJetObjectLoader( m_config ) );
- top::check(m_rcjet_particle->initialize(),"Failed to initialize ParticleLevelRCJetObjectLoader");
- // then create the branches
+ // create the branches
m_particleLevelTreeManager->makeOutputVariable(m_rcjet_pt, "rcjet_pt");
m_particleLevelTreeManager->makeOutputVariable(m_rcjet_eta, "rcjet_eta");
m_particleLevelTreeManager->makeOutputVariable(m_rcjet_phi, "rcjet_phi");
@@ -1217,6 +1199,27 @@ namespace top {
}
+ if (m_makeVarRCJets){
+ std::string VarRC = "vrcjet";
+
+ for (auto& rho : m_VarRCJetRho){
+ for (auto& mass_scale : m_VarRCJetMassScale){
+ std::replace( rho.begin(), rho.end(), '.', '_');
+ std::string name = rho+mass_scale;
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetBranchesParticle[VarRC+"_"+name+"_pt"], VarRC+"_"+name+"_pt");
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetBranchesParticle[VarRC+"_"+name+"_eta"],VarRC+"_"+name+"_eta");
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetBranchesParticle[VarRC+"_"+name+"_phi"],VarRC+"_"+name+"_phi");
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetBranchesParticle[VarRC+"_"+name+"_e"], VarRC+"_"+name+"_e");
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetBranchesParticle[VarRC+"_"+name+"_d12"],VarRC+"_"+name+"_d12"); // requires >= 2 subjets
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetBranchesParticle[VarRC+"_"+name+"_d23"],VarRC+"_"+name+"_d23"); // requires >= 3 subjets
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_pt"], VarRC+"sub_"+name+"_pt"); // vector of vectors for subjet info
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_eta"], VarRC+"sub_"+name+"_eta");
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_phi"], VarRC+"sub_"+name+"_phi");
+ m_particleLevelTreeManager->makeOutputVariable(m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_e"], VarRC+"sub_"+name+"_e");
+ } // end loop over mass parameters
+ } // end loop over multipliers for mass scale
+ } // end if VarRC jets
+
//met
if ( m_config->useTruthMET() ){
m_particleLevelTreeManager->makeOutputVariable(m_met_met, "met_met");
@@ -2125,16 +2128,6 @@ namespace top {
}
if (m_makeRCJets){
- // Execute the re-clustering code
- // - make jet container of small-r jets in the event, put it in TStore, do re-clustering
- top::check(m_rc->execute(event),"Failed to execute RCJetMC15 container");
-
- // Get the name of the container of re-clustered jets in TStore
- m_RCJetContainer = m_rc->rcjetContainerName(event.m_hashValue,event.m_isLoose);
-
- // -- Retrieve the re-clustered jets from TStore & save good re-clustered jets -- //
- const xAOD::JetContainer* rc_jets(nullptr);
- top::check(evtStore()->retrieve(rc_jets,m_RCJetContainer),"Failed to retrieve RC JetContainer");
// re-clustered jet substructure
static SG::AuxElement::ConstAccessor<float> RCSplit12("Split12");
@@ -2177,7 +2170,7 @@ namespace top {
static SG::AuxElement::ConstAccessor<float> RRCJet_e("RRCJet_e");
// Initialize the vectors to be saved as branches
- unsigned int sizeOfRCjets(rc_jets->size());
+ unsigned int sizeOfRCjets(event.m_RCJets.size());
m_rcjet_pt.clear();
m_rcjet_eta.clear();
@@ -2275,11 +2268,10 @@ namespace top {
m_rcjet_L5_clstr.resize(sizeOfRCjets,-999.);
}
unsigned int i = 0;
- for (xAOD::JetContainer::const_iterator jet_itr = rc_jets->begin(); jet_itr != rc_jets->end(); ++jet_itr) {
+ for (auto jet_itr = event.m_RCJets.begin(); jet_itr != event.m_RCJets.end(); ++jet_itr) {
+
const xAOD::Jet* rc_jet = *jet_itr;
- if (!m_rc->passSelection(*rc_jet))
- continue;
m_rcjet_pt[i] = rc_jet->pt();
m_rcjet_eta[i] = rc_jet->eta();
@@ -2357,83 +2349,28 @@ namespace top {
++i;
} // end for-loop over re-clustered jets
- m_rcjet_pt.resize(i);
- m_rcjet_eta.resize(i);
- m_rcjet_phi.resize(i);
- m_rcjet_e.resize(i);
- m_rcjet_d12.resize(i);
- m_rcjet_d23.resize(i);
- m_rcjetsub_pt.resize(i, std::vector<float>());
- m_rcjetsub_eta.resize(i, std::vector<float>());
- m_rcjetsub_phi.resize(i, std::vector<float>());
- m_rcjetsub_e.resize(i, std::vector<float>());
- m_rcjetsub_mv2c10.resize(i, std::vector<float>());
-
- if (m_useRCJSS){
- m_rrcjet_pt.resize(i);
- m_rrcjet_eta.resize(i);
- m_rrcjet_phi.resize(i);
- m_rrcjet_e.resize(i);
-
- m_rcjet_tau21_clstr.resize(i);
- m_rcjet_tau32_clstr.resize(i);
- m_rcjet_tau1_clstr.resize(i);
- m_rcjet_tau2_clstr.resize(i);
- m_rcjet_tau3_clstr.resize(i);
-
- m_rcjet_D2_clstr.resize(i);
- m_rcjet_ECF1_clstr.resize(i);
- m_rcjet_ECF2_clstr.resize(i);
- m_rcjet_ECF3_clstr.resize(i);
-
- m_rcjet_d12_clstr.resize(i);
- m_rcjet_d23_clstr.resize(i);
- m_rcjet_Qw_clstr.resize(i);
-
- m_rcjet_gECF332_clstr.resize(i);
- m_rcjet_gECF461_clstr.resize(i);
- m_rcjet_gECF322_clstr.resize(i);
- m_rcjet_gECF331_clstr.resize(i);
- m_rcjet_gECF422_clstr.resize(i);
- m_rcjet_gECF441_clstr.resize(i);
- m_rcjet_gECF212_clstr.resize(i);
- m_rcjet_gECF321_clstr.resize(i);
- m_rcjet_gECF311_clstr.resize(i);
- m_rcjet_L1_clstr.resize(i);
- m_rcjet_L2_clstr.resize(i);
- m_rcjet_L3_clstr.resize(i);
- m_rcjet_L4_clstr.resize(i);
- m_rcjet_L5_clstr.resize(i);
- }
+
} // end if make rcjets
// end re-clustered jets
/**********************************/
// VarRC jets
if (m_makeVarRCJets){
- // Execute the re-clustering code
- // - make jet container, put it in TStore, do re-clustering
std::string VarRC = "vrcjet";
+ std::unordered_map< std::string,std::shared_ptr<xAOD::JetContainer> > VarRCJets=event.m_VarRCJets;
for (auto& rho : m_VarRCJetRho){
for (auto& mass_scale : m_VarRCJetMassScale){
std::replace( rho.begin(), rho.end(), '.', '_');
std::string name = rho+mass_scale;
- top::check(m_VarRC[name]->execute(event),"Failed to execute RCJetMC15 container");
-
- // Get the name of the container of re-clustered jets in TStore
- m_RCJetContainer = m_VarRC[name]->rcjetContainerName(event.m_hashValue,event.m_isLoose);
-
- // -- Retrieve the re-clustered jets from TStore & save good re-clustered jets -- //
- const xAOD::JetContainer* vrc_jets(nullptr);
- top::check(evtStore()->retrieve(vrc_jets,m_RCJetContainer),"Failed to retrieve RC JetContainer");
-
// re-clustered jet substructure
static SG::AuxElement::ConstAccessor<float> VarRCSplit12("Split12");
static SG::AuxElement::ConstAccessor<float> VarRCSplit23("Split23");
// Initialize the vectors to be saved as branches
- unsigned int sizeOfRCjets(vrc_jets->size());
+
+ xAOD::JetContainer* vrc_jets = VarRCJets[name].get();
+ unsigned int sizeOfRCjets = vrc_jets->size();
m_VarRCjetBranches[VarRC+"_"+name+"_pt"].resize(sizeOfRCjets,-999.);
m_VarRCjetBranches[VarRC+"_"+name+"_eta"].resize(sizeOfRCjets,-999.);
m_VarRCjetBranches[VarRC+"_"+name+"_phi"].resize(sizeOfRCjets,-999.);
@@ -2447,12 +2384,10 @@ namespace top {
m_VarRCjetsubBranches[VarRC+"_"+name+"_sub_mv2c10"].resize(sizeOfRCjets, std::vector<float>());
unsigned int i = 0;
- for (xAOD::JetContainer::const_iterator jet_itr = vrc_jets->begin(); jet_itr != vrc_jets->end(); ++jet_itr) {
- const xAOD::Jet* rc_jet = *jet_itr;
-
- if (!m_VarRC[name]->passSelection(*rc_jet))
- continue;
-
+
+ for (auto jet_ptr : *vrc_jets) {
+ const xAOD::Jet* rc_jet = jet_ptr;
+
m_VarRCjetBranches[VarRC+"_"+name+"_pt"][i] = rc_jet->pt();
m_VarRCjetBranches[VarRC+"_"+name+"_eta"][i] = rc_jet->eta();
m_VarRCjetBranches[VarRC+"_"+name+"_phi"][i] = rc_jet->phi();
@@ -2492,17 +2427,6 @@ namespace top {
} // end for-loop over re-clustered jets
- m_VarRCjetBranches[VarRC+"_"+name+"_pt"].resize(i);
- m_VarRCjetBranches[VarRC+"_"+name+"_eta"].resize(i);
- m_VarRCjetBranches[VarRC+"_"+name+"_phi"].resize(i);
- m_VarRCjetBranches[VarRC+"_"+name+"_e"].resize(i);
- m_VarRCjetBranches[VarRC+"_"+name+"_d12"].resize(i);
- m_VarRCjetBranches[VarRC+"_"+name+"_d23"].resize(i);
- m_VarRCjetsubBranches[VarRC+"_"+name+"_sub_pt"].resize(i, std::vector<float>());
- m_VarRCjetsubBranches[VarRC+"_"+name+"_sub_eta"].resize(i, std::vector<float>());
- m_VarRCjetsubBranches[VarRC+"_"+name+"_sub_phi"].resize(i, std::vector<float>());
- m_VarRCjetsubBranches[VarRC+"_"+name+"_sub_e"].resize(i, std::vector<float>());
- m_VarRCjetsubBranches[VarRC+"_"+name+"_sub_mv2c10"].resize(i, std::vector<float>());
} // end loop over mass parameters
} // end loop over multipliers for mass scale
} // end if make VarRC jets
@@ -3171,17 +3095,7 @@ namespace top {
}
if(m_makeRCJets){
- // Execute the re-clustering code
- // - make jet container of small-r jets in the event, put it in TStore, do re-clustering
- top::check(m_rcjet_particle->execute(plEvent),"Failed to execute ParticleLevelRCJetObjectLoader container");
-
- // Get the name of the container of re-clustered jets
- m_RCJetContainerParticle = m_rcjet_particle->rcjetContainerName();
-
- // -- Retrieve the re-clustered jets from TStore & save good re-clustered jets -- //
- const xAOD::JetContainer* rc_jets_particle(nullptr);
- top::check(evtStore()->retrieve(rc_jets_particle,m_RCJetContainerParticle),"Failed to retrieve particle RC JetContainer");
-
+
// re-clustered jet substructure
static SG::AuxElement::ConstAccessor<float> RCSplit12("Split12");
static SG::AuxElement::ConstAccessor<float> RCSplit23("Split23");
@@ -3221,7 +3135,7 @@ namespace top {
// Initialize the vectors to be saved as branches
- unsigned int sizeOfRCjets(rc_jets_particle->size());
+ unsigned int sizeOfRCjets(plEvent.m_RCJets.size());
m_rcjet_pt.clear();
m_rcjet_eta.clear();
@@ -3317,15 +3231,8 @@ namespace top {
}
unsigned int i = 0;
- std::vector<int> rc_particle_selected_jets;
- rc_particle_selected_jets.clear();
-
- for (xAOD::JetContainer::const_iterator jet_itr = rc_jets_particle->begin(); jet_itr != rc_jets_particle->end(); ++jet_itr) {
- const xAOD::Jet* rc_jet = *jet_itr;
- if (!m_rcjet_particle->passSelection(*rc_jet))
- continue;
-
- rc_particle_selected_jets.push_back(rc_jet->index());
+ for( auto rc_jet : plEvent.m_RCJets){
+
m_rcjet_pt[i] = rc_jet->pt();
m_rcjet_eta[i] = rc_jet->eta();
@@ -3398,20 +3305,6 @@ namespace top {
} // end for-loop over subjets
++i;
} // end for-loop over re-clustered jets
- // we resized earlier to the size of rc_jets_particle container, but then only stored a sub-set
- // so have to resize again to shrink to correct size incase some jets were not stored
- m_rcjet_pt.resize(i,-999);
- m_rcjet_eta.resize(i,-999.);
- m_rcjet_phi.resize(i,-999.);
- m_rcjet_e.resize(i,-999.);
- m_rcjet_d12.resize(i,-999.);
- m_rcjet_d23.resize(i,-999.);
- m_rcjetsub_pt.resize(i, std::vector<float>());
- m_rcjetsub_eta.resize(i, std::vector<float>());
- m_rcjetsub_phi.resize(i, std::vector<float>());
- m_rcjetsub_e.resize(i, std::vector<float>());
- m_rcjetsub_Ghosts_BHadron_Final_Count.resize(i, std::vector<int>());
- m_rcjetsub_Ghosts_CHadron_Final_Count.resize(i, std::vector<int>());
if (m_useRCJSS){
m_rrcjet_pt.resize(i);
@@ -3452,6 +3345,72 @@ namespace top {
}
+ // VarRC jets
+ if (m_makeVarRCJets){
+ std::string VarRC = "vrcjet";
+ std::unordered_map< std::string,std::shared_ptr<xAOD::JetContainer> > VarRCJets=plEvent.m_VarRCJets;
+ for (auto& rho : m_VarRCJetRho){
+ for (auto& mass_scale : m_VarRCJetMassScale){
+ std::replace( rho.begin(), rho.end(), '.', '_');
+ std::string name = rho+mass_scale;
+
+ // re-clustered jet substructure
+ static SG::AuxElement::ConstAccessor<float> VarRCSplit12("Split12");
+ static SG::AuxElement::ConstAccessor<float> VarRCSplit23("Split23");
+
+ // Initialize the vectors to be saved as branches
+
+ xAOD::JetContainer* vrc_jets = VarRCJets[name].get();
+ unsigned int sizeOfRCjets = vrc_jets->size();
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_pt"].resize(sizeOfRCjets,-999.);
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_eta"].resize(sizeOfRCjets,-999.);
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_phi"].resize(sizeOfRCjets,-999.);
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_e"].resize(sizeOfRCjets,-999.);
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_d12"].resize(sizeOfRCjets,-999.);
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_d23"].resize(sizeOfRCjets,-999.);
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_pt"].resize(sizeOfRCjets, std::vector<float>());
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_eta"].resize(sizeOfRCjets, std::vector<float>());
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_phi"].resize(sizeOfRCjets, std::vector<float>());
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_e"].resize(sizeOfRCjets, std::vector<float>());
+
+ unsigned int i = 0;
+
+ for (auto jet_ptr : *vrc_jets) {
+ const xAOD::Jet* rc_jet = jet_ptr;
+
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_pt"][i] = rc_jet->pt();
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_eta"][i] = rc_jet->eta();
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_phi"][i] = rc_jet->phi();
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_e"][i] = rc_jet->e();
+
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_d12"][i] = (VarRCSplit12.isAvailable(*rc_jet)) ? VarRCSplit12(*rc_jet) : -999.;
+ m_VarRCjetBranchesParticle[VarRC+"_"+name+"_d23"][i] = (VarRCSplit23.isAvailable(*rc_jet)) ? VarRCSplit23(*rc_jet) : -999.;
+
+ // loop over subjets
+ const xAOD::Jet* subjet(nullptr);
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_pt"][i].clear(); // clear the vector size (otherwise it grows out of control!)
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_eta"][i].clear();
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_phi"][i].clear();
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_e"][i].clear();
+ for(auto rc_jet_subjet : rc_jet->getConstituents()){
+ subjet = static_cast<const xAOD::Jet*>(rc_jet_subjet->rawConstituent());
+
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_pt"][i].push_back(subjet->pt());
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_eta"][i].push_back(subjet->eta());
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_phi"][i].push_back(subjet->phi());
+ m_VarRCjetsubBranchesParticle[VarRC+"_"+name+"_sub_e"][i].push_back(subjet->e());
+ } // end for-loop over subjets
+ ++i;
+
+ } // end for-loop over re-clustered jets
+
+ } // end loop over mass parameters
+ } // end loop over multipliers for mass scale
+ } // end if make VarRC jets
+ // end VarRC jets
+
+
+
//met
if ( m_config->useTruthMET() ){
m_met_met = plEvent.m_met->met();
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/TopAnalysis/EventSaverFlatNtuple.h b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/TopAnalysis/EventSaverFlatNtuple.h
index f0290e1688afbc0ed91c76ae3bcc0a00d2668626..7ca7186296649c63fe7ff0a283726ca72358f08d 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/TopAnalysis/EventSaverFlatNtuple.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/TopAnalysis/EventSaverFlatNtuple.h
@@ -8,7 +8,6 @@
#include "TopAnalysis/EventSaverBase.h"
#include "TopCorrections/ScaleFactorRetriever.h"
#include "TopEventSelectionTools/TreeManager.h"
-#include "TopObjectSelectionTools/RCJetMC15.h"
// Framework include(s):
#include "AsgTools/AsgTool.h"
@@ -16,8 +15,6 @@
#include "AthContainers/DataVector.h"
#include <unordered_map>
-// fwd declare particle-level class
-class ParticleLevelRCJetObjectLoader;
namespace top {
@@ -583,14 +580,13 @@ private:
std::string m_RCJetContainer; // name for RC jets container in TStore
std::vector<std::string> m_VarRCJetRho;
std::vector<std::string> m_VarRCJetMassScale;
- std::unique_ptr<RCJetMC15> m_rc;
- std::map<std::string,std::unique_ptr<RCJetMC15> > m_VarRC;
std::string m_egamma; // egamma systematic naming scheme
std::string m_muonsyst; // muon systematic naming scheme
std::string m_jetsyst; // jet systematic naming scheme
- std::unordered_map<std::size_t, JetReclusteringTool*> m_jetReclusteringTool;
std::map<std::string,std::vector<float>> m_VarRCjetBranches;
std::map<std::string,std::vector<std::vector<float>>> m_VarRCjetsubBranches;
+ std::map<std::string,std::vector<float>> m_VarRCjetBranchesParticle;
+ std::map<std::string,std::vector<std::vector<float>>> m_VarRCjetsubBranchesParticle;
std::vector<int> m_rcjet_nsub;
std::vector<float> m_rcjet_pt;
std::vector<float> m_rcjet_eta;
@@ -844,11 +840,6 @@ private:
std::vector<std::vector<int>> m_rcjetsub_Ghosts_BHadron_Final_Count;
std::vector<std::vector<int>> m_rcjetsub_Ghosts_CHadron_Final_Count;
- // for rc jets at particle level
- std::unique_ptr<ParticleLevelRCJetObjectLoader> m_rcjet_particle;
- std::string m_RCJetContainerParticle; // name for RC jets container
-
-
// Truth tree inserted variables
// This can be expanded as required
// This is just a first pass at doing this sort of thing
@@ -1161,12 +1152,9 @@ protected:
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;}
- const std::unique_ptr<RCJetMC15>& rc() const { return m_rc;}
- const std::map<std::string,std::unique_ptr<RCJetMC15> >& VarRC() const { return m_VarRC;}
const std::string& egamma() const { return m_egamma;} // egamma systematic naming scheme
const std::string& muonsyst() const { return m_muonsyst;} // muon systematic naming scheme
const std::string& jetsyst() const { return m_jetsyst;} // jet systematic naming scheme
- const std::unordered_map<std::size_t, JetReclusteringTool*>& jetReclusteringTool() const { return m_jetReclusteringTool;}
const std::map<std::string,std::vector<float>>& VarRCjetBranches() const { return m_VarRCjetBranches;}
const std::map<std::string,std::vector<std::vector<float>>>& VarRCjetsubBranches() const { return m_VarRCjetsubBranches;}
const std::vector<int>& rcjet_nsub() const { return m_rcjet_nsub;}
@@ -1378,10 +1366,6 @@ protected:
const std::vector<std::vector<int>>& rcjetsub_Ghosts_BHadron_Final_Count() const { return m_rcjetsub_Ghosts_BHadron_Final_Count;}
const std::vector<std::vector<int>>& rcjetsub_Ghosts_CHadron_Final_Count() const { return m_rcjetsub_Ghosts_CHadron_Final_Count;}
- // for rc jets at particle level
- const std::unique_ptr<ParticleLevelRCJetObjectLoader>& rcjet_particle() const { return m_rcjet_particle;}
- const std::string& RCJetContainerParticle() const { return m_RCJetContainerParticle;} // name for RC jets container
-
// Truth tree inserted variables
// This can be expanded as required
// This is just a first pass at doing this sort of thing
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx
index 6d860f2e0167591df5c1b9420509850bbaf5aac3..425bdd2f6ad3f12c5b022a81efa929f8465d4fba 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx
@@ -336,6 +336,7 @@ int main(int argc, char** argv) {
// make top::Event objects
std::unique_ptr<top::TopEventMaker> topEventMaker( new top::TopEventMaker( "top::TopEventMaker" ) );
top::check(topEventMaker->setProperty( "config" , topConfig ) , "Failed to setProperty of top::TopEventMaker");
+ top::check(topEventMaker->initialize(),"Failed to initialize top::TopEventMaker");
// Debug messages?
// topEventMaker.msg().setLevel(MSG::DEBUG);
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/CMakeLists.txt b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/CMakeLists.txt
index 867e742c8c54e5d74640ced11785ace31ce8419e..320a69fa43dd44606ab78bf7f98276b1a200ed82 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/CMakeLists.txt
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/CMakeLists.txt
@@ -19,11 +19,19 @@ atlas_depends_on_subdirs( PUBLIC
xAODMissingET
xAODTracking
TopConfiguration
- TopPartons )
+ TopPartons
+ JetAnalysisInterfaces
+ JetSubStructureMomentTools
+ JetSubStructureUtils
+ JetReclustering )
# This package uses ROOT:
find_package( ROOT REQUIRED COMPONENTS Core Gpad Tree Hist RIO MathCore Graf )
+# Need fast jet for the RC jet substructure code
+find_package( FastJet COMPONENTS fastjetplugins fastjettools )
+find_package( FastJetContrib COMPONENTS EnergyCorrelator Nsubjettiness )
+
# Generate a CINT dictionary source file:
atlas_add_root_dictionary( TopEvent _cintDictSource
ROOT_HEADERS Root/LinkDef.h
@@ -49,6 +57,12 @@ atlas_add_library( TopEvent Root/*.cxx Root/*.h Root/*.icc
xAODTracking
TopConfiguration
TopPartons
+ JetAnalysisInterfacesLib
+ JetSubStructureMomentToolsLib
+ JetSubStructureUtils
+ JetReclusteringLib
+ ${FASTJET_LIBRARIES}
+ ${FASTJETCONTRIB_LIBRARIES}
${ROOT_LIBRARIES}
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} )
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/LinkDef.h b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/LinkDef.h
index 8c4eeedc0f27b2b94c79b389255b6604fe0304d9..5ec3eac5955fef21a8f87c8d192cdfa00d8d006d 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/LinkDef.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/LinkDef.h
@@ -5,11 +5,13 @@
#include "TopEvent/SystematicEvent.h"
#include "TopEvent/KLFitterResult.h"
#include "TopEvent/PseudoTopResult.h"
+#include "TopEvent/RCJetMC15.h"
#ifdef __CINT__
#pragma extra_include "TopEvent/SystematicEvent.h";
#pragma extra_include "TopEvent/KLFitterResult.h";
#pragma extra_include "TopEvent/PseudoTopResult.h";
+#pragma extra_include "TopEvent/RCJetMC15.h";
#pragma link off all globals;
#pragma link off all classes;
@@ -28,4 +30,6 @@
#pragma link C++ class xAOD::PseudoTopResultContainer+;
#pragma link C++ class xAOD::PseudoTopResultAuxContainer+;
+#pragma link C++ class RCJetMC15+;
+
#endif
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/RCJetMC15.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/RCJetMC15.cxx
similarity index 65%
rename from PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/RCJetMC15.cxx
rename to PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/RCJetMC15.cxx
index 1d72a5b623c10292174cd667372577cf9563cce8..e58de4fe73d20aaa2608cd039a812088693bd911 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/RCJetMC15.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/RCJetMC15.cxx
@@ -14,7 +14,7 @@
// File for initializing and making re-clustered jets.
//
***************************************************************/
-#include "TopObjectSelectionTools/RCJetMC15.h"
+#include "TopEvent/RCJetMC15.h"
#include "TopConfiguration/TopConfig.h"
#include "AsgTools/AsgTool.h"
@@ -174,7 +174,7 @@ StatusCode RCJetMC15::initialize(){
if (treeName.second.compare("nominal")!=0) hash_name = treeName.second; // no extra strings for nominal (so all other non-unique systs have same name as nominal)
- m_InputJetContainer = m_InJetContainerBase+hash_name+m_name;
+ m_InputJetContainer = m_InJetContainerBase+hash_name;
m_OutputJetContainer = m_OutJetContainerBase+hash_name+m_name;
// build a jet re-clustering tool for each case
@@ -216,7 +216,7 @@ StatusCode RCJetMC15::initialize(){
else{
- m_InputJetContainer = m_InJetContainerBase+m_name;
+ m_InputJetContainer = m_InJetContainerBase;
m_OutputJetContainer = m_OutJetContainerBase+m_name;
// map of container names to access in event saver
@@ -256,7 +256,6 @@ StatusCode RCJetMC15::execute(const top::Event& event) {
// 22 Feb 2016:
// Code significantly shortened to make this container
// thanks to email exchange between Davide Gerbaudo & Attila Krasznahorkay
-
typedef ConstDataVector< xAOD::JetContainer > CJets;
std::unique_ptr< CJets > rcjets( new CJets( event.m_jets.begin(), event.m_jets.end(), SG::VIEW_ELEMENTS ) );
top::check( evtStore()->tds()->record( std::move( rcjets ), m_InputJetContainer ), "Failed to put jets in TStore for re-clustering" );
@@ -276,211 +275,138 @@ StatusCode RCJetMC15::execute(const top::Event& event) {
tool_iter->second->execute();
else
m_jetReclusteringTool.at(hash_factor*m_config->nominalHashValue())->execute();
+
+ xAOD::JetContainer* myJets(nullptr);
+ top::check(evtStore()->retrieve(myJets,m_OutputJetContainer),"Failed to retrieve RC JetContainer");
+ for (auto rcjet : *myJets){
+ rcjet->auxdecor<bool>("PassedSelection") = passSelection(*rcjet);
+ }
+
+
+ if (m_useJSS){
+
+ // get the subjets and clusters of the rcjets
+ for (auto rcjet : *myJets){
+
+ std::vector<fastjet::PseudoJet> clusters;
+ // // //EMTOPO CLUSTERS
+ getEMTopoClusters(clusters,rcjet);
+ // // // LCTOPO CLUSTERS
+ //getLCTopoClusters(clusters,rcjet);
+
+ // 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!
+ double tau32 = -1, tau21 = -1, D2 = -1;
+
+ double tau1 = m_nSub1_beta1->result(correctedJet);
+ double tau2 = m_nSub2_beta1->result(correctedJet);
+ double tau3 = m_nSub3_beta1->result(correctedJet);
+
+ if(fabs(tau1) > 1e-8)
+ tau21 = tau2/tau1;
+ else
+ tau21 = -999.0;
+ if(fabs(tau2) > 1e-8)
+ tau32 = tau3/tau2;
+ else
+ tau32 = -999.0;
+
+ double vECF1 = m_ECF1->result(correctedJet);
+ double vECF2 = m_ECF2->result(correctedJet);
+ double vECF3 = m_ECF3->result(correctedJet);
+ if(fabs(vECF2) > 1e-8)
+ D2 = vECF3 * pow(vECF1,3) / pow(vECF2,3);
+ else
+ D2 = -999.0;
+
+ double split12 = m_split12->result(correctedJet);
+ double split23 = m_split23->result(correctedJet);
+ double qw = m_qw->result(correctedJet);
+
+
+ // 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 (fabs(gECF212) > 1e-12){
+ L1 = gECF321 / (pow(gECF212,(1.0)));
+ L2 = gECF331 / (pow(gECF212,(3.0/2.0)));
+ }
+ if (fabs(gECF331) > 1e-12){
+ L3 = gECF311 / (pow(gECF331,(1.0/3.0) ));
+ L4 = gECF322 / (pow(gECF331,(4.0/3.0) ));
+ }
+ if (fabs(gECF441) > 1e-12){
+ L5 = gECF422 / (pow(gECF441, (1.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>("D2_clstr") = D2;
+ rcjet->auxdecor<float>("ECF1_clstr") = vECF1;
+ rcjet->auxdecor<float>("ECF2_clstr") = vECF2;
+ rcjet->auxdecor<float>("ECF3_clstr") = vECF3;
+
+ rcjet->auxdecor<float>("d12_clstr") = split12;
+ rcjet->auxdecor<float>("d23_clstr") = split23;
+ rcjet->auxdecor<float>("Qw_clstr") = qw;
+
+ 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 rcjet loop
+ }// end of if useJSS
+
+
}
- if (m_useJSS){
- // get the rcjets
- xAOD::JetContainer* myJets(nullptr);
- top::check(evtStore()->retrieve(myJets,m_OutputJetContainer),"Failed to retrieve RC JetContainer");
- // // get the clusters (directly we so can try using the LCTopo clusters)
- //const xAOD::CaloClusterContainer* myClusters(nullptr);
- //top::check(evtStore()->retrieve(myClusters, "CaloCalTopoClusters"), "Failed to retrieve CaloCalTopoClusters");
-
- // get the subjets and clusters of the rcjets
- for (auto rcjet : *myJets){
- std::vector<fastjet::PseudoJet> clusters;
- clusters.clear();
- // for (auto cluster : *myClusters){
- // //const xAOD::CaloCluster* cluster_raw = static_cast<const xAOD::CaloCluster*>(cluster->rawConstituents());
-
- // // std::cout << "Cluster E default = " << cluster->e() << std::endl;
- // // std::cout << "Cluster Eta default = " << cluster->eta() << std::endl;
-
- // // std::cout <<" Cluster E getCalE = " << cluster->calE() << std::endl;
- // // std::cout <<" Cluster E getCalEta = " << cluster->calEta() << std::endl;
-
- // // std::cout <<"Cluster cal pt(E) = " << cluster->e(xAOD::CaloCluster_v1::State(1)) << std::endl;
-
-
- // LCTOPO CLUSTERS
-
- // for (auto cluster : *myClusters){
- // for (auto subjet : rcjet->getConstituents() ){
- // const xAOD::Jet* subjet_raw = static_cast<const xAOD::Jet*>(subjet->rawConstituent());
-
- // // //const xAOD::CaloCluster* cluster_raw = static_cast<const xAOD::CaloCluster*>(cluster->rawConstituents());
- // float dR = subjet_raw->p4().DeltaR(cluster->p4());
- // if (dR < 0.4){
- // TLorentzVector temp_p4;
- // temp_p4.SetPtEtaPhiE(cluster->pt((xAOD::CaloCluster_v1::State(1))), cluster->eta((xAOD::CaloCluster_v1::State(1))), cluster->phi((xAOD::CaloCluster_v1::State(1))), cluster->e((xAOD::CaloCluster_v1::State(1))));
- // clusters.push_back(fastjet::PseudoJet(temp_p4.Px(),temp_p4.Py(),temp_p4.Pz(),temp_p4.E()));
- // break;
- // }
- // }
- // }
-
-
- // // //SCALED CLUSTERS
- for (auto subjet : rcjet->getConstituents() ){
- const xAOD::Jet* subjet_raw = static_cast<const xAOD::Jet*>(subjet->rawConstituent());
-
- // Make sure we don't try to access jets that have had the clusters thinned
- bool hasConstituents=true;
- auto links = subjet_raw->constituentLinks();
- for(auto link : links) {
- if(!link.isValid()) {
- ATH_MSG_WARNING("Some of the RC Jet Constituents have been thinned - will not be included in RCJet JSS calculation");
- hasConstituents=false;
- break;
- }
- }
- if (!hasConstituents) {continue;}
+
- for (auto clus_itr : subjet_raw->getConstituents() ){
- if(clus_itr->e() > 0){
- TLorentzVector temp_p4;
-
- //std::cout << "SubJetRaw ConstitScale pt = " << subjet_raw->jetP4(xAOD::JetEMScaleMomentum).pt() << std::endl;
- //std::cout << "SubJetRaw ConstitScale pt = " << subjet_raw->jetP4(xAOD::JetConstitScaleMomentum).pt() << std::endl;
- //std::cout << "SubJetRaw calibrated pt = " << subjet_raw->jetP4().pt() << std::endl;
-
- // const xAOD::CaloCluster* clus_raw = static_cast<const xAOD::CaloCluster*>(clus_itr->rawConstituent());
-
- // std::cout << "Cluster is a " << typeid(clus_itr).name() << std::endl;
- // std::cout <<"Cluster E = " << clus_itr->e() << std::endl;
- //std::cout <<"Cluster raw type = " << clus_raw->type() << std::endl;
-
-
- //std::cout <<"Cluster EM E = " << clus_raw->getRawE() << std::endl;
- //std::cout <<"Cluster LC E = " << clus_raw->calE() << std::endl;
-
- //std::cout <<"Cluster raw pt = " << clus_itr->pt(xAOD::CaloCluster_v1::State(0)) << std::endl;
- //std::cout <<"Cluster cal pt = " << clus_raw->pt(xAOD::CaloCluster_v1::State(1)) << std::endl;
-
- //std::cout <<"Cluster EM pT = " << clus_itr->pt() << std::endl;
- //std::cout <<"Cluster LC pT = " << clus_itr->pt() << std::endl;
-
- // double sf = subjet_raw->jetP4().pt() / subjet_raw->jetP4(xAOD::JetEMScaleMomentum).pt();
- double sf = 1.0;
- //std::cout << "SF = " << sf << std::endl;
- temp_p4.SetPtEtaPhiM(clus_itr->pt()*sf, clus_itr->eta(), clus_itr->phi(), clus_itr->m());
-
- clusters.push_back(fastjet::PseudoJet(temp_p4.Px(),temp_p4.Py(),temp_p4.Pz(),temp_p4.E()));
- }
- }
- }
-
-
-
-
-
- // 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!
- double tau32 = -1, tau21 = -1, D2 = -1;
-
- double tau1 = m_nSub1_beta1->result(correctedJet);
- double tau2 = m_nSub2_beta1->result(correctedJet);
- double tau3 = m_nSub3_beta1->result(correctedJet);
-
- if(fabs(tau1) > 1e-8)
- tau21 = tau2/tau1;
- else
- tau21 = -999.0;
- if(fabs(tau2) > 1e-8)
- tau32 = tau3/tau2;
- else
- tau32 = -999.0;
-
- double vECF1 = m_ECF1->result(correctedJet);
- double vECF2 = m_ECF2->result(correctedJet);
- double vECF3 = m_ECF3->result(correctedJet);
- if(fabs(vECF2) > 1e-8)
- D2 = vECF3 * pow(vECF1,3) / pow(vECF2,3);
- else
- D2 = -999.0;
-
- double split12 = m_split12->result(correctedJet);
- double split23 = m_split23->result(correctedJet);
- double qw = m_qw->result(correctedJet);
-
-
- // 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 (fabs(gECF212) > 1e-12){
- L1 = gECF321 / (pow(gECF212,(1.0)));
- L2 = gECF331 / (pow(gECF212,(3.0/2.0)));
- }
- if (fabs(gECF331) > 1e-12){
- L3 = gECF311 / (pow(gECF331,(1.0/3.0) ));
- L4 = gECF322 / (pow(gECF331,(4.0/3.0) ));
- }
- if (fabs(gECF441) > 1e-12){
- L5 = gECF422 / (pow(gECF441, (1.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>("D2_clstr") = D2;
- rcjet->auxdecor<float>("ECF1_clstr") = vECF1;
- rcjet->auxdecor<float>("ECF2_clstr") = vECF2;
- rcjet->auxdecor<float>("ECF3_clstr") = vECF3;
-
- rcjet->auxdecor<float>("d12_clstr") = split12;
- rcjet->auxdecor<float>("d23_clstr") = split23;
- rcjet->auxdecor<float>("Qw_clstr") = qw;
-
- 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 rcjet loop
- }// end of if useJSS
@@ -589,3 +515,100 @@ bool RCJetMC15::passSelection(const xAOD::Jet& jet) const {
return true;
}
+
+void RCJetMC15::getEMTopoClusters(std::vector<fastjet::PseudoJet>& clusters,const xAOD::Jet* rcjet){
+
+ clusters.clear();
+
+ for (auto subjet : rcjet->getConstituents() ){
+ const xAOD::Jet* subjet_raw = static_cast<const xAOD::Jet*>(subjet->rawConstituent());
+
+ // Make sure we don't try to access jets that have had the clusters thinned
+ bool hasConstituents=true;
+ auto links = subjet_raw->constituentLinks();
+ for(auto link : links) {
+ if(!link.isValid()) {
+ ATH_MSG_WARNING("Some of the RC Jet Constituents have been thinned - will not be included in RCJet JSS calculation");
+ hasConstituents=false;
+ break;
+ }
+ }
+ if (!hasConstituents) {continue;}
+
+ for (auto clus_itr : subjet_raw->getConstituents() ){
+ if(clus_itr->e() > 0){
+ TLorentzVector temp_p4;
+
+ //std::cout << "SubJetRaw ConstitScale pt = " << subjet_raw->jetP4(xAOD::JetEMScaleMomentum).pt() << std::endl;
+ //std::cout << "SubJetRaw ConstitScale pt = " << subjet_raw->jetP4(xAOD::JetConstitScaleMomentum).pt() << std::endl;
+ //std::cout << "SubJetRaw calibrated pt = " << subjet_raw->jetP4().pt() << std::endl;
+
+ // const xAOD::CaloCluster* clus_raw = static_cast<const xAOD::CaloCluster*>(clus_itr->rawConstituent());
+
+ // std::cout << "Cluster is a " << typeid(clus_itr).name() << std::endl;
+ // std::cout <<"Cluster E = " << clus_itr->e() << std::endl;
+ //std::cout <<"Cluster raw type = " << clus_raw->type() << std::endl;
+
+
+ //std::cout <<"Cluster EM E = " << clus_raw->getRawE() << std::endl;
+ //std::cout <<"Cluster LC E = " << clus_raw->calE() << std::endl;
+
+ //std::cout <<"Cluster raw pt = " << clus_itr->pt(xAOD::CaloCluster_v1::State(0)) << std::endl;
+ //std::cout <<"Cluster cal pt = " << clus_raw->pt(xAOD::CaloCluster_v1::State(1)) << std::endl;
+
+ //std::cout <<"Cluster EM pT = " << clus_itr->pt() << std::endl;
+ //std::cout <<"Cluster LC pT = " << clus_itr->pt() << std::endl;
+
+ // double sf = subjet_raw->jetP4().pt() / subjet_raw->jetP4(xAOD::JetEMScaleMomentum).pt();
+ double sf = 1.0;
+ //std::cout << "SF = " << sf << std::endl;
+ temp_p4.SetPtEtaPhiM(clus_itr->pt()*sf, clus_itr->eta(), clus_itr->phi(), clus_itr->m());
+
+ clusters.push_back(fastjet::PseudoJet(temp_p4.Px(),temp_p4.Py(),temp_p4.Pz(),temp_p4.E()));
+ }
+ }
+ }
+
+}
+
+void RCJetMC15::getLCTopoClusters(std::vector<fastjet::PseudoJet>& clusters,const xAOD::Jet* rcjet){
+
+ //LCTOPO CLUSTERS
+ clusters.clear();
+
+ // get the clusters (directly we so can try using the LCTopo clusters)
+ const xAOD::CaloClusterContainer* myClusters(nullptr);
+ top::check(evtStore()->retrieve(myClusters, "CaloCalTopoClusters"), "Failed to retrieve CaloCalTopoClusters");
+
+
+
+ for (auto cluster : *myClusters){
+
+
+ //const xAOD::CaloCluster* cluster_raw = static_cast<const xAOD::CaloCluster*>(cluster->rawConstituents());
+
+ // std::cout << "Cluster E default = " << cluster->e() << std::endl;
+ // std::cout << "Cluster Eta default = " << cluster->eta() << std::endl;
+
+ // std::cout <<" Cluster E getCalE = " << cluster->calE() << std::endl;
+ // std::cout <<" Cluster E getCalEta = " << cluster->calEta() << std::endl;
+
+ // std::cout <<"Cluster cal pt(E) = " << cluster->e(xAOD::CaloCluster_v1::State(1)) << std::endl;
+
+ for (auto subjet : rcjet->getConstituents() ){
+ const xAOD::Jet* subjet_raw = static_cast<const xAOD::Jet*>(subjet->rawConstituent());
+
+ // //const xAOD::CaloCluster* cluster_raw = static_cast<const xAOD::CaloCluster*>(cluster->rawConstituents());
+ float dR = subjet_raw->p4().DeltaR(cluster->p4());
+ if (dR < 0.4){
+ TLorentzVector temp_p4;
+ temp_p4.SetPtEtaPhiE(cluster->pt((xAOD::CaloCluster_v1::State(1))), cluster->eta((xAOD::CaloCluster_v1::State(1))), cluster->phi((xAOD::CaloCluster_v1::State(1))), cluster->e((xAOD::CaloCluster_v1::State(1))));
+ clusters.push_back(fastjet::PseudoJet(temp_p4.Px(),temp_p4.Py(),temp_p4.Pz(),temp_p4.E()));
+ break;
+ }
+ }
+ }
+
+}
+
+
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/TopEventMaker.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/TopEventMaker.cxx
index c2f5c60cf706f32925433d9fcc9a491c64d89e1b..f7bfff733c047532c73757f74fcabb26287f33aa 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/TopEventMaker.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/Root/TopEventMaker.cxx
@@ -23,6 +23,10 @@
#include "TopPartons/PartonHistory.h"
+#include <boost/algorithm/string.hpp>
+
+#include <iostream>
+
namespace top {
TopEventMaker::TopEventMaker( const std::string& name ) :
@@ -32,6 +36,34 @@ namespace top {
declareProperty( "config" , m_config );
}
+ StatusCode TopEventMaker::initialize(){
+ if (m_config->useRCJets() == true){
+ m_rc=std::unique_ptr<RCJetMC15> ( new RCJetMC15( "RCJetMC15" ) );
+ top::check(m_rc->setProperty( "config" , m_config ) , "Failed to set config property of RCJetMC15");
+ top::check(m_rc->initialize(),"Failed to initialize RCJetMC15");
+ }
+
+ if (m_config->useVarRCJets() == true){
+ boost::split(m_VarRCJetRho, m_config->VarRCJetRho(), boost::is_any_of(","));
+ boost::split(m_VarRCJetMassScale, m_config->VarRCJetMassScale(), boost::is_any_of(","));
+
+ for (auto& rho : m_VarRCJetRho){
+ for (auto& mass_scale : m_VarRCJetMassScale){
+ std::replace( rho.begin(), rho.end(), '.', '_');
+ std::string name = rho+mass_scale;
+ m_VarRC[name] = std::unique_ptr<RCJetMC15> ( new RCJetMC15( "VarRCJetMC15_"+name ) );
+ top::check(m_VarRC[name]->setProperty( "config" , m_config ) , "Failed to set config property of VarRCJetMC15");
+ top::check(m_VarRC[name]->setProperty( "VarRCjets", true ) , "Failed to set VarRCjets property of VarRCJetMC15");
+ top::check(m_VarRC[name]->setProperty( "VarRCjets_rho", rho ) , "Failed to set VarRCjets rho property of VarRCJetMC15");
+ top::check(m_VarRC[name]->setProperty( "VarRCjets_mass_scale", mass_scale ) , "Failed to set VarRCjets mass scale property of VarRCJetMC15");
+ top::check(m_VarRC[name]->initialize(),"Failed to initialize VarRCJetMC15");
+ } // end loop over mass scale parameters (e.g., top mass, w mass, etc.)
+ } // end loop over mass scale multiplies (e.g., 1.,2.,etc.)
+ }
+
+ return StatusCode::SUCCESS;
+ }
+
/// As top-xaod isn't an asg::AsgTool, it doesn't have access to all the information
/// Very annoying, as it's actually quite simple
const xAOD::SystematicEventContainer* TopEventMaker::systematicEvents(const std::string& sgKey) const
@@ -249,6 +281,7 @@ namespace top {
xAOD::JetContainer* calibratedJetsTDS(nullptr);
top::check(evtStore()->retrieve(calibratedJetsTDS, m_config->sgKeyJetsTDS(hash,looseJets) ), "Failed to retrieve taus");
+
for (auto index : currentSystematic.goodJets()) {
auto jet = calibratedJetsTDS->at(index);
@@ -275,6 +308,44 @@ namespace top {
}
+ // Reclustered jets
+ if (m_config->useRCJets()){
+ top::check(m_rc->execute(event),"Failed to execute RCJetMC15 container");
+ std::string rcJetContainerName = m_rc->rcjetContainerName(event.m_hashValue,event.m_isLoose);
+ const xAOD::JetContainer* rc_jets(nullptr);
+ top::check(evtStore()->retrieve(rc_jets,rcJetContainerName),"Failed to retrieve RC JetContainer");
+ //Object selection
+ for (auto rcjet : *rc_jets){
+ top::check( rcjet->isAvailable<bool>("PassedSelection") , " Can't find jet decoration \"PassedSelection\" - we need it to decide if we should keep the reclustered jet in the top::Event instance or not!");
+ if(rcjet->auxdataConst<bool>("PassedSelection"))event.m_RCJets.push_back((xAOD::Jet*)rcjet);
+ }
+ }
+ // Variable-R reclustered jets
+ if (m_config->useVarRCJets()){
+ for (auto& rho : m_VarRCJetRho){
+ for (auto& mass_scale : m_VarRCJetMassScale){
+ std::replace( rho.begin(), rho.end(), '.', '_');
+ std::string name = rho+mass_scale;
+ top::check(m_VarRC[name]->execute(event),"Failed to execute RCJetMC15 container");
+
+ // Get the name of the container of re-clustered jets in TStore
+ std::string varRCJetContainerName = m_VarRC[name]->rcjetContainerName(event.m_hashValue,event.m_isLoose);
+
+ // -- Retrieve the re-clustered jets from TStore & save good re-clustered jets -- //
+ const xAOD::JetContainer* vrc_jets(nullptr);
+ top::check(evtStore()->retrieve(vrc_jets,varRCJetContainerName),"Failed to retrieve RC JetContainer");
+
+ event.m_VarRCJets[name] = std::make_shared<xAOD::JetContainer>(SG::VIEW_ELEMENTS);
+ for (auto vrcjet : *vrc_jets){
+ top::check( vrcjet->isAvailable<bool>("PassedSelection") , " Can't find jet decoration \"PassedSelection\" - we need it to decide if we should keep the variable-R reclustered jet in the top::Event instance or not!");
+ if(vrcjet->auxdataConst<bool>("PassedSelection"))event.m_VarRCJets[name]->push_back((xAOD::Jet*)vrcjet);
+ }
+
+ }
+ }
+ }
+
+
//large-R jets
if (m_config->useLargeRJets()) {
///-- Need to read const collections for mini-xaod read back --///
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/Event.h b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/Event.h
index 6216e0afec9220123e6cfef6b5cdcf7eca21b981..44e6b1ea348d1d434a4de1ce169a95af8d70fac1 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/Event.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/Event.h
@@ -54,6 +54,7 @@ public:
m_jets(SG::VIEW_ELEMENTS),
m_photons(SG::VIEW_ELEMENTS),
m_largeJets(SG::VIEW_ELEMENTS),
+ m_RCJets(SG::VIEW_ELEMENTS),
m_trackJets(SG::VIEW_ELEMENTS),
m_tauJets(SG::VIEW_ELEMENTS),
m_met(nullptr),
@@ -96,6 +97,12 @@ public:
///Container of large jets (can be sorted)
xAOD::JetContainer m_largeJets;
+ ///Container of recluster jets (can be sorted)
+ xAOD::JetContainer m_RCJets;
+
+ /// Containers of variable-R reclustered jets (can be sorted)
+ std::unordered_map< std::string,std::shared_ptr<xAOD::JetContainer> > m_VarRCJets;
+
///Container of track jets (can be sorted)
xAOD::JetContainer m_trackJets;
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/TopObjectSelectionTools/RCJetMC15.h b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/RCJetMC15.h
similarity index 95%
rename from PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/TopObjectSelectionTools/RCJetMC15.h
rename to PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/RCJetMC15.h
index 3b790eebe74e1612d6984ddf374e4682db8bbd58..caef7a1c7aaa938892563fae5d8a0f63ab1e5196 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/TopObjectSelectionTools/RCJetMC15.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/RCJetMC15.h
@@ -82,6 +82,9 @@ public:
std::string rcjetContainerName( std::size_t hash_value, bool isLooseEvent );
bool passSelection(const xAOD::Jet& jet) const;
+
+ void getEMTopoClusters(std::vector<fastjet::PseudoJet>& clusters,const xAOD::Jet* rcjet);
+ void getLCTopoClusters(std::vector<fastjet::PseudoJet>& clusters,const xAOD::Jet* rcjet);
private:
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/TopEventMaker.h b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/TopEventMaker.h
index fb9c853f660733c4080ef5a5c308430e291a0687..77c01d47d56aba52f61850e6f497c361dca3ce52 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/TopEventMaker.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEvent/TopEvent/TopEventMaker.h
@@ -27,15 +27,21 @@
#include "TopEvent/Event.h"
#include "TopEvent/SystematicEvent.h"
+#include "TopEvent/RCJetMC15.h"
+
+
+class RCJetMC15;
namespace top {
class TopConfig;
+
class TopEventMaker final : public asg::AsgTool {
public:
explicit TopEventMaker( const std::string& name );
virtual ~TopEventMaker(){}
+ StatusCode initialize();
// Delete Standard constructors
TopEventMaker(const TopEventMaker& rhs) = delete;
@@ -66,6 +72,10 @@ namespace top {
private:
std::shared_ptr<top::TopConfig> m_config;
+ std::unique_ptr<RCJetMC15> m_rc;
+ std::map<std::string,std::unique_ptr<RCJetMC15> > m_VarRC;
+ std::vector<std::string> m_VarRCJetRho;
+ std::vector<std::string> m_VarRCJetMassScale;
};
}
#endif
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/NRCJetSelector.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/NRCJetSelector.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5cf907a53dc48b1872d4ecf4f6ddea6de157d6f8
--- /dev/null
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/NRCJetSelector.cxx
@@ -0,0 +1,26 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "TopEventSelectionTools/NRCJetSelector.h"
+
+namespace top{
+
+ NRCJetSelector::NRCJetSelector(const std::string& params) :
+ SignValueSelector("RCJET_N", params, true) {
+ checkMultiplicityIsInteger();
+ }
+
+ bool NRCJetSelector::apply(const top::Event& event) const {
+ auto func = [&](const xAOD::Jet* jetPtr){return jetPtr->pt() > value();};
+ auto count = std::count_if(event.m_RCJets.begin(), event.m_RCJets.end(), func);
+ return checkInt(count, multiplicity());
+ }
+
+ bool NRCJetSelector::applyParticleLevel(const top::ParticleLevelEvent& event) const {
+ auto func = [&](const xAOD::Jet* jetPtr){return jetPtr->pt() > value();};
+ auto count = std::count_if(event.m_RCJets.begin(), event.m_RCJets.end(), func);
+ return checkInt(count, multiplicity());
+ }
+
+}
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/NVarRCJetSelector.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/NVarRCJetSelector.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b783f61da80b09ecc195db220313ede2d463c169
--- /dev/null
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/NVarRCJetSelector.cxx
@@ -0,0 +1,34 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "TopEventSelectionTools/NVarRCJetSelector.h"
+#include "TopEvent/EventTools.h"
+#include <iostream>
+
+namespace top{
+
+ NVarRCJetSelector::NVarRCJetSelector(const std::string& name,const std::string& params) :
+ SignValueSelector("VRCJET_N " + name, params, true),
+ m_name(name){
+ checkMultiplicityIsInteger();
+ std::cout << "NVarRCJetSelector::NVarRCJetSelector: m_name = " << m_name << std::endl;
+ }
+
+ bool NVarRCJetSelector::apply(const top::Event& event) const {
+ auto func = [&](const xAOD::Jet* jetPtr){return jetPtr->pt() > value();};
+ top::check(event.m_VarRCJets.find(m_name)!=event.m_VarRCJets.end(),"Error in NVarRCJetSelector: Variable-R reclustered jets with parameter " + m_name + " not defined!" );
+ std::unordered_map< std::string,std::shared_ptr<xAOD::JetContainer> > VarRCJets=event.m_VarRCJets;
+ auto count = std::count_if(VarRCJets[m_name]->begin(), VarRCJets[m_name]->end(), func);
+ return checkInt(count, multiplicity());
+ }
+
+ bool NVarRCJetSelector::applyParticleLevel(const top::ParticleLevelEvent& event) const {
+ auto func = [&](const xAOD::Jet* jetPtr){return jetPtr->pt() > value();};
+ top::check(event.m_VarRCJets.find(m_name)!=event.m_VarRCJets.end(),"Error in NVarRCJetSelector: Variable-R reclustered jets with parameter " + m_name + " not defined!" );
+ std::unordered_map< std::string,std::shared_ptr<xAOD::JetContainer> > VarRCJets=event.m_VarRCJets;
+ auto count = std::count_if(VarRCJets[m_name]->begin(), VarRCJets[m_name]->end(), func);
+ return checkInt(count, multiplicity());
+ }
+
+}
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/TopEventSelectionToolsLoader.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/TopEventSelectionToolsLoader.cxx
index cf430d859bf039e639967b6e1059e052427d2d88..5cfcb28a8495dbcafef1a34809b0f96577430840 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/TopEventSelectionToolsLoader.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/Root/TopEventSelectionToolsLoader.cxx
@@ -51,6 +51,8 @@
#include "TopEventSelectionTools/TrigMatchSelector.h"
#include "TopEventSelectionTools/RunNumberSelector.h"
#include "TopEventSelectionTools/NLargeJetSelector.h"
+#include "TopEventSelectionTools/NRCJetSelector.h"
+#include "TopEventSelectionTools/NVarRCJetSelector.h"
#include "TopConfiguration/TopConfig.h"
@@ -77,6 +79,14 @@ namespace top {
return new top::NJetBtagSelector(param,config);
else if (toolname == "LJET_N")
return new top::NLargeJetSelector(param);
+ else if (toolname == "RCJET_N")
+ return new top::NRCJetSelector(param);
+ else if (toolname == "VRCJET_N"){
+ std::istringstream is(param);
+ std::string name;
+ getline(is,name,' ' );
+ return new top::NVarRCJetSelector(name,param.substr(name.size()+1));
+ }
else if (toolname == "MV1_N")
return new top::MV1Selector(param);
else if (toolname == "MV2C20_N")
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/TopEventSelectionTools/NRCJetSelector.h b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/TopEventSelectionTools/NRCJetSelector.h
new file mode 100644
index 0000000000000000000000000000000000000000..d336648f4c738606655cf222b155c0d3e2a65070
--- /dev/null
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/TopEventSelectionTools/NRCJetSelector.h
@@ -0,0 +1,25 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef NRCJETSELECTOR_H_
+#define NRCJETSELECTOR_H_
+
+#include "TopEventSelectionTools/SignValueSelector.h"
+
+
+namespace top {
+
+ class NRCJetSelector : public top::SignValueSelector {
+
+ public:
+
+ explicit NRCJetSelector(const std::string& params);
+
+ bool apply(const top::Event& event) const override;
+ bool applyParticleLevel( const top::ParticleLevelEvent & ) const override;
+ };
+
+}
+
+#endif
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/TopEventSelectionTools/NVarRCJetSelector.h b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/TopEventSelectionTools/NVarRCJetSelector.h
new file mode 100644
index 0000000000000000000000000000000000000000..9cc1efbae30fdd82e09847a82716726867e5bc85
--- /dev/null
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopEventSelectionTools/TopEventSelectionTools/NVarRCJetSelector.h
@@ -0,0 +1,25 @@
+/*
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef NVARRCJETSELECTOR_H_
+#define NVARRCJETSELECTOR_H_
+
+#include "TopEventSelectionTools/SignValueSelector.h"
+
+
+namespace top {
+
+ class NVarRCJetSelector : public top::SignValueSelector {
+ std::string m_name;
+ public:
+
+ explicit NVarRCJetSelector(const std::string& name,const std::string& params);
+
+ bool apply(const top::Event& event) const override;
+ bool applyParticleLevel( const top::ParticleLevelEvent & ) const override;
+ };
+
+}
+
+#endif
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/CMakeLists.txt b/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/CMakeLists.txt
index 91647f83456fc5a0c46d6909ac5073cf1a34a68d..f457b7a388f96957b40a45f61b26418eab649408 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/CMakeLists.txt
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/CMakeLists.txt
@@ -32,11 +32,6 @@ atlas_depends_on_subdirs( PUBLIC
# This package uses ROOT:
find_package( ROOT REQUIRED COMPONENTS Core Gpad Tree Hist RIO MathCore Graf )
-
-# Need fast jet for the RC jet substructure code
-find_package( FastJet COMPONENTS fastjetplugins fastjettools )
-find_package( FastJetContrib COMPONENTS EnergyCorrelator Nsubjettiness )
-
# Custom definitions needed for this package:
add_definitions( -g )
@@ -74,8 +69,6 @@ atlas_add_library( TopObjectSelectionTools Root/*.cxx Root/*.h Root/*.icc
MuonAnalysisInterfacesLib
FTagAnalysisInterfacesLib
${ROOT_LIBRARIES}
- ${FASTJET_LIBRARIES}
- ${FASTJETCONTRIB_LIBRARIES}
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} )
# Install data files from the package:
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/LinkDef.h b/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/LinkDef.h
index 16b7731924d4c6b3964e22dd6122e13f19858f0b..fd41112fe26d0c513ae10eec69df2e9d9dc7625a 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/LinkDef.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopObjectSelectionTools/Root/LinkDef.h
@@ -2,17 +2,11 @@
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
-#include "TopObjectSelectionTools/RCJetMC15.h"
-
#ifdef __CINT__
-#pragma extra_include "TopObjectSelectionTools/RCJetMC15.h";
-
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclass;
-#pragma link C++ class RCJetMC15+;
-
#endif
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelLoader.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelLoader.cxx
index fd3d3c54690dbe4749a62106910d58eea80b4bff..fd3384a21d128d6e362fbba219c4047813499935 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelLoader.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelLoader.cxx
@@ -24,6 +24,8 @@
#include "TopConfiguration/TopConfig.h"
+#include <boost/algorithm/string.hpp>
+
// #define TOP_PARTICLE_LEVEL_DEBUG_OVERLAP_REMOVAL 1
namespace top {
@@ -152,7 +154,33 @@ namespace top {
std::cout << " " << std::setw( 20 ) << "DoOverlapRemoval Mu-Jet? " << std::setw( 5 ) << std::boolalpha << m_config->doParticleLevelOverlapRemovalMuJet() << '\n';
std::cout << " " << std::setw( 20 ) << "DoOverlapRemoval El-Jet? " << std::setw( 5 ) << std::boolalpha << m_config->doParticleLevelOverlapRemovalElJet() << '\n';
std::cout << " " << std::setw( 20 ) << "DoOverlapRemoval Jet-Photon? " << std::setw( 5 ) << std::boolalpha << m_config->doParticleLevelOverlapRemovalJetPhoton() << '\n';
-
+
+
+ if ( m_config->useRCJets()){
+ m_particleLevelRCJetObjectLoader = std::unique_ptr<ParticleLevelRCJetObjectLoader> ( new ParticleLevelRCJetObjectLoader( m_config ) );
+ top::check(m_particleLevelRCJetObjectLoader->initialize(),"Failed to initialize ParticleLevelRCJetObjectLoader");
+
+ }
+
+ if (m_config->useVarRCJets() == true){
+ boost::split(m_VarRCJetRho, m_config->VarRCJetRho(), boost::is_any_of(","));
+ boost::split(m_VarRCJetMassScale, m_config->VarRCJetMassScale(), boost::is_any_of(","));
+
+ for (auto& rho : m_VarRCJetRho){
+ for (auto& mass_scale : m_VarRCJetMassScale){
+ std::replace( rho.begin(), rho.end(), '.', '_');
+ std::string name = rho+mass_scale;
+ m_particleLevelVarRCJetObjectLoader[name] = std::unique_ptr<ParticleLevelRCJetObjectLoader> ( new ParticleLevelRCJetObjectLoader( m_config ) );
+ top::check(m_particleLevelVarRCJetObjectLoader[name]->setProperty( "VarRCjets", true ) , "Failed to set VarRCjets property of VarRCJetMC15");
+ top::check(m_particleLevelVarRCJetObjectLoader[name]->setProperty( "VarRCjets_rho", rho ) , "Failed to set VarRCjets rho property of VarRCJetMC15");
+ top::check(m_particleLevelVarRCJetObjectLoader[name]->setProperty( "VarRCjets_mass_scale", mass_scale ) , "Failed to set VarRCjets mass scale property of VarRCJetMC15");
+ top::check(m_particleLevelVarRCJetObjectLoader[name]->initialize(),"Failed to initialize VarRCJetMC15");
+ } // end loop over mass scale parameters (e.g., top mass, w mass, etc.)
+ } // end loop over mass scale multiplies (e.g., 1.,2.,etc.)
+ }
+
+
+
} else {
std::cout << "Particle level reconstruction is disabled." << '\n';
}
@@ -507,6 +535,49 @@ namespace top {
plEvent.m_largeRJets = m_config->useTruthLargeRJets() ? m_goodLargeRJets.get() : nullptr;
plEvent.m_met = m_config->useTruthMET() ? (* mets)[ "NonInt" ] : nullptr;
+ // Reclustered jets
+ if ( m_config->useRCJets() ){
+ top::check(m_particleLevelRCJetObjectLoader->execute(plEvent),"Failed to execute ParticleLevelRCJetObjectLoader container");
+ // Get the name of the container of re-clustered jets
+ std::string RCJetContainerNane = m_particleLevelRCJetObjectLoader->rcjetContainerName();
+
+ // -- Retrieve the re-clustered jets from TStore & save good re-clustered jets -- //
+ const xAOD::JetContainer* rc_jets(nullptr);
+ top::check(evtStore()->retrieve(rc_jets,RCJetContainerNane),"Failed to retrieve particle RC JetContainer");
+
+ for( auto rcjet : *rc_jets){
+ top::check( rcjet->isAvailable<bool>("PassedSelection") , " Can't find reclustered jet decoration \"PassedSelection\" - we need it to decide if we should keep the reclustered jet in the top::Event instance or not!");
+ if(rcjet->auxdataConst<bool>("PassedSelection"))plEvent.m_RCJets.push_back((xAOD::Jet*)rcjet);
+ }
+
+ }
+ // Variable-R reclustered jets
+ if (m_config->useVarRCJets()){
+ for (auto& rho : m_VarRCJetRho){
+ for (auto& mass_scale : m_VarRCJetMassScale){
+ std::replace( rho.begin(), rho.end(), '.', '_');
+ std::string name = rho+mass_scale;
+ top::check(m_particleLevelVarRCJetObjectLoader[name]->execute(plEvent),"Failed to execute RCJetMC15 container");
+
+ // Get the name of the container of re-clustered jets in TStore
+ std::string varRCJetContainerName = m_particleLevelVarRCJetObjectLoader[name]->rcjetContainerName();
+
+ // -- Retrieve the re-clustered jets from TStore & save good re-clustered jets -- //
+ const xAOD::JetContainer* vrc_jets(nullptr);
+ top::check(evtStore()->retrieve(vrc_jets,varRCJetContainerName),"Failed to retrieve RC JetContainer");
+
+ plEvent.m_VarRCJets[name] = std::make_shared<xAOD::JetContainer>(SG::VIEW_ELEMENTS);
+ for (auto vrcjet : *vrc_jets){
+ top::check( vrcjet->isAvailable<bool>("PassedSelection") , " Can't find jet decoration \"PassedSelection\" - we need it to decide if we should keep the variable-R reclustered jet in the top::Event instance or not!");
+ if(vrcjet->auxdataConst<bool>("PassedSelection"))plEvent.m_VarRCJets[name]->push_back((xAOD::Jet*)vrcjet);
+ }
+
+ }
+ }
+ }
+
+
+
return plEvent;
}
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelRCJetObjectLoader.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelRCJetObjectLoader.cxx
index beb4f794b91da5147ecc0a3d93f9e547de96a17d..3055352e437aa2f63d78aae596e52929825ebdcc 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelRCJetObjectLoader.cxx
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/Root/ParticleLevelRCJetObjectLoader.cxx
@@ -23,38 +23,44 @@
#include "JetSubStructureUtils/EnergyCorrelatorGeneralized.h"
#include "JetSubStructureUtils/EnergyCorrelator.h"
-
-ParticleLevelRCJetObjectLoader::ParticleLevelRCJetObjectLoader( const std::shared_ptr<top::TopConfig> & cfg ): asg::AsgTool( "ParticleLevelRCJetObjectLoader" ),
- m_config(cfg),
- m_ptcut(0.),
- m_etamax(0.),
- m_trim(0.),
- m_radius(0.),
- m_treeName("particleLevel"),
-
- m_InJetContainerBase( "AntiKt4TruthWZJets_RC"),
- m_OutJetContainerBase("AntiKt10RCTrim"),
- m_InputJetContainer("AntiKt4TruthWZJets_RC"),
- m_OutputJetContainer("AntiKt10RCTrim"),
- m_jet_def_rebuild(nullptr),
- m_nSub1_beta1(nullptr),
- m_nSub2_beta1(nullptr),
- m_nSub3_beta1(nullptr),
- m_ECF1(nullptr),
- m_ECF2(nullptr),
- m_ECF3(nullptr),
- m_split12(nullptr),
- m_split23(nullptr),
- m_qw(nullptr),
- m_gECF332(nullptr),
- m_gECF461(nullptr),
- m_gECF322(nullptr),
- m_gECF331(nullptr),
- m_gECF422(nullptr),
- m_gECF441(nullptr),
- m_gECF212(nullptr),
- m_gECF321(nullptr),
- m_gECF311(nullptr){}
+ParticleLevelRCJetObjectLoader::ParticleLevelRCJetObjectLoader( const std::shared_ptr<top::TopConfig> & cfg ): asg::AsgTool( "ParticleLevelRCJetObjectLoader" ),m_config(cfg)
+{
+ m_ptcut = 0.;
+ m_etamax=0.;
+ m_trim=0.;
+ m_radius=0.;
+ m_treeName = "particleLevel";
+ m_InJetContainerBase = "AntiKt4TruthWZJets_RC";
+ m_OutJetContainerBase = "AntiKt10RCTrim";
+ m_InputJetContainer = "AntiKt4TruthWZJets_RC";
+ m_OutputJetContainer = "AntiKt10RCTrim";
+ m_name = "";
+ m_useJSS = false;
+ m_jet_def_rebuild = nullptr;
+ m_nSub1_beta1 = nullptr;
+ m_nSub2_beta1 = nullptr;
+ m_nSub3_beta1 = nullptr;
+ m_ECF1 = nullptr;
+ m_ECF2 = nullptr;
+ m_ECF3 = nullptr;
+ m_split12 = nullptr;
+ m_split23 = nullptr;
+ m_qw = nullptr;
+ m_gECF332 = nullptr;
+ m_gECF461 = nullptr;
+ m_gECF322 = nullptr;
+ m_gECF331 = nullptr;
+ m_gECF422 = nullptr;
+ m_gECF441 = nullptr;
+ m_gECF212 = nullptr;
+ m_gECF321 = nullptr;
+ m_gECF311 = nullptr;
+
+
+ declareProperty( "VarRCjets", m_VarRCjets=false);
+ declareProperty( "VarRCjets_rho", m_VarRCjets_rho="");
+ declareProperty( "VarRCjets_mass_scale",m_VarRCjets_mass_scale="");
+}
ParticleLevelRCJetObjectLoader::~ParticleLevelRCJetObjectLoader() {}
@@ -63,22 +69,46 @@ StatusCode ParticleLevelRCJetObjectLoader::initialize(){
/* Initialize the re-clustered jets */
ATH_MSG_INFO(" Initializing particle level Re-clustered jets ");
- m_ptcut = m_config->RCJetPtcut() ; // for initialize [GeV] & passSelection
- m_etamax = m_config->RCJetEtacut() ; // for passSelection
- m_trim = m_config->RCJetTrimcut() ; // for initialize
- m_radius = m_config->RCJetRadius() ; // for initialize
+ m_name = m_VarRCjets_rho+m_VarRCjets_mass_scale;
+
+
+ if(m_VarRCjets){
+ m_ptcut = m_config->VarRCJetPtcut(); // 100 GeV
+ m_etamax = m_config->VarRCJetEtacut(); // 2.5
+ m_trim = m_config->VarRCJetTrimcut(); // 0.05 (5% jet pT)
+ m_radius = m_config->VarRCJetMaxRadius(); // 1.2 (min=0.4)
+ m_minradius = 0.4; // 0.4 default (until we have smaller jets!)
+ std::string original_rho(m_VarRCjets_rho);
+ std::replace( original_rho.begin(), original_rho.end(), '_', '.');
+ float rho = std::stof(original_rho);
+ float m_scale = mass_scales.at(m_VarRCjets_mass_scale);
+ m_massscale = rho*m_scale*1e-3;
+ }
+ else {
+ m_ptcut = m_config->RCJetPtcut() ; // for initialize [GeV] & passSelection
+ m_etamax = m_config->RCJetEtacut() ; // for passSelection
+ m_trim = m_config->RCJetTrimcut() ; // for initialize
+ m_radius = m_config->RCJetRadius() ; // for initialize
+ m_minradius = -1.0;
+ m_massscale = -1.0;
+ }
+
+
+
m_InputJetContainer = m_InJetContainerBase;
- m_OutputJetContainer = m_OutJetContainerBase;
+ m_OutputJetContainer = m_OutJetContainerBase + m_name;
// build a jet re-clustering tool for each case
- m_jetReclusteringTool = new JetReclusteringTool(m_treeName);
+ m_jetReclusteringTool = new JetReclusteringTool(m_treeName+m_name);
top::check(m_jetReclusteringTool->setProperty("InputJetContainer", m_InputJetContainer),"Failed inputjetcontainer initialize reclustering tool");
top::check(m_jetReclusteringTool->setProperty("OutputJetContainer", m_OutputJetContainer),"Failed outputjetcontainer initialize reclustering tool");
top::check(m_jetReclusteringTool->setProperty("ReclusterRadius", m_radius),"Failed re-clustering radius initialize reclustering tool");
top::check(m_jetReclusteringTool->setProperty("RCJetPtMin", m_ptcut*1e-3),"Failed ptmin [GeV] initialize reclustering tool");
top::check(m_jetReclusteringTool->setProperty("RCJetPtFrac", m_trim),"Failed pT fraction initialize reclustering tool");
+ top::check(m_jetReclusteringTool->setProperty("VariableRMinRadius", m_minradius),"Failed VarRC min radius initialize reclustering tool");
+ top::check(m_jetReclusteringTool->setProperty("VariableRMassScale", m_massscale),"Failed VarRC mass scale initialize reclustering tool");
top::check(m_jetReclusteringTool->initialize(),"Failed to initialize reclustering tool");
if (m_config->useRCJetSubstructure()){
@@ -157,137 +187,147 @@ StatusCode ParticleLevelRCJetObjectLoader::execute(const top::ParticleLevelEvent
// (do not re-make the 'nominal' jet container over & over again!)
if (!evtStore()->contains<xAOD::JetContainer>(m_OutputJetContainer)) {
m_jetReclusteringTool->execute();
- }
-
- if (m_useJSS){
- xAOD::JetContainer* myJets(nullptr);
- top::check(evtStore()->retrieve(myJets,m_OutputJetContainer),"Failed to retrieve Particle Level RC JetContainer");
- for (auto rcjet : *myJets){
- // std::cout <<"retrieved PL rcjets" << std::endl;
- std::vector<fastjet::PseudoJet> clusters;
- clusters.clear();
- for (auto subjet : rcjet->getConstituents() ){
- //std::cout <<"retrieved PL subjets" << std::endl;
- const xAOD::Jet* subjet_raw = static_cast<const xAOD::Jet*>(subjet->rawConstituent());
- //std::cout <<"retrieved PL rawjets" << std::endl;
-
- for (auto clus_itr : subjet_raw->getConstituents() ){
- // std::cout <<"retrieved PL clusters" << std::endl;
- TLorentzVector temp_p4;
- temp_p4.SetPtEtaPhiM(clus_itr->pt(), clus_itr->eta(), clus_itr->phi(), clus_itr->m());
- clusters.push_back(fastjet::PseudoJet(temp_p4.Px(),temp_p4.Py(),temp_p4.Pz(),temp_p4.E()));
- }
- }
-
- 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!
- double tau32 = -1, tau21 = -1, D2 = -1;
-
- double tau1 = m_nSub1_beta1->result(correctedJet);
- double tau2 = m_nSub2_beta1->result(correctedJet);
- double tau3 = m_nSub3_beta1->result(correctedJet);
-
- if(fabs(tau1) > 1e-8)
- tau21 = tau2/tau1;
- else
- tau21 = -999.0;
- if(fabs(tau2) > 1e-8)
- tau32 = tau3/tau2;
- else
- tau32 = -999.0;
-
- double vECF1 = m_ECF1->result(correctedJet);
- double vECF2 = m_ECF2->result(correctedJet);
- double vECF3 = m_ECF3->result(correctedJet);
- if(fabs(vECF2) > 1e-8)
- D2 = vECF3 * pow(vECF1,3) / pow(vECF2,3);
- else
- D2 = -999.0;
-
- double split12 = m_split12->result(correctedJet);
- double split23 = m_split23->result(correctedJet);
- double qw = m_qw->result(correctedJet);
-
- // 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 (fabs(gECF212) > 1e-8){
- L1 = gECF321 / (pow(gECF212,(1.0)));
- L2 = gECF331 / (pow(gECF212,(3.0/2.0)));
- }
- if (fabs(gECF331) > 1e-8){
- L3 = gECF311 / (pow(gECF331,(1.0/3.0) ));
- L4 = gECF322 / (pow(gECF331,(4.0/3.0) ));
- }
- if (fabs(gECF441) > 1e-8){
- L5 = gECF422 / (pow(gECF441, (1.0)));
+ const xAOD::JetContainer* myJets(nullptr);
+ top::check(evtStore()->retrieve(myJets,m_OutputJetContainer),"Failed to retrieve particle RC JetContainer");
+
+ for( auto rcjet : *myJets ){
+ rcjet->auxdecor<bool>("PassedSelection") = passSelection(*rcjet);
}
-
- // 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>("D2_clstr") = D2;
- rcjet->auxdecor<float>("ECF1_clstr") = vECF1;
- rcjet->auxdecor<float>("ECF2_clstr") = vECF2;
- rcjet->auxdecor<float>("ECF3_clstr") = vECF3;
-
- rcjet->auxdecor<float>("d12_clstr") = split12;
- rcjet->auxdecor<float>("d23_clstr") = split23;
- rcjet->auxdecor<float>("Qw_clstr") = qw;
-
- 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();
-
-
- // std::cout << "RCJet pT = " << rcjet->p4().Pt() << std::endl;
- // std::cout << "RRCJet pT= " << correctedJet.pt() << std::endl;
-
- } // end rcjet loop
- } // end if useJSS
+ if (m_useJSS){
+
+ for (auto rcjet : *myJets){
+ // std::cout <<"retrieved PL rcjets" << std::endl;
+ std::vector<fastjet::PseudoJet> clusters;
+ clusters.clear();
+
+ for (auto subjet : rcjet->getConstituents() ){
+ //std::cout <<"retrieved PL subjets" << std::endl;
+ const xAOD::Jet* subjet_raw = static_cast<const xAOD::Jet*>(subjet->rawConstituent());
+ //std::cout <<"retrieved PL rawjets" << std::endl;
+
+ for (auto clus_itr : subjet_raw->getConstituents() ){
+ // std::cout <<"retrieved PL clusters" << std::endl;
+ TLorentzVector temp_p4;
+ temp_p4.SetPtEtaPhiM(clus_itr->pt(), clus_itr->eta(), clus_itr->phi(), clus_itr->m());
+ clusters.push_back(fastjet::PseudoJet(temp_p4.Px(),temp_p4.Py(),temp_p4.Pz(),temp_p4.E()));
+ }
+ }
+
+ 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!
+ double tau32 = -1, tau21 = -1, D2 = -1;
+
+ double tau1 = m_nSub1_beta1->result(correctedJet);
+ double tau2 = m_nSub2_beta1->result(correctedJet);
+ double tau3 = m_nSub3_beta1->result(correctedJet);
+
+ if(fabs(tau1) > 1e-8)
+ tau21 = tau2/tau1;
+ else
+ tau21 = -999.0;
+ if(fabs(tau2) > 1e-8)
+ tau32 = tau3/tau2;
+ else
+ tau32 = -999.0;
+
+ double vECF1 = m_ECF1->result(correctedJet);
+ double vECF2 = m_ECF2->result(correctedJet);
+ double vECF3 = m_ECF3->result(correctedJet);
+ if(fabs(vECF2) > 1e-8)
+ D2 = vECF3 * pow(vECF1,3) / pow(vECF2,3);
+ else
+ D2 = -999.0;
+
+ double split12 = m_split12->result(correctedJet);
+ double split23 = m_split23->result(correctedJet);
+ double qw = m_qw->result(correctedJet);
+
+ // 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 (fabs(gECF212) > 1e-8){
+ L1 = gECF321 / (pow(gECF212,(1.0)));
+ L2 = gECF331 / (pow(gECF212,(3.0/2.0)));
+ }
+ if (fabs(gECF331) > 1e-8){
+ L3 = gECF311 / (pow(gECF331,(1.0/3.0) ));
+ L4 = gECF322 / (pow(gECF331,(4.0/3.0) ));
+ }
+ if (fabs(gECF441) > 1e-8){
+ L5 = gECF422 / (pow(gECF441, (1.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>("D2_clstr") = D2;
+ rcjet->auxdecor<float>("ECF1_clstr") = vECF1;
+ rcjet->auxdecor<float>("ECF2_clstr") = vECF2;
+ rcjet->auxdecor<float>("ECF3_clstr") = vECF3;
+
+ rcjet->auxdecor<float>("d12_clstr") = split12;
+ rcjet->auxdecor<float>("d23_clstr") = split23;
+ rcjet->auxdecor<float>("Qw_clstr") = qw;
+
+ 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();
+
+
+ // std::cout << "RCJet pT = " << rcjet->p4().Pt() << std::endl;
+ // std::cout << "RRCJet pT= " << correctedJet.pt() << std::endl;
+
+
+ } // end rcjet loop
+ } // end if useJSS
+
+
+ }
+
+
// // }
return StatusCode::SUCCESS;
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelEvent.h b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelEvent.h
index 2bf81aa2d37cb161a3ed58e85f2e248c4b3d6dc5..100645fa59c50fa8db9d4123acbc00c6629a3dd0 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelEvent.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelEvent.h
@@ -30,7 +30,8 @@ namespace top {
m_muons( nullptr ),
m_photons( nullptr ),
m_jets( nullptr ),
- m_largeRJets( nullptr ),
+ m_RCJets( SG::VIEW_ELEMENTS ),
+ m_largeRJets( nullptr ),
m_met( nullptr ),
m_selectionDecisions(){}
@@ -49,6 +50,12 @@ namespace top {
/// Pointer to truth level jets
const xAOD::JetContainer * m_jets;
+
+ ///Container of recluster jets (can be sorted)
+ xAOD::JetContainer m_RCJets;
+
+ /// Containers of variable-R reclustered jets (can be sorted)
+ std::unordered_map< std::string,std::shared_ptr<xAOD::JetContainer> > m_VarRCJets;
/// Pointer to the truth level large R jets.
const xAOD::JetContainer * m_largeRJets;
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelLoader.h b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelLoader.h
index bc4ad6980be6fbb122e2bd6711c852c2bc79fff4..ecc2d61b897c03b7e85653608be3527dba47dfb5 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelLoader.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelLoader.h
@@ -13,6 +13,7 @@
#include <memory>
#include "TopParticleLevel/ObjectSelectorBase.h"
+#include "TopParticleLevel/ParticleLevelRCJetObjectLoader.h"
#include "AsgTools/AsgTool.h"
@@ -100,6 +101,11 @@ namespace top {
std::unique_ptr<ObjectSelectorBase<xAOD::TruthParticle> > m_objectSelector_Photon;
std::unique_ptr<ObjectSelectorBase<xAOD::Jet> > m_objectSelector_Jet;
std::unique_ptr<ObjectSelectorBase<xAOD::Jet> > m_objectSelector_LargeRJet;
+ std::unique_ptr<ParticleLevelRCJetObjectLoader> m_particleLevelRCJetObjectLoader;
+ std::map<std::string,std::unique_ptr<ParticleLevelRCJetObjectLoader> > m_particleLevelVarRCJetObjectLoader;
+ std::vector<std::string> m_VarRCJetRho;
+ std::vector<std::string> m_VarRCJetMassScale;
+
private:
// The dressed leptons (shallow copies of the input containers)
std::unique_ptr<xAOD::TruthParticleContainer> m_electronsDressed;
@@ -121,7 +127,7 @@ namespace top {
std::unique_ptr<xAOD::JetContainer> m_goodJets;
std::unique_ptr<xAOD::JetAuxContainer> m_goodJetsAux;
- std::unique_ptr<xAOD::JetContainer> m_goodLargeRJets;
+ std::unique_ptr<xAOD::JetContainer> m_goodLargeRJets;
std::unique_ptr<xAOD::JetAuxContainer> m_goodLargeRJetsAux;
// Flag denoting whether the loader tool is active. Will be set by the
diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelRCJetObjectLoader.h b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelRCJetObjectLoader.h
index 5352c5e4d29bc6e6664d3313f9e49f9e7d90ef1c..26fb730073864dcd594acab8671036d15e96fd3f 100644
--- a/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelRCJetObjectLoader.h
+++ b/PhysicsAnalysis/TopPhys/xAOD/TopParticleLevel/TopParticleLevel/ParticleLevelRCJetObjectLoader.h
@@ -76,12 +76,21 @@ public:
std::string rcjetContainerName(){return m_OutputJetContainer; };
private:
+
+ std::string m_name;
const std::shared_ptr<top::TopConfig> & m_config;
+ bool m_VarRCjets;
+ std::string m_VarRCjets_rho;
+ std::string m_VarRCjets_mass_scale;
+
float m_ptcut; // in GeV
float m_etamax;
float m_trim;
float m_radius;
+ float m_minradius;
+ float m_massscale;
+
bool m_useJSS;
std::string m_treeName;
@@ -91,6 +100,9 @@ private:
std::string m_InputJetContainer;
std::string m_OutputJetContainer;
+
+
+
//Substructure tool definitions
fastjet::JetDefinition* m_jet_def_rebuild;
fastjet::contrib::Nsubjettiness* m_nSub1_beta1;
@@ -112,6 +124,13 @@ private:
JetSubStructureUtils::EnergyCorrelatorGeneralized* m_gECF321;
JetSubStructureUtils::EnergyCorrelatorGeneralized* m_gECF311;
+ std::map<std::string,float> mass_scales = {
+ {"m_t",172500.},
+ {"m_w",80385.},
+ {"m_z",91188.},
+ {"m_h",125090.}};
+
+
//re-clustered jets
// -> need unordered map for systematics
JetReclusteringTool* m_jetReclusteringTool;