Index: branches/fitting/src/InclusiveDiffractiveCrossSections_dz.cpp =================================================================== --- branches/fitting/src/InclusiveDiffractiveCrossSections_dz.cpp (revision 379) +++ branches/fitting/src/InclusiveDiffractiveCrossSections_dz.cpp (revision 380) @@ -1,469 +1,243 @@ #include "Math/SpecFunc.h" #include "TFile.h" #include #include #include #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" #include "DipoleModelParameters.h" +#define PR(x) cout << #x << " = " << (x) << endl; + void InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQ(double Q2, double W2, double beta, double z){ if(beta1) dipoleModel()->createSigma_ep_LookupTable(kinematicPoint[3]); - + double resultT=0; double resultL=0; for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) - integrate=false; - } - if(integrate){ - if(z0){ + if(qMass[i-1]==qMass[i]) + integrate=false; + } + if(integrate){ + if(z0 - cout<<"Warning, InclusiveDiffractiveCrossSections::calculate_dsigmadbetadz_QQG invalid kinematic point!"<1) dipoleModel()->createSigma_ep_LookupTable(xpom); double result=0; - double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius(); - + double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius(); + for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) - integrate=false; - } - if(integrate){ - //beta is the lower limit of z in this treatment of QQG: - if(z0){ + if(qMass[i-1]==qMass[i]) + integrate=false; + } + if(integrate){ + //beta is the lower limit of z in this treatment of QQG: + if(z1) dipoleModel()->createSigma_ep_LookupTable(xpom); - + double result=0; - double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius(); - + double rbhigh=upperIntegrationLimit*dipoleModel()->nucleus()->radius(); + for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - 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); + if(!mWhichFlavours[i]) continue; + 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.; + N_of_r=mDipoleModel->dsigmadb2ep(r, b, xpom)/2.; else - N_of_r=mDipoleModel->coherentDsigmadb2(r, b, xpom)/2.; - + 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.; + N_of_rp=mDipoleModel->dsigmadb2ep(rp, b, xpom)/2.; else - N_of_rp=mDipoleModel->coherentDsigmadb2(rp, b, xpom)/2.; - + 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.; + N_of_rMinusrp=mDipoleModel->dsigmadb2ep(rMinusrp, b, xpom)/2.; else - N_of_rMinusrp=mDipoleModel->coherentDsigmadb2(rMinusrp, b, xpom)/2.; - + 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!"<1) dipoleModel()->createSigma_ep_LookupTable(x); - - // dipoleModel()->setMu02(W2+Q2); - // dipoleModel()->setMu02(C2*Q2+C3); - - double resultT=0; - double resultL=0; - mTotalCST=0; - mTotalCSL=0; - for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - setQuarkIndex(i); - - double zlow=0; - double zhigh=1-zlow; - - mTotalT->SetParameter(0, Q2); - mTotalT->SetParameter(1, W2); - resultT=mGITotalT.Integral(zlow, zhigh); //GeV^2*fm^4 - - mTotalL->SetParameter(0, Q2); - mTotalL->SetParameter(1, W2); - resultL=mGITotalL.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::uiTotalDIST(const double* arg){ - double r=arg[0]; - double b=arg[1]; - - double Q2=mQ2; - double W2=mW2; - double z=mZ; - - double xbj=Kinematics::x(Q2, W2); - if(quarkIndex()>2) xbj=xbj*(1.+4.*pow(mParameters->quarkMass(quarkIndex()), 2)/Q2); //as in RSVV - - double waveOverlap=waveOverlapT(r, z, Q2); //GeV^2 - - double dsigmadb2=0; - if(mA==1) - dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xbj); - else - dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xbj); - - double result=waveOverlap*dsigmadb2/(4*M_PI); //GeV^2 - - //Do angular parts of r- and b-integrals: - result*=(2*M_PI)*b*(2*M_PI)*r; //fm^2*GeV^2 - - return result; -} - -double InclusiveDiffractiveCrossSections::uiTotalDISL(const double* arg){ - double r=arg[0]; - double b=arg[1]; - - double Q2=mQ2; - double W2=mW2; - double z=mZ; - - double xbj=Kinematics::x(Q2, W2); - if(quarkIndex()>2) xbj=xbj*(1.+4.*pow(mParameters->quarkMass(quarkIndex()), 2)/Q2); //as in RSVV - - double waveOverlap=waveOverlapL(r, z, Q2); //GeV^2 - - double dsigmadb2=0; - if(mA==1) - dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xbj); - else - dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xbj); - - double result=waveOverlap*dsigmadb2/(4*M_PI); //GeV^2 - result*=(2*M_PI)*b*(2*M_PI)*r; //fm^2*GeV^2 - - return result; -} - - -double InclusiveDiffractiveCrossSections::uiDsigmadzT(double* x, double* par){ - double z=*x; - double Q2=par[0]; - double W2=par[1]; - mQ2=Q2; - mW2=W2; - mZ=z; - - double epsilon2=z*(1-z)*Q2+pow(mParameters->quarkMass(quarkIndex()), 2); //GeV^2 - double epsilon=sqrt(epsilon2)/hbarc; //fm^-1 - - double rlow=1e-5; - double rhigh=mUpperRinIntegrationDIS/epsilon; //fm - double blow=0; - double bhigh=mUpperBinIntegration; - - double a[2]={rlow, blow}; //r, b - double b[2]={rhigh, bhigh}; - - - ROOT::Math::Functor wfT(this, &InclusiveDiffractiveCrossSections::uiTotalDIST, 2); - ROOT::Math::IntegratorMultiDim igT(ROOT::Math::IntegrationMultiDim::kADAPTIVE, 0., mRelativePrecisionOfIntegration, 1e8); - igT.SetFunction(wfT); - double resultT=igT.Integral(a, b); //fm^4*GeV^2 - - return resultT; -} - -double InclusiveDiffractiveCrossSections::uiDsigmadzL(double* x, double* par){ - double z=*x; - double Q2=par[0]; - double W2=par[1]; - mQ2=Q2; - mW2=W2; - mZ=z; - - double epsilon2=z*(1-z)*Q2+pow(mParameters->quarkMass(quarkIndex()), 2); //GeV^2 - double epsilon=sqrt(epsilon2)/hbarc; //fm^-1 - - double rlow=1e-5; - double rhigh=mUpperRinIntegrationDIS/epsilon; //fm - double blow=0; - double bhigh=mUpperBinIntegration; - - double a[2]={rlow, blow}; //r, b - double b[2]={rhigh, bhigh}; - - - ROOT::Math::Functor wfL(this, &InclusiveDiffractiveCrossSections::uiTotalDISL, 2); - ROOT::Math::IntegratorMultiDim igL(ROOT::Math::IntegrationMultiDim::kADAPTIVE, 0., mRelativePrecisionOfIntegration, 1e8); - igL.SetFunction(wfL); - double resultL=igL.Integral(a, b); //fm^4*GeV^2 - - return resultL; -} - -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); - if(quarkIndex()>2) xbj=xbj*(1.+4.*pow(mParameters->quarkMass(quarkIndex()), 2)/Q2); //as in RSVV - double blow=0; - double bhigh=mUpperBinIntegration; - mDipole->SetParameter(0, r); - mDipole->SetParameter(1, xbj); - double dipoleCS=mGIDipole.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); - if(quarkIndex()>2) xbj=xbj*(1.+4.*pow(mParameters->quarkMass(quarkIndex()), 2)/Q2); //as in RSVV - double blow=0; - double bhigh=mUpperBinIntegration; - - mDipole->SetParameter(0, r); - mDipole->SetParameter(1, xbj); - double dipoleCS=mGIDipole.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=mParameters->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=mParameters->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/fitting/src/DipoleModel.cpp =================================================================== --- branches/fitting/src/DipoleModel.cpp (revision 379) +++ branches/fitting/src/DipoleModel.cpp (revision 380) @@ -1,369 +1,751 @@ //============================================================================== // DipoleModel.cpp // -// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich +// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich // -// This file is part of Sartre. +// This file is part of Sartre version: 1.1 // // 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 . // // Author: Tobias Toll // Last update: // $Date$ // $Author$ //============================================================================== #include #include #include #include #include "DipoleModel.h" #include "TableGeneratorSettings.h" #include "DglapEvolution.h" +#include "Parameters_Minuit2.h" #include "Constants.h" #include "TFile.h" #include "TVector3.h" #include "TMath.h" #include "TH2F.h" #include "TF1.h" #include "Math/WrappedTF1.h" #include "Math/GaussIntegrator.h" - +#include "TObject.h" +#include "Minuit2/MnUserParameters.h" +//#include "LHAPDF/LHAPDF.h" + + #define PR(x) cout << #x << " = " << (x) << endl; - + using namespace std; - + DipoleModel::DipoleModel() { - mConfigurationExists=false; - TableGeneratorSettings* settings = TableGeneratorSettings::instance(); - unsigned int A = settings->A(); - mNucleus.init(A); - mIsInitialized=true; - mParameters = 0; + mConfigurationExists=false; + TableGeneratorSettings* settings = TableGeneratorSettings::instance(); + unsigned int A = settings->A(); + mNucleus.init(A); + mIsInitialized=true; + //comparisons with previous version in ep yields + //best agreement with as(MZ)=0.1196 at NLO: + // mAs.init(0.1196, 2); + // mAs.init(0.13939, 1); //Default Pythia 8 Too large + // mAs.init(params_coupling, params_asOrder); + // mAs.init(0.126, 1);// + + // + // Initialise DGLAP evolution + // + cout << "DipoleModel(): 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.); + mUseSkewednessCorrection=false; } - -DipoleModel::~DipoleModel() -{ - delete mParameters; -} - + +DipoleModel::~DipoleModel() { /* no op */} + const TableGeneratorNucleus* DipoleModel::nucleus() const { return &mNucleus; } - + bool DipoleModel::configurationExists() const { return mConfigurationExists; } - + double DipoleModel::bDependence(double) { return 0; } double DipoleModel::bDependence(double, double) { return 0; } double DipoleModel::dsigmadb2ep(double, double, double) { return 0;} - -void DipoleModel::createCompressedExponentThickness(){ - double zbhigh=30;//30.*(1/muCE)*pow(hbarc, 1); //fm - double rhigh=sqrt(1/1.5)/hbarc; - int numBinsb=10000; - int numBinsr=1000; - - if(mCE_lu) delete mCE_lu; - mCE_lu=new TH2D("CE_lu", "CE_lu", numBinsb, 0, zbhigh, numBinsr, 0, rhigh); //r Upper limit can be changed - double binWidthb=zbhigh/numBinsb; //fm - double norm=0; - for(int i=1; i<=numBinsb; i++){ - for(int j=1; j<=numBinsr; j++){ - double b=mCE_lu->GetXaxis()->GetBinCenter(i); //fm - double r=mCE_lu->GetYaxis()->GetBinCenter(j); //fm - TF1 CE("CE", this, &DipoleModel::uiCompressedExponent, -zbhigh, zbhigh, 4); - CE.SetParameter(0, b); - CE.SetParameter(1, r); - CE.SetParameter(2, mParameters->BG()); - CE.SetParameter(3, mParameters->v()); //exponent parameter - ROOT::Math::WrappedTF1 wf(CE); - ROOT::Math::GaussIntegrator ig; - ig.SetFunction(wf); - ig.SetRelTolerance(1e-6); - double CE_of_b=2.*ig.IntegralUp(0.)/hbarc;//(-zbhigh, zbhigh)/hbarc; //GeV^2 - mCE_lu->SetBinContent(i, j, CE_of_b); - // - // Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db) - // ~2pi*sum(T'(bi)*bi*Delta b) - // - norm+=2.*M_PI*b/hbarc*CE_of_b*binWidthb/hbarc; //GeV^0 - } - } - PR(norm); - mCE_lu->Scale(1./norm); //normalise, GeV^2 -}//BS -double DipoleModel::uiCompressedExponent(double* x, double* par){ - //All in fm - double z=*x; //fm - double b=par[0]; - double r=par[1]; //fm - double u=par[2]; //Radius GeV^-1 - double exponentparameter=par[3]; - //PR(exponentparameter); - //double v=2-exp(-exponentparameter*1e-2/(pow(r,6.))); - //double v=1+exp(-exponentparameter*r*r); - double v=2-exp(-exponentparameter*1e-4/(pow(r,4.))); - //double v=1+exp(-exponentparameter*pow(r,4.)); - //double v=max(2.95673-exponentparameter*r, 1.); - //double v=1+exp(-exponentparameter*pow(r,9)); - //double v=par[2]; - //double v=exponentparameter; - // double arg=u*(pow((sqrt(z*z+b*b)/hbarc),v)); - double arg=u*(pow((z*z+b*b)/hbarc2, 0.5*v)); - return pow(u,3./v)*v / (4.*M_PI*TMath::Gamma(3./v)) * exp(-arg); //GeV^ +double DipoleModel::ddxdsigmadb2(double, double, double) { return 0;} + +//void DipoleModel::setParams(const std::vector& par) {/*this needs to return nothing*/} + +//***********bSat:***************************************************** + +void DipoleModel::setParams(const std::vector& par) +{ + + mBg = par[0]; //par[0];//BS For MZ, BG is kept fixed + mLmass = par[1]; + mCmass = par[2]; + mBmass = par[3]; + mC = par[4]; + mMu_02 = par[5]; //par[5];//BS For MZ, Mu02 is kept fixed + mAg = par[6]; + mLamg = par[7]; + mRmax = par[8]; + mAs_qqg = par[9]; + mRhs = par[10]; + mnLP = par[11]; + mR0WS = par[12]; + mdWS = par[13]; + muCE = par[14]; + mEP = par[15]; + BV = par[16]; + RV = par[17]; + wV = par[18]; + + qMass[0] = mLmass; + qMass[1] = mLmass; + qMass[2] = mLmass; + qMass[3] = mCmass; + qMass[4] = mBmass; + qMass[5] = params_t_mass; + //cout << "DipoleModel mAg = " + std::to_string(mAg) << endl; + + } -//***********bSat:***************************************************** DipoleModel_bSat::DipoleModel_bSat() { - mBDependence = 0; + // mPDF=LHAPDF::mkPDF("HERAPDF20_LO_EIG", 0); //LHAPDF + mBDependence = 0; mSigma_ep_LookupTable = 0; - - // - // Set the parameters. Note that we enforce here the bSat model - // independent of what the settings say. - // - TableGeneratorSettings* settings = TableGeneratorSettings::instance(); - mParameters = new DipoleModelParameters(bSat, settings->dipoleModelParameterSet()); -} - + mThicknessNormTable = 0; +} + DipoleModel_bSat::~DipoleModel_bSat() -{ - delete mSigma_ep_LookupTable; - delete mBDependence; +{ + // delete mPDF; + delete mSigma_ep_LookupTable; + delete mThicknessNormTable; + delete mBDependence; } - + DipoleModel_bSat& DipoleModel_bSat::operator=(const DipoleModel_bSat& dp) { - if (this != &dp) { - delete mBDependence; - delete mSigma_ep_LookupTable; - - DipoleModel::operator=(dp); - mBDependence = new TH2F(*(dp.mBDependence)); - mBDependence->SetDirectory(0); - mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable)); - mSigma_ep_LookupTable->SetDirectory(0); - } - return *this; + if (this != &dp) { + delete mBDependence; + delete mSigma_ep_LookupTable; + delete mThicknessNormTable; + + DipoleModel::operator=(dp); + mBDependence = new TH2F(*(dp.mBDependence)); + mBDependence->SetDirectory(0); + mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable)); + mThicknessNormTable = new TH1F(*(dp.mThicknessNormTable)); + mSigma_ep_LookupTable->SetDirectory(0); + mThicknessNormTable->SetDirectory(0); + + } + return *this; } - + DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp) { - mBDependence = new TH2F(*(dp.mBDependence)); - mBDependence->SetDirectory(0); + mBDependence = new TH2F(*(dp.mBDependence)); + mBDependence->SetDirectory(0); } - + +void DipoleModel::createLaplaceThickness(){ + double zbhigh=30.*mnLP*hbarc; //fm + int numBins=100000; + PR(zbhigh); + mLP_lu=new TH1D("LP_lu", "LP_lu", numBins, 0, zbhigh); + double binWidth=zbhigh/numBins; //fm + double norm=0; + for(int i=1; i<=numBins; i++){ + double b=mLP_lu->GetBinCenter(i); //fm + TF1 LP("LP", this, &DipoleModel_bSat::uiLaplace, -zbhigh, zbhigh, 3); + LP.SetParameter(0, b); + LP.SetParameter(1, mnLP); + //WS.SetParameter(2, params_dWS); + ROOT::Math::WrappedTF1 wf(LP); + ROOT::Math::GaussIntegrator ig; + ig.SetFunction(wf); + ig.SetRelTolerance(1e-3); + double LP_of_b=ig.Integral(-zbhigh, zbhigh)/hbarc; //GeV^2 + mLP_lu->SetBinContent(i, LP_of_b); + // + // Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db) + // ~2pi*sum(T'(bi)*bi*Delta b) + // + norm+=2.*M_PI*b/hbarc*LP_of_b*binWidth/hbarc; //GeV^0 + } + PR(norm); + mLP_lu->Scale(1./norm); //normalise, GeV^2 + // + // Calculate =integral(b^2 T(b)d^2b)=integral(2pi b^3T(b)db) + // + // PR(mLP_lu->Interpolate(0.5)); + // PR(mLP_lu->Interpolate(1)); + // PR(mLP_lu->Interpolate(2)); + + // PR((1/(4*M_PI*mnLP*mnLP))*0.5/(mnLP*hbarc)*TMath::BesselK1(0.5/(mnLP*hbarc))); + // PR((1/(4*M_PI*mnLP*mnLP))*1/(mnLP*hbarc)*TMath::BesselK1(1/(mnLP*hbarc))); + // PR((1/(4*M_PI*mnLP*mnLP))*2/(mnLP*hbarc)*TMath::BesselK1(2/(mnLP*hbarc))); + + double b2_mean=0.; + for(int i=1; i<=numBins; i++){ + double b=mLP_lu->GetBinCenter(i)/hbarc; + double T_of_b=mLP_lu->GetBinContent(i); + b2_mean+=b*b*b*T_of_b; //GeV-1 + } + b2_mean*=2.*M_PI*binWidth/hbarc; //GeV^-2 + PR(sqrt(b2_mean*hbarc2)); +}//BS + + +double DipoleModel::uiLaplace(double* x, double* par){ + //All in fm + double z=*x; //fm + double b=par[0]; //fm + double r=par[1]; //Radius GeV^-1 + //double d=par[2]; //thickness GeV^-1 + + double arg=(sqrt((z*z+b*b))/(hbarc*r)); + return (1/(8*M_PI*r*r*r))*exp(-arg); //GeV^3 +} + +void DipoleModel::createCompressedExponentThickness(){ + double zbhigh=30;//30.*(1/muCE)*pow(hbarc, 1); //fm + double rhigh=sqrt(1/1.5)/hbarc; + int numBinsb=10000; + int numBinsr=1000; + // int numBinsb=100; + // int numBinsr=10; + PR(zbhigh); + PR(rhigh); + mCE_lu=new TH2D("CE_lu", "CE_lu", numBinsb, 0, zbhigh, numBinsr, 0, rhigh); //r Upper limit can be changed + //mCE_lu=new TH1D("CE_lu", "CE_lu", numBinsb, 0, zbhigh);//, numBinsr, 0, rhigh); //r Upper limit can be changed + double binWidthb=zbhigh/numBinsb; //fm + double binWidthr=rhigh/numBinsr; + double norm=0; + PR(mEP); + PR(muCE); + for(int i=1; i<=numBinsb; i++){ + for(int j=1; j<=numBinsr; j++){ + double b=mCE_lu->GetXaxis()->GetBinCenter(i); //fm + double r=mCE_lu->GetYaxis()->GetBinCenter(j); //fm + TF1 CE("CE", this, &DipoleModel::uiCompressedExponent, -zbhigh, zbhigh, 4); + CE.SetParameter(0, b); + CE.SetParameter(1, r); + CE.SetParameter(2, muCE); + CE.SetParameter(3, mEP); + //if(r>0.8) PR(r); + ROOT::Math::WrappedTF1 wf(CE); + ROOT::Math::GaussIntegrator ig; + ig.SetFunction(wf); + ig.SetRelTolerance(1e-6); + double CE_of_b=2.*ig.IntegralUp(0.)/hbarc;//(-zbhigh, zbhigh)/hbarc; //GeV^2 + //PR(CE_of_b); + // PR(b); + // PR(r); + // PR(muCE); + // PR(mEP); + // PR(zbhigh); + // PR(rhigh); + // PR(CE_of_b); + mCE_lu->SetBinContent(i, j, CE_of_b); + // + // Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db) + // ~2pi*sum(T'(bi)*bi*Delta b) + // + norm+=2.*M_PI*b/hbarc*CE_of_b*binWidthb/hbarc; //GeV^0 + } + } + // PR(norm); + // mCE_lu->Scale(1./norm); //normalise, GeV^2 + // + // Calculate =integral(b^2 T(b)d^2b)=integral(2pi b^3T(b)db) + // // + // double* z=0; + // double parameters[3]={1, muCE, mEP}; + + // PR(DipoleModel::uiCompressedExponent(z, parameters)); + // PR(mCE_lu->Interpolate(0.5, 0));//, 1.)); + // PR(mCE_lu->Interpolate(1, 0));//, 2.)); + // PR(mCE_lu->Interpolate(2, 0));//, 4.)); + + // PR(((pow(muCE,3/mEP)*mEP)/(4*M_PI*TMath::Gamma(3/mEP)))); + + // PR((1/pow(2*M_PI*mBg, 1.5))); + + // // double u=muCE; + // // double v=mEP; + // // double z=0; + // // double b=0.5; + // // double arg=u*(pow((sqrt(z*z+b*b)/hbarc),v)); + // // PR(((pow(u,3/v)*v)/(4*M_PI*TMath::Gamma(3/v)))*exp(-arg)); + // // b=1; + // // arg=u*(pow((sqrt(z*z+b*b)/hbarc),v)); + // // PR(((pow(u,3/v)*v)/(4*M_PI*TMath::Gamma(3/v)))*exp(-arg)); + // // b=2; + // // arg=u*(pow((sqrt(z*z+b*b)/hbarc),v)); + // // PR(((pow(u,3/v)*v)/(4*M_PI*TMath::Gamma(3/v)))*exp(-arg)); + + // PR((1/pow(2*M_PI*mBg, 1.5)) * exp(-(0.5*0.5)/(hbarc2*2*mBg))); + // PR((1/pow(2*M_PI*mBg, 1.5)) * exp(-(1*1)/(hbarc2*2*mBg))); + // PR((1/pow(2*M_PI*mBg, 1.5)) * exp(-(2*2)/(hbarc2*2*mBg))); + + // double zvalue=0; + // double values1[4]={0.5, 0., muCE, mEP}; + // double values2[4]={1., 0., muCE, mEP}; + // double values3[4]={2., 0., muCE, mEP}; + // PR(uiCompressedExponent(&zvalue, values1)); + // PR(uiCompressedExponent(&zvalue, values2)); + // PR(uiCompressedExponent(&zvalue, values3)); + // // PR((1/(4*M_PI*mnLP*mnLP))*0.5/(mnLP*hbarc)*TMath::BesselK1(0.5/(mnLP*hbarc))); + // // PR((1/(4*M_PI*mnLP*mnLP))*1/(mnLP*hbarc)*TMath::BesselK1(1/(mnLP*hbarc))); + // // PR((1/(4*M_PI*mnLP*mnLP))*2/(mnLP*hbarc)*TMath::BesselK1(2/(mnLP*hbarc))); + + // // PR(1/(2*M_PI*mBg) * exp(-0*0/((hbarc2)*(2*mBg)))); + // PR((1/pow(2*M_PI*mBg, 1)) * exp(-(0.5*0.5)/(hbarc2*2*mBg))); + // PR((1/pow(2*M_PI*mBg, 1)) * exp(-(1*1)/(hbarc2*2*mBg))); + // PR((1/pow(2*M_PI*mBg, 1)) * exp(-(2*2)/(hbarc2*2*mBg))); + // // /double b2_mean=0.; + // // for(int i=1; i<=numBins; i++){ + // // double b=mCE_lu->GetBinCenter(i)/hbarc; + // // double T_of_b=mCE_lu->GetBinContent(i); + // // b2_mean+=b*b*b*T_of_b; //GeV-1 + // // } + // // b2_mean*=2.*M_PI*binWidth/hbarc; //GeV^-2 + // // PR(sqrt(b2_mean*hbarc2)); +}//BS + + +double DipoleModel::uiCompressedExponent(double* x, double* par){ + //All in fm + double z=*x; //fm + double b=par[0]; + double r=par[1]; //fm + double u=par[2]; //Radius GeV^-1 + double exponentparameter=par[3]; + //PR(exponentparameter); + //double v=2-exp(-exponentparameter*1e-2/(pow(r,6.))); + //double v=1+exp(-exponentparameter*r*r); + double v=2-exp(-exponentparameter*1e-4/(pow(r,4.))); + //double v=1+exp(-exponentparameter*pow(r,4.)); + //double v=max(2.95673-exponentparameter*r, 1.); + //double v=1+exp(-exponentparameter*pow(r,9)); + //double v=par[2]; + //double v=exponentparameter; + // double arg=u*(pow((sqrt(z*z+b*b)/hbarc),v)); + double arg=u*(pow((z*z+b*b)/hbarc2, 0.5*v)); + return pow(u,3./v)*v / (4.*M_PI*TMath::Gamma(3./v)) * exp(-arg); //GeV^ +} + +void DipoleModel::createWoodsSaxonThickness(){ + double zbhigh=3.*params_R0WS*hbarc; //fm + int numBins=100000; + + mWS_lu=new TH1D("WS_lu", "WS_lu", numBins, 0, zbhigh); + double binWidth=zbhigh/numBins; //fm + double norm=0; + for(int i=1; i<=numBins; i++){ + double b=mWS_lu->GetBinCenter(i); //fm + TF1 WS("WS", this, &DipoleModel_bSat::uiWoodsSaxon, -zbhigh, zbhigh, 3); + WS.SetParameter(0, b); + WS.SetParameter(1, params_R0WS); + WS.SetParameter(2, params_dWS); + ROOT::Math::WrappedTF1 wf(WS); + ROOT::Math::GaussIntegrator ig; + ig.SetFunction(wf); + ig.SetRelTolerance(1e-2); + double WS_of_b=ig.Integral(-zbhigh, zbhigh)/hbarc; //GeV^-1 + mWS_lu->SetBinContent(i, WS_of_b); + // + // Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db) + // ~2pi*sum(T'(bi)*bi*Delta b) + // + norm+=2.*M_PI*b/hbarc*WS_of_b*binWidth/hbarc; //GeV^-3 + } + PR(norm); + mWS_lu->Scale(1./norm); //normalise, GeV^2 + // + // Calculate =integral(b^2 T(b)d^2b)=integral(2pi b^3T(b)db) + // + double b2_mean=0.; + for(int i=1; i<=numBins; i++){ + double b=mWS_lu->GetBinCenter(i)/hbarc; + double T_of_b=mWS_lu->GetBinContent(i); + b2_mean+=b*b*b*T_of_b; //GeV-1 + } + b2_mean*=2.*M_PI*binWidth/hbarc; //GeV^-2 + PR(sqrt(b2_mean)); +} + + +double DipoleModel::uiWoodsSaxon(double* x, double* par){ + //All in fm + double z=*x; //fm + double b=par[0]; //fm + double R0=par[1]; //Radius GeV^-1 + double d=par[2]; //thickness GeV^-1 + + double arg=(sqrt((z*z+b*b))/hbarc-R0)/d; + return 1./(1+exp(arg)); //GeV^0 +} + void DipoleModel_bSat::createConfiguration(int iConfiguration) { - if (!mIsInitialized) { - cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl; - exit(1); - } - TableGeneratorSettings* settings = TableGeneratorSettings::instance(); - unsigned int A = mNucleus.A(); - string path=settings->bSatLookupPath(); - ostringstream filename; - filename.str(""); - filename << path << "/bSat_bDependence_A" << A <<".root"; - ifstream ifs(filename.str().c_str()); - if (!ifs) { - cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl; - cout << "Stopping." << endl; - exit(1); - } - TFile* lufile= new TFile(filename.str().c_str()); - ostringstream histoName; - histoName.str( "" ); - histoName << "Configuration_" << iConfiguration; - lufile->GetObject( histoName.str().c_str(), mBDependence ); - mBDependence->SetDirectory(0); - lufile->Close(); - mConfigurationExists=true; + if (!mIsInitialized) { + cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl; + exit(1); + } + TableGeneratorSettings* settings = TableGeneratorSettings::instance(); + unsigned int A = mNucleus.A(); + string path=settings->bSatLookupPath(); + ostringstream filename; + filename.str(""); + filename << path << "/bSat_bDependence_A" << A <<".root"; + ifstream ifs(filename.str().c_str()); + if (!ifs) { + cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl; + cout << "Stopping." << endl; + exit(1); + } + TFile* lufile= new TFile(filename.str().c_str()); + ostringstream histoName; + histoName.str( "" ); + histoName << "Configuration_" << iConfiguration; + lufile->GetObject( histoName.str().c_str(), mBDependence ); + mBDependence->SetDirectory(0); + lufile->Close(); + mConfigurationExists=true; } - + + double DipoleModel_bSat::dsigmadb2(double r, double b, double phi, double xprobe) { - double bDep=bDependence(b, phi); - double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02(); + double bDep=bDependence(b, phi); + double muQ2 = mC/(r*r/hbarc2) + mMu_02; double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; - double result = 2.*(1. - exp(-omega/2)); - return result; + double result = 2.*(1. - exp(-omega/2)); + return result; } - + +double DipoleModel_bSat::ddxdsigmadb2(double r, double b, double xprobe){ + /*#TT + const double BG = mBg; // GeV^-2 + double arg = b*b/(2*BG); + arg /= hbarc2; + double bDep= 1/(2*M_PI*BG) * exp(-arg); + double muQ2 = mC/(r*r/hbarc2) + mMu_02; + double alphas=mAs.alphaS(muQ2); + double Gofx=DglapEvolution::instance().G(xprobe, muQ2); + // double Gofx=mPDF->xfxQ2(21, xprobe, muQ2)/xprobe; + double omegaFactor = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*alphas*bDep; + double h=xprobe*1e-3; + double Gofxp=DglapEvolution::instance().G(xprobe+h, muQ2); + double Gofxm=DglapEvolution::instance().G(xprobe-h, muQ2); + // double Gofxp=mPDF->xfxQ2(21, xprobe+h, muQ2)/(xprobe+h); + // double Gofxm=mPDF->xfxQ2(21, xprobe-h, muQ2)/(xprobe-h); + double ddx_xg=((xprobe+h)*Gofxp-(xprobe-h)*Gofxm)/(2*h); + double result = omegaFactor*ddx_xg*exp(-omegaFactor*xprobe*Gofx/2); + return result; + */ + return 0.; +} + double DipoleModel_bSat::bDependence(double b, double phi) -{ - double result = mBDependence->Interpolate(b, phi); - return result; -} - -double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe) -{ - const double BG = mParameters->BG(); // GeV^-2 - double arg = b*b/(2*BG); - arg /= hbarc2; - double bDep= 1/(2*M_PI*BG) * exp(-arg); - double Mu02 = mParameters->mu02(); // GeV^2 - double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02; - double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); - double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; - double result = 2.*(1. - exp(-omega/2)); - - return result; -} - -double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe){ - xprobe*=1.; - double sigmap = mSigma_ep_LookupTable->Interpolate(r); - int A=nucleus()->A(); - double TA=nucleus()->T(b)/A; - double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) ); - return result; -} - -void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe) -{ - double rbRange=3.*nucleus()->radius(); - TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2); - mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange); - dsigmaForIntegration->SetNpx(1000); - for (int iR=1; iR<=1000; iR++) { - double r=mSigma_ep_LookupTable->GetBinCenter(iR); - dsigmaForIntegration->SetParameter(0, r); - dsigmaForIntegration->SetParameter(1, xprobe); - double result=dsigmaForIntegration->Integral(0, rbRange); - mSigma_ep_LookupTable->SetBinContent(iR, result); - } - delete dsigmaForIntegration; -} - -double DipoleModel_bSat::dsigmadb2epForIntegration(double *x, double* par) -{ - double r=par[0]; - double xprobe=par[1]; - double b = *x; - return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe); +{ + double result = mBDependence->Interpolate(b, phi); + return result; } - - -//***********bNonSat:************************************************* -DipoleModel_bNonSat::DipoleModel_bNonSat() + + +double DipoleModel::protonThickness(double b, double r=0){ + double result=0; + switch(params_Shape){ + case HardSphere: + // + // T(b)=3/2pi*1/R^3*sqrt(R^2-b^2) + // + // =0, =2/5*R^2 + // sigma(b)=0.63*R + // + if(b/hbarcInterpolate(b, r); //GeV^2 + break; + case Gauss: + default: + // + // T(b)=1/(2piBG)exp(-b2/(2BG)) + // + // =0, =BG + // sigma(b)=sqrt(BG) + // + const double BG = mBg; // GeV^-2 + double arg = b*b/(2*BG); //fm2*GeV^2 + arg /= hbarc2; //GeV^0 + result= 1/(2*M_PI*BG) * exp(-arg); //GeV^2 + break; + } + return result; +} + +void DipoleModel_bSat::setMu02(double val){ mMu02=val; } + +void DipoleModel_bSat::useSkewednessCorrection(){ mUseSkewednessCorrection=true; } + +double DipoleModel_bSat::calculateNorm(double w){ + double norm=0; + double step=1e-2; + mBg=4; + mnLP=1.8282; + + // double T_old=0; + for(double b=step; b<3.5; b+=step){ + double T_gauss=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg); + double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*TMath::BesselK1(b/(mnLP*hbarc)); + double T=(1-w)*T_gauss+w*T_laplace; + // T=0.5*(T+T_old); + // T_old=T; + norm+=2*M_PI*b*T/hbarc2*step; + } + return norm; +} + +void DipoleModel_bSat::createNormTable(){ + mBg=4; //#TT + mnLP=1.8282; //#TT + mThicknessNormTable = new TH1F("", "", 1000, -5., 5.); + for(int i=1; i<=1000; i++){ + //Omega = X(r, b_ * T(b), this is a function of X(r, xprobe) + double X=pow(10, mThicknessNormTable->GetBinCenter(i)); + double norm=0; + double db=1e-3; + for(double b=1e-4; b<4.; b+=db){ + //find weight: + double T_gauss=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg); + double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*TMath::BesselK1(b/(mnLP*hbarc)); + double weight=0.5; + double dsigdb2=0; + int count=0; + int countMax=1e4; + double bDep=0; + while(true){ + double weight_old=weight; + bDep= (1-weight)*T_gauss+weight*T_laplace; + double omega = X*bDep; + dsigdb2 = 2.*(1. - exp(-omega/2)); + weight=0.5*dsigdb2; + if(abs(weight_old-weight)/weight < 1e-5 or weight==0 or count>countMax) break; + } + norm+=2*M_PI*b*bDep*db/hbarc2; + } + mThicknessNormTable->SetBinContent(i, norm); + } +} + +// +// Version which interpolates between different shapes: +// +//double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe) +//{ +// // +// //modify dipole size as http://arxiv.org/abs/0807.0325v1 +// // +// //r=mRmax*sqrt(log(1+r*r/(mRmax*mRmax))); +// mBg=4; +// mnLP=1.8282; +// double T_gauss=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg); +// double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*TMath::BesselK1(b/(mnLP*hbarc)); +// double muQ2 = mC/(r*r/hbarc2) + mMu_02; //GeV2 +// double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); +// +// //Use an iterative approach to find the thickness +// double weight=0.5; +// double dsigdb2=0; +// int count=0; +// int countMax=1e4; +// double bDep=0; +// while(true){ +// double weight_old=weight; +// bDep= (1-weight)*T_gauss+weight*T_laplace; +// double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; +// dsigdb2 = 2.*(1. - exp(-omega/2)); +// weight=0.5*dsigdb2; +// if(abs(weight_old-weight)/weight < 1e-5 or weight==0 or count>countMax) break; +// } +// //We need to normalise the thickness +// // double norm=calculateNorm(weight); +// double X=((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg; +// // bDep/=mThicknessNormTable->Interpolate(log10(X)); +// dsigdb2 = 2.*(1. - exp(-X*bDep/2)); +// +// return bDep; +// return dsigdb2; +//} + +//original version: +double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe) { // - // Set the parameters. Note that we enforce here the bNonSat model - // independent of what the settings say. + //modify dipole size as http://arxiv.org/abs/0807.0325v1 // - TableGeneratorSettings* settings = TableGeneratorSettings::instance(); - mParameters = new DipoleModelParameters(bNonSat, settings->dipoleModelParameterSet()); -} - -DipoleModel_bNonSat::~DipoleModel_bNonSat(){} - - -double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe) -{ - const double BG = mParameters->BG(); // GeV^-2 - double arg = b*b/(2*BG); - arg /= hbarc2; - double bDep= 1/(2*M_PI*BG) * exp(-arg); - double Mu02 = mParameters->mu02(); // GeV^2 - double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02; - double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); - double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; - - return omega; -} - -double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe) -{ - double bDep=bDependence(b, phi); - double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02(); + //r=mRmax*sqrt(log(1+r*r/(mRmax*mRmax))); + double bDep=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg); + double muQ2 = mC/(r*r/hbarc2) + mMu_02; //GeV2 double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; - - return omega; -} - -double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){ - int A=nucleus()->A(); - double TA=nucleus()->T(b)/A; - double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02(); - double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); - double result=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg; - return result; -} - - -//***********bCGC:***************************************************** + double dsigdb2 = 2.*(1. - exp(-omega/2)); + return dsigdb2; +} + +//Skewedness: x1+x2=xp +// sigma(xp)^2 = \int_0^xp dx1\int_0^xp dx2 sigma(x1)sigma(x2)delta(xp-x1-x2)= +// = \int_0^xp sigma(x1)sigma(xp-x1) + + +double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe){ + xprobe*=1.; + double sigmap = mSigma_ep_LookupTable->Interpolate(r); + int A=nucleus()->A(); + double TA=nucleus()->T(b)/A; + double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) ); + return result; +} + +void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe) +{ + double rbRange=3.*nucleus()->radius(); + TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2); + mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange); + dsigmaForIntegration->SetNpx(1000); + for (int iR=1; iR<=1000; iR++) { + double r=mSigma_ep_LookupTable->GetBinCenter(iR); + dsigmaForIntegration->SetParameter(0, r); + dsigmaForIntegration->SetParameter(1, xprobe); + double result=dsigmaForIntegration->Integral(0, rbRange); + mSigma_ep_LookupTable->SetBinContent(iR, result); + } + delete dsigmaForIntegration; +} + +double DipoleModel_bSat::dsigmadb2epForIntegration(double *x, double* par) +{ + double r=par[0]; + double xprobe=par[1]; + double b = *x; + return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe); +} + + +//***********bNonSat:************************************************* +DipoleModel_bNonSat::DipoleModel_bNonSat(){} + +DipoleModel_bNonSat::~DipoleModel_bNonSat(){} -DipoleModel_bCGC::DipoleModel_bCGC() +double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe) { // - // Set the parameters. Note that we enforce here the bNonSat model - // independent of what the settings say. + //modify dipole size as http://arxiv.org/abs/0807.0325v1 // - TableGeneratorSettings* settings = TableGeneratorSettings::instance(); - mParameters = new DipoleModelParameters(bCGC, settings->dipoleModelParameterSet()); + r=mRmax*sqrt(log(1+r*r/(mRmax*mRmax))); + //const double BG = bNonSat_BG; // GeV^-2 + //double arg = b*b/(2*BG); + //arg /= hbarc2; + // double bDep= 1/(2*M_PI*BG) * exp(-arg); + double bDep= protonThickness(b, r); + double Mu02 = mMu_02; // GeV^2 + double muQ2 = mC/(r*r/hbarc2) + Mu02; + double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); + double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; + + return omega; +} + +double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe) +{ + double bDep=bDependence(b, phi); + double muQ2 = params_C/(r*r/hbarc2) + mMu_02; + double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); + double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep; + return omega; } +double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){ + int A=nucleus()->A(); + double TA=nucleus()->T(b)/A; + double muQ2 = params_C/(r*r/hbarc2) + mMu_02; + double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2); + double result=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg; + return result; +} -void DipoleModel_bCGC::createConfiguration(int iConfiguration) -{ - if (!mIsInitialized) { - cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl; - exit(1); - } - //iConfiguration is a dummy for bCGC. Do this to avoid compiler warnings: - iConfiguration++; - mNucleus.generate(); - mConfigurationExists=true; -} - -double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x) -{ - double result=1; - for (unsigned int iA=0; iAkappa(); - double N0 = mParameters->N0(); - double x0 = mParameters->x0(); - double lambda = mParameters->lambda(); - double gammas = mParameters->gammas(); - double A = -N0*N0*gammas*gammas/((1-N0)*(1-N0)*log(1-N0)); - double B = 0.5*pow(1-N0,-(1-N0)/(N0*gammas)); - double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b)); - double rQs = r*Qs/hbarc; - double result=0; - if (rQs <= 2) - result = 2*N0*pow(0.5*rQs, 2*(gammas+(1/(kappa*lambda*Y))*log(2/rQs))); - else - result = 2*(1 - exp(-A*log(B*rQs)*log(B*rQs))); - return result; -} - -double DipoleModel_bCGC::bDependence(double b) -{ - double gammas = mParameters->gammas(); - double Bcgc = mParameters->Bcgc(); - return exp(-0.5*b*b/Bcgc/gammas/hbarc2); -} - +// +////***********bCGC:***************************************************** +// +//void DipoleModel_bCGC::createConfiguration(int iConfiguration) +//{ +// if (!mIsInitialized) { +// cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl; +// exit(1); +// } +// //iConfiguration is a dummy for bCGC. Do this to avoid compiler warnings: +// iConfiguration++; +// mNucleus.generate(); +// mConfigurationExists=true; +//} +// +//double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x) +//{ +// double result=1; +// for (unsigned int iA=0; iA. // // Author: Thomas Ullrich // $Date$ // $Author: ullrich $ -//============================================================================== +//==============================================================================DipoleModelParameters::quarkMass(3) #ifndef DipoleModelParameters_h #define DipoleModelParameters_h #include "Enumerations.h" #include #include using namespace std; class Settings; class DipoleModelParameters { public: DipoleModelParameters(Settings*); DipoleModelParameters(DipoleModelType, DipoleModelParameterSet); void setDipoleModelType(DipoleModelType); void setDipoleModelParameterSet(DipoleModelParameterSet); DipoleModelType dipoleModelType() const; DipoleModelParameterSet dipoleModelParameterSet() const; // bSat, bNonSat double mu02() const; double lambdaG() const; double Ag() const; double C() const; double Rshrink() const; //proton parameters double BG() const; //Size parameter (u in compr. exp.) ProtonShape shape() const; double v() const; //2 for Gauss, 1 for Laplace // bCGC double kappa() const; double N0() const; double x0() const; double lambda() const; double gammas() const; double Bcgc() const; double quarkMass(unsigned int) const; double boostedGaussianR2(int vmID); double boostedGaussianNL(int vmID); double boostedGaussianNT(int vmID); double boostedGaussianQuarkMass(int vmID); void setBG(double); void setMu02(double); void setLambdaG(double); void setAg(double); void setC(double); void setRshrink(double); void setV(double); void setQuarkMass(int, double); bool list(ostream& = cout); private: void setupParameters(); void setup_bSat(); void setup_bNonSat(); void setup_bCGC(); void setup_boostedGaussiansWaveFunction(); private: DipoleModelType mDipoleModelType; string mDipoleModelParameterSetName; DipoleModelParameterSet mDipoleModelParameterSet; double mQuarkMass[6]; // bSat, bNonSat double mBG; double mMu02; double mLambdaG; double mAg; double mC; ProtonShape mShape; double mV; double mRshrink; // bCGC double mKappa; double mN0; double mX0; double mLambda; double mGammas; double mBcgc; // boosted Gaussian wave function parameters double mBoostedGaussianR2_rho; double mBoostedGaussianNL_rho; double mBoostedGaussianNT_rho; double mBoostedGaussianQuarkMass_rho; double mBoostedGaussianR2_phi; double mBoostedGaussianNL_phi; double mBoostedGaussianNT_phi; double mBoostedGaussianQuarkMass_phi; double mBoostedGaussianR2_jpsi; double mBoostedGaussianNL_jpsi; double mBoostedGaussianNT_jpsi; double mBoostedGaussianQuarkMass_jpsi; // hold the custom parameter (internal only) vector mCustomParameters; }; inline double DipoleModelParameters::BG() const {return mBG;} inline double DipoleModelParameters::mu02() const {return mMu02;} inline double DipoleModelParameters::C() const {return mC;} inline double DipoleModelParameters::lambdaG() const {return mLambdaG;} inline double DipoleModelParameters::Ag() const {return mAg;} inline double DipoleModelParameters::kappa() const {return mKappa;} inline double DipoleModelParameters::N0() const {return mN0;} inline double DipoleModelParameters::x0() const {return mX0;} inline double DipoleModelParameters::lambda() const {return mLambda;} inline double DipoleModelParameters::gammas() const {return mGammas;} inline double DipoleModelParameters::Bcgc() const {return mBcgc;} inline double DipoleModelParameters::quarkMass(unsigned int i) const { if (i < 6) return mQuarkMass[i]; else return 0; } inline double DipoleModelParameters::Rshrink() const {return mRshrink;} inline ProtonShape DipoleModelParameters::shape() const {return mShape;} inline double DipoleModelParameters::v() const {return mV;} #endif Index: branches/fitting/src/CMakeLists.txt =================================================================== --- branches/fitting/src/CMakeLists.txt (revision 379) +++ branches/fitting/src/CMakeLists.txt (revision 380) @@ -1,159 +1,161 @@ #=============================================================================== # CMakeLists.txt (src) # # Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich # # This file is part of Sartre. # # 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 . # # Author: Thomas Ullrich # Last update: # $Date$ # $Author$ #=============================================================================== include(ExternalProject) cmake_minimum_required (VERSION 3.1) # # Compiler flags for release and debug version # set(CMAKE_C_FLAGS_DEBUG " -g -W") set(CMAKE_CXX_FLAGS_DEBUG " -g -W -Wall -Wextra -Wno-potentially-evaluated-expression -pedantic -Wno-long-long -std=c++11") set(CMAKE_C_FLAGS_RELEASE " -O -W") set(CMAKE_CXX_FLAGS_RELEASE " -O -W -Wall -Wextra -Wno-potentially-evaluated-expression -pedantic -Wno-long-long -std=c++11") #set(CMAKE_C_FLAGS " -O -W") #set(CMAKE_CXX_FLAGS " -O -W -Wall -Wextra -pedantic -Wno-long-long") #set(CMAKE_CXX_STANDARD 11) #set(CMAKE_CXX_STANDARD_REQUIRED on) #set(CMAKE_CXX_EXTENSIONS OFF) # # Find external required packages # (see also FindGSL.cmake and FindROOT.cmke in cmake/modules) # Herer we need only the root library to create the table # tools. # # GSL find_package(GSL REQUIRED) include_directories(${GSL_INCLUDE_DIR}) # ROOT find_package(ROOT REQUIRED) include_directories(${ROOT_INCLUDE_DIR}) set(LIBS ${LIBS} ${ROOT_LIBRARIES}) #BOOST if (MULTITHREADED) set(Boost_USE_STATIC_LIBS ON) set(Boost_USE_MULTITHREADED ON) find_package(Boost 1.39 COMPONENTS thread REQUIRED) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIR}) endif(Boost_FOUND) endif (MULTITHREADED) # # Include files from sartre package # include_directories(${PROJECT_SOURCE_DIR}/src) include_directories(${PROJECT_SOURCE_DIR}/gemini) include_directories(${PROJECT_SOURCE_DIR}/cuba) # # Cuba is built using the autoconf shipped with it # ExternalProject_Add( cuba DOWNLOAD_COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/cuba ${PROJECT_BINARY_DIR}/cuba SOURCE_DIR ${PROJECT_BINARY_DIR}/cuba PREFIX ${PROJECT_BINARY_DIR}/cuba CONFIGURE_COMMAND ${PROJECT_BINARY_DIR}/cuba/configure --prefix=${PROJECT_BINARY_DIR}/cuba/build BUILD_COMMAND make lib BUILD_IN_SOURCE 1 ) # Cuba builds with an external configure/make # Copy the cuba directory from the source to the build directory so we can build it there add_custom_command( OUTPUT cuba COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/cuba ${PROJECT_BINARY_DIR}/cuba ) # # Defines source files for sartre library # set(SARTRE_SRC "AlphaStrong.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Amplitudes.cpp") set(SARTRE_SRC ${SARTRE_SRC} "BreakupProduct.cpp") set(SARTRE_SRC ${SARTRE_SRC} "CrossSection.cpp") set(SARTRE_SRC ${SARTRE_SRC} "DglapEvolution.cpp") set(SARTRE_SRC ${SARTRE_SRC} "DipoleModel.cpp") set(SARTRE_SRC ${SARTRE_SRC} "DipoleModelParameters.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Event.cpp") set(SARTRE_SRC ${SARTRE_SRC} "ExclusiveFinalStateGenerator.cpp") set(SARTRE_SRC ${SARTRE_SRC} "FinalStateGenerator.cpp") set(SARTRE_SRC ${SARTRE_SRC} "FrangibleNucleus.cpp") set(SARTRE_SRC ${SARTRE_SRC} "EventGeneratorSettings.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Integrals.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Kinematics.cpp") set(SARTRE_SRC ${SARTRE_SRC} "ModeFinderFunctor.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Nucleon.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Nucleus.cpp") set(SARTRE_SRC ${SARTRE_SRC} "PhotonFlux.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Sartre.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Settings.cpp") set(SARTRE_SRC ${SARTRE_SRC} "Table.cpp") set(SARTRE_SRC ${SARTRE_SRC} "TableCollection.cpp") set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorNucleus.cpp") set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorSettings.cpp") set(SARTRE_SRC ${SARTRE_SRC} "WaveOverlap.cpp") set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractiveCrossSections.cpp") set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractiveCrossSections_dz.cpp") +set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractiveCrossSections_DIS.cpp") + add_library(sartre ${SARTRE_SRC}) # # Table tools (all stand alone programs) # add_executable(tableInspector tableInspector.cpp) add_executable(tableMerger tableMerger.cpp) add_executable(tableQuery tableQuery.cpp) add_executable(tableDumper tableDumper.cpp) add_executable(tableVarianceMaker tableVarianceMaker.cpp) target_link_libraries(tableInspector ${CMAKE_CURRENT_BINARY_DIR}/libsartre.a ${LIBS}) target_link_libraries(tableMerger ${CMAKE_CURRENT_BINARY_DIR}/libsartre.a ${LIBS}) target_link_libraries(tableQuery ${CMAKE_CURRENT_BINARY_DIR}/libsartre.a ${LIBS}) target_link_libraries(tableDumper ${CMAKE_CURRENT_BINARY_DIR}/libsartre.a ${LIBS}) target_link_libraries(tableVarianceMaker ${CMAKE_CURRENT_BINARY_DIR}/libsartre.a ${LIBS}) add_dependencies(tableInspector sartre) add_dependencies(tableMerger sartre) add_dependencies(tableQuery sartre) add_dependencies(tableDumper sartre) add_dependencies(tableVarianceMaker sartre) # # Install library and include files (make install) within # the distribution tree. Top level CMakeLists.txt will install # Sartre in final destination. # install(TARGETS sartre DESTINATION sartre/lib) install(FILES "${PROJECT_BINARY_DIR}/cuba/build/lib/libcuba.a" DESTINATION sartre/lib) FILE(GLOB AllIncludeFiles *.h) install(FILES ${AllIncludeFiles} DESTINATION sartre/include) install(FILES "${PROJECT_BINARY_DIR}/cuba/build/include/cuba.h" DESTINATION sartre/include) install(TARGETS tableInspector DESTINATION sartre/bin) install(TARGETS tableMerger DESTINATION sartre/bin) install(TARGETS tableQuery DESTINATION sartre/bin) install(TARGETS tableDumper DESTINATION sartre/bin) install(TARGETS tableVarianceMaker DESTINATION sartre/bin) Index: branches/fitting/src/DipoleModel.h =================================================================== --- branches/fitting/src/DipoleModel.h (revision 379) +++ branches/fitting/src/DipoleModel.h (revision 380) @@ -1,110 +1,158 @@ //============================================================================== // DipoleModel.h // -// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich +// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich // -// This file is part of Sartre. +// This file is part of Sartre version: 1.1 // // 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 . // // Author: Tobias Toll // Last update: // $Date$ // $Author$ //============================================================================== #ifndef DipoleModel_h #define DipoleModel_h -#include "AlphaStrong.h" +#include "AlphaStrong.h" +#include "TH2D.h" #include "TableGeneratorNucleus.h" -#include "DipoleModelParameters.h" +#include "Parameters_Minuit2.h" +//#include "LHAPDF/LHAPDF.h" //LHAPDF -class TH2F; -class TH2D; +class TH2F; class TH1F; class DipoleModel { public: DipoleModel(); virtual ~DipoleModel(); const TableGeneratorNucleus* nucleus() const; bool configurationExists() const; virtual void createConfiguration(int)=0; virtual double dsigmadb2(double, double, double, double)=0; virtual double bDependence(double); virtual double bDependence(double, double); - virtual double dsigmadb2ep(double, double, double); - virtual double coherentDsigmadb2(double, double, double) {return 0;}; - virtual void createSigma_ep_LookupTable(double) {/* nothing */}; - - void createCompressedExponentThickness(); + virtual double dsigmadb2ep(double, double, double); +// virtual double calculateNorm(double); + virtual double ddxdsigmadb2(double, double, double); + virtual double coherentDsigmadb2(double, double, double){ return 0;}; + virtual void createSigma_ep_LookupTable(double){/* nothing */}; + virtual void createNormTable(){/* nothing*/}; + virtual void setMu02(double){/* nothing */}; + virtual void useSkewednessCorrection(){/* nothing */}; + virtual void setParams(const std::vector&); + TH1D* mLP_lu; + TH1D* mWS_lu; + TH2D* mCE_lu; + void createWoodsSaxonThickness(); + void createLaplaceThickness(); //BS + void createCompressedExponentThickness(); + + + double mBg=params_B_G; + double mC=params_C; + double mMu_02=params_mu02; + double mAg=params_A_g; + double mLamg=params_lambda_g; + double mRmax=params_r_max; + double mLmass=params_l_mass; + double mCmass=params_c_mass; + double mBmass=params_b_mass; + double qMass[6]; + double mAs_qqg; + double mRhs; + double mnLP=params_nLP; + double mR0WS=params_R0WS; + double mdWS=params_dWS; + double muCE=params_uCE; + //double mvCE=params_vCE; + double mEP=params_EP; + double BV=params_BV; + double RV=params_RV; + double wV=params_wV; + protected: + double uiWoodsSaxon(double*, double*); + double uiLaplace(double*, double*); //BS + double uiCompressedExponent(double*, double*); TableGeneratorNucleus mNucleus; - DipoleModelParameters *mParameters; - - AlphaStrong mAs; + AlphaStrong mAs; bool mConfigurationExists; bool mIsInitialized; + bool mUseSkewednessCorrection; - double uiCompressedExponent(double*, double*); - TH2D* mCE_lu; -}; + double protonThickness(double, double); + +}; class DipoleModel_bSat : public DipoleModel { public: DipoleModel_bSat(); - DipoleModel_bSat(const DipoleModel_bSat&); + ~DipoleModel_bSat(); + DipoleModel_bSat(const DipoleModel_bSat&); DipoleModel_bSat& operator=(const DipoleModel_bSat&); - ~DipoleModel_bSat(); - void createSigma_ep_LookupTable(double); + void createNormTable(); double dsigmadb2ep(double, double, double); - -protected: + void setMu02(double); + void useSkewednessCorrection(); + double coherentDsigmadb2(double, double, double); //BS MOVED ExclusiveCSMain Privacy + +protected: void createConfiguration(int); - double dsigmadb2(double, double, double, double); - double bDependence(double, double); + double dsigmadb2(double, double, double, double); + double bDependence(double, double); + double ddxdsigmadb2(double, double, double); protected: - TH2F* mBDependence; + TH2F* mBDependence; + // HERAPDF20_NNLO_ALPHAS_118 + // LHAPDF::PDF* mPDF;//=LHAPDF::mkPDF("HERAPDF20_NNLO_ALPHAS_118"); //LHAPDF private: + double mMu02; + double calculateNorm(double); + double dsigmadb2epForIntegration(double*, double*); - TH1F* mSigma_ep_LookupTable; - double coherentDsigmadb2(double, double, double); + TH1F* mSigma_ep_LookupTable; + TH1F* mThicknessNormTable; + }; class DipoleModel_bNonSat : public DipoleModel_bSat{ public: DipoleModel_bNonSat(); ~DipoleModel_bNonSat(); double dsigmadb2ep(double, double, double); private: double dsigmadb2(double, double, double, double); double coherentDsigmadb2(double, double, double); -}; -class DipoleModel_bCGC : public DipoleModel { -public: - DipoleModel_bCGC(); +protected: + // LHAPDF::PDF* mPDF; -private: - void createConfiguration(int); - double dsigmadb2(double, double, double, double); - double dsigmadb2ep(double, double, double); - double bDependence(double); }; +//class DipoleModel_bCGC : public DipoleModel { +//private: +//// void createConfiguration(int); +//// double dsigmadb2(double, double, double, double); +//// double dsigmadb2ep(double, double, double); +//// double bDependence(double); +//}; + #endif Index: branches/fitting/src/InclusiveDiffractiveCrossSections_DIS.cpp =================================================================== --- branches/fitting/src/InclusiveDiffractiveCrossSections_DIS.cpp (revision 0) +++ branches/fitting/src/InclusiveDiffractiveCrossSections_DIS.cpp (revision 380) @@ -0,0 +1,192 @@ +#include "Math/SpecFunc.h" +#include "TFile.h" +#include +#include +#include +#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" +#include "DipoleModelParameters.h" + +#define PR(x) cout << #x << " = " << (x) << endl; + +int integrandWrapperDIS( const int *, const double x[], const int *, double fval[], void* ptr2integrals) +{ + InclusiveDiffractiveCrossSections* ip=static_cast( ptr2integrals ); + + // double rlow=1e-4; //fm + // + // double rthigh=exp(-rlow); + // cout<<"#TT 1"<uiDIS(b, r, z, ip->mIsTransverse)*jacobian; + // cout<<"#TT 4"<uiDIS(b, r, z, ip->mIsTransverse)*jacobian; + double z = x[2]; + double mf=ip->qMass[ip->quarkIndex()]; + double eps2=z*(1-z)*ip->mQ2+mf*mf; + double r_max=8./sqrt(eps2); //fm (The factor is units fmGeV) + + double low[2]={0., 1e-4}; //b, r + double high[2]={2.2325002, r_max}; +// double high[2]={2.2325, 12.}; + // double high[3]={2.5, 3., 1.};//{3., 12., 1. + + double b = low[0] + (high[0]-low[0]) * x[0];//fm + double r = low[1] + (high[1]-low[1]) * x[1];//fm + + double jacobian = (high[0]-low[0])*(high[1]-low[1]); //fm2 + + fval[0]=ip->uiDIS(b, r, z)*jacobian; //fm2 + return 0; +} + +void InclusiveDiffractiveCrossSections::operator()(double Q2, double W2){ + mQ2=Q2; + mW2=W2; + + calculateTotal(); +} + +void InclusiveDiffractiveCrossSections::calculateTotal(){ + + const double epsrel=1.e-2, epsabs=1e-12; + const int flags=0, mineval=1e3, maxeval=1e8, key=0; + int nregionsT, nevalT, failT; + int nregionsL, nevalL, failL; + double valT=0, errT=0, probT=0; + double valL=0, errL=0, probL=0; + const char* statefile=0; + const int nvec=1; + + if(!mIsInitialized){ + cout<<"InclusiveDiffractiveCrossSections::dsigmadbeta Warning: InclusiveDiffractiveCrossSections not initialised!"<1) dipoleModel()->createSigma_ep_LookupTable(x); + + double resultT=0; + double resultL=0; + mTotalCST=0; + mTotalCSL=0; + for(unsigned int i=0; i<5; i++){ + if(!mWhichFlavours[i]) continue; + setQuarkIndex(i); + + double zlow=0; + double zhigh=1-zlow; + + // mTotalT->SetParameter(0, Q2); + // mTotalT->SetParameter(1, W2); + // resultT=mGITotalT.Integral(zlow, zhigh); //GeV^2*fm^4 + + // mTotalL->SetParameter(0, Q2); + // mTotalL->SetParameter(1, W2); + // resultL=mGITotalL.Integral(zlow, zhigh); //GeV^2*fm^4 + mIsTransverse=true; + Cuhre(3, 1, integrandWrapperDIS, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT); + resultT=valT; + + mIsTransverse=false; + Cuhre(3, 1, integrandWrapperDIS, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL); + resultL=valL; + + //Store the results: + mTotalCST+=resultT/hbarc2/hbarc2; //GeV^-2 + mTotalCSL+=resultL/hbarc2/hbarc2; //GeV^-2 + } //flavour +} + +double InclusiveDiffractiveCrossSections::uiDIS(double b, double r, double z){ + double Q2=mQ2; + double W2=mW2; + + double xbj=Kinematics::x(Q2, W2); + if(quarkIndex()>2) xbj=xbj*(1.+4.*pow(qMass[quarkIndex()], 2)/Q2); //as in RSVV + double waveOverlap=0.; + if(mIsTransverse) + waveOverlap=waveOverlapT(r, z, Q2); //GeV^2 + else + waveOverlap=waveOverlapL(r, z, Q2); //GeV^2 + double dsigmadb2=0; + if(mA==1) + dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xbj); + else + dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xbj); + + double result=waveOverlap*dsigmadb2/(4*M_PI); //GeV^2 + + //Do angular parts of r- and b-integrals: + result*=(2*M_PI)*b*(2*M_PI)*r; //GeV^0 + + return result; +} +double InclusiveDiffractiveCrossSections::waveOverlapT(double r, double z, double Q2){ + double ef=quarkCharge[quarkIndex()]; + double mf=qMass[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=qMass[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 +} Index: branches/fitting/src/InclusiveDiffractiveCrossSections.cpp =================================================================== --- branches/fitting/src/InclusiveDiffractiveCrossSections.cpp (revision 379) +++ branches/fitting/src/InclusiveDiffractiveCrossSections.cpp (revision 380) @@ -1,832 +1,857 @@ //The scales for alpha_s // // For GBW, the scale z(1-z)k^2 corresponds to the gg-dipole. We want the scale // for the qq-dipole only. The best choice may be to simply take Q2. // // For GBW(beta->0) there is no clear scale choice. // The large Q2 limit suggests small qq-dipole, // while the small beta limit suggest large qq-dipole // // Instead, one may use a different interpolation between the two schemes: // Norm*(Q2*beta*FT_GBW + (s-Q2)(1-beta)*FT_MS) // // Q2/(xpom*s)*FT_GBW + (1-Q2/(xpom*s))*FT_MS // // Q2=y*x*s=y*xpom*beta*s => xpom*s=Q2/beta/y // // Interpolation: beta*y*FT_GBW + (1-beta*y)*FT_MS // // GBW: 4*kappa^2+\mu_0^2 // MS: 4/r^2+\mu_0^2 // GBW(beta=0): 4*kappa^2+\mu_0^2 // THIS FAILS SPECTACULARLY!!! // --Because it assumes the wrong qq-dipole size (the gg-dipole size) // Could the correct interpolation be GBW(beta=0, alpha_s=alpha_s(MS))??? // // TODO: The number of parameters of dPhi0db, and dPhi1db can be reduced to one: z, // by using the kinematicPoint vector instead... // //Notes on Scale in Alpha_s: // In the small Q2 approximation, there is a small qqbar dipole, which splits into // a g-qqbar~gg dipole. Here it is reasonable to assume that the dipole size // felt by the guon being emitted is smaller than the subsequent tripole. // Here the scale should be set by r2~1/Q2, and the natural scale for this is // k2 which is integrated from 0 to Q2. // // For the small beta approximation, the dipole is assumed to be large to begin with // and split into two smaller dipoles vi a gluon emission "inbetween" // Here, the natural scale for the dipole is therefore~r, i.e. the size of the tripole // (the approximation of course being that the tripole is not expected to change size // in the time between the qq->qqg splitting and the interaction with the pomeron.) // // KLMV are using a constant alpha_s, which essentially is a free parameter // (within reasonable values) // They found a different value for as in bSat and bCGC, indicating that there // are other factors in diffraction affecting the normalisation. // (the qq-qqg split should be independent on dipole model) // Therefore, one can can expect a more consistent treatment of the strong coupling // to yield a *worse* description of the data. However, up to normalization, // the description should be better. // // However, this constant only absorbs normalisation for the QQG part of the cross-section. // Furhter more, with these definitions of alpha_s there are two different // scales for the two different approximations, changing that part of the calculation... // //============================================================================== // 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 . // // 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 "TFile.h" #include #include #include #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 "DipoleModelParameters.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::setParamsDCS(const std::vector& par) +{ + // qMass[6] = {par[0], par[0], par[0], par[1], par[2], params_t_mass}; + mBg = par[0]; //par[0];//BS For MZ, BG is kept fixed + mLmass = par[1]; + mCmass = par[2]; + mBmass = par[3]; + mC = par[4]; + mMu_02 = par[5]; //par[5];//BS For MZ, Mu02 is kept fixed + mAg = par[6]; + mLamg = par[7]; + mRmax = par[8]; + mAs_qqg = par[9]; + mRhs = par[10]; + mnLP = par[11]; + mR0WS = par[12]; + mdWS = par[13]; + muCE = par[14]; + mvCE = par[15]; + //mR_L2 = par[16]; + //mN_L = par[17]; + //mN_T = par[18]; + + qMass[0] = mLmass; + qMass[1] = mLmass; + qMass[2] = mLmass; + qMass[3] = mCmass; + qMass[4] = mBmass; + qMass[5] = params_t_mass; + double b_rms=sqrt(2*mBg)*hbarc; + mUpperBinIntegration=4.*b_rms;//20.*sqrt(params_BG)*hbarc; //fm + //PR(mUpperBinIntegration); + + //cout<< "InclusiveDiffractiveCrossSections mAg = " + std::to_string(mAg) << endl; +} void InclusiveDiffractiveCrossSections::init() { cout<<"init InclusiveDiffractiveCrossSections"<dipoleModelType(); - 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); - } - + + DipoleModelType model=settings->dipoleModelType(); + if (model==bSat) { + mDipoleModel = new DipoleModel_bSat; + } + else if(model==bNonSat){ + mDipoleModel = new DipoleModel_bNonSat; + } + else { + cout << "Integrals::init(): Error, model not implemented: "<< model << endl; + exit(1); + } + //Find upper integration limit in b and r: mUpperRinIntegrationDIS=8.;//4*hbarc/quarkMass[0];//fm //8.; mUpperRinIntegrationQQ=20.;//9.*bSat_rmax;//bSat_rmax; //fm mUpperRinIntegrationQQG=13.;//9.*bSat_rmax;//bSat_rmax; //fm if(dipoleModel()->nucleus()->A() > 1) - mUpperBinIntegration=mUpperRinIntegrationDIS; + mUpperBinIntegration=mUpperRinIntegrationDIS; else{ - switch (mParameters->shape()){ - case HardSphere: - cout<<"Proton is a Hard Sphere."<BG()*hbarc; //fm - break; - case Laplace: - cout<<"Proton is a Laplace-distribution."<BG())*hbarc;//fm - cout<<"b_rms="<shape()){ + case HardSphere: + cout<<"Proton is a Hard Sphere."<BG()*hbarc; //fm + break; + case Laplace: + cout<<"Proton is a Laplace-distribution."<qqg splitting + // Setup the strong coupling for the qq->qqg splitting // (should be the same as in DipoleModel, if not default): // // mAs.init(0.1183, 1); - + mA=settings->A(); - + // // Set up all integration functions: // // //Inclusive Difrfaction Integrals: // double precisionQQ=1e-3; double precisionQQG=1e-2; - double precisionDIS=1e-5; + double precisionDIS=1e-5; setRelativePrecisionOfIntegration(precisionQQG); - - mIntegralT=new TF1("IntegralT", this, &InclusiveDiffractiveCrossSections::uiIntegralT, 0., 0.5, 3); + + mIntegralT=new TF1("IntegralT", this, &InclusiveDiffractiveCrossSections::uiIntegralT, 0., 0.5, 3); mWFT=new ROOT::Math::WrappedTF1(*mIntegralT); mGIintegralT.SetFunction(*mWFT); mGIintegralT.SetRelTolerance(precisionQQ); - mIntegralL=new TF1("IntegralL", this, &InclusiveDiffractiveCrossSections::uiIntegralL, 0., 0.5, 3); + mIntegralL=new TF1("IntegralL", this, &InclusiveDiffractiveCrossSections::uiIntegralL, 0., 0.5, 3); mWFL=new ROOT::Math::WrappedTF1(*mIntegralL); mGIintegralL.SetFunction(*mWFL); mGIintegralL.SetRelTolerance(precisionQQ); - + mPhi0=new TF1("Phi0", this, &InclusiveDiffractiveCrossSections::dPhi0db, 0, mUpperBinIntegration, 4); mWFPhi0=new ROOT::Math::WrappedTF1(*mPhi0); mGIPhi0.SetFunction(*mWFPhi0); // mGIPhi0.SetRelTolerance(precisionQQ/2.); mGIPhi0.SetRelTolerance(precisionQQ/100.); - + mPhi1=new TF1("Phi1", this, &InclusiveDiffractiveCrossSections::dPhi1db, 0, mUpperBinIntegration, 4); mWFPhi1=new ROOT::Math::WrappedTF1(*mPhi1); mGIPhi1.SetFunction(*mWFPhi1); // mGIPhi1.SetRelTolerance(precisionQQ/2.); mGIPhi1.SetRelTolerance(precisionQQ/100.); mPhi0r=new TF1("Phi0", this, &InclusiveDiffractiveCrossSections::uiPhi0, 0, 100., 5); mWFPhi0r=new ROOT::Math::WrappedTF1(*mPhi0r); mGIPhi0r.SetFunction(*mWFPhi0r); // mGIPhi0r.SetRelTolerance(precisionQQ/4.); mGIPhi0r.SetRelTolerance(precisionQQ/200.); - + mPhi1r=new TF1("Phi1", this, &InclusiveDiffractiveCrossSections::uiPhi1, 0, 100., 5); mWFPhi1r=new ROOT::Math::WrappedTF1(*mPhi1r); mGIPhi1r.SetFunction(*mWFPhi1r); // mGIPhi1r.SetRelTolerance(precisionQQ/4.); mGIPhi1r.SetRelTolerance(precisionQQ/200.); - + mRGBW=new TF1("frIntegral", this, &InclusiveDiffractiveCrossSections::uiR_GBW, 0, 100., 4); mWFRGBW=new ROOT::Math::WrappedTF1(*mRGBW); mGIRGBW.SetFunction(*mWFRGBW); mGIRGBW.SetRelTolerance(precisionQQG); mRGBWbeta0=new TF1("frIntegral", this, &InclusiveDiffractiveCrossSections::uiR_GBWbeta0, 0, 100, 3); mWFRGBWbeta0=new ROOT::Math::WrappedTF1(*mRGBWbeta0); mGIRGBWbeta0.SetFunction(*mWFRGBWbeta0); mGIRGBWbeta0.SetRelTolerance(precisionQQG); - + // // Inclusive DIS integrals: // //#TT CHECK PRECISION HERE! - mTotalT=new TF1("IntegralT", this, &InclusiveDiffractiveCrossSections::uiDsigmadzT, 0., 1., 2); - mWFTotalT=new ROOT::Math::WrappedTF1(*mTotalT); - mGITotalT.SetFunction(*mWFTotalT); - // mGITotalT.SetRelTolerance(precisionDIS); - mGITotalT.SetRelTolerance(1e-2); - - mTotalL=new TF1("IntegralL", this, &InclusiveDiffractiveCrossSections::uiDsigmadzL, 0., 1., 2); - mWFTotalL=new ROOT::Math::WrappedTF1(*mTotalL); - mGITotalL.SetFunction(*mWFTotalL); - //mGITotalL.SetRelTolerance(precisionDIS); - mGITotalL.SetRelTolerance(1e-2); - - mSigmaT=new TF1("sigmaT", this, &InclusiveDiffractiveCrossSections::uiDsigmadrT, 0., 100., 3); - mWFSigmaT=new ROOT::Math::WrappedTF1(*mSigmaT); - mGISigmaT.SetFunction(*mWFSigmaT); - mGISigmaT.SetRelTolerance(precisionDIS/2.); - - mSigmaL=new TF1("sigmaL", this, &InclusiveDiffractiveCrossSections::uiDsigmadrL, 0., 100., 3); - mWFSigmaL=new ROOT::Math::WrappedTF1(*mSigmaL); - mGISigmaL.SetFunction(*mWFSigmaL); - mGISigmaL.SetRelTolerance(precisionDIS/2.); - - mDipole= new TF1("dp", this, &InclusiveDiffractiveCrossSections::uiDipoleSigma, 0., 100., 2); - mWFDipole=new ROOT::Math::WrappedTF1(*mDipole); - mGIDipole.SetFunction(*mWFDipole); - mGIDipole.SetRelTolerance(precisionDIS/4.); +// mTotalT=new TF1("IntegralT", this, &InclusiveDiffractiveCrossSections::uiDsigmadzT, 0., 1., 2); +// mWFTotalT=new ROOT::Math::WrappedTF1(*mTotalT); +// mGITotalT.SetFunction(*mWFTotalT); +// // mGITotalT.SetRelTolerance(precisionDIS); +// mGITotalT.SetRelTolerance(1e-2); +// +// mTotalL=new TF1("IntegralL", this, &InclusiveDiffractiveCrossSections::uiDsigmadzL, 0., 1., 2); +// mWFTotalL=new ROOT::Math::WrappedTF1(*mTotalL); +// mGITotalL.SetFunction(*mWFTotalL); +// //mGITotalL.SetRelTolerance(precisionDIS); +// mGITotalL.SetRelTolerance(1e-2); +// +// mSigmaT=new TF1("sigmaT", this, &InclusiveDiffractiveCrossSections::uiDsigmadrT, 0., 100., 3); +// mWFSigmaT=new ROOT::Math::WrappedTF1(*mSigmaT); +// mGISigmaT.SetFunction(*mWFSigmaT); +// mGISigmaT.SetRelTolerance(precisionDIS/2.); +// +// mSigmaL=new TF1("sigmaL", this, &InclusiveDiffractiveCrossSections::uiDsigmadrL, 0., 100., 3); +// mWFSigmaL=new ROOT::Math::WrappedTF1(*mSigmaL); +// mGISigmaL.SetFunction(*mWFSigmaL); +// mGISigmaL.SetRelTolerance(precisionDIS/2.); +// +// mDipole= new TF1("dp", this, &InclusiveDiffractiveCrossSections::uiDipoleSigma, 0., 100., 2); +// mWFDipole=new ROOT::Math::WrappedTF1(*mDipole); +// mGIDipole.SetFunction(*mWFDipole); +// mGIDipole.SetRelTolerance(precisionDIS/4.); + + mIsInitialized=true; - 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); -} -*/ + 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=mParameters->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); + 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); + 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; + if (xpom<0 || xpom>1) + result = false; + kinematicPoint[3]=xpom; kinematicPoint[4]=beta; - return result; + return result; } void InclusiveDiffractiveCrossSections::setNumberOfFlavours(unsigned int val){ mNumberOfFlavours=val; } void InclusiveDiffractiveCrossSections::setWhichFlavours(bool val0, bool val1, bool val2, bool val3, bool val4){ - mWhichFlavours[0]=val0;//u - mWhichFlavours[1]=val1;//d - mWhichFlavours[2]=val2;//s - mWhichFlavours[3]=val3;//c - mWhichFlavours[4]=val4;//b + mWhichFlavours[0]=val0;//u + mWhichFlavours[1]=val1;//d + mWhichFlavours[2]=val2;//s + mWhichFlavours[3]=val3;//c + mWhichFlavours[4]=val4;//b } /////////////Cross-Sections////////////////////////////////////////////////////////// void InclusiveDiffractiveCrossSections::calculate_dsigmadbeta_QQ(double Q2, double W2, double beta){ if(beta1) dipoleModel()->createSigma_ep_LookupTable(xpom); - + double resultT=0; double resultL=0; for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) - integrate=false; - } - if(integrate){ - double zlow=zlower(); - double zhigh=0.5; - mIntegralT->SetParameter(0, beta); - mIntegralT->SetParameter(1, xpom); - mIntegralT->SetParameter(2, Q2); - resultT=mGIintegralT.Integral(zlow, zhigh); //GeV^-4 - - mIntegralL->SetParameter(0, beta); - mIntegralL->SetParameter(1, xpom); - mIntegralL->SetParameter(2, Q2); - resultL=mGIintegralL.Integral(zlow, zhigh); //GeV^-6 - } - // - // 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 - //Store the results: - mInclusiveT[i]=resultT*ef*ef * Nc*Q2*Q2 / (16*pow(M_PI, 3)*beta) * xF_to_dsigmadbeta; - mInclusiveL[i]=resultL*ef*ef * Nc*Q2*Q2*Q2 / (4*pow(M_PI, 3)*beta) * xF_to_dsigmadbeta; + if(!mWhichFlavours[i]) continue; + 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) + integrate=false; + } + if(integrate){ + double zlow=zlower(); + double zhigh=0.5; + mIntegralT->SetParameter(0, beta); + mIntegralT->SetParameter(1, xpom); + mIntegralT->SetParameter(2, Q2); + resultT=mGIintegralT.Integral(zlow, zhigh); //GeV^-4 + + mIntegralL->SetParameter(0, beta); + mIntegralL->SetParameter(1, xpom); + mIntegralL->SetParameter(2, Q2); + resultL=mGIintegralL.Integral(zlow, zhigh); //GeV^-6 + } + // + // 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 + //Store the results: + mInclusiveT[i]=resultT*ef*ef * Nc*Q2*Q2 / (16*pow(M_PI, 3)*beta) * xF_to_dsigmadbeta; + mInclusiveL[i]=resultL*ef*ef * Nc*Q2*Q2*Q2 / (4*pow(M_PI, 3)*beta) * xF_to_dsigmadbeta; } } 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; - - mPhi0->SetParameter(0, beta); - mPhi0->SetParameter(1, xpom); - mPhi0->SetParameter(2, z); - mPhi0->SetParameter(3, Q2); - double Phi0=mGIPhi0.IntegralUp(blow); //fm^6 - Phi0 = Phi0 /(hbarc2*hbarc2*hbarc2); //GeV^-6 - - return z*z*z*(1-z)*(1-z)*(1-z)*Phi0; //GeV-6 + double beta=par[0]; + double xpom=par[1]; + double Q2=par[2]; + double z=*x; + double blow=0; + + mPhi0->SetParameter(0, beta); + mPhi0->SetParameter(1, xpom); + mPhi0->SetParameter(2, z); + mPhi0->SetParameter(3, Q2); + double Phi0=mGIPhi0.IntegralUp(blow); //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; - - mPhi0->SetParameter(0, beta); - mPhi0->SetParameter(1, xpom); - mPhi0->SetParameter(2, z); - mPhi0->SetParameter(3, Q2); - double Phi0=mGIPhi0.IntegralUp(blow); //fm^6 - Phi0 = Phi0/hbarc2/hbarc2/hbarc2; //GeV^-6 - - mPhi1->SetParameter(0, beta); - mPhi1->SetParameter(1, xpom); - mPhi1->SetParameter(2, z); - mPhi1->SetParameter(3, Q2); - double Phi1=mGIPhi1.IntegralUp(blow); //fm^6 - Phi1 = Phi1/hbarc2/hbarc2/hbarc2; //GeV^-6 - - double mf=mParameters->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 beta=par[0]; + double xpom=par[1]; + double Q2=par[2]; + double z=*x; + double blow=0; + + mPhi0->SetParameter(0, beta); + mPhi0->SetParameter(1, xpom); + mPhi0->SetParameter(2, z); + mPhi0->SetParameter(3, Q2); + double Phi0=mGIPhi0.IntegralUp(blow); //fm^6 + Phi0 = Phi0/hbarc2/hbarc2/hbarc2; //GeV^-6 + + mPhi1->SetParameter(0, beta); + mPhi1->SetParameter(1, xpom); + mPhi1->SetParameter(2, z); + mPhi1->SetParameter(3, Q2); + double Phi1=mGIPhi1.IntegralUp(blow); //fm^6 + Phi1 = Phi1/hbarc2/hbarc2/hbarc2; //GeV^-6 + + double mf=mParameters->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-4; - - mPhi0r->SetParameter(0, beta); - mPhi0r->SetParameter(1, xpom); - mPhi0r->SetParameter(2, z); - mPhi0r->SetParameter(3, Q2); - mPhi0r->SetParameter(4, b); - double integral=mGIPhi0r.IntegralUp(rlow); //fm^2 - - double result = integral*integral*2*M_PI*b; //fm^5 - return result; + double beta=par[0]; + double xpom=par[1]; + double z=par[2]; + double Q2=par[3]; + double b=*x; + double rlow=1e-4; + + mPhi0r->SetParameter(0, beta); + mPhi0r->SetParameter(1, xpom); + mPhi0r->SetParameter(2, z); + mPhi0r->SetParameter(3, Q2); + mPhi0r->SetParameter(4, b); + double integral=mGIPhi0r.IntegralUp(rlow); //fm^2 + + double result = integral*integral*2*M_PI*b; //fm^5 + return result; } 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-4;; - - mPhi1r->SetParameter(0, beta); - mPhi1r->SetParameter(1, xpom); - mPhi1r->SetParameter(2, z); - mPhi1r->SetParameter(3, Q2); - mPhi1r->SetParameter(4, b); - double integral=mGIPhi1r.IntegralUp(rlow); //fm^2 - double result = integral*integral*2*M_PI*b; //fm^5 - return result; + double beta=par[0]; + double xpom=par[1]; + double z=par[2]; + double Q2=par[3]; + double b=*x; + double rlow=1e-4;; + + mPhi1r->SetParameter(0, beta); + mPhi1r->SetParameter(1, xpom); + mPhi1r->SetParameter(2, z); + mPhi1r->SetParameter(3, Q2); + mPhi1r->SetParameter(4, b); + double integral=mGIPhi1r.IntegralUp(rlow); //fm^2 + double result = integral*integral*2*M_PI*b; //fm^5 + return result; } double InclusiveDiffractiveCrossSections::uiPhi0(double* x, double* par){ double r=*x; double beta=par[0]; double xpom=par[1]; - double z=par[2]; + double z=par[2]; double Q2=par[3]; double b=par[4]; double mf=mParameters->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? + 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); //GeV double besselK0=TMath::BesselK0(epsilon*r/hbarc); double besselJ0=TMath::BesselJ0(kappa*r/hbarc); double dsigmadb2=0; if(mA==1) - dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom); + dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom); else - dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom); + 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 z=par[2]; double Q2=par[3]; double b=par[4]; double mf=mParameters->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? + 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 besselJ1=TMath::BesselJ1(kappa*r/hbarc); double dsigmadb2=0; if(mA==1) - dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom); + dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom); else - dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom); + dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom); return r*besselK1*besselJ1*dsigmadb2; //fm } /////////////////////////////////////////////////////////////////////////////////////////// ////QQG//////////////////////////////////////////////////////////////////////////////////// void InclusiveDiffractiveCrossSections::putAllQQGtoZero(){ - for(int i=0; i<6; i++){ + 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!"<q+qbar+g splitting // If this coupling is 0, there is no QQG and we the QQG part remains 0. // double alpha_s=as_qqg_fixed; if(alpha_s == 0) return; - + // - // dsigma/dbeta = 4*M_PI*M_PI*alpha_em/(Q2*beta) * xpom * F + // 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 + 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!"<3){ - smallBetaExtrapolation=1; - } - mInclusiveQQG[i]=alpha_s*mInclusiveQQG_GBW[i]*smallBetaExtrapolation; - mInclusiveQQG[i]*=xF_to_dsigmadbeta; + if(!mWhichFlavours[i]) continue; + double smallBetaExtrapolation=mInclusiveQQG_MS[i]/mInclusiveQQG_GBWbeta0[i]; + if(i>3){ + smallBetaExtrapolation=1; + } + mInclusiveQQG[i]=alpha_s*mInclusiveQQG_GBW[i]*smallBetaExtrapolation; + mInclusiveQQG[i]*=xF_to_dsigmadbeta; } } ////QQG_GBW//////////////////////////////////////////////////////////////////////////////////// void InclusiveDiffractiveCrossSections::FTQQG_GBW(double Q2, double beta){ mBeta=beta; - double xpom=kinematicPoint[3]; + double xpom=kinematicPoint[3]; // Make the coherent lookup-table if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom); double result=0; double bhigh=mUpperBinIntegration; - + for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) - integrate=false; - } - if(integrate){ - double a[3]={0., 0., beta}; //b, k2, z - double b[3]={bhigh, 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 - } - double ef=quarkCharge[quarkIndex()]; - //Store the result: - mInclusiveQQG_GBW[i]=result*beta/(8*pow(M_PI, 4))*ef*ef; //GeV^0 + if(!mWhichFlavours[i]) continue; + 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) + integrate=false; + } + if(integrate){ + double a[3]={0., 0., beta}; //b, k2, z + double b[3]={bhigh, 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 + } + double ef=quarkCharge[quarkIndex()]; + //Store the result: + mInclusiveQQG_GBW[i]=result*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 xpom=kinematicPoint[3]; - + // double rhigh=min(1./sqrt(z*k2), 1./(sqrt((1-z)*k2)))*hbarc*mUpperRinIntegration; double rhigh=mUpperRinIntegrationQQG; double rlow=0.;; // double rhigh=mUpperRinIntegration; - + mRGBW->SetParameter(0, z); mRGBW->SetParameter(1, b); mRGBW->SetParameter(2, k2); mRGBW->SetParameter(3, xpom); - + double rIntegral=mGIRGBW.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 - + 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=par[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); + dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom); else - dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom); - + 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]; + double xpom=kinematicPoint[3]; // Make the coherent lookup-table if(mA>1) dipoleModel()->createSigma_ep_LookupTable(xpom); double result=0; double bhigh=mUpperBinIntegration; - + for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) - integrate=false; - } - if(integrate){ - double a[3]={0., 0.}; //b, k2 - double b[3]={bhigh, 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 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/(2*pow(M_PI, 4))*ef*ef*mysteryFactor; //GeV^0 - + if(!mWhichFlavours[i]) continue; + 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(mParameters->quarkMass(i-1)==mParameters->quarkMass(i)) + integrate=false; + } + if(integrate){ + double a[3]={0., 0.}; //b, k2 + double b[3]={bhigh, 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 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/(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 xpom=kinematicPoint[3]; - + double rlow=0.; // double rhigh=1./sqrt(k2)*hbarc*mUpperRinIntegration; double rhigh=mUpperRinIntegrationQQG; - + mRGBWbeta0->SetParameter(0, b); mRGBWbeta0->SetParameter(1, k2); mRGBWbeta0->SetParameter(2, xpom); - + double rIntegral=mGIRGBWbeta0.Integral(rlow, rhigh)/hbarc; //GeV^0 double result=2*M_PI*b/hbarc* //GeV^-1 - log(Q2/k2)*//GeV^0 - rIntegral*rIntegral; //GeV^0 - - + 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=par[2]; - + double k=sqrt(k2); - + double besselJ2=0; if(r!=0) - besselJ2=ROOT::Math::cyl_bessel_j(2, k*r/hbarc)/(r/hbarc); - + besselJ2=ROOT::Math::cyl_bessel_j(2, k*r/hbarc)/(r/hbarc); + double result=dsigmadb2tilde(r, b, xpom)*besselJ2; - + return result; //GeV^1 } /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// ////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 rhigh=mUpperRinIntegrationQQG; double bhigh=mUpperBinIntegration; for(unsigned int i=0; i<5; i++){ - if(!mWhichFlavours[i]) continue; - setQuarkIndex(i); - double a[5]={0., 0, 0, 1e-3, 0.}; //r, z, b, rp, phi - double b[5]={rhigh, 1, bhigh, rhigh, 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); - result=ig.Integral(a, b)/hbarc/hbarc/hbarc; //GeV*GeV^-3=GeV^-2 - - //Store the result: - double Cf=4/3; - mInclusiveQQG_MS[i]=result*Cf*Q2/(4*pow(M_PI, 4)*alpha_em); //GeV^0 + if(!mWhichFlavours[i]) continue; + setQuarkIndex(i); + double a[5]={0., 0, 0, 1e-3, 0.}; //r, z, b, rp, phi + double b[5]={rhigh, 1, bhigh, rhigh, 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); + result=ig.Integral(a, b)/hbarc/hbarc/hbarc; //GeV*GeV^-3=GeV^-2 + + //Store the result: + double Cf=4/3; + mInclusiveQQG_MS[i]=result*Cf*Q2/(4*pow(M_PI, 4)*alpha_em); //GeV^0 } } double InclusiveDiffractiveCrossSections::uiQQG_MS(const double* arg){ - //KLMV eq2. (12) and (13) - // From r and b integrals we get factors 4*pi^4*r*b - // From rp integral we get factor rp - + //KLMV eq2. (12) and (13) + // From r and b integrals we get factors 4*pi^4*r*b + // From rp integral we get factor rp + 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 rMinusrp=sqrt(r*r+rp*rp-2*r*rp*cos(phi)); double rprefactor=r*r/(rp*rp*rMinusrp*rMinusrp)*hbarc2; //GeV2 double Nfactor=QQG_MS_Nfactors(r, rp, rMinusrp, xpom, b); double A=rprefactor*Nfactor*Nfactor; //GeV2 - + //WaveOverlap (different prefactor in MS and KMW): double waveOverlap=waveOverlapT(r, z, Q2)/2.; //GeV^2 - + double result= waveOverlap*A //GeV^4 - *2*M_PI*b/hbarc //GeV^-1 - *2*M_PI*r/hbarc //GeV^-1 - *rp; // GeV^-1 = GeV - + *2*M_PI*b/hbarc //GeV^-1 + *2*M_PI*r/hbarc //GeV^-1 + *rp; // GeV^-1 = GeV + return result; //GeV } double InclusiveDiffractiveCrossSections::QQG_MS_Nfactors(double r, double rp, double rMinusrp, double xpom, double b){ double N_of_r=0; if(mA==1) - N_of_r=mDipoleModel->dsigmadb2ep(r, b, xpom)/2.; + N_of_r=mDipoleModel->dsigmadb2ep(r, b, xpom)/2.; else - N_of_r=mDipoleModel->coherentDsigmadb2(r, b, xpom)/2.; - + 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.; + N_of_rp=mDipoleModel->dsigmadb2ep(rp, b, xpom)/2.; else - N_of_rp=mDipoleModel->coherentDsigmadb2(rp, b, xpom)/2.; - + N_of_rp=mDipoleModel->coherentDsigmadb2(rp, b, xpom)/2.; + double N_of_rMinusrp=0; if(mA==1) - N_of_rMinusrp=mDipoleModel->dsigmadb2ep(rMinusrp, b, xpom)/2.; + N_of_rMinusrp=mDipoleModel->dsigmadb2ep(rMinusrp, b, xpom)/2.; else - N_of_rMinusrp=mDipoleModel->coherentDsigmadb2(rMinusrp, b, xpom)/2.; - + N_of_rMinusrp=mDipoleModel->coherentDsigmadb2(rMinusrp, b, xpom)/2.; + double Nfactor=N_of_rp+N_of_rMinusrp-N_of_r-N_of_rp*N_of_rMinusrp; //GeV^0 - + return Nfactor; } /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// Index: branches/fitting/src/DglapEvolution.cpp =================================================================== --- branches/fitting/src/DglapEvolution.cpp (revision 379) +++ branches/fitting/src/DglapEvolution.cpp (revision 380) @@ -1,519 +1,544 @@ -//============================================================================== +//============================================================================== // DglapEvolution.h -// +// // Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich -// +// // This file is part of Sartre. -// +// // 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 . -// +// // Author: Thomas Ullrich // $Date$ // $Author$ -//============================================================================== -#include "DglapEvolution.h" -#include "TableGeneratorSettings.h" -#include "DipoleModelParameters.h" -#include -#include -#include - -using namespace std; - -#define PR(x) cout << #x << " = " << (x) << endl; - -DglapEvolution* DglapEvolution::mInstance = 0; - -DglapEvolution::DglapEvolution() -{ - DipoleModelParameters parameters(TableGeneratorSettings::instance()); - - // +//============================================================================== +#include "DglapEvolution.h" +#include "TableGeneratorSettings.h" +#include "DipoleModelParameters.h" +#include +#include +#include +#include "Parameters_Minuit2.h" + + +using namespace std; + +#define PR(x) cout << #x << " = " << (x) << endl; + +DglapEvolution* DglapEvolution::mInstance = 0; + +//double DglapEvolution::m_Bg = params_B_G; +//double DglapEvolution::m_C = params_C; +//double DglapEvolution::mMu02 = params_mu02; +//double DglapEvolution::mAg = params_A_g; +//double DglapEvolution::mLamdaG = params_lambda_g; +//double DglapEvolution::mMC = params_c_mass; +//double DglapEvolution::mMB = params_b_mass; + +DglapEvolution::DglapEvolution() +{ + DipoleModelParameters parameters(TableGeneratorSettings::instance()); + + // // For speed purposes the key parameters are held as data member. // Here we assign the proper ones depending on the model and the // parameters set choosen. - // - mAg = parameters.Ag(); - mLambdaG = parameters.lambdaG(); - mMu02 = parameters.mu02(); - - // + // +// mAg = parameters.Ag(); +// mLambdaG = parameters.lambdaG(); +// mMu02 = parameters.mu02(); + + mMu02 = params_mu02; + mAg = params_A_g; + mLambdaG = params_lambda_g; + mMC = params_c_mass; + mMB = params_b_mass; + + PR(mMu02); + PR(mAg); + PR(mLambdaG); + PR(mMC); + PR(mMB); + // // Define alpha_s functor - // - mAlphaStrong = new AlphaStrong; - - mAlphaStrong->setMassCharm(parameters.quarkMass(3)); - mAlphaStrong->setMassBottom(parameters.quarkMass(4)); - - // - // Set alpha_s at c, b, t mass - // - mMC = mAlphaStrong->massCharm(); - mMB = mAlphaStrong->massBottom(); - mMT = mAlphaStrong->massTop(); + // + mAlphaStrong = new AlphaStrong; + + mAlphaStrong->setMassCharm(mMC); + mAlphaStrong->setMassBottom(mMB); + + // + // Set alpha_s at c, b, t mass + // +// mMC = mAlphaStrong->massCharm(); +// mMB = mAlphaStrong->massBottom(); + mMT = mAlphaStrong->massTop(); const double FourPi = 4.*M_PI; mALPC = mAlphaStrong->at(mMC*mMC)/FourPi; mALPB = mAlphaStrong->at(mMB*mMB)/FourPi; mALPT = mAlphaStrong->at(mMT*mMT)/FourPi; mALPS = mAlphaStrong->at(mMu02)/FourPi; // // Initialization of support points in n-space and weights for the // Gauss quadrature and of the anomalous dimensions for the RG // evolution at these n-values. // // // Weights and support points for nomalized 8 point gauss quadrature // double wz[9] = {0, 0.101228536290376,0.222381034453374,0.313706645877887, 0.362683783378362,0.362683783378362,0.313706645877887, 0.222381034453374,0.101228536290376}; double zs[9] = {0, -0.960289856497536,-0.796666477413627,-0.525532409916329, -0.183434642495650,0.183434642495650,0.525532409916329, 0.796666477413627,0.960289856497536}; // // Integration contour parameters // double down[18] = {0, 0., 0.5, 1., 2., 3., 4., 6., 8., 1.e1, 1.2e1, 1.5e1, 1.8e1, 2.1e1, 2.4e1, 2.7e1, 3.e1, 3.3e1}; double up[18]; mC = 0.8; double phi = M_PI * 3./4.; double co = cos(phi); double si = sin(phi); mCC = complex(co, si); for (int i=1; i <=16; i++) up[i] = down[i+1]; up[17] = 36.; // // Support points and weights for the gauss integration // (the factor (up-down)/2 is included in the weights) // int k = 0; double sum, diff, z; for (int i=1; i <=17; i++) { sum = up[i] + down[i]; diff = up[i] - down[i]; for (int j=1; j <=8; j++) { k++; z = 0.5 * (sum + diff * zs[j]); mWN[k] = diff / 2.* wz[j]; mN[k] = complex(mC+co*z+1.,si*z); } } - anom(); - - // - // Defaults for lookup table - // - mTableMinX = 1.e-10; - mTableMaxX = 1; - mTableMinQ2 = 1.; - mTableMaxQ2 = 1e6; - mLookupTableIsFilled = false; - mUseLookupTable = false; - mNumberOfNodesInX = mNumberOfNodesInQ2 = 0; + anom(); + + // + // Defaults for lookup table + // + mTableMinX = 1.e-10; + mTableMaxX = 1; + mTableMinQ2 = 1.; + mTableMaxQ2 = 1e7; + mLookupTableIsFilled = false; + mUseLookupTable = false; + mNumberOfNodesInX = mNumberOfNodesInQ2 = 0; } DglapEvolution::~DglapEvolution() { - if (mAlphaStrong) delete mAlphaStrong; - if (mLookupTableIsFilled) { - for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i) - delete [] mLookupTable[i]; - delete [] mLookupTable; - } + if (mAlphaStrong) delete mAlphaStrong; + if (mLookupTableIsFilled) { + for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i) + delete [] mLookupTable[i]; + delete [] mLookupTable; + } } -void DglapEvolution::reinit(){ - DglapEvolution(); -} +void DglapEvolution::reinit(){ + DglapEvolution(); +} double DglapEvolution::alphaSxG(double x, double Q2) { double val = xG(x, Q2); return val*mAlphaStrong->at(Q2); } double DglapEvolution::xG(double x, double Q2) -{ - double result = 0; - - if (!mUseLookupTable) { - result = xG_Engine(x, Q2); - } - else if (mUseLookupTable && !mLookupTableIsFilled) { - cout << "DglapEvolution::xG(): Warning, use of lookup table requested\n" - << " but table is not setup. First use \n" - << " generateLookupTable() to fill table.\n" - << " Will fall back to numeric calculation." - << endl; - result = xG_Engine(x, Q2); - } - else - result = xG_Interpolator(x, Q2); +{ + double result = 0; + + if (!mUseLookupTable) { + result = xG_Engine(x, Q2); + } + else if (mUseLookupTable && !mLookupTableIsFilled) { + cout << "DglapEvolution::xG(): Warning, use of lookup table requested\n" + << " but table is not setup. First use \n" + << " generateLookupTable() to fill table.\n" + << " Will fall back to numeric calculation." + << endl; + result = xG_Engine(x, Q2); + } + else + result = xG_Interpolator(x, Q2); return result; } double DglapEvolution::xG_Engine(double x, double Q2) { // // These are provided by AlphaStrong ensuring // consistency between evolution and alpha_s. // if (mAlphaStrong->order() != 0) { cout << "DglapEvolution::xG_Engine(): Fatal error, alpha_s is in order " << mAlphaStrong->order() << " but DglapEvolution only support order = 0. Stop here." << endl; exit(1); } double alpq = mAlphaStrong->at(Q2)/(4*M_PI); // // Q2 and x dependent quantities // double ax = log(x); double ex = exp(-mC * ax); // // integration length parameter for the mellin inversion // const int nmax = 136; // // Gluon density and output // complex fn[137]; reno(fn, alpq, nmax, mAg, mLambdaG); long double fun = 0; // long double is needed long double fz; complex xnm,cex; for (int i=1; i <= nmax; i++) { xnm = (mC - mN[i]+1.) * ax; cex = exp(xnm) / M_PI * mCC; fz = imag(fn[i]*cex); fun = fun + mWN[i] * fz; } double pa = fun * ex; return pa; } -void DglapEvolution::generateLookupTable(int nx, int nq2) -{ - // - // Delete old lookup table (if it exists) - // - if (mLookupTableIsFilled) { - for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i) - delete [] mLookupTable[i]; - delete [] mLookupTable; - } - - mNumberOfNodesInX = nx; - mNumberOfNodesInQ2 = nq2; - - // - // Create new table - // - mLookupTable = new double*[mNumberOfNodesInQ2]; - for(unsigned int i = 0; i < mNumberOfNodesInQ2; ++i) - mLookupTable[i] = new double[mNumberOfNodesInX]; - - // - // Fill lookup table - // - int printEach = mNumberOfNodesInQ2*mNumberOfNodesInX/10; - int nCount = 0; - cout << "DglapEvolution: generating lookup table "; cout.flush(); - for (unsigned int i = 0; i < mNumberOfNodesInQ2; i++) { - double Q2 = pow(mTableMaxQ2/mTableMinQ2,i*1./(mNumberOfNodesInQ2-1.))*mTableMinQ2; - for (unsigned int j = 0; j < mNumberOfNodesInX; j++) { - double x = pow(mTableMaxX/mTableMinX,j*1./(mNumberOfNodesInX-1.))*mTableMinX; - mLookupTable[i][j] = xG_Engine(x, Q2); - nCount++; - if (nCount%printEach == 0) cout << '.'; cout.flush(); - } - } - cout << " done" << endl; - mLookupTableIsFilled = true; +void DglapEvolution::generateLookupTable(int nx, int nq2) +{ + // + // Delete old lookup table (if it exists) + // + if (mLookupTableIsFilled) { + for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i) + delete [] mLookupTable[i]; + delete [] mLookupTable; + } + + mNumberOfNodesInX = nx; + mNumberOfNodesInQ2 = nq2; + + // + // Create new table + // + mLookupTable = new double*[mNumberOfNodesInQ2]; + for(unsigned int i = 0; i < mNumberOfNodesInQ2; ++i) + mLookupTable[i] = new double[mNumberOfNodesInX]; + + // + // Fill lookup table + // + int printEach = mNumberOfNodesInQ2*mNumberOfNodesInX/10; + int nCount = 0; + cout << "DglapEvolution: generating lookup table "; cout.flush(); + for (unsigned int i = 0; i < mNumberOfNodesInQ2; i++) { + double Q2 = pow(mTableMaxQ2/mTableMinQ2,i*1./(mNumberOfNodesInQ2-1.))*mTableMinQ2; + for (unsigned int j = 0; j < mNumberOfNodesInX; j++) { + double x = pow(mTableMaxX/mTableMinX,j*1./(mNumberOfNodesInX-1.))*mTableMinX; + mLookupTable[i][j] = xG_Engine(x, Q2); + nCount++; + if (nCount%printEach == 0) cout << '.'; cout.flush(); + } + } + cout << " done" << endl; + mLookupTableIsFilled = true; } double DglapEvolution::luovi(double f[4], double arg[4], double z) { - double cof[4]; + double cof[4]; for (int i=0; i < 4; i++ ) cof[i]=f[i]; for (int i=0; i < 3 ; i++) { for (int j=i; j < 3 ; j++) { int jndex = 2 - j; int index = jndex + 1 + i; cof[index]= (cof[index]-cof[index-1])/(arg[index]-arg[jndex]); } } double sum=cof[3]; for (int i=0; i < 3; i++ ) { int index = 2 - i; sum = (z-arg[index])*sum + cof[index]; } return sum; } double DglapEvolution::xG_Interpolator(double x, double Q2) -{ - int q2steps = mNumberOfNodesInQ2-1; - int xsteps = mNumberOfNodesInX-1; - - // - // Position in the Q2 grid +{ + int q2steps = mNumberOfNodesInQ2-1; + int xsteps = mNumberOfNodesInX-1; + + if(Q2>mTableMaxQ2){ + cout<<"DglapEvolution::xg_Interpolator, error: Interpolator called with Q2="<(realq); if (n_q2 <= 0) {n_q2 = 1;} if (n_q2 > q2steps-2) {n_q2 = n_q2-2;} - + // - // Position in the x grid + // Position in the x grid // double realx = xsteps * log(x/mTableMinX)/log(mTableMaxX/mTableMinX); int n_x = static_cast(realx); if (n_x <= 0) { n_x = 1;} if (n_x > xsteps-2){ n_x = n_x-2;} - // + // // Starting the interpolation - // + // double arg[4]; for (int i=0; i < 4; i++) { arg[i] = n_x-1+i;} - + double fu[4], fg[4]; for (int i=0; i < 4; i++) { fu[0] = mLookupTable[n_q2-1+i][n_x-1]; fu[1] = mLookupTable[n_q2-1+i][n_x]; fu[2] = mLookupTable[n_q2-1+i][n_x+1]; fu[3] = mLookupTable[n_q2-1+i][n_x+2]; fg[i] = luovi(fu,arg,realx); } for (int i=0; i < 4; i++) { arg[i] = n_q2-1+i;} return luovi(fg, arg, realq); } void DglapEvolution::reno(complex *fn, double alpq, int nmax, double ag, double lambdag) { // // Mellin-n space Q**2 - evolution of the gluon at LO - // + // // The moments are calculated on an array of moments, mN, suitable // for a (non-adaptive) Gauss quadrature. - // + // // Currently this takes the simplest possible fit form: // xg = A_g x^(-lambdag) (1-x)^(5.6), following Amir&Raju - // + // for (int k1 = 1; k1 <= nmax; k1++) { // // Input moments of the parton densities // at the low scale // complex xn = mN[k1]; complex zero; complex one(1.,0.); complex two(2.,0.); complex gln = ag * (1.0 / (xn + 5.0 - lambdag) - 6.0 / (xn + 4.0 - lambdag) + 15.0 / (xn + 3.0 - lambdag) - 20.0 / (xn + 2.0 -lambdag) + 15.0 / (xn + 1.0 - lambdag) - 6.0 / (xn - lambdag) + 1.0 / (xn - lambdag - 1.0)); int f; double xl, s, alp; complex ep, gl; if (alpq >= mALPC) { // evolution below the charm threshold f = 3; xl = mALPS / alpq; s = log(xl); ep = exp(-mAP[k1][f]*s); gl = gln; gln = ep * gl; } else if ((alpq < mALPC) && (alpq >= mALPB)) { // between thresholds f = 3; xl = mALPS / mALPC; s = log(xl); ep = exp(-mAP[k1][f]*s); gl = gln; gln = ep * gl; f = 4; xl = mALPC / alpq; s = log(xl); ep = exp(-mAP[k1][f]*s); gl = gln; gln = ep * gl; } else if (alpq < mALPB) { // above bottom threshold f = 3; xl = mALPS / mALPC; s = log (xl); ep = exp (- mAP[k1][f]*s); gl = gln; gln = ep * gl; f = 4; alp = mALPB; xl = mALPC / mALPB; s = log (xl); ep = exp(-mAP[k1][f]*s); gl = gln; gln = ep * gl; f = 5; xl = mALPB / alpq; s = log(xl); ep = exp(-mAP[k1][f]*s); gl = gln; gln = ep * gl; } fn[k1] = gln; } } void DglapEvolution::anom() { // // Anomalous dimensions for leading order evolution of parton densities. // The moments are calculated on an externally given array of mellin // moments, mN, suitable for a (non-adaptive) quadrature. // // Present version: the number of active flavours in the factorization // is fixed to ff=3, in the beta-function it varies between f=3 and // f=5. The dimension of the moment array is 136. // double b0, b02; complex ggi, ggf; complex xn, xap; complex gg; for (int k1=1; k1 <= 136; k1++) { xn = mN[k1]; anCalc(ggi, ggf, xn); for (int k2=3; k2 <= 5; k2++) { double f = k2; // anomalous dimensions and related quantities in leading order b0 = 11.- 2./3.* f; b02 = 2.* b0; gg = ggi + f * ggf; xap = gg / b02; mAP[k1][k2] = xap; } } } void DglapEvolution::anCalc(complex& ggi, complex& ggf, complex& xn) { complex xn1 = xn + 1.; complex xn2 = xn + 2.; complex xnm = xn - 1.; // // Leading order // complex cpsi = psiFunction(xn1) + 0.577216; ggi = -22.- 24./(xn * xnm) - 24./(xn1 * xn2) + 24.* cpsi; ggf = 4./3.; } complex DglapEvolution::psiFunction(complex z) { // // psi - function for complex argument // complex sub; while (real(z) < 10) { sub = sub - 1./z; z += 1.; } complex rz = 1./z; complex dz = rz * rz; complex result = sub + log(z) - rz/2.- dz/2520. * ( 210.+ dz * (-21.+10.*dz )); return result; -} - -double DglapEvolution::logDerivative(double x, double Q2) -{ - double alpq = mAlphaStrong->at(Q2)/(4*M_PI); - - // - // Q2 and x dependent quantities - // - double ax = log(x); - double ex = exp(-mC * ax); - - // - // integration length parameter for the mellin inversion - // - const int nmax = 136; - - // - // Gluon density and output - // - complex fn[137]; - reno(fn, alpq, nmax, mAg, mLambdaG); - - long double fun = 0; // long double is needed - long double fz; - complex xnm,cex; - - for (int i=1; i <= nmax; i++) { - xnm = (mC - mN[i]+1.) * ax; - cex = exp(xnm) / M_PI * mCC; - fz = imag(fn[i]*cex); - fun = fun + mWN[i] * fz; - } - double glue = fun * ex; - - fun = 0; - for (int i=1; i <= nmax; i++) { - xnm = (mC - mN[i]+1.) * ax; - cex = exp(xnm) / M_PI * mCC; - fz = imag(fn[i]*cex*mN[i]); - fun = fun + mWN[i] * fz; - } - double pa = fun * ex; - - double lambda = -(1-pa/glue); - - return lambda; -} +} + +double DglapEvolution::logDerivative(double x, double Q2) +{ + double alpq = mAlphaStrong->at(Q2)/(4*M_PI); + + // + // Q2 and x dependent quantities + // + double ax = log(x); + double ex = exp(-mC * ax); + + // + // integration length parameter for the mellin inversion + // + const int nmax = 136; + + // + // Gluon density and output + // + complex fn[137]; + reno(fn, alpq, nmax, mAg, mLambdaG); + + long double fun = 0; // long double is needed + long double fz; + complex xnm,cex; + + for (int i=1; i <= nmax; i++) { + xnm = (mC - mN[i]+1.) * ax; + cex = exp(xnm) / M_PI * mCC; + fz = imag(fn[i]*cex); + fun = fun + mWN[i] * fz; + } + double glue = fun * ex; + + fun = 0; + for (int i=1; i <= nmax; i++) { + xnm = (mC - mN[i]+1.) * ax; + cex = exp(xnm) / M_PI * mCC; + fz = imag(fn[i]*cex*mN[i]); + fun = fun + mWN[i] * fz; + } + double pa = fun * ex; + + double lambda = -(1-pa/glue); + + return lambda; +} Index: branches/fitting/src/Constants.h =================================================================== --- branches/fitting/src/Constants.h (revision 379) +++ branches/fitting/src/Constants.h (revision 380) @@ -1,49 +1,49 @@ //============================================================================== // Constants.h // // Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich // // This file is part of Sartre. // // 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 . // // Author: Thomas Ullrich // Last update: // $Date$ // $Author$ //============================================================================== #ifndef Constants_h #define Constants_h // // General constants // const double electronMass = 0.510998902E-3; // GeV const double electronMass2 = electronMass*electronMass; // GeV^2 const double protonMass = 0.9382700; // GeV const double protonMass2 = protonMass*protonMass; // GeV^2 const double alpha_em = 1/137.036; const double hbarc = 0.197327; // GeV*fm const double hbarc2 = hbarc*hbarc; // GeV^2*fm^2 // // Constants used in dipole model // const double Nc=3.; const double quarkCharge[6] = {-1./3., 2./3., -1./3., 2./3., -1./3., 2./3.};//d, u, s, c, b, t - + // // Constants used in integartion and other numerical operations // const double upperIntegrationLimit=2.5; // factor: nuclear radius -> integration limits (b, r) #endif