Index: trunk/src/readinluminosity.cpp
===================================================================
--- trunk/src/readinluminosity.cpp (revision 8)
+++ trunk/src/readinluminosity.cpp (revision 9)
@@ -1,346 +1,346 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: 265 $: revision of last commit
// $Author:: butter $: author of last commit
// $Date:: 2016-06-06 21:37:26 +0100 #$: date of last commit
//
// Description:
// Added 18->19 for reading in the luminosity table
// Incoherent factor added to table --Joey
//
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include "readinluminosity.h"
#include "starlightconstants.h"
#include "inputParameters.h"
using namespace std;
//______________________________________________________________________________
readLuminosity::readLuminosity(const inputParameters& inputParametersInstance)
: _Warray(0), _Yarray(0), _Farray(0), _Farray1(0), _Farray2(0),
_f_WYarray(0), _g_Earray(0)/*, _g_EQ2array(0)*/
, _ReadInputNPT(inputParametersInstance.nmbPtBinsInterference())
, _ReadInputnumy(inputParametersInstance.nmbRapidityBins())
, _ReadInputnumw(inputParametersInstance.nmbWBins())
, _ReadInputnumega(inputParametersInstance.nmbEnergyBins())
, _ReadInputnumQ2(inputParametersInstance.nmbGammaQ2Bins())
, _ReadInputgg_or_gP(inputParametersInstance.productionMode())
, _ReadInputinterferencemode(inputParametersInstance.interferenceEnabled())
, _baseFileName(inputParametersInstance.baseFileName())
{
}
//______________________________________________________________________________
readLuminosity::~readLuminosity()
{
if(_Warray) delete [] _Warray;
if(_Yarray) delete [] _Yarray;
if(_Farray) delete [] _Farray;
if(_Farray1) delete [] _Farray1;
if(_Farray2) delete [] _Farray2;
// For eSTARLIGHT
if(_f_WYarray) delete [] _f_WYarray;
if(_g_Earray) delete [] _g_Earray;
if(_g_EQ2array) delete _g_EQ2array;
}
//______________________________________________________________________________
void readLuminosity::read()
{
if(!_Warray) _Warray = new double[_ReadInputnumw];
if(!_Yarray) _Yarray = new double[_ReadInputnumy];
if(!_Farray)
{
_Farray = new double*[_ReadInputnumw];
for(int i = 0; i < _ReadInputnumw; i++)
{
_Farray[i] = new double[_ReadInputnumy];
}
}
if(!_Farray1)
{
_Farray1 = new double*[_ReadInputnumw];
for(int i = 0; i < _ReadInputnumw; i++)
{
_Farray1[i] = new double[_ReadInputnumy];
}
}
if(!_Farray2)
{
_Farray2 = new double*[_ReadInputnumw];
for(int i = 0; i < _ReadInputnumw; i++)
{
_Farray2[i] = new double[_ReadInputnumy];
}
}
double dummy[17]; //number of lines used to read in input parameters saved to lookup table[slight.txt].
std::string wyFileName;
wyFileName = _baseFileName +".txt";
// cout << "wyFileName being read in" << wyFileName << endl;
double fpart =0.;
double fptsum=0.;
ifstream wylumfile;
_f_max=0.0;
_f_max1=0.0;
_f_max2=0.0;
wylumfile.open(wyFileName.c_str());
for(int i=0;i < 17;i++){
wylumfile >> dummy[i];
}
int A_1 = dummy[1];
int A_2 = dummy[3];
for(int i=0;i<_ReadInputnumw;i++){
wylumfile >> _Warray[i];
}
for(int i=0;i<_ReadInputnumy;i++){
wylumfile >> _Yarray[i];
}
if( (_ReadInputgg_or_gP == 1) || (A_2 == 1 && A_1 != 1) || (A_1 ==1 && A_2 != 1) ){
for(int i=0;i<_ReadInputnumw;i++){
for(int j=0;j<_ReadInputnumy;j++){
wylumfile >> _Farray[i][j];
if( _Farray[i][j] > _f_max ) _f_max=_Farray[i][j];
}
}
//Normalize farray
for(int i=0;i<_ReadInputnumw;i++){
for(int j=0;j<_ReadInputnumy;j++){
_Farray[i][j]=_Farray[i][j]/_f_max;
}
}
} else {
for(int i=0;i<_ReadInputnumw;i++){
for(int j=0;j<_ReadInputnumy;j++){
wylumfile >> _Farray1[i][j];
}
}
for(int i=0;i<_ReadInputnumw;i++){
for(int j=0;j<_ReadInputnumy;j++){
wylumfile >> _Farray2[i][j];
if( _Farray1[i][j] + _Farray2[i][j] > _f_max ) _f_max=(_Farray1[i][j] + _Farray2[i][j]);
}
}
//Normalize farray, farray1, farray2
for(int i=0;i<_ReadInputnumw;i++){
for(int j=0;j<_ReadInputnumy;j++){
_Farray1[i][j]=_Farray1[i][j]/_f_max;
_Farray2[i][j]=_Farray2[i][j]/_f_max;
_Farray[i][j]=_Farray1[i][j]+_Farray2[i][j];
}
}
}
wylumfile >> _bwnormsave;
if (_ReadInputgg_or_gP != 1 && _ReadInputinterferencemode != 0) {
// only numy/2 y bins here, from 0 (not -ymax) to ymax
double **finterm = new double*[starlightLimits::MAXWBINS];
for (int i = 0; i < starlightLimits::MAXWBINS; i++) finterm[i] = new double[starlightLimits::MAXYBINS];
for (int i=0;i<_ReadInputnumy;i++) {
//fmax=0;
//we want to convert _fptarray to an integral array where fpt(i,j) is near 0, and fpt(j,NPT) ~1. This will facilitate a simple table loookup
fptsum=0.;
for (int j=0;j<_ReadInputNPT;j++) {
wylumfile >> fpart;
finterm[i][j] = fpart;
_fptarray[i][j]=0.;
fptsum=fptsum+fpart;
}
//convert array to integral
_fptarray[i][0]=finterm[i][0]/fptsum;
for (int j=1;j<_ReadInputNPT;j++) {
for (int k=0;k<j;k++) {
_fptarray[i][j]=_fptarray[i][j]+finterm[i][k];
}
_fptarray[i][j]=_fptarray[i][j]/fptsum;
}
}
delete [] finterm;
}
wylumfile.close();
return;
}
//______________________________________________________________________________
void readLuminosity::e_read()
{
if(!_Warray) _Warray = new double[_ReadInputnumw];
if(!_Yarray) _Yarray = new double[_ReadInputnumy];
if(!_BWarray) _BWarray = new double[_ReadInputnumw];
if(!_f_WYmax)
{
_f_WYarray = new double*[_ReadInputnumega];
for(int i = 0; i < _ReadInputnumega; i++)
{
_f_WYarray[i] = new double[_ReadInputnumQ2];
}
}
if(!_g_Earray)
{
_g_Earray = new double*[_ReadInputnumega];
for(int i = 0; i < _ReadInputnumega; i++)
{
_g_Earray[i] = new double[_ReadInputnumQ2];
}
}
double dummy[13]; //number of lines used to read in input parameters saved to lookup table[slight.txt].
std::string wyFileName;
wyFileName = _baseFileName +".txt";
// cout << "wyFileName being read in" << wyFileName << endl;
ifstream wylumfile;
_f_WYmax=0.0;
_g_Emax=0.0;
wylumfile.open(wyFileName.c_str());
for(int i=0;i < 13;i++){
wylumfile >> dummy[i];
}
int A_1 = dummy[1];
int A_2 = dummy[3];
for(int i=0;i<_ReadInputnumw;i++){
wylumfile >> _Warray[i];
}
//No longer needed, can be removed later
for(int i=0;i<_ReadInputnumy;i++){
wylumfile >> _Yarray[i];
}
//New table saving the Breit-Wigner function for form factor
double bw_max = 0 ;
for(int i=0;i<_ReadInputnumw;i++){
wylumfile >> _BWarray[i];
if( _BWarray[i] > bw_max )
bw_max = _BWarray[i];
}
for(int i=0;i<_ReadInputnumw;i++){
_BWarray[i] = _BWarray[i]/bw_max;
}
if( (A_2 == 0 && A_1 >= 1) || (A_1 ==0 && A_2 >= 1) ){
for(int i=0;i<_ReadInputnumega;i++){
for(int j=0;j<_ReadInputnumQ2;j++){
wylumfile >> _f_WYarray[i][j];
if( _f_WYarray[i][j] > _f_WYmax ) _f_WYmax=_f_WYarray[i][j];
//
//wylumfile >> _g_Earray[i][j];
//if( _g_Earray[i][j] > _g_Emax ) _g_Emax = _g_Earray[i][j];
}
}
//Normalize f_WY array, g does not need to be normalized, it is used for normalization
for(int i=0;i<_ReadInputnumega;i++){
for(int j=0;j<_ReadInputnumQ2;j++){
_f_WYarray[i][j] = _f_WYarray[i][j]/( _f_WYmax );
//_g_Earray[i][j] = _g_Earray[i][j]/_g_Emax;
}
}
}
wylumfile >> _bwnormsave;
wylumfile.close();
cout<<"Done reading wylumi file"<<endl;
std::string EQ2FileName;
EQ2FileName = "e_"+_baseFileName+".txt";
ifstream EQlumfile;
EQlumfile.open(EQ2FileName.c_str());
int n_steps;
EQlumfile >> n_steps;
double integrated_max = 0 ;
// _g_EQ2array = new map<string,std::vector<double> >();
_g_EQ2array = new vector< std::pair< double, vector<double> > > ();
while( !EQlumfile.eof() ){
_g_EQ2max = 0 ;
double integral;
std::vector<double> p;
for( int iQ2 = 0 ; iQ2 < _ReadInputnumQ2+3 ; ++iQ2){
if(iQ2 == 0 ){
EQlumfile >> integral;
if( integral > integrated_max)
integrated_max = integral;
}
else{
double temp;
EQlumfile >> temp;
p.push_back(temp);
//cout<<p[iQ2-1]<<endl;
}
if( iQ2 > 2 && p[iQ2-1] > _g_EQ2max){
_g_EQ2max = p[iQ2-1];
}
}
- for( uint iQ2=2; iQ2 < p.size(); ++iQ2)
+ for( unsigned int iQ2=2; iQ2 < p.size(); ++iQ2)
p[iQ2] = p[iQ2]/_g_EQ2max;
_g_EQ2array->push_back(std::pair< double, std::vector<double> >(integral,p));
}
for( std::vector< std::pair<double, std::vector<double> > >::iterator it =_g_EQ2array->begin() ; it != _g_EQ2array->end(); ++it){
//cout<<it->first<<endl;
it->first = it->first/integrated_max;
//cout<<"Now"<<it->first<<endl;
}
//normalize
/*for( std::map<string,std::vector<double> >::iterator it =_g_EQ2array->begin() ; it != _g_EQ2array->end(); ++it){
std::string key = it->first;
if(key=="")
continue;
- for( uint iQ2=2; iQ2 < it->second.size(); ++iQ2)
+ for(unsigned int iQ2=2; iQ2 < it->second.size(); ++iQ2)
it->second[iQ2] = it->second[iQ2]/_g_EQ2max;
}*/
EQlumfile.close();
return;
}
Index: trunk/src/eventfilewriter.cpp
===================================================================
--- trunk/src/eventfilewriter.cpp (revision 8)
+++ trunk/src/eventfilewriter.cpp (revision 9)
@@ -1,189 +1,189 @@
///////////////////////////////////////////////////////////////////////////
//
// Copyright 2010
//
// This file is part of starlight.
//
// starlight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// starlight is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 228 $: revision of last commit
// $Author:: butter $: author of last commit
// $Date:: 2016-01-18 17:10:17 +0000 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#include "eventfilewriter.h"
#include "starlightparticlecodes.h"
//______________________________________________________________________________
eventFileWriter::eventFileWriter()
: fileWriter()
,_writeFullPythia(false)
{ }
//______________________________________________________________________________
eventFileWriter::eventFileWriter(std::string filename)
: fileWriter(filename)
{ }
//______________________________________________________________________________
int eventFileWriter::writeInit(inputParameters &_p)
{
_fileStream<<"CONFIG_OPT: "<<_p.productionMode()<<" "<<_p.prodParticleId()<<" "<<_p.nmbEvents()
<<" "<<_p.quantumGlauber()<<" "<<_p.impulseVM()<<" "<<_p.randomSeed()<<std::endl;
_fileStream<<"BEAM_1: "<<_p.beam1Z()<<" "<<_p.beam1A()<<" "<<_p.beam1LorentzGamma()<<std::endl;
_fileStream<<"BEAM_2: "<<_p.beam2Z()<<" "<<_p.beam2A()<<" "<<_p.beam2LorentzGamma()<<std::endl;
_fileStream<<"PHOTON: "<<_p.nmbEnergyBins()<<" "<<_p.fixedQ2Range()<<" "<<_p.nmbGammaQ2Bins()
<<" "<<_p.minGammaQ2()<<" "<<_p.maxGammaQ2()<<std::endl;
return 0;
}
//______________________________________________________________________________
int eventFileWriter::writeEvent(upcEvent &event, int eventnumber)
{
int numberoftracks = 0;
if(_writeFullPythia)
{
numberoftracks = event.getParticles()->size();
}
else
{
for(unsigned int i = 0; i<event.getParticles()->size(); ++i)
{
if(event.getParticles()->at(i).getStatus() >= 0) numberoftracks++;
}
}
// sometimes we don't have tracks due to wrongly picked W , check it
if(numberoftracks){
eventnumber++;
_fileStream << "EVENT: " << eventnumber << " " << numberoftracks << " " << 1 << std::endl;
if(event.getGammaEnergies()->size()) _fileStream << "GAMMAENERGIES:";
for(unsigned int n = 0; n < event.getGammaEnergies()->size(); n++)
{
_fileStream << " " << event.getGammaEnergies()->at(n);
}
if(event.getGammaEnergies()->size()) _fileStream<< std::endl;
_fileStream <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<std::endl;
int ipart = 0;
std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
{
if(!_writeFullPythia)
{
if((*part).getStatus() < 0) continue;
}
_fileStream << "TRACK: " << " " << starlightParticleCodes::jetsetToGeant((*part).getPdgCode()) <<" "<< (*part).GetPx() << " " << (*part).GetPy()
<< " "<< (*part).GetPz() << " " << eventnumber << " " << ipart << " " << 0 << " "
<< (*part).getPdgCode();
if(_writeFullPythia)
{
lorentzVector vtx = (*part).getVertex();
_fileStream << " " << vtx.GetPx() << " " << vtx.GetPy() << " " << vtx.GetPz() << " " << vtx.GetE();
_fileStream << " " << (*part).getFirstParent() << " " << (*part).getLastParent() << " " << (*part).getFirstDaughter() << " " << (*part).getLastDaughter() << " " << (*part).getStatus();
}
_fileStream <<std::endl;
}
}
return 0;
}
//______________________________________________________________________________
int eventFileWriter::writeEvent(eXEvent &event, int eventnumber)
{
int numberoftracks = 0;
if(_writeFullPythia)
{
numberoftracks = event.getParticles()->size();
}
else
{
for(unsigned int i = 0; i<event.getParticles()->size(); ++i)
{
if(event.getParticles()->at(i).getStatus() >= 0) numberoftracks++;
}
}
// sometimes we don't have tracks due to wrongly picked W , check it
if(numberoftracks){
eventnumber++;
_fileStream << "EVENT: " << eventnumber << " " << numberoftracks << " " << 1 << std::endl;
_fileStream <<"VERTEX: "<<0.<<" "<<0.<<" "<<0.<<" "<<0.<<" "<<1<<" "<<0<<" "<<0<<" "<<numberoftracks<<std::endl;
- for( uint igam = 0 ; igam < event.getGammaEnergies()->size(); ++igam){
+ for(unsigned int igam = 0 ; igam < event.getGammaEnergies()->size(); ++igam){
_fileStream <<"GAMMA: "<<event.getGammaEnergies()->at(igam)<<" "<<event.getGammaMasses()->at(igam)<<std::endl;
lorentzVector gam = event.getGamma()->at(igam);
// Write the photon 4-vector out to file. Might be used in later iterations, so I keep it here
//_fileStream <<"GAMMA VECTOR: "<<gam.GetPx()<<" "<<gam.GetPy()<<" "<<gam.GetPz()<<" "<<gam.GetE()<<" "<<-temp<<std::endl;
}
- for( uint itarget = 0 ; itarget < event.getTarget()->size(); ++itarget){
+ for(unsigned int itarget = 0 ; itarget < event.getTarget()->size(); ++itarget){
lorentzVector target = event.getTarget()->at(itarget);
_fileStream <<"t: "<<event.getVertext()->at(itarget)<<std::endl;
_fileStream <<"TARGET: "<<target.GetPx()<<" "<<target.GetPy()<<" "<<target.GetPz()<<" "<<target.GetE()<<std::endl;
}
- for( uint iel = 0 ; iel < event.getSources()->size(); ++iel){
+ for(unsigned int iel = 0 ; iel < event.getSources()->size(); ++iel){
lorentzVector el = event.getSources()->at(iel);
_fileStream <<"SOURCE: "<<el.GetPx()<<" "<<el.GetPy()<<" "<<el.GetPz()<<" "<<el.GetE()<<std::endl;
}
int ipart = 0;
std::vector<starlightParticle>::const_iterator part = (event.getParticles())->begin();
for (part = event.getParticles()->begin(); part != event.getParticles()->end(); part++, ipart++)
{
if(!_writeFullPythia)
{
if((*part).getStatus() < 0) continue;
}
_fileStream << "TRACK: " << " "<< starlightParticleCodes::jetsetToGeant((*part).getPdgCode()) <<" "<< (*part).GetPx() << " " << (*part).GetPy()
<< " "<< (*part).GetPz() << " " << eventnumber << " " << ipart << " " << 0 << " "
<< (*part).getPdgCode();
if(_writeFullPythia)
{
lorentzVector vtx = (*part).getVertex();
_fileStream << " " << vtx.GetPx() << " " << vtx.GetPy() << " " << vtx.GetPz() << " " << vtx.GetE();
_fileStream << " " << (*part).getFirstParent() << " " << (*part).getLastParent() << " " << (*part).getFirstDaughter() << " " << (*part).getLastDaughter() << " " << (*part).getStatus();
}
_fileStream <<std::endl;
}
}
return 0;
}
Index: trunk/src/gammaeluminosity.cpp
===================================================================
--- trunk/src/gammaeluminosity.cpp (revision 8)
+++ trunk/src/gammaeluminosity.cpp (revision 9)
@@ -1,209 +1,209 @@
///////////////////////////////////////////////////////////////////////////
//
// 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:: 264 $: revision of last commit
// $Author:: jseger $: author of last commit
// $Date:: 2016-06-06 21:05:12 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <fstream>
#include <cmath>
#include "inputParameters.h"
#include "beambeamsystem.h"
#include "beam.h"
#include "starlightconstants.h"
#include "nucleus.h"
#include "bessel.h"
#include "gammaeluminosity.h"
using namespace std;
using namespace starlightConstants;
//______________________________________________________________________________
photonElectronLuminosity::photonElectronLuminosity(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem)
: photonNucleusCrossSection(inputParametersInstance, bbsystem)
,_protonEnergy(inputParametersInstance.protonEnergy())
,_electronEnergy(inputParametersInstance.electronEnergy())
,_beamLorentzGamma(inputParametersInstance.beamLorentzGamma())
,_baseFileName(inputParametersInstance.baseFileName())
,_maxW(inputParametersInstance.maxW())
,_minW(inputParametersInstance.minW())
,_nmbWBins(inputParametersInstance.nmbWBins())
,_maxRapidity(inputParametersInstance.maxRapidity())
,_nmbRapidityBins(inputParametersInstance.nmbRapidityBins())
,_nEBins(inputParametersInstance.nmbEnergyBins())
,_minGammaQ2(inputParametersInstance.minGammaQ2())
,_maxGammaQ2(inputParametersInstance.maxGammaQ2())
,_nmbGammaQ2Bins(inputParametersInstance.nmbGammaQ2Bins())
,_cmsMaxPhotonEnergy(inputParametersInstance.cmsMaxPhotonEnergy())
,_cmsMinPhotonEnergy(inputParametersInstance.cmsMinPhotonEnergy())
,_targetMaxPhotonEnergy(inputParametersInstance.targetMaxPhotonEnergy())
,_targetMinPhotonEnergy(inputParametersInstance.targetMinPhotonEnergy())
,_productionMode(inputParametersInstance.productionMode())
,_beamBreakupMode(inputParametersInstance.beamBreakupMode())
{
cout <<"Creating Luminosity Tables."<<endl;
photonNucleusDifferentialLuminosity();
cout <<"Luminosity Tables created."<<endl;
}
//______________________________________________________________________________
photonElectronLuminosity::~photonElectronLuminosity()
{ }
//______________________________________________________________________________
void photonElectronLuminosity::photonNucleusDifferentialLuminosity()
{
double W,dW, dY;
double Egamma,Y;
//
double testint;
//
double g_E;
double csgA;
double C;
int beam;
//
//int nQ2steps =100;
std::string wyFileName;
wyFileName = _baseFileName +".txt";
ofstream wylumfile;
wylumfile.precision(15);
std::string EQ2FileName;
EQ2FileName = "e_"+_baseFileName +".txt";
ofstream EQ2lumfile;
EQ2lumfile.precision(15);
double bwnorm;
dW = (_wMax-_wMin)/_nWbins;
dY = (_yMax-(-1.0)*_yMax)/_nYbins;
// Look-up table storing double and single photon flux differentials
EQ2lumfile.open(EQ2FileName.c_str());
EQ2lumfile << _nmbGammaQ2Bins <<endl;
// Look-up table storing photo-nuclear cross section
// Write the values of W used in the calculation to slight.txt.
wylumfile.open(wyFileName.c_str());
wylumfile << getbbs().beam1().Z() <<endl;
wylumfile << getbbs().beam1().A() <<endl;
wylumfile << getbbs().beam2().Z() <<endl;
wylumfile << getbbs().beam2().A() <<endl;
wylumfile << _beamLorentzGamma <<endl;
wylumfile << _maxW <<endl;
wylumfile << _minW <<endl;
wylumfile << _nmbWBins <<endl;
wylumfile << _maxRapidity <<endl;
wylumfile << _nmbRapidityBins <<endl;
wylumfile << _productionMode <<endl;
wylumfile << _beamBreakupMode <<endl;
wylumfile << starlightConstants::deuteronSlopePar <<endl;
// Normalize the Breit-Wigner Distribution and write values of W to slight.txt
testint=0.0;
// Grabbing default value for C in the breit-wigner calculation
C=getDefaultC();
for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
W = _wMin + double(i)*dW + 0.5*dW;
testint = testint + breitWigner(W,C)*dW;
wylumfile << W << endl;
}
bwnorm = 1./testint;
// Write the values of Y used in the calculation to slight.txt.
for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
wylumfile << Y << endl;
}
int A_1 = getbbs().beam1().A();
int A_2 = getbbs().beam2().A();
if( A_2 == 0 && A_1 != 0 ){
beam = 1;
} else if( A_1 ==0 && A_2 != 0){
beam = 2;
} else{
beam = 2;
}
// Exponential steps for both tables: Target frame for photon generation and CMS frame for final state generation
double dE_target = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/_nEBins;
// - - - - Calculations in Target frame
for(unsigned int i = 0; i < _nWbins; ++i) {
W = _wMin + double(i)*dW + 0.5*dW;
wylumfile << breitWigner(W,bwnorm) <<endl;
}
for(int j = 0; j < _nEBins ; ++j) {
// Used for calculation of g(Egamma) which must be done in the ion (target) rest frame
Egamma = std::exp( std::log(_targetMinPhotonEnergy)+j*dE_target);
g_E = 0;
// Photon energy limits are determnined in target frame - multiply by Egamma to speed up event generation
std::pair< double, double >* this_energy = Q2arraylimits(Egamma);
double Q2min = this_energy->first;
double Q2max = this_energy->second;
if( Q2min > Q2max)
continue;
//Accounts for the phase space factor
g_E = Egamma*integrated_Q2_dep(Egamma, Q2min, Q2max);
EQ2lumfile << g_E<<endl;
// Write out Q2 range used for generation
EQ2lumfile << Q2min << endl;
EQ2lumfile << Q2max << endl;
- for( uint iQ2 =0 ;iQ2<_nmbGammaQ2Bins; ++iQ2){
+ for( unsigned int iQ2 =0 ;iQ2<_nmbGammaQ2Bins; ++iQ2){
double Q2 = std::exp( std::log(Q2min)+iQ2*std::log(Q2max/Q2min)/100 );
EQ2lumfile<< g(Egamma,Q2) <<endl;
csgA=getcsgA(Egamma,Q2,beam);
wylumfile << csgA << endl;
}
}
EQ2lumfile.close();
wylumfile << bwnorm << endl;
wylumfile.close();
}
string photonElectronLuminosity::gammaTableParse(int ii, int jj)
{
ostringstream tag1, tag2;
tag1<<ii;
tag2<<jj;
string to_ret = tag1.str()+","+tag2.str();
return to_ret;
}