]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/.svn/text-base/narrowResonanceCrossSection.cpp.svn-base
STARLIGHT code and interface
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / .svn / text-base / narrowResonanceCrossSection.cpp.svn-base
CommitLineData
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:: $: revision of last commit
24// $Author:: $: author of last commit
25// $Date:: $: 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 "starlightconstants.h"
39#include "narrowResonanceCrossSection.h"
40
41
42using namespace std;
43using namespace starlightConstants;
44
45
46//______________________________________________________________________________
47narrowResonanceCrossSection::narrowResonanceCrossSection(const beamBeamSystem& bbsystem)
48 :photonNucleusCrossSection(bbsystem)
49{
50 _narrowYmax = inputParametersInstance.maxRapidity();
51 _narrowYmin = -1.0*_narrowYmax;
52 _narrowNumY = inputParametersInstance.nmbRapidityBins();
53 _Ep = inputParametersInstance.protonEnergy();
54}
55
56
57//______________________________________________________________________________
58narrowResonanceCrossSection::~narrowResonanceCrossSection()
59{ }
60
61
62//______________________________________________________________________________
63void
64narrowResonanceCrossSection::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 Av,Wgp,cs,cvma;
70 double W,dY;
71 double y1,y2,y12,ega1,ega2,ega12;
72 // double t,tmin,tmax;
73 double csgA1,csgA2,csgA12,int_r,dR,rate;
74 double tmp;
75 // double ax,bx;
76 double Eth;
77 int J,NY;
78 // int K,NGAUSS;
79
80 NY = _narrowNumY;
81 dY = (_narrowYmax-_narrowYmin)/double(NY);
82
83 cout<<" Using Narrow Resonance ..."<<endl;
84
85 W = getChannelMass();
86 Eth=0.5*(((W+protonMass)*(W+protonMass)-
87 protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
88
89 cout<<" gamma+nucleon Threshold: "<<Eth<<endl;
90 int_r=0.;
91
92 tmp = 0.0;
93
94 for(J=0;J<=(NY-1);J++){
95
96 y1 = _narrowYmin + double(J)*dY;
97 y2 = _narrowYmin + double(J+1)*dY;
98 y12 = 0.5*(y1+y2);
99
100 ega1 = 0.5*W*exp(y1);
101 ega2 = 0.5*W*exp(y2);
102 ega12 = 0.5*W*exp(y12);
103
104 if(ega1 < Eth)
105 continue;
106 if(ega2 > maxPhotonEnergy())
107 continue;
108
109 csgA1=getcsgA(ega1,W);
110
111 // Middle Point =====>>>
112 csgA12=getcsgA(ega12,W);
113
114 // Second Point =====>>>
115 csgA2=getcsgA(ega2,W);
116
117 // Sum the contribution for this W,Y.
118 dR = ega1*photonFlux(ega1)*csgA1;
119 dR = dR + 4.*ega12*photonFlux(ega12)*csgA12;
120 dR = dR + ega2*photonFlux(ega2)*csgA2;
121 tmp = tmp+2.*dR*(dY/6.);
122 dR = dR*(dY/6.);
123
124 // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
125
126 // The 2 accounts for the 2 beams
127 if(getbbs().beam1().A()==getbbs().beam2().A()){
128 dR = 2.*dR;
129 }
130 int_r = int_r+dR;
131 }
132 rate=luminosity()*int_r;
133 cout<<" Cross section (mb): " <<10.*int_r<<endl;
134}