include/Geant4/G4StatAnalysis.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4StatAnalysis inline methods implementation
0027 //
0028 // Author: J.Madsen, 25.10.2018
0029 // --------------------------------------------------------------------
0030
0031 #include <cmath>
0032 #include <fstream>
0033 #include <iomanip>
0034 #include <iostream>
0035 #include <limits>
0036
0037 #include "G4Timer.hh"
0038 #include "G4ios.hh"
0039 #include "globals.hh"
0040 #include "tls.hh"
0041
0042 #include "G4Allocator.hh"
0043 #include "G4Types.hh"
0044
0045 //----------------------------------------------------------------------------//
0046
0047 G4StatAnalysis::G4StatAnalysis() {}
0048
0049 //----------------------------------------------------------------------------//
0050
0051 G4double G4StatAnalysis::GetMean() const
0052 {
0053   return (fHits > 0) ? fSum1 / ((G4double) fHits) : 0.;
0054 }
0055
0056 //----------------------------------------------------------------------------//
0057
0058 const G4double& G4StatAnalysis::GetSum() const { return fSum1; }
0059
0060 //----------------------------------------------------------------------------//
0061
0062 const G4double& G4StatAnalysis::GetSumSquared() const { return fSum2; }
0063
0064 //----------------------------------------------------------------------------//
0065
0066 const G4double& G4StatAnalysis::GetSum1() const { return fSum1; }
0067
0068 //----------------------------------------------------------------------------//
0069
0070 const G4double& G4StatAnalysis::GetSum2() const { return fSum2; }
0071
0072 //----------------------------------------------------------------------------//
0073
0074 const G4int& G4StatAnalysis::GetHits() const { return fHits; }
0075
0076 //----------------------------------------------------------------------------//
0077
0078 G4int G4StatAnalysis::GetNumNonZero() const { return fHits - fZero; }
0079
0080 //----------------------------------------------------------------------------//
0081
0082 G4int G4StatAnalysis::GetNumZero() const { return fZero; }
0083
0084 //----------------------------------------------------------------------------//
0085
0086 void G4StatAnalysis::SetSum(const G4double& val) { fSum1 = val; }
0087
0088 //----------------------------------------------------------------------------//
0089
0090 void G4StatAnalysis::SetSumSquared(const G4double& val) { fSum2 = val; }
0091
0092 //----------------------------------------------------------------------------//
0093
0094 void G4StatAnalysis::SetSum1(const G4double& val) { fSum1 = val; }
0095
0096 //----------------------------------------------------------------------------//
0097
0098 void G4StatAnalysis::SetSum2(const G4double& val) { fSum2 = val; }
0099
0100 //----------------------------------------------------------------------------//
0101
0102 void G4StatAnalysis::SetHits(const G4int& val) { fHits = val; }
0103
0104 //----------------------------------------------------------------------------//
0105
0106 void G4StatAnalysis::SetZero(const G4int& val) { fZero = val; }
0107
0108 //----------------------------------------------------------------------------//
0109
0110 G4double G4StatAnalysis::GetFOM() const
0111 {
0112   G4double elapsed_time = this->GetCpuTime();
0113   G4double relative_err = this->GetRelativeError();
0114   // lambda for equation clarity (will be inlined)
0115   auto compute_figure_of_merit = [&]() {
0116     return (1.0 / (relative_err * relative_err) / elapsed_time);
0117   };
0118   return (std::fabs(relative_err) > 0.0 && elapsed_time > 0.0)
0119            ? compute_figure_of_merit()
0120            : ((fHits > 0) ? 1.0 : 0.0);
0121 }
0122
0123 //----------------------------------------------------------------------------//
0124
0125 G4double G4StatAnalysis::GetRelativeError() const
0126 {
0127   // lambda for equation clarity (will be inlined)
0128   auto compute_relative_error = [&]() {
0129     return (GetStdDev() / GetMean() / std::sqrt((G4double) fHits));
0130   };
0131   return (std::fabs(GetMean()) > 0 && fHits > 0) ? compute_relative_error()
0132                                                  : ((fHits > 0) ? 1.0 : 0.0);
0133 }
0134
0135 //----------------------------------------------------------------------------//
0136
0137 G4double G4StatAnalysis::GetStdDev() const
0138 {
0139   return ::sqrt(std::fabs(GetVariance()));
0140 }
0141
0142 //----------------------------------------------------------------------------//
0143
0144 G4double G4StatAnalysis::GetVariance() const
0145 {
0146   // lambda for equation clarity (will be inlined)
0147   auto compute_variance = [&]() {
0148     return ((fSum2 - (std::pow(fSum1, 2.0) / fHits)) /
0149             (((G4double) fHits) - 1.0));
0150   };
0151   return (fHits > 1) ? compute_variance() : 0.0;
0152 }
0153
0154 //----------------------------------------------------------------------------//
0155
0156 G4double G4StatAnalysis::GetCoeffVariation() const
0157 {
0158   // lambda for equation clarity (will be inlined)
0159   auto coefficient_of_variation = [&]() {
0160     G4double hits = fHits;
0161     return ::sqrt((hits / (hits - 1.0)) *
0162                   ((fSum2 / (fSum1 * fSum1)) - (1.0 / hits)));
0163   };
0164   return (fHits > 1) ? coefficient_of_variation() : 0.0;
0165   // return (fHits > 0 && fabs(fSum1) > 0.0)
0166   //        ? (100.0*GetStdDev()/GetMean()) : 0.0;
0167 }
0168
0169 //----------------------------------------------------------------------------//
0170
0171 G4double G4StatAnalysis::GetEfficiency() const
0172 {
0173   G4double hits  = fHits;
0174   G4double nzero = fHits - fZero;
0175   return (fHits > 0) ? (nzero / hits) : 0.0;
0176 }
0177
0178 //----------------------------------------------------------------------------//
0179
0180 G4double G4StatAnalysis::GetR2Int() const
0181 {
0182   G4double hits = fHits;
0183   return (fHits > 0)
0184            ? (fSum2 / (fSum1 * fSum1)) - 1.0 / (GetEfficiency() * hits)
0185            : 0.0;
0186 }
0187
0188 //----------------------------------------------------------------------------//
0189
0190 G4double G4StatAnalysis::GetR2Eff() const
0191 {
0192   G4double hits = fHits;
0193   return (fHits > 0) ? (1.0 - GetEfficiency()) / (GetEfficiency() * hits) : 0.0;
0194 }
0195
0196 //----------------------------------------------------------------------------//
0197
0198 G4StatAnalysis::operator G4double() const { return this->GetSum(); }
0199
0200 //----------------------------------------------------------------------------//
0201
0202 void G4StatAnalysis::Reset()
0203 {
0204   fHits = 0;
0205   fZero = 0;
0206   fSum1 = 0.0;
0207   fSum2 = 0.0;
0208 }
0209
0210 //----------------------------------------------------------------------------//
0211
0212 void G4StatAnalysis::Add(const G4double& val, const G4double& weight)
0213 {
0214   fHits += 1;
0215   fSum1 += val * weight;
0216   fSum2 += val * val * weight;
0217   if(std::fabs(val * weight) <
0218      std::fabs(GetMean() * std::numeric_limits<double>::epsilon()))
0219     fZero += 1;
0220 }
0221
0222 //----------------------------------------------------------------------------//
0223
0224 void G4StatAnalysis::Rescale(const G4double& factor)
0225 {
0226   fSum1 *= factor;
0227   fSum2 *= factor * factor;
0228 }
0229
0230 //----------------------------------------------------------------------------//
0231
0232 void G4StatAnalysis::PrintInfo(std::ostream& os, const std::string& tab) const
0233 {
0234   G4int _hits     = this->GetHits();
0235   G4double _sum   = this->GetSum();
0236   G4double _sigma = this->GetStdDev();
0237   G4double _coeff = this->GetCoeffVariation();
0238   G4double _error = this->GetRelativeError();
0239   G4double _eff   = this->GetEfficiency();
0240   G4double _fom   = this->GetFOM();
0241   G4double _r2int = this->GetR2Int();
0242   G4double _r2eff = this->GetR2Eff();
0243
0244   using std::fixed;
0245   using std::ios;
0246   using std::left;
0247   using std::right;
0248   using std::scientific;
0249   using std::setprecision;
0250   using std::setw;
0251
0252   std::stringstream ss;
0253   ss << tab  //<< scientific
0254      << setprecision((G4int)os.precision()) << right << _sum << left
0255      << " [sigma: " << right << _sigma << left << " | error: " << right
0256      << _error << left << " | coeff: " << right << _coeff << left
0257      << " | eff: " << right << _eff << left << " | fom: " << right << _fom
0258      << left << " | r2int: " << right << _r2int << left << " | r2eff: " << right
0259      << _r2eff << left << " | hits: " << right << _hits << left << " ]";
0260
0261   os << ss.str();
0262 }
0263
0264 //----------------------------------------------------------------------------//
0265
0266 G4double G4StatAnalysis::GetCpuTime() const
0267 {
0268   tms* startTime = GetCpuClock();
0269   tms endTime;
0270   times(&endTime);
0271   return ((endTime.tms_stime - startTime->tms_stime) +
0272           (endTime.tms_utime - startTime->tms_utime)) /
0273          sysconf(_SC_CLK_TCK);
0274 }
0275
0276 //----------------------------------------------------------------------------//
0277
0278 G4StatAnalysis& G4StatAnalysis::operator+=(const G4double& _val)
0279 {
0280   this->Add(_val);
0281   return *this;
0282 }
0283
0284 //----------------------------------------------------------------------------//
0285
0286 G4StatAnalysis& G4StatAnalysis::operator/=(const G4double& _val)
0287 {
0288   fSum1 /= _val;
0289   fSum2 /= (_val * _val);
0290   return *this;
0291 }
0292
0293 //----------------------------------------------------------------------------//
0294
0295 G4StatAnalysis& G4StatAnalysis::operator+=(const G4StatAnalysis& rhs)
0296 {
0297   fHits += rhs.fHits;
0298   fSum1 += rhs.fSum1;
0299   fSum2 += rhs.fSum2;
0300   fZero += rhs.fZero;
0301   return *this;
0302 }
0303
0304 //----------------------------------------------------------------------------//
0305
0306 G4StatAnalysis& G4StatAnalysis::operator-=(const G4StatAnalysis& rhs)
0307 {
0308   fHits -= rhs.fHits;
0309   fSum1 -= rhs.fSum1;
0310   fSum2 -= rhs.fSum2;
0311   fZero -= rhs.fZero;
0312   return *this;
0313 }
0314
0315 //----------------------------------------------------------------------------//
0316
0317 inline G4Allocator<G4StatAnalysis>*& _aStatAnalysisAllocator_G4MT_TLS_()
0318 {
0319   G4ThreadLocalStatic G4Allocator<G4StatAnalysis>* _instance =
0320     new G4Allocator<G4StatAnalysis>();
0321   return _instance;
0322 }
0323
0324 //----------------------------------------------------------------------------//
0325
0326 void* G4StatAnalysis::operator new(std::size_t)
0327 {
0328   G4Allocator<G4StatAnalysis>& _allocator =
0329     *_aStatAnalysisAllocator_G4MT_TLS_();
0330   return (void*) _allocator.MallocSingle();
0331 }
0332
0333 //----------------------------------------------------------------------------//
0334
0335 void G4StatAnalysis::operator delete(void* _ptr)
0336 {
0337   G4Allocator<G4StatAnalysis>& _allocator =
0338     *_aStatAnalysisAllocator_G4MT_TLS_();
0339   _allocator.FreeSingle((G4StatAnalysis*) _ptr);
0340 }
0341
0342 //----------------------------------------------------------------------------//