No OneTemporary
Actions

Size

48 KB

Subscribers

None

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.0starlightConstants::pibint;
	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.0starlightConstants::pibint;
	_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 + eFrac2ptemp2;

	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 (14 h, 43 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023602
Default Alt Text: (48 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions