PFClusterTools/src/PFResolutionMap.cc

0001 #include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"
0002 #include <fstream>
0003 #include <sstream>
0004 #include <algorithm>
0005 #include <functional>
0006 #include <stdexcept>
0007 #include <TFile.h>
0008 #include <TMath.h>
0009 #include <vector>
0010 using namespace std;
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 const unsigned PFResolutionMap::lineSize_ = 10000;
0014
0015 PFResolutionMap::PFResolutionMap(const char* name,
0016                                  unsigned nbinseta,
0017                                  double mineta,
0018                                  double maxeta,
0019                                  unsigned nbinse,
0020                                  double mine,
0021                                  double maxe,
0022                                  double value)
0023     : TH2D(name, name, nbinseta, mineta, maxeta, nbinse, mine, maxe) {
0024   GetXaxis()->SetTitle("#eta");
0025   GetYaxis()->SetTitle("E");
0026
0027   if (value > 0) {
0028     for (int ie = 1; ie <= GetNbinsY(); ie++) {
0029       for (int ieta = 1; ieta <= GetNbinsX(); ieta++) {
0030         SetBinContent(ieta, ie, value);
0031       }
0032     }
0033   }
0034   // cout<<"creating resolution map "<<endl;
0035   // Print("all");
0036 }
0037
0038 PFResolutionMap::PFResolutionMap(const char* name, const char* mapfile) {
0039   SetTitle(name);
0040   GetXaxis()->SetTitle("#eta");
0041   GetYaxis()->SetTitle("E");
0042   if (!ReadMapFile(mapfile)) {
0043     string err = "PFResolutionMap::PFResolutionMap : cannot read file ";
0044     err += mapfile;
0045     throw invalid_argument(err);
0046   }
0047 }
0048
0049 bool PFResolutionMap::WriteMapFile(const char* mapfile) {
0050   //  assert(fData.size() == fNbinsEta*fNbinsE);
0051
0052   // open the file
0053
0054   std::ofstream outf(mapfile);
0055   if (!outf.good()) {
0056     edm::LogWarning("PFResolutionMap::Write") << " : cannot open file " << mapfile;
0057     return false;
0058   }
0059
0060   outf << (*this) << endl;
0061   if (!outf.good()) {
0062     edm::LogError("PFResolutionMap::Write") << " : corrupted file " << mapfile;
0063     return false;
0064   } else {
0065     mapFile_ = mapfile;
0066     return true;
0067   }
0068 }
0069
0070 bool PFResolutionMap::ReadMapFile(const char* mapfile) {
0071   // open the file
0072
0073   std::ifstream inf(mapfile);
0074   if (!inf.good()) {
0075     // cout<<"PFResolutionMap::Read : cannot open file "<<mapfile<<endl;
0076     return false;
0077   }
0078
0079   // first data describes the map
0080
0081   int nbinseta = 0;
0082   double mineta = 0;
0083   double maxeta = 0;
0084
0085   int nbinse = 0;
0086   double mine = 0;
0087   double maxe = 0;
0088
0089   inf >> nbinseta;
0090   inf >> mineta;
0091   inf >> maxeta;
0092
0093   inf >> nbinse;
0094   inf >> mine;
0095   inf >> maxe;
0096
0097   SetBins(nbinseta, mineta, maxeta, nbinse, mine, maxe);
0098
0099   char s[lineSize_];
0100   int pos = inf.tellg();
0101
0102   // parse map data
0103   int i = 1;
0104   do {
0105     inf.seekg(pos);
0106     inf.getline(s, lineSize_);
0107
0108     pos = inf.tellg();
0109
0110     if (string(s).empty()) {
0111       continue;  // remove empty lines
0112     }
0113
0114     istringstream lin(s);
0115
0116     double dataw;
0117     int j = 1;
0118     do {
0119       lin >> dataw;
0120       SetBinContent(j, i, dataw);
0121       j++;
0122     } while (lin.good());
0123     i++;
0124   } while (inf.good());
0125
0126   if (inf.eof()) {
0127     mapFile_ = mapfile;
0128     return true;
0129   } else {
0130     // cout<<"error"<<endl;
0131     return false;
0132   }
0133
0134   mapFile_ = mapfile;
0135   return true;
0136 }
0137
0138 double PFResolutionMap::getRes(double eta, double phi, double e, int MapEta) {
0139   constexpr double fMinEta = -2.95;
0140   constexpr double fMaxEta = 2.95;
0141   constexpr double fMinE = 0;
0142   constexpr double fMaxE = 100;
0143
0144   if (eta < fMinEta)
0145     eta = fMinEta + 0.001;
0146   if (eta > fMaxEta)
0147     eta = fMaxEta - 0.001;
0148
0149   if (e < fMinE)
0150     e = fMinE + 0.001;
0151   if (e > fMaxE)
0152     e = fMaxE - 0.001;
0153
0154   unsigned bin = FindBin(TMath::Abs(eta), e);
0155
0156   double res = GetBinContent(bin);
0157   if (MapEta > -1) {
0158     if ((eta < 1.48) && IsInAPhiCrack(phi, eta)) {
0159       if (MapEta == 1)
0160         res *= 1.88;
0161       else
0162         res *= 1.65;
0163     }
0164   }
0165   return res;
0166 }
0167
0168 int PFResolutionMap::FindBin(double eta, double e, double z) {
0169   if (e >= GetYaxis()->GetXmax())
0170     e = GetYaxis()->GetXmax() - 0.001;
0171
0172   return TH2D::FindBin(eta, e);
0173 }
0174
0175 ostream& operator<<(ostream& outf, const PFResolutionMap& rm) {
0176   if (!outf.good())
0177     return outf;
0178
0179   // first data describes the map
0180   outf << rm.GetNbinsX() << endl;
0181   outf << rm.GetXaxis()->GetXmin() << endl;
0182   outf << rm.GetXaxis()->GetXmax() << endl;
0183
0184   outf << rm.GetNbinsY() << endl;
0185   outf << rm.GetYaxis()->GetXmin() << endl;
0186   outf << rm.GetYaxis()->GetXmax() << endl;
0187
0188   for (int ie = 1; ie <= rm.GetNbinsY(); ie++) {
0189     for (int ieta = 1; ieta <= rm.GetNbinsX(); ieta++) {
0190       outf << rm.GetBinContent(ieta, ie) << "\t";
0191     }
0192     outf << endl;
0193   }
0194
0195   return outf;
0196 }
0197
0198 bool PFResolutionMap::IsInAPhiCrack(double phi, double eta) {
0199   double dminPhi = dCrackPhi(phi, eta);
0200   bool Is = (TMath::Abs(dminPhi) < 0.005);
0201   return Is;
0202 }
0203
0204 //useful to compute the signed distance to the closest crack in the barrel
0205 double PFResolutionMap::minimum(double a, double b) {
0206   if (TMath::Abs(b) < TMath::Abs(a))
0207     a = b;
0208   return a;
0209 }
0210
0211 #ifndef M_PI
0212 #define M_PI 3.14159265358979323846
0213 #endif
0214
0215 //compute the unsigned distance to the closest phi-crack in the barrel
0216 double PFResolutionMap::dCrackPhi(double phi, double eta) {
0217   constexpr double pi = M_PI;  // 3.14159265358979323846;
0218   constexpr double twopi = 2. * pi;
0219   constexpr double twopiO18 = pi / 9;
0220
0221   //Location of the 18 phi-cracks
0222   constexpr double c0 = 2.97025;
0223   constexpr std::array<double, 18> cPhi{{c0,
0224                                          c0 - twopiO18,
0225                                          c0 - 2 * twopiO18,
0226                                          c0 - 3 * twopiO18,
0227                                          c0 - 4 * twopiO18,
0228                                          c0 - 5 * twopiO18,
0229                                          c0 - 6 * twopiO18,
0230                                          c0 - 7 * twopiO18,
0231                                          c0 - 8 * twopiO18,
0232                                          c0 - 9 * twopiO18,
0233                                          c0 - 10 * twopiO18,
0234                                          c0 - 11 * twopiO18,
0235                                          c0 - 12 * twopiO18,
0236                                          c0 - 13 * twopiO18,
0237                                          c0 - 14 * twopiO18,
0238                                          c0 - 15 * twopiO18,
0239                                          c0 - 16 * twopiO18,
0240                                          c0 - 17 * twopiO18}};
0241
0242   //Shift of this location if eta<0
0243   constexpr double delta_cPhi = 0.00638;
0244
0245   double m;  //the result
0246
0247   //the location is shifted
0248   if (eta < 0) {
0249     phi += delta_cPhi;
0250     if (phi > pi)
0251       phi -= twopi;
0252   }
0253   if (phi >= -pi && phi <= pi) {
0254     //the problem of the extrema
0255     if (phi < cPhi[17] || phi >= cPhi[0]) {
0256       if (phi < 0)
0257         phi += twopi;
0258       m = minimum(phi - cPhi[0], phi - cPhi[17] - twopi);
0259     }
0260
0261     //between these extrema...
0262     else {
0263       bool OK = false;
0264       unsigned i = 16;
0265       while (!OK) {
0266         if (phi < cPhi[i]) {
0267           m = minimum(phi - cPhi[i + 1], phi - cPhi[i]);
0268           OK = true;
0269         } else
0270           i -= 1;
0271       }
0272     }
0273   } else {
0274     m = 0.;  //if there is a problem, we assum that we are in a crack
0275     edm::LogWarning("PFResolutionMap:Problem") << "Problem in dminphi";
0276   }
0277   if (eta < 0)
0278     m = -m;  //because of the disymetry
0279   return m;
0280 }