Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

use deterministic seed for time smearing producers #21354

Merged
merged 2 commits into from Nov 25, 2017
Merged

use deterministic seed for time smearing producers #21354

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
99d612c
use deterministic seed for TrackTimeValueMapProducer
kpedro88 Nov 17, 2017
afb40b3
use deterministic seed for EcalBarrelClusterFastTimer
kpedro88 Nov 22, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
25 changes: 10 additions & 15 deletions RecoParticleFlow/PFSimProducer/plugins/EcalBarrelClusterFastTimer.cc
View file
Open in desktop
Expand Up @@ -25,14 +25,12 @@
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.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"

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

void EcalBarrelClusterFastTimer::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));

edm::Handle<std::vector<reco::PFCluster> > clustersH;
edm::Handle<EcalRecHitCollection> timehitsH;

Expand All @@ -116,6 +103,13 @@ void EcalBarrelClusterFastTimer::produce(edm::StreamID sid, edm::Event& evt, con
const auto& clusters = *clustersH;
const auto& timehits = *timehitsH;

// 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());
std::uint32_t seed = (lumiNum_uint<<10) + (runNum_uint<<20) + evNum_uint;
std::mt19937 rng(seed);

std::vector<std::pair<float,DetId> > times; // perfect times keyed to cluster index
times.reserve(clusters.size());

Expand All @@ -133,8 +127,9 @@ void EcalBarrelClusterFastTimer::produce(edm::StreamID sid, edm::Event& evt, con
// smear once then correct to multiple vertices
for( unsigned i = 0 ; i < clusters.size(); ++i ) {
const float theresolution = reso->getTimeResolution(clusters[i]);
std::normal_distribution<float> gausTime(times[i].first, theresolution);

smeared_times.emplace_back( CLHEP::RandGauss::shoot(rng_engine, times[i].first, theresolution) );
smeared_times.emplace_back( gausTime(rng) );
resolutions.push_back( theresolution );
}

Expand Down
33 changes: 13 additions & 20 deletions SimTracker/TrackAssociation/plugins/TrackTimeValueMapProducer.cc
View file
Open in desktop
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