MdLDFFinder.cc

Go to the documentation of this file.

#include "MdLDFFinder.h"
#include "KascadeGrandeLDF.h"
#include "Likelihood.h"
#include "Likelihood2.h"
#include "Likelihood3.h"
//
#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerMRecData.h>
//
#include <mevt/MEvent.h>
#include <mevt/Counter.h>
#include <mevt/Module.h>
#include <mevt/ModuleRecData.h>
//
#include <fwk/CentralConfig.h>
#include <utl/ConsecutiveEnumFactory.h>
#include <utl/AugerUnits.h>
#include <utl/Branch.h>
#include <utl/ConfigUtils.h>
#include <fwk/RunController.h>
//
#include <fwk/LocalCoordinateSystem.h>
#include <utl/ErrorLogger.h>
#include <utl/CoordinateSystem.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/AxialVector.h>
#include <utl/Math.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Line.h>
//
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameterState.h>
#include <Minuit2/MnUserCovariance.h>
#include <boost/iterator/transform_iterator.hpp>
#include <string>
#include <iostream>
#include <vector>
#include <sstream>

using namespace fwk;
using namespace evt;
using std::cout;
using std::endl;


namespace MdLDFFinderAG {

const char* const MdLDFFinder::kLDFTags[] = { "KascadeGrande", "ToBeDef_I", "ToBeDef_II" };

const char* const MdLDFFinder::kMinMethodTags[] = { "Likelihood", "Likelihood2", "Likelihood3" };


MdLDFFinder::MdLDFFinder()
{ }


MdLDFFinder::~MdLDFFinder()
{ }


VModule::ResultFlag
MdLDFFinder::Init()
{
  // Configuration reading.
  utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdLDFFinder");
  std::string tag;

  LoadConfig(config, "ldfType", tag, kLDFTags[0]);//default LDF
  //Creates enumeration from tags
  fLDFType = utl::ConsecutiveEnumFactory<LDFType, eUndef_II, kLDFTags>::Create(tag);

  //LDF Parameter
  LoadConfig(config, "referenceDistance", fReferenceDistance, 450*utl::m );

  LoadConfig(config, "a0", fA0,  3.00 );
  LoadConfig(config, "a1", fA1, -1.00 );
  LoadConfig(config, "a2", fA2,  0.00 );

  //Fitting and minimization methods
  LoadConfig(config, "minMethodType", tag, kMinMethodTags[0]);//default MinMethod
  //Creates enumeration from tags
  fMinMethodType = utl::ConsecutiveEnumFactory<MinMethodType, eLikelihood3, kMinMethodTags>::Create(tag);

  LoadConfig(config, "useSilentCounters",fUseSilent, true);
  LoadConfig(config, "silentLimit", fSilentLimit, 2 );

  LoadConfig(config, "useSaturatedCounters",fUseSaturated, false);
  LoadConfig(config, "saturationLimit", fSaturationLimit, 64);

  //Minimum Number of good stations to fix beta in the  MLDF fit
  LoadConfig(config, "fixBeta", fFixBeta, false );
  LoadConfig(config, "countersToFixBeta", fCountersToFixBeta, 5);

  //Criteria to fix or relax beta
  //(half those corresponding to regular array)
  //Range where counters are considered
  LoadConfig(config, "LowLimToRelaxBeta", fLowLimToRelaxBeta, 500*utl::m/2);
  LoadConfig(config, "UppLimToRelaxBeta", fUppLimToRelaxBeta, 1500*utl::m/2);
  //Distances between counters (depending on the number of counters)
  LoadConfig(config, "CounterSpacings1", fCounterSpacings1, 500*utl::m/2 );
  LoadConfig(config, "CounterSpacings2", fCounterSpacings2, 400*utl::m/2 );
  LoadConfig(config, "CounterSpacings3", fCounterSpacings3, 300*utl::m/2 );

  LoadConfig(config, "fixCoreFromSd", fFixCoreFromSd, 0 );

  // Init according to selected LDF.
  switch (fLDFType) {
    case eKascadeGrande:
      fLDF.reset( new KascadeGrandeLDF(fReferenceDistance) );
      break;
    default:
      ERROR("LDF UNDEFINED");
      return eFailure;
      //fLDF.reset( new Undefined_I() );
      break;
  }

  // Init according to selected MinMethod.
  switch (fMinMethodType) {
    case eLikelihood:
      INFO("Original likelihood will be used");
      fMinimizer.reset( new Likelihood( fLDF.get(),  fUseSilent, fSilentLimit, fUseSaturated, fSaturationLimit) );
      break;
    case eLikelihood2:
      INFO("Profile likelihood will be used");
      fMinimizer.reset( new Likelihood2(fLDF.get(), fUseSilent, fSilentLimit) );
      break;
    case eLikelihood3:
      INFO("Infinite window likelihood will be used");
      fMinimizer.reset( new Likelihood3(fLDF.get(), fUseSilent, fSilentLimit,  fUseSaturated) );
      break;
  }

  mDetector = &(det::Detector::GetInstance().GetMDetector());
  siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();

  return eSuccess;
}


VModule::ResultFlag
MdLDFFinder::Run(evt::Event& event)
{
  std::ostringstream info("");

  INFO("Starting the muon LDF reconstruction");

  // nothing to do if there is no RecShower
  if (!(event.HasRecShower() &&
        event.GetRecShower().HasSRecShower()) ) {
    ERROR("Current event does not have a recontructed shower.");
    return eContinueLoop;
  }

  if (!event.HasMEvent()) {
    ERROR("Current event does not have an MD event.");
    return eContinueLoop;
  }

  // Pass core position and shower axis from SD reconstruction to minimizer
  const ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower() ;
  const utl::Point& sdCore = sRecShower.GetCorePosition();
  const utl::Vector& sdAxis = sRecShower.GetAxis();

  info.str("");
  info << "Core = (" << sdCore.GetX(siteCS) << ", " << sdCore.GetY(siteCS) << ", " << sdCore.GetZ(siteCS) << ")";
  DEBUGLOG(info.str());

  // Select counters with data
  mevt::MEvent& me = event.GetMEvent();
  const mdet::MDetector& mDetector = det::Detector::GetInstance().GetMDetector();

  // Filling modules SP distances
  FillModulesShowerPlaneDistances(sRecShower, mDetector, me);

  CounterList fCandidateCounters = SelectCandidateCounters(me);

  info.str("");
  info << "candidate counters = " << fCandidateCounters.size();
  DEBUGLOG(info.str());

  info.str("");
  info << "Candidate counters: " << fCandidateCounters;
  DEBUGLOG(info.str());

  // Select silent counters
  CounterList fSilentCounters = SelectSilentCounters(me);

  // Estimate initial values of rho0 and beta
  double rho0Seed = CalculateRho0Seed( fCandidateCounters, sdCore, sdAxis );
  double rho0SeedError = .5*rho0Seed; // using some gross estimate...

  // Define a coordinate system centered in the core reconstructed by the SD
  utl::CoordinateSystemPtr sdCoreCS = fwk::LocalCoordinateSystem::Create(sdCore);

  info.str("");
  info << "Core = (" << sdCore.GetX(sdCoreCS) << ", " << sdCore.GetY(sdCoreCS) << ", " << sdCore.GetZ(sdCoreCS) << ")" << std::endl;
  info << " theta = " << sdAxis.GetTheta(sdCoreCS)/utl::deg << " [deg]\n";
  DEBUGLOG(info.str());

  size_t numberOfFitParameters = 4; // rho0, beta, xcore, ycore

  std::vector<double> fPars(numberOfFitParameters);  // Estimators
  std::vector<double> fErrPars(numberOfFitParameters); // Errors
  std::vector<bool> fFlagPars(numberOfFitParameters); // Fix/free parameter flags

  fPars.at(VMinMethodFunctor::eShowerSize) = rho0Seed;
  fErrPars.at(VMinMethodFunctor::eShowerSize) = rho0SeedError;
  fFlagPars.at(VMinMethodFunctor::eShowerSize) = false;

  double betaSeed = CalculateBeta( sdAxis.GetTheta(sdCoreCS) );
  //Criteria to fix beta in this event
  bool fFixBetaEvent = FixBetaFlag( fCandidateCounters, sdCore, sdAxis );

  fPars.at(VMinMethodFunctor::eBeta)    = betaSeed;
  fErrPars.at(VMinMethodFunctor::eBeta) = 0;  //TO-DO: add error in axis (by means of cos(theta)) recontruction
  fFlagPars.at(VMinMethodFunctor::eBeta) = fFixBetaEvent;

  // Set the initical core as the SD core
  double xcore = sdCore.GetX(siteCS); // referred to the siteCS as the counters
  double ycore = sdCore.GetY(siteCS);
  double zcore = sdCore.GetZ(siteCS);

  fPars.at(VMinMethodFunctor::eCoreX)    = xcore;
  fErrPars.at(VMinMethodFunctor::eCoreX) = 100;
  fFlagPars.at(VMinMethodFunctor::eCoreX) = fFixCoreFromSd;

  ycore = sdCore.GetY(siteCS);
  fPars.at(VMinMethodFunctor::eCoreY)    = ycore;
  fErrPars.at(VMinMethodFunctor::eCoreY) = 100;
  fFlagPars.at(VMinMethodFunctor::eCoreY) = fFixCoreFromSd;

  // Configure the minimizer
  fMinimizer->SetZcore(zcore);
  fMinimizer->SetAxis(sdAxis);
  fMinimizer->SetCandidateCounterList(fCandidateCounters);
  fMinimizer->SetSilentCounterList(fSilentCounters);

  // Run the minimization
  ROOT::Minuit2::FunctionMinimum functionMinimum = fMinimizer->Minimize(fPars, fErrPars, fFlagPars);

  if (!functionMinimum.IsValid()) {
      ERROR("MD reconstruction failed");
      std::cerr << functionMinimum << endl;
      ++RunController::GetInstance().GetRunData().GetNamedCounters()["MdLDFFinder/Failed"];
      return eContinueLoop;
  }

  if (!functionMinimum.HasValidCovariance()) {
      ERROR("MD reconstruction: invalid covariance matrix");
      std::cerr << functionMinimum << endl;
      ++RunController::GetInstance().GetRunData().GetNamedCounters()["MdLDFFinder/Failed"];
      return eContinueLoop;
  }

  ++RunController::GetInstance().GetRunData().GetNamedCounters()["MdLDFFinder/Reconstructed"];

  // Fill the event with reconstructed data
  if (!(event.GetRecShower().HasMRecShower()) )
    event.GetRecShower().MakeMRecShower();

  ShowerMRecData& mRecShower = event.GetRecShower().GetMRecShower();

  if (!mRecShower.HasMLDF())
    mRecShower.MakeMLDF();

  FillReconstructedParameters(mRecShower, functionMinimum.UserState(), zcore);

  mRecShower.SetLdfReconstructed(true);
  mRecShower.SetReferenceDistance( fReferenceDistance );
  //mRecShower.SetAxis(sdAxis);
  mRecShower.SetReferenceAxis(sdAxis);
  mRecShower.SetReferenceCorePosition(sdCore);


  // Set the residuals and calculate the chi2
  SetLdfResiduals(fCandidateCounters, fPars, sdCore, sdAxis);
  double chi2 = CalculateChi2(fCandidateCounters, fPars, sdCore, sdAxis);
  size_t ndof = CalculateDegreesOfFreedom(fCandidateCounters, fFlagPars);
  mRecShower.SetMLDFChi2(chi2, ndof);

  // Set the likelihood (using the -log likelihood!)
  const double nlnLikelihood = fMinimizer->operator()(fPars);
  mRecShower.SetMLDFLikelihood(nlnLikelihood);

  FillTabulatedFunction( mRecShower, fPars, fErrPars );

  return eSuccess;
}


VModule::ResultFlag
  MdLDFFinder::Finish()
{
  return eSuccess;
}


// Private methods

bool
MdLDFFinder::FixBetaFlag( const CounterList counters, const utl::Point& rCore, const utl::Vector& rAxis )
  const
{
  std::ostringstream info;

  if (fFixBeta || counters.size() < fCountersToFixBeta) {  // if fix beta specified in the input
    info.str("");
    info << "Fixing Beta (" << counters.size() << " counters when at least "
         << fCountersToFixBeta << " are required)";
    INFO(info);
    return true;
  }

  typedef CounterList::const_iterator CIt;

  int nCountersInRange = 0;
  double maxDistance   = 0;

  info.str("");
  info << "Core = (" << rCore.GetX(siteCS) << ", " << rCore.GetY(siteCS) << ", " << rCore.GetZ(siteCS) << ')';
  DEBUGLOG(info);

  for (CIt ci = counters.begin(), cend = counters.end(); ci != cend; ++ci) {

    const unsigned int counterId_i = (*ci)->GetId();
    const utl::Point& position_i = mDetector->GetCounter(counterId_i).GetPosition();

    const double rho_i = utl::Cross(rAxis, position_i - rCore).GetMag();

    info.str("");
    info << "counter " << (*ci)->GetId() << ", position = ("
         << position_i.GetX(siteCS) / utl::m << ", "
         << position_i.GetY(siteCS) / utl::m << ", "
         << position_i.GetZ(siteCS) / utl::m << ") m"
         << ", distance to core = " << rho_i;
    DEBUGLOG(info);

    if (fLowLimToRelaxBeta < rho_i && rho_i < fUppLimToRelaxBeta) {
      ++nCountersInRange;
      for (CIt cj = ci+1; cj != cend; ++cj) {

        const unsigned int counterId_j = (*cj)->GetId();
        const utl::Point& position_j = mDetector->GetCounter(counterId_j).GetPosition();

        const double rho_j = utl::Cross(rAxis, position_j - rCore).GetMag();

        if (fLowLimToRelaxBeta < rho_j && rho_j < fUppLimToRelaxBeta) {
          const double distance = std::abs(rho_i - rho_j);
          if (distance > maxDistance)
            maxDistance = distance;
        }
      } // end loop cj
    } // end if

  } // end loop ci

  const bool rvalue = !(
    (nCountersInRange > 1 && maxDistance > fCounterSpacings1) ||
    (nCountersInRange > 2 && maxDistance > fCounterSpacings2) ||
    (nCountersInRange > 3 && maxDistance > fCounterSpacings3)
  );

  info.str("");
  info << "Using " << (rvalue?"fixed":"free") << " beta in the LDF fit "
       << "nCountersInRange " << nCountersInRange << " (out of " <<  counters.size() << ")\n"
       << "maxDistance " << maxDistance/utl::m << " m";
  INFO(info);

  return rvalue;
}


void
MdLDFFinder::FillTabulatedFunction(evt::ShowerMRecData& mRecShower,
                                   const std::vector<double>& fPars,
                                   const std::vector<double>& /*fErrPars*/)
{
  std::ostringstream info("");
/*
  info << "Pars -- Errors\n";
  std::vector<double>::const_iterator it1, it2;
  for ( it1=fPars.begin(), it2=fErrPars.begin(); it1!=fPars.end() && it2!=fErrPars.end(); ++it1, ++it2 ) {
    info << *it1 << " +/- " << *it2 << std::endl;
  }
  DEBUGLOG(info.str());
*/
  utl::TabulatedFunctionErrors& tabulatedMLDF = mRecShower.GetMLDF();
  tabulatedMLDF.Clear();

  const double minRho = 50*utl::m;
  const double maxRho = 3500*utl::m;

  const int nMLDFPoints = 250;
  const double dRho = (maxRho - minRho) / (nMLDFPoints - 1);

  for (int i = 0; i < nMLDFPoints; ++i) {
    const double rho = minRho + i*dRho;
    const double Nmu_r =(*fLDF)(rho, &fPars[0]);
    const double sigma = sqrt(Nmu_r);
    tabulatedMLDF.PushBack(rho, 0, Nmu_r, sigma);
  }
}


void
MdLDFFinder::FillReconstructedParameters(evt::ShowerMRecData& mRecShower,
                                         const ROOT::Minuit2::MnUserParameterState& MnParState,
                                         double zcore)
{
  //DEBUGLOG("");
  //std::cout << MnParState << std::endl;

  // Retrieve parameter and errors from Minuit
  const std::vector<ROOT::Minuit2::MinuitParameter> minuitParameters(MnParState.MinuitParameters());
  //std::vector<double> fPars(MnParState.Params());
  //std::vector<double> fErrPars(MnParState.Errors());

  // Muon density
  double muonDensity0 = minuitParameters[VMinMethodFunctor::eShowerSize].Value();
  double muonDensity0Error = minuitParameters[VMinMethodFunctor::eShowerSize].Error();
  mRecShower.SetNMuRef(muonDensity0, muonDensity0Error);

  // LDF slope
  const double beta = minuitParameters[VMinMethodFunctor::eBeta].Value();
  if (minuitParameters[VMinMethodFunctor::eBeta].IsFixed()) {
    mRecShower.SetBetaFixed(true);
    mRecShower.SetBeta(beta);
  } else {
    mRecShower.SetBetaFixed(false);
    const double betaError = minuitParameters[VMinMethodFunctor::eBeta].Error();
    mRecShower.SetBeta(beta,betaError);
  }

  // Core position
  const double xcore = minuitParameters[VMinMethodFunctor::eCoreX].Value();
  const double ycore = minuitParameters[VMinMethodFunctor::eCoreY].Value();
  const utl::Point core = utl::Point(xcore, ycore, zcore, siteCS);
 // mRecShower.SetCorePosition(core);

  // Filling the core position from the MLDF minimizer only if it was left as free paramerter (ie if a geometry reconstruction was done with UMD)
  // To acces the sd core position use fReferenceCore
  if (mRecShower.IsGeometryReconstructed()) {
   mRecShower.SetCorePosition(core);
  }

  if (minuitParameters[VMinMethodFunctor::eCoreX].IsFixed()) {

    mRecShower.SetCoreFixedLdf(true);

  } else {

    mRecShower.SetCoreFixedLdf(false);

    // Core errors
    const double xcoreError = minuitParameters[VMinMethodFunctor::eCoreX].Error();
    const double ycoreError = minuitParameters[VMinMethodFunctor::eCoreY].Error();
    const double zcoreError = 0;        // approximate core altitude by barycenter
    const utl::Vector coreError = utl::Vector(xcoreError, ycoreError, zcoreError, siteCS);
    mRecShower.SetCoreError(coreError);

    // Core xy correlation
    const ROOT::Minuit2::MnUserCovariance covarianceMatrix = MnParState.Covariance();
    unsigned int xcoreIndex = GetFreeParIndex(minuitParameters, VMinMethodFunctor::eCoreX);
    unsigned int ycoreIndex = GetFreeParIndex(minuitParameters, VMinMethodFunctor::eCoreY);
    const double covaxx = covarianceMatrix(xcoreIndex,xcoreIndex);
    const double covayy = covarianceMatrix(ycoreIndex,ycoreIndex);
    const double covaxy = covarianceMatrix(xcoreIndex,ycoreIndex);
    const double correlationXY =  covaxy / sqrt(covaxx*covayy);
    mRecShower.SetCorrelationXY(correlationXY);

  } // end if

}


unsigned int
MdLDFFinder::GetFreeParIndex(const std::vector<ROOT::Minuit2::MinuitParameter> minuitParameters,
                             unsigned int globalIndex)
  const
{
  // check boundaries
  assert(globalIndex<minuitParameters.size());

  int freeParIndex = -1;
  for (unsigned int i=0; i<=globalIndex; ++i) {
    if(!minuitParameters[i].IsFixed()) {
       freeParIndex++;
    }
  }

  assert(freeParIndex>=0);

  return freeParIndex;
}


double
MdLDFFinder::CalculateBeta(double theta)
  const
{
  return fA0 + fA1/cos(theta) + fA2/pow(cos(theta),2);
}


double
MdLDFFinder::CalculateRho0Seed(const CounterList counters, const utl::Point rCore,
                               const utl::Vector rAxis )
  const
{

  double minDist = std::numeric_limits<double>::max();

  mevt::Counter* nearestCounter = 0;

  // find the counter closest to the reference distance
  for (CounterList::const_iterator ic = counters.begin(); ic != counters.end(); ++ic) {

    mevt::Counter* co = *ic;

    unsigned int counterId = co->GetId();
    utl::Point counterPosition = mDetector->GetCounter(counterId).GetPosition();

    double distanceToCore = cross(rAxis, counterPosition - rCore).GetMag();

    double distanceToR0 = fabs(fReferenceDistance - distanceToCore);

    if (distanceToR0<minDist) {
      nearestCounter = co;
    }

  }

  double density = nearestCounter->GetMuonDensity();

  // minimal density to be used as a seed - avoids reconstruction issues with seed density=0
  const double minSeedDensity = 0.05;

  density = ( density>minSeedDensity? density : minSeedDensity);


  return density;
}


CounterList
MdLDFFinder::SelectSilentCounters(mevt::MEvent& me)
{
  CounterList fSilentCounters;

  // Populate the list of candidate modules and silent counters
  for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd();
     ic != ec; ++ic) {

    if (ic->IsSilent()) {
      fSilentCounters.push_back(&(*ic));
    }

  }  // end counter loop

