Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-12 02:41:38

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::map<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::map<int, Amplitude> SignalMapType;
0093   typedef std::map<uint32_t, SignalMapType> signalMaps;
0094 
0095   const SignalMapType* 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   DeadAPVList.clear();
0386 
0387   // First, have to convert all ADC counts to raw pulse heights so that values can be added properly
0388   // In PreMixing, pulse heights are saved with ADC = sqrt(9.0*PulseHeight) - have to undo.
0389 
0390   // This is done here because it's the only place we have access to EventSetup
0391   // Simultaneously, merge lists of hit channels in each DetId.
0392   // Signal Digis are in the list first, have to merge lists of hit strips on the fly,
0393   // add signals on duplicates later
0394 
0395   OneDetectorRawMap LocalRawMap;
0396 
0397   // Now, loop over hits and add them to the map in the proper sorted order
0398   // Note: We are assuming that the hits from the Signal events have been created in
0399   // "PreMix" mode, rather than in the standard ADC conversion routines.  If not, this
0400   // doesn't work at all.
0401 
0402   // At the moment, both Signal and Reconstituted PU hits have the same compression algorithm.
0403   // If this were different, and one needed gains, the conversion back to pulse height can only
0404   // be done in this routine.  So, yes, there is an extra loop over hits here in the current code,
0405   // because, in principle, one could convert to pulse height during the read/store phase.
0406 
0407   for (SiGlobalIndex::const_iterator IDet = SiHitStorage_.begin(); IDet != SiHitStorage_.end(); IDet++) {
0408     uint32_t detID = IDet->first;
0409 
0410     OneDetectorMap LocalMap = IDet->second;
0411 
0412     //loop over hit strips for this DetId, do conversion to pulse height, store.
0413 
0414     LocalRawMap.clear();
0415 
0416     OneDetectorMap::const_iterator iLocal = LocalMap.begin();
0417     for (; iLocal != LocalMap.end(); ++iLocal) {
0418       uint16_t currentStrip = iLocal->strip();
0419       float signal = float(iLocal->adc());
0420       if (iLocal->adc() == 1022)
0421         signal = 1500.;  // average values for overflows
0422       if (iLocal->adc() == 1023)
0423         signal = 3000.;
0424 
0425       //convert signals back to raw counts
0426 
0427       float ReSignal = signal * signal / 9.0;  // The PreMixing conversion is adc = sqrt(9.0*pulseHeight)
0428 
0429       RawDigi NewRawDigi = std::make_pair(currentStrip, ReSignal);
0430 
0431       LocalRawMap.push_back(NewRawDigi);
0432     }
0433 
0434     // save information for this detiD into global map
0435     SiRawDigis_.insert(SiGlobalRawIndex::value_type(detID, LocalRawMap));
0436   }
0437 
0438   // If we are killing APVs, merge list of dead ones before we digitize
0439 
0440   int NumberOfBxBetweenHIPandEvent = 1e3;
0441 
0442   if (APVSaturationFromHIP_) {
0443     // calculate affected BX parameter
0444 
0445     bool HasAtleastOneAffectedAPV = false;
0446     while (!HasAtleastOneAffectedAPV) {
0447       for (int bx = floor(300.0 / 25.0); bx > 0; bx--) {  //Reminder: make these numbers not hard coded!!
0448         float temp = CLHEP::RandFlat::shoot(engine) < 0.5 ? 1 : 0;
0449         if (temp == 1 && bx < NumberOfBxBetweenHIPandEvent) {
0450           NumberOfBxBetweenHIPandEvent = bx;
0451           HasAtleastOneAffectedAPV = true;
0452         }
0453       }
0454     }
0455 
0456     APVMap::const_iterator iAPVchk;
0457     uint32_t formerID = 0;
0458     uint32_t currentID;
0459     std::bitset<6> NewAPVBits;
0460 
0461     for (APVMap::const_iterator iAPV = theAffectedAPVmap_.begin(); iAPV != theAffectedAPVmap_.end(); ++iAPV) {
0462       currentID = iAPV->first;
0463 
0464       if (currentID == formerID) {  // we have to OR these
0465         for (int ibit = 0; ibit < 6; ++ibit) {
0466           NewAPVBits[ibit] = NewAPVBits[ibit] || (iAPV->second)[ibit];
0467         }
0468       } else {
0469         DeadAPVList[currentID] = NewAPVBits;
0470         //save pointers for next iteration
0471         formerID = currentID;
0472         NewAPVBits = iAPV->second;
0473       }
0474 
0475       iAPVchk = iAPV;
0476       if ((++iAPVchk) == theAffectedAPVmap_.end()) {  //make sure not to lose the last one
0477         DeadAPVList[currentID] = NewAPVBits;
0478       }
0479     }
0480   }
0481   //
0482 
0483   //  Ok, done with merging raw signals and APVs - now add signals on duplicate strips
0484 
0485   // collection of Digis to put in the event
0486   std::vector<edm::DetSet<SiStripDigi>> vSiStripDigi;
0487 
0488   // loop through our collection of detectors, merging hits and making a new list of "signal" digis
0489 
0490   // clear some temporary storage for later digitization:
0491 
0492   signals_.clear();
0493 
0494   // big loop over Detector IDs:
0495   for (SiGlobalRawIndex::const_iterator IDet = SiRawDigis_.begin(); IDet != SiRawDigis_.end(); IDet++) {
0496     uint32_t detID = IDet->first;
0497 
0498     SignalMapType Signals;
0499     Signals.clear();
0500 
0501     OneDetectorRawMap LocalMap = IDet->second;
0502 
0503     //counter variables
0504     int formerStrip = -1;
0505     int currentStrip;
0506     float ADCSum = 0;
0507 
0508     //loop over hit strips for this DetId, add duplicates
0509 
0510     OneDetectorRawMap::const_iterator iLocalchk;
0511     OneDetectorRawMap::const_iterator iLocal = LocalMap.begin();
0512     for (; iLocal != LocalMap.end(); ++iLocal) {
0513       currentStrip = iLocal->first;  // strip is first element
0514 
0515       if (currentStrip == formerStrip) {  // we have to add these digis together
0516 
0517         ADCSum += iLocal->second;  // raw pulse height is second element.
0518       } else {
0519         if (formerStrip != -1) {
0520           Signals.insert(std::make_pair(formerStrip, ADCSum));
0521         }
0522         // save pointers for next iteration
0523         formerStrip = currentStrip;
0524         ADCSum = iLocal->second;  // lone ADC
0525       }
0526 
0527       iLocalchk = iLocal;
0528       if ((++iLocalchk) == LocalMap.end()) {  //make sure not to lose the last one
0529         Signals.insert(std::make_pair(formerStrip, ADCSum));
0530       }
0531     }
0532     // save merged map:
0533     signals_.insert(std::make_pair(detID, Signals));
0534   }
0535 
0536   //Now, do noise, zero suppression, take into account bad channels, etc.
0537   // This section stolen from SiStripDigitizerAlgorithm
0538   // must loop over all detIds in the tracker to get all of the noise added properly.
0539   for (const auto& iu : pDD->detUnits()) {
0540     const StripGeomDetUnit* sgd = dynamic_cast<const StripGeomDetUnit*>(iu);
0541     if (sgd != nullptr) {
0542       uint32_t detID = sgd->geographicalId().rawId();
0543 
0544       edm::DetSet<SiStripDigi> SSD(detID);  // Make empty collection with this detector ID
0545 
0546       int numStrips = (sgd->specificTopology()).nstrips();
0547 
0548       // see if there is some signal on this detector
0549 
0550       const SignalMapType* theSignal(getSignal(detID));
0551 
0552       std::vector<float> detAmpl(numStrips, 0.);
0553       if (theSignal) {
0554         for (const auto& amp : *theSignal) {
0555           detAmpl[amp.first] = amp.second;
0556         }
0557       }
0558 
0559       //removing signal from the dead (and HIP effected) strips
0560       std::vector<bool>& badChannels = allBadChannels[detID];
0561 
0562       for (int strip = 0; strip < numStrips; ++strip) {
0563         if (badChannels[strip])
0564           detAmpl[strip] = 0.;
0565       }
0566 
0567       if (simulateAPVInThisEvent) {
0568         // Get index in apv baseline distributions corresponding to z of detSet and PU
0569         const StripTopology* topol = dynamic_cast<const StripTopology*>(&(sgd->specificTopology()));
0570         LocalPoint localPos = topol->localPosition(0);
0571         GlobalPoint globalPos = sgd->surface().toGlobal(Local3DPoint(localPos.x(), localPos.y(), localPos.z()));
0572         float detSet_z = fabs(globalPos.z());
0573         float detSet_r = globalPos.perp();
0574 
0575         const uint32_t SubDet = DetId(detID).subdetId();
0576         // Simulate APV response for each strip
0577         for (int strip = 0; strip < numStrips; ++strip) {
0578           if (detAmpl[strip] > 0) {
0579             // Convert charge from electrons to fC
0580             double stripCharge = detAmpl[strip] * apv_fCPerElectron_;
0581 
0582             // Get APV baseline
0583             double baselineV = 0;
0584             if (SubDet == SiStripSubdetector::TIB) {
0585               baselineV = apvSimulationParameters->sampleTIB(tTopo->tibLayer(detID), detSet_z, nTruePU, engine);
0586             } else if (SubDet == SiStripSubdetector::TOB) {
0587               baselineV = apvSimulationParameters->sampleTOB(tTopo->tobLayer(detID), detSet_z, nTruePU, engine);
0588             } else if (SubDet == SiStripSubdetector::TID) {
0589               baselineV = apvSimulationParameters->sampleTID(tTopo->tidWheel(detID), detSet_r, nTruePU, engine);
0590             } else if (SubDet == SiStripSubdetector::TEC) {
0591               baselineV = apvSimulationParameters->sampleTEC(tTopo->tecWheel(detID), detSet_r, nTruePU, engine);
0592             }
0593             // Fitted parameters from G Hall/M Raymond
0594             double maxResponse = apv_maxResponse_;
0595             double rate = apv_rate_;
0596 
0597             double outputChargeInADC = 0;
0598             if (baselineV < apv_maxResponse_) {
0599               // Convert V0 into baseline charge
0600               double baselineQ = -1.0 * rate * log(2 * maxResponse / (baselineV + maxResponse) - 1);
0601 
0602               // Add charge deposited in this BX
0603               double newStripCharge = baselineQ + stripCharge;
0604 
0605               // Apply APV response
0606               double signalV = 2 * maxResponse / (1 + exp(-1.0 * newStripCharge / rate)) - maxResponse;
0607               double gain = signalV - baselineV;
0608 
0609               // Convert gain (mV) to charge (assuming linear region of APV) and then to electrons
0610               double outputCharge = gain / apv_mVPerQ_;
0611               outputChargeInADC = outputCharge / apv_fCPerElectron_;
0612             }
0613 
0614             // Output charge back to original container
0615             detAmpl[strip] = outputChargeInADC;
0616           }
0617         }
0618       }
0619 
0620       if (APVSaturationFromHIP_) {
0621         std::bitset<6>& bs = DeadAPVList[detID];
0622 
0623         if (bs.any()) {
0624           // Here below is the scaling function which describes the evolution of the baseline (i.e. how the charge is suppressed).
0625           // This must be replaced as soon as we have a proper modeling of the baseline evolution from VR runs
0626           float Shift =
0627               1 - NumberOfBxBetweenHIPandEvent / floor(300.0 / 25.0);  //Reminder: make these numbers not hardcoded!!
0628           float randomX = CLHEP::RandFlat::shoot(engine);
0629           float scalingValue = (randomX - Shift) * 10.0 / 7.0 - 3.0 / 7.0;
0630 
0631           for (int strip = 0; strip < numStrips; ++strip) {
0632             if (!badChannels[strip] && bs[strip / 128] == 1) {
0633               detAmpl[strip] *= scalingValue > 0 ? scalingValue : 0.0;
0634             }
0635           }
0636         }
0637       }
0638 
0639       SiStripNoises::Range detNoiseRange = noise.getRange(detID);
0640       SiStripApvGain::Range detGainRange = gain.getRange(detID);
0641 
0642       // Gain conversion is already done during signal adding
0643       //convert our signals back to raw counts so that we can add noise properly:
0644 
0645       /*
0646         if(theSignal) {
0647         for(unsigned int iv = 0; iv!=detAmpl.size(); iv++) {
0648         float signal = detAmpl[iv];
0649         if(signal > 0) {
0650         float gainValue = gain.getStripGain(iv, detGainRange);
0651         signal *= theElectronPerADC/gainValue;
0652         detAmpl[iv] = signal;
0653         }
0654         }
0655         }
0656       */
0657 
0658       //SiStripPedestals::Range detPedestalRange = pedestal.getRange(detID);
0659 
0660       // -----------------------------------------------------------
0661 
0662       size_t firstChannelWithSignal = 0;
0663       size_t lastChannelWithSignal = numStrips;
0664 
0665       if (SingleStripNoise) {
0666         std::vector<float> noiseRMSv;
0667         noiseRMSv.clear();
0668         noiseRMSv.insert(noiseRMSv.begin(), numStrips, 0.);
0669         for (int strip = 0; strip < numStrips; ++strip) {
0670           if (!badChannels[strip]) {
0671             float gainValue = gain.getStripGain(strip, detGainRange);
0672             noiseRMSv[strip] = (noise.getNoise(strip, detNoiseRange)) * theElectronPerADC / gainValue;
0673           }
0674         }
0675         theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv, engine);
0676       } else {
0677         int RefStrip = int(numStrips / 2.);
0678         while (RefStrip < numStrips &&
0679                badChannels[RefStrip]) {  //if the refstrip is bad, I move up to when I don't find it
0680           RefStrip++;
0681         }
0682         if (RefStrip < numStrips) {
0683           float RefgainValue = gain.getStripGain(RefStrip, detGainRange);
0684           float RefnoiseRMS = noise.getNoise(RefStrip, detNoiseRange) * theElectronPerADC / RefgainValue;
0685 
0686           theSiNoiseAdder->addNoise(
0687               detAmpl, firstChannelWithSignal, lastChannelWithSignal, numStrips, RefnoiseRMS, engine);
0688         }
0689       }
0690 
0691       DigitalVecType digis;
0692       theSiZeroSuppress->suppress(
0693           theSiDigitalConverter->convert(detAmpl, &gain, detID), digis, detID, noise, threshold);
0694 
0695       SSD.data = digis;
0696 
0697       // stick this into the global vector of detector info
0698       vSiStripDigi.push_back(SSD);
0699 
0700     }  // end of loop over one detector
0701 
0702   }  // end of big loop over all detector IDs
0703 
0704   // put the collection of digis in the event
0705   edm::LogInfo("PreMixingSiStripWorker") << "total # Merged strips: " << vSiStripDigi.size();
0706 
0707   // make new digi collection
0708 
0709   std::unique_ptr<edm::DetSetVector<SiStripDigi>> MySiStripDigis(new edm::DetSetVector<SiStripDigi>(vSiStripDigi));
0710 
0711   // put collection
0712 
0713   e.put(std::move(MySiStripDigis), SiStripDigiCollectionDM_);
0714   e.put(std::make_unique<bool>(simulateAPVInThisEvent), SiStripDigiCollectionDM_ + "SimulatedAPVDynamicGain");
0715 
0716   // clear local storage for this event
0717   SiHitStorage_.clear();
0718   SiRawDigis_.clear();
0719   signals_.clear();
0720 }
0721 
0722 DEFINE_PREMIXING_WORKER(PreMixingSiStripWorker);