FdDoubleBumpFinderKG/ProfileFitter.cc

Go to the documentation of this file.

#include "ProfileFitter.h"
using namespace profileFit;

#include <TMinuit.h>

#include <stdexcept>
#include <iostream>
#include <sstream>
#include <algorithm>
using namespace std;

unsigned int ProfileFitter::fNdof=0;
double ProfileFitter::fVariance=-1;
double ProfileFitter::fThreshold=0;
vector<double> ProfileFitter::fDepth;
vector<double> ProfileFitter::fSize;
vector<double> ProfileFitter::fSizeError;
EFitFunctionType ProfileFitter::fFunctionType=eGaisserHillas;
bool ProfileFitter::fIsConstrained=false;

ProfileFitter::ProfileFitter():
  fVerbosity(-1)
{

}

ProfileFitter::~ProfileFitter()
{

}
  
bool 
ProfileFitter::Fit() 
{

  if (fSizeError.empty() && fVariance < 0)
    {
      ostringstream errorMsg;

      if (fVerbosity > -1)
        {
          errorMsg << " ProfileFitter::Fit() " 
                   << " ERROR: variance is not initialized, exit!"; 
          cerr << errorMsg.str() << endl;
        }
      throw std::runtime_error(errorMsg.str());
    }
    
  fIsConstrained=true;

  switch ( fFunctionType ) {
  case eGaisserHillas:
    return GHFit();
  case eGaisserHillasUnconstrained:
    fIsConstrained=false;
    return GHFit();
  case eDoubleGaisserHillas:
    return DoubleGHFit();
  case eGaussInAge:
    return GaussInAgeFit();
  case eDoubleGaussInAge:
    return GaussInAgeFit();
  }

  fStartParameters.clear();
  fStepParameters.clear();
  
  return false;

}

bool 
ProfileFitter::GHFit() 
{
  
  const unsigned int nPar = eNGHpars;
  TMinuit theMinuit(nPar); 
  theMinuit.SetFCN(ProfileFitter::GHFitFunction);
  InitMinuit(theMinuit);

  if ( fStartParameters.size() != nPar ) {
    pair<double, double> theMax = FindMaximum();
    const double xMax = theMax.first;
    const double nMax = theMax.second;
    
    int ierflag;
    theMinuit.mnparm(eGHnMax,"nMax", nMax, nMax*5e-2,0,0,ierflag);
    theMinuit.mnparm(eGHXmax,"xMax", xMax,       10.,0,0,ierflag);
    theMinuit.mnparm(eGHX0,  "x0",-100.,       10.,-1000,500.,ierflag);
    theMinuit.mnparm(eGHLambda,"lamb",  60.,        5.,10.,150.,ierflag);
  }
  else
    UserStartInit(theMinuit);

  return Minimize(theMinuit);

}

bool
ProfileFitter::GaussInAgeFit()
{

  const unsigned int nPar=eNGiApars;
  TMinuit theMinuit(nPar); 
  theMinuit.SetFCN(ProfileFitter::GaussInAgeFitFunction);
  InitMinuit(theMinuit);

  if ( fStartParameters.size() != nPar ) {
    pair<double, double> theMax = FindMaximum();
    const double xMax = theMax.first;
    const double nMax = theMax.second;
    
    int ierflag;
    theMinuit.mnparm(eGiAnMax,"nMax", nMax,nMax*5e-2,0,0,ierflag);
    theMinuit.mnparm(eGiAXmax,"xMax", xMax,      10.,0,0,ierflag);
    theMinuit.mnparm(eGiASigma1,"sig1",  0.2,       .05,0.,.5,ierflag);
    theMinuit.mnparm(eGiASigma2,"sig2",  0.2,       .05,0.,.5,ierflag);
    
    if ( fFunctionType == eGaussInAge )
      theMinuit.FixParameter(eGiASigma2);

  }
  else 
    UserStartInit(theMinuit);
      
  return Minimize(theMinuit);

}