  return fSilentCounters;
}


CounterList
MdLDFFinder::SelectCandidateCounters(mevt::MEvent& me)
{
  CounterList fCandidateCounters;

  // Populate the list of candidate modules and silent counters
  for (mevt::MEvent::CounterIterator ic = me.CountersBegin(), ec = me.CountersEnd();
     ic != ec; ++ic) {

    if (ic->IsCandidate()) {
      if (fMinMethodType == eLikelihood ){//This check is not necessary in ProfileLikelihood case
        // Loop over the modules to check user defined saturation limit (different from hardaware imposed limit by segmentation)
        for (mevt::Counter::ModuleIterator im = ic->ModulesBegin(); im != ic->ModulesEnd(); ++im) {
          if(!im->IsCandidate()) {// skip non candidate modules
            continue;
          }
          size_t maxWindowsOn = im->GetRecData().GetMaxChannelsOn();
          im->GetRecData().SetSaturationFlag( ((im->GetRecData().IsSaturated()||maxWindowsOn>=fSaturationLimit) ? true : false) );
        }//end if on modules
    }//end if likelihood
    fCandidateCounters.push_back(&(*ic));
    }

  }  // end counter loop

  return fCandidateCounters;
}


void
MdLDFFinder::SetLdfResiduals(CounterList counters, const std::vector<double> parameters,
                             const utl::Point rCore, const utl::Vector rAxis)
  const
{
  // Get the cosine of the zenith angle in the core coordinate system
  const utl::CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(rCore);
  const double cosTheta = rAxis.GetCosTheta(coreCS);

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

    mevt::Counter* co = *ic;

    unsigned int counterId = co->GetId();
    const mdet::Counter& mdetCounter = mDetector->GetCounter(counterId);

    // Loop over the modules
    for (mevt::Counter::ModuleIterator im = co->ModulesBegin(); im != co->ModulesEnd(); ++im) {
      if(im->IsCandidate() && !im->GetRecData().IsSaturated()) {

        unsigned int moduleId = im->GetId();
        const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

        const  utl::Point r_i = mdetModule.GetPosition();
        double distanceToCore = cross( rAxis, r_i - rCore ).GetMag();

        //density of muons from the ldf
        double ldfDensity = (*fLDF)( distanceToCore, &parameters[0] );

       // density of muons from the module data corrected by zenith angle
        double dataDensity = im->GetRecData().GetMuonDensity() / cosTheta;

        // simple chi2 like residual
        double residual = (dataDensity - ldfDensity);

        im->GetRecData().SetLDFResidual(residual);
      }
    } // end module loop

  } // end counter loop
}


