IAP GITLAB

Commit b60221ad authored by Ralf Ulrich's avatar Ralf Ulrich

This is a very intrusive commit.

- rename a lot of classes
- split classes where appropriate
- create new CRMC main class
- no (almost) global variables any more
- boost program_options
- policy based particle output



git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@2888 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 7218d442
......@@ -54,8 +54,8 @@ FIND_PACKAGE (HepMC REQUIRED COMPONENTS HepMC HepMCfio)
## configure a header file to pass some of the CMake settings to the source code
CONFIGURE_FILE (
"${PROJECT_SOURCE_DIR}/src/crmcconfig.h.in"
"${PROJECT_BINARY_DIR}/src/crmcconfig.h"
"${PROJECT_SOURCE_DIR}/src/CRMCconfig.h.in"
"${PROJECT_BINARY_DIR}/src/CRMCconfig.h"
)
#add epos (must!)
......@@ -107,7 +107,7 @@ INCLUDE_DIRECTORIES ("${PROJECT_BINARY_DIR}/src") #cmake config file
INCLUDE_DIRECTORIES ("${HepMC_INCLUDE_DIRS}")
INCLUDE_DIRECTORIES ("${Boost_INCLUDE_DIRS}")
INCLUDE_DIRECTORIES ("${ROOT_INCLUDE_DIR}")
ADD_EXECUTABLE(crmc src/crmc.cc src/crmctrapfpe.c src/crmc-aaa.f src/models.F)
ADD_EXECUTABLE(crmc src/crmc.cc src/CRMC.cc src/CRMCinterface.cc src/CRMCoptions.cc src/OutputPolicyROOT.cc src/OutputPolicyHepMC.cc src/CRMCtrapfpe.c src/crmc-aaa.f src/models.F)
# linking of modules
TARGET_LINK_LIBRARIES (crmc ${HepMC_LIBRARIES})
......
#include <CRMC.h>
#include <CRMCinterface.h>
#include <CRMCoptions.h>
#include <OutputPolicyROOT.h>
#include <OutputPolicyHepMC.h>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <stdlib.h>
#include <string>
#include <stdio.h>
#warning remove these;
#include <TSystem.h>
#include <TStopwatch.h>
#include <CRMCconfig.h> //cmake generated
using namespace std;
template<class OutputPolicy>
CRMC<OutputPolicy>::CRMC(const CRMCoptions& cfg)
: fCfg(cfg) {
}
template<class OutputPolicy>
bool
CRMC<OutputPolicy>::init()
{
setbuf(stdout, 0); // set output to unbuffered
// open FORTRAN IO at first call
crmc_init_f_(fCfg.fSeed,
fCfg.fProjectileMomentum,
fCfg.fTargetMomentum,
fCfg.fProjectileId,
fCfg.fTargetId,
fCfg.fHEModel,
fCfg.fParamFileName.c_str());
OutputPolicy::InitOutput(fCfg);
return true;
}
template<class OutputPolicy>
bool
CRMC<OutputPolicy>::run()
{
TStopwatch watch;
std::cout.precision(10);
for (int iColl = 0; iColl < fCfg.fNCollision; ++iColl){
// cleanup vectors
gCRMC_data.Clean();
if (iColl % 1000 == 0 || (fCfg.fProjectileId+fCfg.fTargetId>400 && iColl %10))
cout << " ==[crmc]==> Collision number " << iColl+1 << endl;
// loop over collisions
crmc_f_(gCRMC_data.fNParticles,
gCRMC_data.fImpactParameter,
gCRMC_data.fPartId[0],
gCRMC_data.fPartPx[0],
gCRMC_data.fPartPy[0],
gCRMC_data.fPartPz[0],
gCRMC_data.fPartEnergy[0],
gCRMC_data.fPartMass[0]);
OutputPolicy::FillEvent(iColl);
}
std::cout.precision(2);
cout << " \n succesfully processed " << fCfg.fNCollision << " collision"
<< ( fCfg.fNCollision > 1 ? "s \n" : " \n" ) ;
watch.Stop();
const double realTime = watch.RealTime();
const double cpuTime = watch.CpuTime();
cout << " in " << realTime << " sec, "
<< ((double)realTime/fCfg.fNCollision) << " sec/collision, with "
<< (cpuTime/realTime*100) << "% cpu usage. \n" << endl;
return true;
}
template<class OutputPolicy>
bool
CRMC<OutputPolicy>::finish()
{
OutputPolicy::CloseOutput();
return true;
}
template class CRMC<OutputPolicyROOT>;
template class CRMC<OutputPolicyHepMC>;
#ifndef __CRMC_H
#define __CRMC_H
// //////////
// //////////
class CRMCoptions;
/*
class VCRMC {
VCRMC();
public:
static VCRMC* Create(const CRMCoptions& cfg);
bool init();
bool run();
bool finish();
void CleanVector();
};
*/
template<class OutputPolicy>
class CRMC : public OutputPolicy {
CRMC();
public:
CRMC(const CRMCoptions& cfg);// : OutputPolicy(cfg), fCfg(cfg) {}
bool init();
bool run();
bool finish();
private:
const CRMCoptions& fCfg;
};
#include <OutputPolicyROOT.h>
#include <OutputPolicyHepMC.h>
typedef CRMC<OutputPolicyROOT> CRMC2ROOT;
typedef CRMC<OutputPolicyHepMC> CRMC2HepMC;
#endif
#include <CRMCinterface.h>
CRMCdata gCRMC_data;
#ifndef __CRMC_INTERFACE_H
#define __CRMC_INTERFACE_H
// to exchange data between fortran and c++
struct CRMCdata {
CRMCdata() : fNParticles(0), fImpactParameter(0) {}
void Clean() { fNParticles = 0; }
// fortran output
const static unsigned int fMaxParticles = 9990; // HEP max number of particles
int fNParticles;
double fImpactParameter;
int fPartId[fMaxParticles];
double fPartPx[fMaxParticles];
double fPartPy[fMaxParticles];
double fPartPz[fMaxParticles];
double fPartEnergy[fMaxParticles];
double fPartMass[fMaxParticles];
};
extern CRMCdata gCRMC_data;
// -------------- forward declarations
extern "C"
......@@ -11,7 +33,7 @@ extern "C"
extern "C"
{
void crmc_init_f_( int&, double&, double&, int&, int&, int&, const char*);
void crmc_init_f_( const int&, const double&, const double&, const int&, const int&, const int&, const char*);
}
extern "C"
......
#include <CRMCoptions.h>
#include <CRMCconfig.h>
#include <iostream>
#include <iomanip>
#warning remove this;
#include <TSystem.h>
#include <boost/program_options.hpp>
namespace po = boost::program_options;
using namespace std;
CRMCoptions::CRMCoptions(int argc, char** argv)
: fError(false),
fNCollision(1),
fSeed(123),
fProjectileId(120),
fTargetId(12),
fHEModel(0),
fProjectileMomentum(300),
fTargetMomentum(0),
fParamFileName("crmc.param"),
fPrefix ( "crmc")
{
CheckEnvironment();
ParseOptions(argc, argv);
}
void
CRMCoptions::CheckEnvironment()
{
const char* conexRoot = gSystem->Getenv("CONEX_ROOT");
if ( conexRoot == 0 ) {
conexRoot = gSystem->Getenv("PWD");
}
}
void
CRMCoptions::ParseOptions(int argc, char** argv)
{
po::options_description desc("Options of CRMC");
ostringstream model_desc;
model_desc << "model [0=EPOS 1.99LHC, 1=EPOS 1.99"
#ifdef __QGSJET01__
<< ", 2=QGSJET01"
#endif
#ifdef __GHEISHA__
<< ", 3=Gheisha"
#endif
#ifdef __PYTHIA__
<< ", 4=Pythia 6.115"
#endif
#ifdef __HIJING__
<< ", 5=Hijing 1.38"
#endif
#ifdef __SIBYLL__
<< ", 6=Sibyll 2.1"
#endif
#ifdef __QGSJETII__
<< ", 7=QGSJETII-03"
#endif
#ifdef __Phojet__
<< ", 8=Phojet"
#endif
<< "]";
desc.add_options()
("help,h", "description of options")
("version,v", "show version and exis")
("output,o", po::value<string>(), "output mode: hepmc, hepmcgz, root, lhe, lhegz, ndst")
// ("verbosity,v", po::value<int>(), "verbosity (default: 0)")
("seed,s", po::value<int>(), "random seed (default: random)")
("number,n", po::value<int>(), "number of collisions")
("model,m", po::value<int>(), model_desc.str().c_str())
("projectile-momentum,p", po::value<double>(), "momentum/(GeV/c)")
("target-momentum,P", po::value<double>(), "momentum/(GeV/c)")
("projectile-id,i", po::value<int>(), "PDG or A*1000")
("target-id,I", po::value<int>(), "PDG or A*1000")
("config,c", po::value<string>(), "config file")
("prefix,x", po::value<string>(), "output prefix")
;
po::variables_map opt;
po::store(po::parse_command_line(argc, argv, desc), opt);
po::notify(opt);
if (opt.count("version")) {
cout << "crmc v" << CRMC_VERSION_MAJOR << "." << CRMC_VERSION_MINOR << endl;
exit(0);
}
// check if user needs help-printout
if (opt.count("help") || opt.count("output")==0) {
cout << endl << desc << endl;
exit(1);
}
const string om = opt["output"].as<string>();
if (om=="hepmc") fOutputMode = eHepMC;
else if (om=="hepmcgz") fOutputMode = eHepMCGZ;
else if (om=="lhe") fOutputMode = eLHE;
else if (om=="lhegz") fOutputMode = eLHEGZ;
else if (om=="root") fOutputMode = eROOT;
else if (om=="ndst") fOutputMode = eNDST;
else {
cout << " Wrong output type: " << om << endl;
cout << endl << desc << endl;
exit(1);
}
// parameter readout
if (opt.count("seed"))
fSeed = opt["seed"].as<int>();
if (opt.count("number"))
fNCollision = opt["number"].as<int>();
if (opt.count("model"))
fHEModel = opt["model"].as<int>();
if (opt.count("projectile-momentum"))
fProjectileMomentum = opt["projectile-momentum"].as<double>();
if (opt.count("target-momentum"))
fTargetMomentum = opt["target-momentum"].as<double>();
if (opt.count("projectile-id"))
fProjectileId = opt["projectile-id"].as<int>();
if (opt.count("target-id"))
fTargetId = opt["target-id"].as<int>();
if (opt.count("config"))
fParamFileName = opt["config"].as<string>();
if (opt.count("prefix"))
fPrefix = opt["prefix"].as<string>();
}
void
CRMCoptions::DumpConfig() const
{
cout << "\n >> crmc <<\n\n"
<< " prefix: " << fPrefix << "\n"
<< " beam'id: " << fProjectileId << "\n"
<< " beam's momentum: " << fProjectileMomentum << "\n"
<< " target's id: " << fTargetId << "\n"
<< " target's momentum: " << fTargetMomentum << "\n\n";
cout << " number of collisions: " << fNCollision << "\n"
<< " parameter file name: " << fParamFileName << "\n"
<< " HE model: " << fHEModel;
if ( fHEModel == 0 )
cout << " (EPOS LHC) \n";
else if ( fHEModel == 1 )
cout << " (EPOS) \n";
else if ( fHEModel == 2 )
cout << " (QGSJET01) \n";
else if ( fHEModel == 3 )
cout << " (Gheisha)\n ";
else if ( fHEModel == 4 )
cout << " (Pythia)\n ";
else if ( fHEModel == 5 )
cout << " (Hijing)\n ";
else if ( fHEModel == 6 )
cout << " (Sibyll)\n ";
else if ( fHEModel == 7 )
cout << " (QGSJETII) \n";
else if ( fHEModel == 8 )
cout << " (Phojet) \n";
else
cout << " (unknown) \n";
cout << endl;
cout.setf(ios::showpoint);
cout.setf(ios::fixed);
cout.precision(3);
}
string
CRMCoptions::GetOutputTypeEnding() const
{
switch (fOutputMode) {
case eHepMC:
return ".hepmc";
break;
case eHepMCGZ:
return ".hepmc.gz";
break;
case eLHE:
return ".lhe";
break;
case eLHEGZ:
return ".lhe.gz";
break;
case eROOT:
return ".root";
break;
case eNDST:
return ".ndst.root";
break;
}
return ".unknown";
}
string
CRMCoptions::GetOutputFileName() const
{
// open output file and connect tree
// check if ROOT_OUT is set (if not assume PWD)
#warning replace TSystem with boost !!!
const char* rootOutDir = gSystem->Getenv("ROOT_OUT");
if (rootOutDir == 0 ) {
rootOutDir = gSystem->Getenv("PWD");
}
#warning rename, by removing ROOT...
// open ROOT output file
ostringstream rootFileName;
rootFileName << rootOutDir << "/" << fPrefix << "_";
switch (fHEModel) {
case 0: rootFileName << "eposlhc"; break;
case 1: rootFileName << "epos"; break;
case 2: rootFileName << "qgsjet"; break;
case 3: rootFileName << "gheisha"; break;
case 4: rootFileName << "pythia"; break;
case 5: rootFileName << "hijing"; break;
case 6: rootFileName << "sibyll"; break;
case 7: rootFileName << "qgsjetII"; break;
case 8: rootFileName << "phojet"; break;
default:
cout << " rootOut: error - unknown model " << fHEModel << endl;
cout << " exit ..." << endl;
exit(1);
break;
}
rootFileName << "_" << fSeed << "_";
switch (fProjectileId) {
case 120 : rootFileName << "pi+"; break;
case -120 : rootFileName << "pi+"; break;
case 1 : rootFileName << "p+"; break;
case -1 : rootFileName << "p-"; break;
case 12 : rootFileName << "C"; break;
case 207 : rootFileName << "Pb"; break;
default:
cout << " rootOut: error - unkown beam's Id" << fProjectileId << endl;
cout << " exit ..." << endl;
exit(1);
break;
}
switch (fTargetId) {
case 1 : rootFileName << "p+"; break;
case -1 : rootFileName << "p-"; break;
case 12 : rootFileName << "C"; break;
case 207 : rootFileName << "Pb"; break;
default:
cout << " rootOut: error - unkown target's Id" << fTargetId << endl;
cout << " exit ..." << endl;
exit(1);
break;
}
rootFileName << "_" << fProjectileMomentum
<< GetOutputTypeEnding();
return rootFileName.str();
}
#ifndef _CRMCoptions_h_
#define _CRMCoptions_h_
#include <string>
class CRMCoptions {
template<typename> friend class CRMC;
friend class OutputPolicyROOT;
friend class OutputPolicyHepMC;
CRMCoptions();
public:
enum EOutputMode {
eHepMC,
eHepMCGZ,
eLHE,
eLHEGZ,
eROOT,
eNDST
};
CRMCoptions(int argc, char** argv);
bool OptionsError() const { return fError; }
void DumpConfig() const;
EOutputMode GetOutputMode() const { return fOutputMode; }
std::string GetOutputTypeEnding() const;
std::string GetOutputFileName() const;
protected:
bool fError;
EOutputMode fOutputMode;
// real data members
int fNCollision;
int fSeed;
int fProjectileId;
int fTargetId;
int fHEModel;
double fProjectileMomentum;
double fTargetMomentum;
std::string fParamFileName;
std::string fPrefix;
private:
void CheckEnvironment();
void ParseOptions(int argc, char** argv);
};
#endif
#include <OutputPolicyHepMC.h>
#include <CRMCoptions.h>
#include <CRMCinterface.h>
// TODO check for has hepmc
#include <HepMC/GenEvent.h>
#include <HepMC/HeavyIon.h>
#include <HepMC/HEPEVT_Wrapper.h>
#include <HepMC/IO_GenEvent.h>
#include <HepMC/IO_HEPEVT.h>
#include <HepMC/PdfInfo.h>
#ifdef HEPMC_HAS_CROSS_SECTION
#include <HepMC/GenCrossSection.h>
#endif
#ifdef HEPMC_HAS_UNITS
#include <HepMC/Units.h>
#endif
#include <iomanip>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <stdlib.h>
#include <string>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#include <boost/iostreams/filter/zlib.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/path.hpp>
#include <CRMCconfig.h> //cmake generated
namespace io = boost::iostreams;
using namespace std;
OutputPolicyHepMC::OutputPolicyHepMC()
{
}
void
OutputPolicyHepMC::InitOutput(const CRMCoptions& cfg)
{
boost::filesystem::path oldFile(cfg.GetOutputFileName());
boost::filesystem::remove(oldFile); //FIXME currently truncate does not seem to work properly in boost
//io::filtering_ostream out; //top to bottom order
fOut = new io::filtering_ostream();
if (cfg.fOutputMode==CRMCoptions::eHepMCGZ)
fOut->push(io::gzip_compressor(io::zlib::best_compression));
fOut->push(io::file_descriptor_sink(cfg.GetOutputFileName()), ios_base::trunc);
// Instantiate an IO strategy for reading from HEPEVT.
hepevtio = new HepMC::IO_HEPEVT();
// Instantiate an IO strategy to write the data to file
ascii_out = new HepMC::IO_GenEvent(*fOut);
//We need to explicitly pass this information to the
// HEPEVT_Wrapper.
HepMC::HEPEVT_Wrapper::set_max_number_entries(gCRMC_data.fMaxParticles); //as used in crmc-aaa.f!!!
HepMC::HEPEVT_Wrapper::set_sizeof_real(8); //as used in crmc-aaa.f!!!
}
void