Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:30:59

0001 // File: SiStripDigitizerAlgorithm.cc
0002 // Description:  Class for digitization.
0003 
0004 // Modified 15/May/2013 mark.grimes@bristol.ac.uk - Modified so that the digi-sim link has the correct
0005 // index for the sim hits stored. It was previously always set to zero (I won't mention that it was
0006 // me who originally wrote that).
0007 
0008 // system include files
0009 #include <memory>
0010 
0011 #include "SimTracker/Common/interface/SimHitSelectorFromDB.h"
0012 
0013 #include "SiStripDigitizer.h"
0014 #include "SiStripDigitizerAlgorithm.h"
0015 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0016 
0017 #include "DataFormats/Common/interface/DetSetVector.h"
0018 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
0019 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
0020 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0021 
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0030 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0031 
0032 #include "FWCore/Framework/interface/ConsumesCollector.h"
0033 //needed for the geometry:
0034 #include "FWCore/Framework/interface/EventSetup.h"
0035 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0036 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0037 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
0038 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0039 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0040 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0041 //needed for the magnetic field:
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 
0044 //Data Base infromations
0045 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
0046 
0047 //Random Number
0048 #include "FWCore/ServiceRegistry/interface/Service.h"
0049 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0050 #include "FWCore/Utilities/interface/Exception.h"
0051 #include "CLHEP/Random/RandFlat.h"
0052 
0053 SiStripDigitizer::SiStripDigitizer(const edm::ParameterSet& conf,
0054                                    edm::ProducesCollector producesCollector,
0055                                    edm::ConsumesCollector& iC)
0056     : hitsProducer(conf.getParameter<std::string>("hitsProducer")),
0057       trackerContainers(conf.getParameter<std::vector<std::string>>("ROUList")),
0058       ZSDigi(conf.getParameter<edm::ParameterSet>("DigiModeList").getParameter<std::string>("ZSDigi")),
0059       SCDigi(conf.getParameter<edm::ParameterSet>("DigiModeList").getParameter<std::string>("SCDigi")),
0060       VRDigi(conf.getParameter<edm::ParameterSet>("DigiModeList").getParameter<std::string>("VRDigi")),
0061       PRDigi(conf.getParameter<edm::ParameterSet>("DigiModeList").getParameter<std::string>("PRDigi")),
0062       useConfFromDB(conf.getParameter<bool>("TrackerConfigurationFromDB")),
0063       zeroSuppression(conf.getParameter<bool>("ZeroSuppression")),
0064       makeDigiSimLinks_(conf.getUntrackedParameter<bool>("makeDigiSimLinks", false)),
0065       includeAPVSimulation_(conf.getParameter<bool>("includeAPVSimulation")),
0066       fracOfEventsToSimAPV_(conf.getParameter<double>("fracOfEventsToSimAPV")),
0067       tTopoToken_(iC.esConsumes()),
0068       pDDToken_(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("GeometryType")))),
0069       pSetupToken_(iC.esConsumes()),
0070       gainToken_(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("Gain")))),
0071       noiseToken_(iC.esConsumes()),
0072       thresholdToken_(iC.esConsumes()),
0073       pedestalToken_(iC.esConsumes()),
0074       deadChannelToken_(iC.esConsumes()) {
0075   if (useConfFromDB) {
0076     detCablingToken_ = iC.esConsumes();
0077   }
0078   if (includeAPVSimulation_) {
0079     apvSimulationParametersToken_ = iC.esConsumes();
0080   }
0081 
0082   const std::string alias("simSiStripDigis");
0083 
0084   producesCollector.produces<edm::DetSetVector<SiStripDigi>>(ZSDigi).setBranchAlias(ZSDigi);
0085   producesCollector.produces<edm::DetSetVector<SiStripRawDigi>>(SCDigi).setBranchAlias(alias + SCDigi);
0086   producesCollector.produces<edm::DetSetVector<SiStripRawDigi>>(VRDigi).setBranchAlias(alias + VRDigi);
0087   producesCollector.produces<edm::DetSetVector<SiStripRawDigi>>("StripAmplitudes")
0088       .setBranchAlias(alias + "StripAmplitudes");
0089   producesCollector.produces<edm::DetSetVector<SiStripRawDigi>>("StripAmplitudesPostAPV")
0090       .setBranchAlias(alias + "StripAmplitudesPostAPV");
0091   producesCollector.produces<edm::DetSetVector<SiStripRawDigi>>("StripAPVBaselines")
0092       .setBranchAlias(alias + "StripAPVBaselines");
0093   producesCollector.produces<edm::DetSetVector<SiStripRawDigi>>(PRDigi).setBranchAlias(alias + PRDigi);
0094   producesCollector.produces<edm::DetSetVector<StripDigiSimLink>>().setBranchAlias(alias + "siStripDigiSimLink");
0095   producesCollector.produces<bool>("SimulatedAPVDynamicGain").setBranchAlias(alias + "SimulatedAPVDynamicGain");
0096   producesCollector.produces<std::vector<std::pair<int, std::bitset<6>>>>("AffectedAPVList")
0097       .setBranchAlias(alias + "AffectedAPV");
0098   for (auto const& trackerContainer : trackerContainers) {
0099     edm::InputTag tag(hitsProducer, trackerContainer);
0100     iC.consumes<std::vector<PSimHit>>(edm::InputTag(hitsProducer, trackerContainer));
0101   }
0102   edm::Service<edm::RandomNumberGenerator> rng;
0103   if (!rng.isAvailable()) {
0104     throw cms::Exception("Configuration")
0105         << "SiStripDigitizer requires the RandomNumberGeneratorService\n"
0106            "which is not present in the configuration file.  You must add the service\n"
0107            "in the configuration file or remove the modules that require it.";
0108   }
0109   theDigiAlgo = std::make_unique<SiStripDigitizerAlgorithm>(conf, iC);
0110 }
0111 
0112 // Virtual destructor needed.
0113 SiStripDigitizer::~SiStripDigitizer() {}
0114 
0115 void SiStripDigitizer::accumulateStripHits(edm::Handle<std::vector<PSimHit>> hSimHits,
0116                                            const TrackerTopology* tTopo,
0117                                            size_t globalSimHitIndex,
0118                                            const unsigned int tofBin) {
0119   // globalSimHitIndex is the index the sim hit will have when it is put in a collection
0120   // of sim hits for all crossings. This is only used when creating digi-sim links if
0121   // configured to do so.
0122 
0123   if (hSimHits.isValid()) {
0124     std::set<unsigned int> detIds;
0125     std::vector<PSimHit> const& simHits = *hSimHits.product();
0126     for (std::vector<PSimHit>::const_iterator it = simHits.begin(), itEnd = simHits.end(); it != itEnd;
0127          ++it, ++globalSimHitIndex) {
0128       unsigned int detId = (*it).detUnitId();
0129       if (detIds.insert(detId).second) {
0130         // The insert succeeded, so this detector element has not yet been processed.
0131         assert(detectorUnits[detId]);
0132         if (detectorUnits[detId]->type().isTrackerStrip()) {  // this test can be removed and replaced by stripdet!=0
0133           auto stripdet = detectorUnits[detId];
0134           //access to magnetic field in global coordinates
0135           GlobalVector bfield = pSetup->inTesla(stripdet->surface().position());
0136           LogDebug("Digitizer ") << "B-field(T) at " << stripdet->surface().position()
0137                                  << "(cm): " << pSetup->inTesla(stripdet->surface().position());
0138           theDigiAlgo->accumulateSimHits(it, itEnd, globalSimHitIndex, tofBin, stripdet, bfield, tTopo, randomEngine_);
0139         }
0140       }
0141     }  // end of loop over sim hits
0142   }
0143 }
0144 
0145 // Functions that gets called by framework every event
0146 void SiStripDigitizer::accumulate(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0147   //Retrieve tracker topology from geometry
0148   const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0149 
0150   // Step A: Get Inputs
0151   for (auto const& trackerContainer : trackerContainers) {
0152     edm::Handle<std::vector<PSimHit>> simHits;
0153     edm::InputTag tag(hitsProducer, trackerContainer);
0154     unsigned int tofBin = StripDigiSimLink::LowTof;
0155     if (trackerContainer.find(std::string("HighTof")) != std::string::npos)
0156       tofBin = StripDigiSimLink::HighTof;
0157 
0158     iEvent.getByLabel(tag, simHits);
0159     accumulateStripHits(simHits, tTopo, crossingSimHitIndexOffset_[tag.encode()], tofBin);
0160     // Now that the hits have been processed, I'll add the amount of hits in this crossing on to
0161     // the global counter. Next time accumulateStripHits() is called it will count the sim hits
0162     // as though they were on the end of this collection.
0163     // Note that this is only used for creating digi-sim links (if configured to do so).
0164     if (simHits.isValid())
0165       crossingSimHitIndexOffset_[tag.encode()] += simHits->size();
0166   }
0167 }
0168 
0169 void SiStripDigitizer::accumulate(PileUpEventPrincipal const& iEvent,
0170                                   edm::EventSetup const& iSetup,
0171                                   edm::StreamID const& streamID) {
0172   const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0173 
0174   //Re-compute luminosity for accumulation for HIP effects
0175   theDigiAlgo->calculateInstlumiScale(PileupInfo_.get());
0176 
0177   // Step A: Get Inputs
0178   for (auto const& trackerContainer : trackerContainers) {
0179     edm::Handle<std::vector<PSimHit>> simHits;
0180     edm::InputTag tag(hitsProducer, trackerContainer);
0181     unsigned int tofBin = StripDigiSimLink::LowTof;
0182     if (trackerContainer.find(std::string("HighTof")) != std::string::npos)
0183       tofBin = StripDigiSimLink::HighTof;
0184 
0185     iEvent.getByLabel(tag, simHits);
0186     accumulateStripHits(simHits, tTopo, crossingSimHitIndexOffset_[tag.encode()], tofBin);
0187     // Now that the hits have been processed, I'll add the amount of hits in this crossing on to
0188     // the global counter. Next time accumulateStripHits() is called it will count the sim hits
0189     // as though they were on the end of this collection.
0190     // Note that this is only used for creating digi-sim links (if configured to do so).
0191     if (simHits.isValid())
0192       crossingSimHitIndexOffset_[tag.encode()] += simHits->size();
0193   }
0194 }
0195 
0196 void SiStripDigitizer::initializeEvent(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0197   // Make sure that the first crossing processed starts indexing the sim hits from zero.
0198   // This variable is used so that the sim hits from all crossing frames have sequential
0199   // indices used to create the digi-sim link (if configured to do so) rather than starting
0200   // from zero for each crossing.
0201   crossingSimHitIndexOffset_.clear();
0202   theAffectedAPVvector.clear();
0203   // Step A: Get Inputs
0204 
0205   if (useConfFromDB) {
0206     iSetup.getData(detCablingToken_).addConnected(theDetIdList);
0207   }
0208 
0209   // Cache random number engine
0210   edm::Service<edm::RandomNumberGenerator> rng;
0211   randomEngine_ = &rng->getEngine(iEvent.streamID());
0212 
0213   theDigiAlgo->initializeEvent(iSetup);
0214 
0215   pDD = &iSetup.getData(pDDToken_);
0216   pSetup = &iSetup.getData(pSetupToken_);
0217 
0218   // We only need to clear and (re)fill detectorUnits when the geometry type IOV changes.
0219   auto ddCache = iSetup.get<TrackerDigiGeometryRecord>().cacheIdentifier();
0220   auto deadChannelCache = iSetup.get<SiStripBadChannelRcd>().cacheIdentifier();
0221   if (ddCache != ddCacheID_ or deadChannelCache != deadChannelCacheID_) {
0222     ddCacheID_ = ddCache;
0223     deadChannelCacheID_ = deadChannelCache;
0224     detectorUnits.clear();
0225 
0226     auto const& deadChannel = iSetup.getData(deadChannelToken_);
0227     for (const auto& iu : pDD->detUnits()) {
0228       unsigned int detId = iu->geographicalId().rawId();
0229       if (iu->type().isTrackerStrip()) {
0230         auto stripdet = dynamic_cast<StripGeomDetUnit const*>(iu);
0231         assert(stripdet != nullptr);
0232         detectorUnits.insert(std::make_pair(detId, stripdet));
0233         theDigiAlgo->initializeDetUnit(stripdet, deadChannel);
0234       }
0235     }
0236   }
0237 }
0238 
0239 void SiStripDigitizer::finalizeEvent(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0240   auto const& gain = iSetup.getData(gainToken_);
0241   auto const& noise = iSetup.getData(noiseToken_);
0242   auto const& threshold = iSetup.getData(thresholdToken_);
0243   auto const& pedestal = iSetup.getData(pedestalToken_);
0244   SiStripApvSimulationParameters const* apvSimulationParameters = nullptr;
0245 
0246   std::unique_ptr<bool> simulateAPVInThisEvent = std::make_unique<bool>(false);
0247   if (includeAPVSimulation_) {
0248     if (CLHEP::RandFlat::shoot(randomEngine_) < fracOfEventsToSimAPV_) {
0249       *simulateAPVInThisEvent = true;
0250       apvSimulationParameters = &iSetup.getData(apvSimulationParametersToken_);
0251     }
0252   }
0253   std::vector<edm::DetSet<SiStripDigi>> theDigiVector;
0254   std::vector<edm::DetSet<SiStripRawDigi>> theRawDigiVector;
0255   std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> theStripAmplitudeVector(new edm::DetSetVector<SiStripRawDigi>());
0256   std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> theStripAmplitudeVectorPostAPV(
0257       new edm::DetSetVector<SiStripRawDigi>());
0258   std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> theStripAPVBaselines(new edm::DetSetVector<SiStripRawDigi>());
0259   std::unique_ptr<edm::DetSetVector<StripDigiSimLink>> pOutputDigiSimLink(new edm::DetSetVector<StripDigiSimLink>);
0260 
0261   const TrackerTopology* tTopo = &iSetup.getData(tTopoToken_);
0262 
0263   // Step B: LOOP on StripGeomDetUnit
0264   theDigiVector.reserve(10000);
0265   theDigiVector.clear();
0266 
0267   for (const auto& iu : pDD->detUnits()) {
0268     if (useConfFromDB) {
0269       //apply the cable map _before_ digitization: consider only the detis that are connected
0270       if (theDetIdList.find(iu->geographicalId().rawId()) == theDetIdList.end())
0271         continue;
0272     }
0273     auto sgd = dynamic_cast<StripGeomDetUnit const*>(iu);
0274     if (sgd != nullptr) {
0275       edm::DetSet<SiStripDigi> collectorZS(iu->geographicalId().rawId());
0276       edm::DetSet<SiStripRawDigi> collectorRaw(iu->geographicalId().rawId());
0277       edm::DetSet<SiStripRawDigi> collectorStripAmplitudes(iu->geographicalId().rawId());
0278       edm::DetSet<SiStripRawDigi> collectorStripAmplitudesPostAPV(iu->geographicalId().rawId());
0279       edm::DetSet<SiStripRawDigi> collectorStripAPVBaselines(iu->geographicalId().rawId());
0280       edm::DetSet<StripDigiSimLink> collectorLink(iu->geographicalId().rawId());
0281 
0282       unsigned int detID = sgd->geographicalId().rawId();
0283       DetId detId(detID);
0284 
0285       theDigiAlgo->digitize(collectorZS,
0286                             collectorRaw,
0287                             collectorStripAmplitudes,
0288                             collectorStripAmplitudesPostAPV,
0289                             collectorStripAPVBaselines,
0290                             collectorLink,
0291                             sgd,
0292                             gain,
0293                             threshold,
0294                             noise,
0295                             pedestal,
0296                             *simulateAPVInThisEvent,
0297                             apvSimulationParameters,
0298                             theAffectedAPVvector,
0299                             randomEngine_,
0300                             tTopo);
0301 
0302       if (!collectorStripAmplitudes.data.empty())
0303         theStripAmplitudeVector->insert(collectorStripAmplitudes);
0304       if (!collectorStripAmplitudesPostAPV.data.empty())
0305         theStripAmplitudeVectorPostAPV->insert(collectorStripAmplitudesPostAPV);
0306       if (!collectorStripAPVBaselines.data.empty())
0307         theStripAPVBaselines->insert(collectorStripAPVBaselines);
0308 
0309       if (zeroSuppression) {
0310         if (!collectorZS.data.empty()) {
0311           theDigiVector.push_back(collectorZS);
0312           if (!collectorLink.data.empty())
0313             pOutputDigiSimLink->insert(collectorLink);
0314         }
0315       } else {
0316         if (!collectorRaw.data.empty()) {
0317           theRawDigiVector.push_back(collectorRaw);
0318           if (!collectorLink.data.empty())
0319             pOutputDigiSimLink->insert(collectorLink);
0320         }
0321       }
0322     }
0323   }
0324   if (zeroSuppression) {
0325     // Step C: create output collection
0326     std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> output_virginraw(new edm::DetSetVector<SiStripRawDigi>());
0327     std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> output_scopemode(new edm::DetSetVector<SiStripRawDigi>());
0328     std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> output_processedraw(new edm::DetSetVector<SiStripRawDigi>());
0329     std::unique_ptr<edm::DetSetVector<SiStripDigi>> output(new edm::DetSetVector<SiStripDigi>(theDigiVector));
0330     std::unique_ptr<std::vector<std::pair<int, std::bitset<6>>>> AffectedAPVList(
0331         new std::vector<std::pair<int, std::bitset<6>>>(theAffectedAPVvector));
0332 
0333     // Step D: write output to file
0334     iEvent.put(std::move(output), ZSDigi);
0335     iEvent.put(std::move(output_scopemode), SCDigi);
0336     iEvent.put(std::move(output_virginraw), VRDigi);
0337     iEvent.put(std::move(theStripAmplitudeVector), "StripAmplitudes");
0338     iEvent.put(std::move(theStripAmplitudeVectorPostAPV), "StripAmplitudesPostAPV");
0339     iEvent.put(std::move(theStripAPVBaselines), "StripAPVBaselines");
0340     iEvent.put(std::move(output_processedraw), PRDigi);
0341     iEvent.put(std::move(AffectedAPVList), "AffectedAPVList");
0342     iEvent.put(std::move(simulateAPVInThisEvent), "SimulatedAPVDynamicGain");
0343     if (makeDigiSimLinks_)
0344       iEvent.put(
0345           std::move(pOutputDigiSimLink));  // The previous EDProducer didn't name this collection so I won't either
0346   } else {
0347     // Step C: create output collection
0348     std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> output_virginraw(
0349         new edm::DetSetVector<SiStripRawDigi>(theRawDigiVector));
0350     std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> output_scopemode(new edm::DetSetVector<SiStripRawDigi>());
0351     std::unique_ptr<edm::DetSetVector<SiStripRawDigi>> output_processedraw(new edm::DetSetVector<SiStripRawDigi>());
0352     std::unique_ptr<edm::DetSetVector<SiStripDigi>> output(new edm::DetSetVector<SiStripDigi>());
0353 
0354     // Step D: write output to file
0355     iEvent.put(std::move(output), ZSDigi);
0356     iEvent.put(std::move(output_scopemode), SCDigi);
0357     iEvent.put(std::move(output_virginraw), VRDigi);
0358     iEvent.put(std::move(theStripAmplitudeVector), "StripAmplitudes");
0359     iEvent.put(std::move(theStripAmplitudeVectorPostAPV), "StripAmplitudesPostAPV");
0360     iEvent.put(std::move(theStripAPVBaselines), "StripAPVBaselines");
0361     iEvent.put(std::move(output_processedraw), PRDigi);
0362     iEvent.put(std::move(simulateAPVInThisEvent), "SimulatedAPVDynamicGain");
0363     if (makeDigiSimLinks_)
0364       iEvent.put(
0365           std::move(pOutputDigiSimLink));  // The previous EDProducer didn't name this collection so I won't either
0366   }
0367   randomEngine_ = nullptr;  // to prevent access outside event
0368 }