Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:25:42 

0001 // File: SiStripDigitizerAlgorithm.cc
0002 // Description:  Steering class for digitization.
0003 // Modified 15/May/2013 mark.grimes@bristol.ac.uk - Modified so that the digi-sim link has the correct
0004 // index for the sim hits stored. It was previously always set to zero (I won't mention that it was
0005 // me who originally wrote that).
0006 // Modified on Feb 11, 2015: prolay.kumar.mal@cern.ch & Jean-Laurent.Agram@cern.ch
0007 //                           Added/Modified the individual strip noise in zero suppression
0008 //                           mode from the conditions DB; previously, the digitizer used to
0009 //                           consider the noise value for individual strips inside a module from
0010 //                           the central strip noise value.
0011 //////////////////////////////////////////////////////////////////////////////////////////////////////////
0012 #include <vector>
0013 #include <algorithm>
0014 #include <iostream>
0015 #include "FWCore/Framework/interface/ConsumesCollector.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "SiStripDigitizerAlgorithm.h"
0020 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0021 #include "SimDataFormats/TrackerDigiSimLink/interface/StripDigiSimLink.h"
0022 #include "DataFormats/Common/interface/DetSetVector.h"
0023 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
0024 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0025 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
0026 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
0027 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0028 #include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
0029 #include "CondFormats/SiStripObjects/interface/SiStripPedestals.h"
0030 #include "CondFormats/SiStripObjects/interface/SiStripBadStrip.h"
0031 #include "CLHEP/Random/RandFlat.h"
0032 
0033 #include <cstring>
0034 #include <sstream>
0035 
0036 #include <boost/algorithm/string.hpp>
0037 
0038 SiStripDigitizerAlgorithm::SiStripDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0039     : theThreshold(conf.getParameter<double>("NoiseSigmaThreshold")),
0040       cmnRMStib(conf.getParameter<double>("cmnRMStib")),
0041       cmnRMStob(conf.getParameter<double>("cmnRMStob")),
0042       cmnRMStid(conf.getParameter<double>("cmnRMStid")),
0043       cmnRMStec(conf.getParameter<double>("cmnRMStec")),
0044       APVSaturationProbScaling_(conf.getParameter<double>("APVSaturationProbScaling")),
0045       makeDigiSimLinks_(conf.getUntrackedParameter<bool>("makeDigiSimLinks", false)),
0046       peakMode(conf.getParameter<bool>("APVpeakmode")),
0047       noise(conf.getParameter<bool>("Noise")),
0048       RealPedestals(conf.getParameter<bool>("RealPedestals")),
0049       SingleStripNoise(conf.getParameter<bool>("SingleStripNoise")),
0050       CommonModeNoise(conf.getParameter<bool>("CommonModeNoise")),
0051       BaselineShift(conf.getParameter<bool>("BaselineShift")),
0052       APVSaturationFromHIP(conf.getParameter<bool>("APVSaturationFromHIP")),
0053       theFedAlgo(conf.getParameter<int>("FedAlgorithm")),
0054       zeroSuppression(conf.getParameter<bool>("ZeroSuppression")),
0055       theElectronPerADC(conf.getParameter<double>(peakMode ? "electronPerAdcPeak" : "electronPerAdcDec")),
0056       theTOFCutForPeak(conf.getParameter<double>("TOFCutForPeak")),
0057       theTOFCutForDeconvolution(conf.getParameter<double>("TOFCutForDeconvolution")),
0058       tofCut(peakMode ? theTOFCutForPeak : theTOFCutForDeconvolution),
0059       cosmicShift(conf.getUntrackedParameter<double>("CosmicDelayShift")),
0060       inefficiency(conf.getParameter<double>("Inefficiency")),
0061       pedOffset((unsigned int)conf.getParameter<double>("PedestalsOffset")),
0062       PreMixing_(conf.getParameter<bool>("PreMixingMode")),
0063       pdtToken_(iC.esConsumes()),
0064       lorentzAngleToken_(iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("LorentzAngle")))),
0065       theSiHitDigitizer(new SiHitDigitizer(conf)),
0066       theSiPileUpSignals(new SiPileUpSignals()),
0067       theSiNoiseAdder(new SiGaussianTailNoiseAdder(theThreshold)),
0068       theSiDigitalConverter(new SiTrivialDigitalConverter(theElectronPerADC, PreMixing_)),
0069       theSiZeroSuppress(new SiStripFedZeroSuppression(theFedAlgo, &iC)),
0070       APVProbabilityFile(conf.getParameter<edm::FileInPath>("APVProbabilityFile")),
0071       includeAPVSimulation_(conf.getParameter<bool>("includeAPVSimulation")),
0072       apv_maxResponse_(conf.getParameter<double>("apv_maxResponse")),
0073       apv_rate_(conf.getParameter<double>("apv_rate")),
0074       apv_mVPerQ_(conf.getParameter<double>("apv_mVPerQ")),
0075       apv_fCPerElectron_(conf.getParameter<double>("apvfCPerElectron")) {
0076   if (peakMode) {
0077     LogDebug("StripDigiInfo") << "APVs running in peak mode (poor time resolution)";
0078   } else {
0079     LogDebug("StripDigiInfo") << "APVs running in deconvolution mode (good time resolution)";
0080   };
0081   if (SingleStripNoise)
0082     LogDebug("SiStripDigitizerAlgorithm") << " SingleStripNoise: ON";
0083   else
0084     LogDebug("SiStripDigitizerAlgorithm") << " SingleStripNoise: OFF";
0085   if (CommonModeNoise)
0086     LogDebug("SiStripDigitizerAlgorithm") << " CommonModeNoise: ON";
0087   else
0088     LogDebug("SiStripDigitizerAlgorithm") << " CommonModeNoise: OFF";
0089   if (PreMixing_ && APVSaturationFromHIP)
0090     throw cms::Exception("PreMixing does not work with HIP loss simulation yet");
0091   if (APVSaturationFromHIP) {
0092     std::string line;
0093     APVProbaFile.open((APVProbabilityFile.fullPath()).c_str());
0094     if (APVProbaFile.is_open()) {
0095       while (getline(APVProbaFile, line)) {
0096         std::vector<std::string> strs;
0097         boost::split(strs, line, boost::is_any_of(" "));
0098         if (strs.size() == 2) {
0099           mapOfAPVprobabilities[std::stoi(strs.at(0))] = std::stof(strs.at(1));
0100         }
0101       }
0102       APVProbaFile.close();
0103     } else
0104       throw cms::Exception("MissingInput") << "It seems that the APV probability list is missing\n";
0105   }
0106 }
0107 
0108 SiStripDigitizerAlgorithm::~SiStripDigitizerAlgorithm() {}
0109 
0110 void SiStripDigitizerAlgorithm::initializeDetUnit(StripGeomDetUnit const* det, const SiStripBadStrip& deadChannel) {
0111   unsigned int detId = det->geographicalId().rawId();
0112   int numStrips = (det->specificTopology()).nstrips();
0113 
0114   SiStripBadStrip::Range detBadStripRange = deadChannel.getRange(detId);
0115   //storing the bad strip of the the module. the module is not removed but just signal put to 0
0116   std::vector<bool>& badChannels = allBadChannels[detId];
0117   badChannels.clear();
0118   badChannels.insert(badChannels.begin(), numStrips, false);
0119   for (SiStripBadStrip::ContainerIterator it = detBadStripRange.first; it != detBadStripRange.second; ++it) {
0120     SiStripBadStrip::data fs = deadChannel.decode(*it);
0121     for (int strip = fs.firstStrip; strip < fs.firstStrip + fs.range; ++strip) {
0122       badChannels[strip] = true;
0123     }
0124   }
0125   firstChannelsWithSignal[detId] = numStrips;
0126   lastChannelsWithSignal[detId] = 0;
0127 
0128   //  if(APVSaturationFromHIP){
0129   //  std::bitset<6> &bs=SiStripTrackerAffectedAPVMap[detId];
0130   //  if(bs.any())theAffectedAPVvector.push_back(std::make_pair(detId,bs));
0131   //}
0132 }
0133 
0134 void SiStripDigitizerAlgorithm::initializeEvent(const edm::EventSetup& iSetup) {
0135   theSiPileUpSignals->reset();
0136   // This should be clear by after all calls to digitize(), but I might as well make sure
0137   associationInfoForDetId_.clear();
0138 
0139   APVSaturationProb_ = APVSaturationProbScaling_;  // reset probability
0140   SiStripTrackerAffectedAPVMap.clear();
0141   FirstLumiCalc_ = true;
0142   FirstDigitize_ = true;
0143 
0144   nTruePU_ = 0;
0145 
0146   //get gain noise pedestal lorentzAngle from ES handle
0147   setParticleDataTable(&iSetup.getData(pdtToken_));
0148   lorentzAngleHandle = iSetup.getHandle(lorentzAngleToken_);
0149 }
0150 
0151 //  Run the algorithm for a given module
0152 //  ------------------------------------
0153 
0154 void SiStripDigitizerAlgorithm::accumulateSimHits(std::vector<PSimHit>::const_iterator inputBegin,
0155                                                   std::vector<PSimHit>::const_iterator inputEnd,
0156                                                   size_t inputBeginGlobalIndex,
0157                                                   unsigned int tofBin,
0158                                                   const StripGeomDetUnit* det,
0159                                                   const GlobalVector& bfield,
0160                                                   const TrackerTopology* tTopo,
0161                                                   CLHEP::HepRandomEngine* engine) {
0162   // produce SignalPoints for all SimHits in detector
0163   unsigned int detID = det->geographicalId().rawId();
0164   int numStrips = (det->specificTopology()).nstrips();
0165 
0166   size_t thisFirstChannelWithSignal = numStrips;
0167   size_t thisLastChannelWithSignal = 0;
0168 
0169   float langle = (lorentzAngleHandle.isValid()) ? lorentzAngleHandle->getLorentzAngle(detID) : 0.;
0170 
0171   std::vector<float> locAmpl(numStrips, 0.);
0172 
0173   // Loop over hits
0174 
0175   uint32_t detId = det->geographicalId().rawId();
0176   // First: loop on the SimHits
0177   if (CLHEP::RandFlat::shoot(engine) > inefficiency) {
0178     AssociationInfoForChannel* pDetIDAssociationInfo;  // I only need this if makeDigiSimLinks_ is true...
0179     if (makeDigiSimLinks_)
0180       pDetIDAssociationInfo = &(associationInfoForDetId_[detId]);  // ...so only search the map if that is the case
0181     std::vector<float>
0182         previousLocalAmplitude;  // Only used if makeDigiSimLinks_ is true. Needed to work out the change in amplitude.
0183 
0184     size_t simHitGlobalIndex = inputBeginGlobalIndex;  // This needs to stored to create the digi-sim link later
0185     for (std::vector<PSimHit>::const_iterator simHitIter = inputBegin; simHitIter != inputEnd;
0186          ++simHitIter, ++simHitGlobalIndex) {
0187       // skip hits not in this detector.
0188       if ((*simHitIter).detUnitId() != detId) {
0189         continue;
0190       }
0191       // check TOF
0192       if (std::fabs(simHitIter->tof() - cosmicShift -
0193                     det->surface().toGlobal(simHitIter->localPosition()).mag() / 30.) < tofCut &&
0194           simHitIter->energyLoss() > 0) {
0195         if (makeDigiSimLinks_)
0196           previousLocalAmplitude = locAmpl;  // Not needed except to make the sim link association.
0197         size_t localFirstChannel = numStrips;
0198         size_t localLastChannel = 0;
0199         // process the hit
0200         theSiHitDigitizer->processHit(
0201             &*simHitIter, *det, bfield, langle, locAmpl, localFirstChannel, localLastChannel, tTopo, engine);
0202 
0203         if (thisFirstChannelWithSignal > localFirstChannel)
0204           thisFirstChannelWithSignal = localFirstChannel;
0205         if (thisLastChannelWithSignal < localLastChannel)
0206           thisLastChannelWithSignal = localLastChannel;
0207 
0208         if (makeDigiSimLinks_) {  // No need to do any of this if truth association was turned off in the configuration
0209           for (size_t stripIndex = 0; stripIndex < locAmpl.size(); ++stripIndex) {
0210             // Work out the amplitude from this SimHit from the difference of what it was before and what it is now
0211             float signalFromThisSimHit = locAmpl[stripIndex] - previousLocalAmplitude[stripIndex];
0212             if (signalFromThisSimHit != 0) {  // If this SimHit had any contribution I need to record it.
0213               auto& associationVector = (*pDetIDAssociationInfo)[stripIndex];
0214               bool addNewEntry = true;
0215               // Make sure the hit isn't in already. I've seen this a few times, it always seems to happen in pairs so I think
0216               // it's something to do with the stereo strips.
0217               for (auto& associationInfo : associationVector) {
0218                 if (associationInfo.trackID == simHitIter->trackId() &&
0219                     associationInfo.eventID == simHitIter->eventId()) {
0220                   // The hit is already in, so add this second contribution and move on
0221                   associationInfo.contributionToADC += signalFromThisSimHit;
0222                   addNewEntry = false;
0223                   break;
0224                 }
0225               }  // end of loop over associationVector
0226               // If the hit wasn't already in create a new association info structure.
0227               if (addNewEntry)
0228                 associationVector.push_back(AssociationInfo{
0229                     simHitIter->trackId(), simHitIter->eventId(), signalFromThisSimHit, simHitGlobalIndex, tofBin});
0230             }  // end of "if( signalFromThisSimHit!=0 )"
0231           }    // end of loop over locAmpl strips
0232         }      // end of "if( makeDigiSimLinks_ )"
0233       }        // end of TOF check
0234     }          // end for
0235   }
0236   theSiPileUpSignals->add(detID, locAmpl, thisFirstChannelWithSignal, thisLastChannelWithSignal);
0237 
0238   if (firstChannelsWithSignal[detID] > thisFirstChannelWithSignal)
0239     firstChannelsWithSignal[detID] = thisFirstChannelWithSignal;
0240   if (lastChannelsWithSignal[detID] < thisLastChannelWithSignal)
0241     lastChannelsWithSignal[detID] = thisLastChannelWithSignal;
0242 }
0243 
0244 //============================================================================
0245 void SiStripDigitizerAlgorithm::calculateInstlumiScale(PileupMixingContent* puInfo) {
0246   //Instlumi scalefactor calculating for dynamic inefficiency
0247 
0248   if (puInfo && FirstLumiCalc_) {
0249     const std::vector<int>& bunchCrossing = puInfo->getMix_bunchCrossing();
0250     const std::vector<float>& TrueInteractionList = puInfo->getMix_TrueInteractions();
0251     const int bunchSpacing = puInfo->getMix_bunchSpacing();
0252 
0253     double RevFreq = 11245.;
0254     double minBXsec = 70.0E-27;  // use 70mb as an approximation
0255     double Bunch = 2100.;        // 2016 value
0256     if (bunchSpacing == 50)
0257       Bunch = Bunch / 2.;
0258 
0259     int pui = 0, p = 0;
0260     std::vector<int>::const_iterator pu;
0261     std::vector<int>::const_iterator pu0 = bunchCrossing.end();
0262 
0263     for (pu = bunchCrossing.begin(); pu != bunchCrossing.end(); ++pu) {
0264       if (*pu == 0) {
0265         pu0 = pu;
0266         p = pui;
0267       }
0268       pui++;
0269     }
0270     if (pu0 != bunchCrossing.end()) {  // found the in-time interaction
0271       nTruePU_ = TrueInteractionList.at(p);
0272       double instLumi = Bunch * nTruePU_ * RevFreq / minBXsec;
0273       APVSaturationProb_ = instLumi / 6.0E33;
0274     }
0275     FirstLumiCalc_ = false;
0276   }
0277 }
0278 
0279 //============================================================================
0280 
0281 void SiStripDigitizerAlgorithm::digitize(edm::DetSet<SiStripDigi>& outdigi,
0282                                          edm::DetSet<SiStripRawDigi>& outrawdigi,
0283                                          edm::DetSet<SiStripRawDigi>& outStripAmplitudes,
0284                                          edm::DetSet<SiStripRawDigi>& outStripAmplitudesPostAPV,
0285                                          edm::DetSet<SiStripRawDigi>& outStripAPVBaselines,
0286                                          edm::DetSet<StripDigiSimLink>& outLink,
0287                                          const StripGeomDetUnit* det,
0288                                          const SiStripGain& gain,
0289                                          const SiStripThreshold& threshold,
0290                                          const SiStripNoises& noiseObj,
0291                                          const SiStripPedestals& pedestal,
0292                                          bool simulateAPVInThisEvent,
0293                                          const SiStripApvSimulationParameters* apvSimulationParameters,
0294                                          std::vector<std::pair<int, std::bitset<6>>>& theAffectedAPVvector,
0295                                          CLHEP::HepRandomEngine* engine,
0296                                          const TrackerTopology* tTopo) {
0297   unsigned int detID = det->geographicalId().rawId();
0298   int numStrips = (det->specificTopology()).nstrips();
0299 
0300   DetId detId(detID);
0301   uint32_t SubDet = detId.subdetId();
0302 
0303   const SiPileUpSignals::SignalMapType* theSignal(theSiPileUpSignals->getSignal(detID));
0304 
0305   std::vector<float> detAmpl(numStrips, 0.);
0306   if (theSignal) {
0307     for (const auto& amp : *theSignal) {
0308       detAmpl[amp.first] = amp.second;
0309     }
0310   }
0311 
0312   //removing signal from the dead (and HIP effected) strips
0313   std::vector<bool>& badChannels = allBadChannels[detID];
0314   for (int strip = 0; strip < numStrips; ++strip) {
0315     if (badChannels[strip]) {
0316       detAmpl[strip] = 0.;
0317     }
0318   }
0319 
0320   if (includeAPVSimulation_ && simulateAPVInThisEvent) {
0321     // Get index in apv baseline distributions corresponding to z of detSet and PU
0322     const StripTopology* topol = dynamic_cast<const StripTopology*>(&(det->specificTopology()));
0323     LocalPoint localPos = topol->localPosition(0);
0324     GlobalPoint globalPos = det->surface().toGlobal(Local3DPoint(localPos.x(), localPos.y(), localPos.z()));
0325     float detSet_z = fabs(globalPos.z());
0326     float detSet_r = globalPos.perp();
0327 
0328     // Store SCD, before APV sim
0329     outStripAmplitudes.reserve(numStrips);
0330     for (int strip = 0; strip < numStrips; ++strip) {
0331       outStripAmplitudes.emplace_back(SiStripRawDigi(detAmpl[strip] / theElectronPerADC));
0332     }
0333 
0334     // Simulate APV response for each strip
0335     for (int strip = 0; strip < numStrips; ++strip) {
0336       if (detAmpl[strip] > 0) {
0337         // Convert charge from electrons to fC
0338         double stripCharge = detAmpl[strip] * apv_fCPerElectron_;
0339 
0340         // Get APV baseline
0341         double baselineV = 0;
0342         if (SubDet == SiStripSubdetector::TIB) {
0343           baselineV = apvSimulationParameters->sampleTIB(tTopo->tibLayer(detId), detSet_z, nTruePU_, engine);
0344         } else if (SubDet == SiStripSubdetector::TOB) {
0345           baselineV = apvSimulationParameters->sampleTOB(tTopo->tobLayer(detId), detSet_z, nTruePU_, engine);
0346         } else if (SubDet == SiStripSubdetector::TID) {
0347           baselineV = apvSimulationParameters->sampleTID(tTopo->tidWheel(detId), detSet_r, nTruePU_, engine);
0348         } else if (SubDet == SiStripSubdetector::TEC) {
0349           baselineV = apvSimulationParameters->sampleTEC(tTopo->tecWheel(detId), detSet_r, nTruePU_, engine);
0350         }
0351         // Store APV baseline for this strip
0352         outStripAPVBaselines.emplace_back(SiStripRawDigi(baselineV));
0353 
0354         // Fitted parameters from G Hall/M Raymond
0355         double maxResponse = apv_maxResponse_;
0356         double rate = apv_rate_;
0357 
0358         double outputChargeInADC = 0;
0359         if (baselineV < apv_maxResponse_) {
0360           // Convert V0 into baseline charge
0361           double baselineQ = -1.0 * rate * log(2 * maxResponse / (baselineV + maxResponse) - 1);
0362 
0363           // Add charge deposited in this BX
0364           double newStripCharge = baselineQ + stripCharge;
0365 
0366           // Apply APV response
0367           double signalV = 2 * maxResponse / (1 + exp(-1.0 * newStripCharge / rate)) - maxResponse;
0368           double gain = signalV - baselineV;
0369 
0370           // Convert gain (mV) to charge (assuming linear region of APV) and then to electrons
0371           double outputCharge = gain / apv_mVPerQ_;
0372           outputChargeInADC = outputCharge / apv_fCPerElectron_;
0373         }
0374 
0375         // Output charge back to original container
0376         detAmpl[strip] = outputChargeInADC;
0377       }
0378     }
0379 
0380     // Store SCD, after APV sim
0381     outStripAmplitudesPostAPV.reserve(numStrips);
0382     for (int strip = 0; strip < numStrips; ++strip)
0383       outStripAmplitudesPostAPV.emplace_back(SiStripRawDigi(detAmpl[strip] / theElectronPerADC));
0384   }
0385 
0386   if (APVSaturationFromHIP) {
0387     //Implementation of the proper charge scaling function. Need consider resaturation effect:
0388     //The probability map gives  the probability that at least one HIP happened during the last N bunch crossings (cfr APV recovery time).
0389     //The impact on the charge depends on the clostest HIP occurance (in terms of bunch crossing).
0390     //The function discribing the APV recovery is therefore the weighted average function which takes into account all possibilities of HIP occurances across the last bx's.
0391 
0392     // do this step here because we now have access to luminosity information
0393     if (FirstDigitize_) {
0394       for (std::map<int, float>::iterator iter = mapOfAPVprobabilities.begin(); iter != mapOfAPVprobabilities.end();
0395            ++iter) {
0396         std::bitset<6> bs;
0397         for (int Napv = 0; Napv < 6; Napv++) {
0398           float cursor = CLHEP::RandFlat::shoot(engine);
0399           bs[Napv] = cursor < iter->second * APVSaturationProb_
0400                          ? true
0401                          : false;  //APVSaturationProb has been scaled by PU luminosity
0402         }
0403         SiStripTrackerAffectedAPVMap[iter->first] = bs;
0404       }
0405 
0406       NumberOfBxBetweenHIPandEvent = 1e3;
0407       bool HasAtleastOneAffectedAPV = false;
0408       while (!HasAtleastOneAffectedAPV) {
0409         for (int bx = floor(300.0 / 25.0); bx > 0; bx--) {  //Reminder: make these numbers not hard coded!!
0410           float temp = CLHEP::RandFlat::shoot(engine) < 0.5 ? 1 : 0;
0411           if (temp == 1 && bx < NumberOfBxBetweenHIPandEvent) {
0412             NumberOfBxBetweenHIPandEvent = bx;
0413             HasAtleastOneAffectedAPV = true;
0414           }
0415         }
0416       }
0417 
0418       FirstDigitize_ = false;
0419     }
0420 
0421     std::bitset<6>& bs = SiStripTrackerAffectedAPVMap[detID];
0422 
0423     if (bs.any()) {
0424       // store this information so it can be saved to the event later
0425       theAffectedAPVvector.push_back(std::make_pair(detID, bs));
0426 
0427       if (!PreMixing_) {
0428         // Here below is the scaling function which describes the evolution of the baseline (i.e. how the charge is suppressed).
0429         // This must be replaced as soon as we have a proper modeling of the baseline evolution from VR runs
0430         float Shift =
0431             1 - NumberOfBxBetweenHIPandEvent / floor(300.0 / 25.0);  //Reminder: make these numbers not hardcoded!!
0432         float randomX = CLHEP::RandFlat::shoot(engine);
0433         float scalingValue = (randomX - Shift) * 10.0 / 7.0 - 3.0 / 7.0;
0434 
0435         for (int strip = 0; strip < numStrips; ++strip) {
0436           if (!badChannels[strip] && bs[strip / 128] == 1) {
0437             detAmpl[strip] *= scalingValue > 0 ? scalingValue : 0.0;
0438           }
0439         }
0440       }
0441     }
0442   }
0443 
0444   SiStripNoises::Range detNoiseRange = noiseObj.getRange(detID);
0445   SiStripApvGain::Range detGainRange = gain.getRange(detID);
0446   SiStripPedestals::Range detPedestalRange = pedestal.getRange(detID);
0447 
0448   // -----------------------------------------------------------
0449 
0450   auto& firstChannelWithSignal = firstChannelsWithSignal[detID];
0451   auto& lastChannelWithSignal = lastChannelsWithSignal[detID];
0452   auto iAssociationInfoByChannel =
0453       associationInfoForDetId_.find(detID);  // Use an iterator so that I can easily remove it once finished
0454 
0455   if (zeroSuppression) {
0456     //Adding the strip noise
0457     //------------------------------------------------------
0458     if (noise) {
0459       if (SingleStripNoise) {
0460         //      std::cout<<"In SSN, detId="<<detID<<std::endl;
0461         std::vector<float> noiseRMSv;
0462         noiseRMSv.clear();
0463         noiseRMSv.insert(noiseRMSv.begin(), numStrips, 0.);
0464         for (int strip = 0; strip < numStrips; ++strip) {
0465           if (!badChannels[strip]) {
0466             float gainValue = gain.getStripGain(strip, detGainRange);
0467             noiseRMSv[strip] = (noiseObj.getNoise(strip, detNoiseRange)) * theElectronPerADC / gainValue;
0468             //std::cout<<"<SiStripDigitizerAlgorithm::digitize>: gainValue: "<<gainValue<<"\tnoiseRMSv["<<strip<<"]: "<<noiseRMSv[strip]<<std::endl;
0469           }
0470         }
0471         theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv, engine);
0472       } else {
0473         int RefStrip = int(numStrips / 2.);
0474         while (RefStrip < numStrips &&
0475                badChannels[RefStrip]) {  //if the refstrip is bad, I move up to when I don't find it
0476           RefStrip++;
0477         }
0478         if (RefStrip < numStrips) {
0479           float RefgainValue = gain.getStripGain(RefStrip, detGainRange);
0480           float RefnoiseRMS = noiseObj.getNoise(RefStrip, detNoiseRange) * theElectronPerADC / RefgainValue;
0481 
0482           theSiNoiseAdder->addNoise(
0483               detAmpl, firstChannelWithSignal, lastChannelWithSignal, numStrips, RefnoiseRMS, engine);
0484           //std::cout<<"<SiStripDigitizerAlgorithm::digitize>: RefgainValue: "<<RefgainValue<<"\tRefnoiseRMS: "<<RefnoiseRMS<<std::endl;
0485         }
0486       }
0487     }  //if noise
0488 
0489     DigitalVecType digis;
0490     theSiZeroSuppress->suppress(
0491         theSiDigitalConverter->convert(detAmpl, &gain, detID), digis, detID, noiseObj, threshold);
0492     // Now do the association to truth. Note that if truth association was turned off in the configuration this map
0493     // will be empty and the iterator will always equal associationInfoForDetId_.end().
0494     if (iAssociationInfoByChannel !=
0495         associationInfoForDetId_.end()) {  // make sure the readings for this DetID aren't completely from noise
0496       for (const auto& iDigi : digis) {
0497         auto& associationInfoByChannel = iAssociationInfoByChannel->second;
0498         const std::vector<AssociationInfo>& associationInfo = associationInfoByChannel[iDigi.channel()];
0499 
0500         // Need to find the total from all sim hits, because this might not be the same as the total
0501         // digitised due to noise or whatever.
0502         float totalSimADC = 0;
0503         for (const auto& iAssociationInfo : associationInfo)
0504           totalSimADC += iAssociationInfo.contributionToADC;
0505         // Now I know that I can loop again and create the links
0506         for (const auto& iAssociationInfo : associationInfo) {
0507           // Note simHitGlobalIndex used to have +1 because TrackerHitAssociator (the only place I can find this value being used)
0508           // expected counting to start at 1, not 0.  Now changed.
0509           outLink.push_back(StripDigiSimLink(iDigi.channel(),
0510                                              iAssociationInfo.trackID,
0511                                              iAssociationInfo.simHitGlobalIndex,
0512                                              iAssociationInfo.tofBin,
0513                                              iAssociationInfo.eventID,
0514                                              iAssociationInfo.contributionToADC / totalSimADC));
0515         }  // end of loop over associationInfo
0516       }    // end of loop over the digis
0517     }      // end of check that iAssociationInfoByChannel is a valid iterator
0518     outdigi.data = digis;
0519   }  //if zeroSuppression
0520 
0521   if (!zeroSuppression) {
0522     //if(noise){
0523     // the constant pedestal offset is needed because
0524     //   negative adc counts are not allowed in case
0525     //   Pedestal and CMN subtraction is performed.
0526     //   The pedestal value read from the conditions
0527     //   is pedValue and after the pedestal subtraction
0528     //   the baseline is zero. The Common Mode Noise
0529     //   is not subtracted from the negative adc counts
0530     //   channels. Adding pedOffset the baseline is set
0531     //   to pedOffset after pedestal subtraction and CMN
0532     //   is subtracted to all the channels since none of
0533     //   them has negative adc value. The pedOffset is
0534     //   treated as a constant component in the CMN
0535     //   estimation and subtracted as CMN.
0536 
0537     //calculating the charge deposited on each APV and subtracting the shift
0538     //------------------------------------------------------
0539     if (BaselineShift) {
0540       theSiNoiseAdder->addBaselineShift(detAmpl, badChannels);
0541     }
0542 
0543     //Adding the strip noise
0544     //------------------------------------------------------
0545     if (noise) {
0546       std::vector<float> noiseRMSv;
0547       noiseRMSv.clear();
0548       noiseRMSv.insert(noiseRMSv.begin(), numStrips, 0.);
0549 
0550       if (SingleStripNoise) {
0551         for (int strip = 0; strip < numStrips; ++strip) {
0552           if (!badChannels[strip])
0553             noiseRMSv[strip] = (noiseObj.getNoise(strip, detNoiseRange)) * theElectronPerADC;
0554         }
0555 
0556       } else {
0557         int RefStrip = 0;  //int(numStrips/2.);
0558         while (RefStrip < numStrips &&
0559                badChannels[RefStrip]) {  //if the refstrip is bad, I move up to when I don't find it
0560           RefStrip++;
0561         }
0562         if (RefStrip < numStrips) {
0563           float noiseRMS = noiseObj.getNoise(RefStrip, detNoiseRange) * theElectronPerADC;
0564           for (int strip = 0; strip < numStrips; ++strip) {
0565             if (!badChannels[strip])
0566               noiseRMSv[strip] = noiseRMS;
0567           }
0568         }
0569       }
0570 
0571       theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv, engine);
0572     }
0573 
0574     //adding the CMN
0575     //------------------------------------------------------
0576     if (CommonModeNoise) {
0577       float cmnRMS = 0.;
0578       DetId detId(detID);
0579       switch (detId.subdetId()) {
0580         case SiStripSubdetector::TIB:
0581           cmnRMS = cmnRMStib;
0582           break;
0583         case SiStripSubdetector::TID:
0584           cmnRMS = cmnRMStid;
0585           break;
0586         case SiStripSubdetector::TOB:
0587           cmnRMS = cmnRMStob;
0588           break;
0589         case SiStripSubdetector::TEC:
0590           cmnRMS = cmnRMStec;
0591           break;
0592       }
0593       cmnRMS *= theElectronPerADC;
0594       theSiNoiseAdder->addCMNoise(detAmpl, cmnRMS, badChannels, engine);
0595     }
0596 
0597     //Adding the pedestals
0598     //------------------------------------------------------
0599 
0600     std::vector<float> vPeds;
0601     vPeds.clear();
0602     vPeds.insert(vPeds.begin(), numStrips, 0.);
0603 
0604     if (RealPedestals) {
0605       for (int strip = 0; strip < numStrips; ++strip) {
0606         if (!badChannels[strip])
0607           vPeds[strip] = (pedestal.getPed(strip, detPedestalRange) + pedOffset) * theElectronPerADC;
0608       }
0609     } else {
0610       for (int strip = 0; strip < numStrips; ++strip) {
0611         if (!badChannels[strip])
0612           vPeds[strip] = pedOffset * theElectronPerADC;
0613       }
0614     }
0615 
0616     theSiNoiseAdder->addPedestals(detAmpl, vPeds);
0617 
0618     //if(!RealPedestals&&!CommonModeNoise&&!noise&&!BaselineShift&&!APVSaturationFromHIP){
0619 
0620     //  edm::LogWarning("SiStripDigitizer")<<"You are running the digitizer without Noise generation and without applying Zero Suppression. ARE YOU SURE???";
0621     //}else{
0622 
0623     DigitalRawVecType rawdigis = theSiDigitalConverter->convertRaw(detAmpl, &gain, detID);
0624 
0625     // Now do the association to truth. Note that if truth association was turned off in the configuration this map
0626     // will be empty and the iterator will always equal associationInfoForDetId_.end().
0627     if (iAssociationInfoByChannel !=
0628         associationInfoForDetId_.end()) {  // make sure the readings for this DetID aren't completely from noise
0629       // N.B. For the raw digis the channel is inferred from the position in the vector.
0630       // I'VE NOT TESTED THIS YET!!!!!
0631       // ToDo Test this properly.
0632       for (size_t channel = 0; channel < rawdigis.size(); ++channel) {
0633         auto& associationInfoByChannel = iAssociationInfoByChannel->second;
0634         const auto iAssociationInfo = associationInfoByChannel.find(channel);
0635         if (iAssociationInfo == associationInfoByChannel.end())
0636           continue;  // Skip if there is no sim information for this channel (i.e. it's all noise)
0637         const std::vector<AssociationInfo>& associationInfo = iAssociationInfo->second;
0638 
0639         // Need to find the total from all sim hits, because this might not be the same as the total
0640         // digitised due to noise or whatever.
0641         float totalSimADC = 0;
0642         for (const auto& iAssociationInfo : associationInfo)
0643           totalSimADC += iAssociationInfo.contributionToADC;
0644         // Now I know that I can loop again and create the links
0645         for (const auto& iAssociationInfo : associationInfo) {
0646           // Note simHitGlobalIndex used to have +1 because TrackerHitAssociator (the only place I can find this value being used)
0647           // expected counting to start at 1, not 0.  Now changed.
0648           outLink.push_back(StripDigiSimLink(channel,
0649                                              iAssociationInfo.trackID,
0650                                              iAssociationInfo.simHitGlobalIndex,
0651                                              iAssociationInfo.tofBin,
0652                                              iAssociationInfo.eventID,
0653                                              iAssociationInfo.contributionToADC / totalSimADC));
0654         }  // end of loop over associationInfo
0655       }    // end of loop over the digis
0656     }      // end of check that iAssociationInfoByChannel is a valid iterator
0657 
0658     outrawdigi.data = rawdigis;
0659 
0660     //}
0661   }
0662 
0663   // Now that I've finished with this entry in the map of associations, I can remove it.
0664   // Note that there might not be an association if the ADC reading is from noise in which
0665   // case associationIsValid will be false.
0666   if (iAssociationInfoByChannel != associationInfoForDetId_.end())
0667     associationInfoForDetId_.erase(iAssociationInfoByChannel);
0668 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!