EmcRecParameter.cxx

Go to the documentation of this file.

//
//  Parameters for Emc Reconstruction
//  
//  No Parameter is allowed to be hard coded into code!
//  
//  Created by Zhe Wang, May 31, 2004
//
#include "EmcRec/EmcRecParameter.h"
#include "TGraph2DErrors.h"
#include <fstream>
#include <sstream>
#include <assert.h>
#include <cstdlib> 
 
//Initialize static data member
EmcRecParameter* EmcRecParameter::fpInstance=0;
pthread_mutex_t EmcRecParameter::m_pthread_lock = PTHREAD_MUTEX_INITIALIZER;
 
// Constructors and destructors 
EmcRecParameter::EmcRecParameter()
{
  string paraPath = getenv("EMCRECROOT");
  if(paraPath==""){} //cout<<"EmcRec environment not set!";
  string paraPath1(paraPath);
  string paraPath2(paraPath);
  paraPath += "/share/EmcRecPara.dat";
  //cout<<paraPath<<endl;
  ifstream in;
  in.open(paraPath.c_str());
  assert(in);
 
  const int maxCharOfOneLine=255;
  char temp[maxCharOfOneLine];
  int inputNo=0;
  while(in.peek()!=EOF) {
    in.getline(temp,maxCharOfOneLine);
    if(temp[0]=='#') continue;
    inputNo++;
    switch(inputNo) {
      case 1:
        istringstream(temp)>>fElectronicsNoiseLevel>>fEThresholdSeed>>fEThresholdCluster>>fLogPosOffset;
        break;
      case 2:
        istringstream(temp)>>fMoliereRadius>>fLateralProfile;
        break;
      case 3:
        istringstream(temp)>>eCorr[0]>>eCorr[1]>>eCorr[2]>>eCorr[3];
        break;
      case 4:
        istringstream(temp)>>sigE[0]>>sigE[1]>>sigE[2]
          >>sigTheta[0]>>sigTheta[1]>>sigPhi[0]>>sigPhi[1];
        break;
      case 5:
        istringstream(temp)>>hitNb[0]>>hitNb[1]>>hitNb[2];
        break;
      case 6:
        istringstream(temp)>>elecBias[0]>>elecBias[1]>>elecBias[2]
          >>elecBias[3]>>elecBias[4];
        break;
      case 7:
        istringstream(temp)>>smCut[0]>>smCut[1]>>smCut[2]>>smCut[3];
        break;
      default:
        break;
    }
  }
  in.close();
  
  paraPath1 += "/share/Peak1.843_10.12calib.dat";
  ifstream in1;
  in1.open(paraPath1.c_str());
  assert(in1);
  int ntheta;
  double sigma,sigmaerr,peakerr;
  for(int i=0;i<56;i++) {
    in1>>ntheta>>sigma>>sigmaerr>>peak[i]>>peakerr;
  }
  in1.close();
 
  // Read energy correction parameters from PhotonCor/McCor
  paraPath2 += "/share/evset.txt";
  ifstream in2;
  in2.open(paraPath2.c_str());
  assert(in2);
  double energy,thetaid,peak1,peakerr1,res,reserr;
  //double ep[18]={0.03,0.04,0.05,0.075,0.1,0.125,0.15,0.2,0.25,0.3,0.4,0.5,0.75,1.0,1.25,1.5,1.75,2.0};
  dt = new TGraph2DErrors();
  dtErr = new TGraph2DErrors();
  for(int i=0;i<504;i++){
    in2>>energy;
    in2>>thetaid;
    in2>>peak1;
    in2>>peakerr1;
    in2>>res;
    in2>>reserr;
    dt->SetPoint(i,energy,thetaid,peak1);
    dt->SetPointError(i,0,0,peakerr1);
    dtErr->SetPoint(i,energy,thetaid,res);
    dtErr->SetPointError(i,0,0,reserr);
  }
  in2.close();
 
}
 
 
EmcRecParameter::~EmcRecParameter()
{
  //   cout<<"deconstructed"<<endl;
  delete dt;
  delete dtErr;
}
 
// static method
//Access to an instance
EmcRecParameter& EmcRecParameter::GetInstance()
{
  if(!Exist()) {
    fpInstance=new EmcRecParameter;
  }
  return *fpInstance;
}
 
bool EmcRecParameter::Exist()
{
  return fpInstance!=0;
}
 
void EmcRecParameter::Kill()
{
  if(Exist()) {
    delete fpInstance;
    fpInstance=0;
  }
}
 
//access to each parameter
double EmcRecParameter::ElectronicsNoiseLevel() const
{
  return fElectronicsNoiseLevel;
}
 
double EmcRecParameter::EThresholdSeed() const
{
  return fEThresholdSeed;
}
 
double EmcRecParameter::EThresholdCluster() const
{
  return fEThresholdCluster;
}
 
double EmcRecParameter::LogPosOffset() const
{
  return fLogPosOffset;
}
 
double EmcRecParameter::TimeMin() const
{
  return fTimeMin;
}
 
double EmcRecParameter::TimeMax() const
{
  return fTimeMax;
}
 
double EmcRecParameter::MoliereRadius() const
{
  return fMoliereRadius;
}
 
double EmcRecParameter::LateralProfile() const
{
  return fLateralProfile;
}
 
double EmcRecParameter::ECorr(int n) const
{
  return eCorr[n];
}
 
double EmcRecParameter::SigE(int n) const
{
  return sigE[n];
}
 
double EmcRecParameter::SigTheta(int n) const
{
  return sigTheta[n];
}
 
double EmcRecParameter::SigPhi(int n) const
{
  return sigPhi[n];
}
 
double EmcRecParameter::HitNb(int n) const
{
  return hitNb[n];
}
 
double EmcRecParameter::ElecBias(int n) const
{
  return elecBias[n];
}
 
double EmcRecParameter::SmCut(int n) const
{
  return smCut[n];
}
 
double EmcRecParameter::Peak(int n) const
{
  return peak[n];
}
 
void EmcRecParameter::SetPositionMode(std::vector<std::string>& mode)
{
  if(mode.size()==2) {
    positionMode1=mode[0];
    positionMode2=mode[1];
  }
}
 
//The following function is copied from PhotonCor/McCor
double EmcRecParameter::ECorrMC(double eg, double theid) const
{
  double Energy5x5=eg;
  if(eg<0.02)eg=0.02;
  if(eg>1.74)eg=1.74;
  if(theid<=0)theid=0.001;
  if(theid>=27)theid=26.999;
  Float_t einter = eg + 0.00001;
  Float_t tinter = theid+0.0001;
  double ecor=dt->Interpolate(einter,tinter);
  if(!(ecor))return Energy5x5;
  if(ecor<0.5)return Energy5x5;
  double EnergyCor=Energy5x5/ecor;
  return EnergyCor;
}
 
// Get energy error
double EmcRecParameter::ErrMC(double eg, double theid) const
{
  if(eg<0.02)eg=0.02;
  if(eg>1.74)eg=1.74;
  if(theid<=0)theid=0.001;
  if(theid>=27)theid=26.999;
  Float_t einter = eg + 0.00001;
  Float_t tinter = theid+0.0001;
  double err=dtErr->Interpolate(einter,tinter);
  return err;
}