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.

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

Use the threaded version of random number generator #2403

Merged
merged 1 commit into from Feb 12, 2014
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
Expand Up @@ -24,6 +24,7 @@
*/

#include "DataFormats/Candidate/interface/Candidate.h"
#include "FWCore/Utilities/interface/StreamID.h"
#include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
#include "GeneratorInterface/PhotosInterface/interface/PhotosInterfaceBase.h"
#include "GeneratorInterface/PhotosInterface/interface/PhotosFactory.h"
Expand All @@ -38,7 +39,7 @@ class GenMuonRadiationAlgorithm
explicit GenMuonRadiationAlgorithm(const edm::ParameterSet&);
~GenMuonRadiationAlgorithm();

reco::Candidate::LorentzVector compFSR(const reco::Candidate::LorentzVector&, int, const reco::Candidate::LorentzVector&, int&);
reco::Candidate::LorentzVector compFSR(const edm::StreamID& streamID, const reco::Candidate::LorentzVector&, int, const reco::Candidate::LorentzVector&, int&);

private:
double beamEnergy_;
Expand All @@ -53,8 +54,6 @@ class GenMuonRadiationAlgorithm
static bool pythia_isInitialized_;

int verbosity_;

static CLHEP::HepRandomEngine* decayRandomEngine;
};