double
MdLDFFinder::CalculateChi2(const CounterList counters, const std::vector<double> parameters,
                           const utl::Point rCore, const utl::Vector rAxis)
  const
{
  double chi2(0);

  // Get the cosine of the zenith angle in the core coordinate system
  const utl::CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(rCore);
  const double cosTheta = rAxis.GetCosTheta(coreCS);

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

    const mevt::Counter* co = *ic;

    unsigned int counterId = co->GetId();
    const mdet::Counter& mdetCounter = mDetector->GetCounter(counterId);

    // Loop over the modules
    for (mevt::Counter::ModuleConstIterator im = co->ModulesBegin(); im != co->ModulesEnd(); ++im) {
      if(im->IsCandidate() && !im->GetRecData().IsSaturated()) {

        unsigned int moduleId = im->GetId();
        const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

        const  utl::Point r_i = mdetModule.GetPosition();
        double distanceToCore = cross( rAxis, r_i - rCore ).GetMag();

        //number of muons from the ldf
        double ldfDensity = (*fLDF)( distanceToCore, &parameters[0] );
        double moduleArea = im->GetRecData().GetActiveArea();
        double ldfMuons = ldfDensity*moduleArea*cosTheta;

        //number of muons from the module data
        double measuredMuons = im->GetRecData().GetNumberOfEstimatedMuons();

        double dMuons = measuredMuons-ldfMuons;
        chi2 += dMuons*dMuons/ldfMuons;

//       Debug
//       DEBUGLOG("");
//       cout << "Module " << co->GetId() << "-" << im->GetId() << ", measured = " << measuredMuons;
//       cout << ", expected = " <<  ldfMuons << ", chi2 = " << dMuons*dMuons/ldfMuons << endl;

      }
    } // end module loop

  } // end counter loop

  return chi2;
}


