VMinMethodFunctor.cc

Go to the documentation of this file.

#include "VMinMethodFunctor.h"

#include <TMath.h>
#include <iostream>
#include <cassert>

using std::cout;
using std::endl;

namespace MdLDFFinderAG {

VMinMethodFunctor::~VMinMethodFunctor()
{
};

ROOT::Minuit2::FunctionMinimum
  VMinMethodFunctor::Minimize(std::vector<double>& pars, std::vector<double>& epars, const std::vector<bool>& fpars)
const
{
  using namespace ROOT::Minuit2;

  DEBUGLOG("Starting");

  std::ostringstream info("");

  info.str("");
  info << "Candidate counters: " << fCandidateCounters << endl;
  info << "Silent counters: " << fSilentCounters;
  INFO(info.str());

  // Check the LDF
  assert ((*fLDFunction).GetName() == "KascadeGrandeLDF");

  DEBUGLOG("Start Minuit2 minimization process (LDF \"a la\" KascadeGrande)");
  MnUserParameters mnpars;
  mnpars.Add("Nmu"  , pars.at(eShowerSize), epars.at(eShowerSize) );
  mnpars.Add("Beta" , pars.at(eBeta), epars.at(eBeta), 0.0, 4.0 );
  mnpars.Add("xcore" , pars.at(eCoreX), epars.at(eCoreX) );
  mnpars.Add("ycore" , pars.at(eCoreY), epars.at(eCoreY) );

/* NOTE from http://root.cern.ch/root/html/TMinuit.html

The result of all the above is that for parameters with limits,
the error reported by Minuit is not exactly equal to the square
root of the diagonal element of the error matrix. The difference
is a measure of how much the limits deform the problem.
If possible, it is suggested not to use limits on parameters,
and the problem goes away.
If for some reason limits are necessary, and you are sensitive to
the difference between the two ways of calculating the errors,
it is suggested to use Minos errors which take into account the non-linearities much more precisely.

*/
  //Limits and fixed parameters
  mnpars.SetLowerLimit("Nmu", 0); // Number of muons must be positive
//  mnpars.SetLowerLimit("Beta",0.0);
//  mnpars.SetUpperLimit("Beta",4.0);  // Disallow steep LDFs

  //Fix criteria
  if (fpars.at(eBeta)) {
    mnpars.Fix("Beta");
  }
  if (fpars.at(eCoreX)) {
    mnpars.Fix("xcore");
  }
  if (fpars.at(eCoreY)) {
    mnpars.Fix("ycore");
  }
  DEBUGLOG("After fixing the fit parameters");

  //Construct from FCNBase + MnUserParameters
  MnMinimize minuit_mn( static_cast<const FCNBase&>((*this)), mnpars );
  DEBUGLOG("After creating the minimizer");

  // Run the minimization
  //FunctionMinimum is a class holding the full result of the minimization;
  //both internal and external (MnUserParameterState) representation available
  //for the parameters at the Minimum
  FunctionMinimum fnc_min = minuit_mn(10000);
  DEBUGLOG("After running the minimization");

  if(!fnc_min.IsValid()) {
    return fnc_min;
  }

/*
  if (fnc_min.IsAboveMaxEdm() || fnc_min.HasReachedCallLimit() || std::isnan(fnc_min.UserParameters().Params()[0])) {
    std::ostringstream info("");
    info << "\n**************** Minimize FAILED ****************\n"
      << fnc_min
      << "**************** Minimize FAILED ****************\n\n";
    INFO(info.str());
    return false;
  }
*/


  info.str("");
  info << "Found minimum after " << fnc_min.NFcn() << " evaluations";
  INFO(info.str());

  // Calculate the covariance matrix
  MnHesse hesseMtrx(1);
  hesseMtrx(minuit_mn.Fcnbase(), fnc_min);

  info.str("");
  info << fnc_min;
  INFO(info.str());

  //Get fit results
  pars  = fnc_min.UserParameters().Params();
  epars = fnc_min.UserParameters().Errors();

  DEBUGLOG("Finishing");

  return fnc_min;

}

void
  VMinMethodFunctor::SetAxis(utl::Vector axis)
{
  rAxis=axis;
  const utl::CoordinateSystemPtr siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
  cosTheta = axis.GetCosTheta(siteCS);
}

double
  VMinMethodFunctor::CalculateSilentLikelihood(const std::vector<double>& par)
const
{

  double nLogLikelihood = 0;

        const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();

  for (CounterList::const_iterator ic = fSilentCounters.begin();
       ic != fSilentCounters.end(); ++ic) {

      mevt::Counter* counter = *ic;
                        unsigned int counterId = counter->GetId();

      // calculate the number of muons expected in the SD station

      // core position
      const double xcore = par[2];
      const double ycore = par[3];
      const utl::Point core = utl::Point(xcore, ycore, zcore, siteCS);

                        utl::Point counterPosition = mDetector.GetCounter(counterId).GetPosition();

      const double coreDistance = DistanceInShowerPlane(counterPosition-core);

      // muon density in the shower plane
      const double muonDensity = (*fLDFunction)( coreDistance, &par[0] );

      // area of the SD station in m2
      const double sdArea = 10;

      // average muon signal in the SD, indedpendent of zenith angle!
      const double muonsInSd = muonDensity * sdArea;

      // calculate the likelihood of 1 counter
      double x = 1;  // argument of ln
      for(unsigned int i = 1; i <= fSilentLimit; ++i) {
        x += TMath::Power(muonsInSd,int(i)) / TMath::Factorial(i);
      }
      const double nLogLikelihood1 =  muonsInSd - TMath::Log(x);

      nLogLikelihood += nLogLikelihood1;
  }

  return nLogLikelihood;
}

// CounterList methods

std::ostream& operator<<(std::ostream& os, const CounterList& cl)
{

    for (CounterList::const_iterator ic = cl.begin(); ic != cl.end(); ++ic) {
      os << (*ic)->GetId() << " ";
    }

    return os;
}





} // end namespace