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 . // // Author: Thomas Ullrich // $Date$ // $Author: ullrich $ //============================================================================== #ifndef DipoleModelParameters_h #define DipoleModelParameters_h #include "Enumerations.h" #include #include 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 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 . // // 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(); 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 . // // Author: Tobias Toll // Last update: // $Date$ // $Author$ //============================================================================== #include "WaveOverlap.h" #include "Constants.h" #include "TableGeneratorSettings.h" #include "DipoleModelParameters.h" #include "TMath.h" #include #include 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/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 . // // 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; - + 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; } 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 . // // Author: Thomas Ullrich // Last update: // $Date$ // $Author$ //============================================================================== #include "Settings.h" #include #include #include #include #include #include "TParticlePDG.h" #include #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(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(time(0))); registerParameter(&mUPC, "UPC", false); registerParameter(&mUPCA, "UPCA", static_cast(1)); } Settings::~Settings() { for (unsigned int k=0; kSetSeed(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> 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>)) { // test SettingsParameter> *var = dynamic_cast>*> (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)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); if (var->name == name) { modeData >> (*(var->address)); isRegistered = true; } } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); if (var->name == name) { modeData >> (*(var->address)); isRegistered = true; } } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); if (var->name == name) { modeData >> (*(var->address)); isRegistered = true; } } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); if (var->name == name) { modeData >> (*(var->address)); isRegistered = true; } } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); if (var->name == name) { *(var->address) = valueString; isRegistered = true; } } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (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)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl; } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl; } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl; } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl; } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]); os << setw(fieldWidth) << var->name.c_str() << "\t" << var->address->c_str() << endl; } else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) { SettingsParameter *var = dynamic_cast*> (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; }