Index: trunk/src/DipoleModelParameters.cpp
===================================================================
--- trunk/src/DipoleModelParameters.cpp (revision 358)
+++ trunk/src/DipoleModelParameters.cpp (revision 359)
@@ -1,468 +1,469 @@
//==============================================================================
// DipoleModelParameters.cpp
//
// Copyright (C) 2016-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: ullrich $
//==============================================================================
#include "DipoleModelParameters.h"
#include "Settings.h"
#include "TableGeneratorSettings.h"
#include
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
DipoleModelParameters::DipoleModelParameters(Settings* settings)
{
mDipoleModelType = settings->dipoleModelType();
mDipoleModelParameterSetName = settings->dipoleModelParameterSetName();
mDipoleModelParameterSet = settings->dipoleModelParameterSet();
setupParameters();
}
DipoleModelParameters::DipoleModelParameters(DipoleModelType mtype, DipoleModelParameterSet pset) :
mDipoleModelType(mtype), mDipoleModelParameterSet(pset)
{
setupParameters();
}
void DipoleModelParameters::setDipoleModelType(DipoleModelType val)
{
mDipoleModelType = val;
setupParameters();
}
void DipoleModelParameters::setDipoleModelParameterSet(DipoleModelParameterSet val)
{
mDipoleModelParameterSet = val;
setupParameters();
}
DipoleModelType DipoleModelParameters::dipoleModelType() const {return mDipoleModelType;}
DipoleModelParameterSet DipoleModelParameters::dipoleModelParameterSet() const {return mDipoleModelParameterSet;}
void DipoleModelParameters::setup_bSat()
{
if (mDipoleModelParameterSet == KMW) {
// KMW paper (arXiv:hep-ph/0606272), Table 3
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
mQuarkMass[3] = 1.4; // c quark
mBG = 4.;
mMu02 = 1.17; // Gev^2
mLambdaG = 0.02;
mAg = 2.55;
mC = 4;
}
else if (mDipoleModelParameterSet == HMPZ) {
// Heikki Mantysaari an Pia Zurita, Phys.Rev. D98 (2018) 036002 (arXiv:1804.05311)
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.03; // u,d,s quarks
mQuarkMass[3] = 1.3528; // c quark
mBG = 4.;
mMu02 = 1.1; // Gev^2
mLambdaG = 0.08289;
mAg = 2.1953;
mC = 2.2894;
}
else if (mDipoleModelParameterSet == CUSTOM) {
if (mCustomParameters.size() < 10) {
cout << "DipoleModelParameters::setupParameters(): Error, require 10 custom parameters for bSAT when" << endl;
cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
exit(1);
}
mQuarkMass[0] = mCustomParameters[0];
mQuarkMass[1] = mCustomParameters[1];
mQuarkMass[2] = mCustomParameters[2];
mQuarkMass[3] = mCustomParameters[3];
mQuarkMass[4] = mCustomParameters[4];
mBG = mCustomParameters[5];
mMu02 = mCustomParameters[6];
mLambdaG = mCustomParameters[7];
mAg = mCustomParameters[8];
mC = mCustomParameters[9];
}
else {
cout << "DipoleModelParameters::setup_bSat(): Error, no known parameters for given dipole model" << endl;
cout << " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
void DipoleModelParameters::setup_bNonSat()
{
if (mDipoleModelParameterSet == KMW) {
// KT paper (arXiv:hep-ph/0304189v3), page 11
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
mQuarkMass[3] = 1.4;
mBG = 4.;
mMu02 = 0.8;
mLambdaG = -0.13;
mAg = 3.5;
mC = 4;
}
else if (mDipoleModelParameterSet == HMPZ) {
// Internal note by Heikki Mantysaari an Pia Zurita, arXiv pending
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.1516; // u,d,s quarks
mQuarkMass[3] = 1.3504;
mBG = 4.;
mMu02 = 1.1;
mLambdaG = -0.006657;
mAg = 3.0391;
mC = 4.2974;
}
else if (mDipoleModelParameterSet == CUSTOM) {
if (mCustomParameters.size() < 10) {
cout << "DipoleModelParameters::setupParameters(): Error, require 10 custom parameters for bNonSAT when" << endl;
cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
exit(1);
}
mQuarkMass[0] = mCustomParameters[0];
mQuarkMass[1] = mCustomParameters[1];
mQuarkMass[2] = mCustomParameters[2];
mQuarkMass[3] = mCustomParameters[3];
mQuarkMass[4] = mCustomParameters[4];
mBG = mCustomParameters[5];
mMu02 = mCustomParameters[6];
mLambdaG = mCustomParameters[7];
mAg = mCustomParameters[8];
mC = mCustomParameters[9];
}
else {
cout << "DipoleModelParameters::setup_bNonSat(): Error, no known parameters for given dipole model" << endl;
cout << " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
void DipoleModelParameters::setup_bCGC()
{
if (mDipoleModelParameterSet == KMW) {
// WK paper (arXiv:0712.2670), Table II
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
mQuarkMass[3] = 1.4;
mKappa = 9.9;
mN0 = 0.558;
mX0 = 1.84e-6;
mLambda = 0.119;
mGammas = 0.46;
mBcgc = 7.5;
}
else if (mDipoleModelParameterSet == CUSTOM) {
if (mCustomParameters.size() < 10) {
cout << "DipoleModelParameters::setup_bCGC(): Error, require 10 custom parameters for bCGC when" << endl;
cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
exit(1);
}
mQuarkMass[0] = mCustomParameters[0];
mQuarkMass[1] = mCustomParameters[1];
mQuarkMass[2] = mCustomParameters[2];
mQuarkMass[3] = mCustomParameters[3];
mKappa = mCustomParameters[4];
mN0 = mCustomParameters[5];
mX0 = mCustomParameters[6];
mLambda = mCustomParameters[7];
mGammas = mCustomParameters[8];
mBcgc = mCustomParameters[9];
}
else {
cout << "DipoleModelParameters::setup_bCGC(): Error, no known parameters for given dipole model"
<< " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
void DipoleModelParameters::setupParameters()
{
TableGeneratorSettings *settings = TableGeneratorSettings::instance();
if (mDipoleModelParameterSet == CUSTOM) mCustomParameters = settings->dipoleModelCustomParameters();
//
// Init
//
mKappa = 0;
mN0 = 0;
mX0 = 0;
mLambda = 0;
mGammas = 0;
mBcgc = 0;
mBG = 0;
mMu02 = 0;
mLambdaG = 0;
mAg = 0;
mC = 0;
mBoostedGaussianR2_rho = 0;
mBoostedGaussianNL_rho = 0;
mBoostedGaussianNT_rho = 0;
mBoostedGaussianQuarkMass_rho = 0;
mBoostedGaussianR2_phi = 0;
mBoostedGaussianNL_phi = 0;
mBoostedGaussianNT_phi = 0;
mBoostedGaussianQuarkMass_phi = 0;
mBoostedGaussianR2_jpsi = 0;
mBoostedGaussianNL_jpsi = 0;
mBoostedGaussianNT_jpsi = 0;
mBoostedGaussianQuarkMass_jpsi = 0;
mBoostedGaussianR2_ups = 0;
mBoostedGaussianNL_ups = 0;
mBoostedGaussianNT_ups = 0;
mBoostedGaussianQuarkMass_ups = 0;
//
// b and t masses (not used, just for completeness)
//
- mQuarkMass[4] = 4.75; // b quark consistent with HMPZ
+ // mQuarkMass[4] = 4.75; // b quark consistent with HMPZ
+ mQuarkMass[4] = 4.2; // b quark consistent with Upsilon wave function.
mQuarkMass[5] = 175.; // t quark consistent with HMPZ
//
// Parameters for boosted Gaussian wave function
//
setup_boostedGaussiansWaveFunction();
//
// Model parameters
//
if (mDipoleModelType == bSat) {
setup_bSat();
}
else if (mDipoleModelType == bNonSat) {
setup_bNonSat();
}
else if (mDipoleModelType == bCGC) {
setup_bCGC();
}
else {
cout << "DipoleModelParameters::setupParameters(): Error, no known parameters for given dipole model" << endl;
cout << " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianR2(int vm)
{
if (vm == 113)
return mBoostedGaussianR2_rho;
else if (vm == 333)
return mBoostedGaussianR2_phi;
else if (vm == 443)
return mBoostedGaussianR2_jpsi;
else if (vm == 553)
return mBoostedGaussianR2_ups;
else {
cout << "DipoleModelParameters::boostedGaussianR2(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianNL(int vm)
{
if (vm == 113)
return mBoostedGaussianNL_rho;
else if (vm == 333)
return mBoostedGaussianNL_phi;
else if (vm == 443)
return mBoostedGaussianNL_jpsi;
else if (vm == 553)
return mBoostedGaussianNL_ups;
else {
cout << "DipoleModelParameters::boostedGaussianNL(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianNT(int vm)
{
if (vm == 113)
return mBoostedGaussianNT_rho;
else if (vm == 333)
return mBoostedGaussianNT_phi;
else if (vm == 443)
return mBoostedGaussianNT_jpsi;
else if (vm == 553)
return mBoostedGaussianNT_ups;
else {
cout << "DipoleModelParameters::boostedGaussianNT(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
void DipoleModelParameters::setup_boostedGaussiansWaveFunction()
{
//
// Technical note:
// The Upsilon is a late addition with parameters coming from
// DKMM (arXiv:hep-ph/1610.06647). Heikki provided the more
// precise normalization constants for N_T and N_L.
// Neither KMW nor HMPZ provided Upsilon parameters, so
// results need to be verified with HERA data first.
- //
+ // Note: All that is taken from DKMM is b-quark mass, the rest
+ // is determined by norm and decay width.
if (mDipoleModelParameterSet == KMW || mDipoleModelParameterSet == HMPZ) {
mBoostedGaussianR2_ups = 0.567;
- mBoostedGaussianNT_ups = 0.481493;
- mBoostedGaussianNL_ups = 0.480264 ; // from Heikki, not provided in DKMM
- mBoostedGaussianQuarkMass_ups = 4.2;
+ mBoostedGaussianNT_ups = 0.481493;
+ mBoostedGaussianNL_ups = 0.480264 ;
+ mBoostedGaussianQuarkMass_ups = 4.2;
}
-
if (mDipoleModelParameterSet == KMW) {
//
// KMW: bSat, bNonSat, and bCGC use the same parameters
// and also do not distingiosh between T and L.
//
mBoostedGaussianR2_rho = 12.9;
mBoostedGaussianNL_rho = 0.853;
mBoostedGaussianNT_rho = 0.911;
mBoostedGaussianQuarkMass_rho = 0.14;
mBoostedGaussianR2_phi = 11.2;
mBoostedGaussianNL_phi = 0.825;
mBoostedGaussianNT_phi= 0.919;
mBoostedGaussianQuarkMass_phi = 0.14;
mBoostedGaussianR2_jpsi = 2.3;
mBoostedGaussianNL_jpsi = 0.575;
mBoostedGaussianNT_jpsi = 0.578;
mBoostedGaussianQuarkMass_jpsi = 1.4;
}
else if (mDipoleModelParameterSet == HMPZ) {
if (mDipoleModelType == bSat) {
mBoostedGaussianR2_rho = 3.6376*3.6376;
mBoostedGaussianNL_rho = 0.8926;
mBoostedGaussianNT_rho = 0.9942;
mBoostedGaussianQuarkMass_rho = 0.03;
mBoostedGaussianR2_phi = 3.3922*3.3922;
mBoostedGaussianNL_phi = 0.8400;
mBoostedGaussianNT_phi = 0.9950;
mBoostedGaussianQuarkMass_phi = 0.03;
mBoostedGaussianR2_jpsi = 1.5070*1.5070;
mBoostedGaussianNL_jpsi = 0.5860;
mBoostedGaussianNT_jpsi = 0.5890;
mBoostedGaussianQuarkMass_jpsi = 1.3528;
}
else if (mDipoleModelType == bNonSat) {
mBoostedGaussianR2_rho = 3.5750*3.5750;
mBoostedGaussianNL_rho = 0.8467;
mBoostedGaussianNT_rho = 0.8978;
mBoostedGaussianQuarkMass_rho = 0.1516;
mBoostedGaussianR2_phi = 3.3530*3.3530;
mBoostedGaussianNL_phi = 0.8196;
mBoostedGaussianNT_phi = 0.9072;
mBoostedGaussianQuarkMass_phi = 0.1516;
mBoostedGaussianR2_jpsi = 1.5071*1.5071;
mBoostedGaussianNL_jpsi = 0.5868;
mBoostedGaussianNT_jpsi = 0.5899;
mBoostedGaussianQuarkMass_jpsi = 1.3504;
}
else if (mDipoleModelType == bCGC) {
cout << "DipoleModelParameters::setup_boostedGaussiansWaveFunction(): "
"Error, no HMPZ wave function parameters for CGC model. Stop." << endl;
exit(1);
}
}
else {
cout << "DipoleModelParameters::setup_boostedGaussiansWaveFunction(): Error, no known parameters for given dipole model "
<< " parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianQuarkMass(int vm)
{
if (vm == 113)
return mBoostedGaussianQuarkMass_rho;
else if (vm == 333)
return mBoostedGaussianQuarkMass_phi;
else if (vm == 443)
return mBoostedGaussianQuarkMass_jpsi;
else if (vm == 553)
return mBoostedGaussianQuarkMass_ups;
else {
cout << "DipoleModelParameters::boostedGaussianQuarkMass(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
bool DipoleModelParameters::list(ostream& os)
{
const int fieldWidth = 32;
os << "\nDipole Model Parameters:" << endl;
os << setw(fieldWidth) << "Set: " << mDipoleModelParameterSetName << endl;
os << setw(fieldWidth) << "Quark masses: "
<< "u=" << mQuarkMass[0]
<< ", d=" << mQuarkMass[1]
<< ", s=" << mQuarkMass[2]
<< ", c=" << mQuarkMass[3]
<< ", b=" << mQuarkMass[4]
<< ", t=" << mQuarkMass[5] << endl;
os << setw(fieldWidth) << "BG: " << mBG << endl;
os << setw(fieldWidth) << "Mu02: " << mMu02 << endl;
os << setw(fieldWidth) << "LambdaG: " << mLambdaG << endl;
os << setw(fieldWidth) << "Ag: " << mAg << endl;
os << setw(fieldWidth) << "C: " << mC << endl;
os << setw(fieldWidth) << "Kappa: " << mKappa << endl;
os << setw(fieldWidth) << "N0: " << mN0 << endl;
os << setw(fieldWidth) << "X0: " << mX0 << endl;
os << setw(fieldWidth) << "Lambda: " << mLambda << endl;
os << setw(fieldWidth) << "Gammas: " << mGammas << endl;
os << setw(fieldWidth) << "Bcgc: " << mBcgc << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_rho: " << mBoostedGaussianR2_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_rho: " << mBoostedGaussianNL_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_rho: " << mBoostedGaussianNT_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_rho: " << mBoostedGaussianQuarkMass_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_phi: " << mBoostedGaussianR2_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_phi: " << mBoostedGaussianNL_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_phi: " << mBoostedGaussianNT_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_phi: " << mBoostedGaussianQuarkMass_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_jpsi: " << mBoostedGaussianR2_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_jpsi: " << mBoostedGaussianNL_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_jpsi: " << mBoostedGaussianNT_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_jpsi: " << mBoostedGaussianQuarkMass_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_ups: " << mBoostedGaussianR2_ups << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_ups: " << mBoostedGaussianNL_ups << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_ups: " << mBoostedGaussianNT_ups << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_ups: " << mBoostedGaussianQuarkMass_ups << endl;
os << endl;
return true;
}
Index: trunk/src/WaveOverlap.h
===================================================================
--- trunk/src/WaveOverlap.h (revision 358)
+++ trunk/src/WaveOverlap.h (revision 359)
@@ -1,75 +1,80 @@
//==============================================================================
// WaveOverlap.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$
//==============================================================================
#ifndef WaveOverlap_h
#define WaveOverlap_h
class DipoleModelParameters;
class WaveOverlap {
public:
WaveOverlap();
virtual ~WaveOverlap();
virtual void setWaveOverlapFunctionParameters(int);
virtual double T(double, double, double)=0;
virtual double L(double, double, double)=0;
virtual void setProcess(int){};
-
+ virtual void testBoostedGaussianParameters(int){};
+
protected:
DipoleModelParameters *mParameters;
};
class WaveOverlapVM : public WaveOverlap {
public:
WaveOverlapVM();
void setWaveOverlapFunctionParameters(int);
void setProcess(int);
-
+ void testBoostedGaussianParameters(int);
+
private:
double T(double, double, double);
double L(double, double, double);
double transverseWaveFunction(double, double);
double longitudinalWaveFunction(double, double);
double dDrTransverseWaveFunction(double, double);
double laplaceRLongitudinalWaveFunction(double, double);
-
+ double uiDecayWidth(double*, double*);
+ double uiNormL(const double*);
+ double uiNormT(const double*);
+
private:
double mNT, mRT2;
double mMf; // mass of quarks in vector meson
double mBoostedGaussianMf; // // mass of quarks in vector meson's boosted Gaussian wave fct.
double mMf2;
double mBoostedGaussianMf2;
double mEf;
double mMV;
double mNL, mRL2;
};
class WaveOverlapDVCS : public WaveOverlap {
private:
double T(double, double, double);
double L(double, double, double);
};
inline void WaveOverlap::setWaveOverlapFunctionParameters(int) {/* no op*/}
#endif
Index: trunk/src/Integrals.cpp
===================================================================
--- trunk/src/Integrals.cpp (revision 358)
+++ trunk/src/Integrals.cpp (revision 359)
@@ -1,737 +1,740 @@
//==============================================================================
// 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;
mDipoleModelForSkewednessCorrection = 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;
mIntegralTForSkewedness = 0;
mIntegralLForSkewedness = 0;
mErrorTForSkewedness = 0;
mErrorLForSkewedness = 0;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
+ mVerbose=settings->verbose();
+
int VMId=settings->vectorMesonId();
mMV = settings->lookupPDG(VMId)->Mass();
mIsUPC = settings->UPC();
if (VMId==113 || VMId==333 || VMId == 443 || VMId == 553) {
mWaveOverlap = new WaveOverlapVM;
mWaveOverlap->setProcess(VMId);
mWaveOverlap->setWaveOverlapFunctionParameters(VMId);
+ if(mVerbose)
+ mWaveOverlap->testBoostedGaussianParameters(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 {
cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
exit(1);
}
mCalculateSkewedness=false;
if(settings->A()==1 and settings->modesToCalculate()!=1 and settings->numberOfConfigurations()==1 and model==bSat) {
mCalculateSkewedness=true;
mDipoleModelForSkewednessCorrection = new DipoleModel_bNonSat;
}
- 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;
if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mErrorImT = integrals.mErrorImT;
mErrorImL = integrals.mErrorImL;
mErrorReT = integrals.mErrorReT;
mErrorReL = integrals.mErrorReL;
mProbImT = integrals.mProbImT;
mProbImL = integrals.mProbImL;
mProbReT = integrals.mProbReT;
mProbReL = integrals.mProbReL;
mMV = integrals.mMV;
}
Integrals& Integrals::operator=(const Integrals& integrals)
{
if (this != &integrals) {
delete mWaveOverlap;
delete mDipoleModel;
delete mDipoleModelForSkewednessCorrection;
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;
if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mErrorImT = integrals.mErrorImT;
mErrorImL = integrals.mErrorImL;
mErrorReT = integrals.mErrorReT;
mErrorReL = integrals.mErrorReL;
mProbImT = integrals.mProbImT;
mProbImL = integrals.mProbImL;
mProbReT = integrals.mProbReT;
mProbReL = integrals.mProbReL;
mMV = integrals.mMV;
mIsInitialized = integrals.mIsInitialized;
mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
}
return *this;
}
Integrals::~Integrals()
{
delete mWaveOverlap;
delete mDipoleModel;
if(mDipoleModelForSkewednessCorrection)
delete mDipoleModelForSkewednessCorrection;
}
void Integrals::operator() (double t, double Q2, double W2)
{
unsigned int A=dipoleModel()->nucleus()->A();
//make sure the configurations have been generated:
if (!mDipoleModel->configurationExists() and A!=1) {
// do not use cout
cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
exit(1);
}
if (setKinematicPoint(t, Q2, W2)) {
if(A==1){
calculateEp();
if(mCalculateSkewedness){
calculateSkewedness();
}
}
else
calculate();
}
else {
fillZeroes();
}
}
void Integrals::operator() (double t, double xpom) //UPC
{
unsigned int A=dipoleModel()->nucleus()->A();
//make sure the configurations have been generated:
if (!mDipoleModel->configurationExists() and A!=1) {
// do not use cout
cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
exit(1);
}
if (setKinematicPoint(t, xpom)) {
if(A==1){
calculateEp();
if(mCalculateSkewedness){
calculateSkewedness();
}
}
else
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; //Q2
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;
//
// 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);
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
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);
}
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);
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
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);
}
//
// 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::calculateEp()
{
//
// 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-4, epsabs=1e-12;
const int flags=0, maxeval=1e9, key=0;
const int mineval=3e6;
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;
//
// Do the integrations
//
/*
bool bContinue=true;
bool isFirst=true;
while(bContinue){
double valTOld=valT;
Cuhre(4, 1, integrandWrapperTep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(abs(valT-valTOld)/valT > epsrel){
mineval*=3;
if(isFirst)
mineval*=30;
isFirst=false;
if(mineval>1e4)
PR(mineval);
}
else
bContinue=false;
}
*/
Cuhre(4, 1, integrandWrapperTep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(failT!=0 and mVerbose)
printf("IntegralsExclusive::calculateEp(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
/*
if(errT/valT > epsrel)
PR(errT/valT);
if(errT < epsabs)
PR(errT);
if(probT>0.5)
PR(probT);
PR(nevalT);
PR(nregionsT);
*/
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
Cuhre(4, 1, integrandWrapperLep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
if(failL!=0 and mVerbose)
printf("IntegralsExclusive::calculateEp(): 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;
//
// Store the probabilities:
//
mProbImT=probT;
mProbImL=probL;
}
void IntegralsExclusive::calculateSkewedness()
{
//
// 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-4, epsabs=1e-12;
const int flags=0, maxeval=1e9, key=0;
const int mineval=3e6;
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;
//
// Do the integrations
//
/*
bool bContinue=true;
bool isFirst=true;
while(bContinue){
double valTOld=valT;
Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(abs(valT-valTOld)/valT > epsrel){
mineval*=3;
if(isFirst)
mineval*=30;
isFirst=false;
if(mineval>1e4)
PR(mineval);
}
else
bContinue=false;
}
*/
Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(failT!=0 and mVerbose)
printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
Cuhre(4, 1, integrandWrapperLForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
if(failL!=0 and mVerbose)
printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
}
//
// Store the results:
//
mIntegralTForSkewedness=valT;
mIntegralLForSkewedness=valL;
//
// Store the errors:
//
mErrorTForSkewedness=errT;
mErrorLForSkewedness=errL;
//
// Store the probabilities:
//
mProbImTForSkewedness=probT;
mProbImLForSkewedness=probL;
}
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;
}
double IntegralsExclusive::uiAmplitudeTep(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
double IntegralsExclusive::uiAmplitudeLep(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
//
// Only for calculating the lamdba for Skewedness Corrections:
//
double IntegralsExclusive::uiAmplitudeTForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
double IntegralsExclusive::uiAmplitudeLForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r, b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
Index: trunk/src/WaveOverlap.cpp
===================================================================
--- trunk/src/WaveOverlap.cpp (revision 358)
+++ trunk/src/WaveOverlap.cpp (revision 359)
@@ -1,186 +1,270 @@
//==============================================================================
// WaveOverlap.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 "WaveOverlap.h"
#include "Constants.h"
#include "TableGeneratorSettings.h"
#include "DipoleModelParameters.h"
#include "TMath.h"
#include
#include
+#include "TF1.h"
+#include "Math/Functor.h"
+#include "Math/IntegratorMultiDim.h"
+#include "Math/WrappedTF1.h"
+#include "Math/GaussIntegrator.h"
+
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
WaveOverlap::WaveOverlap()
{
mParameters = new DipoleModelParameters(TableGeneratorSettings::instance());
}
WaveOverlap::~WaveOverlap() {/* no op*/}
//
// VECTOR MESONS
//
WaveOverlapVM::WaveOverlapVM()
{
mNT = mRT2 = 0;
mMf = 0;
mMf2 = 0;
mEf = 0;
mMV = 0;
mNL = mRL2 = 0;
mBoostedGaussianMf = 0;
-}
+}
+double WaveOverlapVM::uiNormT(const double* var){
+ double z=var[0];
+ double r=var[1];
+ double phi=transverseWaveFunction(r, z); //GeV0
+ double dphidr=dDrTransverseWaveFunction(r, z); //GeV1
+ return 3.*r/hbarc/(z*z*(1-z)*(1-z))
+ *(mMf2*phi*phi + (z*z + (1-z)*(1-z))*dphidr*dphidr); //GeV1
+}
+
+double WaveOverlapVM::uiNormL(const double* var){
+ double z=var[0];
+ double r=var[1];
+ double phi=longitudinalWaveFunction(r, z); //GeV0
+ double d2phidr2=laplaceRLongitudinalWaveFunction(r, z); //GeV2
+ double term=mMV*phi + (mMf2*phi - d2phidr2)/(mMV*z*(1-z)); //GeV1
+ return 3.*r/hbarc*term*term; //GeV1
+}
+
+double WaveOverlapVM::uiDecayWidth(double* var, double* par)
+{
+ double z=*var;
+ double phi=longitudinalWaveFunction(0., z); //GeV0
+ double d2phidr2=laplaceRLongitudinalWaveFunction(0., z); //GeV0
+ double result = mEf*3./M_PI*(mMV*phi+(mMf2*phi-d2phidr2)/(mMV*z*(1-z))); //GeV1
+ return result;
+}
+
+void WaveOverlapVM::testBoostedGaussianParameters(int id)
+{
+ //
+ // This function calculates the resulting normalisation
+ // and decay width resulting from the boosted gaussian parameters
+ // and compare with the actual values. This can be used to
+ // test or modify the boosted gaussian parameters.
+ //
+ // KMW paper hep-ph/0606272 Eqs. (24)-(28)
+ //
+ // Start with decay width:
+ TF1 fDW("fDW", this, &WaveOverlapVM::uiDecayWidth, 0., 1., 0.);
+ ROOT::Math::WrappedTF1 wfDW(fDW);
+ ROOT::Math::GaussIntegrator giDW;
+ giDW.SetFunction(wfDW);
+ giDW.SetAbsTolerance(0.);
+ giDW.SetRelTolerance(1e-5);
+ double f_VL=giDW.Integral(0., 1.); //GeV
+ PR(f_VL);
+ cout<<"The e+e- decay width is: "<ee)="<<4*M_PI*alpha_em*alpha_em*f_VL*f_VL/3/mMV*1e6
+ <<" keV"<ee)=1.340 +- 0.018 keV"<ee)=5.55 +- 0.14 +- 0.02 keV"<ee)=1.27 +- 0.04 keV"<ee)=7.04 +- 0.06 keV"<boostedGaussianR2(val);
mRT2 = mRL2;
mNT = mParameters->boostedGaussianNT(val);
mNL = mParameters->boostedGaussianNL(val);
mBoostedGaussianMf = mParameters->boostedGaussianQuarkMass(val);
mBoostedGaussianMf2 = mBoostedGaussianMf*mBoostedGaussianMf;
}
void WaveOverlapVM::setProcess(int val)
{
switch (val) {
case 113:
mMf = mParameters->quarkMass(1);
mMV = 0.776;
mEf = 1./sqrt(2.);
break;
case 333:
mMf = mParameters->quarkMass(2);
mMV = 1.019;
mEf = 1./3.;
break;
case 443:
mMf = mParameters->quarkMass(3);
mMV = 3.096916;
mEf = 2./3.;
break;
case 553:
mMf = mParameters->quarkMass(4);
mMV = 9.46;
mEf = -1./3.;
break;
default:
cerr << "WaveOverlap::setProcess(): error no such type: " << val << endl;
break;
}
mMf2 = mMf*mMf;
}
//
// DVCS
//
double WaveOverlapDVCS::T(double z, double Q2, double r)
{
// KMW paper hep-ph/0606272 Eq. 17
double term0, term1, term2;
double result=0;
for (int iFlav=0; iFlav<4; iFlav++) {
double mf=mParameters->quarkMass(iFlav);
double ef=quarkCharge[iFlav];
double eps2 = z*(1-z)*Q2 + mf*mf;
double eps = sqrt(eps2);
term0=2.*Nc/M_PI*alpha_em*ef*ef;
term1=( z*z+(1-z)*(1-z) )*
eps*TMath::BesselK1(eps*r/hbarc)*mf*TMath::BesselK1(mf*r/hbarc);
term2=mf*mf*TMath::BesselK0(eps*r/hbarc)*TMath::BesselK0(mf*r/hbarc);
result += term0*(term1+term2);
}
return result;
}
double WaveOverlapDVCS::L(double, double, double) {return 0;}