SiStripAPVAnalysis/src/TT6PedestalCalculator.cc

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