Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:32

0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalTimeSlewSim.h"
0002 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0003 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
0004 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0007 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include "CLHEP/Random/RandGaussQ.h"
0012 
0013 HcalTimeSlewSim::HcalTimeSlewSim(const CaloVSimParameterMap* parameterMap, double minFCToDelay)
0014     : theParameterMap(parameterMap), minFCToDelay_(minFCToDelay) {}
0015 
0016 // not quite adequate to 25ns high-PU regime
0017 double HcalTimeSlewSim::charge(const CaloSamples& samples) const {
0018   double totalCharge = 0.;
0019   for (int i = 0; i < 4; ++i) {
0020     int bin = i + samples.presamples();
0021     if (bin < samples.size()) {
0022       totalCharge += samples[bin];
0023     }
0024   }
0025   return totalCharge;
0026 }
0027 
0028 void HcalTimeSlewSim::delay(CaloSamples& cs,
0029                             CLHEP::HepRandomEngine* engine,
0030                             const HcalTimeSlew* hcalTimeSlew_delay) const {
0031   // HO goes slow, HF shouldn't be used at all
0032   //ZDC not used for the moment
0033 
0034   DetId detId(cs.id());
0035   if (detId.det() == DetId::Calo && detId.subdetId() == HcalZDCDetId::SubdetectorId)
0036     return;
0037   HcalDetId hcalDetId(detId);
0038 
0039   if (hcalDetId.subdet() == HcalBarrel || hcalDetId.subdet() == HcalEndcap || hcalDetId.subdet() == HcalOuter) {
0040     HcalTimeSlew::BiasSetting biasSetting =
0041         (hcalDetId.subdet() == HcalOuter) ? HcalTimeSlew::Slow : HcalTimeSlew::Medium;
0042 
0043     // double totalCharge = charge(cs); // old TS...
0044 
0045     int maxbin = cs.size();
0046     CaloSamples data(detId, maxbin);  // for a temporary copy
0047     data = cs;
0048 
0049     // smearing
0050     int soi = cs.presamples();
0051     double eps = 1.e-6;
0052     double scale_factor = 0.6;
0053     double scale = cs[soi] / scale_factor;
0054     double smearns = 0.;
0055     double cut = minFCToDelay_;  //5. fC (above mean) for signal to be delayed
0056 
0057     const HcalSimParameters& params = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(detId));
0058     if (params.doTimeSmear()) {
0059       double rms = params.timeSmearRMS(scale);
0060       smearns = CLHEP::RandGaussQ::shoot(engine) * rms;
0061       LogDebug("HcalTimeSlewSim") << "TimeSmear charge " << scale << " rms " << rms << " smearns " << smearns;
0062     }
0063 
0064     // cycle over TS',  it - current TS index
0065     for (int it = 0; it < maxbin - 1;) {
0066       double datai = cs[it];
0067       int nts = 0;
0068       double tshift = smearns;
0069       double totalCharge = datai / scale_factor;
0070 
0071       // until we get more precise/reliable QIE8 simulation...
0072       if (totalCharge <= 0.)
0073         totalCharge = eps;  // protecion against negaive v.
0074       tshift += hcalTimeSlew_delay->delay(totalCharge, biasSetting);
0075       if (tshift <= 0.)
0076         tshift = eps;
0077 
0078       if (cut > -999.) {  //preserve compatibility
0079         if ((datai > cut) && (it == 0 || (datai > cs[it - 1]))) {
0080           // number of TS affected by current move depends on the signal shape:
0081           // rising or peaking
0082           nts = 2;  // rising
0083           if (datai > cs[it + 1])
0084             nts = 3;  // peaking
0085 
0086           // 1 or 2 TS to move from here,
0087           // 2d or 3d TS gets leftovers to preserve the sum
0088           for (int j = it; j < it + nts && j < maxbin; ++j) {
0089             // snippet borrowed from  CaloSamples::offsetTime(offset)
0090             // CalibFormats/CaloObjects/src/CaloSamples.cc
0091             double t = j * 25. - tshift;
0092             int firstbin = floor(t / 25.);
0093             double f = t / 25. - firstbin;
0094             int nextbin = firstbin + 1;
0095             double v1 = (firstbin < 0 || firstbin >= maxbin) ? 0. : cs[firstbin];
0096             double v2 = (nextbin < 0 || nextbin >= maxbin) ? 0. : cs[nextbin];
0097 
0098             // Keeping the overal sum intact
0099             if (nts == 2) {
0100               if (j == it)
0101                 data[j] = v2 * f;
0102               else
0103                 data[j] = v1 * (1. - f) + v2;
0104             } else {  // nts = 3
0105               if (j == it)
0106                 data[j] = v2 * f;
0107               if (j == it + 1)
0108                 data[j] = v1 * (1. - f) + v2 * f;
0109               if (j == it + nts - 1)
0110                 data[j] = v1 * (1. - f) + v2;
0111             }
0112 
0113           }  // end of local move of TS', now update...
0114           cs = data;
0115 
0116         }  // end of rising edge or peak finding
0117       } else {
0118         double t = it * 25. - tshift;
0119         int firstbin = floor(t / 25.);
0120         double f = t / 25. - firstbin;
0121         int nextbin = firstbin + 1;
0122         double v2 = (nextbin < 0 || nextbin >= maxbin) ? 0. : data[nextbin];
0123         data[it] = v2 * f;
0124         data[it + 1] += (v2 - data[it]);
0125         cs = data;
0126       }
0127 
0128       if (nts < 3)
0129         it++;
0130       else
0131         it += 2;
0132     }
0133 
0134     // final update of the shifted TS array
0135     cs = data;
0136   }
0137 }