Page MenuHomeHEPForge
Files F8310061

No OneTemporary
Actions

View Options

Index: trunk/include/spectrum.h
===================================================================
--- trunk/include/spectrum.h (revision 293)
+++ trunk/include/spectrum.h (revision 294)
@@ -1,158 +1,158 @@
/*
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
This program 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.
This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SPECTRUM_H
#define SPECTRUM_H
#include <vector>
#include "randomgenerator.h"
class beamBeamSystem;
//class randomGenerator;
class spectrum
{
public:
/** Spectrum must be constructed with beam-beam system, default constructor disallowed */
- spectrum(const randomGenerator &randy, beamBeamSystem *bbs);
+ spectrum(randomGenerator* randy, beamBeamSystem *bbs);
/**
* Generate a table of photon energy probabilities
* Use NK+1 logarithmic steps between Et_min and Eg_max
*/
int generateKsingle();
/**
* Generate a 2-D table of photon energy probabilities
* Use NK+1 x NK+1 logarithmic steps between Et_min and Eg_max
*/
int generateKdouble();
/**
* Get the energy of a single gamma
* @return energy of the gamma
*/
double drawKsingle();
/**
* Get the energy of a single gamma
* @param egamma1 variable passed by reference to get the energy of the frst gamma
* @param egamma2 variable passed by reference to get the energy of the second gamma
* @return energy of the gamma
*/
void drawKdouble(float &egamma1, float &egamma2);
/** Set the beam beam system */
void setBeamBeamSystem(beamBeamSystem *bbs) {
_beamBeamSystem = bbs;
}
/** Set the minimum gamma energy */
void setMinGammaEnergy(double energy) { _eGammaMin = energy; }
/** Set the maximum gamma energy */
void setMaxGammaEnergy(double energy) { _eGammaMax = energy; }
/** Set minimum impact parameter */
void setBmin(double bmin) { _bMin = bmin; }
/** Set maximum impact parameter */
void setBmax(double bmax) { _bMax = bmax; }
protected:
/** Generate the hadron breakup probability table */
virtual bool generateBreakupProbabilities();
/** Needs some explanation */
virtual double getSigma(double /*egamma*/) const {
return 1.05;
}
virtual double getTransformedNofe(double egamma, double b);
/** Minimum impact parameter */
double _bMin;
/** Maximum impact parameter */
double _bMax;
/** Number of bins in impact parameter */
int _nBbins;
/** Vector containing the probability of breakup */
std::vector<double> _probOfBreakup;
/** Beam beam system */
beamBeamSystem *_beamBeamSystem;
private:
double getFnSingle(double egamma) const;
double getFnDouble(double egamma1, double egamma2) const;
/** NK */
int _nK;
/** Contains the 1 photon probabilities */
std::vector<double> _fnSingle;
/** Contains the 2 photon probabilities */
std::vector<std::vector<double> > _fnDouble;
/** Contains the cumulative distribution */
std::vector<double> _fnSingleCumulative;
/** Contains the cumulative distribution */
std::vector<std::vector<double> > _fnDoubleCumulative;
/** */
std::vector<double> _fnDoubleInt;
/** */
std::vector<double> _fnDoubleIntCumulative;
/** Vecotr of gamma energies */
std::vector<double> _eGamma;
/** Min gamma energy */
double _eGammaMin;
/** Max gamma energy */
double _eGammaMax;
/** Z of target */
int _zTarget;
/** A of target */
int _aTarget;
/** Hadron breakup probability is calculated */
bool _hadBreakProbCalculated;
/** random number generator **/
- randomGenerator _randy;
+ randomGenerator* _randy;
/** Default constructed disallowed (not implemented) */
spectrum();
};
#endif // SPECTRUM_H
Index: trunk/include/starlightdpmjet.h
===================================================================
--- trunk/include/starlightdpmjet.h (revision 293)
+++ trunk/include/starlightdpmjet.h (revision 294)
@@ -1,79 +1,79 @@
/*
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
This program 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 2 of the License, or
(at your option) any later version.
This program 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 this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
#ifndef STARLIGHTDPMJET_H
#define STARLIGHTDPMJET_H
#include <eventchannel.h>
#include <spectrum.h>
class starlightDpmJet : public eventChannel
{
public:
- starlightDpmJet(const inputParameters& input,beamBeamSystem& beamsystem);
+ starlightDpmJet(const inputParameters& input, randomGenerator* randy, beamBeamSystem& beamsystem);
int init();
virtual upcEvent produceEvent();
virtual upcEvent produceSingleEvent(int zdirection, float egamma);
virtual upcEvent produceDoubleEvent();
virtual starlightConstants::event produceEvent(int& /*ievent*/) { return starlightConstants::event(); }
void setSingleMode() { _doDoubleEvent = false; }
void setDoubleMode() { _doDoubleEvent = true; }
void setMinGammaEnergy(double energy) { _minGammaEnergy = energy; }
void setMaxGammaEnergy(double energy) { _maxGammaEnergy = energy; }
void setProtonMode(bool v = true) { _protonMode = v; }
private:
/** Contains the photon spectrum */
spectrum *_spectrum;
/** Should we produce a double event? */
bool _doDoubleEvent;
/** Min gamma energy */
double _minGammaEnergy;
/** Max gamma energy */
double _maxGammaEnergy;
/** Proton-nucleus mode */
double _protonMode;
/** Default constructor not implemented */
starlightDpmJet();
};
#endif // STARLIGHTDPMJET_H
Index: trunk/include/spectrumprotonnucleus.h
===================================================================
--- trunk/include/spectrumprotonnucleus.h (revision 293)
+++ trunk/include/spectrumprotonnucleus.h (revision 294)
@@ -1,39 +1,39 @@
/*
<one line to give the library's name and an idea of what it does.>
Copyright (C) 2011 Oystein Djuvsland <oystein.djuvsland@gmail.com>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#ifndef SPECTRUMPROTONLEAD_H
#define SPECTRUMPROTONLEAD_H
#include "spectrum.h"
class beamBeamSystem;
class spectrumProtonNucleus : public spectrum
{
public:
- spectrumProtonNucleus(const randomGenerator &randy, beamBeamSystem *bb);
+ spectrumProtonNucleus(randomGenerator* randy, beamBeamSystem *bb);
virtual double getNucleonNucleonSigma() const { return 7.35; }
protected:
virtual bool generateBreakupProbabilities();
virtual double getSigma(double ) const;
};
#endif // SPECTRUMPROTONLEAD_H
Index: trunk/include/starlightpythia.h
===================================================================
--- trunk/include/starlightpythia.h (revision 293)
+++ trunk/include/starlightpythia.h (revision 294)
@@ -1,92 +1,92 @@
/*
<one line to give the library's name and an idea of what it does.>
Copyright (C) 2011 Oystein Djuvsland <oystein.djuvsland@gmail.com>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#ifndef STARLIGHTPYTHIA_H
#define STARLIGHTPYTHIA_H
#include "upcevent.h"
#include "inputParameters.h"
#include "beambeamsystem.h"
#include "eventchannel.h"
class spectrum;
class starlightPythia : public eventChannel
{
public:
- starlightPythia(const inputParameters& input, beamBeamSystem& bbsystem);
+ starlightPythia(const inputParameters& input,randomGenerator* randy,beamBeamSystem& bbsystem);
virtual ~starlightPythia();
int init(std::string pythiaParams, bool fullEventRecord = false);
virtual upcEvent produceEvent();
virtual upcEvent produceSingleEvent(int /*zdirection*/, float /*egamma*/){return upcEvent();}
virtual upcEvent produceDoubleEvent(){return upcEvent();}
virtual starlightConstants::event produceEvent(int& /*ievent*/){ return starlightConstants::event();}
void setSingleMode() {
_doDoubleEvent = false;
}
void setDoubleMode() {
_doDoubleEvent = true;
}
void setMinGammaEnergy(double energy) {
_minGammaEnergy = energy;
}
void setMaxGammaEnergy(double energy) {
_maxGammaEnergy = energy;
}
void setFullEventRecord(bool fer = true) { _fullEventRecord = fer; }
private:
/** Contains the photon spectrum */
spectrum *_spectrum;
/** Should we produce a double event? */
bool _doDoubleEvent;
/** Min gamma energy */
double _minGammaEnergy;
/** Max gamma energy */
double _maxGammaEnergy;
/** Full event record or not */
bool _fullEventRecord;
/** Prohibited */
starlightPythia();
starlightPythia(const starlightPythia& other);
starlightPythia& operator=(const starlightPythia& other);
bool operator==(const starlightPythia& other) const;
};
#endif // STARLIGHTPYTHIA_H
Index: trunk/src/spectrum.cpp
===================================================================
--- trunk/src/spectrum.cpp (revision 293)
+++ trunk/src/spectrum.cpp (revision 294)
@@ -1,465 +1,465 @@
/*
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
p This program 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.
This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "spectrum.h"
#include <cmath>
#include "beambeamsystem.h"
#include <randomgenerator.h>
#include <iostream>
-spectrum::spectrum(const randomGenerator &randy, beamBeamSystem *bbs) :
+spectrum::spectrum(randomGenerator* randy, beamBeamSystem *bbs) :
_bMin(5.0)
,_bMax(128000.0)
,_nBbins(6400)
,_probOfBreakup(_nBbins)
,_beamBeamSystem(bbs)
,_nK(10000)
,_fnSingle(_nK)
,_fnDouble(_nK)
,_fnSingleCumulative(_nK+1)
,_fnDoubleCumulative(_nK+1)
,_fnDoubleInt(_nK)
,_fnDoubleIntCumulative(_nK+1)
,_eGamma(_nK+1)
,_eGammaMin(6.0)
,_eGammaMax(600000.0)
,_zTarget(82)
,_aTarget(278)
,_hadBreakProbCalculated(false)
,_randy(randy)
{
_eGamma.resize(_nK+1);
_probOfBreakup.resize(_nBbins);
}
int spectrum::generateKsingle()
{
_fnSingle.resize(_nK);
_fnSingleCumulative.resize(_nK+1);
double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
double egamma = _eGammaMin;
for (int i = 0; i < _nK+1; i++)
{
_eGamma[i] = egamma;
egamma = egamma * eg_inc;
}
egamma = _eGammaMin;
double fnorm = 0;
if (_hadBreakProbCalculated == false)
{
_hadBreakProbCalculated = generateBreakupProbabilities();
}
double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
for (int i = 0; i < _nK; i++)
{
double b = _bMin;
double bint = 0.0;
double f1 = 0;
double f2 = 0;
for (int j = 0; j < _nBbins - 1; j++)
{
double bold = b;
if (j == 0)
{
f1 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j]*b;
//std::cout << fProbOfBreakup[j] << std::endl;
}
else
{
f1 = f2;
}
b = b*binc;
f2 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j+1]*b;;
bint = bint + 0.5*(f1+f2)*(b-bold);
}
bint = 2.0*starlightConstants::pi*bint;
if (i == 0)
{
fnorm = 1.0/bint;
}
_fnSingle[i] = bint*(_eGamma[i+1]-_eGamma[i]);
egamma = egamma*eg_inc;
}
_fnSingleCumulative[0] = 0.00;
for (int i = 0; i < _nK; i++)
{
_fnSingleCumulative[i+1] = _fnSingleCumulative[i]+_fnSingle[i];
}
double fnormfactor = 1.00/_fnSingleCumulative[_nK];
for (int i = 0; i < _nK; i++)
{
_fnSingleCumulative[i+1] = fnormfactor*_fnSingleCumulative[i+1];
}
return 0;
}
int spectrum::generateKdouble()
{
//Quick fix for now TODO: Fix it!
_nK = 100;
_fnDouble.resize(_nK);
_fnDoubleInt.resize(_nK);
_fnDoubleIntCumulative.resize(_nK+1);
_fnDoubleCumulative.resize(_nK+1);
for (int i = 0; i < _nK; i++)
{
_fnDouble[i].resize(_nK);
_fnDoubleCumulative[i].resize(_nK+1);
}
_fnDoubleCumulative[_nK].resize(_nK+1);
double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK);
double egamma1 = _eGammaMin;
double egamma2 = _eGammaMin;
for (int i = 0; i < _nK+1; i++)
{
_eGamma[i] = egamma1;
egamma1 = egamma1 * eg_inc;
}
egamma1 = _eGammaMin;
double fnorm = 0;
if (_hadBreakProbCalculated == false)
{
_hadBreakProbCalculated = generateBreakupProbabilities();
}
double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
int nbbins = _nBbins;
for (int i = 0; i < _nK; i++)
{
egamma2 = _eGammaMin;
for (int j = 0; j < _nK; j++)
{
double bint = 0.0;
double b = _bMin;
double f1 = 0;
double f2 = 0;
for (int k = 0; k < nbbins - 1; k++)
{
double bold = b;
if (k == 0)
{
f1 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
* getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k]*b; }
else
{
f1 = f2;
}
b = b*binc;
f2 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b)
* getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k+1]*b;
bint = bint + 0.5*(f1+f2)*(b-bold);
}
bint = 2.0*starlightConstants::pi*bint;
_fnDouble[i][j] = bint * (_eGamma[i+1] - _eGamma[i]) * (_eGamma[j+1] - _eGamma[j]);
egamma2 = egamma2 * eg_inc;
}
egamma1 = egamma1 * eg_inc;
}
for (int i = 0; i < _nK; i++)
{
_fnDoubleInt[i] = 0.0;
for (int j = 0; j < _nK; j++)
{
_fnDoubleInt[i] = _fnDoubleInt[i] + _fnDouble[i][j];
}
}
_fnDoubleIntCumulative[0] = 0.0;
for (int i = 1; i < _nK+1; i++)
{
_fnDoubleIntCumulative[i] = _fnDoubleIntCumulative[i-1] + _fnDoubleInt[i-1];
}
fnorm = 1.0/_fnDoubleIntCumulative[_nK];
for (int i = 0; i < _nK+1; i++)
{
_fnDoubleIntCumulative[i] = fnorm * _fnDoubleIntCumulative[i];
}
return 0;
}
double spectrum::drawKsingle()
{
double xtest = 0;
int itest = 0;
double egamma = 0.0;
- xtest = _randy.Rndom();
+ xtest = _randy->Rndom();
while (xtest > _fnSingleCumulative[itest])
{
itest++;
}
itest = itest - 1;
if (itest >= _nK || itest < 0)
{
std::cerr << "ERROR: itest: " << itest << std::endl;
}
double delta_f = xtest - _fnSingleCumulative[itest];
if (delta_f <= 0.0)
{
std::cout << "WARNING: delta_f: " << delta_f << std::endl;
return -1;
}
double dE = _eGamma[itest+1] - _eGamma[itest];
double dF = _fnSingleCumulative[itest+1] - _fnSingleCumulative[itest];
double delta_e = delta_f*dE/dF;
if (delta_e > (_eGamma[itest+1] - _eGamma[itest]))
{
std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
}
egamma = _eGamma[itest] + delta_e;
return egamma;
}
void spectrum::drawKdouble(float& egamma1, float& egamma2)
{
double xtest1 = 0.0;
double xtest2 = 0.0;
int itest1 = 0;
int itest2 = 0;
- xtest1 = _randy.Rndom();
+ xtest1 = _randy->Rndom();
while (xtest1 > _fnDoubleIntCumulative[itest1])
{
itest1++;
}
itest1 = itest1 - 1;
if (itest1 >= _nK || itest1 < 0)
{
std::cerr << "ERROR: itest1: " << itest1 << std::endl;
}
double delta_f = xtest1 - _fnDoubleIntCumulative[itest1];
if (delta_f <= 0.0)
{
std::cout << "WARNING: delta_f: " << delta_f << std::endl;
}
double dE = _eGamma[itest1+1] - _eGamma[itest1];
double dF = _fnDoubleIntCumulative[itest1+1] - _fnDoubleIntCumulative[itest1];
double delta_e = delta_f*dE/dF;
if (delta_e > (_eGamma[itest1+1] - _eGamma[itest1]))
{
std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
}
egamma1 = _eGamma[itest1] + delta_e;
// Second gamma
std::vector<double> fn_second_cumulative(_nK+1);
fn_second_cumulative[0] = 0.0;
for(int i = 1; i < _nK+1; i++)
{
fn_second_cumulative[i] = fn_second_cumulative[i-1] + _fnDouble[itest1][i-1];
}
double norm_factor = 1.0/fn_second_cumulative[_nK];
for(int i = 0; i < _nK+1; i++)
{
fn_second_cumulative[i] = norm_factor*fn_second_cumulative[i];
}
- xtest2 = _randy.Rndom();
+ xtest2 = _randy->Rndom();
while (xtest2 > fn_second_cumulative[itest2])
{
itest2++;
}
itest2 = itest2 - 1;
if (itest2 >= _nK || itest2 < 0)
{
std::cerr << "ERROR: itest2: " << itest2 << std::endl;
}
delta_f = xtest2 - fn_second_cumulative[itest2];
if (delta_f <= 0.0)
{
std::cout << "WARNING: delta_f: " << delta_f << std::endl;
}
dE = _eGamma[itest2+1] - _eGamma[itest2];
dF = fn_second_cumulative[itest2+1] - fn_second_cumulative[itest2];
delta_e = delta_f*dE/dF;
if (delta_e > (_eGamma[itest2+1] - _eGamma[itest2]))
{
std::cerr << "ERROR: delta_E: " << delta_e << std::endl;
}
egamma2 = _eGamma[itest2] + delta_e;
return;
}
bool spectrum::generateBreakupProbabilities()
{
int nbbins = _nBbins;
double b_min = _bMin;
double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
double b = b_min;
_probOfBreakup.resize(nbbins);
for (int i = 0; i < nbbins; i++)
{
double bimp = b;
double rhad = 0;
rhad = _beamBeamSystem->probabilityOfBreakup(bimp);
_probOfBreakup[i] = exp(-rhad);
b = b*binc;
}
return true;
}
double spectrum::getFnSingle(double egamma) const
{
double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
int i1 = log(egamma/_eGammaMin)/log(eginc);
int i2 = i1 + 1;
double fnSingle = 0.0;
if (i1 < 0 || i2 > _nK)
{
std::cout << "I1, I2 out of bounds. Egamma = " << egamma << std::endl;
std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
fnSingle = 0.0;
}
else
{
double dE = _eGamma[i2] - _eGamma[i1];
double eFrac = (egamma - _eGamma[i1])/dE;
if (eFrac < 0.0 || eFrac > 1.0)
{
std::cout << "WARNING: Efrac = " << eFrac << std::endl;
}
fnSingle = (1.0 - eFrac)*_fnSingle[i1] + eFrac*_fnSingle[i2];
}
return fnSingle;
}
double spectrum::getFnDouble(double egamma1, double egamma2) const
{
double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK));
int i1 = log(egamma1/_eGammaMin)/log(eginc);
int i2 = i1 + 1;
double fnDouble = 0.0;
if (i1 < 0 || i2 > _nK)
{
std::cout << "I1, I2 out of bounds. Egamma1 = " << egamma1 << std::endl;
std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl;
fnDouble = 0.0;
return fnDouble;
}
int j1 = log(egamma2/_eGammaMin)/log(eginc);
int j2 = j1 + 1;
if (j1 < 0 || j2 > _nK)
{
std::cout << "J1, J2 out of bounds. Egamma2 = " << egamma2 << std::endl;
std::cout << "J1, J2 = " << j1 << ", " << j2 << std::endl;
fnDouble = 0.0;
return fnDouble;
}
double dE1 = _eGamma[i2] - _eGamma[i1];
double eFrac1 = (egamma1 - _eGamma[i1])/dE1;
if (eFrac1 < 0.0 || eFrac1 > 1.0)
{
std::cout << "WARNING: Efrac1 = " << eFrac1 << std::endl;
}
double ptemp1 = (1.0 - eFrac1)*_fnDouble[i1][j1] + eFrac1*_fnDouble[i2][j1];
double ptemp2 = (1.0 - eFrac1)*_fnDouble[i1][j2] + eFrac1*_fnDouble[i2][j2];
double dE2 = _eGamma[j2] - _eGamma[j1];
double eFrac2 = (egamma2 - _eGamma[j1])/dE2;
if (eFrac2 < 0.0 || eFrac2 > 1.0)
{
std::cout << "WARNING: Efrac2 = " << eFrac2 << std::endl;
}
fnDouble = (1.0 - eFrac2)*ptemp1 + eFrac2*ptemp2;
return fnDouble;
}
double spectrum::getTransformedNofe(double egamma, double b)
{
double factor = 1.0/(2.0*_beamBeamSystem->beamLorentzGamma());
double res = factor * _beamBeamSystem->beam1().photonDensity(b, egamma*factor);
return res;
}
Index: trunk/src/spectrumprotonnucleus.cpp
===================================================================
--- trunk/src/spectrumprotonnucleus.cpp (revision 293)
+++ trunk/src/spectrumprotonnucleus.cpp (revision 294)
@@ -1,60 +1,60 @@
/*
<one line to give the library's name and an idea of what it does.>
Copyright (C) 2011 Oystein Djuvsland <oystein.djuvsland@gmail.com>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <cmath>
#include "spectrumprotonnucleus.h"
#include "beambeamsystem.h"
#include "beam.h"
#include <iostream>
-spectrumProtonNucleus::spectrumProtonNucleus(const randomGenerator &randy, beamBeamSystem *b) : spectrum(randy,b)
+spectrumProtonNucleus::spectrumProtonNucleus(randomGenerator* randy, beamBeamSystem *b) : spectrum(randy,b)
{
_bMin = 4.0;
}
bool spectrumProtonNucleus::generateBreakupProbabilities()
{
int nbbins = _nBbins;
double b_min = _bMin;
double binc = exp((log(_bMax/_bMin))/(double)_nBbins);
double b = b_min;
_probOfBreakup.resize(nbbins);
for (int i = 0; i < nbbins; i++)
{
_beamBeamSystem->beam1().Z() > 1 ? _probOfBreakup[i] = exp(-getNucleonNucleonSigma()*_beamBeamSystem->beam1().thickness(b)) :
_beamBeamSystem->beam2().Z() > 1 ? _probOfBreakup[i] = exp(-getNucleonNucleonSigma()*_beamBeamSystem->beam2().thickness(b)) :
b < 7.76 ? _probOfBreakup[i] = 0 : _probOfBreakup[i] = 1; // Should always be true though
b = b*binc;
}
return true;
}
double spectrumProtonNucleus::getSigma(double ) const
{
return 0.11;
}
Index: trunk/src/starlightpythia.cpp
===================================================================
--- trunk/src/starlightpythia.cpp (revision 293)
+++ trunk/src/starlightpythia.cpp (revision 294)
@@ -1,130 +1,130 @@
/*
<one line to give the library's name and an idea of what it does.>
Copyright (C) 2011 Oystein Djuvsland <oystein.djuvsland@gmail.com>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <cstdio>
#include <iostream>
#include "starlightpythia.h"
#include "pythiaInterface.h"
#include "spectrumprotonnucleus.h"
#include <cmath>
#include <sstream>
-starlightPythia::starlightPythia(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem) : eventChannel(inputParametersInstance, bbsystem)
+starlightPythia::starlightPythia(const inputParameters& inputParametersInstance,randomGenerator* randy,beamBeamSystem& bbsystem) : eventChannel(inputParametersInstance,randy, bbsystem)
,_spectrum(0)
,_doDoubleEvent(false)
,_minGammaEnergy(inputParametersInstance.minGammaEnergy())
,_maxGammaEnergy(inputParametersInstance.maxGammaEnergy())
,_fullEventRecord(false)
{
}
starlightPythia::~starlightPythia()
{
}
int starlightPythia::init(std::string pythiaParams, bool fullEventRecord)
{
_fullEventRecord = fullEventRecord;
_spectrum = new spectrumProtonNucleus(_randy,&_bbs);
_spectrum->setMinGammaEnergy(_minGammaEnergy);
_spectrum->setMaxGammaEnergy(_maxGammaEnergy);
if(!_doDoubleEvent)
{
_spectrum->generateKsingle();
}
else
{
_spectrum->generateKdouble();
}
// Set to run with varying energies
pythiaInterface::pygive("mstp(171)=1"); // Varying energies
pythiaInterface::pygive("mstp(172)=1"); // Set the energy before generation
std::stringstream ss(pythiaParams);
std::string p;
while(std::getline(ss, p, ';'))
{
if(p.size()>1)
{
pythiaInterface::pygive(p.c_str());
}
}
pythiaInterface::pyinit("FIXT", "gamma", "p", _maxGammaEnergy); // Fixed target, beam, target, beam momentum (GeV/c)
return 0;
}
upcEvent starlightPythia::produceEvent()
{
upcEvent event;
if (!_doDoubleEvent)
{
double gammaE = 0;
do
{
gammaE = _spectrum->drawKsingle();
} while(isnan(gammaE));
event.addGamma(gammaE);
char opt[32];
std::sprintf(opt, "parp(171)=%f", gammaE/_maxGammaEnergy);
pythiaInterface::pygive(opt); // Set the energy of the photon beam (gammaE/1000 * 1000.0);
pythiaInterface::pyevnt(); // Generate event
int zdirection = (_bbs.beam1().Z()==1 ? 1 : -1);
double boost = acosh(_bbs.beamLorentzGamma())*zdirection;
vector3 boostVector(0, 0, tanh(boost));
for(int idx = 0; idx < pyjets_.n; idx++)
{
if(pyjets_.k[0][idx] > 10 && _fullEventRecord==false) continue;
int pdgCode = pyjets_.k[1][idx];
int charge = 0;
if( pdgCode == 12 || pdgCode == 14 || pdgCode == 16 || pdgCode == 22 || pdgCode == 111 || pdgCode == 130 || pdgCode == 321 || pdgCode == 2112)
{
charge = 0;
}
else
{
charge = (pdgCode > 0) - (pdgCode < 0);
}
starlightParticle particle(pyjets_.p[0][idx], pyjets_.p[1][idx], -zdirection*pyjets_.p[2][idx], pyjets_.p[3][idx], pyjets_.p[4][idx], pyjets_.k[1][idx], charge);
if(_fullEventRecord)
{
particle.setFirstDaughter(pyjets_.k[3][idx]);
particle.setLastDaughter(pyjets_.k[4][idx]);
particle.setStatus(pyjets_.k[0][idx]);
}
particle.Boost(boostVector);
event.addParticle(particle);
}
}
return event;
}
Index: trunk/src/starlight.cpp
===================================================================
--- trunk/src/starlight.cpp (revision 293)
+++ trunk/src/starlight.cpp (revision 294)
@@ -1,386 +1,386 @@
///////////////////////////////////////////////////////////////////////////
//
// 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, *_beamSystem);
+ _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, *_beamSystem);
+ _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, *_beamSystem);
+ _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, *_beamSystem);
+ _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 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/starlightdpmjet.cpp
===================================================================
--- trunk/src/starlightdpmjet.cpp (revision 293)
+++ trunk/src/starlightdpmjet.cpp (revision 294)
@@ -1,154 +1,154 @@
/*
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
This program 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 2 of the License, or
(at your option) any later version.
This program 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 this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
#include "starlightdpmjet.h"
#include "spectrum.h"
#include <iostream>
#include <spectrumprotonnucleus.h>
#include "starlightconfig.h"
extern "C"
{
extern struct
{
double slpx[200000];
double slpy[200000];
double slpz[200000];
double sle[200000];
double slm[200000];
int slpid[200000];
int slcharge[200000];
} dpmjetparticle_;
void dt_produceevent_(float* gammaE, int* nparticles);
void dt_getparticle_(int *ipart, int *res);
void dt_initialise_();
}
-starlightDpmJet::starlightDpmJet(const inputParameters& inputParametersInstance, beamBeamSystem& beamsystem ) : eventChannel(inputParametersInstance, beamsystem)
+starlightDpmJet::starlightDpmJet(const inputParameters& inputParametersInstance,randomGenerator* randy,beamBeamSystem& beamsystem ) : eventChannel(inputParametersInstance,randy,beamsystem)
,_spectrum(0)
,_doDoubleEvent(true)
,_minGammaEnergy(6.0)
,_maxGammaEnergy(600000.0)
,_protonMode(false)
{
}
int starlightDpmJet::init()
{
if(_protonMode)
{
_spectrum = new spectrumProtonNucleus(_randy,&_bbs);
}
else
{
_spectrum = new spectrum(_randy,&_bbs);
}
_spectrum->setMinGammaEnergy(_minGammaEnergy);
_spectrum->setMaxGammaEnergy(_maxGammaEnergy);
if(!_doDoubleEvent)
{
_spectrum->generateKsingle();
}
else
{
_spectrum->generateKdouble();
}
return 0;
}
upcEvent starlightDpmJet::produceEvent()
{
upcEvent event;
if (!_doDoubleEvent)
{
int zdirection = 1;
float gammaE = _spectrum->drawKsingle();
event = produceSingleEvent(zdirection, gammaE);
// std::cout << "Gamma energy: " << gammaE << std::endl;
}
else
{
event = produceDoubleEvent();
}
return event;
}
upcEvent starlightDpmJet::produceSingleEvent(int zdirection, float gammaE)
{
upcEvent event;
event.addGamma(gammaE);
int nParticles = 0;
dt_produceevent_(&gammaE, &nParticles);
//In which direction do we go?
double rapidity = _bbs.beam1().rapidity()*zdirection;
for (int i = 0; i < nParticles; i++)
{
starlightParticle particle(dpmjetparticle_.slpx[i], dpmjetparticle_.slpy[i], zdirection*dpmjetparticle_.slpz[i], dpmjetparticle_.sle[i], dpmjetparticle_.slm[i], dpmjetparticle_.slpid[i], dpmjetparticle_.slcharge[i]);
vector3 boostVector(0, 0, tanh(-rapidity));
particle.Boost(boostVector);
event.addParticle(particle);
}
return event;
}
upcEvent starlightDpmJet::produceDoubleEvent()
{
upcEvent event;
float gammaE1 = 0.0;
float gammaE2 = 0.0;
_spectrum->drawKdouble(gammaE1, gammaE2);
// std::cout << "Gamma1 energy: " << gammaE1 << std::endl;
//std::cout << "Gamma2 energy: " << gammaE2 << std::endl;
//In which direction do we go?
- int zdirection = (_randy.Rndom()) < 0.5 ? -1 : 1;
+ int zdirection = (_randy->Rndom()) < 0.5 ? -1 : 1;
event = produceSingleEvent(zdirection, gammaE1);
zdirection = zdirection *-1;
event = event + produceSingleEvent(zdirection, gammaE2);
return event;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:15 PM (12 h, 11 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023602
Default Alt Text
(48 KB)

Event Timeline