Likelihood3.cc

Go to the documentation of this file.

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

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


namespace MdLDFFinderAG {

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


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

    nLogLikelihood += CalculateCandidateLikelihood(par);

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

    return nLogLikelihood;
  }


  double
  Likelihood3::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;

        const unsigned int chOn = im->GetRecData().GetTotalNumberOfChannelsOn();

        const unsigned int segmentation = im->GetRecData().GetSegmentation();

        // check for saturation // TODO  use this to check if the user wants to use saturated modules or not
        //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

        //const double measuredMuons = CornerClippLike::GetNumberOfEstimatedMuonsCorrected(*im, mDetModule, rAxis, siteCS)

        const double ratio = double(chOn) / double(segmentation);

        const double pcc = im->GetRecData().GetCornerClippingProbability();

        nLogLikelihood1 = predictedMuons*(1+pcc)*(1. - ratio) - chOn * std::log(1 - std::exp(-predictedMuons*(1+pcc) / segmentation));
        //nLogLikelihood1 = predictedMuons - chOn * std::log(std::exp(predictedMuons / segmentation) -1);
        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;
  }

}