Page MenuHomeHEPForge
Files F7878816

No OneTemporary
Actions

View Options

Index: trunk/include/inputParameters.h
===================================================================
--- trunk/include/inputParameters.h (revision 308)
+++ trunk/include/inputParameters.h (revision 309)
@@ -1,539 +1,543 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: $: revision of last commit
// $Author:: $: author of last commit
// $Date:: $: 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>
#include <cstdlib>
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(); }
unsigned 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
unsigned 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 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 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
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]
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
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 printVM () const { return _printVM.value(); } ///< returns the printVM value
int impulseVM () const { return _impulseVM.value(); } ///< returns the impulseVM value
int quantumGlauber () const { return _quantumGlauber.value(); } ///< returns the quantum glauber value
double bmin () const { return _bmin.value(); } // returns the minimum impact parameter for BREAKUP_MODE==6
double bmax () const { return _bmax.value(); } // returns the maximum impact parameter for BREAKUP_MODE==6
bool outputHeader () const { return _outputHeader.value(); } // returns whether the output should have parameters at the start
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 inputBranchingRatio () const { return _inputBranchingRatio; }
double deuteronSlopePar () const {return _deuteronSlopePar .value();}
double protonMass () const {return _protonMass .value();}
double pionChargedMass () const {return _pionChargedMass .value();}
double pionNeutralMass () const {return _pionNeutralMass .value();}
double kaonChargedMass () const {return _kaonChargedMass .value();}
double mel () const {return _mel .value();}
double muonMass () const {return _muonMass .value();}
double tauMass () const {return _tauMass .value();}
double f0Mass () const {return _f0Mass .value();}
double f0Width () const {return _f0Width .value();}
double f0BrPiPi () const {return _f0BrPiPi .value();}
double etaMass () const {return _etaMass .value();}
double etaWidth () const {return _etaWidth .value();}
double etaPrimeMass () const {return _etaPrimeMass .value();}
double etaPrimeWidth () const {return _etaPrimeWidth .value();}
double etaCMass () const {return _etaCMass .value();}
double etaCWidth () const {return _etaCWidth .value();}
double f2Mass () const {return _f2Mass .value();}
double f2Width () const {return _f2Width .value();}
double f2BrPiPi () const {return _f2BrPiPi .value();}
double a2Mass () const {return _a2Mass .value();}
double a2Width () const {return _a2Width .value();}
double f2PrimeMass () const {return _f2PrimeMass .value();}
double f2PrimeWidth () const {return _f2PrimeWidth .value();}
double f2PrimeBrKK () const {return _f2PrimeBrKK .value();}
double zoverz03Mass () const {return _zoverz03Mass .value();}
double f0PartialggWidth () const {return _f0PartialggWidth .value();}
double etaPartialggWidth () const {return _etaPartialggWidth .value();}
double etaPrimePartialggWidth() const {return _etaPrimePartialggWidth.value();}
double etaCPartialggWidth () const {return _etaCPartialggWidth .value();}
double f2PartialggWidth () const {return _f2PartialggWidth .value();}
double a2PartialggWidth () const {return _a2PartialggWidth .value();}
double f2PrimePartialggWidth () const {return _f2PrimePartialggWidth .value();}
double zoverz03PartialggWidth() const {return _zoverz03PartialggWidth.value();}
double f0Spin () const {return _f0Spin .value();}
double etaSpin () const {return _etaSpin .value();}
double etaPrimeSpin () const {return _etaPrimeSpin .value();}
double etaCSpin () const {return _etaCSpin .value();}
double f2Spin () const {return _f2Spin .value();}
double a2Spin () const {return _a2Spin .value();}
double f2PrimeSpin () const {return _f2PrimeSpin .value();}
double zoverz03Spin () const {return _zoverz03Spin .value();}
double axionSpin () const {return _axionSpin .value();}
double rho0Mass () const {return _rho0Mass .value();}
double rho0Width () const {return _rho0Width .value();}
double rho0BrPiPi () const {return _rho0BrPiPi .value();}
+ double rho0Bree () const {return _rho0Bree .value();}
+ double rho0Brmumu () const {return _rho0Brmumu .value();}
double rho0PrimeMass () const {return _rho0PrimeMass .value();}
double rho0PrimeWidth () const {return _rho0PrimeWidth .value();}
double rho0PrimeBrPiPi () const {return _rho0PrimeBrPiPi .value();}
double OmegaMass () const {return _OmegaMass .value();}
double OmegaWidth () const {return _OmegaWidth .value();}
double OmegaBrPiPi () const {return _OmegaBrPiPi .value();}
double PhiMass () const {return _PhiMass .value();}
double PhiWidth () const {return _PhiWidth .value();}
double PhiBrKK () const {return _PhiBrKK .value();}
double JpsiMass () const {return _JpsiMass .value();}
double JpsiWidth () const {return _JpsiWidth .value();}
double JpsiBree () const {return _JpsiBree .value();}
double JpsiBrmumu () const {return _JpsiBrmumu .value();}
double JpsiBrppbar () const {return _JpsiBrppbar .value();}
double Psi2SMass () const {return _Psi2SMass .value();}
double Psi2SWidth () const {return _Psi2SWidth .value();}
double Psi2SBree () const {return _Psi2SBree .value();}
double Psi2SBrmumu () const {return _Psi2SBrmumu .value();}
double Upsilon1SMass () const {return _Upsilon1SMass .value();}
double Upsilon1SWidth () const {return _Upsilon1SWidth .value();}
double Upsilon1SBree () const {return _Upsilon1SBree .value();}
double Upsilon1SBrmumu () const {return _Upsilon1SBrmumu .value();}
double Upsilon2SMass () const {return _Upsilon2SMass .value();}
double Upsilon2SWidth () const {return _Upsilon2SWidth .value();}
double Upsilon2SBree () const {return _Upsilon2SBree .value();}
double Upsilon2SBrmumu () const {return _Upsilon2SBrmumu .value();}
double Upsilon3SMass () const {return _Upsilon3SMass .value();}
double Upsilon3SWidth () const {return _Upsilon3SWidth .value();}
double Upsilon3SBree () const {return _Upsilon3SBree .value();}
double Upsilon3SBrmumu () const {return _Upsilon3SBrmumu .value();}
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 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 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 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 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 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 setprintVM (int v) { _printVM = v; } ///< sets the value of _printVM
void setimpulseVM (int v) { _impulseVM = v; } ///< sets the value of _impulseVM
void setquantumGlauber (int v) { _quantumGlauber = v; } ///< sets the value of quantum_glauber
void setbmin (double v) { _bmin=v; } ///< sets the minimum impact parameter (for BREAKUP_MODE==6
void setbmax (double v) { _bmax=v; } ///< sets the minimum impact parameter (for BREAKUP_MODE==6
void setoutputHeader (bool v) { _outputHeader = v; } ///< sets whether the output should have parameters at the start
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
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<unsigned int,VALIDITY_CHECK> _beam1Z; ///< atomic number of beam particle 1
parameter<unsigned int,VALIDITY_CHECK> _beam1A; ///< atomic mass number of beam particle 1
parameter<unsigned 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<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)
/// 6 = set impact parameter range via bmax and bmin
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> _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<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<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> _printVM; ///< Optional parameter to set printing options for VM cross section
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
parameter<double, VALIDITY_CHECK> _bmin; ///< Optional parameter minimum impact parameter for b-range calculation
parameter<double, VALIDITY_CHECK> _bmax; /// < Optional parameter maximum impact parameter for b-range calculation
parameter<bool, NO_VALIDITY_CHECK> _outputHeader; /// < Optional parameter puts parameter information at the top of the output file
parameter<double, VALIDITY_CHECK> _deuteronSlopePar ; ///< deuteron slope parameter (effective temperature) [(GeV/c)^-2]
parameter<double, VALIDITY_CHECK> _protonMass ; ///< mass of the proton [GeV/c^2]
parameter<double, VALIDITY_CHECK> _pionChargedMass ; ///< mass of the pi^+/- [GeV/c^2]
parameter<double, VALIDITY_CHECK> _pionNeutralMass ; ///< mass of the pi^0 [GeV/c^2]
parameter<double, VALIDITY_CHECK> _kaonChargedMass ; ///< mass of the K^+/- [GeV/c^2]
parameter<double, VALIDITY_CHECK> _mel ; ///< mass of the e^+/- [GeV/c^2]
parameter<double, VALIDITY_CHECK> _muonMass ; ///< mass of the mu^+/- [GeV/c^2]
parameter<double, VALIDITY_CHECK> _tauMass ; ///< mass of the tau^+/- [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f0Mass ; ///< mass of the f_0(980) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f0Width ; ///< width of the f_0(980) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f0BrPiPi ; ///< branching ratio f_0(980) -> pi^+ pi^- and pi^0 pi^0
parameter<double, VALIDITY_CHECK> _etaMass ; ///< mass of the eta [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaWidth ; ///< width of the eta [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaPrimeMass ; ///< mass of the eta' [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaPrimeWidth ; ///< width of the eta' [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaCMass ; ///< mass of the eta_c [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaCWidth ; ///< width of the eta_c [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2Mass ; ///< mass of the f_2(1270) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2Width ; ///< width of the f_2(1270) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2BrPiPi ; ///< branching ratio f_2(1270) -> pi^+ pi^-
parameter<double, VALIDITY_CHECK> _a2Mass ; ///< mass of the a_2(1320) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _a2Width ; ///< width of the a_2(1320) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2PrimeMass ; ///< mass of the f'_2(1525) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2PrimeWidth ; ///< width of the f'_2(1525) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2PrimeBrKK ; ///< branching ratio f'_2(1525) -> K^+ K^- and K^0 K^0bar
parameter<double, VALIDITY_CHECK> _zoverz03Mass ; ///< mass of four-quark resonance (rho^0 pair production) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f0PartialggWidth ; ///< partial width f_0(980) -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaPartialggWidth ; ///< partial width eta -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaPrimePartialggWidth; ///< partial width eta' -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _etaCPartialggWidth ; ///< partial width eta_c -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2PartialggWidth ; ///< partial width f_2(1270) -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _a2PartialggWidth ; ///< partial width a_2(1320) -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f2PrimePartialggWidth ; ///< partial width f'_2(1525) -> g g [GeV/c^2]
parameter<double, VALIDITY_CHECK> _zoverz03PartialggWidth; ///< partial width four-quark resonance -> g g (rho^0 pair production) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _f0Spin ; ///< spin of the f_0(980)
parameter<double, VALIDITY_CHECK> _etaSpin ; ///< spin of the eta
parameter<double, VALIDITY_CHECK> _etaPrimeSpin ; ///< spin of the eta'
parameter<double, VALIDITY_CHECK> _etaCSpin ; ///< spin of the eta_c
parameter<double, VALIDITY_CHECK> _f2Spin ; ///< spin of the f_2(1270)
parameter<double, VALIDITY_CHECK> _a2Spin ; ///< spin of the a_2(1320)
parameter<double, VALIDITY_CHECK> _f2PrimeSpin ; ///< spin of the f'_2(1525)
parameter<double, VALIDITY_CHECK> _zoverz03Spin ; ///< spin of the four-quark resonance -> g g (rho^0 pair production)
parameter<double, VALIDITY_CHECK> _axionSpin ; ///< spin of the axion
parameter<double, VALIDITY_CHECK> _rho0Mass ; ///< mass of the rho^0 [GeV/c^2]
parameter<double, VALIDITY_CHECK> _rho0Width ; ///< width of the rho^0 [GeV/c^2]
parameter<double, VALIDITY_CHECK> _rho0BrPiPi ; ///< branching ratio rho^0 -> pi^+ pi^-
+ parameter<double, VALIDITY_CHECK> _rho0Bree ; ///< branching ratio rho^0 -> e^+ e^-
+ parameter<double, VALIDITY_CHECK> _rho0Brmumu ; ///< branching ratio rho^0 -> mu^+ mu^-
parameter<double, VALIDITY_CHECK> _rho0PrimeMass ; ///< mass of the rho'^0 (4 pi^+/- final state) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _rho0PrimeWidth ; ///< width of the rho'^0 (4 pi^+/- final state) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _rho0PrimeBrPiPi ; ///< branching ratio rho'^0 -> pi^+ pi^-
parameter<double, VALIDITY_CHECK> _OmegaMass ; ///< mass of the omega [GeV/c^2]
parameter<double, VALIDITY_CHECK> _OmegaWidth ; ///< width of the omega [GeV/c^2]
parameter<double, VALIDITY_CHECK> _OmegaBrPiPi ; ///< branching ratio omega -> pi^+ pi^-
parameter<double, VALIDITY_CHECK> _PhiMass ; ///< mass of the phi [GeV/c^2]
parameter<double, VALIDITY_CHECK> _PhiWidth ; ///< width of the phi [GeV/c^2]
parameter<double, VALIDITY_CHECK> _PhiBrKK ; ///< branching ratio phi -> K^+ K^-
parameter<double, VALIDITY_CHECK> _JpsiMass ; ///< mass of the J/psi [GeV/c^2]
parameter<double, VALIDITY_CHECK> _JpsiWidth ; ///< width of the J/psi [GeV/c^2]
parameter<double, VALIDITY_CHECK> _JpsiBree ; ///< branching ratio J/psi -> e^+ e^-
parameter<double, VALIDITY_CHECK> _JpsiBrmumu ; ///< branching ratio J/psi -> mu^+ mu^-
parameter<double, VALIDITY_CHECK> _JpsiBrppbar ; ///< branching ratio J/psi -> p pbar
parameter<double, VALIDITY_CHECK> _Psi2SMass ; ///< mass of the psi(2S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Psi2SWidth ; ///< width of the psi(2S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Psi2SBree ; ///< branching ratio psi(2S) -> e^+ e^-
parameter<double, VALIDITY_CHECK> _Psi2SBrmumu ; ///< branching ratio psi(2S) -> mu^+ mu^-
parameter<double, VALIDITY_CHECK> _Upsilon1SMass ; ///< mass of the Upsilon(1S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Upsilon1SWidth ; ///< width of the Upsilon(1S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Upsilon1SBree ; ///< branching ratio Upsilon(1S) -> e^+ e^-
parameter<double, VALIDITY_CHECK> _Upsilon1SBrmumu ; ///< branching ratio Upsilon(1S) -> mu^+ mu^-
parameter<double, VALIDITY_CHECK> _Upsilon2SMass ; ///< mass of the Upsilon(2S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Upsilon2SWidth ; ///< width of the Upsilon(2S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Upsilon2SBree ; ///< branching ratio Upsilon(2S) -> e^+ e^-
parameter<double, VALIDITY_CHECK> _Upsilon2SBrmumu ; ///< branching ratio Upsilon(2S) -> mu^+ mu^-
parameter<double, VALIDITY_CHECK> _Upsilon3SMass ; ///< mass of the Upsilon(3S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Upsilon3SWidth ; ///< width of the Upsilon(3S) [GeV/c^2]
parameter<double, VALIDITY_CHECK> _Upsilon3SBree ; ///< branching ratio Upsilon(3S) -> e^+ e^-
parameter<double, VALIDITY_CHECK> _Upsilon3SBrmumu ; ///< branching ratio Upsilon(3S) -> mu^+ mu^-
starlightConstants::particleTypeEnum _particleType;
starlightConstants::decayTypeEnum _decayType;
starlightConstants::interactionTypeEnum _interactionType;
double _beamLorentzGamma; ///< Lorentz gamma factor of the beams in CMS frame, not an input parameter
double _inputBranchingRatio; ///< Branching ratio defined for each channel
inputParser _ip;
};
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/starlightconstants.h
===================================================================
--- trunk/include/starlightconstants.h (revision 308)
+++ trunk/include/starlightconstants.h (revision 309)
@@ -1,155 +1,157 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: $: revision of last commit
// $Author:: $: author of last commit
// $Date:: $: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#ifndef STARLIGHTCONSTANTS_H_INCLUDE
#define STARLIGHTCONSTANTS_H_INCLUDE
/*
* Constants are set here
*/
namespace starlightConstants
{
// constants
static const double hbarc = 0.1973269718;
static const double hbarcmev = hbarc*1000.;
static const double pi = 3.141592654;
static const double twoPi = 2 * pi;
static const double alpha = 1/137.035999074;
enum particleTypeEnum {
UNKNOWN = 0,
ELECTRON = 11,
MUON = 13,
TAUON = 15,
TAUONDECAY = 10015,
PROTON = 2212,
PION = 211,
KAONCHARGE = 321,
KAONNEUTRAL = 310,
A2 = 115,
ETA = 221,
F2 = 225,
ETAPRIME = 331,
F2PRIME = 335,
ETAC = 441,
F0 = 9010221,
ZOVERZ03 = 33,
RHO = 113,
+ RHO_ee = 113011,
+ RHO_mumu = 113013,
RHOZEUS = 913,
FOURPRONG = 999,
OMEGA = 223,
PHI = 333,
JPSI = 443,
JPSI_ee = 443011,
JPSI_mumu = 443013,
JPSI_ppbar = 4432212,
JPSI2S = 444,
JPSI2S_ee = 444011,
JPSI2S_mumu = 444013,
UPSILON = 553,
UPSILON_ee = 553011,
UPSILON_mumu = 553013,
UPSILON2S = 554,
UPSILON2S_ee = 554011,
UPSILON2S_mumu = 554013,
UPSILON3S = 555,
UPSILON3S_ee = 555011,
UPSILON3S_mumu = 555013,
AXION = 88, //AXION HACK
PHOTON = 22
};
enum decayTypeEnum {
NOTKNOWN = 0,
NARROWVMDEFAULT = 1,
WIDEVMDEFAULT = 2,
PSIFAMILY = 3,
LEPTONPAIR = 4,
SINGLEMESON = 5
};
enum interactionTypeEnum {
UNSPECIFIED = 0,
PHOTONPHOTON = 1,
PHOTONPOMERONNARROW = 2,
PHOTONPOMERONWIDE = 3,
PHOTONPOMERONINCOHERENT = 4,
PHOTONUCLEARSINGLE = 5,
PHOTONUCLEARDOUBLE = 6,
PHOTONUCLEARSINGLEPA = 7,
PHOTONUCLEARSINGLEPAPY = 8
};
enum systemTypeEnum{
NONSTANDARD = 0,
PP = 1,
PA = 2,
AA = 3
};
enum targetTypeEnum {
NOTHADRON = 0,
NUCLEUS = 1, //coherent gamma+A
NUCLEON = 2 //gamma+p or incoherent gamma+A
};
//Structure for each event's set of tracks.
struct event{
public:
int _numberOfTracks;
double px[30],py[30],pz[30];
int _fsParticle[30];
int _charge[30];
//To help track mothers and daughters produced through pythia.
int _mother1[30];
int _mother2[30];
int _daughter1[30];
int _daughter2[30];
//Normally we just set vertices to 0
//But for pythia, we decay additional states
int _numberOfVertices;
double _vertx[10],_verty[10],_vertz[10];
};
} // starlightConstants
#endif // STARLIGHTCONSTANTS_H_INCLUDE
Index: trunk/src/gammaavm.cpp
===================================================================
--- trunk/src/gammaavm.cpp (revision 308)
+++ trunk/src/gammaavm.cpp (revision 309)
@@ -1,949 +1,951 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: $: revision of last commit
// $Author:: $: author of last commit
// $Date:: $: date of last commit
//
// Description:
// Added incoherent t2-> pt2 selection. Following pp selection scheme
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include <cassert>
#include <cmath>
#include "gammaavm.h"
#include "photonNucleusCrossSection.h"
#include "wideResonanceCrossSection.h"
#include "narrowResonanceCrossSection.h"
#include "incoherentVMCrossSection.h"
using namespace std;
//______________________________________________________________________________
Gammaavectormeson::Gammaavectormeson(const inputParameters& inputParametersInstance, randomGenerator* randy, beamBeamSystem& bbsystem):eventChannel(inputParametersInstance, randy, bbsystem), _phaseSpaceGen(0)
{
_VMNPT=inputParametersInstance.nmbPtBinsInterference();
_VMWmax=inputParametersInstance.maxW();
_VMWmin=inputParametersInstance.minW();
_VMYmax=inputParametersInstance.maxRapidity();
_VMYmin=-1.*_VMYmax;
_VMnumw=inputParametersInstance.nmbWBins();
_VMnumy=inputParametersInstance.nmbRapidityBins();
_VMgamma_em=inputParametersInstance.beamLorentzGamma();
_VMinterferencemode=inputParametersInstance.interferenceEnabled();
_VMbslope=0.;//Will define in wide/narrow constructor
_bslopeDef=inputParametersInstance.bslopeDefinition();
_bslopeVal=inputParametersInstance.bslopeValue();
_pEnergy= inputParametersInstance.protonEnergy();
_VMpidtest=inputParametersInstance.prodParticleType();
_VMptmax=inputParametersInstance.maxPtInterference();
_VMdpt=inputParametersInstance.ptBinWidthInterference();
_ProductionMode=inputParametersInstance.productionMode();
N0 = 0; N1 = 0; N2 = 0;
if (_VMpidtest == starlightConstants::FOURPRONG){
// create n-body phase-spage generator
_phaseSpaceGen = new nBodyPhaseSpaceGen(_randy);
}
}
//______________________________________________________________________________
Gammaavectormeson::~Gammaavectormeson()
{
if (_phaseSpaceGen)
delete _phaseSpaceGen;
}
//______________________________________________________________________________
void Gammaavectormeson::pickwy(double &W, double &Y)
{
double dW, dY, xw,xy,xtest,btest;
int IW,IY;
dW = (_VMWmax-_VMWmin)/double(_VMnumw);
dY = (_VMYmax-_VMYmin)/double(_VMnumy);
L201pwy:
xw = _randy->Rndom();
W = _VMWmin + xw*(_VMWmax-_VMWmin);
if (W < 2 * _ip->pionChargedMass())
goto L201pwy;
IW = int((W-_VMWmin)/dW);
xy = _randy->Rndom();
Y = _VMYmin + xy*(_VMYmax-_VMYmin);
IY = int((Y-_VMYmin)/dY);
xtest = _randy->Rndom();
if( xtest > _Farray[IW][IY] )
goto L201pwy;
N0++;
// Determine the target nucleus
// For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
if( _ProductionMode == 2 || _ProductionMode ==3){
_TargetBeam = 2;
} else {
_TargetBeam = 1;
}
} else if( _bbs.beam1().A() != 1 && _bbs.beam2().A()==1 ){
if( _ProductionMode == 2 || _ProductionMode ==3){
_TargetBeam = 1;
} else {
_TargetBeam = 2;
}
} else {
btest = _randy->Rndom();
if ( btest < _Farray1[IW][IY]/_Farray[IW][IY] ){
_TargetBeam = 2;
N2++;
} else {
_TargetBeam = 1;
N1++;
}
}
}
//______________________________________________________________________________
void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
double W,
double px0, double py0, double pz0,
double& px1, double& py1, double& pz1,
double& px2, double& py2, double& pz2,
int& iFbadevent)
{
// This routine decays a particle into two particles of mass mdec,
// taking spin into account
double pmag;
double phi,theta,Ecm;
double betax,betay,betaz;
double mdec=0.0;
double E1=0.0,E2=0.0;
// set the mass of the daughter particles
mdec=getDaughterMass(ipid);
// calculate the magnitude of the momenta
if(W < 2*mdec){
cout<<" ERROR: W="<<W<<endl;
iFbadevent = 1;
return;
}
pmag = sqrt(W*W/4. - mdec*mdec);
// pick an orientation, based on the spin
// phi has a flat distribution in 2*pi
phi = _randy->Rndom()*2.*starlightConstants::pi;
// find theta, the angle between one of the outgoing particles and
// the beamline, in the outgoing particles' rest frame
theta=getTheta(ipid);
// compute unboosted momenta
px1 = sin(theta)*cos(phi)*pmag;
py1 = sin(theta)*sin(phi)*pmag;
pz1 = cos(theta)*pmag;
px2 = -px1;
py2 = -py1;
pz2 = -pz1;
Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
betax = -(px0/Ecm);
betay = -(py0/Ecm);
betaz = -(pz0/Ecm);
transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
if(iFbadevent == 1)
return;
}
//______________________________________________________________________________
// decays a particle into four particles with isotropic angular distribution
bool Gammaavectormeson::fourBodyDecay
(starlightConstants::particleTypeEnum& ipid,
const double , // E (unused)
const double W, // mass of produced particle
const double* p, // momentum of produced particle; expected to have size 3
lorentzVector* decayVecs, // array of Lorentz vectors of daughter particles; expected to have size 4
int& iFbadevent)
{
const double parentMass = W;
// set the mass of the daughter particles
const double daughterMass = getDaughterMass(ipid);
if (parentMass < 4 * daughterMass){
cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
iFbadevent = 1;
return false;
}
// construct parent four-vector
const double parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
+ parentMass * parentMass);
const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
// setup n-body phase-space generator
assert(_phaseSpaceGen);
static bool firstCall = true;
if (firstCall) {
const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
_phaseSpaceGen->setDecay(4, m);
// estimate maximum phase-space weight
_phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
firstCall = false;
}
// generate phase-space event
if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
return false;
// set Lorentzvectors of decay daughters
for (unsigned int i = 0; i < 4; ++i)
decayVecs[i] = _phaseSpaceGen->daughter(i);
return true;
}
//______________________________________________________________________________
double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
{
//This will return the daughter particles mass, and the final particles outputed id...
double mdec=0.;
switch(_VMpidtest){
case starlightConstants::RHO:
case starlightConstants::RHOZEUS:
case starlightConstants::FOURPRONG:
case starlightConstants::OMEGA:
mdec = _ip->pionChargedMass();
ipid = starlightConstants::PION;
break;
case starlightConstants::PHI:
mdec = _ip->kaonChargedMass();
ipid = starlightConstants::KAONCHARGE;
break;
case starlightConstants::JPSI:
mdec = _ip->mel();
ipid = starlightConstants::ELECTRON;
break;
+ case starlightConstants::RHO_ee:
case starlightConstants::JPSI_ee:
mdec = _ip->mel();
ipid = starlightConstants::ELECTRON;
break;
+ case starlightConstants::RHO_mumu:
case starlightConstants::JPSI_mumu:
mdec = _ip->muonMass();
ipid = starlightConstants::MUON;
break;
case starlightConstants::JPSI_ppbar:
mdec = _ip->protonMass();
ipid = starlightConstants::PROTON;
break;
case starlightConstants::JPSI2S_ee:
mdec = _ip->mel();
ipid = starlightConstants::ELECTRON;
break;
case starlightConstants::JPSI2S_mumu:
mdec = _ip->muonMass();
ipid = starlightConstants::MUON;
break;
case starlightConstants::JPSI2S:
case starlightConstants::UPSILON:
case starlightConstants::UPSILON2S:
case starlightConstants::UPSILON3S:
mdec = _ip->muonMass();
ipid = starlightConstants::MUON;
break;
case starlightConstants::UPSILON_ee:
case starlightConstants::UPSILON2S_ee:
case starlightConstants::UPSILON3S_ee:
mdec = _ip->mel();
ipid = starlightConstants::ELECTRON;
break;
case starlightConstants::UPSILON_mumu:
case starlightConstants::UPSILON2S_mumu:
case starlightConstants::UPSILON3S_mumu:
mdec = _ip->muonMass();
ipid = starlightConstants::MUON;
break;
default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
}
return mdec;
}
//______________________________________________________________________________
double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid)
{
//This depends on the decay angular distribution
//Valid for rho, phi, omega.
double theta=0.;
double xtest=0.;
double dndtheta=0.;
L200td:
theta = starlightConstants::pi*_randy->Rndom();
xtest = _randy->Rndom();
// Follow distribution for helicity +/-1
// Eq. 19 of J. Breitweg et al., Eur. Phys. J. C2, 247 (1998)
// SRK 11/14/2000
switch(ipid){
case starlightConstants::MUON:
case starlightConstants::ELECTRON:
//VM->ee/mumu
dndtheta = sin(theta)*(1.+(cos(theta)*cos(theta)));
break;
case starlightConstants::PROTON:
//Pick this angular distribution for J/psi --> ppbar
dndtheta = sin(theta)*(1.+(0.605*cos(theta)*cos(theta)));
break;
case starlightConstants::PION:
case starlightConstants::KAONCHARGE:
//rhos etc
dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
break;
default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
}//end of switch
if(xtest > dndtheta)
goto L200td;
return theta;
}
//______________________________________________________________________________
double Gammaavectormeson::getWidth()
{
return _width;
}
//______________________________________________________________________________
double Gammaavectormeson::getMass()
{
return _mass;
}
//______________________________________________________________________________
double Gammaavectormeson::getSpin()
{
return 1.0; //VM spins are the same
}
//______________________________________________________________________________
void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)
{
// This subroutine calculates momentum and energy of vector meson
// given W and Y, without interference. Subroutine vmpt handles
// production with interference
double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
double px1,py1,px2,py2;
double pt,xt,xtest,ytest;
double t2;
//Find Egam,Epom in CM frame
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
// This is pA
if( _ProductionMode == 2 || _ProductionMode ==3 ){
Egam = 0.5*W*exp(Y);
Epom = 0.5*W*exp(-Y);
}else{
Egam = 0.5*W*exp(-Y);
Epom = 0.5*W*exp(Y);
}
} else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
// This is Ap
if( _ProductionMode == 2 || _ProductionMode == 3 ){
Egam = 0.5*W*exp(-Y);
Epom = 0.5*W*exp(Y);
}else{
Egam = 0.5*W*exp(Y);
Epom = 0.5*W*exp(-Y);
}
} else {
// This is pp or AA
if( _TargetBeam == 1 ){
Egam = 0.5*W*exp(-Y);
Epom = 0.5*W*exp(Y);
}
else {
Egam = 0.5*W*exp(Y);
Epom = 0.5*W*exp(-Y);
}
}
// } else if( _ProductionMode == 2 || _ProductionMode==3){
// Egam = 0.5*W*exp(-Y);
// Epom = 0.5*W*exp(Y);
// } else {
// Egam = 0.5*W*exp(Y);
// Epom = 0.5*W*exp(-Y);
// }
pt1 = pTgamma(Egam);
phi1 = 2.*starlightConstants::pi*_randy->Rndom();
if( (_bbs.beam1().A()==1 && _bbs.beam2().A()==1) ||
(_ProductionMode == 4) ) {
if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
// Use dipole form factor for light VM
L555vm:
xtest = 2.0*_randy->Rndom();
double ttest = xtest*xtest;
ytest = _randy->Rndom();
double t0 = 1./2.23;
double yprob = xtest*_bbs.beam1().dipoleFormFactor(ttest,t0)*_bbs.beam1().dipoleFormFactor(ttest,t0);
if( ytest > yprob ) goto L555vm;
t2 = ttest;
pt2 = xtest;
}else{
//Use dsig/dt= exp(-_VMbslope*t) for heavy VM
double bslope_tdist = _VMbslope;
double Wgammap = 0.0;
switch(_bslopeDef){
case 0:
//This is the default, as before
bslope_tdist = _VMbslope;
break;
case 1:
//User defined value of bslope. BSLOPE_VALUE default is 4.0 if not set.
bslope_tdist = _bslopeVal;
if( N0 <= 1 )cout<<" ATTENTION: Using user defined value of bslope = "<<_bslopeVal<<endl;
break;
case 2:
//This is Wgammap dependence of b from H1 (Eur. Phys. J. C 46 (2006) 585)
Wgammap = sqrt(4.*Egam*_pEnergy);
bslope_tdist = 4.63 + 4.*0.164*log(Wgammap/90.0);
if( N0 <= 1 )cout<<" ATTENTION: Using energy dependent value of bslope!"<<endl;
break;
default:
cout<<" Undefined setting for BSLOPE_DEFINITION "<<endl;
}
xtest = _randy->Rndom();
// t2 = (-1./_VMbslope)*log(xtest);
t2 = (-1./bslope_tdist)*log(xtest);
pt2 = sqrt(1.*t2);
}
} else {
// >> Check tmin
tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
if(tmin > 0.5){
cout<<" WARNING: tmin= "<<tmin<<endl;
cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl;
cout<<" Will pick a new W,Y "<<endl;
tcheck = 1;
return;
}
L203vm:
xt = _randy->Rndom();
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
if( _ProductionMode == 2 || _ProductionMode ==3){
// Changed '8' to '32' 6 times below to extend the region of the p_T calculation up to 1 GeV.c SRK May 28, 2019
// 'pt2' is the maximum vector meson momentum. For heavy nuclei, the '32'coefficient corresonds to about 1 GeV/c
// The downside of the larger coefficient is that the sampling runs more slowly. This could be made into a parameter
pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
}else{
pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
}
} else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
if( _ProductionMode == 2 || _ProductionMode ==3){
pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
}else{
pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
}
} else if (_TargetBeam==1) {
pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
} else {
pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
}
xtest = _randy->Rndom();
t2 = tmin + pt2*pt2;
double comp=0.0;
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
if( _ProductionMode == 2 || _ProductionMode ==3){
comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
}else{
comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
}
} else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
if( _ProductionMode == 2 || _ProductionMode ==3){
comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
}else{
comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
}
} else if (_TargetBeam==1) {
comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
} else {
comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
}
if( xtest > comp ) goto L203vm;
}//else end from pp
phi2 = 2.*starlightConstants::pi*_randy->Rndom();
px1 = pt1*cos(phi1);
py1 = pt1*sin(phi1);
px2 = pt2*cos(phi2);
py2 = pt2*sin(phi2);
// Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
px = px1 + px2;
py = py1 + py2;
pt = sqrt( px*px + py*py );
E = sqrt(W*W+pt*pt)*cosh(Y);
pz = sqrt(W*W+pt*pt)*sinh(Y);
}
//______________________________________________________________________________
double Gammaavectormeson::pTgamma(double E)
{
// returns on random draw from pp(E) distribution
double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
int satisfy =0;
ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
//sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
Cm = sqrt(3.)*E/_VMgamma_em;
// If E is very small, the drawing of a pT below is extre_ip->mel()y slow.
// ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E.
// Should have no observable consequences (JN, SRK 11-Sep-2014)
if( E < 0.0005 )return Cm;
//the amplitude of the p_t spectrum at the maximum
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
if( _ProductionMode == 2 || _ProductionMode ==3 ){
singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
}else{
singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
}
} else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
if( _ProductionMode == 2 || _ProductionMode ==3){
singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
}else{
singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
}
} else if (_TargetBeam == 1) {
singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
} else {
singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
}
Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
//pick a test value pp, and find the amplitude there
x = _randy->Rndom();
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
if( _ProductionMode == 2 || _ProductionMode ==3){
pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
}else{
pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
}
} else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
if( _ProductionMode == 2 || _ProductionMode ==3){
pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
}else{
pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
}
} else if (_TargetBeam == 1) {
pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
} else {
pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
}
test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
while(satisfy==0){
u = _randy->Rndom();
if(u*Coef <= test)
{
satisfy =1;
}
else{
x =_randy->Rndom();
if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
if( _ProductionMode == 2 || _ProductionMode ==3){
pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
}else{
pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
}
} else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
if( _ProductionMode == 2 || _ProductionMode ==3){
pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
}else{
pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
}
} else if (_TargetBeam == 1) {
pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
} else {
pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
}
test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
}
}
return pp;
}
//______________________________________________________________________________
void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
int&) // tcheck (unused)
{
// This function calculates momentum and energy of vector meson
// given W and Y, including interference.
// It gets the pt distribution from a lookup table.
double dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
int IY=0,j=0;
dY = (_VMYmax-_VMYmin)/double(_VMnumy);
// Y is already fixed; choose a pt
// Follow the approach in pickwy
// in _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
// Changed, now works -y to +y.
IY=int((Y-_VMYmin)/dY);
if (IY > (_VMnumy)-1){
IY=(_VMnumy)-1;
}
yleft=(Y-_VMYmin)-(IY)*dY;
yfract=yleft*dY;
xpt=_randy->Rndom();
for(j=0;j<_VMNPT;j++){
if (xpt < _fptarray[IY][j]) goto L60;
}
if(j == _VMNPT) j = _VMNPT-1;
L60:
// now do linear interpolation - start with extremes
if (j == 0){
pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
goto L80;
}
if (j == _VMNPT-1){
pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
goto L80;
}
// we're in the middle
ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
// at an extreme in y?
if (IY == (_VMnumy)-1){
pt=pt1;
goto L120;
}
L80:
// interpolate in y repeat for next fractional y bin
for(j=0;j<_VMNPT;j++){
if (xpt < _fptarray[IY+1][j]) goto L90;
}
if(j == _VMNPT) j = _VMNPT-1;
L90:
// now do linear interpolation - start with extremes
if (j == 0){
pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
goto L100;
}
if (j == _VMNPT-1){
pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
goto L100;
}
// we're in the middle
ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
L100:
// now interpolate in y
pt=yfract*pt2+(1-yfract)*pt1;
L120:
// we have a pt
theta=2.*starlightConstants::pi*_randy->Rndom();
px=pt*cos(theta);
py=pt*sin(theta);
E = sqrt(W*W+pt*pt)*cosh(Y);
pz = sqrt(W*W+pt*pt)*sinh(Y);
// randomly choose to make pz negative 50% of the time
if(_randy->Rndom()>=0.5) pz = -pz;
}
//______________________________________________________________________________
starlightConstants::event Gammaavectormeson::produceEvent(int&)
{
//Note used; return default event
return starlightConstants::event();
}
//______________________________________________________________________________
upcEvent Gammaavectormeson::produceEvent()
{
// The new event type
upcEvent event;
int iFbadevent=0;
int tcheck=0;
starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
if (_VMpidtest == starlightConstants::FOURPRONG) {
double comenergy = 0;
double mom[3] = {0, 0, 0};
double E = 0;
lorentzVector decayVecs[4];
do {
double rapidity = 0;
pickwy(comenergy, rapidity);
if (_VMinterferencemode == 0)
momenta(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
else if (_VMinterferencemode==1)
vmpt(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
} while (!fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent));
if ((iFbadevent == 0) and (tcheck == 0))
for (unsigned int i = 0; i < 4; ++i) {
starlightParticle daughter(decayVecs[i].GetPx(),
decayVecs[i].GetPy(),
decayVecs[i].GetPz(),
sqrt(decayVecs[i].GetPx()*decayVecs[i].GetPx()+decayVecs[i].GetPy()*decayVecs[i].GetPy()+decayVecs[i].GetPz()*decayVecs[i].GetPz()+0.139*0.139),//energy
starlightConstants::UNKNOWN, // _mass
ipid*(2*(i/2)-1), // make half of the particles pi^+, half pi^-
i/2);
event.addParticle(daughter);
}
} else {
double comenergy = 0.;
double rapidity = 0.;
double E = 0.;
double momx=0.,momy=0.,momz=0.;
double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
bool accepted = false;
do{
pickwy(comenergy,rapidity);
if (_VMinterferencemode==0){
momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
} else if (_VMinterferencemode==1){
vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck);
}
_nmbAttempts++;
vmpid = ipid;
twoBodyDecay(ipid,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
double pt1chk = sqrt(px1*px1+py1*py1);
double pt2chk = sqrt(px2*px2+py2*py2);
double eta1 = pseudoRapidity(px1, py1, pz1);
double eta2 = pseudoRapidity(px2, py2, pz2);
if(_ptCutEnabled && !_etaCutEnabled){
if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
accepted = true;
_nmbAccepted++;
}
}
else if(!_ptCutEnabled && _etaCutEnabled){
if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
accepted = true;
_nmbAccepted++;
}
}
else if(_ptCutEnabled && _etaCutEnabled){
if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
accepted = true;
_nmbAccepted++;
}
}
}
else if(!_ptCutEnabled && !_etaCutEnabled)
_nmbAccepted++;
}while((_ptCutEnabled || _etaCutEnabled) && !accepted);
if (iFbadevent==0&&tcheck==0) {
int q1=0,q2=0;
int ipid1,ipid2=0;
double xtest = _randy->Rndom();
if (xtest<0.5)
{
q1=1;
q2=-1;
}
else {
q1=-1;
q2=1;
}
if ( ipid == 11 || ipid == 13 ){
ipid1 = -q1*ipid;
ipid2 = -q2*ipid;
} else {
ipid1 = q1*ipid;
ipid2 = q2*ipid;
}
double md = getDaughterMass(vmpid);
double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1);
starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
event.addParticle(particle1);
double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2);
starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
event.addParticle(particle2);
}
}
return event;
}
double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
{
double pT = sqrt(px*px + py*py);
double p = sqrt(pz*pz + pT*pT);
double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
return eta;
}
//______________________________________________________________________________
Gammaanarrowvm::Gammaanarrowvm(const inputParameters& input, randomGenerator* randy, beamBeamSystem& bbsystem):Gammaavectormeson(input, randy, bbsystem)
{
cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
read();
cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
narrowResonanceCrossSection sigma(input, bbsystem);
sigma.crossSectionCalculation(_bwnormsave);
setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
_VMbslope=sigma.slopeParameter();
}
//______________________________________________________________________________
Gammaanarrowvm::~Gammaanarrowvm()
{ }
//______________________________________________________________________________
Gammaaincoherentvm::Gammaaincoherentvm(const inputParameters& input, randomGenerator* randy, beamBeamSystem& bbsystem):Gammaavectormeson(input, randy, bbsystem)
{
cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
read();
cout<<"Creating and calculating crosssection. Gammaaincoherentvm()"<<endl;
incoherentVMCrossSection sigma(input, bbsystem);
sigma.crossSectionCalculation(_bwnormsave);
setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
_VMbslope=sigma.slopeParameter();
}
//______________________________________________________________________________
Gammaaincoherentvm::~Gammaaincoherentvm()
{ }
//______________________________________________________________________________
Gammaawidevm::Gammaawidevm(const inputParameters& input, randomGenerator* randy, beamBeamSystem& bbsystem):Gammaavectormeson(input, randy, bbsystem)
{
cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
read();
cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
wideResonanceCrossSection sigma(input, bbsystem);
sigma.crossSectionCalculation(_bwnormsave);
setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
_VMbslope=sigma.slopeParameter();
}
//______________________________________________________________________________
Gammaawidevm::~Gammaawidevm()
{ }
Index: trunk/src/starlight.cpp
===================================================================
--- trunk/src/starlight.cpp (revision 308)
+++ trunk/src/starlight.cpp (revision 309)
@@ -1,386 +1,388 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: $: revision of last commit
// $Author:: $: author of last commit
// $Date:: $: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#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 "incoherentPhotonNucleusLuminosity.h"
#include "upcevent.h"
#include "eventfilewriter.h"
#include "starlight.h"
using namespace std;
using namespace starlightConstants;
starlight::starlight() :
_beamSystem (0),
_eventChannel (0),
_nmbEventsPerFile (100),
_isInitialised (false),
_inputParameters (0),
_randomGenerator (0)
{ }
starlight::~starlight()
{ }
bool
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->beam1Z() > _inputParameters->beam1A() ){
printErr << endl << "A must be >= Z; A beam1 = "<<_inputParameters->beam1A()<<", Z beam1 = "<<_inputParameters->beam1Z()<<". Terminating."<<endl ;
return false;
}
if( _inputParameters->beam2Z() > _inputParameters->beam2A() ){
printErr << endl << "A must be >= Z; A beam2 = "<<_inputParameters->beam2A()<<", Z beam2 = "<<_inputParameters->beam2Z()<<". Terminating."<<endl ;
return false;
}
if( _inputParameters->interactionType() == PHOTONPOMERONINCOHERENT && _inputParameters->beam1A() == 1 &&
_inputParameters->beam1Z() == 1 && _inputParameters->beam2A() == 1 && _inputParameters->beam2Z() ){
printErr << endl << " Do not use PROD_MODE = 4 for pp collisions. Use PROD_MODE = 2 or 3 instead. Terminating."<<endl;
return false;
}
bool createChannel = true;
switch (_inputParameters->interactionType()) {
case PHOTONPHOTON:
if (!lumTableIsValid) {
printInfo << "creating luminosity table for photon-photon channel" << endl;
twoPhotonLuminosity(*_inputParameters, _beamSystem->beam1(), _beamSystem->beam2());
}
break;
case PHOTONPOMERONNARROW: // narrow and wide resonances use
case PHOTONPOMERONWIDE: // the same luminosity function
if (!lumTableIsValid) {
printInfo << "creating luminosity table for coherent photon-Pomeron channel" <<endl;
photonNucleusLuminosity lum(*_inputParameters, *_beamSystem);
}
break;
case PHOTONPOMERONINCOHERENT: // narrow and wide resonances use
if (!lumTableIsValid) {
printInfo << "creating luminosity table for incoherent photon-Pomeron channel" << endl;
incoherentPhotonNucleusLuminosity lum(*_inputParameters, *_beamSystem);
}
break;
#ifdef ENABLE_DPMJET
case PHOTONUCLEARSINGLE:
createChannel = false;
_eventChannel = new starlightDpmJet(*_inputParameters,_randomGenerator,*_beamSystem);
std::cout << "CREATING PHOTONUCLEAR/DPMJET SINGLE" << std::endl;
dynamic_cast<starlightDpmJet*>(_eventChannel)->setSingleMode();
dynamic_cast<starlightDpmJet*>(_eventChannel)->setMinGammaEnergy(_inputParameters->minGammaEnergy());
dynamic_cast<starlightDpmJet*>(_eventChannel)->setMaxGammaEnergy(_inputParameters->maxGammaEnergy());
dynamic_cast<starlightDpmJet*>(_eventChannel)->init();
break;
case PHOTONUCLEARDOUBLE:
createChannel = false;
_eventChannel = new starlightDpmJet(*_inputParameters,_randomGenerator,*_beamSystem);
std::cout << "CREATING PHOTONUCLEAR/DPMJET DOUBLE" << std::endl;
dynamic_cast<starlightDpmJet*>(_eventChannel)->setDoubleMode();
dynamic_cast<starlightDpmJet*>(_eventChannel)->setMinGammaEnergy(_inputParameters->minGammaEnergy());
dynamic_cast<starlightDpmJet*>(_eventChannel)->setMaxGammaEnergy(_inputParameters->maxGammaEnergy());
dynamic_cast<starlightDpmJet*>(_eventChannel)->init();
break;
case PHOTONUCLEARSINGLEPA:
createChannel = false;
_eventChannel = new starlightDpmJet(*_inputParameters,_randomGenerator,*_beamSystem);
std::cout << "CREATING PHOTONUCLEAR/DPMJET SINGLE" << std::endl;
dynamic_cast<starlightDpmJet*>(_eventChannel)->setSingleMode();
dynamic_cast<starlightDpmJet*>(_eventChannel)->setProtonMode();
dynamic_cast<starlightDpmJet*>(_eventChannel)->setMinGammaEnergy(_inputParameters->minGammaEnergy());
dynamic_cast<starlightDpmJet*>(_eventChannel)->setMaxGammaEnergy(_inputParameters->maxGammaEnergy());
dynamic_cast<starlightDpmJet*>(_eventChannel)->init();
break;
#endif
#ifdef ENABLE_PYTHIA6
case PHOTONUCLEARSINGLEPAPY:
createChannel = false;
_eventChannel = new starlightPythia(*_inputParameters,randomGenerator,*_beamSystem);
std::cout << "CREATING PHOTONUCLEAR/PYTHIA SINGLE" << std::endl;
dynamic_cast<starlightPythia*>(_eventChannel)->setSingleMode();
dynamic_cast<starlightPythia*>(_eventChannel)->setMinGammaEnergy(_inputParameters->minGammaEnergy());
dynamic_cast<starlightPythia*>(_eventChannel)->setMaxGammaEnergy(_inputParameters->maxGammaEnergy());
dynamic_cast<starlightPythia*>(_eventChannel)->init(_inputParameters->pythiaParams(), _inputParameters->pythiaFullEventRecord());
break;
#endif
default:
{
printWarn << "unknown interaction type '" << _inputParameters->interactionType() << "'."
<< " cannot initialize starlight." << endl;
return false;
}
}
if(createChannel)
{
if (!createEventChannel())
return false;
}
_isInitialised = true;
return true;
}
upcEvent
starlight::produceEvent()
{
if (!_isInitialised) {
printErr << "trying to generate event but Starlight is not initialised. aborting." << endl;
exit(-1);
}
return _eventChannel->produceEvent();
}
bool
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;
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
>> interferenceEnabled >> interferenceStrength
>> 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
starlight::createEventChannel()
{
switch (_inputParameters->prodParticleType()) {
case ELECTRON:
case MUON:
case TAUON:
case TAUONDECAY:
{
_eventChannel = new Gammagammaleptonpair(*_inputParameters, _randomGenerator, *_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
// PythiaOutput = true;
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
}
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, _randomGenerator, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammagammasingle event channel." << endl;
return false;
}
}
case RHO:
+ case RHO_ee:
+ case RHO_mumu:
case RHOZEUS:
case FOURPRONG:
case OMEGA:
case PHI:
case JPSI:
case JPSI_ee:
case JPSI_mumu:
case JPSI_ppbar:
case JPSI2S:
case JPSI2S_ee:
case JPSI2S_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() == PHOTONPOMERONNARROW) {
_eventChannel = new Gammaanarrowvm(*_inputParameters, _randomGenerator, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammaanarrowvm event channel." << endl;
return false;
}
}
if (_inputParameters->interactionType() == PHOTONPOMERONWIDE) {
_eventChannel = new Gammaawidevm(*_inputParameters, _randomGenerator, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammaawidevm event channel." << endl;
return false;
}
}
if (_inputParameters->interactionType() == PHOTONPOMERONINCOHERENT) {
_eventChannel = new Gammaaincoherentvm(*_inputParameters, _randomGenerator, *_beamSystem);
if (_eventChannel)
return true;
else {
printWarn << "cannot construct Gammaanarrowvm 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/photonNucleusCrossSection.cpp
===================================================================
--- trunk/src/photonNucleusCrossSection.cpp (revision 308)
+++ trunk/src/photonNucleusCrossSection.cpp (revision 309)
@@ -1,926 +1,950 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: $: revision of last commit
// $Author:: $: author of last commit
// $Date:: $: date of last commit
//
// Description:
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include <cmath>
#include "reportingUtils.h"
#include "starlightconstants.h"
#include "bessel.h"
#include "photonNucleusCrossSection.h"
using namespace std;
using namespace starlightConstants;
//______________________________________________________________________________
photonNucleusCrossSection::photonNucleusCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
: _ip(&inputParametersInstance),
_nWbins (inputParametersInstance.nmbWBins() ),
_nYbins (inputParametersInstance.nmbRapidityBins() ),
_wMin (inputParametersInstance.minW() ),
_wMax (inputParametersInstance.maxW() ),
_yMax (inputParametersInstance.maxRapidity() ),
_beamLorentzGamma (inputParametersInstance.beamLorentzGamma() ),
_bbs (bbsystem ),
_protonEnergy (inputParametersInstance.protonEnergy() ),
_particleType (inputParametersInstance.prodParticleType() ),
_beamBreakupMode (inputParametersInstance.beamBreakupMode() ),
_productionMode (inputParametersInstance.productionMode() ),
_sigmaNucleus (_bbs.beam2().A() )
{
// new options - impulse aproximation (per Joakim) and Quantum Glauber (per SK) SKQG
_impulseSelected = inputParametersInstance.impulseVM();
_quantumGlauber = inputParametersInstance.quantumGlauber();
switch(_particleType) {
case RHO:
+ case RHO_ee:
+ case RHO_mumu:
_slopeParameter = 11.0; // [(GeV/c)^{-2}]
_vmPhotonCoupling = 2.02;
_ANORM = -2.75;
_BNORM = 0.0;
_defaultC = 1.0;
_channelMass = _ip->rho0Mass();
_width = _ip->rho0Width();
break;
case RHOZEUS:
_slopeParameter =11.0;
_vmPhotonCoupling=2.02;
_ANORM=-2.75;
_BNORM=1.84;
_defaultC=1.0;
_channelMass = _ip->rho0Mass();
_width = _ip->rho0Width();
break;
case FOURPRONG:
_slopeParameter = 11.0;
_vmPhotonCoupling = 2.02;
_ANORM = -2.75;
_BNORM = 0;
_defaultC = 1.0;
_channelMass = _ip->rho0PrimeMass();
_width = _ip->rho0PrimeWidth();
break;
case OMEGA:
_slopeParameter=10.0;
_vmPhotonCoupling=23.13;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->OmegaMass();
_width = _ip->OmegaWidth();
break;
case PHI:
_slopeParameter=7.0;
_vmPhotonCoupling=13.71;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->PhiMass();
_width = _ip->PhiWidth();
break;
case JPSI:
case JPSI_ee:
case JPSI_mumu:
case JPSI_ppbar:
_slopeParameter=4.0;
_vmPhotonCoupling=10.45;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->JpsiMass();
_width = _ip->JpsiWidth();
break;
case JPSI2S:
case JPSI2S_ee:
case JPSI2S_mumu:
_slopeParameter=4.3;
_vmPhotonCoupling=26.39;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->Psi2SMass();
_width = _ip->Psi2SWidth();
break;
case UPSILON:
case UPSILON_ee:
case UPSILON_mumu:
_slopeParameter=4.0;
_vmPhotonCoupling=125.37;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->Upsilon1SMass();
_width = _ip->Upsilon1SWidth();
break;
case UPSILON2S:
case UPSILON2S_ee:
case UPSILON2S_mumu:
_slopeParameter=4.0;
_vmPhotonCoupling=290.84;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->Upsilon2SMass();
_width = _ip->Upsilon2SWidth();
break;
case UPSILON3S:
case UPSILON3S_ee:
case UPSILON3S_mumu:
_slopeParameter=4.0;
_vmPhotonCoupling=415.10;
_ANORM=-2.75;
_BNORM=0.0;
_defaultC=1.0;
_channelMass = _ip->Upsilon3SMass();
_width = _ip->Upsilon3SWidth();
break;
default:
cout <<"No sigma constants parameterized for pid: "<<_particleType
<<" GammaAcrosssection"<<endl;
}
_maxPhotonEnergy = 12. * _beamLorentzGamma * hbarc/(_bbs.beam1().nuclearRadius()+_bbs.beam2().nuclearRadius());
}
//______________________________________________________________________________
photonNucleusCrossSection::~photonNucleusCrossSection()
{ }
//______________________________________________________________________________
void
photonNucleusCrossSection::crossSectionCalculation(const double)
{
cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::getcsgA(const double Egamma,
const double W,
const int beam)
{
//This function returns the cross-section for photon-nucleus interaction
//producing vectormesons
double Av,Wgp,cs,cvma;
double t,tmin,tmax;
double csgA,ax,bx;
int NGAUSS;
// DATA FOR GAUSS INTEGRATION
double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
NGAUSS = 6;
// Find gamma-proton CM energy
Wgp = sqrt(2. * Egamma * (_protonEnergy
+ sqrt(_protonEnergy * _protonEnergy - _ip->protonMass() * _ip->protonMass()))
+ _ip->protonMass() * _ip->protonMass());
//Used for A-A
tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
if ((_bbs.beam1().A() == 1) && (_bbs.beam2().A() == 1)){
// proton-proton, no scaling needed
csgA = sigmagp(Wgp);
} else {
// coherent AA interactions
int A_1 = _bbs.beam1().A();
int A_2 = _bbs.beam2().A();
// Calculate V.M.+proton cross section
// cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
cs = sigma_N(Wgp); //Use member function instead
// Calculate V.M.+nucleus cross section
cvma = sigma_A(cs,beam);
// Do impulse approximation here
if( _impulseSelected == 1){
if( beam == 1 ){
cvma = A_1*cs;
} else if ( beam == 2 ){
cvma = A_2*cs;
}
}
// Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
tmax = tmin + 0.25;
ax = 0.5 * (tmax - tmin);
bx = 0.5 * (tmax + tmin);
csgA = 0.;
for (int k = 1; k < NGAUSS; ++k) {
t = ax * xg[k] + bx;
if( A_1 == 1 && A_2 != 1){
csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
}else if(A_2 ==1 && A_1 != 1){
csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
}else{
if( beam==1 ){
csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
}else if(beam==2){
csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
}else{
cout<<"Something went wrong here, beam= "<<beam<<endl;
}
}
t = ax * (-xg[k]) + bx;
if( A_1 == 1 && A_2 != 1){
csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
}else if(A_2 ==1 && A_1 != 1){
csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
}else{
if( beam==1 ){
csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
}else if(beam==2){
csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
}else{
cout<<"Something went wrong here, beam= "<<beam<<endl;
}
}
}
csgA = 0.5 * (tmax - tmin) * csgA;
csgA = Av * csgA;
}
return csgA;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::photonFlux(const double Egamma, const int beam)
{
// This routine gives the photon flux as a function of energy Egamma
// It works for arbitrary nuclei and gamma; the first time it is
// called, it calculates a lookup table which is used on
// subsequent calls.
// It returns dN_gamma/dE (dimensions 1/E), not dI/dE
// energies are in GeV, in the lab frame
// rewritten 4/25/2001 by SRK
// NOTE: beam (=1,2) defines the photon TARGET
double lEgamma,Emin,Emax;
static double lnEmax, lnEmin, dlnE;
double stepmult,energy,rZ;
int nbstep,nrstep,nphistep,nstep;
double bmin,bmax,bmult,biter,bold,integratedflux;
double fluxelement,deltar,riter;
double deltaphi,phiiter,dist;
static double dide[401];
double lnElt;
double flux_r;
double Xvar;
int Ilt;
double RNuc=0.,RSum=0.;
RSum=_bbs.beam1().nuclearRadius()+_bbs.beam2().nuclearRadius();
if( beam == 1){
rZ=double(_bbs.beam2().Z());
RNuc = _bbs.beam1().nuclearRadius();
} else {
rZ=double(_bbs.beam1().Z());
RNuc = _bbs.beam2().nuclearRadius();
}
static int Icheck = 0;
static int Ibeam = 0;
//Check first to see if pp
if( _bbs.beam1().A()==1 && _bbs.beam2().A()==1 ){
int nbsteps = 400;
double bmin = 0.5;
double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
double local_sum=0.0;
// Impact parameter loop
for (int i = 0; i<=nbsteps;i++){
double bnn0 = bmin*exp(i*dlnb);
double bnn1 = bmin*exp((i+1)*dlnb);
double db = bnn1-bnn0;
double ppslope = 19.0;
double GammaProfile = exp(-bnn0*bnn0/(2.*hbarc*hbarc*ppslope));
double PofB0 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
GammaProfile = exp(-bnn1*bnn1/(2.*hbarc*hbarc*ppslope));
double PofB1 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
double loc_nofe0 = _bbs.beam1().photonDensity(bnn0,Egamma);
double loc_nofe1 = _bbs.beam2().photonDensity(bnn1,Egamma);
local_sum += 0.5*loc_nofe0*(1. - PofB0)*2.*starlightConstants::pi*bnn0*db;
local_sum += 0.5*loc_nofe1*(1. - PofB1)*2.*starlightConstants::pi*bnn1*db;
}
// End Impact parameter loop
return local_sum;
}
// first call or new beam? - initialize - calculate photon flux
Icheck=Icheck+1;
// Do the numerical integration only once for symmetric systems.
if( Icheck > 1 && _bbs.beam1().A() == _bbs.beam2().A() && _bbs.beam1().Z() == _bbs.beam2().Z() ) goto L1000f;
// For asymmetric systems check if we have another beam
if( Icheck > 1 && beam == Ibeam ) goto L1000f;
Ibeam = beam;
// Nuclear breakup is done by PofB
// collect number of integration steps here, in one place
nbstep=1200;
nrstep=60;
nphistep=40;
// this last one is the number of energy steps
nstep=100;
// following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
Emin=1.E-5;
if (_beamLorentzGamma < 500)
Emin=1.E-3;
Emax=_maxPhotonEnergy;
// Emax=12.*hbarc*_beamLorentzGamma/RSum;
// >> lnEmin <-> ln(Egamma) for the 0th bin
// >> lnEmax <-> ln(Egamma) for the last bin
lnEmin=log(Emin);
lnEmax=log(Emax);
dlnE=(lnEmax-lnEmin)/nstep;
printf("Calculating photon flux from Emin = %e GeV to Emax = %e GeV (CM frame) for source with Z = %3.0f \n", Emin, Emax, rZ);
stepmult= exp(log(Emax/Emin)/double(nstep));
energy=Emin;
for (int j = 1; j<=nstep;j++){
energy=energy*stepmult;
// integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
// use exponential steps
bmin=0.8*RSum; //Start slightly below 2*Radius
bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
bmult=exp(log(bmax/bmin)/double(nbstep));
biter=bmin;
integratedflux=0.;
if( (_bbs.beam1().A() == 1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A() == 1 && _bbs.beam1().A() != 1) ){
// This is pA
if( _productionMode == PHOTONPOMERONINCOHERENT ){
// This pA incoherent, proton is the target
int nbsteps = 400;
double bmin = 0.7*RSum;
double bmax = 2.0*RSum + (8.0*_beamLorentzGamma*hbarc/energy);
double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
double local_sum=0.0;
// Impact parameter loop
for (int i = 0; i<=nbsteps; i++){
double bnn0 = bmin*exp(i*dlnb);
double bnn1 = bmin*exp((i+1)*dlnb);
double db = bnn1-bnn0;
double PofB0 = _bbs.probabilityOfBreakup(bnn0);
double PofB1 = _bbs.probabilityOfBreakup(bnn1);
double loc_nofe0 = 0.0;
double loc_nofe1 = 0.0;
if( _bbs.beam1().A() == 1 ){
loc_nofe0 = _bbs.beam2().photonDensity(bnn0,energy);
loc_nofe1 = _bbs.beam2().photonDensity(bnn1,energy);
}
else if( _bbs.beam2().A() == 1 ){
loc_nofe0 = _bbs.beam1().photonDensity(bnn0,energy);
loc_nofe1 = _bbs.beam1().photonDensity(bnn1,energy);
}
// cout<<" i: "<<i<<" bnn0: "<<bnn0<<" PofB0: "<<PofB0<<" loc_nofe0: "<<loc_nofe0<<endl;
local_sum += 0.5*loc_nofe0*PofB0*2.*starlightConstants::pi*bnn0*db;
local_sum += 0.5*loc_nofe1*PofB1*2.*starlightConstants::pi*bnn1*db;
} // End Impact parameter loop
integratedflux = local_sum;
} else if ( _productionMode == PHOTONPOMERONNARROW || _productionMode == PHOTONPOMERONWIDE ){
// This is pA coherent, nucleus is the target
double localbmin = 0.0;
if( _bbs.beam1().A() == 1 ){
localbmin = _bbs.beam2().nuclearRadius() + 0.7;
}
if( _bbs.beam2().A() == 1 ){
localbmin = _bbs.beam1().nuclearRadius() + 0.7;
}
integratedflux = nepoint(energy,localbmin);
}
}else{
// This is AA
for (int jb = 1; jb<=nbstep;jb++){
bold=biter;
biter=biter*bmult;
// When we get to b>20R_A change methods - just take the photon flux
// at the center of the nucleus.
if (biter > (10.*RNuc)){
// if there is no nuclear breakup or only hadronic breakup, which only
// occurs at smaller b, we can analytically integrate the flux from b~20R_A
// to infinity, following Jackson (2nd edition), Eq. 15.54
Xvar=energy*biter/(hbarc*_beamLorentzGamma);
// Here, there is nuclear breakup. So, we can't use the integrated flux
// However, we can do a single flux calculation, at the center of the nucleus
// Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
// this is the flux per unit area
fluxelement = (rZ*rZ*alpha*energy)*
(bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
((pi*_beamLorentzGamma*hbarc)*
(pi*_beamLorentzGamma*hbarc));
} else {
// integrate over nuclear surface. n.b. this assumes total shadowing -
// treat photons hitting the nucleus the same no matter where they strike
fluxelement=0.;
deltar=RNuc/double(nrstep);
riter=-deltar/2.;
for (int jr =1; jr<=nrstep;jr++){
riter=riter+deltar;
// use symmetry; only integrate from 0 to pi (half circle)
deltaphi=pi/double(nphistep);
phiiter=0.;
for( int jphi=1;jphi<= nphistep;jphi++){
phiiter=(double(jphi)-0.5)*deltaphi;
// dist is the distance from the center of the emitting nucleus
// to the point in question
dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
Xvar=energy*dist/(hbarc*_beamLorentzGamma);
flux_r = (rZ*rZ*alpha*energy)*
(bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
((pi*_beamLorentzGamma*hbarc)*
(pi*_beamLorentzGamma*hbarc));
// The surface element is 2.* delta phi* r * delta r
// The '2' is because the phi integral only goes from 0 to pi
fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
// end phi and r integrations
}//for(jphi)
}//for(jr)
// average fluxelement over the nuclear surface
fluxelement=fluxelement/(pi*RNuc*RNuc);
}//else
// multiply by volume element to get total flux in the volume element
fluxelement=fluxelement*2.*pi*biter*(biter-bold);
// modulate by the probability of nuclear breakup as f(biter)
// cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl;
if (_beamBreakupMode > 1){
fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
}
// cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl;
integratedflux=integratedflux+fluxelement;
} //end loop over impact parameter
} //end of else (pp, pA, AA)
// In lookup table, store k*dN/dk because it changes less
// so the interpolation should be better
dide[j]=integratedflux*energy;
}//end loop over photon energy
// for 2nd and subsequent calls, use lookup table immediately
L1000f:
lEgamma=log(Egamma);
if (lEgamma < (lnEmin+dlnE) || lEgamma > lnEmax){
flux_r=0.0;
// cout<<" WARNING: Egamma outside defined range. Egamma= "<<Egamma
// <<" "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl;
}
else{
// >> Egamma between Ilt and Ilt+1
Ilt = int((lEgamma-lnEmin)/dlnE);
// >> ln(Egamma) for first point
lnElt = lnEmin + Ilt*dlnE;
// >> Interpolate
flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
flux_r = flux_r/Egamma;
}
return flux_r;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::nepoint(const double Egamma,
const double bmin)
{
// Function for the spectrum of virtual photons,
// dn/dEgamma, for a point charge q=Ze sweeping
// past the origin with velocity gamma
// (=1/SQRT(1-(V/c)**2)) integrated over impact
// parameter from bmin to infinity
// See Jackson eq15.54 Classical Electrodynamics
// Declare Local Variables
double beta,X,C1,bracket,nepoint_r;
beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma)));
X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc);
bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X)
-bessel::dbesk0(X)*bessel::dbesk0(X));
bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X);
// Note: NO Z*Z!!
C1=(2.*alpha)/pi;
nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket;
return nepoint_r;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::sigmagp(const double Wgp)
{
// Function for the gamma-proton --> VectorMeson
// cross section. Wgp is the gamma-proton CM energy.
// Unit for cross section: fm**2
double sigmagp_r=0.;
switch(_particleType)
{
case RHO:
+ case RHO_ee:
+ case RHO_mumu:
case RHOZEUS:
case FOURPRONG:
sigmagp_r=1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
break;
case OMEGA:
sigmagp_r=1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
break;
case PHI:
sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
break;
case JPSI:
case JPSI_ee:
case JPSI_mumu:
case JPSI_ppbar:
sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
// sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
break;
case JPSI2S:
case JPSI2S_ee:
case JPSI2S_mumu:
sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
sigmagp_r*=0.166;
// sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
break;
case UPSILON:
case UPSILON_ee:
case UPSILON_mumu:
// >> This is W**1.7 dependence from QCD calculations
// sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-10*6.4*exp(0.74*log(Wgp));
break;
case UPSILON2S:
case UPSILON2S_ee:
case UPSILON2S_mumu:
// sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-10*2.9*exp(0.74*log(Wgp));
break;
case UPSILON3S:
case UPSILON3S_ee:
case UPSILON3S_mumu:
// sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-10*2.1*exp(0.74*log(Wgp));
break;
default: cout<< "!!! ERROR: Unidentified Vector Meson: "<< _particleType <<endl;
}
return sigmagp_r;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::sigma_A(const double sig_N, const int beam)
{
// Nuclear Cross Section
// sig_N,sigma_A in (fm**2)
double sum;
double b,bmax,Pint,arg,sigma_A_r;
int NGAUSS;
double xg[17]=
{.0,
.0483076656877383162,.144471961582796493,
.239287362252137075, .331868602282127650,
.421351276130635345, .506899908932229390,
.587715757240762329, .663044266930215201,
.732182118740289680, .794483795967942407,
.849367613732569970, .896321155766052124,
.934906075937739689, .964762255587506430,
.985611511545268335, .997263861849481564
};
double ag[17]=
{.0,
.0965400885147278006, .0956387200792748594,
.0938443990808045654, .0911738786957638847,
.0876520930044038111, .0833119242269467552,
.0781938957870703065, .0723457941088485062,
.0658222227763618468, .0586840934785355471,
.0509980592623761762, .0428358980222266807,
.0342738629130214331, .0253920653092620595,
.0162743947309056706, .00701861000947009660
};
NGAUSS=16;
// Check if one or both beams are nuclei
int A_1 = _bbs.beam1().A();
int A_2 = _bbs.beam2().A();
if( A_1 == 1 && A_2 == 1)cout<<" This is pp, you should not be here..."<<endl;
// CALCULATE P(int) FOR b=0.0 - bmax (fm)
bmax = 25.0;
sum = 0.;
for(int IB=1;IB<=NGAUSS;IB++){
b = 0.5*bmax*xg[IB]+0.5*bmax;
if( A_1 == 1 && A_2 != 1){
arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
}else if(A_2 == 1 && A_1 != 1){
arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
}else{
// Check which beam is target
if( beam == 1 ){
arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
}else if( beam==2 ){
arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
}else{
cout<<" Something went wrong here, beam= "<<beam<<endl;
}
}
Pint=1.0-exp(arg);
// If this is a quantum Glauber calculation, use the quantum Glauber formula
if (_quantumGlauber == 1){Pint=2.0*(1.0-exp(arg/2.0));}
sum=sum+2.*pi*b*Pint*ag[IB];
b = 0.5*bmax*(-xg[IB])+0.5*bmax;
if( A_1 == 1 && A_2 != 1){
arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
}else if(A_2 == 1 && A_1 != 1){
arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
}else{
// Check which beam is target
if( beam == 1 ){
arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
}else if(beam==2){
arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
}else{
cout<<" Something went wrong here, beam= "<<beam<<endl;
}
}
Pint=1.0-exp(arg);
// If this is a quantum Glauber calculation, use the quantum Glauber formula
if (_quantumGlauber == 1){Pint=2.0*(1.0-exp(arg/2.0));}
sum=sum+2.*pi*b*Pint*ag[IB];
}
sum=0.5*bmax*sum;
sigma_A_r=sum;
return sigma_A_r;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::sigma_N(const double Wgp)
{
// Nucleon Cross Section in (fm**2)
double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
return cs;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::breitWigner(const double W,
const double C)
{
// use simple fixed-width s-wave Breit-Wigner without coherent background for rho'
// (PDG '08 eq. 38.56)
if(_particleType==FOURPRONG) {
if (W < 4.01 * _ip->pionChargedMass())
return 0;
const double termA = _channelMass * _width;
const double termA2 = termA * termA;
const double termB = W * W - _channelMass * _channelMass;
double C = _ANORM * _ANORM * termA2 / (termB * termB + termA2);
return C/W; // this is dsigma/dW, not dsigma/ds need to convert?
}
// Relativistic Breit-Wigner according to J.D. Jackson,
// Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength
// of the resonant term and b the strength of the non-resonant
// term. C is an overall normalization.
double ppi=0.,ppi0=0.,GammaPrim,rat;
double aa,bb,cc;
double nrbw_r;
// width depends on energy - Jackson Eq. A.2
// if below threshold, then return 0. Added 5/3/2001 SRK
// 0.5% extra added for safety margin
// omega added here 10/26/2014 SRK
if( _particleType==RHO ||_particleType==RHOZEUS || _particleType==OMEGA){
if (W < 2.01*_ip->pionChargedMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt( ((W/2.)*(W/2.)) - _ip->pionChargedMass() * _ip->pionChargedMass());
ppi0=0.358;
}
+
+ // handle rho-->e+e- properly
+ if (_particleType==RHO_ee){
+ if(W<2.*_ip->mel()){
+ nrbw_r=0.;
+ return nrbw_r;
+ }
+ ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
+ ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
+ }
+
+ // handle rho-->mu+mu- properly
+ if (_particleType==RHO_mumu){
+ if(W<2.*_ip->muonMass()){
+ nrbw_r=0.;
+ return nrbw_r;
+ }
+ ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
+ ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
+ }
// handle phi-->K+K- properly
if (_particleType == PHI){
if (W < 2.*_ip->kaonChargedMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt( ((W/2.)*(W/2.))- _ip->kaonChargedMass()*_ip->kaonChargedMass());
ppi0=sqrt( ((_channelMass/2.)*(_channelMass/2.))-_ip->kaonChargedMass()*_ip->kaonChargedMass());
}
//handle J/Psi-->e+e- properly
if (_particleType==JPSI || _particleType==JPSI2S){
if(W<2.*_ip->mel()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
}
if (_particleType==JPSI_ee){
if(W<2.*_ip->mel()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
}
if (_particleType==JPSI_mumu){
if(W<2.*_ip->muonMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
}
if (_particleType==JPSI_ppbar){
if(W<2.*_ip->protonMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->protonMass()*_ip->protonMass());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->protonMass()*_ip->protonMass());
}
if (_particleType==JPSI2S_ee){
if(W<2.*_ip->mel()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
}
if (_particleType==JPSI2S_mumu){
if(W<2.*_ip->muonMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
}
if(_particleType==UPSILON || _particleType==UPSILON2S ||_particleType==UPSILON3S ){
if (W<2.*_ip->muonMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
}
if(_particleType==UPSILON_mumu || _particleType==UPSILON2S_mumu ||_particleType==UPSILON3S_mumu ){
if (W<2.*_ip->muonMass()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
}
if(_particleType==UPSILON_ee || _particleType==UPSILON2S_ee ||_particleType==UPSILON3S_ee ){
if (W<2.*_ip->mel()){
nrbw_r=0.;
return nrbw_r;
}
ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
}
if(ppi==0.&&ppi0==0.)
cout<<"Improper Gammaacrosssection::breitwigner, ppi&ppi0=0."<<endl;
rat=ppi/ppi0;
GammaPrim=_width*(_channelMass/W)*rat*rat*rat;
aa=_ANORM*sqrt(GammaPrim*_channelMass*W);
bb=W*W-_channelMass*_channelMass;
cc=_channelMass*GammaPrim;
// First real part squared
nrbw_r = (( (aa*bb)/(bb*bb+cc*cc) + _BNORM)*( (aa*bb)/(bb*bb+cc*cc) + _BNORM));
// Then imaginary part squared
nrbw_r = nrbw_r + (( (aa*cc)/(bb*bb+cc*cc) )*( (aa*cc)/(bb*bb+cc*cc) ));
// Alternative, a simple, no-background BW, following J. Breitweg et al.
// Eq. 15 of Eur. Phys. J. C2, 247 (1998). SRK 11/10/2000
// nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2
nrbw_r = C*nrbw_r;
return nrbw_r;
}
Index: trunk/src/inputParameters.cpp
===================================================================
--- trunk/src/inputParameters.cpp (revision 308)
+++ trunk/src/inputParameters.cpp (revision 309)
@@ -1,844 +1,874 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: $: revision of last commit // $Author:: $: author of last commit
// $Date:: $: 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",0),
_nmbRapidityBins ("RAP_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),
_interferenceStrength ("IF_STRENGTH",0),
_maxPtInterference ("INT_PT_MAX",0),
_nmbPtBinsInterference ("INT_PT_N_BINS",0),
_ptBinWidthInterference("INT_PT_WIDTH",0),
_protonEnergy ("PROTON_ENERGY",0, NOT_REQUIRED),
_minGammaEnergy ("MIN_GAMMA_ENERGY",6.0, NOT_REQUIRED),
_maxGammaEnergy ("MAX_GAMMA_ENERGY",600000.0, NOT_REQUIRED),
_pythiaParams ("PYTHIA_PARAMS","", NOT_REQUIRED),
_pythiaFullEventRecord ("PYTHIA_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),
_printVM ("PRINT_VM",0,NOT_REQUIRED),
_impulseVM ("SELECT_IMPULSE_VM",0,NOT_REQUIRED),
_quantumGlauber ("QUANTUM_GLAUBER",0,NOT_REQUIRED),
_bmin ("BMIN",0,NOT_REQUIRED),
_bmax ("BMAX",0,NOT_REQUIRED),
_outputHeader ("OUTPUT_HEADER", false, NOT_REQUIRED),
_deuteronSlopePar ("deuteronSlopePar" , 9.5 , NOT_REQUIRED),
_protonMass ("protonMass" , 0.938272046 , NOT_REQUIRED),
_pionChargedMass ("pionChargedMass" , 0.13957018 , NOT_REQUIRED),
_pionNeutralMass ("pionNeutralMass" , 0.1349766 , NOT_REQUIRED),
_kaonChargedMass ("kaonChargedMass" , 0.493677 , NOT_REQUIRED),
_mel ("mel" , 0.000510998928, NOT_REQUIRED),
_muonMass ("muonMass" , 0.1056583715 , NOT_REQUIRED),
_tauMass ("tauMass" , 1.77682 , NOT_REQUIRED),
_f0Mass ("f0Mass" , 0.990 , NOT_REQUIRED),
_f0Width ("f0Width" , 0.100 , NOT_REQUIRED),
_f0BrPiPi ("f0BrPiPi" , 1.0 , NOT_REQUIRED),
_etaMass ("etaMass" , 0.547862 , NOT_REQUIRED),
_etaWidth ("etaWidth" , 0.00000131 , NOT_REQUIRED),
_etaPrimeMass ("etaPrimeMass" , 0.95778 , NOT_REQUIRED),
_etaPrimeWidth ("etaPrimeWidth" , 0.000198 , NOT_REQUIRED),
_etaCMass ("etaCMass" , 2.9836 , NOT_REQUIRED),
_etaCWidth ("etaCWidth" , 0.0322 , NOT_REQUIRED),
_f2Mass ("f2Mass" , 1.2751 , NOT_REQUIRED),
_f2Width ("f2Width" , 0.1851 , NOT_REQUIRED),
_f2BrPiPi ("f2BrPiPi" , 0.561 , NOT_REQUIRED),
_a2Mass ("a2Mass" , 1.3183 , NOT_REQUIRED),
_a2Width ("a2Width" , 0.105 , NOT_REQUIRED),
_f2PrimeMass ("f2PrimeMass" , 1.525 , NOT_REQUIRED),
_f2PrimeWidth ("f2PrimeWidth" , 0.073 , NOT_REQUIRED),
_f2PrimeBrKK ("f2PrimeBrKK" , 0.887 , NOT_REQUIRED),
_zoverz03Mass ("zoverz03Mass" , 1.540 , NOT_REQUIRED),
_f0PartialggWidth ("f0PartialggWidth" , 0.29E-6 , NOT_REQUIRED),
_etaPartialggWidth ("etaPartialggWidth" , 0.516E-6 , NOT_REQUIRED),
_etaPrimePartialggWidth("etaPrimePartialggWidth", 4.35E-6 , NOT_REQUIRED),
_etaCPartialggWidth ("etaCPartialggWidth" , 5.0E-6 , NOT_REQUIRED),
_f2PartialggWidth ("f2PartialggWidth" , 3.03E-6 , NOT_REQUIRED),
_a2PartialggWidth ("a2PartialggWidth" , 1.0E-6 , NOT_REQUIRED),
_f2PrimePartialggWidth ("f2PrimePartialggWidth" , 0.081E-6 , NOT_REQUIRED),
_zoverz03PartialggWidth("zoverz03PartialggWidth", 0.1E-6 , NOT_REQUIRED),
_f0Spin ("f0Spin" , 0.0 , NOT_REQUIRED),
_etaSpin ("etaSpin" , 0.0 , NOT_REQUIRED),
_etaPrimeSpin ("etaPrimeSpin" , 0.0 , NOT_REQUIRED),
_etaCSpin ("etaCSpin" , 0.0 , NOT_REQUIRED),
_f2Spin ("f2Spin" , 2.0 , NOT_REQUIRED),
_a2Spin ("a2Spin" , 2.0 , NOT_REQUIRED),
_f2PrimeSpin ("f2PrimeSpin" , 2.0 , NOT_REQUIRED),
_zoverz03Spin ("zoverz03Spin" , 2.0 , NOT_REQUIRED),
_axionSpin ("axionSpin" , 0.0 , NOT_REQUIRED),
_rho0Mass ("rho0Mass" , 0.769 , NOT_REQUIRED),
_rho0Width ("rho0Width" , 0.1517 , NOT_REQUIRED),
_rho0BrPiPi ("rho0BrPiPi" , 1.0 , NOT_REQUIRED),
+ _rho0Bree ("rho0Bree" , 0.0000472 , NOT_REQUIRED), // added from PDGlive (25 Jun 2019)
+ _rho0Brmumu ("rho0Brmumu" , 0.0000455 , NOT_REQUIRED), // added from PDGlive (25 Jun 2019)
_rho0PrimeMass ("rho0PrimeMass" , 1.540 , NOT_REQUIRED),
_rho0PrimeWidth ("rho0PrimeWidth" , 0.570 , NOT_REQUIRED),
_rho0PrimeBrPiPi ("rho0PrimeBrPiPi" , 1.0 , NOT_REQUIRED),
_OmegaMass ("OmegaMass" , 0.78265 , NOT_REQUIRED),
_OmegaWidth ("OmegaWidth" , 0.00849 , NOT_REQUIRED),
_OmegaBrPiPi ("OmegaBrPiPi" , 0.0153 , NOT_REQUIRED),
_PhiMass ("PhiMass" , 1.019461 , NOT_REQUIRED),
_PhiWidth ("PhiWidth" , 0.004266 , NOT_REQUIRED),
_PhiBrKK ("PhiBrKK" , 0.489 , NOT_REQUIRED),
_JpsiMass ("JpsiMass" , 3.096916 , NOT_REQUIRED),
_JpsiWidth ("JpsiWidth" , 0.0000929 , NOT_REQUIRED),
_JpsiBree ("JpsiBree" , 0.05971 , NOT_REQUIRED),
_JpsiBrmumu ("JpsiBrmumu" , 0.05961 , NOT_REQUIRED),
_JpsiBrppbar ("JpsiBrppbar" , 0.002120 , NOT_REQUIRED),
_Psi2SMass ("Psi2SMass" , 3.686109 , NOT_REQUIRED),
_Psi2SWidth ("Psi2SWidth" , 0.000299 , NOT_REQUIRED),
_Psi2SBree ("Psi2SBree" , 0.00789 , NOT_REQUIRED),
_Psi2SBrmumu ("Psi2SBrmumu" , 0.0079 , NOT_REQUIRED),
_Upsilon1SMass ("Upsilon1SMass" , 9.46030 , NOT_REQUIRED),
_Upsilon1SWidth ("Upsilon1SWidth" , 0.00005402 , NOT_REQUIRED),
_Upsilon1SBree ("Upsilon1SBree" , 0.0238 , NOT_REQUIRED),
_Upsilon1SBrmumu ("Upsilon1SBrmumu" , 0.0248 , NOT_REQUIRED),
_Upsilon2SMass ("Upsilon2SMass" , 10.02326 , NOT_REQUIRED),
_Upsilon2SWidth ("Upsilon2SWidth" , 0.00003198 , NOT_REQUIRED),
_Upsilon2SBree ("Upsilon2SBree" , 0.0191 , NOT_REQUIRED),
_Upsilon2SBrmumu ("Upsilon2SBrmumu" , 0.0193 , NOT_REQUIRED),
_Upsilon3SMass ("Upsilon3SMass" , 10.3552 , NOT_REQUIRED),
_Upsilon3SWidth ("Upsilon3SWidth" , 0.00002032 , NOT_REQUIRED),
_Upsilon3SBree ("Upsilon3SBree" , 0.0218 , NOT_REQUIRED),
_Upsilon3SBrmumu ("Upsilon3SBrmumu" , 0.0218 , 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(_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(_pythiaParams);
_ip.addParameter(_pythiaFullEventRecord);
_ip.addParameter(_xsecCalcMethod);
_ip.addParameter(_axionMass); // AXION HACK
_ip.addParameter(_bslopeDefinition);
_ip.addParameter(_bslopeValue);
_ip.addParameter(_printVM);
_ip.addParameter(_impulseVM);
_ip.addParameter(_quantumGlauber);
_ip.addParameter(_bmin);
_ip.addParameter(_bmax);
_ip.addParameter(_outputHeader);
_ip.addParameter(_deuteronSlopePar );
_ip.addParameter(_protonMass );
_ip.addParameter(_pionChargedMass );
_ip.addParameter(_pionNeutralMass );
_ip.addParameter(_kaonChargedMass );
_ip.addParameter(_mel );
_ip.addParameter(_muonMass );
_ip.addParameter(_tauMass );
_ip.addParameter(_f0Mass );
_ip.addParameter(_f0Width );
_ip.addParameter(_f0BrPiPi );
_ip.addParameter(_etaMass );
_ip.addParameter(_etaWidth );
_ip.addParameter(_etaPrimeMass );
_ip.addParameter(_etaPrimeWidth );
_ip.addParameter(_etaCMass );
_ip.addParameter(_etaCWidth );
_ip.addParameter(_f2Mass );
_ip.addParameter(_f2Width );
_ip.addParameter(_f2BrPiPi );
_ip.addParameter(_a2Mass );
_ip.addParameter(_a2Width );
_ip.addParameter(_f2PrimeMass );
_ip.addParameter(_f2PrimeWidth );
_ip.addParameter(_f2PrimeBrKK );
_ip.addParameter(_zoverz03Mass );
_ip.addParameter(_f0PartialggWidth );
_ip.addParameter(_etaPartialggWidth );
_ip.addParameter(_etaPrimePartialggWidth);
_ip.addParameter(_etaCPartialggWidth );
_ip.addParameter(_f2PartialggWidth );
_ip.addParameter(_a2PartialggWidth );
_ip.addParameter(_f2PrimePartialggWidth );
_ip.addParameter(_zoverz03PartialggWidth);
_ip.addParameter(_f0Spin );
_ip.addParameter(_etaSpin );
_ip.addParameter(_etaPrimeSpin );
_ip.addParameter(_etaCSpin );
_ip.addParameter(_f2Spin );
_ip.addParameter(_a2Spin );
_ip.addParameter(_f2PrimeSpin );
_ip.addParameter(_zoverz03Spin );
_ip.addParameter(_axionSpin );
_ip.addParameter(_rho0Mass );
_ip.addParameter(_rho0Width );
_ip.addParameter(_rho0BrPiPi );
+ _ip.addParameter(_rho0Bree );
+ _ip.addParameter(_rho0Brmumu );
_ip.addParameter(_rho0PrimeMass );
_ip.addParameter(_rho0PrimeWidth );
_ip.addParameter(_rho0PrimeBrPiPi );
_ip.addParameter(_OmegaMass );
_ip.addParameter(_OmegaWidth );
_ip.addParameter(_OmegaBrPiPi );
_ip.addParameter(_PhiMass );
_ip.addParameter(_PhiWidth );
_ip.addParameter(_PhiBrKK );
_ip.addParameter(_JpsiMass );
_ip.addParameter(_JpsiWidth );
_ip.addParameter(_JpsiBree );
_ip.addParameter(_JpsiBrmumu );
_ip.addParameter(_JpsiBrppbar );
_ip.addParameter(_Psi2SMass );
_ip.addParameter(_Psi2SWidth );
_ip.addParameter(_Psi2SBree );
_ip.addParameter(_Psi2SBrmumu );
_ip.addParameter(_Upsilon1SMass );
_ip.addParameter(_Upsilon1SWidth );
_ip.addParameter(_Upsilon1SBree );
_ip.addParameter(_Upsilon1SBrmumu );
_ip.addParameter(_Upsilon2SMass );
_ip.addParameter(_Upsilon2SWidth );
_ip.addParameter(_Upsilon2SBree );
_ip.addParameter(_Upsilon2SBrmumu );
_ip.addParameter(_Upsilon3SMass );
_ip.addParameter(_Upsilon3SWidth );
_ip.addParameter(_Upsilon3SBree );
_ip.addParameter(_Upsilon3SBrmumu );
}
//______________________________________________________________________________
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;
}
/// temporary output test
printWarn<< * this<<endl;
return true;
}
bool inputParameters::init()
{
// Calculate beam gamma in CMS frame
double rap1 = acosh(beam1LorentzGamma());
double rap2 = -acosh(beam2LorentzGamma());
_beamLorentzGamma = cosh((rap1-rap2)/2);
std::cout << "Rapidity beam 1: " << rap1 << ", rapidity beam 2: " << rap2 << ", rapidity CMS system: " << (rap1+rap2)/2 << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
_ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
_protonEnergy = _beamLorentzGamma * protonMass();
// check for deuteron or tritium - these must be the second beam
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;}
+ switch(_prodParticleId.value()) {
+ case 113:
+ case 113011:
+ case 113013:
+ if(_productionMode.value()==2){
+ printWarn << endl<< "For rho meson production, you should choose the wide resonance option (production mode = 3)" << endl;
+ return false;
+ }
+ default:
+ break;
+ }
// 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;
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 = 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 = a2Mass();
width = 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 = etaMass();
width = 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 = f2Mass();
width = 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 = f2BrPiPi();
break;
case 331: // eta'(958)
_particleType = ETAPRIME;
_decayType = SINGLEMESON;
mass = etaPrimeMass();
width = 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 = f2PrimeMass();
width = 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 = f2PrimeBrKK();
break;
case 441: // eta_c(1s)
_particleType = ETAC;
_decayType = SINGLEMESON;
mass = etaCMass();
width = 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 = f0Mass();
width = 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 = 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 = rho0Mass();
width = rho0Width();
defaultMinW = 2 * pionChargedMass();
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = rho0BrPiPi();
break;
+ case 113011: // rho(770) -> e+ e-
+ _particleType = RHO_ee;
+ _decayType = WIDEVMDEFAULT;
+ mass = rho0Mass();
+ width = rho0Width();
+ defaultMinW = 2 * pionChargedMass();
+ defaultMaxW = mass + 5 * width;
+ _inputBranchingRatio = rho0Bree();
+ break;
+ case 113013: // rho(770) -> mu+ mu -
+ _particleType = RHO_mumu;
+ _decayType = WIDEVMDEFAULT;
+ mass = rho0Mass();
+ width = rho0Width();
+ defaultMinW = 2 * pionChargedMass();
+ defaultMaxW = mass + 5 * width;
+ _inputBranchingRatio = rho0Brmumu();
+ break;
case 913: // rho(770) with direct pi+pi- decay, interference given by ZEUS data
_particleType = RHOZEUS;
_decayType = WIDEVMDEFAULT;
mass = rho0Mass();
width = rho0Width();
defaultMinW = 2 * pionChargedMass();
defaultMaxW = mass + 5 * width; // use the same 1.5GeV max mass as ZEUS
_inputBranchingRatio = rho0BrPiPi();
break;
case 999: // pi+pi-pi+pi- phase space decay
_particleType = FOURPRONG;
_decayType = WIDEVMDEFAULT;
mass = rho0PrimeMass();
width = 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
if (defaultMaxW>10.){defaultMaxW=10.;} //SRK May 29, 2019, to avoid an unduly high defaultMaxW
_inputBranchingRatio = 1.0;
break;
case 223: // omega(782)
_particleType = OMEGA;
_decayType = NARROWVMDEFAULT;
mass = OmegaMass();
width = OmegaWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = OmegaBrPiPi();
break;
case 333: // phi(1020)
_particleType = PHI;
_decayType = NARROWVMDEFAULT;
mass = PhiMass();
width = PhiWidth();
defaultMinW = 2 * kaonChargedMass();
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = PhiBrKK();
break;
case 443: // J/psi
_particleType = JPSI;
_decayType = NARROWVMDEFAULT;
mass = JpsiMass();
width = JpsiWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (JpsiBree() + JpsiBrmumu())/2.;
break;
case 443011: // J/psi
_particleType = JPSI_ee;
_decayType = NARROWVMDEFAULT;
mass = JpsiMass();
width = JpsiWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = JpsiBree();
break;
case 443013: // J/psi
cout<<"In inputParameters setting J/psi mass!"<<endl;
_particleType = JPSI_mumu;
_decayType = NARROWVMDEFAULT;
mass = JpsiMass();
width = JpsiWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = JpsiBrmumu();
break;
case 4432212: // J/psi
cout<<"In inputParameters setting J/psi mass!"<<endl;
_particleType = JPSI_ppbar;
_decayType = NARROWVMDEFAULT;
mass = JpsiMass();
width = JpsiWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = JpsiBrppbar();
break;
case 444: // psi(2S)
_particleType = JPSI2S;
_decayType = NARROWVMDEFAULT;
mass = Psi2SMass();
width = Psi2SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (Psi2SBree() + Psi2SBrmumu())/2.;
break;
case 444011: // psi(2S)
_particleType = JPSI2S_ee;
_decayType = NARROWVMDEFAULT;
mass = Psi2SMass();
width = Psi2SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Psi2SBree();
break;
case 444013: // psi(2S)
_particleType = JPSI2S_mumu;
_decayType = NARROWVMDEFAULT;
mass = Psi2SMass();
width = Psi2SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Psi2SBrmumu();
break;
case 553: // Upsilon(1S)
_particleType = UPSILON;
_decayType = NARROWVMDEFAULT;
mass = Upsilon1SMass();
width = Upsilon1SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (Upsilon1SBree() + Upsilon1SBrmumu())/2.;
break;
case 553011: // Upsilon
_particleType = UPSILON_ee;
_decayType = NARROWVMDEFAULT;
mass = Upsilon1SMass();
width = Upsilon1SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Upsilon1SBree();
break;
case 553013: // Upsilon
_particleType = UPSILON_mumu;
_decayType = NARROWVMDEFAULT;
mass = Upsilon1SMass();
width = Upsilon1SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Upsilon1SBrmumu();
break;
case 554: // Upsilon(2S)
_particleType = UPSILON2S;
_decayType = NARROWVMDEFAULT;
mass = Upsilon2SMass();
width = Upsilon2SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (Upsilon2SBree() + Upsilon2SBrmumu())/2.;
break;
case 554011: // Upsilon(2S)
_particleType = UPSILON2S_ee;
_decayType = NARROWVMDEFAULT;
mass = Upsilon2SMass();
width = Upsilon2SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Upsilon2SBree();
break;
case 554013: // Upsilon(2S)
_particleType = UPSILON2S_mumu;
_decayType = NARROWVMDEFAULT;
mass = Upsilon2SMass();
width = Upsilon2SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Upsilon2SBrmumu();
break;
case 555: // Upsilon(3S)
_particleType = UPSILON3S;
_decayType = NARROWVMDEFAULT;
mass = Upsilon3SMass();
width = Upsilon3SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = (Upsilon3SBree() + Upsilon3SBrmumu())/2.;
break;
case 555011: // Upsilon(3S)
_particleType = UPSILON3S_ee;
_decayType = NARROWVMDEFAULT;
mass = Upsilon3SMass();
width = Upsilon3SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = Upsilon3SBree();
break;
case 555013: // Upsilon(3S)
_particleType = UPSILON3S_mumu;
_decayType = NARROWVMDEFAULT;
mass = Upsilon3SMass();
width = Upsilon3SWidth();
defaultMinW = mass - 5 * width;
defaultMaxW = mass + 5 * width;
_inputBranchingRatio = 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;
}
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
<< " 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;
}
}
out <<" Quantum Glauber parameter............... "<<_quantumGlauber.value()<<endl;
out <<" Impulse VM parameter.................... "<<_impulseVM.value()<<endl;
// if (_quantumGlauber.value()==1) {out << " Quantum Glauber calculation being used"<< endl;}
//if (_quantumGlauber.value()!=1) {out << " Classical Glauber calculation being used"<< endl;}
if (_beamBreakupMode.value()==8) {
out <<" Minimum impact parameter.................."<<_bmin.value()<<" fm"<<endl;
out <<" Maximum impact parameter.................."<<_bmax.value()<<" fm"<<endl;
}
// Add some checks here SRK September, 2017
if (_beamBreakupMode.value()==8 && _bmin.value() > _bmax.value()) {
out <<"****Error bmin > bmax"<<endl;
exit(-1);
}
if (_beamBreakupMode.value() == 8 && _productionMode.value() !=1) {
out <<"****Error. Cannot use beam breakup mode 8 (with bmin and bmax) with photonuclear interactions"<<endl;
exit(-1);
}
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
<< "QUANTUM_GLAUBER"<<quantumGlauber () <<endl
<< "BMIN" <<bmin ()<<endl
<< "BMAX" <<bmax ()<<endl;
return out;
}
std::string
inputParameters::parameterValueKey() const
{
return parameterbase::_parameters.validationKey();
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:58 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805712
Default Alt Text
(162 KB)

Event Timeline