Likelihood.cc

Go to the documentation of this file.

#include "Likelihood.h"
#include <sstream>

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


namespace MdLDFFinderAG {

  Likelihood::Likelihood(const VLDFFunctor* const ldfFnc, const bool useSil, const unsigned int silLim, const bool useSat, const size_t satLim)
  {
    fLDFunction = ldfFnc;
    fUseSilent = useSil;
    fSilentLimit = silLim;
    fUseSaturated = useSat;
    fSaturationLimit = satLim;
    siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
  }


  double
  Likelihood::operator()(const std::vector<double>& par)
    const
  {
    double nLogLikelihood = 0;

    nLogLikelihood += CalculateCandidateLikelihood(par);

    if (fUseSilent)
      nLogLikelihood += CalculateSilentLikelihood(par);

    return nLogLikelihood;
  }


  double
  Likelihood::CalculateCandidateLikelihood(const std::vector<double>& par)
    const
  {
    //DEBUGLOG( "Starting Likelihood evaluation");

    double nLogLikelihood = 0;
    std::ostringstream os;
    os.str("");
       
    os << "Nmu = " << par[0] << endl;
    //DEBUGLOG ( os.str() );

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

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

    for (auto ic = fCandidateCounters.begin(), icEnd = fCandidateCounters.end(); ic != icEnd; ++ic) {

      const mevt::Counter* co = *ic;

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

      // Loop over the modules
      for (auto im = co->ModulesBegin(); im != co->ModulesEnd(); ++im) {

        if (!im->IsCandidate()) // skip non candidate modules
          continue;

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

        const double moduleArea = im->GetRecData().GetActiveArea();

        const utl::Point modulePosition = mdetModule.GetPosition();

        const double coreDistance = DistanceInShowerPlane(modulePosition - core);

        // muon density in the shower plane
        const double muonDensity = (*fLDFunction)(coreDistance, &par[0]);
        // number of muons expected in the module
        const double predictedMuons = muonDensity * moduleArea * cosTheta;

        // check for saturation
        bool isSaturated = im->GetRecData().SaturationFlag();

        // DEBUG
        /*
        os.str("");
        os << "Processing counter " << co->GetId() << " - module " << im->GetId() << std::endl;
        os << "Segmentation " << im->GetRecData().GetSegmentation()  << std::endl;
        os << "position = (" << modulePosition.GetX(siteCS) << ", " << modulePosition.GetY(siteCS)
           << ", " << modulePosition.GetZ(siteCS) << ")" << std::endl
           << "Distance to the core = " << coreDistance/utl::m << std::endl;
        os << "Predicted muons = " << muonDensity << " * " <<  moduleArea << " * " << cosTheta
            << " = " << predictedMuons << std::endl;
        DEBUGLOG( os.str() );
        */

        // calculate the negative of the log likelihood
        double nLogLikelihood1; // likelihood of a counter
        if (!isSaturated) {  // if not saturared

          const double measuredMuons = im->GetRecData().GetNumberOfEstimatedMuons();
          nLogLikelihood1 = predictedMuons - measuredMuons * std::log(predictedMuons);

        } else {  // if saturated

          const unsigned int windowsOn = im->GetRecData().GetNumberOfChannelsOn();
          nLogLikelihood1 = -utl::LogarithmOfNormalComplementCDF(windowsOn, predictedMuons, std::sqrt(predictedMuons));
        }

        nLogLikelihood += nLogLikelihood1; // add the likelihood of this module

      }

    }
    /*
    os.str("");
    os << "Total -logL = " << nLogLikelihood << std::endl;
    os << "Finishing Likelihood evaluation" << std::endl ;
    DEBUGLOG( os.str() );
    */
    return nLogLikelihood;
  }

}