CMS 3D CMS Logo

RPCSimTriv.cc

Go to the documentation of this file.
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       // Here I hould check if the RPC are up side down;
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 }

Generated on Tue Jun 9 17:47:50 2009 for CMSSW by  doxygen 1.5.4