InvMatrixUtils.cc

Go to the documentation of this file.

 
#include "Calibration/Tools/interface/InvMatrixUtils.h"
#include "Calibration/Tools/interface/InvMatrixCommonDefs.h"
#include "TStyle.h"
#include "TROOT.h"
#include "CLHEP/Geometry/Point3D.h"
//#include "ConfigParser.h"
#include "Calibration/Tools/interface/matrixSaver.h"
 
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <cassert>
 
void setStyle() {
  gROOT->SetStyle("Plain");
  gStyle->SetTextSize(0.5);
  //gStyle->SetOptStat (1111111) ;
  gStyle->SetOptStat(0);
  //gStyle->SetOptFit (1111) ;
  gStyle->SetOptFit(0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleX(0.08);
  gStyle->SetTitleY(0.97);
  gStyle->SetPalette(1, nullptr);
  gStyle->SetStatColor(0);
  gStyle->SetFrameFillStyle(0);
  gStyle->SetFrameFillColor(0);
  return;
}
 
// -----------------------------------------------------------------
 
TCanvas *getGlobalCanvas(std::string name) {
  setStyle();
  TCanvas *globalCanvas = static_cast<TCanvas *>(gROOT->FindObject(name.c_str()));
  if (globalCanvas) {
    globalCanvas->Clear();
    globalCanvas->UseCurrentStyle();
    globalCanvas->SetWindowSize(700, 600);
 
  } else {
    globalCanvas = new TCanvas(name.c_str(), name.c_str(), 700, 600);
  }
  return globalCanvas;
}
 
// -----------------------------------------------------------------
 
TFile *getGlobalTFile(std::string name) {
  //    std::cout << "writing " << name << std::endl ;
  //    setStyle () ;
  TFile *globalTFile = (TFile *)gROOT->FindObject(name.c_str());
  if (!globalTFile) {
    //        std::cout << "does not exist. creating it " << std::endl;
    globalTFile = new TFile(name.c_str(), "RECREATE");
  }
 
  return globalTFile;
}
 
// -----------------------------------------------------------------
 
int saveGlobalTFile(std::string name) {
  TFile *globalTFile = static_cast<TFile *>(gROOT->FindObject(name.c_str()));
  if (!globalTFile)
    return 1;
  globalTFile->Write();
  globalTFile->Close();
  delete globalTFile;
  return 0;
}
 
// -----------------------------------------------------------------
 
CLHEP::HepMatrix *getSavedMatrix(const std::string &name) {
  matrixSaver reader;
  CLHEP::HepMatrix *savedMatrix;
  if (reader.touch(name)) {
    savedMatrix = static_cast<CLHEP::HepMatrix *>(reader.getMatrix(name));
  } else {
    savedMatrix = new CLHEP::HepMatrix(SCMaxEta, SCMaxPhi, 0);
  }
 
  return savedMatrix;
}
 
//========================================================================
//da usare
//HepGeom::Point3D<double> TBimpactPoint = TBposition (amplitude,m_beamEnergy) ;
 
/* 
dove trovare il codice di Chiara in CMSSW, per la ricostruzione
della posizione:
 
http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/RecoTBCalo/EcalTBAnalysisCoreTools/src/TBPositionCalc.cc?rev=1.1&cvsroot=CMSSW&content-type=text/vnd.viewcvs-markup
 
*/
 
HepGeom::Point3D<Float_t> TBposition(const Float_t amplit[7][7],
                                     const Float_t beamEne,
                                     const Float_t w0,
                                     const Float_t x0,
                                     const Float_t a0,
                                     const Float_t sideX,  // crystal geometry, in mm
                                     const Float_t sideY) {
  // variables
  Float_t caloX = 0.;
  Float_t caloY = 0.;
  Float_t sumWeight = 0.;
  Float_t depth = x0 * (log(beamEne) + a0);  // shower depthh, in mm
  Float_t sin3 = 0.052335956;                // sin (3 degrees) , sin3 = sin(3.*3.141592654/180.)
 
  Float_t invE3x3 = 1. / get3x3(amplit);
 
  // loop over 3x3 crystals
  for (int eta = 2; eta <= 4; eta++) {
    for (int phi = 2; phi <= 4; phi++) {
      Float_t weight = log(amplit[eta][phi] * invE3x3) + w0;
      if (weight > 0) {
        caloX += (eta - 3) * sideX * weight;
        caloY -= (phi - 3) * sideY * weight;
        sumWeight += weight;
      }
    }
  }
 
  caloX /= sumWeight;
  caloY /= sumWeight;
 
  // correction for shower depthh
  caloX -= depth * sin3;
  caloY -= depth * sin3;
 
  // FIXME the z set to zero
  HepGeom::Point3D<Float_t> TBposition(caloX, caloY, 0);
 
  return TBposition;
}
 
// -----------------------------------------------------------------
 
double get5x5(const Float_t energy[7][7]) {
  double total = 0.;
 
  for (int eta = 1; eta < 6; ++eta)
    for (int phi = 1; phi < 6; ++phi)
      total += energy[eta][phi];
 
  return total;
}
 
// -----------------------------------------------------------------
 
double get3x3(const Float_t energy[7][7]) {
  double total = 0.;
 
  for (int eta = 2; eta < 5; ++eta)
    for (int phi = 2; phi < 5; ++phi)
      total += energy[eta][phi];
 
  return total;
}
 
// -----------------------------------------------------------------
 
/* int parseConfigFile (const TString& config)
{
    if (gConfigParser) return 1 ;
    
    std::cout << "parsing "
        << config << " file"
        << std::endl ;
    
    
    gConfigParser = new ConfigParser () ;
    if ( !(gConfigParser->init(config)) )
        {
        std::cout << "Analysis::parseConfigFile: Could not open configuration file "
        << config << std::endl;
        perror ("Analysis::parseConfigFile: ") ;
        exit (-1) ;
        }
    gConfigParser->print () ;
    return 0 ;
}*/
 
// -----------------------------------------------------------------
 
// per un certo eta, il cristallo puo' essere qualsiasi intero
// Calibrationtra xmin e xmax
// lo posso fissare solo sfruttando anche phi
int xtalFromEtaPhi(const int &myEta, const int &myPhi) {
  int xMin = 20 * myEta + 1;
  int xMax = 20 * (myEta + 1) + 1;
 
  int myCryst = 999999;
 
  for (int x = xMin; x < xMax; x++) {
    if (phiFromXtal(x) == myPhi)
      myCryst = x;
  }
  return myCryst;
}
 
// -----------------------------------------------------------------
 
int xtalFromiEtaiPhi(const int &iEta, const int &iPhi) {
  assert(iEta >= 1);
  assert(iEta <= 85);
  assert(iPhi >= 1);
  assert(iPhi <= 20);
  return 20 * (iEta - 1) + 21 - iPhi;
}
 
// -----------------------------------------------------------------
 
int etaFromXtal(const int &xtal) {
  //  return floor (static_cast<double> ((xtal-1) / 20)) ;
  return int(floor((xtal - 1) / 20));
}
 
// -----------------------------------------------------------------
 
int phiFromXtal(const int &xtal) {
  int phi = (xtal - 1) - 20 * etaFromXtal(xtal);
  return (20 - phi - 1);
}
 
// -----------------------------------------------------------------
 
int ietaFromXtal(const int &xtal) { return etaFromXtal(xtal) + 1; }
 
// -----------------------------------------------------------------
 
int iphiFromXtal(const int &xtal) { return phiFromXtal(xtal) + 1; }
 
// -----------------------------------------------------------------
 
int extract(std::vector<int> *output, const std::string &dati) {
  std::ifstream _dati(dati.c_str());
  // loop over the file
  while (!_dati.eof()) {
    // get the line
    std::string dataline;
    do {
      getline(_dati, dataline, '\n');
    } while (*dataline.begin() == '#');
    std::stringstream linea(dataline);
    // loop over the line
    while (!linea.eof()) {
      int buffer = -1;
      linea >> buffer;
      if (buffer != -1)
        output->push_back(buffer);
    }  // loop over the line
  }    // loop over the file
  return output->size();
}
 
// -----------------------------------------------------------------
 
// FIXME questi eta, phi sono quelli della matrice CLHEP,
// FIXME non quelli del super-modulo, giusto?
int writeCalibTxt(const CLHEP::HepMatrix &AmplitudeMatrix,
                  const CLHEP::HepMatrix &SigmaMatrix,
                  const CLHEP::HepMatrix &StatisticMatrix,
                  std::string fileName) {
  // look for the reference crystal
  double reference = 0.;
  for (int eta = 0; eta < SCMaxEta; ++eta)
    for (int phi = 0; phi < SCMaxPhi; ++phi) {
      if (AmplitudeMatrix[eta][phi] && SigmaMatrix[eta][phi] < 100 /*FIXME sigmaCut*/) {
        reference = AmplitudeMatrix[eta][phi];
        std::cout << "[InvMatrixUtils][writeCalibTxt] reference crystal: "
                  << "(" << eta << "," << phi << ") -> " << reference << "\n";
        break;
      }
    }
  if (!reference) {
    std::cerr << "ERROR: no calibration coefficients found" << std::endl;
    return 1;
  }
 
  // open the file for output
  std::ofstream txt_outfile;
  txt_outfile.open(fileName.c_str());
  txt_outfile << "# xtal\tcoeff\tsigma\tevt\tisGood\n";
 
  // loop over the crystals
  for (int eta = 0; eta < SCMaxEta; ++eta)
    for (int phi = 0; phi < SCMaxPhi; ++phi) {
      int isGood = 1;
      if (AmplitudeMatrix[eta][phi] == 0)
        isGood = 0;
      if (SigmaMatrix[eta][phi] > 100 /*FIXME sigmaCut*/)
        isGood = 0;
      txt_outfile << xtalFromEtaPhi(eta, phi) << "\t" << AmplitudeMatrix[eta][phi] / reference << "\t"
                  << SigmaMatrix[eta][phi] << "\t" << StatisticMatrix[eta][phi] << "\t" << isGood << "\n";
    }
 
  // save and close the file
  txt_outfile.close();
  return 0;
}
 
// -----------------------------------------------------------------
 
int writeCMSSWCoeff(const CLHEP::HepMatrix &amplMatrix,
                    double calibThres,
                    float ERef,
                    const CLHEP::HepMatrix &sigmaMatrix,
                    const CLHEP::HepMatrix &statisticMatrix,
                    std::string fileName,
                    std::string genTag,
                    std::string method,
                    std::string version,
                    std::string type) {
  // open the file for output
  std::ofstream txt_outfile;
  txt_outfile.open(fileName.c_str());
  txt_outfile << "1\n";   // super-module number
  txt_outfile << "-1\n";  // number of events
  txt_outfile << genTag << "\n";
  txt_outfile << method << "\n";
  txt_outfile << version << "\n";
  txt_outfile << type << "\n";
 
  double reference = ERef;
 
  // loop over crystals
  for (int eta = 0; eta < SCMaxEta; ++eta)
    for (int phi = 0; phi < SCMaxPhi; ++phi) {
      if (amplMatrix[eta][phi] <= calibThres)
        txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
      else
        txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << reference / amplMatrix[eta][phi] << "\t"
                    << sigmaMatrix[eta][phi] << "\t" << statisticMatrix[eta][phi] << "\t" << 1 << "\n";
    }  // loop over crystals
 
  // save and close the file
  txt_outfile.close();
  return 0;
}
 
// -----------------------------------------------------------------
 
int writeCMSSWCoeff(const CLHEP::HepMatrix &amplMatrix,
                    double calibThres,
                    int etaRef,
                    int phiRef,
                    const CLHEP::HepMatrix &sigmaMatrix,
                    const CLHEP::HepMatrix &statisticMatrix,
                    std::string fileName,
                    std::string genTag,
                    std::string method,
                    std::string version,
                    std::string type) {
  // open the file for output
  std::ofstream txt_outfile;
  txt_outfile.open(fileName.c_str());
  txt_outfile << "1\n";   // super-module number
  txt_outfile << "-1\n";  // number of events
  txt_outfile << genTag << "\n";
  txt_outfile << method << "\n";
  txt_outfile << version << "\n";
  txt_outfile << type << "\n";
 
  if (amplMatrix[etaRef - 1][phiRef - 1] == 0) {
    std::cerr << "The reference crystal: (" << etaRef << "," << phiRef << ") is out of range\n";
    return 1;
  }
  double reference = amplMatrix[etaRef - 1][phiRef - 1];
 
  // loop over crystals
  for (int eta = 0; eta < SCMaxEta; ++eta)
    for (int phi = 0; phi < SCMaxPhi; ++phi) {
      if (amplMatrix[eta][phi] <= calibThres)
        txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
      else
        txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << reference / amplMatrix[eta][phi] << "\t"
                    << sigmaMatrix[eta][phi] << "\t" << statisticMatrix[eta][phi] << "\t" << 1 << "\n";
    }  // loop over crystals
 
  // save and close the file
  txt_outfile.close();
  return 0;
}
 
// -----------------------------------------------------------------
 
int translateCoeff(const CLHEP::HepMatrix &calibcoeff,
                   const CLHEP::HepMatrix &sigmaMatrix,
                   const CLHEP::HepMatrix &statisticMatrix,
                   std::string SMnumber,
                   double calibThres,
                   std::string fileName,
                   std::string genTag,
                   std::string method,
                   std::string version,
                   std::string type) {
  // open the file for output
  std::ofstream txt_outfile;
  txt_outfile.open(fileName.c_str());
  txt_outfile << SMnumber << "\n";  // super-module number
  txt_outfile << "-1\n";            // number of events
  txt_outfile << genTag << "\n";
  txt_outfile << method << "\n";
  txt_outfile << version << "\n";
  txt_outfile << type << "\n";
 
  // loop over crystals
  for (int eta = 0; eta < SCMaxEta; ++eta)
    for (int phi = 0; phi < SCMaxPhi; ++phi) {
      if (calibcoeff[eta][phi] < calibThres) {
        txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
        std::cout << "[translateCoefff][" << SMnumber << "]\t WARNING crystal " << xtalFromiEtaiPhi(eta + 1, phi + 1)
                  << " calib coeff below threshold: "
                  << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
      } else
        txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << calibcoeff[eta][phi] << "\t"
                    << sigmaMatrix[eta][phi] << "\t" << statisticMatrix[eta][phi] << "\t" << 1 << "\n";
    }  // loop over crystals
 
  // save and close the file
  txt_outfile.close();
  return 0;
}
 
// -----------------------------------------------------------------
 
int readCMSSWcoeff(CLHEP::HepMatrix &calibcoeff, const std::string &inputFileName, double defaultVal) {
  std::ifstream CMSSWfile;
  CMSSWfile.open(inputFileName.c_str());
  std::string buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  while (!CMSSWfile.eof()) {
    int xtalnum;
    CMSSWfile >> xtalnum;
    double coeff;
    CMSSWfile >> coeff;
    double buffer;
    CMSSWfile >> buffer;
    int good;
    CMSSWfile >> good;
    CMSSWfile >> good;
    if (!good)
      coeff = defaultVal;  //FIXME 0 o 1?
    calibcoeff[etaFromXtal(xtalnum)][phiFromXtal(xtalnum)] = coeff;
  }
  return 0;
}
 
// -----------------------------------------------------------------
 
int readCMSSWcoeffForComparison(CLHEP::HepMatrix &calibcoeff, const std::string &inputFileName) {
  std::ifstream CMSSWfile;
  CMSSWfile.open(inputFileName.c_str());
  std::string buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  CMSSWfile >> buffer;
  while (!CMSSWfile.eof()) {
    int xtalnum;
    CMSSWfile >> xtalnum;
    double coeff;
    CMSSWfile >> coeff;
    double buffer;
    CMSSWfile >> buffer;
    int good;
    CMSSWfile >> good;
    CMSSWfile >> good;
    if (!good)
      coeff = 0.;  //FIXME 0 o 1?
    calibcoeff[etaFromXtal(xtalnum)][phiFromXtal(xtalnum)] = coeff;
  }
  return 0;
}
 
// -----------------------------------------------------------------
 
TH1D *smartProfile(TH2F *strip, double width) {
  TProfile *stripProfile = strip->ProfileX();
 
  // (from FitSlices of TH2.h)
 
  double xmin = stripProfile->GetXaxis()->GetXmin();
  double xmax = stripProfile->GetXaxis()->GetXmax();
  int profileBins = stripProfile->GetNbinsX();
 
  std::string name = strip->GetName();
  name += "_smart";
  TH1D *prof = new TH1D(name.c_str(), strip->GetTitle(), profileBins, xmin, xmax);
 
  int cut = 0;  // minimum number of entries per fitted bin
  int nbins = strip->GetXaxis()->GetNbins();
  int binmin = 1;
  int ngroup = 1;  // bins per step
  int binmax = nbins;
 
  // loop over the strip bins
  for (int bin = binmin; bin <= binmax; bin += ngroup) {
    TH1D *hpy = strip->ProjectionY("_temp", bin, bin + ngroup - 1, "e");
    if (hpy == nullptr)
      continue;
    int nentries = Int_t(hpy->GetEntries());
    if (nentries == 0 || nentries < cut) {
      delete hpy;
      continue;
    }
 
    Int_t biny = bin + ngroup / 2;
 
    hpy->GetXaxis()->SetRangeUser(hpy->GetMean() - width * hpy->GetRMS(), hpy->GetMean() + width * hpy->GetRMS());
    prof->Fill(strip->GetXaxis()->GetBinCenter(biny), hpy->GetMean());
    prof->SetBinError(biny, hpy->GetRMS());
 
    delete hpy;
  }  // loop over the bins
 
  delete stripProfile;
  return prof;
}
 
// -----------------------------------------------------------------
 
TH1D *smartGausProfile(TH2F *strip, double width) {
  TProfile *stripProfile = strip->ProfileX();
 
  // (from FitSlices of TH2.h)
 
  double xmin = stripProfile->GetXaxis()->GetXmin();
  double xmax = stripProfile->GetXaxis()->GetXmax();
  int profileBins = stripProfile->GetNbinsX();
 
  std::string name = strip->GetName();
  name += "_smartGaus";
  TH1D *prof = new TH1D(name.c_str(), strip->GetTitle(), profileBins, xmin, xmax);
 
  int cut = 0;  // minimum number of entries per fitted bin
  int nbins = strip->GetXaxis()->GetNbins();
  int binmin = 1;
  int ngroup = 1;  // bins per step
  int binmax = nbins;
 
  // loop over the strip bins
  for (int bin = binmin; bin <= binmax; bin += ngroup) {
    TH1D *hpy = strip->ProjectionY("_temp", bin, bin + ngroup - 1, "e");
    if (hpy == nullptr)
      continue;
    int nentries = Int_t(hpy->GetEntries());
    if (nentries == 0 || nentries < cut) {
      delete hpy;
      continue;
    }
 
    Int_t biny = bin + ngroup / 2;
 
    TF1 *gaussian =
        new TF1("gaussian", "gaus", hpy->GetMean() - width * hpy->GetRMS(), hpy->GetMean() + width * hpy->GetRMS());
    gaussian->SetParameter(1, hpy->GetMean());
    gaussian->SetParameter(2, hpy->GetRMS());
    hpy->Fit("gaussian", "RQL");
 
    hpy->GetXaxis()->SetRangeUser(hpy->GetMean() - width * hpy->GetRMS(), hpy->GetMean() + width * hpy->GetRMS());
    prof->Fill(strip->GetXaxis()->GetBinCenter(biny), gaussian->GetParameter(1));
    prof->SetBinError(biny, gaussian->GetParameter(2));
 
    delete gaussian;
    delete hpy;
  }  // loop over the bins
 
  delete stripProfile;
  return prof;
}
 
// -----------------------------------------------------------------
 
TH1D *smartError(TH1D *strip) {
  double xmin = strip->GetXaxis()->GetXmin();
  double xmax = strip->GetXaxis()->GetXmax();
  int stripsBins = strip->GetNbinsX();
 
  std::string name = strip->GetName();
  name += "_error";
  TH1D *error = new TH1D(name.c_str(), strip->GetTitle(), stripsBins, xmin, xmax);
 
  int binmin = 1;
  int ngroup = 1;  // bins per step
  int binmax = stripsBins;
  for (int bin = binmin; bin <= binmax; bin += ngroup) {
    double dummyError = strip->GetBinError(bin);
    error->SetBinContent(bin, dummyError);
  }
  return error;
}
 
// -----------------------------------------------------------------
 
double effectiveSigma(TH1F &histogram, int vSteps) {
  double totInt = histogram.Integral();
  int maxBin = histogram.GetMaximumBin();
  int maxBinVal = int(histogram.GetBinContent(maxBin));
  int totBins = histogram.GetNbinsX();
  double area = totInt;
  double threshold = 0;
  double vStep = maxBinVal / vSteps;
  int leftBin = 1;
  int rightBin = totBins - 1;
  //loop over the vertical range
  while (area / totInt > 0.683) {
    threshold += vStep;
    // loop toward the left
    for (int back = maxBin; back > 0; --back) {
      if (histogram.GetBinContent(back) < threshold) {
        leftBin = back;
        break;
      }
    }  // loop toward the left
 
    // loop toward the right
    for (int fwd = maxBin; fwd < totBins; ++fwd) {
      if (histogram.GetBinContent(fwd) < threshold) {
        rightBin = fwd;
        break;
      }
    }  // loop toward the right
    area = histogram.Integral(leftBin, rightBin);
  }  //loop over the vertical range
 
  histogram.GetXaxis()->SetRange(leftBin, rightBin);
  // double sigmaEff = histogram.GetRMS () ;
  double halfWidthRange = 0.5 * (histogram.GetBinCenter(rightBin) - histogram.GetBinCenter(leftBin));
  return halfWidthRange;
}
 
// -----------------------------------------------------------------
 
std::pair<int, int> findSupport(TH1F &histogram, double thres) {
  int totBins = histogram.GetNbinsX();
  if (thres >= histogram.GetMaximum())
    return std::pair<int, int>(0, totBins);
 
  int leftBin = totBins - 1;
  // search from left for the minimum
  for (int bin = 1; bin < totBins; ++bin) {
    if (histogram.GetBinContent(bin) > thres) {
      leftBin = bin;
      break;
    }
  }  // search from left for the minimum
  int rightBin = 1;
  // search from right for the maximum
  for (int bin = totBins - 1; bin > 0; --bin) {
    if (histogram.GetBinContent(bin) > thres) {
      rightBin = bin;
      break;
    }
  }  // search from right for the maximum
  return std::pair<int, int>(leftBin, rightBin);
}
 
// -----------------------------------------------------------------
 
void mtrTransfer(double output[SCMaxEta][SCMaxPhi], CLHEP::HepMatrix *input, double Default) {
  for (int eta = 0; eta < SCMaxEta; ++eta)
    for (int phi = 0; phi < SCMaxPhi; ++phi) {
      if ((*input)[eta][phi])
        output[eta][phi] = (*input)[eta][phi];
      else
        output[eta][phi] = Default;
    }
  return;
}
 
// -----------------------------------------------------------------
 
double etaCorrE1E25(int eta) {
  double p0 = 0.807883;
  double p1 = 0.000182551;
  double p2 = -5.76961e-06;
  double p3 = 7.41903e-08;
  double p4 = -2.25384e-10;
 
  double corr;
  if (eta < 6)
    corr = p0;
  else
    corr = p0 + p1 * eta + p2 * eta * eta + p3 * eta * eta * eta + p4 * eta * eta * eta * eta;
  return corr / p0;
}
// -----------------------------------------------------------------
 
double etaCorrE1E49(int eta) {
  double p0 = 0.799895;
  double p1 = 0.000235487;
  double p2 = -8.26496e-06;
  double p3 = 1.21564e-07;
  double p4 = -4.83286e-10;
 
  double corr;
  if (eta < 8)
    corr = p0;
  else
    corr = p0 + p1 * eta + p2 * eta * eta + p3 * eta * eta * eta + p4 * eta * eta * eta * eta;
  return corr / p0;
}
// -----------------------------------------------------------------
 
double etaCorrE1E9(int eta) {
  if (eta < 4)
    return 1.0;
  // grazie Paolo
  double p0 = 0.834629;
  double p1 = 0.00015254;
  double p2 = -4.91784e-06;
  double p3 = 6.54652e-08;
  double p4 = -2.4894e-10;
 
  double corr;
  if (eta < 6)
    corr = p0;
  else
    corr = p0 + p1 * eta + p2 * eta * eta + p3 * eta * eta * eta + p4 * eta * eta * eta * eta;
  return corr / p0;
}