Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <getopt.h>
#include <fstream>
#include <string>
#include "PwaUtils/AbsLh.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "ErrLogger/ErrLogger.hh"
AbsLh::AbsLh(boost::shared_ptr<const EvtDataBaseList> theEvtList) :
_evtListPtr(theEvtList)
{
_evtDataVec=_evtListPtr->getDataVecs();
_evtMCVec=_evtListPtr->getMcVecs();
}
AbsLh::AbsLh(boost::shared_ptr<AbsLh> theAbsLhPtr):
_evtListPtr(theAbsLhPtr->getEventList())
{
_evtDataVec=_evtListPtr->getDataVecs();
_evtMCVec=_evtListPtr->getMcVecs();
}
AbsLh::~AbsLh()
{
}
double AbsLh::calcLogLh(fitParams& theParamVal){
double logLH=0.;
double logLH_data=0.;
double weightSum=0.;
std::vector<EvtData*>::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<EvtData*>::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=weightSum *(LH_mc/_evtMCVec.size()-1)*(LH_mc/_evtMCVec.size()-1)
-2.*logLH_data
+2.*weightSum*logLH_mc_Norm;
Info << "current LH = " << logLH << endmsg;
return logLH;
}