Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // File: DigiSimLinkAlgorithm.cc
0002 // Description:  Steering class for digitization.
0003 
0004 #include <vector>
0005 #include <algorithm>
0006 #include <iostream>
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "DigiSimLinkAlgorithm.h"
0010 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0011 #include "DataFormats/Common/interface/DetSetVector.h"
0012 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
0013 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0014 
0015 #include "CLHEP/Random/RandFlat.h"
0016 
0017 #define CBOLTZ (1.38E-23)
0018 #define e_SI (1.6E-19)
0019 
0020 DigiSimLinkAlgorithm::DigiSimLinkAlgorithm(const edm::ParameterSet& conf) : conf_(conf) {
0021   theThreshold = conf_.getParameter<double>("NoiseSigmaThreshold");
0022   theFedAlgo = conf_.getParameter<int>("FedAlgorithm");
0023   peakMode = conf_.getParameter<bool>("APVpeakmode");
0024   theElectronPerADC = conf_.getParameter<double>(peakMode ? "electronPerAdcPeak" : "electronPerAdcDec");
0025   noise = conf_.getParameter<bool>("Noise");
0026   zeroSuppression = conf_.getParameter<bool>("ZeroSuppression");
0027   theTOFCutForPeak = conf_.getParameter<double>("TOFCutForPeak");
0028   theTOFCutForDeconvolution = conf_.getParameter<double>("TOFCutForDeconvolution");
0029   cosmicShift = conf_.getUntrackedParameter<double>("CosmicDelayShift");
0030   inefficiency = conf_.getParameter<double>("Inefficiency");
0031   RealPedestals = conf_.getParameter<bool>("RealPedestals");
0032   SingleStripNoise = conf_.getParameter<bool>("SingleStripNoise");
0033   CommonModeNoise = conf_.getParameter<bool>("CommonModeNoise");
0034   BaselineShift = conf_.getParameter<bool>("BaselineShift");
0035   APVSaturationFromHIP = conf_.getParameter<bool>("APVSaturationFromHIP");
0036   APVSaturationProb = conf_.getParameter<double>("APVSaturationProb");
0037   cmnRMStib = conf_.getParameter<double>("cmnRMStib");
0038   cmnRMStob = conf_.getParameter<double>("cmnRMStob");
0039   cmnRMStid = conf_.getParameter<double>("cmnRMStid");
0040   cmnRMStec = conf_.getParameter<double>("cmnRMStec");
0041   pedOffset = (unsigned int)conf_.getParameter<double>("PedestalsOffset");
0042   PreMixing_ = conf_.getParameter<bool>("PreMixingMode");
0043   if (peakMode) {
0044     tofCut = theTOFCutForPeak;
0045     LogDebug("StripDigiInfo") << "APVs running in peak mode (poor time resolution)";
0046   } else {
0047     tofCut = theTOFCutForDeconvolution;
0048     LogDebug("StripDigiInfo") << "APVs running in deconvolution mode (good time resolution)";
0049   };
0050 
0051   theSiHitDigitizer = new SiHitDigitizer(conf_);
0052   theDigiSimLinkPileUpSignals = new DigiSimLinkPileUpSignals();
0053   theSiNoiseAdder = new SiGaussianTailNoiseAdder(theThreshold);
0054   theSiDigitalConverter = new SiTrivialDigitalConverter(theElectronPerADC, PreMixing_);
0055   theSiZeroSuppress = new SiStripFedZeroSuppression(theFedAlgo, nullptr);
0056 }
0057 
0058 DigiSimLinkAlgorithm::~DigiSimLinkAlgorithm() {
0059   delete theSiHitDigitizer;
0060   delete theDigiSimLinkPileUpSignals;
0061   delete theSiNoiseAdder;
0062   delete theSiDigitalConverter;
0063   delete theSiZeroSuppress;
0064   //delete particle;
0065   //delete pdt;
0066 }
0067 
0068 //  Run the algorithm for a given module
0069 //  ------------------------------------
0070 
0071 void DigiSimLinkAlgorithm::run(edm::DetSet<SiStripDigi>& outdigi,
0072                                edm::DetSet<SiStripRawDigi>& outrawdigi,
0073                                const std::vector<std::pair<const PSimHit*, int> >& input,
0074                                StripGeomDetUnit const* det,
0075                                GlobalVector bfield,
0076                                float langle,
0077                                edm::ESHandle<SiStripGain>& gainHandle,
0078                                edm::ESHandle<SiStripThreshold>& thresholdHandle,
0079                                edm::ESHandle<SiStripNoises>& noiseHandle,
0080                                edm::ESHandle<SiStripPedestals>& pedestalHandle,
0081                                edm::ESHandle<SiStripBadStrip>& deadChannelHandle,
0082                                const TrackerTopology* tTopo,
0083                                CLHEP::HepRandomEngine* engine) {
0084   theDigiSimLinkPileUpSignals->reset();
0085   unsigned int detID = det->geographicalId().rawId();
0086   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
0087   SiStripApvGain::Range detGainRange = gainHandle->getRange(detID);
0088   SiStripPedestals::Range detPedestalRange = pedestalHandle->getRange(detID);
0089   SiStripBadStrip::Range detBadStripRange = deadChannelHandle->getRange(detID);
0090   numStrips = (det->specificTopology()).nstrips();
0091 
0092   //stroing the bad stip of the the module. the module is not removed but just signal put to 0
0093   std::vector<bool> badChannels;
0094   badChannels.clear();
0095   badChannels.insert(badChannels.begin(), numStrips, false);
0096   SiStripBadStrip::data fs;
0097   for (SiStripBadStrip::ContainerIterator it = detBadStripRange.first; it != detBadStripRange.second; ++it) {
0098     fs = deadChannelHandle->decode(*it);
0099     for (int strip = fs.firstStrip; strip < fs.firstStrip + fs.range; ++strip)
0100       badChannels[strip] = true;
0101   }
0102 
0103   // local amplitude of detector channels (from processed PSimHit)
0104   //  locAmpl.clear();
0105   detAmpl.clear();
0106   //  locAmpl.insert(locAmpl.begin(),numStrips,0.);
0107   // total amplitude of detector channels
0108   detAmpl.insert(detAmpl.begin(), numStrips, 0.);
0109 
0110   firstChannelWithSignal = numStrips;
0111   lastChannelWithSignal = 0;
0112 
0113   // First: loop on the SimHits
0114   std::vector<std::pair<const PSimHit*, int> >::const_iterator simHitIter = input.begin();
0115   std::vector<std::pair<const PSimHit*, int> >::const_iterator simHitIterEnd = input.end();
0116   if (CLHEP::RandFlat::shoot(engine) > inefficiency) {
0117     for (; simHitIter != simHitIterEnd; ++simHitIter) {
0118       locAmpl.clear();
0119       locAmpl.insert(locAmpl.begin(), numStrips, 0.);
0120       // check TOF
0121       if (std::fabs(((*simHitIter).first)->tof() - cosmicShift -
0122                     det->surface().toGlobal(((*simHitIter).first)->localPosition()).mag() / 30.) < tofCut &&
0123           ((*simHitIter).first)->energyLoss() > 0) {
0124         localFirstChannel = numStrips;
0125         localLastChannel = 0;
0126         // process the hit
0127         theSiHitDigitizer->processHit(
0128             ((*simHitIter).first), *det, bfield, langle, locAmpl, localFirstChannel, localLastChannel, tTopo, engine);
0129 
0130         //APV Killer to simulate HIP effect
0131         //------------------------------------------------------
0132 
0133         if (APVSaturationFromHIP && !zeroSuppression) {
0134           int pdg_id = ((*simHitIter).first)->particleType();
0135           particle = pdt->particle(pdg_id);
0136           if (particle != nullptr) {
0137             float charge = particle->charge();
0138             bool isHadron = particle->isHadron();
0139             if (charge != 0 && isHadron) {
0140               if (CLHEP::RandFlat::shoot(engine) < APVSaturationProb) {
0141                 int FirstAPV = localFirstChannel / 128;
0142                 int LastAPV = localLastChannel / 128;
0143                 std::cout << "-------------------HIP--------------" << std::endl;
0144                 std::cout << "Killing APVs " << FirstAPV << " - " << LastAPV << " " << detID << std::endl;
0145                 for (int strip = FirstAPV * 128; strip < LastAPV * 128 + 128; ++strip)
0146                   badChannels[strip] = true;  //doing like that I remove the signal information only after the
0147                                               //stip that got the HIP but it remains the signal of the previous
0148                                               //one. I'll make a further loop to remove all signal
0149               }
0150             }
0151           }
0152         }
0153 
0154         theDigiSimLinkPileUpSignals->add(
0155             locAmpl, localFirstChannel, localLastChannel, ((*simHitIter).first), (*simHitIter).second);
0156 
0157         // sum signal on strips
0158         for (size_t iChannel = localFirstChannel; iChannel < localLastChannel; iChannel++) {
0159           if (locAmpl[iChannel] > 0.) {
0160             //if(!badChannels[iChannel]) detAmpl[iChannel]+=locAmpl[iChannel];
0161             //locAmpl[iChannel]=0;
0162             detAmpl[iChannel] += locAmpl[iChannel];
0163           }
0164         }
0165         if (firstChannelWithSignal > localFirstChannel)
0166           firstChannelWithSignal = localFirstChannel;
0167         if (lastChannelWithSignal < localLastChannel)
0168           lastChannelWithSignal = localLastChannel;
0169       }
0170     }
0171   }
0172 
0173   //removing signal from the dead (and HIP effected) strips
0174   for (int strip = 0; strip < numStrips; ++strip)
0175     if (badChannels[strip])
0176       detAmpl[strip] = 0;
0177 
0178   const DigiSimLinkPileUpSignals::HitToDigisMapType& theLink(theDigiSimLinkPileUpSignals->dumpLink());
0179   const DigiSimLinkPileUpSignals::HitCounterToDigisMapType& theCounterLink(
0180       theDigiSimLinkPileUpSignals->dumpCounterLink());
0181 
0182   if (zeroSuppression) {
0183     if (noise) {
0184       int RefStrip = int(numStrips / 2.);
0185       while (RefStrip < numStrips &&
0186              badChannels[RefStrip]) {  //if the refstrip is bad, I move up to when I don't find it
0187         RefStrip++;
0188       }
0189       if (RefStrip < numStrips) {
0190         float noiseRMS = noiseHandle->getNoise(RefStrip, detNoiseRange);
0191         float gainValue = gainHandle->getStripGain(RefStrip, detGainRange);
0192         theSiNoiseAdder->addNoise(detAmpl,
0193                                   firstChannelWithSignal,
0194                                   lastChannelWithSignal,
0195                                   numStrips,
0196                                   noiseRMS * theElectronPerADC / gainValue,
0197                                   engine);
0198       }
0199     }
0200     digis.clear();
0201     theSiZeroSuppress->suppress(theSiDigitalConverter->convert(detAmpl, gainHandle.product(), detID),
0202                                 digis,
0203                                 detID,
0204                                 *noiseHandle,
0205                                 *thresholdHandle);
0206     push_link(digis, theLink, theCounterLink, detAmpl, detID);
0207     outdigi.data = digis;
0208   }
0209 
0210   if (!zeroSuppression) {
0211     //if(noise){
0212     // the constant pedestal offset is needed because
0213     //   negative adc counts are not allowed in case
0214     //   Pedestal and CMN subtraction is performed.
0215     //   The pedestal value read from the conditions
0216     //   is pedValue and after the pedestal subtraction
0217     //   the baseline is zero. The Common Mode Noise
0218     //   is not subtracted from the negative adc counts
0219     //   channels. Adding pedOffset the baseline is set
0220     //   to pedOffset after pedestal subtraction and CMN
0221     //   is subtracted to all the channels since none of
0222     //   them has negative adc value. The pedOffset is
0223     //   treated as a constant component in the CMN
0224     //   estimation and subtracted as CMN.
0225 
0226     //calculating the charge deposited on each APV and subtracting the shift
0227     //------------------------------------------------------
0228     if (BaselineShift) {
0229       theSiNoiseAdder->addBaselineShift(detAmpl, badChannels);
0230     }
0231 
0232     //Adding the strip noise
0233     //------------------------------------------------------
0234     if (noise) {
0235       std::vector<float> noiseRMSv;
0236       noiseRMSv.clear();
0237       noiseRMSv.insert(noiseRMSv.begin(), numStrips, 0.);
0238 
0239       if (SingleStripNoise) {
0240         for (int strip = 0; strip < numStrips; ++strip) {
0241           if (!badChannels[strip])
0242             noiseRMSv[strip] = (noiseHandle->getNoise(strip, detNoiseRange)) * theElectronPerADC;
0243         }
0244 
0245       } else {
0246         int RefStrip = 0;  //int(numStrips/2.);
0247         while (RefStrip < numStrips &&
0248                badChannels[RefStrip]) {  //if the refstrip is bad, I move up to when I don't find it
0249           RefStrip++;
0250         }
0251         if (RefStrip < numStrips) {
0252           float noiseRMS = noiseHandle->getNoise(RefStrip, detNoiseRange) * theElectronPerADC;
0253           for (int strip = 0; strip < numStrips; ++strip) {
0254             if (!badChannels[strip])
0255               noiseRMSv[strip] = noiseRMS;
0256           }
0257         }
0258       }
0259 
0260       theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv, engine);
0261     }
0262 
0263     //adding the CMN
0264     //------------------------------------------------------
0265     if (CommonModeNoise) {
0266       float cmnRMS = 0.;
0267       DetId detId(detID);
0268       uint32_t SubDet = detId.subdetId();
0269       if (SubDet == 3) {
0270         cmnRMS = cmnRMStib;
0271       } else if (SubDet == 4) {
0272         cmnRMS = cmnRMStid;
0273       } else if (SubDet == 5) {
0274         cmnRMS = cmnRMStob;
0275       } else if (SubDet == 6) {
0276         cmnRMS = cmnRMStec;
0277       }
0278       cmnRMS *= theElectronPerADC;
0279       theSiNoiseAdder->addCMNoise(detAmpl, cmnRMS, badChannels, engine);
0280     }
0281 
0282     //Adding the pedestals
0283     //------------------------------------------------------
0284 
0285     std::vector<float> vPeds;
0286     vPeds.clear();
0287     vPeds.insert(vPeds.begin(), numStrips, 0.);
0288 
0289     if (RealPedestals) {
0290       for (int strip = 0; strip < numStrips; ++strip) {
0291         if (!badChannels[strip])
0292           vPeds[strip] = (pedestalHandle->getPed(strip, detPedestalRange) + pedOffset) * theElectronPerADC;
0293       }
0294     } else {
0295       for (int strip = 0; strip < numStrips; ++strip) {
0296         if (!badChannels[strip])
0297           vPeds[strip] = pedOffset * theElectronPerADC;
0298       }
0299     }
0300 
0301     theSiNoiseAdder->addPedestals(detAmpl, vPeds);
0302 
0303     //if(!RealPedestals&&!CommonModeNoise&&!noise&&!BaselineShift&&!APVSaturationFromHIP){
0304     //  edm::LogWarning("DigiSimLinkAlgorithm")<<"You are running the digitizer without Noise generation and without applying Zero Suppression. ARE YOU SURE???";
0305     //}else{
0306 
0307     rawdigis.clear();
0308     rawdigis = theSiDigitalConverter->convertRaw(detAmpl, gainHandle.product(), detID);
0309     push_link_raw(rawdigis, theLink, theCounterLink, detAmpl, detID);
0310     outrawdigi.data = rawdigis;
0311 
0312     //}
0313   }
0314 }
0315 
0316 void DigiSimLinkAlgorithm::push_link(const DigitalVecType& digis,
0317                                      const HitToDigisMapType& htd,
0318                                      const HitCounterToDigisMapType& hctd,
0319                                      const std::vector<float>& afterNoise,
0320                                      unsigned int detID) {
0321   link_coll.clear();
0322   for (DigitalVecType::const_iterator i = digis.begin(); i != digis.end(); i++) {
0323     // Instead of checking the validity of the links against the digis,
0324     //  let's loop over digis and push the corresponding link
0325     HitToDigisMapType::const_iterator mi(htd.find(i->strip()));
0326     if (mi == htd.end())
0327       continue;
0328     HitCounterToDigisMapType::const_iterator cmi(hctd.find(i->strip()));
0329     std::map<const PSimHit*, Amplitude> totalAmplitudePerSimHit;
0330     for (std::vector<std::pair<const PSimHit*, Amplitude> >::const_iterator simul = (*mi).second.begin();
0331          simul != (*mi).second.end();
0332          simul++) {
0333       totalAmplitudePerSimHit[(*simul).first] += (*simul).second;
0334     }
0335 
0336     //--- include the noise as well
0337     double totalAmplitude1 = afterNoise[(*mi).first];
0338 
0339     //--- digisimlink
0340     int sim_counter = 0;
0341     for (std::map<const PSimHit*, Amplitude>::const_iterator iter = totalAmplitudePerSimHit.begin();
0342          iter != totalAmplitudePerSimHit.end();
0343          iter++) {
0344       float threshold = 0.;
0345       float fraction = (*iter).second / totalAmplitude1;
0346       if (fraction >= threshold) {
0347         // Noise fluctuation could make fraction>1. Unphysical, set it by hand = 1.
0348         if (fraction > 1.)
0349           fraction = 1.;
0350         for (std::vector<std::pair<const PSimHit*, int> >::const_iterator simcount = (*cmi).second.begin();
0351              simcount != (*cmi).second.end();
0352              ++simcount) {
0353           if ((*iter).first == (*simcount).first)
0354             sim_counter = (*simcount).second;
0355         }
0356         link_coll.push_back(StripDigiSimLink((*mi).first,                 //channel
0357                                              ((*iter).first)->trackId(),  //simhit trackId
0358                                              sim_counter,                 //simhit counter
0359                                              ((*iter).first)->eventId(),  //simhit eventId
0360                                              fraction));                  //fraction
0361       }
0362     }
0363   }
0364 }
0365 
0366 void DigiSimLinkAlgorithm::push_link_raw(const DigitalRawVecType& digis,
0367                                          const HitToDigisMapType& htd,
0368                                          const HitCounterToDigisMapType& hctd,
0369                                          const std::vector<float>& afterNoise,
0370                                          unsigned int detID) {
0371   link_coll.clear();
0372   int nstrip = -1;
0373   for (DigitalRawVecType::const_iterator i = digis.begin(); i != digis.end(); i++) {
0374     nstrip++;
0375     // Instead of checking the validity of the links against the digis,
0376     //  let's loop over digis and push the corresponding link
0377     HitToDigisMapType::const_iterator mi(htd.find(nstrip));
0378     HitCounterToDigisMapType::const_iterator cmi(hctd.find(nstrip));
0379     if (mi == htd.end())
0380       continue;
0381     std::map<const PSimHit*, Amplitude> totalAmplitudePerSimHit;
0382     for (std::vector<std::pair<const PSimHit*, Amplitude> >::const_iterator simul = (*mi).second.begin();
0383          simul != (*mi).second.end();
0384          simul++) {
0385       totalAmplitudePerSimHit[(*simul).first] += (*simul).second;
0386     }
0387 
0388     //--- include the noise as well
0389     double totalAmplitude1 = afterNoise[(*mi).first];
0390 
0391     //--- digisimlink
0392     int sim_counter_raw = 0;
0393     for (std::map<const PSimHit*, Amplitude>::const_iterator iter = totalAmplitudePerSimHit.begin();
0394          iter != totalAmplitudePerSimHit.end();
0395          iter++) {
0396       float threshold = 0.;
0397       float fraction = (*iter).second / totalAmplitude1;
0398       if (fraction >= threshold) {
0399         //Noise fluctuation could make fraction>1. Unphysical, set it by hand.
0400         if (fraction > 1.)
0401           fraction = 1.;
0402         //add counter information
0403         for (std::vector<std::pair<const PSimHit*, int> >::const_iterator simcount = (*cmi).second.begin();
0404              simcount != (*cmi).second.end();
0405              ++simcount) {
0406           if ((*iter).first == (*simcount).first)
0407             sim_counter_raw = (*simcount).second;
0408         }
0409         link_coll.push_back(StripDigiSimLink((*mi).first,                 //channel
0410                                              ((*iter).first)->trackId(),  //simhit trackId
0411                                              sim_counter_raw,             //simhit counter
0412                                              ((*iter).first)->eventId(),  //simhit eventId
0413                                              fraction));                  //fraction
0414       }
0415     }
0416   }
0417 }
[ 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!