size_t
MdLDFFinder::CalculateDegreesOfFreedom(const CounterList counters, const std::vector<bool> parameterFlags)
  const
{
  size_t ndof(0);

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

    const mevt::Counter* co = *ic;

    // Loop over the modules
    for (mevt::Counter::ModuleConstIterator im = co->ModulesBegin(); im != co->ModulesEnd(); ++im)
      if (im->IsCandidate() && !im->GetRecData().IsSaturated()) {
        ++ndof;
      } // end module loop
  } // end counter loop

  // count the number of free parameters
  const size_t nFitParameters = std::count(parameterFlags.begin(),parameterFlags.end(), false);

  // substract the number of fit parameters to the number of degree of freedom
  ndof -= nFitParameters;

  return ndof;
}

void
MdLDFFinder::FillModulesShowerPlaneDistances(const ShowerSRecData& sRecShower, const mdet::MDetector& mdet, mevt::MEvent& me)
{
  const utl::CoordinateSystemPtr sdShowerCS = sRecShower.GetShowerCoordinateSystem();

  for (mevt::MEvent::CounterIterator cIt = me.CountersBegin(); cIt != me.CountersEnd(); ++cIt) {
    int cId = cIt->GetId();
    const mdet::Counter& cDet = mdet.GetCounter(cId);
    for (mevt::Counter::ModuleIterator mIt = cIt->ModulesBegin(); mIt != cIt->ModulesEnd(); ++mIt) {
      if (!mIt->HasRecData())
        continue;
      int mId = mIt->GetId();
      const mdet::Module& mDet = cDet.GetModule(mId);
      const utl::Point& mPosition = mDet.GetPosition();
      const double r = mPosition.GetRho(sdShowerCS);
      mevt::ModuleRecData& mRecData = mIt->GetRecData();
      mRecData.SetSPDistance(r);

    }
  }

}


} // end namespace MdLDFFinderAG