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;
-}
+}