Index: branches/fitting/src/Integrals.cpp =================================================================== --- branches/fitting/src/Integrals.cpp (revision 380) +++ branches/fitting/src/Integrals.cpp (revision 381) @@ -1,494 +1,498 @@ //============================================================================== // Integrals.cpp // // 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: Tobias Toll // Last update: // $Date$ // $Author$ //============================================================================== #include #include #include #include "Integrals.h" #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 "IntegrandWrappers.h" #include "TF1.h" #include "TH1F.h" #include "cuba.h" #define PRf(x) printf(#x); printf("=%g \n", (x)); #define PRi(x) printf(#x); printf("=%d \n", (x)); #define PR(x) cout << #x << " = " << (x) << endl; using namespace std; Integrals::Integrals() { mIsInitialized=false; mRelativePrecisionOfIntegration = 0; mWaveOverlap = 0; mDipoleModel = 0; mIntegralImT = 0; mIntegralImL = 0; mIntegralReT = 0; mIntegralReL = 0; mErrorImT = 0; mErrorImL = 0; mErrorReT = 0; mErrorReL = 0; mProbImT = 0; mProbImL = 0; mProbReT = 0; mProbReL = 0; TableGeneratorSettings* settings = TableGeneratorSettings::instance(); int VMId=settings->vectorMesonId(); mMV = settings->lookupPDG(VMId)->Mass(); if (VMId==113 || VMId==333 || VMId == 443) { mWaveOverlap = new WaveOverlapVM; mWaveOverlap->setProcess(VMId); mWaveOverlap->setWaveOverlapFunctionParameters(VMId); } else if (VMId==22) { mWaveOverlap = new WaveOverlapDVCS; } else { cout << "Integrals::init(): Error, no exclusive production implemented for: "<< VMId << endl; exit(1); } DipoleModelType model=settings->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 if (model==bCGC) { +// mDipoleModel = new DipoleModel_bCGC; +// } else { cout << "Integrals::init(): Error, model not implemented: "<< model << endl; exit(1); } mVerbose=settings->verbose(); mIsInitialized=true; } Integrals::Integrals(const Integrals& integrals) { mIsInitialized = integrals.mIsInitialized; mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration; if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS)) mWaveOverlap = new WaveOverlapDVCS; else mWaveOverlap = new WaveOverlapVM; if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)) mDipoleModel = new DipoleModel_bSat; else if(typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat)) mDipoleModel = new DipoleModel_bNonSat; - else - mDipoleModel = new DipoleModel_bCGC; +// else +// mDipoleModel = new DipoleModel_bCGC; + else{ + cout<<"Wrong Model!"<configurationExists()) { // do not use cout cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl; exit(1); } //store present kinematic point: if (setKinematicPoint(t, Q2, W2)) { calculate(); } else { fillZeroes(); } } void Integrals::operator() (double t, double xpom) //UPC { //make sure the configurations have been generated: if (!mDipoleModel->configurationExists()) { // do not use cout cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl; exit(1); } //store present kinematic point: if (setKinematicPoint(t, xpom)) { calculate(); } else { fillZeroes(); } } void Integrals::fillZeroes(){ //Store the results mIntegralImT=0; mIntegralReT=0; mIntegralImL=0; mIntegralReL=0; //Store the errors: mErrorImT=0; mErrorImL=0; mErrorReT=0; mErrorReL=0; //Store the probabilities: mProbImT=0; mProbImL=0; mProbReT=0; mProbReL=0; } //*********EXCLUSIVE VECTOR MESONS OR DVCS: ******************************** void IntegralsExclusive::coherentIntegrals(double t, double Q2, double W2) { if(setKinematicPoint(t, Q2, W2)){ if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){ //store present kinematic point: double xprobe=kinematicPoint[3]; dipoleModel()->createSigma_ep_LookupTable(xprobe); } calculateCoherent(); } else fillZeroes(); } void IntegralsExclusive::coherentIntegrals(double t, double xpom) { if(setKinematicPoint(t, xpom)){ if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){ dipoleModel()->createSigma_ep_LookupTable(xpom); } calculateCoherent(); } else fillZeroes(); } IntegralsExclusive::IntegralsExclusive() { } IntegralsExclusive& IntegralsExclusive::operator=(const IntegralsExclusive& cobj) { if (this != &cobj) { Integrals::operator=(cobj); copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint); } return *this; } IntegralsExclusive::IntegralsExclusive(const IntegralsExclusive& cobj) : Integrals(cobj) { copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint); } bool IntegralsExclusive::setKinematicPoint(double t, double xpom) //UPC { bool result = true; kinematicPoint[0]=t; kinematicPoint[1]=0; kinematicPoint[2]=0; //W2 is not used kinematicPoint[3]=xpom; return result; } bool IntegralsExclusive::setKinematicPoint(double t, double Q2, double W2) { bool result = true; kinematicPoint[0]=t; kinematicPoint[1]=Q2; kinematicPoint[2]=W2; double xprobe=Kinematics::xpomeron(t, Q2, W2, mMV); if (xprobe<0 || xprobe>1) result = false; kinematicPoint[3]=xprobe; return result; } void IntegralsExclusive::calculate() { // // This function calls a wrapper from where the // integral is calculated with the Cuhre method. // Pass this Integrals object as the fourth (void*) argument of the Cuhre function. // const double epsrel=1.e-3, epsabs=1e-12; const int flags=0, mineval=1e4, maxeval=1e8, key=0; int nregionsTIm, nevalTIm, failTIm; int nregionsTRe, nevalTRe, failTRe; int nregionsLIm, nevalLIm, failLIm; int nregionsLRe, nevalLRe, failLRe; double valTIm=0, errTIm=0, probTIm=0; double valLIm=0, errLIm=0, probLIm=0; double valTRe=0, errTRe=0, probTRe=0; double valLRe=0, errLRe=0, probLRe=0; const char* statefile=0; const int nvec=1; // double probabilityCutOff=1e-6; unsigned int A=dipoleModel()->nucleus()->A(); // // Do the integrations // Cuhre(4, 1, integrandWrapperTIm, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsTIm, &nevalTIm, &failTIm, &valTIm, &errTIm, &probTIm); if(failTIm!=0 and mVerbose) printf("IntegralsExclusive::calculate(): Warning: Integration TIm did not reach desired precision! Error code=%d \n", failTIm); Cuhre(4, 1, integrandWrapperLIm, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsLIm, &nevalLIm, &failLIm, &valLIm, &errLIm, &probLIm); if(failLIm!=0 and mVerbose) printf("IntegralsExclusive::calculate(): Warning: Integration LIm did not reach desired precision! Error code=%d \n", failLIm); // // In the case of proton, there is angular symmetry, which makes the real part zero. // if(A>1){ Cuhre(4, 1, integrandWrapperTRe, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsTRe, &nevalTRe, &failTRe, &valTRe, &errTRe, &probTRe); if(failTRe!=0 and mVerbose) printf("IntegralsExclusive::calculate(): Warning: Integration TRe did not reach desired precision! Error code=%d \n", failTRe); Cuhre(4, 1, integrandWrapperLRe, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsLRe, &nevalLRe, &failLRe, &valLRe, &errLRe, &probLRe); if(failLRe!=0 and mVerbose) printf("IntegralsExclusive::calculate(): Warning: Integration LRe did not reach desired precision! Error code=%d \n", failLRe); }//if A // // Store the results: // mIntegralImT=valTIm; mIntegralReT=valTRe; mIntegralImL=valLIm; mIntegralReL=valLRe; // // Store the errors: // mErrorImT=errTIm; mErrorImL=errLIm; mErrorReT=errTRe; mErrorReL=errLRe; // // Store the probabilities: // mProbImT=probTIm; mProbImL=probLIm; mProbReT=probTRe; mProbReL=probLRe; } void IntegralsExclusive::calculateCoherent() { // // As calculate() but for the coherent case // const double epsrel=1e-6, epsabs=1e-12; const int flags=0, mineval=1e1, maxeval=1e8, key=0; int nregionsT, nevalT, failT; int nregionsL, nevalL, failL; double valT, errT, probT; double valL, errL, probL; const int nvec=1; const char* statefile=0; // // Do the integrations // Cuhre(3, 1, integrandWrapperCoherentAmplitudeT, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT); if(failT!=0 and mVerbose) printf("IntegralsExclusive::calculateCoherent(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT); Cuhre(3, 1, integrandWrapperCoherentAmplitudeL, this, nvec, epsrel, epsabs, flags, mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL); if(failL!=0 and mVerbose) printf("IntegralsExclusive::calculateCoherent(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL); // // Store the results // mIntegralImT=valT; mIntegralImL=valL; // // Store the errors // mErrorImT=errT; mErrorImL=errL; } // // The following functions are the Integrands in the Amplitudes: // double IntegralsExclusive::uiAmplitudeTIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta) { double cosArg = (b/hbarc)*Delta*cos(phi); double waveOverlap = mWaveOverlap->T(z, Q2, r); double dsigdb2 = dipoleModel()->dsigmadb2(r , b , phi, xprobe); double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc); double result=0.5*r/hbarc2*waveOverlap* BesselJ0*b*cos(cosArg)*dsigdb2; return result; } double IntegralsExclusive::uiAmplitudeTRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta) { double sinArg = b*Delta*cos(phi)/hbarc; double waveOverlap = mWaveOverlap->T(z, Q2, r); double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe); double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc); double result=0.5*r/hbarc2*waveOverlap* BesselJ0*b*sin(sinArg)*dsigdb2; return result; } double IntegralsExclusive::uiAmplitudeLIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta) { double waveOverlap = mWaveOverlap->L(z, Q2, r); double cosArg = b*Delta*cos(phi)/hbarc; double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe); double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc); double result=0.5*r/hbarc2*waveOverlap* BesselJ0*b*cos(cosArg)*dsigdb2; return result; } double IntegralsExclusive::uiAmplitudeLRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta) { double waveOverlap = mWaveOverlap->L(z, Q2, r); double sinArg = b*Delta*cos(phi)/hbarc; double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe); double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc); double result=0.5*r/hbarc2*waveOverlap* BesselJ0*b*sin(sinArg)*dsigdb2; return result; } double IntegralsExclusive::uiCoherentAmplitudeT(double b, double z, double r, double Q2, double Delta) { double waveOverlap = mWaveOverlap->T(z, Q2, r); double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc); double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc); double xprobe=kinematicPoint[3]; double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe); double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean; return result; } double IntegralsExclusive::uiCoherentAmplitudeL(double b, double z, double r, double Q2, double Delta) { double waveOverlap = mWaveOverlap->L(z, Q2, r); double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc); double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc); double xprobe=kinematicPoint[3]; double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe); double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean; return result; } Index: branches/fitting/src/Parameters_Minuit2.h =================================================================== --- branches/fitting/src/Parameters_Minuit2.h (revision 0) +++ branches/fitting/src/Parameters_Minuit2.h (revision 381) @@ -0,0 +1,71 @@ +#ifndef parameters_minuit2_h +#define parameters_minuit2_h + + +#include +#include +#include +#include "Enumerations.h" +#include + +using namespace std; +using namespace ROOT::Minuit2; + +struct FitParameters //Temporary location +{ + const vector *values; + const MnUserParameters *parameter; + double alphas_mur; +}; + +const double params_l_mass =0.03;//0.1001; +const double params_c_mass =1.3210; +const double params_b_mass = 4.75; +const double params_t_mass = 172.; +const double params_xtil_c_mass = 1.3210; + +const double params_mu02 =1.1; + +const double params_B_G = 4.; +const double params_A_g =2.0670; +const double params_C =1.8178; +const double params_lambda_g = 0.09575; +const double params_r_max =0.8678; + +const double params_BV =3.; +const double params_RV =0.2; +const double params_wV =9.; + + +const double params_A_s = 0; +const double params_lambda_s = 0; + +const double params_coupling = 0.1184; +const int params_asOrder = 1; + +const double params_as_qqg = 0.; + +const double params_Ctilde = 48.0; +//const double params_quarkMass[6] = {0.0, 0.0, 0.0, 1.27, 4.18, 172.}; // b and t not in fit +//Shapes +//const ProtonShape params_Shape=Gauss; +const ProtonShape params_Shape=Gauss; + + +//Hard Sphere +const double params_Rhs =7.8004; //GeV-1 + + +//Woods-Saxon +const double params_R0WS=7.8; //GeV^-1 +const double params_dWS=0.01; //GeV^-1 + + +//Laplace +const double params_nLP = 1.8282; + +//CompressedExponent +const double params_uCE = 0.7538246273225; +const double params_EP = 1.176352891591; + +#endif Index: branches/fitting/src/Amplitudes.cpp =================================================================== --- branches/fitting/src/Amplitudes.cpp (revision 380) +++ branches/fitting/src/Amplitudes.cpp (revision 381) @@ -1,311 +1,311 @@ //============================================================================== // Amplitudes.cpp // // 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: Tobias Toll // Last update: // $Date$ // $Author$ //============================================================================== -//#define SARTRE_IN_MULTITHREADED_MODE 1 +#define SARTRE_IN_MULTITHREADED_MODE 1 #include #include #include "Amplitudes.h" #include "Constants.h" #include "TableGeneratorSettings.h" #include "DglapEvolution.h" #include "Enumerations.h" #include "Kinematics.h" #include "Integrals.h" #include "DipoleModel.h" #if defined(SARTRE_IN_MULTITHREADED_MODE) #include -#endif +#endif #define PR(x) cout << #x << " = " << (x) << endl; using namespace std; Amplitudes::Amplitudes() { mAmplitudeT = 0; mAmplitudeL = 0; mAmplitudeT2 = 0; mAmplitudeL2 = 0; mNumberOfConfigurations = 0; mTheModes = 0; mA = 0; mErrorT = 0; mErrorL = 0; mErrorT2 = 0; mErrorL2 = 0; TableGeneratorSettings* settings = TableGeneratorSettings::instance(); mNumberOfConfigurations = settings->numberOfConfigurations(); mVerbose = settings->verbose(); // // Create a vector containing instances of the Integrals class // and initialize them: // for (int i=0; i<=mNumberOfConfigurations; i++) { mIntegrals.push_back(new IntegralsExclusive); } mA=settings->A(); mUPC=settings->UPC(); // // Get the modes to calculate: // 0: analytically averaged over configurations // 1: only analytically // 2: Both and averaged over configurations // mTheModes=settings->modesToCalculate(); } Amplitudes& Amplitudes::operator=(const Amplitudes& amp) { if (this != &) { for (unsigned int i=0; idipoleModel()->createConfiguration(i); } //void Amplitudes::calculate(double t, double Q2, double W2) void Amplitudes::calculate(double* kinematicPoint) { double t=0, Q2=0, W2=0, xpom=0; if(!mUPC){ t =kinematicPoint[0]; Q2=kinematicPoint[1]; W2=kinematicPoint[2]; } else{ t =kinematicPoint[0]; xpom=kinematicPoint[1]; } #if defined(SARTRE_IN_MULTITHREADED_MODE) // multithreaded version if (mA == 1) { cout << "Amplitudes::calculate(): Multithreaded mode (SARTRE_IN_MULTITHREADED_MODE)" << endl; cout << " is not supported for ep (A=1). Stopping." << endl; exit(1); } // // Create a vector containing the threads: // std::vector vThreads; vThreads.clear(); // // Create the thread group: // boost::thread_group gThreads; if (mTheModes==0 || mTheModes == 2){ //Start loop over configurations, each calculated on a separate thread: for (int i=0; i1 && (mTheModes==1 || mTheModes == 0)) { if(!mUPC) mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2); else mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom); } if (mTheModes==0 || mTheModes == 2) { // Wait for all threads to finish before continuing main thread: gThreads.join_all(); // Clean up threads vThreads.clear(); } #else // unforked version if ((mTheModes==0 || mTheModes == 2) || mA==1) { //Start loop over configurations: for (int i=0; ioperator()(t, Q2, W2); else mIntegrals.at(i)->operator()(t, xpom); } } // // Calculate coherent cross-section according to eq.(47) in KT arXiv:hep-ph/0304189v3, // (only in eA) // if (mA>1 && (mTheModes==1 || mTheModes == 0)){ if(!mUPC) mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2); else mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom); } #endif // // Calculate the resulting : // double coherentT=0, coherentL=0; double errCoherentT=0, errCoherentL=0; if ((mTheModes==0 || mTheModes == 2) || mA==1) { double totalT = 0; double totalL = 0; double err2TotalT=0, err2TotalL=0; double probabilityCutOff = 1e-6; for (int i=0; iintegralImT(); double valreT=mIntegrals.at(i)->integralReT(); double valimL=mIntegrals.at(i)->integralImL(); double valreL=mIntegrals.at(i)->integralReL(); double errimT=mIntegrals.at(i)->errorImT(); double errimL=mIntegrals.at(i)->errorImL(); double errreT=mIntegrals.at(i)->errorReT(); double errreL=mIntegrals.at(i)->errorReL(); double probimT=mIntegrals.at(i)->probImT(); double probimL=mIntegrals.at(i)->probImL(); double probreT=mIntegrals.at(i)->probReT(); double probreL=mIntegrals.at(i)->probReL(); if (probimT > probabilityCutOff || probimL > probabilityCutOff || probreT > probabilityCutOff || probreL > probabilityCutOff){ if(mVerbose) { cout<< "Amplitudes::calculate(): Warning, Integrals may not have reached desired precision" < probreT and probimT > probimL and probimT > probreL ? cout<< " The probability for this is "< probimT and probreT > probimL and probreT > probreL ? cout<< " The probability for this is "< probimT and probimL > probreT and probimL > probreL ? cout<< " The probability for this is "<1 && (mTheModes==1 || mTheModes == 0)) { double coherentKTT=mIntegrals.at(mNumberOfConfigurations)->integralImT(); double coherentKTL=mIntegrals.at(mNumberOfConfigurations)->integralImL(); double errCoherentKTT=mIntegrals.at(mNumberOfConfigurations)->errorImT(); double errCoherentKTL=mIntegrals.at(mNumberOfConfigurations)->errorImL(); mAmplitudeT=coherentKTT; mAmplitudeL=coherentKTL; mErrorT=errCoherentKTT; mErrorL=errCoherentKTL; } else { mAmplitudeT=coherentT/mNumberOfConfigurations; mAmplitudeL=coherentL/mNumberOfConfigurations; mErrorT=errCoherentT/mNumberOfConfigurations; mErrorL=errCoherentL/mNumberOfConfigurations; } } Index: branches/fitting/src/InclusiveDiffractiveCrossSections.h =================================================================== --- branches/fitting/src/InclusiveDiffractiveCrossSections.h (revision 380) +++ branches/fitting/src/InclusiveDiffractiveCrossSections.h (revision 381) @@ -1,216 +1,243 @@ // // To calculate total cross-sections use // void calculate(double, double, double, double); // This integrates over z, it is a 5D integration. N.B. Slow! // // For exclusive final state use // void calculate(double, double, double, double, double); // This does not inteagrate over z, and is a much faster 2D integration // -#ifndef InclusiveDiffractiveCrossSections_h -#define InclusiveDiffractiveCrossSections_h +#ifndef InclusiveDiffractiveCrossSections_h +#define InclusiveDiffractiveCrossSections_h #include "AlphaStrong.h" #include "TH1.h" -#include "TMath.h" -#include "TF1.h" +#include "TMath.h" +#include "TF1.h" #include "Math/WrappedTF1.h" #include "Math/GaussIntegrator.h" class DipoleModel; class DipoleModelParameters; class InclusiveDiffractiveCrossSections { - - public: - InclusiveDiffractiveCrossSections(); - ~InclusiveDiffractiveCrossSections(); - InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections&); - InclusiveDiffractiveCrossSections& operator=(const InclusiveDiffractiveCrossSections&); - - DipoleModel* dipoleModel() const; - - void init(); - void setZ(double); - void setQuarkIndex(int); - double z(); - int quarkIndex(); - double zlower(); - void setNumberOfFlavours(unsigned int); - void setWhichFlavours(bool, bool, bool, bool, bool); - - void calculate_dsigmadbeta_QQ(double, double, double); - void calculate_dsigmadbetadz_QQ(double, double, double, double); - - double* dsigmadbetaT_QQ(); - double* dsigmadbetaL_QQ(); - - void setBetaMin(double); - double betaMin(); - - AlphaStrong mAs; - - void setRelativePrecisionOfIntegration(double); - - double mUpperBinIntegration; - double mUpperRinIntegrationQQ; - double mUpperRinIntegrationDIS; - double mUpperRinIntegrationQQG; - - double as_qqg_fixed; - - //Total (non-diffractive) cross-sections: - double uiDipoleSigma(double*, double*); - double waveOverlapT(double, double, double); - double waveOverlapL(double, double, double); - double uiDsigmadrT(double*, double*); - double uiDsigmadrL(double*, double*); - double uiDsigmadzT(double*, double*); - double uiDsigmadzL(double*, double*); - void calculateTotal(double, double); - - //QQG - void putAllQQGtoZero(); - void calculate_dsigmadbeta_QQG(double, double, double); - void calculate_dsigmadbetadz_QQG(double, double, double, double); - double* dsigmadbetaT_QQG(); - - //QQG_MS - double uiQQG_MS(const double*); - double QQG_MS_Nfactors(double, double, double, double, double); - double uiQQGdz_MS(const double*); - void FTQQG_MS(double, double); - void FTQQGdz_MS(double, double, double); - double* dsigmadbetaT_QQG_MS(); - - //QQG_GBW: - void FTQQG_GBW(double, double); - void FTQQGdz_GBW(double, double, double); - double dsigmadb2tilde(double, double, double); - double uiR_GBW(double*, double*); - double uiQQG_GBW(const double*); - double uiQQGdz_GBW(const double*); - double* dsigmadbetaT_QQG_GBW(); - //QQG_GBW beta->0: - void FTQQG_GBWbeta0(double); - double uiQQG_GBWbeta0(const double*); - double uiR_GBWbeta0(double*, double*); - double* dsigmadbetaT_QQG_GBWbeta0(); - - double uiTotalDIST(const double*); - double uiTotalDISL(const double*); - - double mTotalCST; - double mTotalCSL; - //End Total (non-diffractive) cross-sections: - - double uiPhi0(double*, double*); - double uiPhi1(double*, double*); - bool setKinematicPoint(double, double, double, double); - private: - // bool setKinematicPoint(double, double, double, double); - double uiIntegralL(double*, double*); - double uiIntegralT(double*, double*); - double dPhi1db(double*, double*); - double dPhi0db(double*, double*); - - double buiA(double*, double*); - double phiuiA(double*, double*); - double ruiA(double*, double*); - double zui(double*, double*); - double rui(double*, double*); - - // double uiPhi0(double*, double*); - // double uiPhi1(double*, double*); - - double kinematicPoint[6]; - double mZ; - double mR; - double mB; - double mBeta; - DipoleModel* mDipoleModel; - int mQuarkIndex; - bool mIsInitialized; - double mInclusiveT[6]; - double mInclusiveL[6]; - double mInclusiveQQG[6]; - double mInclusiveQQG_MS[6]; - double mInclusiveQQG_GBW[6]; - double mInclusiveQQG_GBWbeta0[6]; - unsigned int mA; - unsigned int mNumberOfFlavours; - bool mWhichFlavours[5]; - double mRelativePrecisionOfIntegration; - double mBetaMin; - double mQ2; - double mW2; - - DipoleModelParameters* mParameters; - - TF1* mIntegralT; - ROOT::Math::WrappedTF1* mWFT; - ROOT::Math::GaussIntegrator mGIintegralT; - - TF1* mIntegralL; - ROOT::Math::WrappedTF1* mWFL; - ROOT::Math::GaussIntegrator mGIintegralL; - - TF1* mPhi0; - ROOT::Math::WrappedTF1* mWFPhi0; - ROOT::Math::GaussIntegrator mGIPhi0; - - TF1* mPhi1; - ROOT::Math::WrappedTF1* mWFPhi1; - ROOT::Math::GaussIntegrator mGIPhi1; - - TF1* mPhi0r; - ROOT::Math::WrappedTF1* mWFPhi0r; - ROOT::Math::GaussIntegrator mGIPhi0r; - - TF1* mPhi1r; - ROOT::Math::WrappedTF1* mWFPhi1r; - ROOT::Math::GaussIntegrator mGIPhi1r; - - TF1* mRGBW; - ROOT::Math::WrappedTF1* mWFRGBW; - ROOT::Math::GaussIntegrator mGIRGBW; - - TF1* mRGBWbeta0; - ROOT::Math::WrappedTF1* mWFRGBWbeta0; - ROOT::Math::GaussIntegrator mGIRGBWbeta0; - - TF1* mTotalT; - ROOT::Math::WrappedTF1* mWFTotalT; - ROOT::Math::GaussIntegrator mGITotalT; - - TF1* mTotalL; - ROOT::Math::WrappedTF1* mWFTotalL; - ROOT::Math::GaussIntegrator mGITotalL; - - TF1* mSigmaT; - ROOT::Math::WrappedTF1* mWFSigmaT; - ROOT::Math::GaussIntegrator mGISigmaT; - - TF1* mSigmaL; - ROOT::Math::WrappedTF1* mWFSigmaL; - ROOT::Math::GaussIntegrator mGISigmaL; - - TF1* mDipole; - ROOT::Math::WrappedTF1* mWFDipole; - ROOT::Math::GaussIntegrator mGIDipole; - + +public: + InclusiveDiffractiveCrossSections(); + ~InclusiveDiffractiveCrossSections(); + InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections&); + InclusiveDiffractiveCrossSections& operator=(const InclusiveDiffractiveCrossSections&); + + DipoleModel* dipoleModel() const; + + void init(); + void setZ(double); + void setQuarkIndex(int); + double z(); + int quarkIndex(); + double zlower(); + void setNumberOfFlavours(unsigned int); + void setWhichFlavours(bool, bool, bool, bool, bool); + + void calculate_dsigmadbeta_QQ(double, double, double); + void calculate_dsigmadbetadz_QQ(double, double, double, double); + + double* dsigmadbetaT_QQ(); + double* dsigmadbetaL_QQ(); + + void setBetaMin(double); + double betaMin(); + + AlphaStrong mAs; + + void setRelativePrecisionOfIntegration(double); + + double mUpperBinIntegration; + double mUpperRinIntegrationQQ; + double mUpperRinIntegrationDIS; + double mUpperRinIntegrationQQG; + + double as_qqg_fixed; + + //Total (non-diffractive) cross-sections: + double uiDipoleSigma(double*, double*); + double waveOverlapT(double, double, double); + double waveOverlapL(double, double, double); + void calculateTotal(); + void operator()(double, double); + + //QQG + void putAllQQGtoZero(); + void calculate_dsigmadbeta_QQG(double, double, double); + void calculate_dsigmadbetadz_QQG(double, double, double, double); + double* dsigmadbetaT_QQG(); + + //QQG_MS + double uiQQG_MS(const double*); + double QQG_MS_Nfactors(double, double, double, double, double); + double uiQQGdz_MS(const double*); + void FTQQG_MS(double, double); + void FTQQGdz_MS(double, double, double); + double* dsigmadbetaT_QQG_MS(); + + //QQG_GBW: + void FTQQG_GBW(double, double); + void FTQQGdz_GBW(double, double, double); + double dsigmadb2tilde(double, double, double); + double uiR_GBW(double*, double*); + double uiQQG_GBW(const double*); + double uiQQGdz_GBW(const double*); + double* dsigmadbetaT_QQG_GBW(); + //QQG_GBW beta->0: + void FTQQG_GBWbeta0(double); + double uiQQG_GBWbeta0(const double*); + double uiR_GBWbeta0(double*, double*); + double* dsigmadbetaT_QQG_GBWbeta0(); + + double uiTotalDIST(const double*); + double uiTotalDISL(const double*); + + double mTotalCST; + double mTotalCSL; + //End Total (non-diffractive) cross-sections: + + double uiPhi0(double*, double*); + double uiPhi1(double*, double*); + bool setKinematicPoint(double, double, double, double); + + void setParamsDCS(const std::vector&); + + double mBg; + double mC; + double mMu_02; + double mAg; + double mLamg; + double mRmax; + double mLmass; + double mCmass; + double mBmass; + double qMass[6]; + double mAs_qqg; + double mRhs; + double mnLP; + double mR0WS; + double mdWS; + double muCE; + double mvCE; + double mR_L2; + double mN_L; + double mN_T; + + double mQ2; + + bool mIsTransverse; + + double uiDIS(double b, double r, double z); + + +private: + // bool setKinematicPoint(double, double, double, double); + double uiIntegralL(double*, double*); + double uiIntegralT(double*, double*); + double dPhi1db(double*, double*); + double dPhi0db(double*, double*); + + double buiA(double*, double*); + double phiuiA(double*, double*); + double ruiA(double*, double*); + double zui(double*, double*); + double rui(double*, double*); + + // double uiPhi0(double*, double*); + // double uiPhi1(double*, double*); + + double kinematicPoint[6]; + double mZ; + double mR; + double mB; + double mBeta; + DipoleModel* mDipoleModel; + int mQuarkIndex; + bool mIsInitialized; + double mInclusiveT[6]; + double mInclusiveL[6]; + double mInclusiveQQG[6]; + double mInclusiveQQG_MS[6]; + double mInclusiveQQG_GBW[6]; + double mInclusiveQQG_GBWbeta0[6]; + unsigned int mA; + unsigned int mNumberOfFlavours; + bool mWhichFlavours[5]; + double mRelativePrecisionOfIntegration; + double mBetaMin; + double mW2; + + DipoleModelParameters* mParameters; + + TF1* mIntegralT; + ROOT::Math::WrappedTF1* mWFT; + ROOT::Math::GaussIntegrator mGIintegralT; + + TF1* mIntegralL; + ROOT::Math::WrappedTF1* mWFL; + ROOT::Math::GaussIntegrator mGIintegralL; + + TF1* mPhi0; + ROOT::Math::WrappedTF1* mWFPhi0; + ROOT::Math::GaussIntegrator mGIPhi0; + + TF1* mPhi1; + ROOT::Math::WrappedTF1* mWFPhi1; + ROOT::Math::GaussIntegrator mGIPhi1; + + TF1* mPhi0r; + ROOT::Math::WrappedTF1* mWFPhi0r; + ROOT::Math::GaussIntegrator mGIPhi0r; + + TF1* mPhi1r; + ROOT::Math::WrappedTF1* mWFPhi1r; + ROOT::Math::GaussIntegrator mGIPhi1r; + + TF1* mRGBW; + ROOT::Math::WrappedTF1* mWFRGBW; + ROOT::Math::GaussIntegrator mGIRGBW; + + TF1* mRGBWbeta0; + ROOT::Math::WrappedTF1* mWFRGBWbeta0; + ROOT::Math::GaussIntegrator mGIRGBWbeta0; + + TF1* mTotalT; + ROOT::Math::WrappedTF1* mWFTotalT; + ROOT::Math::GaussIntegrator mGITotalT; + + TF1* mTotalL; + ROOT::Math::WrappedTF1* mWFTotalL; + ROOT::Math::GaussIntegrator mGITotalL; + + TF1* mSigmaT; + ROOT::Math::WrappedTF1* mWFSigmaT; + ROOT::Math::GaussIntegrator mGISigmaT; + + TF1* mSigmaL; + ROOT::Math::WrappedTF1* mWFSigmaL; + ROOT::Math::GaussIntegrator mGISigmaL; + + TF1* mDipole; + ROOT::Math::WrappedTF1* mWFDipole; + ROOT::Math::GaussIntegrator mGIDipole; + }; -inline DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel() const { return mDipoleModel; } +inline DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel() const { return mDipoleModel; } inline double InclusiveDiffractiveCrossSections::z(){ return mZ; } inline int InclusiveDiffractiveCrossSections::quarkIndex(){ return mQuarkIndex; } inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQ(){ return mInclusiveT; } inline double* InclusiveDiffractiveCrossSections::dsigmadbetaL_QQ(){ return mInclusiveL; } inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG(){ return mInclusiveQQG; } inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG_MS(){ return mInclusiveQQG_MS; } inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG_GBW(){ return mInclusiveQQG_GBW; } inline double* InclusiveDiffractiveCrossSections::dsigmadbetaT_QQG_GBWbeta0(){ return mInclusiveQQG_GBWbeta0; } inline double InclusiveDiffractiveCrossSections::betaMin(){ return mBetaMin; } #endif Index: branches/fitting/src/InclusiveDiffractiveCrossSections_DIS.cpp =================================================================== --- branches/fitting/src/InclusiveDiffractiveCrossSections_DIS.cpp (revision 380) +++ branches/fitting/src/InclusiveDiffractiveCrossSections_DIS.cpp (revision 381) @@ -1,192 +1,176 @@ #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]={3., 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; + const double epsrel=1.e-3, epsabs=1e-12; + const int flags=0, mineval=1e2, 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/IntegrandWrappers.h =================================================================== --- branches/fitting/src/IntegrandWrappers.h (revision 380) +++ branches/fitting/src/IntegrandWrappers.h (revision 381) @@ -1,170 +1,170 @@ //============================================================================== // IntegrandWrappers.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: Tobias Toll // Last update: // $Date$ // $Author$ //============================================================================== // // Wrapper functions meeting the definition of an integrand (see cuba.h): // typedef int (*integrand_t)(const int *ndim, const double x[], // const int *ncomp, double f[], void *userdata); // // Use the void* argument of the integrand to pass a pointer to an // integrals object, which does the actual calculation via its // uiAmplitudeXXX method. // // Cuba only allows integrations with limits 0 and 1, therefore all variables // in the integrations have to be scaled to their actual limits, // and compensated by a Jacobian. // // It saves time to do the integration over four separate functions rather than // letting the integrand being 4D, since in the latter case the integration // continues until all 4 integrals have reached the desired precision, // which is not necessary. // A quick test yields that this way is ~20% faster. // //============================================================================== int integrandWrapperTIm( const int*, const double x[], const int *, double fval[], void* ptr2integrals) { IntegralsExclusive* ip = static_cast( ptr2integrals ); double upperBRValue=2.5*ip->dipoleModel()->nucleus()->radius(); double low[4] = {0, 0, 1e-2, 0}; // b, z, r, phi - lower double high[4] = {upperBRValue, 1, upperBRValue, 2.*M_PI}; // b, z, r, phi - upper double b = low[0] + (high[0] - low[0]) * x[0]; //fm double z = low[1] + (high[1] - low[1]) * x[1]; double r = low[2] + (high[2] - low[2]) * x[2]; //fm double phi = low[3] + (high[3] - low[3]) * x[3]; double jacobian = (high[0]-low[0])*(high[1]-low[1])*(high[2]-low[2])*(high[3]-low[3]); double t=ip->kinematicPoint[0]; double Q2=ip->kinematicPoint[1]; double xprobe=ip->kinematicPoint[3]; double Delta = sqrt(fabs(t)); fval[0]=ip->uiAmplitudeTIm(b, z, r, phi, Q2, xprobe, Delta)*jacobian; return 0; } int integrandWrapperTRe( const int *, const double x[], const int *, double fval[], void* ptr2integrals) { IntegralsExclusive* ip = static_cast( ptr2integrals ); double upperBRValue=2.5*ip->dipoleModel()->nucleus()->radius(); double low[4] = {0, 0, 1e-2, 0}; // b, z, r, phi - lower double high[4] = {upperBRValue, 1, upperBRValue, 2.*M_PI}; // b, z, r, phi - upper double b = low[0] + (high[0] - low[0]) * x[0]; //fm double z = low[1] + (high[1] - low[1]) * x[1]; double r = low[2] + (high[2] - low[2]) * x[2]; //fm double phi = low[3] + (high[3] - low[3]) * x[3]; double jacobian = (high[0]-low[0])*(high[1]-low[1])*(high[2]-low[2])*(high[3]-low[3]); double t=ip->kinematicPoint[0]; double Q2=ip->kinematicPoint[1]; double xprobe=ip->kinematicPoint[3]; double Delta = sqrt(fabs(t)); fval[0]=ip->uiAmplitudeTRe(b, z, r, phi, Q2, xprobe, Delta)*jacobian; return 0; } int integrandWrapperLIm( const int *, const double x[], const int *, double fval[], void* ptr2integrals) { IntegralsExclusive* ip = static_cast( ptr2integrals ); double upperBRValue=2.5*ip->dipoleModel()->nucleus()->radius(); double low[4] = {0, 0, 1e-2, 0}; // b, z, r, phi - lower double high[4] = {upperBRValue, 1, upperBRValue, 2.*M_PI}; // b, z, r, phi - upper double b = low[0] + (high[0] - low[0]) * x[0]; //fm double z = low[1] + (high[1] - low[1]) * x[1]; double r = low[2] + (high[2] - low[2]) * x[2]; //fm double phi = low[3] + (high[3] - low[3]) * x[3]; double jacobian = (high[0]-low[0])*(high[1]-low[1])*(high[2]-low[2])*(high[3]-low[3]); double t=ip->kinematicPoint[0]; double Q2=ip->kinematicPoint[1]; double xprobe=ip->kinematicPoint[3]; double Delta = sqrt(fabs(t)); fval[0]=ip->uiAmplitudeLIm(b, z, r, phi, Q2, xprobe, Delta)*jacobian; return 0; } int integrandWrapperLRe( const int *, const double x[], const int *, double fval[], void* ptr2integrals) { IntegralsExclusive* ip = static_cast( ptr2integrals ); double upperBRValue=2.5*ip->dipoleModel()->nucleus()->radius(); double low[4] = {0, 0, 1e-2, 0}; // b, z, r, phi - lower double high[4] = {upperBRValue, 1, upperBRValue, 2.*M_PI}; // b, z, r, phi - upper double b = low[0] + (high[0] - low[0]) * x[0]; //fm double z = low[1] + (high[1] - low[1]) * x[1]; double r = low[2] + (high[2] - low[2]) * x[2]; //fm double phi = low[3] + (high[3] - low[3]) * x[3]; double jacobian = (high[0]-low[0])*(high[1]-low[1])*(high[2]-low[2])*(high[3]-low[3]); double t=ip->kinematicPoint[0]; double Q2=ip->kinematicPoint[1]; double xprobe=ip->kinematicPoint[3]; double Delta = sqrt(fabs(t)); fval[0]=ip->uiAmplitudeLRe(b, z, r, phi, Q2, xprobe, Delta)*jacobian; return 0; } int integrandWrapperCoherentAmplitudeT( const int *, const double x[], const int *, double fval[], void* ptr2integrals) { IntegralsExclusive* ip = static_cast( ptr2integrals ); double upperBRValue=2.5*ip->dipoleModel()->nucleus()->radius(); double low[5] = {0, 0, 1e-2}; // b, z, r - lower double high[5] = {upperBRValue, 1, upperBRValue}; // b, z, r - upper double b = low[0] + (high[0] - low[0]) * x[0]; //fm double z = low[1] + (high[1] - low[1]) * x[1]; double r = low[2] + (high[2] - low[2]) * x[2]; //fm double jacobian = (high[0]-low[0])*(high[1]-low[1])*(high[2]-low[2]); double t=ip->kinematicPoint[0]; double Q2=ip->kinematicPoint[1]; double Delta = sqrt(fabs(t)); fval[0]=ip->uiCoherentAmplitudeT(b, z, r, Q2, Delta)*jacobian; return 0; } int integrandWrapperCoherentAmplitudeL( const int *, const double x[], const int *, double fval[], void* ptr2integrals) { IntegralsExclusive* ip = static_cast( ptr2integrals ); double upperBRValue=2.5*ip->dipoleModel()->nucleus()->radius(); double low[5] = {0, 0, 1e-2}; // b, z, r - lower double high[5] = {upperBRValue, 1, upperBRValue}; // b, z, r - upper double b = low[0] + (high[0] - low[0]) * x[0]; //fm double z = low[1] + (high[1] - low[1]) * x[1]; double r = low[2] + (high[2] - low[2]) * x[2]; //fm double jacobian = (high[0]-low[0])*(high[1]-low[1])*(high[2]-low[2]); double t=ip->kinematicPoint[0]; double Q2=ip->kinematicPoint[1]; double Delta = sqrt(fabs(t)); fval[0]=ip->uiCoherentAmplitudeL(b, z, r, Q2, Delta)*jacobian; return 0; -} +}