diff --git a/Configuration/Eras/python/Modifier_phase2_rpc_devel_cff.py b/Configuration/Eras/python/Modifier_phase2_rpc_devel_cff.py new file mode 100644 index 0000000000000..3c98dc624011d --- /dev/null +++ b/Configuration/Eras/python/Modifier_phase2_rpc_devel_cff.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms + +phase2_rpc_devel = cms.Modifier() + diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index ba76d044bb352..b6edd93f92a45 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -3397,6 +3397,53 @@ def condition(self, fragment, stepList, key, hasHarvest): offset = 0.85, ) +# RPC development WF +class UpgradeWorkflow_rpcDevel(UpgradeWorkflow): + def __init__(self, digi = {}, reco = {}, harvest = {}, **kwargs): + super(UpgradeWorkflow_rpcDevel, self).__init__( + steps = [ + 'DigiTrigger', + 'RecoGlobal', + 'RecoGlobalFakeHLT', + 'HARVESTGlobal', + 'HARVESTGlobalFakeHLT', + 'ALCAPhase2', + ], + PU = [ + 'DigiTrigger', + 'RecoGlobal', + 'RecoGlobalFakeHLT', + 'HARVESTGlobal', + 'HARVESTGlobalFakeHLT', + 'ALCAPhase2', + ], + **kwargs + ) + self.__digi = digi + self.__reco = reco + self.__harvest = harvest + + def setup_(self, step, stepName, stepDict, k, properties): + mods = {'--era': stepDict[step][k]['--era'] + ',phase2_rpc_devel'} + + if 'Digi' in step: + mods |= self.__digi + elif 'Reco' in step: + mods |= self.__reco + elif 'HARVEST' in step: + mods |= self.__harvest + + stepDict[stepName][k] = merge([mods, stepDict[step][k]]) + + def condition(self, fragment, stepList, key, hasHarvest): + return fragment == "TTbar_14TeV" and 'Run4' in key + + +upgradeWFs['rpcDevel'] = UpgradeWorkflow_rpcDevel( + suffix = '_rpcDevel', + offset = 0.62, +) + # check for duplicates in offsets or suffixes offsets = [specialWF.offset for specialType,specialWF in upgradeWFs.items()] suffixes = [specialWF.suffix for specialType,specialWF in upgradeWFs.items()] diff --git a/DataFormats/RPCDigi/interface/IRPCDigi.h b/DataFormats/RPCDigi/interface/IRPCDigi.h new file mode 100644 index 0000000000000..784c8dac0fc59 --- /dev/null +++ b/DataFormats/RPCDigi/interface/IRPCDigi.h @@ -0,0 +1,44 @@ +#ifndef IRPCDigi_IRPCDigi_h +#define IRPCDigi_IRPCDigi_h + +/** \class IRPCDigi + * + * Digi for Improved Resisitive Plate Chamber (IRPC) + * + * \author Borislav Pavlov - University of Sofia +*/ + +#include +#include + +class IRPCDigi { +public: + explicit IRPCDigi(int strip, int bxLR, int bxHR, int sbxLR, int sbxHR, int fineLR, int fineHR); + IRPCDigi(); + + bool operator==(const IRPCDigi& digi) const; + bool operator<(const IRPCDigi& digi) const; + void print() const; + int strip() const { return strip_; } + int bx() const { return bxLR_; } + int sbx() const { return sbxLR_; } + int bxLR() const { return bxLR_; } + int bxHR() const { return bxHR_; } + int sbxLR() const { return sbxLR_; } + int sbxHR() const { return sbxHR_; } + int fineLR() const { return fineLR_; } + int fineHR() const { return fineHR_; } + +private: + uint16_t strip_; + int32_t bxLR_; //BX from low radius FEB + int32_t bxHR_; //BX from high radius FEB + int8_t sbxLR_; //sub-BX from low radius FEB + int8_t sbxHR_; //sub-BX from high radius FEB + int8_t fineLR_; //high resolution time rom low radius FEB + int8_t fineHR_; //high resolution time rom high radius FEB +}; + +std::ostream& operator<<(std::ostream& o, const IRPCDigi& digi); + +#endif diff --git a/DataFormats/RPCDigi/interface/IRPCDigiCollection.h b/DataFormats/RPCDigi/interface/IRPCDigiCollection.h new file mode 100644 index 0000000000000..fa9641db3feaf --- /dev/null +++ b/DataFormats/RPCDigi/interface/IRPCDigiCollection.h @@ -0,0 +1,15 @@ +#ifndef IRPCDigi_IRPCDigiCollection_h +#define IRPCDigi_IRPCDigiCollection_h +/** \class RPCDigiCollection + * + * \author Borislav Pavlov + * \date 14 July 2021 + */ + +#include +#include +#include + +typedef MuonDigiCollection IRPCDigiCollection; + +#endif diff --git a/DataFormats/RPCDigi/interface/IRPCDigiTime.h b/DataFormats/RPCDigi/interface/IRPCDigiTime.h new file mode 100644 index 0000000000000..3f4afe53a2be4 --- /dev/null +++ b/DataFormats/RPCDigi/interface/IRPCDigiTime.h @@ -0,0 +1,14 @@ +#include "DataFormats/RPCDigi/interface/IRPCDigi.h" + +class IRPCDigiTime { +public: + IRPCDigiTime(const IRPCDigi& adigi); + float time(); + float coordinateY(); + float timeLR(); + float timeHR(); + +private: + IRPCDigi theDigi; + float TDC2Time(int BX, int SBX, int FT); +}; diff --git a/DataFormats/RPCDigi/interface/RPCDigiPhase2.h b/DataFormats/RPCDigi/interface/RPCDigiPhase2.h new file mode 100644 index 0000000000000..ca83bbba436c3 --- /dev/null +++ b/DataFormats/RPCDigi/interface/RPCDigiPhase2.h @@ -0,0 +1,34 @@ +#ifndef RPCDigi_RPCDigiPhase2_h +#define RPCDigi_RPCDigiPhase2_h + +/** \class RPCDigiPhase2 + * + * Digi for Resisitive Plate Chamber, after Phase2 upgrade + * + * \author Borislav Pavlov - University of Sofia +*/ + +#include +#include + +class RPCDigiPhase2 { +public: + explicit RPCDigiPhase2(int strip, int bx, int sbx); + RPCDigiPhase2(); + + bool operator==(const RPCDigiPhase2& digi) const; + bool operator<(const RPCDigiPhase2& digi) const; + void print() const; + int strip() const { return strip_; } + int bx() const { return bx_; } + int sbx() const { return sbx_; } + +private: + uint16_t strip_; + int32_t bx_; // for BX + int8_t sbx_; // for sub-BX +}; + +std::ostream& operator<<(std::ostream& o, const RPCDigiPhase2& digi); + +#endif diff --git a/DataFormats/RPCDigi/interface/RPCDigiPhase2Collection.h b/DataFormats/RPCDigi/interface/RPCDigiPhase2Collection.h new file mode 100644 index 0000000000000..33ed49df28e73 --- /dev/null +++ b/DataFormats/RPCDigi/interface/RPCDigiPhase2Collection.h @@ -0,0 +1,15 @@ +#ifndef RPCDigi_RPCDigiPhase2Collection_h +#define RPCDigi_RPCDigiPhase2Collection_h +/** \class RPCDigiCollection + * + * \author Borislav Pavlov + * \date 14 June 2024 + */ + +#include +#include +#include + +typedef MuonDigiCollection RPCDigiPhase2Collection; + +#endif diff --git a/DataFormats/RPCDigi/interface/RPCDigiPhase2Time.h b/DataFormats/RPCDigi/interface/RPCDigiPhase2Time.h new file mode 100644 index 0000000000000..c3128ed8d2cd5 --- /dev/null +++ b/DataFormats/RPCDigi/interface/RPCDigiPhase2Time.h @@ -0,0 +1,10 @@ +#include "DataFormats/RPCDigi/interface/RPCDigiPhase2.h" + +class RPCDigiPhase2Time { +public: + RPCDigiPhase2Time(const RPCDigiPhase2& adigi); + float time(); + +private: + RPCDigiPhase2 theDigi; +}; diff --git a/DataFormats/RPCDigi/src/IRPCDigi.cc b/DataFormats/RPCDigi/src/IRPCDigi.cc new file mode 100644 index 0000000000000..5329220272c4c --- /dev/null +++ b/DataFormats/RPCDigi/src/IRPCDigi.cc @@ -0,0 +1,33 @@ +/** \file + * + * + * \author Borislav Pavlov - University of Sofia + * + */ + +#include "DataFormats/RPCDigi/interface/IRPCDigi.h" +#include + +IRPCDigi::IRPCDigi(int strip, int bxLR, int bxHR, int sbxLR, int sbxHR, int fineLR, int fineHR) + : strip_(strip), bxLR_(bxLR), bxHR_(bxHR), sbxLR_(sbxLR), sbxHR_(sbxHR), fineLR_(fineLR), fineHR_(fineHR) {} + +IRPCDigi::IRPCDigi() : strip_(0), bxLR_(0), sbxLR_(0) {} + +// Comparison +bool IRPCDigi::operator==(const IRPCDigi& digi) const { + if (strip_ != digi.strip() || bxLR_ != digi.bx()) + return false; + return true; +} + +///Precedence operator +bool IRPCDigi::operator<(const IRPCDigi& digi) const { + if (digi.bx() == this->bx()) + return digi.strip() < this->strip(); + else + return digi.bx() < this->bx(); +} + +std::ostream& operator<<(std::ostream& o, const IRPCDigi& digi) { return o << " " << digi.strip() << " " << digi.bx(); } + +void IRPCDigi::print() const { std::cout << "Strip " << strip() << " bx " << bx() << std::endl; } diff --git a/DataFormats/RPCDigi/src/IRPCDigiTime.cc b/DataFormats/RPCDigi/src/IRPCDigiTime.cc new file mode 100644 index 0000000000000..e339fc5b06f2a --- /dev/null +++ b/DataFormats/RPCDigi/src/IRPCDigiTime.cc @@ -0,0 +1,16 @@ +#include "DataFormats/RPCDigi/interface/IRPCDigiTime.h" + +IRPCDigiTime::IRPCDigiTime(const IRPCDigi& adigi) : theDigi(adigi) {} + +float IRPCDigiTime::time() { return (timeLR() + timeHR()) / 2.; } + +float IRPCDigiTime::coordinateY() { + const double signal_speed = 0.66 * 299792458e-7; //signal propagation speed [cm/ns] + return signal_speed * (timeLR() - timeHR()) / 2.; +} + +float IRPCDigiTime::timeLR() { return TDC2Time(theDigi.bxLR(), theDigi.sbxLR(), theDigi.fineLR()); } + +float IRPCDigiTime::timeHR() { return TDC2Time(theDigi.bxHR(), theDigi.sbxHR(), theDigi.fineHR()); } + +float IRPCDigiTime::TDC2Time(int BX, int SBX, int FT) { return 25. * BX + 2.5 * SBX + 0.2 * FT; } diff --git a/DataFormats/RPCDigi/src/RPCDigiPhase2.cc b/DataFormats/RPCDigi/src/RPCDigiPhase2.cc new file mode 100644 index 0000000000000..58d9842cc3ae4 --- /dev/null +++ b/DataFormats/RPCDigi/src/RPCDigiPhase2.cc @@ -0,0 +1,34 @@ +/** \file + * + * + * \author Borislav Pavlov - University of Sofia + * + */ + +#include "DataFormats/RPCDigi/interface/RPCDigiPhase2.h" +#include + +RPCDigiPhase2::RPCDigiPhase2(int strip, int bx, int sbx) : strip_(strip), bx_(bx), sbx_(sbx) {} + +RPCDigiPhase2::RPCDigiPhase2() : strip_(0), bx_(0), sbx_(0) {} + +// Comparison +bool RPCDigiPhase2::operator==(const RPCDigiPhase2& digi) const { + if (strip_ != digi.strip() || bx_ != digi.bx()) + return false; + return true; +} + +///Precedence operator +bool RPCDigiPhase2::operator<(const RPCDigiPhase2& digi) const { + if (digi.bx() == this->bx()) + return digi.strip() < this->strip(); + else + return digi.bx() < this->bx(); +} + +std::ostream& operator<<(std::ostream& o, const RPCDigiPhase2& digi) { + return o << " " << digi.strip() << " " << digi.bx(); +} + +void RPCDigiPhase2::print() const { std::cout << "Strip " << strip() << " bx " << bx() << std::endl; } diff --git a/DataFormats/RPCDigi/src/RPCDigiPhase2Time.cc b/DataFormats/RPCDigi/src/RPCDigiPhase2Time.cc new file mode 100644 index 0000000000000..16006aa4d8060 --- /dev/null +++ b/DataFormats/RPCDigi/src/RPCDigiPhase2Time.cc @@ -0,0 +1,7 @@ +#include "DataFormats/RPCDigi/interface/RPCDigiPhase2Time.h" + +RPCDigiPhase2Time::RPCDigiPhase2Time(const RPCDigiPhase2& adigi) : theDigi(adigi) {} + +float RPCDigiPhase2Time::time() { + return 25. * theDigi.bx() + 1.5625 * theDigi.sbx(); // 25./16. = 1.5625 ns +} diff --git a/DataFormats/RPCDigi/src/classes.h b/DataFormats/RPCDigi/src/classes.h index 8a38045fdb89f..37e295e8f8815 100644 --- a/DataFormats/RPCDigi/src/classes.h +++ b/DataFormats/RPCDigi/src/classes.h @@ -1,5 +1,9 @@ #include #include +#include +#include +#include +#include #include "DataFormats/RPCDigi/interface/RPCRawDataCounts.h" #include "DataFormats/RPCDigi/interface/RPCRawSynchro.h" #include "DataFormats/RPCDigi/interface/RPCDigiL1Link.h" diff --git a/DataFormats/RPCDigi/src/classes_def.xml b/DataFormats/RPCDigi/src/classes_def.xml index 5f5842115f977..be74804b50639 100644 --- a/DataFormats/RPCDigi/src/classes_def.xml +++ b/DataFormats/RPCDigi/src/classes_def.xml @@ -8,6 +8,24 @@ + + + + + + + + + + + + + + + + + + diff --git a/IOMC/RandomEngine/python/IOMC_cff.py b/IOMC/RandomEngine/python/IOMC_cff.py index c8c2b93668ca0..be87675dcadd4 100644 --- a/IOMC/RandomEngine/python/IOMC_cff.py +++ b/IOMC/RandomEngine/python/IOMC_cff.py @@ -53,6 +53,14 @@ initialSeed = cms.untracked.uint32(1234567), engineName = FullSimEngine ), + simMuonIRPCDigis = cms.PSet( + initialSeed = cms.untracked.uint32(1234567), + engineName = FullSimEngine + ), + simMuonRPCDigisPhase2 = cms.PSet( + initialSeed = cms.untracked.uint32(1234567), + engineName = FullSimEngine + ), # # HI generation & simulation is a special processing/chain, # integrated since 330 cycle diff --git a/SimMuon/Configuration/python/SimMuon_EventContent_cff.py b/SimMuon/Configuration/python/SimMuon_EventContent_cff.py index ba911c199eede..e8ea83ff7caa2 100644 --- a/SimMuon/Configuration/python/SimMuon_EventContent_cff.py +++ b/SimMuon/Configuration/python/SimMuon_EventContent_cff.py @@ -74,6 +74,9 @@ phase2_muon.toModify( SimMuonRECO, outputCommands = SimMuonRECO.outputCommands + ['keep *DigiSimLinkedmDetSetVector_simMuonME0Digis_*_*'] ) phase2_muon.toModify( SimMuonPREMIX, outputCommands = SimMuonPREMIX.outputCommands + ['keep *_mix_g4SimHitsMuonME0Hits_*'] ) +from Configuration.Eras.Modifier_phase2_rpc_devel_cff import phase2_rpc_devel +phase2_rpc_devel.toModify( SimMuonFEVTDEBUG, outputCommands = SimMuonFEVTDEBUG.outputCommands + ['keep *_simMuonRPCDigisPhase2_*_*', + 'keep *_simMuonIRPCDigis_*_*'] ) # FastSim uses different naming convention from Configuration.Eras.Modifier_fastSim_cff import fastSim diff --git a/SimMuon/Configuration/python/SimMuon_cff.py b/SimMuon/Configuration/python/SimMuon_cff.py index a39f7f8f103ab..a2a7f370ad0a9 100644 --- a/SimMuon/Configuration/python/SimMuon_cff.py +++ b/SimMuon/Configuration/python/SimMuon_cff.py @@ -1,17 +1,23 @@ import FWCore.ParameterSet.Config as cms # Muon Digitization (CSC, DT, RPC electronics responce) + # CSC digitizer # from SimMuon.CSCDigitizer.muonCSCDigis_cfi import * from CalibMuon.CSCCalibration.CSCChannelMapper_cfi import * from CalibMuon.CSCCalibration.CSCIndexer_cfi import * + # DT digitizer # from SimMuon.DTDigitizer.muondtdigi_cfi import * + # RPC digitizer -# -from SimMuon.RPCDigitizer.muonrpcdigi_cfi import * +# +# If your branch still uses the legacy filename muonrpcdigi_cfi.py, +# change only this import line back to that filename. +from SimMuon.RPCDigitizer.muonRPCDigis_cfi import * + muonDigiTask = cms.Task(simMuonCSCDigis, simMuonDTDigis, simMuonRPCDigis) muonDigi = cms.Sequence(muonDigiTask) @@ -27,6 +33,14 @@ # while GE0 is in development, just turn off ME0 tasks _phase2_ge0 = _phase2_muonDigiTask.copyAndExclude([muonME0DigiTask]) +# phase2_rpc_devel: +# keep legacy simMuonRPCDigis for side-by-side validation +# and add the new Phase-2 RPC + iRPC digis +_phase2_rpc_devel = _phase2_muonDigiTask.copyAndExclude([muonME0DigiTask]) +_phase2_rpc_devel.add(simMuonRPCDigisPhase2) +_phase2_rpc_devel.add(simMuonIRPCDigis) + + from Configuration.Eras.Modifier_run2_GEM_2017_cff import run2_GEM_2017 run2_GEM_2017.toReplaceWith( muonDigiTask, _run3_muonDigiTask ) from Configuration.Eras.Modifier_run3_GEM_cff import run3_GEM @@ -35,4 +49,5 @@ phase2_muon.toReplaceWith( muonDigiTask, _phase2_muonDigiTask ) from Configuration.Eras.Modifier_phase2_GE0_cff import phase2_GE0 phase2_GE0.toReplaceWith( muonDigiTask, _phase2_ge0 ) - +from Configuration.Eras.Modifier_phase2_rpc_devel_cff import phase2_rpc_devel +phase2_rpc_devel.toReplaceWith(muonDigiTask, _phase2_rpc_devel) \ No newline at end of file diff --git a/SimMuon/RPCDigitizer/python/muonRPCDigis_cfi.py b/SimMuon/RPCDigitizer/python/muonRPCDigis_cfi.py index 6b7f809c1f48e..6001b0b04697e 100644 --- a/SimMuon/RPCDigitizer/python/muonRPCDigis_cfi.py +++ b/SimMuon/RPCDigitizer/python/muonRPCDigis_cfi.py @@ -35,7 +35,7 @@ #the digitizer for PhaseII muon upgrade is RPCSimModelTiming and for the moment is based on RPCSimAverageNoiseEffCls from Configuration.Eras.Modifier_fastSim_cff import fastSim -_simMuonRPCDigisPhaseII = cms.EDProducer("RPCandIRPCDigiProducer", +_simMuonRPCDigisUpgradeTDR = cms.EDProducer("RPCandIRPCDigiProducer", Noise = cms.bool(True), digiModelConfig = cms.PSet( signalPropagationSpeed = cms.double(0.66), @@ -89,8 +89,72 @@ digiIRPCModel = cms.string('RPCSimModelTiming') ) + +simMuonRPCDigisPhase2 = cms.EDProducer("RPCDigiPhase2Producer", + Noise = cms.bool(True), + digiModelConfig = cms.PSet( + signalPropagationSpeed = cms.double(0.66), + timingRPCOffset = cms.double(50.0), + Frate = cms.double(1.0), + printOutDigitizer = cms.bool(False), + cosmics = cms.bool(False), + deltatimeAdjacentStrip = cms.double(3.0), + linkGateWidth = cms.double(20.0), + Rate = cms.double(0.0), + timeResolution = cms.double(1.5), + averageClusterSize = cms.double(1.5), + Gate = cms.double(25.0), + averageEfficiency = cms.double(0.95), + Nbxing = cms.int32(9), + BX_range = cms.int32(5), + timeJitter = cms.double(0.1), + sigmaY = cms.double(2.), # resolution of 2 cm + do_Y_coordinate = cms.bool(False), + digitizeElectrons = cms.bool(False), + IRPC_time_resolution = cms.double(1.5),# intrinsic time resolution of 1.5 ns + IRPC_electronics_jitter = cms.double(0.1)# resolution of 100 ps + ), + doBkgNoise = cms.bool(False), #False - no noise and bkg simulation + Signal = cms.bool(True), + mixLabel = cms.string('mix'), + InputCollection = cms.string('g4SimHitsMuonRPCHits'), + digiModel = cms.string('RPCSimModelTimingPhase2') +) + +#the digitizer for IRPC chambers is IRPCSimModelTiming and is based on RPCSimAverageNoiseEffCls +simMuonIRPCDigis = cms.EDProducer("IRPCDigiProducer", + Noise = cms.bool(True), + digiIRPCModelConfig = cms.PSet( + signalPropagationSpeed = cms.double(0.66), + timingRPCOffset = cms.double(50.0), + Frate = cms.double(1.0), + printOutDigitizer = cms.bool(False), + cosmics = cms.bool(False), + deltatimeAdjacentStrip = cms.double(3.0), + linkGateWidth = cms.double(20.0), + Rate = cms.double(0.0), + timeResolution = cms.double(1.0), + averageClusterSize = cms.double(1.5), + Gate = cms.double(25.0), + averageEfficiency = cms.double(0.95), + Nbxing = cms.int32(9), + BX_range = cms.int32(5), + timeJitter = cms.double(0.1), + IRPC_time_resolution = cms.double(1),# resolution of 1 ns + IRPC_electronics_jitter = cms.double(0.1),# resolution of 100 ps + sigmaY = cms.double(2.), # resolution of 2 cm + do_Y_coordinate = cms.bool(True), + digitizeElectrons = cms.bool(False), + ), + doBkgNoise = cms.bool(False), #False - no noise and bkg simulation + Signal = cms.bool(True), + mixLabel = cms.string('mix'), + InputCollection = cms.string('g4SimHitsMuonRPCHits'), + digiIRPCModel = cms.string('IRPCSimModelTiming') +) + from Configuration.Eras.Modifier_phase2_muon_cff import phase2_muon -phase2_muon.toReplaceWith( simMuonRPCDigis, _simMuonRPCDigisPhaseII ) +phase2_muon.toReplaceWith( simMuonRPCDigis, _simMuonRPCDigisUpgradeTDR ) from Configuration.ProcessModifiers.premix_stage2_cff import premix_stage2 premix_stage2.toModify(simMuonRPCDigis, mixLabel = "mixData") diff --git a/SimMuon/RPCDigitizer/src/IRPCDigiProducer.cc b/SimMuon/RPCDigitizer/src/IRPCDigiProducer.cc new file mode 100644 index 0000000000000..f7454fa41cbac --- /dev/null +++ b/SimMuon/RPCDigitizer/src/IRPCDigiProducer.cc @@ -0,0 +1,109 @@ +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" +#include "SimMuon/RPCDigitizer/src/IRPCDigiProducer.h" +#include "SimMuon/RPCDigitizer/src/IRPCDigitizer.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h" +#include +#include + +#include +#include + +#include "FWCore/Framework/interface/MakerMacros.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" + +//Random Number +#include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "CLHEP/Random/RandFlat.h" + +namespace CLHEP { + class HepRandomEngine; +} + +IRPCDigiProducer::IRPCDigiProducer(const edm::ParameterSet& ps) { + produces(); + produces("IRPCDigiSimLink"); + + //Name of Collection used for create the XF + mix_ = ps.getParameter("mixLabel"); + collection_for_XF = ps.getParameter("InputCollection"); + + edm::Service rng; + if (!rng.isAvailable()) { + throw cms::Exception("Configuration") + << "RPCDigitizer requires the RandomNumberGeneratorService\n" + "which is not present in the configuration file. You must add the service\n" + "in the configuration file or remove the modules that require it."; + }; + + theRPCSimSetUpIRPC = new RPCSimSetUp(ps); + theIRPCDigitizer = new IRPCDigitizer(ps); + crossingFrameToken = consumes>(edm::InputTag(mix_, collection_for_XF)); + geomToken = esConsumes(); + noiseToken = esConsumes(); + clsToken = esConsumes(); +} + +IRPCDigiProducer::~IRPCDigiProducer() { + delete theIRPCDigitizer; + delete theRPCSimSetUpIRPC; +} + +void IRPCDigiProducer::beginRun(const edm::Run& r, const edm::EventSetup& eventSetup) { + edm::ESHandle hGeom = eventSetup.getHandle(geomToken); + const RPCGeometry* pGeom = &*hGeom; + _pGeom = &*hGeom; + + edm::ESHandle noiseRcd = eventSetup.getHandle(noiseToken); + + edm::ESHandle clsRcd = eventSetup.getHandle(clsToken); + //eventSetup.get().get(clsRcd); + + //setup the two digi models + theRPCSimSetUpIRPC->setGeometry(pGeom); + theRPCSimSetUpIRPC->setRPCSetUp(noiseRcd->getVNoise(), clsRcd->getCls()); + + //setup the two digitizers + theIRPCDigitizer->setGeometry(pGeom); + theIRPCDigitizer->setRPCSimSetUp(theRPCSimSetUpIRPC); +} + +void IRPCDigiProducer::produce(edm::Event& e, const edm::EventSetup& eventSetup) { + edm::Service rng; + CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID()); + + LogDebug("IRPCDigiProducer") << "[IRPCDigiProducer::produce] got the CLHEP::HepRandomEngine engine from " + "the edm::Event.streamID() and edm::Service"; + LogDebug("IRPCDigiProducer") << "[IRPCDigiProducer::produce] test the CLHEP::HepRandomEngine by firing " + "once RandFlat ---- this must be the first time in SimMuon/RPCDigitizer"; + LogDebug("IRPCDigiProducer") << "[IRPCDigiProducer::produce] to activate the test go in " + "IRPCDigiProducer.cc and uncomment the line below"; + + const edm::Handle>& cf = e.getHandle(crossingFrameToken); + + std::unique_ptr> hits(new MixCollection(cf.product())); + + // Create empty output + std::unique_ptr pDigis(new IRPCDigiCollection()); + std::unique_ptr IRPCDigitSimLink(new IRPCDigitizerSimLinks()); + + theIRPCDigitizer->doAction(*hits, *pDigis, *IRPCDigitSimLink, engine); //make "IRPC" digitizer do the action + + e.put(std::move(pDigis)); + //store the SimDigiLinks in the event + e.put(std::move(IRPCDigitSimLink), "IRPCDigiSimLink"); +} diff --git a/SimMuon/RPCDigitizer/src/IRPCDigiProducer.h b/SimMuon/RPCDigitizer/src/IRPCDigiProducer.h new file mode 100644 index 0000000000000..a5273eef9434f --- /dev/null +++ b/SimMuon/RPCDigitizer/src/IRPCDigiProducer.h @@ -0,0 +1,57 @@ +#ifndef IRPCDigiProducer_h +#define IRPCDigiProducer_h + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimMuon/RPCDigitizer/src/IRPCDigitizer.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "CondFormats/RPCObjects/interface/RPCStripNoises.h" +#include "CondFormats/DataRecord/interface/RPCStripNoisesRcd.h" +#include "CondFormats/RPCObjects/interface/RPCClusterSize.h" +#include "CondFormats/DataRecord/interface/RPCClusterSizeRcd.h" + +class RPCGeometry; +class RPCSimSetUp; +class RPCSynchronizer; + +class IRPCDigiProducer : public edm::stream::EDProducer { +public: + // typedef edm::DetSetVector RPCDigiSimLinks; + typedef IRPCDigitizer::RPCDigiSimLinks IRPCDigitizerSimLinks; + + explicit IRPCDigiProducer(const edm::ParameterSet& ps); + ~IRPCDigiProducer() override; + + void beginRun(const edm::Run&, const edm::EventSetup&) override; + + /**Produces the EDM products,*/ + void produce(edm::Event& e, const edm::EventSetup& c) override; + + void setRPCSetUp(const std::vector&, const std::vector&); + +private: + IRPCDigitizer* theIRPCDigitizer; + RPCSimSetUp* theRPCSimSetUpIRPC; + // RPCSimSetUp* theRPCSimSetUp; + + //Name of Collection used for create the XF + std::string mix_; + std::string collection_for_XF; + + //Token for accessing data + edm::EDGetTokenT> crossingFrameToken; + const RPCGeometry* _pGeom; + + //EventSetup Tokens + edm::ESGetToken geomToken; + edm::ESGetToken noiseToken; + edm::ESGetToken clsToken; +}; + +#endif diff --git a/SimMuon/RPCDigitizer/src/IRPCDigitizer.cc b/SimMuon/RPCDigitizer/src/IRPCDigitizer.cc index abd7ea4258197..5c9f4a3735edf 100644 --- a/SimMuon/RPCDigitizer/src/IRPCDigitizer.cc +++ b/SimMuon/RPCDigitizer/src/IRPCDigitizer.cc @@ -18,7 +18,7 @@ IRPCDigitizer::IRPCDigitizer(const edm::ParameterSet& config) IRPCDigitizer::~IRPCDigitizer() = default; void IRPCDigitizer::doAction(MixCollection& simHits, - RPCDigiCollection& rpcDigis, + IRPCDigiCollection& rpcDigis, RPCDigiSimLinks& rpcDigiSimLink, CLHEP::HepRandomEngine* engine) { theRPCSim->setRPCSimSetUp(theSimSetUp); @@ -49,7 +49,9 @@ void IRPCDigitizer::doAction(MixCollection& simHits, } theRPCSim->fillDigis((*r)->id(), rpcDigis); - rpcDigiSimLink.insert(theRPCSim->rpcDigiSimLinks()); + if (rpcDigiSimLink.find((theRPCSim->rpcDigiSimLinks()).detId()) == rpcDigiSimLink.end()) { + rpcDigiSimLink.insert(theRPCSim->rpcDigiSimLinks()); + } } } diff --git a/SimMuon/RPCDigitizer/src/IRPCDigitizer.h b/SimMuon/RPCDigitizer/src/IRPCDigitizer.h index fb054205b266b..b7d7b748c642b 100644 --- a/SimMuon/RPCDigitizer/src/IRPCDigitizer.h +++ b/SimMuon/RPCDigitizer/src/IRPCDigitizer.h @@ -12,7 +12,7 @@ #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h" #include "SimDataFormats/RPCDigiSimLink/interface/RPCDigiSimLink.h" #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" -#include "DataFormats/RPCDigi/interface/RPCDigiCollection.h" +#include "DataFormats/RPCDigi/interface/IRPCDigiCollection.h" #include "Geometry/RPCGeometry/interface/RPCGeometry.h" #include "SimDataFormats/CrossingFrame/interface/MixCollection.h" #include @@ -38,7 +38,7 @@ class IRPCDigitizer { // *** digitize *** void doAction(MixCollection& simHits, - RPCDigiCollection& rpcDigis, + IRPCDigiCollection& rpcDigis, RPCDigiSimLinks& rpcDigiSimLink, CLHEP::HepRandomEngine*); diff --git a/SimMuon/RPCDigitizer/src/IRPCSimModelTiming.cc b/SimMuon/RPCDigitizer/src/IRPCSimModelTiming.cc new file mode 100644 index 0000000000000..3569647c9f517 --- /dev/null +++ b/SimMuon/RPCDigitizer/src/IRPCSimModelTiming.cc @@ -0,0 +1,280 @@ +#include "Geometry/RPCGeometry/interface/RPCRoll.h" +#include "Geometry/RPCGeometry/interface/RPCRollSpecs.h" +#include "SimMuon/RPCDigitizer/src/IRPCSimModelTiming.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" + +#include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h" +#include "Geometry/CommonTopologies/interface/RectangularStripTopology.h" +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" +#include "Geometry/RPCGeometry/interface/RPCGeomServ.h" + +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" +#include "DataFormats/RPCDigi/interface/IRPCDigiCollection.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CLHEP/Random/RandFlat.h" +#include "CLHEP/Random/RandPoissonQ.h" +#include "CLHEP/Random/RandGaussQ.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +IRPCSimModelTiming::IRPCSimModelTiming(const edm::ParameterSet& config) : RPCSim(config) { + aveEff = config.getParameter("averageEfficiency"); + aveCls = config.getParameter("averageClusterSize"); + resRPC = config.getParameter("timeResolution"); + timOff = config.getParameter("timingRPCOffset"); + dtimCs = config.getParameter("deltatimeAdjacentStrip"); + resEle = config.getParameter("timeJitter"); + sspeed = config.getParameter("signalPropagationSpeed"); + lbGate = config.getParameter("linkGateWidth"); + rpcdigiprint = config.getParameter("printOutDigitizer"); + + rate = config.getParameter("Rate"); + nbxing = config.getParameter("Nbxing"); + gate = config.getParameter("Gate"); + frate = config.getParameter("Frate"); + sigmaY = config.getParameter("sigmaY"); + eledig = config.getParameter("digitizeElectrons"); + + if (rpcdigiprint) { + edm::LogInfo("RPC digitizer parameters") << "Average Efficiency = " << aveEff; + edm::LogInfo("RPC digitizer parameters") << "Average Cluster Size = " << aveCls << " strips"; + edm::LogInfo("RPC digitizer parameters") << "RPC Time Resolution = " << resRPC << " ns"; + edm::LogInfo("RPC digitizer parameters") << "RPC Signal formation time = " << timOff << " ns"; + edm::LogInfo("RPC digitizer parameters") << "RPC adjacent strip delay = " << dtimCs << " ns"; + edm::LogInfo("RPC digitizer parameters") << "Electronic Jitter = " << resEle << " ns"; + edm::LogInfo("RPC digitizer parameters") << "Signal propagation time = " << sspeed << " x c"; + edm::LogInfo("RPC digitizer parameters") << "Link Board Gate Width = " << lbGate << " ns"; + } + + _rpcSync = new RPCSynchronizer(config); +} + +IRPCSimModelTiming::~IRPCSimModelTiming() { delete _rpcSync; } + +void IRPCSimModelTiming::simulate(const RPCRoll* roll, + const edm::PSimHitContainer& rpcHits, + CLHEP::HepRandomEngine* engine) { + _rpcSync->setRPCSimSetUp(getRPCSimSetUp()); + theRpcDigiSimLinks.clear(); + theDetectorHitMap.clear(); + theRpcDigiSimLinks = RPCDigiSimLinks(roll->id().rawId()); + + RPCDetId rpcId = roll->id(); + RPCGeomServ RPCname(rpcId); + + const Topology& topology = roll->specs()->topology(); + + for (edm::PSimHitContainer::const_iterator _hit = rpcHits.begin(); _hit != rpcHits.end(); ++_hit) { + if (!eledig && _hit->particleType() == 11) + continue; + + // Here I would check if the RPC are up side down; + const LocalPoint& entr = _hit->entryPoint(); + const TrapezoidalStripTopology* top_ = dynamic_cast(&(roll->topology())); + float striplength = (top_->stripLength()); + + std::pair TDCs = + _rpcSync->getDoubleTiming(&(*_hit), engine, striplength); // get the timing for the two TDCs + std::tuple tdc1 = + _rpcSync->getBX_SBX_fine_time(TDCs.first); //calculates the BX,subBX and fine_time for the first TDC + std::tuple tdc2 = + _rpcSync->getBX_SBX_fine_time(TDCs.second); //calculates the BX,subBX and fine_time for the second TDC + + float posX = roll->strip(_hit->localPosition()) - static_cast(roll->strip(_hit->localPosition())); + + std::vector veff = (getRPCSimSetUp())->getEff(rpcId.rawId()); + + // Effinciecy + int centralStrip = topology.channel(entr) + 1; + + float fire = CLHEP::RandFlat::shoot(engine); + + if (fire < veff[centralStrip - 1]) { + int fstrip = centralStrip; + int lstrip = centralStrip; + + // Compute the cluster size + int clsize = this->getClSize(rpcId.rawId(), posX, engine); // This is for cluster size chamber by chamber + std::vector cls; + cls.push_back(centralStrip); + if (clsize > 1) { + for (int cl = 0; cl < (clsize - 1) / 2; cl++) { + if (centralStrip - cl - 1 >= 1) { + fstrip = centralStrip - cl - 1; + cls.push_back(fstrip); + } + if (centralStrip + cl + 1 <= roll->nstrips()) { + lstrip = centralStrip + cl + 1; + cls.push_back(lstrip); + } + } + if (clsize % 2 == 0) { + // insert the last strip according to the + // simhit position in the central strip + int lr = LeftRightNeighbour(*roll, entr, centralStrip); + if (lr == 1) { + if (lstrip < roll->nstrips()) { + lstrip++; + cls.push_back(lstrip); + } + } else { + if (fstrip > 1) { + fstrip--; + cls.push_back(fstrip); + } + } + } + } + //digitize all the strips in the cluster + //in the previuos version some strips were dropped + //leading to un-physical "shift" of the cluster + for (std::vector::iterator i = cls.begin(); i != cls.end(); i++) { + std::pair digi(*i, std::get<0>(tdc1)); + IRPCDigi adigi(*i, + std::get<0>(tdc1), + std::get<0>(tdc2), + std::get<1>(tdc1), + std::get<1>(tdc2), + std::get<2>(tdc1), + std::get<2>(tdc2)); + irpc_digis.insert(adigi); + theDetectorHitMap.insert(DetectorHitMap::value_type(digi, &(*_hit))); + } + } + } +} + +void IRPCSimModelTiming::simulateNoise(const RPCRoll* roll, CLHEP::HepRandomEngine* engine) { + RPCDetId rpcId = roll->id(); + RPCGeomServ RPCname(rpcId); + std::vector vnoise = (getRPCSimSetUp())->getNoise(rpcId.rawId()); + std::vector veff = (getRPCSimSetUp())->getEff(rpcId.rawId()); + unsigned int nstrips = roll->nstrips(); + double area = 0.0; + float striplength, xmin, xmax; + if (rpcId.region() == 0) { + const RectangularStripTopology* top_ = dynamic_cast(&(roll->topology())); + xmin = (top_->localPosition(0.)).x(); + xmax = (top_->localPosition((float)roll->nstrips())).x(); + striplength = (top_->stripLength()); + area = striplength * (xmax - xmin); + } else { + const TrapezoidalStripTopology* top_ = dynamic_cast(&(roll->topology())); + xmin = (top_->localPosition(0.)).x(); + xmax = (top_->localPosition((float)roll->nstrips())).x(); + striplength = (top_->stripLength()); + area = striplength * (xmax - xmin); + } + + for (unsigned int j = 0; j < vnoise.size(); ++j) { + if (j >= nstrips) + break; + + double ave = vnoise[j] * nbxing * gate * area * 1.0e-9 * frate / ((float)roll->nstrips()); + + CLHEP::RandPoissonQ randPoissonQ(*engine, ave); + N_hits = randPoissonQ.fire(); + for (int i = 0; i < N_hits; i++) { + double TDC1_time = CLHEP::RandFlat::shoot(engine, (nbxing * gate) / gate); + int TDC1_BX = (static_cast(TDC1_time)) - nbxing / 2; + double TDC2_time = CLHEP::RandFlat::shoot(engine, (nbxing * gate) / gate); + int TDC2_BX = (static_cast(TDC2_time)) - nbxing / 2; + int TDC1_SBX = CLHEP::RandFlat::shootInt(long(0), long(10)); + int TDC2_SBX = CLHEP::RandFlat::shootInt(long(0), long(10)); + int TDC1_fine = CLHEP::RandFlat::shootInt(long(0), long(12)); + int TDC2_fine = CLHEP::RandFlat::shootInt(long(0), long(12)); + + IRPCDigi adigi(j + 1, TDC1_BX, TDC2_BX, TDC1_SBX, TDC2_SBX, TDC1_fine, TDC2_fine); + irpc_digis.insert(adigi); + } + } +} + +int IRPCSimModelTiming::getClSize(uint32_t id, float posX, CLHEP::HepRandomEngine* engine) { + std::vector clsForDetId = getRPCSimSetUp()->getCls(id); + + int cnt = 1; + int min = 1; + double func = 0.0; + std::vector sum_clsize; + + sum_clsize.clear(); + sum_clsize = clsForDetId; + int vectOffset(0); + + double rr_cl = CLHEP::RandFlat::shoot(engine); + + if (0.0 <= posX && posX < 0.2) { + func = clsForDetId[19] * (rr_cl); + vectOffset = 0; + } + if (0.2 <= posX && posX < 0.4) { + func = clsForDetId[39] * (rr_cl); + vectOffset = 20; + } + if (0.4 <= posX && posX < 0.6) { + func = clsForDetId[59] * (rr_cl); + vectOffset = 40; + } + if (0.6 <= posX && posX < 0.8) { + func = clsForDetId[79] * (rr_cl); + vectOffset = 60; + } + if (0.8 <= posX && posX < 1.0) { + func = clsForDetId[89] * (rr_cl); + vectOffset = 80; + } + + for (int i = vectOffset; i < (vectOffset + 20); i++) { + cnt++; + if (func > clsForDetId[i]) { + min = cnt; + } else if (func < clsForDetId[i]) { + break; + } + } + return min; +} + +int IRPCSimModelTiming::LeftRightNeighbour(const RPCRoll& roll, const LocalPoint& hit_pos, int strip) { + //if left return -1 + //if right return +1 + + int leftStrip = strip - 1; + int rightStrip = strip + 1; + + if (leftStrip < 0) + return +1; + if (rightStrip > roll.nstrips()) + return -1; + + double deltawL = fabs((roll.centreOfStrip(leftStrip)).x() - hit_pos.x()); + double deltawR = fabs((roll.centreOfStrip(rightStrip)).x() - hit_pos.x()); + + if (deltawL >= deltawR) { + return +1; + } else { + return -1; + } +} diff --git a/SimMuon/RPCDigitizer/src/IRPCSimModelTiming.h b/SimMuon/RPCDigitizer/src/IRPCSimModelTiming.h new file mode 100644 index 0000000000000..fa344c890968f --- /dev/null +++ b/SimMuon/RPCDigitizer/src/IRPCSimModelTiming.h @@ -0,0 +1,69 @@ +#ifndef RPCDigitizer_IRPCSimModelTiming_h +#define RPCDigitizer_IRPCSimModelTiming_h + +/** \class RPCSimAverage + * Class for the RPC strip response simulation based + * on a parametrized model (ORCA-based) + * + * \author Borislav Pavlov -- University of Sofia + */ +#include "SimMuon/RPCDigitizer/src/RPCSim.h" +#include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h" +#include "SimMuon/RPCDigitizer/src/RPCSimAsymmetricCls.h" +#include "SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEffCls.h" + +#include +#include +#include +#include +#include +#include +#include "FWCore/Framework/interface/EventSetup.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +class RPCGeometry; + +namespace CLHEP { + class HepRandomEngine; +} + +class IRPCSimModelTiming : public RPCSim { +public: + IRPCSimModelTiming(const edm::ParameterSet& config); + ~IRPCSimModelTiming() override; + void simulate(const RPCRoll* roll, const edm::PSimHitContainer& rpcHits, CLHEP::HepRandomEngine*) override; + void simulateNoise(const RPCRoll*, CLHEP::HepRandomEngine*) override; + int getClSize(uint32_t id, float posX, CLHEP::HepRandomEngine*); + int LeftRightNeighbour(const RPCRoll& roll, const LocalPoint& hit_pos, int strip); + +protected: + void init() override {} + + double aveEff; + double aveCls; + double resRPC; + double timOff; + double dtimCs; + double resEle; + double sspeed; + double lbGate; + bool rpcdigiprint; + bool eledig; + + int N_hits; + int nbxing; + double rate; + double gate; + double frate; + bool do_Y; + double sigmaY; + + std::map > clsMap; + std::vector sum_clsize; + std::vector clsForDetId; + std::ifstream* infile; + + RPCSynchronizer* _rpcSync; +}; +#endif diff --git a/SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.cc b/SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.cc new file mode 100644 index 0000000000000..0e9647c003718 --- /dev/null +++ b/SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.cc @@ -0,0 +1,109 @@ +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" +#include "SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.h" +#include "SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h" +#include +#include + +#include +#include + +#include "FWCore/Framework/interface/MakerMacros.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" + +//Random Number +#include "FWCore/AbstractServices/interface/RandomNumberGenerator.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "CLHEP/Random/RandFlat.h" + +namespace CLHEP { + class HepRandomEngine; +} + +RPCDigiPhase2Producer::RPCDigiPhase2Producer(const edm::ParameterSet& ps) { + produces(); + produces("RPCDigiPhase2SimLink"); + + //Name of Collection used for create the XF + mix_ = ps.getParameter("mixLabel"); + collection_for_XF = ps.getParameter("InputCollection"); + + edm::Service rng; + if (!rng.isAvailable()) { + throw cms::Exception("Configuration") + << "RPCDigitizerPhase2 requires the RandomNumberGeneratorService\n" + "which is not present in the configuration file. You must add the service\n" + "in the configuration file or remove the modules that require it."; + }; + + theRPCSimSetUpRPC = new RPCSimSetUp(ps); + theRPCDigitizerPhase2 = new RPCDigitizerPhase2(ps); + crossingFrameToken = consumes>(edm::InputTag(mix_, collection_for_XF)); + geomToken = esConsumes(); + noiseToken = esConsumes(); + clsToken = esConsumes(); +} + +RPCDigiPhase2Producer::~RPCDigiPhase2Producer() { + delete theRPCDigitizerPhase2; + delete theRPCSimSetUpRPC; +} + +void RPCDigiPhase2Producer::beginRun(const edm::Run& r, const edm::EventSetup& eventSetup) { + edm::ESHandle hGeom = eventSetup.getHandle(geomToken); + const RPCGeometry* pGeom = &*hGeom; + _pGeom = &*hGeom; + edm::ESHandle noiseRcd = eventSetup.getHandle(noiseToken); + + edm::ESHandle clsRcd = eventSetup.getHandle(clsToken); + + //setup the two digi models + theRPCSimSetUpRPC->setGeometry(pGeom); + theRPCSimSetUpRPC->setRPCSetUp(noiseRcd->getVNoise(), clsRcd->getCls()); + + //setup the two digitizers + theRPCDigitizerPhase2->setGeometry(pGeom); + theRPCDigitizerPhase2->setRPCSimSetUp(theRPCSimSetUpRPC); +} + +void RPCDigiPhase2Producer::produce(edm::Event& e, const edm::EventSetup& eventSetup) { + edm::Service rng; + CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID()); + + LogDebug("RPCDigiPhase2Producer") << "[RPCDigiPhase2Producer::produce] got the CLHEP::HepRandomEngine engine from " + "the edm::Event.streamID() and edm::Service"; + LogDebug("RPCDigiPhase2Producer") << "[RPCDigiPhase2Producer::produce] test the CLHEP::HepRandomEngine by firing " + "once RandFlat ---- this must be the first time in SimMuon/RPCDigitizerPhase2"; + LogDebug("RPCDigiPhase2Producer") << "[RPCDigiPhase2Producer::produce] to activate the test go in " + "RPCDigiPhase2Producer.cc and uncomment the line below"; + + edm::Handle> cf; + e.getByToken(crossingFrameToken, cf); + + std::unique_ptr> hits(new MixCollection(cf.product())); + + // Create empty output + std::unique_ptr pDigis(new RPCDigiPhase2Collection()); + std::unique_ptr RPCDigitSimLink(new RPCDigitizerPhase2SimLinks()); + + theRPCDigitizerPhase2->doAction( + *hits, *pDigis, *RPCDigitSimLink, engine); //make "bakelite RPC" digitizer do the action + + e.put(std::move(pDigis)); + //store the SimDigiLinks in the event + e.put(std::move(RPCDigitSimLink), "RPCDigiPhase2SimLink"); +} diff --git a/SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.h b/SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.h new file mode 100644 index 0000000000000..4da91a67e1aba --- /dev/null +++ b/SimMuon/RPCDigitizer/src/RPCDigiPhase2Producer.h @@ -0,0 +1,56 @@ +#ifndef RPCDigiPhase2Producer_h +#define RPCDigiPhase2Producer_h + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/ESGetToken.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "CondFormats/RPCObjects/interface/RPCStripNoises.h" +#include "CondFormats/DataRecord/interface/RPCStripNoisesRcd.h" +#include "CondFormats/RPCObjects/interface/RPCClusterSize.h" +#include "CondFormats/DataRecord/interface/RPCClusterSizeRcd.h" + +class RPCGeometry; +class RPCSimSetUp; +class RPCSynchronizer; + +class RPCDigiPhase2Producer : public edm::stream::EDProducer { +public: + // typedef edm::DetSetVector RPCDigiSimLinks; + typedef RPCDigitizerPhase2::RPCDigiSimLinks RPCDigitizerPhase2SimLinks; + + explicit RPCDigiPhase2Producer(const edm::ParameterSet& ps); + ~RPCDigiPhase2Producer() override; + + void beginRun(const edm::Run&, const edm::EventSetup&) override; + + /**Produces the EDM products,*/ + void produce(edm::Event& e, const edm::EventSetup& c) override; + + void setRPCSetUp(const std::vector&, const std::vector&); + +private: + RPCDigitizerPhase2* theRPCDigitizerPhase2; + RPCSimSetUp* theRPCSimSetUpRPC; + + //Name of Collection used for create the XF + std::string mix_; + std::string collection_for_XF; + + //Token for accessing data + edm::EDGetTokenT> crossingFrameToken; + const RPCGeometry* _pGeom; + + //EventSetup Tokens + edm::ESGetToken geomToken; + edm::ESGetToken noiseToken; + edm::ESGetToken clsToken; +}; + +#endif diff --git a/SimMuon/RPCDigitizer/src/RPCDigiProducer.cc b/SimMuon/RPCDigitizer/src/RPCDigiProducer.cc index b322b8733d5fb..f14c61bcfdecf 100644 --- a/SimMuon/RPCDigitizer/src/RPCDigiProducer.cc +++ b/SimMuon/RPCDigitizer/src/RPCDigiProducer.cc @@ -48,7 +48,7 @@ RPCDigiProducer::RPCDigiProducer(const edm::ParameterSet& ps) { "in the configuration file or remove the modules that require it."; } theRPCSimSetUp = new RPCSimSetUp(ps); - theDigitizer = new RPCDigitizer(ps); + theDigitizer = new RPCDigitizer(ps, true); crossingFrameToken = consumes>(edm::InputTag(mix_, collection_for_XF)); geomToken = esConsumes(); noiseToken = esConsumes(); diff --git a/SimMuon/RPCDigitizer/src/RPCDigitizer.cc b/SimMuon/RPCDigitizer/src/RPCDigitizer.cc index fb5c5c0ebfd1f..464429ae1b215 100644 --- a/SimMuon/RPCDigitizer/src/RPCDigitizer.cc +++ b/SimMuon/RPCDigitizer/src/RPCDigitizer.cc @@ -9,10 +9,11 @@ // default constructor allocates default wire and strip digitizers -RPCDigitizer::RPCDigitizer(const edm::ParameterSet& config) +RPCDigitizer::RPCDigitizer(const edm::ParameterSet& config, bool type) : theRPCSim{RPCSimFactory::get()->create(config.getParameter("digiModel"), config.getParameter("digiModelConfig"))}, - theNoise{config.getParameter("doBkgNoise")} {} + theNoise{config.getParameter("doBkgNoise")}, + theType(type) {} RPCDigitizer::~RPCDigitizer() = default; @@ -39,7 +40,7 @@ void RPCDigitizer::doAction(MixCollection& simHits, RPCDetId id = (*r)->id(); const edm::PSimHitContainer& rollSimHits = hitMap[id]; - if (!((*r)->isIRPC())) { + if ((*r)->isIRPC() != theType) { theRPCSim->simulate(*r, rollSimHits, engine); if (theNoise) { @@ -48,7 +49,9 @@ void RPCDigitizer::doAction(MixCollection& simHits, } theRPCSim->fillDigis((*r)->id(), rpcDigis); - rpcDigiSimLink.insert(theRPCSim->rpcDigiSimLinks()); + if (rpcDigiSimLink.find((theRPCSim->rpcDigiSimLinks()).detId()) == rpcDigiSimLink.end()) { + rpcDigiSimLink.insert(theRPCSim->rpcDigiSimLinks()); + } } } diff --git a/SimMuon/RPCDigitizer/src/RPCDigitizer.h b/SimMuon/RPCDigitizer/src/RPCDigitizer.h index eb5e3a5f4672d..1f221c00c8eca 100644 --- a/SimMuon/RPCDigitizer/src/RPCDigitizer.h +++ b/SimMuon/RPCDigitizer/src/RPCDigitizer.h @@ -33,7 +33,8 @@ namespace CLHEP { class RPCDigitizer { public: typedef edm::DetSetVector RPCDigiSimLinks; - RPCDigitizer(const edm::ParameterSet& config); + RPCDigitizer(const edm::ParameterSet& config, bool type); + ~RPCDigitizer(); // *** digitize *** @@ -57,6 +58,7 @@ class RPCDigitizer { std::unique_ptr theRPCSim; RPCSimSetUp* theSimSetUp; bool theNoise; + bool theType; // true for RPC, false for iRPC }; #endif diff --git a/SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.cc b/SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.cc new file mode 100644 index 0000000000000..9be623ab9a1b5 --- /dev/null +++ b/SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.cc @@ -0,0 +1,62 @@ +#include "SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.h" +#include "SimMuon/RPCDigitizer/src/RPCSimFactory.h" +#include "SimMuon/RPCDigitizer/src/RPCSim.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "Geometry/RPCGeometry/interface/RPCRoll.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" + +// default constructor allocates default wire and strip digitizers + +RPCDigitizerPhase2::RPCDigitizerPhase2(const edm::ParameterSet& config) + : theRPCSim{RPCSimFactory::get()->create(config.getParameter("digiModel"), + config.getParameter("digiModelConfig"))}, + theNoise{config.getParameter("doBkgNoise")} {} + +RPCDigitizerPhase2::~RPCDigitizerPhase2() = default; + +void RPCDigitizerPhase2::doAction(MixCollection& simHits, + RPCDigiPhase2Collection& rpcDigis, + RPCDigiSimLinks& rpcDigiSimLink, + CLHEP::HepRandomEngine* engine) { + theRPCSim->setRPCSimSetUp(theSimSetUp); + + // arrange the hits by roll + std::map hitMap; + for (MixCollection::MixItr hitItr = simHits.begin(); hitItr != simHits.end(); ++hitItr) { + hitMap[hitItr->detUnitId()].push_back(*hitItr); + } + + if (!theGeometry) { + throw cms::Exception("Configuration") + << "RPCDigitizerPhase2 requires the RPCGeometry \n which is not present in the configuration file. You must " + "add the " + "service\n in the configuration file or remove the modules that require it."; + } + + const std::vector& rpcRolls = theGeometry->rolls(); + for (auto r = rpcRolls.begin(); r != rpcRolls.end(); r++) { + RPCDetId id = (*r)->id(); + const edm::PSimHitContainer& rollSimHits = hitMap[id]; + + if (!((*r)->isIRPC())) { + theRPCSim->simulate(*r, rollSimHits, engine); + + if (theNoise) { + theRPCSim->simulateNoise(*r, engine); + } + } + + theRPCSim->fillDigis((*r)->id(), rpcDigis); + if (rpcDigiSimLink.find((theRPCSim->rpcDigiSimLinks()).detId()) == rpcDigiSimLink.end()) { + rpcDigiSimLink.insert(theRPCSim->rpcDigiSimLinks()); + } + } +} + +const RPCRoll* RPCDigitizerPhase2::findDet(int detId) const { + assert(theGeometry != nullptr); + const GeomDetUnit* detUnit = theGeometry->idToDetUnit(RPCDetId(detId)); + return dynamic_cast(detUnit); +} diff --git a/SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.h b/SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.h new file mode 100644 index 0000000000000..6810c87020563 --- /dev/null +++ b/SimMuon/RPCDigitizer/src/RPCDigitizerPhase2.h @@ -0,0 +1,62 @@ +#ifndef SimMuon_RPCDigitizerPhase2_h +#define SimMuon_RPCDigitizerPhase2_h +// + +/** \class RPCDigitizerPhase2 + * Digitizer class for RPC Phase2 upgrade + * + * \author Borislav Pavlov -- Sofia University + * + */ +#include "DataFormats/Common/interface/DetSetVector.h" +#include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h" +#include "SimDataFormats/RPCDigiSimLink/interface/RPCDigiSimLink.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "DataFormats/RPCDigi/interface/RPCDigiPhase2Collection.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include +#include "CLHEP/Random/RandomEngine.h" + +namespace edm { + class ParameterSet; +} + +class RPCRoll; +class RPCSim; +class RPCSimSetUp; + +namespace CLHEP { + class HepRandomEngine; +} + +class RPCDigitizerPhase2 { +public: + typedef edm::DetSetVector RPCDigiSimLinks; + RPCDigitizerPhase2(const edm::ParameterSet& config); + ~RPCDigitizerPhase2(); + + // *** digitize *** + void doAction(MixCollection& simHits, + RPCDigiPhase2Collection& rpcDigis, + RPCDigiSimLinks& rpcDigiSimLink, + CLHEP::HepRandomEngine*); + + /// sets geometry + void setGeometry(const RPCGeometry* geom) { theGeometry = geom; } + + void setRPCSimSetUp(RPCSimSetUp* simsetup) { theSimSetUp = simsetup; } + + RPCSimSetUp* getRPCSimSetUp() { return theSimSetUp; } + + /// finds the rpc det unit in the geometry associated with this det ID + const RPCRoll* findDet(int detId) const; + +private: + const RPCGeometry* theGeometry; + std::unique_ptr theRPCSim; + RPCSimSetUp* theSimSetUp; + bool theNoise; +}; + +#endif diff --git a/SimMuon/RPCDigitizer/src/RPCSim.cc b/SimMuon/RPCDigitizer/src/RPCSim.cc index 79a2c8d0d6707..07637a96b75ae 100644 --- a/SimMuon/RPCDigitizer/src/RPCSim.cc +++ b/SimMuon/RPCDigitizer/src/RPCSim.cc @@ -21,6 +21,16 @@ void RPCSim::fillDigis(int rollDetId, RPCDigiCollection& digis) { } strips.clear(); + for (auto it : rpc_digis) { + if (it.bx() != -999) { + digis.insertDigi(RPCDetId(rollDetId), it); + this->addLinks(it.strip(), it.bx()); + } + } + rpc_digis.clear(); +} + +void RPCSim::fillDigis(int rollDetId, IRPCDigiCollection& digis) { for (auto it : irpc_digis) { if (it.bx() != -999) { digis.insertDigi(RPCDetId(rollDetId), it); @@ -30,6 +40,16 @@ void RPCSim::fillDigis(int rollDetId, RPCDigiCollection& digis) { irpc_digis.clear(); } +void RPCSim::fillDigis(int rollDetId, RPCDigiPhase2Collection& digis) { + for (auto it : rpc_digis_phase2) { + if (it.bx() != -999) { + digis.insertDigi(RPCDetId(rollDetId), it); + this->addLinks(it.strip(), it.bx()); + } + } + rpc_digis_phase2.clear(); +} + void RPCSim::addLinks(unsigned int strip, int bx) { std::pair digi(strip, bx); std::pair channelHitItr = theDetectorHitMap.equal_range(digi); diff --git a/SimMuon/RPCDigitizer/src/RPCSim.h b/SimMuon/RPCDigitizer/src/RPCSim.h index ea23a39562c7d..20d3fd449f732 100644 --- a/SimMuon/RPCDigitizer/src/RPCSim.h +++ b/SimMuon/RPCDigitizer/src/RPCSim.h @@ -8,6 +8,8 @@ */ #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" #include "DataFormats/RPCDigi/interface/RPCDigiCollection.h" +#include "DataFormats/RPCDigi/interface/RPCDigiPhase2Collection.h" +#include "DataFormats/RPCDigi/interface/IRPCDigiCollection.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Framework/interface/EventSetup.h" @@ -40,6 +42,10 @@ class RPCSim { virtual void fillDigis(int rollDetId, RPCDigiCollection& digis); + virtual void fillDigis(int rollDetId, IRPCDigiCollection& digis); + + virtual void fillDigis(int rollDetId, RPCDigiPhase2Collection& digis); + void setRPCSimSetUp(RPCSimSetUp* setup) { theSimSetUp = setup; } RPCSimSetUp* getRPCSimSetUp() { return theSimSetUp; } @@ -53,7 +59,9 @@ class RPCSim { protected: std::set > strips; - std::set irpc_digis; + std::set rpc_digis; + std::set rpc_digis_phase2; + std::set irpc_digis; //--------NEW--------------------- diff --git a/SimMuon/RPCDigitizer/src/RPCSimModelTiming.cc b/SimMuon/RPCDigitizer/src/RPCSimModelTiming.cc index 918aedf1c05a0..3dc9aef5f8e13 100644 --- a/SimMuon/RPCDigitizer/src/RPCSimModelTiming.cc +++ b/SimMuon/RPCDigitizer/src/RPCSimModelTiming.cc @@ -155,7 +155,7 @@ void RPCSimModelTiming::simulate(const RPCRoll* roll, adigi.setY(smearedPositionY); adigi.setDeltaY(sigmaY); } - irpc_digis.insert(adigi); + rpc_digis.insert(adigi); theDetectorHitMap.insert(DetectorHitMap::value_type(digi, &(*_hit))); } } @@ -205,7 +205,7 @@ void RPCSimModelTiming::simulateNoise(const RPCRoll* roll, CLHEP::HepRandomEngin adigi.setY(positionY); adigi.setDeltaY(sigmaY); } - irpc_digis.insert(adigi); + rpc_digis.insert(adigi); } } } diff --git a/SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.cc b/SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.cc new file mode 100644 index 0000000000000..9c453bde1611c --- /dev/null +++ b/SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.cc @@ -0,0 +1,267 @@ +#include "Geometry/RPCGeometry/interface/RPCRoll.h" +#include "Geometry/RPCGeometry/interface/RPCRollSpecs.h" +#include "SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" + +#include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h" +#include "Geometry/CommonTopologies/interface/RectangularStripTopology.h" +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" +#include "Geometry/RPCGeometry/interface/RPCGeomServ.h" + +#include + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CLHEP/Random/RandFlat.h" +#include "CLHEP/Random/RandPoissonQ.h" +#include "CLHEP/Random/RandGaussQ.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +RPCSimModelTimingPhase2::RPCSimModelTimingPhase2(const edm::ParameterSet& config) : RPCSim(config) { + aveEff = config.getParameter("averageEfficiency"); + aveCls = config.getParameter("averageClusterSize"); + resRPC = config.getParameter("timeResolution"); + timOff = config.getParameter("timingRPCOffset"); + dtimCs = config.getParameter("deltatimeAdjacentStrip"); + resEle = config.getParameter("timeJitter"); + sspeed = config.getParameter("signalPropagationSpeed"); + lbGate = config.getParameter("linkGateWidth"); + rpcdigiprint = config.getParameter("printOutDigitizer"); + + rate = config.getParameter("Rate"); + nbxing = config.getParameter("Nbxing"); + gate = config.getParameter("Gate"); + frate = config.getParameter("Frate"); + eledig = config.getParameter("digitizeElectrons"); + + if (rpcdigiprint) { + edm::LogInfo("RPC digitizer parameters") << "Average Efficiency = " << aveEff; + edm::LogInfo("RPC digitizer parameters") << "Average Cluster Size = " << aveCls << " strips"; + edm::LogInfo("RPC digitizer parameters") << "RPC Time Resolution = " << resRPC << " ns"; + edm::LogInfo("RPC digitizer parameters") << "RPC Signal formation time = " << timOff << " ns"; + edm::LogInfo("RPC digitizer parameters") << "RPC adjacent strip delay = " << dtimCs << " ns"; + edm::LogInfo("RPC digitizer parameters") << "Electronic Jitter = " << resEle << " ns"; + edm::LogInfo("RPC digitizer parameters") << "Signal propagation time = " << sspeed << " x c"; + edm::LogInfo("RPC digitizer parameters") << "Link Board Gate Width = " << lbGate << " ns"; + } + + _rpcSync = new RPCSynchronizer(config); +} + +RPCSimModelTimingPhase2::~RPCSimModelTimingPhase2() { delete _rpcSync; } + +void RPCSimModelTimingPhase2::simulate(const RPCRoll* roll, + const edm::PSimHitContainer& rpcHits, + CLHEP::HepRandomEngine* engine) { + _rpcSync->setRPCSimSetUp(getRPCSimSetUp()); + theRpcDigiSimLinks.clear(); + theDetectorHitMap.clear(); + theRpcDigiSimLinks = RPCDigiSimLinks(roll->id().rawId()); + + RPCDetId rpcId = roll->id(); + RPCGeomServ RPCname(rpcId); + + const Topology& topology = roll->specs()->topology(); + + for (edm::PSimHitContainer::const_iterator _hit = rpcHits.begin(); _hit != rpcHits.end(); ++_hit) { + if (!eledig && _hit->particleType() == 11) + continue; + // Here I hould check if the RPC are up side down; + const LocalPoint& entr = _hit->entryPoint(); + + float striplength; + if (rpcId.region() == 0) { + const RectangularStripTopology* top_ = dynamic_cast(&(roll->topology())); + striplength = (top_->stripLength()); + + } else { + const TrapezoidalStripTopology* top_ = dynamic_cast(&(roll->topology())); + striplength = (top_->stripLength()); + } + + double precise_time = _rpcSync->getTiming(&(*_hit), engine, striplength); + std::pair tdc = _rpcSync->getBX_SBX(precise_time); + float posX = roll->strip(_hit->localPosition()) - static_cast(roll->strip(_hit->localPosition())); + + std::vector veff = (getRPCSimSetUp())->getEff(rpcId.rawId()); + + // Effinciecy + int centralStrip = topology.channel(entr) + 1; + float fire = CLHEP::RandFlat::shoot(engine); + + if (fire < veff[centralStrip - 1]) { + int fstrip = centralStrip; + int lstrip = centralStrip; + + // Compute the cluster size + int clsize = this->getClSize(rpcId.rawId(), posX, engine); // This is for cluster size chamber by chamber + std::vector cls; + cls.push_back(centralStrip); + if (clsize > 1) { + for (int cl = 0; cl < (clsize - 1) / 2; cl++) { + if (centralStrip - cl - 1 >= 1) { + fstrip = centralStrip - cl - 1; + cls.push_back(fstrip); + } + if (centralStrip + cl + 1 <= roll->nstrips()) { + lstrip = centralStrip + cl + 1; + cls.push_back(lstrip); + } + } + if (clsize % 2 == 0) { + // insert the last strip according to the + // simhit position in the central strip + int lr = LeftRightNeighbour(*roll, entr, centralStrip); + if (lr == 1) { + if (lstrip < roll->nstrips()) { + lstrip++; + cls.push_back(lstrip); + } + } else { + if (fstrip > 1) { + fstrip--; + cls.push_back(fstrip); + } + } + } + } + //digitize all the strips in the cluster + //in the previuos version some strips were dropped + //leading to un-physical "shift" of the cluster + for (std::vector::iterator i = cls.begin(); i != cls.end(); i++) { + std::pair digi(*i, tdc.first); + RPCDigiPhase2 adigi(*i, tdc.first, tdc.second); + rpc_digis_phase2.insert(adigi); + theDetectorHitMap.insert(DetectorHitMap::value_type(digi, &(*_hit))); + } + } + } +} + +void RPCSimModelTimingPhase2::simulateNoise(const RPCRoll* roll, CLHEP::HepRandomEngine* engine) { + RPCDetId rpcId = roll->id(); + RPCGeomServ RPCname(rpcId); + std::vector vnoise = (getRPCSimSetUp())->getNoise(rpcId.rawId()); + std::vector veff = (getRPCSimSetUp())->getEff(rpcId.rawId()); + unsigned int nstrips = roll->nstrips(); + double area = 0.0; + float striplength, xmin, xmax; + if (rpcId.region() == 0) { + const RectangularStripTopology* top_ = dynamic_cast(&(roll->topology())); + xmin = (top_->localPosition(0.)).x(); + xmax = (top_->localPosition((float)roll->nstrips())).x(); + striplength = (top_->stripLength()); + area = striplength * (xmax - xmin); + } else { + const TrapezoidalStripTopology* top_ = dynamic_cast(&(roll->topology())); + xmin = (top_->localPosition(0.)).x(); + xmax = (top_->localPosition((float)roll->nstrips())).x(); + striplength = (top_->stripLength()); + area = striplength * (xmax - xmin); + } + + for (unsigned int j = 0; j < vnoise.size(); ++j) { + if (j >= nstrips) + break; + + double ave = vnoise[j] * nbxing * gate * area * 1.0e-9 * frate / ((float)roll->nstrips()); + + CLHEP::RandPoissonQ randPoissonQ(*engine, ave); + N_hits = randPoissonQ.fire(); + for (int i = 0; i < N_hits; i++) { + double precise_time = CLHEP::RandFlat::shoot(engine, (nbxing * gate) / gate); + int time_hit = (static_cast(precise_time)) - nbxing / 2; + int sbx = CLHEP::RandFlat::shootInt(long(0), long(10)); + RPCDigiPhase2 adigi(j + 1, time_hit, sbx); + rpc_digis_phase2.insert(adigi); + } + } +} + +int RPCSimModelTimingPhase2::getClSize(uint32_t id, float posX, CLHEP::HepRandomEngine* engine) { + std::vector clsForDetId = getRPCSimSetUp()->getCls(id); + + int cnt = 1; + int min = 1; + double func = 0.0; + std::vector sum_clsize; + + sum_clsize.clear(); + sum_clsize = clsForDetId; + int vectOffset(0); + + double rr_cl = CLHEP::RandFlat::shoot(engine); + + if (0.0 <= posX && posX < 0.2) { + func = clsForDetId[19] * (rr_cl); + vectOffset = 0; + } + if (0.2 <= posX && posX < 0.4) { + func = clsForDetId[39] * (rr_cl); + vectOffset = 20; + } + if (0.4 <= posX && posX < 0.6) { + func = clsForDetId[59] * (rr_cl); + vectOffset = 40; + } + if (0.6 <= posX && posX < 0.8) { + func = clsForDetId[79] * (rr_cl); + vectOffset = 60; + } + if (0.8 <= posX && posX < 1.0) { + func = clsForDetId[89] * (rr_cl); + vectOffset = 80; + } + + for (int i = vectOffset; i < (vectOffset + 20); i++) { + cnt++; + if (func > clsForDetId[i]) { + min = cnt; + } else if (func < clsForDetId[i]) { + break; + } + } + return min; +} + +int RPCSimModelTimingPhase2::LeftRightNeighbour(const RPCRoll& roll, const LocalPoint& hit_pos, int strip) { + //if left return -1 + //if right return +1 + + int leftStrip = strip - 1; + int rightStrip = strip + 1; + + if (leftStrip < 0) + return +1; + if (rightStrip > roll.nstrips()) + return -1; + + double deltawL = fabs((roll.centreOfStrip(leftStrip)).x() - hit_pos.x()); + double deltawR = fabs((roll.centreOfStrip(rightStrip)).x() - hit_pos.x()); + + if (deltawL >= deltawR) { + return +1; + } else { + return -1; + } +} diff --git a/SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.h b/SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.h new file mode 100644 index 0000000000000..21e3f619b6c4c --- /dev/null +++ b/SimMuon/RPCDigitizer/src/RPCSimModelTimingPhase2.h @@ -0,0 +1,69 @@ +#ifndef RPCDigitizer_RPCSimModelTimingPhase2_h +#define RPCDigitizer_RPCSimModelTimingPhase2_h + +/** \class RPCSimAverage + * Class for the RPC strip response simulation based + * on a parametrized model (ORCA-based) + * + * \author Borislav Pavlov -- University of Sofia + */ +#include "SimMuon/RPCDigitizer/src/RPCSim.h" +#include "SimMuon/RPCDigitizer/src/RPCSynchronizer.h" +#include "SimMuon/RPCDigitizer/src/RPCSimAsymmetricCls.h" +#include "SimMuon/RPCDigitizer/src/RPCSimAverageNoiseEffCls.h" + +#include +#include +#include +#include +#include +#include +#include "FWCore/Framework/interface/EventSetup.h" +#include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +class RPCGeometry; + +namespace CLHEP { + class HepRandomEngine; +} + +class RPCSimModelTimingPhase2 : public RPCSim { +public: + RPCSimModelTimingPhase2(const edm::ParameterSet& config); + ~RPCSimModelTimingPhase2() override; + void simulate(const RPCRoll* roll, const edm::PSimHitContainer& rpcHits, CLHEP::HepRandomEngine*) override; + void simulateNoise(const RPCRoll*, CLHEP::HepRandomEngine*) override; + int getClSize(uint32_t id, float posX, CLHEP::HepRandomEngine*); + int LeftRightNeighbour(const RPCRoll& roll, const LocalPoint& hit_pos, int strip); + +protected: + void init() override {} + + double aveEff; + double aveCls; + double resRPC; + double timOff; + double dtimCs; + double resEle; + double sspeed; + double lbGate; + bool rpcdigiprint; + bool eledig; + + int N_hits; + int nbxing; + double rate; + double gate; + double frate; + bool do_Y; + double sigmaY; + + std::map > clsMap; + std::vector sum_clsize; + std::vector clsForDetId; + std::ifstream* infile; + + RPCSynchronizer* _rpcSync; +}; +#endif diff --git a/SimMuon/RPCDigitizer/src/RPCSynchronizer.cc b/SimMuon/RPCDigitizer/src/RPCSynchronizer.cc index 8ba5088a7bd05..f564f15906f89 100644 --- a/SimMuon/RPCDigitizer/src/RPCSynchronizer.cc +++ b/SimMuon/RPCDigitizer/src/RPCSynchronizer.cc @@ -131,7 +131,6 @@ int RPCSynchronizer::getSimHitBx(const PSimHit* simhit, CLHEP::HepRandomEngine* } } } - return bx; } @@ -145,8 +144,6 @@ int RPCSynchronizer::getSimHitBxAndTimingForIRPC(const PSimHit* simhit, CLHEP::H float tof = simhit->timeOfFlight(); //automatic variable to prevent memory leak - - // float rr_el = CLHEP::RandGaussQ::shoot(engine, 0.,resEle); float rr_el = CLHEP::RandGaussQ::shoot(engine, 0., irpc_electronics_jitter); RPCDetId SimDetId(simhit->detUnitId()); @@ -183,10 +180,7 @@ int RPCSynchronizer::getSimHitBxAndTimingForIRPC(const PSimHit* simhit, CLHEP::H } float prop_time = distanceFromEdge / sspeed; - - // double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0.,resRPC); double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0., irpc_timing_res); - double total_time = tof + prop_time + timOff + rr_tim1 + rr_el; // Bunch crossing assignment @@ -221,15 +215,140 @@ int RPCSynchronizer::getSimHitBxAndTimingForIRPC(const PSimHit* simhit, CLHEP::H if (inf_time < time_differ && time_differ < sup_time) { bx = n; - the_exact_time = exact_time_differ; - the_smeared_time = time_differ; - - //cout<<"Debug\t"<getRPCSimSetUp(); + float timeref = simsetup->getTime(simhit->detUnitId()); + + LocalPoint simHitPos = simhit->localPosition(); + float tof = simhit->timeOfFlight(); + + //automatic variable to prevent memory leak + float rr_el = CLHEP::RandGaussQ::shoot(engine, 0., irpc_electronics_jitter); + + RPCDetId SimDetId(simhit->detUnitId()); + float distanceFromEdge = 0; + float half_stripL = StripLength / 2.; + + if (SimDetId.region() == 0) { + distanceFromEdge = half_stripL + simHitPos.y(); + } else { + distanceFromEdge = half_stripL - simHitPos.y(); + } + + float prop_time = distanceFromEdge / sspeed; + + // double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0.,resRPC); + double rr_tim1 = CLHEP::RandGaussQ::shoot(engine, 0., irpc_timing_res); + + double total_time = tof + prop_time + timOff + rr_tim1 + rr_el; + + // Bunch crossing assignment + double time_differ = 0.; + + if (cosmics) { + time_differ = (total_time - (timeref + ((half_stripL / sspeed) + timOff))) / cosmicPar; + } else if (!cosmics) { + time_differ = total_time - (timeref + (half_stripL / sspeed) + timOff); + } + + return time_differ; +} + +std::pair RPCSynchronizer::getDoubleTiming(const PSimHit* simhit, + CLHEP::HepRandomEngine* engine, + float StripLength) { + RPCSimSetUp* simsetup = this->getRPCSimSetUp(); + float timeref = simsetup->getTime(simhit->detUnitId()); + + LocalPoint simHitPos = simhit->localPosition(); + float tof = simhit->timeOfFlight(); + RPCDetId SimDetId(simhit->detUnitId()); + + double rpc_resolution = CLHEP::RandGaussQ::shoot(engine, 0., irpc_timing_res); + + //rpc_time simulate the iRPC timing resolution + double rpc_time = tof - timeref + rpc_resolution; + + //First FEB time resolution simulation + float feb_resolution = CLHEP::RandGaussQ::shoot(engine, 0., irpc_electronics_jitter); + + // The correct signal propagation is StripLength/2 + signalSign*simHitPos.y()/sspeed, but + // StripLength/2 sohould be substituted in order to have time = 0 when the particle hits the center of the RPC, + // so signal propagation is StripLength/2 + signalSign*simHitPos.y()/sspeed - StripLength/2 = signalSign*simHitPos.y()/sspeed + double tdc_LR_time = rpc_time + simHitPos.y() / sspeed + feb_resolution; + + //In the similar way for the second FEB TDC + feb_resolution = CLHEP::RandGaussQ::shoot(engine, 0., irpc_electronics_jitter); + double tdc_HR_time = rpc_time - simHitPos.y() / sspeed + feb_resolution; + + if (cosmics) { + tdc_LR_time /= cosmicPar; + tdc_HR_time /= cosmicPar; + } + + std::pair TDCs; + TDCs.first = tdc_LR_time; + TDCs.second = tdc_HR_time; + return TDCs; +} + +int RPCSynchronizer::getBX(float time) { + int bx = -999; + double inf_time = 0; + double sup_time = 0; + + for (int n = -N_BX; n <= N_BX; ++n) { + if (cosmics) { + inf_time = (-lbGate / 2 + n * LHCGate) / cosmicPar; + sup_time = (lbGate / 2 + n * LHCGate) / cosmicPar; + } else if (!cosmics) { + inf_time = -lbGate / 2 + n * LHCGate; + sup_time = lbGate / 2 + n * LHCGate; + } + + if (inf_time < time && time < sup_time) { + bx = n; + break; + } + } + return bx; +} + +std::pair RPCSynchronizer::getBX_SBX(float time) { + const float LB_clock = 25.; // 25 ns + const float LB_precise_clock = 1.5625; // 25./16. = 1.5625 ns + int BX = int(time / LB_clock); + if (time < 0) + BX--; + double dt = time - BX * LB_clock; + int SBX = int(dt / LB_precise_clock); + std::pair tdc; + tdc.first = BX; + tdc.second = SBX; + return tdc; +} + +std::tuple RPCSynchronizer::getBX_SBX_fine_time(float time) { + const float LB_clock = 25.; // 25 ns + const float LB_precise_clock = 2.5; // 2.5 ns + const float LB_fine_clock = 0.2; // 200 ps = 0.2 ns + int BX = int(time / LB_clock); + if (time < 0) + BX--; + double dt = time - BX * LB_clock; + int SBX = int(dt / LB_precise_clock); + dt = time - BX * LB_clock - SBX * LB_precise_clock; + int fine_time = int(dt / LB_fine_clock); + std::tuple tdc; + tdc = std::make_tuple(BX, SBX, fine_time); + return tdc; +} diff --git a/SimMuon/RPCDigitizer/src/RPCSynchronizer.h b/SimMuon/RPCDigitizer/src/RPCSynchronizer.h index 563b878f3cca2..1aa62fb75131b 100644 --- a/SimMuon/RPCDigitizer/src/RPCSynchronizer.h +++ b/SimMuon/RPCDigitizer/src/RPCSynchronizer.h @@ -43,6 +43,12 @@ class RPCSynchronizer { RPCSimSetUp* getRPCSimSetUp() { return theSimSetUp; } double getExactTime() const { return the_exact_time; } double getSmearedTime() const { return the_smeared_time; } + float getTiming(const PSimHit* simhit, CLHEP::HepRandomEngine* engine, float StripLength); + std::pair getDoubleTiming(const PSimHit* simhit, CLHEP::HepRandomEngine* engine, float StripLength); + int getBX(float time); + std::pair getBX_SBX(float time); + //std::pair getFineTime(const PSimHit* simhit, CLHEP::HepRandomEngine* engine,float StripLength); + std::tuple getBX_SBX_fine_time(float time); private: double resRPC; diff --git a/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.cc b/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.cc index 322434f48349f..34f53176ed261 100644 --- a/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.cc +++ b/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.cc @@ -1,7 +1,7 @@ #include "SimMuon/RPCDigitizer/src/RPCSimSetUp.h" #include "SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.h" #include "SimMuon/RPCDigitizer/src/RPCDigitizer.h" -#include "SimMuon/RPCDigitizer/src/IRPCDigitizer.h" +//#include "SimMuon/RPCDigitizer/src/IRPCDigitizer.h" #include "Geometry/Records/interface/MuonGeometryRecord.h" #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h" #include "SimDataFormats/CrossingFrame/interface/MixCollection.h" @@ -51,8 +51,9 @@ RPCandIRPCDigiProducer::RPCandIRPCDigiProducer(const edm::ParameterSet& ps) { theRPCSimSetUpRPC = new RPCSimSetUp(ps); theRPCSimSetUpIRPC = new RPCSimSetUp(ps); - theRPCDigitizer = new RPCDigitizer(ps); - theIRPCDigitizer = new IRPCDigitizer(ps); + theRPCDigitizer = new RPCDigitizer(ps, true); + //theIRPCDigitizer = new IRPCDigitizer(ps); + theIRPCDigitizer = new RPCDigitizer(ps, false); crossingFrameToken = consumes>(edm::InputTag(mix_, collection_for_XF)); geomToken = esConsumes(); noiseToken = esConsumes(); diff --git a/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.h b/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.h index 45d9a872ada96..db9e32913ab0a 100644 --- a/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.h +++ b/SimMuon/RPCDigitizer/src/RPCandIRPCDigiProducer.h @@ -10,7 +10,6 @@ #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" #include "SimMuon/RPCDigitizer/src/RPCDigitizer.h" -#include "SimMuon/RPCDigitizer/src/IRPCDigitizer.h" #include "Geometry/Records/interface/MuonGeometryRecord.h" #include "CondFormats/RPCObjects/interface/RPCStripNoises.h" #include "CondFormats/DataRecord/interface/RPCStripNoisesRcd.h" @@ -38,7 +37,7 @@ class RPCandIRPCDigiProducer : public edm::stream::EDProducer