Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
#include "CalibTracker/SiStripAPVAnalysis/interface/SimpleNoiseCalculator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <cmath>
#include <numeric>
#include <algorithm>
using namespace std;
//
//  Constructors
//
SimpleNoiseCalculator::SimpleNoiseCalculator() : numberOfEvents(0), alreadyUsedEvent(false) {
  if (false)
    cout << "Constructing SimpleNoiseCalculator " << endl;
  init();
}
//
SimpleNoiseCalculator::SimpleNoiseCalculator(int evnt_ini, bool use_DB) : numberOfEvents(0), alreadyUsedEvent(false) {
  if (false)
    cout << "Constructing SimpleNoiseCalculator " << endl;
  useDB_ = use_DB;
  eventsRequiredToCalibrate_ = evnt_ini;
  //  eventsRequiredToUpdate_    = evnt_iter;
  //  cutToAvoidSignal_          = sig_cut;
  init();
}
//
// Initialization.
//
void SimpleNoiseCalculator::init() {
  theCMPSubtractedSignal.clear();
  theNoise.clear();
  theNoiseSum.clear();
  theNoiseSqSum.clear();
  theEventPerStrip.clear();
  // theStatus.setCalibrating();
}
//
//  Destructor
//
SimpleNoiseCalculator::~SimpleNoiseCalculator() {
  if (false)
    cout << "Destructing SimpleNoiseCalculator " << endl;
}
//
// Update the Status of Noise Calculation
//
void SimpleNoiseCalculator::updateStatus() {
  if ((theStatus.isCalibrating() && numberOfEvents >= eventsRequiredToCalibrate_) ||
      (useDB_ == true && numberOfEvents == 1)) {
    theStatus.setUpdating();
  }
}
//
// Calculate and update (when needed) Noise Values
//
void SimpleNoiseCalculator::updateNoise(ApvAnalysis::PedestalType& in) {
  if (alreadyUsedEvent == false) {
    alreadyUsedEvent = true;
    numberOfEvents++;

    if (numberOfEvents == 1 && theNoise.size() != in.size()) {
      edm::LogError("SimpleNoiseCalculator:updateNoise")
          << " You did not initialize the Noise correctly prior to noise calibration.";
    }

    // Initialize sums used for estimating noise.
    if (numberOfEvents == 1) {
      theNoiseSum.clear();
      theNoiseSqSum.clear();
      theEventPerStrip.clear();

      theNoiseSum.resize(in.size(), 0.0);
      theNoiseSqSum.resize(in.size(), 0.0);
      theEventPerStrip.resize(in.size(), 0);
    }

    unsigned int i;

    // At every event Update sums used for estimating noise.
    for (i = 0; i < in.size(); i++) {
      theNoiseSum[i] += in[i];
      theNoiseSqSum[i] += in[i] * in[i];
      theEventPerStrip[i]++;
    }

    // Calculate noise.
    if ((theStatus.isCalibrating() && numberOfEvents == eventsRequiredToCalibrate_) || theStatus.isUpdating()) {
      theCMPSubtractedSignal.clear();
      theNoise.clear();

      for (i = 0; i < in.size(); i++) {
        double avVal = (theEventPerStrip[i]) ? theNoiseSum[i] / (theEventPerStrip[i]) : 0.0;
        double sqAvVal = (theEventPerStrip[i]) ? theNoiseSqSum[i] / (theEventPerStrip[i]) : 0.0;
        double corr_fac = (theEventPerStrip[i] > 1) ? (theEventPerStrip[i] / (theEventPerStrip[i] - 1)) : 1.0;
        double rmsVal = (sqAvVal - avVal * avVal > 0.0) ? sqrt(corr_fac * (sqAvVal - avVal * avVal)) : 0.0;

        theCMPSubtractedSignal.push_back(static_cast<float>(avVal));

        theNoise.push_back(static_cast<float>(rmsVal));

        if (false)
          cout << " SimpleNoiseCalculator::updateNoise " << theNoiseSum[i] << " " << theNoiseSqSum[i] << " "
               << theEventPerStrip[i] << " " << avVal << " " << sqAvVal << " " << (sqAvVal - avVal * avVal) << " "
               << rmsVal << endl;
      }
    }
    updateStatus();
  }
}
//
// Define New Event
//
void SimpleNoiseCalculator::newEvent() { alreadyUsedEvent = false; }