Page MenuHomeHEPForge
Files F9501494

No OneTemporary
Actions

View Options

Index: branches/inclusive-diffraction/src/InclusiveDiffractiveCrossSections.cpp
===================================================================
--- branches/inclusive-diffraction/src/InclusiveDiffractiveCrossSections.cpp (revision 0)
+++ branches/inclusive-diffraction/src/InclusiveDiffractiveCrossSections.cpp (revision 137)
@@ -0,0 +1,1018 @@
+//==============================================================================
+// InclusiveDiffractiveCrossSecitions.cpp
+//
+// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+//
+// This file is part of Sartre version: 2.0
+//
+// This program 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.
+// This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+//
+// Author: Tobias Toll
+// Last update:
+// $Date: 2013-05-29 16:26:19 -0400 (Wed, 29 May 2013) $
+// $Author: ttoll@bnl.gov $
+//==============================================================================
+//How to get exclusive cross-sections:
+// init()
+// calculate_dsigmadbetadz_QQ(double, double, double, double);
+// calculate_dsigmadbetadz_QQG(double, double, double, double);
+// Add all elements in
+// double* dsigmadbetaT_QQ(); and double* dsigmadbetaT_QQG(); and
+// and in double* dsigmadbetaL_QQ();
+
+#include "Math/SpecFunc.h"
+
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+#include "Constants.h"
+#include "Nucleus.h"
+#include "DipoleModel.h"
+#include "AlphaStrong.h"
+#include "Math/IntegratorMultiDim.h"
+#include "Math/Functor.h"
+#include "TMath.h"
+#include "WaveOverlap.h"
+#include "Kinematics.h"
+#include "TableGeneratorSettings.h"
+#include "Enumerations.h"
+#include "TF1.h"
+#include "TH1F.h"
+#include "cuba.h"
+#include "InclusiveDiffractiveCrossSections.h"
+#include "DglapEvolution.h"
+#include "Math/WrappedTF1.h"
+#include "Math/GaussIntegrator.h"
+
+#define PR(x) cout << #x << " = " << (x) << endl;
+
+using namespace std;
+
+InclusiveDiffractiveCrossSections::InclusiveDiffractiveCrossSections(){}
+
+void InclusiveDiffractiveCrossSections::setRelativePrecisionOfIntegration(double val){
+ mRelativePrecisionOfIntegration=val;
+}
+
+void InclusiveDiffractiveCrossSections::init()
+{
+ cout<<"init InclusiveDiffractiveCrossSections"<<endl;
+ mDipoleModel = 0;
+ mQuarkIndex = 0;
+ mBetaMin = 0;
+
+ setRelativePrecisionOfIntegration(1e-3);
+ TableGeneratorSettings* settings = TableGeneratorSettings::instance();
+
+ DipoleModelType model=settings->dipoleModel();
+ if (model==bSat) {
+ mDipoleModel = new DipoleModel_bSat;
+ }
+ else if(model==bNonSat){
+ mDipoleModel = new DipoleModel_bNonSat;
+ }
+ else if (model==bCGC) {
+ mDipoleModel = new DipoleModel_bCGC;
+ }
+ else {
+ cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
+ exit(1);
+ }
+ if (model==bSat || model==bNonSat){
+ cout << "InclusiveDiffractiveCrossSections::init(): Initializing the DGLAP engine:" << endl;
+ //Upper scale in the DGLAP evoultion. Should be high enough to cover LHeC:
+ DglapEvolution::instance().setS(1e10); //this should do it...
+ DglapEvolution::instance().G(0.01, 4.);
+ }
+ mA=settings->A();
+ mIsInitialized=true;
+
+}
+InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){}
+/*
+InclusiveDiffractiveCrossSections& InclusiveDiffractiveCrossSections::operator=(const InclusiveDiffractiveCrossSections& cobj)
+{
+ if (this != &cobj) {
+ Integrals::operator=(cobj);
+ copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
+ }
+ return *this;
+}
+
+InclusiveDiffractiveCrossSections::InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections& cobj) : Integrals(cobj)
+{
+ copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
+}
+*/
+void InclusiveDiffractiveCrossSections::setZ(double val) { mZ = val; }
+void InclusiveDiffractiveCrossSections::setQuarkIndex(int val) { mQuarkIndex = val; }
+
+double InclusiveDiffractiveCrossSections::zlower(){
+ double mf=quarkMass[quarkIndex()];
+ double Q2=kinematicPoint[1];
+ double beta=kinematicPoint[4];
+ double MX2=Q2*(1-beta)/beta;
+ double root2=1.-4*mf*mf/MX2;
+ double root=0;
+ if(root2>0)
+ root=sqrt(root2);
+ return 0.5*(1-root);
+}
+bool InclusiveDiffractiveCrossSections::setKinematicPoint(double t, double Q2, double W2, double beta)
+{
+ bool result=true;
+ kinematicPoint[0]=t;
+ kinematicPoint[1]=Q2;
+ kinematicPoint[2]=W2;
+ double x=Kinematics::x(Q2, W2);
+ double xpom=x/beta;
+ if (xpom<0 || xpom>1)
+ result = false;
+ kinematicPoint[3]=xpom;
+ kinematicPoint[4]=beta;
+ return result;
+}
+
+void InclusiveDiffractiveCrossSections::setNumberOfFlavours(unsigned int val){ mNumberOfFlavours=val; }
+
+/////////////Cross-Sections//////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQ(double Q2, double W2, double beta){
+ if(beta<mBetaMin) return;
+ if(!mIsInitialized){
+ cout<<"InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQ Warning: InclusiveDiffractiveCrossSections not initialised!"<<endl;
+ exit(0);
+ }
+ for(int i=0; i<6; i++){
+ mInclusiveT[i]=0;
+ mInclusiveL[i]=0;
+ }
+ // Set Kinematics and check for validity simultaneously
+ if(!setKinematicPoint(0, Q2, W2, beta)){
+ cout<<"InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQ Warning: invalid kinematic point"<<endl;
+ return;
+ }
+ double xpom=kinematicPoint[3];
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom);
+
+ double resultT=0;
+ double resultL=0;
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ double ef=quarkCharge[i];
+ //
+ // The integration will be identical for equal quark masses
+ // This is checked for in order not to redo the same integral
+ // If integrate=false, the previous results will be reused
+ //
+ bool integrate=true;
+ if(i>0){
+ if(quarkMass[i-1]==quarkMass[i])
+ integrate=false;
+ }
+ if(integrate){
+ double zlow=zlower();
+ double zhigh=0.5;
+ TF1 IntegralT("IntegralT", this, &InclusiveDiffractiveCrossSections::uiIntegralT, zlow, zhigh, 3);
+ IntegralT.SetParameter(0, beta);
+ IntegralT.SetParameter(1, xpom);
+ IntegralT.SetParameter(2, Q2);
+ ROOT::Math::WrappedTF1 wfT(IntegralT);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ resultT=ig.Integral(zlow, zhigh); //GeV^-4
+
+ TF1 IntegralL("IntegralL", this, &InclusiveDiffractiveCrossSections::uiIntegralL, zlow, zhigh, 3);
+ IntegralL.SetParameter(0, beta);
+ IntegralL.SetParameter(1, xpom);
+ IntegralL.SetParameter(2, Q2);
+ ROOT::Math::WrappedTF1 wfL(IntegralL);
+ ROOT::Math::GaussIntegrator ig2;
+ ig2.SetFunction(wfL);
+ ig2.SetRelTolerance(mRelativePrecisionOfIntegration);
+ resultL=ig2.Integral(zlow, zhigh); //GeV^-6
+ }
+ //Store the results:
+ mInclusiveT[i]=resultT*ef*ef*alpha_em*Nc*Q2/4/M_PI/beta/beta;
+ mInclusiveL[i]=resultL*ef*ef*alpha_em*Nc*Q2*Q2/M_PI/beta/beta;
+ }
+}
+
+void InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQ(double Q2, double W2, double beta, double z){
+ if(beta<mBetaMin) return;
+ if(!mIsInitialized){
+ cout<<"InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQ Warning: InclusiveDiffractiveCrossSections not initialised!"<<endl;
+ exit(0);
+ }
+ for(int i=0; i<6; i++){
+ mInclusiveT[i]=0;
+ mInclusiveL[i]=0;
+ }
+ // Set Kinematics and check for validity simultaneously
+ if(!setKinematicPoint(0, Q2, W2, beta)){
+ cout<<"InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQ Warning: invalid kinematic point"<<endl;
+ return;
+ }
+ double xpom=kinematicPoint[3];
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(kinematicPoint[3]);
+
+ double resultT=0;
+ double resultL=0;
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ double ef=quarkCharge[i];
+ //
+ // The integration will be identical for equal quark masses
+ // This is checked for in order not to redo the same integral
+ // If integrate=false, the previous results will be reused
+ //
+ bool integrate=true;
+ if(i>0){
+ if(quarkMass[i-1]==quarkMass[i])
+ integrate=false;
+ }
+ if(integrate){
+ if(z<zlower()){
+ resultT=0;
+ resultL=0;
+ continue;
+ }
+ double par[3]={beta, xpom, Q2};
+ resultT=uiIntegralT(&z, par);
+ resultL=uiIntegralL(&z, par);
+ }
+ //Store the results:
+ mInclusiveT[i]=resultT*ef*ef*alpha_em*Nc*Q2/4/M_PI/beta/beta;
+ mInclusiveL[i]=resultL*ef*ef*alpha_em*Nc*Q2*Q2/M_PI/beta/beta;
+ }
+}
+
+double InclusiveDiffractiveCrossSections::uiIntegralL(double* x, double* par){
+ double beta=par[0];
+ double xpom=par[1];
+ double Q2=par[2];
+ double z=*x;
+ double blow=0;
+ double bhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 phi0("Phi0", this, &InclusiveDiffractiveCrossSections::dPhi0db, blow, bhigh, 4);
+ phi0.SetParameter(0, beta);
+ phi0.SetParameter(1, xpom);
+ phi0.SetParameter(2, z);
+ phi0.SetParameter(3, Q2);
+ ROOT::Math::WrappedTF1 wfT(phi0);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ double Phi0=ig.Integral(blow, bhigh); //fm^6
+ Phi0 = Phi0 / hbarc2/hbarc2/hbarc2; //GeV^-6
+ return z*z*z*(1-z)*(1-z)*(1-z)*Phi0; //GeV-6
+}
+
+double InclusiveDiffractiveCrossSections::uiIntegralT(double* x, double* par){
+ double beta=par[0];
+ double xpom=par[1];
+ double Q2=par[2];
+ double z=*x;
+ double blow=0;
+ double bhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 phi0("Phi0", this, &InclusiveDiffractiveCrossSections::dPhi0db, blow, bhigh, 4);
+ phi0.SetParameter(0, beta);
+ phi0.SetParameter(1, xpom);
+ phi0.SetParameter(2, z);
+ phi0.SetParameter(3, Q2);
+ ROOT::Math::WrappedTF1 wf0(phi0);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wf0);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ double Phi0=ig.Integral(blow, bhigh); //fm^6
+ Phi0 = Phi0/hbarc2/hbarc2/hbarc2; //GeV^-6
+
+ TF1 phi1("Phi1", this, &InclusiveDiffractiveCrossSections::dPhi1db, blow, bhigh, 4);
+ phi1.SetParameter(0, beta);
+ phi1.SetParameter(1, xpom);
+ phi1.SetParameter(2, z);
+ phi1.SetParameter(3, Q2);
+ ROOT::Math::WrappedTF1 wf1(phi1);
+ ROOT::Math::GaussIntegrator ig2;
+ ig2.SetFunction(wf1);
+ ig2.SetRelTolerance(mRelativePrecisionOfIntegration);
+ double Phi1=ig2.Integral(blow, bhigh); //fm^6
+ Phi1 = Phi1/hbarc2/hbarc2/hbarc2; //GeV^-6
+
+ double mf=quarkMass[mQuarkIndex]; //GeV
+ double epsilon2=z*(1-z)*Q2+mf*mf; //GeV^2
+
+ return z*(1-z)*(epsilon2*(z*z+(1-z)*(1-z))*Phi1+mf*mf*Phi0); //GeV^-4
+}
+
+double InclusiveDiffractiveCrossSections::dPhi0db(double* x, double* par){
+ double beta=par[0];
+ double xpom=par[1];
+ double z=par[2];
+ double Q2=par[3];
+ double b=*x;
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 phi0("phi0", this, &InclusiveDiffractiveCrossSections::uiPhi0, rlow, rhigh, 5);
+
+ phi0.SetParameter(0, beta);
+ phi0.SetParameter(1, xpom);
+ phi0.SetParameter(2, z);
+ phi0.SetParameter(3, Q2);
+ phi0.SetParameter(4, b);
+ ROOT::Math::WrappedTF1 wfT(phi0);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ double integral=ig.Integral(rlow, rhigh); //fm^2
+ return integral*integral*2*M_PI*b; //fm^5
+}
+
+double InclusiveDiffractiveCrossSections::dPhi1db(double* x, double* par){
+ double beta=par[0];
+ double xpom=par[1];
+ double z=par[2];
+ double Q2=par[3];
+ double b=*x;
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 phi1("phi1", this, &InclusiveDiffractiveCrossSections::uiPhi1, rlow, rhigh, 5);
+
+ phi1.SetParameter(0, beta);
+ phi1.SetParameter(1, xpom);
+ phi1.SetParameter(2, z);
+ phi1.SetParameter(3, Q2);
+ phi1.SetParameter(4, b);
+ ROOT::Math::WrappedTF1 wfT(phi1);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ double integral=ig.Integral(rlow, rhigh); //fm^2
+ return integral*integral*2*M_PI*b; //fm^5
+}
+
+double InclusiveDiffractiveCrossSections::uiPhi0(double* x, double* par){
+ double r=*x;
+ double beta=par[0];
+ double xpom=par[1];
+ double z=par[2];
+ double Q2=par[3];
+ double b=par[4];
+ double mf=quarkMass[mQuarkIndex];
+ double MX2=Q2*(1-beta)/beta;
+ if(1-4*mf*mf/MX2 <0) return 0; //This should never happen but it does. Edge effect in integration limits?
+ double kappa2=z*(1-z)*MX2-mf*mf;
+ double kappa=sqrt(kappa2);
+ double epsilon2=z*(1-z)*Q2+mf*mf;
+ double epsilon=sqrt(epsilon2);
+ double besselK0=TMath::BesselK0(epsilon*r/hbarc);
+ double besselJ0=TMath::BesselJ0(kappa*r/hbarc);
+ double dsigmadb2=1;
+ if(mA==1)
+ dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
+ else
+ dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
+ return r*besselK0*besselJ0*dsigmadb2; //fm
+}
+double InclusiveDiffractiveCrossSections::uiPhi1(double* x, double* par){
+ double r=*x;
+ double beta=par[0];
+ double xpom=par[1];
+ double z=par[2];
+ double Q2=par[3];
+ double b=par[4];
+ double mf=quarkMass[mQuarkIndex];
+ double MX2=Q2*(1-beta)/beta;
+ if(1-4*mf*mf/MX2 <0) return 0; //This should never happen but it does. Edge effect in integration limits?
+ double kappa2=z*(1-z)*MX2-mf*mf;
+ double kappa=sqrt(kappa2);
+ double epsilon2=z*(1-z)*Q2+mf*mf;
+ double epsilon=sqrt(epsilon2);
+ double besselK1=TMath::BesselK1(epsilon*r/hbarc);
+ double besselJ1=TMath::BesselJ1(kappa*r/hbarc);
+ double dsigmadb2=0;
+ if(mA==1)
+ dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
+ else
+ dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
+ return r*besselK1*besselJ1*dsigmadb2; //fm
+}
+///////////////////////////////////////////////////////////////////////////////////////////
+////QQG////////////////////////////////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::putAllQQGtoZero(){
+ for(int i=0; i<6; i++){
+ mInclusiveQQG[i]=0;
+ mInclusiveQQG_GBW[i]=0;
+ mInclusiveQQG_MS[i]=0;
+ mInclusiveQQG_GBWbeta0[i]=0;
+ }
+}
+void InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQG(double Q2, double W2, double beta){
+ if(!mIsInitialized){
+ cout<<"InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQG Warning: InclusiveDiffractiveCrossSections not initialised!"<<endl;
+ exit(0);
+ }
+ putAllQQGtoZero();
+ // dsigma/dbeta = 4*M_PI*M_PI*alpha_em/(Q2*beta) * xpom * F
+ double xF_to_dsigmadbeta=4*M_PI*M_PI*alpha_em/Q2/beta; //GeV^-2
+ if(!setKinematicPoint(0., Q2, W2, beta)){ //this is calculated in the limit beta->0
+ cout<<"Warning, InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQG invalid kinematic point!"<<endl;
+ return;
+ }
+ double xpom=kinematicPoint[3];
+ FTQQG_MS(Q2, xpom);
+ FTQQG_GBW(Q2, beta);
+ FTQQG_GBWbeta0(Q2);
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ mInclusiveQQG[i]=mInclusiveQQG_GBW[i]*mInclusiveQQG_MS[i]/mInclusiveQQG_GBWbeta0[i];
+ mInclusiveQQG[i]*=xF_to_dsigmadbeta;
+ }
+}
+void InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQG(double Q2, double W2, double beta, double z){
+ if(!mIsInitialized){
+ cout<<"InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQ Warning: InclusiveDiffractiveCrossSections not initialised!"<<endl;
+ exit(0);
+ }
+ putAllQQGtoZero();
+ // dsigma/dbeta = 4*M_PI*M_PI*alpha_em/(Q2*beta) * xpom * F
+ double xF_to_dsigmadbeta=4*M_PI*M_PI*alpha_em/Q2/beta; //GeV^-2
+ if(!setKinematicPoint(0., Q2, W2, beta)){ //this is calculated in the limit beta->0
+ cout<<"Warning, InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQG invalid kinematic point!"<<endl;
+ return;
+ }
+ double xpom=kinematicPoint[3];
+ FTQQGdz_MS(Q2, xpom, z);
+ FTQQGdz_GBW(Q2, beta, z);
+ FTQQG_GBWbeta0(Q2);
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ mInclusiveQQG[i]=mInclusiveQQG_GBW[i]*mInclusiveQQG_MS[i]/mInclusiveQQG_GBWbeta0[i];
+ mInclusiveQQG[i]*=xF_to_dsigmadbeta;
+ }
+
+}
+///////////////////////////////////////////////////////////////////////////////////////////
+////QQG_GBW/dz/////////////////////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::FTQQGdz_GBW(double Q2, double beta, double z){
+ mBeta=beta;
+ double xpom=kinematicPoint[3];
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom);
+ double result=0;
+ double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ //
+ // The integration will be identical for equal quark masses
+ // This is checked for in order not to redo the same integral
+ // If integrate=false, the previous results will be reused
+ //
+ bool integrate=true;
+ if(i>0){
+ if(quarkMass[i-1]==quarkMass[i])
+ integrate=false;
+ }
+ if(integrate){
+ //beta is the lower limit of z in this treatment of QQG:
+ if(z<beta) {
+ cout<<"Warning in InclusiveDiffractiveCrossSections::FTQQGdz_GBW: z<beta, returning"<<endl;
+ return;
+ }
+ mZ=z;
+ double a[2]={0., 0.}; //b, k2
+ double b[2]={rbhigh, Q2};
+ ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiQQG_GBW, 2);
+ ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
+ ig.SetFunction(wf);
+ ig.SetAbsTolerance(0.);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ result=ig.Integral(a, b)/hbarc; //GeV^-1*GeV^-1*GeV^2=GeV^0
+ }
+ //Store the result:
+ double alpha_s=0.14;
+ double ef=quarkCharge[quarkIndex()];
+ mInclusiveQQG_GBW[i]=result*alpha_s*beta/(8*pow(M_PI, 4))*ef*ef; //GeV^0
+ }
+}
+double InclusiveDiffractiveCrossSections::uiQQGdz_GBW(const double* arg){
+ double b=arg[0];
+ double k2=arg[1];
+
+ double z=mZ;
+ double Q2=kinematicPoint[1];
+ double beta=mBeta;
+
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 frIntegral("frIntegral", this, &InclusiveDiffractiveCrossSections::uiR_GBW, rlow, rhigh, 3);
+
+ frIntegral.SetParameter(0, z);
+ frIntegral.SetParameter(1, b);
+ frIntegral.SetParameter(2, k2);
+
+ ROOT::Math::WrappedTF1 wfT(frIntegral);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ double rIntegral=ig.Integral(rlow, rhigh)/hbarc; //GeV^-2
+ double result=2*M_PI*b/hbarc* //GeV^-1
+ k2*k2*log(Q2/k2)*((1-beta/z)*(1-beta/z)+beta*beta/(z*z))* //GeV^4
+ rIntegral*rIntegral; //GeV^-4
+ return result; //GeV^-1
+}
+
+////QQG_GBW////////////////////////////////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::FTQQG_GBW(double Q2, double beta){
+ mBeta=beta;
+ double xpom=kinematicPoint[3];
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom);
+ double result=0;
+ double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ //
+ // The integration will be identical for equal quark masses
+ // This is checked for in order not to redo the same integral
+ // If integrate=false, the previous results will be reused
+ //
+ bool integrate=true;
+ if(i>0){
+ if(quarkMass[i-1]==quarkMass[i])
+ integrate=false;
+ }
+ if(integrate){
+ double a[3]={0., 0., beta}; //b, k2, z
+ double b[3]={rbhigh, Q2, 1.};
+ ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiQQG_GBW, 3);
+ ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
+ ig.SetFunction(wf);
+ ig.SetAbsTolerance(0.);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ result=ig.Integral(a, b)/hbarc; //GeV^-1*GeV^-1*GeV^2=GeV^0
+ }
+ //Store the result:
+ double alpha_s=0.14;
+ double ef=quarkCharge[quarkIndex()];
+ mInclusiveQQG_GBW[i]=result*alpha_s*beta/(8*pow(M_PI, 4))*ef*ef; //GeV^0
+ }
+}
+
+double InclusiveDiffractiveCrossSections::uiQQG_GBW(const double* arg){
+ double b=arg[0];
+ double k2=arg[1];
+ double z=arg[2];
+
+ double Q2=kinematicPoint[1];
+ double beta=mBeta;
+
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 frIntegral("frIntegral", this, &InclusiveDiffractiveCrossSections::uiR_GBW, rlow, rhigh, 3);
+
+ frIntegral.SetParameter(0, z);
+ frIntegral.SetParameter(1, b);
+ frIntegral.SetParameter(2, k2);
+
+ ROOT::Math::WrappedTF1 wfT(frIntegral);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ double rIntegral=ig.Integral(rlow, rhigh)/hbarc; //GeV^-2
+ double result=2*M_PI*b/hbarc* //GeV^-1
+ k2*k2*log(Q2/k2)*((1-beta/z)*(1-beta/z)+beta*beta/(z*z))* //GeV^4
+ rIntegral*rIntegral; //GeV^-4
+ return result; //GeV^-1
+}
+
+double InclusiveDiffractiveCrossSections::uiR_GBW(double* x, double* par){
+ double r=*x;
+ double z=par[0];
+ double b=par[1];
+ double k2=par[2];
+
+ double xpom=kinematicPoint[3];
+
+ double k=sqrt(k2);
+
+ double besselK2=TMath::BesselK(2, sqrt(z)*k*r/hbarc);
+ double besselJ2=ROOT::Math::cyl_bessel_j(2, sqrt(1-z)*k*r/hbarc);
+
+ double result=r/hbarc*dsigmadb2tilde(r, b, xpom)*besselK2*besselJ2;
+
+ return result; //GeV^-1
+}
+
+double InclusiveDiffractiveCrossSections::dsigmadb2tilde(double r, double b, double xpom){
+ double dsigmadb2=0;
+ if(mA==1)
+ dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
+ else
+ dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
+
+ double termToSquare=1-0.5*dsigmadb2;
+ return 2*(1-termToSquare*termToSquare);
+}
+////QQG_GBW beta=0///////////////////////////////////////////////////////////////////////
+// It turns out that it not as simple as just putting beta=0 in the argument...
+void InclusiveDiffractiveCrossSections::FTQQG_GBWbeta0(double Q2){
+ double xpom=kinematicPoint[3];
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom);
+ double result=0;
+ double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ //
+ // The integration will be identical for equal quark masses
+ // This is checked for in order not to redo the same integral
+ // If integrate=false, the previous results will be reused
+ //
+ bool integrate=true;
+ if(i>0){
+ if(quarkMass[i-1]==quarkMass[i])
+ integrate=false;
+ }
+ if(integrate){
+ double a[3]={0., 0.}; //b, k2
+ double b[3]={rbhigh, Q2};
+ ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiQQG_GBWbeta0, 2);
+ ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
+ ig.SetFunction(wf);
+ ig.SetAbsTolerance(0.);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ result=ig.Integral(a, b)/hbarc; //GeV^-1*GeV^-1*GeV^2=GeV^0
+ }
+ //Store the result:
+ double alpha_s=0.14;
+ double ef=quarkCharge[quarkIndex()];
+ //This is both from Cyrille and confirmed by letting beta approach 0 in the other calculation
+ double mysteryFactor=2./3;
+ mInclusiveQQG_GBWbeta0[i]=result*alpha_s/(2*pow(M_PI, 4))*ef*ef*mysteryFactor; //GeV^0
+ }
+}
+
+double InclusiveDiffractiveCrossSections::uiQQG_GBWbeta0(const double* arg){
+ double b=arg[0];
+ double k2=arg[1];
+
+ double Q2=kinematicPoint[1];
+
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ TF1 frIntegral("frIntegral", this, &InclusiveDiffractiveCrossSections::uiR_GBWbeta0, rlow, rhigh, 2);
+ frIntegral.SetParameter(0, b);
+ frIntegral.SetParameter(1, k2);
+
+ ROOT::Math::WrappedTF1 wfT(frIntegral);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ double rIntegral=ig.Integral(rlow, rhigh)/hbarc; //GeV^0
+ double result=2*M_PI*b/hbarc* //GeV^-1
+ log(Q2/k2)*//GeV^0
+ rIntegral*rIntegral; //GeV^0
+ return result; //GeV^-1
+}
+double InclusiveDiffractiveCrossSections::uiR_GBWbeta0(double* x, double* par){
+ double r=*x;
+ double b=par[0];
+ double k2=par[1];
+
+ double xpom=kinematicPoint[3];
+ double k=sqrt(k2);
+
+ double besselJ2=0;
+ if(r!=0)
+ besselJ2=ROOT::Math::cyl_bessel_j(2, k*r/hbarc)/(r/hbarc);
+
+ double result=dsigmadb2tilde(r, b, xpom)*besselJ2;
+ return result; //GeV^1
+}
+
+///////////////////////////////////////////////////////////////////////////////////////////
+////QQGdz_MS////////////////////////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::FTQQGdz_MS(double Q2, double xpom, double z){
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom);
+
+ double result=0;
+ double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ mZ=z;
+ double a[4]={1e-3, 0, 1e-3, 0.}; //r, z, b, rp, phi
+ double b[4]={rbhigh, rbhigh, rbhigh, 2*M_PI};
+ ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiQQG_MS, 4);
+ ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
+ ig.SetFunction(wf);
+ ig.SetAbsTolerance(0.);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ result=ig.Integral(a, b)/hbarc/hbarc/hbarc; //GeV^3*GeV^-3=GeV^0
+
+ //Store the result:
+ double alpha_s=0.14;
+ double Cf=4/3;
+ mInclusiveQQG_MS[i]=result*Cf*alpha_s*Q2/(4*pow(M_PI, 4)*alpha_em);
+ }
+}
+
+double InclusiveDiffractiveCrossSections::uiQQGdz_MS(const double* arg){
+ double r=arg[0];
+ double b=arg[2];
+ double rp=arg[3];
+ double phi=arg[4];
+
+ double z=mZ;
+ double xpom=kinematicPoint[3];
+ double Q2=kinematicPoint[1];
+
+ //Calculate A
+ double N_of_r=0;
+ if(mA==1)
+ N_of_r=mDipoleModel->dsigmadb2ep(r, b, xpom)/2.;
+ else
+ N_of_r=mDipoleModel->coherentDsigmadb2(r, b, xpom)/2.;
+
+ double N_of_rp=0;
+ if(mA==1)
+ N_of_rp=mDipoleModel->dsigmadb2ep(rp, b, xpom)/2.;
+ else
+ N_of_rp=mDipoleModel->coherentDsigmadb2(rp, b, xpom)/2.;
+
+ double rMinusrp=sqrt(r*r+rp*rp-2*r*rp*cos(phi));
+ double N_of_rMinusrp=0;
+ if(mA==1)
+ N_of_rMinusrp=mDipoleModel->dsigmadb2ep(rMinusrp, b, xpom)/2.;
+ else
+ N_of_rMinusrp=mDipoleModel->coherentDsigmadb2(rMinusrp, b, xpom)/2.;
+
+ double rprefactor=rp * r*r/(rp*rp*rMinusrp*rMinusrp)*hbarc; //GeV
+ double Nfactor=N_of_rp+N_of_rMinusrp-N_of_r-N_of_rp*N_of_rMinusrp; //GeV^0
+
+ double A=rprefactor*Nfactor*Nfactor; //GeV
+ //WaveOverlap (different prefactor in MS and KMW):
+ double waveOverlap=waveOverlapT(r, z, Q2)/(2.*M_PI);///(4*M_PI); //GeV^2
+
+ double result= waveOverlap*A*2*M_PI*b/hbarc*2*M_PI*r/hbarc; //GeV^4*GeV*GeV^-2=GeV^3
+ return result; //GeV^3
+}
+///////////////////////////////////////////////////////////////////////////////////////////
+////QQG_MS////////////////////////////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::FTQQG_MS(double Q2, double xpom){
+ // Make the coherent lookup-table
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom);
+
+ double result=0;
+ double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ double a[5]={1e-3, 0, 0, 1e-3, 0.}; //r, z, b, rp, phi
+ double b[5]={rbhigh, 1, rbhigh, rbhigh, 2*M_PI};
+ ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiQQG_MS, 5);
+ ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
+ ig.SetFunction(wf);
+ ig.SetAbsTolerance(0.);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ result=ig.Integral(a, b)/hbarc/hbarc/hbarc; //GeV^3*GeV^-3=GeV^0
+
+ //Store the result:
+ double alpha_s=0.14;
+ double Cf=4/3;
+ mInclusiveQQG_MS[i]=result*Cf*alpha_s*Q2/(4*pow(M_PI, 4)*alpha_em);
+ }
+}
+
+double InclusiveDiffractiveCrossSections::uiQQG_MS(const double* arg){
+ double r=arg[0];
+ double z=arg[1];
+ double b=arg[2];
+ double rp=arg[3];
+ double phi=arg[4];
+
+ double xpom=kinematicPoint[3];
+ double Q2=kinematicPoint[1];
+
+ //Calculate A
+ double N_of_r=0;
+ if(mA==1)
+ N_of_r=mDipoleModel->dsigmadb2ep(r, b, xpom)/2.;
+ else
+ N_of_r=mDipoleModel->coherentDsigmadb2(r, b, xpom)/2.;
+
+ double N_of_rp=0;
+ if(mA==1)
+ N_of_rp=mDipoleModel->dsigmadb2ep(rp, b, xpom)/2.;
+ else
+ N_of_rp=mDipoleModel->coherentDsigmadb2(rp, b, xpom)/2.;
+
+ double rMinusrp=sqrt(r*r+rp*rp-2*r*rp*cos(phi));
+ double N_of_rMinusrp=0;
+ if(mA==1)
+ N_of_rMinusrp=mDipoleModel->dsigmadb2ep(rMinusrp, b, xpom)/2.;
+ else
+ N_of_rMinusrp=mDipoleModel->coherentDsigmadb2(rMinusrp, b, xpom)/2.;
+
+ double rprefactor=rp * r*r/(rp*rp*rMinusrp*rMinusrp)*hbarc; //GeV
+ double Nfactor=N_of_rp+N_of_rMinusrp-N_of_r-N_of_rp*N_of_rMinusrp; //GeV^0
+
+ double A=rprefactor*Nfactor*Nfactor; //GeV
+ //WaveOverlap (different prefactor in MS and KMW):
+ double waveOverlap=waveOverlapT(r, z, Q2)/(2.*M_PI);///(4*M_PI); //GeV^2
+
+ double result= waveOverlap*A*2*M_PI*b/hbarc*2*M_PI*r/hbarc; //GeV^4*GeV*GeV^-2=GeV^3
+ return result; //GeV^3
+}
+///////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////////////////
+/////////////Total Cross-Sections//////////////////////////////////////////////////////////
+void InclusiveDiffractiveCrossSections::calculateTotal(double Q2, double W2){
+ if(!mIsInitialized){
+ cout<<"InclusiveDiffractiveCrossSections::dsigmadbeta Warning: InclusiveDiffractiveCrossSections not initialised!"<<endl;
+ exit(0);
+ }
+ for(int i=0; i<6; i++){
+ mInclusiveT[i]=0;
+ mInclusiveL[i]=0;
+ }
+ // Make the coherent lookup-table
+ double x=Kinematics::x(Q2, W2);
+ if(mA>1) dipoleModel()->createSigma_ep_LookupTable(x);
+
+ double resultT=0;
+ double resultL=0;
+ mTotalCST=0;
+ mTotalCSL=0;
+ for(unsigned int i=0; i<mNumberOfFlavours; i++){
+ setQuarkIndex(i);
+ double zlow=0;
+ double zhigh=1;
+ TF1 IntegralT("IntegralT", this, &InclusiveDiffractiveCrossSections::uiDsigmadzT, zlow, zhigh, 2);
+ IntegralT.SetParameter(0, Q2);
+ IntegralT.SetParameter(1, W2);
+ ROOT::Math::WrappedTF1 wfT(IntegralT);
+ ROOT::Math::GaussIntegrator igT;
+ igT.SetFunction(wfT);
+ igT.SetRelTolerance(mRelativePrecisionOfIntegration);
+ resultT=igT.Integral(zlow, zhigh); //GeV^2*fm^4
+
+ TF1 IntegralL("IntegralL", this, &InclusiveDiffractiveCrossSections::uiDsigmadzL, zlow, zhigh, 2);
+ IntegralL.SetParameter(0, Q2);
+ IntegralL.SetParameter(1, W2);
+ ROOT::Math::WrappedTF1 wfL(IntegralL);
+ ROOT::Math::GaussIntegrator igL;
+ igL.SetFunction(wfL);
+ igL.SetRelTolerance(mRelativePrecisionOfIntegration);
+ resultL=igL.Integral(zlow, zhigh); //GeV^2*fm^4
+ //Store the results:
+ mTotalCST+=resultT/hbarc2/hbarc2; //GeV^-2
+ mTotalCSL+=resultL/hbarc2/hbarc2; //GeV^-2
+ } //flavour
+}
+
+double InclusiveDiffractiveCrossSections::uiDsigmadzT(double* x, double* par){
+ double z=*x;
+ double Q2=par[0];
+ double W2=par[1];
+
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+ TF1 sigma("sigma", this, &InclusiveDiffractiveCrossSections::uiDsigmadrT, rlow, rhigh, 3);
+ sigma.SetParameter(0, z);
+ sigma.SetParameter(1, Q2);
+ sigma.SetParameter(2, W2);
+ ROOT::Math::WrappedTF1 wfT(sigma);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ double dipoleCS=ig.Integral(rlow, rhigh); //GeV^2*fm^4
+
+ return dipoleCS;
+}
+double InclusiveDiffractiveCrossSections::uiDsigmadzL(double* x, double* par){
+ double z=*x;
+ double Q2=par[0];
+ double W2=par[1];
+
+ double rlow=1e-3;
+ double rhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+ TF1 sigma("sigma", this, &InclusiveDiffractiveCrossSections::uiDsigmadrL, rlow, rhigh, 3);
+ sigma.SetParameter(0, z);
+ sigma.SetParameter(1, Q2);
+ sigma.SetParameter(2, W2);
+ ROOT::Math::WrappedTF1 wfT(sigma);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration/10.);
+ double dipoleCS=ig.Integral(rlow, rhigh); //GeV^2*fm^4
+
+ return dipoleCS;
+}
+double InclusiveDiffractiveCrossSections::uiDsigmadrT(double* x, double* par){
+ double r=*x;
+ double z=par[0];
+ double Q2=par[1];
+ double W2=par[2];
+
+ double waveOverlap=waveOverlapT(r, z, Q2);
+ double xbj=Kinematics::x(Q2, W2);
+ double blow=0;
+ double bhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+ TF1 dp("dp", this, &InclusiveDiffractiveCrossSections::uiDipoleSigma, blow, bhigh, 2);
+ dp.SetParameter(0, r);
+ dp.SetParameter(1, xbj);
+ ROOT::Math::WrappedTF1 wfT(dp);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ double dipoleCS=ig.Integral(blow, bhigh); //fm^2
+
+ return waveOverlap*dipoleCS/2.*r; //GeV^2*fm^3
+}
+double InclusiveDiffractiveCrossSections::uiDsigmadrL(double* x, double* par){
+ double r=*x;
+ double z=par[0];
+ double Q2=par[1];
+ double W2=par[2];
+
+ double waveOverlap=waveOverlapL(r, z, Q2);
+ double xbj=Kinematics::x(Q2, W2);
+ double blow=0;
+ double bhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius();
+ TF1 dp("dp", this, &InclusiveDiffractiveCrossSections::uiDipoleSigma, blow, bhigh, 2);
+ dp.SetParameter(0, r);
+ dp.SetParameter(1, xbj);
+ ROOT::Math::WrappedTF1 wfT(dp);
+ ROOT::Math::GaussIntegrator ig;
+ ig.SetFunction(wfT);
+ ig.SetRelTolerance(mRelativePrecisionOfIntegration);
+ double dipoleCS=ig.Integral(blow, bhigh); //fm^2
+
+ return waveOverlap*dipoleCS/2.*r; //GeV^2*fm^3
+}
+
+double InclusiveDiffractiveCrossSections::waveOverlapT(double r, double z, double Q2){
+ double ef=quarkCharge[quarkIndex()];
+ double mf=quarkMass[quarkIndex()];
+ double epsilon2=z*(1-z)*Q2+mf*mf;
+ double epsilon=sqrt(epsilon2);
+ double besselK1=TMath::BesselK1(epsilon*r/hbarc);
+ double besselK0=TMath::BesselK0(epsilon*r/hbarc);
+
+ double term1=2*Nc/M_PI*alpha_em*ef*ef;
+ double term2=(z*z+(1-z)*(1-z))*epsilon2*besselK1*besselK1;
+ double term3=mf*mf*besselK0*besselK0;
+
+ return term1*(term2+term3); //GeV^2
+}
+
+double InclusiveDiffractiveCrossSections::waveOverlapL(double r, double z, double Q2){
+ double ef=quarkCharge[quarkIndex()];
+ double mf=quarkMass[quarkIndex()];
+ double epsilon2=z*(1-z)*Q2+mf*mf;
+ double epsilon=sqrt(epsilon2);
+ double besselK0=TMath::BesselK0(epsilon*r/hbarc);
+
+ double term1=8*Nc/M_PI*alpha_em*ef*ef;
+ double term2=Q2*z*z*(1-z)*(1-z)*besselK0*besselK0;
+
+ return term1*term2; //GeV^2
+}
+
+double InclusiveDiffractiveCrossSections::uiDipoleSigma(double* x, double* par){
+ double b=*x;
+ double r=par[0];
+ double xbj=par[1];
+ double dsigmadb2=0;
+ if(mA==1)
+ dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xbj);
+ else
+ dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xbj);
+ return 2*M_PI*b*dsigmadb2; //fm
+}
+
+void InclusiveDiffractiveCrossSections::setBetaMin(double val){ mBetaMin=val; }
Index: branches/inclusive-diffraction/src/InclusiveDiffractiveCrossSections.h
===================================================================
--- branches/inclusive-diffraction/src/InclusiveDiffractiveCrossSections.h (revision 0)
+++ branches/inclusive-diffraction/src/InclusiveDiffractiveCrossSections.h (revision 137)
@@ -0,0 +1,138 @@
+//
+// To calculate total cross-sections use
+// void calculate(double, double, double, double);
+// This integrates over z, it is a 5D integration. N.B. Slow!
+//
+// For exclusive final state use
+// void calculate(double, double, double, double, double);
+// This does not inteagrate over z, and is a much faster 2D integration
+//
+
+#ifndef InclusiveDiffractiveCrossSections_h
+#define InclusiveDiffractiveCrossSections_h
+
+class DipoleModel;
+class InclusiveDiffractiveCrossSections {
+
+ public:
+ InclusiveDiffractiveCrossSections();
+ ~InclusiveDiffractiveCrossSections();
+ InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections&);
+ InclusiveDiffractiveCrossSections& operator=(const InclusiveDiffractiveCrossSections&);
+
+
+ DipoleModel* dipoleModel() const;
+
+ void init();
+ void setZ(double);
+ void setQuarkIndex(int);
+ double z();
+ int quarkIndex();
+ double zlower();
+ void setNumberOfFlavours(unsigned int);
+
+ void calculate_dsigmadbeta_QQ(double, double, double);
+ void calculate_dsigmadbetadz_QQ(double, double, double, double);
+
+ double* dsigmadbetaT_QQ();
+ double* dsigmadbetaL_QQ();
+
+ void setBetaMin(double);
+ double betaMin();
+
+ void setRelativePrecisionOfIntegration(double);
+
+ //Total (non-diffractive) cross-sections:
+ double uiDipoleSigma(double*, double*);
+ double waveOverlapT(double, double, double);
+ double waveOverlapL(double, double, double);
+ double uiDsigmadrT(double*, double*);
+ double uiDsigmadrL(double*, double*);
+ double uiDsigmadzT(double*, double*);
+ double uiDsigmadzL(double*, double*);
+ void calculateTotal(double, double);
+
+ //QQG
+ void putAllQQGtoZero();
+ void calculate_dsigmadbeta_QQG(double, double, double);
+ void calculate_dsigmadbetadz_QQG(double, double, double, double);
+ double* dsigmadbetaT_QQG();
+
+ //QQG_MS
+ double uiQQG_MS(const double*);
+ double uiQQGdz_MS(const double*);
+ void FTQQG_MS(double, double);
+ void FTQQGdz_MS(double, double, double);
+ double* dsigmadbetaT_QQG_MS();
+
+ //QQG_GBW:
+ void FTQQG_GBW(double, double);
+ void FTQQGdz_GBW(double, double, double);
+ double dsigmadb2tilde(double, double, double);
+ double uiR_GBW(double*, double*);
+ double uiQQG_GBW(const double*);
+ double uiQQGdz_GBW(const double*);
+ double* dsigmadbetaT_QQG_GBW();
+ //QQG_GBW beta->0:
+ void FTQQG_GBWbeta0(double);
+ double uiQQG_GBWbeta0(const double*);
+ double uiR_GBWbeta0(double*, double*);
+ double* dsigmadbetaT_QQG_GBWbeta0();
+
+ double mTotalCST;
+ double mTotalCSL;
+ //End Total (non-diffractive) cross-sections:
+
+ double uiPhi0(double*, double*);
+ double uiPhi1(double*, double*);
+ bool setKinematicPoint(double, double, double, double);
+ private:
+ // bool setKinematicPoint(double, double, double, double);
+ double uiIntegralL(double*, double*);
+ double uiIntegralT(double*, double*);
+ double dPhi1db(double*, double*);
+ double dPhi0db(double*, double*);
+
+ double buiA(double*, double*);
+ double phiuiA(double*, double*);
+ double ruiA(double*, double*);
+ double zui(double*, double*);
+ double rui(double*, double*);
+
+ // double uiPhi0(double*, double*);
+ // double uiPhi1(double*, double*);
+
+ double kinematicPoint[6];
+ double mZ;
+ double mR;
+ double mB;
+ double mBeta;
+ DipoleModel* mDipoleModel;
+ int mQuarkIndex;
+ bool mIsInitialized;
+ double mInclusiveT[6];
+ double mInclusiveL[6];
+ double mInclusiveQQG[6];
+ double mInclusiveQQG_MS[6];
+ double mInclusiveQQG_GBW[6];
+ double mInclusiveQQG_GBWbeta0[6];
+ unsigned int mA;
+ unsigned int mNumberOfFlavours;
+ double mRelativePrecisionOfIntegration;
+
+ double mBetaMin;
+};
+
+inline DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel() const { return mDipoleModel; }
+inline double InclusiveDiffractiveCrossSections::z(){ return mZ; }
+inline int InclusiveDiffractiveCrossSections::quarkIndex(){ return mQuarkIndex; }
+inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQ(){ return mInclusiveT; }
+inline double* InclusiveDiffractiveCrossSections::dsigmadbetaL_QQ(){ return mInclusiveL; }
+
+inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG(){ return mInclusiveQQG; }
+inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG_MS(){ return mInclusiveQQG_MS; }
+inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG_GBW(){ return mInclusiveQQG_GBW; }
+inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG_GBWbeta0(){ return mInclusiveQQG_GBWbeta0; }
+
+inline double InclusiveDiffractiveCrossSections::betaMin(){ return mBetaMin; }
+#endif

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:31 PM (22 h, 26 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486634
Default Alt Text
(42 KB)

Event Timeline