RPCDigitizer/src/RPCSimAverageNoise.cc

0001 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0002 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
0003 #include "SimMuon/RPCDigitizer/src/RPCSimAverageNoise.h"
0004
0005 #include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h"
0006 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0007 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0008
0009 #include <cmath>
0010
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0017 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0018 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0019 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0020 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0021 #include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h"
0022
0023 #include <cstring>
0024 #include <iostream>
0025 #include <fstream>
0026 #include <string>
0027 #include <vector>
0028 #include <cstdlib>
0029 #include <utility>
0030 #include <map>
0031
0032 #include "CLHEP/Random/RandFlat.h"
0033 #include "CLHEP/Random/RandPoissonQ.h"
0034
0035 using namespace std;
0036
0037 RPCSimAverageNoise::RPCSimAverageNoise(const edm::ParameterSet& config) : RPCSim(config) {
0038   aveEff = config.getParameter<double>("averageEfficiency");
0039   aveCls = config.getParameter<double>("averageClusterSize");
0040   resRPC = config.getParameter<double>("timeResolution");
0041   timOff = config.getParameter<double>("timingRPCOffset");
0042   dtimCs = config.getParameter<double>("deltatimeAdjacentStrip");
0043   resEle = config.getParameter<double>("timeJitter");
0044   sspeed = config.getParameter<double>("signalPropagationSpeed");
0045   lbGate = config.getParameter<double>("linkGateWidth");
0046   rpcdigiprint = config.getParameter<bool>("printOutDigitizer");
0047   rate = config.getParameter<double>("Rate");
0048   nbxing = config.getParameter<int>("Nbxing");
0049   gate = config.getParameter<double>("Gate");
0050   frate = config.getParameter<double>("Frate");
0051
0052   if (rpcdigiprint) {
0053     std::cout << "Average Efficiency        = " << aveEff << std::endl;
0054     std::cout << "Average Cluster Size      = " << aveCls << " strips" << std::endl;
0055     std::cout << "RPC Time Resolution       = " << resRPC << " ns" << std::endl;
0056     std::cout << "RPC Signal formation time = " << timOff << " ns" << std::endl;
0057     std::cout << "RPC adjacent strip delay  = " << dtimCs << " ns" << std::endl;
0058     std::cout << "Electronic Jitter         = " << resEle << " ns" << std::endl;
0059     std::cout << "Signal propagation time   = " << sspeed << " x c" << std::endl;
0060     std::cout << "Link Board Gate Width     = " << lbGate << " ns" << std::endl;
0061   }
0062
0063   _rpcSync = new RPCSynchronizer(config);
0064 }
0065
0066 RPCSimAverageNoise::~RPCSimAverageNoise() { delete _rpcSync; }
0067
0068 int RPCSimAverageNoise::getClSize(float posX, CLHEP::HepRandomEngine* engine) {
0069   std::map<int, std::vector<double> > clsMap = getRPCSimSetUp()->getClsMap();
0070
0071   int cnt = 1;
0072   int min = 1;
0073   double func = 0.0;
0074   std::vector<double> sum_clsize;
0075
0076   double rr_cl = CLHEP::RandFlat::shoot(engine);
0077   if (0.0 <= posX && posX < 0.2) {
0078     func = (clsMap[1])[(clsMap[1]).size() - 1] * (rr_cl);
0079     sum_clsize = clsMap[1];
0080   }
0081   if (0.2 <= posX && posX < 0.4) {
0082     func = (clsMap[2])[(clsMap[2]).size() - 1] * (rr_cl);
0083     sum_clsize = clsMap[2];
0084   }
0085   if (0.4 <= posX && posX < 0.6) {
0086     func = (clsMap[3])[(clsMap[3]).size() - 1] * (rr_cl);
0087     sum_clsize = clsMap[3];
0088   }
0089   if (0.6 <= posX && posX < 0.8) {
0090     func = (clsMap[4])[(clsMap[4]).size() - 1] * (rr_cl);
0091     sum_clsize = clsMap[4];
0092   }
0093   if (0.8 <= posX && posX < 1.0) {
0094     func = (clsMap[5])[(clsMap[5]).size() - 1] * (rr_cl);
0095     sum_clsize = clsMap[5];
0096   }
0097
0098   for (vector<double>::iterator iter = sum_clsize.begin(); iter != sum_clsize.end(); ++iter) {
0099     cnt++;
0100     if (func > (*iter)) {
0101       min = cnt;
0102     } else if (func < (*iter)) {
0103       break;
0104     }
0105   }
0106   return min;
0107 }
0108
0109 void RPCSimAverageNoise::simulate(const RPCRoll* roll,
0110                                   const edm::PSimHitContainer& rpcHits,
0111                                   CLHEP::HepRandomEngine* engine) {
0112   _rpcSync->setRPCSimSetUp(getRPCSimSetUp());
0113   theRpcDigiSimLinks.clear();
0114   theDetectorHitMap.clear();
0115   theRpcDigiSimLinks = RPCDigiSimLinks(roll->id().rawId());
0116
0117   const Topology& topology = roll->specs()->topology();
0118
0119   for (edm::PSimHitContainer::const_iterator _hit = rpcHits.begin(); _hit != rpcHits.end(); ++_hit) {
0120     // Here I hould check if the RPC are up side down;
0121     const LocalPoint& entr = _hit->entryPoint();
0122     int time_hit = _rpcSync->getSimHitBx(&(*_hit), engine);
0123     float posX = roll->strip(_hit->localPosition()) - static_cast<int>(roll->strip(_hit->localPosition()));
0124
0125     // Effinciecy
0126
0127     if (CLHEP::RandFlat::shoot(engine) < aveEff) {
0128       int centralStrip = topology.channel(entr) + 1;
0129       int fstrip = centralStrip;
0130       int lstrip = centralStrip;
0131       // Compute the cluster size
0132       //double w = CLHEP::RandFlat::shoot(engine);
0133       //if (w < 1.e-10) w=1.e-10;
0134       int clsize = this->getClSize(posX, engine);
0135
0136       std::vector<int> cls;
0137       cls.push_back(centralStrip);
0138       if (clsize > 1) {
0139         for (int cl = 0; cl < (clsize - 1) / 2; cl++) {
0140           if (centralStrip - cl - 1 >= 1) {
0141             fstrip = centralStrip - cl - 1;
0142             cls.push_back(fstrip);
0143           }
0144           if (centralStrip + cl + 1 <= roll->nstrips()) {
0145             lstrip = centralStrip + cl + 1;
0146             cls.push_back(lstrip);
0147           }
0148         }
0149         if (clsize % 2 == 0) {
0150           // insert the last strip according to the
0151           // simhit position in the central strip
0152           double deltaw = roll->centreOfStrip(centralStrip).x() - entr.x();
0153           if (deltaw < 0.) {
0154             if (lstrip < roll->nstrips()) {
0155               lstrip++;
0156               cls.push_back(lstrip);
0157             }
0158           } else {
0159             if (fstrip > 1) {
0160               fstrip--;
0161               cls.push_back(fstrip);
0162             }
0163           }
0164         }
0165       }
0166
0167       for (std::vector<int>::iterator i = cls.begin(); i != cls.end(); i++) {
0168         // Check the timing of the adjacent strip
0169         std::pair<int, int> digi(*i, time_hit);
0170
0171         theDetectorHitMap.insert(DetectorHitMap::value_type(digi, &(*_hit)));
0172         strips.insert(digi);
0173       }
0174     }
0175   }
0176 }
0177
0178 void RPCSimAverageNoise::simulateNoise(const RPCRoll* roll, CLHEP::HepRandomEngine* engine) {
0179   RPCDetId rpcId = roll->id();
0180   std::vector<float> vnoise = (getRPCSimSetUp())->getNoise(rpcId.rawId());
0181   unsigned int nstrips = roll->nstrips();
0182
0183   double area = 0.0;
0184
0185   if (rpcId.region() == 0) {
0186     const RectangularStripTopology* top_ = dynamic_cast<const RectangularStripTopology*>(&(roll->topology()));
0187     float xmin = (top_->localPosition(0.)).x();
0188     float xmax = (top_->localPosition((float)roll->nstrips())).x();
0189     float striplength = (top_->stripLength());
0190     area = striplength * (xmax - xmin);
0191   } else {
0192     const TrapezoidalStripTopology* top_ = dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology()));
0193     float xmin = (top_->localPosition(0.)).x();
0194     float xmax = (top_->localPosition((float)roll->nstrips())).x();
0195     float striplength = (top_->stripLength());
0196     area = striplength * (xmax - xmin);
0197   }
0198   for (unsigned int j = 0; j < vnoise.size(); ++j) {
0199     if (j >= nstrips)
0200       break;
0201
0202     double ave = frate * vnoise[j] * nbxing * gate * area * 1.0e-9;
0203     CLHEP::RandPoissonQ randPoissonQ(*engine, ave);
0204     N_hits = randPoissonQ.fire();
0205
0206     for (int i = 0; i < N_hits; i++) {
0207       int time_hit = (static_cast<int>(CLHEP::RandFlat::shoot(engine, (nbxing * gate) / gate))) - nbxing / 2;
0208       std::pair<int, int> digi(j + 1, time_hit);
0209       strips.insert(digi);
0210     }
0211   }
0212 }