CSCDigitizer/src/CSCWireElectronicsSim.cc

0001 #include "DataFormats/CSCDigi/interface/CSCWireDigi.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CSCGeometry/interface/CSCChamberSpecs.h"
0004 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0005 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0006 #include "SimMuon/CSCDigitizer/src/CSCAnalogSignal.h"
0007 #include "SimMuon/CSCDigitizer/src/CSCWireElectronicsSim.h"
0008
0009 #include <CLHEP/Random/RandGaussQ.h>
0010 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0011 #include <CLHEP/Units/SystemOfUnits.h>
0012
0013 #include <iostream>
0014
0015 CSCWireElectronicsSim::CSCWireElectronicsSim(const edm::ParameterSet &p)
0016     : CSCBaseElectronicsSim(p), theFraction(0.5), theWireNoise(0.0), theWireThreshold(0.) {
0017   fillAmpResponse();
0018 }
0019
0020 void CSCWireElectronicsSim::initParameters() {
0021   nElements = theLayerGeometry->numberOfWireGroups();
0022   theWireNoise = theSpecs->wireNoise(theShapingTime) * CLHEP::e_SI * pow(10.0, 15);
0023   theWireThreshold = theWireNoise * 8;
0024 }
0025
0026 int CSCWireElectronicsSim::readoutElement(int element) const { return theLayerGeometry->wireGroup(element); }
0027
0028 void CSCWireElectronicsSim::fillDigis(CSCWireDigiCollection &digis, CLHEP::HepRandomEngine *engine) {
0029   if (theSignalMap.empty()) {
0030     return;
0031   }
0032
0033   // Loop over analog signals, run the fractional discriminator on each one,
0034   // and save the DIGI in the layer.
0035   for (CSCSignalMap::iterator mapI = theSignalMap.begin(), lastSignal = theSignalMap.end(); mapI != lastSignal;
0036        ++mapI) {
0037     int wireGroup = (*mapI).first;
0038     const CSCAnalogSignal &signal = (*mapI).second;
0039     LogTrace("CSCWireElectronicsSim") << "CSCWireElectronicsSim: dump of wire signal follows... " << signal;
0040     int signalSize = signal.getSize();
0041
0042     int timeWord = 0;  // and this will remain if too early or late (<bx-6 or >bx+9)
0043
0044     // the way we handle noise in this chamber is by randomly varying
0045     // the threshold
0046     float threshold = theWireThreshold;
0047     if (doNoise_) {
0048       threshold += CLHEP::RandGaussQ::shoot(engine) * theWireNoise;
0049     }
0050     for (int ibin = 0; ibin < signalSize; ++ibin) {
0051       if (signal.getBinValue(ibin) > threshold) {
0052         // jackpot.  Now define this signal as everything up until
0053         // the signal goes below zero.
0054         int lastbin = signalSize;
0055         int i;
0056         for (i = ibin; i < signalSize; ++i) {
0057           if (signal.getBinValue(i) < 0.) {
0058             lastbin = i;
0059             break;
0060           }
0061         }
0062
0063         float qMax = 0.0;
0064         // in this loop, find the max charge and the 'fifth' electron arrival
0065         for (i = ibin; i < lastbin; ++i) {
0066           float next_charge = signal.getBinValue(i);
0067           if (next_charge > qMax) {
0068             qMax = next_charge;
0069           }
0070         }
0071
0072         int bin_firing_FD = 0;
0073         for (i = ibin; i < lastbin; ++i) {
0074           if (signal.getBinValue(i) >= qMax * theFraction) {
0075             bin_firing_FD = i;
0076             //@@ Long-standing but unlikely minor bug, I (Tim) think - following
0077             //'break' was missing...
0078             //@@ ... So if both ibins 0 and 1 could fire FD, we'd flag the
0079             // firing bin as 1 not 0
0080             //@@ (since the above test was restricted to bin_firing_FD==0 too).
0081             break;
0082           }
0083         }
0084
0085         float tofOffset = timeOfFlightCalibration(wireGroup);
0086         int chamberType = theSpecs->chamberType();
0087
0088         // Note that CSCAnalogSignal::superimpose does not reset theTimeOffset
0089         // to the earliest of the two signal's time offsets. If it did then we
0090         // could handle signals from any time range e.g. form pileup events many
0091         // bx's from the signal bx (bx=0). But then we would be wastefully
0092         // storing signals over times which we can never see in the real
0093         // detector, because only hits within a few bx's of bx=0 are read out.
0094         // Instead, the working time range for wire hits is always started from
0095         // theSignalStartTime, set as a parameter in the config file.
0096         // On the other hand, if any of the overlapped CSCAnalogSignals happens
0097         // to have a timeOffset earlier than theSignalStartTime (which is
0098         // currently set to -100 ns) then we're in trouble. For pileup events
0099         // this would mean events from collisions earlier than 4 bx before the
0100         // signal bx.
0101
0102         float fdTime = theSignalStartTime + theSamplingTime * bin_firing_FD;
0103         if (doNoise_) {
0104           fdTime += theTimingCalibrationError[chamberType] * CLHEP::RandGaussQ::shoot(engine);
0105         }
0106
0107         float bxFloat = (fdTime - tofOffset - theBunchTimingOffsets[chamberType]) / theBunchSpacing + theOffsetOfBxZero;
0108         int bxInt = static_cast<int>(bxFloat);
0109         if (bxFloat >= 0 && bxFloat < 16) {
0110           timeWord |= (1 << bxInt);
0111           // discriminator stays high for 35 ns
0112           if (bxFloat - bxInt > 0.6) {
0113             timeWord |= (1 << (bxInt + 1));
0114           }
0115         }
0116
0117         // Wire digi as of Oct-2006 adapted to real data: time word has 16 bits
0118         // with set bit flagging appropriate bunch crossing, and bx 0
0119         // corresponding to the 7th bit, 'bit 6':
0120
0121         //      1st bit set (bit 0) <-> bx -6
0122         //      2nd              1  <-> bx -5
0123         //      ...           ...       ....
0124         //      7th              6  <-> bx  0
0125         //      8th              7  <-> bx +1
0126         //      ...           ...       ....
0127         //     16th             15  <-> bx +9
0128
0129         // skip over all the time bins used for this digi
0130         ibin = lastbin;
0131       }  // if over threshold
0132     }  // loop over time bins in signal
0133
0134     // Only create a wire digi if there is a wire hit within [-6 bx, +9 bx]
0135     if (timeWord != 0) {
0136       CSCWireDigi newDigi(wireGroup, timeWord);
0137       LogTrace("CSCWireElectronicsSim") << "CSCWireElectronicsSim: " << newDigi;
0138       digis.insertDigi(layerId(), newDigi);
0139       addLinks(channelIndex(wireGroup));
0140     }
0141   }  // loop over wire signals
0142 }
0143
0144 float CSCWireElectronicsSim::calculateAmpResponse(float t) const {
0145   static const float fC_by_ns = 1000000;
0146   static const float resistor = 20000;
0147   static const float amplifier_pole = 1 / 7.5;
0148   static const float fastest_chamber_exp_risetime = 10.;
0149   static const float p0 = amplifier_pole;
0150   static const float p1 = 1 / fastest_chamber_exp_risetime;
0151
0152   static const float dp = p0 - p1;
0153
0154   // ENABLE DISC:
0155
0156   static const float norm = -12 * resistor * p1 * pow(p0 / dp, 4) / fC_by_ns;
0157
0158   float enable_disc_volts =
0159       norm * (exp(-p0 * t) * (1 + t * dp + pow(t * dp, 2) / 2 + pow(t * dp, 3) / 6) - exp(-p1 * t));
0160   static const float collectionFraction = 0.12;
0161   static const float igain = 1. / 0.005;  // volts per fC
0162   return enable_disc_volts * igain * collectionFraction;
0163 }
0164
0165 float CSCWireElectronicsSim::timeOfFlightCalibration(int wireGroup) const {
0166   // calibration is done for groups of 8 wire groups, facetiously
0167   // called wireGroupGroups
0168   int middleWireGroup = wireGroup - wireGroup % 8 + 4;
0169   int numberOfWireGroups = theLayerGeometry->numberOfWireGroups();
0170   if (middleWireGroup > numberOfWireGroups)
0171     middleWireGroup = numberOfWireGroups;
0172
0173   GlobalPoint centerOfGroupGroup = theLayer->centerOfWireGroup(middleWireGroup);
0174   float averageDist = centerOfGroupGroup.mag();
0175   float averageTOF = averageDist * CLHEP::cm / c_light;  // Units of c_light: mm/ns
0176
0177   LogTrace("CSCWireElectronicsSim") << "CSCWireElectronicsSim: TofCalib  wg = " << wireGroup
0178                                     << " mid wg = " << middleWireGroup << " av dist = " << averageDist
0179                                     << " av tof = " << averageTOF;
0180
0181   return averageTOF;
0182 }