--- /dev/null
+///////////////////////////////////////////////////////////////////////////
+//
+// 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:: 45 $: revision of last commit
+// $Author:: bgrube $: author of last commit
+// $Date:: 2011-02-27 20:52:35 +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 "incoherentPhotonNucleusLuminosity.h"
+
+
+using namespace std;
+
+
+//______________________________________________________________________________
+incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusLuminosity(beamBeamSystem& bbsystem)
+ : photonNucleusCrossSection(bbsystem)
+{
+ cout <<"Creating Luminosity Tables for incoherent vector meson production."<<endl;
+ incoherentPhotonNucleusDifferentialLuminosity();
+ cout <<"Luminosity Tables created."<<endl;
+}
+
+
+//______________________________________________________________________________
+incoherentPhotonNucleusLuminosity::~incoherentPhotonNucleusLuminosity()
+{ }
+
+
+//______________________________________________________________________________
+void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLuminosity()
+{
+ // double Av,Wgp,cs,cvma;
+ double W,dW,dY;
+ double Egamma,Y;
+ // double t,tmin,tmax;
+ double testint,dndWdY;
+ // double ax,bx;
+ double C;
+
+ ofstream wylumfile;
+ wylumfile.precision(15);
+
+ double bwnorm,Eth;
+
+ dW = (_wMax - _wMin)/_nWbins;
+ dY = (_yMax-(-1.0)*_yMax)/_nYbins;
+
+ // Write the values of W used in the calculation to slight.txt.
+ wylumfile.open("slight.txt");
+ wylumfile << inputParametersInstance.parameterValueKey() << endl;
+ wylumfile << getbbs().beam1().Z() <<endl;
+ wylumfile << getbbs().beam1().A() <<endl;
+ wylumfile << getbbs().beam2().Z() <<endl;
+ wylumfile << getbbs().beam2().A() <<endl;
+ wylumfile << inputParametersInstance.beamLorentzGamma() <<endl;
+ wylumfile << inputParametersInstance.maxW() <<endl;
+ wylumfile << inputParametersInstance.minW() <<endl;
+ wylumfile << inputParametersInstance.nmbWBins() <<endl;
+ wylumfile << inputParametersInstance.maxRapidity() <<endl;
+ wylumfile << inputParametersInstance.nmbRapidityBins() <<endl;
+ wylumfile << inputParametersInstance.productionMode() <<endl;
+ wylumfile << inputParametersInstance.beamBreakupMode() <<endl;
+ wylumfile << inputParametersInstance.interferenceEnabled() <<endl;
+ wylumfile << inputParametersInstance.interferenceStrength() <<endl;
+ wylumfile << inputParametersInstance.coherentProduction() <<endl;
+ wylumfile << inputParametersInstance.incoherentFactor() <<endl;
+ wylumfile << inputParametersInstance.deuteronSlopePar() <<endl;
+ wylumfile << inputParametersInstance.maxPtInterference() <<endl;
+ wylumfile << inputParametersInstance.nmbPtBinsInterference() <<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;
+ }
+
+ // Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin
+ // +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
+ // (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
+
+ for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
+
+ W = _wMin + double(i)*dW + 0.5*dW;
+
+ double Ep = inputParametersInstance.protonEnergy();
+
+ Eth=0.5*(((W+starlightConstants::protonMass)*(W+starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
+ (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
+
+ for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
+
+ Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
+
+ int A_1 = getbbs().beam1().A();
+ int A_2 = getbbs().beam2().A();
+ if( A_2 == 1 && A_1 != 1 ){
+ // pA, first beam is the nucleus
+ Egamma = 0.5*W*exp(Y);
+ } else if( A_1 ==1 && A_2 != 1){
+ // pA, second beam is the nucleus
+ Egamma = 0.5*W*exp(-Y);
+ } else {
+ Egamma = 0.5*W*exp(Y);
+ }
+
+ dndWdY = 0.;
+
+ if(Egamma > Eth){
+ if(Egamma > maxPhotonEnergy())Egamma = maxPhotonEnergy();
+ double Wgp = sqrt(2.*Egamma*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
+ starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
+
+ double localsig = sigmagp(Wgp);
+ // int localz = 0;
+ // double localbmin = 0;
+ if( A_1 == 1 && A_2 != 1 ){
+ // localbmin = getbbs().beam2().nuclearRadius() + 0.7;
+ // localz = getbbs().beam2().Z();
+ // dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm);
+ dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm);
+ }else if (A_2 ==1 && A_1 !=1){
+ // localbmin = getbbs().beam1().nuclearRadius() + 0.7;
+ // localz = getbbs().beam1().Z();
+ // dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm);
+ dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm);
+ }else{
+ double csVN = sigma_N(Wgp);
+ double csVA = sigma_A(csVN);
+ double csgA= (csVA/csVN)*sigmagp(Wgp);
+ dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm);
+ }
+ }
+
+ wylumfile << dndWdY << endl;
+ }
+ }
+
+ wylumfile << bwnorm << endl;
+ wylumfile << inputParametersInstance.parameterValueKey() << endl;
+ wylumfile.close();
+
+// wylumfile.open("slight.txt",ios::app);
+ cout << "bwnorm: "<< bwnorm <<endl;
+// wylumfile << bwnorm << endl;
+// wylumfile << inputParametersInstance.parameterValueKey() << endl;
+// wylumfile.close();
+}
+
+
+//______________________________________________________________________________
+double incoherentPhotonNucleusLuminosity::nofe(double Egamma, double bimp)
+{
+ //Function for the calculation of the "photon density".
+ //nofe=numberofgammas/(energy*area)
+ //Assume beta=1.0 and gamma>>1, i.e. neglect the (1/gamma^2)*K0(x) term
+
+ double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.;
+
+ X = (bimp*Egamma)/(_beamLorentzGamma*starlightConstants::hbarc);
+
+ if( X <= 0.0)
+ cout<<"In nofe, X= "<<X<<endl;
+
+ factor1 = (double(getbbs().beam1().Z()*getbbs().beam1().Z())
+ *starlightConstants::alpha)/(starlightConstants::pi*starlightConstants::pi);
+
+ factor2 = 1./(Egamma*bimp*bimp);
+ factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X));
+ nofex = factor1*factor2*factor3;
+ return nofex;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff