Merge pull request #13 from SUEPPhysics/addRho

Add rho

.gitignore

@@ -0,0 +1,2 @@
+__init__.py
+*.pyc

README.md

@@ -2,9 +2,7 @@
 ```
 cmsrel CMSSW_10_6_0 #you can use CMSSW_11_1_0
 cd CMSSW_10_6_0/src
-mkdir PhysicsTools
-cd PhysicsTools
-git clone -b mods https://github.com/SUEPPhysics/SUEPScouting.git 
+git clone -b mods https://github.com/SUEPPhysics/SUEPScouting.git PhysicsTools/SUEPScouting
 ```
 
 # To setup and compile

plugins/BuildFile.xml

@@ -6,6 +6,8 @@
 <use name="FWCore/ServiceRegistry"/>
 <use name="JetMETCorrections/Algorithms"/>
 <use name="HLTrigger/HLTcore"/>
+<use name="DataFormats/Math"/>
+<use name="DataFormats/Scouting"/>
 <use name="fastjet"/>
 <use name="fastjet-contrib"/>
 <use name="root"/>

plugins/FixedGridRhoProducerFastjetScouting.cc

@@ -0,0 +1,44 @@
+#include "PhysicsTools/SUEPScouting/plugins/FixedGridRhoProducerFastjetScouting.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "DataFormats/Common/interface/View.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+
+using namespace std;
+
+FixedGridRhoProducerFastjetScouting::FixedGridRhoProducerFastjetScouting(const edm::ParameterSet& iConfig)
+    : bge_(iConfig.getParameter<double>("maxRapidity"), iConfig.getParameter<double>("gridSpacing")) {
+  pfCandidatesTag_ = iConfig.getParameter<edm::InputTag>("pfCandidatesTag");
+  electronsTag_ = iConfig.getParameter<edm::InputTag>("electronsTag");
+  produces<double>();
+
+  input_pfcoll_token_ = consumes<std::vector<ScoutingParticle>>(pfCandidatesTag_);
+  input_elecoll_token_ = consumes<std::vector<ScoutingElectron>>(electronsTag_);
+}
+
+FixedGridRhoProducerFastjetScouting::~FixedGridRhoProducerFastjetScouting() {}
+
+void FixedGridRhoProducerFastjetScouting::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  std::vector<fastjet::PseudoJet> inputs;
+
+  edm::Handle<std::vector<ScoutingParticle>> pfColl;
+  iEvent.getByToken(input_pfcoll_token_, pfColl);
+  for (std::vector<ScoutingParticle>::const_iterator ibegin = pfColl->begin(), iend = pfColl->end(), i = ibegin; i != iend;
+       ++i) {
+    math::PtEtaPhiMLorentzVector v{i->pt(), i->eta(), i->phi(), i->m()};
+    inputs.push_back(fastjet::PseudoJet(v.px(), v.py(), v.pz(), v.energy()));
+  }
+
+  //electrons are not in PF candidate collection
+  edm::Handle<std::vector<ScoutingElectron>> eleColl;
+  iEvent.getByToken(input_elecoll_token_, eleColl);
+  for (std::vector<ScoutingElectron>::const_iterator ibegin = eleColl->begin(), iend = eleColl->end(), i = ibegin; i != iend;
+       ++i) {
+    math::PtEtaPhiMLorentzVector v{i->pt(), i->eta(), i->phi(), i->m()};
+    inputs.push_back(fastjet::PseudoJet(v.px(), v.py(), v.pz(), v.energy()));
+  }
+  bge_.set_particles(inputs);
+  iEvent.put(std::make_unique<double>(bge_.rho()));
+}
+
+DEFINE_FWK_MODULE(FixedGridRhoProducerFastjetScouting);

plugins/FixedGridRhoProducerFastjetScouting.h

@@ -0,0 +1,27 @@
+#ifndef PhysicsTools_SUEPScouting_plugins_FixedGridRhoProducerFastjetScouting_h
+#define PhysicsTools_SUEPScouting_plugins_FixedGridRhoProducerFastjetScouting_h
+
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/Scouting/interface/ScoutingElectron.h"
+#include "DataFormats/Scouting/interface/ScoutingParticle.h"
+#include "fastjet/tools/GridMedianBackgroundEstimator.hh"
+
+class FixedGridRhoProducerFastjetScouting : public edm::stream::EDProducer<> {
+public:
+  explicit FixedGridRhoProducerFastjetScouting(const edm::ParameterSet& iConfig);
+  ~FixedGridRhoProducerFastjetScouting() override;
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  edm::InputTag pfCandidatesTag_;
+  edm::InputTag electronsTag_;
+  fastjet::GridMedianBackgroundEstimator bge_;
+
+  edm::EDGetTokenT<std::vector<ScoutingParticle>> input_pfcoll_token_;
+  edm::EDGetTokenT<std::vector<ScoutingElectron>> input_elecoll_token_;
+};
+
+#endif

plugins/ScoutingNanoAOD.cc

