RPCDigitizer/src/RPCSimAverage.cc

0001 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0002 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
0003 #include "SimMuon/RPCDigitizer/src/RPCSimAverage.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 RPCSimAverage::RPCSimAverage(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
0051   if (rpcdigiprint) {
0052     std::cout << "Average Efficiency        = " << aveEff << std::endl;
0053     std::cout << "Average Cluster Size      = " << aveCls << " strips" << std::endl;
0054     std::cout << "RPC Time Resolution       = " << resRPC << " ns" << std::endl;
0055     std::cout << "RPC Signal formation time = " << timOff << " ns" << std::endl;
0056     std::cout << "RPC adjacent strip delay  = " << dtimCs << " ns" << std::endl;
0057     std::cout << "Electronic Jitter         = " << resEle << " ns" << std::endl;
0058     std::cout << "Signal propagation time   = " << sspeed << " x c" << std::endl;
0059     std::cout << "Link Board Gate Width     = " << lbGate << " ns" << std::endl;
0060   }
0061
0062   _rpcSync = new RPCSynchronizer(config);
0063 }
0064
0065 RPCSimAverage::~RPCSimAverage() { delete _rpcSync; }
0066
0067 int RPCSimAverage::getClSize(float posX, CLHEP::HepRandomEngine* engine) {
0068   std::map<int, std::vector<double> > clsMap = getRPCSimSetUp()->getClsMap();
0069
0070   int cnt = 1;
0071   int min = 1;
0072   double func = 0.0;
0073   std::vector<double> sum_clsize;
0074
0075   double rr_cl = CLHEP::RandFlat::shoot(engine);
0076   if (0.0 <= posX && posX < 0.2) {
0077     func = (clsMap[1])[(clsMap[1]).size() - 1] * (rr_cl);
0078     sum_clsize = clsMap[1];
0079   }
0080   if (0.2 <= posX && posX < 0.4) {
0081     func = (clsMap[2])[(clsMap[2]).size() - 1] * (rr_cl);
0082     sum_clsize = clsMap[2];
0083   }
0084   if (0.4 <= posX && posX < 0.6) {
0085     func = (clsMap[3])[(clsMap[3]).size() - 1] * (rr_cl);
0086     sum_clsize = clsMap[3];
0087   }
0088   if (0.6 <= posX && posX < 0.8) {
0089     func = (clsMap[4])[(clsMap[4]).size() - 1] * (rr_cl);
0090     sum_clsize = clsMap[4];
0091   }
0092   if (0.8 <= posX && posX < 1.0) {
0093     func = (clsMap[5])[(clsMap[5]).size() - 1] * (rr_cl);
0094     sum_clsize = clsMap[5];
0095   }
0096
0097   for (vector<double>::iterator iter = sum_clsize.begin(); iter != sum_clsize.end(); ++iter) {
0098     cnt++;
0099     if (func > (*iter)) {
0100       min = cnt;
0101     } else if (func < (*iter)) {
0102       break;
0103     }
0104   }
0105   return min;
0106 }
0107
0108 void RPCSimAverage::simulate(const RPCRoll* roll,
0109                              const edm::PSimHitContainer& rpcHits,
0110                              CLHEP::HepRandomEngine* engine) {
0111   _rpcSync->setRPCSimSetUp(getRPCSimSetUp());
0112   theRpcDigiSimLinks.clear();
0113   theDetectorHitMap.clear();
0114   theRpcDigiSimLinks = RPCDigiSimLinks(roll->id().rawId());
0115
0116   const Topology& topology = roll->specs()->topology();
0117
0118   for (edm::PSimHitContainer::const_iterator _hit = rpcHits.begin(); _hit != rpcHits.end(); ++_hit) {
0119     // Here I hould check if the RPC are up side down;
0120     const LocalPoint& entr = _hit->entryPoint();
0121
0122     //    const LocalPoint& exit=_hit->exitPoint();
0123
0124     float posX = roll->strip(_hit->localPosition()) - static_cast<int>(roll->strip(_hit->localPosition()));
0125     int time_hit = _rpcSync->getSimHitBx(&(*_hit), engine);
0126
0127     // Effinciecy
0128
0129     if (CLHEP::RandFlat::shoot(engine) < aveEff) {
0130       int centralStrip = topology.channel(entr) + 1;
0131       int fstrip = centralStrip;
0132       int lstrip = centralStrip;
0133       // Compute the cluster size
0134       //double w = CLHEP::RandFlat::shoot(engine);
0135       //if (w < 1.e-10) w=1.e-10;
0136       int clsize = this->getClSize(posX, engine);
0137
0138       std::vector<int> cls;
0139       cls.push_back(centralStrip);
0140       if (clsize > 1) {
0141         for (int cl = 0; cl < (clsize - 1) / 2; cl++)
0142           if (centralStrip - cl - 1 >= 1) {
0143             fstrip = centralStrip - cl - 1;
0144             cls.push_back(fstrip);
0145           }
0146         for (int cl = 0; cl < (clsize - 1) / 2; cl++)
0147           if (centralStrip + cl + 1 <= roll->nstrips()) {
0148             lstrip = centralStrip + cl + 1;
0149             cls.push_back(lstrip);
0150           }
0151         if (clsize % 2 == 0) {
0152           // insert the last strip according to the
0153           // simhit position in the central strip
0154           double deltaw = roll->centreOfStrip(centralStrip).x() - entr.x();
0155           if (deltaw < 0.) {
0156             if (lstrip < roll->nstrips()) {
0157               lstrip++;
0158               cls.push_back(lstrip);
0159             }
0160           } else {
0161             if (fstrip > 1) {
0162               fstrip--;
0163               cls.push_back(fstrip);
0164             }
0165           }
0166         }
0167       }
0168
0169       for (std::vector<int>::iterator i = cls.begin(); i != cls.end(); i++) {
0170         // Check the timing of the adjacent strip
0171         std::pair<int, int> digi(*i, time_hit);
0172         theDetectorHitMap.insert(DetectorHitMap::value_type(digi, &(*_hit)));
0173         strips.insert(digi);
0174       }
0175     }
0176   }
0177 }
0178
0179 void RPCSimAverage::simulateNoise(const RPCRoll* roll, CLHEP::HepRandomEngine* engine) {
0180   RPCDetId rpcId = roll->id();
0181   int nstrips = roll->nstrips();
0182   double area = 0.0;
0183
0184   if (rpcId.region() == 0) {
0185     const RectangularStripTopology* top_ = dynamic_cast<const RectangularStripTopology*>(&(roll->topology()));
0186     float xmin = (top_->localPosition(0.)).x();
0187     float xmax = (top_->localPosition((float)roll->nstrips())).x();
0188     float striplength = (top_->stripLength());
0189     area = striplength * (xmax - xmin);
0190   } else {
0191     const TrapezoidalStripTopology* top_ = dynamic_cast<const TrapezoidalStripTopology*>(&(roll->topology()));
0192     float xmin = (top_->localPosition(0.)).x();
0193     float xmax = (top_->localPosition((float)roll->nstrips())).x();
0194     float striplength = (top_->stripLength());
0195     area = striplength * (xmax - xmin);
0196   }
0197
0198   double ave = rate * nbxing * gate * area * 1.0e-9;
0199
0200   CLHEP::RandPoissonQ randPoissonQ(*engine, ave);
0201   N_hits = randPoissonQ.fire();
0202
0203   for (int i = 0; i < N_hits; i++) {
0204     int strip = static_cast<int>(CLHEP::RandFlat::shoot(engine, 1, nstrips));
0205     int time_hit;
0206     time_hit = (static_cast<int>(CLHEP::RandFlat::shoot(engine, (nbxing * gate) / gate))) - nbxing / 2;
0207     std::pair<int, int> digi(strip, time_hit);
0208     strips.insert(digi);
0209   }
0210 }