RPCDigitizer/src/RPCSynchronizer.cc

0001 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0002 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
0003 #include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h"
0004 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0005 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0013 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0014 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0015 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0016 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0017 #include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h"
0018 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0019
0020 #include "CLHEP/Random/RandGaussQ.h"
0021
0022 #include <cstring>
0023 #include <iostream>
0024 #include <fstream>
0025 #include <string>
0026 #include <vector>
0027 #include <cstdlib>
0028 #include <cmath>
0029
0030 using namespace std;
0031
0032 RPCSynchronizer::RPCSynchronizer(const edm::ParameterSet& config) {
0033   resRPC = config.getParameter<double>("timeResolution");
0034   timOff = config.getParameter<double>("timingRPCOffset");
0035   dtimCs = config.getParameter<double>("deltatimeAdjacentStrip");
0036   resEle = config.getParameter<double>("timeJitter");
0037   sspeed = config.getParameter<double>("signalPropagationSpeed");
0038   lbGate = config.getParameter<double>("linkGateWidth");
0039   LHCGate = config.getParameter<double>("Gate");
0040   cosmics = config.getParameter<bool>("cosmics");
0041   irpc_timing_res = config.getParameter<double>("IRPC_time_resolution");
0042   irpc_electronics_jitter = config.getParameter<double>("IRPC_electronics_jitter");
0043   N_BX = config.getParameter<int>("BX_range");
0044   //"magic" parameter for cosmics
0045   cosmicPar = 37.62;
0046
0047   double c = 299792458;  // [m/s]
0048   //light speed in [cm/ns]
0049   cspeed = c * 1e+2 * 1e-9;
0050   //signal propagation speed [cm/ns]
0051   sspeed = sspeed * cspeed;
0052 }
0053
0054 RPCSynchronizer::~RPCSynchronizer() {}
0055
0056 int RPCSynchronizer::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* engine) {
0057   RPCSimSetUp* simsetup = this->getRPCSimSetUp();
0058   const RPCGeometry* geometry = simsetup->getGeometry();
0059   float timeref = simsetup->getTime(simhit->detUnitId());
0060
0061   int bx = -999;
0062   LocalPoint simHitPos = simhit->localPosition();
0063   float tof = simhit->timeOfFlight();
0064
0065   //automatic variable to prevent memory leak
0066
0067   float rr_el = CLHEP::RandGaussQ::shoot(engine, 0., resEle);
0068
0069   RPCDetId SimDetId(simhit->detUnitId());
0070
0071   const RPCRoll* SimRoll = nullptr;
0072
0073   for (TrackingGeometry::DetContainer::const_iterator it = geometry->dets().begin(); it != geometry->dets().end();
0074        it++) {
0075     if (dynamic_cast<const RPCChamber*>(*it) != nullptr) {
0076       auto ch = dynamic_cast<const RPCChamber*>(*it);
0077
0078       std::vector<const RPCRoll*> rollsRaf = (ch->rolls());
0079       for (std::vector<const RPCRoll*>::iterator r = rollsRaf.begin(); r != rollsRaf.end(); ++r) {
0080         if ((*r)->id() == SimDetId) {
0081           SimRoll = &(*(*r));
0082           break;
0083         }
0084       }
0085     }
0086   }
0087
0088   if (SimRoll != nullptr) {
0089     float distanceFromEdge = 0;
0090     float half_stripL = 0.;
0091
0092     if (SimRoll->id().region() == 0) {
0093       const RectangularStripTopology* top_ = dynamic_cast<const RectangularStripTopology*>(&(SimRoll->topology()));
0094       half_stripL = top_->stripLength() / 2;
0095       distanceFromEdge = half_stripL + simHitPos.y();
0096     } else {
0097       const TrapezoidalStripTopology* top_ = dynamic_cast<const TrapezoidalStripTopology*>(&(SimRoll->topology()));
0098       half_stripL = top_->stripLength() / 2;
0099       distanceFromEdge = half_stripL - simHitPos.y();
0100     }
0101
0102     float prop_time = distanceFromEdge / sspeed;
0103
0104     double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0., resRPC);
0105     double total_time = tof + prop_time + timOff + rr_tim1 + rr_el;
0106
0107     // Bunch crossing assignment
0108     double time_differ = 0.;
0109
0110     if (cosmics) {
0111       time_differ = (total_time - (timeref + ((half_stripL / sspeed) + timOff))) / cosmicPar;
0112     } else if (!cosmics) {
0113       time_differ = total_time - (timeref + (half_stripL / sspeed) + timOff);
0114     }
0115
0116     double inf_time = 0;
0117     double sup_time = 0;
0118
0119     for (int n = -N_BX; n <= N_BX; ++n) {
0120       if (cosmics) {
0121         inf_time = (-lbGate / 2 + n * LHCGate) / cosmicPar;
0122         sup_time = (lbGate / 2 + n * LHCGate) / cosmicPar;
0123       } else if (!cosmics) {
0124         inf_time = -lbGate / 2 + n * LHCGate;
0125         sup_time = lbGate / 2 + n * LHCGate;
0126       }
0127
0128       if (inf_time < time_differ && time_differ < sup_time) {
0129         bx = n;
0130         break;
0131       }
0132     }
0133   }
0134
0135   return bx;
0136 }
0137
0138 int RPCSynchronizer::getSimHitBxAndTimingForIRPC(const PSimHit* simhit, CLHEP::HepRandomEngine* engine) {
0139   RPCSimSetUp* simsetup = this->getRPCSimSetUp();
0140   const RPCGeometry* geometry = simsetup->getGeometry();
0141   float timeref = simsetup->getTime(simhit->detUnitId());
0142
0143   int bx = -999;
0144   LocalPoint simHitPos = simhit->localPosition();
0145   float tof = simhit->timeOfFlight();
0146
0147   //automatic variable to prevent memory leak
0148
0149   //  float rr_el = CLHEP::RandGaussQ::shoot(engine, 0.,resEle);
0150   float rr_el = CLHEP::RandGaussQ::shoot(engine, 0., irpc_electronics_jitter);
0151
0152   RPCDetId SimDetId(simhit->detUnitId());
0153
0154   const RPCRoll* SimRoll = nullptr;
0155
0156   for (TrackingGeometry::DetContainer::const_iterator it = geometry->dets().begin(); it != geometry->dets().end();
0157        it++) {
0158     if (dynamic_cast<const RPCChamber*>(*it) != nullptr) {
0159       auto ch = dynamic_cast<const RPCChamber*>(*it);
0160
0161       std::vector<const RPCRoll*> rollsRaf = (ch->rolls());
0162       for (std::vector<const RPCRoll*>::iterator r = rollsRaf.begin(); r != rollsRaf.end(); ++r) {
0163         if ((*r)->id() == SimDetId) {
0164           SimRoll = &(*(*r));
0165           break;
0166         }
0167       }
0168     }
0169   }
0170
0171   if (SimRoll != nullptr) {
0172     float distanceFromEdge = 0;
0173     float half_stripL = 0.;
0174
0175     if (SimRoll->id().region() == 0) {
0176       const RectangularStripTopology* top_ = dynamic_cast<const RectangularStripTopology*>(&(SimRoll->topology()));
0177       half_stripL = top_->stripLength() / 2;
0178       distanceFromEdge = half_stripL + simHitPos.y();
0179     } else {
0180       const TrapezoidalStripTopology* top_ = dynamic_cast<const TrapezoidalStripTopology*>(&(SimRoll->topology()));
0181       half_stripL = top_->stripLength() / 2;
0182       distanceFromEdge = half_stripL - simHitPos.y();
0183     }
0184
0185     float prop_time = distanceFromEdge / sspeed;
0186
0187     //    double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0.,resRPC);
0188     double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0., irpc_timing_res);
0189
0190     double total_time = tof + prop_time + timOff + rr_tim1 + rr_el;
0191
0192     // Bunch crossing assignment
0193     double time_differ = 0.;
0194
0195     if (cosmics) {
0196       time_differ = (total_time - (timeref + ((half_stripL / sspeed) + timOff))) / cosmicPar;
0197     } else if (!cosmics) {
0198       time_differ = total_time - (timeref + (half_stripL / sspeed) + timOff);
0199     }
0200
0201     double exact_total_time = tof + prop_time + timOff;
0202     double exact_time_differ = 0.;
0203
0204     if (cosmics) {
0205       exact_time_differ = (exact_total_time - (timeref + ((half_stripL / sspeed) + timOff))) / cosmicPar;
0206     } else if (!cosmics) {
0207       exact_time_differ = exact_total_time - (timeref + (half_stripL / sspeed) + timOff);
0208     }
0209
0210     double inf_time = 0;
0211     double sup_time = 0;
0212
0213     for (int n = -N_BX; n <= N_BX; ++n) {
0214       if (cosmics) {
0215         inf_time = (-lbGate / 2 + n * LHCGate) / cosmicPar;
0216         sup_time = (lbGate / 2 + n * LHCGate) / cosmicPar;
0217       } else if (!cosmics) {
0218         inf_time = -lbGate / 2 + n * LHCGate;
0219         sup_time = lbGate / 2 + n * LHCGate;
0220       }
0221
0222       if (inf_time < time_differ && time_differ < sup_time) {
0223         bx = n;
0224         the_exact_time = exact_time_differ;
0225         the_smeared_time = time_differ;
0226
0227         //cout<<"Debug\t"<<inf_time<<'\t'<<sup_time<<endl;
0228         ////if(bx)
0229         //  cout<<"Bingo\t"<<time_differ<<'\t'<<bx<<'\t'<<exact_time_differ<<'\t'<<exact_time_differ-time_differ<<'\t'<<exact_time_differ-bx*25.<<endl;
0230         break;
0231       }
0232     }
0233   }
0234   return bx;
0235 }