SiStripDigitizer/plugins/PreMixingSiStripWorker.cc

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "FWCore/Framework/interface/ProducesCollector.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0009 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0010 #include "CondFormats/SiStripObjects/interface/SiStripBadStrip.h"
0011 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
0012
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0015
0016 //Data Formats
0017 #include "DataFormats/Common/interface/DetSetVector.h"
0018 #include "DataFormats/Common/interface/DetSet.h"
0019 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
0020 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0021
0022 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0023 #include "CondFormats/SiStripObjects/interface/SiStripPedestals.h"
0024 #include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
0025 #include "CondFormats/SiStripObjects/interface/SiStripApvSimulationParameters.h"
0026 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
0027 #include "SimTracker/SiStripDigitizer/interface/SiTrivialDigitalConverter.h"
0028 #include "SimTracker/SiStripDigitizer/interface/SiGaussianTailNoiseAdder.h"
0029 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripFedZeroSuppression.h"
0030
0031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0032 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0034 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0035 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0036
0037 #include "CLHEP/Random/RandFlat.h"
0038
0039 #include "SimGeneral/PreMixingModule/interface/PreMixingWorker.h"
0040 #include "SimGeneral/PreMixingModule/interface/PreMixingWorkerFactory.h"
0041
0042 #include <map>
0043 #include <memory>
0044
0045 class PreMixingSiStripWorker : public PreMixingWorker {
0046 public:
0047   PreMixingSiStripWorker(const edm::ParameterSet& ps, edm::ProducesCollector, edm::ConsumesCollector&& iC);
0048   ~PreMixingSiStripWorker() override = default;
0049
0050   void initializeEvent(edm::Event const& e, edm::EventSetup const& c) override;
0051   void addSignals(edm::Event const& e, edm::EventSetup const& es) override;
0052   void addPileups(PileUpEventPrincipal const& pep, edm::EventSetup const& es) override;
0053   void put(edm::Event& e, edm::EventSetup const& iSetup, std::vector<PileupSummaryInfo> const& ps, int bs) override;
0054
0055 private:
0056   void DMinitializeDetUnit(StripGeomDetUnit const* det, const edm::EventSetup& iSetup);
0057
0058   // data specifiers
0059
0060   edm::InputTag SistripLabelSig_;        // name given to collection of SiStrip digis
0061   edm::InputTag SiStripPileInputTag_;    // InputTag for pileup strips
0062   std::string SiStripDigiCollectionDM_;  // secondary name to be given to new SiStrip digis
0063
0064   edm::InputTag SistripAPVLabelSig_;  // where to find vector of dead APVs
0065   edm::InputTag SiStripAPVPileInputTag_;
0066   std::string SistripAPVListDM_;  // output tag
0067
0068   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> const pDDToken_;
0069   edm::ESGetToken<SiStripBadStrip, SiStripBadChannelRcd> const deadChannelToken_;
0070   edm::ESGetToken<SiStripGain, SiStripGainSimRcd> const gainToken_;
0071   edm::ESGetToken<SiStripNoises, SiStripNoisesRcd> const noiseToken_;
0072   edm::ESGetToken<SiStripThreshold, SiStripThresholdRcd> const thresholdToken_;
0073   //edm::ESGetToken<SiStripPedestals, SiStripPedestalsRcd> const pedestalToken_;
0074   edm::ESGetToken<SiStripApvSimulationParameters, SiStripApvSimulationParametersRcd> apvSimulationParametersToken_;
0075   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0076
0077   //
0078
0079   typedef float Amplitude;
0080   typedef std::pair<uint16_t, Amplitude> RawDigi;  // Replacement for SiStripDigi with pulse height instead of integer ADC
0081   typedef std::vector<SiStripDigi> OneDetectorMap;  // maps by strip ID for later combination - can have duplicate strips
0082   typedef std::vector<RawDigi> OneDetectorRawMap;  // maps by strip ID for later combination - can have duplicate strips
0083   typedef std::map<uint32_t, OneDetectorMap> SiGlobalIndex;                      // map to all data for each detector ID
0084   typedef std::vector<std::pair<uint32_t, OneDetectorRawMap>> SiGlobalRawIndex;  // map to all data for each detector ID
0085
0086   typedef SiDigitalConverter::DigitalVecType DigitalVecType;
0087
0088   SiGlobalIndex SiHitStorage_;
0089   SiGlobalRawIndex SiRawDigis_;
0090
0091   // variables for temporary storage of mixed hits:
0092   typedef std::vector<std::pair<int, Amplitude>> SignalMapType;
0093   typedef std::map<uint32_t, SignalMapType> signalMaps;
0094
0095   SignalMapType const* getSignal(uint32_t detID) const {
0096     auto where = signals_.find(detID);
0097     if (where == signals_.end()) {
0098       return nullptr;
0099     }
0100     return &where->second;
0101   }
0102
0103   signalMaps signals_;
0104
0105   // to keep track of dead APVs from HIP interactions
0106   typedef std::multimap<uint32_t, std::bitset<6>> APVMap;
0107
0108   APVMap theAffectedAPVmap_;
0109
0110   // for noise adding:
0111
0112   bool SingleStripNoise;
0113   bool peakMode;
0114   double theThreshold;
0115   double theElectronPerADC;
0116   bool APVSaturationFromHIP_;
0117   int theFedAlgo;
0118
0119   std::unique_ptr<SiGaussianTailNoiseAdder> theSiNoiseAdder;
0120   std::unique_ptr<SiStripFedZeroSuppression> theSiZeroSuppress;
0121   std::unique_ptr<SiTrivialDigitalConverter> theSiDigitalConverter;
0122
0123   const TrackerGeometry* pDD;
0124
0125   // bad channels for each detector ID
0126   std::map<unsigned int, std::vector<bool>> allBadChannels;
0127   // channels killed by HIP interactions for each detector ID
0128   std::map<unsigned int, std::vector<bool>> allHIPChannels;
0129   // first and last channel wit signal for each detector ID
0130   std::map<unsigned int, size_t> firstChannelsWithSignal;
0131   std::map<unsigned int, size_t> lastChannelsWithSignal;
0132
0133   bool includeAPVSimulation_;
0134   const double fracOfEventsToSimAPV_;
0135   const double apv_maxResponse_;
0136   const double apv_rate_;
0137   const double apv_mVPerQ_;
0138   const double apv_fCPerElectron_;
0139
0140   //----------------------------
0141
0142   class StrictWeakOrdering {
0143   public:
0144     bool operator()(SiStripDigi i, SiStripDigi j) const { return i.strip() < j.strip(); }
0145   };
0146
0147   class StrictWeakRawOrdering {
0148   public:
0149     bool operator()(RawDigi i, RawDigi j) const { return i.first < j.first; }
0150   };
0151 };
0152
0153 PreMixingSiStripWorker::PreMixingSiStripWorker(const edm::ParameterSet& ps,
0154                                                edm::ProducesCollector producesCollector,
0155                                                edm::ConsumesCollector&& iC)
0156     : pDDToken_(iC.esConsumes(edm::ESInputTag("", ps.getParameter<std::string>("GeometryType")))),
0157       deadChannelToken_(iC.esConsumes()),
0158       gainToken_(iC.esConsumes(edm::ESInputTag("", ps.getParameter<std::string>("Gain")))),
0159       noiseToken_(iC.esConsumes()),
0160       thresholdToken_(iC.esConsumes()),
0161       SingleStripNoise(ps.getParameter<bool>("SingleStripNoise")),
0162       peakMode(ps.getParameter<bool>("APVpeakmode")),
0163       theThreshold(ps.getParameter<double>("NoiseSigmaThreshold")),
0164       theElectronPerADC(ps.getParameter<double>(peakMode ? "electronPerAdcPeak" : "electronPerAdcDec")),
0165       APVSaturationFromHIP_(ps.getParameter<bool>("APVSaturationFromHIP")),
0166       theFedAlgo(ps.getParameter<int>("FedAlgorithm_PM")),
0167       theSiZeroSuppress(new SiStripFedZeroSuppression(theFedAlgo, &iC)),
0168       theSiDigitalConverter(new SiTrivialDigitalConverter(theElectronPerADC, false)),  // no premixing
0169       includeAPVSimulation_(ps.getParameter<bool>("includeAPVSimulation")),
0170       fracOfEventsToSimAPV_(ps.getParameter<double>("fracOfEventsToSimAPV")),
0171       apv_maxResponse_(ps.getParameter<double>("apv_maxResponse")),
0172       apv_rate_(ps.getParameter<double>("apv_rate")),
0173       apv_mVPerQ_(ps.getParameter<double>("apv_mVPerQ")),
0174       apv_fCPerElectron_(ps.getParameter<double>("apvfCPerElectron"))
0175
0176 {
0177   // declare the products to produce
0178
0179   SistripLabelSig_ = ps.getParameter<edm::InputTag>("SistripLabelSig");
0180   SiStripPileInputTag_ = ps.getParameter<edm::InputTag>("SiStripPileInputTag");
0181
0182   SiStripDigiCollectionDM_ = ps.getParameter<std::string>("SiStripDigiCollectionDM");
0183   SistripAPVListDM_ = ps.getParameter<std::string>("SiStripAPVListDM");
0184
0185   producesCollector.produces<edm::DetSetVector<SiStripDigi>>(SiStripDigiCollectionDM_);
0186   producesCollector.produces<bool>(SiStripDigiCollectionDM_ + "SimulatedAPVDynamicGain");
0187
0188   if (APVSaturationFromHIP_) {
0189     SistripAPVLabelSig_ = ps.getParameter<edm::InputTag>("SistripAPVLabelSig");
0190     SiStripAPVPileInputTag_ = ps.getParameter<edm::InputTag>("SistripAPVPileInputTag");
0191     iC.consumes<std::vector<std::pair<int, std::bitset<6>>>>(SistripAPVLabelSig_);
0192   }
0193   iC.consumes<edm::DetSetVector<SiStripDigi>>(SistripLabelSig_);
0194
0195   if (includeAPVSimulation_) {
0196     tTopoToken_ = iC.esConsumes();
0197     apvSimulationParametersToken_ = iC.esConsumes();
0198   }
0199   // clear local storage for this event
0200   SiHitStorage_.clear();
0201
0202   edm::Service<edm::RandomNumberGenerator> rng;
0203   if (!rng.isAvailable()) {
0204     throw cms::Exception("Psiguration") << "SiStripDigitizer requires the RandomNumberGeneratorService\n"
0205                                            "which is not present in the psiguration file.  You must add the service\n"
0206                                            "in the configuration file or remove the modules that require it.";
0207   }
0208
0209   theSiNoiseAdder = std::make_unique<SiGaussianTailNoiseAdder>(theThreshold);
0210 }
0211
0212 void PreMixingSiStripWorker::initializeEvent(const edm::Event& e, edm::EventSetup const& iSetup) {
0213   // initialize individual detectors so we can copy real digitization code:
0214
0215   pDD = &iSetup.getData(pDDToken_);
0216
0217   for (auto iu = pDD->detUnits().begin(); iu != pDD->detUnits().end(); ++iu) {
0218     unsigned int detId = (*iu)->geographicalId().rawId();
0219     DetId idet = DetId(detId);
0220     unsigned int isub = idet.subdetId();
0221     if ((isub == StripSubdetector::TIB) || (isub == StripSubdetector::TID) || (isub == StripSubdetector::TOB) ||
0222         (isub == StripSubdetector::TEC)) {
0223       auto stripdet = dynamic_cast<StripGeomDetUnit const*>((*iu));
0224       assert(stripdet != nullptr);
0225       DMinitializeDetUnit(stripdet, iSetup);
0226     }
0227   }
0228 }
0229
0230 void PreMixingSiStripWorker::DMinitializeDetUnit(StripGeomDetUnit const* det, const edm::EventSetup& iSetup) {
0231   SiStripBadStrip const& deadChannel = iSetup.getData(deadChannelToken_);
0232
0233   unsigned int detId = det->geographicalId().rawId();
0234   int numStrips = (det->specificTopology()).nstrips();
0235
0236   SiStripBadStrip::Range detBadStripRange = deadChannel.getRange(detId);
0237   //storing the bad strip of the the module. the module is not removed but just signal put to 0
0238   std::vector<bool>& badChannels = allBadChannels[detId];
0239   std::vector<bool>& hipChannels = allHIPChannels[detId];
0240   badChannels.clear();
0241   badChannels.insert(badChannels.begin(), numStrips, false);
0242   hipChannels.clear();
0243   hipChannels.insert(hipChannels.begin(), numStrips, false);
0244
0245   for (SiStripBadStrip::ContainerIterator it = detBadStripRange.first; it != detBadStripRange.second; ++it) {
0246     SiStripBadStrip::data fs = deadChannel.decode(*it);
0247     for (int strip = fs.firstStrip; strip < fs.firstStrip + fs.range; ++strip)
0248       badChannels[strip] = true;
0249   }
0250   firstChannelsWithSignal[detId] = numStrips;
0251   lastChannelsWithSignal[detId] = 0;
0252 }
0253
0254 void PreMixingSiStripWorker::addSignals(const edm::Event& e, edm::EventSetup const& es) {
0255   // fill in maps of hits
0256
0257   edm::Handle<edm::DetSetVector<SiStripDigi>> input;
0258
0259   if (e.getByLabel(SistripLabelSig_, input)) {
0260     OneDetectorMap LocalMap;
0261
0262     //loop on all detsets (detectorIDs) inside the input collection
0263     edm::DetSetVector<SiStripDigi>::const_iterator DSViter = input->begin();
0264     for (; DSViter != input->end(); DSViter++) {
0265 #ifdef DEBUG
0266       LogDebug("PreMixingSiStripWorker") << "Processing DetID " << DSViter->id;
0267 #endif
0268
0269       LocalMap.clear();
0270       LocalMap.reserve((DSViter->data).size());
0271       LocalMap.insert(LocalMap.end(), (DSViter->data).begin(), (DSViter->data).end());
0272
0273       SiHitStorage_.insert(SiGlobalIndex::value_type(DSViter->id, LocalMap));
0274     }
0275   }
0276
0277   // keep here for future reference.  In current implementation, HIP killing is done once in PU file
0278   /*  if(APVSaturationFromHIP_) {
0279       edm::Handle<std::vector<std::pair<int,std::bitset<6>> > >  APVinput;
0280
0281       if( e.getByLabel(SistripAPVLabelSig_,APVinput) ) {
0282
0283       std::vector<std::pair<int,std::bitset<6>> >::const_iterator entry = APVinput->begin();
0284       for( ; entry != APVinput->end(); entry++) {
0285       theAffectedAPVmap_.insert(APVMap::value_type(entry->first, entry->second));
0286       }
0287       }
0288       } */
0289
0290 }  // end of addSiStripSignals
0291
0292 void PreMixingSiStripWorker::addPileups(PileUpEventPrincipal const& pep, edm::EventSetup const& es) {
0293   LogDebug("PreMixingSiStripWorker") << "\n===============> adding pileups from event  " << pep.principal().id()
0294                                      << " for bunchcrossing " << pep.bunchCrossing();
0295
0296   // fill in maps of hits; same code as addSignals, except now applied to the pileup events
0297
0298   edm::Handle<edm::DetSetVector<SiStripDigi>> inputHandle;
0299   pep.getByLabel(SiStripPileInputTag_, inputHandle);
0300
0301   if (inputHandle.isValid()) {
0302     const auto& input = *inputHandle;
0303
0304     OneDetectorMap LocalMap;
0305
0306     //loop on all detsets (detectorIDs) inside the input collection
0307     edm::DetSetVector<SiStripDigi>::const_iterator DSViter = input.begin();
0308     for (; DSViter != input.end(); DSViter++) {
0309 #ifdef DEBUG
0310       LogDebug("PreMixingSiStripWorker") << "Pileups: Processing DetID " << DSViter->id;
0311 #endif
0312
0313       // find correct local map (or new one) for this detector ID
0314
0315       SiGlobalIndex::const_iterator itest;
0316
0317       itest = SiHitStorage_.find(DSViter->id);
0318
0319       if (itest != SiHitStorage_.end()) {  // this detID already has hits, add to existing map
0320
0321         LocalMap = itest->second;
0322
0323         // fill in local map with extra channels
0324         LocalMap.insert(LocalMap.end(), (DSViter->data).begin(), (DSViter->data).end());
0325         std::stable_sort(LocalMap.begin(), LocalMap.end(), PreMixingSiStripWorker::StrictWeakOrdering());
0326         SiHitStorage_[DSViter->id] = LocalMap;
0327
0328       } else {  // fill local storage with this information, put in global collection
0329
0330         LocalMap.clear();
0331         LocalMap.reserve((DSViter->data).size());
0332         LocalMap.insert(LocalMap.end(), (DSViter->data).begin(), (DSViter->data).end());
0333
0334         SiHitStorage_.insert(SiGlobalIndex::value_type(DSViter->id, LocalMap));
0335       }
0336     }
0337
0338     if (APVSaturationFromHIP_) {
0339       edm::Handle<std::vector<std::pair<int, std::bitset<6>>>> inputAPVHandle;
0340       pep.getByLabel(SiStripAPVPileInputTag_, inputAPVHandle);
0341
0342       if (inputAPVHandle.isValid()) {
0343         const auto& APVinput = inputAPVHandle;
0344
0345         std::vector<std::pair<int, std::bitset<6>>>::const_iterator entry = APVinput->begin();
0346         for (; entry != APVinput->end(); entry++) {
0347           theAffectedAPVmap_.insert(APVMap::value_type(entry->first, entry->second));
0348         }
0349       }
0350     }
0351   }
0352 }
0353
0354 void PreMixingSiStripWorker::put(edm::Event& e,
0355                                  edm::EventSetup const& iSetup,
0356                                  std::vector<PileupSummaryInfo> const& ps,
0357                                  int bs) {
0358   // set up machinery to do proper noise adding:
0359   SiStripGain const& gain = iSetup.getData(gainToken_);
0360   SiStripNoises const& noise = iSetup.getData(noiseToken_);
0361   SiStripThreshold const& threshold = iSetup.getData(thresholdToken_);
0362   //SiStripPedestals const& pedestal = iSetup.getData(pedestalToken_);
0363   SiStripApvSimulationParameters const* apvSimulationParameters = nullptr;
0364   TrackerTopology const* tTopo = nullptr;
0365
0366   edm::Service<edm::RandomNumberGenerator> rng;
0367   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0368
0369   const bool simulateAPVInThisEvent = includeAPVSimulation_ && (CLHEP::RandFlat::shoot(engine) < fracOfEventsToSimAPV_);
0370   float nTruePU = 0.;  // = ps.getTrueNumInteractions();
0371   if (simulateAPVInThisEvent) {
0372     tTopo = &iSetup.getData(tTopoToken_);
0373     apvSimulationParameters = &iSetup.getData(apvSimulationParametersToken_);
0374     const auto it = std::find_if(
0375         std::begin(ps), std::end(ps), [](const PileupSummaryInfo& bxps) { return bxps.getBunchCrossing() == 0; });
0376     if (it != std::begin(ps)) {
0377       nTruePU = it->getTrueNumInteractions();
0378     } else {
0379       edm::LogWarning("PreMixingSiStripWorker") << "Could not find PileupSummaryInfo for current bunch crossing";
0380       nTruePU = std::begin(ps)->getTrueNumInteractions();
0381     }
0382   }
0383
0384   std::map<int, std::bitset<6>> DeadAPVList;
0385
0386   // First, have to convert all ADC counts to raw pulse heights so that values can be added properly
0387   // In PreMixing, pulse heights are saved with ADC = sqrt(9.0*PulseHeight) - have to undo.
0388
0389   // This is done here because it's the only place we have access to EventSetup
0390   // Simultaneously, merge lists of hit channels in each DetId.
0391   // Signal Digis are in the list first, have to merge lists of hit strips on the fly,
0392   // add signals on duplicates later
0393
0394   // Now, loop over hits and add them to the map in the proper sorted order
0395   // Note: We are assuming that the hits from the Signal events have been created in
0396   // "PreMix" mode, rather than in the standard ADC conversion routines.  If not, this
0397   // doesn't work at all.
0398
0399   // At the moment, both Signal and Reconstituted PU hits have the same compression algorithm.
0400   // If this were different, and one needed gains, the conversion back to pulse height can only
0401   // be done in this routine.  So, yes, there is an extra loop over hits here in the current code,
0402   // because, in principle, one could convert to pulse height during the read/store phase.
0403
0404   for (SiGlobalIndex::const_iterator IDet = SiHitStorage_.begin(); IDet != SiHitStorage_.end(); IDet++) {
0405     uint32_t detID = IDet->first;
0406
0407     auto const& LocalMap = IDet->second;
0408
0409     //loop over hit strips for this DetId, do conversion to pulse height, store.
0410     SiRawDigis_.emplace_back(detID, OneDetectorRawMap());
0411     auto& localRawMap = SiRawDigis_.back().second;
0412     localRawMap.reserve(LocalMap.size());
0413     for (auto const& iLocal : LocalMap) {
0414       uint16_t currentStrip = iLocal.strip();
0415       float signal = float(iLocal.adc());
0416       if (iLocal.adc() == 1022)
0417         signal = 1500.f;  // average values for overflows
0418       if (iLocal.adc() == 1023)
0419         signal = 3000.f;
0420
0421       //convert signals back to raw counts
0422       constexpr float inv9 = 1.f / 9.0f;
0423       float ReSignal = signal * signal * inv9;  // The PreMixing conversion is adc = sqrt(9.0*pulseHeight)
0424
0425       // RawDigi NewRawDigi = std::make_pair(currentStrip, ReSignal);
0426
0427       localRawMap.emplace_back(currentStrip, ReSignal);
0428     }
0429   }
0430
0431   // If we are killing APVs, merge list of dead ones before we digitize
0432
0433   int NumberOfBxBetweenHIPandEvent = 1e3;
0434
0435   if (APVSaturationFromHIP_) {
0436     // calculate affected BX parameter
0437
0438     bool HasAtleastOneAffectedAPV = false;
0439     while (!HasAtleastOneAffectedAPV) {
0440       for (int bx = floor(300.0 / 25.0); bx > 0; bx--) {  //Reminder: make these numbers not hard coded!!
0441         float temp = CLHEP::RandFlat::shoot(engine) < 0.5 ? 1 : 0;
0442         if (temp == 1 && bx < NumberOfBxBetweenHIPandEvent) {
0443           NumberOfBxBetweenHIPandEvent = bx;
0444           HasAtleastOneAffectedAPV = true;
0445         }
0446       }
0447     }
0448
0449     APVMap::const_iterator iAPVchk;
0450     uint32_t formerID = 0;
0451     uint32_t currentID;
0452     std::bitset<6> NewAPVBits;
0453
0454     for (APVMap::const_iterator iAPV = theAffectedAPVmap_.begin(); iAPV != theAffectedAPVmap_.end(); ++iAPV) {
0455       currentID = iAPV->first;
0456
0457       if (currentID == formerID) {  // we have to OR these
0458         for (int ibit = 0; ibit < 6; ++ibit) {
0459           NewAPVBits[ibit] = NewAPVBits[ibit] || (iAPV->second)[ibit];
0460         }
0461       } else {
0462         DeadAPVList[currentID] = NewAPVBits;
0463         //save pointers for next iteration
0464         formerID = currentID;
0465         NewAPVBits = iAPV->second;
0466       }
0467
0468       iAPVchk = iAPV;
0469       if ((++iAPVchk) == theAffectedAPVmap_.end()) {  //make sure not to lose the last one
0470         DeadAPVList[currentID] = NewAPVBits;
0471       }
0472     }
0473   }
0474   //
0475
0476   //  Ok, done with merging raw signals and APVs - now add signals on duplicate strips
0477
0478   // collection of Digis to put in the event
0479   std::vector<edm::DetSet<SiStripDigi>> vSiStripDigi;
0480
0481   // loop through our collection of detectors, merging hits and making a new list of "signal" digis
0482
0483   // clear some temporary storage for later digitization:
0484
0485   signals_.clear();
0486
0487   // big loop over Detector IDs:
0488   for (auto const& IDet : SiRawDigis_) {
0489     uint32_t detID = IDet.first;
0490
0491     // create merged map:
0492     auto& lSignals = signals_[detID];
0493
0494     auto const& LocalMap = IDet.second;
0495
0496     // guess max occupancy...
0497     lSignals.reserve(std::min(LocalMap.size(), 512ul));
0498
0499     //counter variables
0500     int formerStrip = -1;
0501     int currentStrip;
0502     float ADCSum = 0;
0503
0504     //loop over hit strips for this DetId, add duplicates
0505
0506     OneDetectorRawMap::const_iterator iLocalchk;
0507     OneDetectorRawMap::const_iterator iLocal = LocalMap.begin();
0508     for (; iLocal != LocalMap.end(); ++iLocal) {
0509       currentStrip = iLocal->first;  // strip is first element
0510
0511       if (currentStrip == formerStrip) {  // we have to add these digis together
0512
0513         ADCSum += iLocal->second;  // raw pulse height is second element.
0514       } else {
0515         if (formerStrip != -1) {
0516           lSignals.emplace_back(formerStrip, ADCSum);
0517         }
0518         // save pointers for next iteration
0519         formerStrip = currentStrip;
0520         ADCSum = iLocal->second;  // lone ADC
0521       }
0522
0523       iLocalchk = iLocal;
0524       if ((++iLocalchk) == LocalMap.end()) {  //make sure not to lose the last one
0525         lSignals.emplace_back(formerStrip, ADCSum);
0526       }
0527     }
0528   }
0529
0530   //Now, do noise, zero suppression, take into account bad channels, etc.
0531   // This section stolen from SiStripDigitizerAlgorithm
0532   // must loop over all detIds in the tracker to get all of the noise added properly.
0533   for (const auto& iu : pDD->detUnits()) {
0534     const StripGeomDetUnit* sgd = dynamic_cast<const StripGeomDetUnit*>(iu);
0535     if (sgd != nullptr) {
0536       uint32_t detID = sgd->geographicalId().rawId();
0537
0538       edm::DetSet<SiStripDigi> SSD(detID);  // Make empty collection with this detector ID
0539
0540       int numStrips = (sgd->specificTopology()).nstrips();
0541
0542       // see if there is some signal on this detector
0543
0544       auto const* lSignal = getSignal(detID);
0545
0546       std::vector<float> detAmpl(numStrips, 0.);
0547       if (lSignal) {
0548         for (const auto& amp : *lSignal) {
0549           detAmpl[amp.first] = amp.second;
0550         }
0551       }
0552
0553       //removing signal from the dead (and HIP effected) strips
0554       std::vector<bool>& badChannels = allBadChannels[detID];
0555
0556       for (int strip = 0; strip < numStrips; ++strip) {
0557         if (badChannels[strip])
0558           detAmpl[strip] = 0.;
0559       }
0560
0561       if (simulateAPVInThisEvent) {
0562         // Get index in apv baseline distributions corresponding to z of detSet and PU
0563         const StripTopology* topol = dynamic_cast<const StripTopology*>(&(sgd->specificTopology()));
0564         LocalPoint localPos = topol->localPosition(0);
0565         GlobalPoint globalPos = sgd->surface().toGlobal(Local3DPoint(localPos.x(), localPos.y(), localPos.z()));
0566         float detSet_z = fabs(globalPos.z());
0567         float detSet_r = globalPos.perp();
0568
0569         const uint32_t SubDet = DetId(detID).subdetId();
0570         // Simulate APV response for each strip
0571         for (int strip = 0; strip < numStrips; ++strip) {
0572           if (detAmpl[strip] > 0) {
0573             // Convert charge from electrons to fC
0574             double stripCharge = detAmpl[strip] * apv_fCPerElectron_;
0575
0576             // Get APV baseline
0577             double baselineV = 0;
0578             if (SubDet == SiStripSubdetector::TIB) {
0579               baselineV = apvSimulationParameters->sampleTIB(tTopo->tibLayer(detID), detSet_z, nTruePU, engine);
0580             } else if (SubDet == SiStripSubdetector::TOB) {
0581               baselineV = apvSimulationParameters->sampleTOB(tTopo->tobLayer(detID), detSet_z, nTruePU, engine);
0582             } else if (SubDet == SiStripSubdetector::TID) {
0583               baselineV = apvSimulationParameters->sampleTID(tTopo->tidWheel(detID), detSet_r, nTruePU, engine);
0584             } else if (SubDet == SiStripSubdetector::TEC) {
0585               baselineV = apvSimulationParameters->sampleTEC(tTopo->tecWheel(detID), detSet_r, nTruePU, engine);
0586             }
0587             // Fitted parameters from G Hall/M Raymond
0588             double maxResponse = apv_maxResponse_;
0589             double rate = apv_rate_;
0590
0591             double outputChargeInADC = 0;
0592             if (baselineV < apv_maxResponse_) {
0593               // Convert V0 into baseline charge
0594               double baselineQ = -1.0 * rate * log(2 * maxResponse / (baselineV + maxResponse) - 1);
0595
0596               // Add charge deposited in this BX
0597               double newStripCharge = baselineQ + stripCharge;
0598
0599               // Apply APV response
0600               double signalV = 2 * maxResponse / (1 + exp(-1.0 * newStripCharge / rate)) - maxResponse;
0601               double gain = signalV - baselineV;
0602
0603               // Convert gain (mV) to charge (assuming linear region of APV) and then to electrons
0604               double outputCharge = gain / apv_mVPerQ_;
0605               outputChargeInADC = outputCharge / apv_fCPerElectron_;
0606             }
0607
0608             // Output charge back to original container
0609             detAmpl[strip] = outputChargeInADC;
0610           }
0611         }
0612       }
0613
0614       if (APVSaturationFromHIP_) {
0615         std::bitset<6>& bs = DeadAPVList[detID];
0616
0617         if (bs.any()) {
0618           // Here below is the scaling function which describes the evolution of the baseline (i.e. how the charge is suppressed).
0619           // This must be replaced as soon as we have a proper modeling of the baseline evolution from VR runs
0620           float Shift =
0621               1 - NumberOfBxBetweenHIPandEvent / floor(300.0 / 25.0);  //Reminder: make these numbers not hardcoded!!
0622           float randomX = CLHEP::RandFlat::shoot(engine);
0623           float scalingValue = (randomX - Shift) * 10.0 / 7.0 - 3.0 / 7.0;
0624
0625           for (int strip = 0; strip < numStrips; ++strip) {
0626             if (!badChannels[strip] && bs[strip / 128] == 1) {
0627               detAmpl[strip] *= scalingValue > 0 ? scalingValue : 0.0;
0628             }
0629           }
0630         }
0631       }
0632
0633       SiStripNoises::Range detNoiseRange = noise.getRange(detID);
0634       SiStripApvGain::Range detGainRange = gain.getRange(detID);
0635
0636       // Gain conversion is already done during signal adding
0637       //convert our signals back to raw counts so that we can add noise properly:
0638
0639       /*
0640         if(theSignal) {
0641         for(unsigned int iv = 0; iv!=detAmpl.size(); iv++) {
0642         float signal = detAmpl[iv];
0643         if(signal > 0) {
0644         float gainValue = gain.getStripGain(iv, detGainRange);
0645         signal *= theElectronPerADC/gainValue;
0646         detAmpl[iv] = signal;
0647         }
0648         }
0649         }
0650       */
0651
0652       //SiStripPedestals::Range detPedestalRange = pedestal.getRange(detID);
0653
0654       // -----------------------------------------------------------
0655
0656       size_t firstChannelWithSignal = 0;
0657       size_t lastChannelWithSignal = numStrips;
0658
0659       if (SingleStripNoise) {
0660         std::vector<float> noiseRMSv;
0661         noiseRMSv.clear();
0662         noiseRMSv.insert(noiseRMSv.begin(), numStrips, 0.);
0663         for (int strip = 0; strip < numStrips; ++strip) {
0664           if (!badChannels[strip]) {
0665             float gainValue = gain.getStripGain(strip, detGainRange);
0666             noiseRMSv[strip] = (noise.getNoise(strip, detNoiseRange)) * theElectronPerADC / gainValue;
0667           }
0668         }
0669         theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv, engine);
0670       } else {
0671         int RefStrip = int(numStrips / 2.);
0672         while (RefStrip < numStrips &&
0673                badChannels[RefStrip]) {  //if the refstrip is bad, I move up to when I don't find it
0674           RefStrip++;
0675         }
0676         if (RefStrip < numStrips) {
0677           float RefgainValue = gain.getStripGain(RefStrip, detGainRange);
0678           float RefnoiseRMS = noise.getNoise(RefStrip, detNoiseRange) * theElectronPerADC / RefgainValue;
0679
0680           theSiNoiseAdder->addNoise(
0681               detAmpl, firstChannelWithSignal, lastChannelWithSignal, numStrips, RefnoiseRMS, engine);
0682         }
0683       }
0684
0685       auto const& data = theSiDigitalConverter->convert(detAmpl, &gain, detID);
0686       // if suppression is trival just copy data
0687       if (5 == theFedAlgo)
0688         SSD.data = data;
0689       else
0690         theSiZeroSuppress->suppress(data, SSD.data, detID, noise, threshold);
0691
0692       LogDebug("PreMixingSiStripWorker") << detID << " " << data.size() << ' ' << SSD.data.size();
0693
0694       // stick this into the global vector of detector info
0695       vSiStripDigi.push_back(SSD);
0696
0697     }  // end of loop over one detector
0698
0699   }  // end of big loop over all detector IDs
0700
0701   // put the collection of digis in the event
0702   edm::LogInfo("PreMixingSiStripWorker") << "total # Merged strips: " << vSiStripDigi.size();
0703
0704   // make new digi collection
0705
0706   std::unique_ptr<edm::DetSetVector<SiStripDigi>> MySiStripDigis(new edm::DetSetVector<SiStripDigi>(vSiStripDigi));
0707
0708   // put collection
0709
0710   e.put(std::move(MySiStripDigis), SiStripDigiCollectionDM_);
0711   e.put(std::make_unique<bool>(simulateAPVInThisEvent), SiStripDigiCollectionDM_ + "SimulatedAPVDynamicGain");
0712
0713   // clear local storage for this event
0714   SiHitStorage_.clear();
0715   SiRawDigis_.clear();
0716   signals_.clear();
0717 }
0718
0719 DEFINE_PREMIXING_WORKER(PreMixingSiStripWorker);