EcalLaserAnalyzer/src/TMom.cc

0001 /*
0002  *  \class TMom
0003  *
0004  *  \author: Julie Malcles - CEA/Saclay
0005  */
0006
0007 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
0008 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMarkov.h>
0009 #include <TMath.h>
0010
0011 #include <cassert>
0012
0013 using namespace std;
0014
0015 //ClassImp(TMom)
0016
0017 // Default Constructor...
0018 TMom::TMom() { init(0.0, 10.0e9); }
0019
0020 // Constructor...
0021 TMom::TMom(double cutlow, double cuthigh) { init(cutlow, cuthigh); }
0022 TMom::TMom(const std::vector<double>& cutlow, const std::vector<double>& cuthigh) { init(cutlow, cuthigh); }
0023
0024 // Destructor
0025 TMom::~TMom() {}
0026
0027 void TMom::init(double cutlow, double cuthigh) {
0028   nevt = 0;
0029   mean = 0;
0030   mean2 = 0;
0031   mean3 = 0;
0032   sum = 0;
0033   sum2 = 0;
0034   sum3 = 0;
0035   rms = 0;
0036   M3 = 0;
0037   peak = 0;
0038   min = 10.0e9;
0039   max = 0.;
0040   _cutLow.clear();
0041   _cutHigh.clear();
0042   _dimCut = 1;
0043   _cutLow.push_back(cutlow);
0044   _cutHigh.push_back(cuthigh);
0045   for (int i = 0; i < 101; i++) {
0046     bing[i] = 0;
0047   }
0048 }
0049 void TMom::init(const std::vector<double>& cutlow, const std::vector<double>& cuthigh) {
0050   nevt = 0;
0051   mean = 0;
0052   mean2 = 0;
0053   mean3 = 0;
0054   sum = 0;
0055   sum2 = 0;
0056   sum3 = 0;
0057   rms = 0;
0058   M3 = 0;
0059   peak = 0;
0060   min = 10.0e9;
0061   max = 0.;
0062   assert(cutlow.size() == cuthigh.size());
0063   _cutLow.clear();
0064   _cutHigh.clear();
0065   _dimCut = cutlow.size();
0066   _cutLow = cutlow;
0067   _cutHigh = cuthigh;
0068   for (int i = 0; i < 101; i++) {
0069     bing[i] = 0;
0070   }
0071 }
0072 void TMom::setCut(double cutlow, double cuthigh) {
0073   _cutLow.clear();
0074   _cutHigh.clear();
0075   _dimCut = 1;
0076   _cutLow.push_back(cutlow);
0077   _cutHigh.push_back(cuthigh);
0078 }
0079 void TMom::setCut(const std::vector<double>& cutlow, const std::vector<double>& cuthigh) {
0080   assert(cutlow.size() == cuthigh.size());
0081   _cutLow.clear();
0082   _cutHigh.clear();
0083   _dimCut = cutlow.size();
0084   _cutLow = cutlow;
0085   _cutHigh = cuthigh;
0086 }
0087
0088 void TMom::addEntry(double val) {
0089   std::vector<double> dumb;
0090   dumb.push_back(val);
0091   addEntry(val, dumb);
0092 }
0093
0094 void TMom::addEntry(double val, const std::vector<double>& valcut) {
0095   int passingAllCuts = 1;
0096
0097   for (int iCut = 0; iCut < _dimCut; iCut++) {
0098     int passing;
0099     if (valcut.at(iCut) > _cutLow.at(iCut) && valcut.at(iCut) <= _cutHigh.at(iCut)) {
0100       passing = 1;
0101     } else
0102       passing = 0;
0103     passingAllCuts *= passing;
0104   }
0105
0106   if (passingAllCuts == 1) {
0107     nevt += 1;
0108     sum += val;
0109     sum2 += val * val;
0110     sum3 += val * val * val;
0111     if (val > max)
0112       max = val;
0113     if (val < min)
0114       min = val;
0115
0116     // for peak stuff
0117     _ampl.push_back(val);
0118   }
0119 }
0120
0121 double TMom::getMean() {
0122   if (nevt != 0)
0123     mean = sum / nevt;
0124   else
0125     mean = 0.0;
0126   return mean;
0127 }
0128
0129 double TMom::getMean2() {
0130   if (nevt != 0)
0131     mean2 = sum2 / nevt;
0132   else
0133     mean2 = 0.0;
0134   return mean2;
0135 }
0136 double TMom::getMean3() {
0137   if (nevt != 0)
0138     mean3 = sum3 / nevt;
0139   else
0140     mean3 = 0.0;
0141   return mean3;
0142 }
0143
0144 int TMom::getNevt() { return nevt; }
0145
0146 double TMom::getRMS() {
0147   double m = getMean();
0148   double m2 = getMean2();
0149   if (nevt != 0)
0150     rms = TMath::Sqrt(m2 - m * m);
0151   else
0152     rms = 0.0;
0153   return rms;
0154 }
0155
0156 double TMom::getM3() {
0157   double m = getMean();
0158   double m2 = getMean2();
0159   double m3 = getMean3();
0160   double sig = getRMS();
0161
0162   if (nevt != 0 && sig != 0)
0163     M3 = (m3 - 3.0 * m * (m2 - m * m) - m * m * m) / (sig * sig * sig);
0164   else
0165     M3 = 0.0;
0166   return M3;
0167 }
0168
0169 double TMom::getMin() { return min; }
0170 double TMom::getMax() { return max; }
0171
0172 std::vector<double> TMom::getPeak() {
0173   std::vector<double> p;
0174   double wbin = (max - min) / 100.;
0175   int bung;
0176
0177   for (unsigned int i = 0; i < _ampl.size(); i++) {
0178     if (wbin <= 0.0)
0179       bung = 1;
0180     else
0181       bung = (int)((_ampl.at(i) - min) / wbin) + 1;
0182     if (1 <= bung && bung <= 100)
0183       bing[bung]++;
0184   }
0185
0186   TMarkov* peakM = new TMarkov();
0187
0188   int popmax = 0;
0189
0190   for (int k = 1; k < 101; k++) {
0191     if (bing[k] > popmax) {
0192       popmax = bing[k];
0193     }
0194   }
0195
0196   peakM->peakFinder(&bing[0]);
0197   p.push_back(peakM->getPeakValue(0));
0198   p.push_back(peakM->getPeakValue(1));
0199
0200   return p;
0201 }