Merge pull request #20878 from gpetruc/egsmearer_semideterministic

introduce semi-deterministic smearing for e/gamma objects
cms-sw · Oct 12, 2017 · b652b5c · b652b5c
2 parents 561fe85 + 7a44c12
commit b652b5c
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 0 deletions.
diff --git a/EgammaAnalysis/ElectronTools/plugins/CalibratedElectronProducersRun2.cc b/EgammaAnalysis/ElectronTools/plugins/CalibratedElectronProducersRun2.cc
@@ -17,6 +17,8 @@
 #include "EgammaAnalysis/ElectronTools/interface/ElectronEnergyCalibratorRun2.h"
 
 #include <vector>
+#include <random>
+#include <TRandom2.h>
 
 template<typename T>
 class CalibratedElectronProducerRun2T: public edm::stream::EDProducer<>
@@ -33,6 +35,8 @@ class CalibratedElectronProducerRun2T: public edm::stream::EDProducer<>
 
         EpCombinationTool theEpCombinationTool;
         ElectronEnergyCalibratorRun2 theEnCorrectorRun2;
+        std::unique_ptr<TRandom> theSemiDeterministicRng;
+
 };
 
 template<typename T>
@@ -42,6 +46,10 @@ CalibratedElectronProducerRun2T<T>::CalibratedElectronProducerRun2T( const edm::
   theEpCombinationTool(),
   theEnCorrectorRun2(theEpCombinationTool, conf.getParameter<bool>("isMC"), conf.getParameter<bool>("isSynchronization"), conf.getParameter<std::string>("correctionFile"))
 {
+  if (conf.existsAs<bool>("semiDeterministic") && conf.getParameter<bool>("semiDeterministic")) {
+    theSemiDeterministicRng.reset(new TRandom2());
+    theEnCorrectorRun2.initPrivateRng(theSemiDeterministicRng.get());
+  }
   produces<std::vector<T> >();
 }
 
@@ -63,6 +71,17 @@ CalibratedElectronProducerRun2T<T>::produce( edm::Event & iEvent, const edm::Eve
     std::unique_ptr<std::vector<T> > out(new std::vector<T>());
     out->reserve(in->size());   
 
+    if (theSemiDeterministicRng && !in->empty()) { // no need to set a seed if in is empty
+        const auto & first = in->front();
+        std::seed_seq seeder = {int(iEvent.id().event()), int(iEvent.id().luminosityBlock()), int(iEvent.id().run()),
+          int(in->size()), int(std::numeric_limits<int>::max()*first.phi()/M_PI) & 0xFFF, int(first.pdgId())};
+        uint32_t seed = 0, tries = 10;
+        do {
+            seeder.generate(&seed,&seed+1); tries++;
+        } while (seed == 0 && tries < 10);
+        theSemiDeterministicRng->SetSeed(seed ? seed : iEvent.id().event());
+    }
+
     for (const T &ele : *in) {
         out->push_back(ele);
         theEnCorrectorRun2.calibrate(out->back(), iEvent.id().run(), iEvent.streamID());

diff --git a/EgammaAnalysis/ElectronTools/plugins/CalibratedPhotonProducersRun2.cc b/EgammaAnalysis/ElectronTools/plugins/CalibratedPhotonProducersRun2.cc
@@ -13,6 +13,8 @@
 #include "EgammaAnalysis/ElectronTools/interface/PhotonEnergyCalibratorRun2.h"
 
 #include <vector>
+#include <random>
+#include <TRandom2.h>
 
 template<typename T>
 class CalibratedPhotonProducerRun2T: public edm::stream::EDProducer<> {
@@ -24,13 +26,18 @@ class CalibratedPhotonProducerRun2T: public edm::stream::EDProducer<> {
 private:
   edm::EDGetTokenT<edm::View<T> > thePhotonToken;
   PhotonEnergyCalibratorRun2 theEnCorrectorRun2;
+  std::unique_ptr<TRandom> theSemiDeterministicRng;
 };
 
 template<typename T>
 CalibratedPhotonProducerRun2T<T>::CalibratedPhotonProducerRun2T( const edm::ParameterSet & conf ) :
   thePhotonToken(consumes<edm::View<T> >(conf.getParameter<edm::InputTag>("photons"))),
   theEnCorrectorRun2(conf.getParameter<bool>("isMC"), conf.getParameter<bool>("isSynchronization"), conf.getParameter<std::string >("correctionFile")) {
 
+  if (conf.existsAs<bool>("semiDeterministic") && conf.getParameter<bool>("semiDeterministic")) {
+    theSemiDeterministicRng.reset(new TRandom2());
+    theEnCorrectorRun2.initPrivateRng(theSemiDeterministicRng.get());
+  }
   produces<std::vector<T> >();
 }
 
@@ -45,6 +52,17 @@ CalibratedPhotonProducerRun2T<T>::produce( edm::Event & iEvent, const edm::Event
   edm::Handle<edm::View<T> > in;
   iEvent.getByToken(thePhotonToken, in);
 
+  if (theSemiDeterministicRng && !in->empty()) { // no need to set a seed if in is empty
+      const auto & first = in->front();
+      std::seed_seq seeder = {int(iEvent.id().event()), int(iEvent.id().luminosityBlock()), int(iEvent.id().run()),
+          int(in->size()), int(std::numeric_limits<int>::max()*first.phi()/M_PI) & 0xFFF, int(first.pdgId())};
+      uint32_t seed = 0, tries = 10;
+      do {
+          seeder.generate(&seed,&seed+1); tries++;
+      } while (seed == 0 && tries < 10);
+      theSemiDeterministicRng->SetSeed(seed ? seed : iEvent.id().event());
+  }
+
   std::unique_ptr<std::vector<T> > out(new std::vector<T>());
   out->reserve(in->size());