+++ /dev/null
-///////////////////////////////////////////////////////////////////////////
-//
-// Copyright 2010
-//
-// This file is part of starlight.
-//
-// starlight is free software: you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation, either version 3 of the License, or
-// (at your option) any later version.
-//
-// starlight is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 164 $: revision of last commit
-// $Author:: odjuvsla $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-#include <cmath>
-
-#include "starlightconstants.h"
-#include "wideResonanceCrossSection.h"
-
-
-using namespace std;
-using namespace starlightConstants;
-
-
-//______________________________________________________________________________
-wideResonanceCrossSection::wideResonanceCrossSection(const beamBeamSystem& bbsystem)
- : photonNucleusCrossSection(bbsystem)//hrm
-{
- _wideWmax = _wMax;
- _wideWmin = _wMin;
- _wideYmax = _yMax;
- _wideYmin = -1.0 * _wideYmax;
- _Ep = inputParametersInstance.protonEnergy();
-}
-
-
-//______________________________________________________________________________
-wideResonanceCrossSection::~wideResonanceCrossSection()
-{
-
-}
-
-
-//______________________________________________________________________________
-void
-wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
-{
- // This subroutine calculates the cross-section assuming a wide
- // (Breit-Wigner) resonance.
-
- // double Av,Wgp,cs,cvma;
- double W,dW,dY;
- double y1,y2,y12,ega1,ega2,ega12;
- // double t,tmin,tmax;
- double csgA1,csgA2,csgA12,int_r,dR,rate;
- double dsigdW,dsigdWalt,dndW,tmp;
- double dsigdW2;
- // double ax,bx;
- double Eth;
- int I,J,NW,NY;
- // int K,NGAUSS;
-
- // ----------------- !!!!!!!!!!!!!!!!!!!! -----------------------------
-
- double bwnorm =bwnormsave;//used to transfer the bwnorm from the luminosity tables
-
- // --------------------------------------------------------------------
- //gamma+nucleon threshold.
-
- Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
- -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
-
- NW = 100;
- dW = (_wideWmax-_wideWmin)/double(NW);
-
- NY = 1200;
- dY = (_wideYmax-_wideYmin)/double(NY);
-
- if (getBNORM() == 0.){
- cout<<" Using Breit-Wigner Resonance Profile."<<endl;
- }
- else{
- cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
- }
-
- cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
-
- int_r=0.;
-
- for(I=0;I<=NW-1;I++){
-
- W = _wideWmin + double(I)*dW + 0.5*dW;
-
- tmp = 0.0;
- dsigdW=0.0;
- dsigdW2=0.0;
- dsigdWalt=0.0;
- dndW=0.0;
-
- for(J=0;J<=NY-1;J++){
-
- y1 = _wideYmin + double(J)*dY;
- y2 = _wideYmin + double(J+1)*dY;
- y12 = 0.5*(y1+y2);
-
- ega1 = 0.5*W*exp(y1);
- ega2 = 0.5*W*exp(y2);
- ega12 = 0.5*W*exp(y12);
-
- if(ega1 < Eth) continue;
- if(ega2 > maxPhotonEnergy()) continue;
- // check it !!
-
- if(J == 0){
- // >> 1st Point (Calculated only the first time) =====>>>
- //ega1 used.
- csgA1=getcsgA(ega1,W);
- }
- else{
- csgA1 = csgA2;
- }
-
- // >> Middle Point =====>>>
- csgA12=getcsgA(ega12,W);
-
- // >> Second Point =====>>>
- csgA2=getcsgA(ega2,W);
-
- //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
- dR = ega1*photonFlux(ega1)*csgA1;
- dR = dR + 4.*ega12*photonFlux(ega12)*csgA12;
- dR = dR + ega2*photonFlux(ega2)*csgA2;
- tmp = tmp+2.*dR*(dY/6.);
- dR = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
-
- //For identical beams, we double. Either may emit photon/scatter
- //For large differences in Z, we approx, that only beam1 emits photon
- //and beam2 scatters, eg d-Au where beam1=au and beam2=d
- if(getbbs().beam1().A()==getbbs().beam2().A()){
- dR = 2.*dR;
- }
- int_r = int_r+dR;
- }
- }
-
- rate=luminosity()*int_r;
-
- cout<<" Cross section (mb): "<<10.*int_r<<endl;
- cout<<" Production rate : "<<rate<<" Hz"<<endl;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff