Skip to content
Snippets Groups Projects
Commit bff3832e authored by Julian Pychy's avatar Julian Pychy
Browse files

omega pi0 pwa: added event weight, new histograms

parent 7a43c14f
No related branches found
No related tags found
No related merge requests found
......@@ -7,9 +7,10 @@
CBElsaReader::CBElsaReader()
{}
CBElsaReader::CBElsaReader(const std::vector<std::string>& files, int particles, int skip):
CBElsaReader::CBElsaReader(const std::vector<std::string>& files, int particles, int skip, bool useWeight):
numParticles(particles),
linesToSkip(skip)
linesToSkip(skip),
_useWeight(useWeight)
{
if (0 == files.size()) {
Alert << "empty list of event files" << endmsg;
......@@ -38,10 +39,19 @@ bool CBElsaReader::fillAll(EventList& evtList)
double e,px,py,pz;
Event* newEvent = new Event();
int parts;
if(_useWeight)
{
double weight;
currentStream >> weight;
newEvent->addWeight(weight);
}
for (parts = 0; parts < numParticles; parts++) {
currentStream >> e >> px >> py >> pz;
newEvent->addParticle(e*1.e-3,px*1.e-3,py*1.e-3,pz*1.e-3);
}
if (!currentStream.fail()) {
evtList.add(newEvent);
for (parts = 0; parts < linesToSkip; parts++)
......
......@@ -16,7 +16,7 @@ class CBElsaReader : public EventReader
{
public:
CBElsaReader();
CBElsaReader(const std::vector<std::string>& files, int particles, int skip);
CBElsaReader(const std::vector<std::string>& files, int particles, int skip, bool useWeight=false);
virtual ~CBElsaReader();
virtual bool fillAll(EventList& evtList);
......@@ -27,6 +27,7 @@ private:
std::ifstream currentStream;
int numParticles;
int linesToSkip;
bool _useWeight;
};
#endif
#include "Event/Event.hh"
#include <vector>
Event::Event()
Event::Event() : weight(1.)
{}
Event::~Event()
......@@ -43,6 +43,11 @@ void Event::addParticle(double e, double px, double py, double pz)
return;
}
void Event::addWeight(double theWeight)
{
weight = theWeight;
}
int Event::size()
{
return particles.size();
......
......@@ -26,12 +26,14 @@ public:
~Event();
void addParticle(double e, double px, double py, double pz);
void addWeight(double);
Vector4<float>* p4(unsigned int i);
float* pid(unsigned int i);
int size();
double Weight(){ return weight; }
private:
std::vector<EvtPartData*> particles;
double weight;
......
......@@ -50,11 +50,13 @@ double AbsOmegaPiLhLS::calcLogLh(const OmegaPiDataLS::fitParamVal& theParamVal){
double logLH=0.;
double logLH_data=0.;
double weightSum=0.;
std::vector<OmegaPiDataLS::OmPiEvtDataLS*>::iterator iterd;
for (iterd=_evtDataVec.begin(); iterd!=_evtDataVec.end(); ++iterd){
double intensity=calcEvtIntensity((*iterd), theParamVal);
if (intensity>0.) logLH_data+=log10(intensity);
if (intensity>0.) logLH_data+= ((*iterd)->eventWeight) * log10(intensity);
weightSum+= (*iterd)->eventWeight;
}
double LH_mc=0.;
......@@ -68,9 +70,9 @@ double AbsOmegaPiLhLS::calcLogLh(const OmegaPiDataLS::fitParamVal& theParamVal){
double logLH_mc_Norm=0.;
if (LH_mc>0.) logLH_mc_Norm=log10(LH_mc/_evtMCVec.size());
logLH=_evtDataVec.size()/2.*(LH_mc/_evtMCVec.size()-1)*(LH_mc/_evtMCVec.size()-1)
-logLH_data
+_evtDataVec.size()*logLH_mc_Norm;
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;
......
......@@ -41,7 +41,7 @@ namespace OmegaPiDataLS {
Vector4<float> pi0HeliOmega4Vec;
Vector4<float> pi0HeliOmega4Vec2;
float cosPi0HeliOmega4Vec;
float eventWeight;
map<Spin,map<Spin,map<Spin,complex<double> > > > Dfp; //Wigner D functions for omega pi0 production
map<Spin,map<Spin,map<Spin,complex<double> > > > Dfd; //Wigner D functions for omega decay to pi0 gamma
......
......@@ -41,8 +41,8 @@ void OmegaPiEventListLS::read4Vecs(EventList& evtList, std::vector<OmPiEvtDataLS
for(int i=0; i<2; ++i){
if ( fabs(anEvent->p4(i)->Mass()-0.78195)<0.01 ) omega_4V=*(anEvent->p4(i));
else if ( fabs(anEvent->p4(i)->Mass()-0.13497)<0.01 ) piRec_4V=*(anEvent->p4(i));
if ( fabs(anEvent->p4(i)->Mass()-0.78195)<0.041 ) omega_4V=*(anEvent->p4(i));
else if ( fabs(anEvent->p4(i)->Mass()-0.13497)<0.041 ) piRec_4V=*(anEvent->p4(i));
else {
Alert <<"this is neither an omega nor a pi0 particle!!!" << endmsg;
exit(1);
......@@ -68,7 +68,7 @@ void OmegaPiEventListLS::read4Vecs(EventList& evtList, std::vector<OmPiEvtDataLS
}
Vector4<float> pi0FromOmega4V=*(anEvent->p4(2));
if ( fabs(pi0FromOmega4V.M()-0.13497)>0.01 ) {
if ( fabs(pi0FromOmega4V.M()-0.13497)>0.041 ) {
Alert <<"the third particle is not the pi0 from the omega decay" << endmsg;
exit(1);
}
......@@ -95,6 +95,8 @@ void OmegaPiEventListLS::read4Vecs(EventList& evtList, std::vector<OmPiEvtDataLS
theOmPiEvtData->pi0HeliOmega4Vec=pi0HeliOmega4V;
theOmPiEvtData->pi0HeliOmega4Vec2=pi0HeliOmega4V2;
theOmPiEvtData->cosPi0HeliOmega4Vec=costDecHeli(cm_4V, omega_4V, pi0FromOmega4V);
theOmPiEvtData->eventWeight = anEvent->Weight();
for (Spin j=0; j<=_jmax; j++){
for (Spin M=-1; M<=1; M++){
......
......@@ -17,14 +17,25 @@ OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<const AbsOmegaPiEventListLS> theE
_cosOmegaHeliMcHist(0),
_cosOmegaHeliFittedHist(0),
_cosOmegaAccCorHist(0),
_cosOmegaFittedAccCorHist(0),
_cosPi0FromOmegaDataHeli(0),
_cosPi0FromOmegaMcHeli(0),
_cosPi0FromOmegaFittedHeli(0),
_cosPi0FromOmegaAccCorHeli(0),
_cosPi0FromOmegaFittedAccCorHeli(0),
_treimanYangDataHist(0),
_treimanYangMcHist(0),
_treimanYangFittedHist(0),
_cosPi0FromOmegaDataHeli1(0)
_treimanYangAccCorHist(0),
_cosPi0FromOmegaDataHeli1(0),
_thetaPhiPi0FromOmegaDataHeli(0),
_thetaPhiPi0FromOmegaMcHeli(0),
_thetaPhiPi0FromOmegaFittedHeli(0),
_thetaPhiPi0FromOmegaAccCorHeli(0),
_thetaPhiPi0FromOmegaFittedAccCorHeli(0),
_prodDecThetaPhiPi0FromOmegaDataHeli(0),
_prodDecThetaPhiPi0FromOmegaMcHeli(0),
_prodDecThetaPhiPi0FromOmegaFittedHeli(0)
{
if(0==theEvtList){
Alert <<"OmegaPiEventList* theEvtList is a 0 pointer !!!!" << endmsg;
......@@ -36,40 +47,47 @@ OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<const AbsOmegaPiEventListLS> theE
initRootStuff(thePathToRootFile);
const std::vector<OmPiEvtDataLS*> dataList=theEvtList->getDataVecs();
_cosOmegaHeliDataHist->Sumw2();
_cosOmegaHeliMcHist->Sumw2();
_cosPi0FromOmegaDataHeli->Sumw2();
_cosPi0FromOmegaMcHeli->Sumw2();
_treimanYangDataHist->Sumw2();
_thetaPhiPi0FromOmegaDataHeli->Sumw2();
_thetaPhiPi0FromOmegaMcHeli->Sumw2();
_prodDecThetaPhiPi0FromOmegaDataHeli->Sumw2();
_prodDecThetaPhiPi0FromOmegaMcHeli->Sumw2();
const std::vector<OmPiEvtDataLS*> dataList=theEvtList->getDataVecs();
std::vector<OmPiEvtDataLS*>::const_iterator it=dataList.begin();
while(it!=dataList.end())
{
plotCosOmegaHeli(_cosOmegaHeliDataHist, (*it), 1.);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaDataHeli, (*it), 1.);
plotTreimanYang(_treimanYangDataHist, (*it), 1.);
{
double dataEventWeight = ((*it)->eventWeight);
plotCosOmegaHeli(_cosOmegaHeliDataHist, (*it), dataEventWeight);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaDataHeli, (*it), dataEventWeight);
plotTreimanYang(_treimanYangDataHist, (*it), dataEventWeight);
plotCosPi0FromOmegaHeli1(_cosPi0FromOmegaDataHeli1, (*it), 1.);
plotThetaPhiPi0FromOmegaHeli(_thetaPhiPi0FromOmegaDataHeli,(*it), dataEventWeight);
plotProdDecThetaPhiPi0FromOmegaHeli(_prodDecThetaPhiPi0FromOmegaDataHeli, (*it), dataEventWeight);
++it;
}
}
const std::vector<OmPiEvtDataLS*> mcList=theEvtList->getMcVecs();
it=mcList.begin();
while(it!=mcList.end())
{
{
plotCosOmegaHeli(_cosOmegaHeliMcHist, (*it), 1.);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaMcHeli, (*it), 1.);
plotTreimanYang(_treimanYangMcHist, (*it), 1.);
plotThetaPhiPi0FromOmegaHeli(_thetaPhiPi0FromOmegaMcHeli,(*it),1.);
plotProdDecThetaPhiPi0FromOmegaHeli(_prodDecThetaPhiPi0FromOmegaMcHeli, (*it), 1.);
++it;
}
_cosOmegaHeliDataHist->Sumw2();
_cosOmegaHeliMcHist->Sumw2();
_cosOmegaAccCorHist=(TH1F*) _cosOmegaHeliDataHist->Clone();
_cosOmegaAccCorHist->Divide(_cosOmegaHeliMcHist);
_cosOmegaAccCorHist->SetName("_cosOmegaAccCorHist");
_cosOmegaAccCorHist->SetTitle("cos(#Theta) #omega heli acc cor");
}
_cosPi0FromOmegaDataHeli->Sumw2();
_cosPi0FromOmegaMcHeli->Sumw2();
_cosPi0FromOmegaAccCorHeli=(TH1F*) _cosPi0FromOmegaDataHeli->Clone();
_cosPi0FromOmegaAccCorHeli->Divide(_cosPi0FromOmegaMcHeli);
_cosPi0FromOmegaAccCorHeli->SetName("_cosPi0FromOmegaAccCorHeli");
_cosPi0FromOmegaAccCorHeli->SetTitle("cos(#Theta) #pi^0 heli (#omega decay) acc cor");
_cosOmegaAccCorHist=doAccCor(_cosOmegaHeliDataHist, _cosOmegaHeliMcHist, "_cosOmegaAccCorHist", "cos(#Theta) #omega heli acc cor");
_cosPi0FromOmegaAccCorHeli=doAccCor(_cosPi0FromOmegaDataHeli, _cosPi0FromOmegaMcHeli, "_cosPi0FromOmegaAccCorHeli", "cos(#Theta) #pi^0 heli (#omega decay) acc cor");
_treimanYangAccCorHist=doAccCor(_treimanYangDataHist, _treimanYangMcHist, "_treimanYangAccCorHist", "_treimanYangAccCorHist");
_thetaPhiPi0FromOmegaAccCorHeli=doAccCor(_thetaPhiPi0FromOmegaDataHeli, _thetaPhiPi0FromOmegaMcHeli,"_thetaPhiPi0FromOmegaAccCorHeli", "_thetaPhiPi0FromOmegaAccCorHeli");
}
OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<AbsOmegaPiLhLS> omegaPiLh, OmegaPiDataLS::fitParamVal& fitParam, const std::string &thePathToRootFile) :
......@@ -78,14 +96,25 @@ OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<AbsOmegaPiLhLS> omegaPiLh, OmegaP
_cosOmegaHeliMcHist(0),
_cosOmegaHeliFittedHist(0),
_cosOmegaAccCorHist(0),
_cosOmegaFittedAccCorHist(0),
_cosPi0FromOmegaDataHeli(0),
_cosPi0FromOmegaMcHeli(0),
_cosPi0FromOmegaFittedHeli(0),
_cosPi0FromOmegaAccCorHeli(0),
_cosPi0FromOmegaFittedAccCorHeli(0),
_treimanYangDataHist(0),
_treimanYangMcHist(0),
_treimanYangFittedHist(0),
_cosPi0FromOmegaDataHeli1(0)
_treimanYangAccCorHist(0),
_cosPi0FromOmegaDataHeli1(0),
_thetaPhiPi0FromOmegaDataHeli(0),
_thetaPhiPi0FromOmegaMcHeli(0),
_thetaPhiPi0FromOmegaFittedHeli(0),
_thetaPhiPi0FromOmegaAccCorHeli(0),
_thetaPhiPi0FromOmegaFittedAccCorHeli(0),
_prodDecThetaPhiPi0FromOmegaDataHeli(0),
_prodDecThetaPhiPi0FromOmegaMcHeli(0),
_prodDecThetaPhiPi0FromOmegaFittedHeli(0)
{
if(0==omegaPiLh){
Alert <<"OmegaPiLh* omegaPiLh is a 0 pointer !!!!" << endmsg;
......@@ -103,15 +132,29 @@ OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<AbsOmegaPiLhLS> omegaPiLh, OmegaP
initRootStuff(thePathToRootFile);
_cosOmegaHeliDataHist->Sumw2();
_cosOmegaHeliMcHist->Sumw2();
_cosPi0FromOmegaDataHeli->Sumw2();
_cosPi0FromOmegaMcHeli->Sumw2();
_treimanYangDataHist->Sumw2();
_thetaPhiPi0FromOmegaDataHeli->Sumw2();
_thetaPhiPi0FromOmegaMcHeli->Sumw2();
_prodDecThetaPhiPi0FromOmegaDataHeli->Sumw2();
_prodDecThetaPhiPi0FromOmegaMcHeli->Sumw2();
std::vector<OmPiEvtDataLS*> dataList=theEvtList->getDataVecs();
std::vector<OmPiEvtDataLS*>::iterator it=dataList.begin();
while(it!=dataList.end())
{
plotCosOmegaHeli(_cosOmegaHeliDataHist, (*it), 1.);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaDataHeli, (*it), 1.);
plotTreimanYang(_treimanYangDataHist, (*it), 1.);
double dataEventWeight = ((*it)->eventWeight);
plotCosOmegaHeli(_cosOmegaHeliDataHist, (*it), dataEventWeight);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaDataHeli, (*it), dataEventWeight);
plotTreimanYang(_treimanYangDataHist, (*it), dataEventWeight);
plotCosPi0FromOmegaHeli1(_cosPi0FromOmegaDataHeli1, (*it), 1.);
plotThetaPhiPi0FromOmegaHeli(_thetaPhiPi0FromOmegaDataHeli,(*it), dataEventWeight);
plotProdDecThetaPhiPi0FromOmegaHeli(_prodDecThetaPhiPi0FromOmegaDataHeli, (*it), dataEventWeight);
++it;
}
......@@ -122,25 +165,21 @@ OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<AbsOmegaPiLhLS> omegaPiLh, OmegaP
plotCosOmegaHeli(_cosOmegaHeliMcHist, (*it), 1.);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaMcHeli, (*it), 1.);
plotTreimanYang(_treimanYangMcHist, (*it), 1.);
plotThetaPhiPi0FromOmegaHeli(_thetaPhiPi0FromOmegaMcHeli,(*it),1.);
plotProdDecThetaPhiPi0FromOmegaHeli(_prodDecThetaPhiPi0FromOmegaMcHeli, (*it), 1.);
double evtWeight= omegaPiLh->calcEvtIntensity((*it), fitParam);
plotCosOmegaHeli(_cosOmegaHeliFittedHist, (*it), evtWeight);
plotCosPi0FromOmegaHeli(_cosPi0FromOmegaFittedHeli, (*it), evtWeight);
plotTreimanYang(_treimanYangFittedHist, (*it), evtWeight);
plotThetaPhiPi0FromOmegaHeli(_thetaPhiPi0FromOmegaFittedHeli, (*it), evtWeight);
plotProdDecThetaPhiPi0FromOmegaHeli(_prodDecThetaPhiPi0FromOmegaFittedHeli, (*it), evtWeight);
++it;
}
_cosOmegaHeliDataHist->Sumw2();
_cosOmegaHeliMcHist->Sumw2();
_cosOmegaAccCorHist=(TH1F*) _cosOmegaHeliDataHist->Clone();
_cosOmegaAccCorHist->Divide(_cosOmegaHeliMcHist);
_cosOmegaAccCorHist->SetName("_cosOmegaAccCorHist");
_cosOmegaAccCorHist->SetTitle("cos(#Theta) #omega heli acc cor");
_cosPi0FromOmegaDataHeli->Sumw2();
_cosPi0FromOmegaMcHeli->Sumw2();
_cosPi0FromOmegaAccCorHeli=(TH1F*) _cosPi0FromOmegaDataHeli->Clone();
_cosPi0FromOmegaAccCorHeli->Divide(_cosPi0FromOmegaMcHeli);
_cosPi0FromOmegaAccCorHeli->SetName("_cosPi0FromOmegaAccCorHeli");
_cosPi0FromOmegaAccCorHeli->SetTitle("cos(#Theta) #pi^0 heli (#omega decay) acc cor");
_cosOmegaAccCorHist=doAccCor(_cosOmegaHeliDataHist, _cosOmegaHeliMcHist, "_cosOmegaAccCorHist", "cos(#Theta) #omega heli acc cor");
_cosPi0FromOmegaAccCorHeli=doAccCor(_cosPi0FromOmegaDataHeli, _cosPi0FromOmegaMcHeli, "_cosPi0FromOmegaAccCorHeli", "cos(#Theta) #pi^0 heli (#omega decay) acc cor");
_treimanYangAccCorHist=doAccCor(_treimanYangDataHist, _treimanYangMcHist, "_treimanYangAccCorHist", "_treimanYangAccCorHist");
_thetaPhiPi0FromOmegaAccCorHeli=doAccCor(_thetaPhiPi0FromOmegaDataHeli, _thetaPhiPi0FromOmegaMcHeli,"_thetaPhiPi0FromOmegaAccCorHeli", "_thetaPhiPi0FromOmegaAccCorHeli");
double integralData=_cosOmegaHeliDataHist->Integral();
......@@ -152,6 +191,13 @@ OmegaPiHistLS::OmegaPiHistLS(boost::shared_ptr<AbsOmegaPiLhLS> omegaPiLh, OmegaP
_cosOmegaHeliFittedHist->Scale(integralData/integralFitted);
_cosPi0FromOmegaFittedHeli->Scale(integralData/integralFitted);
_treimanYangFittedHist->Scale(integralData/integralFitted);
_thetaPhiPi0FromOmegaFittedHeli->Scale(integralData/integralFitted);
_prodDecThetaPhiPi0FromOmegaFittedHeli->Scale(integralData/integralFitted);
_cosOmegaFittedAccCorHist=doAccCor(_cosOmegaHeliFittedHist, _cosOmegaHeliMcHist, "_cosOmegaFittedAccCorHist", "cos(#Theta) #omega heli fitted acc cor");
_cosPi0FromOmegaFittedAccCorHeli=doAccCor(_cosPi0FromOmegaFittedHeli, _cosPi0FromOmegaMcHeli, "_cosPi0FromOmegaFittedAccCorHeli", "cos(#Theta) #pi^0 heli (#omega decay) fitted acc cor");
_treimanYangFittedAccCorHist=doAccCor(_treimanYangFittedHist, _treimanYangMcHist, "_treimanYangFittedAccCorHist", "_treimanYangFittedAccCorHist");
_thetaPhiPi0FromOmegaFittedAccCorHeli=doAccCor(_thetaPhiPi0FromOmegaFittedHeli, _thetaPhiPi0FromOmegaMcHeli,"_thetaPhiPi0FromOmegaFittedAccCorHeli", "_thetaPhiPi0FromOmegaFittedAccCorHeli");
}
OmegaPiHistLS::~OmegaPiHistLS()
......@@ -181,6 +227,36 @@ void OmegaPiHistLS::initRootStuff(const std::string &thePathToRootFile)
_treimanYangMcHist= new TH1F("_treimanYangMcHist","Treiman Yang angle MC",101, -TMath::Pi(), TMath::Pi());
_treimanYangFittedHist= new TH1F("_treimanYangFittedHist","Treiman Yang angle fit result",101, -TMath::Pi(), TMath::Pi());
_cosPi0FromOmegaDataHeli1= new TH1F("_cosPi0FromOmegaDataHeli1","cos(#Theta) #pi^{0} heli (#omega decay) data 1",101, -1.0, 1.0);
_thetaPhiPi0FromOmegaDataHeli = new TH2F("_thetaPhiPi0FromOmegaDataHeli","_thetaPhiPi0FromOmegaDataHeli",11,-1,1,11,-TMath::Pi(),TMath::Pi());
_thetaPhiPi0FromOmegaMcHeli = new TH2F("_thetaPhiPi0FromOmegaMcHeli","_thetaPhiPi0FromOmegaMcHeli",11,-1,1,11,-TMath::Pi(),TMath::Pi());
_thetaPhiPi0FromOmegaFittedHeli = new TH2F("_thetaPhiPi0FromOmegaFittedHeli","_thetaPhiPi0FromOmegaFittedHeli",11,-1,1,11,-TMath::Pi(),TMath::Pi());
_prodDecThetaPhiPi0FromOmegaDataHeli = new TH3F("__prodDecThetaPhiPi0FromOmegaDataHeli",
"__prodDecThetaPhiPi0FromOmegaDataHeli", 8,-1,1,8,-1,1,8,-TMath::Pi(),TMath::Pi());
_prodDecThetaPhiPi0FromOmegaMcHeli = new TH3F("__prodDecThetaPhiPi0FromOmegaMcHeli",
"__prodDecThetaPhiPi0FromOmegaMcHeli", 8,-1,1,8,-1,1,8,-TMath::Pi(),TMath::Pi());
_prodDecThetaPhiPi0FromOmegaFittedHeli = new TH3F("__prodDecThetaPhiPi0FromOmegaFittedHeli",
"__prodDecThetaPhiPi0FromOmegaFittedHeli", 8,-1,1,8,-1,1,8,-TMath::Pi(),TMath::Pi());
}
TH1F* OmegaPiHistLS::doAccCor(TH1F* dataHist, TH1F* mcHist, char* name, char* title)
{
TH1F* accCorHist = (TH1F*)(dataHist->Clone());
accCorHist->Divide(mcHist);
accCorHist->SetName(name);
accCorHist->SetTitle(title);
return accCorHist;
}
TH2F* OmegaPiHistLS::doAccCor(TH2F* dataHist, TH2F* mcHist, char* name, char* title)
{
TH2F* accCorHist = (TH2F*)(dataHist->Clone());
accCorHist->Divide(mcHist);
accCorHist->SetName(name);
accCorHist->SetTitle(title);
return accCorHist;
}
void OmegaPiHistLS::plotCosOmegaHeli(TH1F* theHisto, const OmPiEvtDataLS* theEvtData, double weight)
......@@ -205,3 +281,16 @@ void OmegaPiHistLS::plotTreimanYang(TH1F* theHisto, const OmPiEvtDataLS* theEvtD
Vector4<float> piFromOmegaHeli4V=theEvtData->pi0HeliOmega4Vec;
theHisto->Fill(piFromOmegaHeli4V.Phi(), weight);
}
void OmegaPiHistLS::plotThetaPhiPi0FromOmegaHeli(TH2F* theHisto, const OmPiEvtDataLS* theEvtData, double weight)
{
Vector4<float> piFromOmegaHeli4V=theEvtData->pi0HeliOmega4Vec;
theHisto->Fill(piFromOmegaHeli4V.CosTheta(), piFromOmegaHeli4V.Phi(), weight);
}
void OmegaPiHistLS::plotProdDecThetaPhiPi0FromOmegaHeli(TH3F* theHisto, const OmPiEvtDataLS* theEvtData, double weight)
{
Vector4<float> piFromOmegaHeli4V=theEvtData->pi0HeliOmega4Vec;
Vector4<float> omegaHeli4V=theEvtData->omegaHeliCm4Vec;
theHisto->Fill(omegaHeli4V.CosTheta(),piFromOmegaHeli4V.CosTheta(), piFromOmegaHeli4V.Phi(), weight);
}
......@@ -14,6 +14,7 @@
// #include <TSystem.h>
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "Examples/pbarpToOmegaPiLS/OmegaPiDataLS.hh"
#include "TH3F.h"
// #include "Examples/pbarpToOmegaPiLS/AbsOmegaPiEventListLS.hh"
using OmegaPiDataLS::OmPiEvtDataLS;
......@@ -49,14 +50,32 @@ private:
TH1F* _cosOmegaHeliMcHist;
TH1F* _cosOmegaHeliFittedHist;
TH1F* _cosOmegaAccCorHist;
TH1F* _cosOmegaFittedAccCorHist;
TH1F* _cosPi0FromOmegaDataHeli;
TH1F* _cosPi0FromOmegaMcHeli;
TH1F* _cosPi0FromOmegaFittedHeli;
TH1F* _cosPi0FromOmegaAccCorHeli;
TH1F* _cosPi0FromOmegaFittedAccCorHeli;
TH1F* _treimanYangDataHist;
TH1F* _treimanYangMcHist;
TH1F* _treimanYangFittedHist;
TH1F* _treimanYangAccCorHist;
TH1F* _treimanYangFittedAccCorHist;
TH1F* _cosPi0FromOmegaDataHeli1;
TH2F* _thetaPhiPi0FromOmegaDataHeli;
TH2F* _thetaPhiPi0FromOmegaMcHeli;
TH2F* _thetaPhiPi0FromOmegaFittedHeli;
TH2F* _thetaPhiPi0FromOmegaAccCorHeli;
TH2F* _thetaPhiPi0FromOmegaFittedAccCorHeli;
TH3F* _prodDecThetaPhiPi0FromOmegaDataHeli;
TH3F* _prodDecThetaPhiPi0FromOmegaMcHeli;
TH3F* _prodDecThetaPhiPi0FromOmegaFittedHeli;
unsigned _lmax;
unsigned _pbarmom;
......@@ -65,6 +84,10 @@ private:
void plotCosPi0FromOmegaHeli(TH1F* theHisto, const OmPiEvtDataLS* theEvtData, double weight);
void plotCosPi0FromOmegaHeli1(TH1F* theHisto, const OmPiEvtDataLS* theEvtData, double weight);
void plotTreimanYang(TH1F* theHisto, const OmPiEvtDataLS* theEvtData, double weight);
void plotThetaPhiPi0FromOmegaHeli(TH2F* theHisto, const OmPiEvtDataLS* theEvtData, double weight);
void plotProdDecThetaPhiPi0FromOmegaHeli(TH3F* theHisto, const OmPiEvtDataLS* theEvtData, double weight);
TH1F* doAccCor(TH1F* dataHist, TH1F* mcHist, char* name, char* title);
TH2F* doAccCor(TH2F* dataHist, TH2F* mcHist, char* name, char* title);
};
#endif
......@@ -617,13 +617,20 @@ int main(int __argc,char *__argv[]){
std::string piomegaDatFile;
std::string piomegaMcFile;
int nParticlesPerEvt=0;
bool readWeightData=true;
if (theAppParams.getLhMode()=="OmegaPiLhGamma" || (theAppParams.getLhMode()=="OmegaPiLhGammaBw") ){
nParticlesPerEvt=3;
constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/510_",theAppParams.getPbarMom(),piomegaDatFile);
constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile);
if(theAppParams.getAppExecMode() == 3){ // Alternative EvtGen-Input for spin density calculation
readWeightData=false;
constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/SDM/SDM_",theAppParams.getPbarMom(),piomegaDatFile);
constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/SDM/SDM_",theAppParams.getPbarMom(),piomegaMcFile);
}
else{
constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/newselection/510_",theAppParams.getPbarMom(),piomegaDatFile);
constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/newselection/mc510_",theAppParams.getPbarMom(),piomegaMcFile);
}
}
else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){
nParticlesPerEvt=4;
constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/om3pi_",theAppParams.getPbarMom(),piomegaDatFile);
......@@ -636,29 +643,25 @@ int main(int __argc,char *__argv[]){
}
// constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/510_",theAppParams.getPbarMom(),piomegaDatFile);
if (checkFileExist(piomegaDatFile))
{
// constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/510_",theAppParams.getPbarMom(),piomegaDatFile);
if (checkFileExist(piomegaDatFile)){
Info << "Using Data file " << piomegaDatFile << endmsg;
}
else
{
}
else{
Alert <<"Data file for pbarMom= " << theAppParams.getPbarMom() << " not available !" << endmsg;
Alert << "File " << piomegaDatFile << " is missing !" << endmsg;
exit(1);
}
}
// constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile);
if (checkFileExist(piomegaMcFile))
{
// constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile);
if (checkFileExist(piomegaMcFile)){
Info << "Using Monte Carlo file " << piomegaMcFile << endmsg;
}
else
{
}
else{
Alert <<"Monte Carlo file for pbarMom= " << theAppParams.getPbarMom() << "not available !" << endmsg;
Alert << "File " << piomegaMcFile << " is missing !" << endmsg;
exit(1);
}
}
Info << "data file: " << piomegaDatFile << endmsg;
Info << "mc file: " << piomegaMcFile << endmsg;
......@@ -673,7 +676,6 @@ int main(int __argc,char *__argv[]){
std::vector<std::string> fileNames;
fileNames.push_back(piomegaDatFile);
// CBElsaReader eventReader(fileNames, 3, 0);
EventList theDataEventList;
boost::shared_ptr<CBElsaReader> eventReaderPtr;
......@@ -684,7 +686,7 @@ int main(int __argc,char *__argv[]){
boost::shared_ptr<CBElsaReader> eventReaderMCPtr;
EventList theMcEventList;
eventReaderPtr = boost::shared_ptr<CBElsaReader>(new CBElsaReader(fileNames,nParticlesPerEvt , 0));
eventReaderPtr = boost::shared_ptr<CBElsaReader>(new CBElsaReader(fileNames,nParticlesPerEvt, 0, readWeightData));
eventReaderPtr->fillAll(theDataEventList);
Info << "\nFile has " << theDataEventList.size() << " events. Each event has "
......@@ -712,11 +714,11 @@ int main(int __argc,char *__argv[]){
if (theAppParams.getLhMode()=="OmegaPiLhGamma" || (theAppParams.getLhMode()=="OmegaPiLhGammaBw") ){
theOmegaPiEventPtr = boost::shared_ptr<const AbsOmegaPiEventListLS>(new OmegaPiEventListLS (theDataEventList, theMcEventList, theAppParams.getLMax()+1, theAppParams.getPbarMom() ) );
theRootFilePath << "./" << theAppParams.getName() << "OmegaPi0Fit_lmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".root";
theRootFilePath << "./" << theAppParams.getName() << "OmegaPi0Fit_Lmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".root";
}
else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){
theOmegaPiEventPtr = boost::shared_ptr<const AbsOmegaPiEventListLS>(new OmegaTo3PiEventListLS (theDataEventList, theMcEventList, theAppParams.getLMax()+1, theAppParams.getPbarMom() ) );
theRootFilePath << "./" << theAppParams.getName() << "OmegaTo3PiFit_lmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".root";
theRootFilePath << "./" << theAppParams.getName() << "OmegaTo3PiFit_Lmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".root";
}
if (theAppParams.getAppExecMode()==ApplicationParameterLS::HistTest){
......@@ -760,6 +762,13 @@ int main(int __argc,char *__argv[]){
case ApplicationParameterLS::Minuit:
{
bExecFinish = Minuit(theAppParams, theParamVal, theOmegaPiLh, true);
//now dump the final fit result
std::ostringstream theParamFilePathMin;
theParamFilePathMin << "./" << theAppParams.getName() << "LastFitParamOmegaPi0Fit_Mmax" << theAppParams.getLMax() << "_mom" << theAppParams.getPbarMom() << ".txt";
ofstream outfileMin (theParamFilePathMin.str().c_str());
theOmegaPiLh->dumpCurrentResult(outfileMin, theParamVal);
outfileMin.close();
break;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment