Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-02 08:36:15

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/GlobalSystemOfUnits.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 }