Likelihood2.cc

Go to the documentation of this file.

#include "Likelihood2.h"

#include <utl/ErrorLogger.h>
#include <utl/Math.h>
#include <utl/Trace.h>
//
#include <vector>
#include <iostream>
#include <limits>


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

namespace MdLDFFinderAG 
{

Likelihood2::Likelihood2(const VLDFFunctor* ldfFnc, bool useSil, unsigned int silLim) 
{
  fLDFunction = ldfFnc;
  fUseSilent = useSil;
  fSilentLimit = silLim;
  siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();
}

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

  nLogLikelihood += CalculateCandidateLikelihood(par);

  if (fUseSilent==true) {
    nLogLikelihood += CalculateSilentLikelihood(par);
  }

  return nLogLikelihood;

}

double
  Likelihood2::CalculateCandidateLikelihood(const std::vector<double>& par)
const
{
//  INFO("Starting");

  std::ostringstream info("");

//  cout << "Starting Likelihood2 evaluation" << endl ;

  double nLogLikelihood = 0;

  // Debug
/*  INFO("Parameters");
  cout << "Nmu = " << par[0] << endl;
  cout << "Beta = " << par[1] << endl;
  cout << "xcore = " << par[2] << endl;
  cout << "ycore = " << par[3] << endl;    
*/

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

//  info.str("");
//  info << "Core = (" << core.GetX(siteCS) << ", " << core.GetY(siteCS) << ", " << core.GetZ(siteCS) << ")" << endl;
//  INFO(info.str());

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

  for (CounterList::const_iterator ic = fCandidateCounters.begin(); 
     ic != fCandidateCounters.end(); ++ic) {
        
    const mevt::Counter* co = *ic;  

    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;  
      } 

      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);
      size_t activeChannels = im->GetNumberOfActiveChannels(); 
      const utl::TraceUI& channelsOn = im->GetRecData().GetChannelsOn();

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

      // DEBUG
/*      INFO("");
      cout << "Processing counter " << co->GetId() << " - module " << im->GetId() << endl
           << "position = (" << modulePosition.GetX(siteCS) << ", " << modulePosition.GetY(siteCS)
           << ", " << modulePosition.GetZ(siteCS) << ")" << endl 
           << "Distance to the core = " << coreDistance << endl;
      cout << "LDF(r) = " << muonDensity << endl;
      cout << "Predicted muons = " << predictedMuons << endl;
      cout << "channelsOn = ";       
      for(std::vector<size_t>::iterator it=channelsOn.begin(); it!=channelsOn.end(); ++it) {
         cout << *it << " ";
      }
      cout <<  endl;
*/

      ProfLike nLogLike1(activeChannels,channelsOn);
//      INFO("After creating the likelihood");

      const double nLogLikelihood1 = nLogLike1(predictedMuons); // calculate the negative of the log likelihood
//      INFO("After evaluating the likelihood");

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

      // DEBUG
/*      INFO(""); 
      cout << "nLogLikelihood1 = " << nLogLikelihood1;
      cout << ", nLogLikelihood = " <<  nLogLikelihood << endl;       
*/
    } 

  }  // end counter loop

//    cout << "Total -logL = " << nLogLikelihood << endl;
//    cout << "Finishing Likelihood2 evaluation" << endl ;

//  INFO("Finishing");

  return nLogLikelihood;
}
  
}