SiStripAPVAnalysis/src/SimplePedestalCalculator.cc

0001 #include "CalibTracker/SiStripAPVAnalysis/interface/SimplePedestalCalculator.h"
0002
0003 #include <cmath>
0004 #include <numeric>
0005 #include <algorithm>
0006
0007 using namespace std;
0008 SimplePedestalCalculator::SimplePedestalCalculator(int evnt_ini) : numberOfEvents(0), alreadyUsedEvent(false) {
0009   if (false)
0010     cout << "Constructing SimplePedestalCalculator " << endl;
0011   eventsRequiredToCalibrate = evnt_ini;
0012   //  eventsRequiredToUpdate    = evnt_iter;
0013   //  cutToAvoidSignal          = sig_cut;
0014   init();
0015 }
0016 //
0017 // Initialization.
0018 //
0019 void SimplePedestalCalculator::init() {
0020   theRawNoise.clear();
0021   thePedestal.clear();
0022   thePedSum.clear();
0023   thePedSqSum.clear();
0024   theEventPerStrip.clear();
0025   theStatus.setCalibrating();
0026 }
0027 //
0028 //  -- Destructor
0029 //
0030 SimplePedestalCalculator::~SimplePedestalCalculator() {
0031   if (false)
0032     cout << "Destructing SimplePedestalCalculator " << endl;
0033 }
0034
0035 //
0036 // -- Set Pedestal Update Status
0037 //
0038 void SimplePedestalCalculator::updateStatus() {
0039   if (theStatus.isCalibrating() && numberOfEvents >= eventsRequiredToCalibrate) {
0040     theStatus.setUpdating();
0041   }
0042 }
0043
0044 //
0045 // -- Initialize or Update (when needed) Pedestal Values
0046 //
0047 void SimplePedestalCalculator::updatePedestal(ApvAnalysis::RawSignalType& in) {
0048   if (alreadyUsedEvent == false) {
0049     alreadyUsedEvent = true;
0050     numberOfEvents++;
0051     if (theStatus.isCalibrating()) {
0052       initializePedestal(in);
0053     } else if (theStatus.isUpdating()) {
0054       refinePedestal(in);
0055     }
0056     updateStatus();
0057   }
0058 }
0059
0060 //
0061 // -- Initialize Pedestal Values using a set of events (eventsRequiredToCalibrate)
0062 //
0063 void SimplePedestalCalculator::initializePedestal(ApvAnalysis::RawSignalType& in) {
0064   if (numberOfEvents == 1) {
0065     thePedSum.clear();
0066     thePedSqSum.clear();
0067     theEventPerStrip.clear();
0068
0069     thePedSum.reserve(128);
0070     thePedSqSum.reserve(128);
0071     theEventPerStrip.reserve(128);
0072
0073     thePedSum.resize(in.data.size(), 0);
0074     thePedSqSum.resize(in.data.size(), 0);
0075     theEventPerStrip.resize(in.data.size(), 0);
0076   }
0077
0078   //eventsRequiredToCalibrate is considered the minimum number of events to be used
0079
0080   if (numberOfEvents <= eventsRequiredToCalibrate) {
0081     edm::DetSet<SiStripRawDigi>::const_iterator i = in.data.begin();
0082     int ii = 0;
0083     for (; i != in.data.end(); i++) {
0084       thePedSum[ii] += (*i).adc();
0085       thePedSqSum[ii] += ((*i).adc()) * ((*i).adc());
0086       theEventPerStrip[ii]++;
0087       ii++;
0088     }
0089   }
0090   if (numberOfEvents == eventsRequiredToCalibrate) {
0091     thePedestal.clear();
0092     theRawNoise.clear();
0093     edm::DetSet<SiStripRawDigi>::const_iterator i = in.data.begin();
0094     int ii = 0;
0095     for (; i != in.data.end(); i++) {
0096       // the pedestal is calculated as int, as required by FED.
0097       int avVal = (theEventPerStrip[ii]) ? thePedSum[ii] / theEventPerStrip[ii] : 0;
0098
0099       double sqAvVal = (theEventPerStrip[ii]) ? thePedSqSum[ii] / theEventPerStrip[ii] : 0.0;
0100       double corr_fac = (theEventPerStrip[ii] > 1) ? (theEventPerStrip[ii] / (theEventPerStrip[ii] - 1)) : 1.0;
0101       double rmsVal = (sqAvVal - avVal * avVal > 0.0) ? sqrt(corr_fac * (sqAvVal - avVal * avVal)) : 0.0;
0102       thePedestal.push_back(static_cast<float>(avVal));
0103       theRawNoise.push_back(static_cast<float>(rmsVal));
0104       ii++;
0105     }
0106   }
0107 }
0108
0109 //
0110 // -- Update Pedestal Values when needed.
0111 //
0112
0113 void SimplePedestalCalculator::refinePedestal(ApvAnalysis::RawSignalType& in) {
0114   // keep adding th adc count for any events
0115
0116   unsigned int ii = 0;
0117   ApvAnalysis::RawSignalType::const_iterator i = in.data.begin();
0118   for (; i < in.data.end(); i++) {
0119     thePedSum[ii] += (*i).adc();
0120     thePedSqSum[ii] += ((*i).adc()) * ((*i).adc());
0121     theEventPerStrip[ii]++;
0122
0123     ii++;
0124   }
0125
0126   // calculate a new pedestal any events, so it will come for free when for the last event
0127
0128   for (unsigned int iii = 0; iii < in.data.size(); iii++) {
0129     if (theEventPerStrip[iii] > 10) {
0130       int avVal = (theEventPerStrip[iii]) ? thePedSum[iii] / theEventPerStrip[iii] : 0;
0131
0132       double sqAvVal = (theEventPerStrip[iii]) ? thePedSqSum[iii] / theEventPerStrip[iii] : 0.0;
0133
0134       double rmsVal = (sqAvVal - avVal * avVal > 0.0) ? sqrt(sqAvVal - avVal * avVal) : 0.0;
0135
0136       if (avVal != 0) {
0137         thePedestal[iii] = static_cast<float>(avVal);
0138         theRawNoise[iii] = static_cast<float>(rmsVal);
0139       }
0140     }
0141   }
0142 }
0143
0144 //
0145 // Define New Event
0146 //
0147
0148 void SimplePedestalCalculator::newEvent() { alreadyUsedEvent = false; }