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 <