Index: trunk/CMakeLists.txt
===================================================================
--- trunk/CMakeLists.txt (revision 14)
+++ trunk/CMakeLists.txt (revision 15)
@@ -1,331 +1,349 @@
###########################################################################
#
# Copyright 2010
#
# This file is part of Starlight.
#
# Starlight is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# Starlight is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Starlight. If not, see <http://www.gnu.org/licenses/>.
#
###########################################################################
#
# File and Version Information:
# $Rev:: 247 $: revision of last commit
# $Author:: butter $: author of last commit
# $Date:: 2016-03-05 00:59:24 +0000 #$: date of last commit
#
# Description:
# Starlight build file
#
#
###########################################################################
# check if cmake has the required version
cmake_minimum_required(VERSION 2.6.0 FATAL_ERROR)
-
# set verbosity
set(CMAKE_VERBOSE_MAKEFILE 0) # if set to 1 compile and link commands are displayed during build
# the same effect can be achieved by calling 'make VERBOSE=1'
# The version number. 9999 indicates trunk
set (Starlight_VERSION_MAJOR 9999)
set (Starlight_VERSION_MINOR 1)
set (Starlight_VERSION_MINOR_MINOR 0)
# define project
project(starlight)
#find_package (Threads)
# load some common cmake macros
# set path, where to look first for cmake modules, before ${CMAKE_ROOT}/Modules/ is checked
set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake_modules")
message(STATUS "Using cmake module path '${CMAKE_MODULE_PATH}'")
include(CommonMacros)
# force out-of-source builds.
enforce_out_of_source_build()
# warn user if system is not UNIX
if(NOT UNIX)
message(FATAL_ERROR "This is an unsupported system.")
endif()
message(STATUS "Detected host system '${CMAKE_HOST_SYSTEM_NAME}' version '${CMAKE_HOST_SYSTEM_VERSION}' architecture '${CMAKE_HOST_SYSTEM_PROCESSOR}'")
message(STATUS "Compiling for system '${CMAKE_SYSTEM_NAME}' version '${CMAKE_SYSTEM_VERSION}' architecture '${CMAKE_SYSTEM_PROCESSOR}'")
-option (CPP11 "Enable compilation with C++11 features" OFF)
+option (CPP11 "Enable compilation with C++11 features" ON)
# define build types
# set a default build type for single-configuration CMake generators, if no build type is set.
set(CMAKE_BUILD_TYPE Debug)
if(NOT CMAKE_CONFIGURATION_TYPES AND NOT CMAKE_BUILD_TYPE)
message(STATUS "No build type was specified. Setting build type to 'Release'.")
set(CMAKE_BUILD_TYPE Release)
endif()
# common compiler flags
if (CMAKE_COMPILER_IS_GNUCC)
execute_process(COMMAND ${CMAKE_CXX_COMPILER} -dumpversion OUTPUT_VARIABLE GCC_VERSION)
message(STATUS "GCC_VERSTION")
message(STATUS ${GCC_VERSION})
if (GCC_VERSION VERSION_GREATER 4.6 OR GCC_VERSION VERSION_EQUAL 4.6)
message(STATUS "GCC_VERSION>=4.6")
if(CPP11)
set(CMAKE_CXX_FLAGS "-Wall -Wextra -Werror -Wno-error=unused-but-set-variable -Wno-error=unused-but-set-parameter -std=c++11")
message(STATUS "Enabling usage of C++11 features")
else()
set(CMAKE_CXX_FLAGS "-Wall -Wextra -Werror -Wno-error=unused-but-set-variable -Wno-error=unused-but-set-parameter")
endif()
else()
message(STATUS "GCC_VERSION<4.6")
set(CMAKE_CXX_FLAGS "-Wall -Wextra -Werror -Wno-error=unused-but-set-variable -Wno-error=unused-but-set-parameter")
if(CPP11)
message(WARNING "C++11 feautures not supported for your compiler")
endif()
endif()
else()
message(STATUS "Not GCC")
set(CMAKE_CXX_FLAGS "-Wall -Wextra -Werror")
if(CPP11)
message(WARNING "C++11 feautures not supported for your compiler")
endif()
endif()
# flags for specific build types
set(CMAKE_CXX_FLAGS_DEBUG "-g")
set(CMAKE_CXX_FLAGS_RELEASE "-O3")
set(CMAKE_CXX_LDFLAGS_DEBUG "-g")
# report global build settings
message(STATUS "Using CXX compiler '${CMAKE_CXX_COMPILER}'")
message(STATUS "Using CXX general compiler flags '${CMAKE_CXX_FLAGS}'")
foreach(_BUILD_TYPE "DEBUG" "MINSIZEREL" "RELEASE" "RELWITHDEBINFO")
message(STATUS "Using CXX compiler flags '${CMAKE_CXX_FLAGS_${_BUILD_TYPE}}' for build type ${_BUILD_TYPE}")
endforeach()
message(STATUS "Build type is '${CMAKE_BUILD_TYPE}'")
# redirect output files
#set(LIBRARY_OUTPUT_PATH "${CMAKE_SOURCE_DIR}/lib")
message(STATUS "Using library output path '${LIBRARY_OUTPUT_PATH}'")
#set(EXECUTABLE_OUTPUT_PATH "${CMAKE_SOURCE_DIR}/bin")
message(STATUS "Using executable output path '${EXECUTABLE_OUTPUT_PATH}'")
# make CMAKE_SOURCE_DIR accessible in source code via predefined macro CMAKE_SOURCE_DIR
if(CMAKE_SOURCE_DIR)
add_definitions(-D'CMAKE_SOURCE_DIR=\"${CMAKE_SOURCE_DIR}\"')
else()
add_definitions(-D'CMAKE_SOURCE_DIR=\"\"')
endif()
# make SVN version string accessible in source code via predefined macro SVN_VERSION
find_package(Subversion)
if(Subversion_FOUND)
# unfortunately CMAKE only parses 'svn info'
find_program(SVNVERSION_EXECUTABLE
svnversion
)
if(NOT SVNVERSION_EXECUTABLE)
message(STATUS "Could not find subversion command 'svnversion'. Repository version unknown.")
else()
execute_process(
COMMAND ${SVNVERSION_EXECUTABLE} "${CMAKE_SOURCE_DIR}"
OUTPUT_VARIABLE SVN_VERSION
RESULT_VARIABLE _SVNVERSION_RETURN
OUTPUT_STRIP_TRAILING_WHITESPACE)
if(NOT ${_SVNVERSION_RETURN})
message(STATUS "Subversion repository revision is '${SVN_VERSION}'")
else()
message(STATUS "Error running 'svnversion'. Repository version unknown.")
set(SVN_VERSION "")
endif()
endif()
else()
message(STATUS "Could not find subversion installation. Repository version unknown.")
endif()
if(SVN_VERSION)
add_definitions(-D'SVN_VERSION=\"${SVN_VERSION}\"')
else()
add_definitions(-D'SVN_VERSION=\"\"')
endif()
# setup doxygen
find_package(Doxygen)
if(NOT DOXYGEN_FOUND)
message(WARNING "Cannot find Doxygen. No HTML documentation will be generated.")
else()
set(DOXYGEN_TARGET "doxygen")
set(DOXYGEN_DOC_DIR "${CMAKE_SOURCE_DIR}/doxygen")
set(DOXYGEN_CONF "${CMAKE_SOURCE_DIR}/starlightDoxyfile.conf")
message(STATUS "Run 'make ${DOXYGEN_TARGET}' to create Doxygen documentation files in '${DOXYGEN_DOC_DIR}'")
add_custom_target(${DOXYGEN_TARGET}
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_CONF}
DEPENDS ${DOXYGEN_CONF}
WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
)
endif()
+option (ENABLE_HEPMC3 "Enable compilation against hepmc3 (necessary for hepmc3 output)" OFF)
+if(ENABLE_HEPMC3)
+ find_package(HepMC3)
+ if(HEPMC3_FOUND)
+ set(optionalLibs ${optionalLibs} ${HEPMC3_LIBRARY})
+ option(ENABLE_HEPMC3 "Should we use the HepMC3 library" ON)
+ else()
+ option(ENABLE_HEPMC3 "Should we use the HepMC3 library" OFF)
+ endif()
+endif()
# setup Pythia 8
option (ENABLE_PYTHIA "Enable compilation against pythia (necessary for certain processes)" OFF)
if(ENABLE_PYTHIA)
find_package(Pythia8)
if(PYTHIA8_FOUND)
set(optionalLibs ${optionalLibs} ${PYTHIA8_LIBRARY})
# find_package(LHAPDF REQUIRED) # implemented for dummy version in Pythia8
# set(optionalLibs ${optionalLibs} ${LHAPDF_LIBRARIES})#liblhapdfdummy
# removed from v8.2, if you want to add your own lhapdf set, uncomment
# find_package + set - jb 05192015
option(ENABLE_PYTHIA "Should we use the Pythia8 library" ON)
else()
option(ENABLE_PYTHIA "Should we use the Pythia8 library" OFF)
endif()
endif()
# setup Pythia 6
option (ENABLE_PYTHIA6 "Enable compilation against pythia 6 (necessary for certain processes)" OFF)
if(ENABLE_PYTHIA6)
find_package(Pythia6 REQUIRED)
if(PYTHIA6_FOUND)
set(optionalLibs ${optionalLibs} ${PYTHIA6_LIBRARY})
option (ENABLE_PYTHIA6 "Enable compilation against pythia 6 (necessary for certain processes)" ON)
include_directories(pythia6)
else(PYTHIA6_FOUND)
option (ENABLE_PYTHIA6 "Enable compilation against pythia 6 (necessary for certain processes)" OFF)
endif(PYTHIA6_FOUND)
endif()
# setup DPMJET
option (ENABLE_DPMJET "Enable compilation against DPMJet" OFF)
if(ENABLE_DPMJET)
find_package(DPMJet REQUIRED)
if(DPMJET_FOUND)
option (ENABLE_DPMJET "Enable compilation against DPMJet" ON)
else(DPMJET_FOUND)
option (ENABLE_DPMJET "Enable compilation against DPMJet" OFF)
endif(DPMJET_FOUND)
endif(ENABLE_DPMJET)
# set include directories
set(INCLUDE_DIRECTORIES
${CMAKE_SOURCE_DIR}/include
${PROJECT_BINARY_DIR}
${PYTHIA8_INCLUDE_DIR}#uncommented 05192015
)
include_directories(${INCLUDE_DIRECTORIES})
# Set our source files, include the generated dictionary
set(SOURCES
src/bessel.cpp
src/beam.cpp
src/inputParameters.cpp
src/beambeamsystem.cpp
src/starlightparticle.cpp
src/gammaaluminosity.cpp
src/randomgenerator.cpp
src/nucleus.cpp
src/eventchannel.cpp
src/gammaavm.cpp
src/gammagammasingle.cpp
src/photonNucleusCrossSection.cpp
src/wideResonanceCrossSection.cpp
src/narrowResonanceCrossSection.cpp
src/readinluminosity.cpp
src/twophotonluminosity.cpp
src/gammagammaleptonpair.cpp
src/upcevent.cpp
src/vector3.cpp
src/lorentzvector.cpp
src/filewriter.cpp
src/eventfilewriter.cpp
src/starlightparticlecodes.cpp
src/nBodyPhaseSpaceGen.cpp
src/inputParser.cpp
src/incoherentPhotonNucleusLuminosity.cpp
src/incoherentVMCrossSection.cpp
# eSTARlight
src/eXevent.cpp
src/gammaeluminosity.cpp
src/e_wideResonanceCrossSection.cpp
src/e_narrowResonanceCrossSection.cpp
src/e_starlight.cpp
src/e_starlightStandalone.cpp
)
+
+if (ENABLE_HEPMC3)
+ set (SOURCES
+ ${SOURCES}
+ src/hepmc3writer.cpp
+ )
+ include_directories(${HEPMC3_INCLUDE_DIR})
+ add_definitions( -DHEPMC3_ON )
+endif()
if(ENABLE_PYTHIA)
set (SOURCES
${SOURCES}
#src/PythiaStarlight.cpp
src/pythiadecayer.cpp
)
include_directories(${PYTHIA8_INCLUDE_DIR})
endif()
if(ENABLE_PYTHIA6)
set (SOURCES
${SOURCES}
src/starlightpythia.cpp
src/spectrum.cpp
src/spectrumprotonnucleus.cpp
)
endif()
if(ENABLE_DPMJET)
set (SOURCES
${SOURCES}
src/starlightdpmjet.cpp
src/spectrum.cpp
src/spectrumprotonnucleus.cpp
)
endif()
# add Starlight library to the build system
set(THIS_LIB "Starlib")
add_library(${THIS_LIB} STATIC ${SOURCES})
#make_shared_library("${THIS_LIB}" "${SOURCES}"
# "${PYTHIA8_LIBRARY}"
# "${LHAPDF_LIBRARIES}"
#)
if(ENABLE_DPMJET)
enable_language(Fortran)
set(DPMJET_LIB "DpmJetLib")
message(STATUS "DPMJet objects: ${DPMJET_OBJECTS}")
add_library(${DPMJET_LIB} STATIC dpmjet/dpmjetint.f ${DPMJET_OBJECTS})
set(optionalLibs ${optionalLibs} ${DPMJET_LIB})
set(optionalLibs ${optionalLibs} "gfortran")
endif()
if(ENABLE_PYTHIA6)
enable_language(Fortran)
endif()
# add starlight executable to the build system
add_executable(e_starlight src/e_main.cpp)
target_link_libraries(e_starlight Starlib ${optionalLibs})# ${CMAKE_THREAD_LIBS_INIT})
configure_file (
"${PROJECT_SOURCE_DIR}/starlightconfig.h.in"
"${PROJECT_BINARY_DIR}/starlightconfig.h"
)
# Erase xsec values in case changes in code affects the xsec, executed during make process
add_custom_command (TARGET Starlib POST_BUILD COMMAND touch ARGS slight.txt)
add_custom_command (TARGET Starlib POST_BUILD COMMAND cp ARGS slight.txt slight.txt.bak)
add_custom_command (TARGET Starlib POST_BUILD COMMAND echo ARGS '' > slight.txt )
message(STATUS "Cmake did not find any errors. run 'make' to build the project.")
message(STATUS "On multi-core machines 'make -j#', where # is the number of parallel jobs, can speedup compilation considerably.")
Index: trunk/include/filewriter.h
===================================================================
--- trunk/include/filewriter.h (revision 14)
+++ trunk/include/filewriter.h (revision 15)
@@ -1,78 +1,78 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 28 $: revision of last commit
// $Author:: bgrube $: author of last commit
// $Date:: 2010-12-10 18:30:01 +0000 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#ifndef FILEWRITER_H
#define FILEWRITER_H
#include <string>
#include <fstream>
#include "upcevent.h"
#include "eXevent.h"
class fileWriter
{
public:
/** Default constructor */
fileWriter();
/** Constructor with filename */
fileWriter(const std::string& fileName);
/** Destructor */
virtual ~fileWriter();
/** open the file */
int open();
/** open file with given filename */
int open(const std::string& fileName);
/** close the file */
- int close();
+ virtual int close();
/** Set the filename we're writing to */
void setFileName(const std::string& fileName) { _fileName = fileName; }
protected:
/** The file name */
std::string _fileName;
/** The file stream */
std::ofstream _fileStream;
};
#endif // FILEWRITER_H
Index: trunk/include/hepmc3writer.h
===================================================================
--- trunk/include/hepmc3writer.h (revision 0)
+++ trunk/include/hepmc3writer.h (revision 15)
@@ -0,0 +1,34 @@
+
+#ifndef HEPMC3WRITER_H
+#define HEPMC3WRITER_H
+
+#include "inputParameters.h"
+#include "eXevent.h"
+#include "HepMC3/WriterAscii.h"
+
+using HepMC3::WriterAscii;
+
+class hepMC3Writer
+{
+
+ public:
+ /** Default Constructor **/
+ hepMC3Writer();
+ ~hepMC3Writer();
+
+ /** Init HepMC3 Writer with input parameters **/
+ int initWriter(const inputParameters & param);
+
+ /** Write an eX event to file **/
+ int writeEvent(const eXEvent &event, int eventnumber);
+
+ /* closes file */
+ void close(){ _hepmc3_output->close(); };
+
+ private:
+
+ WriterAscii * _hepmc3_output;
+
+};
+
+#endif // HEPMC3WRITER_H
Index: trunk/include/eventfilewriter.h
===================================================================
--- trunk/include/eventfilewriter.h (revision 14)
+++ trunk/include/eventfilewriter.h (revision 15)
@@ -1,72 +1,86 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 97 $: revision of last commit
// $Author:: odjuvsla $: author of last commit
// $Date:: 2012-10-22 22:25:35 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#ifndef EVENTFILEWRITER_H
#define EVENTFILEWRITER_H
+#ifdef HEPMC3_ON
+class hepMC3Writer;
+#endif
#include <string>
#include "filewriter.h"
#include "inputParameters.h"
class eventFileWriter : public fileWriter
{
public:
/** Default constructor */
eventFileWriter();
/** Constructor with name */
eventFileWriter(std::string filename);
/** Write out simulation set up */
int writeInit(inputParameters ¶m );
/** Write an UPC event to file */
int writeEvent(upcEvent &event, int eventnumber);
/** Write an eX event to file */
int writeEvent(eXEvent &event, int eventnumber);
/** Set if we want to write full pythia information */
void writeFullPythiaInfo(bool v) { _writeFullPythia = v; }
+
+ /** Set if we want to write full pythia information */
+ void writeFullHepMC3Info(bool v) { _writeFullHepMC3 = v; }
+
+ /** close the file */
+ int close();
private:
bool _writeFullPythia;
-
+ bool _writeFullHepMC3;
+
+#ifdef HEPMC3_ON
+ hepMC3Writer * _hepmc3writer;
+#endif
+
};
#endif // EVENTFILEWRITER_H
Index: trunk/include/inputParameters.h
===================================================================
--- trunk/include/inputParameters.h (revision 14)
+++ trunk/include/inputParameters.h (revision 15)
@@ -1,405 +1,408 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 276 $: revision of last commit
// $Author:: jnystrand $: author of last commit
// $Date:: 2016-09-13 19:54:42 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#ifndef INPUTPARAMETERS_H
#define INPUTPARAMETERS_H
#include "starlightconstants.h"
#include "inputParser.h"
#include <string>
#include <ostream>
#include <vector>
#include <sstream>
class parameterbase;
class parameterlist
{
public:
parameterlist() : _parameters(0) {}
void add(parameterbase* p) {
_parameters.push_back(p);
}
// Returns a string with a key of the current state of the parameter list
// only
inline std::string validationKey();
private:
std::vector<parameterbase*> _parameters;
};
// Base class for parameters, needed to keep a list of parameters
class parameterbase
{
public:
// Add this to parameter list
parameterbase()
{
_parameters.add(this);
}
virtual std::string validationkey() = 0;
template<typename T>
std::string toString(T v)
{
std::stringstream s;
s << v;
return s.str();
}
inline friend std::ostream& operator<<(std::ostream& os, const parameterbase& par);
// List of all parameters
static parameterlist _parameters;
};
// Need to init the static variable
// parameterlist parameterbase::_parameters;
// The actual parameter class
// validate parameter specifies if the parameter should be a part of the validity check of the current parameters
template<typename T, bool validate>
class parameter : public parameterbase
{
public:
// Constructor
parameter(const std::string &name, T value, bool required = true) :parameterbase(),_name(name), _value(value), _validate(validate), _required(required) {}
parameter &operator=(T v) { _value = v; return *this;}
T* ptr() const {
return const_cast<T*>(&_value);
}
T value() const { return _value; }
std::string name() const { return _name;}
bool required() const { return _required; }
void setValue(T v) { _value = v; }
void setName(std::string name) { _name = name; }
void setRequired(bool r) { _required = r; }
// Validation key for this parameter
std::string validationkey()
{
return (_validate ? _name + ":" + toString(_value) + "-" : std::string(""));
}
template<typename S, bool v>
inline friend std::ostream& operator<<(std::ostream& os, const parameter<S,v>& par);
private:
std::string _name;
T _value; // Value
bool _validate; // true if a change in the parameter invalidates x-sec tables
bool _required; // true if this is required option.
parameter();
};
template<typename S, bool v>
inline std::ostream& operator<<(std::ostream& os, const parameter<S,v>& par)
{
os << par._value;
return os;
}
std::ostream& operator<<(std::ostream& os, const parameterbase& par)
{
os << par._parameters.validationKey();
return os;
}
std::string parameterlist::validationKey()
{
std::stringstream s;
for(unsigned int i = 0; i < _parameters.size(); ++i)
{
s << _parameters[i]->validationkey(); // Will print names and values of validation parameters
}
return s.str();
}
class inputParameters {
public:
inputParameters();
~inputParameters();
bool init();
bool configureFromFile(const std::string &configFileName = "./config/slight.in");
std::string baseFileName () const { return _baseFileName.value(); }
int beam1Z () const { return _beam1Z.value(); } ///< returns atomic number of beam particle 1
unsigned int beam1A () const { return _beam1A.value(); } ///< returns atomic mass number of beam particle 1
int beam2Z () const { return _beam2Z.value(); } ///< returns atomic number of beam particle 2
unsigned int beam2A () const { return _beam2A.value(); } ///< returns atomic mass number of beam particle 2
double targetLorentzGamma () const { return _targetLorentzGamma; } ///< returns Lorentz gamma factor of source in target frame
double beamLorentzGamma () const { return _beamLorentzGamma; } ///< returns Lorentz gamma factor of both beams in beam CMS frame
double beam1LorentzGamma () const { return _beam1LorentzGamma.value(); } ///< returns Lorentz gamma factor of beam 1 in collider frame
double beam2LorentzGamma () const { return _beam2LorentzGamma.value(); } ///< returns Lorentz gamma factor of beam 2 in collider frame
double targetMaxPhotonEnergy () const { return _targetMaxPhotonEnergy; } ///< returns maximum photon energy
double cmsMaxPhotonEnergy () const { return _cmsMaxPhotonEnergy; } ///< returns maximum photon energy
double targetMinPhotonEnergy () const { return _targetMinPhotonEnergy; } ///< returns maximum photon energy
double cmsMinPhotonEnergy () const { return _cmsMinPhotonEnergy; } ///< returns maximum photon energy
double maxW () const { return _maxW.value(); } ///< returns maximum mass W of produced hadronic system [GeV/c^2]
double minW () const { return _minW.value(); } ///< returns minimum mass W of produced hadronic system [GeV/c^2]
unsigned int nmbWBins () const { return _nmbWBins.value(); } ///< returns number of W bins in lookup table
double maxRapidity () const { return _maxRapidity.value(); } ///< returns maximum absolute value of rapidity
unsigned int nmbRapidityBins () const { return _nmbRapidityBins.value(); } ///< returns number of rapidity bins in lookup table
unsigned int nmbEnergyBins () const { return _nmbEnergyBins.value(); } ///< return the number of Egamma bins for eSTARlight lookup
bool ptCutEnabled () const { return _ptCutEnabled.value(); } ///< returns cut in pt
double ptCutMin () const { return _ptCutMin.value(); } ///< returns minimum pt
double ptCutMax () const { return _ptCutMax.value(); } ///< returns maximum pt
bool etaCutEnabled () const { return _etaCutEnabled.value(); } ///< returns cut in eta
double etaCutMin () const { return _etaCutMin.value(); } ///< returns minimum eta
double etaCutMax () const { return _etaCutMax.value(); } ///< returns maximum eta
int productionMode () const { return _productionMode.value(); } ///< returns production mode
unsigned int nmbEvents () const { return _nmbEventsTot.value(); } ///< returns total number of events to generate
int prodParticleId () const { return _prodParticleId.value(); } ///< returns PDG particle ID of produced particle
int randomSeed () const { return _randomSeed.value(); } ///< returns seed for random number generator
int beamBreakupMode () const { return _beamBreakupMode.value(); } ///< returns breakup mode for beam particles
bool interferenceEnabled () const { return _interferenceEnabled.value(); } ///< returns whether interference is taken into account
double interferenceStrength () const { return _interferenceStrength.value(); } ///< returns percentage of interference
double maxPtInterference () const { return _maxPtInterference.value(); } ///< returns maximum p_T for interference calculation [GeV/c]
int nmbPtBinsInterference () const { return _nmbPtBinsInterference.value(); } ///< returns number of p_T bins for interference calculation
double ptBinWidthInterference() const { return _ptBinWidthInterference.value(); } ///< returns width of p_T bins for interference calculation [GeV/c]
double minGammaEnergy () const { return _minGammaEnergy.value(); } ///< returns minimum gamma energy in case of photo nuclear processes [GeV]
double maxGammaEnergy () const { return _maxGammaEnergy.value(); } ///< returns maximum gamma energy in case of photo nuclear processes [GeV]
double minGammaQ2 () const { return _minGammaQ2.value(); } ///< return minimum gamma virtuality
double maxGammaQ2 () const { return _maxGammaQ2.value(); } ///< return maximum gamma virtuality
bool fixedQ2Range () const { return _fixedQ2Range; } ///< return state of Q2 range
unsigned int nmbGammaQ2Bins () const { return _nmbGammaQ2Bins.value(); } ///< return number of gamma q2 bins
std::string pythiaParams () const { return _pythiaParams.value(); } ///< returns parameters to be passed to pythia
bool pythiaFullEventRecord () const { return _pythiaFullEventRecord.value(); } ///< returns if the full pythia event record should be printed
+ bool hepmc3FullEventRecord () const { return _hepmc3FullEventRecord.value(); } ///< returns if the full hepmc3 event record should be printed
int xsecCalcMethod () const { return _xsecCalcMethod.value(); } ///< returns the method used for the x-sec calculation
double axionMass () const { return _axionMass.value(); } ///< returns axion mass //AXION HACK
int bslopeDefinition () const { return _bslopeDefinition.value(); } ///< returns the definition of b-slope
double bslopeValue () const { return _bslopeValue.value(); } ///< returns the value of b-slope
int impulseVM () const { return _impulseVM.value(); } ///< returns the impulseVM value
int quantumGlauber () const { return _quantumGlauber.value(); } ///< returns the quantum glauber value
starlightConstants::particleTypeEnum prodParticleType () const { return _particleType; } ///< returns type of produced particle
starlightConstants::decayTypeEnum prodParticleDecayType() const { return _decayType; } ///< returns decay type of produced particle
starlightConstants::interactionTypeEnum interactionType () const { return _interactionType; } ///< returns interaction type
double protonEnergy () const { return _protonEnergy.value(); }
double electronEnergy () const { return _electronEnergy.value(); }
double inputBranchingRatio () const { return _inputBranchingRatio; }
double targetRadius () const { return _targetR; }
void setBaseFileName (std::string v ) { _baseFileName = v; }
void setBeam1Z (unsigned int v) { _beam1Z = v; } ///< sets atomic number of beam particle 1
void setBeam1A (unsigned int v) { _beam1A = v; } ///< sets atomic mass number of beam particle 1
void setBeam2Z (unsigned int v) { _beam2Z = v; } ///< sets atomic number of beam particle 2
void setBeam2A (unsigned int v) { _beam2A = v; } ///< sets atomic mass number of beam particle 2
void setTargetLorentzGamma (double v) { _targetLorentzGamma = v; } ///< sets Lorentz gamma factor of both beams in beam CMS frame
void setBeamLorentzGamma (double v) { _beamLorentzGamma = v; } ///< sets Lorentz gamma factor of both beams in beam CMS frame
void setBeam1LorentzGamma (double v) { _beam1LorentzGamma = v; } ///< sets Lorentz gamma factor of beam 1 in collider frame
void setBeam2LorentzGamma (double v) { _beam2LorentzGamma = v; } ///< sets Lorentz gamma factor of beam 2 in collider frame
//void setMaxPhotonEnergy (double v) { _maxPhotonEnergy = v ; } ///< sets maximim photon energy
void setMaxW (double v) { _maxW = v; } ///< sets maximum mass W of produced hadronic system [GeV/c^2]
void setMinW (double v) { _minW = v; } ///< sets minimum mass W of produced hadronic system [GeV/c^2]
void setNmbWBins (unsigned int v) { _nmbWBins = v; } ///< sets number of W bins in lookup table
void setMaxRapidity (double v) { _maxRapidity = v; } ///< sets maximum absolute value of rapidity
void setNmbRapidityBins (unsigned int v) { _nmbRapidityBins = v; } ///< sets number of rapidity bins in lookup table
void setNmbEgaBins (unsigned int v) { _nmbEnergyBins = v; } ///< sets number of Ega bins for eSTARlight
void setPtCutEnabled (bool v) { _ptCutEnabled = v; } ///< sets cut in pt
void setPtCutMin (double v) { _ptCutMin = v; } ///< sets minimum pt
void setPtCutMax (double v) { _ptCutMax = v; } ///< sets maximum pt
void setEtaCutEnabled (bool v) { _etaCutEnabled = v; } ///< sets cut in eta
void setEtaCutMin (double v) { _etaCutMin = v; } ///< sets minimum eta
void setEtaCutMax (double v) { _etaCutMax = v; } ///< sets maximum eta
void setProductionMode (int v) { _productionMode = v; } ///< sets production mode
void setNmbEvents (unsigned int v) { _nmbEventsTot = v; } ///< sets total number of events to generate
void setProdParticleId (int v) { _prodParticleId = v; } ///< sets PDG particle ID of produced particle
void setRandomSeed (int v) { _randomSeed = v; } ///< sets seed for random number generator
void setBeamBreakupMode (int v) { _beamBreakupMode = v; } ///< sets breakup mode for beam particles
void setInterferenceEnabled (bool v) { _interferenceEnabled = v; } ///< sets whether interference is taken into account
void setInterferenceStrength (double v) { _interferenceStrength = v; } ///< sets percentage of interference
void setMaxPtInterference (double v) { _maxPtInterference = v; } ///< sets maximum p_T for voiderference calculation [GeV/c]
void setNmbPtBinsInterference (int v) { _nmbPtBinsInterference = v; } ///< sets number of p_T bins for interference calculation
void setPtBinWidthInterference(double v) { _ptBinWidthInterference = v; } ///< sets width of p_T bins for voiderference calculation [GeV/c]
void setMinGammaEnergy (double v) { _minGammaEnergy = v; } ///< sets minimum gamma energy in case of photo nuclear processes [GeV]
void setMaxGammaEnergy (double v) { _maxGammaEnergy = v; } ///< sets maximum gamma energy in case of photo nuclear processes [GeV]
void setMinGammaQ2 (double v) { _minGammaQ2 = v; } ///< sets minimum gamma virtuality in case of photo nuclear processes [GeV]
void setMaxGammaQ2 (double v) { _maxGammaQ2 = v; } ///< sets maximum gamma virtuality in case of photo nuclear processes [GeV]
void setPythiaParams (std::string v) { _pythiaParams = v; } ///< sets parameters to be passed to pythia
void setPythiaFullEventRecord (bool v) { _pythiaFullEventRecord = v; } ///< sets if the full pythia event record should be prvoided
+ void setHepMC3FullEventRecord (bool v) { _hepmc3FullEventRecord = v; } ///< sets if the full hepmc3 event record should be prvoided
void setXsecCalcMethod (int v) { _xsecCalcMethod = v; } ///< sets the method used for the x-sec calculation
void setAxionMass (double v) { _axionMass = v; } ///< sets axion mass //AXION HACK
void setbslopeDefinition (int v) { _bslopeDefinition = v; } ///< sets the definition of b slope
void setbslopeValue (double v) { _bslopeValue = v; } ///< sets the value of b slope
void setimpulseVM (int v) { _impulseVM = v; } ///< sets the value of _impulseVM
void setquantumGlauber (int v) { _quantumGlauber = v; } ///< sets the value of quantum_glauber
void setProdParticleType (starlightConstants::particleTypeEnum v) { _particleType = v; } ///< sets type of produced particle
void setProdParticleDecayType (starlightConstants::decayTypeEnum v) { _decayType = v; } ///< sets decay type of produced particle
void setInteractionType (starlightConstants::interactionTypeEnum v) { _interactionType = v; } ///< sets interaction type
void setProtonEnergy (double v) { _protonEnergy = v; } ///< sets proton energy
void setElectronEnergy (double v) { _electronEnergy = v; } ///< sets electron energy
template<typename T>
inline bool setParameter(std::string expression);
std::ostream& print(std::ostream& out) const; ///< prints parameter summary
std::ostream& write(std::ostream& out) const; ///< writes parameters back to an ostream
std::string parameterValueKey() const; ///< Generates key for the current parameters
private:
// To indicate if the crossection table should be re-calculated if parameter changes
#define VALIDITY_CHECK true
#define NO_VALIDITY_CHECK false
std::string _configFileName; ///< path to configuration file (default = ./config/slight.in)
// config file parameters
parameter<std::string,NO_VALIDITY_CHECK> _baseFileName;
parameter<int,VALIDITY_CHECK> _beam1Z; ///< atomic number of beam particle 1
parameter<unsigned int,VALIDITY_CHECK> _beam1A; ///< atomic mass number of beam particle 1
parameter<int,VALIDITY_CHECK> _beam2Z; ///< atomic number of beam particle 2
parameter<unsigned int,VALIDITY_CHECK> _beam2A; ///< atomic mass number of beam particle 2
parameter<double, VALIDITY_CHECK> _beam1LorentzGamma; ///< Lorentz gamma factor of beam 1 in collider frame
parameter<double, VALIDITY_CHECK> _beam2LorentzGamma; ///< Lorentz gamma factor of beam 2 in collider frame
parameter<double, VALIDITY_CHECK> _maxW; ///< maximum mass W of produced hadronic system [GeV/c^2]
parameter<double, VALIDITY_CHECK> _minW; ///< minimum mass W of produced hadronic system; if set to -1 default value is taken [GeV/c^2]
parameter<unsigned int, VALIDITY_CHECK> _nmbWBins; ///< number of W bins in lookup table
parameter<double, VALIDITY_CHECK> _maxRapidity; ///< maximum absolute value of rapidity
parameter<unsigned int, VALIDITY_CHECK> _nmbRapidityBins; ///< number of rapidity bins in lookup table
parameter<unsigned int, VALIDITY_CHECK> _nmbEnergyBins; ///< number of Egamma bins in lookup table
parameter<bool, VALIDITY_CHECK> _ptCutEnabled; ///< en/disables cut in pt
parameter<double, VALIDITY_CHECK> _ptCutMin; ///< minimum pt, if cut is enabled
parameter<double, VALIDITY_CHECK> _ptCutMax; ///< maximum pt, if cut is enabled
parameter<bool, VALIDITY_CHECK> _etaCutEnabled; ///< en/disables cut in eta
parameter<double, VALIDITY_CHECK> _etaCutMin; ///< minimum eta, if cut is enabled
parameter<double, VALIDITY_CHECK> _etaCutMax; ///< maximum eta, if cut is enabled
parameter<unsigned int, VALIDITY_CHECK> _productionMode; ///< \brief production mode
///<
///< 1 = photon-photon fusion,
///< 2 = narrow vector meson resonance in photon-Pomeron fusion,
///< 3 = Breit-Wigner vector meson resonance in photon-Pomeron fusion
parameter<unsigned int, VALIDITY_CHECK> _nmbEventsTot; ///< total number of events to generate
parameter<unsigned int, VALIDITY_CHECK> _prodParticleId; ///< PDG particle ID of produced particle
parameter<unsigned int, VALIDITY_CHECK> _randomSeed; ///< seed for random number generator
///<
///< 1 = ASCII
///< 2 = GSTARtext,
///< 3 = PAW ntuple (not working)
parameter<unsigned int, VALIDITY_CHECK> _beamBreakupMode; ///< \brief breakup mode for beam particles
///<
///< 1 = hard sphere nuclei (b > 2R),
///< 2 = both nuclei break up (XnXn),
///< 3 = a single neutron from each nucleus (1n1n),
///< 4 = neither nucleon breaks up (with b > 2R),
///< 5 = no hadronic break up (similar to option 1, but with the actual hadronic interaction)
parameter<bool, VALIDITY_CHECK> _interferenceEnabled; ///< if VALIDITY_CHECK, interference is taken into account
parameter<double, VALIDITY_CHECK> _interferenceStrength; ///< percentage of interference: from 0 = none to 1 = full
parameter<double, VALIDITY_CHECK> _maxPtInterference; ///< maximum p_T for interference calculation [GeV/c]
parameter<unsigned int, VALIDITY_CHECK> _nmbPtBinsInterference; ///< number of p_T bins for interference calculation
parameter<double, VALIDITY_CHECK> _ptBinWidthInterference; ///< width of p_T bins for interference calculation [GeV/c]
parameter<double, VALIDITY_CHECK> _protonEnergy;
parameter<double, VALIDITY_CHECK> _electronEnergy;
parameter<double, VALIDITY_CHECK> _minGammaEnergy; ///< minimum gamma energy in case of photo nuclear processes [GeV]
parameter<double, VALIDITY_CHECK> _maxGammaEnergy; ///< maximum gamma energy in case of photo nuclear processes [GeV]
parameter<double, VALIDITY_CHECK> _minGammaQ2; ///< minimum gamma Q2 in case of photo nuclear processes
parameter<double, VALIDITY_CHECK> _maxGammaQ2; ///< maximum gamma Q2 in case of photo nuclear processes
parameter<unsigned int, VALIDITY_CHECK> _nmbGammaQ2Bins; ///< number of gamma q2 bins
parameter<std::string,NO_VALIDITY_CHECK> _pythiaParams; ///< semi-colon separated parameters to pass to pythia, e.g. "mstj(1)=0;paru(13)=0.1"
parameter<bool, NO_VALIDITY_CHECK> _pythiaFullEventRecord; ///< if the full pythia event record should be in the output
+ parameter<bool, NO_VALIDITY_CHECK> _hepmc3FullEventRecord; ///< if the full hepmc3 event record should be in the output
parameter<unsigned int, VALIDITY_CHECK> _xsecCalcMethod; ///< Select x-sec calc method. (0 is standard starlight method, 1 must be used for assym. collisions (e.g. p-A), but is slow)
parameter<double, VALIDITY_CHECK> _axionMass; ///Axion mass//AXION HACK
parameter<unsigned int, VALIDITY_CHECK> _bslopeDefinition; ///< Optional parameter to set different values of slope parameter
parameter<double, VALIDITY_CHECK> _bslopeValue; ///< Value of slope parameter when _bslopeDefinition is set to 1
parameter<unsigned int, VALIDITY_CHECK> _impulseVM; ///< Optional parameter to use impulse approximation (no nuclear effects) for VM cross section.
parameter<unsigned int, VALIDITY_CHECK> _quantumGlauber; ///< Optional parameter. Set = 1 to use Quantum Glauber calculation, rather than Classical Glauber
starlightConstants::particleTypeEnum _particleType;
starlightConstants::decayTypeEnum _decayType;
starlightConstants::interactionTypeEnum _interactionType;
double _targetLorentzGamma; ///< Lorentz gamma factor of the source in target frame, not an input parameter
double _beamLorentzGamma; ///< Lorentz gamma factor of the beams in CMS frame, not an input parameter
double _targetR;
double _cmsMaxPhotonEnergy;
double _cmsMinPhotonEnergy;
double _targetMaxPhotonEnergy;
double _targetMinPhotonEnergy;
double _inputBranchingRatio; ///< Branching ratio defined for each channel
bool _fixedQ2Range;
inputParser _ip;
};
template<typename T>
inline
bool inputParameters::setParameter(std::string expression)
{
return _ip.parseString(expression);
}
inline
std::ostream&
operator <<(std::ostream& out,
const inputParameters& par)
{
return par.print(out);
}
#endif // INPUTPARAMETERS_H
Index: trunk/include/e_starlightStandalone.h
===================================================================
--- trunk/include/e_starlightStandalone.h (revision 14)
+++ trunk/include/e_starlightStandalone.h (revision 15)
@@ -1,80 +1,81 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 263 $: revision of last commit
// $Author:: butter $: author of last commit
// $Date:: 2016-06-04 23:03:58 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#ifndef E_STARLIGHTSTANDALONE_H
#define E_STARLIGHTSTANDALONE_H
#include <string>
class upcEvent;
class e_starlight;
class inputParameters;
class e_starlightStandalone {
public:
e_starlightStandalone();
~e_starlightStandalone();
bool init(); ///< reads configuration file and initializes startlight
bool run (); ///< creates output file and runs starlight
const std::string& baseFileName () const { return _baseFileName; }
const std::string& configFileName () const { return _configFileName; } ///< returns path to config file
const std::string& eventDataFileName() const { return _eventDataFileName; } ///< returns path to output file
void setBaseFileName (const std::string& baseFileName ) { _baseFileName = baseFileName; } ///< sets path to base file
void setConfigFileName (const std::string& configFileName ) { _configFileName = configFileName; } ///< sets path to config file
void setEventDataFileName(const std::string& eventDataFileName) { _eventDataFileName = eventDataFileName; } ///< sets path to output file
void boostEvent(eXEvent &e); ///< Boost event from beam CMS to lab system
+ void flipZEvent(eXEvent &e); ///< flip z to -z for all output particles
private:
std::string _baseFileName; ///< path to base filename
std::string _configFileName; ///< path to configuration file
std::string _eventDataFileName; ///< path to output file
e_starlight* _starlight; ///< pointer to starlight instance
inputParameters* _inputParameters; ///< pointer to parameter instance
unsigned int _nmbEventsTot; ///< total number of events to generate (taken from configuration file)
unsigned int _nmbEventsPerFile; ///< maximum number of events written to a single file (not yet implemented)
};
#endif // E_STARLIGHTSTANDALONE_H
Index: trunk/src/e_starlight.cpp
===================================================================
--- trunk/src/e_starlight.cpp (revision 14)
+++ trunk/src/e_starlight.cpp (revision 15)
@@ -1,343 +1,344 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2017
//
// This file is part of estarlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 263 $: revision of last commit
// $Author:: mlomnitz $: author of last commit
// $Date:: 02/28/2017 #$: date of last commit
//
// Description:
//
// Main interface for eSTARlight code inherited from STARlight
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include <cstdlib>
#include "starlightconfig.h"
#ifdef ENABLE_PYTHIA
#include "PythiaStarlight.h"
#endif
#ifdef ENABLE_DPMJET
#include "starlightdpmjet.h"
#endif
#ifdef ENABLE_PYTHIA6
#include "starlightpythia.h"
#endif
#include "reportingUtils.h"
#include "inputParameters.h"
#include "eventchannel.h"
#include "gammagammaleptonpair.h"
#include "gammagammasingle.h"
#include "gammaavm.h"
#include "twophotonluminosity.h"
#include "gammaaluminosity.h"
#include "gammaeluminosity.h"
#include "incoherentPhotonNucleusLuminosity.h"
#include "upcevent.h"
#include "e_starlight.h"
using namespace std;
using namespace starlightConstants;
e_starlight::e_starlight() :
_beamSystem (0),
_eventChannel (0),
_nmbEventsPerFile (100),
_isInitialised (false),
_inputParameters (0)
{ }
e_starlight::~e_starlight()
{ }
bool
e_starlight::init()
{
// ???
if(Starlight_VERSION_MAJOR == 9999)
{
cout << endl << "#########################################" << endl
<< " Initialising Starlight version: trunk..." << endl
<< "#########################################" << endl << endl;
}
else
{
cout << endl << "#########################################" << endl
<< " Initialising Starlight version: " << Starlight_VERSION_MAJOR << "."
<< Starlight_VERSION_MINOR << "." << Starlight_VERSION_MINOR_MINOR << "..." << endl
<< "#########################################" << endl << endl;
}
_nmbEventsPerFile = _inputParameters->nmbEvents(); // for now we write only one file...
_beamSystem = new beamBeamSystem(*_inputParameters);
// This sets the precision used by cout for printing
cout.setf(ios_base::fixed,ios_base::floatfield);
cout.precision(3);
std::string _baseFileName;
_baseFileName = _inputParameters->baseFileName();
_lumLookUpTableFileName = _baseFileName + ".txt";
const bool lumTableIsValid = luminosityTableIsValid();
// Do some sanity checks of the input parameters here.
if( _inputParameters->beam1A() != 0 && _inputParameters->beam2A() != 0 ){
- printErr << endl << "One of the two beams must be an electron in eSTARlight" << endl;
+ // printErr << endl << "One of the two beams must be an electron in eSTARlight" << endl;
return false;
}
if( _inputParameters->beam1A() == 0 && _inputParameters->beam2A() == 0 ){
- printErr << endl << "Only one of the two beams can be electron in eSTARlight"<< endl;
+ // printErr << endl << "Only one of the two beams can be electron in eSTARlight"<< endl;
return false;
}
if( _inputParameters->beam1A() == 0 && std::abs(_inputParameters->beam1Z()) != 1 ){
- printErr << endl << "Beam 1 should be electron, but has wrong Z"<< endl;
+ // printErr << endl << "Beam 1 should be electron, but has wrong Z"<< endl;
return false;
}
if( _inputParameters->beam2A() == 0 && _inputParameters->beam2Z() != std::abs(1) ){
- printErr << endl << "Beam 1 should be electron, but has wrong Z"<< endl;
+ // printErr << endl << "Beam 1 should be electron, but has wrong Z"<< endl;
return false;
}
bool createChannel = true;
switch (_inputParameters->interactionType()) {
// Photon-photon scattering is not implemented in eSTARlight as of yet
/*
case PHOTONPHOTON:
if (!lumTableIsValid) {
printInfo << "creating luminosity table for photon-photon channel" << endl;
twoPhotonLuminosity(*_inputParameters, _beamSystem->beam1(), _beamSystem->beam2());
}
break;
*/
case E_PHOTONPOMERONNARROW: // narrow and wide resonances use
case E_PHOTONPOMERONWIDE: // the same luminosity function
if (!lumTableIsValid) {
printInfo << "creating luminosity table for coherent photon-Pomeron channel" <<endl;
photonElectronLuminosity lum(*_inputParameters, *_beamSystem);
}
break;
// Incoherent electro-production is still pending
/* case PHOTONPOMERONINCOHERENT: // narrow and wide resonances use
if (!lumTableIsValid) {
printInfo << "creating luminosity table for incoherent photon-Pomeron channel" << endl;
incoherentPhotonNucleusLuminosity lum(*_inputParameters, *_beamSystem);
}
break;*/
default:
{
printWarn << "unknown interaction type '" << _inputParameters->interactionType() << "'."
<< " cannot initialize eSTARlight." << endl;
return false;
}
}
if(createChannel)
{
if (!createEventChannel())
return false;
}
_isInitialised = true;
return true;
}
eXEvent
e_starlight::produceEvent()
{
if (!_isInitialised) {
- printErr << "trying to generate event but Starlight is not initialised. aborting." << endl;
+ // printErr << "trying to generate event but Starlight is not initialised. aborting." << endl;
exit(-1);
}
return _eventChannel->e_produceEvent();
}
bool
e_starlight::luminosityTableIsValid() const
{
printInfo << "using random seed = " << _inputParameters->randomSeed() << endl;
ifstream lumLookUpTableFile(_lumLookUpTableFileName.c_str());
lumLookUpTableFile.precision(15);
if ((!lumLookUpTableFile) || (!lumLookUpTableFile.good())) {
return false;
}
unsigned int beam1Z, beam1A, beam2Z, beam2A;
double beamLorentzGamma = 0;
double maxW = 0, minW = 0;
unsigned int nmbWBins;
double maxRapidity = 0;
unsigned int nmbRapidityBins;
int productionMode, beamBreakupMode;
//Remove interference options for eX collisions (no logner symmetric)
// bool interferenceEnabled = false;
// double interferenceStrength = 0;
bool coherentProduction = false;
double incoherentFactor = 0, deuteronSlopePar = 0, maxPtInterference = 0;
int nmbPtBinsInterference;
std::string validationKey;
if (!(lumLookUpTableFile
>> validationKey
>> beam1Z >> beam1A
>> beam2Z >> beam2A
>> beamLorentzGamma
>> maxW >> minW >> nmbWBins
>> maxRapidity >> nmbRapidityBins
>> productionMode
>> beamBreakupMode
>> maxPtInterference
>> nmbPtBinsInterference
>> coherentProduction >> incoherentFactor
>> deuteronSlopePar))
// cannot read parameters from lookup table file
return false;
std::string validationKeyEnd;
while(!lumLookUpTableFile.eof())
{
lumLookUpTableFile >> validationKeyEnd;
}
lumLookUpTableFile.close();
return (validationKey == _inputParameters->parameterValueKey() && validationKeyEnd == validationKey);
return true;
}
bool
e_starlight::createEventChannel()
{
switch (_inputParameters->prodParticleType()) {
case ELECTRON:
case MUON:
case TAUON:
case TAUONDECAY:
{
_eventChannel = new Gammagammaleptonpair(*_inputParameters, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammagammaleptonpair event channel." << endl;
return false;
}
}
case A2: // jetset/pythia
case ETA: // jetset/pythia
case ETAPRIME: // jetset/pythia
case ETAC: // jetset/pythia
case F0: // jetset/pythia
{
-#ifdef ENABLE_PYTHIA
+
+ //#ifdef ENABLE_PYTHIA
// PythiaOutput = true;
- cout<<"Pythia is enabled!"<<endl;
+ // cout<<"Pythia is enabled!"<<endl;
// return true;
-#else
- printWarn << "Starlight is not compiled against Pythia8; "
- << "jetset event channel cannot be used." << endl;
- return false;
-#endif
+//#else
+// printWarn << "Starlight is not compiled against Pythia8; "
+// << "jetset event channel cannot be used." << endl;
+// return false;
+//#endif
}
case F2:
case F2PRIME:
case ZOVERZ03:
case AXION: // AXION HACK
{
// #ifdef ENABLE_PYTHIA
cout<<" This is f2, f2prim, rho^0 rho^0, or axion "<<endl;
_eventChannel= new Gammagammasingle(*_inputParameters, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammagammasingle event channel." << endl;
return false;
}
}
case RHO:
case RHOZEUS:
case FOURPRONG:
case OMEGA:
case PHI:
case JPSI:
case JPSI2S:
case JPSI2S_ee:
case JPSI2S_mumu:
case JPSI_ee:
case JPSI_mumu:
case UPSILON:
case UPSILON_ee:
case UPSILON_mumu:
case UPSILON2S:
case UPSILON2S_ee:
case UPSILON2S_mumu:
case UPSILON3S:
case UPSILON3S_ee:
case UPSILON3S_mumu:
{
if (_inputParameters->interactionType() == E_PHOTONPOMERONNARROW) {
_eventChannel = new e_Gammaanarrowvm(*_inputParameters, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammaanarrowvm event channel." << endl;
return false;
}
}
if (_inputParameters->interactionType() == E_PHOTONPOMERONWIDE) {
_eventChannel = new e_Gammaawidevm(*_inputParameters, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammaawidevm event channel." << endl;
return false;
}
}
printWarn << "interaction type '" << _inputParameters->interactionType() << "' "
<< "cannot be used with particle type '" << _inputParameters->prodParticleType() << "'. "
<< "cannot create event channel." << endl;
return false;
}
default:
{
printWarn << "unknown event channel '" << _inputParameters->prodParticleType() << "'."
<< " cannot create event channel." << endl;
return false;
}
}
}
Index: trunk/src/e_starlightStandalone.cpp
===================================================================
--- trunk/src/e_starlightStandalone.cpp (revision 14)
+++ trunk/src/e_starlightStandalone.cpp (revision 15)
@@ -1,187 +1,188 @@
//////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 270 $: revision of last commit
// $Author:: jnystrand $: author of last commit
// $Date:: 2016-07-08 16:31:51 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <chrono>
#include "reportingUtils.h"
#include "e_starlight.h"
#include "inputParameters.h"
#include "eventfilewriter.h"
#include "e_starlightStandalone.h"
using namespace std;
e_starlightStandalone::e_starlightStandalone()
: _configFileName ("slight.in"),
_starlight (0),
_nmbEventsTot (1),
_nmbEventsPerFile (_nmbEventsTot)
{ }
e_starlightStandalone::~e_starlightStandalone()
{ }
bool
e_starlightStandalone::init()
{
_inputParameters = new inputParameters();
// read input parameters from config file
_inputParameters->configureFromFile(_configFileName);
if (!_inputParameters->init()) {
printWarn << "problems initializing input parameters. cannot initialize starlight." << endl;
return false;
}
// copy input file to one with baseFileName naming scheme
std::string inputCopyName, _baseFileName;
_baseFileName = _inputParameters->baseFileName();
if (_baseFileName != "slight") {
inputCopyName = _baseFileName +".in";
ofstream inputCopyFile;
inputCopyFile.open(inputCopyName.c_str());
std::ifstream infile(_configFileName.c_str());
if ((!infile) || (!infile.good()))
{
return -1;
}
int lineSize = 256;
char tmp[lineSize];
while (!infile.getline(tmp, lineSize).eof())
{
cout << tmp << endl;
inputCopyFile << tmp << endl;
}
inputCopyFile.close();
}
// get the number of events
// for now we write everything to one file
_nmbEventsTot = _inputParameters->nmbEvents();
_nmbEventsPerFile = _nmbEventsTot;
// create the starlight object
_starlight = new e_starlight();
// give starlight the input parameters
_starlight->setInputParameters(_inputParameters);
// initialize starlight
return _starlight->init();
}
bool
e_starlightStandalone::run()
{
if (!_starlight) {
printWarn << "null pointer to starlight object. make sure that init() was called. "
<< "cannot generate events." << endl;
return false;
}
// open output file
eventFileWriter fileWriter;
fileWriter.writeFullPythiaInfo(_inputParameters->pythiaFullEventRecord());
+ fileWriter.writeFullHepMC3Info(_inputParameters->hepmc3FullEventRecord());
_baseFileName = _inputParameters->baseFileName();
_eventDataFileName = _baseFileName +".out";
fileWriter.open(_eventDataFileName);
//
fileWriter.writeInit(*_inputParameters);
//
printInfo << "generating events:" << endl;
unsigned int nmbEvents = 0;
std::chrono::steady_clock::time_point begin= std::chrono::steady_clock::now();
while (nmbEvents < _nmbEventsTot) {
for (unsigned int iEvent = 0; (iEvent < _nmbEventsPerFile) && (nmbEvents < _nmbEventsTot);
++iEvent, ++nmbEvents) {
progressIndicator(iEvent, _nmbEventsTot, true, 4);
eXEvent event = _starlight->produceEvent();
// Boost event from back to lab reference frame
boostEvent(event);
- // flipZEvent(event);
+ flipZEvent(event);
fileWriter.writeEvent(event, iEvent);
}
}
std::chrono::steady_clock::time_point end= std::chrono::steady_clock::now();
float running_total = 1E-3*std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
cout<<"Total time "<<running_total<<" s ("<<1E3*running_total/nmbEvents<<" ms/ev)"<<endl;
fileWriter.close();
if( _starlight->nmbAttempts() == 0 )return true;
double _branchingRatio = _inputParameters->inputBranchingRatio();
printInfo << "number of attempts = " << _starlight->nmbAttempts() << ", "
<< "number of accepted events = " << _starlight->nmbAccepted() << endl;
double selectedCrossSection =
((double)_starlight->nmbAccepted()/_starlight->nmbAttempts())*_branchingRatio*_starlight->getTotalCrossSection();
if (selectedCrossSection > 1.){
cout<< " The cross section of the generated sample is "<<selectedCrossSection<<" barn."<<endl;
} else if (1.E3*selectedCrossSection > 1.){
cout<< " The cross section of the generated sample is "<<1.E3*selectedCrossSection<<" mb."<<endl;
} else if (1.E6*selectedCrossSection > 1.){
cout<< " The cross section of the generated sample is "<<1.E6*selectedCrossSection<<" microbarn."<<endl;
} else if (1.E9*selectedCrossSection > 1.){
cout<< " The cross section of the generated sample is "<<1.E9*selectedCrossSection<<" nanobarn."<<endl;
} else if (1.E12*selectedCrossSection > 1.){
cout<< " The cross section of the generated sample is "<<1.E12*selectedCrossSection<<" picobarn."<<endl;
} else {
cout<< " The cross section of the generated sample is " <<1.E15*selectedCrossSection<<" femtob."<<endl;
}
return true;
}
void e_starlightStandalone::boostEvent(eXEvent &event)
{
// Calculate CMS boost
double rap1 = acosh(_inputParameters->beam1LorentzGamma());
double rap2 = -acosh(_inputParameters->beam2LorentzGamma());
double boost = (rap1+rap2)/2.;
event.boost(boost, rap2); // Boost back to laboratory reference frame. Electron initially in target frame and V.M. in CMS
// Assuming electron is beam1
}
void e_starlightStandalone::flipZEvent(eXEvent &event)
{
event.flipZ(); //Change all z to -z
}
Index: trunk/src/hepmc3writer.cpp
===================================================================
--- trunk/src/hepmc3writer.cpp (revision 0)
+++ trunk/src/hepmc3writer.cpp (revision 15)
@@ -0,0 +1,73 @@
+#include "HepMC3/GenVertex.h"
+#include "HepMC3/GenVertex_fwd.h"
+#include "HepMC3/FourVector.h"
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/GenParticle.h"
+#include "HepMC3/WriterAscii.h"
+#include "HepMC3/Print.h"
+
+#include "inputParameters.h"
+#include "eXevent.h"
+#include "hepmc3writer.h"
+
+using HepMC3::FourVector;
+using HepMC3::Print;
+
+hepMC3Writer::hepMC3Writer()
+{
+}
+
+int hepMC3Writer::initWriter(const inputParameters ¶m)
+{
+
+ std::string hepmc3_filename = param.baseFileName() + ".hepmc";
+
+ /** need to pass the parameters to HepMC3 **/
+ param.beam1Z(); //just do something for now
+ _hepmc3_output = new HepMC3::WriterAscii(hepmc3_filename);
+ return 0;
+}
+
+int hepMC3Writer::writeEvent(const eXEvent &event, int eventnumber)
+{
+
+ /** Make HepMC3 Event and Vertex ... Currently only 1 Vertex per Event **/
+ HepMC3::GenEvent hepmc3_evt(HepMC3::Units::GEV , HepMC3::Units::MM);
+ HepMC3::GenVertexPtr hepmc3_vertex_to_write = std::make_shared<HepMC3::GenVertex>(FourVector(0,0,0,0));
+
+ /** Get starlight particle vector **/
+ const std::vector<starlightParticle> * particle_vector = event.getParticles();
+ lorentzVector gamma_lorentzVec = (*event.getGamma())[0];
+ FourVector hepmc3_gamma_four_vector = FourVector(gamma_lorentzVec.GetPx(),
+ gamma_lorentzVec.GetPy(),
+ gamma_lorentzVec.GetPz(),
+ gamma_lorentzVec.GetE());
+ HepMC3::GenParticlePtr hepmc3_gamma = std::make_shared<HepMC3::GenParticle>( hepmc3_gamma_four_vector, 22, 0 ); // status currently set to 0... need to double check
+ hepmc3_vertex_to_write->add_particle_in( hepmc3_gamma );
+
+ /** Takes e_starlight events and converts to hepmc3 format **/
+ for ( std::vector<starlightParticle>::const_iterator particle_iter = (*particle_vector).begin();
+ particle_iter != (*particle_vector).end();
+ ++particle_iter)
+ {
+ int hepmc3_pid = (*particle_iter).getPdgCode();
+ /** pass to HepMC3 FourVector **/
+ FourVector hepmc3_four_vector = FourVector( (*particle_iter).GetPx(),
+ (*particle_iter).GetPy(),
+ (*particle_iter).GetPz(),
+ (*particle_iter).GetE());
+
+ HepMC3::GenParticlePtr hepmc3_particle = std::make_shared<HepMC3::GenParticle>( hepmc3_four_vector, hepmc3_pid, 0 ); // status currently set to 0... need to double check
+ hepmc3_vertex_to_write->add_particle_out( hepmc3_particle );
+ }
+
+ /** add vertex to HepMC3 Event**/
+ hepmc3_evt.add_vertex( hepmc3_vertex_to_write );
+ _hepmc3_output->write_event(hepmc3_evt);
+ // Print::listing(hepmc3_evt);
+ // Print::content(hepmc3_evt);
+ hepmc3_evt.clear();
+
+ return eventnumber;
+}
+
Index: trunk/src/filewriter.cpp
===================================================================
--- trunk/src/filewriter.cpp (revision 14)
+++ trunk/src/filewriter.cpp (revision 15)
@@ -1,96 +1,94 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 28 $: revision of last commit
// $Author:: bgrube $: author of last commit
// $Date:: 2010-12-10 18:30:01 +0000 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <exception>
#include <cstdlib>
#include "filewriter.h"
using namespace std;
fileWriter::fileWriter()
: _fileName(""),
_fileStream()
{ }
fileWriter::fileWriter(const string& fileName) :
_fileName(fileName)
,_fileStream(fileName.c_str())
{ }
-
fileWriter::~fileWriter()
{ }
-
int fileWriter::open()
{
try
{
_fileStream.open(_fileName.c_str());
_fileStream.precision(15);
_fileStream.setf(ios::fixed);
}
catch (const ios::failure & error)
{
cerr << "I/O exception: " << error.what() << endl;
return EXIT_FAILURE;
}
return 0;
}
int fileWriter::open(const string& fileName)
{
_fileName = fileName;
return open();
}
int fileWriter::close()
{
try
{
_fileStream.close();
}
catch (const ios::failure & error)
{
cerr << "I/O exception: " << error.what() << endl;
return EXIT_FAILURE;
}
return 0;
}
Index: trunk/src/eventfilewriter.cpp
===================================================================
--- trunk/src/eventfilewriter.cpp (revision 14)
+++ trunk/src/eventfilewriter.cpp (revision 15)
@@ -1,189 +1,234 @@
-///////////////////////////////////////////////////////////////////////////
+
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 228 $: revision of last commit
// $Author:: butter $: author of last commit
// $Date:: 2016-01-18 17:10:17 +0000 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
+#ifdef HEPMC3_ON
+#include "hepmc3writer.h"
+#endif
+
+#include <iostream>
+#include <exception>
+#include <cstdlib>
#include "eventfilewriter.h"
#include "starlightparticlecodes.h"
+using namespace std;
+
//______________________________________________________________________________
eventFileWriter::eventFileWriter()
: fileWriter()
,_writeFullPythia(false)
+,_writeFullHepMC3(false)
{ }
-
//______________________________________________________________________________
eventFileWriter::eventFileWriter(std::string filename)
: fileWriter(filename)
{ }
//______________________________________________________________________________
int eventFileWriter::writeInit(inputParameters &_p)
{
_fileStream<<"CONFIG_OPT: "<<_p.productionMode()<<" "<<_p.prodParticleId()<<" "<<_p.nmbEvents()
<<" "<<_p.quantumGlauber()<<" "<<_p.impulseVM()<<" "<<_p.randomSeed()<<std::endl;
_fileStream<<"BEAM_1: "<<_p.beam1Z()<<" "<<_p.beam1A()<<" "<<_p.beam1LorentzGamma()<<std::endl;
_fileStream<<"BEAM_2: "<<_p.beam2Z()<<" "<<_p.beam2A()<<" "<<_p.beam2LorentzGamma()<<std::endl;
_fileStream<<"PHOTON: "<<_p.nmbEnergyBins()<<" "<<_p.fixedQ2Range()<<" "<<_p.nmbGammaQ2Bins()
<<" "<<_p.minGammaQ2()<<" "<<_p.maxGammaQ2()<<std::endl;
+ if ( _writeFullHepMC3 )
+ {
+#ifdef HEPMC3_ON
+ _hepmc3writer = new hepMC3Writer();
+ _hepmc3writer->initWriter(_p);
+#else
+ std::cout << "***HepMC3 file format not written --- eStarlight not compiled with HepMC3 Enabled***" << std::endl;
+ _writeFullHepMC3 = false;
+#endif
+ }
+
return 0;
}
//______________________________________________________________________________
int eventFileWriter::writeEvent(upcEvent &event, int eventnumber)
{
int numberoftracks = 0;
if(_writeFullPythia)
{
numberoftracks = event.getParticles()->size();
}
else
{
for(unsigned int i = 0; i<event.getParticles()->size(); ++i)
{
if(event.getParticles()->at(i).getStatus() >= 0) numberoftracks++;
}
}
-
// sometimes we don't have tracks due to wrongly picked W , check it
if(numberoftracks){
eventnumber++;
_fileStream << "EVENT: " << eventnumber << " " << numberoftracks << " " << 1 << std::endl;
if(event.getGammaEnergies()->size()) _fileStream << "GAMMAENERGIES:";
for(unsigned int n = 0; n < event.getGammaEnergies()->size(); n++)
{
_fileStream << " " << event.getGammaEnergies()->at(n);
}
if(event.getGammaEnergies()->size()) _fileStream<< std::endl;
_fileStream <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<std::endl;
int ipart = 0;
std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
{
if(!_writeFullPythia)
{
if((*part).getStatus() < 0) continue;
}
_fileStream << "TRACK: " << " " << starlightParticleCodes::jetsetToGeant((*part).getPdgCode()) <<" "<< (*part).GetPx() << " " << (*part).GetPy()
<< " "<< (*part).GetPz() << " " << eventnumber << " " << ipart << " " << 0 << " "
<< (*part).getPdgCode();
if(_writeFullPythia)
{
lorentzVector vtx = (*part).getVertex();
_fileStream << " " << vtx.GetPx() << " " << vtx.GetPy() << " " << vtx.GetPz() << " " << vtx.GetE();
_fileStream << " " << (*part).getFirstParent() << " " << (*part).getLastParent() << " " << (*part).getFirstDaughter() << " " << (*part).getLastDaughter() << " " << (*part).getStatus();
}
_fileStream <<std::endl;
}
}
-
+
return 0;
}
//______________________________________________________________________________
int eventFileWriter::writeEvent(eXEvent &event, int eventnumber)
{
int numberoftracks = 0;
if(_writeFullPythia)
{
numberoftracks = event.getParticles()->size();
}
else
{
for(unsigned int i = 0; i<event.getParticles()->size(); ++i)
{
if(event.getParticles()->at(i).getStatus() >= 0) numberoftracks++;
}
}
// sometimes we don't have tracks due to wrongly picked W , check it
if(numberoftracks){
eventnumber++;
_fileStream << "EVENT: " << eventnumber << " " << numberoftracks << " " << 1 << std::endl;
_fileStream <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<std::endl;
- for(unsigned int igam = 0 ; igam < event.getGammaEnergies()->size(); ++igam){
+ for( uint igam = 0 ; igam < event.getGammaEnergies()->size(); ++igam){
_fileStream <<"GAMMA: "<<event.getGammaEnergies()->at(igam)<<" "<<event.getGammaMasses()->at(igam)<<std::endl;
lorentzVector gam = event.getGamma()->at(igam);
// Write the photon 4-vector out to file. Might be used in later iterations, so I keep it here
//_fileStream <<"GAMMA VECTOR: "<<gam.GetPx()<<" "<<gam.GetPy()<<" "<<gam.GetPz()<<" "<<gam.GetE()<<" "<<-temp<<std::endl;
}
- for(unsigned int itarget = 0 ; itarget < event.getTarget()->size(); ++itarget){
+ for( uint itarget = 0 ; itarget < event.getTarget()->size(); ++itarget){
lorentzVector target = event.getTarget()->at(itarget);
_fileStream <<"t: "<<event.getVertext()->at(itarget)<<std::endl;
_fileStream <<"TARGET: "<<target.GetPx()<<" "<<target.GetPy()<<" "<<target.GetPz()<<" "<<target.GetE()<<std::endl;
}
- for(unsigned int iel = 0 ; iel < event.getSources()->size(); ++iel){
+ for( uint iel = 0 ; iel < event.getSources()->size(); ++iel){
lorentzVector el = event.getSources()->at(iel);
_fileStream <<"SOURCE: "<<el.GetPx()<<" "<<el.GetPy()<<" "<<el.GetPz()<<" "<<el.GetE()<<std::endl;
}
int ipart = 0;
std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
{
if(!_writeFullPythia)
{
if((*part).getStatus() < 0) continue;
}
_fileStream << "TRACK: " << " "<< starlightParticleCodes::jetsetToGeant((*part).getPdgCode()) <<" "<< (*part).GetPx() << " " << (*part).GetPy()
<< " "<< (*part).GetPz() << " " << eventnumber << " " << ipart << " " << 0 << " "
<< (*part).getPdgCode();
if(_writeFullPythia)
{
lorentzVector vtx = (*part).getVertex();
_fileStream << " " << vtx.GetPx() << " " << vtx.GetPy() << " " << vtx.GetPz() << " " << vtx.GetE();
_fileStream << " " << (*part).getFirstParent() << " " << (*part).getLastParent() << " " << (*part).getFirstDaughter() << " " << (*part).getLastDaughter() << " " << (*part).getStatus();
}
_fileStream <<std::endl;
}
}
+#ifdef HEPMC3_ON
+ if ( _writeFullHepMC3 ){
+ _hepmc3writer->writeEvent(event,eventnumber);
+ }
+#endif
+
return 0;
}
+int eventFileWriter::close()
+{
+
+#ifdef HEPMC3_ON
+ if ( _writeFullHepMC3 ){
+ _hepmc3writer->close();
+ }
+#endif
+
+ try
+ {
+ _fileStream.close();
+ }
+ catch (const ios::failure & error)
+ {
+ cerr << "I/O exception: " << error.what() << endl;
+ return EXIT_FAILURE;
+ }
+ return 0;
+}
Index: trunk/src/inputParameters.cpp
===================================================================
--- trunk/src/inputParameters.cpp (revision 14)
+++ trunk/src/inputParameters.cpp (revision 15)
@@ -1,703 +1,705 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 276 $: revision of last commit // $Author:: jnystra#$: author of last commit
// $Date:: 2016-09-13 19:54:42 +0100 #$: date of last commit //
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include "reportingUtils.h"
#include "starlightconstants.h"
#include "inputParameters.h"
#include "inputParser.h"
#include "starlightconfig.h"
#include <cmath>
#include <cstring>
#include "randomgenerator.h"
using namespace std;
using namespace starlightConstants;
parameterlist parameterbase::_parameters;
#define REQUIRED true
#define NOT_REQUIRED false
//______________________________________________________________________________
inputParameters::inputParameters()
: _baseFileName ("baseFileName","slight"),
_beam1Z ("BEAM_1_Z",0),
_beam1A ("BEAM_1_A",0),
_beam2Z ("BEAM_2_Z",0),
_beam2A ("BEAM_2_A",0),
_beam1LorentzGamma ("BEAM_1_GAMMA",0),
_beam2LorentzGamma ("BEAM_2_GAMMA",0),
_maxW ("W_MAX",0),
_minW ("W_MIN",0),
_nmbWBins ("W_N_BINS",0),
_maxRapidity ("RAP_MAX",9,NOT_REQUIRED),
_nmbRapidityBins ("RAP_N_BINS",100,NOT_REQUIRED),
_nmbEnergyBins ("EGA_N_BINS",0),
_ptCutEnabled ("CUT_PT",false, NOT_REQUIRED),
_ptCutMin ("PT_MIN",0, NOT_REQUIRED),
_ptCutMax ("PT_MAX",0, NOT_REQUIRED),
_etaCutEnabled ("CUT_ETA",false, NOT_REQUIRED),
_etaCutMin ("ETA_MIN",0, NOT_REQUIRED),
_etaCutMax ("ETA_MAX",0, NOT_REQUIRED),
_productionMode ("PROD_MODE",0),
_nmbEventsTot ("N_EVENTS",0),
_prodParticleId ("PROD_PID",0),
_randomSeed ("RND_SEED",0),
_beamBreakupMode ("BREAKUP_MODE",0),
_interferenceEnabled ("INTERFERENCE",false,NOT_REQUIRED),
_interferenceStrength ("IF_STRENGTH",0,NOT_REQUIRED),
_maxPtInterference ("INT_PT_MAX",0,NOT_REQUIRED),
_nmbPtBinsInterference ("INT_PT_N_BINS",0),
_ptBinWidthInterference("INT_PT_WIDTH",0),
_protonEnergy ("PROTON_ENERGY",0, NOT_REQUIRED),
_electronEnergy ("ELECTRON_ENERGY",0, NOT_REQUIRED),
_minGammaEnergy ("MIN_GAMMA_ENERGY",6.0, NOT_REQUIRED),
_maxGammaEnergy ("MAX_GAMMA_ENERGY",600000.0, NOT_REQUIRED),
_minGammaQ2 ("MIN_GAMMA_Q2",0,NOT_REQUIRED),
_maxGammaQ2 ("MAX_GAMMA_Q2",0,NOT_REQUIRED),
_nmbGammaQ2Bins ("INT_GAMMA_Q2_BINS",400,NOT_REQUIRED),
_pythiaParams ("PYTHIA_PARAMS","", NOT_REQUIRED),
_pythiaFullEventRecord ("PYTHIA_FULL_EVENTRECORD",false, NOT_REQUIRED),
+ _hepmc3FullEventRecord ("HEPMC3_FULL_EVENTRECORD",false, NOT_REQUIRED),
_xsecCalcMethod ("XSEC_METHOD",0, NOT_REQUIRED),
_axionMass ("AXION_MASS",50, NOT_REQUIRED), // AXION HACK
_bslopeDefinition ("BSLOPE_DEFINITION",0, NOT_REQUIRED),
_bslopeValue ("BSLOPE_VALUE",4.0,NOT_REQUIRED),
_impulseVM ("SELECT_IMPULSE_VM",0,NOT_REQUIRED),
_quantumGlauber ("QUANTUM_GLAUBER",0,NOT_REQUIRED)
{
// All parameters must be initialised in initialisation list!
// If not: error: 'parameter<T, validate>::parameter() [with T = unsigned int, bool validate = true]' is private
// or similar
_ip.addParameter(_baseFileName);
_ip.addParameter(_beam1Z);
_ip.addParameter(_beam2Z);
_ip.addParameter(_beam1A);
_ip.addParameter(_beam2A);
_ip.addParameter(_beam1LorentzGamma);
_ip.addParameter(_beam2LorentzGamma);
_ip.addParameter(_maxW);
_ip.addParameter(_minW);
_ip.addParameter(_nmbWBins);
_ip.addParameter(_maxRapidity);
_ip.addParameter(_nmbRapidityBins);
//
_ip.addParameter(_nmbEnergyBins);
_ip.addParameter(_ptCutEnabled);
_ip.addParameter(_ptCutMin);
_ip.addParameter(_ptCutMax);
_ip.addParameter(_etaCutEnabled);
_ip.addParameter(_etaCutMax);
_ip.addParameter(_etaCutMin);
_ip.addParameter(_productionMode);
_ip.addParameter(_nmbEventsTot);
_ip.addParameter(_prodParticleId);
_ip.addParameter(_randomSeed);
_ip.addParameter(_beamBreakupMode);
_ip.addParameter(_interferenceEnabled);
_ip.addParameter(_interferenceStrength);
_ip.addParameter(_maxPtInterference);
_ip.addParameter(_nmbPtBinsInterference);
_ip.addParameter(_minGammaEnergy);
_ip.addParameter(_maxGammaEnergy);
_ip.addParameter(_minGammaQ2);
_ip.addParameter(_maxGammaQ2);
_ip.addParameter(_nmbGammaQ2Bins);
_ip.addParameter(_pythiaParams);
_ip.addParameter(_pythiaFullEventRecord);
+ _ip.addParameter(_hepmc3FullEventRecord);
_ip.addParameter(_xsecCalcMethod);
_ip.addParameter(_axionMass); // AXION HACK
_ip.addParameter(_bslopeDefinition);
_ip.addParameter(_bslopeValue);
_ip.addParameter(_impulseVM);
_ip.addParameter(_quantumGlauber);
}
//______________________________________________________________________________
inputParameters::~inputParameters()
{ }
//______________________________________________________________________________
bool
inputParameters::configureFromFile(const std::string &_configFileName)
{
int nParameters = _ip.parseFile(_configFileName);
if(nParameters == -1)
{
printWarn << "could not open file '" << _configFileName << "'" << endl;
return false;
}
if(_ip.validateParameters(cerr))
printInfo << "successfully read input parameters from '" << _configFileName << "'" << endl;
else {
printWarn << "problems reading input parameters from '" << _configFileName << "'" << endl
<< *this;
return false;
}
return true;
}
bool inputParameters::init()
{
// Calculate beam gamma in CMS frame
double rap1 = acosh(beam1LorentzGamma());
double rap2 = -acosh(beam2LorentzGamma());
_beamLorentzGamma = cosh((rap1-rap2)/2);
_targetLorentzGamma = cosh(rap1-rap2);
if( beam2A() == 1) //proton case 0.87 fm = 4.4 GeV^{-1}
_targetR = 4.4;
else
_targetR = 6.1 * pow(beam2A(), 1./3.);
//_targetR = 4.4/2.;
_fixedQ2Range = false;
std::cout << "Rapidity beam 1: " << rap1 << ", rapidity beam 2: " << rap2 << ", rapidity CMS system: " << (rap1+rap2)/2 << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
std::cout << "Rapidity beam 1 in beam 2 frame: " << rap1-rap2 << ", beam 1 gamma in beam 2 frame: " << _targetLorentzGamma<< std::endl;
_ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
_protonEnergy = _beamLorentzGamma * protonMass;
// Electron energy in the target frame of reference
_electronEnergy = _targetLorentzGamma * starlightConstants::mel;
if((beam1Z()==1) && (beam1A()==2)){
if((beam2Z()==1) && (beam2A()==2)){
printWarn << "deuteron-deuteron collisions are not supported" << endl;
return false;}
printWarn << "deuterium must be listed as the second nucleus" << endl;
return false;}
if( ((beam1Z()==1) && (beam1A()==3)) || ((beam2Z()==1) && (beam2A()==3)) ){
printWarn << "tritium is not currently supported" << endl;
return false;}
// check that rho production uses wide resonance option
if(_prodParticleId.value()==113 && _productionMode.value()==2){
printWarn << endl<< "For rho meson production, you should choose the wide resonance option (production mode = 3)" << endl;
return false;}
// define interaction type
switch (productionMode()) {
case 1:
_interactionType = PHOTONPHOTON;
break;
case 2:
_interactionType = PHOTONPOMERONNARROW;
break;
case 3:
_interactionType = PHOTONPOMERONWIDE;
break;
case 4:
_interactionType = PHOTONPOMERONINCOHERENT;
break;
case 5:
_interactionType = PHOTONUCLEARSINGLE;
break;
case 6:
_interactionType = PHOTONUCLEARDOUBLE;
break;
case 7:
_interactionType = PHOTONUCLEARSINGLEPA;
break;
case 8:
_interactionType = PHOTONUCLEARSINGLEPAPY;
break;
// case 9:
// _interactionType = PHOTONPHOTONKNIEHL;
// break;
// case 10:
// _interactionType = PHOTONPHOTONKNIEHLMODIFIED;
// break;
case 12:
_interactionType = E_PHOTONPOMERONNARROW;
break;
case 13:
_interactionType = E_PHOTONPOMERONWIDE;
break;
default:
printWarn << "unknown production mode '" << _productionMode << "'" << endl;
return false;
}
if( (_productionMode.value() ==4) && (_interferenceEnabled.value())) {
printWarn << " cannot enable interference for incoherent production " << endl;
return false;
}
double mass = 0;
double width = 0;
double defaultMinW = 0; // default for _minW, unless it is defined later [GeV/c^2]
double defaultMaxW = 0; // default for _maxW, unless it is defined later [GeV/c^2]
switch (prodParticleId()) {
case 11: // e+e- pair
_particleType = ELECTRON;
_decayType = LEPTONPAIR;
mass = starlightConstants::mel;
defaultMinW = 2*mass; // default is 0.01; up to 0.15 is safe for Summer 2000 triggering for e+e- pairs
defaultMaxW = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 13: // mu+mu- pair
_particleType = MUON;
_decayType = LEPTONPAIR;
defaultMinW = 2 * muonMass;
defaultMaxW = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 15: // tau+tau- pair
_particleType = TAUON;
_decayType = LEPTONPAIR;
defaultMinW = 2 * tauMass;
defaultMaxW = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 10015: // tau+tau- pair
_particleType = TAUONDECAY;
_decayType = LEPTONPAIR;
defaultMinW = 2 * tauMass;
defaultMaxW = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
// case 24: // W+W- pair
// _particleType = W;
// _decayType = WW;
// defaultMinW = 2 * muonMass;
// break;
case 115: // a_2(1320)
_particleType = A2;
_decayType = SINGLEMESON;
mass = starlightConstants::a2Mass;
width = starlightConstants::a2Width;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 221: // eta
_particleType = ETA;
_decayType = SINGLEMESON;
mass = starlightConstants::etaMass;
width = starlightConstants::etaWidth;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 225: // f_2(1270)
_particleType = F2;
_decayType = SINGLEMESON;
mass = starlightConstants::f2Mass;
width = starlightConstants::f2Width;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = starlightConstants::f2BrPiPi;
break;
case 331: // eta'(958)
_particleType = ETAPRIME;
_decayType = SINGLEMESON;
mass = starlightConstants::etaPrimeMass;
width = starlightConstants::etaPrimeWidth;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 335: // f_2'(1525)
_particleType = F2PRIME;
_decayType = SINGLEMESON;
mass = starlightConstants::f2PrimeMass;
width = starlightConstants::f2PrimeWidth;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = starlightConstants::f2PrimeBrKK;
break;
case 441: // eta_c(1s)
_particleType = ETAC;
_decayType = SINGLEMESON;
mass = starlightConstants::etaCMass;
width = starlightConstants::etaCWidth;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 9010221: // f_0(980), was orginally called 10221? updated to standard number
_particleType = F0;
_decayType = SINGLEMESON;
mass = starlightConstants::f0Mass;
width = starlightConstants::f0Width;
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = starlightConstants::f0BrPiPi;
break;
case 33: // Z"/Z0 This is the rho^0 rho^0 final state SRK
_particleType = ZOVERZ03;
_decayType = SINGLEMESON;
defaultMinW = 4*pionChargedMass;
defaultMaxW = 1.6; // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 88: // axion// AXION HACK, till break statement
_particleType = AXION;
_decayType = SINGLEMESON;
mass = _axionMass.value();
width = 1/(64*starlightConstants::pi)*mass*mass*mass/(1000*1000);//Fix Lambda=1000 GeV, rescaling is trivial.
defaultMinW = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
defaultMaxW = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
break; // AXION HACK, end
// case 25: // Higgs
// _particleType = HIGGS;
// _decayType = SINGLEMESON;
// break;
case 113: // rho(770)
_particleType = RHO;
_decayType = WIDEVMDEFAULT;
mass = starlightConstants::rho0Mass;
width = starlightConstants::rho0Width;
defaultMinW = 2 * pionChargedMass;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::rho0BrPiPi;
break;
case 913: // rho(770) with direct pi+pi- decay, interference given by ZEUS data
_particleType = RHOZEUS;
_decayType = WIDEVMDEFAULT;
mass = starlightConstants::rho0Mass;
width = starlightConstants::rho0Width;
defaultMinW = 2 * pionChargedMass;
defaultMaxW = mass + 5 * width; // use the same 1.5GeV max mass as ZEUS
_inputBranchingRatio = starlightConstants::rho0BrPiPi;
break;
case 999: // pi+pi-pi+pi- phase space decay
_particleType = FOURPRONG;
_decayType = WIDEVMDEFAULT;
mass = starlightConstants::rho0PrimeMass;
width = starlightConstants::rho0PrimeWidth;
defaultMinW = 4 * pionChargedMass;
defaultMaxW = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
_inputBranchingRatio = 1.0;
break;
case 223: // omega(782)
_particleType = OMEGA;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::OmegaMass;
width = starlightConstants::OmegaWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::OmegaBrPiPi;
break;
case 333: // phi(1020)
_particleType = PHI;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::PhiMass;
width = starlightConstants::PhiWidth;
defaultMinW = 2 * kaonChargedMass;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::PhiBrKK;
break;
case 443: // J/psi
_particleType = JPSI;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::JpsiMass;
width = starlightConstants::JpsiWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (starlightConstants::JpsiBree + starlightConstants::JpsiBrmumu)/2.;
break;
case 443011: // J/psi
_particleType = JPSI_ee;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::JpsiMass;
width = starlightConstants::JpsiWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::JpsiBree;
break;
case 443013: // J/psi
cout<<"In inputParameters setting J/psi mass!"<<endl;
_particleType = JPSI_mumu;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::JpsiMass;
width = starlightConstants::JpsiWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::JpsiBrmumu;
break;
case 444: // psi(2S)
_particleType = JPSI2S;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Psi2SMass;
width = starlightConstants::Psi2SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (starlightConstants::Psi2SBree + starlightConstants::Psi2SBrmumu)/2.;
break;
case 444011: // psi(2S)
_particleType = JPSI2S_ee;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Psi2SMass;
width = starlightConstants::Psi2SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Psi2SBree;
break;
case 444013: // psi(2S)
_particleType = JPSI2S_mumu;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Psi2SMass;
width = starlightConstants::Psi2SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Psi2SBrmumu;
break;
case 553: // Upsilon(1S)
_particleType = UPSILON;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon1SMass;
width = starlightConstants::Upsilon1SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (starlightConstants::Upsilon1SBree + starlightConstants::Upsilon1SBrmumu)/2.;
break;
case 553011: // Upsilon
_particleType = UPSILON_ee;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon1SMass;
width = starlightConstants::Upsilon1SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Upsilon1SBree;
break;
case 553013: // Upsilon
_particleType = UPSILON_mumu;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon1SMass;
width = starlightConstants::Upsilon1SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Upsilon1SBrmumu;
break;
case 554: // Upsilon(2S)
_particleType = UPSILON2S;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon2SMass;
width = starlightConstants::Upsilon2SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (starlightConstants::Upsilon2SBree + starlightConstants::Upsilon2SBrmumu)/2.;
break;
case 554011: // Upsilon(2S)
_particleType = UPSILON2S_ee;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon2SMass;
width = starlightConstants::Upsilon2SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Upsilon2SBree;
break;
case 554013: // Upsilon(2S)
_particleType = UPSILON2S_mumu;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon2SMass;
width = starlightConstants::Upsilon2SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Upsilon2SBrmumu;
break;
case 555: // Upsilon(3S)
_particleType = UPSILON3S;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon3SMass;
width = starlightConstants::Upsilon3SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (starlightConstants::Upsilon3SBree + starlightConstants::Upsilon3SBrmumu)/2.;
break;
case 555011: // Upsilon(3S)
_particleType = UPSILON3S_ee;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon3SMass;
width = starlightConstants::Upsilon3SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Upsilon3SBree;
break;
case 555013: // Upsilon(3S)
_particleType = UPSILON3S_mumu;
_decayType = NARROWVMDEFAULT;
mass = starlightConstants::Upsilon3SMass;
width = starlightConstants::Upsilon3SWidth;
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = starlightConstants::Upsilon3SBrmumu;
break;
default:
printWarn << "unknown particle ID " << _prodParticleId << endl;
return false;
} // _prodParticleId
if (_minW.value() == -1)
_minW = defaultMinW;
if (_maxW.value() == -1)
_maxW = defaultMaxW;
if ( _maxW.value() <= _minW.value() ) {
printWarn << "maxW must be greater than minW" << endl;
return false;
}
// Sanity check on Q2 range in case it is specified by user
if( _minGammaQ2.value() != 0 || _maxGammaQ2.value() != 0){
if( _minGammaQ2.value() <0 || _maxGammaQ2.value() <=_minGammaQ2.value() )
printWarn << "Input values for min and max Q2 are inconsistent: "<<_minGammaQ2.value()<<","<<_maxGammaQ2.value()
<<". Continuing with default "<<endl;
else
_fixedQ2Range = true;
}
//Photon energy limits in C.M.S frame, used for some basic safety chekcs
_cmsMinPhotonEnergy = 0.5*(((mass+protonMass)*(mass+protonMass)-
protonMass*protonMass)/(_protonEnergy.value()+sqrt(_protonEnergy.value()*_protonEnergy.value()-protonMass*protonMass)));
_cmsMaxPhotonEnergy = 0.5*mass*exp(9);
// Photon limits in target frame: this is where the photon flux is well described
_targetMaxPhotonEnergy = (_targetLorentzGamma - 10. ) *starlightConstants::mel;
_targetMinPhotonEnergy = mass;
printInfo << "using the following " << *this;
return true;
}
//______________________________________________________________________________
ostream&
inputParameters::print(ostream& out) const
{
out << "starlight parameters:" << endl
<< " base file name ...................... '" << _baseFileName.value() << "'" << endl
<< " beam 1 atomic number ................... " << _beam1Z.value() << endl
<< " beam 1 atomic mass number .............. " << _beam1A.value() << endl
<< " beam 2 atomic number ................... " << _beam2Z.value() << endl
<< " beam 2 atomic mass number .............. " << _beam2A.value() << endl
<< " Lorentz gamma of beams in CM frame ..... " << _beamLorentzGamma << endl
<< " mass W of produced hadronic system ..... " << _minW.value() << " < W < " << _maxW.value() << " GeV/c^2" << endl
<< " # of W bins ............................ " << _nmbWBins.value() << endl
<< " maximum absolute value for rapidity .... " << _maxRapidity.value() << endl
<< " # of rapidity bins ..................... " << _nmbRapidityBins.value() << endl
<< " # of Egamma bins ....................... " << _nmbEnergyBins.value() << endl
<< " cut in pT............................... " << yesNo(_ptCutEnabled.value()) << endl;
if (_ptCutEnabled.value()) {
out << " minumum pT.......................... " << _ptCutMin.value() << " GeV/c" << endl
<< " maximum pT.......................... " << _ptCutMax.value() << " GeV/c" << endl;}
out << " cut in eta.............................. " << yesNo(_etaCutEnabled.value()) << endl;
if (_etaCutEnabled.value()) {
out << " minumum eta......................... " << _etaCutMin.value() << endl
<< " maximum eta......................... " << _etaCutMax.value() << endl;}
out << " production mode ........................ " << _productionMode.value() << endl
<< " number of events to generate ........... " << _nmbEventsTot.value() << endl
<< " PDG ID of produced particle ............ " << _prodParticleId.value() << endl
<< " seed for random generator .............. " << _randomSeed.value() << endl
<< " breakup mode for beam particles ........ " << _beamBreakupMode.value() << endl
<< " interference enabled ................... " << yesNo(_interferenceEnabled.value()) << endl;
if (_interferenceEnabled.value()) {
out << " interference strength .................. " << _interferenceStrength.value() << endl
<< " maximum p_T for interference calc. ..... " << _maxPtInterference.value() << " GeV/c" << endl
<< " # of p_T bins for interference calc. ... " << _nmbPtBinsInterference.value() << endl;}
if (_productionMode.value()!=1) {
if (_productionMode.value()==4) {
out << " coherent scattering off nucleus ........ no" << endl;}
else {
out << " coherent scattering off nucleus ........ yes" << endl;
}
}
if (_fixedQ2Range == true ){
out << " fixed photon Q2 range .................. " << _minGammaQ2.value() << " < Q2 < "
<<_maxGammaQ2.value() << " GeV/c^2 "<<endl;
}
out<< " Q2_BINS" <<nmbGammaQ2Bins () <<endl;
out <<" Quantum Glauber parameter............... "<<_quantumGlauber.value()<<endl;
out <<" Impulse VM parameter.................... "<<_impulseVM.value()<<endl;
return out;
}
//______________________________________________________________________________
ostream&
inputParameters::write(ostream& out) const
{
out << "baseFileName" << baseFileName () <<endl
<< "BEAM_1_Z" << beam1Z () <<endl
<< "BEAM_2_Z" << beam1A () <<endl
<< "BEAM_1_A" << beam2Z () <<endl
<< "BEAM_2_A" << beam2A () <<endl
<< "BEAM_GAMMA" << beamLorentzGamma () <<endl
<< "W_MAX" << maxW () <<endl
<< "W_MIN" << minW () <<endl
<< "W_N_BINS" << nmbWBins () <<endl
<< "RAP_MAX" << maxRapidity () <<endl
<< "RAP_N_BINS" << nmbRapidityBins () <<endl
<< "CUT_PT" << ptCutEnabled () <<endl
<< "PT_MIN" << ptCutMin () <<endl
<< "PT_MAX" << ptCutMax () <<endl
<< "CUT_ETA" << etaCutEnabled () <<endl
<< "ETA_MIN" << etaCutMin () <<endl
<< "ETA_MAX" << etaCutMax () <<endl
<< "PROD_MODE" << productionMode () <<endl
<< "N_EVENTS" << nmbEvents () <<endl
<< "PROD_PID" << prodParticleId () <<endl
<< "RND_SEED" << randomSeed () <<endl
<< "BREAKUP_MODE" << beamBreakupMode () <<endl
<< "INTERFERENCE" << interferenceEnabled () <<endl
<< "IF_STRENGTH" << interferenceStrength () <<endl
<< "INT_PT_MAX" << maxPtInterference () <<endl
<< "INT_PT_N_BINS" << nmbPtBinsInterference() <<endl;
out<< "USING FIXED RANGE"<<_fixedQ2Range<<endl;
if( _fixedQ2Range == true ){
out << "Q2_MIN" << minGammaQ2 () <<endl
<< "Q2_MAX" << maxGammaQ2 () <<endl;
}
out<< " Q2_BINS" <<nmbGammaQ2Bins () <<endl;
return out;
}
std::string
inputParameters::parameterValueKey() const
{
return parameterbase::_parameters.validationKey();
}
Index: trunk/cmake_modules/FindHepMC3.cmake
===================================================================
--- trunk/cmake_modules/FindHepMC3.cmake (revision 0)
+++ trunk/cmake_modules/FindHepMC3.cmake (revision 15)
@@ -0,0 +1,27 @@
+
+
+FIND_PATH(HEPMC3_INCLUDE_DIR HepMC3/GenEvent.h HINTS /usr/local/include $ENV{HEPMC3DIR}/include)
+#FIND_LIBRARY(HEPMC3_LIBRARY NAMES HepMC3 PATHS $ENV{HEPMC3_ROOT_DIR}/lib)
+FIND_LIBRARY(HEPMC3_LIBRARY NAMES HepMC3 PATHS /usr/lib /usr/lib/HepMC3 /usr/local/lib $ENV{HEPMC3DIR}/lib)
+
+IF (HEPMC3_INCLUDE_DIR AND HEPMC3_LIBRARY)
+ SET(HEPMC3_FOUND TRUE)
+ENDIF (HEPMC3_INCLUDE_DIR AND HEPMC3_LIBRARY)
+
+
+IF (HEPMC3_FOUND)
+ IF (NOT HEPMC3_FIND_QUIETLY)
+ MESSAGE(STATUS "Found HepMC3: ${HEPMC3_LIBRARY}")
+ MESSAGE(STATUS "Found HepMC3 include: ${HEPMC3_INCLUDE_DIR}")
+ ENDIF (NOT HEPMC3_FIND_QUIETLY)
+ELSE (HEPMC3_FOUND)
+ IF (HEPMC3_FIND_REQUIRED)
+ MESSAGE(FATAL_ERROR "Could not find HepMC3. We search first in the normal library paths, then in $HEPMC3IR")
+ ELSE(HEPMC3_FIND_REQUIRED)
+ IF(NOT HEPMC3_FIND_QUIETLY)
+ MESSAGE(STATUS "Could not find HepMC3. We search first in the normal library paths, then in $HEPMC3DIR")
+ ENDIF(NOT HEPMC3_FIND_QUIETLY)
+ ENDIF (HEPMC3_FIND_REQUIRED)
+
+ENDIF (HEPMC3_FOUND)
+