#endif
Expand Down
Expand Up @@ -142,8 +142,8 @@ void GenMuonRadCorrAnalyzer::analyze(const edm::Event& evt, const edm::EventSetu
double muonMinusRad = 0.;
int muonMinusRad_error = 0;
if ( muonRadiationAlgo_ ) {
muonPlusRad = muonRadiationAlgo_->compFSR(genMuonPlusP4_beforeRad, +1, genMuonMinusP4_beforeRad, muonPlusRad_error).E();
muonMinusRad = muonRadiationAlgo_->compFSR(genMuonMinusP4_beforeRad, -1, genMuonPlusP4_beforeRad, muonMinusRad_error).E();
muonPlusRad = muonRadiationAlgo_->compFSR(evt.streamID(), genMuonPlusP4_beforeRad, +1, genMuonMinusP4_beforeRad, muonPlusRad_error).E();
muonMinusRad = muonRadiationAlgo_->compFSR(evt.streamID(), genMuonMinusP4_beforeRad, -1, genMuonPlusP4_beforeRad, muonMinusRad_error).E();
} else {
muonPlusRad = genMuonPlusP4_beforeRad.E() - genMuonPlusP4_afterRad.E();
muonMinusRad = genMuonMinusP4_beforeRad.E() - genMuonMinusP4_afterRad.E();
Expand Down
2 changes: 2 additions & 0 deletions TauAnalysis/MCEmbeddingTools/plugins/MCParticleReplacer.h
Expand Up @@ -56,6 +56,8 @@ class MCParticleReplacer : public edm::EDProducer
evt_->put(product, instanceName);
}

edm::StreamID getStreamID() const { assert(evt_ != nullptr); return evt_->streamID(); }

private:
enum HepMcMode { kInvalid = 0, kNew, kReplace };
static HepMcMode stringToHepMcMode(const std::string& name);
Expand Down
34 changes: 15 additions & 19 deletions TauAnalysis/MCEmbeddingTools/plugins/ParticleReplacerZtautau.cc
Expand Up @@ -36,8 +36,6 @@ const double breitWignerWidthW = 2.141;
const double nomMassZ = 91.1876;
const double breitWignerWidthZ = 2.4952;

CLHEP::HepRandomEngine* ParticleReplacerZtautau::decayRandomEngine = nullptr;

bool ParticleReplacerZtautau::tauola_isInitialized_ = false;

typedef std::vector<reco::Particle> ParticleCollection;
Expand All @@ -53,7 +51,6 @@ ParticleReplacerZtautau::ParticleReplacerZtautau(const edm::ParameterSet& cfg)
tauola_ = (gen::TauolaInterfaceBase*)(TauolaFactory::get()->create("Tauolapp113a",cfg.getParameter<edm::ParameterSet>("TauolaOptions")));
// settings?
// usesResource(edm::SharedResourceNames::kTauola);
// you must call tauola_->setRandomEngine(decayRandomEngine); every event to properly pass the random number with the multi-threading will add below by Rnd-gen
maxNumberOfAttempts_ = ( cfg.exists("maxNumberOfAttempts") ) ?
cfg.getParameter<int>("maxNumberOfAttempts") : 10000;

Expand Down Expand Up @@ -147,17 +144,6 @@ ParticleReplacerZtautau::ParticleReplacerZtautau(const edm::ParameterSet& cfg)

rfRotationAngle_ = cfg.getParameter<double>("rfRotationAngle")*TMath::Pi()/180.;

edm::Service<edm::RandomNumberGenerator> rng;
if ( !rng.isAvailable() )
throw cms::Exception("Configuration")
<< "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
<< "which appears to be absent. Please add that service to your configuration\n"
<< "or remove the modules that require it.\n";

// this is a global variable defined in GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc
decayRandomEngine = &rng->getEngine();
tauola_->setRandomEngine(decayRandomEngine); // you must call tauola_->setRandomEngine(decayRandomEngine); every event to properly pass the random number with the multi-threading

std::string applyMuonRadiationCorrection_string = cfg.getParameter<std::string>("applyMuonRadiationCorrection");
if ( applyMuonRadiationCorrection_string != "" ) {
edm::ParameterSet cfgMuonRadiationAlgo;
Expand Down Expand Up @@ -246,6 +232,16 @@ namespace

std::auto_ptr<HepMC::GenEvent> ParticleReplacerZtautau::produce(const std::vector<reco::Particle>& muons, const reco::Vertex* evtVtx, const HepMC::GenEvent* genEvt, MCParticleReplacer* producer)
{
edm::Service<edm::RandomNumberGenerator> rng;
if ( !rng.isAvailable() )
throw cms::Exception("Configuration")
<< "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
<< "which appears to be absent. Please add that service to your configuration\n"
<< "or remove the modules that require it.\n";

CLHEP::HepRandomEngine& decayRandomEngine = rng->getEngine(producer->getStreamID());
tauola_->setRandomEngine(&decayRandomEngine);

if ( evtVtx != 0 ) throw cms::Exception("Configuration")
<< "ParticleReplacerZtautau does NOT support using primary vertex as the origin for taus !!\n";

Expand Down Expand Up @@ -286,7 +282,7 @@ std::auto_ptr<HepMC::GenEvent> ParticleReplacerZtautau::produce(const std::vecto
reco::Particle embedLepton1(muon1.charge(), muon1.p4(), muon1.vertex(), muon1.pdgId(), 0, true);
reco::Particle embedLepton2(muon2.charge(), muon2.p4(), muon2.vertex(), muon2.pdgId(), 0, true);
if ( targetParticle1AbsPdgID_ != 13 ) {
transformMuMu2LepLep(&embedLepton1, &embedLepton2);
transformMuMu2LepLep(decayRandomEngine, &embedLepton1, &embedLepton2);
}
embedParticles.push_back(embedLepton1);
embedParticles.push_back(embedLepton2);
Expand Down Expand Up @@ -403,7 +399,7 @@ std::auto_ptr<HepMC::GenEvent> ParticleReplacerZtautau::produce(const std::vecto
if ( deltaR(otherParticle->p4(), embedParticle->p4()) > 1.e-3 ) sumOtherParticlesP4 += otherParticle->p4();
}
int errorFlag = 0;
reco::Candidate::LorentzVector muonFSR = muonRadiationAlgo_->compFSR(embedParticle->p4(), embedParticle->charge(), sumOtherParticlesP4, errorFlag);
reco::Candidate::LorentzVector muonFSR = muonRadiationAlgo_->compFSR(producer->getStreamID(), embedParticle->p4(), embedParticle->charge(), sumOtherParticlesP4, errorFlag);
double embedParticlePx_corrected = embedParticle->px() + muonFSR.px();
double embedParticlePy_corrected = embedParticle->py() + muonFSR.py();
double embedParticlePz_corrected = embedParticle->pz() + muonFSR.pz();
Expand Down Expand Up @@ -481,7 +477,7 @@ std::auto_ptr<HepMC::GenEvent> ParticleReplacerZtautau::produce(const std::vecto
// << " Py = " << embedParticle->py() << ","
// << " Pz = " << embedParticle->pz() << std::endl;
int errorFlag = 0;
reco::Candidate::LorentzVector muonFSR = muonRadiationAlgo_->compFSR(embedParticle->p4(), embedParticle->charge(), sumOtherParticlesP4, errorFlag);
reco::Candidate::LorentzVector muonFSR = muonRadiationAlgo_->compFSR(producer->getStreamID(), embedParticle->p4(), embedParticle->charge(), sumOtherParticlesP4, errorFlag);
//std::cout << "muonFSR:"
// << " En = " << muonFSR.E() << ","
// << " Px = " << muonFSR.px() << ","
Expand Down Expand Up @@ -888,7 +884,7 @@ namespace
}
}

void ParticleReplacerZtautau::transformMuMu2LepLep(reco::Particle* muon1, reco::Particle* muon2)
void ParticleReplacerZtautau::transformMuMu2LepLep(CLHEP::HepRandomEngine& randomEngine, reco::Particle* muon1, reco::Particle* muon2)
{
//--- transform a muon pair into an electron/tau pair,
// taking into account the difference between muon and electron/tau mass
Expand Down Expand Up @@ -918,7 +914,7 @@ void ParticleReplacerZtautau::transformMuMu2LepLep(reco::Particle* muon1, reco::
if ( rfRotationAngle_ != 0. ) {
double rfRotationAngle_value = rfRotationAngle_;
if ( rfRotationAngle_ == -1. ) {
double u = decayRandomEngine->flat();
double u = randomEngine.flat();
rfRotationAngle_value = 2.*TMath::Pi()*u;
}

Expand Down
Expand Up @@ -51,7 +51,7 @@ class ParticleReplacerZtautau : public ParticleReplacerBase
virtual void endJob();

private:
void transformMuMu2LepLep(reco::Particle*, reco::Particle*);
void transformMuMu2LepLep(CLHEP::HepRandomEngine& randomEngine, reco::Particle*, reco::Particle*);
void transformMuMu2TauNu(reco::Particle*, reco::Particle*);

HepMC::GenEvent* processEventWithTauola(HepMC::GenEvent*);
Expand Down Expand Up @@ -133,8 +133,6 @@ class ParticleReplacerZtautau : public ParticleReplacerBase

int maxNumberOfAttempts_;

static CLHEP::HepRandomEngine* decayRandomEngine;

};

#endif
Expand Down
28 changes: 14 additions & 14 deletions TauAnalysis/MCEmbeddingTools/src/GenMuonRadiationAlgorithm.cc
Expand Up @@ -28,7 +28,6 @@ const double protonMass = 0.938272;

bool GenMuonRadiationAlgorithm::photos_isInitialized_ = false;
bool GenMuonRadiationAlgorithm::pythia_isInitialized_ = false;
CLHEP::HepRandomEngine* GenMuonRadiationAlgorithm::decayRandomEngine = nullptr;

class myPythia6ServiceWithCallback : public gen::Pythia6Service
{
Expand Down Expand Up @@ -166,16 +165,6 @@ GenMuonRadiationAlgorithm::GenMuonRadiationAlgorithm(const edm::ParameterSet& cf
photos_(0),
pythia_(0)
{
edm::Service<edm::RandomNumberGenerator> rng;
if ( !rng.isAvailable() )
throw cms::Exception("Configuration")
<< "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
<< "which appears to be absent. Please add that service to your configuration\n"
<< "or remove the modules that require it.\n";

// this is a global variable defined in GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc
decayRandomEngine = &rng->getEngine();

std::string mode_string = cfg.getParameter<std::string>("mode");
if ( mode_string == "pythia" ) mode_ = kPYTHIA;
else if ( mode_string == "photos" ) mode_ = kPHOTOS;
Expand All @@ -187,7 +176,6 @@ GenMuonRadiationAlgorithm::GenMuonRadiationAlgorithm(const edm::ParameterSet& cf
photos_ = (gen::PhotosInterfaceBase*)(PhotosFactory::get()->create("Photos2155", cfg.getParameter<edm::ParameterSet>("PhotosOptions")));
//settings?
//usesResource(edm::SharedResourceNames::kPhotos);
//you must call photos_->setRandomEngine(decayRandomEngine); every event to properly pass the random number with the multi-threading will add below by Rnd-gen
}

verbosity_ = ( cfg.exists("verbosity") ) ?
Expand Down Expand Up @@ -236,9 +224,21 @@ namespace
}
}

reco::Candidate::LorentzVector GenMuonRadiationAlgorithm::compFSR(const reco::Candidate::LorentzVector& muonP4, int muonCharge,
const reco::Candidate::LorentzVector& otherP4, int& errorFlag)
reco::Candidate::LorentzVector GenMuonRadiationAlgorithm::compFSR(const edm::StreamID& streamID,
const reco::Candidate::LorentzVector& muonP4, int muonCharge,
const reco::Candidate::LorentzVector& otherP4, int& errorFlag)
{
edm::Service<edm::RandomNumberGenerator> rng;
if ( !rng.isAvailable() )
throw cms::Exception("Configuration")
<< "The GenMuonRadiationAlgorithm module requires the RandomNumberGeneratorService\n"
<< "which appears to be absent. Please add that service to your configuration\n"
<< "or remove the modules that require it.\n";

// this is a global variable defined in GeneratorInterface/ExternalDecays/src/ExternalDecayDriver.cc
CLHEP::HepRandomEngine& decayRandomEngine = rng->getEngine(streamID);
photos_->setRandomEngine(&decayRandomEngine);

if ( verbosity_ ) {
std::cout << "<GenMuonRadiationAlgorithm::compMuonFSR>:" << std::endl;
std::cout << " muon: En = " << muonP4.E() << ", Pt = " << muonP4.pt() << ", eta = " << muonP4.eta() << ", phi = " << muonP4.phi() << ", charge = " << muonCharge << std::endl;
Expand Down