ProfileChi2.cc

Go to the documentation of this file.

#include "ProfileChi2.h"

#include <fevt/EyeRecData.h>

#include <evt/GaisserHillas4Parameter.h>

#include "../FdEnergyDepositFinderKG/ProfileFitFunctions.h"
#include "../FdEnergyDepositFinderKG/GHShapeParameters.h"


#include <utl/ErrorLogger.h>
#include <limits>

#include <fwk/CentralConfig.h>
#include <det/Detector.h>

using namespace fevt;
using namespace utl;
using namespace std;

namespace FdProfileConstrainedGeometryFit {

void ProfileChi2::Init() {
  Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch("FdProfileConstrainedGeometryFit");
  topBranch.GetChild("leavingAtmoIsError").GetData(fLeavingAtmoIsError);
  topBranch.GetChild("antiAliasingFilterCorrection").GetData(fAafCorrection); //VN: fitter is static
  topBranch.GetChild("useNoiseBins").GetData(fUseNoiseBins);
  ostringstream info;
  info<<"\n\t leavingAtmoIsError is " << (fLeavingAtmoIsError?"on":"off");
  info<<"\n\t antiAliasingFilterCorrection is " << (fAafCorrection?"on":"off");
  info<<"\n\t useNoiseBins is " << (fUseNoiseBins?"on":"off");
  INFO(info);
  topBranch.GetChild("onlyDirectLight").GetData(fOnlyDirectLight);
  
  const Branch profileBranch = topBranch.GetChild("profile");
  //VN to manage static fitter
  //FdEnergyDepositFinderKG::GHShapeParameters ghPars(profileBranch);
  
  //fProfileFitter.SetShapePars(ghPars);
  ghPars = FdEnergyDepositFinderKG::GHShapeParameters(profileBranch);
  
  Branch kUnivBranch = profileBranch.GetChild("kUnivConstrained");

  if (kUnivBranch.GetChild("constrained").Get<bool>()) {
    fIsUnivConstrained = true;
    kUnivBranch.GetChild("function").GetData(fUnivFunction);
    kUnivBranch.GetChild("ksigma").GetData(fksigma);
  }
}

double
ProfileChi2::operator()(Eye &eyeCopy)
{
  fCFMatrixCalculator.SetMethod(FdEnergyDepositFinderKG::CFMatrixCalculator::eFast);
  fCFMatrixCalculator.SetOnlyDirect(fOnlyDirectLight);
  fCFMatrixCalculator.SetVerbosity(-1);

  try {
    fCFMatrixCalculator.BuildMatrix(eyeCopy, fLeavingAtmoIsError);
  } catch (OutOfBoundException&) {
    return numeric_limits<double>::infinity();
  }
  if (!fCFMatrixCalculator.HasMatrix()
    || fCFMatrixCalculator.GetSize() < 2)
    {
      ERROR( "No Matrix built.");
      return numeric_limits<double>::infinity();
    }

  fProfileCalculator.CalculateProfile(eyeCopy, fCFMatrixCalculator, true);

  if (!fUseLightFlux) {
    if (FitProfile(eyeCopy,FdEnergyDepositFinderKG::eGH4dEdXChi2)
        /*&&  FitProfile(eyeCopy,FdEnergyDepositFinderKG::eGH2dEdXChi2)*/)
      {
        const double chi2pfl =
          eyeCopy.GetRecData().GetFRecShower().GetGHParameters().GetChiSquare();
        // /eyeCopy.GetRecData().GetFRecShower().GetGHParameters().GetNdof();
        return chi2pfl;
      }
    else return numeric_limits<double>::infinity();
  } else {
    if (/*FitProfile(eyeCopy,FdEnergyDepositFinderKG::eGH2dEdXChi2)
      &&*/ FitProfile(eyeCopy,FdEnergyDepositFinderKG::eGH4dEdXChi2) //VN GH2 not precise enough -> Xmax bias
      && FitProfile(eyeCopy,FdEnergyDepositFinderKG::eGH4LightLogLike))
      {
        const double chi2pfl = fProfileFitter.GetLHChi2();
        //VN from Chi2 to Likelihood
          //eyeCopy.GetRecData().GetFRecShower().GetGHParameters().GetChiSquare();
        // /eyeCopy.GetRecData().GetFRecShower().GetGHParameters().GetNdof();
        return chi2pfl;
      }
    else return numeric_limits<double>::infinity();
  }
}


bool
ProfileChi2::FitProfile(Eye& eye,
                        const FdEnergyDepositFinderKG::EFitFunctionType type)
{
  using namespace FdEnergyDepositFinderKG;
  using namespace evt;

  ShowerFRecData& shower = eye.GetRecData().GetFRecShower();
  //shower.SetEnergyCutoff(fElectronEnergyThreshold);
  fProfileFitter.ResetStartParameters();

  if (!shower.HasGHParameters()) {
    shower.MakeGHParameters(GaisserHillas4Parameter(ghPars.Type()));
    if (!GuessGHParameters(shower))
      return false;
    //DumpCurrentParameters(shower, false, eUndefined);
  }
  GaisserHillas4Parameter& gh4Pars = shower.GetGHParameters<GaisserHillas4Parameter>();

cout<<"Type: "<< GetFunctionTypeName(ghPars.Type())<<endl;
cout<<"L: "<< gh4Pars.GetShapeParameter(gh::eUSP_L)/(g/cm2)<< endl;
cout<<"R: "<< gh4Pars.GetShapeParameter(gh::eUSP_R)<< endl;


  //if (!shower.HasGHParameters<GaisserHillas4Parameter>()) {
  // return false; // we really need GH4 parameters
  //}

  
  fProfileFitter.SetFitMode(ProfileFitter::eNMax);
  {
    const GaisserHillas4Parameter& gh4Pars = shower.GetGHParameters<GaisserHillas4Parameter>();
    fProfileFitter.SetStartParameters(gh4Pars);
  }
  fProfileFitter.SetLightFluxData(fCFMatrixCalculator);

  fProfileFitter.SetProfileData(shower.GetEnergyDeposit());
  fProfileFitter.SetFunctionType(type);

  //VN: in separate inicialization function because fitter is static
  //fProfileFitter.SetUnivConstrained(fIsUnivConstrained, fA1, fA2, fksigma);

  if(fProfileFitter.Fit()){
    GaisserHillas4Parameter& gh4Pars = shower.GetGHParameters<GaisserHillas4Parameter>();
    fProfileFitter.FillGHParameters(gh4Pars);

    //    return true;
    double calorEnergy = gh4Pars.GetIntegral();
    double calorEnergyErr = calorEnergy * gh4Pars.GetIntegralError();

    // refit at final step using energy as parameter instead of Xmax
    if (type == eFinalFitStep) {
      fProfileFitter.SetFitMode(ProfileFitter::eIntegral);
      fProfileFitter.SetStartParameters(gh4Pars);
      if (fProfileFitter.Fit())
        fProfileFitter.GetEnergy(calorEnergy, calorEnergyErr);
      else
        return false;
    }
    shower.SetEmEnergy(calorEnergy, calorEnergyErr, ShowerFRecData::eStatistical);
    shower.SetEmEnergyError(0, ShowerFRecData::eAtmospheric);

    // rebuilt matrix if Xmax changed too much and not in final round
    const double currXmax = gh4Pars.GetXMax();
    const double maxDeltaXmax = 10*g/cm2;
    const bool updateCFM =
      fabs(currXmax-fCFMatrixCalculator.GetXmax()) > maxDeltaXmax &&
      type != eFinalFitStep;

    if (updateCFM)
      fCFMatrixCalculator.BuildMatrix(eye, fLeavingAtmoIsError);

    if (fCFMatrixCalculator.HasMatrix()) {

      fProfileCalculator.CalculateProfile(eye, fCFMatrixCalculator);
      //DumpCurrentParameters(shower, updateCFM, type);

      return true;
    }
  }

  return false;
}

bool
ProfileChi2::GuessGHParameters(evt::ShowerFRecData& shower)
  const
{
  using namespace FdEnergyDepositFinderKG;
  using namespace evt;
  
  GaisserHillas4Parameter& gh4Pars = shower.GetGHParameters<GaisserHillas4Parameter>();

  // running average

  const TabulatedFunctionErrors& energyDeposit = shower.GetEnergyDeposit();
  const unsigned int nPoints = energyDeposit.GetNPoints();
  const double depthBin = 50*g/cm2;

  bool firstTime = true;
  double xMax = 0;
  double dEdXmax = 0;
  for (unsigned int i = 0; i < nPoints; ++i) {
    const double firstDepth = energyDeposit.GetX(i);
    unsigned int nAverage = 0;
    double wSum = 0;
    double ywSum = 0;
    double xSum = 0;
    for (unsigned int j = i; j < nPoints; ++j) {
      const double thisDepth = energyDeposit.GetX(j);
      ++nAverage;
      xSum += thisDepth;
      const double w = 1 / pow(energyDeposit.GetYErr(j), 2);
      wSum += w;
      ywSum += energyDeposit.GetY(j)*w;

      if (thisDepth - firstDepth > depthBin) {
        const double yAverage = ywSum/wSum;
        if (firstTime || yAverage > dEdXmax) {
          dEdXmax = yAverage;
          xMax = xSum / nAverage;
          firstTime = false;
        }
        break;
      }
    }
  }

  if (ghPars.Type() == gh::eClassic) {
    //VN: values better for low energies? 
    gh4Pars.SetShapeParameter(gh::eX0, /*-120*/1*g/cm2, 0);
    gh4Pars.SetShapeParameter(gh::eLambda, /*61*/86*g/cm2, 0);
  }

 if (ghPars.Type() == gh::eUSP) {
       
   gh4Pars.SetShapeParameter(ghPars.Type(eOne),
                             ghPars.Mean(eOne,1e18*eV), ghPars.Sigma(eOne,1e18*eV));
   gh4Pars.SetShapeParameter(ghPars.Type(eTwo),
                             ghPars.Mean(eTwo,1e18*eV), ghPars.Sigma(eTwo,1e18*eV));   
  }

  if (firstTime) {
    WARNING("can not estimate start parameters");
    gh4Pars.SetXMax(700*g/cm2, 0);
    gh4Pars.SetNMax(1*PeV/g*cm2, 0, true);
    //gh4Pars.SetNMax(energyDeposit.GetY(0), 0, true);
    return false;
  } else {
    gh4Pars.SetXMax(xMax, 0);
    gh4Pars.SetNMax(dEdXmax, 0, true);
  }
  return true;
}

} //namespace