STARLIGHT/starlight/src/wideResonanceCrossSection.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:: 164                         $: revision of last commit
  24 // $Author:: odjuvsla                 $: author of last commit
  25 // $Date:: 2013-10-06 16:18:08 +0200 #$: 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 "wideResonanceCrossSection.h"
  40
  41
  42 using namespace std;
  43 using namespace starlightConstants;
  44
  45
  46 //______________________________________________________________________________
  47 wideResonanceCrossSection::wideResonanceCrossSection(const beamBeamSystem&  bbsystem)
  48         : photonNucleusCrossSection(bbsystem)//hrm
  49 {
  50         _wideWmax = _wMax;
  51         _wideWmin = _wMin;
  52         _wideYmax = _yMax;
  53         _wideYmin = -1.0 * _wideYmax;
  54         _Ep       = inputParametersInstance.protonEnergy();
  55 }
  56
  57
  58 //______________________________________________________________________________
  59 wideResonanceCrossSection::~wideResonanceCrossSection()
  60 {
  61
  62 }
  63
  64
  65 //______________________________________________________________________________
  66 void
  67 wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
  68 {
  69         //     This subroutine calculates the cross-section assuming a wide
  70         //     (Breit-Wigner) resonance.
  71
  72         // double Av,Wgp,cs,cvma;
  73         double W,dW,dY;
  74         double y1,y2,y12,ega1,ega2,ega12;
  75         // double t,tmin,tmax;
  76         double csgA1,csgA2,csgA12,int_r,dR,rate;
  77         double dsigdW,dsigdWalt,dndW,tmp;
  78         double dsigdW2;
  79         // double ax,bx;
  80         double Eth;
  81         int    I,J,NW,NY;
  82         // int    K,NGAUSS;
  83
  84         // ----------------- !!!!!!!!!!!!!!!!!!!! -----------------------------
  85
  86         double bwnorm =bwnormsave;//used to transfer the bwnorm from the luminosity tables
  87
  88         // --------------------------------------------------------------------
  89         //gamma+nucleon threshold.
  90
  91         Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
  92                   -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
  93
  94         NW   = 100;
  95         dW   = (_wideWmax-_wideWmin)/double(NW);
  96
  97         NY   =  1200;
  98         dY   = (_wideYmax-_wideYmin)/double(NY);
  99
 100         if (getBNORM()  ==  0.){
 101                 cout<<" Using Breit-Wigner Resonance Profile."<<endl;
 102         }
 103         else{
 104                 cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
 105         }
 106
 107         cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
 108
 109         int_r=0.;
 110
 111         for(I=0;I<=NW-1;I++){
 112
 113                 W = _wideWmin + double(I)*dW + 0.5*dW;
 114
 115                 tmp = 0.0;
 116                 dsigdW=0.0;
 117                 dsigdW2=0.0;
 118                 dsigdWalt=0.0;
 119                 dndW=0.0;
 120
 121                 for(J=0;J<=NY-1;J++){
 122
 123                         y1  = _wideYmin + double(J)*dY;
 124                         y2  = _wideYmin + double(J+1)*dY;
 125                         y12 = 0.5*(y1+y2);
 126
 127                         ega1  = 0.5*W*exp(y1);
 128                         ega2  = 0.5*W*exp(y2);
 129                         ega12 = 0.5*W*exp(y12);
 130
 131                         if(ega1 < Eth) continue;
 132                         if(ega2 > maxPhotonEnergy()) continue;
 133                         // check it !!
 134
 135                         if(J == 0){
 136                                 // >> 1st Point (Calculated only the first time)     =====>>>
 137                                 //ega1 used.
 138                                 csgA1=getcsgA(ega1,W);
 139                         }
 140                         else{
 141                                 csgA1 = csgA2;
 142                         }
 143
 144                         //         >> Middle Point                      =====>>>
 145                         csgA12=getcsgA(ega12,W);
 146
 147                         //         >> Second Point                      =====>>>
 148                         csgA2=getcsgA(ega2,W);
 149
 150                         //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
 151                         dR  = ega1*photonFlux(ega1)*csgA1;
 152                         dR  = dR + 4.*ega12*photonFlux(ega12)*csgA12;
 153                         dR  = dR + ega2*photonFlux(ega2)*csgA2;
 154                         tmp = tmp+2.*dR*(dY/6.);
 155                         dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
 156
 157                         //For identical beams, we double.  Either may emit photon/scatter
 158                         //For large differences in Z, we approx, that only beam1 emits photon
 159                         //and beam2 scatters, eg d-Au where beam1=au and beam2=d
 160                         if(getbbs().beam1().A()==getbbs().beam2().A()){
 161                                 dR  = 2.*dR;
 162                         }
 163                         int_r = int_r+dR;
 164                 }
 165         }
 166
 167         rate=luminosity()*int_r;
 168
 169         cout<<" Cross section (mb): "<<10.*int_r<<endl;
 170         cout<<" Production rate   : "<<rate<<" Hz"<<endl;
 171 }