bool 
ProfileFitter::DoubleGHFit() 
{

  const unsigned int nPar=eNDGHpars;
  TMinuit theMinuit(nPar); 
  theMinuit.SetFCN(ProfileFitter::DoubleGHFitFunction);
  InitMinuit(theMinuit);


  if ( fStartParameters.size() != nPar ) {

    pair<double, double> theMax = FindMaximum();
    const double nMax = theMax.second;
    const double xMax = theMax.first;
    
    double leftSum=0.;
    double rightSum=0.;
    for ( unsigned int i=0;i<fDepth.size(); i++ ) {
      if ( fDepth[i] < xMax )
        leftSum+=fSize[i];
      else
        rightSum+=fSize[i];
    }
    
    double deltaX = 50.;
    if ( leftSum > rightSum )
      deltaX*= (-1);

    int ierflag;
    theMinuit.mnparm(eDGHnMax1,   "nMax1", nMax*0.5,     nMax*1e-1, 0.,     0.,    ierflag);
    theMinuit.mnparm(eDGHXmax1,   "xMax1", xMax-deltaX,  10.,       0.,     0.,    ierflag);
    theMinuit.mnparm(eDGHX01,     "x01",   -100.,        10.,       -1000., 500.,  ierflag);
    theMinuit.mnparm(eDGHLambda1, "lamb1", 61.,          5.,        10.,    150.,  ierflag);
    
    theMinuit.mnparm(eDGHnMax2,   "nMax2", nMax*0.5,     nMax*1e-1, 0.,     0.,    ierflag);
    theMinuit.mnparm(eDGHXmax2,   "xMax2", xMax+deltaX,  10.,       0.,     0.,    ierflag);
    theMinuit.mnparm(eDGHX02,     "x02",   -100.+deltaX, 10.,       -1000., 500.,  ierflag);
    theMinuit.mnparm(eDGHLambda2, "lamb2", 61.,          5.,        10.,    150.,  ierflag);
  }
  else
    UserStartInit(theMinuit);
  
  return Minimize(theMinuit);

}


std::pair<double, double> 
ProfileFitter::FindMaximum() 
  const
{

  vector<double>::const_iterator maximum = max_element(fSize.begin(),
                                                       fSize.end());

  const double nMax = *maximum;
  const int iMax    = (maximum-fSize.begin());
  const double xMax = fDepth[iMax];

  return pair<double,double>(xMax,nMax);

}


void 
ProfileFitter::GHFitFunction(int& /*npar*/, double* const /*gin*/, 
                             double& value, double* const par,
                             const int /*iFlag*/) 
{
  
  double chi2 = 0.;
  const double nMax = par[eGHnMax];
  const double xMax = par[eGHXmax];
  const double x0   = par[eGHX0];
  const double lambda = par[eGHLambda];

  fNdof=0;

  for ( unsigned int i=0;i<fDepth.size();i++ ) {
    
    const double y  = fSize[i];

    if ( fSizeError.size() || y > fThreshold ) {

      const double x  = fDepth[i];
      const double fx = (x>x0)?GHFunction(x,x0,lambda,xMax,nMax):0.;
      const double vy = fSizeError.empty()?fVariance:pow(fSizeError[i],2);

      chi2 += pow(y-fx,2)/vy;
      fNdof++;

    }
  }

  // constraints a la FdProfileReconstructor
  if ( fIsConstrained ) {
    const double fdefX0 = -121.;
    const double rmsX0  =  172.;
    chi2+=(x0-fdefX0)*(x0-fdefX0)/rmsX0/rmsX0;
    const double fdefLambda = 61.;
    const double rmsLambda  =  13.;
    chi2+=(lambda-fdefLambda)*(lambda-fdefLambda)/rmsLambda/rmsLambda;
    
    fNdof-=2;
  }

  value = chi2;
}

void 
ProfileFitter::GaussInAgeFitFunction(int& /*npar*/, double* const /*gin*/, 
                                     double& value, double* const par,
                                     const int /*iFlag*/) 
{
  
  double chi2 = 0.;
  const double nMax   = par[eGiAnMax];
  const double xMax   = par[eGiAXmax];
  const double sigma1 = par[eGiASigma1];
  const double sigma2 = (fFunctionType == eGaussInAge )?
    //    par[eGiASigma1]:
    // Giller ICRC 2005
    0.0340 + 0.7682*par[eGiASigma1]:par[eGiASigma2];

  fNdof=0;

  for ( unsigned int i=0;i<fDepth.size();i++ ) {
    
    const double y  = fSize[i];

    if ( y > fThreshold ) {

      const double x  = fDepth[i];
      const double fx = GaussInAgeFunction(x,nMax,xMax,sigma1,sigma2);
      const double vy = fSizeError.empty()?fVariance:pow(fSizeError[i],2);
      chi2 += pow(y-fx,2)/vy;
      fNdof++;

    }
  }

  if ( fFunctionType == eDoubleGaussInAge ) {
    // constraints a la FdProfileReconstructor
    const double meanSigma1  = 0.2;
    const double vSigma1 =  pow(0.03,2);
    const double meanSigma2  = 0.2;
    const double vSigma2 =  pow(0.03,2);
    chi2+=pow(sigma1-meanSigma1,2)/vSigma1;
    chi2+=pow(sigma2-meanSigma2,2)/vSigma2;
  }

  fNdof-=2;

  value = chi2;

}

