RPDigiProducer/plugins/RPDigiProducer.cc

0001 // user include files
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012
0013 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0014 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0015 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0016 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0017
0018 #include "DataFormats/CTPPSDigi/interface/TotemRPDigi.h"
0019 #include "DataFormats/CTPPSDigi/interface/TotemRPDigi.h"
0020 #include "DataFormats/Common/interface/DetSetVector.h"
0021 #include "DataFormats/CTPPSDetId/interface/TotemRPDetId.h"
0022 #include "DataFormats/Common/interface/DetSet.h"
0023
0024 //Random Number
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0027
0028 #include "CondFormats/PPSObjects/interface/TotemAnalysisMask.h"
0029 #include "CondFormats/DataRecord/interface/TotemReadoutRcd.h"
0030 #include "SimPPS/RPDigiProducer/interface/RPSimTypes.h"
0031 #include "SimPPS/RPDigiProducer/plugins/RPDetDigitizer.h"
0032 #include "SimPPS/RPDigiProducer/plugins/DeadChannelsManager.h"
0033
0034 #include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h"
0035 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0036
0037 // system include files
0038 #include <memory>
0039 #include <vector>
0040 #include <map>
0041 #include <string>
0042 #include <iostream>
0043 #include <cstdlib>  // I need it for random numbers
0044
0045 // user include files
0046
0047 //
0048 // class decleration
0049 //
0050
0051 namespace CLHEP {
0052   class HepRandomEngine;
0053 }
0054
0055 class RPDigiProducer : public edm::stream::EDProducer<> {
0056 public:
0057   explicit RPDigiProducer(const edm::ParameterSet&);
0058   ~RPDigiProducer() override = default;
0059
0060   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0061
0062 private:
0063   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0064   void produce(edm::Event&, const edm::EventSetup&) override;
0065
0066   edm::DetSet<TotemRPDigi> convertRPStripDetSet(const edm::DetSet<TotemRPDigi>&);
0067
0068   // ----------member data ---------------------------
0069   std::vector<std::string> RP_hit_containers_;
0070   typedef std::map<unsigned int, std::vector<PSimHit>> simhit_map;
0071   typedef simhit_map::iterator simhit_map_iterator;
0072
0073   edm::ParameterSet conf_;
0074   std::map<RPDetId, std::unique_ptr<RPDetDigitizer>> theAlgoMap;
0075
0076   CLHEP::HepRandomEngine* rndEngine_ = nullptr;
0077   int verbosity_;
0078
0079   /**
0080        * this variable answers the question whether given channel is dead or not
0081        */
0082   DeadChannelsManager deadChannelsManager;
0083   /**
0084        * this variable indicates whether we take into account dead channels or simulate as if all
0085        * channels work ok (by default we do not simulate dead channels)
0086        */
0087   bool simulateDeadChannels;
0088
0089   edm::EDGetTokenT<CrossingFrame<PSimHit>> tokenCrossingFrameTotemRP;
0090   edm::ESGetToken<TotemAnalysisMask, TotemReadoutRcd> tokenAnalysisMask;
0091   edm::ESGetToken<CTPPSRPAlignmentCorrectionsData, VeryForwardMisalignedGeometryRecord> alignmentToken;
0092   edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> geomToken;
0093 };
0094
0095 RPDigiProducer::RPDigiProducer(const edm::ParameterSet& conf) : conf_(conf) {
0096   //now do what ever other initialization is needed
0097   produces<edm::DetSetVector<TotemRPDigi>>();
0098
0099   // register data to consume
0100   tokenCrossingFrameTotemRP = consumes<CrossingFrame<PSimHit>>(edm::InputTag("mix", "g4SimHitsTotemHitsRP", ""));
0101
0102   RP_hit_containers_ = conf.getParameter<std::vector<std::string>>("ROUList");
0103   verbosity_ = conf.getParameter<int>("RPVerbosity");
0104
0105   simulateDeadChannels = false;
0106   if (conf.exists(
0107           "simulateDeadChannels")) {  //check if "simulateDeadChannels" variable is defined in configuration file
0108     simulateDeadChannels = conf.getParameter<bool>("simulateDeadChannels");
0109   }
0110   if (simulateDeadChannels) {
0111     tokenAnalysisMask = esConsumes();
0112   }
0113   alignmentToken = esConsumes();
0114   geomToken = esConsumes();
0115 }
0116
0117 //
0118 // member functions
0119 //
0120
0121 // ------------ method called to produce the data  ------------
0122 void RPDigiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0123   using namespace edm;
0124
0125   // initialize random engine
0126   if (!rndEngine_) {
0127     Service<RandomNumberGenerator> rng;
0128     if (!rng.isAvailable()) {
0129       throw cms::Exception("Configuration")
0130           << "This class requires the RandomNumberGeneratorService\n"
0131              "which is not present in the configuration file.  You must add the service\n"
0132              "in the configuration file or remove the modules that require it.";
0133     }
0134     rndEngine_ = &(rng->getEngine(iEvent.streamID()));
0135   }
0136
0137   // Step A: Get Inputs
0138   edm::Handle<CrossingFrame<PSimHit>> cf;
0139   iEvent.getByLabel("mix", "g4SimHitsTotemHitsRP", cf);
0140
0141   if (verbosity_) {
0142     edm::LogInfo("RPDigiProducer") << "\n\n=================== Starting SimHit access"
0143                                    << "  ==================="
0144                                    << "\n";
0145
0146     MixCollection<PSimHit> col{cf.product(), std::pair(-0, 0)};
0147     MixCollection<PSimHit>::iterator cfi;
0148     int count = 0;
0149     for (cfi = col.begin(); cfi != col.end(); cfi++) {
0150       edm::LogInfo("RPDigiProducer") << " Hit " << count << " has tof " << cfi->timeOfFlight() << " trackid "
0151                                      << cfi->trackId() << " bunchcr " << cfi.bunch() << " trigger " << cfi.getTrigger()
0152                                      << ", from EncodedEventId: " << cfi->eventId().bunchCrossing() << " "
0153                                      << cfi->eventId().event() << " bcr from MixCol " << cfi.bunch() << "\n";
0154       edm::LogInfo("RPDigiProducer") << " Hit: " << (*cfi) << "\n";
0155       count++;
0156     }
0157   }
0158
0159   MixCollection<PSimHit> allRPHits{cf.product(), std::pair(0, 0)};
0160
0161   if (verbosity_)
0162     edm::LogInfo("RPDigiProducer") << "Input MixCollection size = " << allRPHits.size() << "\n";
0163
0164   //Loop on PSimHit
0165   simhit_map simHitMap_;
0166   simHitMap_.clear();
0167
0168   MixCollection<PSimHit>::iterator isim;
0169   for (isim = allRPHits.begin(); isim != allRPHits.end(); ++isim) {
0170     simHitMap_[(*isim).detUnitId()].push_back((*isim));
0171   }
0172
0173   // Step B: LOOP on hits in event
0174   std::vector<edm::DetSet<TotemRPDigi>> DigiVector;
0175   DigiVector.reserve(400);
0176   DigiVector.clear();
0177
0178   CTPPSRPAlignmentCorrectionsData const* alignments = nullptr;
0179   if (auto rec = iSetup.tryToGet<VeryForwardMisalignedGeometryRecord>()) {
0180     alignments = &iSetup.getData(alignmentToken);
0181   }
0182   auto const& geom = iSetup.getData(geomToken);
0183
0184   for (simhit_map_iterator it = simHitMap_.begin(); it != simHitMap_.end(); ++it) {
0185     edm::DetSet<TotemRPDigi> digi_collector(it->first);
0186
0187     if (theAlgoMap.find(it->first) == theAlgoMap.end()) {
0188       theAlgoMap[it->first] = std::make_unique<RPDetDigitizer>(conf_, *rndEngine_, it->first, alignments, geom);
0189     }
0190
0191     std::vector<int> input_links;
0192     simromanpot::DigiPrimaryMapType output_digi_links;
0193
0194     (theAlgoMap.find(it->first)->second)
0195         ->run(simHitMap_[it->first], input_links, digi_collector.data, output_digi_links);
0196
0197     if (!digi_collector.data.empty()) {
0198       DigiVector.push_back(convertRPStripDetSet(digi_collector));
0199     }
0200   }
0201
0202   // Step C: create empty output collection
0203   std::unique_ptr<edm::DetSetVector<TotemRPDigi>> digi_output(new edm::DetSetVector<TotemRPDigi>(DigiVector));
0204
0205   if (verbosity_) {
0206     edm::LogInfo("RPDigiProducer") << "digi_output->size()=" << digi_output->size() << "\n";
0207   }
0208   // Step D: write output to file
0209   iEvent.put(std::move(digi_output));
0210 }
0211
0212 // ------------ method called once each job just before starting event loop  ------------
0213 void RPDigiProducer::beginRun(const edm::Run& beginrun, const edm::EventSetup& es) {
0214   // get analysis mask to mask channels
0215   if (simulateDeadChannels) {
0216     //set analysisMask in deadChannelsManager
0217     deadChannelsManager = DeadChannelsManager(&es.getData(tokenAnalysisMask));
0218   }
0219 }
0220
0221 edm::DetSet<TotemRPDigi> RPDigiProducer::convertRPStripDetSet(const edm::DetSet<TotemRPDigi>& rpstrip_detset) {
0222   edm::DetSet<TotemRPDigi> rpdigi_detset(rpstrip_detset.detId());
0223   rpdigi_detset.reserve(rpstrip_detset.size());
0224
0225   for (std::vector<TotemRPDigi>::const_iterator stripIterator = rpstrip_detset.data.begin();
0226        stripIterator < rpstrip_detset.data.end();
0227        ++stripIterator) {
0228     rpdigi_detset.push_back(TotemRPDigi(stripIterator->stripNumber()));
0229   }
0230
0231   return rpdigi_detset;
0232 }
0233
0234 void RPDigiProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0235   //RPSiDetDigitizer
0236   //all distances in [mm]
0237   edm::ParameterSetDescription desc;
0238   desc.add<bool>("RPLandauFluctuations", true);
0239   desc.add<bool>("RPDisplacementOn", false);
0240   desc.add<int>("RPVerbosity", 0);
0241   desc.add<double>("RPVFATThreshold", 9000.0);
0242   desc.add<double>("RPTopEdgePosition", 1.5);
0243   desc.add<double>("RPActiveEdgeSmearing", 0.013);
0244   desc.add<double>("RPEquivalentNoiseCharge300um", 1000.0);
0245   desc.add<int>("RPVFATTriggerMode", 2);
0246   desc.add<std::vector<double>>("RPInterStripSmearing",
0247                                 {
0248                                     0.011,
0249                                 });
0250   desc.add<double>("RPSharingSigmas", 5.0);  //how many sigmas taken into account for the edges and inter strips
0251   desc.add<double>("RPGeVPerElectron", 3.61e-09);
0252   desc.add<double>("RPActiveEdgePosition", 0.034);  //from the physical edge
0253   desc.add<bool>("RPDeadStripSimulationOn", false);
0254   desc.add<std::vector<std::string>>("ROUList",
0255                                      {
0256                                          "TotemHitsRP",
0257                                      });
0258   desc.add<bool>("RPNoNoise", false);
0259   desc.add<bool>("RPDigiSimHitRelationsPresistence", false);  //save links betweend digi, clusters and OSCAR/Geant4 hits
0260   desc.add<std::string>("mixLabel", "mix");
0261   desc.add<int>("RPChargeDivisionsPerThickness", 5);
0262   desc.add<double>("RPDeltaProductionCut", 0.120425);  //[MeV]
0263   desc.add<double>("RPBottomEdgePosition", 1.5);
0264   desc.add<double>("RPBottomEdgeSmearing", 0.011);
0265   desc.add<double>("RPTopEdgeSmearing", 0.011);
0266   desc.add<std::string>("InputCollection", "g4SimHitsTotemHitsRP");
0267   desc.add<double>("RPInterStripCoupling",
0268                    1.0);  //fraction of charge going to the strip, the missing part is taken by its neighbours
0269   desc.add<double>("RPDeadStripProbability", 0.001);
0270   desc.add<int>("RPChargeDivisionsPerStrip", 15);
0271   descriptions.add("RPSiDetDigitizer", desc);
0272 }
0273
0274 DEFINE_FWK_MODULE(RPDigiProducer);