SiStripAPVAnalysis/src/TT6NoiseCalculator.cc

0001 #include "CalibTracker/SiStripAPVAnalysis/interface/TT6NoiseCalculator.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include <cmath>
0004 #include <numeric>
0005 #include <algorithm>
0006 using namespace std;
0007 //
0008 //  Constructors
0009 //
0010 TT6NoiseCalculator::TT6NoiseCalculator() : numberOfEvents(0), alreadyUsedEvent(false) {
0011   if (false)
0012     cout << "Constructing TT6NoiseCalculator " << endl;
0013   init();
0014 }
0015 //
0016 TT6NoiseCalculator::TT6NoiseCalculator(int evnt_ini, int evnt_iter, float sig_cut)
0017     : numberOfEvents(0), alreadyUsedEvent(false) {
0018   if (false)
0019     cout << "Constructing TT6NoiseCalculator " << endl;
0020   eventsRequiredToCalibrate_ = evnt_ini;
0021   eventsRequiredToUpdate_ = evnt_iter;
0022   cutToAvoidSignal_ = sig_cut;
0023   init();
0024 }
0025 //
0026 // Initialization.
0027 //
0028 void TT6NoiseCalculator::init() {
0029   theCMPSubtractedSignal.clear();
0030   theNoise.clear();
0031   theNoiseSum.clear();
0032   theNoiseSqSum.clear();
0033   theEventPerStrip.clear();
0034   theStatus.setCalibrating();
0035 }
0036 //
0037 //  Destructor
0038 //
0039 TT6NoiseCalculator::~TT6NoiseCalculator() {
0040   if (false)
0041     cout << "Destructing TT6NoiseCalculator " << endl;
0042 }
0043 //
0044 // Update the Status of Noise Calculation
0045 //
0046 void TT6NoiseCalculator::updateStatus() {
0047   if (theStatus.isCalibrating() && numberOfEvents >= eventsRequiredToCalibrate_) {
0048     theStatus.setUpdating();
0049   }
0050 }
0051 //
0052 // Calculate and update (when needed) Noise Values
0053 //
0054 void TT6NoiseCalculator::updateNoise(ApvAnalysis::PedestalType& in) {
0055   if (alreadyUsedEvent == false) {
0056     alreadyUsedEvent = true;
0057     numberOfEvents++;
0058
0059     if (numberOfEvents == 1 && theNoise.size() != in.size()) {
0060       edm::LogError("TT6NoiseCalculator:updateNoise")
0061           << " You did not initialize the Noise correctly prior to noise calibration.";
0062     }
0063
0064     // Initialize sums used for estimating noise.
0065     if ((theStatus.isCalibrating() && numberOfEvents == 1) ||
0066         (theStatus.isUpdating() && (numberOfEvents - eventsRequiredToCalibrate_) % eventsRequiredToUpdate_ == 1)) {
0067       theNoiseSum.clear();
0068       theNoiseSqSum.clear();
0069       theEventPerStrip.clear();
0070
0071       theNoiseSum.resize(in.size(), 0.0);
0072       theNoiseSqSum.resize(in.size(), 0.0);
0073       theEventPerStrip.resize(in.size(), 0);
0074     }
0075
0076     unsigned int i;
0077
0078     // Update sums used for estimating noise.
0079     for (i = 0; i < in.size(); i++) {
0080       if (fabs(in[i]) < cutToAvoidSignal_ * theNoise[i]) {
0081         theNoiseSum[i] += in[i];
0082         theNoiseSqSum[i] += in[i] * in[i];
0083         theEventPerStrip[i]++;
0084       }
0085     }
0086
0087     // Calculate noise.
0088     if ((theStatus.isCalibrating() && numberOfEvents == eventsRequiredToCalibrate_) ||
0089         (theStatus.isUpdating() && (numberOfEvents - eventsRequiredToCalibrate_) % eventsRequiredToUpdate_ == 0)) {
0090       theCMPSubtractedSignal.clear();
0091       theNoise.clear();
0092
0093       for (i = 0; i < in.size(); i++) {
0094         double avVal = (theEventPerStrip[i]) ? theNoiseSum[i] / (theEventPerStrip[i]) : 0.0;
0095         double sqAvVal = (theEventPerStrip[i]) ? theNoiseSqSum[i] / (theEventPerStrip[i]) : 0.0;
0096         double corr_fac = (theEventPerStrip[i] > 1) ? (theEventPerStrip[i] / (theEventPerStrip[i] - 1)) : 1.0;
0097         double rmsVal = (sqAvVal - avVal * avVal > 0.0) ? sqrt(corr_fac * (sqAvVal - avVal * avVal)) : 0.0;
0098
0099         theCMPSubtractedSignal.push_back(static_cast<float>(avVal));
0100
0101         theNoise.push_back(static_cast<float>(rmsVal));
0102
0103         if (false)
0104           cout << " TT6NoiseCalculator::updateNoise " << theNoiseSum[i] << " " << theNoiseSqSum[i] << " "
0105                << theEventPerStrip[i] << " " << avVal << " " << sqAvVal << " " << (sqAvVal - avVal * avVal) << " "
0106                << rmsVal << endl;
0107       }
0108     }
0109     updateStatus();
0110   }
0111 }
0112 //
0113 // Define New Event
0114 //
0115 void TT6NoiseCalculator::newEvent() { alreadyUsedEvent = false; }