diff --git a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc index cf105ee62a0b5e8d0ef313007548bea525b02aca..5c6f96be7db0de751ce24bd889aaa29143d82d2c 100644 --- a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc +++ b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.cc @@ -79,6 +79,8 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) po::options_description common("Common Options"); common.add_options() + ("datFile",po::value<string>(&_dataFile), "full path of data file") + ("mcFile",po::value<string>(&_mcFile), "full path of Monte Carlo file") ("lmax", po::value<unsigned>(&lMax)->default_value(lMax),"choose lmax.") ("pbarmom", po::value<unsigned>(&pbarMom)->default_value(pbarMom),"choose pbar momentum.") ("errLogMode,e", po::value<string>(&strErrLogMode)->default_value(strErrLogMode), @@ -88,7 +90,7 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) ("name,n", po::value<string>(&strName)->default_value("myApp"), "Name that is attached to all otuput file names to be able to run multiple fits in parallel.") ("LhMode", po::value<std::string>(&theLhMode)->default_value(theLhMode), - "Specifies the likelihood mode: OmegaPiLhGamma or OmegaTo3PiLhGamma") + "Specifies the likelihood mode: OmegaPiLhGamma, OmegaTo3PiLhGamma or OmegaTo3PiLhProd") ; po::options_description config("Configuration file options"); @@ -121,8 +123,6 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) "Specifies whether results should be returned even if they are not better than before") ("waitFactor", po::value<boost::uint32_t>(&waitFactor)->default_value(waitFactor), "Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individuum") - ("SourcePath", po::value<std::string>(&theSourcePath)->default_value(theSourcePath), - "Specifies the path to root directory of the source") ; po::options_description cmdline_options; @@ -334,23 +334,24 @@ bool ApplicationParameterLS::parseCommandLine(int argc, char **argv) if(verbose){ - std::cout << "\nRunning with the following options using " << configFile << ":\n" - << "nProducerThreads = " << (boost::uint16_t)nProducerThreads << "\n" // boost::uint8_t not printable on gcc ??? - << "populationSize = " << populationSize << "\n" - << "nParents = " << nParents << "\n" - << "maxIterations = " << maxIterations << "\n" - << "maxMinutes = " << maxMinutes << "\n" - << "reportIteration = " << reportIteration << "\n" - << "rScheme = " << (boost::uint16_t)rScheme << "\n" - << "sortingScheme = " << smode << "\n" - << "arraySize = " << arraySize << "\n" - << "processingCycles = " << processingCycles << "\n" - << "returnRegardless = " << (returnRegardless?"true":"false") << "\n" - << "lMax = " << lMax << "\n" - << "pbarMom = " << pbarMom << "\n" - << "errLogMode = " << strErrLogMode << "\n" - << "SourcePath = " << theSourcePath << "\n" - << "LhMode = " << theLhMode << "\n" + std::cout << "\nRunning with the following options using " << configFile << ":\n" + << "data file: " << _dataFile <<"\n" + << "mc file: " << _mcFile <<"\n" + << "nProducerThreads = " << (boost::uint16_t)nProducerThreads << "\n" // boost::uint8_t not printable on gcc ??? + << "populationSize = " << populationSize << "\n" + << "nParents = " << nParents << "\n" + << "maxIterations = " << maxIterations << "\n" + << "maxMinutes = " << maxMinutes << "\n" + << "reportIteration = " << reportIteration << "\n" + << "rScheme = " << (boost::uint16_t)rScheme << "\n" + << "sortingScheme = " << smode << "\n" + << "arraySize = " << arraySize << "\n" + << "processingCycles = " << processingCycles << "\n" + << "returnRegardless = " << (returnRegardless?"true":"false") << "\n" + << "lMax = " << lMax << "\n" + << "pbarMom = " << pbarMom << "\n" + << "errLogMode = " << strErrLogMode << "\n" + << "LhMode = " << theLhMode << "\n" << endl; } diff --git a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh index 26d642b4ecbc2e2e16faafbe6b338d92f57c341d..71659d23d6f713770256a54654ce3fe843a28a93 100644 --- a/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh +++ b/Examples/pbarpToOmegaPiLS/GArgumentParserLS.hh @@ -71,6 +71,8 @@ class ApplicationParameterLS public: ApplicationParameterLS(int argc,char **argv) : configFile("./GOmegaPiProject.cfg") + , _dataFile("/data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/PWA_0900.dat") + , _mcFile("/data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/mcPWA_0900.dat") , parallelizationMode(1) , ip("localhost") , port(10000) @@ -97,13 +99,14 @@ class ApplicationParameterLS , lMax(3) , pbarMom(600) , errLogMode(debug) - , theSourcePath("./") - , theLhMode("OmegaPiLhGamma") + , theLhMode("OmegaTo3PiLhGamma") { if (!parseCommandLine(argc, argv)) throw false; } const std::string& getConfigFile() const { return configFile;} + const std::string dataFile() const {return _dataFile;} + const std::string mcFile() const {return _mcFile;} const boost::uint16_t& getParallelizationMode() const {return parallelizationMode;} const bool& getServerMode() const { return serverMode; } const std::string& getIp() const { return ip; } @@ -126,7 +129,6 @@ class ApplicationParameterLS const enErrLogMode& getErrLogMode() const { return errLogMode; } const Gem::Common::serializationMode& getSerMode() const { return serMode; } const enExecMode& getAppExecMode() const { return enAppExecMode; } - const std::string& getSourcePath() const { return theSourcePath; } const std::string& getLhMode() const { return theLhMode; } const std::string& getPathStartParamFile() const { return strPathStartParamFile; } const std::string& getMinuitFixParamFile() const { return strMinuitFixParamFile; } @@ -146,7 +148,9 @@ protected: std::string strName; std::string strPathStartParamFile; std::string strMinuitFixParamFile; - std::string configFile; + std::string configFile; + std::string _dataFile; + std::string _mcFile; boost::uint16_t parallelizationMode; bool serverMode; std::string ip; @@ -176,7 +180,6 @@ protected: unsigned lMax; unsigned pbarMom; enErrLogMode errLogMode; - std::string theSourcePath; std::string theLhMode; }; diff --git a/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg b/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg index bc0f456aa6442c74c22a64c3d088c0017004e47e..c9b7e8f406a840de2b518b15c39914e5eb2380e6 100644 --- a/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg +++ b/Examples/pbarpToOmegaPiLS/GOmegaPiProjectLS.cfg @@ -3,7 +3,8 @@ # # # Authors: Ruediger Berlich # ################################################################################# - +datFile = /data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/PWA_0900.dat +mcFile = /data/puru2/julian/Pawian_0811/Pawian/Examples/pbarpToOmegaPiLS/data/mcPWA_0900.dat # The amount of random number producer threads nProducerThreads = 16 #nProducerThreads = 4 @@ -38,7 +39,8 @@ returnRegardless = 1 # Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individual waitFactor = 2 # Influences the maximum waiting time of the GBrokerEA after the arrival of the first evaluated individual -lmax = 2 -pbarmom = 600 +lmax = 4 +pbarmom = 900 errLogMode = debug -SourcePath = /data/sleipnir1/bertram/PawianGit1107/Pawian/ +appMode = 1 +LhMode = OmegaTo3PiLhProd \ No newline at end of file diff --git a/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc b/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc index 733849ae035b95ac51d19687587f94c81add8596..509bddaa577bcb09f387023f2ed021000852ea17 100644 --- a/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc +++ b/Examples/pbarpToOmegaPiLS/OmegaPiLSTestApp.cc @@ -18,6 +18,7 @@ #include "Examples/pbarpToOmegaPiLS/AbsOmegaPiLhLS.hh" #include "Examples/pbarpToOmegaPiLS/OmegaPiLhPi0GammaLS.hh" #include "Examples/pbarpToOmegaPiLS/OmegaTo3PiLhPi0GammaLS.hh" +#include "Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh" #include "Examples/pbarpToOmegaPiLS/MOmegaPiFcnLS.hh" #include "Examples/pbarpToOmegaPiLS/SpinDensityHistLS.hh" @@ -614,21 +615,20 @@ int main(int __argc,char *__argv[]){ exit(1); } - std::string piomegaDatFile; - std::string piomegaMcFile; + std::string piomegaDatFile=theAppParams.dataFile(); + std::string piomegaMcFile=theAppParams.mcFile(); + Info << "data file: " << piomegaDatFile << endmsg; + Info << "mc file: " << piomegaMcFile << endmsg; + 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); - } - else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){ - nParticlesPerEvt=4; - constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/om3pi_",theAppParams.getPbarMom(),piomegaDatFile); - constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mcom3pi_",theAppParams.getPbarMom(),piomegaMcFile); - } + nParticlesPerEvt=3; + } + else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhProd") || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){ + nParticlesPerEvt=4; + } else { Alert << "LhMode: " << theAppParams.getLhMode() << " does not exist!!!" << endmsg; @@ -636,7 +636,6 @@ int main(int __argc,char *__argv[]){ } - // constructPath(theAppParams.getSourcePath()+"/Examples/pbarpToOmegaPiLS/data/510_",theAppParams.getPbarMom(),piomegaDatFile); if (checkFileExist(piomegaDatFile)){ Info << "Using Data file " << piomegaDatFile << endmsg; } @@ -646,7 +645,6 @@ int main(int __argc,char *__argv[]){ exit(1); } - // constructPath(theAppParams.getSourcePath()+ "/Examples/pbarpToOmegaPiLS/data/mc510_",theAppParams.getPbarMom(),piomegaMcFile); if (checkFileExist(piomegaMcFile)){ Info << "Using Monte Carlo file " << piomegaMcFile << endmsg; } @@ -661,7 +659,9 @@ int main(int __argc,char *__argv[]){ ParticleTable pTable; PdtParser parser; - std::string pdtFile(theAppParams.getSourcePath()+"/Particle/pdt.table"); + + std::string theSourcePath=getenv("CMAKE_SOURCE_DIR"); + std::string pdtFile(theSourcePath+"/Particle/pdt.table"); if (!parser.parse(pdtFile, pTable)) { Alert << "Error: could not parse " << pdtFile << endmsg; exit(1); @@ -709,11 +709,10 @@ int main(int __argc,char *__argv[]){ 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"; } - else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma" || (theAppParams.getLhMode()=="OmegaTo3PiLhGammaBw")){ + else if ( (theAppParams.getLhMode()=="OmegaTo3PiLhGamma") || (theAppParams.getLhMode()=="OmegaTo3PiLhProd") || (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"; } - if (theAppParams.getAppExecMode()==ApplicationParameterLS::HistTest){ histTest(theOmegaPiEventPtr, theRootFilePath.str()); exit(1); @@ -737,6 +736,9 @@ int main(int __argc,char *__argv[]){ else if (theAppParams.getLhMode()=="OmegaTo3PiLhGamma"){ theOmegaPiLh = boost::shared_ptr<AbsOmegaPiLhLS>(new OmegaTo3PiLhPi0GammaLS(theOmegaPiEventPtr, pbarpToOmegaPi0StatesPtr)); } + else if (theAppParams.getLhMode()=="OmegaTo3PiLhProd"){ + theOmegaPiLh = boost::shared_ptr<AbsOmegaPiLhLS>(new OmegaTo3PiLhProdLS(theOmegaPiEventPtr, pbarpToOmegaPi0StatesPtr)); + } else{ Alert <<"LhMode " << theAppParams.getLhMode() << " doesn't exist !!! " << endmsg; exit(1); diff --git a/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.cc b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.cc new file mode 100644 index 0000000000000000000000000000000000000000..0a425aba9e092b8527bd9ecafe89a5741a67d870 --- /dev/null +++ b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.cc @@ -0,0 +1,84 @@ +#include <getopt.h> + +#include "Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh" +#include "Examples/pbarpToOmegaPiLS/AbsOmegaPiEventListLS.hh" +#include "PwaUtils/pbarpStatesLS.hh" +#include "Examples/pbarpToOmegaPiLS/pbarpToOmegaPi0StatesLS.hh" +#include "ErrLogger/ErrLogger.hh" + +// #include <geneva/GConstrainedDoubleObject.hpp> + + +OmegaTo3PiLhProdLS::OmegaTo3PiLhProdLS(boost::shared_ptr<const AbsOmegaPiEventListLS> theEvtList, boost::shared_ptr<const pbarpToOmegaPi0StatesLS> theStates) : + AbsOmegaPiLhLS(theEvtList, theStates, "OmegaTo3PiLhProd") +{ +} + +OmegaTo3PiLhProdLS::OmegaTo3PiLhProdLS(boost::shared_ptr<OmegaTo3PiLhProdLS> theOmegaTo3PiLhProdLSPtr): + AbsOmegaPiLhLS(theOmegaTo3PiLhProdLSPtr, "OmegaTo3PiLhProd") +{ +} + +OmegaTo3PiLhProdLS::~OmegaTo3PiLhProdLS() +{ +} + +double OmegaTo3PiLhProdLS::calcLogLh(const OmegaPiDataLS::fitParamVal& theParamVal){ + + double result=AbsOmegaPiLhLS::calcLogLh(theParamVal); + if (_globalItCounter%10 == 0) printFitParams(std::cout, theParamVal); + +// if (_globalItCounter%1000 == 0){ +// std::ofstream theStream ("currentResult.dat"); +// dumpCurrentResult(theStream, theParamVal); +// theStream.close(); +// } + + return result; +} + + +double OmegaTo3PiLhProdLS::calcEvtIntensity(OmegaPiDataLS::OmPiEvtDataLS* theData, const OmegaPiDataLS::fitParamVal& theParamVal){ + + double result=0.; + + for (Spin lamOmega=-1; lamOmega<=1; ++lamOmega){ + result+= norm(calcCoherentAmp(lamOmega,0, theParamVal, _singlet_JPCLS_States, theData)); //singlet amp + result+= norm(calcCoherentAmp(lamOmega,0, theParamVal, _triplet0_JPCLS_States, theData)); //triplet0 + result+= norm(calcCoherentAmp(lamOmega,1, theParamVal, _tripletp1_JPCLS_States, theData)); //triplet+1 + result+= norm(calcCoherentAmp(lamOmega,-1, theParamVal, _tripletm1_JPCLS_States, theData)); //triplet-1 + } + + return result; +} + + +complex<double> OmegaTo3PiLhProdLS::calcCoherentAmp(Spin lamOm, Spin Minit, const OmegaPiDataLS::fitParamVal& theParamVal, std::vector< boost::shared_ptr<const JPCLSls> >& theJPCLSlsStates, OmegaPiDataLS::OmPiEvtDataLS* theData){ + + complex<double> result(0.,0.); + + std::vector< boost::shared_ptr<const JPCLSls> >::const_iterator it; + for ( it=theJPCLSlsStates.begin(); it!=theJPCLSlsStates.end(); ++it){ + + + std::map< boost::shared_ptr<const JPCLSls>, pair<double, double>, pawian::Collection::SharedPtrLess >::const_iterator itmap; + itmap = theParamVal.lsParam.find((*it)); + pair<double, double> fitPair= (*itmap).second; + + if (fabs(lamOm)> (*it)->J ) continue; + + complex<double> omegaProdAmp=calcOmegaProdPartAmp(Minit, lamOm, (*it), fitPair, theData); + + result += omegaProdAmp; + + } + + + + return result; + } + + +void OmegaTo3PiLhProdLS::print(std::ostream& os) const{ + os << "OmegaTo3PiLhProdLS\n"; +} diff --git a/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh new file mode 100644 index 0000000000000000000000000000000000000000..f0f0aad8e81d01f5de67aaec056713969ca63df3 --- /dev/null +++ b/Examples/pbarpToOmegaPiLS/OmegaTo3PiLhProdLS.hh @@ -0,0 +1,63 @@ +#ifndef _OmegaTo3PiLhProdLS_H +#define _OmegaTo3PiLhProdLS_H + +#include <iostream> +#include <fstream> +#include <string> +#include <vector> +#include <complex> + +#include <cassert> +#include <boost/shared_ptr.hpp> + +#include "TROOT.h" +// #include <TSystem.h> +#include "qft++/topincludes/relativistic-quantum-mechanics.hh" +#include "Examples/pbarpToOmegaPiLS/OmegaPiDataLS.hh" +#include "Examples/pbarpToOmegaPiLS/AbsOmegaPiLhLS.hh" + + +#include "PwaUtils/DataUtils.hh" + +// using namespace std; + + +class AbsOmegaPiEventListLS; +class pbarpToOmegaPi0StatesLS; + +class OmegaTo3PiLhProdLS : public AbsOmegaPiLhLS{ + +public: + + // create/copy/destroy: + + ///Constructor + OmegaTo3PiLhProdLS(boost::shared_ptr<const AbsOmegaPiEventListLS>, boost::shared_ptr<const pbarpToOmegaPi0StatesLS>); + OmegaTo3PiLhProdLS(boost::shared_ptr<OmegaTo3PiLhProdLS>); + + /** Destructor */ + virtual ~OmegaTo3PiLhProdLS(); + + virtual AbsOmegaPiLhLS* clone_() const{ + return new OmegaTo3PiLhProdLS(_omegaPiEventListPtr, _omegaPi0StatesPtr); + } + + + // Getters: + virtual double calcLogLh(const OmegaPiDataLS::fitParamVal& theParamVal); + virtual double calcEvtIntensity(OmegaPiDataLS::OmPiEvtDataLS* theData, const OmegaPiDataLS::fitParamVal& theParamVal); + + + virtual void print(std::ostream& os) const; + +protected: + virtual complex<double> calcCoherentAmp(Spin lamOm, Spin Minit, const OmegaPiDataLS::fitParamVal& theParamVal, std::vector< boost::shared_ptr<const JPCLSls> >& theJPCLSlsStates, OmegaPiDataLS::OmPiEvtDataLS* theData); + + + +private: + + +}; + +#endif