Skip to content

Commit

Permalink
use deterministic seed for TrackTimeValueMapProducer
Browse files Browse the repository at this point in the history
  • Loading branch information
kpedro88 committed Nov 17, 2017
1 parent 6a9e6e4 commit d68a901
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 28 deletions.
8 changes: 0 additions & 8 deletions IOMC/RandomEngine/python/IOMC_cff.py
Expand Up @@ -205,14 +205,6 @@
from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing
phase2_timing.toModify(
RandomNumberGeneratorService,
trackTimeValueMapProducer = cms.PSet(
initialSeed = cms.untracked.uint32(1234567),
engineName = cms.untracked.string('HepJamesRandom')
),
gsfTrackTimeValueMapProducer = cms.PSet(
initialSeed = cms.untracked.uint32(1234567),
engineName = cms.untracked.string('HepJamesRandom')
),
ecalBarrelClusterFastTimer = cms.PSet(
initialSeed = cms.untracked.uint32(1234567),
engineName = cms.untracked.string('HepJamesRandom')
Expand Down
33 changes: 13 additions & 20 deletions SimTracker/TrackAssociation/plugins/TrackTimeValueMapProducer.cc
Expand Up @@ -30,16 +30,11 @@
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"


#include <random>
#include <memory>

#include "SimTracker/TrackAssociation/interface/ResolutionModel.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "CLHEP/Random/RandGauss.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"

#include "FWCore/Utilities/interface/isFinite.h"
#include "FWCore/Utilities/interface/transform.h"

Expand Down Expand Up @@ -110,21 +105,9 @@ TrackTimeValueMapProducer::TrackTimeValueMapProducer(const edm::ParameterSet& co
produces<edm::ValueMap<float> >(tracksName_+name);
produces<edm::ValueMap<float> >(tracksName_+name+resolution);
}
// get RNG engine
edm::Service<edm::RandomNumberGenerator> rng;
if (!rng.isAvailable()){
throw cms::Exception("Configuration")
<< "TrackTimeValueMapProducer::TrackTimeValueMapProducer() - RandomNumberGeneratorService is not present in configuration file.\n"
<< "Add the service in the configuration file or remove the modules that require it.";
}

}

void TrackTimeValueMapProducer::produce(edm::StreamID sid, edm::Event& evt, const edm::EventSetup& es) const {
// get RNG engine
edm::Service<edm::RandomNumberGenerator> rng;
auto rng_engine = &(rng->getEngine(sid));

// get sim track associators
std::vector<edm::Handle<reco::TrackToTrackingParticleAssociator> > associators;
for( const auto& token : associators_ ) {
Expand Down Expand Up @@ -176,6 +159,15 @@ void TrackTimeValueMapProducer::produce(edm::StreamID sid, edm::Event& evt, cons
double meanSimTime = sumSimTime/double(nsim);
double varSimTime = sumSimTimeSq/double(nsim) - meanSimTime*meanSimTime;
double rmsSimTime = std::sqrt(std::max(0.1*0.1,varSimTime));
std::normal_distribution<float> gausSimTime(meanSimTime, rmsSimTime);

// get event-based seed for RNG
unsigned int runNum_uint = static_cast <unsigned int> (evt.id().run());
unsigned int lumiNum_uint = static_cast <unsigned int> (evt.id().luminosityBlock());
unsigned int evNum_uint = static_cast <unsigned int> (evt.id().event());
unsigned int tkChi2_uint = uint32_t(TrackCollection.empty() ? 0 : TrackCollection.refAt(0)->chi2()/0.01);
std::uint32_t seed = tkChi2_uint + (lumiNum_uint<<10) + (runNum_uint<<20) + evNum_uint;
std::mt19937 rng(seed);

for( unsigned itk = 0; itk < TrackCollection.size(); ++itk ) {
const auto tkref = TrackCollection.refAt(itk);
Expand All @@ -192,7 +184,7 @@ void TrackTimeValueMapProducer::produce(edm::StreamID sid, edm::Event& evt, cons
generalTrackTimes.push_back(time);
}
else {
float rndtime = CLHEP::RandGauss::shoot(rng_engine, meanSimTime, rmsSimTime);
float rndtime = gausSimTime(rng);
generalTrackTimes.push_back(rndtime);
if (track_tps != associatedTracks.back().end() && track_tps->val.size() > 1) {
LogDebug("TooManyTracks") << "track matched to " << track_tps->val.size() << " tracking particles!" << std::endl;
Expand All @@ -213,7 +205,8 @@ void TrackTimeValueMapProducer::produce(edm::StreamID sid, edm::Event& evt, cons
bool inAcceptance = absEta < etaMax_ && absEta >= etaMin_ && tk.p()>pMin_ && (absEta>etaMaxForPtThreshold_ || tk.pt()>ptMin_);
if (inAcceptance) {
const float resolution = reso->getTimeResolution(tk);
times.push_back( CLHEP::RandGauss::shoot(rng_engine, generalTrackTimes[i], resolution) );
std::normal_distribution<float> gausGeneralTime(generalTrackTimes[i], resolution);
times.push_back( gausGeneralTime(rng) );
resos.push_back( resolution );
}
else {
Expand Down

0 comments on commit d68a901

Please sign in to comment.