TMom.cc

Go to the documentation of this file.

/* 
 *  \class TMom
 *
 *  \author: Julie Malcles - CEA/Saclay
 */

#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMarkov.h>
#include <TMath.h>

#include <cassert>

using namespace std;

//ClassImp(TMom)

// Default Constructor...
TMom::TMom() { init(0.0, 10.0e9); }

// Constructor...
TMom::TMom(double cutlow, double cuthigh) { init(cutlow, cuthigh); }
TMom::TMom(const std::vector<double>& cutlow, const std::vector<double>& cuthigh) { init(cutlow, cuthigh); }

// Destructor
TMom::~TMom() {}

void TMom::init(double cutlow, double cuthigh) {
  nevt = 0;
  mean = 0;
  mean2 = 0;
  mean3 = 0;
  sum = 0;
  sum2 = 0;
  sum3 = 0;
  rms = 0;
  M3 = 0;
  peak = 0;
  min = 10.0e9;
  max = 0.;
  _cutLow.clear();
  _cutHigh.clear();
  _dimCut = 1;
  _cutLow.push_back(cutlow);
  _cutHigh.push_back(cuthigh);
  for (int i = 0; i < 101; i++) {
    bing[i] = 0;
  }
}
void TMom::init(const std::vector<double>& cutlow, const std::vector<double>& cuthigh) {
  nevt = 0;
  mean = 0;
  mean2 = 0;
  mean3 = 0;
  sum = 0;
  sum2 = 0;
  sum3 = 0;
  rms = 0;
  M3 = 0;
  peak = 0;
  min = 10.0e9;
  max = 0.;
  assert(cutlow.size() == cuthigh.size());
  _cutLow.clear();
  _cutHigh.clear();
  _dimCut = cutlow.size();
  _cutLow = cutlow;
  _cutHigh = cuthigh;
  for (int i = 0; i < 101; i++) {
    bing[i] = 0;
  }
}
void TMom::setCut(double cutlow, double cuthigh) {
  _cutLow.clear();
  _cutHigh.clear();
  _dimCut = 1;
  _cutLow.push_back(cutlow);
  _cutHigh.push_back(cuthigh);
}
void TMom::setCut(const std::vector<double>& cutlow, const std::vector<double>& cuthigh) {
  assert(cutlow.size() == cuthigh.size());
  _cutLow.clear();
  _cutHigh.clear();
  _dimCut = cutlow.size();
  _cutLow = cutlow;
  _cutHigh = cuthigh;
}

void TMom::addEntry(double val) {
  std::vector<double> dumb;
  dumb.push_back(val);
  addEntry(val, dumb);
}

void TMom::addEntry(double val, const std::vector<double>& valcut) {
  int passingAllCuts = 1;

  for (int iCut = 0; iCut < _dimCut; iCut++) {
    int passing;
    if (valcut.at(iCut) > _cutLow.at(iCut) && valcut.at(iCut) <= _cutHigh.at(iCut)) {
      passing = 1;
    } else
      passing = 0;
    passingAllCuts *= passing;
  }

  if (passingAllCuts == 1) {
    nevt += 1;
    sum += val;
    sum2 += val * val;
    sum3 += val * val * val;
    if (val > max)
      max = val;
    if (val < min)
      min = val;

    // for peak stuff
    _ampl.push_back(val);
  }
}

double TMom::getMean() {
  if (nevt != 0)
    mean = sum / nevt;
  else
    mean = 0.0;
  return mean;
}

double TMom::getMean2() {
  if (nevt != 0)
    mean2 = sum2 / nevt;
  else
    mean2 = 0.0;
  return mean2;
}
double TMom::getMean3() {
  if (nevt != 0)
    mean3 = sum3 / nevt;
  else
    mean3 = 0.0;
  return mean3;
}

int TMom::getNevt() { return nevt; }

double TMom::getRMS() {
  double m = getMean();
  double m2 = getMean2();
  if (nevt != 0)
    rms = TMath::Sqrt(m2 - m * m);
  else
    rms = 0.0;
  return rms;
}

double TMom::getM3() {
  double m = getMean();
  double m2 = getMean2();
  double m3 = getMean3();
  double sig = getRMS();

  if (nevt != 0 && sig != 0)
    M3 = (m3 - 3.0 * m * (m2 - m * m) - m * m * m) / (sig * sig * sig);
  else
    M3 = 0.0;
  return M3;
}

double TMom::getMin() { return min; }
double TMom::getMax() { return max; }

std::vector<double> TMom::getPeak() {
  std::vector<double> p;
  double wbin = (max - min) / 100.;
  int bung;

  for (unsigned int i = 0; i < _ampl.size(); i++) {
    if (wbin <= 0.0)
      bung = 1;
    else
      bung = (int)((_ampl.at(i) - min) / wbin) + 1;
    if (1 <= bung && bung <= 100)
      bing[bung]++;
  }

  TMarkov* peakM = new TMarkov();

  int popmax = 0;

  for (int k = 1; k < 101; k++) {
    if (bing[k] > popmax) {
      popmax = bing[k];
    }
  }

  peakM->peakFinder(&bing[0]);
  p.push_back(peakM->getPeakValue(0));
  p.push_back(peakM->getPeakValue(1));

  return p;
}