MuonIsolation/src/NominalEfficiencyThresholds.cc

0001 #include "RecoMuon/MuonIsolation/src/NominalEfficiencyThresholds.h"
0002
0003 #include <cmath>
0004 #include <iostream>
0005 #include <cstring>
0006 #include <cstdlib>
0007 #include <cstdio>
0008 using namespace muonisolation;
0009 using namespace std;
0010
0011 NominalEfficiencyThresholds::EtaBounds::EtaBounds() {
0012   float BaseEtaBin = 0.087;
0013   theBounds[0] = 0.0;
0014   for (int it = 1; it <= 20; it++)
0015     theBounds[it] = it * BaseEtaBin;
0016   theBounds[21] = 1.83;
0017   theBounds[22] = 1.93;
0018   theBounds[23] = 2.043;
0019   theBounds[24] = 2.172;
0020   theBounds[25] = 2.322;
0021   theBounds[26] = 2.5;
0022   theBounds[27] = 2.65;
0023   theBounds[28] = 3.0;
0024   theBounds[29] = 3.13;
0025   theBounds[30] = 3.305;
0026   theBounds[31] = 3.48;
0027   theBounds[32] = 3.655;
0028 }
0029
0030 int NominalEfficiencyThresholds::EtaBounds::towerFromEta(double eta) const {
0031   int number = 0;
0032   for (int num = 1; num <= NumberOfTowers; num++) {
0033     if (fabs(eta) > theBounds[num - 1]) {
0034       number = num;
0035       if (eta < 0)
0036         number = -num;
0037     }
0038   }
0039   if (fabs(eta) >= theBounds[NumberOfTowers - 1])
0040     number = 0;
0041   return number;
0042 }
0043
0044 vector<double> NominalEfficiencyThresholds::bins() const {
0045   vector<double> result;
0046   for (unsigned int i = 1; i <= EtaBounds::NumberOfTowers; i++)
0047     result.push_back(etabounds(i));
0048   return result;
0049 }
0050
0051 bool NominalEfficiencyThresholds::EfficiencyBin::operator()(const EfficiencyBin &e1, const EfficiencyBin &e2) const {
0052   return e1.eff <= e2.eff_previous;
0053 }
0054
0055 bool NominalEfficiencyThresholds::locless::operator()(const ThresholdLocation &l1, const ThresholdLocation &l2) const {
0056   int itow1 = abs(etabounds.towerFromEta(l1.eta));
0057   int itow2 = abs(etabounds.towerFromEta(l2.eta));
0058   if (itow1 < itow2)
0059     return true;
0060   if (itow1 == itow2 && l1.cone < l2.cone)
0061     return true;
0062   return false;
0063 }
0064
0065 NominalEfficiencyThresholds::NominalEfficiencyThresholds(const string &infile) {
0066   FILE *INFILE;
0067   char buffer[81];
0068   char tag[5];
0069   if ((INFILE = fopen(infile.c_str(), "r")) == nullptr) {
0070     cout << "Cannot open input file " << infile << endl;
0071     return;
0072   }
0073   ThresholdLocation location;
0074   EfficiencyBin eb;
0075   eb.eff = 0;
0076   float thr_val;
0077   while (fgets(buffer, 80, INFILE)) {
0078     sscanf(buffer, "%4s", tag);
0079     if (strcmp(tag, "ver:") == 0) {
0080       cout << " NominalEfficiencyThresholds: " << infile << " comment: " << buffer << endl;
0081       thresholds.clear();
0082     }
0083     if (strcmp(tag, "loc:") == 0) {
0084       sscanf(buffer, "%*s %f %*s %d", &location.eta, &location.cone);
0085       eb.eff = 0.;
0086     }
0087     if (strcmp(tag, "thr:") == 0) {
0088       eb.eff_previous = eb.eff;
0089       sscanf(buffer, "%*s %f %f", &eb.eff, &thr_val);
0090       add(location, ThresholdConstituent(eb, thr_val));
0091     }
0092   }
0093   fclose(INFILE);
0094   cout << "... done" << endl;
0095   //dump();
0096 }
0097
0098 void NominalEfficiencyThresholds::add(ThresholdLocation location, ThresholdConstituent threshold) {
0099   MapType::iterator ploc = thresholds.find(location);
0100   if (ploc == thresholds.end()) {
0101     ThresholdConstituents mt;
0102     mt.insert(threshold);
0103     thresholds[location] = mt;
0104   } else {
0105     //    cout << "insert element ("<<threshold.first.eff<<","
0106     //                            <<threshold.first.eff_previous<<") to map "<<endl;
0107     (*ploc).second.insert(threshold);
0108     //    cout << "new size is:"<< (*ploc).second.size() <<endl;
0109   }
0110 }
0111
0112 void NominalEfficiencyThresholds::dump() {
0113   MapType::iterator ploc;
0114   for (ploc = thresholds.begin(); ploc != thresholds.end(); ploc++) {
0115     cout << "eta: " << (*ploc).first.eta << " icone: " << (*ploc).first.cone << endl;
0116     ThresholdConstituents::iterator it;
0117     for (it = (*ploc).second.begin(); it != (*ploc).second.end(); it++) {
0118       cout << " eff: " << (*it).first.eff << " eff_previous: " << (*it).first.eff_previous
0119            << " cut value: " << (*it).second << endl;
0120     }
0121   }
0122 }
0123
0124 float NominalEfficiencyThresholds::thresholdValueForEfficiency(ThresholdLocation location, float eff_thr) const {
0125   MapType::const_iterator ploc = thresholds.find(location);
0126   if (ploc == thresholds.end()) {
0127     cout << "NominalEfficiencyThresholds: Problem:can't find location in the map :( " << location.eta << " "
0128          << location.cone << " " << eff_thr << endl;
0129     return -1;
0130   }
0131
0132   const float epsilon = 1.e-6;
0133   EfficiencyBin eb = {eff_thr, eff_thr - epsilon};
0134   ThresholdConstituents::const_iterator it = (*ploc).second.find(eb);
0135   if (it == (*ploc).second.end()) {
0136     cout << "NominalEfficiencyThresholds: Problem:can't find threshold in the map :(" << endl;
0137     return -1;
0138   }
0139
0140   return (*it).second;
0141 }