void 
ProfileFitter::DoubleGHFitFunction(int& /*npar*/, double* const /*gin*/, 
                                   double& value, double* const par,
                                   const int /*iFlag*/) 
{
  double chi2 = 0.;
  const double nMax1   = par[eDGHnMax1];
  const double xMax1   = par[eDGHXmax1];
  const double x01     = par[eDGHX01];
  const double lambda1 = par[eDGHLambda1];
  const double nMax2   = par[eDGHnMax2];
  const double xMax2   = par[eDGHXmax2];
  const double x02     = par[eDGHX02];
  const double lambda2 = par[eDGHLambda2];

  fNdof=0;

  for (unsigned int i = 0; i < fDepth.size(); i++) {
    
    const double y = fSize[i];

    if ( y > fThreshold ) {

      const double x   = fDepth[i];
      const double fx1 = (x>x01)?GHFunction(x,x01,lambda1,xMax1,nMax1):0.;
      const double fx2 = (x>x02)?GHFunction(x,x02,lambda2,xMax2,nMax2):0.;
      const double vy  = fSizeError.empty()?fVariance:pow(fSizeError[i],2);

      chi2 += pow(y-fx1-fx2,2)/vy;
      fNdof++;

    }
  }

  // constraints a la FdProfileReconstructor but with constraint to X0-Xmax stead of X0 to comply with double GH fit

  const double fdefX0_minus_Xmax = -816.;
  const double rmsX0_minus_Xmax  =  183.;
  chi2+=(x01-xMax1-fdefX0_minus_Xmax)*(x01-xMax1-fdefX0_minus_Xmax)/rmsX0_minus_Xmax/rmsX0_minus_Xmax;
  chi2+=(x02-xMax2-fdefX0_minus_Xmax)*(x02-xMax2-fdefX0_minus_Xmax)/rmsX0_minus_Xmax/rmsX0_minus_Xmax;

  const double fdefLambda = 61.;
  const double rmsLambda  =  13.;
  chi2+=(lambda1-fdefLambda)*(lambda1-fdefLambda)/rmsLambda/rmsLambda;
  chi2+=(lambda2-fdefLambda)*(lambda2-fdefLambda)/rmsLambda/rmsLambda;

  fNdof-=4;

  //only positive dEdX values. This circumvents problem with limits in TMinuit
  if ( nMax1<0. || nMax2<0.)
    chi2+=1e99;

  value = chi2;
}



bool 
ProfileFitter::Minimize(TMinuit& theMinuit) 
{

  int ierflag;
  double arglist[2]={100000, 1.};
  theMinuit.mnexcm("MIGRAD", arglist, 2, ierflag);
  
  if (ierflag) {
    if (fVerbosity > -1)
      cerr << " MIGRAD failed. Error: " << ierflag << endl;
    return false;
  }

  theMinuit.mnexcm("HESSE", arglist, 2, ierflag);

  double amin, edm, errdef;
  int nvpar, nparx, icstat;
  theMinuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
  fStatus = icstat;
  fChi2 = amin;

  fFitParameters.resize(nparx);
  fFitErrors.resize(nparx);
  for ( int i=0;i<nparx;i++ ) 
    theMinuit.GetParameter(i, fFitParameters[i], fFitErrors[i]);

  return true;
}


void 
ProfileFitter::InitMinuit(TMinuit& theMinuit) 
  const 
{

  int ierflag;
  if(fVerbosity==2)
    theMinuit.SetPrintLevel(0); 
  else if (fVerbosity>=3)
    theMinuit.SetPrintLevel(fVerbosity); 
  else {
    theMinuit.SetPrintLevel(-1); 
    theMinuit.mnexcm("SET NOW", 0 ,0,ierflag);
  }

  double arglist=1.;
  theMinuit.mnexcm("SET ERR",&arglist,1,ierflag);

}

void 
ProfileFitter::UserStartInit(TMinuit& theMinuit) {

  int ierflag;
  for ( unsigned int i=0;i<fStartParameters.size();i++ ) {
    ostringstream name;
    name << "p" << i;
    if(fLimitsUp.size()!=fStartParameters.size() || fLimitsDown.size()!=fStartParameters.size())
      {
        theMinuit.mnparm(i,name.str().c_str(),
                         fStartParameters[i],
                         fStepParameters[i],0,0,ierflag);
        if(fVerbosity>0) cout << "no limits set" << fLimitsUp.size()<< " " << fStartParameters.size()<<endl;
      }
    else
      theMinuit.mnparm(i,name.str().c_str(),
                       fStartParameters[i],
                       fStepParameters[i],fLimitsDown[i],fLimitsUp[i],ierflag);
  }
}