Index: trunk/src/DipoleModelParameters.h
===================================================================
--- trunk/src/DipoleModelParameters.h (revision 342)
+++ trunk/src/DipoleModelParameters.h (revision 343)
@@ -1,133 +1,137 @@
//==============================================================================
// DipoleModelParameters.h
//
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// $Date$
// $Author: ullrich $
//==============================================================================
#ifndef DipoleModelParameters_h
#define DipoleModelParameters_h
#include "Enumerations.h"
#include <iostream>
#include <vector>
using namespace std;
class Settings;
class DipoleModelParameters {
public:
DipoleModelParameters(Settings*);
DipoleModelParameters(DipoleModelType, DipoleModelParameterSet);
void setDipoleModelType(DipoleModelType);
void setDipoleModelParameterSet(DipoleModelParameterSet);
DipoleModelType dipoleModelType() const;
DipoleModelParameterSet dipoleModelParameterSet() const;
// bSat, bNonSat
double BG() const;
double mu02() const;
double lambdaG() const;
double Ag() const;
double C() const;
// bCGC
double kappa() const;
double N0() const;
double x0() const;
double lambda() const;
double gammas() const;
double Bcgc() const;
double quarkMass(unsigned int) const;
double boostedGaussianR2(int vmID);
double boostedGaussianNL(int vmID);
double boostedGaussianNT(int vmID);
double boostedGaussianQuarkMass(int vmID);
bool list(ostream& = cout);
private:
void setupParameters();
void setup_bSat();
void setup_bNonSat();
void setup_bCGC();
void setup_boostedGaussiansWaveFunction();
private:
DipoleModelType mDipoleModelType;
string mDipoleModelParameterSetName;
DipoleModelParameterSet mDipoleModelParameterSet;
double mQuarkMass[6];
// bSat, bNonSat
double mBG;
double mMu02;
double mLambdaG;
double mAg;
double mC;
// bCGC
double mKappa;
double mN0;
double mX0;
double mLambda;
double mGammas;
double mBcgc;
// boosted Gaussian wave function parameters
double mBoostedGaussianR2_rho;
double mBoostedGaussianNL_rho;
double mBoostedGaussianNT_rho;
double mBoostedGaussianQuarkMass_rho;
double mBoostedGaussianR2_phi;
double mBoostedGaussianNL_phi;
double mBoostedGaussianNT_phi;
double mBoostedGaussianQuarkMass_phi;
double mBoostedGaussianR2_jpsi;
double mBoostedGaussianNL_jpsi;
double mBoostedGaussianNT_jpsi;
double mBoostedGaussianQuarkMass_jpsi;
-
+ double mBoostedGaussianR2_ups;
+ double mBoostedGaussianNL_ups;
+ double mBoostedGaussianNT_ups;
+ double mBoostedGaussianQuarkMass_ups;
+
// hold the custom parameter (internal only)
vector<double> mCustomParameters;
};
inline double DipoleModelParameters::BG() const {return mBG;}
inline double DipoleModelParameters::mu02() const {return mMu02;}
inline double DipoleModelParameters::C() const {return mC;}
inline double DipoleModelParameters::lambdaG() const {return mLambdaG;}
inline double DipoleModelParameters::Ag() const {return mAg;}
inline double DipoleModelParameters::kappa() const {return mKappa;}
inline double DipoleModelParameters::N0() const {return mN0;}
inline double DipoleModelParameters::x0() const {return mX0;}
inline double DipoleModelParameters::lambda() const {return mLambda;}
inline double DipoleModelParameters::gammas() const {return mGammas;}
inline double DipoleModelParameters::Bcgc() const {return mBcgc;}
inline double DipoleModelParameters::quarkMass(unsigned int i) const
{
if (i < 6)
return mQuarkMass[i];
else
return 0;
}
#endif
Index: trunk/src/Integrals.cpp
===================================================================
--- trunk/src/Integrals.cpp (revision 342)
+++ trunk/src/Integrals.cpp (revision 343)
@@ -1,737 +1,737 @@
//==============================================================================
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include <iostream>
#include <cmath>
#include <algorithm>
#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();
int VMId=settings->vectorMesonId();
mMV = settings->lookupPDG(VMId)->Mass();
mIsUPC = settings->UPC();
- if (VMId==113 || VMId==333 || VMId == 443) {
+ if (VMId==113 || VMId==333 || VMId == 443 || VMId == 553) {
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 {
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 342)
+++ trunk/src/WaveOverlap.cpp (revision 343)
@@ -1,181 +1,186 @@
//==============================================================================
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "WaveOverlap.h"
#include "Constants.h"
#include "TableGeneratorSettings.h"
#include "DipoleModelParameters.h"
#include "TMath.h"
#include <iostream>
#include <cmath>
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::transverseWaveFunction(double r, double z)
{
// KMW paper hep-ph/0606272 Eq. 33
return mNT*z*(1-z)*exp(-(mBoostedGaussianMf2*mRT2)/(8*z*(1-z)) -
(2*z*(1-z)*r*r)/mRT2/hbarc2 + (mBoostedGaussianMf2*mRT2)/2);
}
double WaveOverlapVM::dDrTransverseWaveFunction(double r, double z)
{
return transverseWaveFunction(r,z) * (-4*z*(1-z)*r/mRT2/hbarc);
}
double WaveOverlapVM::T(double z, double Q2, double r)
{
// KMW paper hep-ph/0606272 Eq. 21
// Units:
// Q2 in GeV^2
// r in fm
const double e = sqrt(4*M_PI*alpha_em);
double eps2 = z*(1-z)*Q2 + mMf2;
double eps = sqrt(eps2);
double term0 = mEf*e*(Nc/(M_PI*z*(1-z)));
double term1 = mMf2*TMath::BesselK0(r*eps/hbarc)*transverseWaveFunction(r,z);
double term2 = (z*z + (1-z)*(1-z))*eps*TMath::BesselK1(r*eps/hbarc)*dDrTransverseWaveFunction(r,z);
return term0*(term1-term2);
}
double WaveOverlapVM::L(double z, double Q2, double r)
{
// KMW paper hep-ph/0606272 Eq. 22
// Units:
// Q2 in GeV^2
// r in fm
const double e = sqrt(4*M_PI*alpha_em);
double eps2 = z*(1-z)*Q2 + mMf2;
double eps = sqrt(eps2);
double result = mEf*e*(Nc/M_PI)*2*sqrt(Q2)*z*(1-z);
result *= TMath::BesselK0(r*eps/hbarc);
double term1 = mMV*longitudinalWaveFunction(r,z);
double term2 = mMf2*longitudinalWaveFunction(r,z) - laplaceRLongitudinalWaveFunction(r,z);
term2 /= (mMV*z*(1-z));
return result*(term1+term2);
}
double WaveOverlapVM::longitudinalWaveFunction(double r, double z)
{
// KMW paper hep-ph/0606272 Eq. 33
return mNL*z*(1-z)*exp(-(mBoostedGaussianMf2*mRL2)/(8*z*(1-z)) -
(2*z*(1-z)*r*r)/mRL2/hbarc2 + (mBoostedGaussianMf2*mRL2)/2);
}
double WaveOverlapVM::laplaceRLongitudinalWaveFunction(double r, double z)
{
double t = 4*z*(1-z)/mRL2;
return longitudinalWaveFunction(r,z)* (r*r*t*t/hbarc2 - 2*t);
}
void WaveOverlapVM::setWaveOverlapFunctionParameters(int val)
{
// KMW paper hep-ph/0606272 Table II
//
// mRL2 (GeV^-2)
//
mRL2 = mParameters->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;}
Index: trunk/src/Settings.cpp
===================================================================
--- trunk/src/Settings.cpp (revision 342)
+++ trunk/src/Settings.cpp (revision 343)
@@ -1,580 +1,580 @@
//==============================================================================
// Settings.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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "Settings.h"
#include <typeinfo>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <ctype.h>
#include "TParticlePDG.h"
#include <cstdlib>
#define PR(x) cout << #x << " = " << (x) << endl;
TRandom3 Settings::mRandomGenerator;
Settings::Settings()
{
mPDG = new TDatabasePDG;
mPDG->ReadPDGTable();
//
// Setup name table (map) for nuclei
//
mPeriodicTable[1] = "H";
mPeriodicTable[2] = "He";
mPeriodicTable[3] = "Li";
mPeriodicTable[4] = "Be";
mPeriodicTable[5] = "B";
mPeriodicTable[6] = "C";
mPeriodicTable[7] = "N";
mPeriodicTable[8] = "O";
mPeriodicTable[9] = "F";
mPeriodicTable[10] = "Ne";
mPeriodicTable[11] = "Na";
mPeriodicTable[12] = "Mg";
mPeriodicTable[13] = "Al";
mPeriodicTable[14] = "Si";
mPeriodicTable[15] = "P";
mPeriodicTable[16] = "S";
mPeriodicTable[17] = "Cl";
mPeriodicTable[18] = "Ar";
mPeriodicTable[19] = "K";
mPeriodicTable[20] = "Ca";
mPeriodicTable[21] = "Sc";
mPeriodicTable[22] = "Ti";
mPeriodicTable[23] = "V";
mPeriodicTable[24] = "Cr";
mPeriodicTable[25] = "Mn";
mPeriodicTable[26] = "Fe";
mPeriodicTable[27] = "Co";
mPeriodicTable[28] = "Ni";
mPeriodicTable[29] = "Cu";
mPeriodicTable[30] = "Zn";
mPeriodicTable[31] = "Ga";
mPeriodicTable[32] = "Ge";
mPeriodicTable[33] = "As";
mPeriodicTable[34] = "Se";
mPeriodicTable[35] = "Br";
mPeriodicTable[36] = "Kr";
mPeriodicTable[37] = "Rb";
mPeriodicTable[38] = "Sr";
mPeriodicTable[39] = "Y";
mPeriodicTable[40] = "Zr";
mPeriodicTable[41] = "Nb";
mPeriodicTable[42] = "Mo";
mPeriodicTable[43] = "Tc";
mPeriodicTable[44] = "Ru";
mPeriodicTable[45] = "Rh";
mPeriodicTable[46] = "Pd";
mPeriodicTable[47] = "Ag";
mPeriodicTable[48] = "Cd";
mPeriodicTable[49] = "In";
mPeriodicTable[50] = "Sn";
mPeriodicTable[51] = "Sb";
mPeriodicTable[52] = "Te";
mPeriodicTable[53] = "I";
mPeriodicTable[54] = "Xe";
mPeriodicTable[55] = "Cs";
mPeriodicTable[56] = "Ba";
mPeriodicTable[57] = "La";
mPeriodicTable[58] = "Ce";
mPeriodicTable[59] = "Pr";
mPeriodicTable[60] = "Nd";
mPeriodicTable[61] = "Pm";
mPeriodicTable[62] = "Sm";
mPeriodicTable[63] = "Eu";
mPeriodicTable[64] = "Gd";
mPeriodicTable[65] = "Tb";
mPeriodicTable[66] = "Dy";
mPeriodicTable[67] = "Ho";
mPeriodicTable[68] = "Er";
mPeriodicTable[69] = "Tm";
mPeriodicTable[70] = "Yb";
mPeriodicTable[71] = "Lu";
mPeriodicTable[72] = "Hf";
mPeriodicTable[73] = "Ta";
mPeriodicTable[74] = "W";
mPeriodicTable[75] = "Re";
mPeriodicTable[76] = "Os";
mPeriodicTable[77] = "Ir";
mPeriodicTable[78] = "Pt";
mPeriodicTable[79] = "Au";
mPeriodicTable[80] = "Hg";
mPeriodicTable[81] = "Tl";
mPeriodicTable[82] = "Pb";
mPeriodicTable[83] = "Bi";
mPeriodicTable[84] = "Po";
mPeriodicTable[85] = "At";
mPeriodicTable[86] = "Rn";
mPeriodicTable[87] = "Fr";
mPeriodicTable[88] = "Ra";
mPeriodicTable[89] = "Ac";
mPeriodicTable[90] = "Th";
mPeriodicTable[91] = "Pa";
mPeriodicTable[92] = "U";
mPeriodicTable[93] = "Np";
mPeriodicTable[94] = "Pu";
mPeriodicTable[95] = "Am";
mPeriodicTable[96] = "Cm";
mPeriodicTable[97] = "Bk";
mPeriodicTable[251] = "Cf";
mPeriodicTable[252] = "Es";
mPeriodicTable[100] = "Fm";
mPeriodicTable[258] = "Md";
mPeriodicTable[102] = "No";
mPeriodicTable[103] = "Lr";
mPeriodicTable[261] = "Rf";
mPeriodicTable[105] = "Db";
mPeriodicTable[106] = "Sg";
mPeriodicTable[107] = "Bh";
mPeriodicTable[108] = "Hs";
mPeriodicTable[109] = "Mt";
//
- // Register parameters
+ // Register parameters (ptr, name, default)
//
registerParameter(&mUserInt, "userInt", 0);
registerParameter(&mUserDouble, "userDouble", 0.);
registerParameter(&mUserString, "userString", string(""));
registerParameter(&mA, "A", static_cast<unsigned int>(1));
registerParameter(&mVectorMesonId, "vectorMesonId", 443); // J/psi
registerParameter(&mDipoleModelName, "dipoleModel", string("bSat"));
registerParameter(&mDipoleModelParameterSetName, "dipoleModelParameterSet", string("KMW"));
registerParameter(&mTableSetTypeName, "tableSetType", string("total_and_coherent"));
registerParameter(&mQ2min, "Q2min", 1000.); // no limits if max <= min
registerParameter(&mQ2max, "Q2max", 0.);
registerParameter(&mWmin, "Wmin", 1000.);
registerParameter(&mWmax, "Wmax", 0.);
registerParameter(&mXpomMin, "xpomMin", 1e-8);
registerParameter(&mXpomMax, "xpomMax", 1.);
registerParameter(&mVerbose, "verbose", false);
registerParameter(&mVerboseLevel, "verboseLevel", 0);
registerParameter(&mRootfile, "rootfile", string("sartre.root"));
registerParameter(&mSeed, "seed", static_cast<unsigned int>(time(0)));
registerParameter(&mUPC, "UPC", false);
registerParameter(&mUPCA, "UPCA", static_cast<unsigned int>(1));
}
Settings::~Settings()
{
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
delete mRegisteredParameters[k];
}
mRegisteredParameters.clear();
}
int Settings::userInt() const {return mUserInt;}
double Settings::userDouble() const {return mUserDouble;}
string Settings::userString() const {return mUserString;}
void Settings::setUserInt(int val) {mUserInt = val;}
void Settings::setUserDouble(double val) {mUserDouble = val;}
void Settings::setUserString(const string& val) {mUserString = val;}
unsigned int Settings::seed() const {return mSeed;}
void Settings::setSeed(unsigned int val)
{
mSeed = val;
mRandomGenerator.SetSeed(mSeed);
gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()
}
bool Settings::readSettingsFromFile(const char *file)
{
if (!file) return false; // nothing to do
mRuncard = string(file);
mLines.clear();
//
// Open file
//
ifstream ifs(file);
if (!ifs) {
cout << "Settings::readSettingsFromFile(): error, cannot open file '"
<< file << "'." << endl;
return false;
}
//
// Read file into vector of strings, skip comments and empty lines
//
while (ifs.good() && !ifs.eof()) {
string line;
getline (ifs, line);
if (ifs.eof() && line.empty()) break;
// empty line
if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) continue;
// if first character is not a letter/digit, then taken to be a comment.
int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
if (!isalnum(line[firstChar])) continue;
mLines.push_back(line);
}
ifs.close(); // done with I/O
//
// Process vector of strings one at a time and use
// it to set registered variables.
//
for (unsigned int i=0; i<mLines.size(); i++) {
// replace '=' by blank to make parsing simpler.
while (mLines[i].find("=") != string::npos) {
int firstEqual = mLines[i].find_first_of("=");
mLines[i].replace(firstEqual, 1, " ");
}
// replace ':' by blank to make parsing simpler.
while (mLines[i].find(":") != string::npos) {
int firstEqual = mLines[i].find_first_of(":");
mLines[i].replace(firstEqual, 1, " ");
}
// get identifier string
istringstream splitLine(mLines[i]);
string name;
splitLine >> name;
// find value string
string valueString;
splitLine >> valueString;
if (!splitLine) {
cout << "Settings::readSettingsFromFile(): error, value of variable '"
<< name.c_str() << "' not recognized." << endl;
}
istringstream modeData(valueString);
//
// Loop over registered variables and see which fits.
// Not particular elegant but does the job and saves
// a lot of programming in derived classes.
//
bool isRegistered = false;
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<vector<double>>)) { // test
SettingsParameter<vector<double>> *var = dynamic_cast<SettingsParameter<vector<double>>*> (mRegisteredParameters[k]);
if (var->name == name) {
var->address->push_back(atof(valueString.c_str())); // first value
while (splitLine.good() && !splitLine.eof()) { // get remaining
string nextValue;
splitLine >> nextValue;
var->address->push_back(atof(nextValue.c_str()));
if (splitLine.eof()) break;
}
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<double>)) {
SettingsParameter<double> *var = dynamic_cast<SettingsParameter<double>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<int>)) {
SettingsParameter<int> *var = dynamic_cast<SettingsParameter<int>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned int>)) {
SettingsParameter<unsigned int> *var = dynamic_cast<SettingsParameter<unsigned int>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned long>)) {
SettingsParameter<unsigned long> *var = dynamic_cast<SettingsParameter<unsigned long>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<string>)) {
SettingsParameter<string> *var = dynamic_cast<SettingsParameter<string>*> (mRegisteredParameters[k]);
if (var->name == name) {
*(var->address) = valueString;
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<bool>)) {
SettingsParameter<bool> *var = dynamic_cast<SettingsParameter<bool>*> (mRegisteredParameters[k]);
if (var->name == name) {
isRegistered = true;
if (valueString == string("true") ||
valueString == string("True") ||
valueString == string("TRUE") ||
valueString == string("on") ||
valueString == string("On") ||
valueString == string("ON") ||
valueString == string("Yes") ||
valueString == string("yes") ||
valueString == string("YES") ||
valueString == string("T") ||
valueString == string("t") ||
valueString == string("1") ) {
*(var->address) = true;
}
else {
*(var->address) = false;
}
}
}
}
if (!isRegistered) {
cout << "Settings::readSettingsFromFile(): error, parameter identifier '" <<
name.c_str() << "' not recognized." << endl;
}
}
//
// Consolidate input (after burner)
//
consolidateCommonSettings();
consolidateSettings(); // overloaded
return true;
}
bool Settings::list(ostream& os)
{
const int fieldWidth = 28;
os << "\nRun Settings:" << endl;
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<double>)) {
SettingsParameter<double> *var = dynamic_cast<SettingsParameter<double>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<int>)) {
SettingsParameter<int> *var = dynamic_cast<SettingsParameter<int>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned int>)) {
SettingsParameter<unsigned int> *var = dynamic_cast<SettingsParameter<unsigned int>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned long>)) {
SettingsParameter<unsigned long> *var = dynamic_cast<SettingsParameter<unsigned long>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<string>)) {
SettingsParameter<string> *var = dynamic_cast<SettingsParameter<string>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << var->address->c_str() << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<bool>)) {
SettingsParameter<bool> *var = dynamic_cast<SettingsParameter<bool>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << (*(var->address) ? "true" : "false") << endl;
}
}
os << endl;
return true;
}
TParticlePDG* Settings::lookupPDG(int id) const
{
if (mPDG)
return mPDG->GetParticle(id);
else
return 0;
}
string Settings::particleName(int pdgID)
{
string name("unknown");
if (abs(pdgID) < 1000000000) { // particle
if (mPDG) {
TParticlePDG *part = lookupPDG(pdgID);
if (part) name = part->GetName();
}
if (abs(pdgID) == 990) name = "pomeron";
}
else { // nucleus in 10LZZZAAAI PDG format
int id = pdgID;
// int iso = id%10;
id /= 10;
int A = id%1000;
id /= 1000;
int Z = id%1000;
stringstream namestream;
namestream << mPeriodicTable[Z] << "(" << A << ")";
name = namestream.str();
}
return name;
}
void Settings::setVerbose(bool val) {
mVerbose = val;
if (mVerbose && mVerboseLevel == 0) mVerboseLevel = 1;
if (!mVerbose && mVerboseLevel != 0) mVerboseLevel = 0;
}
bool Settings::verbose() const {return mVerbose;}
void Settings::setVerboseLevel(int val) {
mVerboseLevel = val;
if (!mVerbose && mVerboseLevel != 0) mVerbose = true;
if (mVerbose && mVerboseLevel == 0) mVerbose = false;
}
int Settings::verboseLevel() const {return mVerboseLevel;}
void Settings::setQ2min(double val) { mQ2min = val;}
double Settings::Q2min() const {return mQ2min;}
double Settings::Qmin() const {return sqrt(mQ2min);}
void Settings::setQ2max(double val) { mQ2max = val;}
double Settings::Q2max() const {return mQ2max;}
double Settings::Qmax() const {return sqrt(mQ2max);}
void Settings::setW2min(double val) { mWmin = sqrt(val);}
void Settings::setWmin(double val) { mWmin = val;}
double Settings::Wmin() const {return mWmin;}
double Settings::W2min() const {return mWmin*mWmin;}
void Settings::setW2max(double val) { mWmax = sqrt(val);}
void Settings::setWmax(double val) { mWmax = val;}
double Settings::Wmax() const {return mWmax;}
double Settings::W2max() const {return mWmax*mWmax;}
void Settings::setXpomMin(double val) {mXpomMin = val;}
void Settings::setXpomMax(double val) {mXpomMax = val;}
double Settings::xpomMin() const {return mXpomMin;}
double Settings::xpomMax() const {return mXpomMax;}
int Settings::vectorMesonId() const {return mVectorMesonId;}
void Settings::setVectorMesonId(int val) {mVectorMesonId = val;}
string Settings::dipoleModelName() const {return mDipoleModelName;}
DipoleModelType Settings::dipoleModelType() const {return mDipoleModelType;}
void Settings::setDipoleModelType(DipoleModelType val)
{
mDipoleModelType = val;
if (mDipoleModelType == bSat)
mDipoleModelName = string("bSat");
else if (mDipoleModelType == bNonSat)
mDipoleModelName = string("bNonSat");
else if (mDipoleModelType == bCGC)
mDipoleModelName = string("bCGC");
}
unsigned int Settings::A() const {return mA;}
void Settings::setA(unsigned int val) {mA = val;}
void Settings::setRootfile(const char* val){ mRootfile = val; }
string Settings::rootfile() const { return mRootfile; }
string Settings::dipoleModelParameterSetName() const {return mDipoleModelParameterSetName;}
DipoleModelParameterSet Settings::dipoleModelParameterSet() const {return mDipoleModelParameterSet;}
void Settings::setDipoleModelParameterSet(DipoleModelParameterSet val)
{
mDipoleModelParameterSet = val;
if (mDipoleModelParameterSet == KMW)
mDipoleModelParameterSetName = string("KMW");
else if (mDipoleModelParameterSet == HMPZ)
mDipoleModelParameterSetName = string("HMPZ");
else if (mDipoleModelParameterSet == CUSTOM)
mDipoleModelParameterSetName = string("CUSTOM");
}
string Settings::tableSetTypeName() const {return mTableSetTypeName;}
TableSetType Settings::tableSetType() const {return mTableSetType;}
void Settings::setTableSetType(TableSetType val)
{
mTableSetType = val;
if (mTableSetType == total_and_coherent)
mTableSetTypeName = string("total_and_coherent");
else if (mTableSetType == coherent_and_incoherent)
mTableSetTypeName = string("coherent_and_incoherent");
}
void Settings::consolidateCommonSettings()
{
//
// Check if verbose levels and flags are consistent
// The verbose flag superseeds the verboseLevel.
//
if (mVerbose && mVerboseLevel == 0) mVerboseLevel = 1;
if (mVerboseLevel != 0 && !mVerbose) mVerboseLevel = 0;
//
// Set random generator seed
//
mRandomGenerator.SetSeed(mSeed);
gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()
//
// Dipole Model
//
if (mDipoleModelName == string("bSat"))
mDipoleModelType = bSat;
else if (mDipoleModelName == string("bNonSat"))
mDipoleModelType = bNonSat;
else if (mDipoleModelName == string("bCGC"))
mDipoleModelType = bCGC;
else {
cout << "Settings::consolidateCommonSettings(): Error, dipole model '"
<< mDipoleModelName.c_str() << "' is not defined." << endl;
exit(1);
}
//
// Dipole Model Parameter Set
//
if (mDipoleModelParameterSetName == string("KMW"))
mDipoleModelParameterSet = KMW;
else if (mDipoleModelParameterSetName == string("HMPZ"))
mDipoleModelParameterSet = HMPZ;
else if (mDipoleModelParameterSetName == string("CUSTOM"))
mDipoleModelParameterSet = CUSTOM;
else {
cout << "Settings::consolidateCommonSettings(): Error, dipole model parameter set'"
<< mDipoleModelParameterSetName.c_str() << "' is not defined." << endl;
exit(1);
}
//
// Table Set Type
//
if (mTableSetTypeName == string("total_and_coherent"))
mTableSetType = total_and_coherent;
else if (mTableSetTypeName == string("coherent_and_incoherent"))
mTableSetType = coherent_and_incoherent;
else {
cout << "Settings::consolidateCommonSettings(): Error, table set type '"
<< mTableSetTypeName.c_str() << "' is not defined." << endl;
exit(1);
}
}
void Settings::setUPC(bool val){ mUPC = val; }
bool Settings::UPC() const { return mUPC; }
void Settings::setUPCA(unsigned int val){ mUPCA = val; }
unsigned int Settings::UPCA() const { return mUPCA; }
Index: trunk/src/DipoleModelParameters.cpp
===================================================================
--- trunk/src/DipoleModelParameters.cpp (revision 342)
+++ trunk/src/DipoleModelParameters.cpp (revision 343)
@@ -1,435 +1,468 @@
//==============================================================================
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// $Date$
// $Author: ullrich $
//==============================================================================
#include "DipoleModelParameters.h"
#include "Settings.h"
#include "TableGeneratorSettings.h"
#include <iomanip>
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;
-
+ 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) {
- // Internal note by Heikki Mantysaari an Pia Zurita, arXiv pending
+ // 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;
-
+ 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() < 9) {
- cout << "DipoleModelParameters::setupParameters(): Error, require 8 custom parameters for bSAT when" << endl;
+ 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];
-
- mBG = mCustomParameters[4];
- mMu02 = mCustomParameters[5];
- mLambdaG = mCustomParameters[6];
- mAg = mCustomParameters[7];
- mC = mCustomParameters[8];
+ 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() < 9) {
- cout << "DipoleModelParameters::setupParameters(): Error, require 8 custom parameters for bNonSAT when" << endl;
+ 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[4];
- mMu02 = mCustomParameters[5];
- mLambdaG = mCustomParameters[6];
- mAg = mCustomParameters[7];
- mC = mCustomParameters[8];
+ 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[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.
+ //
+ 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;
+ }
+
if (mDipoleModelParameterSet == KMW) {
//
// KMW: bSat, bNonSat, and bCGC use the same parameters
- // and also does not distingiosh between T and L.
+ // 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;
}