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