Egamma/src/EffectiveAreas.cc

0001 #include "CommonTools/Egamma/interface/EffectiveAreas.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 #include <cmath>
0005 #include <fstream>
0006 #include <string>
0007 #include <sstream>
0008
0009 EffectiveAreas::EffectiveAreas(const std::string& filename, const bool quadraticEAflag)
0010     : filename_(filename), quadraticEAflag_(quadraticEAflag) {
0011   // Open the file with the effective area constants
0012   std::ifstream inputFile;
0013   inputFile.open(filename_.c_str());
0014   if (!inputFile.is_open())
0015     throw cms::Exception("EffectiveAreas config failure") << "failed to open the file " << filename_ << std::endl;
0016
0017   // Read file line by line
0018   std::string line;
0019   const float undef = -999;
0020
0021   if (!quadraticEAflag_) {
0022     while (getline(inputFile, line)) {
0023       if (line[0] == '#')
0024         continue;  // skip the comments lines
0025       float etaMin = undef, etaMax = undef, effArea = undef;
0026       std::stringstream ss(line);
0027       ss >> etaMin >> etaMax >> effArea;
0028       // In case if the format is messed up, there are letters
0029       // instead of numbers, or not exactly three numbers in the line,
0030       // it is likely that one or more of these vars never changed
0031       // the original "undef" value:
0032       if (etaMin == undef || etaMax == undef || effArea == undef)
0033         throw cms::Exception("EffectiveAreas config failure")
0034             << "wrong file format, file name " << filename_ << std::endl;
0035
0036       absEtaMin_.push_back(etaMin);
0037       absEtaMax_.push_back(etaMax);
0038       effectiveAreaValues_.push_back(effArea);
0039     }
0040
0041   }
0042
0043   else {
0044     while (getline(inputFile, line)) {
0045       if (line[0] == '#')
0046         continue;  // skip the comments lines
0047
0048       float etaMin = undef, etaMax = undef;
0049       float linEffArea = undef, quadEffArea = undef;
0050
0051       std::stringstream ss(line);
0052
0053       ss >> etaMin >> etaMax >> linEffArea >> quadEffArea;
0054
0055       // In case if the format is messed up, there are letters
0056       // instead of numbers, or not exactly three numbers in the line,
0057       // it is likely that one or more of these vars never changed
0058       // the original "undef" value:
0059       if (etaMin == undef || etaMax == undef || linEffArea == undef || quadEffArea == undef)
0060         throw cms::Exception("EffectiveAreas config failure")
0061             << "wrong file format, file name " << filename_ << std::endl;
0062
0063       absEtaMin_.push_back(etaMin);
0064       absEtaMax_.push_back(etaMax);
0065       linearEffectiveAreaValues_.push_back(linEffArea);
0066       quadraticEffectiveAreaValues_.push_back(quadEffArea);
0067     }
0068   }
0069
0070   // Extra consistency checks are in the function below.
0071   // If any of them fail, an exception is thrown.
0072   checkConsistency();
0073 }
0074
0075 // Return effective area for given eta
0076 const float EffectiveAreas::getEffectiveArea(float eta) const {
0077   float effArea = 0;
0078   uint nEtaBins = absEtaMin_.size();
0079   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0080     if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
0081       effArea = effectiveAreaValues_[iEta];
0082       break;
0083     }
0084   }
0085
0086   return effArea;
0087 }
0088
0089 // Return linear term of EA for given eta
0090 const float EffectiveAreas::getLinearEA(float eta) const {
0091   float linEffArea = 0;
0092   uint nEtaBins = absEtaMin_.size();
0093
0094   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0095     if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
0096       linEffArea = linearEffectiveAreaValues_[iEta];
0097       break;
0098     }
0099   }
0100
0101   return linEffArea;
0102 }
0103
0104 // Return quadratic term of EA for given eta
0105 const float EffectiveAreas::getQuadraticEA(float eta) const {
0106   float quadEffArea = 0;
0107   uint nEtaBins = absEtaMin_.size();
0108
0109   for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0110     if (std::abs(eta) >= absEtaMin_[iEta] && std::abs(eta) < absEtaMax_[iEta]) {
0111       quadEffArea = quadraticEffectiveAreaValues_[iEta];
0112       break;
0113     }
0114   }
0115
0116   return quadEffArea;
0117 }
0118
0119 void EffectiveAreas::printEffectiveAreas() const {
0120   printf("EffectiveAreas: source file %s\n", filename_.c_str());
0121
0122   if (!quadraticEAflag_) {
0123     printf("  eta_min   eta_max    effective area\n");
0124     uint nEtaBins = absEtaMin_.size();
0125
0126     for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0127       printf("  %8.4f    %8.4f   %8.5f\n", absEtaMin_[iEta], absEtaMax_[iEta], effectiveAreaValues_[iEta]);
0128     }
0129   }
0130
0131   else {
0132     printf("  eta_min   eta_max    EA linear term    EA quadratic term\n");
0133     uint nEtaBins = absEtaMin_.size();
0134     for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0135       printf("  %8.4f    %8.4f   %8.5f    %8.5f\n",
0136              absEtaMin_[iEta],
0137              absEtaMax_[iEta],
0138              linearEffectiveAreaValues_[iEta],
0139              quadraticEffectiveAreaValues_[iEta]);
0140     }
0141   }
0142 }
0143
0144 // Basic common sense checks
0145 void EffectiveAreas::checkConsistency() const {
0146   // There should be at least one eta range with one constant
0147
0148   if (!quadraticEAflag_) {
0149     if (effectiveAreaValues_.empty())
0150       throw cms::Exception("EffectiveAreas config failure")
0151           << "found no effective area constants in the file " << filename_ << std::endl;
0152
0153     uint nEtaBins = absEtaMin_.size();
0154
0155     for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0156       // The low limit should be lower than the upper limit
0157       if (!(absEtaMin_[iEta] < absEtaMax_[iEta]))
0158         throw cms::Exception("EffectiveAreas config failure")
0159             << "eta ranges improperly defined (min>max) in the file" << filename_ << std::endl;
0160
0161       // The low limit of the next range should be (near) equal to the
0162       // upper limit of the previous range
0163       if (iEta != nEtaBins - 1)  // don't do the check for the last bin
0164         if (!(absEtaMin_[iEta + 1] - absEtaMax_[iEta] < 0.0001))
0165           throw cms::Exception("EffectiveAreas config failure")
0166               << "eta ranges improperly defined (disjointed) in the file " << filename_ << std::endl;
0167
0168       // The effective area should be non-negative number,
0169       // and should be less than the whole calorimeter area
0170       // eta range -2.5 to 2.5, phi 0 to 2pi => Amax = 5*2*pi ~= 31.4
0171       if (!(effectiveAreaValues_[iEta] >= 0 && effectiveAreaValues_[iEta] < 31.4))
0172         throw cms::Exception("EffectiveAreas config failure")
0173             << "effective area values are too large or negative in the file" << filename_ << std::endl;
0174     }
0175   }
0176
0177   else {
0178     if (linearEffectiveAreaValues_.empty() || quadraticEffectiveAreaValues_.empty())
0179       throw cms::Exception("EffectiveAreas config failure")
0180           << "found no effective area constants (linear and/or quadratic) in the file " << filename_ << std::endl;
0181
0182     uint nEtaBins = absEtaMin_.size();
0183
0184     for (uint iEta = 0; iEta < nEtaBins; iEta++) {
0185       // The low limit should be lower than the upper limit
0186       if (!(absEtaMin_[iEta] < absEtaMax_[iEta]))
0187         throw cms::Exception("EffectiveAreas config failure")
0188             << "eta ranges improperly defined (min>max) in the file" << filename_ << std::endl;
0189
0190       // The low limit of the next range should be (near) equal to the
0191       // upper limit of the previous range
0192       if (iEta != nEtaBins - 1)  // don't do the check for the last bin
0193         if (!(absEtaMin_[iEta + 1] - absEtaMax_[iEta] < 0.0001))
0194           throw cms::Exception("EffectiveAreas config failure")
0195               << "eta ranges improperly defined (disjointed) in the file " << filename_ << std::endl;
0196
0197       // The linear effective area should be non-negative number,
0198       // and should be less than the whole calorimeter area
0199       // eta range -2.5 to 2.5, phi 0 to 2pi => Amax = 5*2*pi ~= 31.4
0200       if (!(linearEffectiveAreaValues_[iEta] >= 0 && linearEffectiveAreaValues_[iEta] < 31.4))
0201         throw cms::Exception("EffectiveAreas config failure")
0202             << "effective area values are too large or negative in the file" << filename_ << std::endl;
0203     }
0204   }
0205 }