File indexing completed on 2024-05-10 02:20:57
0001 #include "IOMC/EventVertexGenerators/interface/GaussianZBeamSpotFilter.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008
0009 #include "HepMC/GenRanges.h"
0010 #include <CLHEP/Units/SystemOfUnits.h>
0011 #include <CLHEP/Random/RandomEngine.h>
0012
0013 GaussianZBeamSpotFilter::GaussianZBeamSpotFilter(const edm::ParameterSet& iPSet)
0014 : src_(iPSet.getParameter<edm::InputTag>("src")),
0015 baseSZ_(iPSet.getParameter<double>("baseSZ") * CLHEP::cm),
0016 baseZ0_(iPSet.getParameter<double>("baseZ0") * CLHEP::cm),
0017 newSZ_(iPSet.getParameter<double>("newSZ") * CLHEP::cm),
0018 newZ0_(iPSet.getParameter<double>("newZ0") * CLHEP::cm),
0019 srcToken_(consumes<edm::HepMCProduct>(src_)) {
0020 edm::Service<edm::RandomNumberGenerator> rng;
0021 if (!rng.isAvailable()) {
0022 throw cms::Exception("Configuration")
0023 << "The GaussianZBeamSpotFilter requires the RandomNumberGeneratorService\n"
0024 "which is not present in the configuration file. You must add the service\n"
0025 "in the configuration file or remove the modules that require it.";
0026 }
0027
0028 if (baseZ0_ != newZ0_) {
0029 edm::LogError("InconsistentZPosition") << "Z0 : old " << baseZ0_ << " and new " << newZ0_ << " are different";
0030 }
0031
0032 if (baseSZ_ < newSZ_) {
0033 edm::LogError("InconsistentZSigma") << "Sigma Z : old " << baseSZ_ << " smaller than new " << newSZ_;
0034 }
0035 }
0036
0037 bool GaussianZBeamSpotFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0038 bool pass = true;
0039
0040 edm::Service<edm::RandomNumberGenerator> rng;
0041 CLHEP::HepRandomEngine& engine = rng->getEngine(iEvent.streamID());
0042
0043 const edm::Handle<edm::HepMCProduct>& HepMCEvt = iEvent.getHandle(srcToken_);
0044
0045 HepMC::GenEvent::vertex_const_iterator vitr = HepMCEvt->GetEvent()->vertex_range().begin();
0046
0047 if (vitr != HepMCEvt->GetEvent()->vertex_range().end()) {
0048 double vtxZ = (*vitr)->point3d().z();
0049 double gaussRatio = std::exp(-(vtxZ - newZ0_) * (vtxZ - newZ0_) / (2.0 * newSZ_ * newSZ_) +
0050 (vtxZ - baseZ0_) * (vtxZ - baseZ0_) / (2.0 * baseSZ_ * baseSZ_));
0051 if (engine.flat() > gaussRatio) {
0052 pass = false;
0053 }
0054 edm::LogVerbatim("GaussianZBeam") << "base sigmaZ = " << baseSZ_ << " new sigmaZ = " << newSZ_ << " vtxZ = " << vtxZ
0055 << " gaussian ratio = " << gaussRatio << " pass = " << pass;
0056 }
0057
0058 return pass;
0059 }
