Index: trunk/src/e_starlightStandalone.cpp
===================================================================
--- trunk/src/e_starlightStandalone.cpp (revision 10)
+++ trunk/src/e_starlightStandalone.cpp (revision 11)
@@ -1,179 +1,180 @@
//////////////////////////////////////////////////////////////////////////
//
// 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 .
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 270 $: revision of last commit
// $Author:: jnystrand $: author of last commit
// $Date:: 2016-07-08 16:31:51 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
#include
+#include
#include "reportingUtils.h"
#include "e_starlight.h"
#include "inputParameters.h"
#include "eventfilewriter.h"
#include "e_starlightStandalone.h"
using namespace std;
e_starlightStandalone::e_starlightStandalone()
: _configFileName ("slight.in"),
_starlight (0),
_nmbEventsTot (1),
_nmbEventsPerFile (_nmbEventsTot)
{ }
e_starlightStandalone::~e_starlightStandalone()
{ }
bool
e_starlightStandalone::init()
{
_inputParameters = new inputParameters();
// read input parameters from config file
_inputParameters->configureFromFile(_configFileName);
if (!_inputParameters->init()) {
printWarn << "problems initializing input parameters. cannot initialize starlight." << endl;
return false;
}
// copy input file to one with baseFileName naming scheme
std::string inputCopyName, _baseFileName;
_baseFileName = _inputParameters->baseFileName();
if (_baseFileName != "slight") {
inputCopyName = _baseFileName +".in";
ofstream inputCopyFile;
inputCopyFile.open(inputCopyName.c_str());
std::ifstream infile(_configFileName.c_str());
if ((!infile) || (!infile.good()))
{
return -1;
}
int lineSize = 256;
char tmp[lineSize];
while (!infile.getline(tmp, lineSize).eof())
{
cout << tmp << endl;
inputCopyFile << tmp << endl;
}
inputCopyFile.close();
}
// get the number of events
// for now we write everything to one file
_nmbEventsTot = _inputParameters->nmbEvents();
_nmbEventsPerFile = _nmbEventsTot;
// create the starlight object
_starlight = new e_starlight();
// give starlight the input parameters
_starlight->setInputParameters(_inputParameters);
// initialize starlight
return _starlight->init();
}
bool
e_starlightStandalone::run()
{
if (!_starlight) {
printWarn << "null pointer to starlight object. make sure that init() was called. "
<< "cannot generate events." << endl;
return false;
}
// open output file
eventFileWriter fileWriter;
fileWriter.writeFullPythiaInfo(_inputParameters->pythiaFullEventRecord());
_baseFileName = _inputParameters->baseFileName();
_eventDataFileName = _baseFileName +".out";
fileWriter.open(_eventDataFileName);
//
fileWriter.writeInit(*_inputParameters);
//
printInfo << "generating events:" << endl;
unsigned int nmbEvents = 0;
std::chrono::steady_clock::time_point begin= std::chrono::steady_clock::now();
while (nmbEvents < _nmbEventsTot) {
for (unsigned int iEvent = 0; (iEvent < _nmbEventsPerFile) && (nmbEvents < _nmbEventsTot);
++iEvent, ++nmbEvents) {
progressIndicator(iEvent, _nmbEventsTot, true, 4);
eXEvent event = _starlight->produceEvent();
// Boost event from back to lab reference frame
boostEvent(event);
fileWriter.writeEvent(event, iEvent);
}
}
std::chrono::steady_clock::time_point end= std::chrono::steady_clock::now();
float running_total = 1E-3*std::chrono::duration_cast(end - begin).count();
cout<<"Total time "<nmbAttempts() == 0 )return true;
double _branchingRatio = _inputParameters->inputBranchingRatio();
printInfo << "number of attempts = " << _starlight->nmbAttempts() << ", "
<< "number of accepted events = " << _starlight->nmbAccepted() << endl;
double selectedCrossSection =
((double)_starlight->nmbAccepted()/_starlight->nmbAttempts())*_branchingRatio*_starlight->getTotalCrossSection();
if (selectedCrossSection > 1.){
cout<< " The cross section of the generated sample is "< 1.){
cout<< " The cross section of the generated sample is "<<1.E3*selectedCrossSection<<" mb."< 1.){
cout<< " The cross section of the generated sample is "<<1.E6*selectedCrossSection<<" microbarn."< 1.){
cout<< " The cross section of the generated sample is "<<1.E9*selectedCrossSection<<" nanobarn."< 1.){
cout<< " The cross section of the generated sample is "<<1.E12*selectedCrossSection<<" picobarn."<beam1LorentzGamma());
double rap2 = -acosh(_inputParameters->beam2LorentzGamma());
double boost = (rap1+rap2)/2.;
event.boost(boost, rap2); // Boost back to laboratory reference frame. Electron initially in target frame and V.M. in CMS
// Assuming electron is beam1
}
Index: trunk/src/e_wideResonanceCrossSection.cpp
===================================================================
--- trunk/src/e_wideResonanceCrossSection.cpp (revision 10)
+++ trunk/src/e_wideResonanceCrossSection.cpp (revision 11)
@@ -1,351 +1,344 @@
///////////////////////////////////////////////////////////////////////////
//
// 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 .
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 264 $: revision of last commit
// $Author:: mlomnitz $: author of last commit
// $Date:: 2017-03-14 21:05:12 +0100 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////
//#define _makeGammaPQ2_
#include
#include
#include
#include "starlightconstants.h"
#include "e_wideResonanceCrossSection.h"
using namespace std;
using namespace starlightConstants;
//______________________________________________________________________________
e_wideResonanceCrossSection::e_wideResonanceCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
: photonNucleusCrossSection(inputParametersInstance, bbsystem)//hrm
{
_wideWmax = _wMax;
_wideWmin = _wMin;
_targetMaxPhotonEnergy=inputParametersInstance.targetMaxPhotonEnergy();
_targetMinPhotonEnergy=inputParametersInstance.targetMinPhotonEnergy();
_Ep = inputParametersInstance.protonEnergy();
_electronEnergy = inputParametersInstance.electronEnergy();
_target_beamLorentz = inputParametersInstance.beamLorentzGamma();
_VMnumEgamma = inputParametersInstance.nmbEnergyBins();
_useFixedRange = inputParametersInstance.fixedQ2Range();
_gammaMinQ2 = inputParametersInstance.minGammaQ2();
_gammaMaxQ2 = inputParametersInstance.maxGammaQ2();
_targetRadii = inputParametersInstance.targetRadius();
}
//______________________________________________________________________________
e_wideResonanceCrossSection::~e_wideResonanceCrossSection()
{
}
//______________________________________________________________________________
void
e_wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
{
// This subroutine calculates the cross-section assuming a wide
// (Breit-Wigner) resonance.
double W,dW, dEgamma, minEgamma;
double ega[3] = {0};
double int_r,dR;
double int_r2,dR2;
- double Eth;
int iW,nW,iEgamma,nEgamma,beam;
double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
-
- //gamma+nucleon threshold.
- Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
- -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
-
+
// For W integration
nW = 100;
dW = (_wideWmax-_wideWmin)/double(nW);
// For Egamma integration
nEgamma = 1000;
dEgamma = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/nEgamma;
minEgamma = std::log(_targetMinPhotonEnergy);
if (getBNORM() == 0.){
cout<<" Using Breit-Wigner Resonance Profile."<= 1 ){
// eA, first beam is the nucleus and is in this case the target
beam = 1;
} else if( A_1 ==0 && A_2 >= 1){
// eA, second beam is the nucleus and is in this case the target
beam = 2;
} else {
beam = 2;
}
int_r=0.;
int_r2 = 0.;
// Do this first for the case when the first beam is the photon emitter
// The variable beam (=1,2) defines which nucleus is the target
// Integration done using Simpson's rule
for(iW=0;iW<=nW-1;iW++){
W = _wideWmin + double(iW)*dW + 0.5*dW;
int nQ2 = 1000;
for(iEgamma = 0 ; iEgamma < nEgamma; ++iEgamma){ // Integral over photon energy
// Target frame photon energies
ega[0] = exp(minEgamma + iEgamma*dEgamma );
ega[1] = exp(minEgamma + (iEgamma+1)*dEgamma );
ega[2] = 0.5*(ega[0]+ega[1]);
// Integral over Q2
double full_int[3] = {0}; // Full e+X --> e+X+V.M. cross section
double dndE[3] = {0}; // Full e+X --> e+X+V.M. cross section
//
for( int iEgaInt = 0 ; iEgaInt < 3; ++iEgaInt){ // Loop over the energies for the three points to integrate over Q2
//These are the physical limits
double Q2_min = std::pow(starlightConstants::mel*ega[iEgaInt],2.0)/(_electronEnergy*(_electronEnergy-ega[iEgaInt]));
double Q2_max = 4.*_electronEnergy*(_electronEnergy-ega[iEgaInt]);
if(_useFixedRange == true){
if( Q2_min < _gammaMinQ2 )
Q2_min = _gammaMinQ2;
if( Q2_max > _gammaMaxQ2 )
Q2_max = _gammaMaxQ2;
}
double lnQ2ratio = std::log(Q2_max/Q2_min)/nQ2;
double lnQ2_min = std::log(Q2_min);
//
- int q2end = -99;
for( int iQ2 = 0 ; iQ2 < nQ2; ++iQ2){ // Integral over photon virtuality
//
double q2_1 = exp( lnQ2_min + iQ2*lnQ2ratio);
double q2_2 = exp( lnQ2_min + (iQ2+1)*lnQ2ratio);
double q2_12 = (q2_2+q2_1)/2.;
//
// Integrating cross section
full_int[iEgaInt] += (q2_2-q2_1)*( g(ega[iEgaInt],q2_1)*getcsgA(ega[iEgaInt],q2_1,beam)
+ g(ega[iEgaInt],q2_2)*getcsgA(ega[iEgaInt],q2_2,beam)
+ 4.*g(ega[iEgaInt],q2_12)*getcsgA(ega[iEgaInt],q2_12,beam) );
// Effective flux
dndE[iEgaInt] +=(q2_2-q2_1)*( getcsgA_Q2_dep(q2_1)*photonFlux(ega[iEgaInt],q2_1)
+getcsgA_Q2_dep(q2_2)*photonFlux(ega[iEgaInt],q2_2)
+4.*getcsgA_Q2_dep(q2_12)*photonFlux(ega[iEgaInt],q2_12) );
}
- q2end = -99 ;
full_int[iEgaInt] = full_int[iEgaInt]/6.;
dndE[iEgaInt] = dndE[iEgaInt]/6.;
}
// Finishing cross-section integral
dR = full_int[0];
dR += full_int[1];
dR += 4.*full_int[2];
dR = dR*(ega[1]-ega[0])/6.;
int_r = int_r + dR*breitWigner(W,bwnorm)*dW;
// Finishing effective flux integral
// Finishing integral over the effective photon flux
dR2 = dndE[0];
dR2 += dndE[1];
dR2 += 4.*dndE[2];
dR2 = dR2*(ega[1]-ega[0])/6.;
int_r2 = int_r2 + dR2*breitWigner(W,bwnorm)*dW;
}
}
cout< VM+X ", int_r/int_r2);
setPhotonNucleusSigma(0.01*int_r);
//
#ifdef _makeGammaPQ2_
makeGammaPQ2dependence(bwnormsave);
#endif
}
#ifdef _makeGammaPQ2_
//______________________________________________________________________________
void
e_wideResonanceCrossSection::makeGammaPQ2dependence( double bwnormsave)
{
// This subroutine calculates the Q2 dependence of
// gamma+X -> VM + X cross section for a narrow resonance
int const nQ2bins = 19;
double const q2Edge[nQ2bins+1] = { 0.,1.,2.,3., 4.,5.,
6.,7.,8.,9.,10.,
11.,12.,13.,14.,15.,
20.,30.,40.,50.};
double W = 0,dW,dY;
double y1,y2,y12,ega1,ega2,ega12;
double int_r,dR,dR2;
double csgA1, csgA2, csgA12;
double Eth;
int I,J,NW,NY,beam;
double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
NW = 100;
dW = (_wideWmax-_wideWmin)/double(NW);
if (getBNORM() == 0.){
cout<<" Using Breit-Wigner Resonance Profile."< VM_p "< VM_p (nanob)"<= 1 ){
// pA, first beam is the nucleus and is in this case the target
ega1 = 0.5*W*exp(-y1);
ega2 = 0.5*W*exp(-y2);
ega12 = 0.5*W*exp(-y12);
beam = 1;
} else if( A_1 ==0 && A_2 >= 1){
// pA, second beam is the nucleus and is in this case the target
ega1 = 0.5*W*exp(y1);
ega2 = 0.5*W*exp(y2);
ega12 = 0.5*W*exp(y12);
// photon energy in Target frame
beam = 2;
} else {
ega1 = 0.5*W*exp(y1);
ega2 = 0.5*W*exp(y2);
ega12 = 0.5*W*exp(y12);
// photon energy in Target frame
beam = 2;
}
//
if(ega1 < Eth || ega2 < Eth)
continue;
if(ega2 > maxPhotonEnergy() || ega1 > maxPhotonEnergy() )
continue;
target_ega1 = ega1*exp_target_cm;
target_ega12 = ega12*exp_target_cm;
target_ega2 = ega2*exp_target_cm;
//cout<<"Nortmalizations "<> Middle Point =====>>>
csgA12=getcsgA(ega12,W,beam);
double full_range_12 = integrated_x_section(target_ega12);
// >> Second Point =====>>>
csgA2=getcsgA(ega2,W,beam);
double full_range_2 = integrated_x_section(target_ega2);
//
//>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
dR = q2_cor*csgA1;
dR = dR + 4.*q2_cor*csgA12;
dR = dR + q2_cor*csgA2;
dR = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
//
dR2 = full_range_1*csgA1;
dR2 = dR2 + 4.*full_range_12*csgA12;
dR2 = dR2 + full_range_2*csgA2;
dR2 = dR2*(dY/6.)*breitWigner(W,bwnorm)*dW;
//
int_r = int_r+dR;
int_r2 = int_r2 +dR2;
}
//
}
if( iQ2 ==0 )
cout<<"Full range "< 1.){
cout<< name.c_str() <<0.01*x_section<<" barn."< 1.){
cout<< name.c_str() <<10.*x_section<<" mb."< 1.){
cout<< name.c_str() <<10000.*x_section<<" microb."< 1.){
cout<< name.c_str() <<10000000.*x_section<<" nanob."< 1.){
cout<< name.c_str() <<1.E10*x_section<<" picob."<.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 276 $: revision of last commit
// $Author:: jnystrand $: author of last commit
// $Date:: 2016-09-13 19:54:42 +0100 #$: date of last commit
//
// Description:
//
//
///////////////////////////////////////////////////////////////////////////
#include
#include
#include
#include "reportingUtils.h"
#include "starlightconstants.h"
#include "bessel.h"
#include "photonNucleusCrossSection.h"
using namespace std;
using namespace starlightConstants;
//______________________________________________________________________________
photonNucleusCrossSection::photonNucleusCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
: _nWbins (inputParametersInstance.nmbWBins() ),
_nYbins (inputParametersInstance.nmbRapidityBins() ),
_wMin (inputParametersInstance.minW() ),
_wMax (inputParametersInstance.maxW() ),
_yMax (inputParametersInstance.maxRapidity() ),
_beamLorentzGamma (inputParametersInstance.beamLorentzGamma() ),
_bbs (bbsystem ),
_protonEnergy (inputParametersInstance.protonEnergy() ),
_electronEnergy (inputParametersInstance.electronEnergy() ),
_particleType (inputParametersInstance.prodParticleType() ),
_beamBreakupMode (inputParametersInstance.beamBreakupMode() ),
_productionMode (inputParametersInstance.productionMode() ),
_sigmaNucleus (_bbs.beam2().A() ),
_fixedQ2range (inputParametersInstance.fixedQ2Range() ),
_minQ2 (inputParametersInstance.minGammaQ2() ),
_maxQ2 (inputParametersInstance.maxGammaQ2() ),
_maxPhotonEnergy (inputParametersInstance.cmsMaxPhotonEnergy()),
_cmsMinPhotonEnergy(inputParametersInstance.cmsMinPhotonEnergy()),
_targetRadii (inputParametersInstance.targetRadius() )
{
// new options - impulse aproximation (per Joakim) and Quantum Glauber (per SK) SKQG
_impulseSelected = inputParametersInstance.impulseVM();
_quantumGlauber = inputParametersInstance.quantumGlauber();
switch(_particleType) {
case RHO:
- _slopeParameter = 11.0; // [(GeV/c)^{-2}]
- _vmPhotonCoupling = 2.02;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [ GeV^{-2}]
- _ANORM = -2.75;
- _BNORM = 0.0;
- _defaultC = 1.0;
- _channelMass = starlightConstants::rho0Mass;
- _width = starlightConstants::rho0Width;
- break;
+ _slopeParameter = 11.0; // [(GeV/c)^{-2}]
+ _vmPhotonCoupling = 2.02;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [ GeV^{-2}]
+ _ANORM = -2.75;
+ _BNORM = 0.0;
+ _defaultC = 1.0;
+ _channelMass = starlightConstants::rho0Mass;
+ _width = starlightConstants::rho0Width;
+ break;
case RHOZEUS:
- _slopeParameter =11.0;
- _vmPhotonCoupling=2.02;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=1.84;
- _defaultC=1.0;
- _channelMass = starlightConstants::rho0Mass;
- _width = starlightConstants::rho0Width;
- break;
+ _slopeParameter =11.0;
+ _vmPhotonCoupling=2.02;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=1.84;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::rho0Mass;
+ _width = starlightConstants::rho0Width;
+ break;
case FOURPRONG:
- _slopeParameter = 11.0;
- _vmPhotonCoupling = 2.02;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM = -2.75;
- _BNORM = 0;
- _defaultC = 11.0;
- _channelMass = starlightConstants::rho0PrimeMass;
- _width = starlightConstants::rho0PrimeWidth;
- break;
+ _slopeParameter = 11.0;
+ _vmPhotonCoupling = 2.02;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM = -2.75;
+ _BNORM = 0;
+ _defaultC = 11.0;
+ _channelMass = starlightConstants::rho0PrimeMass;
+ _width = starlightConstants::rho0PrimeWidth;
+ break;
case OMEGA:
- _slopeParameter=10.0;
- _vmPhotonCoupling=23.13;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::OmegaMass;
- _width = starlightConstants::OmegaWidth;
- break;
+ _slopeParameter=10.0;
+ _vmPhotonCoupling=23.13;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::OmegaMass;
+ _width = starlightConstants::OmegaWidth;
+ break;
case PHI:
- _slopeParameter=7.0;
- _vmPhotonCoupling=13.71;
- _vmQ2Power_c1 = 2.15;
- _vmQ2Power_c2 = 0.0074; // [GeV^{-2}]
- //_vmQ2Power_c1 = 2.45;
- //_vmQ2Power_c2 = 0.0011; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::PhiMass;
- _width = starlightConstants::PhiWidth;
- break;
+ _slopeParameter=7.0;
+ _vmPhotonCoupling=13.71;
+ _vmQ2Power_c1 = 2.15;
+ _vmQ2Power_c2 = 0.0074; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::PhiMass;
+ _width = starlightConstants::PhiWidth;
+ break;
case JPSI:
case JPSI_ee:
- //_vmQ2Power_c1 = 2.44;
- // _vmQ2Power_c2 = 0.0074; // [GeV^{-2}]
+ _slopeParameter=4.0;
+ _vmPhotonCoupling=10.45;
_vmQ2Power_c1 = 2.45;
_vmQ2Power_c2 = 0.00084; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::JpsiMass;
+ _width = starlightConstants::JpsiWidth;
+ break;
case JPSI_mumu:
- _slopeParameter=4.0;
- _vmPhotonCoupling=10.45;
- // Best so far
- //_vmQ2Power_c1 = 2.60;
- //_vmQ2Power_c2 = 0.0057; // [GeV^{-2}]
- //Clsoe second
- //_vmQ2Power_c1 = 2.30;
- //_vmQ2Power_c2 = 0.0098; // [GeV^{-2}]
- //New winner
- _vmQ2Power_c1 = 2.36;
- _vmQ2Power_c2 = 0.0029; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::JpsiMass;
- _width = starlightConstants::JpsiWidth;
- break;
+ _slopeParameter=4.0;
+ _vmPhotonCoupling=10.45;
+ _vmQ2Power_c1 = 2.36;
+ _vmQ2Power_c2 = 0.0029; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::JpsiMass;
+ _width = starlightConstants::JpsiWidth;
+ break;
case JPSI2S:
case JPSI2S_ee:
case JPSI2S_mumu:
- _slopeParameter=4.3;
- _vmPhotonCoupling=26.39;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::Psi2SMass;
- _width = starlightConstants::Psi2SWidth;
- break;
+ _slopeParameter=4.3;
+ _vmPhotonCoupling=26.39;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::Psi2SMass;
+ _width = starlightConstants::Psi2SWidth;
+ break;
case UPSILON:
case UPSILON_ee:
case UPSILON_mumu:
- _slopeParameter=4.0;
- _vmPhotonCoupling=125.37;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::Upsilon1SMass;
- _width = starlightConstants::Upsilon1SWidth;
- break;
+ _slopeParameter=4.0;
+ _vmPhotonCoupling=125.37;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::Upsilon1SMass;
+ _width = starlightConstants::Upsilon1SWidth;
+ break;
case UPSILON2S:
case UPSILON2S_ee:
case UPSILON2S_mumu:
- _slopeParameter=4.0;
- _vmPhotonCoupling=290.84;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::Upsilon2SMass;
- _width = starlightConstants::Upsilon2SWidth;
- break;
+ _slopeParameter=4.0;
+ _vmPhotonCoupling=290.84;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::Upsilon2SMass;
+ _width = starlightConstants::Upsilon2SWidth;
+ break;
case UPSILON3S:
case UPSILON3S_ee:
case UPSILON3S_mumu:
- _slopeParameter=4.0;
- _vmPhotonCoupling=415.10;
- _vmQ2Power_c1 = 2.09;
- _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
- _ANORM=-2.75;
- _BNORM=0.0;
- _defaultC=1.0;
- _channelMass = starlightConstants::Upsilon3SMass;
- _width = starlightConstants::Upsilon3SWidth;
- break;
+ _slopeParameter=4.0;
+ _vmPhotonCoupling=415.10;
+ _vmQ2Power_c1 = 2.09;
+ _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
+ _ANORM=-2.75;
+ _BNORM=0.0;
+ _defaultC=1.0;
+ _channelMass = starlightConstants::Upsilon3SMass;
+ _width = starlightConstants::Upsilon3SWidth;
+ break;
default:
cout <<"No sigma constants parameterized for pid: "<<_particleType
<<" GammaAcrosssection"< 1 || 2.*targetEgamma/(Q2+W*W) < _targetRadii)
return 0;
double temp_pz = sqrt(pz_squared);
// Now boost to CM frame
double Egamma = targetEgamma*cosh(beam_y) - temp_pz*sinh(beam_y);
if( Egamma < _cmsMinPhotonEnergy || Egamma > _maxPhotonEnergy){
return 0;
}
//cout<<" ::: Lomnitz test in photonNucleus ::: pz^2 = "<20R_A change methods - just take the photon flux
// at the center of the nucleus.
if (biter > (10.*RNuc)){
// if there is no nuclear breakup or only hadronic breakup, which only
// occurs at smaller b, we can analytically integrate the flux from b~20R_A
// to infinity, following Jackson (2nd edition), Eq. 15.54
Xvar=energy*biter/(hbarc*_beamLorentzGamma);
// Here, there is nuclear breakup. So, we can't use the integrated flux
// However, we can do a single flux calculation, at the center of the nucleus
// Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
// this is the flux per unit area
fluxelement = (rZ*rZ*alpha*energy)*
(bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
((pi*_beamLorentzGamma*hbarc)*
(pi*_beamLorentzGamma*hbarc));
} else {
// integrate over nuclear surface. n.b. this assumes total shadowing -
// treat photons hitting the nucleus the same no matter where they strike
fluxelement=0.;
deltar=RNuc/double(nrstep);
riter=-deltar/2.;
for (int jr =1; jr<=nrstep;jr++){
riter=riter+deltar;
// use symmetry; only integrate from 0 to pi (half circle)
deltaphi=pi/double(nphistep);
phiiter=0.;
for( int jphi=1;jphi<= nphistep;jphi++){
phiiter=(double(jphi)-0.5)*deltaphi;
// dist is the distance from the center of the emitting nucleus
// to the point in question
dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
Xvar=energy*dist/(hbarc*_beamLorentzGamma);
flux_r = (rZ*rZ*alpha*energy)*
(bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
((pi*_beamLorentzGamma*hbarc)*
(pi*_beamLorentzGamma*hbarc));
// The surface element is 2.* delta phi* r * delta r
// The '2' is because the phi integral only goes from 0 to pi
fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
// end phi and r integrations
}//for(jphi)
}//for(jr)
// average fluxelement over the nuclear surface
fluxelement=fluxelement/(pi*RNuc*RNuc);
}//else
// multiply by volume element to get total flux in the volume element
fluxelement=fluxelement*2.*pi*biter*(biter-bold);
// modulate by the probability of nuclear breakup as f(biter)
// cout<<" jb: "< 1){
fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
}
// cout<<" jb: "< lnEmax){
flux_r=0.0;
cout<<" ERROR: Egamma outside defined range. Egamma= "<> Egamma between Ilt and Ilt+1
Ilt = int((lEgamma-lnEmin)/dlnE);
// >> ln(Egamma) for first point
lnElt = lnEmin + Ilt*dlnE;
// >> Interpolate
flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
flux_r = flux_r/Egamma;
}
return flux_r;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::integrated_Q2_dep(double const Egamma, double const _min, double const _max)
{
//Integration over full limits gives more accurate result
double Q2_min = std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
double Q2_max = 4.*_electronEnergy*(_electronEnergy-Egamma);
//double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
if( _min != 0 || _max !=0){
if( _min > Q2_min )
Q2_min = _min;
if( _max < Q2_max )
Q2_max = _max;
}
// Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
// nstep = 10,000 100,000 & 10,000,0000
int nstep = 1000;
double ln_min = std::log(Q2_min);
double ratio = std::log(Q2_max/Q2_min)/nstep;
double g_int = 0;
double g_int2 = 0 ;
double g_int3 = 0;
//cout<<"*** Lomnitz **** Energy "< Q2_min)
Q2_min = _min;
if( _max < Q2_max)
Q2_max = _max;
}
// Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
// nstep = 10,000 100,000 & 10,000,0000
int nstep = 1000;
double ln_min = std::log(Q2_min);
double ratio = std::log(Q2_max/Q2_min)/nstep;
double g_int = 0;
for ( int ii = 0 ; ii< nstep; ++ii){
double x1 = std::exp(ln_min+(double)ii*ratio);
double x3 = std::exp(ln_min+(double)(ii+1)*ratio);
double x2 = (x3+x1)/2.;
//Tests from HERA https://arxiv.org/pdf/hep-ex/9601009.pdf
// g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
}
//return g_int2*g_int3/36.;
return g_int/6.;
}
//______________________________________________________________________________
pair< double, double >*photonNucleusCrossSection::Q2arraylimits(double const Egamma)
{
//double Q2max = 2.*Egamma/_targetRadii - _wMax*_wMax;
double Q2max = 4.*_electronEnergy*(_electronEnergy-Egamma);
double Q2min= std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
if( _fixedQ2range == true){
if( Q2min < _minQ2 )
Q2min = _minQ2;
if( Q2max > _maxQ2 )
Q2max = _maxQ2;
//cout<<" Imposed limits "<* to_ret = new std::pair(Q2min,Q2max);
return to_ret;
}
int Nstep = 1000;
//
double ratio = std::log(Q2max/Q2min)/(double)Nstep;
double ln_min = std::log(Q2min);
// -- - -
const double limit = 1E-9;
std::vectorQ2_array;
int iNstep = 0;
double g_step = 1.;
while( g_step>limit ){
double Q2 = std::exp(ln_min+iNstep*ratio);
if(Q2>Q2max)
break;
g_step = g(Egamma,Q2);
Q2_array.push_back(g_step);
iNstep++;
}
if( std::exp(ln_min+iNstep*ratio) < Q2max)
Q2max = std::exp(ln_min+iNstep*ratio);
//cout<* to_ret = new std::pair(Q2min,Q2max);
return to_ret;
}
//______________________________________________________________________________
double
photonNucleusCrossSection::g(double const Egamma,
double const Q2)
{
//return photonFlux(Egamma,Q2)*getcsgA_Q2_dep(Q2);
return photonFlux(Egamma,Q2); //This could be done more elegantly in the future, but this one change should account for everything
}
//______________________________________________________________________________
double
photonNucleusCrossSection::photonFlux(double const Egamma,
double const Q2)
{
//Need to check units later
//double const hbar = starlightConstants::hbarc / 2.99*pow(10,14); // 6.582 * pow (10,-16) [eVs]
//double omega = Egamma/ hbar;
//Do we even need a lookup table for this case? This should return N(E,Q2) from dn = N(E,Q2) dE dQ2
double const ratio = Egamma/_electronEnergy;
double const minQ2 = std::pow( starlightConstants::mel*Egamma,2.0) / (_electronEnergy*(_electronEnergy - Egamma));
double to_ret = alpha/(pi) *( 1- ratio + ratio*ratio/2. - (1-ratio)*( fabs(minQ2/Q2)) );
//Comparisons:
// double temp = pow(2.*_electronEnergy-Egamma,2.)/(Egamma*Egamma + Q2) + 1. - 4.*starlightConstants::mel*starlightConstants::mel/Q2;
//temp = alpha*temp*Egamma/(4.*Q2*pi*_electronEnergy*_electronEnergy);
// cout<<" *** Lomnitz *** Testing photon flux approximation for electron energy "<<_electronEnergy<<" gamma "< VectorMeson
// cross section. Wgp is the gamma-proton CM energy.
// Unit for cross section: fm**2
double sigmagp_r=0.;
switch(_particleType)
{
case RHO:
case RHOZEUS:
case FOURPRONG:
sigmagp_r=1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
break;
case OMEGA:
sigmagp_r=1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
break;
case PHI:
sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
break;
case JPSI:
case JPSI_ee:
case JPSI_mumu:
sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
// sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
break;
case JPSI2S:
case JPSI2S_ee:
case JPSI2S_mumu:
sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
sigmagp_r*=0.166;
// sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
break;
case UPSILON:
case UPSILON_ee:
case UPSILON_mumu:
// >> This is W**1.7 dependence from QCD calculations
// sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-10*6.4*exp(0.74*log(Wgp));
break;
case UPSILON2S:
case UPSILON2S_ee:
case UPSILON2S_mumu:
// sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-10*2.9*exp(0.74*log(Wgp));
break;
case UPSILON3S:
case UPSILON3S_ee:
case UPSILON3S_mumu:
// sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
sigmagp_r*=sigmagp_r;
sigmagp_r*=1.E-10*2.1*exp(0.74*log(Wgp));
break;
default: cout<< "!!! ERROR: Unidentified Vector Meson: "<< _particleType <