]>
Commit | Line | Data |
---|---|---|
da32329d AM |
1 | /////////////////////////////////////////////////////////////////////////// |
2 | // | |
3 | // Copyright 2010 | |
4 | // | |
5 | // This file is part of starlight. | |
6 | // | |
7 | // starlight is free software: you can redistribute it and/or modify | |
8 | // it under the terms of the GNU General Public License as published by | |
9 | // the Free Software Foundation, either version 3 of the License, or | |
10 | // (at your option) any later version. | |
11 | // | |
12 | // starlight is distributed in the hope that it will be useful, | |
13 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
15 | // GNU General Public License for more details. | |
16 | // | |
17 | // You should have received a copy of the GNU General Public License | |
18 | // along with starlight. If not, see <http://www.gnu.org/licenses/>. | |
19 | // | |
20 | /////////////////////////////////////////////////////////////////////////// | |
21 | // | |
22 | // File and Version Information: | |
23 | // $Rev:: 45 $: revision of last commit | |
24 | // $Author:: bgrube $: author of last commit | |
25 | // $Date:: 2011-02-27 20:52:35 +0100 #$: date of last commit | |
26 | // | |
27 | // Description: | |
28 | // | |
29 | // | |
30 | // | |
31 | /////////////////////////////////////////////////////////////////////////// | |
32 | ||
33 | ||
34 | #include <iostream> | |
35 | #include <fstream> | |
36 | #include <cmath> | |
37 | ||
38 | #include "inputParameters.h" | |
39 | #include "beambeamsystem.h" | |
40 | #include "beam.h" | |
41 | #include "starlightconstants.h" | |
42 | #include "nucleus.h" | |
43 | #include "bessel.h" | |
44 | #include "incoherentPhotonNucleusLuminosity.h" | |
45 | ||
46 | ||
47 | using namespace std; | |
45d54d9a | 48 | using namespace starlightConstants; |
da32329d AM |
49 | |
50 | ||
51 | //______________________________________________________________________________ | |
52 | incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusLuminosity(beamBeamSystem& bbsystem) | |
53 | : photonNucleusCrossSection(bbsystem) | |
54 | { | |
55 | cout <<"Creating Luminosity Tables for incoherent vector meson production."<<endl; | |
56 | incoherentPhotonNucleusDifferentialLuminosity(); | |
57 | cout <<"Luminosity Tables created."<<endl; | |
58 | } | |
59 | ||
60 | ||
61 | //______________________________________________________________________________ | |
62 | incoherentPhotonNucleusLuminosity::~incoherentPhotonNucleusLuminosity() | |
63 | { } | |
64 | ||
65 | ||
66 | //______________________________________________________________________________ | |
67 | void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLuminosity() | |
68 | { | |
69 | // double Av,Wgp,cs,cvma; | |
70 | double W,dW,dY; | |
71 | double Egamma,Y; | |
72 | // double t,tmin,tmax; | |
73 | double testint,dndWdY; | |
74 | // double ax,bx; | |
75 | double C; | |
76 | ||
77 | ofstream wylumfile; | |
78 | wylumfile.precision(15); | |
79 | ||
80 | double bwnorm,Eth; | |
81 | ||
82 | dW = (_wMax - _wMin)/_nWbins; | |
83 | dY = (_yMax-(-1.0)*_yMax)/_nYbins; | |
84 | ||
85 | // Write the values of W used in the calculation to slight.txt. | |
86 | wylumfile.open("slight.txt"); | |
be31281a | 87 | // wylumfile << inputParametersInstance.parameterValueKey() << endl; |
da32329d AM |
88 | wylumfile << getbbs().beam1().Z() <<endl; |
89 | wylumfile << getbbs().beam1().A() <<endl; | |
90 | wylumfile << getbbs().beam2().Z() <<endl; | |
91 | wylumfile << getbbs().beam2().A() <<endl; | |
92 | wylumfile << inputParametersInstance.beamLorentzGamma() <<endl; | |
93 | wylumfile << inputParametersInstance.maxW() <<endl; | |
94 | wylumfile << inputParametersInstance.minW() <<endl; | |
95 | wylumfile << inputParametersInstance.nmbWBins() <<endl; | |
96 | wylumfile << inputParametersInstance.maxRapidity() <<endl; | |
97 | wylumfile << inputParametersInstance.nmbRapidityBins() <<endl; | |
98 | wylumfile << inputParametersInstance.productionMode() <<endl; | |
99 | wylumfile << inputParametersInstance.beamBreakupMode() <<endl; | |
100 | wylumfile << inputParametersInstance.interferenceEnabled() <<endl; | |
101 | wylumfile << inputParametersInstance.interferenceStrength() <<endl; | |
102 | wylumfile << inputParametersInstance.coherentProduction() <<endl; | |
103 | wylumfile << inputParametersInstance.incoherentFactor() <<endl; | |
45d54d9a | 104 | wylumfile << starlightConstants::deuteronSlopePar <<endl; |
da32329d AM |
105 | wylumfile << inputParametersInstance.maxPtInterference() <<endl; |
106 | wylumfile << inputParametersInstance.nmbPtBinsInterference() <<endl; | |
107 | ||
108 | // Normalize the Breit-Wigner Distribution and write values of W to slight.txt | |
109 | testint=0.0; | |
110 | //Grabbing default value for C in the breit-wigner calculation | |
111 | C=getDefaultC(); | |
112 | for(unsigned int i = 0; i <= _nWbins - 1; ++i) { | |
113 | W = _wMin + double(i)*dW + 0.5*dW; | |
114 | testint = testint + breitWigner(W,C)*dW; | |
115 | wylumfile << W << endl; | |
116 | } | |
117 | bwnorm = 1./testint; | |
118 | ||
119 | // Write the values of Y used in the calculation to slight.txt. | |
120 | for(unsigned int i = 0; i <= _nYbins - 1; ++i) { | |
121 | Y = -1.0*_yMax + double(i)*dY + 0.5*dY; | |
122 | wylumfile << Y << endl; | |
123 | } | |
124 | ||
125 | // Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin | |
126 | // +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/ | |
127 | // (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass))); | |
128 | ||
129 | for(unsigned int i = 0; i <= _nWbins - 1; ++i) { | |
130 | ||
131 | W = _wMin + double(i)*dW + 0.5*dW; | |
132 | ||
133 | double Ep = inputParametersInstance.protonEnergy(); | |
134 | ||
135 | Eth=0.5*(((W+starlightConstants::protonMass)*(W+starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/ | |
136 | (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass))); | |
137 | ||
138 | for(unsigned int j = 0; j <= _nYbins - 1; ++j) { | |
139 | ||
140 | Y = -1.0*_yMax + double(j)*dY + 0.5*dY; | |
141 | ||
142 | int A_1 = getbbs().beam1().A(); | |
143 | int A_2 = getbbs().beam2().A(); | |
144 | if( A_2 == 1 && A_1 != 1 ){ | |
145 | // pA, first beam is the nucleus | |
146 | Egamma = 0.5*W*exp(Y); | |
147 | } else if( A_1 ==1 && A_2 != 1){ | |
148 | // pA, second beam is the nucleus | |
149 | Egamma = 0.5*W*exp(-Y); | |
150 | } else { | |
151 | Egamma = 0.5*W*exp(Y); | |
152 | } | |
153 | ||
154 | dndWdY = 0.; | |
155 | ||
156 | if(Egamma > Eth){ | |
157 | if(Egamma > maxPhotonEnergy())Egamma = maxPhotonEnergy(); | |
158 | double Wgp = sqrt(2.*Egamma*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass* | |
159 | starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass); | |
160 | ||
161 | double localsig = sigmagp(Wgp); | |
162 | // int localz = 0; | |
163 | // double localbmin = 0; | |
164 | if( A_1 == 1 && A_2 != 1 ){ | |
165 | // localbmin = getbbs().beam2().nuclearRadius() + 0.7; | |
166 | // localz = getbbs().beam2().Z(); | |
167 | // dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm); | |
168 | dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm); | |
169 | }else if (A_2 ==1 && A_1 !=1){ | |
170 | // localbmin = getbbs().beam1().nuclearRadius() + 0.7; | |
171 | // localz = getbbs().beam1().Z(); | |
172 | // dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm); | |
173 | dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm); | |
174 | }else{ | |
175 | double csVN = sigma_N(Wgp); | |
176 | double csVA = sigma_A(csVN); | |
177 | double csgA= (csVA/csVN)*sigmagp(Wgp); | |
178 | dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm); | |
179 | } | |
180 | } | |
181 | ||
182 | wylumfile << dndWdY << endl; | |
be31281a | 183 | |
da32329d AM |
184 | } |
185 | } | |
186 | ||
187 | wylumfile << bwnorm << endl; | |
188 | wylumfile << inputParametersInstance.parameterValueKey() << endl; | |
189 | wylumfile.close(); | |
190 | ||
be31281a | 191 | // cout << "bwnorm: "<< bwnorm <<endl; |
192 | ||
da32329d AM |
193 | } |
194 | ||
195 | ||
196 | //______________________________________________________________________________ | |
197 | double incoherentPhotonNucleusLuminosity::nofe(double Egamma, double bimp) | |
198 | { | |
199 | //Function for the calculation of the "photon density". | |
200 | //nofe=numberofgammas/(energy*area) | |
201 | //Assume beta=1.0 and gamma>>1, i.e. neglect the (1/gamma^2)*K0(x) term | |
202 | ||
203 | double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.; | |
204 | ||
205 | X = (bimp*Egamma)/(_beamLorentzGamma*starlightConstants::hbarc); | |
206 | ||
207 | if( X <= 0.0) | |
208 | cout<<"In nofe, X= "<<X<<endl; | |
209 | ||
210 | factor1 = (double(getbbs().beam1().Z()*getbbs().beam1().Z()) | |
211 | *starlightConstants::alpha)/(starlightConstants::pi*starlightConstants::pi); | |
212 | ||
213 | factor2 = 1./(Egamma*bimp*bimp); | |
214 | factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X)); | |
215 | nofex = factor1*factor2*factor3; | |
216 | return nofex; | |
217 | } |