MuonRPCGeometry/src/PTStatistics.cc

0001 #include "Validation/MuonRPCGeometry/interface/PTStatistics.h"
0002
0003 #include <sstream>
0004 #include <numeric>
0005 #include <iostream>
0006 //-----------------------------------------------------------------------------
0007 //  PTStatistics
0008 //-----------------------------------------------------------------------------
0009 static std::vector<long double> initializeRates() {
0010   std::vector<long double> rates(RPCpg::ptBins_s, 0);
0011   //std::cout << "Initilizing rates" << std::endl;
0012
0013   // Note bin=0 is empty during generation
0014   // bin=0 is used only when calculating efficiencies (for storing muons,that werent found)
0015   for (unsigned int i = 1; i < rates.size(); ++i) {
0016     long double low = RPCpg::pts[i];
0017     long double high = RPCpg::pts[i + 1];
0018     long double rt = RPCpg::rate(low) - RPCpg::rate(high);
0019
0020     /* std::cout << "PtCode " << i
0021               << " " << low
0022               << " " << high
0023               << " " << rt
0024               << std::endl;*/
0025     rates[i] = rt;
0026   }
0027   return rates;
0028 }
0029
0030 const std::vector<long double> PTStatistics::m_rates = initializeRates();
0031
0032 PTStatistics::PTStatistics() { this->assign(RPCpg::ptBins_s, 0); }
0033
0034 void PTStatistics::update(PTStatistics& otherPtStats) {
0035   for (unsigned int i = 0; i < this->size(); ++i) {
0036     //this->at(i)+=otherPtStats.at(i);
0037     (*this)[i] += otherPtStats[i];
0038   }
0039 }
0040 std::string PTStatistics::toString() {
0041   std::stringstream ss;
0042   ss << "PTStats:";
0043   for (unsigned int i = 0; i < this->size(); ++i) {
0044     ss << " " << (*this)[i];
0045   }
0046
0047   return ss.str();
0048 }
0049
0050 long double PTStatistics::eff(int ptCut) {  // ptCut=0 -> total rate
0051   //int eqOrAbovePtCut = 0;
0052   //for(unsigned int i=ptCut;i<this->size();++i) eqOrAbovePtCut += this->at(i);
0053   // return double(eqOrAbovePtCut)/this->sum();
0054   return double(sum(ptCut)) / this->sum();
0055 }
0056
0057 long int PTStatistics::sum(const int& ptCut) const {
0058   //inline int PTStatistics::sum(const int & ptCut) const{
0059   //return std::accumulate(this->begin(),this->end(),0);
0060   long int eqOrAbovePtCut = 0;
0061   unsigned int size = this->size();
0062   //for(unsigned int i=ptCut;i<size;++i) eqOrAbovePtCut += this->at(i);
0063   for (unsigned int i = ptCut; i < size; ++i)
0064     eqOrAbovePtCut += (*this)[i];
0065   return eqOrAbovePtCut;
0066 }
0067
0068 //inline int PTStatistics::sum() const{
0069 long int PTStatistics::sum() const {
0070   //return std::accumulate(this->begin(),this->end(),0);
0071   long int eqOrAbovePtCut = 0;
0072   //unsigned int size = this->size();
0073   //for(unsigned int i=ptCut;i<size;++i) eqOrAbovePtCut += this->at(i);
0074   //for(unsigned int i=ptCut;i<size;++i) eqOrAbovePtCut += (*this)[i];
0075   PTStatistics::const_iterator it = this->begin();
0076   PTStatistics::const_iterator itend = this->end();
0077   for (; it != itend; ++it)
0078     eqOrAbovePtCut += *it;
0079
0080   return eqOrAbovePtCut;
0081 }
0082
0083 long double PTStatistics::sumR(const int& ptCut) const {
0084   //return std::accumulate(this->begin(),this->end(),0);
0085   long double eqOrAbovePtCut = 0;
0086   unsigned int size = this->size();
0087   //for(unsigned int i=ptCut;i<size;++i) eqOrAbovePtCut += this->at(i);
0088   for (unsigned int i = ptCut; i < size; ++i)
0089     eqOrAbovePtCut += (*this)[i] * m_rates[i];
0090   return eqOrAbovePtCut;
0091 }
0092
0093 long double PTStatistics::sumR() const {
0094   //return std::accumulate(this->begin(),this->end(),0);
0095   long double eqOrAbovePtCut = 0;
0096   //unsigned int size = this->size();
0097   //for(unsigned int i=ptCut;i<size;++i) eqOrAbovePtCut += this->at(i);
0098   //for(unsigned int i=ptCut;i<size;++i) eqOrAbovePtCut += (*this)[i];
0099   //PTStatistics::const_iterator it = this->begin();
0100   //PTStatistics::const_iterator itend = this->end();
0101   //PtCode 0 - muons not found
0102   unsigned int size = this->size();
0103   for (unsigned int i = 1; i < size; ++i)
0104     eqOrAbovePtCut += (*this)[i] * m_rates[i];
0105   //for(;it!=itend;++it) eqOrAbovePtCut += *it*m_rates[i];
0106
0107   return eqOrAbovePtCut;
0108 }