STARLIGHT code and interface

[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / wideResonanceCrossSection.cpp

diff --git a/STARLIGHT/starlight/src/wideResonanceCrossSection.cpp b/STARLIGHT/starlight/src/wideResonanceCrossSection.cpp
new file mode 100644 (file)
index 0000000..e538c4d
--- /dev/null
+++ b/STARLIGHT/starlight/src/wideResonanceCrossSection.cpp
@@ -0,0 +1,171 @@
+///////////////////////////////////////////////////////////////////////////
+//
+//    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(0),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;
+}