Skip to content
Snippets Groups Projects
Commit bc68f627 authored by Michael Leyhe's avatar Michael Leyhe
Browse files

added Examples/JpsiGamEtaPiPi

parent df11dbfb
No related branches found
No related tags found
No related merge requests found
Showing
with 1652 additions and 0 deletions
#include "Examples/JpsiGamKsKlKK/FitParamErrorMatrix.hh"
using namespace ROOT::Minuit2;
FitParamErrorMatrix::FitParamErrorMatrix( ROOT::Minuit2::MnUserCovariance covMatrix, ROOT::Minuit2::MnUserParameters usrParams){
MnUserTransformation trafo = usrParams.Trafo();
_size = usrParams.Parameters().size();
_errMatrix = boost_matrix(_size,_size);
for(unsigned int irow=0; irow<_size; irow++){
for(unsigned int icol=0; icol<_size; icol++){
_errMatrix(irow,icol)=0.0;
}
}
for(unsigned int row=0;row<covMatrix.Nrow(); row++){
for(unsigned int col=0;col<covMatrix.Nrow(); col++){
double element = covMatrix(row,col);
int trow = trafo.ExtOfInt(row);
int tcol = trafo.ExtOfInt(col);
_errMatrix(trow,tcol)=element;
}
}
}
FitParamErrorMatrix::~FitParamErrorMatrix( ){
}
FitParamErrorMatrix::FitParamErrorMatrix(std::vector<double> &theData, int size){
_size=size;
_errMatrix = boost_matrix(_size,_size);
unsigned int irow(0), icol(0);
std::vector<double>::iterator dataIter;
for(dataIter=theData.begin(); dataIter!=theData.end(); dataIter++){
double val = *dataIter;
_errMatrix(irow,icol) = val;
if(icol==_size-1){
icol=0;
irow++;
}else{
icol++;
}
}
}
void FitParamErrorMatrix::Print(std::ostream &os){
for(unsigned int irow=0;irow<_size;irow++){
for(unsigned int icol=0;icol<_size;icol++){
os << _errMatrix(irow,icol) << " ";
}
os << std::endl;
}
}
std::vector<double> FitParamErrorMatrix::Data(){
std::vector<double> theData;
for(unsigned int irow=0;irow<_size;irow++){
for(unsigned int icol=0;icol<_size;icol++){
theData.push_back(_errMatrix(irow,icol) );
}
}
return theData;
}
void FitParamErrorMatrix::Write(std::ostream& os){
std::vector<double> data = Data();
std::vector<double>::iterator dataIter = data.begin();
for(; dataIter!=data.end();dataIter++){
os << *dataIter << std::endl;
}
}
#ifndef _FitParamErrorMatrix_H
#define _FitParamErrorMatrix_H
#include <iostream>
#include <fstream>
#include "Minuit2/MnUserCovariance.h"
#include "Minuit2/MnUserParameters.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
typedef boost::numeric::ublas::matrix<double> boost_matrix;
class FitParamErrorMatrix {
public:
FitParamErrorMatrix( ROOT::Minuit2::MnUserCovariance covMatrix, ROOT::Minuit2::MnUserParameters usrParams );
FitParamErrorMatrix( int size){_errMatrix = boost_matrix(size,size);}
FitParamErrorMatrix(std::vector<double> &theData, int size);
virtual ~FitParamErrorMatrix();
void Print(std::ostream& os);
void Write(std::ostream& os);
std::vector<double> Data();
double operator() (int row, int col) {return _errMatrix(row,col);}
private:
boost_matrix _errMatrix;
unsigned int _size;
};
#endif
#include <iostream>
#include <fstream>
#include "FitParamErrorMatrixStreamer.hh"
#include "ErrLogger/ErrLogger.hh"
FitParamErrorMatrixStreamer::FitParamErrorMatrixStreamer( std::string errFile ){
std::ifstream inputStream;
inputStream.open(errFile.c_str());
if (!inputStream) {
Alert << "can not open " << errFile ;
exit(1);
}
while (!inputStream.eof()) {
double err;
inputStream >> err;
if (!inputStream.fail()) {
_matrixData.push_back(err);
}
}
}
#ifndef _FitParamErrorMatrixStreamer_H
#define _FitParamErrorMatrixStreamer_H
#include <string>
#include <vector>
class FitParamErrorMatrixStreamer {
public :
FitParamErrorMatrixStreamer( std::string errFile );
virtual ~FitParamErrorMatrixStreamer(){};
void matrixData( std::vector<double>& theData, int& nrows ){
theData=_matrixData;
int sizeOfData=(int) theData.size();
nrows=1;
while(nrows*nrows!=sizeOfData){
nrows++;
}
}
private :
std::vector<double> _matrixData;
};
#endif
#include "Examples/JpsiGamEtaPiPi/FitParamIndex.hh"
FitParamIndex::FitParamIndex(fitParams& theParams){
int indexCounter=0;
//iterate amplitude magnitueds
std::map<int, std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > >& startMagMap=theParams.Mags;
std::map<int, std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > >::iterator itMagMap;
for (itMagMap=startMagMap.begin(); itMagMap!=startMagMap.end(); ++itMagMap){
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >::map jpclsMap = itMagMap->second;
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >::iterator jpclsMapIter;
for (jpclsMapIter=jpclsMap.begin(); jpclsMapIter!=jpclsMap.end(); ++jpclsMapIter){
_Mags[itMagMap->first][jpclsMapIter->first] = indexCounter;
indexCounter++;
}
}
//loop again and fill phi index
for (itMagMap=startMagMap.begin(); itMagMap!=startMagMap.end(); ++itMagMap){
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >::map jpclsMap = itMagMap->second;
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >::iterator jpclsMapIter;
for (jpclsMapIter=jpclsMap.begin(); jpclsMapIter!=jpclsMap.end(); ++jpclsMapIter){
_Phis[itMagMap->first][jpclsMapIter->first] = indexCounter;
indexCounter++;
}
}
}
#ifndef _FitParamIndex_H
#define _FitParamIndex_H
#include "PwaUtils/FitParamsBase.hh"
class FitParamIndex {
public:
FitParamIndex( fitParams& theParams );
virtual ~FitParamIndex(){};
int Mag(int amp, boost::shared_ptr<const JPCLS> state){ return _Mags[amp][state]; }
int Phi(int amp, boost::shared_ptr<const JPCLS> state){ return _Phis[amp][state]; }
std::map<int, std::map< boost::shared_ptr<const JPCLS>, int, pawian::Collection::SharedPtrLess > > _Mags;
std::map<int, std::map< boost::shared_ptr<const JPCLS>, int, pawian::Collection::SharedPtrLess > > _Phis;
};
#endif
project :
;
lib JpsiGamEtaPiPi :
[ glob *.cc : *App.cc ]
$(TOP)/PwaUtils//PwaUtils
$(TOP)/Event//Event
$(TOP)/qft++//qft++
: <use>$(TOP)//Geneva
:
: <library>$(TOP)//Geneva ;
exe MJpsiGamEtaPiPiApp : MJpsiGamEtaPiPiApp.cc JpsiGamEtaPiPi : ;
exe MJpsiGamEtaPiPiParserTestApp : JpsiGamEtaPiPiParserTestApp.cc JpsiGamEtaPiPi : ;
project :
;
lib JpsiGamKlKsKK :
[ glob *.cc : *App.cc ]
$(TOP)/PwaUtils//PwaUtils
$(TOP)/Event//Event
$(TOP)/qft++//qft++
: <use>$(TOP)//Geneva
:
: <library>$(TOP)//Geneva ;
exe MJpsiGamKlKsKKApp : MJpsiGamKlKsKKApp.cc JpsiGamKlKsKK : ;
exe MJpsiGamKsKlKKParserTestApp : JpsiGamKsKlKKParserTestApp.cc JpsiGamKlKsKK : ;
#################################################################################
# Configuration file for JpsiToEtaPiPi PWA #
# #
# Authors: Bertram Kopf (RUB) #
#################################################################################
errLogMode = debug
datFile = /data/liema/michaell/PwaJpsiGamEtaPiPi/FS4Vectors-Data-JpsiGamEtaPipPim.dat
mcFile = /data/liema/michaell/PwaJpsiGamEtaPiPi/FS4Vectors-PHSPMC-JpsiGamEtaPipPim.dat
paramFile = /data/liema/michaell/Pawian/Examples/JpsiGamEtaPiPi/startParamBase.dat
startHypo = production
enableHyp = f11285Hyp
enableHyp = eta1295Hyp
enableHyp = eta1405Hyp
#enableHyp = f02020FlatteHyp
#enableHyp = f22300Hyp
#enableHyp = eta21870Hyp
#enableHyp = f1Hyp
enableHyp = usePhasespace
mnParFix = J1P-1C-1L1S1PsiToEtacGamPhi
massRangeMin = 1.28
massRangeMax = 1.30
mode = qaMode
# Determines whether information about configuration options should be emitted
verbose = 1
#ifndef _JpsiGamEtaPiPiData_H
#define _JpsiGamEtaPiPiData_H
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include <boost/shared_ptr.hpp>
#include <map>
#include "Utils/PawianCollectionUtils.hh"
#include "PwaUtils/DataUtils.hh"
struct enumJpsiGamEtaPiPiData{
enum {V4_Psi=0, V4_EtaPipPim_HeliPsi, V4_PipPim_HeliPsi, V4_EtaPip_HeliPsi, V4_EtaPim_HeliPsi,
V4_gamma_HeliPsi, V4_Eta_HeliPsi, V4_Pip_HeliPsi, V4_Pim_HeliPsi,
V4_Eta_HeliEtaPipPim, V4_Pip_HeliEtaPipPim, V4_Pim_HeliEtaPipPim, V4_EtaPip_HeliEtaPipPim, V4_EtaPim_HeliEtaPipPim, V4_PipPim_HeliEtaPipPim,
V4_Pip_HeliPipPim, V4_Pim_HeliPipPim, V4_Pip_HeliEtaPip, V4_Pim_HeliEtaPim, V4_Eta_HeliEtaPip, V4_Eta_HeliEtaPim,
V4_normPipPimDecHeliEtaPipPim, V4_normEtaPipDecHeliEtaPipPim,
n4Vecs, Df_Psi=n4Vecs, Df_Spin0, Df_Spin2, Df_KsKl, Df_KpKm, nDfuncts
};
static const std::string& name(unsigned int t)
{
static std::string fitName[enumJpsiGamEtaPiPiData::nDfuncts]
={"Psi", "EtaPipPim_HeliPsi", "PipPim_HeliPsi","EtaPip_HeliPsi","EtaPim_HeliPsi",
"gamma_HeliPsi", "Eta_HeliPsi", "Pip_HeliPsi", "Pim_HeliPsi",
"Eta_HeliEtaPipPim", "Pip_HeliEtaPipPim", "Pim_HeliEtaPipPim", "EtaPip_HeliEtaPipPim", "EtaPim_HeliEtaPipPim", "PipPim_HeliEtaPipPim",
"Pip_HeliPipPim", "Pim_HeliPipPim", "Pip_HeliEtaPip", "Pim_HeliEtaPim", "Eta_HeliEtaPip", "Eta_HeliEtaPim",
"Df_Psi","Df_Spin0","Df_Spin2", "Df_KsKl","Df_KpKm"
};
if (t<0 || t>=enumJpsiGamEtaPiPiData::nDfuncts) assert(0);
return fitName[t];
}
};
#endif /* _JpsiGamEtaPiPiData_H */
#include <utility>
struct dynamicModelParams{
enum enumDynamicModel {BreitWigner, BreitWignerBlattW, Flatte, MassIndependent};
enumDynamicModel dynamicModel;
double mass;
double width;
double gFactor1;
double gFactor2;
pair<double, double> massPair1;
pair<double, double> massPair2;
complex<double> value;
};
#include <getopt.h>
#include <fstream>
#include <string>
#include "Examples/JpsiGamEtaPiPi/JpsiGamEtaPiPiEventList.hh"
#include "Examples/JpsiGamEtaPiPi/JpsiGamEtaPiPiData.hh"
#include "Event/EventList.hh"
#include "PwaUtils/KinUtils.hh"
#include "Event/Event.hh"
#include "ErrLogger/ErrLogger.hh"
JpsiGamEtaPiPiEventList::JpsiGamEtaPiPiEventList(EventList& evtListData, EventList& evtListMc)
{
read4Vecs(evtListData, _evtDataList);
read4Vecs(evtListMc, _mcDataList);
}
JpsiGamEtaPiPiEventList::~JpsiGamEtaPiPiEventList()
{
}
void JpsiGamEtaPiPiEventList::read4Vecs(EventList& evtList, std::vector<EvtData*>& theEvtList){
Event* anEvent;
int evtCount = 0;
while ((anEvent = evtList.nextEvent())){
if (evtCount%10000 == 0) Info << "4vec calculation for event " << evtCount ; // << endmsg;
Vector4<float> gam = *(anEvent->p4(0));
Vector4<float> eta = *(anEvent->p4(1));
Vector4<float> pip = *(anEvent->p4(2));
Vector4<float> pim = *(anEvent->p4(3));
Vector4<float> V4_psi = gam+eta+pip+pim;
if (evtCount%10000 == 0){
Info << "psi 4vec" << "\n"
<< " px: " << V4_psi.Px() <<"\t"
<< " py: " << V4_psi.Py() <<"\t"
<< " pz: " << V4_psi.Pz() <<"\t"
<< " e : " << V4_psi.E() << "\t"
<< " m : " << V4_psi.M() ; // << endmsg;
}
Vector4<float> V4_all_Lab(gam+eta+pip+pim);
Vector4<float> V4_EtaPipPim_Lab(eta+pip+pim);
Vector4<float> V4_EtaPip_Lab(eta+pip);
Vector4<float> V4_EtaPim_Lab(eta+pim);
Vector4<float> V4_PipPim_Lab(pip+pim);
Vector4<float> V4_Pip_Lab(pip);
Vector4<float> V4_Pim_Lab(pim);
Vector4<float> V4_Eta_Lab(eta);
Vector4<float> V4_EtaPipPim_HeliPsi(eta+pip+pim);
V4_EtaPipPim_HeliPsi.Boost(V4_psi);
Vector4<float> V4_gamma_HeliPsi(gam);
V4_gamma_HeliPsi.Boost(V4_psi);
Vector4<float> V4_PipPim_HeliPsi(pip+pim);
V4_PipPim_HeliPsi.Boost(V4_psi);
Vector4<float> V4_EtaPip_HeliPsi(eta+pip);
V4_EtaPip_HeliPsi.Boost(V4_psi);
Vector4<float> V4_EtaPim_HeliPsi(eta+pim);
V4_EtaPim_HeliPsi.Boost(V4_psi);
Vector4<float> V4_Eta_HeliPsi(eta);
V4_Eta_HeliPsi.Boost(V4_psi);
Vector4<float> V4_Pip_HeliPsi(pip);
V4_Pip_HeliPsi.Boost(V4_psi);
Vector4<float> V4_Pim_HeliPsi(pim);
V4_Pim_HeliPsi.Boost(V4_psi);
Vector4<float> V4_Eta_HeliEtaPipPim = helicityVec(V4_all_Lab, V4_EtaPipPim_Lab, V4_Eta_Lab);
Vector4<float> V4_Pip_HeliEtaPipPim = helicityVec(V4_all_Lab, V4_EtaPipPim_Lab, V4_Pip_Lab);
Vector4<float> V4_Pim_HeliEtaPipPim = helicityVec(V4_all_Lab, V4_EtaPipPim_Lab, V4_Pim_Lab);
Vector4<float> V4_EtaPip_HeliEtaPipPim = helicityVec(V4_all_Lab, V4_EtaPipPim_Lab, V4_EtaPip_Lab);
Vector4<float> V4_EtaPim_HeliEtaPipPim = helicityVec(V4_all_Lab, V4_EtaPipPim_Lab, V4_EtaPim_Lab);
Vector4<float> V4_PipPim_HeliEtaPipPim = helicityVec(V4_all_Lab, V4_EtaPipPim_Lab, V4_PipPim_Lab);
Vector4<float> V4_Pip_HeliPipPim = helicityVec(V4_all_Lab, V4_PipPim_Lab, V4_Pip_Lab);
Vector4<float> V4_Pim_HeliPipPim = helicityVec(V4_all_Lab, V4_PipPim_Lab, V4_Pim_Lab);
Vector4<float> V4_Pip_HeliEtaPip = helicityVec(V4_all_Lab, V4_EtaPip_Lab, V4_Pip_Lab);
Vector4<float> V4_Pim_HeliEtaPim = helicityVec(V4_all_Lab, V4_EtaPim_Lab, V4_Pim_Lab);
Vector4<float> V4_Eta_HeliEtaPip = helicityVec(V4_all_Lab, V4_EtaPip_Lab, V4_Eta_Lab);
Vector4<float> V4_Eta_HeliEtaPim = helicityVec(V4_all_Lab, V4_EtaPim_Lab, V4_Eta_Lab);
/*
Vector4<float> V4_EtaPip_HeliEtaPipPim=helicityVec(V4_all_Lab, V4_KsKlKpKm_Lab, V4_KsKl_Lab);
Vector4<float> V4_KpKm_HeliKsKlKpKm=helicityVec(V4_all_Lab, V4_KsKlKpKm_Lab, V4_KpKm_Lab);
Vector4<float> V4_Ks_HeliKsKl=helicityVec(V4_KsKlKpKm_Lab, V4_KsKl_Lab, V4_Ks_Lab);
Vector4<float> V4_Kp_HeliKpKm=helicityVec(V4_KsKlKpKm_Lab, V4_KpKm_Lab, V4_Kp_Lab);
Vector4<float> V4_Kp_HeliKsKlKpKm=helicityVec(V4_all_Lab, V4_KsKlKpKm_Lab, V4_Kp_Lab);
Vector4<float> V4_Km_HeliKsKlKpKm=helicityVec(V4_all_Lab, V4_KsKlKpKm_Lab, V4_Km_Lab);
*/
Vector4<float> V4_normPipPimDecHeliEtaPipPim(0.5*(V4_Pip_HeliEtaPipPim.T()+V4_Pim_HeliEtaPipPim.T()),
V4_Pim_HeliEtaPipPim.Y()*V4_Pip_HeliEtaPipPim.Z()-V4_Pim_HeliEtaPipPim.Z()*V4_Pip_HeliEtaPipPim.Y(),
V4_Pim_HeliEtaPipPim.Z()*V4_Pip_HeliEtaPipPim.X()-V4_Pim_HeliEtaPipPim.X()*V4_Pip_HeliEtaPipPim.Z(),
V4_Pim_HeliEtaPipPim.X()*V4_Pip_HeliEtaPipPim.Y()-V4_Pim_HeliEtaPipPim.Y()*V4_Pip_HeliEtaPipPim.X());
Vector4<float> V4_normEtaPipDecHeliEtaPipPim(0.5*(V4_Eta_HeliEtaPipPim.T()+V4_Pip_HeliEtaPipPim.T()),
V4_Pip_HeliEtaPipPim.Y()*V4_Eta_HeliEtaPipPim.Z()-V4_Pip_HeliEtaPipPim.Z()*V4_Eta_HeliEtaPipPim.Y(),
V4_Pip_HeliEtaPipPim.Z()*V4_Eta_HeliEtaPipPim.X()-V4_Pip_HeliEtaPipPim.X()*V4_Eta_HeliEtaPipPim.Z(),
V4_Pip_HeliEtaPipPim.X()*V4_Eta_HeliEtaPipPim.Y()-V4_Pip_HeliEtaPipPim.Y()*V4_Eta_HeliEtaPipPim.X());
EvtData* evtData=new EvtData();
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Psi] = V4_psi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_EtaPipPim_HeliPsi] = V4_EtaPipPim_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_PipPim_HeliPsi] = V4_PipPim_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_EtaPip_HeliPsi] = V4_EtaPip_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_EtaPim_HeliPsi] = V4_EtaPim_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_gamma_HeliPsi] = V4_gamma_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Eta_HeliPsi] = V4_Eta_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pip_HeliPsi] = V4_Pip_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pim_HeliPsi] = V4_Pim_HeliPsi;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Eta_HeliEtaPipPim] = V4_Eta_HeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pip_HeliEtaPipPim] = V4_Pip_HeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pim_HeliEtaPipPim] = V4_Pim_HeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_EtaPip_HeliEtaPipPim] = V4_EtaPip_HeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_EtaPim_HeliEtaPipPim] = V4_EtaPim_HeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_PipPim_HeliEtaPipPim] = V4_PipPim_HeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pip_HeliPipPim] = V4_Pip_HeliPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pim_HeliPipPim] = V4_Pim_HeliPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pip_HeliEtaPip] = V4_Pip_HeliEtaPip;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Pim_HeliEtaPim] = V4_Pim_HeliEtaPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Eta_HeliEtaPip] = V4_Eta_HeliEtaPip;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_Eta_HeliEtaPim] = V4_Eta_HeliEtaPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_normPipPimDecHeliEtaPipPim] = V4_normPipPimDecHeliEtaPipPim;
evtData->FourVecs[enumJpsiGamEtaPiPiData::V4_normEtaPipDecHeliEtaPipPim] = V4_normEtaPipDecHeliEtaPipPim;
// calculate and store WignerD functions for Psi -> X gamma (JPC =0-=, 0++, 2++)
Spin jPsi=1;
for (Spin M=-1; M<=1; M=M+2){
for (Spin lam=-1; lam<=1; lam++){
evtData->WignerDs[enumJpsiGamEtaPiPiData::Df_Psi][jPsi][M][lam]=Wigner_D(0.,V4_EtaPipPim_HeliPsi.Theta(),0,jPsi,M,lam);
}
}
/* //WignerD function for 2+ -> phi phi
Spin jTensor =2;
for(Spin M=-jTensor; M<=jTensor; M++){
for (Spin lam=-2; lam<=2; lam++){
evtData->WignerDs[enumJpsiGamEtaPiPiData::Df_Spin2][jTensor][M][lam] = Wigner_D(V4_KsKl_HeliKsKlKpKm.Phi(),V4_KsKl_HeliKsKlKpKm.Theta(), 0,jTensor,M,lam);
}
}
//WignerD function for 1+ -> phi phi
Spin jAxialVec =1;
for(Spin M=-jAxialVec; M<=jAxialVec; M++){
for (Spin lam=-1; lam<=1; lam++){
evtData->WignerDs[enumJpsiGamEtaPiPiData::Df_Spin2][jAxialVec][M][lam] = Wigner_D(V4_KsKl_HeliKsKlKpKm.Phi(),V4_KsKl_HeliKsKlKpKm.Theta(), 0,jAxialVec,M,lam); //use Df_Spin2 (!) to be able to use same function for fJ(J=1,2) decay amplitudes
}
}
//WignerD function for 0-/0+ -> phi phi
Spin jScalar =0;
Spin M =0;
Spin lam =0;
evtData->WignerDs[enumJpsiGamEtaPiPiData::Df_Spin0][jScalar][M][lam] = Wigner_D(V4_KsKl_HeliKsKlKpKm.Phi(),V4_KsKl_HeliKsKlKpKm.Theta(), 0,jScalar,M,lam);
//WignerD function for phi -> K+ K- and phi -> KS KL
Spin phiSpin=1;
for(Spin M=-phiSpin; M<=phiSpin; M++){
Spin lam=0;
evtData->WignerDs[enumJpsiGamEtaPiPiData::Df_KsKl][phiSpin][M][lam] = Wigner_D(V4_Ks_HeliKsKl.Phi(),V4_Ks_HeliKsKl.Theta(), 0,phiSpin,M,lam);
evtData->WignerDs[enumJpsiGamEtaPiPiData::Df_KpKm][phiSpin][M][lam] = Wigner_D(V4_Kp_HeliKpKm.Phi(),V4_Kp_HeliKpKm.Theta(), 0,phiSpin,M,lam);
}
*/
evtData->evtWeight=1.;
theEvtList.push_back(evtData);
++evtCount;
}
}
#ifndef _JpsiGamEtaPiPiEventList_H
#define _JpsiGamEtaPiPiEventList_H
#include <iostream>
#include <vector>
#include <cassert>
// #include <TSystem.h>
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "PwaUtils/EvtDataBaseList.hh"
class EventList;
class JpsiGamEtaPiPiEventList : public EvtDataBaseList {
public:
// create/copy/destroy:
///Constructor
JpsiGamEtaPiPiEventList(EventList& evtListData, EventList& evtListMc);
/** Destructor */
virtual ~JpsiGamEtaPiPiEventList();
// Getters:
protected:
void read4Vecs(EventList& evtList, std::vector<EvtData*>& theEvtList);
private:
};
#endif
#include <getopt.h>
#include <fstream>
#include <string>
#include "Examples/JpsiGamEtaPiPi/JpsiGamEtaPiPiFitParams.hh"
#include "Examples/JpsiGamEtaPiPi/JpsiGamEtaPiPiStates.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "ErrLogger/ErrLogger.hh"
JpsiGamEtaPiPiFitParams::JpsiGamEtaPiPiFitParams()
{
JpsiGamEtaPiPiStates theStates;
theStates.print(std::cout);
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToEta1405Gamma]=theStates.PsiToEtaGammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToEta1295Gamma]=theStates.PsiToEtaGammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF11285Gamma]=theStates.PsiToF1GammaStates();
/*
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF22010Gamma]=theStates.PsiToF2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF22300Gamma]=theStates.PsiToF2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF22340Gamma]=theStates.PsiToF2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF1Gamma]=theStates.PsiToF1GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToEta21870Gamma]=theStates.PsiToEta2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::F02020ToPhiPhi]=theStates.F0ToPhiPhiStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::F22300ToPhiPhi]=theStates.F2ToPhiPhiStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::Eta21870ToPhiPhi]=theStates.Eta2ToPhiPhiStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::F1ToPhiPhi]=theStates.F1ToPhiPhiStates();
*/
}
JpsiGamEtaPiPiFitParams::JpsiGamEtaPiPiFitParams(fitParams& theStartparams, fitParams& theErrorparams) :
FitParamsBase(theStartparams, theErrorparams)
{
JpsiGamEtaPiPiStates theStates;
theStates.print(std::cout);
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToEta1405Gamma]=theStates.PsiToEtaGammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToEta1295Gamma]=theStates.PsiToEtaGammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF11285Gamma]=theStates.PsiToF1GammaStates();
/*
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF22010Gamma]=theStates.PsiToF2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF22300Gamma]=theStates.PsiToF2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF22340Gamma]=theStates.PsiToF2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToF1Gamma]=theStates.PsiToF1GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::PsiToEta21870Gamma]=theStates.PsiToEta2GammaStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::F02020ToPhiPhi]=theStates.F0ToPhiPhiStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::F22300ToPhiPhi]=theStates.F2ToPhiPhiStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::Eta21870ToPhiPhi]=theStates.Eta2ToPhiPhiStates();
_jpclsMap[paramEnumJpsiGamEtaPiPi::F1ToPhiPhi]=theStates.F1ToPhiPhiStates();
*/
}
JpsiGamEtaPiPiFitParams::~JpsiGamEtaPiPiFitParams()
{
}
#ifndef _JpsiGamEtaPiPiFitParams_H
#define _JpsiGamEtaPiPiFitParams_H
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <map>
#include <complex>
#include <cassert>
#include <boost/shared_ptr.hpp>
#include "TROOT.h"
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "PwaUtils/DataUtils.hh"
#include "PwaUtils/FitParamsBase.hh"
#include "Utils/PawianCollectionUtils.hh"
#include "Minuit2/MnUserParameters.h"
// using namespace std;
using namespace ROOT::Minuit2;
struct paramEnumJpsiGamEtaPiPi{
enum { PsiToEta1405Gamma=0, PsiToEta1295Gamma, PsiToF11285Gamma, nAmps,
eta1405=nAmps, eta1295, f11285, nMasses, ngFactors, phaseSpace=ngFactors, nOthers };
/* enum { PsiToEtacGamma=0, PsiToEta2225Gamma, PsiToF02020Gamma, PsiToF22010Gamma, PsiToF22300Gamma, PsiToF22340Gamma, PsiToEta21870Gamma, PsiToF1Gamma,
F02020ToPhiPhi, F22300ToPhiPhi, Eta21870ToPhiPhi, F1ToPhiPhi, nAmps,
etac=nAmps, eta2225, f02020, f22010, f22300 ,f22340, eta21870, f1, nMasses, f02020gKK=nMasses, f02020gPhiPhi, ngFactors,
phaseSpace=ngFactors, nOthers };*/
static const std::string& name(unsigned int t)
{
static std::string fitName[paramEnumJpsiGamEtaPiPi::nOthers]
={"PsiToEta1405Gamma", "PsiToEta1295Gamma", "PsiToF11285Gamma",
"eta1405", "eta1295", "f11285", "phaseSpace"};
if (t<0 || t>=paramEnumJpsiGamEtaPiPi::nOthers) assert(0);
return fitName[t];
}
};
class JpsiGamEtaPiPiFitParams : public FitParamsBase {
public:
JpsiGamEtaPiPiFitParams();
JpsiGamEtaPiPiFitParams(fitParams& theStartparams, fitParams& theErrorparams);
virtual ~JpsiGamEtaPiPiFitParams();
virtual const std::string ampName(int index) {return paramEnumJpsiGamEtaPiPi::name(index);}
virtual const std::string massName(int index) {return paramEnumJpsiGamEtaPiPi::name(index);}
virtual const std::string widthName(int index) {return paramEnumJpsiGamEtaPiPi::name(index);}
virtual const std::string gFactorName(int index) {return paramEnumJpsiGamEtaPiPi::name(index);}
virtual const std::string otherName(int index) {return paramEnumJpsiGamEtaPiPi::name(index);}
virtual int ampIdxMin() {return paramEnumJpsiGamEtaPiPi::PsiToEta1295Gamma;}
virtual int ampIdxMax() {return paramEnumJpsiGamEtaPiPi::nAmps-1;}
virtual int massIdxMin() {return paramEnumJpsiGamEtaPiPi::nAmps;}
virtual int massIdxMax() {return paramEnumJpsiGamEtaPiPi::nMasses-1;}
virtual int gFactorIdxMin() {return paramEnumJpsiGamEtaPiPi::nMasses;}
virtual int gFactorIdxMax() {return paramEnumJpsiGamEtaPiPi::ngFactors-1;}
virtual int otherIdxMin() {return paramEnumJpsiGamEtaPiPi::ngFactors;}
virtual int otherIdxMax() {return paramEnumJpsiGamEtaPiPi::nOthers-1;}
protected:
private:
};
#endif
This diff is collapsed.
#ifndef _JpsiGamEtaPiPiHist_H
#define _JpsiGamEtaPiPiHist_H
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <utility>
#include <cassert>
#include <boost/shared_ptr.hpp>
#include "TROOT.h"
#include <TSystem.h>
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "qft++/topincludes/tensor.hh"
#include "PwaUtils/EvtDataBaseList.hh"
#include "PwaUtils/FitParamsBase.hh"
#include "PwaUtils/AbsLh.hh"
#include "Examples/JpsiGamEtaPiPi/FitParamErrorMatrix.hh"
class TFile;
class TH2F;
class TH1F;
class TNtuple;
class JpsiGamEtaPiPiProdLh;
class EvtDataBaseList;
class FitParamErrorMatrix;
class JpsiGamEtaPiPiHist {
public:
// create/copy/destroy:
///Constructor
JpsiGamEtaPiPiHist(boost::shared_ptr<const EvtDataBaseList>);
//JpsiGamEtaPiPiHist(boost::shared_ptr<AbsLh>, fitParams&);
JpsiGamEtaPiPiHist(JpsiGamEtaPiPiProdLh* theJpsiGamEtaPiPiLh, fitParams& fitParam, FitParamErrorMatrix* theErrorMatrix );
void fill();
void setMassRange(std::pair<double, double> theMassRange){ _massRange = theMassRange; }
/** Destructor */
virtual ~JpsiGamEtaPiPiHist();
// Getters:
protected:
private:
TFile* _theTFile;
/*
TH2F* _dalitzDataHist;
TH2F* _dalitzMcHist;
TH2F* _dalitzFittedHist;
*/
TH1F* _EtaPiPiMassDataHist;
TH1F* _EtaPiPiMassMcHist;
TH1F* _EtaPiPiMassFittedHist;
TH1F* _PipPimMassDataHist;
TH1F* _PipPimMassMcHist;
TH1F* _PipPimMassFittedHist;
TH1F* _EtaPiMassDataHist;
TH1F* _EtaPiMassMcHist;
TH1F* _EtaPiMassFittedHist;
TH1F* _costEta_EtaPipHeliDataHist;
TH1F* _costEta_EtaPipHeliMcHist;
TH1F* _costEta_EtaPipHeliFittedHist;
TH1F* _phiEta_EtaPipHeliDataHist;
TH1F* _phiEta_EtaPipHeliMcHist;
TH1F* _phiEta_EtaPipHeliFittedHist;
TH1F* _costPip_PipPimHeliDataHist;
TH1F* _costPip_PipPimHeliMcHist;
TH1F* _costPip_PipPimHeliFittedHist;
TH1F* _phiPip_PipPimHeliDataHist;
TH1F* _phiPip_PipPimHeliMcHist;
TH1F* _phiPip_PipPimHeliFittedHist;
TH1F* _costGamCmDataHist;
TH1F* _costGamCmMcHist;
TH1F* _costGamCmFittedHist;
/*
TH1F* _costPhi_KpKmDataHist;
TH1F* _costPhi_KpKmMcHist;
TH1F* _costPhi_KpKmFittedHist;
TH1F* _phiPhi_KpKmDataHist;
TH1F* _phiPhi_KpKmMcHist;
TH1F* _phiPhi_KpKmFittedHist;
*/
TH1F* _chiDataHist;
TH1F* _chiMcHist;
TH1F* _chiFittedHist;
TNtuple* _dataTuple;
TNtuple* _mcTuple;
TNtuple* _massIndepTuple;
std::pair<double, double> _massRange;
void initRootStuff();
void plotDalitz(TH2F* theHisto, EvtData* theData, double weight);
void plotEtaPipPim(TH1F* theHisto, EvtData* theData, double weight);
void plotEtaPi(TH1F* theHisto, EvtData* theData, double weight);
void plotPipPim(TH1F* theHisto, EvtData* theData, double weight);
void plotCostPhiEta(TH1F* theCostHisto, TH1F* thePhiHisto, EvtData* theData, double weight);
void plotCostPhiPip(TH1F* theCostHisto, TH1F* thePhiHisto, EvtData* theData, double weight);
// void plotCostPhi_PhiPhiHeli(TH1F* theCostHisto, TH1F* thePhiHisto, const Vector4<double>& the4Vec, double weight);
void plotCostGam(TH1F* theCostHisto, EvtData* theData, double weight);
void plotChi(TH1F* theChiHisto, EvtData* theData, double weight);
void fillTuple( TNtuple* theTuple, EvtData* theData, double weight);
double decayAngleChi(const Vector4<double>& v4_p,const Vector4<double>& v4_d1,
const Vector4<double>& v4_d2,const Vector4<double>& v4_h1,
const Vector4<double>& v4_h2 ) ;
JpsiGamEtaPiPiProdLh* _theJpsiGamEtaPiPiLh;
fitParams _fitParam;
FitParamErrorMatrix* _errMatrix;
};
#endif
// Bertram Kopf (RUB)
#include "Examples/JpsiGamEtaPiPi/JpsiGamEtaPiPiParser.hh"
#include "ErrLogger/ErrLogger.hh"
#include <iterator>
#include <iostream>
#include <fstream>
using namespace std;
/************************************************************************************************/
/**
* A function that parses the command line for all required parameters
*/
bool JpsiGamEtaPiPiParser::parseCommandLine(int argc, char **argv)
{
try
{
bool verbose;
string strErrLogMode="debug";
// Check the command line options. Uses the Boost program options library.
string strAppName(argv[0]);
size_t found = strAppName.rfind("/")+1;
if (found != string::npos) strAppName=strAppName.substr(found);
string strDesc="Usage: " + strAppName + " [options]";
po::options_description desc(strDesc);
desc.add_options()
("help,h", "emit help message")
("configFile,c",po::value<std::string>(&_configFile)->default_value(_configFile),
"The name of the configuration file holding further configuration options")
;
po::options_description common("Common Options");
common.add_options()
("errLogMode,e", po::value<string>(&strErrLogMode)->default_value(strErrLogMode),"choose mode for Error logger.")
("datFile",po::value<string>(&_dataFile), "full path of data file")
("mcFile",po::value<string>(&_mcFile), "full path of Monte Carlo file")
("paramFile",po::value<string>(&_paramFile), "file with start parameters for fit or QA (full path)")
("startHypo",po::value<string>(&_startHypo), "choose the hyopthesis to start")
("enableHyp",po::value< vector<string> >(&_enabledHyps), "enable hypotheses")
("mode",po::value<string>(&_mode), "modes are: pwa, dumpDefaulParams, qaMode, plotmode")
("massIndependentFit", po::value<bool>(&_massIndependentFit), "enable/disable mass independence in fit")
("commonProdPhases",po::value<bool>(&_useCommonProductionPhases), "enable/disable common production phases")
;
po::options_description config("Configuration file options");
config.add_options()
("verbose",po::value<bool>(&verbose)->default_value(true), "Determines whether additional information should be emitted")
("mnParFix",po::value< vector<string> >(&_mnParFixs), "minuit parameters can be fixed here")
("massRangeMin",po::value<double>(&_massMin), "min of eta pi pi mass range for mass indep. fit")
("massRangeMax",po::value<double>(&_massMax), "max of eta pi pi mass range for mass indep. fit")
;
po::options_description cmdline_options;
cmdline_options.add(desc).add(common);
po::options_description config_file_options;
config_file_options.add(config).add(common);
po::variables_map vm;
po::store(po::parse_command_line(argc, argv, cmdline_options), vm);
po::notify(vm);
// Check the name of the configuation file
if(_configFile.empty() || _configFile == "empty" || _configFile == "unknown")
{
std::cout << cmdline_options << endl;
stringstream strError;
strError << "Error: Invalid configuration file name given: \"" << _configFile << "\"";
throw runtime_error(strError.str());
}
std::ifstream ifs(_configFile.c_str());
if(!ifs.good())
{
stringstream strError;
strError << "Error accessing configuratiocommonn file " << _configFile;
std::cout << cmdline_options << endl;
throw runtime_error(strError.str());
}
store(po::parse_config_file(ifs, config_file_options), vm);
po::notify(vm);
// Emit a help message, if necessary
if (vm.count("help"))
{
std::cout << cmdline_options << endl;
exit(0);
}
if(strErrLogMode == "debug") _errLogMode = debug;
else if(strErrLogMode == "trace") _errLogMode = trace;
else if(strErrLogMode == "routine") _errLogMode = routine;
else if(strErrLogMode == "warning") _errLogMode = warning;
else if(strErrLogMode == "error") _errLogMode = error;
else if(strErrLogMode == "alert") _errLogMode = alert;
else
{
_errLogMode = debug;
Warning << "ErrorLogger not (properly) set -> Use mode 'DEBUG' " ; // << endmsg;
}
if(verbose){
std::cout << "\nRunning with the following options using " << _configFile << ":\n\n"
<< "Error log mode: " << _errLogMode <<"\n\n"
<< "data file: " << _dataFile <<"\n\n"
<< "mc file: " << _mcFile <<"\n\n"
<< "file with start parameters for fit or qa: " << _paramFile << "\n\n"
<< "startHypo: " << _startHypo << "\n\n"
<< "mode: " << _mode << "\n\n"
<< endl;
std::vector<std::string>::const_iterator it;
for (it=_enabledHyps.begin(); it!=_enabledHyps.end();++it){
std::cout << "hypothesis\t" << (*it) << "\t enabled\n";
}
std::cout << std::endl;
for (it=_mnParFixs.begin(); it!=_mnParFixs.end();++it){
std::cout << "minuit parameter\t" << (*it) << "\t fixed\n";
}
std::cout << std::endl;
}
}
catch( std::exception & e )
{
cerr << "Error parsing the command line:" << endl;
cerr << e.what() << std::endl;
cerr << "You can use -h or --help to obtain the description of the program parameters." << endl;
cerr << "This is the command line options\n" << endl;
return false;
}
catch(...){
std::cerr << "Error parsing the command line. Use -h or --help to see the description of the program paramters." << endl;
return false;
}
return true;
}
#ifndef JpsiGamEtaPiPiParser_HH
#define JpsiGamEtaPiPiParser_HH
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <utility>
// Boost headers go here
#include <boost/version.hpp>
#if BOOST_VERSION < 103600
#error "Error: Boost should at least have version 1.36 !"
#endif /* BOOST_VERSION */
#include <boost/program_options.hpp>
#include <boost/filesystem.hpp>
// For Microsoft-compatible compilers
#if defined(_MSC_VER) && (_MSC_VER >= 1020)
#pragma once
#endif
// GenEvA headers go here
//#include <common/GCommonEnums.hpp>
//#include <common/GSerializationHelperFunctionsT.hpp>
//#include <geneva/GOptimizationEnums.hpp>
namespace po = boost::program_options;
//using namespace Gem::Geneva;
class JpsiGamEtaPiPiParser
{
public:
typedef enum tagerrLogMode {debug,trace,routine,warning,error,alert} enErrLogMode;
public:
JpsiGamEtaPiPiParser(int argc,char **argv)
: _configFile("./JpsiGamEtaPiPi.cfg")
, _errLogMode(debug)
, _dataFile("")
, _mcFile("")
, _paramFile("")
, _startHypo("base")
, _mode("plotmode")
, _massIndependentFit(false)
, _useCommonProductionPhases(false)
, _massMin(0.7)
, _massMax(3.1)
{
if (!parseCommandLine(argc, argv)) throw false;
}
const std::string& getConfigFile() const { return _configFile;}
const enErrLogMode& getErrLogMode() const { return _errLogMode; }
const std::string dataFile() const {return _dataFile;}
const std::string mcFile() const {return _mcFile;}
const std::string fitParamFile() const {return _paramFile;}
const std::vector<std::string>& enabledHyps() const { return _enabledHyps; }
const std::string startHypo() const {return _startHypo;}
const std::string mode() const {return _mode;}
const std::vector<std::string>& fixedParams() const { return _mnParFixs; }
const bool massIndependentFit() const {return _massIndependentFit; }
const bool useCommonProductionPhases() const {return _useCommonProductionPhases; }
const std::pair<double, double> massRange() const { return std::make_pair( _massMin, _massMax ) ; }
protected:
bool parseCommandLine(int argc,char **argv);
protected:
std::string _configFile;
enErrLogMode _errLogMode;
std::string _dataFile;
std::string _mcFile;
std::string _paramFile;
std::string _startHypo;
std::string _mode;
std::vector<std::string> _enabledHyps;
std::vector<std::string> _mnParFixs;
bool _massIndependentFit;
bool _useCommonProductionPhases;
double _massMin;
double _massMax;
};
#endif /* JpsiGamEtaPiPiParser_HH */
#include <iostream>
#include <cstring>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <boost/shared_ptr.hpp>
#include "Examples/JpsiGamEtaPiPi/JpsiGamEtaPiPiParser.hh"
#include "ErrLogger/ErrLogger.hh"
int main(int __argc,char *__argv[]){
// Parse the command line
static JpsiGamEtaPiPiParser theAppParams(__argc, __argv);
// Inform the audience about the execution mode
// emitExecutionMode(theAppParams.getAppExecMode());
return 0;
}
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