AbsLhNew.cc

#include <getopt.h>
#include <fstream>
#include <string>

#include "PwaUtils/AbsLhNew.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "ErrLogger/ErrLogger.hh"

AbsLhNew::AbsLhNew(boost::shared_ptr<const EvtDataBaseListNew> theEvtList) :
  _evtListPtr(theEvtList)
{
  _evtDataVec=_evtListPtr->getDataVecs();
  _evtMCVec=_evtListPtr->getMcVecs();
}

AbsLhNew::AbsLhNew(boost::shared_ptr<AbsLhNew> theAbsLhPtr):
  _evtListPtr(theAbsLhPtr->getEventList())
{
  _evtDataVec=_evtListPtr->getDataVecs();
  _evtMCVec=_evtListPtr->getMcVecs();
}

AbsLhNew::~AbsLhNew()
{
}

double AbsLhNew::calcLogLh(fitParamsNew& theParamVal){
 
  double logLH=0.;
  double logLH_data=0.;
  double weightSum=0.;

  std::vector<EvtDataNew*>::iterator iterd;
  for (iterd=_evtDataVec.begin(); iterd!=_evtDataVec.end(); ++iterd){
    double intensity=calcEvtIntensity((*iterd), theParamVal);
    if (intensity>0.) logLH_data+=((*iterd)->evtWeight)*log(intensity);
    weightSum+= (*iterd)->evtWeight;
  } 

  double LH_mc=0.;
  
  std::vector<EvtDataNew*>::iterator iterm;
  for (iterm=_evtMCVec.begin(); iterm!=_evtMCVec.end(); ++iterm){
           double intensity=calcEvtIntensity((*iterm), theParamVal);
           LH_mc+=intensity;
         }

  double logLH_mc_Norm=0.;
  if (LH_mc>0.) logLH_mc_Norm=log(LH_mc/_evtMCVec.size());

  logLH=0.5*weightSum *(LH_mc/_evtMCVec.size()-1.)*(LH_mc/_evtMCVec.size()-1.)
    -logLH_data
    +weightSum*logLH_mc_Norm;

  Info << "current LH = " << logLH << endmsg;

 return logLH;

}

void AbsLhNew::setHyps( const std::map<const std::string, bool>& theMap, bool& theHyp, std::string& theKey){

  std::map<const std::string, bool>::const_iterator iter= theMap.find(theKey);
  
  if (iter !=theMap.end()){
    theHyp= iter->second;
    DebugMsg<< "hypothesis " << iter->first << "\t" << theHyp <<endmsg;
    _hypMap[iter->first]= iter->second;
  }
  else{
    Alert << theKey << " does not exist!!!" <<endmsg;
    exit(0);
  }
}