//#include <getopt.h> //#include <fstream> //#include <string> #include <math.h> #include <stdio.h> #include "Minuit2/MnUserParameters.h" #include "Examples/MATpbarpToOmegaPi/MOmegaPiFcn.hh" #include "Examples/MATpbarpToOmegaPi/OmegaPiLh.hh" #include "ErrLogger/ErrLogger.hh" using namespace ROOT::Minuit2; MOmegaPiFcn::MOmegaPiFcn(boost::shared_ptr<OmegaPiLh> omegaPiLh) : _omegaPiLhPtr(omegaPiLh), _barpToOmegaPi0States(omegaPiLh->omegaPi0States()) { if (0==_omegaPiLhPtr) { Alert << "OmegaPiLh pointer is 0 !!!!" << endmsg; exit(1); } } MOmegaPiFcn::~MOmegaPiFcn() { } double MOmegaPiFcn::operator()(const std::vector<double>& par) const { OmegaPiData::fitParamVal theFitParmValTmp; setFitParamVal(theFitParmValTmp, par); double result=_omegaPiLhPtr->calcLogLh(theFitParmValTmp); // print fit paramss std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaSinglet=_barpToOmegaPi0States->jpclsSinglet(); std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaTriplet0=_barpToOmegaPi0States->jpclsTriplet0(); std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaTriplet1=_barpToOmegaPi0States->jpclsTriplet1(); std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itJPCLS; DebugMsg << "logLh= " << result <<endmsg; DebugMsg << "***fit parameter singlet states*** " <<endmsg; for ( itJPCLS=JPCLSOmegaSinglet.begin(); itJPCLS!=JPCLSOmegaSinglet.end(); ++itJPCLS){ DebugMsg<< (*itJPCLS)->name()<< endmsg; std::pair<double, double> tmpParam=theFitParmValTmp.omegaProdSinglet[(*itJPCLS)]; DebugMsg <<"\t mag:" << tmpParam.first <<"\t phi:" << tmpParam.second << endmsg; } DebugMsg << "***fit parameter triplet m=0 states*** " <<endmsg; for ( itJPCLS=JPCLSOmegaTriplet0.begin(); itJPCLS!=JPCLSOmegaTriplet0.end(); ++itJPCLS){ DebugMsg<< (*itJPCLS)->name()<< endmsg; std::pair<double, double> tmpParam=theFitParmValTmp.omegaProdTriplet0[(*itJPCLS)]; DebugMsg <<"\t mag:" << tmpParam.first <<"\t phi:" << tmpParam.second << endmsg; } DebugMsg << "***fit parameter triplet m=1 states*** " <<endmsg; for ( itJPCLS=JPCLSOmegaTriplet1.begin(); itJPCLS!=JPCLSOmegaTriplet1.end(); ++itJPCLS){ DebugMsg<< (*itJPCLS)->name()<< endmsg; std::pair<double, double> tmpParam=theFitParmValTmp.omegaProdTriplet1[(*itJPCLS)]; DebugMsg <<"\t mag:" << tmpParam.first <<"\t phi:" << tmpParam.second << endmsg; } DebugMsg << endmsg; return result; } double MOmegaPiFcn::Up() const { return .5; } void MOmegaPiFcn::setMnUsrParams(MnUserParameters& upar){ std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itJPCLS; std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaSinglet=_barpToOmegaPi0States->jpclsSinglet(); int counter=0; for ( itJPCLS=JPCLSOmegaSinglet.begin(); itJPCLS!=JPCLSOmegaSinglet.end(); ++itJPCLS){ //now fill the fitParameterMap std::string magStr=(*itJPCLS)->name()+"S"+"mag"; std::string phiStr=(*itJPCLS)->name()+"S"+"phi"; upar.Add(magStr, 0.5, .1, 0., 1.); if (counter>0) upar.Add(phiStr, 0., .1, -M_PI, M_PI); counter++; } std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaTriplet0=_barpToOmegaPi0States->jpclsTriplet0(); counter=0; for ( itJPCLS=JPCLSOmegaTriplet0.begin(); itJPCLS!=JPCLSOmegaTriplet0.end(); ++itJPCLS){ //now fill the fitParameterMap std::string magStr=(*itJPCLS)->name()+"T0"+"mag"; std::string phiStr=(*itJPCLS)->name()+"T0"+"phi"; upar.Add(magStr, 0.5, .1, 0., 1.); if (counter>0) upar.Add(phiStr, 0., .1, -M_PI, M_PI); counter++; } std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaTriplet1=_barpToOmegaPi0States->jpclsTriplet1(); counter=0; for ( itJPCLS=JPCLSOmegaTriplet1.begin(); itJPCLS!=JPCLSOmegaTriplet1.end(); ++itJPCLS){ //now fill the fitParameterMap std::string magStr=(*itJPCLS)->name()+"T1"+"mag"; std::string phiStr=(*itJPCLS)->name()+"T1"+"phi"; upar.Add(magStr, 0.5, .1, 0., 1.); if (counter>0) upar.Add(phiStr, 0., .1, -M_PI, M_PI); counter++; } } void MOmegaPiFcn::setFitParamVal(OmegaPiData::fitParamVal& theParamVal, const std::vector<double>& par) const{ std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itJPCLS; std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaSinglet=_barpToOmegaPi0States->jpclsSinglet(); std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaTriplet0=_barpToOmegaPi0States->jpclsTriplet0(); std::vector< boost::shared_ptr<const JPCLS> > JPCLSOmegaTriplet1=_barpToOmegaPi0States->jpclsTriplet1(); if (par.size()!= JPCLSOmegaSinglet.size()*2+JPCLSOmegaTriplet0.size()*2+JPCLSOmegaTriplet1.size()*2-3) { Alert << "size of parameters wrong!!! par.size()=" << par.size() << "\tJPCLSOmegaSinglet.size()+JPCLSOmegaTriplet0.size()+JPCLSOmegaTriplet1.size()-3=" << JPCLSOmegaSinglet.size()*2+JPCLSOmegaTriplet0.size()*2+JPCLSOmegaTriplet1.size()*2-3 << endmsg; exit(1); } int counter=0; for ( itJPCLS=JPCLSOmegaSinglet.begin(); itJPCLS!=JPCLSOmegaSinglet.end(); ++itJPCLS){ //now fill the fitParameterMap double mag=par[counter]; counter++; double phi=0.; if (counter>1){ phi=par[counter]; counter++; } std::pair <double,double> tmpParameter=make_pair(mag,phi); theParamVal.omegaProdSinglet[(*itJPCLS)]=tmpParameter; } for ( itJPCLS=JPCLSOmegaTriplet0.begin(); itJPCLS!=JPCLSOmegaTriplet0.end(); ++itJPCLS){ //now fill the fitParameterMap double mag=par[counter]; counter++; double phi=0.; if (counter>JPCLSOmegaSinglet.size()*2){ phi=par[counter]; counter++; } std::pair <double,double> tmpParameter=make_pair(mag,phi); theParamVal.omegaProdTriplet0[(*itJPCLS)]=tmpParameter; } for ( itJPCLS=JPCLSOmegaTriplet1.begin(); itJPCLS!=JPCLSOmegaTriplet1.end(); ++itJPCLS){ //now fill the fitParameterMap double mag=par[counter]; counter++; double phi=0.; if (counter>JPCLSOmegaSinglet.size()*2+JPCLSOmegaTriplet0.size()*2-1){ phi=par[counter]; counter++; } std::pair <double,double> tmpParameter=make_pair(mag,phi); theParamVal.omegaProdTriplet1[(*itJPCLS)]=tmpParameter; } }