Reference Manual

Go to the documentation of this file.
00001 // File: SiStripDigitizerAlgorithm.cc
00002 // Description:  Steering class for digitization.
00003 
00004 #include <vector>
00005 #include <algorithm>
00006 #include <iostream>
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00010 #include "SiStripDigitizerAlgorithm.h"
00011 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00012 #include "DataFormats/Common/interface/DetSetVector.h"
00013 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
00014 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00015 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
00016 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
00017 #include "CondFormats/DataRecord/interface/SiStripCondDataRecords.h"
00018 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
00019 #include "CondFormats/SiStripObjects/interface/SiStripThreshold.h"
00020 #include "CondFormats/SiStripObjects/interface/SiStripPedestals.h"
00021 #include "CondFormats/SiStripObjects/interface/SiStripBadStrip.h"
00022 #include "CLHEP/Random/RandFlat.h"
00023 
00024 SiStripDigitizerAlgorithm::SiStripDigitizerAlgorithm(const edm::ParameterSet& conf, CLHEP::HepRandomEngine& eng):
00025   lorentzAngleName(conf.getParameter<std::string>("LorentzAngle")),
00026   theThreshold(conf.getParameter<double>("NoiseSigmaThreshold")),
00027   cmnRMStib(conf.getParameter<double>("cmnRMStib")),
00028   cmnRMStob(conf.getParameter<double>("cmnRMStob")),
00029   cmnRMStid(conf.getParameter<double>("cmnRMStid")),
00030   cmnRMStec(conf.getParameter<double>("cmnRMStec")),
00031   APVSaturationProb(conf.getParameter<double>("APVSaturationProb")),
00032   makeDigiSimLinks_(conf.getUntrackedParameter<bool>("makeDigiSimLinks", false)),
00033   peakMode(conf.getParameter<bool>("APVpeakmode")),
00034   noise(conf.getParameter<bool>("Noise")),
00035   RealPedestals(conf.getParameter<bool>("RealPedestals")), 
00036   SingleStripNoise(conf.getParameter<bool>("SingleStripNoise")),
00037   CommonModeNoise(conf.getParameter<bool>("CommonModeNoise")),
00038   BaselineShift(conf.getParameter<bool>("BaselineShift")),
00039   APVSaturationFromHIP(conf.getParameter<bool>("APVSaturationFromHIP")),
00040   theFedAlgo(conf.getParameter<int>("FedAlgorithm")),
00041   zeroSuppression(conf.getParameter<bool>("ZeroSuppression")),
00042   theElectronPerADC(conf.getParameter<double>( peakMode ? "electronPerAdcPeak" : "electronPerAdcDec" )),
00043   theTOFCutForPeak(conf.getParameter<double>("TOFCutForPeak")),
00044   theTOFCutForDeconvolution(conf.getParameter<double>("TOFCutForDeconvolution")),
00045   tofCut(peakMode ? theTOFCutForPeak : theTOFCutForDeconvolution),
00046   cosmicShift(conf.getUntrackedParameter<double>("CosmicDelayShift")),
00047   inefficiency(conf.getParameter<double>("Inefficiency")),
00048   pedOffset((unsigned int)conf.getParameter<double>("PedestalsOffset")),
00049   theSiHitDigitizer(new SiHitDigitizer(conf, eng)),
00050   theSiPileUpSignals(new SiPileUpSignals()),
00051   theSiNoiseAdder(new SiGaussianTailNoiseAdder(theThreshold, eng)),
00052   theSiDigitalConverter(new SiTrivialDigitalConverter(theElectronPerADC)),
00053   theSiZeroSuppress(new SiStripFedZeroSuppression(theFedAlgo)),
00054   theFlatDistribution(new CLHEP::RandFlat(eng, 0., 1.)) {
00055 
00056   if (peakMode) {
00057     LogDebug("StripDigiInfo")<<"APVs running in peak mode (poor time resolution)";
00058   } else {
00059     LogDebug("StripDigiInfo")<<"APVs running in deconvolution mode (good time resolution)";
00060   };
00061 }
00062 
00063 SiStripDigitizerAlgorithm::~SiStripDigitizerAlgorithm(){
00064 }
00065 
00066 void
00067 SiStripDigitizerAlgorithm::initializeDetUnit(StripGeomDetUnit* det, const edm::EventSetup& iSetup) {
00068   edm::ESHandle<SiStripBadStrip> deadChannelHandle;
00069   iSetup.get<SiStripBadChannelRcd>().get(deadChannelHandle);
00070 
00071   unsigned int detId = det->geographicalId().rawId();
00072   int numStrips = (det->specificTopology()).nstrips();  
00073 
00074   SiStripBadStrip::Range detBadStripRange = deadChannelHandle->getRange(detId);
00075   //storing the bad strip of the the module. the module is not removed but just signal put to 0
00076   std::vector<bool>& badChannels = allBadChannels[detId];
00077   badChannels.clear();
00078   badChannels.insert(badChannels.begin(), numStrips, false);
00079   for(SiStripBadStrip::ContainerIterator it = detBadStripRange.first; it != detBadStripRange.second; ++it) {
00080     SiStripBadStrip::data fs = deadChannelHandle->decode(*it);
00081     for(int strip = fs.firstStrip; strip < fs.firstStrip + fs.range; ++strip) badChannels[strip] = true;
00082   }
00083   firstChannelsWithSignal[detId] = numStrips;
00084   lastChannelsWithSignal[detId]= 0;
00085 }
00086 
00087 void
00088 SiStripDigitizerAlgorithm::initializeEvent(const edm::EventSetup& iSetup) {
00089   theSiPileUpSignals->reset();
00090 
00091   //get gain noise pedestal lorentzAngle from ES handle
00092   edm::ESHandle<ParticleDataTable> pdt;
00093   iSetup.getData(pdt);
00094   setParticleDataTable(&*pdt);
00095   iSetup.get<SiStripLorentzAngleSimRcd>().get(lorentzAngleName,lorentzAngleHandle);
00096 }
00097 
00098 //  Run the algorithm for a given module
00099 //  ------------------------------------
00100 
00101 void
00102 SiStripDigitizerAlgorithm::accumulateSimHits(std::vector<PSimHit>::const_iterator inputBegin,
00103                                              std::vector<PSimHit>::const_iterator inputEnd,
00104                                              const StripGeomDetUnit* det,
00105                                              const GlobalVector& bfield,
00106                                              const TrackerTopology *tTopo) {
00107   // produce SignalPoints for all SimHits in detector
00108   unsigned int detID = det->geographicalId().rawId();
00109   int numStrips = (det->specificTopology()).nstrips();  
00110 
00111   std::vector<bool>& badChannels = allBadChannels[detID];
00112   size_t thisFirstChannelWithSignal = numStrips;
00113   size_t thisLastChannelWithSignal = numStrips;
00114 
00115   float langle = (lorentzAngleHandle.isValid()) ? lorentzAngleHandle->getLorentzAngle(detID) : 0.;
00116 
00117   std::vector<double> locAmpl(numStrips, 0.);
00118 
00119   // Loop over hits
00120 
00121   uint32_t detId = det->geographicalId().rawId();
00122   // First: loop on the SimHits
00123   if(theFlatDistribution->fire()>inefficiency) {
00124     for (std::vector<PSimHit>::const_iterator simHitIter = inputBegin; simHitIter != inputEnd; ++simHitIter) {
00125       // skip hits not in this detector.
00126       if((*simHitIter).detUnitId() != detId) {
00127         continue;
00128       }
00129       // check TOF
00130       if (std::fabs(simHitIter->tof() - cosmicShift - det->surface().toGlobal(simHitIter->localPosition()).mag()/30.) < tofCut && simHitIter->energyLoss()>0) {
00131         size_t localFirstChannel = numStrips;
00132         size_t localLastChannel  = 0;
00133         // process the hit
00134         theSiHitDigitizer->processHit(&*simHitIter, *det, bfield, langle, locAmpl, localFirstChannel, localLastChannel, tTopo);
00135           
00136                   //APV Killer to simulate HIP effect
00137                   //------------------------------------------------------
00138                   
00139                   if(APVSaturationFromHIP&&!zeroSuppression){
00140                     int pdg_id = simHitIter->particleType();
00141                         particle = pdt->particle(pdg_id);
00142                         if(particle != NULL){
00143                                 float charge = particle->charge();
00144                                 bool isHadron = particle->isHadron();
00145                             if(charge!=0 && isHadron){
00146                                         if(theFlatDistribution->fire()<APVSaturationProb){
00147                                                 int FirstAPV = localFirstChannel/128;
00148                                                 int LastAPV = localLastChannel/128;
00149                                                 std::cout << "-------------------HIP--------------" << std::endl;
00150                                                 std::cout << "Killing APVs " << FirstAPV << " - " <<LastAPV << " " << detID <<std::endl;
00151                                                 for(int strip = FirstAPV*128; strip < LastAPV*128 +128; ++strip) {
00152                                                         badChannels[strip] = true;
00153                                                 }
00154                                                 //doing like that I remove the signal information only after the 
00155                                                 //stip that got the HIP but it remains the signal of the previous
00156                                                 //one. I'll make a further loop to remove all signal
00157                                         }
00158                                 }
00159                         }
00160               }             
00161                 
00162     
00163         if(thisFirstChannelWithSignal > localFirstChannel) thisFirstChannelWithSignal = localFirstChannel;
00164         if(thisLastChannelWithSignal < localLastChannel) thisLastChannelWithSignal = localLastChannel;
00165       }
00166     } // end for
00167   }
00168   theSiPileUpSignals->add(detID, locAmpl, thisFirstChannelWithSignal, thisLastChannelWithSignal);
00169 
00170   if(firstChannelsWithSignal[detID] > thisFirstChannelWithSignal) firstChannelsWithSignal[detID] = thisFirstChannelWithSignal;
00171   if(lastChannelsWithSignal[detID] < thisLastChannelWithSignal) lastChannelsWithSignal[detID] = thisLastChannelWithSignal;
00172 }
00173 
00174 void
00175 SiStripDigitizerAlgorithm::digitize(
00176                            edm::DetSet<SiStripDigi>& outdigi,
00177                            edm::DetSet<SiStripRawDigi>& outrawdigi,
00178                            const StripGeomDetUnit *det,
00179                            edm::ESHandle<SiStripGain> & gainHandle,
00180                            edm::ESHandle<SiStripThreshold> & thresholdHandle,
00181                            edm::ESHandle<SiStripNoises> & noiseHandle,
00182                            edm::ESHandle<SiStripPedestals> & pedestalHandle) {
00183   unsigned int detID = det->geographicalId().rawId();
00184   int numStrips = (det->specificTopology()).nstrips();  
00185 
00186   const SiPileUpSignals::SignalMapType* theSignal(theSiPileUpSignals->getSignal(detID));  
00187 
00188   std::vector<double> detAmpl(numStrips, 0.);
00189   if(theSignal) {
00190     for(const auto& amp : *theSignal) {
00191       detAmpl[amp.first] = amp.second;
00192     }
00193   }
00194 
00195   //removing signal from the dead (and HIP effected) strips
00196   std::vector<bool>& badChannels = allBadChannels[detID];
00197   for(int strip =0; strip < numStrips; ++strip) if(badChannels[strip]) detAmpl[strip] = 0.;
00198 
00199   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
00200   SiStripApvGain::Range detGainRange = gainHandle->getRange(detID);
00201   SiStripPedestals::Range detPedestalRange = pedestalHandle->getRange(detID);
00202 
00203 // -----------------------------------------------------------
00204 
00205   auto& firstChannelWithSignal = firstChannelsWithSignal[detID];
00206   auto& lastChannelWithSignal = lastChannelsWithSignal[detID];
00207 
00208   if(zeroSuppression){
00209     if(noise){
00210           int RefStrip = int(numStrips/2.);
00211           while(RefStrip<numStrips&&badChannels[RefStrip]){ //if the refstrip is bad, I move up to when I don't find it
00212                 RefStrip++;
00213           }
00214           if(RefStrip<numStrips){
00215                 float noiseRMS = noiseHandle->getNoise(RefStrip,detNoiseRange);
00216                 float gainValue = gainHandle->getStripGain(RefStrip, detGainRange);
00217                 theSiNoiseAdder->addNoise(detAmpl,firstChannelWithSignal,lastChannelWithSignal,numStrips,noiseRMS*theElectronPerADC/gainValue);
00218           }
00219         }
00220     DigitalVecType digis;
00221     theSiZeroSuppress->suppress(theSiDigitalConverter->convert(detAmpl, gainHandle, detID), digis, detID,noiseHandle,thresholdHandle);
00222     outdigi.data = digis;
00223   }
00224   
00225   
00226   if(!zeroSuppression){
00227     //if(noise){
00228       // the constant pedestal offset is needed because
00229       //   negative adc counts are not allowed in case
00230       //   Pedestal and CMN subtraction is performed.
00231       //   The pedestal value read from the conditions
00232       //   is pedValue and after the pedestal subtraction
00233       //   the baseline is zero. The Common Mode Noise
00234       //   is not subtracted from the negative adc counts
00235       //   channels. Adding pedOffset the baseline is set
00236       //   to pedOffset after pedestal subtraction and CMN
00237       //   is subtracted to all the channels since none of
00238       //   them has negative adc value. The pedOffset is
00239       //   treated as a constant component in the CMN
00240       //   estimation and subtracted as CMN.
00241       
00242          
00243                 //calculating the charge deposited on each APV and subtracting the shift
00244                 //------------------------------------------------------
00245                 if(BaselineShift){
00246                    theSiNoiseAdder->addBaselineShift(detAmpl, badChannels);
00247                 }
00248                 
00249                 //Adding the strip noise
00250                 //------------------------------------------------------                                                 
00251                 if(noise){
00252                     std::vector<float> noiseRMSv;
00253                         noiseRMSv.clear();
00254                     noiseRMSv.insert(noiseRMSv.begin(),numStrips,0.);
00255                         
00256                     if(SingleStripNoise){
00257                             for(int strip=0; strip< numStrips; ++strip){
00258                                         if(!badChannels[strip]) noiseRMSv[strip] = (noiseHandle->getNoise(strip,detNoiseRange))* theElectronPerADC;
00259                                 }
00260                         
00261                 } else {
00262                             int RefStrip = 0; //int(numStrips/2.);
00263                             while(RefStrip<numStrips&&badChannels[RefStrip]){ //if the refstrip is bad, I move up to when I don't find it
00264                                         RefStrip++;
00265                                 }
00266                                 if(RefStrip<numStrips){
00267                                         float noiseRMS = noiseHandle->getNoise(RefStrip,detNoiseRange) *theElectronPerADC;
00268                                         for(int strip=0; strip< numStrips; ++strip){
00269                                         if(!badChannels[strip]) noiseRMSv[strip] = noiseRMS;
00270                                         }
00271                                 }
00272                         }
00273                         
00274                         theSiNoiseAdder->addNoiseVR(detAmpl, noiseRMSv);
00275                 }                       
00276                 
00277                 //adding the CMN
00278                 //------------------------------------------------------
00279         if(CommonModeNoise){
00280                   float cmnRMS = 0.;
00281                   DetId  detId(detID);
00282                   uint32_t SubDet = detId.subdetId();
00283                   if(SubDet==3){
00284                     cmnRMS = cmnRMStib;
00285                   }else if(SubDet==4){
00286                     cmnRMS = cmnRMStid;
00287                   }else if(SubDet==5){
00288                     cmnRMS = cmnRMStob;
00289                   }else if(SubDet==6){
00290                     cmnRMS = cmnRMStec;
00291                   }
00292                   cmnRMS *= theElectronPerADC;
00293           theSiNoiseAdder->addCMNoise(detAmpl, cmnRMS, badChannels);
00294                 }
00295                 
00296                         
00297                 //Adding the pedestals
00298                 //------------------------------------------------------
00299                 
00300                 std::vector<float> vPeds;
00301                 vPeds.clear();
00302                 vPeds.insert(vPeds.begin(),numStrips,0.);
00303                 
00304                 if(RealPedestals){
00305                     for(int strip=0; strip< numStrips; ++strip){
00306                            if(!badChannels[strip]) vPeds[strip] = (pedestalHandle->getPed(strip,detPedestalRange)+pedOffset)* theElectronPerADC;
00307                     }
00308         } else {
00309                     for(int strip=0; strip< numStrips; ++strip){
00310                           if(!badChannels[strip]) vPeds[strip] = pedOffset* theElectronPerADC;
00311                         }
00312                 }
00313                 
00314                 theSiNoiseAdder->addPedestals(detAmpl, vPeds);  
00315                 
00316                  
00317         //if(!RealPedestals&&!CommonModeNoise&&!noise&&!BaselineShift&&!APVSaturationFromHIP){
00318     //  edm::LogWarning("SiStripDigitizer")<<"You are running the digitizer without Noise generation and without applying Zero Suppression. ARE YOU SURE???";
00319     //}else{                                                     
00320     
00321     DigitalRawVecType rawdigis = theSiDigitalConverter->convertRaw(detAmpl, gainHandle, detID);
00322     outrawdigi.data = rawdigis;
00323         
00324         //}
00325   }
00326 }
/afs/cern.ch/work/c/cnuttens/private/SiStripDev/DocumentationProduction/CMSSW_6_1_0/src/SimTracker/SiStripDigitizer/plugins/SiStripDigitizerAlgorithm.cc
