STARLIGHT code and interface

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

diff --git a/STARLIGHT/starlight/src/narrowResonanceCrossSection.cpp b/STARLIGHT/starlight/src/narrowResonanceCrossSection.cpp
new file mode 100644 (file)
index 0000000..f238c5c
--- /dev/null
+++ b/STARLIGHT/starlight/src/narrowResonanceCrossSection.cpp
@@ -0,0 +1,134 @@
+///////////////////////////////////////////////////////////////////////////
+//
+//    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 "narrowResonanceCrossSection.h"
+
+
+using namespace std;
+using namespace starlightConstants;
+
+
+//______________________________________________________________________________
+narrowResonanceCrossSection::narrowResonanceCrossSection(const beamBeamSystem&  bbsystem)
+       :photonNucleusCrossSection(bbsystem)
+{
+       _narrowYmax = inputParametersInstance.maxRapidity();
+       _narrowYmin = -1.0*_narrowYmax;
+       _narrowNumY = inputParametersInstance.nmbRapidityBins();
+       _Ep         = inputParametersInstance.protonEnergy();   
+}
+
+
+//______________________________________________________________________________
+narrowResonanceCrossSection::~narrowResonanceCrossSection()
+{ }
+
+
+//______________________________________________________________________________
+void
+narrowResonanceCrossSection::crossSectionCalculation(const double)  // _bwnormsave (unused)
+{
+       // This subroutine calculates the vector meson cross section assuming
+       // a narrow resonance.  For reference, see STAR Note 386.
+  
+       // double Av,Wgp,cs,cvma;
+       double W,dY;
+       double y1,y2,y12,ega1,ega2,ega12;
+       // double t,tmin,tmax;
+       double csgA1,csgA2,csgA12,int_r,dR,rate;
+       double tmp;
+       // double ax,bx;
+       double Eth;
+       int    J,NY;
+       // int    K,NGAUSS;
+  
+       NY   =  _narrowNumY;
+       dY   = (_narrowYmax-_narrowYmin)/double(NY);
+  
+       cout<<" Using Narrow Resonance ..."<<endl;
+  
+       W = getChannelMass();
+       Eth=0.5*(((W+protonMass)*(W+protonMass)-
+                 protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
+  
+       cout<<" gamma+nucleon  Threshold: "<<Eth<<endl;
+       int_r=0.;
+  
+       tmp = 0.0;
+  
+       for(J=0;J<=(NY-1);J++){
+    
+               y1  = _narrowYmin + double(J)*dY;
+               y2  = _narrowYmin + 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;
+
+               csgA1=getcsgA(ega1,W);
+    
+               // Middle Point                      =====>>>
+               csgA12=getcsgA(ega12,W); 
+
+               // Second Point                      =====>>>
+               csgA2=getcsgA(ega2,W);
+    
+               // Sum the contribution for this W,Y.
+               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.);
+
+               // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
+
+               // The 2 accounts for the 2 beams
+               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;
+}