diff --git a/Examples/benchmark/coupledChannel/benchmarkCoupledChannel_runscriptGraNum b/Examples/benchmark/coupledChannel/benchmarkCoupledChannel_runscriptGraNum index 78873c62be1effa7907976a9e3bda165e5b4ffbd..c13d9de437fb64fd9dd9d344194ebd023bbe7756 100755 --- a/Examples/benchmark/coupledChannel/benchmarkCoupledChannel_runscriptGraNum +++ b/Examples/benchmark/coupledChannel/benchmarkCoupledChannel_runscriptGraNum @@ -1,7 +1,7 @@ #!/bin/sh echo Starting server for coupled channel pwa ... -$TOP_DIR/bin/coupledChannelApp -c global.cfg --pbarpFiles pbarpToEtaEtaPiReaction.cfg --pbarpFiles pbarpToKKPiReaction.cfg --epemFiles psiToGamKKpi.cfg --resFiles D0ToKmpippi0.cfg --mode serverGradientNum !>& log.server & +$TOP_DIR/bin/coupledChannelApp -c global.cfg --pbarpFiles pbarpToEtaEtaPiReaction.cfg --pbarpFiles pbarpToKKPiReaction.cfg --epemFiles psiToGamKKpi.cfg --resFiles D0ToKmpippi0.cfg --mode serverGradientNum !>& logGraNum.server & sleep 15 diff --git a/Examples/benchmark/coupledChannel/pbarpToEtaEtaPiReaction.cfg b/Examples/benchmark/coupledChannel/pbarpToEtaEtaPiReaction.cfg index 7398989b25f717d6592c6bbdbc20c77ebbdff9d8..48fdc9a5ea2b5e1d8e26eca3ba5a45b431c3330b 100644 --- a/Examples/benchmark/coupledChannel/pbarpToEtaEtaPiReaction.cfg +++ b/Examples/benchmark/coupledChannel/pbarpToEtaEtaPiReaction.cfg @@ -15,7 +15,7 @@ ratioMcToData = 3 useDataEventWeight = false useMCEventWeight = false -paramFile = ./startparamsV4.dat +paramFile = ./startparams.dat mode = pwa verbose = 1 diff --git a/MinFunctions/PwaFcnServer.cc b/MinFunctions/PwaFcnServer.cc index 422301c09877bb91bafa7dc56406b994a2a648a4..35805978bd0465f4537ee0c633aa6b8113931193 100644 --- a/MinFunctions/PwaFcnServer.cc +++ b/MinFunctions/PwaFcnServer.cc @@ -34,6 +34,7 @@ #include "PwaUtils/AbsLh.hh" #include "PwaUtils/DataUtils.hh" #include "PwaUtils/NetworkServer.hh" +#include "Utils/PawianConstants.hh" #include "ConfigParser/ParserBase.hh" #include "ErrLogger/ErrLogger.hh" @@ -43,6 +44,7 @@ template<typename T> PwaFcnServer<T>::PwaFcnServer(std::shared_ptr<NetworkServer> netServer) : AbsFcn<T>() , _networkServerPtr(netServer) + , _numStepSize(std::sqrt(std::numeric_limits<double>::epsilon()*150.)) { this->_defaultPawianParms = GlobalEnv::instance()->defaultPawianParams(); this->_currentPawianParms = GlobalEnv::instance()->startPawianParams(); @@ -60,41 +62,6 @@ double PwaFcnServer<T>::operator()(const std::vector<double>& par) const this->_currentPawianParms->SetAllValues(par); ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); - // std::map<ChannelID, LHData> theLHDataMap; - // std::ostringstream output; - // std::ostringstream outputLHDump; - // _networkServerPtr->BroadcastParams(this->_currentPawianParms->Params()); - // if(!_networkServerPtr->WaitForLH(theLHDataMap)) - // result = 0; - // else{ - // if(theLHDataMap.size() > 1){ - // outputLHDump << result << "\t"; - // } - // // Add LLHs of different channels - // output << "current LH = "; - // for(auto it = theLHDataMap.begin(); it!=theLHDataMap.end();++it){ - // (*it).second.weightSum = _networkServerPtr->weightSum((*it).first); - // (*it).second.squaredWeightSum = _networkServerPtr->squaredWeightSum((*it).first); - // (*it).second.num_mc = _networkServerPtr->numMCs((*it).first); - // double channelLH = AbsLh::mergeLogLhData((*it).second, (*it).first); - // result += channelLH; - // output << std::setprecision(16) << channelLH << "\t"; - // outputLHDump << std::setprecision(16) << channelLH << "\t"; - // } - // if(theLHDataMap.size() > 1){ - // output << "sum = " << result; - // } - // } - - // if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeLhPrint() == 0){ - // InfoMsg << output.str() << endmsg; - // } - // if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeTimer() == 0) this->printTimer(); - // if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeParamsPrint() == 0) this->printFitParams(this->_currentPawianParms); - // if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeParamsDump() == 0) this->dumpFitParams(this->_currentPawianParms); - // if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeLhDump() == 0) this->dumpLhVals(outputLHDump.str()); - // this->_fcnCounter++; - // return result; return collectLH(); } @@ -104,36 +71,42 @@ double PwaFcnServer<T>::collectLH() const{ std::map<ChannelID, LHData> theLHDataMap; std::ostringstream output; std::ostringstream outputLHDump; + bool lhPrint=false; + if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeLhPrint() == 0) lhPrint=true; + bool lhDump=false; + if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeLhDump() == 0) lhDump=true; + _networkServerPtr->BroadcastParams(this->_currentPawianParms->Params()); if(!_networkServerPtr->WaitForLH(theLHDataMap)) result = 0; else{ - if(theLHDataMap.size() > 1){ + if(lhDump && theLHDataMap.size() > 1){ outputLHDump << result << "\t"; } // Add LLHs of different channels - output << "current LH = "; + if(lhPrint) output << "current LH = "; for(auto it = theLHDataMap.begin(); it!=theLHDataMap.end();++it){ (*it).second.weightSum = _networkServerPtr->weightSum((*it).first); (*it).second.squaredWeightSum = _networkServerPtr->squaredWeightSum((*it).first); (*it).second.num_mc = _networkServerPtr->numMCs((*it).first); double channelLH = AbsLh::mergeLogLhData((*it).second, (*it).first); result += channelLH; - output << std::setprecision(16) << channelLH << "\t"; - outputLHDump << std::setprecision(16) << channelLH << "\t"; + if(lhPrint) output << std::setprecision(16) << channelLH << "\t"; + if(lhDump) outputLHDump << std::setprecision(16) << channelLH << "\t"; } - if(theLHDataMap.size() > 1){ + if(lhPrint && theLHDataMap.size() > 1){ output << "sum = " << result; } } - if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeLhPrint() == 0){ + if(lhPrint){ InfoMsg << output.str() << endmsg; } + if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeTimer() == 0) this->printTimer(); if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeParamsPrint() == 0) this->printFitParams(this->_currentPawianParms); if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeParamsDump() == 0) this->dumpFitParams(this->_currentPawianParms); - if(this->_fcnCounter%GlobalEnv::instance()->parser()->stepSizeLhDump() == 0) this->dumpLhVals(outputLHDump.str()); + if(lhDump) this->dumpLhVals(outputLHDump.str()); this->_fcnCounter++; return result; } @@ -145,41 +118,32 @@ std::vector<double> PwaFcnServer<T>::Gradient(const std::vector<double>& par) co this->_currentPawianParms->SetAllValues(par); ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); double LHBase=collectLH(); - double epsilon=1.e-6; + double epsilon; + for(unsigned int i=0; i<par.size(); ++i){ if(this->_currentPawianParms->IsFixed(i)) resultVec.at(i)=0.; else{ double currentVal=this->_currentPawianParms->Value(i); + epsilon=_numStepSize*std::abs(currentVal); + if ((this->_currentPawianParms->GetName(i)).substr( (this->_currentPawianParms->GetName(i)).length() - 3 ) == "Phi"){ + epsilon=_numStepSize*PawianConstants::pi; + } + else if (std::abs(currentVal)<1.e-10) epsilon=_numStepSize*1.e-10; + double dx=(currentVal+epsilon)-currentVal; if(this->_currentPawianParms->HasLimits(i) && std::abs(currentVal-this->_currentPawianParms->UpperLimit(i))<1.e-6){ this->_currentPawianParms->SetValue(i, currentVal-epsilon); ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); double currentLH=collectLH(); - resultVec.at(i)=(LHBase-currentLH)/epsilon; - // InfoMsg << "resultVecLow.at(" << i << ")= " << resultVec.at(i) << endmsg; + resultVec.at(i)=(LHBase-currentLH)/dx; + // InfoMsg << "resultVecLow.at(" << i << ")= " << resultVec.at(i) << endmsg; } else{ this->_currentPawianParms->SetValue(i, currentVal+epsilon); ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); double currentLH=collectLH(); - resultVec.at(i)=(currentLH-LHBase)/epsilon; - // InfoMsg << "resultVecHigh.at(" << i << ")= " << resultVec.at(i) << endmsg; + resultVec.at(i)=(currentLH-LHBase)/dx; + // InfoMsg << "resultVecHigh.at(" << i << ")= " << resultVec.at(i) << endmsg; } - // else{ - // this->_currentPawianParms->SetValue(i, currentVal-1.e-6); - // ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); - // double currentLH=collectLH(); - // resultVec.at(i)=(LHBase-currentLH)/1.e-6; - // InfoMsg << "resultVecLow.at(" << i << ")= " << resultVec.at(i) << endmsg; - // } - // this->_currentPawianParms->SetValue(i, currentVal+1.e-6); - // ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); - // double currentLHh=collectLH(); - // this->_currentPawianParms->SetValue(i, currentVal-1.e-6); - // ParamDepHandler::instance()->ApplyDependencies(this->_currentPawianParms); - // double currentLHl=collectLH(); - // resultVec.at(i)=(currentLHh-currentLHl)/(2.*1.e-6); - //InfoMsg << "resultVec.at(" << i << ")= " << resultVec.at(i) << endmsg; - this->_currentPawianParms->SetValue(i, currentVal); } } diff --git a/MinFunctions/PwaFcnServer.hh b/MinFunctions/PwaFcnServer.hh index c13f033d19268473bffda240591d08bfc2fea4e6..d9435c455f57b876e5a9f72ae87db638a947ae4b 100644 --- a/MinFunctions/PwaFcnServer.hh +++ b/MinFunctions/PwaFcnServer.hh @@ -49,6 +49,7 @@ namespace ROOT { protected: virtual double collectLH() const; std::shared_ptr<NetworkServer> _networkServerPtr; + const double _numStepSize; }; } // namespace Minuit2 } // namespace ROOT