@@ -130,6 +130,7 @@ class ScoutingNanoAOD : public edm::one::EDAnalyzer<edm::one::SharedResources, e
   const edm::EDGetTokenT<std::vector<PileupSummaryInfo> >        pileupInfoToken;
   const edm::EDGetTokenT<std::vector<reco::GenParticle> >  	gensToken;
   //const edm::EDGetTokenT<GenEventInfoProduct>               genEvtInfoToken;
+  const edm::EDGetTokenT<double>  	rhoToken;
 
   std::vector<std::string> triggerPathsVector;
   std::map<std::string, int> triggerPathsMap;
@@ -332,6 +333,8 @@ class ScoutingNanoAOD : public edm::one::EDAnalyzer<edm::one::SharedResources, e
   vector<Float16_t> Vertex_ndof;
   vector<Float16_t> Vertex_isValidVtx;
 
+  float rho;
+
   // Event shape variables
   float event_isotropy;
   float event_circularity;
@@ -367,6 +370,7 @@ ScoutingNanoAOD::ScoutingNanoAOD(const edm::ParameterSet& iConfig):
   verticesToken            (consumes<std::vector<ScoutingVertex> >           (iConfig.getParameter<edm::InputTag>("vertices"))),
   pileupInfoToken          (consumes<std::vector<PileupSummaryInfo> >        (iConfig.getParameter<edm::InputTag>("pileupinfo"))),
   gensToken                (consumes<std::vector<reco::GenParticle> >        (iConfig.getParameter<edm::InputTag>("gens"))),
+  rhoToken                 (consumes<double>                                 (iConfig.getParameter<edm::InputTag>("rho"))),
   //genEvtInfoToken          (consumes<GenEventInfoProduct>                    (iConfig.getParameter<edm::InputTag>("geneventinfo"))),    
   doL1                     (iConfig.existsAs<bool>("doL1")               ?    iConfig.getParameter<bool>  ("doL1")            : false),
   doData                     (iConfig.existsAs<bool>("doData")               ?    iConfig.getParameter<bool>  ("doData")            : false),
@@ -561,6 +565,8 @@ ScoutingNanoAOD::ScoutingNanoAOD(const edm::ParameterSet& iConfig):
   tree->Branch("FatJet_mtrim"       ,&FatJet_mtrim    );
   tree->Branch("FatJet_nconst"      ,&FatJet_nconst   );
 
+  tree->Branch("rho", &rho);
+
   tree->Branch("event_isotropy"        ,&event_isotropy     );
   tree->Branch("event_circularity"     ,&event_circularity  );
   tree->Branch("event_sphericity"      ,&event_sphericity   );
@@ -1257,6 +1263,10 @@ for(int e = 0; e < static_cast<int>(truth_pts.size()); e++){//loop over pf cands
     n_pfcand_tot++;
   }
 
+  Handle<double> rhoH;
+  iEvent.getByToken(rhoToken, rhoH);
+  rho = *rhoH;
+
  // done for all events, no need to reset?
  EventShapeVariables event_algo(event_tracks);
  event_isotropy    = event_algo.isotropy();

test/ScoutingNanoAOD_cfg.py

@@ -186,6 +186,13 @@
     useLHEWeights = cms.bool(params.useWeights)
 )
 
+# get rho producer
+process.fixedGridRhoFastjetAllScouting = cms.EDProducer("FixedGridRhoProducerFastjetScouting",
+    pfCandidatesTag = cms.InputTag("hltScoutingPFPacker"),
+    electronsTag = cms.InputTag("hltScoutingEgammaPacker"),
+    maxRapidity = cms.double(5.0),
+    gridSpacing = cms.double(0.55),
+)
 
 HLTInfo = [
   "DST_DoubleMu1_noVtx_CaloScouting_v*",
@@ -216,6 +223,7 @@
   vertexinfo         = cms.InputTag("hltScoutingPFPacker","") ##Toggle this for 2016 instead of the next line,
 else:
   vertexinfo = cms.InputTag("hltScoutingPrimaryVertexPacker","primaryVtx")
+
 process.mmtree = cms.EDAnalyzer('ScoutingNanoAOD',
 	doL1 = cms.bool(False),
 	doData = cms.bool(params.data),
@@ -253,6 +261,7 @@
     gens         = cms.InputTag("genParticles"),
 	#vertices         = cms.InputTag("hltScoutingMuonPacker","displacedVtx"),
     #geneventinfo     = cms.InputTag("generator"),
+    rho          = cms.InputTag("fixedGridRhoFastjetAllScouting"),
 
     # for JEC corrections eventually
     #L1corrAK4_DATA = cms.FileInPath('CMSDIJET/DijetScoutingRootTreeMaker/data/80X_dataRun2_HLT_v12/80X_dataRun2_HLT_v12_L1FastJet_AK4CaloHLT.txt'),
@@ -261,4 +270,9 @@
 )
 #process.Tracer = cms.Service("Tracer")
 
+# add any intermediate modules to this task list
+# then unscheduled mode will call them automatically when the final module (mmtree) consumes their products
+process.myTask = cms.Task(process.fixedGridRhoFastjetAllScouting)
+
 process.p = cms.Path(                  process.mmtree)
+process.p.associate(process.myTask)