Something went wrong on our end
-
Julian Pychy authoredde5a8cc7
pbarpReactionApp.cc 15.68 KiB
#include <iostream>
#include <string>
#include <cstdlib>
#include <vector>
#include <map>
#include <iterator>
#include <boost/shared_ptr.hpp>
#include "TROOT.h"
#include "pbarpUtils/pbarpParser.hh"
#include "Particle/ParticleTable.hh"
#include "Particle/Particle.hh"
#include "Particle/PdtParser.hh"
#include "ErrLogger/ErrLogger.hh"
#include "pbarpUtils/pbarpStatesLS.hh"
#include "PwaUtils/AbsLh.hh"
#include "PwaUtils/FitParamsBase.hh"
#include "PwaUtils/StreamFitParmsBase.hh"
#include "PwaUtils/PwaFcnBase.hh"
#include "PwaUtils/PwaCovMatrix.hh"
#include "Utils/PawianCollectionUtils.hh"
#include "Utils/ErrLogUtils.hh"
#include "pbarpUtils/pbarpEnv.hh"
#include "pbarpUtils/pbarpReaction.hh"
#include "pbarpUtils/pbarpBaseLh.hh"
#include "pbarpUtils/pbarpHeliLh.hh"
#include "pbarpUtils/pbarpCanoLh.hh"
#include "pbarpUtils/pbarpEvtReader.hh"
#include "PwaUtils/EvtDataBaseList.hh"
//#include "pbarpUtils/pbarpEventList.hh"
#include "pbarpUtils/pbarpHist.hh"
#include "pbarpUtils/spinDensityHist.hh"
#include "Event/Event.hh"
#include "Event/EventList.hh"
#include "PwaUtils/NetworkServer.hh"
#include "PwaUtils/NetworkClient.hh"
#include "PwaUtils/WelcomeScreen.hh"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnStrategy.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnScan.h"
int main(int __argc,char *__argv[]){
clock_t start, end;
start= clock();
Info << welcomeScreen << endmsg;
// Disable output buffering
setvbuf(stdout, NULL, _IONBF, 0);
// Parse the command line
pbarpParser* theAppParams=new pbarpParser(__argc, __argv);
// Set the desired error logging mode
setErrLogMode(theAppParams->getErrLogMode());
pbarpEnv::instance()->setup(theAppParams);
boost::shared_ptr<pbarpReaction> thepbarpReaction=pbarpEnv::instance()->reaction();
thepbarpReaction->print(std::cout);
// boost::shared_ptr<pbarpDataBaseList> thepbarbDataBaseListPtr(new pbarpDataBaseList());
// boost::shared_ptr<AbsLh> theLhPtr(new pbarpBaseLh(thepbarbDataBaseListPtr));
std::string mode=theAppParams->mode();
boost::shared_ptr<FitParamsBase> theFitParamBase=boost::shared_ptr<FitParamsBase>(new FitParamsBase());
const std::string datFile=theAppParams->dataFile();
const std::string mcFile=theAppParams->mcFile();
Info << "data file: " << datFile ; // << endmsg;
Info << "mc file: " << mcFile ; // << endmsg;
std::vector<std::string> dataFileNames;
dataFileNames.push_back(datFile);
std::vector<std::string> mcFileNames;
mcFileNames.push_back(mcFile);
bool withEvtWeight=theAppParams->useEvtWeight();
Info << "EvtWeight: " << withEvtWeight << endmsg;
int noFinalStateParticles=pbarpEnv::instance()->noFinalStateParticles();
pbarpEvtReader eventReaderData(dataFileNames, noFinalStateParticles, 0, withEvtWeight);
EventList eventsData;
eventReaderData.fillAll(eventsData);
Info << "\nFile has " << eventsData.size() << " events. Each event has "
<< eventsData.nextEvent()->size() << " final state particles.\n" ; // << endmsg;
eventsData.rewind();
Event* anEvent;
int evtCount = 0;
while ((anEvent = eventsData.nextEvent()) != 0 && evtCount < 10) {
Info << "\n";
for(int i=0; i<noFinalStateParticles; ++i){
Info << (*anEvent->p4(i)) << "\tm = " << anEvent->p4(i)->Mass() << "\n";
}
Info << "\n" << endmsg;
; // << endmsg;
++evtCount;
}
eventsData.rewind();
pbarpEvtReader eventReaderMc(mcFileNames, noFinalStateParticles, 0, false);
EventList mcData;
eventReaderMc.fillAll(mcData);
Info << "\nFile has " << mcData.size() << " events. Each event has "
<< mcData.nextEvent()->size() << " final state particles.\n" ; // << endmsg;
mcData.rewind();
evtCount = 0;
while ((anEvent = mcData.nextEvent()) != 0 && evtCount < 10) {
Info << "\n";
for(int i=0; i<noFinalStateParticles; ++i){
Info << (*anEvent->p4(i)) << "\tm = " << anEvent->p4(i)->Mass() << "\n";
}
Info << "\n" << endmsg;
; // << endmsg;
++evtCount;
}
mcData.rewind();
// boost::shared_ptr<pbarpEventList> pbarpEventListPtr(new pbarpEventList()); boost::shared_ptr<EvtDataBaseList> pbarpEventListPtr(new EvtDataBaseList(pbarpEnv::instance()));
pbarpEventListPtr->ratioMcToData(theAppParams->ratioMcToData());
pbarpEventListPtr->read(eventsData, mcData);
// pbarpEventListPtr->read4Vecs();
std::string prodFormalism=theAppParams->productionFormalism();
boost::shared_ptr<AbsLh> theLhPtr;
if(prodFormalism=="Cano") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpCanoLh(pbarpEventListPtr));
else if(prodFormalism=="Heli") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpHeliLh(pbarpEventListPtr));
else {
Alert << "prodFormalism\t" << prodFormalism << "\tdoesn't exist!!!" << endmsg;
exit(1);
}
if (mode=="dumpDefaultParams"){
fitParams defaultVal;
fitParams defaultErr;
theLhPtr->getDefaultParams(defaultVal, defaultErr);
std::stringstream defaultparamsname;
defaultparamsname << "defaultparams" << pbarpEnv::instance()->outputFileNameSuffix() << ".dat";
std::ofstream theStreamDefault ( defaultparamsname.str().c_str() );
theFitParamBase->dumpParams(theStreamDefault, defaultVal, defaultErr);
return 0;
}
std::string paramStreamerPath=theAppParams->fitParamFile();
std::string outputFileNameSuffix= pbarpEnv::instance()->outputFileNameSuffix();
StreamFitParmsBase theParamStreamer(paramStreamerPath, theLhPtr);
fitParams theStartparams=theParamStreamer.getFitParamVal();
fitParams theErrorparams=theParamStreamer.getFitParamErr();
PwaFcnBase theFcn(theLhPtr, theFitParamBase, outputFileNameSuffix);
MnUserParameters upar;
theFitParamBase->setMnUsrParams(upar, theStartparams, theErrorparams);
std::cout << "\n\n**************** Minuit Fit parameter **************************" << std::endl;
for (int i=0; i<int(upar.Params().size()); ++i){
std::cout << upar.Name(i) << "\t" << upar.Value(i) << "\t" << upar.Error(i) << std::endl;
}
const std::vector<std::string> fixedParams=theAppParams->fixedParams();
const unsigned int noOfFreeFitParams = upar.Params().size()-fixedParams.size();
if (mode=="qaMode"){
Info << "\nThe parameter values are: " << "\n" << endmsg;
theFitParamBase->printParams(theStartparams);
Info << "\nThe parameter errors are: " << "\n" << endmsg;
theFitParamBase->printParams(theErrorparams);
double theLh=theLhPtr->calcLogLh(theStartparams);
Info <<"theLh = "<< theLh << endmsg;
pbarpHist theHist(theLhPtr, theStartparams);
double evtWeightSumData = pbarpEventListPtr->NoOfWeightedDataEvts();
double BICcriterion=2.*theLh+noOfFreeFitParams*log(evtWeightSumData);
double AICcriterion=2.*theLh+2.*noOfFreeFitParams;
double AICccriterion=AICcriterion+2.*noOfFreeFitParams*(noOfFreeFitParams+1)/(evtWeightSumData-noOfFreeFitParams-1);
Info << "noOfFreeFitParams:\t" <<noOfFreeFitParams;
Info << "evtWeightSumData:\t" <<evtWeightSumData;
Info << "BIC:\t" << BICcriterion << endmsg;
Info << "AIC:\t" << AICcriterion << endmsg;
Info << "AICc:\t" << AICccriterion << endmsg;
std::ostringstream qaSummaryFileName; qaSummaryFileName << "qaSummary" << outputFileNameSuffix << ".dat";
std::ofstream theQaStream ( qaSummaryFileName.str().c_str() );
theQaStream << "BIC\t" << BICcriterion << "\n";
theQaStream << "AICa\t" << AICcriterion << "\n";
theQaStream << "AICc\t" << AICccriterion << "\n";
theQaStream << "logLh\t" << theLh << "\n";
theQaStream << "free parameter\t" << noOfFreeFitParams << "\n";
theQaStream.close();
end= clock();
double cpuTime= (end-start)/ (CLOCKS_PER_SEC);
Info << "cpuTime:\t" << cpuTime << "\tsec" << endmsg;
return 1;
}
if (mode=="pwa"){
bool cacheAmps = theAppParams->cacheAmps();
Info << "caching amplitudes enabled / disabled:\t" << cacheAmps << endmsg;
if (cacheAmps) theLhPtr->cacheAmplitudes();
std::vector<std::string>::const_iterator itFix;
for (itFix=fixedParams.begin(); itFix!=fixedParams.end(); ++itFix){
upar.Fix( (*itFix) );
}
MnMigrad migrad(theFcn, upar);
Info <<"start migrad "<< endmsg;
FunctionMinimum min = migrad();
if(!min.IsValid()) {
//try with higher strategy
Info <<"FM is invalid, try with strategy = 2."<< endmsg;
MnMigrad migrad2(theFcn, min.UserState(), MnStrategy(2));
min = migrad2();
}
MnUserParameters finalUsrParameters=min.UserParameters();
const std::vector<double> finalParamVec=finalUsrParameters.Params();
fitParams finalFitParams=theStartparams;
theFitParamBase->getFitParamVal(finalParamVec, finalFitParams);
theFitParamBase->printParams(finalFitParams);
double theLh=theLhPtr->calcLogLh(finalFitParams);
Info <<"theLh = "<< theLh << endmsg;
const std::vector<double> finalParamErrorVec=finalUsrParameters.Errors();
fitParams finalFitErrs=theErrorparams;
theFitParamBase->getFitParamVal(finalParamErrorVec, finalFitErrs);
std::ostringstream finalResultname;
finalResultname << "finalResult" << outputFileNameSuffix << ".dat";
std::ofstream theStream ( finalResultname.str().c_str() );
theFitParamBase->dumpParams(theStream, finalFitParams, finalFitErrs);
std::cout << min << std::endl;
std::cout << "\n\n**************** Minuit FunctionMinimum information ******************" << std::endl;
if(min.IsValid()) std::cout << "Function minimum is valid." << std::endl;
else std::cout << "*** Function minimum is invalid! ***" << std::endl;
if(!min.HasValidParameters()) std::cout << "hasValidParameters() returned FALSE" << std::endl;
if(!min.HasValidCovariance()) std::cout << "hasValidCovariance() returned FALSE" << std::endl;
if(!min.HasAccurateCovar()) std::cout << "hasAccurateCovar() returned FALSE" << std::endl;
if(!min.HasPosDefCovar()) std::cout << "hasPosDefCovar() returned FALSE" << std::endl;
if(!min.HasMadePosDefCovar()) std::cout << "hasMadePosDefCovar() returned FALSE" << std::endl;
if(!min.HasCovariance()) std::cout << "hasCovariance() returned FALSE" << std::endl;
if(min.HasReachedCallLimit()) std::cout << "hasReachedCallLimit() returned TRUE" << std::endl; if(min.IsAboveMaxEdm()) std::cout << "isAboveMaxEdm() returned TRUE" << std::endl;
if(min.HesseFailed()) std::cout << "hesseFailed() returned TRUE" << std::endl;
MnUserCovariance theCovMatrix = min.UserCovariance();
std::ostringstream serializationFileName;
serializationFileName << "serializedOutput" << outputFileNameSuffix << ".dat";
std::ofstream serializationStream(serializationFileName.str().c_str());
boost::archive::text_oarchive boostOutputArchive(serializationStream);
if(min.IsValid()){
PwaCovMatrix thePwaCovMatrix(theCovMatrix, finalUsrParameters, finalFitParams);
boostOutputArchive << thePwaCovMatrix;
}
return 0;
}
if(mode == "spinDensity"){
std::string serializationFileName = pbarpEnv::instance()->serializationFileName();
std::ifstream serializationStream(serializationFileName.c_str());
if(!serializationStream.is_open()){
Alert << "Could not open serialization file." << endmsg;
return 0;
}
boost::archive::text_iarchive boostInputArchive(serializationStream);
PwaCovMatrix thePwaCovMatrix;
boostInputArchive >> thePwaCovMatrix;
spinDensityHist theSpinDensityHists(theLhPtr, theStartparams, thePwaCovMatrix);
}
if(mode == "server"){
std::vector<std::string>::const_iterator itFix;
for (itFix=fixedParams.begin(); itFix!=fixedParams.end(); ++itFix){
upar.Fix( (*itFix) );
}
boost::shared_ptr<NetworkServer> theServer(new NetworkServer(theAppParams->serverPort(),
theAppParams->noOfClients(),
pbarpEventListPtr->getDataVecs().size(),
pbarpEventListPtr->getMcVecs().size()));
theFcn.SetServerMode(theServer);
theServer->WaitForFirstClientLogin();
MnMigrad migrad(theFcn, upar);
Info <<"start migrad "<< endmsg;
FunctionMinimum min = migrad();
if(!min.IsValid()) {
//try with higher strategy
Info <<"FM is invalid, try with strategy = 2."<< endmsg;
MnMigrad migrad2(theFcn, min.UserState(), MnStrategy(2));
min = migrad2();
}
theServer->SendClosingMessage();
Info << "Closing server." << endmsg;
MnUserParameters finalUsrParameters=min.UserParameters();
const std::vector<double> finalParamVec=finalUsrParameters.Params();
fitParams finalFitParams=theStartparams;
theFitParamBase->getFitParamVal(finalParamVec, finalFitParams);
theFitParamBase->printParams(finalFitParams);
const std::vector<double> finalParamErrorVec=finalUsrParameters.Errors();
fitParams finalFitErrs=theErrorparams;
theFitParamBase->getFitParamVal(finalParamErrorVec, finalFitErrs);
std::ostringstream finalResultname;
finalResultname << "finalResult" << outputFileNameSuffix << ".dat";
std::ofstream theStream ( finalResultname.str().c_str() );
theFitParamBase->dumpParams(theStream, finalFitParams, finalFitErrs);
std::cout << min << std::endl;
std::cout << "\n\n**************** Minuit FunctionMinimum information ******************" << std::endl;
if(min.IsValid()) std::cout << "Function minimum is valid." << std::endl;
else std::cout << "*** Function minimum is invalid! ***" << std::endl;
if(!min.HasValidParameters()) std::cout << "hasValidParameters() returned FALSE" << std::endl;
if(!min.HasValidCovariance()) std::cout << "hasValidCovariance() returned FALSE" << std::endl;
if(!min.HasAccurateCovar()) std::cout << "hasAccurateCovar() returned FALSE" << std::endl;
if(!min.HasPosDefCovar()) std::cout << "hasPosDefCovar() returned FALSE" << std::endl;
if(!min.HasMadePosDefCovar()) std::cout << "hasMadePosDefCovar() returned FALSE" << std::endl;
if(!min.HasCovariance()) std::cout << "hasCovariance() returned FALSE" << std::endl;
if(min.HasReachedCallLimit()) std::cout << "hasReachedCallLimit() returned TRUE" << std::endl;
if(min.IsAboveMaxEdm()) std::cout << "isAboveMaxEdm() returned TRUE" << std::endl;
if(min.HesseFailed()) std::cout << "hesseFailed() returned TRUE" << std::endl;
MnUserCovariance theCovMatrix = min.UserCovariance();
std::ostringstream serializationFileName;
serializationFileName << "serializedOutput" << outputFileNameSuffix << ".dat";
std::ofstream serializationStream(serializationFileName.str().c_str());
boost::archive::text_oarchive boostOutputArchive(serializationStream);
if(min.IsValid()){
PwaCovMatrix thePwaCovMatrix(theCovMatrix, finalUsrParameters, finalFitParams);
boostOutputArchive << thePwaCovMatrix;
}
}
if(mode == "client"){
bool cacheAmps = theAppParams->cacheAmps();
Info << "caching amplitudes enabled / disabled:\t" << cacheAmps << endmsg;
if (cacheAmps) theLhPtr->cacheAmplitudes();
std::ostringstream portStringStream;
portStringStream << theAppParams->serverPort();
NetworkClient theClient(theAppParams->serverAddress(), portStringStream.str());
if(!theClient.Login())
return 0;
while(true){
if(!theClient.WaitForParams())
return 0;
const std::vector<double> currentParamVec=theClient.GetParams();
fitParams currentFitParams=theStartparams;
theFitParamBase->getFitParamVal(currentParamVec, currentFitParams);
LHData theLHData;
theLhPtr->calcLogLhDataClient(currentFitParams, theLHData, theClient.GetEventLimits());
if(!theClient.SendLH(theLHData.logLH_data, theLHData.weightSum, theLHData.LH_mc))
return 0;
}
}
return 0;
}