EcalLaserAnalyzer/src/TAPD.cc

0001 /*
0002  *  \class TAPD
0003  *
0004  *  \author: Julie Malcles  - CEA/Saclay
0005  */
0006
0007 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMom.h>
0008 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TAPD.h>
0009 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMarkov.h>
0010 #include <TMath.h>
0011
0012 using namespace std;
0013 #include <iostream>
0014 #include <cassert>
0015
0016 //ClassImp(TAPD)
0017
0018 // Default Constructor...
0019 TAPD::TAPD() { init(); }
0020
0021 // Destructor
0022 TAPD::~TAPD() {}
0023
0024 void TAPD::init() {
0025   for (int j = 0; j < nOutVar; j++) {
0026     _apdcuts[0][j].clear();
0027     _apdcuts[1][j].clear();
0028     _cutvars[j].clear();
0029
0030     _apdcuts[0][j].push_back(0.0);
0031     _apdcuts[1][j].push_back(10.0e6);
0032     _cutvars[j].push_back(j);
0033
0034     mom[j] = new TMom();
0035   }
0036 }
0037
0038 void TAPD::addEntry(double apd, double pn, double pn0, double pn1, double time) {
0039   addEntry(apd, pn, pn0, pn1, time, 0.0, 0.0);
0040 }
0041
0042 void TAPD::addEntry(double apd, double pn, double pn0, double pn1, double time, double apd0, double apd1) {
0043   double val[nOutVar];
0044   std::vector<double> valcuts[nOutVar];
0045
0046   val[iAPD] = apd;
0047   if (pn != 0)
0048     val[iAPDoPN] = apd / pn;
0049   else
0050     val[iAPDoPN] = 0.0;
0051   if (pn0 != 0)
0052     val[iAPDoPN0] = apd / pn0;
0053   else
0054     val[iAPDoPN0] = 0.0;
0055   if (pn1 != 0)
0056     val[iAPDoPN1] = apd / pn1;
0057   else
0058     val[iAPDoPN1] = 0.0;
0059   val[iTime] = time;
0060   if (apd0 != 0.)
0061     val[iAPDoAPD0] = apd / apd0;
0062   else
0063     val[iAPDoAPD0] = 0.0;
0064   if (apd1 != 0.)
0065     val[iAPDoAPD1] = apd / apd1;
0066   else
0067     val[iAPDoAPD1] = 0.0;
0068
0069   for (int ivar = 0; ivar < nOutVar; ivar++) {
0070     int dimcut = _cutvars[ivar].size();
0071     for (int ic = 0; ic < dimcut; ic++) {
0072       assert(_cutvars[ivar].at(ic) < nOutVar);
0073       valcuts[ivar].push_back(val[_cutvars[ivar].at(ic)]);
0074     }
0075   }
0076
0077   for (int ivar = 0; ivar < nOutVar; ivar++) {
0078     mom[ivar]->addEntry(val[ivar], valcuts[ivar]);
0079     //    std::cout << "addEntry: val[ivar=" << ivar <<"] = "<<val[ivar]<< std::endl;
0080
0081     for (size_t ic = 0; ic < _cutvars[ivar].size(); ic++) {
0082       //      std::cout << "addEntry: valcuts[ivar="<< ivar <<"][ic="<<ic<<"] = "<<valcuts[ivar].at(ic)<< std::endl;
0083       for (size_t iv = 0; iv < _cutvars[ivar].size(); iv++) {
0084         //  std::cout <<"low cut:"<<_apdcuts[0][ivar].at(iv)<<", high cut:"<<_apdcuts[1][ivar].at(iv)<<", cutvar: "<<_cutvars[ivar].at(iv)<< std::endl;
0085       }
0086     }
0087   }
0088 }
0089
0090 void TAPD::setCut(int ivar, double mean, double sig) {
0091   assert(ivar < nOutVar);
0092
0093   std::vector<int> cutvar;
0094   cutvar.push_back(ivar);
0095
0096   std::vector<double> lowcut;
0097   std::vector<double> highcut;
0098
0099   double low = mean - 2.0 * sig;
0100   if (low < 0)
0101     low = 0.0;
0102   double high = mean + 2.0 * sig;
0103
0104   lowcut.push_back(low);
0105   highcut.push_back(high);
0106
0107   setCut(ivar, cutvar, lowcut, highcut);
0108 }
0109
0110 void TAPD::setCut(int ivar,
0111                   const std::vector<int>& cutVars,
0112                   const std::vector<double>& lowCut,
0113                   const std::vector<double>& highCut) {
0114   assert(ivar < nOutVar);
0115   int cutdim = cutVars.size();
0116   assert(cutdim < nOutVar);
0117   assert(cutdim == (int)lowCut.size());
0118   assert(cutdim == (int)highCut.size());
0119
0120   _apdcuts[0][ivar].clear();
0121   _apdcuts[1][ivar].clear();
0122   _cutvars[ivar].clear();
0123
0124   for (int ic = 0; ic < cutdim; ic++) {
0125     // FINISH THIS
0126     if (lowCut.at(ic) > 0) {
0127       _apdcuts[0][ivar].push_back(lowCut.at(ic));
0128     } else
0129       _apdcuts[0][ivar].push_back(0.0);
0130
0131     _apdcuts[1][ivar].push_back(highCut.at(ic));
0132     _cutvars[ivar].push_back(cutVars.at(ic));
0133   }
0134
0135   mom[ivar]->setCut(_apdcuts[0][ivar], _apdcuts[1][ivar]);
0136 }
0137
0138 // Simple 1D cuts on main variable at 2 sigmas
0139 // ===========================================
0140
0141 void TAPD::setAPDCut(double mean, double sig) { setCut(TAPD::iAPD, mean, sig); }
0142 void TAPD::setAPDoPNCut(double mean, double sig) { setCut(TAPD::iAPDoPN, mean, sig); }
0143 void TAPD::setAPDoPN0Cut(double mean, double sig) { setCut(TAPD::iAPDoPN0, mean, sig); }
0144 void TAPD::setAPDoPN1Cut(double mean, double sig) { setCut(TAPD::iAPDoPN1, mean, sig); }
0145 void TAPD::setTimeCut(double mean, double sig) { setCut(TAPD::iTime, mean, sig); }
0146
0147 // More complicated 2D cuts
0148 // =========================
0149
0150 // Cut on main var and Time:
0151 void TAPD::set2DCut(int ivar, const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0152   assert(lowCut.size() == 2);
0153   assert(highCut.size() == 2);
0154   std::vector<int> cutVars;
0155   cutVars.push_back(ivar);
0156   cutVars.push_back(TAPD::iTime);
0157   setCut(ivar, cutVars, lowCut, highCut);
0158 }
0159
0160 void TAPD::set2DAPDCut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0161   set2DCut(TAPD::iAPD, lowCut, highCut);
0162 }
0163 void TAPD::set2DAPDoPNCut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0164   set2DCut(TAPD::iAPDoPN, lowCut, highCut);
0165 }
0166 void TAPD::set2DAPDoPN0Cut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0167   set2DCut(TAPD::iAPDoPN0, lowCut, highCut);
0168 }
0169 void TAPD::set2DAPDoPN1Cut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0170   set2DCut(TAPD::iAPDoPN1, lowCut, highCut);
0171 }
0172
0173 void TAPD::set2DAPDoAPD0Cut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0174   assert(lowCut.size() == 2);
0175   assert(highCut.size() == 2);
0176   std::vector<int> cutVars;
0177   cutVars.push_back(TAPD::iAPD);
0178   cutVars.push_back(TAPD::iTime);
0179   setCut(TAPD::iAPDoAPD0, cutVars, lowCut, highCut);
0180 }
0181 void TAPD::set2DAPDoAPD1Cut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0182   assert(lowCut.size() == 2);
0183   assert(highCut.size() == 2);
0184   std::vector<int> cutVars;
0185   cutVars.push_back(TAPD::iAPD);
0186   cutVars.push_back(TAPD::iTime);
0187   setCut(TAPD::iAPDoAPD1, cutVars, lowCut, highCut);
0188 }
0189
0190 void TAPD::set2DTimeCut(const std::vector<double>& lowCut, const std::vector<double>& highCut) {
0191   assert(lowCut.size() == 2);
0192   assert(highCut.size() == 2);
0193   std::vector<int> cutVars;
0194   cutVars.push_back(TAPD::iAPD);
0195   cutVars.push_back(TAPD::iTime);
0196   setCut(TAPD::iTime, cutVars, lowCut, highCut);
0197 }
0198
0199 std::vector<double> TAPD::get(int ivar) {
0200   std::vector<double> res;
0201
0202   if (ivar < nOutVar) {
0203     res.push_back(mom[ivar]->getMean());
0204     res.push_back(mom[ivar]->getRMS());
0205     res.push_back(mom[ivar]->getM3());
0206     res.push_back(mom[ivar]->getNevt());
0207     res.push_back(mom[ivar]->getMin());
0208     res.push_back(mom[ivar]->getMax());
0209   }
0210
0211   //  std::cout << "In get: ivar="<< ivar << ", mean="<< mom[ivar]->getMean()<<" res size="<< res.size()<< std::endl;
0212
0213   return res;
0214 }
0215
0216 std::vector<double> TAPD::getAPD() {
0217   std::vector<double> x = get(TAPD::iAPD);
0218   return x;
0219 }
0220 std::vector<double> TAPD::getAPDoPN() {
0221   std::vector<double> x = get(TAPD::iAPDoPN);
0222   return x;
0223 }
0224 std::vector<double> TAPD::getAPDoPN0() {
0225   std::vector<double> x = get(TAPD::iAPDoPN0);
0226   return x;
0227 }
0228 std::vector<double> TAPD::getAPDoPN1() {
0229   std::vector<double> x = get(TAPD::iAPDoPN1);
0230   return x;
0231 }
0232 std::vector<double> TAPD::getTime() {
0233   std::vector<double> x = get(TAPD::iTime);
0234   return x;
0235 }
0236 std::vector<double> TAPD::getAPDoAPD0() {
0237   std::vector<double> x = get(TAPD::iAPDoAPD0);
0238   // std::cout<< "In GetAPDoAPD0: x[0]="<< x.at(0) << std::endl;
0239   return x;
0240 }
0241 std::vector<double> TAPD::getAPDoAPD1() {
0242   std::vector<double> x = get(TAPD::iAPDoAPD1);
0243   return x;
0244 }