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:: 193                         $: revision of last commit
  24 // $Author:: jnystrand                $: author of last commit
  25 // $Date:: 2014-12-01 20:39:46 +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 "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 W,dW,dY;
  73         double y1,y2,y12,ega1,ega2,ega12;
  74         double csgA1,csgA2,csgA12,int_r,dR;
  75         double dsigdW,dsigdWalt,dndW;
  76         double dsigdW2;
  77         double Eth;
  78         int    I,J,NW,NY,beam;
  79
  80         double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
  81
  82         cout<<" in wideResonanceCrossSection, bwnorm: "<<bwnorm<<endl;
  83         //gamma+nucleon threshold.
  84         Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
  85                   -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
  86
  87         NW   = 100;
  88         dW   = (_wideWmax-_wideWmin)/double(NW);
  89
  90         NY   =  1200;
  91         dY   = (_wideYmax-_wideYmin)/double(NY);
  92
  93         if (getBNORM()  ==  0.){
  94                 cout<<" Using Breit-Wigner Resonance Profile."<<endl;
  95         }
  96         else{
  97                 cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
  98         }
  99
 100         cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
 101
 102         int A_1 = getbbs().beam1().A();
 103         int A_2 = getbbs().beam2().A();
 104
 105         int_r=0.;
 106
 107         // Do this first for the case when the first beam is the photon emitter
 108         // Treat pA separately with defined beams
 109         // The variable beam (=1,2) defines which nucleus is the target
 110         for(I=0;I<=NW-1;I++){
 111
 112                 W = _wideWmin + double(I)*dW + 0.5*dW;
 113
 114                 dsigdW=0.0;
 115                 dsigdW2=0.0;
 116                 dsigdWalt=0.0;
 117                 dndW=0.0;
 118
 119                 for(J=0;J<=NY-1;J++){
 120
 121                         y1  = _wideYmin + double(J)*dY;
 122                         y2  = _wideYmin + double(J+1)*dY;
 123                         y12 = 0.5*(y1+y2);
 124
 125                         if( A_2 == 1 && A_1 != 1 ){
 126                           // pA, first beam is the nucleus and is in this case the target
 127                           // Egamma = 0.5*W*exp(-Y);
 128                           ega1  = 0.5*W*exp(-y1);
 129                           ega2  = 0.5*W*exp(-y2);
 130                           ega12 = 0.5*W*exp(-y12);
 131                           beam = 1;
 132                         } else if( A_1 ==1 && A_2 != 1){
 133                           // pA, second beam is the nucleus and is in this case the target
 134                           // Egamma = 0.5*W*exp(Y);
 135                           ega1  = 0.5*W*exp(y1);
 136                           ega2  = 0.5*W*exp(y2);
 137                           ega12 = 0.5*W*exp(y12);
 138                           beam = 2;
 139                         } else {
 140                           // Egamma = 0.5*W*exp(Y);
 141                           ega1  = 0.5*W*exp(y1);
 142                           ega2  = 0.5*W*exp(y2);
 143                           ega12 = 0.5*W*exp(y12);
 144                           beam = 2;
 145                         }
 146
 147                         //                      ega1  = 0.5*W*exp(y1);
 148                         //      ega2  = 0.5*W*exp(y2);
 149                         //      ega12 = 0.5*W*exp(y12);
 150
 151                         if(ega1 < Eth) continue;
 152                         if(ega2 > maxPhotonEnergy()) continue;
 153
 154                         csgA1=getcsgA(ega1,W,beam);
 155
 156                         //         >> Middle Point                      =====>>>
 157                         csgA12=getcsgA(ega12,W,beam);
 158
 159                         //         >> Second Point                      =====>>>
 160                         csgA2=getcsgA(ega2,W,beam);
 161
 162                         //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
 163                         dR  = ega1*photonFlux(ega1,beam)*csgA1;
 164                         dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
 165                         dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
 166                         dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
 167
 168                         // cout<<" W: "<<W<<" Y: "<<y12<<" dR: "<<dR<<endl;
 169                         //For identical beams, we double.  Either may emit photon/scatter
 170                         //For large differences in Z, we approx, that only beam1 emits photon
 171                         //and beam2 scatters, eg d-Au where beam1=au and beam2=d
 172                         //if(getbbs().beam1().A()==getbbs().beam2().A()){
 173                         //      dR  = 2.*dR;
 174                         //}
 175                         int_r = int_r+dR;
 176                 }
 177         }
 178
 179         // Repeat the loop for the case when the second beam is the photon emitter.
 180         // Don't repeat for pA
 181         if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
 182           for(I=0;I<=NW-1;I++){
 183
 184                 W = _wideWmin + double(I)*dW + 0.5*dW;
 185
 186                 dsigdW=0.0;
 187                 dsigdW2=0.0;
 188                 dsigdWalt=0.0;
 189                 dndW=0.0;
 190
 191                 for(J=0;J<=NY-1;J++){
 192
 193                         y1  = _wideYmin + double(J)*dY;
 194                         y2  = _wideYmin + double(J+1)*dY;
 195                         y12 = 0.5*(y1+y2);
 196
 197                         beam = 1;
 198                         ega1  = 0.5*W*exp(-y1);
 199                         ega2  = 0.5*W*exp(-y2);
 200                         ega12 = 0.5*W*exp(-y12);
 201
 202                         if(ega2 < Eth) continue;
 203                         if(ega1 > maxPhotonEnergy()) continue;
 204
 205                         csgA1=getcsgA(ega1,W,beam);
 206
 207                         //         >> Middle Point                      =====>>>
 208                         csgA12=getcsgA(ega12,W,beam);
 209
 210                         //         >> Second Point                      =====>>>
 211                         csgA2=getcsgA(ega2,W,beam);
 212
 213                         //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
 214                         dR  = ega1*photonFlux(ega1,beam)*csgA1;
 215                         dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
 216                         dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
 217                         dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
 218
 219                         //For identical beams, we double.  Either may emit photon/scatter
 220                         //For large differences in Z, we approx, that only beam1 emits photon
 221                         //and beam2 scatters, eg d-Au where beam1=au and beam2=d
 222                         // if(getbbs().beam1().A()==getbbs().beam2().A()){
 223                         //      dR  = 2.*dR;
 224                         // }
 225                         int_r = int_r+dR;
 226                 }
 227           }
 228         }
 229         cout<<" Cross section (mb): "<<10.*int_r<<endl;
 230
 231 }