EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
incoherentVMCrossSection.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file incoherentVMCrossSection.cpp
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 //
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 //
32 
33 
34 #include <iostream>
35 #include <fstream>
36 #include <cmath>
37 
38 #include "starlightconstants.h"
39 #include "incoherentVMCrossSection.h"
40 
41 
42 using namespace std;
43 using namespace starlightConstants;
44 
45 
46 //______________________________________________________________________________
47 incoherentVMCrossSection::incoherentVMCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
48  :photonNucleusCrossSection(inputParametersInstance, bbsystem)
49 {
50  _narrowYmax = inputParametersInstance.maxRapidity();
51  _narrowYmin = -1.0*_narrowYmax;
52  _narrowNumY = inputParametersInstance.nmbRapidityBins();
53  _Ep = inputParametersInstance.protonEnergy();
54 }
55 
56 
57 //______________________________________________________________________________
58 incoherentVMCrossSection::~incoherentVMCrossSection()
59 { }
60 
61 
62 //______________________________________________________________________________
63 void
64 incoherentVMCrossSection::crossSectionCalculation(const double) // _bwnormsave (unused)
65 {
66  // This subroutine calculates the vector meson cross section assuming
67  // a narrow resonance. For reference, see STAR Note 386.
68 
69  double W,dY;
70  double y1,y2,y12,ega1,ega2,ega12;
71  double csgA1,csgA2,csgA12,int_r,dR;
72  double Wgp,csVN,csVA;
73  double Eth;
74  int J,NY,beam;
75 
76  NY = _narrowNumY;
77  dY = (_narrowYmax-_narrowYmin)/double(NY);
78 
79  cout<<" Using Narrow Resonance ..."<<endl;
80 
81  W = getChannelMass();
82  Eth=0.5*(((W+protonMass)*(W+protonMass)-
83  protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
84 
85  // cout<<" gamma+nucleon Threshold: "<<Eth<<endl;
86  printf(" gamma+nucleon threshold: %e GeV \n", Eth);
87 
88  int A_1 = getbbs().electronBeam().A();
89  int A_2 = getbbs().targetBeam().A();
90 
91  int_r=0.;
92 
93  // Do this first for the case when the first beam is the photon emitter
94  // Treat pA separately with defined beams
95  // The variable beam (=1,2) defines which nucleus is the target
96  for(J=0;J<=(NY-1);J++){
97 
98  // This is the fdefault
99  y1 = _narrowYmin + double(J)*dY;
100  y2 = _narrowYmin + double(J+1)*dY;
101  y12 = 0.5*(y1+y2);
102 
103  if( A_2 == 1 && A_1 != 1 ){
104  // pA, first beam is the nucleus and photon emitter
105  ega1 = 0.5*W*exp(y1);
106  ega2 = 0.5*W*exp(y2);
107  ega12 = 0.5*W*exp(y12);
108  beam = 2;
109  } else if( A_1 ==1 && A_2 != 1){
110  // pA, second beam is the nucleus and photon emitter
111  ega1 = 0.5*W*exp(-y1);
112  ega2 = 0.5*W*exp(-y2);
113  ega12 = 0.5*W*exp(-y12);
114  beam = 1;
115  } else {
116  ega1 = 0.5*W*exp(y1);
117  ega2 = 0.5*W*exp(y2);
118  ega12 = 0.5*W*exp(y12);
119  beam = 2;
120  }
121 
122  // This is for checking things in the lab frame
123  // y1lab = _narrowYmin + double(J)*dY;
124  // y2lab = _narrowYmin + double(J+1)*dY;
125  // y12lab = 0.5*(y1lab+y2lab);
126 
127  // p+Pb
128  // y1 = y1lab + 0.465;
129  // y2 = y2lab + 0.465;
130  // y12 = y12lab + 0.465;
131  // ega1 = 0.5*W*exp(y1);
132  // ega2 = 0.5*W*exp(y2);
133  // ega12 = 0.5*W*exp(y12);
134 
135  // Pb+p
136  // y1 = y1lab - 0.465;
137  // y2 = y2lab - 0.465;
138  // y12 = y12lab - 0.465;
139  // ega1 = 0.5*W*exp(-y1);
140  // ega2 = 0.5*W*exp(-y2);
141  // ega12 = 0.5*W*exp(-y12);
142 
143 
144  if(ega1 < Eth)
145  continue;
146  if(ega2 > maxPhotonEnergy())
147  continue;
148 
149  // First point
150  Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
151  +starlightConstants::protonMass*starlightConstants::protonMass);
152  csVN = sigma_N(Wgp);
153  csVA = sigma_A(csVN,beam);
154  csgA1 = (csVA/csVN)*sigmagp(Wgp);
155  if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
156  csgA1 = sigmagp(Wgp);
157  }
158 
159  // Middle point
160  Wgp = sqrt(2.*ega12*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
161  +starlightConstants::protonMass*starlightConstants::protonMass);
162  csVN = sigma_N(Wgp);
163  csVA = sigma_A(csVN,beam);
164  csgA12 = (csVA/csVN)*sigmagp(Wgp);
165  if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
166  csgA12 = sigmagp(Wgp);
167  }
168 
169  // Last point
170  Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
171  +starlightConstants::protonMass*starlightConstants::protonMass);
172  csVN = sigma_N(Wgp);
173  csVA = sigma_A(csVN,beam);
174  csgA2 = (csVA/csVN)*sigmagp(Wgp);
175  if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
176  csgA2 = sigmagp(Wgp);
177  }
178 
179  dR = ega1*photonFlux(ega1,beam)*csgA1;
180  dR = dR + 4*ega12*photonFlux(ega12,beam)*csgA12;
181  dR = dR + ega2*photonFlux(ega2,beam)*csgA2;
182  dR = dR*(dY/6.);
183 
184  // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
185  // cout<<" y: "<<y12lab<<" egamma: "<<ega12<<" flux: "<<ega12*photonFlux(ega12)<<" W: "<<Wgpm<<" Wflux: "<<Wgpm*photonFlux(ega12)<<" sigma_gA (nb): "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*ega12*photonFlux(ega12)*csgA12<<endl;
186 
187  int_r = int_r+dR;
188 
189  }
190 
191  // Repeat the loop for the case when the second beam is the photon emitter.
192  // Don't repeat for pA
193  if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
194  for(J=0;J<=(NY-1);J++){
195 
196  // This is the fdefault
197  y1 = _narrowYmin + double(J)*dY;
198  y2 = _narrowYmin + double(J+1)*dY;
199  y12 = 0.5*(y1+y2);
200 
201  beam = 1;
202  ega1 = 0.5*W*exp(-y1);
203  ega2 = 0.5*W*exp(-y2);
204  ega12 = 0.5*W*exp(-y12);
205 
206  // This is for checking things in the lab frame
207  // y1lab = _narrowYmin + double(J)*dY;
208  // y2lab = _narrowYmin + double(J+1)*dY;
209  // y12lab = 0.5*(y1lab+y2lab);
210 
211  // p+Pb
212  // y1 = y1lab + 0.465;
213  // y2 = y2lab + 0.465;
214  // y12 = y12lab + 0.465;
215  // ega1 = 0.5*W*exp(y1);
216  // ega2 = 0.5*W*exp(y2);
217  // ega12 = 0.5*W*exp(y12);
218 
219  // Pb+p
220  // y1 = y1lab - 0.465;
221  // y2 = y2lab - 0.465;
222  // y12 = y12lab - 0.465;
223  // ega1 = 0.5*W*exp(-y1);
224  // ega2 = 0.5*W*exp(-y2);
225  // ega12 = 0.5*W*exp(-y12);
226 
227 
228  if(ega2 < Eth)
229  continue;
230  if(ega1 > maxPhotonEnergy())
231  continue;
232 
233  // First point
234  Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
235  +starlightConstants::protonMass*starlightConstants::protonMass);
236  csVN = sigma_N(Wgp);
237  csVA = sigma_A(csVN,beam);
238  csgA1 = (csVA/csVN)*sigmagp(Wgp);
239  if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
240  csgA1 = sigmagp(Wgp);
241  }
242 
243  // Middle point
244  Wgp = sqrt(2.*ega12*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
245  +starlightConstants::protonMass*starlightConstants::protonMass);
246  csVN = sigma_N(Wgp);
247  csVA = sigma_A(csVN,beam);
248  csgA12 = (csVA/csVN)*sigmagp(Wgp);
249  if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
250  csgA12 = sigmagp(Wgp);
251  }
252 
253  // Last point
254  Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
255  +starlightConstants::protonMass*starlightConstants::protonMass);
256  csVN = sigma_N(Wgp);
257  csVA = sigma_A(csVN,beam);
258  csgA2 = (csVA/csVN)*sigmagp(Wgp);
259  if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
260  csgA2 = sigmagp(Wgp);
261  }
262 
263  dR = ega1*photonFlux(ega1,beam)*csgA1;
264  dR = dR + 4*ega12*photonFlux(ega12,beam)*csgA12;
265  dR = dR + ega2*photonFlux(ega2,beam)*csgA2;
266  dR = dR*(dY/6.);
267 
268  // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
269  // cout<<" y: "<<y12lab<<" egamma: "<<ega12<<" flux: "<<ega12*photonFlux(ega12)<<" W: "<<Wgpm<<" Wflux: "<<Wgpm*photonFlux(ega12)<<" sigma_gA (nb): "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*ega12*photonFlux(ega12)*csgA12<<endl;
270 
271  int_r = int_r+dR;
272 
273  }
274  }
275 
276  cout<<endl;
277  if (0.01*int_r > 1.){
278  cout<< " Total cross section: "<<0.01*int_r<<" barn."<<endl;
279  } else if (10.*int_r > 1.){
280  cout<< " Total cross section: " <<10.*int_r<<" mb."<<endl;
281  } else if (10000.*int_r > 1.){
282  cout<< " Total cross section: " <<10000.*int_r<<" microb."<<endl;
283  } else if (10000000.*int_r > 1.){
284  cout<< " Total cross section: " <<10000000.*int_r<<" nanob."<<endl;
285  } else if (1.E10*int_r > 1.){
286  cout<< " Total cross section: "<<1.E10*int_r<<" picob."<<endl;
287  } else {
288  cout<< " Total cross section: " <<1.E13*int_r<<" femtob."<<endl;
289  }
290  cout<<endl;
291  setPhotonNucleusSigma(0.01*int_r);
292 }