estarlight/src/e_wideResonanceCrossSection.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 264                         $: revision of last commit
0024 // $Author:: mlomnitz                   $: author of last commit
0025 // $Date:: 2017-03-14 21:05:12 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032 //#define _makeGammaPQ2_
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "starlightconstants.h"
0039 #include "e_wideResonanceCrossSection.h"
0040
0041
0042 using namespace std;
0043 using namespace starlightConstants;
0044
0045
0046 //______________________________________________________________________________
0047 e_wideResonanceCrossSection::e_wideResonanceCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem&  bbsystem)
0048     : photonNucleusCrossSection(inputParametersInstance, bbsystem)//hrm
0049 {
0050     _wideWmax = _wMax;
0051     _wideWmin = _wMin;
0052     _targetMaxPhotonEnergy=inputParametersInstance.targetMaxPhotonEnergy();
0053     _targetMinPhotonEnergy=inputParametersInstance.targetMinPhotonEnergy();
0054     _Ep       = inputParametersInstance.protonEnergy();
0055     _electronEnergy = inputParametersInstance.electronEnergy();
0056     _target_beamLorentz = inputParametersInstance.beamLorentzGamma();
0057     _VMnumEgamma = inputParametersInstance.nmbEnergyBins();
0058     _useFixedRange = inputParametersInstance.fixedQ2Range();
0059     _gammaMinQ2 = inputParametersInstance.minGammaQ2();
0060     _gammaMaxQ2 = inputParametersInstance.maxGammaQ2();
0061     _targetRadii = inputParametersInstance.targetRadius();
0062 }
0063
0064
0065 //______________________________________________________________________________
0066 e_wideResonanceCrossSection::~e_wideResonanceCrossSection()
0067 {
0068
0069 }
0070
0071
0072 //______________________________________________________________________________
0073 void
0074 e_wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
0075 {
0076     //     This subroutine calculates the cross-section assuming a wide
0077     //     (Breit-Wigner) resonance.
0078
0079   double W,dW, dEgamma, minEgamma;
0080     double ega[3] = {0};
0081     double int_r,dR;
0082     double int_r2,dR2;
0083     int    iW,nW,iEgamma,nEgamma,beam;
0084
0085     double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
0086
0087     // For W integration
0088     nW   = 100;
0089     dW   = (_wideWmax-_wideWmin)/double(nW);
0090     // For Egamma integration
0091     nEgamma = 1000;
0092     dEgamma = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/nEgamma;
0093     minEgamma = std::log(_targetMinPhotonEnergy);
0094
0095
0096     if (getBNORM()  ==  0.){
0097         cout<<" Using Breit-Wigner Resonance Profile."<<endl;
0098     }
0099     else{
0100         cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
0101     }
0102
0103     cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
0104
0105         int A_1 = getbbs().electronBeam().A();
0106         int A_2 = getbbs().targetBeam().A();
0107
0108     if( A_2 == 0 && A_1 >= 1 ){
0109           // eA, first beam is the nucleus and is in this case the target
0110           beam = 1;
0111         } else if( A_1 ==0 && A_2 >= 1){
0112       // eA, second beam is the nucleus and is in this case the target
0113           beam = 2;
0114         } else {
0115           beam = 2;
0116         }
0117
0118     int_r=0.;
0119     int_r2 = 0.;
0120         // Do this first for the case when the first beam is the photon emitter
0121         // The variable beam (=1,2) defines which nucleus is the target
0122     // Integration done using Simpson's rule
0123     for(iW=0;iW<=nW-1;iW++){
0124
0125         W = _wideWmin + double(iW)*dW + 0.5*dW;
0126         int nQ2 = 1000;
0127         for(iEgamma = 0 ; iEgamma < nEgamma; ++iEgamma){    // Integral over photon energy
0128           // Target frame photon energies
0129           ega[0] = exp(minEgamma + iEgamma*dEgamma );
0130           ega[1] = exp(minEgamma + (iEgamma+1)*dEgamma );
0131           ega[2] = 0.5*(ega[0]+ega[1]);
0132           // Integral over Q2
0133           double full_int[3] = {0}; // Full e+X --> e+X+V.M. cross section
0134           double dndE[3] = {0}; // Full e+X --> e+X+V.M. cross section
0135           //
0136           for( int iEgaInt = 0 ; iEgaInt < 3; ++iEgaInt){    // Loop over the energies for the three points to integrate over Q2
0137             //These are the physical limits
0138             double Q2_min = std::pow(starlightConstants::mel*ega[iEgaInt],2.0)/(_electronEnergy*(_electronEnergy-ega[iEgaInt]));
0139             double Q2_max = 4.*_electronEnergy*(_electronEnergy-ega[iEgaInt]);
0140             if(_useFixedRange == true){
0141               if( Q2_min < _gammaMinQ2 )
0142             Q2_min = _gammaMinQ2;
0143               if( Q2_max > _gammaMaxQ2 )
0144             Q2_max = _gammaMaxQ2;
0145             }
0146             double lnQ2ratio = std::log(Q2_max/Q2_min)/nQ2;
0147             double lnQ2_min = std::log(Q2_min);
0148             //
0149             for( int iQ2 = 0 ; iQ2 < nQ2; ++iQ2){     // Integral over photon virtuality
0150               //
0151               double q2_1 = exp( lnQ2_min + iQ2*lnQ2ratio);
0152               double q2_2 = exp( lnQ2_min + (iQ2+1)*lnQ2ratio);
0153               double q2_12 = (q2_2+q2_1)/2.;
0154               //
0155               // Integrating cross section
0156               full_int[iEgaInt] += (q2_2-q2_1)*( g(ega[iEgaInt],q2_1)*getcsgA(ega[iEgaInt],q2_1,beam)
0157                              + g(ega[iEgaInt],q2_2)*getcsgA(ega[iEgaInt],q2_2,beam)
0158                              + 4.*g(ega[iEgaInt],q2_12)*getcsgA(ega[iEgaInt],q2_12,beam) );
0159               // Effective flux
0160               dndE[iEgaInt] +=(q2_2-q2_1)*( getcsgA_Q2_dep(q2_1)*photonFlux(ega[iEgaInt],q2_1)
0161                             +getcsgA_Q2_dep(q2_2)*photonFlux(ega[iEgaInt],q2_2)
0162                             +4.*getcsgA_Q2_dep(q2_12)*photonFlux(ega[iEgaInt],q2_12) );
0163             }
0164             full_int[iEgaInt] = full_int[iEgaInt]/6.;
0165             dndE[iEgaInt] = dndE[iEgaInt]/6.;
0166           }
0167           // Finishing cross-section integral
0168           dR = full_int[0];
0169           dR += full_int[1];
0170           dR += 4.*full_int[2];
0171           dR = dR*(ega[1]-ega[0])/6.;
0172           int_r = int_r + dR*breitWigner(W,bwnorm)*dW;
0173           // Finishing effective flux integral
0174               // Finishing integral over the effective photon flux
0175           dR2 = dndE[0];
0176           dR2 += dndE[1];
0177           dR2 += 4.*dndE[2];
0178           dR2 = dR2*(ega[1]-ega[0])/6.;
0179           int_r2 = int_r2 + dR2*breitWigner(W,bwnorm)*dW;
0180         }
0181     }
0182     cout<<endl;
0183     if(_useFixedRange == true){
0184       cout<<" Using fixed Q2 range "<<_gammaMinQ2 << " < Q2 < "<<_gammaMaxQ2<<endl;
0185     }
0186     printCrossSection(" Total cross section: ",int_r);
0187     //printCrossSection(" gamma+X --> VM+X ", int_r/int_r2);
0188     setPhotonNucleusSigma(0.01*int_r);
0189     //
0190 #ifdef _makeGammaPQ2_
0191     makeGammaPQ2dependence(bwnormsave);
0192 #endif
0193 }
0194
0195 #ifdef _makeGammaPQ2_
0196 //______________________________________________________________________________
0197 void
0198 e_wideResonanceCrossSection::makeGammaPQ2dependence( double bwnormsave)
0199 {
0200     // This subroutine calculates the Q2 dependence of
0201         // gamma+X -> VM + X cross section for a narrow resonance
0202
0203         int const nQ2bins = 19;
0204     double const q2Edge[nQ2bins+1] = { 0.,1.,2.,3., 4.,5.,
0205                        6.,7.,8.,9.,10.,
0206                        11.,12.,13.,14.,15.,
0207                        20.,30.,40.,50.};
0208     double W = 0,dW,dY;
0209     double y1,y2,y12,ega1,ega2,ega12;
0210     double int_r,dR,dR2;
0211     double csgA1, csgA2, csgA12;
0212     double Eth;
0213     int    I,J,NW,NY,beam;
0214
0215     double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
0216
0217     NW   = 100;
0218     dW   = (_wideWmax-_wideWmin)/double(NW);
0219
0220     if (getBNORM()  ==  0.){
0221         cout<<" Using Breit-Wigner Resonance Profile."<<endl;
0222     }
0223     else{
0224         cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
0225     }
0226
0227     cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
0228
0229     //Lomnitz old used for XX
0230     Eth=0.5*(((W+protonMass)*(W+protonMass)-
0231               protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
0232     // Adapted for eX
0233     //Eth=0.5*(((W+starlightConstants::mel)*(W +starlightConstants::mel)-
0234     //    starlightConstants::mel*starlightConstants::mel)/(_electronEnergy+sqrt(_electronEnergy*_electronEnergy-starlightConstants::mel*starlightConstants::mel)));
0235
0236         printf(" gamma+nucleon threshold: %e GeV \n", Eth);
0237
0238         int A_1 = getbbs().electronBeam().A();
0239         int A_2 = getbbs().targetBeam().A();
0240
0241
0242         // Do this first for the case when the first beam is the photon emitter
0243         // The variable beam (=1,2) defines which nucleus is the target
0244     // Target beam ==2 so repidity is negative. Can generalize later
0245     ///
0246     cout<<" Lomnitz debug :: sigma_gamma_p --> VM_p "<<endl;
0247     cout<<" Q2+MV2 \t \t"<<" sigma_gamma_p --> VM_p (nanob)"<<endl;
0248     double target_cm = acosh(_target_beamLorentz);
0249     // another - sign from subraction in addition rule
0250     double exp_target_cm = exp(-target_cm);
0251     double int_r2;
0252     for( int iQ2 = 0 ; iQ2 < nQ2bins; ++iQ2){
0253       int_r=0.;
0254       int_r2=0.;
0255       //
0256       double q2_cor = getcsgA_Q2_dep( (q2Edge[iQ2+1] + q2Edge[iQ2])/2. );
0257       for(I=0;I<=NW-1;I++){
0258
0259         W = _wideWmin + double(I)*dW + 0.5*dW;
0260         for(J=0;J<=(NY-1);J++){
0261           y1  = _wideYmin + double(J)*dY;
0262           y2  = _wideYmin + double(J+1)*dY;
0263           y12 = 0.5*(y1+y2);
0264           double target_ega1, target_ega2, target_ega12;
0265           if( A_2 == 0 && A_1 >= 1 ){
0266         // pA, first beam is the nucleus and is in this case the target
0267         ega1  = 0.5*W*exp(-y1);
0268         ega2  = 0.5*W*exp(-y2);
0269         ega12 = 0.5*W*exp(-y12);
0270         beam = 1;
0271           } else if( A_1 ==0 && A_2 >= 1){
0272         // pA, second beam is the nucleus and is in this case the target
0273         ega1  = 0.5*W*exp(y1);
0274         ega2  = 0.5*W*exp(y2);
0275         ega12 = 0.5*W*exp(y12);
0276         // photon energy in Target frame
0277         beam = 2;
0278           } else {
0279         ega1  = 0.5*W*exp(y1);
0280         ega2  = 0.5*W*exp(y2);
0281         ega12 = 0.5*W*exp(y12);
0282         // photon energy in Target frame
0283         beam = 2;
0284           }
0285           //
0286           if(ega1 < Eth || ega2 < Eth)
0287         continue;
0288           if(ega2 > maxPhotonEnergy() || ega1 > maxPhotonEnergy() )
0289         continue;
0290           target_ega1 = ega1*exp_target_cm;
0291           target_ega12 = ega12*exp_target_cm;
0292           target_ega2 = ega2*exp_target_cm;
0293           //cout<<"Nortmalizations "<<integrated_x_section(ega1,0,50)<<endl;
0294           //
0295           csgA1=getcsgA(ega1,W,beam);
0296           double full_range_1 = integrated_x_section(target_ega1);
0297           //         >> Middle Point                      =====>>>
0298           csgA12=getcsgA(ega12,W,beam);
0299           double full_range_12 = integrated_x_section(target_ega12);
0300           //         >> Second Point                      =====>>>
0301           csgA2=getcsgA(ega2,W,beam);
0302           double full_range_2 = integrated_x_section(target_ega2);
0303           //
0304
0305           //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
0306           dR  = q2_cor*csgA1;
0307           dR  = dR + 4.*q2_cor*csgA12;
0308           dR  = dR + q2_cor*csgA2;
0309           dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
0310           //
0311           dR2  = full_range_1*csgA1;
0312           dR2  = dR2 + 4.*full_range_12*csgA12;
0313           dR2  = dR2 + full_range_2*csgA2;
0314           dR2  = dR2*(dY/6.)*breitWigner(W,bwnorm)*dW;
0315           //
0316           int_r = int_r+dR;
0317           int_r2 = int_r2 +dR2;
0318         }
0319       //
0320       }
0321       if( iQ2 ==0 )
0322         cout<<"Full range "<<int_r2*10000000<<endl;
0323       cout<<(q2Edge[iQ2+1]+q2Edge[iQ2])/2.+W*W<<" ,  "<<10000000.*int_r<<endl;
0324     }
0325
0326 }
0327 #endif
0328 void e_wideResonanceCrossSection::printCrossSection(const string name, const double x_section)
0329 {
0330   if (0.01*x_section > 1.){
0331     cout<< name.c_str() <<0.01*x_section<<" barn."<<endl;
0332   } else if (10.*x_section > 1.){
0333     cout<< name.c_str() <<10.*x_section<<" mb."<<endl;
0334   } else if (10000.*x_section > 1.){
0335     cout<< name.c_str() <<10000.*x_section<<" microb."<<endl;
0336   } else if (10000000.*x_section > 1.){
0337     cout<< name.c_str() <<10000000.*x_section<<" nanob."<<endl;
0338   } else if (1.E10*x_section > 1.){
0339     cout<< name.c_str() <<1.E10*x_section<<" picob."<<endl;
0340   } else {
0341     cout<< name.c_str() <<1.E13*x_section<<" femtob."<<endl;
0342   }
0343   //cout<<endl;
0344 }