00001 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
00002 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
00003 #include "SimMuon/RPCDigitizer/src/RPCSimTriv.h"
00004 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
00005 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
00006 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
00007 #include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h"
00008
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012 #include "CLHEP/Random/RandomEngine.h"
00013 #include "CLHEP/Random/RandFlat.h"
00014 #include <cmath>
00015 #include <utility>
00016 #include <map>
00017
00018 #include "CLHEP/config/CLHEP.h"
00019 #include "CLHEP/Random/Random.h"
00020 #include "CLHEP/Random/RandFlat.h"
00021 #include "CLHEP/Random/RandPoissonQ.h"
00022
00023
00024 RPCSimTriv::RPCSimTriv(const edm::ParameterSet& config) : RPCSim(config){
00025
00026 rate=config.getParameter<double>("Rate");
00027 nbxing=config.getParameter<int>("Nbxing");
00028 gate=config.getParameter<double>("Gate");
00029
00030 edm::Service<edm::RandomNumberGenerator> rng;
00031 if ( ! rng.isAvailable()) {
00032 throw cms::Exception("Configuration")
00033 << "RPCDigitizer requires the RandomNumberGeneratorService\n"
00034 "which is not present in the configuration file. You must add the service\n"
00035 "in the configuration file or remove the modules that require it.";
00036 }
00037 _rpcSync = new RPCSynchronizer(config);
00038
00039 rndEngine = &(rng->getEngine());
00040 }
00041
00042 RPCSimTriv::~RPCSimTriv(){}
00043
00044 void
00045 RPCSimTriv::simulate(const RPCRoll* roll,
00046 const edm::PSimHitContainer& rpcHits)
00047 {
00048
00049 _rpcSync->setRPCSimSetUp(getRPCSimSetUp());
00050 theRpcDigiSimLinks.clear();
00051 theDetectorHitMap.clear();
00052 theRpcDigiSimLinks = RPCDigiSimLinks(roll->id().rawId());
00053
00054 const Topology& topology=roll->specs()->topology();
00055 for (edm::PSimHitContainer::const_iterator _hit = rpcHits.begin();
00056 _hit != rpcHits.end(); ++_hit){
00057
00058 int type = _hit->particleType();
00059 if (type == 13 || type == -13){
00060
00061 const LocalPoint& entr=_hit->entryPoint();
00062 int time_hit = _rpcSync->getSimHitBx(&(*_hit));
00063 std::pair<int, int> digi(topology.channel(entr)+1,
00064 time_hit);
00065
00066 theDetectorHitMap.insert(DetectorHitMap::value_type(digi,&(*_hit)));
00067 strips.insert(digi);
00068 }
00069 }
00070 }
00071
00072
00073 void RPCSimTriv::simulateNoise(const RPCRoll* roll)
00074 {
00075
00076 RPCDetId rpcId = roll->id();
00077 int nstrips = roll->nstrips();
00078 double area = 0.0;
00079
00080 if ( rpcId.region() == 0 )
00081 {
00082 const RectangularStripTopology* top_ = dynamic_cast<const
00083 RectangularStripTopology*>(&(roll->topology()));
00084 float xmin = (top_->localPosition(0.)).x();
00085 float xmax = (top_->localPosition((float)roll->nstrips())).x();
00086 float striplength = (top_->stripLength());
00087 area = striplength*(xmax-xmin);
00088 }
00089 else
00090 {
00091 const TrapezoidalStripTopology* top_=dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology()));
00092 float xmin = (top_->localPosition(0.)).x();
00093 float xmax = (top_->localPosition((float)roll->nstrips())).x();
00094 float striplength = (top_->stripLength());
00095 area = striplength*(xmax-xmin);
00096 }
00097
00098 double ave = rate*nbxing*gate*area*1.0e-9;
00099 poissonDistribution_ = new CLHEP::RandPoissonQ(rndEngine, ave);
00100 N_hits = poissonDistribution_->fire();
00101
00102 for (int i = 0; i < N_hits; i++ ){
00103
00104 flatDistribution = new CLHEP::RandFlat(rndEngine, 1, nstrips);
00105 int strip = static_cast<int>(flatDistribution->fire());
00106 int time_hit;
00107
00108 flatDistribution = new CLHEP::RandFlat(rndEngine, (nbxing*gate)/gate);
00109 time_hit = (static_cast<int>(flatDistribution->fire())) - nbxing/2;
00110
00111 std::pair<int, int> digi(strip,time_hit);
00112 strips.insert(digi);
00113 }
00114
00115 }