STARLIGHT/starlight/src/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 //
  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 "starlightconstants.h"
  39 #include "incoherentVMCrossSection.h"
  40
  41
  42 using namespace std;
  43 using namespace starlightConstants;
  44
  45
  46 //______________________________________________________________________________
  47 incoherentVMCrossSection::incoherentVMCrossSection(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 //______________________________________________________________________________
  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 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 Wgp,csVN,csVA;
  75         double tmp;
  76         // double ax,bx;
  77         double Eth;
  78         int    J,NY;
  79         // int    K,NGAUSS;
  80
  81         NY   =  _narrowNumY;
  82         dY   = (_narrowYmax-_narrowYmin)/double(NY);
  83
  84         cout<<" Using Narrow Resonance ..."<<endl;
  85
  86         W = getChannelMass();
  87         Eth=0.5*(((W+protonMass)*(W+protonMass)-
  88                   protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
  89
  90         cout<<" gamma+nucleon  Threshold: "<<Eth<<endl;
  91         int_r=0.;
  92
  93         tmp = 0.0;
  94
  95         for(J=0;J<=(NY-1);J++){
  96
  97                 y1  = _narrowYmin + double(J)*dY;
  98                 y2  = _narrowYmin + double(J+1)*dY;
  99                 y12 = 0.5*(y1+y2);
 100
 101                 ega1  = 0.5*W*exp(y1);
 102                 ega2  = 0.5*W*exp(y2);
 103                 ega12 = 0.5*W*exp(y12);
 104
 105                 if(ega1 < Eth)
 106                         continue;
 107                 if(ega2 > maxPhotonEnergy())
 108                         continue;
 109
 110                 // First point
 111                 Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
 112                                   +starlightConstants::protonMass*starlightConstants::protonMass);
 113                 csVN = sigma_N(Wgp);
 114                 csVA = sigma_A(csVN);
 115                 csgA1 = (csVA/csVN)*sigmagp(Wgp);
 116                 if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
 117                   csgA1 = sigmagp(Wgp);
 118                 }
 119
 120                 // Middle point
 121                 Wgp = sqrt(2.*ega12*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
 122                                   +starlightConstants::protonMass*starlightConstants::protonMass);
 123                 csVN = sigma_N(Wgp);
 124                 csVA = sigma_A(csVN);
 125                 csgA12 = (csVA/csVN)*sigmagp(Wgp);
 126                 if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
 127                   csgA12 = sigmagp(Wgp);
 128                 }
 129
 130                 // Last point
 131                 Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
 132                                   +starlightConstants::protonMass*starlightConstants::protonMass);
 133                 csVN = sigma_N(Wgp);
 134                 csVA = sigma_A(csVN);
 135                 csgA2 = (csVA/csVN)*sigmagp(Wgp);
 136                 if( getbbs().beam1().A() == 1 || getbbs().beam2().A()==1 ){
 137                   csgA2 = sigmagp(Wgp);
 138                 }
 139
 140                 dR = ega1*photonFlux(ega1)*csgA1;
 141                 dR = dR + 4*ega12*photonFlux(ega12)*csgA12;
 142                 dR = dR + ega2*photonFlux(ega2)*csgA2;
 143                 dR = dR*(dY/6.);
 144
 145                 // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
 146
 147                 // The 2 accounts for the 2 beams
 148                 if(getbbs().beam1().A()==getbbs().beam2().A()){
 149                         dR  = 2.*dR;
 150                 }
 151
 152                 int_r = int_r+dR;
 153
 154         }
 155         rate=luminosity()*int_r;
 156         cout<<" Cross section (mb): " <<10.*int_r<<endl;
 157 }