Index: trunk/src/Amplitudes.cpp
===================================================================
--- trunk/src/Amplitudes.cpp (revision 377)
+++ trunk/src/Amplitudes.cpp (revision 378)
@@ -1,334 +1,334 @@
//==============================================================================
// Amplitudes.cpp
//
-// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
//#define SARTRE_IN_MULTITHREADED_MODE 1
#include
#include
#include "Amplitudes.h"
#include "Constants.h"
#include "TableGeneratorSettings.h"
#include "DglapEvolution.h"
#include "Enumerations.h"
#include "Kinematics.h"
#include "Integrals.h"
#include "DipoleModel.h"
#if defined(SARTRE_IN_MULTITHREADED_MODE)
#include
#endif
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
Amplitudes::Amplitudes()
{
mAmplitudeT = 0;
mAmplitudeL = 0;
mAmplitudeT2 = 0;
mAmplitudeL2 = 0;
mNumberOfConfigurations = 0;
mTheModes = 0;
mA = 0;
mErrorT = 0;
mErrorL = 0;
mErrorT2 = 0;
mErrorL2 = 0;
mAmplitudeTForSkewednessCorrection = 0;
mAmplitudeLForSkewednessCorrection = 0;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mNumberOfConfigurations = settings->numberOfConfigurations();
mVerbose = settings->verbose();
//
// Create a vector containing instances of the Integrals class
// and initialize them:
//
for (int i=0; i<=mNumberOfConfigurations; i++) {
mIntegrals.push_back(new IntegralsExclusive);
}
mA=settings->A();
mUPC=settings->UPC();
isBNonSat = false;
if(settings->dipoleModelType() == bNonSat)
isBNonSat = true;
//
// Get the modes to calculate:
// 0: analytically averaged over configurations
// 1: only analytically
// 2: Both and averaged over configurations
//
mTheModes=settings->modesToCalculate();
}
Amplitudes& Amplitudes::operator=(const Amplitudes& amp)
{
if (this != &) {
for (unsigned int i=0; idipoleModel()->createConfiguration(i);
}
//void Amplitudes::calculate(double t, double Q2, double W2)
void Amplitudes::calculate(double* kinematicPoint)
{
double t=0, Q2=0, W2=0, xpom=0;
if(!mUPC){
t =kinematicPoint[0];
Q2=kinematicPoint[1];
W2=kinematicPoint[2];
}
else{
t =kinematicPoint[0];
xpom=kinematicPoint[1];
}
#if defined(SARTRE_IN_MULTITHREADED_MODE) // multithreaded version
if (mA == 1) {
cout << "Amplitudes::calculate(): Multithreaded mode (SARTRE_IN_MULTITHREADED_MODE)" << endl;
cout << " is not supported for ep (A=1). Stopping." << endl;
exit(1);
}
//
// Create a vector containing the threads:
//
std::vector vThreads;
vThreads.clear();
//
// Create the thread group:
//
boost::thread_group gThreads;
if (mTheModes==0 || mTheModes == 2){
//Start loop over configurations, each calculated on a separate thread:
for (int i=0; i1 && (mTheModes==1 || mTheModes == 0)) {
if(!mUPC)
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2);
else
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom);
}
if (mTheModes==0 || mTheModes == 2) {
// Wait for all threads to finish before continuing main thread:
gThreads.join_all();
// Clean up threads
vThreads.clear();
}
#else // unforked version
if ((mTheModes==0 || mTheModes == 2) || mA==1) {
//Start loop over configurations:
for (int i=0; ioperator()(t, Q2, W2);
else
mIntegrals.at(i)->operator()(t, xpom);
}
}
//
// Calculate coherent cross-section according to eq.(47) in KT arXiv:hep-ph/0304189v3,
// (only in eA)
//
if (mA>1 && (mTheModes==1 || mTheModes == 0)){
if(!mUPC)
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2);
else
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom);
}
#endif
//
// Calculate the resulting :
//
double coherentT=0, coherentL=0;
double errCoherentT=0, errCoherentL=0;
if ((mTheModes==0 || mTheModes == 2) || mA==1) {
double totalT = 0;
double totalL = 0;
double err2TotalT=0, err2TotalL=0;
double probabilityCutOff = 1e-6;
for (int i=0; iintegralImT();
double valreT=mIntegrals.at(i)->integralReT();
double valimL=mIntegrals.at(i)->integralImL();
double valreL=mIntegrals.at(i)->integralReL();
double errimT=mIntegrals.at(i)->errorImT();
double errimL=mIntegrals.at(i)->errorImL();
double errreT=mIntegrals.at(i)->errorReT();
double errreL=mIntegrals.at(i)->errorReL();
double probimT=mIntegrals.at(i)->probImT();
double probimL=mIntegrals.at(i)->probImL();
double probreT=mIntegrals.at(i)->probReT();
double probreL=mIntegrals.at(i)->probReL();
if (probimT > probabilityCutOff || probimL > probabilityCutOff ||
probreT > probabilityCutOff || probreL > probabilityCutOff){
if(mVerbose) {
cout<< "Amplitudes::calculate(): Warning, Integrals may not have reached desired precision" < probreT and probimT > probimL and probimT > probreL ?
cout<< " The probability for this is "< probimT and probreT > probimL and probreT > probreL ?
cout<< " The probability for this is "< probimT and probimL > probreT and probimL > probreL ?
cout<< " The probability for this is "<1 && (mTheModes==1 || mTheModes == 0)) {
double coherentKTT=mIntegrals.at(mNumberOfConfigurations)->integralImT();
double coherentKTL=mIntegrals.at(mNumberOfConfigurations)->integralImL();
double errCoherentKTT=mIntegrals.at(mNumberOfConfigurations)->errorImT();
double errCoherentKTL=mIntegrals.at(mNumberOfConfigurations)->errorImL();
mAmplitudeT=coherentKTT;
mAmplitudeL=coherentKTL;
mErrorT=errCoherentKTT;
mErrorL=errCoherentKTL;
}
else {
mAmplitudeT=coherentT/mNumberOfConfigurations;
mAmplitudeL=coherentL/mNumberOfConfigurations;
mErrorT=errCoherentT/mNumberOfConfigurations;
mErrorL=errCoherentL/mNumberOfConfigurations;
}
if(mA==1 and mTheModes != 1 and mNumberOfConfigurations == 1){
if(isBNonSat){
mAmplitudeTForSkewednessCorrection = mAmplitudeT;
mAmplitudeLForSkewednessCorrection = mAmplitudeL;
}
else{
mAmplitudeTForSkewednessCorrection = mIntegrals.at(0)->integralTForSkewedness();
mAmplitudeLForSkewednessCorrection = mIntegrals.at(0)->integralLForSkewedness();
}
}
}
Index: trunk/src/EventGeneratorSettings.cpp
===================================================================
--- trunk/src/EventGeneratorSettings.cpp (revision 377)
+++ trunk/src/EventGeneratorSettings.cpp (revision 378)
@@ -1,116 +1,116 @@
//==============================================================================
// EventGeneratorSettings.cpp
//
-// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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 "EventGeneratorSettings.h"
#include "Constants.h"
#include
#include
EventGeneratorSettings* EventGeneratorSettings::mInstance = 0; // initialize static
EventGeneratorSettings* EventGeneratorSettings::instance()
{
if (mInstance == 0)
mInstance = new EventGeneratorSettings;
return mInstance;
}
EventGeneratorSettings::EventGeneratorSettings()
{
//
// Register all parameters.
// list() will print them in the order they are defined here.
//
registerParameter(&mElectronBeamEnergy, "eBeamEnergy", 10.);
registerParameter(&mHadronBeamEnergy, "hBeamEnergy", 100.);
registerParameter(&mNumberOfEvents, "numberOfEvents", static_cast(10000));
registerParameter(&mTimesToShow, "timesToShow", static_cast(100));
registerParameter(&mCorrectForRealAmplitude, "correctForRealAmplitude", true);
registerParameter(&mCorrectSkewedness, "correctSkewedness", true);
registerParameter(&mEnableNuclearBreakup, "enableNuclearBreakup", false);
registerParameter(&mMaxLambdaUsedInCorrections, "maxLambdaUsedInCorrections", 0.65);
registerParameter(&mMaxNuclearExcitationEnergy, "maxNuclearExcitationEnergy", 1.4);
registerParameter(&mApplyPhotonFlux, "applyPhotonFlux", true);
}
void EventGeneratorSettings::consolidateSettings() // called after runcard is read
{
//
// Disable nuclear breakup in UPC mode for now (July 2018)
//
if (mUPC && mEnableNuclearBreakup) {
cout << "EventGeneratorSettings::consolidateSettings(): Nuclear breakup is currently not possible in UPC mode.\n"
<< " Switched off now." << endl;
mEnableNuclearBreakup = false;
}
}
//
// Access functions
//
void EventGeneratorSettings::setNumberOfEvents(unsigned long val) { mNumberOfEvents = val;}
unsigned long EventGeneratorSettings::numberOfEvents() const {return mNumberOfEvents;}
void EventGeneratorSettings::setTimesToShow(unsigned int val) { mTimesToShow = val;}
unsigned int EventGeneratorSettings::timesToShow() const {return mTimesToShow;}
double EventGeneratorSettings::electronBeamEnergy() const {return mElectronBeamEnergy;}
void EventGeneratorSettings::setElectronBeamEnergy(double val) {mElectronBeamEnergy = val;}
double EventGeneratorSettings::hadronBeamEnergy() const {return mHadronBeamEnergy;}
void EventGeneratorSettings::setHadronBeamEnergy(double val) {mHadronBeamEnergy = val;}
TLorentzVector EventGeneratorSettings::eBeam() const
{
double mass2;
if(mUPC)
mass2 = protonMass2; //#TT should be mass/nucleon in UPC...?
else
mass2 = electronMass2;
return TLorentzVector(0, 0, -sqrt(mElectronBeamEnergy*mElectronBeamEnergy-mass2), mElectronBeamEnergy);
}
TLorentzVector EventGeneratorSettings::hBeam() const
{
return TLorentzVector(0, 0, sqrt(mHadronBeamEnergy*mHadronBeamEnergy-protonMass2), mHadronBeamEnergy);
}
bool EventGeneratorSettings::correctForRealAmplitude() const {return mCorrectForRealAmplitude;}
void EventGeneratorSettings::setCorrectForRealAmplitude(bool val) {mCorrectForRealAmplitude = val;}
bool EventGeneratorSettings::correctSkewedness() const {return mCorrectSkewedness;}
void EventGeneratorSettings::setCorrectSkewedness(bool val) {mCorrectSkewedness = val;}
bool EventGeneratorSettings::enableNuclearBreakup() const {return mEnableNuclearBreakup;}
void EventGeneratorSettings::setEnableNuclearBreakup(bool val) {mEnableNuclearBreakup = val;}
double EventGeneratorSettings::maxNuclearExcitationEnergy() const {return mMaxNuclearExcitationEnergy;}
void EventGeneratorSettings::setMaxNuclearExcitationEnergy(double val) {mMaxNuclearExcitationEnergy = val;}
double EventGeneratorSettings::maxLambdaUsedInCorrections() const {return mMaxLambdaUsedInCorrections;}
void EventGeneratorSettings::setMaxLambdaUsedInCorrections(double val) {mMaxLambdaUsedInCorrections = val;}
bool EventGeneratorSettings::applyPhotonFlux() const {return mApplyPhotonFlux;}
void EventGeneratorSettings::setApplyPhotonFlux(bool val) {mApplyPhotonFlux = val;}
Index: trunk/src/Event.cpp
===================================================================
--- trunk/src/Event.cpp (revision 377)
+++ trunk/src/Event.cpp (revision 378)
@@ -1,142 +1,142 @@
//==============================================================================
// Event.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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 "Event.h"
#include "EventGeneratorSettings.h"
#include
#include
void Event::list(ostream& os) const
{
ios::fmtflags fmt = os.flags(); // store io flags
//
// Event characteristics
//
os << setw(6) << "evt = " << setw(10) << left << eventNumber;
os << setw(7) << right << "Q2 = " << setw(8) << left << fixed << setprecision(3) << Q2;
os << setw(10) << right << "x = " << scientific << x << endl;
os << setw(16+7) << right << "W = " << setw(8) << left << fixed << W;
os << setw(10) << right << "y = " << fixed << y << endl;
os << setw(16+7) << right << "t = " << setw(8) << left << fixed << t;
os << setw(10) << right << "xpom = " << scientific << xpom << endl;
os << setw(16+7) << right << "pol = " << setw(8) << left << (polarization == transverse ? 'T' : 'L');
os << setw(10) << right << "diff = " << (diffractiveMode == coherent ? "coherent" : "incoherent") << endl;
os << endl;
//
// Particle record
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
os << right;
os << setw(4) << "#"
<< setw(12) << "id"
<< setw(7) << "name"
<< setw(13) << "status"
<< setw(11) << "parents"
<< setw(14) << "daughters"
<< setw(10) << "px"
<< setw(9) << "py"
<< setw(10) << "pz"
<< setw(10) << "E"
<< setw(13) << "m";
os << endl;
for (unsigned int i=0; iparticleName(particles[i].pdgId);
os << setw(12) << left << name.c_str() << ' ';
// status
os << setw(4) << right << particles[i].status << ' ';
// parents (max 2)
int nparents = particles[i].parents.size();
if (nparents == 0 ) {
os << setw(4) << right << '-' << ' ';
os << ' ';
os << setw(4) << right << '-' << ' ';
}
if (nparents == 1 ) {
os << setw(4) << right << particles[i].parents[0] << ' ';
os << ' ';
os << setw(4) << right << '-' << ' ';
}
if (nparents == 2 ) {
os << setw(4) << right << particles[i].parents[0] << ' ';
os << ' ';
os << setw(4) << right << particles[i].parents[1] << ' ';
}
if (nparents > 2 ) {
os << setw(4) << right << particles[i].parents[0] << ' ';
os << '-';
os << setw(4) << right << particles[i].parents[nparents-1] << ' ';
}
os << " ";
// daughters (max 2)
int ndaughters = particles[i].daughters.size();
if (ndaughters == 0 ) {
os << setw(4) << right << '-' << ' ';
os << ' ';
os << setw(4) << right << '-' << ' ';
}
if (ndaughters == 1 ) {
os << setw(4) << right << particles[i].daughters[0] << ' ';
os << ' ';
os << setw(4) << right << '-' << ' ';
}
if (ndaughters == 2 ) {
os << setw(4) << right << particles[i].daughters[0] << ' ';
os << ' ';
os << setw(4) << right << particles[i].daughters[1] << ' ';
}
if (ndaughters > 2 ) {
os << setw(4) << right << particles[i].daughters[0] << ' ';
os << '-';
os << setw(4) << right << particles[i].daughters[ndaughters-1] << ' ';
}
os << " ";
// 4-momentum
os << fixed;
os << setw(8) << right << setprecision(3) << particles[i].p.Px() << ' ';
os << setw(8) << right << setprecision(3) << particles[i].p.Py() << ' ';
os << setw(9) << right << setprecision(3) << particles[i].p.Pz() << ' ';
os << setw(9) << right << setprecision(3) << particles[i].p.E() << ' ';
os << " ";
// mass
if (fabs(particles[i].p.M()) < 0.01) {
os << setw(10) << right << setprecision(3) << scientific << particles[i].p.M() << ' ';
os << fixed;
}
else
os << setw(10) << right << setprecision(3) << particles[i].p.M() << ' ';
os << endl;
}
os << endl;
os.flags(fmt); // restore io flags
}
Index: trunk/src/Constants.h
===================================================================
--- trunk/src/Constants.h (revision 377)
+++ trunk/src/Constants.h (revision 378)
@@ -1,49 +1,49 @@
//==============================================================================
// Constants.h
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
#ifndef Constants_h
#define Constants_h
//
// General constants
//
const double electronMass = 0.510998902E-3; // GeV
const double electronMass2 = electronMass*electronMass; // GeV^2
const double protonMass = 0.9382700; // GeV
const double protonMass2 = protonMass*protonMass; // GeV^2
const double alpha_em = 1/137.036;
const double hbarc = 0.197327; // GeV*fm
const double hbarc2 = hbarc*hbarc; // GeV^2*fm^2
//
// Constants used in dipole model
//
const double Nc=3.;
const double quarkCharge[6] = {2./3., -1./3., -1./3., 2./3., -1./3., 2./3.}; // u, d, s, c, b, t
//
// Constants used in integartion and other numerical operations
//
const double upperIntegrationLimit=2.5; // factor: nuclear radius -> integration limits (b, r)
#endif
Index: trunk/src/AlphaStrong.cpp
===================================================================
--- trunk/src/AlphaStrong.cpp (revision 377)
+++ trunk/src/AlphaStrong.cpp (revision 378)
@@ -1,449 +1,449 @@
//==============================================================================
// AlphaStrong.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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
#include
#include "AlphaStrong.h"
#include "Math/BrentRootFinder.h"
AlphaStrong::AlphaStrong(int order, double fr2, double mur, double asmur,
double mc, double mb, double mt)
{
//
// order = 0 (LO), 1 (NLO), 2 (NNLO), 3 (NNNLO).
// fr2 = ratio of mu_f^2 to mu_r^2 (must be a fixed value).
// mur = input renormalisation scale (in GeV) for alpha_s.
// asmur = input value of alpha_s at the renormalisation scale mur.
// mc,mb,mt = heavy quark masses in GeV.
// Note: the default parameters set for Sartre in the header file.
//
const double eps = 1e-10;
const int maxf = 10000;
mWasCalled = 0;
double r0, asi, a, b;
mMassC = mc;
mMassB = mb;
mMassT = mt;
mR0 = 1./sqrt(fr2);
mOrder = order;
mFR2 = fr2;
mScale = mur;
mAsScale = asmur;
if (mur*sqrt(fr2) <= mc) {
r0 = mur;
asi = asmur;
}
else { // Solve for alpha_s at R0 = 1 GeV.
a = 0.02; // lower bound for alpha_s(R0)
b = 2.00; // upper bound for alpha_s(R0)
r0 = mR0;
//
// Now get alpha_s(R0) corresponding to alpha_s(MUR)
// using Brent root finder from ROOT library.
//
rootFinderFormula fun;
fun.mParent = this;
ROOT::Math::BrentRootFinder rootfinder;
rootfinder.SetFunction(fun, a, b);
rootfinder.Solve(maxf, 0, eps);
asi = rootfinder.Root();
}
initR0(order, fr2, r0, asi, mc, mb, mt);
}
double AlphaStrong::findR0(double asi)
{
//
// Function for dzerox
//
initR0(mOrder, mFR2, mR0, asi, mMassC, mMassB, mMassT);
return at(mScale*mScale) - mAsScale; // solve equal to zero
}
void AlphaStrong::initR0(int order, double fr2, double R0, double ASI, double MC, double MB, double MT)
{
//
// order = 0 (LO), 1 (NLO), 2 (NNLO), 3 (NNNLO).
// fr2 = ratio of mu_f^2 to mu_r^2 (must be a fixed value).
// R0 = input renormalisation scale (in GeV) for alphas_s.
// ASI = input value of alpha_s at the renormalisation scale R0.
// MC,MB,MT = heavy quark masses in GeV.
// Must have R0*sqrt(fr2) <= MC to call this function.
//
//
// QCD colour factors
//
mCA = 3.;
mCF = 4./3.;
mTR = 0.5;
//
// The lowest integer values of the Zeta function
//
mZeta[1] = 0.57721566490153;
mZeta[2] = 1.644934066848226;
mZeta[3] = 1.202056903159594;
mZeta[4] = 1.082323233711138;
mZeta[5] = 1.036927755143370;
mZeta[6] = 1.017343061984449;
mVarFlavourNumScheme = 1; // variable flavour-number scheme (VFNS)
//mVarFlavourNumScheme = 0; // fixed flavour-number scheme (FFNS)
mNumFlavorsFFNS = 4; // number of flavours for FFNS
mPertubativeOrder = order; // perturbative order of alpha_s
mIntegrationSteps = 20; // num. steps in Runge-Kutta integration
double R20 = R0*R0; // input renormalisation scale
double MC2 = MC*MC; // mu_f^2 for charm threshold
double MB2 = MB*MB; // mu_f^2 for bottom threshold
double MT2 = MT*MT; // mu_f^2 for top threshold
mLogFR = log(fr2); // log of ratio of mu_f^2 to mu_r^2
mM20 = R20 * fr2; // input factorisation scale
//
// Stop some nonsense
//
if ( (mVarFlavourNumScheme == 0) && (mNumFlavorsFFNS < 3) ) {
cout << "AlphaStrong::initR0(): Wrong flavour number for FFNS evolution. STOP." << endl;
exit(1);
}
if ( (mVarFlavourNumScheme == 0) && (mNumFlavorsFFNS > 5) ) {
cout << "AlphaStrong::initR0(): Wrong flavour number for FFNS evolution. STOP" << endl;
exit(1);
}
if ( mPertubativeOrder > 3 ) {
cout << "AlphaStrong::initR0(): Specified order in a_s too high ("
<< "mM20 = " << mM20 << ", MC2 = " << MC2 << "). STOP" << endl;
exit(1);
}
if ( (mVarFlavourNumScheme != 0) && (fr2 > 4.001) ) {
cout << "AlphaStrong::initR0(): Too low mu_r for VFNS evolution. STOP" << endl;
exit(1);
}
if ( (mVarFlavourNumScheme == 1) && (mM20 > MC2) ) {
cout << "AlphaStrong::initR0(): Too high mu_0 for VFNS evolution ("
<< "mM20 = " << mM20 << ", MC2 = " << MC2 << "). STOP" << endl;
exit(1);
}
if ( (ASI > 2.) || (ASI < 2.e-2) ) {
cout << "AlphaStrong::initR0(): alpha_s out of range. STOP" << endl;
exit(1);
}
if ( (mVarFlavourNumScheme == 1) && (MC2 > MB2) ) {
cout << "AlphaStrong::initR0(): Wrong charm-bottom mass hierarchy. STOP" << endl;
exit(1);
}
if ( (mVarFlavourNumScheme == 1) && (MB2 > MT2) ) {
cout << "AlphaStrong::initR0(): Wrong bottom-top mass hierarchy. STOP" << endl;
exit(1);
}
betaFunction();
// Keep a_s = alpha_s(mu_r^2)/(4 pi) at the input scale R0.
mAS0 = ASI / (4*M_PI);
if (mVarFlavourNumScheme != 0) {
evolution (MC2, MB2, MT2);
}
}
double AlphaStrong::at(double R2)
{
double M2 = R2 * exp(mLogFR);
int NF;
double R20, ASI, ASF, R2T, R2B, R2C;
if (mVarFlavourNumScheme == 0) {
//
// Fixed number of flavours
//
NF = mNumFlavorsFFNS;
R20 = mM20 * R2/M2;
ASI = mAS0;
ASF = as(R2, R20, mAS0, NF);
}
else {
//
// Variable number of flavours
//
if (M2 > mM2T) {
NF = 6;
R2T = mM2T * R2/M2;
ASI = mAST;
ASF = as(R2, R2T, mAST, NF);
}
else if (M2 > mM2B) {
NF = 5;
R2B = mM2B * R2/M2;
ASI = mASB;
ASF = as(R2, R2B, mASB, NF);
}
else if (M2 > mM2C) {
NF = 4;
R2C = mM2C * R2/M2;
ASI = mASC;
ASF = as(R2, R2C, mASC, NF);
}
else {
NF = 3;
R20 = mM20 * R2/M2;
ASI = mAS0;
ASF = as(R2, R20, mAS0, NF);
}
}
//
// Final value of alpha_s
//
double result = 4.*M_PI*ASF;
return result;
}
double AlphaStrong::asnf1(double ASNF, double LOGRH, int NF)
{
//
// The threshold matching of the QCD coupling in the MS(bar) scheme,
// a_s = alpha_s(mu_r^2)/(4 pi), for NF -> NF + 1 active flavours
// up to order a_s^4 (NNNLO).
//
// The value ASNF of a_s for NF flavours at the matching scale, the
// logarithm LOGRH = ln (mu_r^2/m_H^2) -- where m_H is the pole mass
// of the heavy quark -- and NF are passed as arguments to the
// function asnf1(). The order of the expansion NAORD (defined as
// the 'n' in N^nLO) is provided as data member(s).
//
// The matching coefficients are inverted from Chetyrkin, Kniehl and
// Steinhauser, Phys. Rev. Lett. 79 (1997) 2184. The QCD colour
// factors have been hard-wired in these results. The lowest integer
// values of the Zeta function are stored as data member.
//
//
// The coupling-constant matching coefficients (CMC's) up to NNNLO
// (calculated and saved in the first call of this routine)
//
if (mWasCalled != 1) {
mCCMCoefficients[1][0] = 0.;
mCCMCoefficients[1][1] = 2./3.;
mCCMCoefficients[2][0] = 14./3.;
mCCMCoefficients[2][1] = 38./3.;
mCCMCoefficients[2][2] = 4./9.;
mCCMCoefficientsI30 = + 80507./432. * mZeta[3] + 58933./1944. + 128./3. * mZeta[2] * (1.+ log(2.)/3.);
mCCMCoefficientsF30 = - 64./9. * (mZeta[2] + 2479./3456.);
mCCMCoefficientsI31 = 8941./27.;
mCCMCoefficientsF31 = - 409./27.;
mCCMCoefficients[3][2] = 511./9.;
mCCMCoefficients[3][3] = 8./27.;
mWasCalled = 1;
}
//
// The N_f dependent CMC's, and the alpha_s matching at order NAORD
//
mCCMCoefficients[3][0] = mCCMCoefficientsI30 + NF * mCCMCoefficientsF30;
mCCMCoefficients[3][1] = mCCMCoefficientsI31 + NF * mCCMCoefficientsF31;
double result = ASNF;
if (mPertubativeOrder == 0) return result;
double ASP = ASNF;
for (int K1 = 1; K1 <= mPertubativeOrder; K1++) {
ASP = ASP * ASNF;
double LRHP = 1.;
for (int K2 = 0; K2 <= K1; K2++) {
result = result + ASP * mCCMCoefficients[K1][K2] * LRHP;
LRHP = LRHP * LOGRH;
}
}
return result;
}
void AlphaStrong::evolution (double MC2, double MB2, double MT2)
{
//
// The function evolution() for the evolution of a_s = alpha_s/(4 pi)
// from a three-flavour initial scale to the four- to six-flavour
// thresholds (identified with the squares of the corresponding quark
// masses). The results are kept as data member.
//
// The input scale mM20 = mu_(f,0)^2 and the corresponding value
// mAS0 of a_s are provided as data member as is the fixed scale
// logarithm mLogFR = ln (mu_f^2/mu_r^2). The alpha_s
// matching is done by the function asnf1.
//
//
// Coupling constants at and evolution distances to/between thresholds
//
double R20 = mM20 * exp(-mLogFR);
//
// Charm
//
mM2C = MC2;
double R2C = mM2C * R20/mM20;
double ASC3 = as(R2C, R20, mAS0, 3);
// double SC = log (mAS0 / ASC3); // not used
mASC = asnf1 (ASC3, -mLogFR, 3);
//
// Bottom
//
mM2B = MB2;
double R2B = mM2B * R20/mM20;
double ASB4 = as(R2B, R2C, mASC, 4);
// double SB = log (mASC / ASB4); // not used
mASB = asnf1 (ASB4, -mLogFR, 4);
//
// Top
//
mM2T = MT2;
double R2T = mM2T * R20/mM20;
double AST5 = as(R2T, R2B, mASB, 5);
// double ST = log (mASB / AST5); // not used
mAST = asnf1 (AST5, -mLogFR, 5);
}
//
// The beta functions funBeta'n' at N^nLO for n = 1, 2, and 3
//
double AlphaStrong::funBeta1(double A, int NF) {return - A*A * (mBeta0[NF] + A * mBeta1[NF]);}
double AlphaStrong::funBeta2(double A, int NF) {return - A*A * (mBeta0[NF] + A * (mBeta1[NF] + A * mBeta2[NF]));}
double AlphaStrong::funBeta3(double A, int NF) {return - A*A * (mBeta0[NF] + A * (mBeta1[NF] + A * (mBeta2[NF] + A * mBeta3[NF])));}
double AlphaStrong::as(double R2, double R20, double AS0, int NF)
{
//
// The running coupling of QCD,
//
// AS = a_s = alpha_s(mu_r^2)/(4 pi),
//
// obtained by integrating the evolution equation for a fixed number
// of massless flavours NF. Except at leading order (LO), AS is
// obtained using a fourth-order Runge-Kutta integration.
//
// The initial and final scales R20 and R2, the value AS0 at
// R20, and NF are passed as function arguments. The coefficients
// of the beta function up to a_s^5 (N^3LO) are provided by data
// member mBetaN. The order of the expansion mPertubativeOrder
// (defined as the 'n' in N^nLO) and the number of steps
// mIntegrationSteps for the integration beyond LO are also kept
// as data member.
//
const double SXTH = 0.166666666666666;
//
// Initial value, evolution distance and step size
//
double result = AS0;
double LRRAT = log (R2/R20);
double DLR = LRRAT / mIntegrationSteps;
//
// Solution of the evolution equation depending on mPertubativeOrder
// (fourth-order Runge-Kutta beyond the leading order)
//
double XK0, XK1, XK2, XK3;
if (mPertubativeOrder == 0) {
result = AS0 / (1 + mBeta0[NF] * AS0 * LRRAT);
}
else if (mPertubativeOrder == 1) {
for (int K1 = 1; K1 <= mIntegrationSteps; K1++) {
XK0 = DLR * funBeta1 (result, NF);
XK1 = DLR * funBeta1 (result + 0.5 * XK0, NF);
XK2 = DLR * funBeta1 (result + 0.5 * XK1, NF);
XK3 = DLR * funBeta1 (result + XK2, NF);
result = result + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3);
}
}
else if (mPertubativeOrder == 2) {
for (int K1 = 1; K1 <= mIntegrationSteps; K1++) {
XK0 = DLR * funBeta2 (result, NF);
XK1 = DLR * funBeta2 (result + 0.5 * XK0, NF);
XK2 = DLR * funBeta2 (result + 0.5 * XK1, NF);
XK3 = DLR * funBeta2 (result + XK2, NF);
result = result + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3);
}
}
else if (mPertubativeOrder == 3) {
for (int K1 = 1; K1 <= mIntegrationSteps; K1++) {
XK0 = DLR * funBeta3 (result, NF);
XK1 = DLR * funBeta3 (result + 0.5 * XK0, NF);
XK2 = DLR * funBeta3 (result + 0.5 * XK1, NF);
XK3 = DLR * funBeta3 (result + XK2, NF);
result = result + SXTH * (XK0 + 2.* XK1 + 2.* XK2 + XK3);
}
}
return result;
}
void AlphaStrong::betaFunction()
{
//
// The function betaFunction() for the coefficients mBeta0...mBeta3 of the
// beta function of QCD up to order alpha_s^5 (N^3LO), normalized by
//
// d a_s / d ln mu_r^2 = - BETA0 a_s^2 - BETA1 a_s^3 - ...
//
// with a_s = alpha_s/(4*pi).
//
// The MSbar coefficients are written to the common-block BETACOM for
// NF = 3...6 (parameters NFMIN, NFMAX) quark flavours.
//
// The factors mCF, mCA and mTF are data members.
// Beyond NLO the QCD colour factors are hard-wired in this routine,
// and the numerical coefficients are truncated to six digits.
//
const int NFMIN = 3;
const int NFMAX = 6;
//
// The full LO and NLO coefficients
//
double B00 = 11./3. * mCA;
double B01 = -4./3. * mTR;
double B10 = 34./3. * mCA*mCA;
double B11 = -20./3. * mCA*mTR - 4.* mCF*mTR;
//
// Flavour-number loop and output to the array
//
for (int NF = NFMIN; NF <= NFMAX; NF++) {
mBeta0[NF] = B00 + B01 * NF;
mBeta1[NF] = B10 + B11 * NF;
mBeta2[NF] = 1428.50 - 279.611 * NF + 6.01852 * NF*NF;
mBeta3[NF] = 29243.0 - 6946.30 * NF + 405.089 * NF*NF + 1.49931 * NF*NF*NF;
}
}
Index: trunk/src/CrossSection.h
===================================================================
--- trunk/src/CrossSection.h (revision 377)
+++ trunk/src/CrossSection.h (revision 378)
@@ -1,104 +1,104 @@
//==============================================================================
// CrossSection.h
//
-// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
//
// Functor class.
//
// operator() returns d^3sig/(dt dQ2 dW2) in nb/GeV^6.
//
//===============================================================================
#ifndef CrossSection_h
#define CrossSection_h
#include "Enumerations.h"
#include "PhotonFlux.h"
class TRandom3;
class EventGeneratorSettings;
class TableCollection;
class CrossSection {
public:
CrossSection(TableCollection* = 0, TableCollection* = 0);
~CrossSection();
double operator()(double t, double Q2, double W2);
double operator()(double t, double xpom); // UPC version
double operator()(const double*); // array of t, Q2, W2 or t, xpom for UPC
double unuranPDF(const double*); // for UNU.RAN using log(Q2) (or log(xpom) for UPC)
// and returning log of cross-section
void setTableCollection(TableCollection*);
void setProtonTableCollection(TableCollection*);
GammaPolarization polarizationOfLastCall() const;
DiffractiveMode diffractiveModeOfLastCall() const;
void setCheckKinematics(bool);
protected:
friend class Sartre; // mostly for debugging and QA
double dsigdtdQ2dW2_total_checked(double t, double Q2, double W2);
double dsigdtdxp_total_checked(double t, double xpom);
double dsigdt_total(double t, double Q2, double W2, GammaPolarization) const; // modified
double dsigdt_coherent(double t, double Q2, double W2, GammaPolarization) const;
double dsigdt_incoherent(double t, double Q2, double W2, GammaPolarization) const; // new
double dsigdtdQ2dW2_total(double t, double Q2, double W2, GammaPolarization) const;
double dsigdtdQ2dW2_coherent(double t, double Q2, double W2, GammaPolarization) const;
double dsigdt_total(double t, double xpom) const; // UPC version
double dsigdt_coherent(double t, double xpom) const; // UPC version
double dsigdt_incoherent(double t, double xpom) const; // UPC version
double dsigdtdxp_total(double t, double xpom) const; // UPC
double dsigdtdxp_coherent(double t, double xpom) const; // UPC
double logDerivateOfAmplitude(double t, double Q2, double W2, GammaPolarization) const;
double logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization) const;
double realAmplitudeCorrection(double t, double Q2, double W2, GammaPolarization pol) const;
double skewednessCorrection(double t, double Q2, double W2, GammaPolarization pol) const;
double logDerivateOfAmplitude(double t, double xpom) const; // UPC version
double logDerivateOfGluonDensity(double t, double xpom) const; // UPC version
double realAmplitudeCorrection(double t, double xpom) const; // UPC version
double skewednessCorrection(double t, double xpom) const; // UPC version
double skewednessCorrection(double lambda) const;
double realAmplitudeCorrection(double lambda) const;
double UPCPhotonFlux(double t, double xpom) const;
private:
TRandom3* mRandom;
GammaPolarization mPolarization;
DiffractiveMode mDiffractiveMode;
PhotonFlux mPhotonFlux;
TableCollection* mTableCollection;
TableCollection* mProtonTableCollection;
EventGeneratorSettings* mSettings;
double mS;
double mVmMass;
bool mCheckKinematics;
};
#endif
Index: trunk/src/DipoleModelParameters.h
===================================================================
--- trunk/src/DipoleModelParameters.h (revision 377)
+++ trunk/src/DipoleModelParameters.h (revision 378)
@@ -1,137 +1,137 @@
//==============================================================================
// DipoleModelParameters.h
//
-// Copyright (C) 2016-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2016-2019 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/CMakeLists.txt
===================================================================
--- trunk/src/CMakeLists.txt (revision 377)
+++ trunk/src/CMakeLists.txt (revision 378)
@@ -1,156 +1,156 @@
#===============================================================================
# CMakeLists.txt (src)
#
-# Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+# Copyright (C) 2010-2019 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(ExternalProject)
cmake_minimum_required (VERSION 3.1)
#
# Compiler flags for release and debug version
#
set(CMAKE_C_FLAGS "-W")
set(CMAKE_C_FLAGS_DEBUG "-g")
set(CMAKE_C_FLAGS_RELEASE "-O")
message("COMPILER = ${CMAKE_CXX_COMPILER_ID}")
set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -pedantic -Wno-long-long")
if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-potentially-evaluated-expression")
endif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS_DEBUG "-g")
set(CMAKE_CXX_FLAGS_RELEASE "-O")
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
#
# Find external required packages
# (see also FindGSL.cmake and FindROOT.cmke in cmake/modules)
# Herer we need only the root library to create the table
# tools.
#
# GSL
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIR})
# ROOT
find_package(ROOT REQUIRED)
include_directories(${ROOT_INCLUDE_DIR})
set(LIBS ${LIBS} ${ROOT_LIBRARIES})
#BOOST
if (MULTITHREADED)
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED ON)
find_package(Boost 1.39 COMPONENTS thread REQUIRED)
if(Boost_FOUND)
include_directories(${Boost_INCLUDE_DIR})
endif(Boost_FOUND)
endif (MULTITHREADED)
#
# Include files from sartre package
#
include_directories(${PROJECT_SOURCE_DIR}/src)
include_directories(${PROJECT_SOURCE_DIR}/gemini)
include_directories(${PROJECT_SOURCE_DIR}/cuba)
#
# Cuba is built using the autoconf shipped with it
#
ExternalProject_Add(
cuba
DOWNLOAD_COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/cuba ${PROJECT_BINARY_DIR}/cuba
SOURCE_DIR ${PROJECT_BINARY_DIR}/cuba
PREFIX ${PROJECT_BINARY_DIR}/cuba
CONFIGURE_COMMAND ${PROJECT_BINARY_DIR}/cuba/configure --prefix=${PROJECT_BINARY_DIR}/cuba/build
BUILD_COMMAND make lib
BUILD_IN_SOURCE 1
)
#
# Defines source files for sartre library
#
set(SARTRE_SRC "AlphaStrong.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Amplitudes.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "BreakupProduct.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "CrossSection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DglapEvolution.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModel.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModelParameters.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Event.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ExclusiveFinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FrangibleNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "EventGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Integrals.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Kinematics.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ModeFinderFunctor.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleon.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "PhotonFlux.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Sartre.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Settings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Table.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableCollection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "WaveOverlap.cpp")
add_library(sartre ${SARTRE_SRC})
#
# Table tools (all stand alone programs)
#
add_executable(tableInspector tableInspector.cpp)
add_executable(tableMerger tableMerger.cpp)
add_executable(tableQuery tableQuery.cpp)
add_executable(tableDumper tableDumper.cpp)
add_executable(tableVarianceMaker tableVarianceMaker.cpp)
target_link_libraries(tableInspector sartre ${LIBS})
target_link_libraries(tableMerger sartre ${LIBS})
target_link_libraries(tableQuery sartre ${LIBS})
target_link_libraries(tableDumper sartre ${LIBS})
target_link_libraries(tableVarianceMaker sartre ${LIBS})
add_dependencies(tableInspector sartre)
add_dependencies(tableMerger sartre)
add_dependencies(tableQuery sartre)
add_dependencies(tableDumper sartre)
add_dependencies(tableVarianceMaker sartre)
#
# Install library and include files (make install) within
# the distribution tree. Top level CMakeLists.txt will install
# Sartre in final destination.
#
install(TARGETS sartre DESTINATION sartre/lib)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/lib/libcuba.a" DESTINATION sartre/lib)
FILE(GLOB AllIncludeFiles *.h)
install(FILES ${AllIncludeFiles} DESTINATION sartre/include)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/include/cuba.h" DESTINATION sartre/include)
install(TARGETS tableInspector DESTINATION sartre/bin)
install(TARGETS tableMerger DESTINATION sartre/bin)
install(TARGETS tableQuery DESTINATION sartre/bin)
install(TARGETS tableDumper DESTINATION sartre/bin)
install(TARGETS tableVarianceMaker DESTINATION sartre/bin)
Index: trunk/src/ExclusiveFinalStateGenerator.h
===================================================================
--- trunk/src/ExclusiveFinalStateGenerator.h (revision 377)
+++ trunk/src/ExclusiveFinalStateGenerator.h (revision 378)
@@ -1,41 +1,41 @@
//==============================================================================
// ExclusiveFinalStateGenerator.h
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
#ifndef ExclusiveFinalStateGenerator_h
#define ExclusiveFinalStateGenerator_h
#include "FinalStateGenerator.h"
class ExclusiveFinalStateGenerator : public FinalStateGenerator {
public:
ExclusiveFinalStateGenerator();
~ExclusiveFinalStateGenerator();
bool generate(int id, double t, double y, double Q2,
bool isIncoherent, int A, Event *event);
//UPC version:
bool generate(int id, double t, double xpom,
bool isIncoherent, int A, Event *event);
double xpomMin(double massVM, double t, TLorentzVector hBeam, TLorentzVector eBeam);
};
#endif
Index: trunk/src/BreakupProduct.h
===================================================================
--- trunk/src/BreakupProduct.h (revision 377)
+++ trunk/src/BreakupProduct.h (revision 378)
@@ -1,42 +1,42 @@
//==============================================================================
// BreakupProduct.h
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
#ifndef BreakupProduct_h
#define BreakupProduct_h
#include "TLorentzVector.h"
#include
#include
using namespace std;
struct BreakupProduct {
double Z;
double A;
double emissionTime; // in units of 1E-21 seconds since the creation of the compound nucleus
long pdgId; // PDG particle ID (for nuclei 10LZZZAAAI)
TLorentzVector p; // GeV units
string name;
};
ostream & operator<<(ostream&, const BreakupProduct&);
#endif
Index: trunk/src/DglapEvolution.h
===================================================================
--- trunk/src/DglapEvolution.h (revision 377)
+++ trunk/src/DglapEvolution.h (revision 378)
@@ -1,145 +1,145 @@
//==============================================================================
// DglapEvolution.h
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
// LO DGLAP evolution
//
// This class is a singleton.
//
// The evolution functions xG() and alphasxG()
// take as arguments:
//
// x : Bjorken momentum fraction
// Q2: Virtuality in GeV^2
//
// DglapEvolution::xG() returns x*G(x, Q^2)
// DglapEvolution::alphaSxG() returns alpha_s(Q^2)*x*G(x,Q^2)
//
// This is a rewrite of code written in F77. The documentation is
// largely left as provided in the Fortran version but includes
// the updated variable and function names. Data originally stored
// in common blocks and several static variables turned into data members.
// Some relics of the original F77 code remained such as array
// indexing starting at 1 (arrays are size +1 here).
// The code was substantially simplified and is limited to work
-// only in LO.
-//
-// Lookup Table:
-// To speed up things the class features a lookup table mechanism.
-// Lookup tables are generated by calling:
-//
-// DglapEvolution::generateLookupTable(int nx = 200, int nq2 = 200);
-//
-// where nx and nq2 is the granularity in x and Q2. Loopkup tables
-// are kept in memory. The method generateLookupTable() can be called
-// multiple times with different sized if needed. The use of the lookup
-// table is switched on or off at any time via:
-//
-// DglapEvolution::useLookupTable(bool val);
-//
-// If switched off the lookup table is *not* deleted. It just instructs
-// xG() to not use the lookup table.
+// only in LO.
+//
+// Lookup Table:
+// To speed up things the class features a lookup table mechanism.
+// Lookup tables are generated by calling:
+//
+// DglapEvolution::generateLookupTable(int nx = 200, int nq2 = 200);
+//
+// where nx and nq2 is the granularity in x and Q2. Loopkup tables
+// are kept in memory. The method generateLookupTable() can be called
+// multiple times with different sized if needed. The use of the lookup
+// table is switched on or off at any time via:
+//
+// DglapEvolution::useLookupTable(bool val);
+//
+// If switched off the lookup table is *not* deleted. It just instructs
+// xG() to not use the lookup table.
// The current state can be checked via lookupTableIsUsed().
//==============================================================================
#ifndef DglapEvolution_h
#define DglapEvolution_h
#include
#include
#include "AlphaStrong.h"
using namespace std;
class DglapEvolution {
public:
static DglapEvolution& instance();
~DglapEvolution();
double xG(double x, double Q2);
- double alphaSxG(double x, double Q2);
+ double alphaSxG(double x, double Q2);
double logDerivative(double x, double Q2); // dlog(xG)/dlog(1/x)
-
- void generateLookupTable(int nx = 200, int nq2 = 200);
- void useLookupTable(bool);
- bool lookupTableIsUsed() const;
+
+ void generateLookupTable(int nx = 200, int nq2 = 200);
+ void useLookupTable(bool);
+ bool lookupTableIsUsed() const;
private:
DglapEvolution();
private:
void anom();
void anCalc(complex& ggi, complex& ggf, complex& xn);
complex psiFunction(complex);
complex lngam(complex X);
complex beta(complex, complex);
- void reno(complex* fn, double alpq, int nmax, double ag, double lambdag);
+ void reno(complex* fn, double alpq, int nmax, double ag, double lambdag);
double xG_Interpolator(double x, double Q2);
double xG_Engine(double x, double Q2);
- double luovi( double f[4], double arg[4], double z);
+ double luovi( double f[4], double arg[4], double z);
private:
static DglapEvolution* mInstance;
double mMu02;
double mAg;
double mLambdaG;
complex mCC;
complex mAP[137][6];
complex mN[137];
double mWN[137];
double mC;
double mALPS;
double mALPC;
double mALPB;
double mALPT;
double mMUR;
double mMC;
double mMB;
- double mMT;
-
- bool mLookupTableIsFilled;
- bool mUseLookupTable;
- unsigned int mNumberOfNodesInX;
- unsigned int mNumberOfNodesInQ2;
- double mTableMinX;
- double mTableMaxX;
- double mTableMinQ2;
- double mTableMaxQ2;
- double **mLookupTable;
+ double mMT;
+
+ bool mLookupTableIsFilled;
+ bool mUseLookupTable;
+ unsigned int mNumberOfNodesInX;
+ unsigned int mNumberOfNodesInQ2;
+ double mTableMinX;
+ double mTableMaxX;
+ double mTableMinQ2;
+ double mTableMaxQ2;
+ double **mLookupTable;
AlphaStrong *mAlphaStrong;
};
-
-inline void DglapEvolution::useLookupTable(bool val) {mUseLookupTable = val;}
-
-inline bool DglapEvolution::lookupTableIsUsed() const
-{
- return mUseLookupTable && mLookupTableIsFilled;
-}
+
+inline void DglapEvolution::useLookupTable(bool val) {mUseLookupTable = val;}
+
+inline bool DglapEvolution::lookupTableIsUsed() const
+{
+ return mUseLookupTable && mLookupTableIsFilled;
+}
inline DglapEvolution& DglapEvolution::instance()
{
if (!mInstance) mInstance = new DglapEvolution;
return *mInstance;
}
#endif
Index: trunk/src/DipoleModel.h
===================================================================
--- trunk/src/DipoleModel.h (revision 377)
+++ trunk/src/DipoleModel.h (revision 378)
@@ -1,104 +1,104 @@
//==============================================================================
// DipoleModel.h
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#ifndef DipoleModel_h
#define DipoleModel_h
#include "AlphaStrong.h"
#include "TableGeneratorNucleus.h"
#include "DipoleModelParameters.h"
class TH2F;
class TH1F;
class DipoleModel {
public:
DipoleModel();
virtual ~DipoleModel();
const TableGeneratorNucleus* nucleus() const;
bool configurationExists() const;
virtual void createConfiguration(int)=0;
virtual double dsigmadb2(double, double, double, double)=0;
virtual double bDependence(double);
virtual double bDependence(double, double);
virtual double dsigmadb2ep(double, double, double);
virtual double coherentDsigmadb2(double, double, double) {return 0;};
virtual void createSigma_ep_LookupTable(double) {/* nothing */};
protected:
TableGeneratorNucleus mNucleus;
DipoleModelParameters *mParameters;
AlphaStrong mAs;
bool mConfigurationExists;
bool mIsInitialized;
};
class DipoleModel_bSat : public DipoleModel {
public:
DipoleModel_bSat();
DipoleModel_bSat(const DipoleModel_bSat&);
DipoleModel_bSat& operator=(const DipoleModel_bSat&);
~DipoleModel_bSat();
void createSigma_ep_LookupTable(double);
protected:
void createConfiguration(int);
double dsigmadb2(double, double, double, double);
double bDependence(double, double);
double dsigmadb2ep(double, double, double);
protected:
TH2F* mBDependence;
private:
double dsigmadb2epForIntegration(double*, double*);
TH1F* mSigma_ep_LookupTable;
double coherentDsigmadb2(double, double, double);
};
class DipoleModel_bNonSat : public DipoleModel_bSat{
public:
DipoleModel_bNonSat();
~DipoleModel_bNonSat();
private:
double dsigmadb2(double, double, double, double);
double dsigmadb2ep(double, double, double);
double coherentDsigmadb2(double, double, double);
};
class DipoleModel_bCGC : public DipoleModel {
public:
DipoleModel_bCGC();
private:
void createConfiguration(int);
double dsigmadb2(double, double, double, double);
double dsigmadb2ep(double, double, double);
double bDependence(double);
};
#endif
Index: trunk/src/Amplitudes.h
===================================================================
--- trunk/src/Amplitudes.h (revision 377)
+++ trunk/src/Amplitudes.h (revision 378)
@@ -1,110 +1,110 @@
//==============================================================================
// Amplitudes.h
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
//
// The Amplitudes class calculate \gamma-A dsigma/dt cross-sections
// for a given t, Q2, W2 for T and L photons.
// It also calculates the respective coherent parts of the cross-sections
//
// The results are accessed by the public functions:
// double amplitudeT() for
// double amplitudeL() for
// double amplitudeT2() for <|A_T|^2>
// double amplitudeL2() for <|A_L|^2>
//
// Units are <|A|^2> in nb/GeV^2 and in sqrt(nb/GeV^2)
//
//===============================================================================
#ifndef Amplitudes_h
#define Amplitudes_h
#include
using namespace std;
class IntegralsExclusive;
class Amplitudes {
public:
Amplitudes();
Amplitudes(const Amplitudes&);
~Amplitudes();
Amplitudes& operator=(const Amplitudes&);
// void calculate(double, double, double);
void calculate(double*);
void generateConfigurations();
double amplitudeT() const;
double amplitudeL() const;
double amplitudeT2() const;
double amplitudeL2() const;
double errorT() const;
double errorL() const;
double errorT2() const;
double errorL2() const;
double amplitudeTForSkewednessCorrection() const;
double amplitudeLForSkewednessCorrection() const;
private:
// vector of instances of the integrals class:
vector mIntegrals;
// these are the results:
double mAmplitudeT;
double mAmplitudeL;
double mAmplitudeT2;
double mAmplitudeL2;
double mAmplitudeTForSkewednessCorrection;
double mAmplitudeLForSkewednessCorrection;
//...and the (absolute) errors:
double mErrorT;
double mErrorL;
double mErrorT2;
double mErrorL2;
int mNumberOfConfigurations;
int mTheModes;
unsigned int mA;
bool mUPC;
bool mVerbose;
bool isBNonSat;
};
inline double Amplitudes::amplitudeT() const {return mAmplitudeT;}
inline double Amplitudes::amplitudeL() const {return mAmplitudeL;}
inline double Amplitudes::amplitudeT2() const {return mAmplitudeT2;}
inline double Amplitudes::amplitudeL2() const {return mAmplitudeL2;}
inline double Amplitudes::amplitudeTForSkewednessCorrection() const {return mAmplitudeTForSkewednessCorrection;}
inline double Amplitudes::amplitudeLForSkewednessCorrection() const {return mAmplitudeLForSkewednessCorrection;}
inline double Amplitudes::errorT() const {return mErrorT;}
inline double Amplitudes::errorL() const {return mErrorL;}
inline double Amplitudes::errorT2() const {return mErrorT2;}
inline double Amplitudes::errorL2() const {return mErrorL2;}
#endif
Index: trunk/src/EventGeneratorSettings.h
===================================================================
--- trunk/src/EventGeneratorSettings.h (revision 377)
+++ trunk/src/EventGeneratorSettings.h (revision 378)
@@ -1,103 +1,103 @@
//==============================================================================
// EventGeneratorSettings.h
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
//
// Singleton class.
//
// Notes on verbose settings:
// verbose = false: most messages are suppressed, only severe errors are
// printed. That is the default mode.
// verbose = true: amount of print out depends on the verbose level.
// verboseLevel: 1 default - useful for production mode, additional info
// 2 more print-out, information of internal processes,
// still OK for production runs.
// 3 even more print-out, OK for test runs but not production
// 10 only for debugging - massive output for every event
//===============================================================================
#ifndef EventGeneratorSettings_h
#define EventGeneratorSettings_h
#include "Settings.h"
#include "TLorentzVector.h"
using namespace std;
class EventGeneratorSettings : public Settings {
public:
static EventGeneratorSettings* instance();
void setNumberOfEvents(unsigned long);
unsigned long numberOfEvents() const;
void setTimesToShow(unsigned int);
unsigned int timesToShow() const;
double electronBeamEnergy() const;
void setElectronBeamEnergy(double);
double hadronBeamEnergy() const;
void setHadronBeamEnergy(double);
TLorentzVector eBeam() const;
TLorentzVector hBeam() const;
bool correctForRealAmplitude() const;
void setCorrectForRealAmplitude(bool);
bool correctSkewedness() const;
void setCorrectSkewedness(bool);
bool enableNuclearBreakup() const;
void setEnableNuclearBreakup(bool);
double maxNuclearExcitationEnergy() const;
void setMaxNuclearExcitationEnergy(double);
double maxLambdaUsedInCorrections() const;
void setMaxLambdaUsedInCorrections(double);
bool applyPhotonFlux() const;
void setApplyPhotonFlux(bool);
private:
EventGeneratorSettings();
void consolidateSettings();
private:
static EventGeneratorSettings* mInstance;
private:
unsigned long mNumberOfEvents;
unsigned int mTimesToShow;
double mElectronBeamEnergy;
double mHadronBeamEnergy;
bool mCorrectForRealAmplitude;
bool mCorrectSkewedness;
bool mEnableNuclearBreakup;
bool mApplyPhotonFlux;
double mMaxLambdaUsedInCorrections;
double mMaxNuclearExcitationEnergy;
};
#endif
Index: trunk/src/Enumerations.h
===================================================================
--- trunk/src/Enumerations.h (revision 377)
+++ trunk/src/Enumerations.h (revision 378)
@@ -1,33 +1,33 @@
//==============================================================================
// Enumerations.h
//
-// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
#ifndef Enumerations_h
#define Enumerations_h
enum DipoleModelType {bSat, bNonSat, bCGC};
enum GammaPolarization {transverse, longitudinal};
enum AmplitudeMoment {mean_A, mean_A2, variance_A, lambda_real, lambda_skew};
enum DiffractiveMode {coherent, incoherent};
enum DipoleModelParameterSet {KMW, HMPZ, CUSTOM};
enum TableSetType {total_and_coherent, coherent_and_incoherent};
#endif
Index: trunk/src/CrossSection.cpp
===================================================================
--- trunk/src/CrossSection.cpp (revision 377)
+++ trunk/src/CrossSection.cpp (revision 378)
@@ -1,921 +1,921 @@
//==============================================================================
// CrossSection.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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 "CrossSection.h"
#include "EventGeneratorSettings.h"
#include "Kinematics.h"
#include "TableCollection.h"
#include "Table.h"
#include "Constants.h"
#include "TH2F.h"
#include "TFile.h"
#include "DglapEvolution.h"
#include
#include
#include
#define PR(x) cout << #x << " = " << (x) << endl;
CrossSection::CrossSection(TableCollection* tc, TableCollection* ptc)
{
mSettings = EventGeneratorSettings::instance();
mRandom = mSettings->randomGenerator();
mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
mPhotonFlux.setS(mS);
mTableCollection = tc;
mProtonTableCollection = ptc;
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mSettings->vectorMesonId());
mVmMass = vectorMesonPDG->Mass();
mCheckKinematics = true;
}
CrossSection::~CrossSection() { /* no op */ }
void CrossSection::setTableCollection(TableCollection* tc) {mTableCollection = tc;}
void CrossSection::setProtonTableCollection(TableCollection* ptc) {mProtonTableCollection = ptc;}
void CrossSection::setCheckKinematics(bool val) {mCheckKinematics = val;}
double CrossSection::operator()(double t, double Q2, double W2)
{
return dsigdtdQ2dW2_total_checked(t, Q2, W2);
}
// UPC Version
double CrossSection::operator()(double t, double xpom)
{
return dsigdtdxp_total_checked(t, xpom);
}
double CrossSection::operator()(const double* array)
{
if (mSettings->UPC())
return dsigdtdxp_total_checked(array[0], array[1]);
else
return dsigdtdQ2dW2_total_checked(array[0], array[1], array[2]);
}
//
// PDF passed to UNURAN
// Array holds: t, log(Q2), W2 for e+p/A running
// t, log(xpom) for UPC
//
double CrossSection::unuranPDF(const double* array) // array is t, log(Q2), W2
{ // or t and log(xpom)
double result = 0;
if (mSettings->UPC()) {
double xpom = exp(array[1]);
result = dsigdtdxp_total_checked(array[0], xpom); // t, xpom
result *= xpom; // Jacobian
}
else {
double Q2 = exp(array[1]);
result = dsigdtdQ2dW2_total_checked(array[0], Q2, array[2]); // t, Q2, W2
result *= Q2; // Jacobian
}
return log(result);
}
double CrossSection::dsigdtdQ2dW2_total_checked(double t, double Q2, double W2)
{
double result = 0;
//
// Check if in kinematically allowed region
//
if (mCheckKinematics && !Kinematics::valid(mS, t, Q2, W2, mVmMass, false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): warning, t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2)
<< " is outside of kinematically allowed region. Return 0." << endl;
return result;
}
//
// Total cross-section dsig2/(dQ2 dW2 dt)
// This is the probability density function needed for UNU.RAN
//
double csT = dsigdtdQ2dW2_total(t, Q2, W2, transverse);
double csL = dsigdtdQ2dW2_total(t, Q2, W2, longitudinal);
result = csT + csL;
//
// Polarization
//
if (mRandom->Uniform(result) <= csT)
mPolarization = transverse;
else
mPolarization = longitudinal;
//
// Diffractive Mode
//
double sampleRange = (mPolarization == transverse ? csT : csL);
double sampleDivider = dsigdtdQ2dW2_coherent(t, Q2, W2, mPolarization);
if (mRandom->Uniform(sampleRange) <= sampleDivider)
mDiffractiveMode = coherent;
else
mDiffractiveMode = incoherent;
//
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): " << result;
cout << " at t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2);
cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal");
cout << " ," << (mDiffractiveMode == coherent ? "coherent" : "incoherent");
cout << ')' << endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, return value = NaN at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, return value = inf at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
if (result < 0) {
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, negative cross-section at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
return result;
}
double CrossSection::dsigdtdxp_total_checked(double t, double xpom)
{
double result = 0;
//
// Check if in kinematically allowed region
//
if (mCheckKinematics && !Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
t, xpom, mVmMass,
(mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "CrossSection::dsigdtdxp_total_checked(): warning, t=" << t << ", xpom=" << xpom
<< " is outside of kinematically allowed region. Return 0." << endl;
return result;
}
//
// Total cross-section dsig2/(dt dxp)
// This is the probability density function needed for UNU.RAN
//
result = dsigdtdxp_total(t, xpom);
//
// Polarization
//
mPolarization = transverse; // always
//
// Diffractive Mode
//
double sampleRange = result;
double sampleDivider = dsigdtdxp_coherent(t, xpom);
if (mRandom->Uniform(sampleRange) <= sampleDivider)
mDiffractiveMode = coherent;
else
mDiffractiveMode = incoherent;
//
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) {
cout << "CrossSection::dsigdtdxp_total_checked(): " << result;
cout << " at t=" << t << ", xp=" << xpom;
cout << " (" << (mDiffractiveMode == coherent ? "coherent" : "incoherent");
cout << ')' << endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "CrossSection::dsigdtdxp_total_checked(): Error, return value = NaN at" << endl;
cout << " t=" << t << ", xp=" << xpom << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "CrossSection::dsigdtdxp_total_checked(): Error, return value = inf at" << endl;
cout << " t=" << t << ", xp=" << xpom << endl;
result = 0;
}
if (result < 0) {
cout << "CrossSection::dsigdtdxp_total_checked(): Error, negative cross-section at" << endl;
cout << " t=" << t << ", xp=" << xpom << endl;
result = 0;
}
return result;
}
double CrossSection::dsigdt_total(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = 0;
if (mSettings->tableSetType() == coherent_and_incoherent) {
result = dsigdt_coherent(t, Q2, W2, pol) + dsigdt_incoherent(t, Q2, W2, pol);
}
else if (mSettings->tableSetType() == total_and_coherent) {
result = mTableCollection->get(Q2, W2, t, pol, mean_A2); // units now are fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()){
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()){
lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
result *= skewednessCorrection(lambda);
}
}
return result;
}
//
// UPC Version
//
double CrossSection::dsigdt_total(double t, double xpom) const
{
double result = 0;
if (mSettings->tableSetType() == coherent_and_incoherent) {
result = dsigdt_coherent(t, xpom) + dsigdt_incoherent(t, xpom);
}
else if (mSettings->tableSetType() == total_and_coherent) {
result = mTableCollection->get(xpom, t, mean_A2); // units now are fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()) {
lambda = logDerivateOfAmplitude(t, xpom);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()) {
lambda = logDerivateOfGluonDensity(t, xpom);
result *= skewednessCorrection(lambda);
}
}
return result;
}
double CrossSection::dsigdt_coherent(double t, double Q2, double W2, GammaPolarization pol) const
{
double val = mTableCollection->get(Q2, W2, t, pol, mean_A); // fm^2
double result = val*val/(16*M_PI); // units now are fm^4
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()) {
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()) {
lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
result *= skewednessCorrection(lambda);
}
return result;
}
//
// UPC Version
//
double CrossSection::dsigdt_coherent(double t, double xpom) const
{
double val = mTableCollection->get(xpom, t, mean_A); // fm^2
double result = val*val/(16*M_PI); // units now are fm^4
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()) {
lambda = logDerivateOfAmplitude(t, xpom);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()) {
lambda = logDerivateOfGluonDensity(t, xpom);
result *= skewednessCorrection(lambda);
}
return result;
}
double CrossSection::dsigdtdQ2dW2_total(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = dsigdt_total(t, Q2, W2, pol);
if (mSettings->applyPhotonFlux()) result *= mPhotonFlux(Q2,W2,pol);
return result;
}
//
// UPC version
//
double CrossSection::dsigdtdxp_total(double t, double xpom) const
{
double result = dsigdt_total(t, xpom);
if (mSettings->applyPhotonFlux()) {
result *= UPCPhotonFlux(t, xpom);
}
return result;
}
double CrossSection::dsigdt_incoherent(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = mTableCollection->get(Q2, W2, t, pol, variance_A); // fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()){
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()){
lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
result *= skewednessCorrection(lambda);
}
return result;
}
//
// UPC Version
//
double CrossSection::dsigdt_incoherent(double t, double xpom) const
{
double result = mTableCollection->get(xpom, t, variance_A); // fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()) {
lambda = logDerivateOfAmplitude(t, xpom);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()) {
lambda = logDerivateOfGluonDensity(t, xpom);
result *= skewednessCorrection(lambda);
}
return result;
}
double CrossSection::dsigdtdQ2dW2_coherent(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = dsigdt_coherent(t, Q2, W2, pol);
if (mSettings->applyPhotonFlux()) result *= mPhotonFlux(Q2,W2,pol);
return result;
}
//
// UPC Version
//
double CrossSection::dsigdtdxp_coherent(double t, double xpom) const
{
double result = dsigdt_coherent(t, xpom);
if (mSettings->applyPhotonFlux()) {
result *= UPCPhotonFlux(t, xpom);
}
return result;
}
double CrossSection::UPCPhotonFlux(double t, double xpom) const{
double Egamma = Kinematics::Egamma(xpom, t, mVmMass,
mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy());
double result = mPhotonFlux(Egamma);
//
// Jacobian = dEgamma/dN, such that dsig/dtdxp = dsigd/tdEgam*dEgam/dxp
//
double h=xpom*1e-3;
double Egamma_p = Kinematics::Egamma(xpom+h, t, mVmMass,
mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy());
double Egamma_m = Kinematics::Egamma(xpom-h, t, mVmMass,
mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy());
double jacobian = (Egamma_p-Egamma_m)/(2*h);
result *= fabs(jacobian);
return result;
}
GammaPolarization CrossSection::polarizationOfLastCall() const {return mPolarization;}
DiffractiveMode CrossSection::diffractiveModeOfLastCall() const {return mDiffractiveMode;}
double CrossSection::logDerivateOfAmplitude(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfAmplitude(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
// Note (TU): if the lambda value from a table is not valid and not > 0,
// get() returns lambda=0 and table=0. This enforces a renewed
// calculation.
//
lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_real, table);
if (!table) { // no lambda value from correct table, use fallback solution
lambdaFromTable = false;
double value = mProtonTableCollection->get(Q2, W2, t, pol, mean_A, table); // use obtained table from here on
if (value <= 0) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << value << '.' << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", pol="
<< (pol == transverse ? 'T' : 'L') << endl;
return 0;
}
if (!table) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid pointer to lookup table." << endl;
return 0;
}
//
// Note: the derivate taken from values in the table is a delicate issue.
// The standard interpolation method(s) used in Table::get() are
// at times not accurate and causes ripples on lambda(W2).
// However, interpolation methods or robust fitting is by far too
// expensive in terms of CPU time per point.
//
//
// Calculate the derivative using simple method.
//
double derivative;
double hplus, hminus;
double dW2 = table->binWidthW2();
hplus = hminus = 0.5*dW2; // half bin width is found to be the best choice after some testing
hplus = min(hplus, table->maxW2()-W2);
hminus = min(hminus, W2-table->minW2());
hminus -= numeric_limits::epsilon();
hplus -= numeric_limits::epsilon();
if (hminus < 0) hminus = 0;
if (hplus < 0) hplus = 0;
if (hplus == 0 && hminus == 0) {
cout << "CrossSection::logDerivateOfAmplitude(): Warning, cannot find derivative." << endl;
return 0;
}
double a = table->get(Q2, W2+hplus, t);
if (a <= 0) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << a << '.' << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2+hplus) << ", pol="
<< (pol == transverse ? 'T' : 'L') << endl;
return 0;
}
double b = table->get(Q2, W2-hminus, t);
if (b <= 0) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << b << '.' << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2-hminus) << ", pol="
<< (pol == transverse ? 'T' : 'L') << endl;
return 0;
}
derivative = (a-b)/(hplus+hminus);
//
// Finally calculate lambda
//
double jacobian = (W2-protonMass2+Q2)/value;
lambda = jacobian*derivative;
} // end fall back solution
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
if (lambdaFromTable)
cout << "Info, lambda taken from table." << endl;
else
cout << "Info, lambda derived numerically from proton amplitude table" << endl;
}
//
// Check lambda value.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
//
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
cout << " Set to 0." << endl;
return 0;
}
if (std::isnan(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
cout << " Set to 0." << endl;
return 0;
}
return lambda;
}
//
// UPC only version
//
double CrossSection::logDerivateOfAmplitude(double t, double xpom) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfAmplitude(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// Usual numeric issues at boundaries.
// Subtracting an eps in log(xpom) does the trick.
//
if (xpom > mProtonTableCollection->maxX()) {
xpom = exp(log(xpom)-numeric_limits::epsilon());
}
if (xpom < mProtonTableCollection->minX()) {
xpom = exp(log(xpom)+numeric_limits::epsilon());
}
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
lambda = mProtonTableCollection->get(xpom, t, lambda_real, table);
if (!table) { // no lambda value from correct table, use fallback solution
lambdaFromTable = false;
(void) mProtonTableCollection->get(xpom, t, mean_A, table); // use obtained table from here on
if (!table) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid pointer to lookup table." << endl;
return 0;
}
//
// Here's the tricky part (see comments in non-UPC version above)
//
double theLogxpom = log(xpom);
double dlogxpom = table->binWidthX(); // assuming table is in log x ???????? FIX later
double maxLogxpom = log(table->maxX());
double derivative;
double hplus, hminus;
hplus = hminus = 0.5*dlogxpom;
hplus = min(hplus, fabs(maxLogxpom-theLogxpom));
hminus = min(hminus, fabs(theLogxpom-log(table->minX())));
hminus -= numeric_limits::epsilon();
hplus -= numeric_limits::epsilon();
if (hminus < 0) hminus = 0;
if (hplus < 0) hplus = 0;
if (hplus == 0 && hminus == 0) {
cout << "CrossSection::logDerivateOfAmplitude(): Warning, cannot find derivative." << endl;
return 0;
}
double a = table->get(exp(theLogxpom+hplus), t);
if (a <= 0) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << a << '.' << endl;
cout << " t=" << t << ", W=" << sqrt(exp(theLogxpom+hplus)) << endl;
return 0;
}
double b = table->get(exp(theLogxpom-hminus), t);
if (a <= 0) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << b << '.' << endl;
cout << " t=" << t << ", W=" << sqrt(exp(theLogxpom-hminus)) << endl;
return 0;
}
derivative = log(a/b)/(hplus+hminus);
//
// Finally calculate lambda
// Directly dlog(A)/-dlog(xpom) here.
//
lambda = -derivative;
}
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
if (lambdaFromTable)
cout << "Info, lambda taken from table." << endl;
else
cout << "Info, lambda derived numerically from proton amplitude table" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
}
//
// Check lambda value.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
//
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "xpom=" << xpom << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda = infinity for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
cout << " Set to 0." << endl;
return 0;
}
if (std::isnan(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
cout << " Set to 0." << endl;
return 0;
}
return lambda;
}
double CrossSection::logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_skew, table);
if (!table) {
//
// no lambda value from correct table, use fallback solution
// lambda_skew ~ lambda_real
//
lambda = logDerivateOfAmplitude(t, Q2, W2, pol) ; // This is an approximation in this case
} // end fall back solution
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
if (lambdaFromTable)
cout << "Info, lambda taken from table." << endl;
else
cout << "Info, lambda taken from logDerivateOfAmplitude as an approximation" << endl;
}
//
// Checking lambda.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
//
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
cout << " Set to 0." << endl;
return 0;
}
if (std::isnan(lambda)) {
cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
cout << " Set to 0." << endl;
return 0;
}
return lambda;
}
//
// UPC version
//
double CrossSection::logDerivateOfGluonDensity(double t, double xpom) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// Usual numeric issues at boundaries.
// Subtracting an eps in log(xpom) does the trick.
//
if (xpom > mProtonTableCollection->maxX()) {
xpom = exp(log(xpom)-numeric_limits::epsilon());
}
if (xpom < mProtonTableCollection->minX()) {
xpom = exp(log(xpom)+numeric_limits::epsilon());
}
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
lambda = mProtonTableCollection->get(xpom, t, lambda_skew, table);
if (!table) {
//
// no lambda value from correct table, use fallback solution
// lambda_skew ~ lambda_real
//
lambdaFromTable = false;
lambda = logDerivateOfAmplitude(t, xpom);
}
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
if (lambdaFromTable)
cout << "Info, lambda taken from table." << endl;
else
cout << "Info, lambda taken from logDerivateOfAmplitude as approximation." << endl;
}
//
// Check lambda value.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
//
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "xpom=" << xpom << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = infinity for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
cout << " Set to 0." << endl;
return 0;
}
if (std::isnan(lambda)) {
cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
cout << " Set to 0." << endl;
return 0;
}
return lambda;
}
double CrossSection::realAmplitudeCorrection(double lambda) const
{
//
// Correction factor for real amplitude contribution
//
double beta = tan(lambda*M_PI/2.);
double correction = 1 + beta*beta;
return correction;
}
double CrossSection::realAmplitudeCorrection(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
double correction = realAmplitudeCorrection(lambda);
// correction *= exp(-10*Kinematics::x(Q2, W2)); // damped
return correction;
}
//
// UPC version
//
double CrossSection::realAmplitudeCorrection(double t, double xpom) const
{
double lambda = logDerivateOfAmplitude(t, xpom);
double correction = realAmplitudeCorrection(lambda);
return correction;
}
double CrossSection::skewednessCorrection(double lambda) const
{
//
// Skewedness correction
//
double R = pow(2.,2*lambda+3)*TMath::Gamma(lambda+2.5)/(sqrt(M_PI)*TMath::Gamma(lambda+4));
double correction = R*R;
return correction;
}
double CrossSection::skewednessCorrection(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
double correction = skewednessCorrection(lambda);
// correction *= exp(-10*Kinematics::x(Q2, W2)); // damped
return correction;
}
//
// UPC version
//
double CrossSection::skewednessCorrection(double t, double xpom) const
{
double lambda = logDerivateOfAmplitude(t, xpom);
double correction = skewednessCorrection(lambda);
return correction;
}
Index: trunk/src/DipoleModelParameters.cpp
===================================================================
--- trunk/src/DipoleModelParameters.cpp (revision 377)
+++ trunk/src/DipoleModelParameters.cpp (revision 378)
@@ -1,474 +1,474 @@
//==============================================================================
// DipoleModelParameters.cpp
//
-// Copyright (C) 2016-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2016-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// $Date$
// $Author: ullrich $
//==============================================================================
#include "DipoleModelParameters.h"
#include "Settings.h"
#include "TableGeneratorSettings.h"
#include
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
DipoleModelParameters::DipoleModelParameters(Settings* settings)
{
mDipoleModelType = settings->dipoleModelType();
mDipoleModelParameterSetName = settings->dipoleModelParameterSetName();
mDipoleModelParameterSet = settings->dipoleModelParameterSet();
setupParameters();
}
DipoleModelParameters::DipoleModelParameters(DipoleModelType mtype, DipoleModelParameterSet pset) :
mDipoleModelType(mtype), mDipoleModelParameterSet(pset)
{
setupParameters();
}
void DipoleModelParameters::setDipoleModelType(DipoleModelType val)
{
mDipoleModelType = val;
setupParameters();
}
void DipoleModelParameters::setDipoleModelParameterSet(DipoleModelParameterSet val)
{
mDipoleModelParameterSet = val;
setupParameters();
}
DipoleModelType DipoleModelParameters::dipoleModelType() const {return mDipoleModelType;}
DipoleModelParameterSet DipoleModelParameters::dipoleModelParameterSet() const {return mDipoleModelParameterSet;}
void DipoleModelParameters::setup_bSat()
{
if (mDipoleModelParameterSet == KMW) {
// KMW paper (arXiv:hep-ph/0606272), Table 3
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
mQuarkMass[3] = 1.4; // c quark
mBG = 4.;
mMu02 = 1.17; // Gev^2
mLambdaG = 0.02;
mAg = 2.55;
mC = 4;
}
else if (mDipoleModelParameterSet == HMPZ) {
// Heikki Mantysaari an Pia Zurita, Phys.Rev. D98 (2018) 036002 (arXiv:1804.05311)
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.03; // u,d,s quarks
mQuarkMass[3] = 1.3528; // c quark
mBG = 4.;
mMu02 = 1.1; // Gev^2
mLambdaG = 0.08289;
mAg = 2.1953;
mC = 2.2894;
}
else if (mDipoleModelParameterSet == CUSTOM) {
if (mCustomParameters.size() < 10) {
cout << "DipoleModelParameters::setupParameters(): Error, require 10 custom parameters for bSAT when" << endl;
cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
exit(1);
}
mQuarkMass[0] = mCustomParameters[0];
mQuarkMass[1] = mCustomParameters[1];
mQuarkMass[2] = mCustomParameters[2];
mQuarkMass[3] = mCustomParameters[3];
mQuarkMass[4] = mCustomParameters[4];
mBG = mCustomParameters[5];
mMu02 = mCustomParameters[6];
mLambdaG = mCustomParameters[7];
mAg = mCustomParameters[8];
mC = mCustomParameters[9];
}
else {
cout << "DipoleModelParameters::setup_bSat(): Error, no known parameters for given dipole model" << endl;
cout << " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
void DipoleModelParameters::setup_bNonSat()
{
if (mDipoleModelParameterSet == KMW) {
// KT paper (arXiv:hep-ph/0304189v3), page 11
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
mQuarkMass[3] = 1.4;
mBG = 4.;
mMu02 = 0.8;
mLambdaG = -0.13;
mAg = 3.5;
mC = 4;
}
else if (mDipoleModelParameterSet == HMPZ) {
// Internal note by Heikki Mantysaari an Pia Zurita, arXiv pending
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.1516; // u,d,s quarks
mQuarkMass[3] = 1.3504;
mBG = 4.;
mMu02 = 1.1;
mLambdaG = -0.006657;
mAg = 3.0391;
mC = 4.2974;
}
else if (mDipoleModelParameterSet == CUSTOM) {
if (mCustomParameters.size() < 10) {
cout << "DipoleModelParameters::setupParameters(): Error, require 10 custom parameters for bNonSAT when" << endl;
cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
exit(1);
}
mQuarkMass[0] = mCustomParameters[0];
mQuarkMass[1] = mCustomParameters[1];
mQuarkMass[2] = mCustomParameters[2];
mQuarkMass[3] = mCustomParameters[3];
mQuarkMass[4] = mCustomParameters[4];
mBG = mCustomParameters[5];
mMu02 = mCustomParameters[6];
mLambdaG = mCustomParameters[7];
mAg = mCustomParameters[8];
mC = mCustomParameters[9];
}
else {
cout << "DipoleModelParameters::setup_bNonSat(): Error, no known parameters for given dipole model" << endl;
cout << " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
void DipoleModelParameters::setup_bCGC()
{
if (mDipoleModelParameterSet == KMW) {
// WK paper (arXiv:0712.2670), Table II
mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
mQuarkMass[3] = 1.4;
mKappa = 9.9;
mN0 = 0.558;
mX0 = 1.84e-6;
mLambda = 0.119;
mGammas = 0.46;
mBcgc = 7.5;
}
else if (mDipoleModelParameterSet == CUSTOM) {
if (mCustomParameters.size() < 10) {
cout << "DipoleModelParameters::setup_bCGC(): Error, require 10 custom parameters for bCGC when" << endl;
cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
exit(1);
}
mQuarkMass[0] = mCustomParameters[0];
mQuarkMass[1] = mCustomParameters[1];
mQuarkMass[2] = mCustomParameters[2];
mQuarkMass[3] = mCustomParameters[3];
mKappa = mCustomParameters[4];
mN0 = mCustomParameters[5];
mX0 = mCustomParameters[6];
mLambda = mCustomParameters[7];
mGammas = mCustomParameters[8];
mBcgc = mCustomParameters[9];
}
else {
cout << "DipoleModelParameters::setup_bCGC(): Error, no known parameters for given dipole model"
<< " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
void DipoleModelParameters::setupParameters()
{
TableGeneratorSettings *settings = TableGeneratorSettings::instance();
if (mDipoleModelParameterSet == CUSTOM) mCustomParameters = settings->dipoleModelCustomParameters();
//
// Init
//
mKappa = 0;
mN0 = 0;
mX0 = 0;
mLambda = 0;
mGammas = 0;
mBcgc = 0;
mBG = 0;
mMu02 = 0;
mLambdaG = 0;
mAg = 0;
mC = 0;
mBoostedGaussianR2_rho = 0;
mBoostedGaussianNL_rho = 0;
mBoostedGaussianNT_rho = 0;
mBoostedGaussianQuarkMass_rho = 0;
mBoostedGaussianR2_phi = 0;
mBoostedGaussianNL_phi = 0;
mBoostedGaussianNT_phi = 0;
mBoostedGaussianQuarkMass_phi = 0;
mBoostedGaussianR2_jpsi = 0;
mBoostedGaussianNL_jpsi = 0;
mBoostedGaussianNT_jpsi = 0;
mBoostedGaussianQuarkMass_jpsi = 0;
mBoostedGaussianR2_ups = 0;
mBoostedGaussianNL_ups = 0;
mBoostedGaussianNT_ups = 0;
mBoostedGaussianQuarkMass_ups = 0;
//
// b and t masses (not used, just for completeness)
//
// mQuarkMass[4] = 4.75; // b quark consistent with HMPZ
mQuarkMass[4] = 4.2; // b quark consistent with Upsilon wave function.
mQuarkMass[5] = 175.; // t quark consistent with HMPZ
//
// Parameters for boosted Gaussian wave function
//
setup_boostedGaussiansWaveFunction();
//
// Model parameters
//
if (mDipoleModelType == bSat) {
setup_bSat();
}
else if (mDipoleModelType == bNonSat) {
setup_bNonSat();
}
else if (mDipoleModelType == bCGC) {
setup_bCGC();
}
else {
cout << "DipoleModelParameters::setupParameters(): Error, no known parameters for given dipole model" << endl;
cout << " and requested parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianR2(int vm)
{
if (vm == 113)
return mBoostedGaussianR2_rho;
else if (vm == 333)
return mBoostedGaussianR2_phi;
else if (vm == 443)
return mBoostedGaussianR2_jpsi;
else if (vm == 553)
return mBoostedGaussianR2_ups;
else {
cout << "DipoleModelParameters::boostedGaussianR2(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianNL(int vm)
{
if (vm == 113)
return mBoostedGaussianNL_rho;
else if (vm == 333)
return mBoostedGaussianNL_phi;
else if (vm == 443)
return mBoostedGaussianNL_jpsi;
else if (vm == 553)
return mBoostedGaussianNL_ups;
else {
cout << "DipoleModelParameters::boostedGaussianNL(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianNT(int vm)
{
if (vm == 113)
return mBoostedGaussianNT_rho;
else if (vm == 333)
return mBoostedGaussianNT_phi;
else if (vm == 443)
return mBoostedGaussianNT_jpsi;
else if (vm == 553)
return mBoostedGaussianNT_ups;
else {
cout << "DipoleModelParameters::boostedGaussianNT(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
void DipoleModelParameters::setup_boostedGaussiansWaveFunction()
{
//
// Technical note:
// The Upsilon is a late addition with parameters coming from
// DKMM (arXiv:hep-ph/1610.06647). Heikki provided the more
// precise normalization constants for N_T and N_L.
// Neither KMW nor HMPZ provided Upsilon parameters, so
// results need to be verified with HERA data first.
// Note: All that is taken from DKMM is b-quark mass, the rest
// is determined by norm and decay width.
//
if (mDipoleModelParameterSet == KMW || mDipoleModelParameterSet == HMPZ) {
mBoostedGaussianR2_ups = 0.567;
mBoostedGaussianNT_ups = 0.481493;
mBoostedGaussianNL_ups = 0.480264 ;
mBoostedGaussianQuarkMass_ups = 4.2;
}
//
// rho, phi, and J/psi
//
if (mDipoleModelParameterSet == KMW) {
//
// KMW: bSat, bNonSat, and bCGC use the same parameters
// and also do not distingiosh between T and L.
//
mBoostedGaussianR2_rho = 12.9;
mBoostedGaussianNL_rho = 0.853;
mBoostedGaussianNT_rho = 0.911;
mBoostedGaussianQuarkMass_rho = 0.14;
mBoostedGaussianR2_phi = 11.2;
mBoostedGaussianNL_phi = 0.825;
mBoostedGaussianNT_phi= 0.919;
mBoostedGaussianQuarkMass_phi = 0.14;
mBoostedGaussianR2_jpsi = 2.3;
mBoostedGaussianNL_jpsi = 0.575;
mBoostedGaussianNT_jpsi = 0.578;
mBoostedGaussianQuarkMass_jpsi = 1.4;
}
else if (mDipoleModelParameterSet == HMPZ) {
if (mDipoleModelType == bSat) {
mBoostedGaussianR2_rho = 3.6376*3.6376;
mBoostedGaussianNL_rho = 0.8926;
mBoostedGaussianNT_rho = 0.9942;
mBoostedGaussianQuarkMass_rho = 0.03;
mBoostedGaussianR2_phi = 3.3922*3.3922;
mBoostedGaussianNL_phi = 0.8400;
mBoostedGaussianNT_phi = 0.9950;
mBoostedGaussianQuarkMass_phi = 0.03;
mBoostedGaussianR2_jpsi = 1.5070*1.5070;
mBoostedGaussianNL_jpsi = 0.5860;
mBoostedGaussianNT_jpsi = 0.5890;
mBoostedGaussianQuarkMass_jpsi = 1.3528;
}
else if (mDipoleModelType == bNonSat) {
mBoostedGaussianR2_rho = 3.5750*3.5750;
mBoostedGaussianNL_rho = 0.8467;
mBoostedGaussianNT_rho = 0.8978;
mBoostedGaussianQuarkMass_rho = 0.1516;
mBoostedGaussianR2_phi = 3.3530*3.3530;
mBoostedGaussianNL_phi = 0.8196;
mBoostedGaussianNT_phi = 0.9072;
mBoostedGaussianQuarkMass_phi = 0.1516;
mBoostedGaussianR2_jpsi = 1.5071*1.5071;
mBoostedGaussianNL_jpsi = 0.5868;
mBoostedGaussianNT_jpsi = 0.5899;
mBoostedGaussianQuarkMass_jpsi = 1.3504;
}
else if (mDipoleModelType == bCGC) {
cout << "DipoleModelParameters::setup_boostedGaussiansWaveFunction(): "
"Error, no HMPZ wave function parameters for CGC model. Stop." << endl;
exit(1);
}
}
else {
cout << "DipoleModelParameters::setup_boostedGaussiansWaveFunction(): Error, no known parameters for given dipole model "
<< " parmeter set "
<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
exit(1);
}
}
double DipoleModelParameters::boostedGaussianQuarkMass(int vm)
{
if (vm == 113)
return mBoostedGaussianQuarkMass_rho;
else if (vm == 333)
return mBoostedGaussianQuarkMass_phi;
else if (vm == 443)
return mBoostedGaussianQuarkMass_jpsi;
else if (vm == 553)
return mBoostedGaussianQuarkMass_ups;
else {
cout << "DipoleModelParameters::boostedGaussianQuarkMass(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
exit(1);
}
}
bool DipoleModelParameters::list(ostream& os)
{
const int fieldWidth = 32;
os << "\nDipole Model Parameters:" << endl;
os << setw(fieldWidth) << "Set: " << mDipoleModelParameterSetName << endl;
os << setw(fieldWidth) << "Quark masses: "
<< "u=" << mQuarkMass[0]
<< ", d=" << mQuarkMass[1]
<< ", s=" << mQuarkMass[2]
<< ", c=" << mQuarkMass[3]
<< ", b=" << mQuarkMass[4]
<< ", t=" << mQuarkMass[5] << endl;
os << setw(fieldWidth) << "BG: " << mBG << endl;
os << setw(fieldWidth) << "Mu02: " << mMu02 << endl;
os << setw(fieldWidth) << "LambdaG: " << mLambdaG << endl;
os << setw(fieldWidth) << "Ag: " << mAg << endl;
os << setw(fieldWidth) << "C: " << mC << endl;
os << setw(fieldWidth) << "Kappa: " << mKappa << endl;
os << setw(fieldWidth) << "N0: " << mN0 << endl;
os << setw(fieldWidth) << "X0: " << mX0 << endl;
os << setw(fieldWidth) << "Lambda: " << mLambda << endl;
os << setw(fieldWidth) << "Gammas: " << mGammas << endl;
os << setw(fieldWidth) << "Bcgc: " << mBcgc << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_rho: " << mBoostedGaussianR2_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_rho: " << mBoostedGaussianNL_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_rho: " << mBoostedGaussianNT_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_rho: " << mBoostedGaussianQuarkMass_rho << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_phi: " << mBoostedGaussianR2_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_phi: " << mBoostedGaussianNL_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_phi: " << mBoostedGaussianNT_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_phi: " << mBoostedGaussianQuarkMass_phi << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_jpsi: " << mBoostedGaussianR2_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_jpsi: " << mBoostedGaussianNL_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_jpsi: " << mBoostedGaussianNT_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_jpsi: " << mBoostedGaussianQuarkMass_jpsi << endl;
os << setw(fieldWidth) << "BoostedGaussianR2_ups: " << mBoostedGaussianR2_ups << endl;
os << setw(fieldWidth) << "BoostedGaussianNL_ups: " << mBoostedGaussianNL_ups << endl;
os << setw(fieldWidth) << "BoostedGaussianNT_ups: " << mBoostedGaussianNT_ups << endl;
os << setw(fieldWidth) << "BoostedGaussianQuarkMass_ups: " << mBoostedGaussianQuarkMass_ups << endl;
os << endl;
return true;
}
Index: trunk/src/Event.h
===================================================================
--- trunk/src/Event.h (revision 377)
+++ trunk/src/Event.h (revision 378)
@@ -1,77 +1,77 @@
//==============================================================================
// Event.h
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
#ifndef Event_h
#define Event_h
#include
#include
#include
#include "TLorentzVector.h"
#include "Enumerations.h"
using namespace std;
class Particle {
public:
int index; // starts at 0 equals index in particle vector
int pdgId; // particle ID according to PDG scheme
int status; // 1 = not decayed/final, 2 = decayed
TLorentzVector p;
vector parents;
vector daughters;
};
class Event {
public:
unsigned long eventNumber;
//
// Event kinematics
//
double Q2;
double W;
double t;
double x;
double s;
double y;
double xpom;
double beta;
//
// Event traits
//
GammaPolarization polarization;
DiffractiveMode diffractiveMode;
//
// List of particles in event.
// First two are always beam particles.
//
vector particles;
public:
void list(ostream& = cout) const;
};
#endif
Index: trunk/src/AlphaStrong.h
===================================================================
--- trunk/src/AlphaStrong.h (revision 377)
+++ trunk/src/AlphaStrong.h (revision 378)
@@ -1,166 +1,166 @@
//==============================================================================
// AlphaStrong.h
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
//
// Stand-alone code for alpha_s cannibalised (with permission)
// from Andreas Vogt's QCD-PEGASUS package (hep-ph/0408244).
// The running coupling alpha_s is obtained at N^mLO (m = 0,1,2,3)
// by solving the renormalisation group equation in the MSbar scheme
// by a fourth-order Runge-Kutta integration. Transitions from
// n_f to n_f+1 flavours are made when the factorisation scale
// mu_f equals the pole masses m_h (h = c,b,t). At exactly
// the thresholds m_{c,b,t}, the number of flavours n_f = {3,4,5}.
// The top quark mass should be set to be very large to evolve with
// a maximum of five flavours. The factorisation scale mu_f may be
// a constant multiple of the renormalisation scale mu_r. The input
// factorisation scale mu_(f,0) should be less than or equal to
// the charm quark mass. However, if it is greater than the
// charm quark mass, the value of alpha_s at mu_(f,0) = 1 GeV will
// be found using a root-finding algorithm.
//
// This is a rewrite of the original alphaS.f code originally
// written in F77. The documentation is largely left as provided
// in the Fortran version but includes the updated variable and
// function names. Data originally stored in common
// blocks and several static variables turned into data members.
// Some relics of the original F77 code remained such as array
// indexing starting at 1 (arrays are size +1 here).
// The old root finder used was replaced by a call to a root
// finder from the ROOT Math library.
//
// Original author: Graeme Watt
// C++ version: Thomas Ullrich (BNL)
//
//
// Example of usage.
// The constructor takes the following arguments:
//
// order // perturbative order (N^mLO,m=0,1,2,3) - always 0 in Sartre
// fr2 // ratio of mu_f^2 to mu_r^2
// mur // input mu_r in GeV
// asmur // input value of alpha_s at mu_r
// mc // charm quark mass
// mb // bottom quark mass
// mt 10 // top quark mass
//
// AlphaStrong myAlphaS(order, fr2, mur, asmur, mc, mb, mt);
//
// Then get alpha_s at a renormalisation scale mu_r with:
//
// mu_r = 10.; // scale in GeV
// double result = myAlphaS.at(mu_r*mu_r); // argument in GeV^2
//
// Note that the argument takes the square of the renormilaztion scale
// in units of GeV^2.
//
// The default arguments of the constructor are tailored for
// Sartre and work equally well with KMW and HMPZ parameters.
//==============================================================================
#ifndef AlphaStrong_h
#define AlphaStrong_h
#include
#include "Math/IFunction.h"
using namespace std;
class AlphaStrong {
public:
AlphaStrong(int order = 0, double fr2 = 1,
double mur = 91.1876, double asmur = 0.1183,
double mc = 1.3528, double mb = 4.75, double mt = 175);
double at(double Q2); // scale in GeV^2
int order() const; // used order
double massCharm() const; // used charm quark mass
double massBottom() const; // used bottom quark mass
double massTop() const; // used bottom quark mass
double ratioFR2() const; // used ratio of mu_f^2 to mu_r^2
double alphasAtMuR() const; // used input value of alpha_s at mu_r
private:
class rootFinderFormula : public ROOT::Math::IBaseFunctionOneDim {
public:
double DoEval(double v) const {return mParent->findR0(v);}
ROOT::Math::IBaseFunctionOneDim* Clone() const {return new rootFinderFormula();}
AlphaStrong* mParent;
};
private:
double findR0(double asi);
double asnf1 (double asnf, double logrh, int nf);
double as(double r2, double r20, double as0, int nf);
double funBeta1(double, int);
double funBeta2(double, int);
double funBeta3(double, int);
double sign(double, double);
void evolution (double mc2, double mb2, double mt2);
void initR0(int order, double fr2, double r0, double asi, double mc, double mb, double mt);
void betaFunction();
private:
double mFR2;
double mScale;
double mAsScale;
double mMassC;
double mMassB;
double mMassT;
double mR0;
double mOrder;
double mZeta[7];
double mCF;
double mCA;
double mTR;
double mAS0;
double mM20;
double mLogFR;
double mASC;
double mM2C;
double mASB;
double mM2B;
double mAST;
double mM2T;
double mBeta0[7];
double mBeta1[7];
double mBeta2[7];
double mBeta3[7];
double mCCMCoefficients[4][4];
double mCCMCoefficientsI30;
double mCCMCoefficientsI31;
double mCCMCoefficientsF30;
double mCCMCoefficientsF31;
int mWasCalled;
int mPertubativeOrder;
int mIntegrationSteps;
int mVarFlavourNumScheme;
int mNumFlavorsFFNS;
};
inline int AlphaStrong::order() const {return mOrder;}
inline double AlphaStrong::massCharm() const {return mMassC;}
inline double AlphaStrong::massBottom() const {return mMassB;}
inline double AlphaStrong::massTop() const {return mMassT;}
inline double AlphaStrong::ratioFR2() const {return mFR2;}
inline double AlphaStrong::alphasAtMuR() const {return mAsScale;}
inline double AlphaStrong::sign(double a, double b) {return ((b < 0) ? -fabs(a) : fabs(a));}
#endif
Index: trunk/src/ExclusiveFinalStateGenerator.cpp
===================================================================
--- trunk/src/ExclusiveFinalStateGenerator.cpp (revision 377)
+++ trunk/src/ExclusiveFinalStateGenerator.cpp (revision 378)
@@ -1,607 +1,607 @@
//==============================================================================
// ExclusiveFinalStateGenerator.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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 "ExclusiveFinalStateGenerator.h"
#include "EventGeneratorSettings.h"
#include "Kinematics.h"
#include "Event.h"
#include "Math/BrentRootFinder.h"
#include "Math/GSLRootFinder.h"
#include "Math/RootFinderAlgorithms.h"
#include "Math/IFunction.h"
#include
#include
#include
#include
#define PR(x) cout << #x << " = " << (x) << endl;
//-------------------------------------------------------------------------------
//
// Helper class needed to find root in
// ExclusiveFinalStateGenerator::generate()
//
//-------------------------------------------------------------------------------
class ScatteredProtonEnergyFormula : public ROOT::Math::IBaseFunctionOneDim
{
public:
double DoEval(double) const;
ROOT::Math::IBaseFunctionOneDim* Clone() const;
void calculateValidRange(double&, double&);
public:
double mT;
double mVmMass2;
double mMY2;
double mPhi; // azimuthal angle for scattered proton
TLorentzVector mProtonIn;
TLorentzVector mVirtualPhoton;
};
ROOT::Math::IBaseFunctionOneDim* ScatteredProtonEnergyFormula::Clone() const
{
return new ScatteredProtonEnergyFormula();
}
double ScatteredProtonEnergyFormula::DoEval(double Ep) const
{
double m2 = mProtonIn.M2();
double pzp = (mT - m2 - mMY2 + 2*mProtonIn.E() * Ep)/(2*mProtonIn.Pz());
double term = Ep*Ep-pzp*pzp-mMY2;
if (term < 0) return 99999; // out of kinematically allowed range
double ptp = sqrt(term);
TLorentzVector p_out(ptp*cos(mPhi), ptp*sin(mPhi), pzp, Ep);
double f = (mVirtualPhoton + mProtonIn - p_out)*(mVirtualPhoton + mProtonIn - p_out) - mVmMass2;
return f;
}
void ScatteredProtonEnergyFormula::calculateValidRange(double& lower, double& upper)
{
double m2 = mProtonIn.M2();
double term1 = mT-m2-mMY2;
double termA = mProtonIn.E()*term1;
double termB = sqrt(mProtonIn.Pz()*mProtonIn.Pz()*(term1*term1-4*m2*mMY2));
double termC = -2*m2;
lower = (termA+termB)/termC;
upper = (termA-termB)/termC;
if (lower > upper) swap(lower, upper);
lower += numeric_limits::epsilon();
upper -= numeric_limits::epsilon();
}
//-------------------------------------------------------------------------------
//
// Implementation of ExclusiveFinalStateGenerator
//
//-------------------------------------------------------------------------------
ExclusiveFinalStateGenerator::ExclusiveFinalStateGenerator() {/* no op */}
ExclusiveFinalStateGenerator::~ExclusiveFinalStateGenerator() {/* no op */}
bool ExclusiveFinalStateGenerator::generate(int id, double t, double y, double Q2,
bool isIncoherent, int A, Event *event)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
//
// The beam particles must be present in the event list
//
int ePos = -1;
int hPos = -1;
bool parentsOK = true;
if (event->particles.size() == 2) {
if (abs(event->particles[0].pdgId) == 11) {
ePos = 0;
hPos = 1;
}
else if (abs(event->particles[1].pdgId) == 11) {
ePos = 1;
hPos = 0;
}
else
parentsOK = false;
}
else
parentsOK = false;
if (!parentsOK) {
cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
return false;
}
//
// Store arguments locally
// (Some could also be obtained from the event structure)
//
mA = A;
mT = t;
if (mT > 0) mT = -mT; // ensure t<0
mQ2 = Q2;
mY = y;
mIsIncoherent = isIncoherent;
mElectronBeam = event->particles[ePos].p;
mHadronBeam = event->particles[hPos].p;
mMassVM = settings->lookupPDG(id)->Mass();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
//
// Constants
//
double const twopi = 2*M_PI;
double const hMass2 = mHadronBeam.M2();
//
// Incoherent diffarction
//
// Generate hadron dissociation mass according to
// dN/dM2 ~ 1/M2. Lower bound is of course the hadron
// mass and upper bound is some arbitrary value (for now).
//
// Note that we calculate and quote eA kinematics always in
// units of 'per nucleon'. Our model of incoherence is that the
// difference of the diffractive mass of one (1) proton out
// of the nucleus gives the final excitation energy E*.
// Hence we have to calculate E* and divide it by A to keep
// the kinematic consistent.
//
if (mIsIncoherent && mA > 1) {
const double lower = hMass2;
const double upper = 9; // GeV2
mMY2 = lower*upper/(upper - rndm->Uniform()*(upper-lower));
double MY_per_nucleon = (sqrt(hMass2)*(mA-1) + sqrt(mMY2))/mA;
mMY2 = MY_per_nucleon*MY_per_nucleon;
if (mMY2 < hMass2) mMY2 = hMass2;
}
else {
mMY2 = hMass2;
}
//
// Re-engineer scattered electron
//
// e'=(E', pt', pz') -> 3 unknowns
//
// Three equations:
// 1: me*me=E'*E'-pt'*pt'-pz'*pz'
// 2: Q2=-(e-e')^2=-2*me*me + 2*(E*E'-pz*pz')
// 3: W2=(P+e-e')^2=mp2+2*me2+2*(Ep*E-Pz*pz)-2*(Ep*E'-Pz*pz')-2*(E*E'-pz*pz')
//
double Ee=mElectronBeam.E();
double Pe=mElectronBeam.Pz();
double Ep=mHadronBeam.E();
double Pp=mHadronBeam.Pz();
double W=event->W;
double W2=W*W;
// Take masses from the beams in case they are not actually electrons or protons
double me2=mElectronBeam.M2();
double mp2=mHadronBeam.M2();
//
// What we want for each particle:
//
double E, pz, pt, px, py, phi;
//
// Equations 2 and 3 yield:
//
E = Pe*(W2-mp2-2*Ee*Ep) + (Pp+Pe)*Q2 + 2*Pe*Pe*Pp + 2*me2*Pp;
E /= 2*(Ee*Pp-Ep*Pe);
pz = Ee*(W2-mp2) + (Ep+Ee)*Q2 + 2*Ee*Pe*Pp + 2*Ep*me2 - 2*Ee*Ee*Ep;
pz /= 2*(Ee*Pp-Ep*Pe);
//
// Equation 1:
//
pt = sqrt(E*E-pz*pz-me2);
phi = rndm->Uniform(twopi);
TLorentzVector theScatteredElectron(pt*sin(phi), pt*cos(phi), pz, E);
//
// Re-engineer virtual photon
//
// gamma=E-E'
E=mElectronBeam.E()-theScatteredElectron.E();
pz=mElectronBeam.Pz()-theScatteredElectron.Pz();
px=mElectronBeam.Px()-theScatteredElectron.Px();
py=mElectronBeam.Py()-theScatteredElectron.Py();
TLorentzVector theVirtualPhoton = TLorentzVector(px, py, pz, E);
//
// Re-engineer scattered proton/dissociated proton
//
// No analytic solution. Need to run a root finder that does
// not need derivates but uses a bracketing algorithm (Brent).
// Correct brackets are crucial since ScatteredProtonEnergyFormula
// produces sqrt(-x) if outside the kinematically allowed range (it
// actually catches it and returns a large positive number, 0 doesn't
// work).
//
//
// Setup formula to solve root
//
phi = rndm->Uniform(twopi);
ScatteredProtonEnergyFormula formula;
formula.mT = mT;
formula.mVmMass2 = mMassVM*mMassVM;
formula.mPhi = phi;
formula.mProtonIn = mHadronBeam;
formula.mVirtualPhoton = theVirtualPhoton;
formula.mMY2 = mMY2;
//
// Find correct brackets to start with
//
double lower, upper;
formula.calculateValidRange(lower, upper);
if (upper > mHadronBeam.E() + theVirtualPhoton.E()) // limit excessive values
upper = mHadronBeam.E() + theVirtualPhoton.E(); // make it easier for Brent
//
// Run root finder
//
ROOT::Math::BrentRootFinder rootfinder;
rootfinder.SetFunction(formula, lower, upper);
rootfinder.Solve(10000, 0, 1.e-12);
E = rootfinder.Root();
if (/* rootfinder.Status() || */ fabs(formula(E)) > 1e-6) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, cannot find root. No final state defined." << endl;
return false;
}
//
// Outgoing proton (hadron) system
//
pz = (mT- hMass2 - mMY2 + 2*mHadronBeam.E()*E)/(2*mHadronBeam.Pz());
pt = sqrt(E*E-pz*pz-mMY2);
px = pt*cos(phi);
py = pt*sin(phi);
TLorentzVector theScatteredProton(px, py, pz, E);
//
// Finally the vector meson
//
TLorentzVector theVectorMeson((mHadronBeam + mElectronBeam) - (theScatteredElectron + theScatteredProton));
//
// Check for numerical glitches
//
if (!isValid(theScatteredElectron)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered electron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theScatteredProton)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered hadron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theVectorMeson)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, vector meson 4-vector is invalid." << endl;
return false;
}
//
// Add particles to event record
//
event->particles.resize(2+5);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int vm = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[vm].index = vm;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[vm].p = theVectorMeson;
event->particles[pomeron].p = theScatteredProton - mHadronBeam;
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[vm].pdgId = id;
event->particles[pomeron].pdgId = 990;
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[vm].status = 1;
event->particles[pomeron].status = 2;
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[vm].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(vm);
event->particles[gamma].daughters.push_back(pomeron);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
return true;
}
//
// UPC version
//
// Note: We call the beam particle that emits the photon "electron"
// even if its a proton/nucleus
//
bool ExclusiveFinalStateGenerator::generate(int id, double t, double xpom,
bool isIncoherent, int A, Event *event)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
//
// The beam particles must be present in the event list
//
int ePos = 0;
int hPos = 1;
bool parentsOK = true;
if (event->particles.size() != 2)
parentsOK = false;
if (!parentsOK) {
cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
return false;
}
//
// Store arguments locally
// (Some could also be obtained from the event structure)
//
mA = A;
mT = t;
if (mT > 0) mT = -mT; // ensure t<0
mIsIncoherent = isIncoherent;
mElectronBeam = event->particles[ePos].p;
mHadronBeam = event->particles[hPos].p;
mMassVM = settings->lookupPDG(id)->Mass();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mXp=xpom;
//
// Constants
//
double const twopi = 2*M_PI;
double const hMass2 = mHadronBeam.M2();
//
// Incoherent diffraction
//
// Generate hadron dissociation mass according to
// dN/dM2 ~ 1/M2. Lower bound is of course the hadron
// mass and upper bound is some arbitrary value (for now).
//
// Note that we calculate and quote eA kinematics always in
// units of 'per nucleon'. Our model of incoherence is that the
// difference of the diffractive mass of one (1) proton out
// of the nucleus gives the final excitation energy E*.
// Hence we have to calculate E* and divide it by A to keep
// the kinematic consistent.
//
if (mIsIncoherent && mA > 1) {
const double lower = hMass2;
const double upper = 9; // GeV2
mMY2 = lower*upper/(upper - rndm->Uniform()*(upper-lower));
double MY_per_nucleon = (sqrt(hMass2)*(mA-1) + sqrt(mMY2))/mA;
mMY2 = MY_per_nucleon*MY_per_nucleon;
if (mMY2 < hMass2) mMY2 = hMass2;
}
else {
mMY2 = hMass2;
}
//
// Check for smallest allowed xpom
//
double xpom_min=Kinematics::xpomMin(mMassVM, mT, mHadronBeam, mElectronBeam, mMY2-hMass2);
if(xpom < xpom_min){
if (settings->verboseLevel() > 2) cout<<"xpom = "<verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered hadron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theVectorMeson)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, vector meson 4-vector is invalid." << endl;
return false;
}
//
// Add particles to event record
//
event->particles.resize(2+5);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int vm = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[vm].index = vm;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[vm].p = theVectorMeson;
event->particles[pomeron].p = thePomeron;
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[vm].pdgId = id;
event->particles[pomeron].pdgId = 990;
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[vm].status = 1;
event->particles[pomeron].status = 2;
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[vm].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(vm);
event->particles[gamma].daughters.push_back(pomeron);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
//fill event structure
double y=mHadronBeam*theVirtualPhoton/(mHadronBeam*mElectronBeam);
double W2=(theVirtualPhoton+mHadronBeam).M2();
mQ2=-theVirtualPhoton.M2();
event->Q2=mQ2;
event->W=sqrt(W2);
event->y=y;
return true;
}
Index: trunk/src/BreakupProduct.cpp
===================================================================
--- trunk/src/BreakupProduct.cpp (revision 377)
+++ trunk/src/BreakupProduct.cpp (revision 378)
@@ -1,37 +1,37 @@
//==============================================================================
// BreakupProduct.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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 "BreakupProduct.h"
#include
ostream & operator<<(ostream& os, const BreakupProduct& p)
{
ios::fmtflags fmt = os.flags(); // store io flags
os << setw(5) << right << p.name << " (A=" << p.A << ",Z=" << p.Z << ") \tid="
<< setw(11) << left << p.pdgId << " time=" << setprecision(3) << setw(10) << left << p.emissionTime
<< "\t p=(" << p.p.Px() << ", " << p.p.Py() << ", " << p.p.Pz() << ", " << p.p.E() << ')';
os.flags(fmt); // restore io flags
return os;
}
Index: trunk/src/DglapEvolution.cpp
===================================================================
--- trunk/src/DglapEvolution.cpp (revision 377)
+++ trunk/src/DglapEvolution.cpp (revision 378)
@@ -1,510 +1,510 @@
//==============================================================================
// DglapEvolution.h
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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$
//==============================================================================
#include "DglapEvolution.h"
#include "TableGeneratorSettings.h"
#include "DipoleModelParameters.h"
#include
#include
#include
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
DglapEvolution* DglapEvolution::mInstance = 0;
DglapEvolution::DglapEvolution()
{
DipoleModelParameters parameters(TableGeneratorSettings::instance());
//
// For speed purposes the key parameters are held as data member.
// Here we assign the proper ones depending on the model and the
// parameters set choosen.
//
mAg = parameters.Ag();
mLambdaG = parameters.lambdaG();
mMu02 = parameters.mu02();
//
// Define alpha_s functor
//
mAlphaStrong = new AlphaStrong;
//
// Set alpha_s at c, b, t mass
//
mMC = mAlphaStrong->massCharm();
mMB = mAlphaStrong->massBottom();
mMT = mAlphaStrong->massTop();
const double FourPi = 4.*M_PI;
mALPC = mAlphaStrong->at(mMC*mMC)/FourPi;
mALPB = mAlphaStrong->at(mMB*mMB)/FourPi;
mALPT = mAlphaStrong->at(mMT*mMT)/FourPi;
mALPS = mAlphaStrong->at(mMu02)/FourPi;
//
// Initialization of support points in n-space and weights for the
// Gauss quadrature and of the anomalous dimensions for the RG
// evolution at these n-values.
//
//
// Weights and support points for nomalized 8 point gauss quadrature
//
double wz[9] = {0, 0.101228536290376,0.222381034453374,0.313706645877887,
0.362683783378362,0.362683783378362,0.313706645877887,
0.222381034453374,0.101228536290376};
double zs[9] = {0, -0.960289856497536,-0.796666477413627,-0.525532409916329,
-0.183434642495650,0.183434642495650,0.525532409916329,
0.796666477413627,0.960289856497536};
//
// Integration contour parameters
//
double down[18] = {0, 0., 0.5, 1., 2., 3., 4., 6., 8.,
1.e1, 1.2e1, 1.5e1, 1.8e1, 2.1e1, 2.4e1, 2.7e1, 3.e1, 3.3e1};
double up[18];
mC = 0.8;
double phi = M_PI * 3./4.;
double co = cos(phi);
double si = sin(phi);
mCC = complex(co, si);
for (int i=1; i <=16; i++) up[i] = down[i+1];
up[17] = 36.;
//
// Support points and weights for the gauss integration
// (the factor (up-down)/2 is included in the weights)
//
int k = 0;
double sum, diff, z;
for (int i=1; i <=17; i++) {
sum = up[i] + down[i];
diff = up[i] - down[i];
for (int j=1; j <=8; j++) {
k++;
z = 0.5 * (sum + diff * zs[j]);
mWN[k] = diff / 2.* wz[j];
mN[k] = complex(mC+co*z+1.,si*z);
}
}
anom();
//
// Defaults for lookup table
//
mTableMinX = 1.e-10;
mTableMaxX = 1;
mTableMinQ2 = 1.;
mTableMaxQ2 = 1e6;
mLookupTableIsFilled = false;
mUseLookupTable = false;
mNumberOfNodesInX = mNumberOfNodesInQ2 = 0;
}
DglapEvolution::~DglapEvolution()
{
if (mAlphaStrong) delete mAlphaStrong;
if (mLookupTableIsFilled) {
for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
delete [] mLookupTable[i];
delete [] mLookupTable;
}
}
double DglapEvolution::alphaSxG(double x, double Q2)
{
double val = xG(x, Q2);
return val*mAlphaStrong->at(Q2);
}
double DglapEvolution::xG(double x, double Q2)
{
double result = 0;
if (!mUseLookupTable) {
result = xG_Engine(x, Q2);
}
else if (mUseLookupTable && !mLookupTableIsFilled) {
cout << "DglapEvolution::xG(): Warning, use of lookup table requested\n"
<< " but table is not setup. First use \n"
<< " generateLookupTable() to fill table.\n"
<< " Will fall back to numeric calculation."
<< endl;
result = xG_Engine(x, Q2);
}
else
result = xG_Interpolator(x, Q2);
return result;
}
double DglapEvolution::xG_Engine(double x, double Q2)
{
//
// These are provided by AlphaStrong ensuring
// consistency between evolution and alpha_s.
//
if (mAlphaStrong->order() != 0) {
cout << "DglapEvolution::xG_Engine(): Fatal error, alpha_s is in order "
<< mAlphaStrong->order()
<< " but DglapEvolution only support order = 0. Stop here." << endl;
exit(1);
}
double alpq = mAlphaStrong->at(Q2)/(4*M_PI);
//
// Q2 and x dependent quantities
//
double ax = log(x);
double ex = exp(-mC * ax);
//
// integration length parameter for the mellin inversion
//
const int nmax = 136;
//
// Gluon density and output
//
complex fn[137];
reno(fn, alpq, nmax, mAg, mLambdaG);
long double fun = 0; // long double is needed
long double fz;
complex xnm,cex;
for (int i=1; i <= nmax; i++) {
xnm = (mC - mN[i]+1.) * ax;
cex = exp(xnm) / M_PI * mCC;
fz = imag(fn[i]*cex);
fun = fun + mWN[i] * fz;
}
double pa = fun * ex;
return pa;
}
void DglapEvolution::generateLookupTable(int nx, int nq2)
{
//
// Delete old lookup table (if it exists)
//
if (mLookupTableIsFilled) {
for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
delete [] mLookupTable[i];
delete [] mLookupTable;
}
mNumberOfNodesInX = nx;
mNumberOfNodesInQ2 = nq2;
//
// Create new table
//
mLookupTable = new double*[mNumberOfNodesInQ2];
for(unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
mLookupTable[i] = new double[mNumberOfNodesInX];
//
// Fill lookup table
//
int printEach = mNumberOfNodesInQ2*mNumberOfNodesInX/10;
int nCount = 0;
cout << "DglapEvolution: generating lookup table "; cout.flush();
for (unsigned int i = 0; i < mNumberOfNodesInQ2; i++) {
double Q2 = pow(mTableMaxQ2/mTableMinQ2,i*1./(mNumberOfNodesInQ2-1.))*mTableMinQ2;
for (unsigned int j = 0; j < mNumberOfNodesInX; j++) {
double x = pow(mTableMaxX/mTableMinX,j*1./(mNumberOfNodesInX-1.))*mTableMinX;
mLookupTable[i][j] = xG_Engine(x, Q2);
nCount++;
if (nCount%printEach == 0) cout << '.'; cout.flush();
}
}
cout << " done" << endl;
mLookupTableIsFilled = true;
}
double DglapEvolution::luovi(double f[4], double arg[4], double z)
{
double cof[4];
for (int i=0; i < 4; i++ ) cof[i]=f[i];
for (int i=0; i < 3 ; i++) {
for (int j=i; j < 3 ; j++) {
int jndex = 2 - j;
int index = jndex + 1 + i;
cof[index]= (cof[index]-cof[index-1])/(arg[index]-arg[jndex]);
}
}
double sum=cof[3];
for (int i=0; i < 3; i++ ) {
int index = 2 - i;
sum = (z-arg[index])*sum + cof[index];
}
return sum;
}
double DglapEvolution::xG_Interpolator(double x, double Q2)
{
int q2steps = mNumberOfNodesInQ2-1;
int xsteps = mNumberOfNodesInX-1;
//
// Position in the Q2 grid
//
double realq = q2steps * log(Q2/mTableMinQ2)/log(mTableMaxQ2/mTableMinQ2);
int n_q2 = static_cast(realq);
if (n_q2 <= 0) {n_q2 = 1;}
if (n_q2 > q2steps-2) {n_q2 = n_q2-2;}
//
// Position in the x grid
//
double realx = xsteps * log(x/mTableMinX)/log(mTableMaxX/mTableMinX);
int n_x = static_cast(realx);
if (n_x <= 0) { n_x = 1;}
if (n_x > xsteps-2){ n_x = n_x-2;}
//
// Starting the interpolation
//
double arg[4];
for (int i=0; i < 4; i++) { arg[i] = n_x-1+i;}
double fu[4], fg[4];
for (int i=0; i < 4; i++) {
fu[0] = mLookupTable[n_q2-1+i][n_x-1];
fu[1] = mLookupTable[n_q2-1+i][n_x];
fu[2] = mLookupTable[n_q2-1+i][n_x+1];
fu[3] = mLookupTable[n_q2-1+i][n_x+2];
fg[i] = luovi(fu,arg,realx);
}
for (int i=0; i < 4; i++) { arg[i] = n_q2-1+i;}
return luovi(fg, arg, realq);
}
void DglapEvolution::reno(complex *fn, double alpq, int nmax, double ag, double lambdag)
{
//
// Mellin-n space Q**2 - evolution of the gluon at LO
//
// The moments are calculated on an array of moments, mN, suitable
// for a (non-adaptive) Gauss quadrature.
//
// Currently this takes the simplest possible fit form:
// xg = A_g x^(-lambdag) (1-x)^(5.6), following Amir&Raju
//
for (int k1 = 1; k1 <= nmax; k1++) {
//
// Input moments of the parton densities
// at the low scale
//
complex xn = mN[k1];
complex gln = ag * (1.0 / (xn + 5.0 - lambdag)
- 6.0 / (xn + 4.0 - lambdag)
+ 15.0 / (xn + 3.0 - lambdag)
- 20.0 / (xn + 2.0 -lambdag)
+ 15.0 / (xn + 1.0 - lambdag)
- 6.0 / (xn - lambdag)
+ 1.0 / (xn - lambdag - 1.0));
int f;
double xl, s, alp;
complex ep, gl;
if (alpq >= mALPC) { // evolution below the charm threshold
f = 3;
xl = mALPS / alpq;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
}
else if ((alpq < mALPC) && (alpq >= mALPB)) { // between thresholds
f = 3;
xl = mALPS / mALPC;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
f = 4;
xl = mALPC / alpq;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
}
else if (alpq < mALPB) { // above bottom threshold
f = 3;
xl = mALPS / mALPC;
s = log (xl);
ep = exp (- mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
f = 4;
alp = mALPB;
xl = mALPC / mALPB;
s = log (xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
f = 5;
xl = mALPB / alpq;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
}
fn[k1] = gln;
}
}
void DglapEvolution::anom() {
//
// Anomalous dimensions for leading order evolution of parton densities.
// The moments are calculated on an externally given array of mellin
// moments, mN, suitable for a (non-adaptive) quadrature.
//
// Present version: the number of active flavours in the factorization
// is fixed to ff=3, in the beta-function it varies between f=3 and
// f=5. The dimension of the moment array is 136.
//
double b0, b02;
complex ggi, ggf;
complex xn, xap;
complex gg;
for (int k1=1; k1 <= 136; k1++) {
xn = mN[k1];
anCalc(ggi, ggf, xn);
for (int k2=3; k2 <= 5; k2++) {
double f = k2;
// anomalous dimensions and related quantities in leading order
b0 = 11.- 2./3.* f;
b02 = 2.* b0;
gg = ggi + f * ggf;
xap = gg / b02;
mAP[k1][k2] = xap;
}
}
}
void DglapEvolution::anCalc(complex& ggi,
complex& ggf, complex& xn)
{
complex xn1 = xn + 1.;
complex xn2 = xn + 2.;
complex xnm = xn - 1.;
//
// Leading order
//
complex cpsi = psiFunction(xn1) + 0.577216;
ggi = -22.- 24./(xn * xnm) - 24./(xn1 * xn2) + 24.* cpsi;
ggf = 4./3.;
}
complex DglapEvolution::psiFunction(complex z)
{
//
// psi - function for complex argument
//
complex sub;
while (real(z) < 10) {
sub = sub - 1./z;
z += 1.;
}
complex rz = 1./z;
complex dz = rz * rz;
complex result = sub + log(z) - rz/2.- dz/2520. * ( 210.+ dz * (-21.+10.*dz ));
return result;
}
double DglapEvolution::logDerivative(double x, double Q2)
{
double alpq = mAlphaStrong->at(Q2)/(4*M_PI);
//
// Q2 and x dependent quantities
//
double ax = log(x);
double ex = exp(-mC * ax);
//
// integration length parameter for the mellin inversion
//
const int nmax = 136;
//
// Gluon density and output
//
complex fn[137];
reno(fn, alpq, nmax, mAg, mLambdaG);
long double fun = 0; // long double is needed
long double fz;
complex xnm,cex;
for (int i=1; i <= nmax; i++) {
xnm = (mC - mN[i]+1.) * ax;
cex = exp(xnm) / M_PI * mCC;
fz = imag(fn[i]*cex);
fun = fun + mWN[i] * fz;
}
double glue = fun * ex;
fun = 0;
for (int i=1; i <= nmax; i++) {
xnm = (mC - mN[i]+1.) * ax;
cex = exp(xnm) / M_PI * mCC;
fz = imag(fn[i]*cex*mN[i]);
fun = fun + mWN[i] * fz;
}
double pa = fun * ex;
double lambda = -(1-pa/glue);
return lambda;
}
Index: trunk/src/DipoleModel.cpp
===================================================================
--- trunk/src/DipoleModel.cpp (revision 377)
+++ trunk/src/DipoleModel.cpp (revision 378)
@@ -1,313 +1,313 @@
//==============================================================================
// DipoleModel.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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
#include "DipoleModel.h"
#include "TableGeneratorSettings.h"
#include "DglapEvolution.h"
#include "Constants.h"
#include "TFile.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH2F.h"
#include "TF1.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
DipoleModel::DipoleModel()
{
mConfigurationExists=false;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = settings->A();
mNucleus.init(A);
mIsInitialized=true;
mParameters = 0;
}
DipoleModel::~DipoleModel()
{
delete mParameters;
}
const TableGeneratorNucleus* DipoleModel::nucleus() const { return &mNucleus; }
bool DipoleModel::configurationExists() const { return mConfigurationExists; }
double DipoleModel::bDependence(double) { return 0; }
double DipoleModel::bDependence(double, double) { return 0; }
double DipoleModel::dsigmadb2ep(double, double, double) { return 0;}
//***********bSat:*****************************************************
DipoleModel_bSat::DipoleModel_bSat()
{
mBDependence = 0;
mSigma_ep_LookupTable = 0;
//
// Set the parameters. Note that we enforce here the bSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bSat, settings->dipoleModelParameterSet());
}
DipoleModel_bSat::~DipoleModel_bSat()
{
delete mSigma_ep_LookupTable;
delete mBDependence;
}
DipoleModel_bSat& DipoleModel_bSat::operator=(const DipoleModel_bSat& dp)
{
if (this != &dp) {
delete mBDependence;
delete mSigma_ep_LookupTable;
DipoleModel::operator=(dp);
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable));
mSigma_ep_LookupTable->SetDirectory(0);
}
return *this;
}
DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp)
{
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
}
void DipoleModel_bSat::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = mNucleus.A();
string path=settings->bSatLookupPath();
ostringstream filename;
filename.str("");
filename << path << "/bSat_bDependence_A" << A <<".root";
ifstream ifs(filename.str().c_str());
if (!ifs) {
cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl;
cout << "Stopping." << endl;
exit(1);
}
TFile* lufile= new TFile(filename.str().c_str());
ostringstream histoName;
histoName.str( "" );
histoName << "Configuration_" << iConfiguration;
lufile->GetObject( histoName.str().c_str(), mBDependence );
mBDependence->SetDirectory(0);
lufile->Close();
mConfigurationExists=true;
}
double DipoleModel_bSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
double bDep=bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
double DipoleModel_bSat::bDependence(double b, double phi)
{
double result = mBDependence->Interpolate(b, phi);
return result;
}
double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
{
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe){
xprobe*=1.;
double sigmap = mSigma_ep_LookupTable->Interpolate(r);
int A=nucleus()->A();
double TA=nucleus()->T(b)/A;
double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) );
return result;
}
void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe)
{
double rbRange=3.*nucleus()->radius();
TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2);
mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange);
dsigmaForIntegration->SetNpx(1000);
for (int iR=1; iR<=1000; iR++) {
double r=mSigma_ep_LookupTable->GetBinCenter(iR);
dsigmaForIntegration->SetParameter(0, r);
dsigmaForIntegration->SetParameter(1, xprobe);
double result=dsigmaForIntegration->Integral(0, rbRange);
mSigma_ep_LookupTable->SetBinContent(iR, result);
}
delete dsigmaForIntegration;
}
double DipoleModel_bSat::dsigmadb2epForIntegration(double *x, double* par)
{
double r=par[0];
double xprobe=par[1];
double b = *x;
return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe);
}
//***********bNonSat:*************************************************
DipoleModel_bNonSat::DipoleModel_bNonSat()
{
//
// Set the parameters. Note that we need bNonSat to calculate
// the skewedness correction for bSat.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(settings->dipoleModelType(), settings->dipoleModelParameterSet());
}
DipoleModel_bNonSat::~DipoleModel_bNonSat(){}
double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe)
{
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
double bDep=bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){
int A=nucleus()->A();
double TA=nucleus()->T(b)/A;
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double result=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
return result;
}
//***********bCGC:*****************************************************
DipoleModel_bCGC::DipoleModel_bCGC()
{
//
// Set the parameters. Note that we enforce here the bNonSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bCGC, settings->dipoleModelParameterSet());
}
void DipoleModel_bCGC::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
//iConfiguration is a dummy for bCGC. Do this to avoid compiler warnings:
iConfiguration++;
mNucleus.generate();
mConfigurationExists=true;
}
double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x)
{
double result=1;
for (unsigned int iA=0; iAkappa();
double N0 = mParameters->N0();
double x0 = mParameters->x0();
double lambda = mParameters->lambda();
double gammas = mParameters->gammas();
double A = -N0*N0*gammas*gammas/((1-N0)*(1-N0)*log(1-N0));
double B = 0.5*pow(1-N0,-(1-N0)/(N0*gammas));
double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b));
double rQs = r*Qs/hbarc;
double result=0;
if (rQs <= 2)
result = 2*N0*pow(0.5*rQs, 2*(gammas+(1/(kappa*lambda*Y))*log(2/rQs)));
else
result = 2*(1 - exp(-A*log(B*rQs)*log(B*rQs)));
return result;
}
double DipoleModel_bCGC::bDependence(double b)
{
double gammas = mParameters->gammas();
double Bcgc = mParameters->Bcgc();
return exp(-0.5*b*b/Bcgc/gammas/hbarc2);
}
Index: trunk/src/FinalStateGenerator.cpp
===================================================================
--- trunk/src/FinalStateGenerator.cpp (revision 377)
+++ trunk/src/FinalStateGenerator.cpp (revision 378)
@@ -1,50 +1,50 @@
//==============================================================================
// FinalStateGenerator.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 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 "FinalStateGenerator.h"
#include "Event.h"
#include
using namespace std;
FinalStateGenerator::FinalStateGenerator()
{
mT = 0;
mQ2 = 0;
mY = 0;
mS = 0;
mMY2 = 0;
mMassVM = 0;
mA = 0;
mIsIncoherent = false;
}
FinalStateGenerator::~FinalStateGenerator() {/* no op */}
bool FinalStateGenerator::isValid(TLorentzVector & v) const
{
for (int i=0; i<4; i++) {
if (std::isnan(v[0])) return false;
if (std::isinf(v[0])) return false;
}
return true;
}