From ca9c7d1503021242e78fa0b7680d80db78969cd0 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Thu, 25 Jul 2019 15:39:09 +0200 Subject: [PATCH 1/8] Initial commit for SimPPS/RPDigiProducer, cleaning the code, BuildFiles, cms code conventions adaptations, etc --- SimPPS/RPDigiProducer/BuildFile.xml | 2 + .../interface/RPEnergyDepositUnit.h | 30 ++ .../RPDigiProducer/interface/RPSignalPoint.h | 32 +++ SimPPS/RPDigiProducer/interface/RPSimTypes.h | 33 +++ SimPPS/RPDigiProducer/plugins/BuildFile.xml | 7 + .../plugins/DeadChannelsManager.cc | 54 ++++ .../plugins/DeadChannelsManager.h | 43 +++ .../RPDigiProducer/plugins/RPDetDigitizer.cc | 65 +++++ .../RPDigiProducer/plugins/RPDetDigitizer.h | 53 ++++ .../RPDigiProducer/plugins/RPDigiProducer.cc | 264 ++++++++++++++++++ .../plugins/RPDisplacementGenerator.cc | 97 +++++++ .../plugins/RPDisplacementGenerator.h | 54 ++++ .../plugins/RPGaussianTailNoiseAdder.cc | 46 +++ .../plugins/RPGaussianTailNoiseAdder.h | 20 ++ .../plugins/RPHitChargeConverter.cc | 24 ++ .../plugins/RPHitChargeConverter.h | 29 ++ .../RPLinearChargeCollectionDrifter.cc | 45 +++ .../plugins/RPLinearChargeCollectionDrifter.h | 24 ++ .../plugins/RPLinearChargeDivider.cc | 134 +++++++++ .../plugins/RPLinearChargeDivider.h | 42 +++ .../plugins/RPLinearInduceChargeOnStrips.cc | 51 ++++ .../plugins/RPLinearInduceChargeOnStrips.h | 23 ++ .../RPDigiProducer/plugins/RPPileUpSignals.cc | 27 ++ .../RPDigiProducer/plugins/RPPileUpSignals.h | 24 ++ .../RPDigiProducer/plugins/RPVFATSimulator.cc | 45 +++ .../RPDigiProducer/plugins/RPVFATSimulator.h | 33 +++ .../test/RPDigiProducerTest_cfg.py | 101 +++++++ 27 files changed, 1402 insertions(+) create mode 100755 SimPPS/RPDigiProducer/BuildFile.xml create mode 100644 SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h create mode 100644 SimPPS/RPDigiProducer/interface/RPSignalPoint.h create mode 100644 SimPPS/RPDigiProducer/interface/RPSimTypes.h create mode 100755 SimPPS/RPDigiProducer/plugins/BuildFile.xml create mode 100644 SimPPS/RPDigiProducer/plugins/DeadChannelsManager.cc create mode 100644 SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPDigiProducer.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPPileUpSignals.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPPileUpSignals.h create mode 100644 SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc create mode 100644 SimPPS/RPDigiProducer/plugins/RPVFATSimulator.h create mode 100644 SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py diff --git a/SimPPS/RPDigiProducer/BuildFile.xml b/SimPPS/RPDigiProducer/BuildFile.xml new file mode 100755 index 0000000000000..3b471750b7917 --- /dev/null +++ b/SimPPS/RPDigiProducer/BuildFile.xml @@ -0,0 +1,2 @@ + + diff --git a/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h b/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h new file mode 100644 index 0000000000000..c29cb5682884e --- /dev/null +++ b/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h @@ -0,0 +1,30 @@ +#ifndef SimPPS_RPDigiProducer_RP_ENERGY_DEPOSIT_UNIT_H +#define SimPPS_RPDigiProducer_RP_ENERGY_DEPOSIT_UNIT_H + +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalVector.h" + +/** + * Class which allows to "follow" an elementary charge in the silicon. + * It basically defines a quantum of energy, with a position. + */ +class RPEnergyDepositUnit { +public: + RPEnergyDepositUnit() : energy_(0), position_(0, 0, 0) {} + RPEnergyDepositUnit(double energy, double x, double y, double z) : energy_(energy), position_(x, y, z) {} + RPEnergyDepositUnit(double energy, const Local3DPoint& position) : energy_(energy), position_(position) {} + inline double X() const { return position_.x(); } + inline double Y() const { return position_.y(); } + inline double Z() const { return position_.z(); } + inline double Energy() const { return energy_; } + inline Local3DPoint Position() const { return position_; } + + inline void setEnergy(double e) { energy_ = e; } + inline void setPosition(Local3DPoint p) { position_ = p; } + +private: + double energy_; + Local3DPoint position_; +}; + +#endif //SimPPS_RPDigiProducer_RP_ENERGY_DEPOSIT_UNIT_H diff --git a/SimPPS/RPDigiProducer/interface/RPSignalPoint.h b/SimPPS/RPDigiProducer/interface/RPSignalPoint.h new file mode 100644 index 0000000000000..87a17cd408dcb --- /dev/null +++ b/SimPPS/RPDigiProducer/interface/RPSignalPoint.h @@ -0,0 +1,32 @@ +#ifndef SimPPS_RPDigiProducer_RP_SignalPoint_H +#define SimPPS_RPDigiProducer_RP_SignalPoint_H + +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalVector.h" + +/** + * An elementar charge point, with position, sigma from diffusion and tof. + */ +class RPSignalPoint { +public: + RPSignalPoint() : pos_(0, 0), sigma_(0), charge_(0) {} + + RPSignalPoint(double x, double y, double s, double charge) : pos_(x, y), sigma_(s), charge_(charge) {} + + inline LocalPoint Position() const { return pos_; } + inline double X() const { return pos_.x(); } + inline double Y() const { return pos_.y(); } + inline double Sigma() const { return sigma_; } + inline double Charge() const { return charge_; } + + inline void setSigma(double s) { sigma_ = s; } + inline void setPosition(LocalPoint p) { pos_ = p; } + inline void setCharge(double charge) { charge_ = charge; } + +private: + LocalPoint pos_; + double sigma_; + double charge_; +}; + +#endif //SimPPS_RPDigiProducer_RP_SignalPoint_H diff --git a/SimPPS/RPDigiProducer/interface/RPSimTypes.h b/SimPPS/RPDigiProducer/interface/RPSimTypes.h new file mode 100644 index 0000000000000..a765073c497d1 --- /dev/null +++ b/SimPPS/RPDigiProducer/interface/RPSimTypes.h @@ -0,0 +1,33 @@ +#ifndef SimPPS_RPDigiProducer_RPSimTypes_H +#define SimPPS_RPDigiProducer_RPSimTypes_H + +#include +#include +#include + +#include "SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h" +#include "SimPPS/RPDigiProducer/interface/RPSignalPoint.h" + +typedef uint32_t RPDetId; + +namespace simromanpot { + typedef std::map strip_charge_map; + typedef std::vector charge_induced_on_surface; + typedef std::vector energy_path_distribution; + typedef std::set > HitsContainer; + typedef std::set >::const_iterator HitsContainerIter; + typedef std::set > TriggerContainer; + typedef std::set >::const_iterator TriggerContainerIter; + typedef std::vector > > + DigiPrimaryMapType; //for each digi in the output the vector of the number of PSimHit and its weight + typedef std::vector > + SingleDigiPrimaryMapType; //for one digi in the output the vector of the number of PSimHit and its weight + typedef std::vector > > + TriggerPrimaryMapType; //for each digi in the output the vector of the number of PSimHit and its weight + typedef std::map > > + strip_charge_map_links_type; //for each strip the indeces of PSimHit and amounts of charge generated by it is given + typedef std::map > + TriggerContainerLinkMap; //for each strip the indeces of PSimHit and amounts of charge generated by it is given +} // namespace simromanpot + +#endif //SimPPS_RPDigiProducer_RPSimTypes_H diff --git a/SimPPS/RPDigiProducer/plugins/BuildFile.xml b/SimPPS/RPDigiProducer/plugins/BuildFile.xml new file mode 100755 index 0000000000000..4d308dd305519 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/BuildFile.xml @@ -0,0 +1,7 @@ + + + + + + + diff --git a/SimPPS/RPDigiProducer/plugins/DeadChannelsManager.cc b/SimPPS/RPDigiProducer/plugins/DeadChannelsManager.cc new file mode 100644 index 0000000000000..5f25150fd7779 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/DeadChannelsManager.cc @@ -0,0 +1,54 @@ +#include "SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h" +#include "CondFormats/CTPPSReadoutObjects/interface/TotemSymbId.h" +#include "CondFormats/CTPPSReadoutObjects/interface/TotemDAQMapping.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" +#include "SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h" +#include +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +DeadChannelsManager::DeadChannelsManager() { analysisMaskPresent = false; } + +DeadChannelsManager::DeadChannelsManager(edm::ESHandle _analysisMask) { + analysisMask = _analysisMask; + analysisMaskPresent = true; +} + +bool DeadChannelsManager::isChannelDead(RPDetId detectorId, unsigned short stripNumber) { + unsigned int symbolicId = RPDisplacementGenerator::rawToDecId(detectorId) * 10; //convert to symbolic ID + + unsigned int vfat = stripNumber / 128; + symbolicId += vfat; //add vfatID to symbolic ID + stripNumber = stripNumber - vfat * 128; //convert strip number to a number from range <0; 127> + TotemSymbID totemSymbolicId; + totemSymbolicId.symbolicID = symbolicId; + if (analysisMaskPresent) { + std::map::const_iterator vfatIter = + analysisMask->analysisMask.find(totemSymbolicId); + if (vfatIter != analysisMask->analysisMask.end()) { + TotemVFATAnalysisMask vfatMask = vfatIter->second; + //if channel is dead return true + if (vfatMask.fullMask || vfatMask.maskedChannels.find(stripNumber) != vfatMask.maskedChannels.end()) { + return true; + } + } + } + return false; +} + +void DeadChannelsManager::displayMap() { + if (analysisMaskPresent) { + std::map::const_iterator vfatIter; + for (vfatIter = analysisMask->analysisMask.begin(); vfatIter != analysisMask->analysisMask.end(); vfatIter++) { + LogDebug("PPSDigiProducer::DeadChannelsManager") << vfatIter->first.symbolicID << "\n"; + TotemVFATAnalysisMask am = vfatIter->second; + if (am.fullMask) { + LogDebug("PPSDigiProducer::DeadChannelsManager") << " full mask\n"; + } else { + std::set::iterator setIterator; + for (setIterator = am.maskedChannels.begin(); setIterator != am.maskedChannels.end(); setIterator++) { + LogDebug("PPSDigiProducer::DeadChannelsManager") << " " << (int)(*setIterator) << "\n"; + } + } + } + } +} diff --git a/SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h b/SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h new file mode 100644 index 0000000000000..b2d6a08502920 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h @@ -0,0 +1,43 @@ +#ifndef SimPPS_RPDigiProducer_DEAD_CHANNELS_MANAGER +#define SimPPS_RPDigiProducer_DEAD_CHANNELS_MANAGER + +#include "FWCore/Framework/interface/ESHandle.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" +#include "CondFormats/CTPPSReadoutObjects/interface/TotemAnalysisMask.h" + +/* + * This purpose of this class is to answer the question whether a channel (given by detectorId + * and stripNumber) is dead or not. This class uses analysisMask which is provided + * by DAQMappingSourceXML. + * @author Jakub Smajek + */ +class DeadChannelsManager { +private: + edm::ESHandle analysisMask; + bool analysisMaskPresent; //this variable indicates whether analysisMask is present or not + +public: + /** + * This constructor allows us to set analysisMask. The analysisMask can be read from + * EventSetup. + */ + DeadChannelsManager(edm::ESHandle analysisMask); + DeadChannelsManager(); + /** + * This function answers the question whether given channel is dead or not. + * RPDetId - detector ID given in raw form, this function has to convert raw ID to symbolic + * stripNumber - ID of the strip, it is a number from range <0; 511>, this function has to convert + * it into a vfat ID and a number from range <0; 127> + * + * It is assumed that: + * channels 0 - 127 are in vfat number 0 + * channels 128 - 255 are in vfat number 1 + * channels 256 - 383 are in vfat number 2 + * channels 384 - 511 are in vfat number 3 + */ + bool isChannelDead(RPDetId detectorId, unsigned short stripNumber); + void displayMap(); + static uint32_t rawToDecId(uint32_t raw); +}; + +#endif diff --git a/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc new file mode 100644 index 0000000000000..f0d3b3150d6a9 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc @@ -0,0 +1,65 @@ +#include +#include + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h" +#include "Geometry/VeryForwardRPTopology/interface/RPTopology.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +RPDetDigitizer::RPDetDigitizer(const edm::ParameterSet ¶ms, + CLHEP::HepRandomEngine &eng, + RPDetId det_id, + const edm::EventSetup &iSetup) + : det_id_(det_id) { + verbosity_ = params.getParameter("RPVerbosity"); + numStrips_ = RPTopology().DetStripNo(); + theNoiseInElectrons = params.getParameter("RPEquivalentNoiseCharge300um"); + theStripThresholdInE = params.getParameter("RPVFATThreshold"); + noNoise_ = params.getParameter("RPNoNoise"); + misalignment_simulation_on_ = params.getParameter("RPDisplacementOn"); + links_persistence_ = params.getParameter("RPDigiSimHitRelationsPresistence"); + + theRPGaussianTailNoiseAdder = + std::make_unique(numStrips_, theNoiseInElectrons, theStripThresholdInE, verbosity_); + theRPPileUpSignals = std::make_unique(params, det_id_); + theRPVFATSimulator = std::make_unique(params, det_id_); + theRPHitChargeConverter = std::make_unique(params, eng, det_id_); + theRPDisplacementGenerator = std::make_unique(params, det_id_, iSetup); +} + +void RPDetDigitizer::run(const std::vector &input, + const std::vector &input_links, + std::vector &output_digi, + simromanpot::DigiPrimaryMapType &output_digi_links) { + if (verbosity_) + LogDebug("RPDetDigitizer ") << det_id_ << " received input.size()=" << input.size() << "\n"; + theRPPileUpSignals->reset(); + + bool links_persistence_checked = links_persistence_ && input_links.size() == input.size(); + + int input_size = input.size(); + for (int i = 0; i < input_size; ++i) { + simromanpot::strip_charge_map the_strip_charge_map; + if (misalignment_simulation_on_) + the_strip_charge_map = theRPHitChargeConverter->processHit(theRPDisplacementGenerator->displace(input[i])); + else + the_strip_charge_map = theRPHitChargeConverter->processHit(input[i]); + + if (verbosity_) + LogDebug("RPHitChargeConverter ") << det_id_ << " returned hits=" << the_strip_charge_map.size() << "\n"; + if (links_persistence_checked) + theRPPileUpSignals->add(the_strip_charge_map, input_links[i]); + else + theRPPileUpSignals->add(the_strip_charge_map, 0); + } + + const simromanpot::strip_charge_map &theSignal = theRPPileUpSignals->dumpSignal(); + simromanpot::strip_charge_map_links_type &theSignalProvenance = theRPPileUpSignals->dumpLinks(); + simromanpot::strip_charge_map afterNoise; + if (noNoise_) + afterNoise = theSignal; + else + afterNoise = theRPGaussianTailNoiseAdder->addNoise(theSignal); + + theRPVFATSimulator->ConvertChargeToHits(afterNoise, theSignalProvenance, output_digi, output_digi_links); +} diff --git a/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h new file mode 100644 index 0000000000000..7073b60208282 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h @@ -0,0 +1,53 @@ +#ifndef SimPPS_RPDigiProducer_RP_DET_DIGITIZER_H +#define SimPPS_RPDigiProducer_RP_DET_DIGITIZER_H + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/CTPPSDigi/interface/TotemRPDigi.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" + +#include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h" + +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" +#include "SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h" +#include "SimPPS/RPDigiProducer/plugins/RPVFATSimulator.h" +#include "SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h" +#include "SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h" +#include "SimPPS/RPDigiProducer/plugins/RPPileUpSignals.h" + +#include +#include + +namespace CLHEP { + class HepRandomEngine; +} + +class RPDetDigitizer { +public: + RPDetDigitizer(const edm::ParameterSet ¶ms, + CLHEP::HepRandomEngine &eng, + RPDetId det_id, + const edm::EventSetup &iSetup); + void run(const std::vector &input, + const std::vector &input_links, + std::vector &output_digi, + simromanpot::DigiPrimaryMapType &output_digi_links); + +private: + std::unique_ptr theRPGaussianTailNoiseAdder; + std::unique_ptr theRPPileUpSignals; + std::unique_ptr theRPHitChargeConverter; + std::unique_ptr theRPVFATSimulator; + std::unique_ptr theRPDisplacementGenerator; + +private: + int numStrips_; + double theNoiseInElectrons; // Noise (RMS) in units of electrons. + double theStripThresholdInE; // Strip noise treshold in electorns. + bool noNoise_; //if the nos is included + RPDetId det_id_; + bool misalignment_simulation_on_; + int verbosity_; + bool links_persistence_; +}; + +#endif //SimCTPPS_RPDigiProducer_RP_DET_DIGITIZER_H diff --git a/SimPPS/RPDigiProducer/plugins/RPDigiProducer.cc b/SimPPS/RPDigiProducer/plugins/RPDigiProducer.cc new file mode 100644 index 0000000000000..2c5bbd8b275ab --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPDigiProducer.cc @@ -0,0 +1,264 @@ +// user include files +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" + +#include "DataFormats/CTPPSDigi/interface/TotemRPDigi.h" +#include "DataFormats/CTPPSDigi/interface/TotemRPDigi.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" +#include "DataFormats/Common/interface/DetSet.h" + +//Random Number +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/RandomNumberGenerator.h" + +#include "CondFormats/CTPPSReadoutObjects/interface/TotemAnalysisMask.h" +#include "CondFormats/DataRecord/interface/TotemReadoutRcd.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" +#include "SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h" +#include "SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h" + +// system include files +#include +#include +#include +#include +#include +#include // I need it for random numbers + +// user include files + +// +// class decleration +// + +namespace CLHEP { + class HepRandomEngine; +} + +class RPDigiProducer : public edm::EDProducer { +public: + explicit RPDigiProducer(const edm::ParameterSet&); + ~RPDigiProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void beginRun(const edm::Run&, const edm::EventSetup&) override; + void produce(edm::Event&, const edm::EventSetup&) override; + + edm::DetSet convertRPStripDetSet(const edm::DetSet&); + + // ----------member data --------------------------- + std::vector RP_hit_containers_; + typedef std::map> simhit_map; + typedef simhit_map::iterator simhit_map_iterator; + + edm::ParameterSet conf_; + std::map> theAlgoMap; + + CLHEP::HepRandomEngine* rndEngine_ = nullptr; + int verbosity_; + + /** + * this variable answers the question whether given channel is dead or not + */ + DeadChannelsManager deadChannelsManager; + /** + * this variable indicates whether we take into account dead channels or simulate as if all + * channels work ok (by default we do not simulate dead channels) + */ + bool simulateDeadChannels; + + edm::EDGetTokenT> tokenCrossingFrameTotemRP; +}; + +RPDigiProducer::RPDigiProducer(const edm::ParameterSet& conf) : conf_(conf) { + //now do what ever other initialization is needed + produces>(); + + // register data to consume + tokenCrossingFrameTotemRP = consumes>(edm::InputTag("mix", "g4SimHitsTotemHitsRP", "")); + + RP_hit_containers_ = conf.getParameter>("ROUList"); + verbosity_ = conf.getParameter("RPVerbosity"); + + simulateDeadChannels = false; + if (conf.exists( + "simulateDeadChannels")) { //check if "simulateDeadChannels" variable is defined in configuration file + simulateDeadChannels = conf.getParameter("simulateDeadChannels"); + } +} + +RPDigiProducer::~RPDigiProducer() { + // do anything here that needs to be done at desctruction time + // (e.g. close files, deallocate resources etc.) +} +// +// member functions +// + +// ------------ method called to produce the data ------------ +void RPDigiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + using namespace edm; + + // initialize random engine + if (!rndEngine_) { + Service rng; + if (!rng.isAvailable()) { + throw cms::Exception("Configuration") + << "This class 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."; + } + rndEngine_ = &(rng->getEngine(iEvent.streamID())); + } + + // Step A: Get Inputs + edm::Handle> cf; + iEvent.getByLabel("mix", "g4SimHitsTotemHitsRP", cf); + + if (verbosity_) { + edm::LogInfo("RPDigiProducer") << "\n\n=================== Starting SimHit access" + << " ===================" + << "\n"; + + MixCollection col{cf.product(), std::pair(-0, 0)}; + MixCollection::iterator cfi; + int count = 0; + for (cfi = col.begin(); cfi != col.end(); cfi++) { + edm::LogInfo("RPDigiProducer") << " Hit " << count << " has tof " << cfi->timeOfFlight() << " trackid " + << cfi->trackId() << " bunchcr " << cfi.bunch() << " trigger " << cfi.getTrigger() + << ", from EncodedEventId: " << cfi->eventId().bunchCrossing() << " " + << cfi->eventId().event() << " bcr from MixCol " << cfi.bunch() << "\n"; + edm::LogInfo("RPDigiProducer") << " Hit: " << (*cfi) << "\n"; + count++; + } + } + + MixCollection allRPHits{cf.product(), std::pair(0, 0)}; + + if (verbosity_) + edm::LogInfo("RPDigiProducer") << "Input MixCollection size = " << allRPHits.size() << "\n"; + + //Loop on PSimHit + simhit_map simHitMap_; + simHitMap_.clear(); + + MixCollection::iterator isim; + for (isim = allRPHits.begin(); isim != allRPHits.end(); ++isim) { + simHitMap_[(*isim).detUnitId()].push_back((*isim)); + } + + // Step B: LOOP on hits in event + std::vector> DigiVector; + DigiVector.reserve(400); + DigiVector.clear(); + + for (simhit_map_iterator it = simHitMap_.begin(); it != simHitMap_.end(); ++it) { + edm::DetSet digi_collector(it->first); + + if (theAlgoMap.find(it->first) == theAlgoMap.end()) { + theAlgoMap[it->first] = + std::unique_ptr(new RPDetDigitizer(conf_, *rndEngine_, it->first, iSetup)); + } + + std::vector input_links; + simromanpot::DigiPrimaryMapType output_digi_links; + + (theAlgoMap.find(it->first)->second) + ->run(simHitMap_[it->first], input_links, digi_collector.data, output_digi_links); + + if (!digi_collector.data.empty()) { + DigiVector.push_back(convertRPStripDetSet(digi_collector)); + } + } + + // Step C: create empty output collection + std::unique_ptr> digi_output(new edm::DetSetVector(DigiVector)); + + if (verbosity_) { + edm::LogInfo("RPDigiProducer") << "digi_output->size()=" << digi_output->size() << "\n"; + } + // Step D: write output to file + iEvent.put(std::move(digi_output)); +} + +// ------------ method called once each job just before starting event loop ------------ +void RPDigiProducer::beginRun(const edm::Run& beginrun, const edm::EventSetup& es) { + // get analysis mask to mask channels + if (simulateDeadChannels) { + edm::ESHandle analysisMask; + es.get().get(analysisMask); + deadChannelsManager = DeadChannelsManager(analysisMask); //set analysisMask in deadChannelsManager + } +} + +edm::DetSet RPDigiProducer::convertRPStripDetSet(const edm::DetSet& rpstrip_detset) { + edm::DetSet rpdigi_detset(rpstrip_detset.detId()); + rpdigi_detset.reserve(rpstrip_detset.size()); + + for (std::vector::const_iterator stripIterator = rpstrip_detset.data.begin(); + stripIterator < rpstrip_detset.data.end(); + ++stripIterator) { + rpdigi_detset.push_back(TotemRPDigi(stripIterator->getStripNumber())); + } + + return rpdigi_detset; +} + +void RPDigiProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //RPSiDetDigitizer + //all distances in [mm] + edm::ParameterSetDescription desc; + desc.add("RPLandauFluctuations", true); + desc.add("RPDisplacementOn", false); + desc.add("RPVerbosity", 0); + desc.add("RPVFATThreshold", 9000.0); + desc.add("RPTopEdgePosition", 1.5); + desc.add("RPActiveEdgeSmearing", 0.013); + desc.add("RPEquivalentNoiseCharge300um", 1000.0); + desc.add("RPVFATTriggerMode", 2); + desc.add>("RPInterStripSmearing", + { + 0.011, + }); + desc.add("RPSharingSigmas", 5.0); //how many sigmas taken into account for the edges and inter strips + desc.add("RPGeVPerElectron", 3.61e-09); + desc.add("RPActiveEdgePosition", 0.034); //from the physical edge + desc.add("RPDeadStripSimulationOn", false); + desc.add>("ROUList", + { + "TotemHitsRP", + }); + desc.add("RPNoNoise", false); + desc.add("RPDigiSimHitRelationsPresistence", false); //save links betweend digi, clusters and OSCAR/Geant4 hits + desc.add("mixLabel", "mix"); + desc.add("RPChargeDivisionsPerThickness", 5); + desc.add("RPDeltaProductionCut", 0.120425); //[MeV] + desc.add("RPBottomEdgePosition", 1.5); + desc.add("RPBottomEdgeSmearing", 0.011); + desc.add("RPTopEdgeSmearing", 0.011); + desc.add("InputCollection", "g4SimHitsTotemHitsRP"); + desc.add("RPInterStripCoupling", + 1.0); //fraction of charge going to the strip, the missing part is taken by its neighbours + desc.add("RPDeadStripProbability", 0.001); + desc.add("RPChargeDivisionsPerStrip", 15); + descriptions.add("RPSiDetDigitizer", desc); +} + +DEFINE_FWK_MODULE(RPDigiProducer); diff --git a/SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.cc b/SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.cc new file mode 100644 index 0000000000000..31414d36049b4 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.cc @@ -0,0 +1,97 @@ +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h" +#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" +#include "CondFormats/CTPPSReadoutObjects/interface/CTPPSRPAlignmentCorrectionsData.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h" + +#include +#include + +using namespace std; +using namespace edm; + +RPDisplacementGenerator::RPDisplacementGenerator(const edm::ParameterSet &ps, + RPDetId _detId, + const edm::EventSetup &iSetup) + : detId_(_detId) { + isOn_ = ps.getParameter("RPDisplacementOn"); + + // read the alignment correction + ESHandle alignments; + if (auto rec = iSetup.tryToGet()) { + iSetup.get().get(alignments); + } + + unsigned int decId = rawToDecId(detId_); + + math::XYZVectorD S_m; + RotationMatrix R_m; + + if (alignments.isValid()) { + const CTPPSRPAlignmentCorrectionData &ac = alignments->getFullSensorCorrection(decId); + S_m = ac.getTranslation(); + R_m = ac.getRotationMatrix(); + } else + isOn_ = false; + + // transform shift and rotation to the local coordinate frame + ESHandle geom; + iSetup.get().get(geom); + const DetGeomDesc *g = geom->getSensor(detId_); + const DDRotationMatrix &R_l = g->rotation(); + rotation_ = R_l.Inverse() * R_m.Inverse() * R_l; + shift_ = R_l.Inverse() * R_m.Inverse() * S_m; + + LogDebug("RPDisplacementGenerator").log([&](auto &log) { + log << " det id = " << decId << ", isOn = " << isOn_ << "\n"; + if (isOn_) { + log << " shift = " << shift_ << "\n"; + log << " rotation = " << rotation_ << "\n"; + } + }); +} + +Local3DPoint RPDisplacementGenerator::displacePoint(const Local3DPoint &p) { + /// input is in mm, shifts are in mm too + + DDTranslation v(p.x(), p.y(), p.z()); + v = rotation_ * v - shift_; + + return Local3DPoint(v.x(), v.y(), v.z()); +} + +PSimHit RPDisplacementGenerator::displace(const PSimHit &input) { + if (!isOn_) + return input; + + const Local3DPoint &ep = input.entryPoint(), &xp = input.exitPoint(); + const Local3DPoint &dep = displacePoint(ep), &dxp = displacePoint(xp); + + LogDebug("RPDisplacementGenerator::displace\n") << " entry point: " << ep << " -> " << dep << "\n" + << " exit point : " << xp << " -> " << dxp << "\n"; + + return PSimHit(dep, + dxp, + input.pabs(), + input.tof(), + input.energyLoss(), + input.particleType(), + input.detUnitId(), + input.trackId(), + input.thetaAtEntry(), + input.phiAtEntry(), + input.processType()); +} + +uint32_t RPDisplacementGenerator::rawToDecId(uint32_t raw) { + return ((raw >> CTPPSDetId::startArmBit) & CTPPSDetId::maskArm) * 1000 + + ((raw >> CTPPSDetId::startStationBit) & CTPPSDetId::maskStation) * 100 + + ((raw >> CTPPSDetId::startRPBit) & CTPPSDetId::maskRP) * 10 + + ((raw >> TotemRPDetId::startPlaneBit) & TotemRPDetId::maskPlane); +} diff --git a/SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h b/SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h new file mode 100644 index 0000000000000..9911935dc0d2e --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPDisplacementGenerator.h @@ -0,0 +1,54 @@ +#ifndef SimPPS_RPDigiProducer_RP_DISPLACEMENT_GENERATOR_H +#define SimPPS_RPDigiProducer_RP_DISPLACEMENT_GENERATOR_H + +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" +#include "DetectorDescription/Core/interface/DDRotationMatrix.h" +#include "DetectorDescription/Core/interface/DDTranslation.h" +#include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h" +#include +#include + +namespace edm { + class ParameterSet; + class EventSetup; +} // namespace edm + +class PSimHit; + +/** + * \ingroup TotemDigiProduction + * \brief This class introduces displacements of RP. + * It actually shifts and rotates PSimHit positions. It doesn't test whether the displaced + * hit is still on the detector's surface. This check takes place later in the process. It + * is done via edge effectivity. + * + * PSimHit points are given in the "local Det frame" (PSimHit.h) + */ + +class RPDisplacementGenerator { +public: + typedef ROOT::Math::Rotation3D RotationMatrix; + + RPDisplacementGenerator(const edm::ParameterSet &, RPDetId, const edm::EventSetup &); + + /// returns displaced PSimHit + PSimHit displace(const PSimHit &); + + static uint32_t rawToDecId(uint32_t raw); + +private: + /// ID of the detector + RPDetId detId_; + + /// displacement + DDTranslation shift_; + DDRotationMatrix rotation_; + + /// set to false to bypass displacements + bool isOn_; + + /// displaces a point + Local3DPoint displacePoint(const Local3DPoint &); +}; + +#endif diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc new file mode 100644 index 0000000000000..1d74a4ba16972 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc @@ -0,0 +1,46 @@ +#include "SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include +#include "TMath.h" +#include "TRandom.h" +#include "CLHEP/Random/RandGauss.h" + +using namespace std; + +RPGaussianTailNoiseAdder::RPGaussianTailNoiseAdder(int numStrips, + double theNoiseInElectrons, + double theStripThresholdInE, + int verbosity) + : numStrips_(numStrips), theNoiseInElectrons(theNoiseInElectrons), theStripThresholdInE(theStripThresholdInE) { + verbosity_ = verbosity; + strips_above_threshold_prob_ = TMath::Erfc(theStripThresholdInE / sqrt(2.0) / theNoiseInElectrons) / 2; +} + +simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanpot::strip_charge_map &theSignal) { + simromanpot::strip_charge_map the_strip_charge_map; + + // noise on strips with signal: + for (simromanpot::strip_charge_map::const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) { + double noise = CLHEP::RandGauss::shoot(0.0, theNoiseInElectrons); + the_strip_charge_map[i->first] = i->second + noise; + if (verbosity_) + edm::LogInfo("RPDigiProducer") << "noise added to signal strips: " << noise << "\n"; + } + + // noise on the other strips + int strips_no_above_threshold = gRandom->Binomial(numStrips_, strips_above_threshold_prob_); + + for (int j = 0; j < strips_no_above_threshold; j++) { + int strip = gRandom->Integer(numStrips_); + if (the_strip_charge_map[strip] == 0) { + the_strip_charge_map[strip] = 2 * theStripThresholdInE; + //only binary decision later, no need to simulate the noise precisely, + //enough to know that it exceeds the threshold + if (verbosity_) + edm::LogInfo("RPDigiProducer") << "nonsignal strips noise :" << strip << " " << the_strip_charge_map[strip] + << "\n"; + } + } + + return the_strip_charge_map; +} diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h new file mode 100644 index 0000000000000..31a5a7c385ab4 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h @@ -0,0 +1,20 @@ +#ifndef SimPPS_RPDigiProducer_RP_GAUSSIAN_TAIL_NOISE_ADDER_H +#define SimPPS_RPDigiProducer_RP_GAUSSIAN_TAIL_NOISE_ADDER_H + +#include "SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" + +class RPGaussianTailNoiseAdder { +public: + RPGaussianTailNoiseAdder(int numStrips, double theNoiseInElectrons, double theStripThresholdInE, int verbosity); + simromanpot::strip_charge_map addNoise(const simromanpot::strip_charge_map &theSignal); + +private: + int numStrips_; + double theNoiseInElectrons; + double theStripThresholdInE; + double strips_above_threshold_prob_; + int verbosity_; +}; + +#endif //SimPPS_RPDigiProducer_RP_GAUSSIAN_TAIL_NOISE_ADDER_H diff --git a/SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.cc b/SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.cc new file mode 100644 index 0000000000000..f74f8d80c43b3 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.cc @@ -0,0 +1,24 @@ +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h" +#include "SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h" +#include "SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h" +#include "SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h" + +RPHitChargeConverter::RPHitChargeConverter(const edm::ParameterSet ¶ms, CLHEP::HepRandomEngine &eng, RPDetId det_id) + : det_id_(det_id) { + verbosity_ = params.getParameter("RPVerbosity"); + theRPChargeDivider = std::make_unique(params, eng, det_id); + theRPChargeCollectionDrifter = std::make_unique(params, det_id); + theRPInduceChargeOnStrips = std::make_unique(params, det_id); +} + +RPHitChargeConverter::~RPHitChargeConverter() {} + +simromanpot::strip_charge_map RPHitChargeConverter::processHit(const PSimHit &hit) { + simromanpot::energy_path_distribution ions_along_path = theRPChargeDivider->divide(hit); + if (verbosity_) + edm::LogInfo("HitChargeConverter") << det_id_ << " clouds no generated on the path=" << ions_along_path.size() + << "\n"; + return theRPInduceChargeOnStrips->Induce(theRPChargeCollectionDrifter->Drift(ions_along_path)); +} diff --git a/SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h b/SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h new file mode 100644 index 0000000000000..183d5fd339c5f --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPHitChargeConverter.h @@ -0,0 +1,29 @@ +#ifndef SimPPS_RPDigiProducer_RP_HIT_CHARGE_CONVERTER_H +#define SimPPS_RPDigiProducer_RP_HIT_CHARGE_CONVERTER_H + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h" +#include "SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h" +#include "SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" + +#include + +class RPHitChargeConverter { +public: + RPHitChargeConverter(const edm::ParameterSet ¶ms_, CLHEP::HepRandomEngine &eng, RPDetId det_id); + ~RPHitChargeConverter(); + + simromanpot::strip_charge_map processHit(const PSimHit &hit); + +private: + const RPDetId det_id_; + + std::unique_ptr theRPChargeDivider; + std::unique_ptr theRPChargeCollectionDrifter; + std::unique_ptr theRPInduceChargeOnStrips; + int verbosity_; +}; + +#endif //SimPPS_RPDigiProducer_RP_HIT_CHARGE_CONVERTER_H diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc new file mode 100644 index 0000000000000..ad0f8b70809a6 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc @@ -0,0 +1,45 @@ +#include "SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h" +#include "Geometry/VeryForwardRPTopology/interface/RPTopology.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include +#include + +RPLinearChargeCollectionDrifter::RPLinearChargeCollectionDrifter(const edm::ParameterSet ¶ms, RPDetId det_id) { + verbosity_ = params.getParameter("RPVerbosity"); + GeV_per_electron_ = params.getParameter("RPGeVPerElectron"); + charge_cloud_sigmas_vect_ = params.getParameter >("RPInterStripSmearing"); + det_thickness_ = RPTopology().DetThickness(); + det_id_ = det_id; +} + +simromanpot::charge_induced_on_surface RPLinearChargeCollectionDrifter::Drift( + const simromanpot::energy_path_distribution &energy_deposition) { + simromanpot::charge_induced_on_surface temp_; + temp_.resize(energy_deposition.size()); + for (unsigned int i = 0; i < energy_deposition.size(); i++) { + temp_[i].setPosition(LocalPoint(energy_deposition[i].X(), energy_deposition[i].Y())); + temp_[i].setSigma(getSigma(energy_deposition[i].Z())); //befor charge_cloud_sigma_ used, now a vector of sigmas; + temp_[i].setCharge(energy_deposition[i].Energy() / GeV_per_electron_); + if (verbosity_) { + edm::LogInfo("RPLinearChargeCollectionDrifter") + << det_id_ << " :" << temp_[i].Position() << " " << temp_[i].Sigma() << " " << temp_[i].Charge() << "\n"; + } + } + return temp_; +} +double RPLinearChargeCollectionDrifter::getSigma(double z) { + if (charge_cloud_sigmas_vect_.size() == 1) + return charge_cloud_sigmas_vect_[0]; + + double factor = (z / det_thickness_) * (charge_cloud_sigmas_vect_.size() - 1); + double lo_i = floor(factor); + double hi_i = ceil(factor); + if (lo_i == hi_i) { + return charge_cloud_sigmas_vect_[(int)factor]; + } else { + double lo_weight = hi_i - factor; + double hi_weight = factor - lo_i; + + return charge_cloud_sigmas_vect_[(int)lo_i] * lo_weight + charge_cloud_sigmas_vect_[(int)hi_i] * hi_weight; + } +} diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h new file mode 100644 index 0000000000000..a452091e303d9 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.h @@ -0,0 +1,24 @@ +#ifndef SimPPS_RPDigiProducer_RP_LINEAR_CHARGE_COLLECTION_DRIFTER_H +#define SimPPS_RPDigiProducer_RP_LINEAR_CHARGE_COLLECTION_DRIFTER_H + +#include +#include +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" + +class RPLinearChargeCollectionDrifter { +public: + RPLinearChargeCollectionDrifter(const edm::ParameterSet ¶ms, RPDetId det_id); + simromanpot::charge_induced_on_surface Drift(const simromanpot::energy_path_distribution &energy_deposition); + +private: + std::vector charge_cloud_sigmas_vect_; + double GeV_per_electron_; + int verbosity_; + double det_thickness_; + RPDetId det_id_; + + double getSigma(double z); //z - z position +}; + +#endif //SimPPS_RPDigiProducer_RP_LINEAR_CHARGE_COLLECTION_DRIFTER_H diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc new file mode 100644 index 0000000000000..d1d4db43fe3f6 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc @@ -0,0 +1,134 @@ +#include "SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalVector.h" +#include "Geometry/VeryForwardRPTopology/interface/RPTopology.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +RPLinearChargeDivider::RPLinearChargeDivider(const edm::ParameterSet& params, + CLHEP::HepRandomEngine& eng, + RPDetId det_id) + : params_(params), rndEngine_(eng), det_id_(det_id) { + verbosity_ = params.getParameter("RPVerbosity"); + + fluctuate_ = std::make_unique(); + + // To Run APV in peak instead of deconvolution mode, which degrades the time resolution. + //use: SimpleConfigurable SiLinearChargeDivider::peakMode(false,"SiStripDigitizer:APVpeakmode"); + + // To Enable interstrip Landau fluctuations within a cluster. + //use: SimpleConfigurable SiLinearChargeDivider::fluctuateCharge(true,"SiStripDigitizer:LandauFluctuations"); + fluctuateCharge_ = params.getParameter("RPLandauFluctuations"); + + // Number of segments per strip into which charge is divided during + // simulation. If large, precision of simulation improves. + //to do so: SimpleConfigurable SiLinearChargeDivider::chargeDivisionsPerStrip(10,"SiStripDigitizer:chargeDivisionsPerStrip"); + chargedivisionsPerStrip_ = params.getParameter("RPChargeDivisionsPerStrip"); + chargedivisionsPerThickness_ = params.getParameter("RPChargeDivisionsPerThickness"); + + // delta cutoff in MeV, has to be same as in OSCAR (0.120425 MeV corresponding // to 100um range for electrons) + // SimpleConfigurable SiLinearChargeDivider::deltaCut(0.120425, + deltaCut_ = params.getParameter("RPDeltaProductionCut"); + + RPTopology rp_det_topol; + pitch_ = rp_det_topol.DetPitch(); + thickness_ = rp_det_topol.DetThickness(); +} + +RPLinearChargeDivider::~RPLinearChargeDivider() {} + +simromanpot::energy_path_distribution RPLinearChargeDivider::divide(const PSimHit& hit) { + LocalVector direction = hit.exitPoint() - hit.entryPoint(); + if (direction.z() > 10 || direction.x() > 200 || direction.y() > 200) { + the_energy_path_distribution_.clear(); + return the_energy_path_distribution_; + } + + int NumberOfSegmentation_y = (int)(1 + chargedivisionsPerStrip_ * fabs(direction.y()) / pitch_); + int NumberOfSegmentation_z = (int)(1 + chargedivisionsPerThickness_ * fabs(direction.z()) / thickness_); + int NumberOfSegmentation = std::max(NumberOfSegmentation_y, NumberOfSegmentation_z); + + double eLoss = hit.energyLoss(); // Eloss in GeV + + the_energy_path_distribution_.resize(NumberOfSegmentation); + + if (fluctuateCharge_) { + int pid = hit.particleType(); + double momentum = hit.pabs(); + double length = direction.mag(); // Track length in Silicon + FluctuateEloss(pid, momentum, eLoss, length, NumberOfSegmentation, the_energy_path_distribution_); + for (int i = 0; i < NumberOfSegmentation; i++) { + the_energy_path_distribution_[i].setPosition(hit.entryPoint() + + double((i + 0.5) / NumberOfSegmentation) * direction); + } + } else { + for (int i = 0; i < NumberOfSegmentation; i++) { + the_energy_path_distribution_[i].setPosition(hit.entryPoint() + + double((i + 0.5) / NumberOfSegmentation) * direction); + the_energy_path_distribution_[i].setEnergy(eLoss / (double)NumberOfSegmentation); + } + } + + if (verbosity_) { + edm::LogInfo("RPLinearChargeDivider") << det_id_ << " charge along the track:\n"; + double sum = 0; + for (unsigned int i = 0; i < the_energy_path_distribution_.size(); i++) { + edm::LogInfo("RPLinearChargeDivider") + << the_energy_path_distribution_[i].X() << " " << the_energy_path_distribution_[i].Y() << " " + << the_energy_path_distribution_[i].Z() << " " << the_energy_path_distribution_[i].Energy() << "\n"; + sum += the_energy_path_distribution_[i].Energy(); + } + edm::LogInfo("RPLinearChargeDivider") << "energy dep. sum=" << sum << "\n"; + } + + return the_energy_path_distribution_; +} + +void RPLinearChargeDivider::FluctuateEloss(int pid, + double particleMomentum, + double eloss, + double length, + int NumberOfSegs, + simromanpot::energy_path_distribution& elossVector) { + double particleMass = 139.6; // Mass in MeV, Assume pion + pid = std::abs(pid); + if (pid != 211) { // Mass in MeV + if (pid == 11) + particleMass = 0.511; + else if (pid == 13) + particleMass = 105.7; + else if (pid == 321) + particleMass = 493.7; + else if (pid == 2212) + particleMass = 938.3; + } + + double segmentLength = length / NumberOfSegs; + + // Generate charge fluctuations. + double de = 0.; + double sum = 0.; + double segmentEloss = (eloss * 1000) / NumberOfSegs; //eloss in MeV + for (int i = 0; i < NumberOfSegs; i++) { + // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV, + // track segment length in mm, segment eloss in MeV + // Returns fluctuated eloss in MeV + double deltaCutoff = deltaCut_; + de = fluctuate_->SampleFluctuations( + particleMomentum * 1000, particleMass, deltaCutoff, segmentLength, segmentEloss, &(rndEngine_)) / + 1000; //convert to GeV + elossVector[i].setEnergy(de); + sum += de; + } + + if (sum > 0.) { // If fluctuations give eloss>0. + // Rescale to the same total eloss + double ratio = eloss / sum; + for (int ii = 0; ii < NumberOfSegs; ii++) + elossVector[ii].setEnergy(ratio * elossVector[ii].Energy()); + } else { // If fluctuations gives 0 eloss + double averageEloss = eloss / NumberOfSegs; + for (int ii = 0; ii < NumberOfSegs; ii++) + elossVector[ii].setEnergy(averageEloss); + } + return; +} diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h new file mode 100644 index 0000000000000..998470a86760d --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.h @@ -0,0 +1,42 @@ +#ifndef SimPPS_RPDigiProducer_RP_LINEAR_CHARGE_DIVIDER_H +#define SimPPS_RPDigiProducer_RP_LINEAR_CHARGE_DIVIDER_H + +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimTracker/Common/interface/SiG4UniversalFluctuation.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +namespace CLHEP { + class HepRandomEngine; +} + +class RPLinearChargeDivider { +public: + RPLinearChargeDivider(const edm::ParameterSet& params, CLHEP::HepRandomEngine& eng, RPDetId det_id); + ~RPLinearChargeDivider(); + simromanpot::energy_path_distribution divide(const PSimHit& hit); + +private: + const edm::ParameterSet& params_; + CLHEP::HepRandomEngine& rndEngine_; + RPDetId det_id_; + + bool fluctuateCharge_; + int chargedivisionsPerStrip_; + int chargedivisionsPerThickness_; + double deltaCut_; + double pitch_; + double thickness_; + simromanpot::energy_path_distribution the_energy_path_distribution_; + std::unique_ptr fluctuate_; + int verbosity_; + + void FluctuateEloss(int pid, + double particleMomentum, + double eloss, + double length, + int NumberOfSegs, + simromanpot::energy_path_distribution& elossVector); +}; + +#endif //SimPPS_RPDigiProducer_RP_LINEAR_CHARGE_DIVIDER_H diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc b/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc new file mode 100644 index 0000000000000..6a49a44f12162 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc @@ -0,0 +1,51 @@ +#include "SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include + +RPLinearInduceChargeOnStrips::RPLinearInduceChargeOnStrips(const edm::ParameterSet ¶ms, RPDetId det_id) + : det_id_(det_id), theRPDetTopology(params) { + verbosity_ = params.getParameter("RPVerbosity"); + signalCoupling_.clear(); + double coupling_constant_ = params.getParameter("RPInterStripCoupling"); + signalCoupling_.push_back(coupling_constant_); + signalCoupling_.push_back((1.0 - coupling_constant_) / 2.); + + no_of_strips_ = theRPDetTopology.DetStripNo(); +} + +simromanpot::strip_charge_map RPLinearInduceChargeOnStrips::Induce( + const simromanpot::charge_induced_on_surface &charge_map) { + theStripChargeMap.clear(); + const double sqrt_2 = sqrt(2.0); + if (verbosity_) + edm::LogInfo("RPLinearInduceChargeOnStrips ") << det_id_ << " : Clouds to be induced:" << charge_map.size() << "\n"; + for (simromanpot::charge_induced_on_surface::const_iterator i = charge_map.begin(); i != charge_map.end(); ++i) { + double hit_pos; + std::vector relevant_strips = + theRPDetTopology.GetStripsInvolved((*i).X(), (*i).Y(), (*i).Sigma(), hit_pos); + if (verbosity_) { + edm::LogInfo("RPLinearInduceChargeOnStrips ") + << det_id_ << " : relevant_strips" << relevant_strips.size() << "\n"; + } + for (std::vector::const_iterator j = relevant_strips.begin(); j != relevant_strips.end(); ++j) { + double strip_begin = (*j).LowerBoarder(); + double strip_end = (*j).HigherBoarder(); + double effic = (*j).EffFactor(); + double sigma = (*i).Sigma(); + unsigned short str_no = (*j).StripNo(); + + double charge_on_strip = (TMath::Erfc((strip_begin - hit_pos) / sqrt_2 / sigma) / 2.0 - + TMath::Erfc((strip_end - hit_pos) / sqrt_2 / sigma) / 2.0) * + (*i).Charge() * effic; + if (verbosity_) + edm::LogInfo("RPLinearInduceChargeOnStrips") << "Efficiency " << det_id_ << " :" << effic << "\n"; + + for (int k = -signalCoupling_.size() + 1; k < (int)signalCoupling_.size(); ++k) { + if ((str_no + k) >= 0 && (str_no + k) < no_of_strips_) + theStripChargeMap[str_no + k] += charge_on_strip * signalCoupling_[abs(k)]; + } + } + } + + return theStripChargeMap; +} diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h b/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h new file mode 100644 index 0000000000000..fd7a1e2639f21 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.h @@ -0,0 +1,23 @@ +#ifndef SimPPS_RPDigiProducer_RP_LINEAR_INDUCE_CHARGE_ON_STRIPS_H +#define SimPPS_RPDigiProducer_RP_LINEAR_INDUCE_CHARGE_ON_STRIPS_H + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" +#include "Geometry/VeryForwardRPTopology/interface/RPSimTopology.h" + +class RPLinearInduceChargeOnStrips { +public: + RPLinearInduceChargeOnStrips(const edm::ParameterSet ¶ms, RPDetId det_id); + simromanpot::strip_charge_map Induce(const simromanpot::charge_induced_on_surface &charge_map); + +private: + RPDetId det_id_; + std::vector signalCoupling_; + simromanpot::strip_charge_map theStripChargeMap; + RPSimTopology theRPDetTopology; + int no_of_strips_; + int verbosity_; +}; + +#endif diff --git a/SimPPS/RPDigiProducer/plugins/RPPileUpSignals.cc b/SimPPS/RPDigiProducer/plugins/RPPileUpSignals.cc new file mode 100644 index 0000000000000..f36d2a53af5e7 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPPileUpSignals.cc @@ -0,0 +1,27 @@ + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "SimPPS/RPDigiProducer/plugins/RPPileUpSignals.h" +#include + +RPPileUpSignals::RPPileUpSignals(const edm::ParameterSet ¶ms, RPDetId det_id) : det_id_(det_id) { + links_persistence_ = params.getParameter("RPDigiSimHitRelationsPresistence"); + verbosity_ = params.getParameter("RPVerbosity"); +} + +void RPPileUpSignals::reset() { + the_strip_charge_piled_up_map_.clear(); + the_strip_charge_piled_up_map_links_.clear(); +} + +void RPPileUpSignals::add(const simromanpot::strip_charge_map &charge_induced, int PSimHitIndex) { + for (simromanpot::strip_charge_map::const_iterator i = charge_induced.begin(); i != charge_induced.end(); ++i) { + the_strip_charge_piled_up_map_[i->first] += i->second; + if (links_persistence_ && i->second > 0) { + the_strip_charge_piled_up_map_links_[i->first].push_back(std::pair(PSimHitIndex, i->second)); + if (verbosity_) { + edm::LogInfo("RPPileUpSignals") << "Det id=" << det_id_ << " strip=" << i->first << " charge=" << i->second + << "\n"; + } + } + } +} diff --git a/SimPPS/RPDigiProducer/plugins/RPPileUpSignals.h b/SimPPS/RPDigiProducer/plugins/RPPileUpSignals.h new file mode 100644 index 0000000000000..343cb66f2c7aa --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPPileUpSignals.h @@ -0,0 +1,24 @@ +#ifndef SimPPS_RPDigiProducer_RP_PILE_UP_SIGNALS_H +#define SimPPS_RPDigiProducer_RP_PILE_UP_SIGNALS_H + +#include +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" + +class RPPileUpSignals { +public: + RPPileUpSignals(const edm::ParameterSet ¶ms, RPDetId det_id); + void reset(); + void add(const simromanpot::strip_charge_map &charge_induced, int PSimHitIndex); + inline const simromanpot::strip_charge_map &dumpSignal() { return the_strip_charge_piled_up_map_; } + inline simromanpot::strip_charge_map_links_type &dumpLinks() { return the_strip_charge_piled_up_map_links_; } + +private: + simromanpot::strip_charge_map the_strip_charge_piled_up_map_; + simromanpot::strip_charge_map_links_type the_strip_charge_piled_up_map_links_; + bool links_persistence_; + RPDetId det_id_; + bool verbosity_; +}; + +#endif //SimPPS_RPDigiProducer_RP_PILE_UP_SIGNALS_H diff --git a/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc b/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc new file mode 100644 index 0000000000000..6eb7ef2863fad --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc @@ -0,0 +1,45 @@ +#include "SimPPS/RPDigiProducer/plugins/RPVFATSimulator.h" +#include "Geometry/VeryForwardRPTopology/interface/RPTopology.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include +#include "TRandom.h" +#include + +RPVFATSimulator::RPVFATSimulator(const edm::ParameterSet ¶ms, RPDetId det_id) : params_(params), det_id_(det_id) { + threshold_ = params.getParameter("RPVFATThreshold"); + dead_strip_probability_ = params.getParameter("RPDeadStripProbability"); + dead_strips_simulation_on_ = params.getParameter("RPDeadStripSimulationOn"); + strips_no_ = RPTopology().DetStripNo(); + verbosity_ = params.getParameter("RPVerbosity"); + links_persistence_ = params.getParameter("RPDigiSimHitRelationsPresistence"); +} + +void RPVFATSimulator::ConvertChargeToHits(const simromanpot::strip_charge_map &signals, + simromanpot::strip_charge_map_links_type &theSignalProvenance, + std::vector &output_digi, + simromanpot::DigiPrimaryMapType &output_digi_links) { + for (simromanpot::strip_charge_map::const_iterator i = signals.begin(); i != signals.end(); ++i) { + //one threshold per hybrid + unsigned short strip_no = i->first; + if (i->second > threshold_ && (!dead_strips_simulation_on_ || dead_strips_.find(strip_no) == dead_strips_.end())) { + output_digi.push_back(TotemRPDigi(strip_no)); + if (links_persistence_) { + output_digi_links.push_back(theSignalProvenance[strip_no]); + if (verbosity_) { + edm::LogInfo("RPVFatSimulator") << " digi links size=" << theSignalProvenance[strip_no].size() << "\n"; + for (unsigned int u = 0; u < theSignalProvenance[strip_no].size(); ++u) { + edm::LogInfo("RPVFatSimulator") + << " digi: particle=" << theSignalProvenance[strip_no][u].first + << " energy [electrons]=" << theSignalProvenance[strip_no][u].second << "\n"; + } + } + } + } + } + + if (verbosity_) { + for (unsigned int i = 0; i < output_digi.size(); ++i) { + edm::LogInfo("RPVFATSimulator") << output_digi[i].getStripNumber() << "\n"; + } + } +} diff --git a/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.h b/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.h new file mode 100644 index 0000000000000..1aa7281d9b379 --- /dev/null +++ b/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.h @@ -0,0 +1,33 @@ +#ifndef RP_VFAT_SIMULATION_H +#define RP_VFAT_SIMULATION_H + +#include + +#include "DataFormats/CTPPSDigi/interface/TotemRPDigi.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "SimPPS/RPDigiProducer/interface/RPSimTypes.h" + +class RPVFATSimulator { +public: + RPVFATSimulator(const edm::ParameterSet ¶ms, RPDetId det_id); + void ConvertChargeToHits(const simromanpot::strip_charge_map &signals, + simromanpot::strip_charge_map_links_type &theSignalProvenance, + std::vector &output_digi, + simromanpot::DigiPrimaryMapType &output_digi_links); + +private: + typedef std::set > dead_strip_set; + const edm::ParameterSet ¶ms_; + RPDetId det_id_; + double dead_strip_probability_; + bool dead_strips_simulation_on_; + dead_strip_set dead_strips_; + int verbosity_; + + unsigned short strips_no_; + double threshold_; + bool links_persistence_; +}; + +#endif diff --git a/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py b/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py new file mode 100644 index 0000000000000..9723c0acfdd70 --- /dev/null +++ b/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py @@ -0,0 +1,101 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("RPDigiProducerTest") + +# Specify the maximum events to simulate +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(10) +) + +# Configure the output module (save the result in a file) +# Configure the output module +process.o1 = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string('file:RPDigiProducerTest_output.root') +) + + +process.load("SimGeneral.HepPDTESSource.pdt_cfi") + + +# Configure if you want to detail or simple log information. +# LoggerMax -- detail log info output including: errors.log, warnings.log, infos.log, debugs.log +# LoggerMin -- simple log info output to the standard output (e.g. screen) +process.load("Configuration.TotemCommon.LoggerMin_cfi") + + +################## STEP 1 +process.source = cms.Source("EmptySource") + +################## STEP 2 - process.generator + +# Use random number generator service +process.load("Configuration.TotemCommon.RandomNumbers_cfi") + +# Monte Carlo gun - elastic specific +energy = "7000" +import IOMC.Elegent.ElegentSource_cfi +process.generator = IOMC.Elegent.ElegentSource_cfi.generator +process.generator.fileName = IOMC.Elegent.ElegentSource_cfi.ElegentDefaultFileName(energy) + +# particle generator paramteres +process.generator.t_min = '6E-2' # beta* specific +process.generator.t_max = '6E-1' # beta* specific +energy = "1180" + + +################## STEP 3 process.SmearingGenerator + +# declare optics parameters +# process.load("Configuration.TotemOpticsConfiguration.OpticsConfig_7000GeV_90_cfi") + +# Smearing +process.load("IOMC.SmearingGenerator.SmearingGenerator_cfi") + +################## STEP 4 process.OptInfo + +process.OptInfo = cms.EDAnalyzer("OpticsInformation") + +process.load("Configuration.TotemOpticsConfiguration.OpticsConfig_1180GeV_11_cfi") + +################## STEP 5 process.*process.g4SimHits + +# Geometry - beta* specific +# process.load("Configuration.TotemCommon.geometryRP_cfi") +# process.XMLIdealGeometryESSource.geomXMLFiles.append('Geometry/TotemRPData/data/RP_Beta_90_150_out/RP_Dist_Beam_Cent.xml') + +# Magnetic Field, by default we have 3.8T +process.load("Configuration.StandardSequences.MagneticField_cff") + +# G4 simulation & proton transport +process.load("Configuration.TotemCommon.g4SimHits_cfi") +#process.g4SimHits.Physics.BeamProtTransportSetup = process.BeamProtTransportSetup +process.g4SimHits.Generator.HepMCProductLabel = 'generator' # The input source for G4 module is connected to "process.source". + +process.load("Configuration.TotemCommon.geometryRP_cfi") +process.XMLIdealGeometryESSource.geomXMLFiles.append('Geometry/TotemRPData/data/RP_1180_Beta_11_220/RP_Dist_Beam_Cent.xml') + +process.g4SimHits.Physics.BeamProtTransportSetup = process.BeamProtTransportSetup + + + +################## STEP 6 process.mix*process.RPSiDetDigitizer + +# No pile up for the mixing module +process.load("SimGeneral.MixingModule.mixNoPU_cfi") + +########################### DIGI+RECO RP ########################################## + +# process.load("SimPPS.RPDigiProducer.RPSiDetConf_cfi") + + + +process.p1 = cms.Path(process.generator + *process.SmearingGenerator + *process.g4SimHits +# *process.mix +# *process.RPSiDetDigitizer + ) + +process.outpath = cms.EndPath(process.o1) + +# print process.dumpConfig() From e24963c7cd4d8e89d0e77157881c091caee3cb51 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Tue, 30 Jul 2019 13:39:27 +0200 Subject: [PATCH 2/8] Moving to LocalPoint acessor methods --- SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h | 5 +---- SimPPS/RPDigiProducer/interface/RPSignalPoint.h | 4 +--- .../plugins/RPLinearChargeCollectionDrifter.cc | 4 ++-- SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc | 4 ++-- .../RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc | 2 +- 5 files changed, 7 insertions(+), 12 deletions(-) diff --git a/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h b/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h index c29cb5682884e..6e00ab8af8394 100644 --- a/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h +++ b/SimPPS/RPDigiProducer/interface/RPEnergyDepositUnit.h @@ -13,11 +13,8 @@ class RPEnergyDepositUnit { RPEnergyDepositUnit() : energy_(0), position_(0, 0, 0) {} RPEnergyDepositUnit(double energy, double x, double y, double z) : energy_(energy), position_(x, y, z) {} RPEnergyDepositUnit(double energy, const Local3DPoint& position) : energy_(energy), position_(position) {} - inline double X() const { return position_.x(); } - inline double Y() const { return position_.y(); } - inline double Z() const { return position_.z(); } inline double Energy() const { return energy_; } - inline Local3DPoint Position() const { return position_; } + inline const Local3DPoint& Position() const { return position_; } inline void setEnergy(double e) { energy_ = e; } inline void setPosition(Local3DPoint p) { position_ = p; } diff --git a/SimPPS/RPDigiProducer/interface/RPSignalPoint.h b/SimPPS/RPDigiProducer/interface/RPSignalPoint.h index 87a17cd408dcb..53262653fb7ca 100644 --- a/SimPPS/RPDigiProducer/interface/RPSignalPoint.h +++ b/SimPPS/RPDigiProducer/interface/RPSignalPoint.h @@ -13,9 +13,7 @@ class RPSignalPoint { RPSignalPoint(double x, double y, double s, double charge) : pos_(x, y), sigma_(s), charge_(charge) {} - inline LocalPoint Position() const { return pos_; } - inline double X() const { return pos_.x(); } - inline double Y() const { return pos_.y(); } + inline const LocalPoint& Position() const { return pos_; } inline double Sigma() const { return sigma_; } inline double Charge() const { return charge_; } diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc index ad0f8b70809a6..f5ba5ac6afcb8 100644 --- a/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc @@ -17,8 +17,8 @@ simromanpot::charge_induced_on_surface RPLinearChargeCollectionDrifter::Drift( simromanpot::charge_induced_on_surface temp_; temp_.resize(energy_deposition.size()); for (unsigned int i = 0; i < energy_deposition.size(); i++) { - temp_[i].setPosition(LocalPoint(energy_deposition[i].X(), energy_deposition[i].Y())); - temp_[i].setSigma(getSigma(energy_deposition[i].Z())); //befor charge_cloud_sigma_ used, now a vector of sigmas; + temp_[i].setPosition(LocalPoint(energy_deposition[i].Position().x(), energy_deposition[i].Position().y())); + temp_[i].setSigma(getSigma(energy_deposition[i].Position().z())); //befor charge_cloud_sigma_ used, now a vector of sigmas; temp_[i].setCharge(energy_deposition[i].Energy() / GeV_per_electron_); if (verbosity_) { edm::LogInfo("RPLinearChargeCollectionDrifter") diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc index d1d4db43fe3f6..d35266f97cc81 100644 --- a/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc @@ -73,8 +73,8 @@ simromanpot::energy_path_distribution RPLinearChargeDivider::divide(const PSimHi double sum = 0; for (unsigned int i = 0; i < the_energy_path_distribution_.size(); i++) { edm::LogInfo("RPLinearChargeDivider") - << the_energy_path_distribution_[i].X() << " " << the_energy_path_distribution_[i].Y() << " " - << the_energy_path_distribution_[i].Z() << " " << the_energy_path_distribution_[i].Energy() << "\n"; + << the_energy_path_distribution_[i].Position().x() << " " << the_energy_path_distribution_[i].Position().y() << " " + << the_energy_path_distribution_[i].Position().z() << " " << the_energy_path_distribution_[i].Energy() << "\n"; sum += the_energy_path_distribution_[i].Energy(); } edm::LogInfo("RPLinearChargeDivider") << "energy dep. sum=" << sum << "\n"; diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc b/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc index 6a49a44f12162..05fe899a85ab5 100644 --- a/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc +++ b/SimPPS/RPDigiProducer/plugins/RPLinearInduceChargeOnStrips.cc @@ -22,7 +22,7 @@ simromanpot::strip_charge_map RPLinearInduceChargeOnStrips::Induce( for (simromanpot::charge_induced_on_surface::const_iterator i = charge_map.begin(); i != charge_map.end(); ++i) { double hit_pos; std::vector relevant_strips = - theRPDetTopology.GetStripsInvolved((*i).X(), (*i).Y(), (*i).Sigma(), hit_pos); + theRPDetTopology.GetStripsInvolved((*i).Position().x(), (*i).Position().y(), (*i).Sigma(), hit_pos); if (verbosity_) { edm::LogInfo("RPLinearInduceChargeOnStrips ") << det_id_ << " : relevant_strips" << relevant_strips.size() << "\n"; From 6a4b97018139d5aa7ddfaf6d1d2afe97c91581ef Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Tue, 30 Jul 2019 13:58:09 +0200 Subject: [PATCH 3/8] implementing a simple line break triggered by code-format --- .../plugins/RPLinearChargeCollectionDrifter.cc | 3 ++- SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc | 5 +++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc index f5ba5ac6afcb8..7e770a8ea6ec2 100644 --- a/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeCollectionDrifter.cc @@ -18,7 +18,8 @@ simromanpot::charge_induced_on_surface RPLinearChargeCollectionDrifter::Drift( temp_.resize(energy_deposition.size()); for (unsigned int i = 0; i < energy_deposition.size(); i++) { temp_[i].setPosition(LocalPoint(energy_deposition[i].Position().x(), energy_deposition[i].Position().y())); - temp_[i].setSigma(getSigma(energy_deposition[i].Position().z())); //befor charge_cloud_sigma_ used, now a vector of sigmas; + temp_[i].setSigma( + getSigma(energy_deposition[i].Position().z())); //befor charge_cloud_sigma_ used, now a vector of sigmas; temp_[i].setCharge(energy_deposition[i].Energy() / GeV_per_electron_); if (verbosity_) { edm::LogInfo("RPLinearChargeCollectionDrifter") diff --git a/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc index d35266f97cc81..05e42fee63ebd 100644 --- a/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc +++ b/SimPPS/RPDigiProducer/plugins/RPLinearChargeDivider.cc @@ -73,8 +73,9 @@ simromanpot::energy_path_distribution RPLinearChargeDivider::divide(const PSimHi double sum = 0; for (unsigned int i = 0; i < the_energy_path_distribution_.size(); i++) { edm::LogInfo("RPLinearChargeDivider") - << the_energy_path_distribution_[i].Position().x() << " " << the_energy_path_distribution_[i].Position().y() << " " - << the_energy_path_distribution_[i].Position().z() << " " << the_energy_path_distribution_[i].Energy() << "\n"; + << the_energy_path_distribution_[i].Position().x() << " " << the_energy_path_distribution_[i].Position().y() + << " " << the_energy_path_distribution_[i].Position().z() << " " << the_energy_path_distribution_[i].Energy() + << "\n"; sum += the_energy_path_distribution_[i].Energy(); } edm::LogInfo("RPLinearChargeDivider") << "energy dep. sum=" << sum << "\n"; From a0fd21157692bfdd031b39e83b9510be1d4e2186 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Wed, 7 Aug 2019 19:48:39 +0200 Subject: [PATCH 4/8] Moving to CLHEP random generator only and changing TMath::Erfc to std::erfc --- .../RPDigiProducer/plugins/RPDetDigitizer.cc | 2 +- .../plugins/RPGaussianTailNoiseAdder.cc | 23 ++++++++++++------- .../plugins/RPGaussianTailNoiseAdder.h | 3 ++- 3 files changed, 18 insertions(+), 10 deletions(-) diff --git a/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc index f0d3b3150d6a9..039f61adee81a 100644 --- a/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc +++ b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc @@ -20,7 +20,7 @@ RPDetDigitizer::RPDetDigitizer(const edm::ParameterSet ¶ms, links_persistence_ = params.getParameter("RPDigiSimHitRelationsPresistence"); theRPGaussianTailNoiseAdder = - std::make_unique(numStrips_, theNoiseInElectrons, theStripThresholdInE, verbosity_); + std::make_unique(numStrips_, theNoiseInElectrons, theStripThresholdInE, eng, verbosity_); theRPPileUpSignals = std::make_unique(params, det_id_); theRPVFATSimulator = std::make_unique(params, det_id_); theRPHitChargeConverter = std::make_unique(params, eng, det_id_); diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc index 1d74a4ba16972..b7c67b7eb83b3 100644 --- a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc @@ -1,19 +1,26 @@ #include "SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include -#include "TMath.h" -#include "TRandom.h" #include "CLHEP/Random/RandGauss.h" +#include "CLHEP/Random/RandFlat.h" +#include "CLHEP/Random/RandBinomial.h" + +#include using namespace std; RPGaussianTailNoiseAdder::RPGaussianTailNoiseAdder(int numStrips, double theNoiseInElectrons, double theStripThresholdInE, + CLHEP::HepRandomEngine &eng, int verbosity) - : numStrips_(numStrips), theNoiseInElectrons(theNoiseInElectrons), theStripThresholdInE(theStripThresholdInE) { - verbosity_ = verbosity; - strips_above_threshold_prob_ = TMath::Erfc(theStripThresholdInE / sqrt(2.0) / theNoiseInElectrons) / 2; + : numStrips_(numStrips), + theNoiseInElectrons(theNoiseInElectrons), + theStripThresholdInE(theStripThresholdInE), + rndEngine_(eng), + verbosity_(verbosity) +{ + strips_above_threshold_prob_ = std::erfc(theStripThresholdInE / sqrt(2.0) / theNoiseInElectrons) / 2; } simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanpot::strip_charge_map &theSignal) { @@ -21,17 +28,17 @@ simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanp // noise on strips with signal: for (simromanpot::strip_charge_map::const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) { - double noise = CLHEP::RandGauss::shoot(0.0, theNoiseInElectrons); + double noise = CLHEP::RandGauss::shoot(&(rndEngine_),0.0, theNoiseInElectrons); the_strip_charge_map[i->first] = i->second + noise; if (verbosity_) edm::LogInfo("RPDigiProducer") << "noise added to signal strips: " << noise << "\n"; } // noise on the other strips - int strips_no_above_threshold = gRandom->Binomial(numStrips_, strips_above_threshold_prob_); + int strips_no_above_threshold = CLHEP::RandBinomial::shoot(&(rndEngine_),(long)numStrips_, strips_above_threshold_prob_); for (int j = 0; j < strips_no_above_threshold; j++) { - int strip = gRandom->Integer(numStrips_); + int strip = CLHEP::RandFlat::shootInt(numStrips_); if (the_strip_charge_map[strip] == 0) { the_strip_charge_map[strip] = 2 * theStripThresholdInE; //only binary decision later, no need to simulate the noise precisely, diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h index 31a5a7c385ab4..b844d6b615ee8 100644 --- a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h @@ -6,7 +6,7 @@ class RPGaussianTailNoiseAdder { public: - RPGaussianTailNoiseAdder(int numStrips, double theNoiseInElectrons, double theStripThresholdInE, int verbosity); + RPGaussianTailNoiseAdder(int numStrips, double theNoiseInElectrons, double theStripThresholdInE, CLHEP::HepRandomEngine &eng, int verbosity); simromanpot::strip_charge_map addNoise(const simromanpot::strip_charge_map &theSignal); private: @@ -14,6 +14,7 @@ class RPGaussianTailNoiseAdder { double theNoiseInElectrons; double theStripThresholdInE; double strips_above_threshold_prob_; + CLHEP::HepRandomEngine& rndEngine_; int verbosity_; }; From ef1d8333a3cb79931545d2eb8bcd5781a67b25a7 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Wed, 7 Aug 2019 20:26:54 +0200 Subject: [PATCH 5/8] applying scram code-format changes - identation --- .../RPDigiProducer/plugins/RPDetDigitizer.cc | 4 ++-- .../plugins/RPGaussianTailNoiseAdder.cc | 23 ++++++++----------- .../plugins/RPGaussianTailNoiseAdder.h | 8 +++++-- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc index 039f61adee81a..b751da0e6fc24 100644 --- a/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc +++ b/SimPPS/RPDigiProducer/plugins/RPDetDigitizer.cc @@ -19,8 +19,8 @@ RPDetDigitizer::RPDetDigitizer(const edm::ParameterSet ¶ms, misalignment_simulation_on_ = params.getParameter("RPDisplacementOn"); links_persistence_ = params.getParameter("RPDigiSimHitRelationsPresistence"); - theRPGaussianTailNoiseAdder = - std::make_unique(numStrips_, theNoiseInElectrons, theStripThresholdInE, eng, verbosity_); + theRPGaussianTailNoiseAdder = std::make_unique( + numStrips_, theNoiseInElectrons, theStripThresholdInE, eng, verbosity_); theRPPileUpSignals = std::make_unique(params, det_id_); theRPVFATSimulator = std::make_unique(params, det_id_); theRPHitChargeConverter = std::make_unique(params, eng, det_id_); diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc index b7c67b7eb83b3..6c82e89159b38 100644 --- a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc @@ -9,17 +9,13 @@ using namespace std; -RPGaussianTailNoiseAdder::RPGaussianTailNoiseAdder(int numStrips, - double theNoiseInElectrons, - double theStripThresholdInE, - CLHEP::HepRandomEngine &eng, - int verbosity) - : numStrips_(numStrips), - theNoiseInElectrons(theNoiseInElectrons), - theStripThresholdInE(theStripThresholdInE), - rndEngine_(eng), - verbosity_(verbosity) -{ +RPGaussianTailNoiseAdder::RPGaussianTailNoiseAdder( + int numStrips, double theNoiseInElectrons, double theStripThresholdInE, CLHEP::HepRandomEngine &eng, int verbosity) + : numStrips_(numStrips), + theNoiseInElectrons(theNoiseInElectrons), + theStripThresholdInE(theStripThresholdInE), + rndEngine_(eng), + verbosity_(verbosity) { strips_above_threshold_prob_ = std::erfc(theStripThresholdInE / sqrt(2.0) / theNoiseInElectrons) / 2; } @@ -28,14 +24,15 @@ simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanp // noise on strips with signal: for (simromanpot::strip_charge_map::const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) { - double noise = CLHEP::RandGauss::shoot(&(rndEngine_),0.0, theNoiseInElectrons); + double noise = CLHEP::RandGauss::shoot(&(rndEngine_), 0.0, theNoiseInElectrons); the_strip_charge_map[i->first] = i->second + noise; if (verbosity_) edm::LogInfo("RPDigiProducer") << "noise added to signal strips: " << noise << "\n"; } // noise on the other strips - int strips_no_above_threshold = CLHEP::RandBinomial::shoot(&(rndEngine_),(long)numStrips_, strips_above_threshold_prob_); + int strips_no_above_threshold = + CLHEP::RandBinomial::shoot(&(rndEngine_), (long)numStrips_, strips_above_threshold_prob_); for (int j = 0; j < strips_no_above_threshold; j++) { int strip = CLHEP::RandFlat::shootInt(numStrips_); diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h index b844d6b615ee8..c75217c209fb5 100644 --- a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h @@ -6,7 +6,11 @@ class RPGaussianTailNoiseAdder { public: - RPGaussianTailNoiseAdder(int numStrips, double theNoiseInElectrons, double theStripThresholdInE, CLHEP::HepRandomEngine &eng, int verbosity); + RPGaussianTailNoiseAdder(int numStrips, + double theNoiseInElectrons, + double theStripThresholdInE, + CLHEP::HepRandomEngine &eng, + int verbosity); simromanpot::strip_charge_map addNoise(const simromanpot::strip_charge_map &theSignal); private: @@ -14,7 +18,7 @@ class RPGaussianTailNoiseAdder { double theNoiseInElectrons; double theStripThresholdInE; double strips_above_threshold_prob_; - CLHEP::HepRandomEngine& rndEngine_; + CLHEP::HepRandomEngine &rndEngine_; int verbosity_; }; From ffe4ffa9d449c5dbbaf0e4351119d846f38e2533 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Fri, 9 Aug 2019 14:59:58 +0200 Subject: [PATCH 6/8] Removing BuildFile dependency with Reco objetcs, replacing by the needed pacakges --- SimPPS/RPDigiProducer/plugins/BuildFile.xml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SimPPS/RPDigiProducer/plugins/BuildFile.xml b/SimPPS/RPDigiProducer/plugins/BuildFile.xml index 4d308dd305519..5f117c4ade35e 100755 --- a/SimPPS/RPDigiProducer/plugins/BuildFile.xml +++ b/SimPPS/RPDigiProducer/plugins/BuildFile.xml @@ -1,6 +1,8 @@ - - + + + + From 2c003a12e920ff7a0c0b60309f874883ebea1542 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Mon, 12 Aug 2019 19:17:54 +0200 Subject: [PATCH 7/8] Removing a call to the random generator without the initialized engine that was left behind; adapting the test config to the CMS standard --- .../plugins/RPGaussianTailNoiseAdder.cc | 2 +- .../RPDigiProducer/plugins/RPVFATSimulator.cc | 1 - .../test/RPDigiProducerTest_cfg.py | 104 ++++++++++-------- 3 files changed, 61 insertions(+), 46 deletions(-) diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc index 6c82e89159b38..78aa1f88a2495 100644 --- a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc @@ -35,7 +35,7 @@ simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanp CLHEP::RandBinomial::shoot(&(rndEngine_), (long)numStrips_, strips_above_threshold_prob_); for (int j = 0; j < strips_no_above_threshold; j++) { - int strip = CLHEP::RandFlat::shootInt(numStrips_); + int strip = CLHEP::RandFlat::shootInt(&(rndEngine_),numStrips_); if (the_strip_charge_map[strip] == 0) { the_strip_charge_map[strip] = 2 * theStripThresholdInE; //only binary decision later, no need to simulate the noise precisely, diff --git a/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc b/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc index 6eb7ef2863fad..24244a3e83b72 100644 --- a/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc +++ b/SimPPS/RPDigiProducer/plugins/RPVFATSimulator.cc @@ -2,7 +2,6 @@ #include "Geometry/VeryForwardRPTopology/interface/RPTopology.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include -#include "TRandom.h" #include RPVFATSimulator::RPVFATSimulator(const edm::ParameterSet ¶ms, RPDetId det_id) : params_(params), det_id_(det_id) { diff --git a/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py b/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py index 9723c0acfdd70..6f3a834d5845d 100644 --- a/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py +++ b/SimPPS/RPDigiProducer/test/RPDigiProducerTest_cfg.py @@ -1,4 +1,5 @@ import FWCore.ParameterSet.Config as cms +import math process = cms.Process("RPDigiProducerTest") @@ -17,66 +18,72 @@ process.load("SimGeneral.HepPDTESSource.pdt_cfi") -# Configure if you want to detail or simple log information. -# LoggerMax -- detail log info output including: errors.log, warnings.log, infos.log, debugs.log -# LoggerMin -- simple log info output to the standard output (e.g. screen) -process.load("Configuration.TotemCommon.LoggerMin_cfi") - - ################## STEP 1 process.source = cms.Source("EmptySource") ################## STEP 2 - process.generator - -# Use random number generator service -process.load("Configuration.TotemCommon.RandomNumbers_cfi") - -# Monte Carlo gun - elastic specific -energy = "7000" -import IOMC.Elegent.ElegentSource_cfi -process.generator = IOMC.Elegent.ElegentSource_cfi.generator -process.generator.fileName = IOMC.Elegent.ElegentSource_cfi.ElegentDefaultFileName(energy) - -# particle generator paramteres -process.generator.t_min = '6E-2' # beta* specific -process.generator.t_max = '6E-1' # beta* specific -energy = "1180" - +process.load('Configuration.StandardSequences.Generator_cff') +process.load('GeneratorInterface.Core.genFilterSummary_cff') +process.load("IOMC.RandomEngine.IOMC_cff") +process.RandomNumberGeneratorService.generator.initialSeed = 456789 +process.RandomNumberGeneratorService.g4SimHits.initialSeed = 9876 +process.RandomNumberGeneratorService.VtxSmeared.initialSeed = 123456789 +process.RandomNumberGeneratorService.LHCTransport.engineName = cms.untracked.string('TRandom3') +phi_min = -math.pi +phi_max = math.pi +t_min = 0. +t_max = 2.0 +xi_min = 0.01 +xi_max = 0.2 +ecms = 13000. +process.generator = cms.EDProducer("RandomtXiGunProducer", + PGunParameters = cms.PSet( + PartID = cms.vint32(2212), + MinPhi = cms.double(phi_min), + MaxPhi = cms.double(phi_max), + ECMS = cms.double(ecms), + Mint = cms.double(t_min), + Maxt = cms.double(t_max), + MinXi = cms.double(xi_min), + MaxXi = cms.double(xi_max) + ), + Verbosity = cms.untracked.int32(0), + psethack = cms.string('single protons'), + FireBackward = cms.bool(True), + FireForward = cms.bool(True), + firstRun = cms.untracked.uint32(1), + ) + +process.ProductionFilterSequence = cms.Sequence(process.generator) ################## STEP 3 process.SmearingGenerator # declare optics parameters -# process.load("Configuration.TotemOpticsConfiguration.OpticsConfig_7000GeV_90_cfi") # Smearing -process.load("IOMC.SmearingGenerator.SmearingGenerator_cfi") +process.load('IOMC.EventVertexGenerators.VtxSmearedRealistic25ns13TeVEarly2017Collision_cfi') ################## STEP 4 process.OptInfo -process.OptInfo = cms.EDAnalyzer("OpticsInformation") - -process.load("Configuration.TotemOpticsConfiguration.OpticsConfig_1180GeV_11_cfi") ################## STEP 5 process.*process.g4SimHits -# Geometry - beta* specific -# process.load("Configuration.TotemCommon.geometryRP_cfi") -# process.XMLIdealGeometryESSource.geomXMLFiles.append('Geometry/TotemRPData/data/RP_Beta_90_150_out/RP_Dist_Beam_Cent.xml') - # Magnetic Field, by default we have 3.8T +process.load('Configuration.StandardSequences.SimIdeal_cff') process.load("Configuration.StandardSequences.MagneticField_cff") +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase1_2017_realistic', '') +process.load('PhysicsTools.HepMCCandAlgos.genParticles_cfi') # G4 simulation & proton transport -process.load("Configuration.TotemCommon.g4SimHits_cfi") -#process.g4SimHits.Physics.BeamProtTransportSetup = process.BeamProtTransportSetup -process.g4SimHits.Generator.HepMCProductLabel = 'generator' # The input source for G4 module is connected to "process.source". - -process.load("Configuration.TotemCommon.geometryRP_cfi") -process.XMLIdealGeometryESSource.geomXMLFiles.append('Geometry/TotemRPData/data/RP_1180_Beta_11_220/RP_Dist_Beam_Cent.xml') -process.g4SimHits.Physics.BeamProtTransportSetup = process.BeamProtTransportSetup +from Geometry.VeryForwardGeometry.geometryPPS_CMSxz_fromDD_2017_cfi import XMLIdealGeometryESSource_CTPPS +process.XMLIdealGeometryESSource = XMLIdealGeometryESSource_CTPPS.clone() +process.load("SimG4Core.Application.g4SimHits_cfi") +process.g4SimHits.Generator.HepMCProductLabel = 'LHCTransport' # The input source for G4 module is connected to "process.source". ################## STEP 6 process.mix*process.RPSiDetDigitizer @@ -86,16 +93,25 @@ ########################### DIGI+RECO RP ########################################## # process.load("SimPPS.RPDigiProducer.RPSiDetConf_cfi") +process.generation_step = cms.Path(process.pgen) +process.simulation_step = cms.Path(process.psim) +process.g4Simhits_step = cms.Path(process.g4SimHits) - -process.p1 = cms.Path(process.generator - *process.SmearingGenerator - *process.g4SimHits -# *process.mix -# *process.RPSiDetDigitizer - ) +process.simulation_step = cms.Path(process.psim) +process.g4Simhits_step = cms.Path(process.g4SimHits) +process.genfiltersummary_step = cms.EndPath(process.genFilterSummary) process.outpath = cms.EndPath(process.o1) +process.schedule = cms.Schedule(process.generation_step,process.genfiltersummary_step,process.g4Simhits_step,process.outpath) + +# filter all path with the production filter sequence +for path in process.paths: + getattr(process,path)._seq = process.ProductionFilterSequence * getattr(process,path)._seq + + # print process.dumpConfig() +from SimPPS.PPSSimTrackProducer.SimTrackProducerForFullSim_cff import customise +process = customise(process) + From be4b814d3d68c2f50356992c8910c3559380fc06 Mon Sep 17 00:00:00 2001 From: Luiz Mundim Date: Mon, 12 Aug 2019 19:44:57 +0200 Subject: [PATCH 8/8] Adding a space before an argument since code-format failed --- SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc index 78aa1f88a2495..12820a193b6b1 100644 --- a/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc +++ b/SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc @@ -35,7 +35,7 @@ simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanp CLHEP::RandBinomial::shoot(&(rndEngine_), (long)numStrips_, strips_above_threshold_prob_); for (int j = 0; j < strips_no_above_threshold; j++) { - int strip = CLHEP::RandFlat::shootInt(&(rndEngine_),numStrips_); + int strip = CLHEP::RandFlat::shootInt(&(rndEngine_), numStrips_); if (the_strip_charge_map[strip] == 0) { the_strip_charge_map[strip] = 2 * theStripThresholdInE; //only binary decision later, no need to simulate the noise precisely,