LikeFCN.cc

Go to the documentation of this file.

#include "LikeFCN.h"

#include <cmath>
#include <iostream>
#include <iomanip>

namespace MdLDFFinderAG {

  const double LikeFCN::EPSILON = 1e-4;


  LikeFCN::LikeFCN(const unsigned int n0, const SegmentsOnMulti& km0, const double mu0) :
    fNumberSegments(n0),
    fSegmentsOnMulti(km0),
    fMuons(mu0)
  { }


  LikeFCN::LikeFCN(const unsigned int n0, const std::vector<unsigned int>& k0, const double mu0) :
    fNumberSegments(n0),
    fMuons(mu0)
  {
    for (std::vector<unsigned int>::const_iterator it1 = k0.begin(); it1 != k0.end(); ++it1)
      fSegmentsOnMulti.push_back(std::make_pair(*it1, 1));
  }


  // Evaluates the negative of the logarithm of the likelihood (which is always positive!)
  // of mu_i expected muons in each time bin with the condition that the sum of the mu_i's
  // is a fixed mu (profile likelihood)
  // Input: par vector are mu_i with i = 1,...,w-1 and w the number of time bins
  // Output: negative of the log likelihood for mu

  double
  LikeFCN::operator()(const std::vector<double>& par)
    const
  {
    //  Build a vector holding the mu_i's (i=1,..,w), the last element is
    // (mu minus the sum of the mu_i's)

    std::vector<double> vmu(par);
    double parWsum = 0;  // sum vmu weighted with k multiplicy

    std::vector<double>::const_iterator it1;
    SegmentsOnMulti::const_iterator it2;

    for (it1 = par.begin(), it2 = fSegmentsOnMulti.begin();
         it1 != par.end(); ++it1, ++it2) {

      //std::cout << "mu_i = " << *it1 << ", m = " << (*it2).second << std::endl;

      parWsum += *it1 * (it2->second);

    }

    // calculate mulast using the multiplicity of the last window
    double mulast = (fMuons - parWsum) / fSegmentsOnMulti.back().second;

    if (mulast < EPSILON)
      mulast = EPSILON;  // trick to avoid divergences

    // Debug
    //std::cout << "In LikeFCN::operator(). mulast =  " << mulast << std::endl;

    vmu.push_back(mulast);

    // Evaluate likelihood
    double nlnLike = 0;      // negative log of the likelihood (>0!)

    for (it1 = vmu.begin(), it2 = fSegmentsOnMulti.begin();
         it2 != fSegmentsOnMulti.end() && it1 != vmu.end(); ++it1, ++it2) {

      const double k = it2->first;
      const double m = it2->second;
      const double p = 1 - exp(-(*it1) / fNumberSegments);

      nlnLike -= m*(k*log(p) + (fNumberSegments - k)*log(1 - p));

    }

    //std::cout << "In LikeFCN::operator(). nlnLike =  ";
    //std::cout << std::setprecision(9) << nlnLike << std::endl;

    return nlnLike;
  }


  double
  LikeFCN::operator()(const double* const mup, const double* const)
    const
  {
    std::vector<double> k;
    k.push_back(mup[0]);
    k.push_back(mup[1]);
    return this->operator()(k);
  }

}