Index: trunk/src/DipoleModelParameters.h
===================================================================
--- trunk/src/DipoleModelParameters.h (revision 392)
+++ trunk/src/DipoleModelParameters.h (revision 393)
@@ -1,137 +1,140 @@
//==============================================================================
// DipoleModelParameters.h
//
// 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;
+ double rMax() const; // BSTT only
// 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;
+ double mRMax; // BSTT only
// 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::rMax() const {return mRMax;}
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/Table.cpp
===================================================================
--- trunk/src/Table.cpp (revision 392)
+++ trunk/src/Table.cpp (revision 393)
@@ -1,1453 +1,1455 @@
//==============================================================================
// Table.cpp
//
// 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$
//==============================================================================
//
// Table data is stored in a TH3F.
//
// All information is stored in the histogram title.
//
// Table ID is encoded as follows:
//
// histogram title is a number that is to be interpreted
// as an uint64_t with the bits set as follows:
//
// bit 0: content type: 0 for , 1 for (see also bits 33 and 45)
// bit 1: polarization: T=1, L=0
// bit 2: t encoding: 0 for |t|, 1 for log(|t|)
// bit 3: W2 encoding: 0 for lin, 1 for log
// bit 4: Q2 encoding: 0 for lin, 1 for log
// bit 5-7: dipole model type (Enumerations.h)
// bit 8-15: mass number A
// bit 16-31: vector meson ID (PDG)
// bit 32: content encoding: 0 in lin, 1 in log
// bit 33: content type is lambda_ (bits 0 and 45 are zero in this case)
// bit 34-41: priority
// bit 42-44: dipole model parameter set (Enumerations.h)
// bit 45: content type is variance_A (bit 33 and 45 are zero in this case)
// bit 46: upc only table (tables with one Q2 bin only)
// bit 47: xpomeron encoding: 0 for lin, 1 for log (UPC only)
// bit 48: content type is lambda_skewedness (bits 0, 33, and 45 are zero in this case).
// bit 49-63: not used
//
// log implies ln
// bit 0 is the LSB
//
// Internal histogram: x = Q2, y = W2, z = t (or the logs)
// x = xpomeron, y = t (or the logs) // UPC only
//
// Actual value is taken at the center of the bin.
// min/max functions (e.g. minQ2()) return the value of the first/last
// bin center.
//==============================================================================
#include "Table.h"
#include "TFile.h"
#include "TH2F.h"
#include "TH3F.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "GridSpline.h"
#define PR(x) cout << #x << " = " << (x) << endl;
Table::Table()
{
mID = 0;
m3DHisto = 0;
m2DHisto = 0;
mFillCounter = 0;
mBackupFrequence = 0;
}
Table& Table::operator=(const Table& tab)
{
if (this != &tab) {
delete m3DHisto;
delete m2DHisto;
mID = tab.mID;
if (tab.m3DHisto) {
m3DHisto = new TH3F(*(tab.m3DHisto));
m3DHisto->SetDirectory(0);
}
if (tab.m2DHisto) {
m2DHisto = new TH2F(*(tab.m2DHisto));
m2DHisto->SetDirectory(0);
}
mFilename = tab.mID;
mFillCounter = tab.mID;
mBackupFrequence = tab.mID;
mBackupPrefix = tab.mID;
mLastBackupFilename = tab.mID;
}
return *this;
}
Table::Table(const Table& tab)
{
mID = tab.mID;
if (tab.m3DHisto) {
m3DHisto = new TH3F(*(tab.m3DHisto));
m3DHisto->SetDirectory(0);
}
if (tab.m2DHisto) {
m2DHisto = new TH2F(*(tab.m2DHisto));
m2DHisto->SetDirectory(0);
}
mFilename = tab.mID;
mFillCounter = tab.mID;
mBackupFrequence = tab.mID;
mBackupPrefix = tab.mID;
mLastBackupFilename = tab.mID;
}
Table::~Table()
{
delete m2DHisto;
delete m3DHisto;
}
unsigned int Table::create(int nbinsQ2, double Q2min, double Q2max,
int nbinsW2, double W2min, double W2max,
int nbinsT, double tmin, double tmax,
bool logQ2, bool logW2, bool logt, bool logC,
AmplitudeMoment mom, GammaPolarization pol,
unsigned int A, int vm,
DipoleModelType model, DipoleModelParameterSet pset,
const char* filename, unsigned char priority)
{
if (m3DHisto != 0) {
cout << "Table:create(): cannot create, table already exists." << endl;
return 0;
}
if (nbinsQ2 < 2) {
cout << "Table:create(): number of bins in Q2 must be at least 2." << endl;
return 0;
}
if (nbinsW2 < 2) {
cout << "Table:create(): number of bins in W2 must be at least 2." << endl;
return 0;
}
if (nbinsT < 2) {
cout << "Table:create(): number of bins in t must be at least 2." << endl;
return 0;
}
if (filename) mFilename = string(filename); // name of file where table is written to
//
// Create table ID first
//
mID = (vm << 16);
mID |= (A << 8);
mID |= (model << 5);
mID |= (static_cast(pset) << 42);
if (mom == mean_A2) mID |= 1;
if (pol == transverse) mID |= (1 << 1);
if (logt) mID |= (1 << 2);
if (logW2) mID |= (1 << 3);
if (logQ2) mID |= (1 << 4);
if (logC) mID |= (static_cast(1) << 32);
if (mom == lambda_real) mID |= (static_cast(1) << 33);
if (mom == lambda_skew) mID |= (static_cast(1) << 48);
mID |= (static_cast(priority) << 34);
if (mom == variance_A) mID |= (static_cast(1) << 45);
ostringstream titlestream;
titlestream << mID;
string title = titlestream.str();
//
// Binning
//
// Interpolate() only works up to the bin center
// of the first and last bin. To guarantee that
// the full range is accessible we shift the min
// and max of each axis.
//
tmin = fabs(tmin);
tmax = fabs(tmax);
if (tmin > tmax) swap(tmin, tmax);
if (logQ2) Q2min = log(Q2min);
if (logQ2) Q2max = log(Q2max);
if (logW2) W2min = log(W2min);
if (logW2) W2max = log(W2max);
if (logt) tmin = log(tmin);
if (logt) tmax = log(tmax);
double binWidthQ2 = (Q2max - Q2min)/(nbinsQ2-1);
double binWidthW2 = (W2max - W2min)/(nbinsW2-1);
double binWidthT = (tmax - tmin)/(nbinsT-1);
//
// Book 3D histo to hold table
//
m3DHisto = new TH3F("table", title.c_str(),
nbinsQ2, Q2min-binWidthQ2/2., Q2max+binWidthQ2/2.,
nbinsW2, W2min-binWidthW2/2., W2max+binWidthW2/2.,
nbinsT, tmin-binWidthT/2., tmax+binWidthT/2.);
m3DHisto->SetDirectory(0);
return nbinsQ2*nbinsW2*nbinsT; // return total number of bins
}
//
// Overloaded function for UPC only
//
unsigned int Table::create(int nbinsX, double Xmin, double Xmax,
int nbinsT, double tmin, double tmax,
bool logX, bool logt, bool logC,
AmplitudeMoment mom, GammaPolarization pol,
unsigned int A, int vm,
DipoleModelType model, DipoleModelParameterSet pset,
const char* filename, unsigned char priority)
{
if (m2DHisto != 0) {
cout << "Table:create(): cannot create, table already exists." << endl;
return 0;
}
if (nbinsX < 2) {
cout << "Table:create(): number of bins in x must be at least 2." << endl;
return 0;
}
if (nbinsT < 2) {
cout << "Table:create(): number of bins in t must be at least 2." << endl;
return 0;
}
if (filename) mFilename = string(filename); // name of file where table is written to
//
// Create table ID first
//
mID = (vm << 16);
mID |= (A << 8);
mID |= (model << 5);
mID |= (static_cast(pset) << 42);
if (mom == mean_A2) mID |= 1;
if (pol == transverse) mID |= (1 << 1);
if (logt) mID |= (1 << 2);
if (logC) mID |= (static_cast(1) << 32);
if (mom == lambda_real) mID |= (static_cast(1) << 33);
if (mom == lambda_skew) mID |= (static_cast(1) << 48);
mID |= (static_cast(priority) << 34);
if (mom == variance_A) mID |= (static_cast(1) << 45);
mID |= (static_cast(1) << 46); // flag UPC
if (logX) mID |= (static_cast(1) << 47);
ostringstream titlestream;
titlestream << mID;
string title = titlestream.str();
//
// Binning
//
// Interpolate() only works up to the bin center
// of the first and last bin. To guarantee that
// the full range is accessible we shift the min
// and max of each axis.
//
tmin = fabs(tmin);
tmax = fabs(tmax);
if (tmin > tmax) swap(tmin, tmax);
if (logX) Xmin = log(Xmin);
if (logX) Xmax = log(Xmax);
if (logt) tmin = log(tmin);
if (logt) tmax = log(tmax);
double binWidthX = (Xmax - Xmin)/(nbinsX-1);
double binWidthT = (tmax - tmin)/(nbinsT-1);
//
// Book 2D histo to hold table
//
m2DHisto = new TH2F("table", title.c_str(),
nbinsX, Xmin-binWidthX/2., Xmax+binWidthX/2.,
nbinsT, tmin-binWidthT/2., tmax+binWidthT/2.);
m2DHisto->SetDirectory(0);
return nbinsX*nbinsT; // return total number of bins
}
int Table::numberOfEntries() const
{
if (m3DHisto) {
int nx = m3DHisto->GetXaxis()->GetNbins();
int ny = m3DHisto->GetYaxis()->GetNbins();
int nz = m3DHisto->GetZaxis()->GetNbins();
return nx*ny*nz;
}
else if (m2DHisto) { // UPC case
int nx = m2DHisto->GetXaxis()->GetNbins();
int ny = m2DHisto->GetYaxis()->GetNbins();
return nx*ny;
}
else
return 0;
}
void Table::binXYZ(int globalBin, int& binx, int& biny, int& binz) const
{
//
// Find correct bins for each axis.
// Don't use ROOT global bins here,
// they are a mess since they cross
// over underflow and overflow bins.
//
// All bins returned count starting
// at 1. The global bin starts at
// 0 as usual in C++.
//
int nx = m3DHisto->GetXaxis()->GetNbins();
int ny = m3DHisto->GetYaxis()->GetNbins();
binx = globalBin%nx;
biny = ((globalBin-binx)/nx)%ny;
binz = ((globalBin-binx)/nx -biny)/ny;
binx++; biny++; binz++;
}
//
// UPC version of binXYZ()
//
void Table::binXY(int globalBin, int& binx, int& biny) const
{
//
// Find correct bins for each axis.
// Don't use ROOT global bins here,
// they are a mess since they cross
// over underflow and overflow bins.
//
// All bins returned count starting
// at 1. The global bin starts at
// 0 as usual in C++.
//
int nx = m2DHisto->GetXaxis()->GetNbins();
binx = globalBin%nx;
biny = (globalBin-binx)/nx;
binx++; biny++;
}
void Table::binCenter(int i, double& Q2, double& W2, double& t) const
{
int binx, biny, binz;
binXYZ(i, binx, biny, binz);
// cout << i << '\t' << binx << '\t' << biny << '\t' << binz << endl;
double x = m3DHisto->GetXaxis()->GetBinCenter(binx);
double y = m3DHisto->GetYaxis()->GetBinCenter(biny);
double z = m3DHisto->GetZaxis()->GetBinCenter(binz);
Q2 = isLogQ2() ? exp(x) : x;
W2 = isLogW2() ? exp(y) : y;
t = isLogT() ? -exp(z) : -z;
if (t > 0) t = 0; // avoid numeric glitch
}
//
// UPC overload
//
void Table::binCenter(int i, double& xpom, double& t) const
{
int binx, biny;
binXY(i, binx, biny);
// cout << i << '\t' << binx << '\t' << biny << endl;
double x = m2DHisto->GetXaxis()->GetBinCenter(binx);
double y = m2DHisto->GetYaxis()->GetBinCenter(biny);
xpom = isLogX() ? exp(x) : x;
t = isLogT() ? -exp(y) : -y;
if (t > 0) t = 0; // avoid numeric glitch
}
void Table::fill(int i, double val, double err)
{
int binx, biny, binz;
if (isUPC())
binXY(i, binx, biny);
else
binXYZ(i, binx, biny, binz);
if (isLogContent()) {
if (val == 0) {
// cout << "Table::fill(): warning, log scale requested but value = 0." << endl;
val = numeric_limits::min();
}
val = log(fabs(val));
}
if (isUPC())
m2DHisto->SetBinContent(binx, biny, val);
else
m3DHisto->SetBinContent(binx, biny, binz, val);
if(err!=0.) {
if (isUPC())
m2DHisto->SetBinError(binx, biny, err);
else
m3DHisto->SetBinError(binx, biny, binz, err);
}
mFillCounter++;
//
// Check if backup is due
//
if (mBackupFrequence)
if (mFillCounter%mBackupFrequence == 0) backup(i);
}
bool Table::fexists(const char* filename) const
{
ifstream ifs(filename);
return !ifs.fail();
}
void Table::write(const char* filename)
{
//
// 'filename' is optional argument. Null value implies
// that one was already passed via create(). Check
// that this is the case. If one is given, it is used
// no matter if it already was defined or not.
//
if (filename)
mFilename = string(filename);
else {
if (mFilename.empty()) { // should be defined but is not
cout << "Table::write(): Warning, no filename supplied. Will use 'sartre_table.root'." << endl;
mFilename = string("sartre_table.root");
}
}
//
// Check if file exist. If so alter
// the filename. Creating tables is
// a time consuming business so we do
// not want to lose precious data if
// something goes wrong here and we do
// want to prevent accidents.
//
while (fexists(mFilename.c_str())) {
mFilename += string(".save");
}
TFile hfile(mFilename.c_str(),"RECREATE");
if (isUPC())
m2DHisto->Write();
else
m3DHisto->Write();
hfile.Close();
cout << "Table::write(): table stored in file '" << mFilename.c_str()<< "'." << endl;
if (mBackupFrequence && !mLastBackupFilename.empty())
remove(mLastBackupFilename.c_str());
}
bool Table::read(const char* filename)
{
//
// Read histogram into memory
//
TFile *file = TFile::Open(filename,"READ");
if (file == 0) {
cout << "Table::read(): failed opening file '" << filename << "'." << endl;
return false;
}
//
// Clean up
//
delete m3DHisto; m3DHisto = 0;
delete m2DHisto; m2DHisto = 0;
//
// Read plain object first and check title which
// tells us if it is an UPC table or a std one.
// Then cast into the proper object.
//
auto ptr = file->Get("table");
if (ptr == 0) {
cout << "Table::read(): failed retrieving table from file '" << filename << "'." << endl;
return false;
}
mID = atoll(ptr->GetTitle()); // for isUPC() to work
if (isUPC()) {
m2DHisto = reinterpret_cast(ptr);
m2DHisto->SetDirectory(0);
}
else {
m3DHisto = reinterpret_cast(ptr);
m3DHisto->SetDirectory(0);
}
if (m2DHisto == 0 && m3DHisto == 0) {
cout << "Table::read(): failed retrieving table from file '" << filename << "'." << endl;
return false;
}
file->Close();
mFilename = string(filename);
return true;
}
int Table::vectorMesonId() const {return ((mID >> 16) & 0xFFFF);}
unsigned int Table::dipoleModelType() const {return ((mID >> 5) & 0x7);}
unsigned int Table::dipoleModelParameterSet() const {return ((mID >> 42) & 0x7);}
unsigned int Table::A() const {return ((mID >> 8) & 0xFF);}
bool Table::isLongitudinal() const {return !isTransverse();}
bool Table::isTransverse() const {return (mID & (1 << 1));}
GammaPolarization Table::polarization() const {return (isTransverse() ? transverse : longitudinal);}
AmplitudeMoment Table::moment() const {
if (isLambdaA())
return lambda_real;
else if (isMeanA())
return mean_A;
else if (isVarianceA())
return variance_A;
else if (isLambdaSkew())
return lambda_skew;
else
return mean_A2;
}
bool Table::isMeanA() const {
// bits 0, 33, 45 and 48 are not set
return !(mID & 1) && !(mID & (static_cast(1) << 33)) && !(mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isMeanA2() const {
// bit 0 is set, bits 33, 45 and 48 are not set
return (mID & 1) && !(mID & (static_cast(1) << 33)) && !(mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isLambdaA() const {
// bit 33 is set, bits 0 and 45 are not set
return !(mID & 1) && (mID & (static_cast(1) << 33)) && !(mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isLambdaSkew() const {
// bit 48 is set, bits 0, 33, and 45 are not set
return !(mID & 1) && (mID & (static_cast(1) << 48)) && !(mID & (static_cast(1) << 33))
&& !(mID & (static_cast(1) << 45));
}
bool Table::isVarianceA() const {
// bit 45 is set, bitss 0, 33, and 48 are not set
return !(mID & 1) && !(mID & (static_cast(1) << 33)) && (mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isLogQ2() const {return (mID & (1 << 4));}
bool Table::isLogW2() const {return (mID & (1 << 3));}
bool Table::isLogT() const {return (mID & (1 << 2));}
bool Table::isLogX() const {return (mID & (static_cast(1) << 47));}
bool Table::isLogContent() const {return (mID & (static_cast(1) << 32));}
unsigned int Table::priority() const {return ((mID >> 34) & 0xFF);}
bool Table::isUPC() const {return (mID & (static_cast(1) << 46));}
uint64_t Table::id() const {return mID;}
double Table::binWidthQ2() const
{
if (!m3DHisto) return 0;
return m3DHisto->GetXaxis()->GetBinWidth(1);
}
double Table::binWidthW2() const
{
if (!m3DHisto) return 0;
return m3DHisto->GetYaxis()->GetBinWidth(1);
}
double Table::binWidthT() const
{
if (isUPC()) {
return m2DHisto->GetYaxis()->GetBinWidth(1);
}
else
return m3DHisto->GetZaxis()->GetBinWidth(1);
}
double Table::binWidthX() const
{
if (!m2DHisto) return 0;
return m2DHisto->GetXaxis()->GetBinWidth(1);
}
double Table::get(double Q2, double W2, double t) const
{
if (m3DHisto == 0) return 0;
//
// Transform variables to how they are stored in the table
//
double x = isLogQ2() ? log(Q2) : Q2;
double y = isLogW2() ? log(W2) : W2;
t = fabs(t);
double z = isLogT() ? log(t) : t;
//
// Tiny numerical glitches will cause TH3F::Interpolate() to fail
// since it requires that the variables lie between the first
// and last bin center, excluding the centers.
// Here we enforce that this is the case. The downside of this
// is that *all* values of Q2, W2, t will be forced to lie within.
// Hence the user has to make sure that the values are within
// the boundaries (see minQ2(), maxQ2(), minW2() etc.) before
// calling get().
// In this case the corrections below are tiny and are only
// applied to avoid minor glitches that do not affect the results.
//
TAxis *axis = m3DHisto->GetXaxis();
x = max(x, axis->GetBinCenter(1)+numeric_limits::epsilon());
x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
axis = m3DHisto->GetYaxis();
y = max(y, axis->GetBinCenter(1)+numeric_limits::epsilon());
y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
axis = m3DHisto->GetZaxis();
z = max(z, axis->GetBinCenter(1)+numeric_limits::epsilon());
z = min(z, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
// double result = InterpolateGridSpline(x, y, z); // tmp uncommented until 0's in tables are cleared
// double result = m3DHisto->Interpolate(x, y, z);
double result = modInterpolation(x, y, z);
if (result == 0 && isLogContent()) {
cout << "Table::get(): warning, 0 is not a valid table content when working in log scale." << endl;
}
if (isLogContent()) result = exp(result);
return result;
}
//
// UPC version
//
double Table::get(double xpom, double t) const
{
if (m2DHisto == 0) return 0;
//
// Transform variables to how they are stored in the table
//
double x = isLogX() ? log(xpom) : xpom;
t = fabs(t);
double y = isLogT() ? log(t) : t;
//
// Tiny numerical glitches will cause TH2F::Interpolate() to fail
// since it requires that the variables lie between the first
// and last bin center, excluding the centers.
// Here we enforce that this is the case. The downside of this
// is that *all* values of x, t will be forced to lie within.
// Hence the user has to make sure that the values are within
// the boundaries (see minX(), maxX(), minT() etc.) before
// calling get().
// In this case the corrections below are tiny and are only
// applied to avoid minor glitches that do not affect the results.
//
TAxis *axis = m2DHisto->GetXaxis();
x = max(x, axis->GetBinCenter(1)+numeric_limits::epsilon());
x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
axis = m2DHisto->GetYaxis();
y = max(y, axis->GetBinCenter(1)+numeric_limits::epsilon());
y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
double result = modInterpolation(x, y); // xpom, t
if (result == 0 && isLogContent()) {
cout << "Table::get(): warning, 0 is not a valid table content when working in log scale." << endl;
}
if (isLogContent()) result = exp(result);
return result;
}
double Table::minQ2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetXaxis();
double val = axis->GetBinCenter(1);
return isLogQ2() ? exp(val) : val;
}
double Table::maxQ2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetXaxis();
double val = axis->GetBinCenter(axis->GetNbins());
return isLogQ2() ? exp(val) : val;
}
double Table::minW2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetYaxis();
double val = axis->GetBinCenter(1);
return isLogW2() ? exp(val) : val;
}
double Table::maxW2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetYaxis();
double val = axis->GetBinCenter(axis->GetNbins());
return isLogW2() ? exp(val) : val;
}
double Table::minT() const
{
if (m2DHisto == 0 && m3DHisto == 0) return 0;
TAxis *axis = 0;
if (isUPC())
axis = m2DHisto->GetYaxis();
else
axis = m3DHisto->GetZaxis();
double val = axis->GetBinCenter(axis->GetNbins()); // t always as |t| in table
return isLogT() ? -exp(val) : -val;
}
double Table::maxT() const
{
if (m2DHisto == 0 && m3DHisto == 0) return 0;
TAxis *axis = 0;
if (isUPC())
axis = m2DHisto->GetYaxis();
else
axis = m3DHisto->GetZaxis();
double val = axis->GetBinCenter(1); // t always as |t| in table
double result = isLogT() ? -exp(val) : -val;
if (result > 0) {
// cout << "Table::maxT(): warning, t > 0, t (" << result << ") set to 0 now." << endl;
result = 0;
}
return result;
}
double Table::minX() const
{
if (m2DHisto == 0) return 0;
TAxis *axis = m2DHisto->GetXaxis();
double val = axis->GetBinCenter(1);
return isLogX() ? exp(val) : val;
}
double Table::maxX() const
{
if (m2DHisto == 0) return 0;
TAxis *axis = m2DHisto->GetXaxis();
double val = axis->GetBinCenter(axis->GetNbins());
return isLogX() ? exp(val) : val;
}
string Table::filename() const {return mFilename;}
bool Table::writeToBinaryFile(const string& filename, bool verbose) const
{
//
// Open binary file
//
ofstream binaryFile;
binaryFile.open (filename, ios::out | ios::binary);
if (binaryFile.fail()) {
cout << "Table::writeToBinaryFile: failed to open file '" << filename << "'." << endl;
return false;
}
//
// Loop over all entries
//
double Q2, W2, xpom, t;
int binx, biny, binz;
double rawContent;
for (int i=0; iGetBinContent(binx, biny);
if (isLogContent()) {
rawContent = exp(rawContent);
}
binaryFile.write(reinterpret_cast(&i), sizeof(i));
binaryFile.write(reinterpret_cast(&xpom), sizeof(xpom));
binaryFile.write(reinterpret_cast(&t), sizeof(t));
binaryFile.write(reinterpret_cast(&rawContent), sizeof(rawContent));
if (verbose)
cout << "i=" << i << "\txpom=" << xpom << "\tt=" << t << "\tvalue=" << rawContent << endl;
}
else {
binXYZ(i, binx, biny, binz);
binCenter(i, Q2, W2, t);
rawContent = m3DHisto->GetBinContent(binx, biny, binz);
if (isLogContent()) {
rawContent = exp(rawContent);
}
binaryFile.write(reinterpret_cast(&i), sizeof(i));
binaryFile.write(reinterpret_cast(&Q2), sizeof(Q2));
binaryFile.write(reinterpret_cast(&W2), sizeof(W2));
binaryFile.write(reinterpret_cast(&t), sizeof(t));
binaryFile.write(reinterpret_cast(&rawContent), sizeof(rawContent));
if (verbose)
cout << "i=" << i << "\tQ2=" << Q2 << "\tW2=" << W2 << "\tt=" << t << "\tvalue=" << rawContent << endl;
}
}
binaryFile.close();
return true;
}
void Table::list(ostream& os, bool printContent, bool printStatistics, int startBin, int endBin) const
{
//
// List table info as compact as possible
//
ios::fmtflags fmt = os.flags(); // store current I/O flags
if ((isUPC() && !m2DHisto) || (!isUPC() && !m3DHisto)) {
os << "Table::list(): table is undefined." << endl;
return;
}
//
// First row: table id, A, and content
//
os << setw(11) << "Table ID = " << setw(16) << left << mID;
os << setw(7) << right << "A = " << setw(4) << left << A();
os << setw(20) << right << "content = " << (isLogContent() ? "log(" : "");
if (isMeanA())
os << "";
else if (isMeanA2())
os << "";
else if (isLambdaA())
os << "lambda_";
else if (isVarianceA())
os << "variance_A";
else if (isLambdaSkew())
os << "lambda_skew";
else
os << "unknown";
os << (isLogContent() ? ")" : "") << endl;
//
// Second row: vm and polarization
//
os << setw(34) << right << "vmId = " << setw(4) << left << vectorMesonId();
os << setw(20) << right << "polarization = " << (isTransverse() ? 'T' : 'L') << endl;
//
// Third row: dipole model and parameter set
//
os << setw(34) << right << "model = " << setw(7) << left;
if (dipoleModelType() == bSat)
os << "bSat ";
else if (dipoleModelType() == bNonSat)
os << "bNonSat";
else
os << "bCGC ";
os << setw(17) << right << "parameter set = ";
if (dipoleModelParameterSet() == KMW)
os << "KMW ";
else if (dipoleModelParameterSet() == HMPZ)
os << "HMPZ";
+ else if (dipoleModelParameterSet() == BSTT)
+ os << "BSTT";
else if (dipoleModelParameterSet() == CUSTOM)
os << "CUSTOM";
else
os << "? ";
os << endl;
//
// Third row: priority and UPC
//
os << setw(34) << right << "priority = " << setw(4) << left << priority();
os << setw(20) << right << "UPC = " << (isUPC() ? "yes" : "no") << endl;
//
// Next three rows: Q2, W2, t bins, range and log/linear
// or for UPC: x, t
//
if (isUPC()) {
os << setw(35) << right << "xp = [" << minX() << ", " << maxX()
<< "; bins = " << m2DHisto->GetXaxis()->GetNbins() << (isLogX() ? "; log]" : "; lin]") << endl;
os << setw(35) << right << "t = [" << minT() << ", " << maxT()
<< "; bins = " << m2DHisto->GetYaxis()->GetNbins() << (isLogT() ? "; log]" : "; lin]") << endl;
}
else {
os << setw(35) << right << "Q2 = [" << minQ2() << ", " << maxQ2()
<< "; bins = " << m3DHisto->GetXaxis()->GetNbins() << (isLogQ2() ? "; log]" : "; lin]") << endl;
os << setw(35) << right << "W2 = [" << minW2() << ", " << maxW2()
<< "; bins = " << m3DHisto->GetYaxis()->GetNbins() << (isLogW2() ? "; log]" : "; lin]") << endl;
os << setw(35) << right << "t = [" << minT() << ", " << maxT()
<< "; bins = " << m3DHisto->GetZaxis()->GetNbins() << (isLogT() ? "; log]" : "; lin]") << endl;
}
//
// Filename at the end
//
os << setw(34) << right << "file = " << mFilename.c_str() << endl;
os << endl;
//
// Print content (optional)
//
double Q2, W2, t, xpom;
int binx, biny, binz;
if (printContent) {
int thePrecision = os.precision();
double rawContent = 0;
if(endBin==0) endBin=numberOfEntries();
for (int i=startBin; iGetBinContent(binx, biny);
double error = m2DHisto->GetBinError(binx, biny);
os << "bin = " << setw(8) << left << i;
os << "xp = " << setw(13) << left << fixed << setprecision(5) << xpom;
os << "t = " << setw(16) << left << scientific << t;
os << "value = " << setw(15) << left << scientific << value;
os << "(binx = " << setw(5) << left << binx;
os << "biny = " << setw(5) << left << biny;
os << "content = " << setw(14) << left << scientific << rawContent;
os << "error = " << left << error << ')';
}
else {
binXYZ(i, binx, biny, binz);
binCenter(i, Q2, W2, t);
double value = get(Q2, W2, t);
rawContent = m3DHisto->GetBinContent(binx, biny, binz);
double error = m3DHisto->GetBinError(binx, biny, binz);
os << "bin = " << setw(8) << left << i;
os << "Q2 = " << setw(13) << left << fixed << setprecision(5) << Q2;
os << "W2 = " << setw(12) << left << fixed << W2;
os << "t = " << setw(16) << left << scientific << t;
os << "value = " << setw(15) << left << scientific << value;
os << "(binx = " << setw(5) << left << binx;
os << "biny = " << setw(5) << left << biny;
os << "binz = " << setw(5) << left << binz;
os << "content = " << setw(14) << left << scientific << rawContent;
os << "error = " << left << error << ')';
}
if ( (isLogContent() && rawContent <= log(numeric_limits::min()*2)) || (rawContent == 0)) cout << " Z";
cout << endl;
}
os << endl;
os.precision(thePrecision);
}
//
// Print statistics (optional)
//
if (printStatistics) {
int nEmpty = 0;
int nInvalid = 0;
int nNegative = 0;
double sum = 0;
double maxContent = 0;
double minContent = numeric_limits::max();
int minBin, maxBin;
minBin = maxBin = 0;
double c = 0;
for (int i=0; iGetBinContent(binx, biny);
}
else {
binXYZ(i, binx, biny, binz);
c = m3DHisto->GetBinContent(binx, biny, binz);
}
if (c == 0) nEmpty++;
if ( !isLogContent() && c < 0) nNegative++;
if (std::isnan(c) || std::isinf(c) || !std::isfinite(c)) nInvalid++;
if (isLogContent()) c = exp(c);
if (c > maxContent) {maxContent = c; maxBin = i;}
if (c >= 0 && c < minContent) {minContent = c; minBin = i;}
if (c > 0) sum += c;
}
os << setw(34) << right << "total number of cells = " << numberOfEntries() << endl;
os << setw(34) << right << "cells with negative content = " << nNegative << endl;
os << setw(34) << right << "cells with no (0) content = " << nEmpty << endl;
os << setw(34) << right << "cells with invalid content = " << nInvalid << endl;
if (isUPC()) {
binCenter(minBin, xpom, t);
os << setw(34) << right << "minimum content = " << minContent << " at xp = " << xpom << ", t = " << t << endl;
binCenter(maxBin, xpom, t);
os << setw(34) << right << "maximum content = " << maxContent << " at xp = " << xpom << ", t = " << t << endl;
}
else {
binCenter(minBin, Q2, W2, t);
os << setw(34) << right << "minimum content = " << minContent << " at Q2 = " << Q2 << ", W2 = " << W2 << ", t = " << t << endl;
binCenter(maxBin, Q2, W2, t);
os << setw(34) << right << "maximum content = " << maxContent << " at Q2 = " << Q2 << ", W2 = " << W2 << ", t = " << t << endl;
}
os << setw(34) << right << "sum = " << sum << endl;
os << endl;
}
os.flags(fmt); // restore I/O flags
}
double Table::InterpolateGridSpline(double x, double y, double z) const // Q2, W2, t
{
//
// The algorithm was taken from Cristian Constantin Lalescu.
// See http://arxiv.org/abs/0905.3564 for details.
//
// Grid points: 4, 6, or 8
// Spline order: 3,5 for order 4
// 3,5,7,9 for order 6
// 3,5,7,9,11,13 for order 8
//
//
// Find cell that refer to the position
//
int nx = m3DHisto->GetXaxis()->FindBin(x);
int ny = m3DHisto->GetYaxis()->FindBin(y);
int nz = m3DHisto->GetZaxis()->FindBin(z);
//
// Define the # of grid points depending on how far we are away
// from the edge. If at the edge we fall back to linear interpolation
// as provided by TH3.
// There must be a smarter way doing this than the code below
//
int gridPoints = 0;
if (nx-3 >= 1 && nx+4 <= m3DHisto->GetXaxis()->GetNbins() &&
ny-3 >= 1 && ny+4 <= m3DHisto->GetYaxis()->GetNbins() &&
nz-3 >= 1 && nz+4 <= m3DHisto->GetZaxis()->GetNbins()) {
gridPoints = 8;
}
else if (nx-2 >= 1 && nx+3 <= m3DHisto->GetXaxis()->GetNbins() &&
ny-2 >= 1 && ny+3 <= m3DHisto->GetYaxis()->GetNbins() &&
nz-2 >= 1 && nz+3 <= m3DHisto->GetZaxis()->GetNbins()) {
gridPoints = 6;
}
else if (nx-1 >= 1 && nx+2 <= m3DHisto->GetXaxis()->GetNbins() &&
ny-1 >= 1 && ny+2 <= m3DHisto->GetYaxis()->GetNbins() &&
nz-1 >= 1 && nz+2 <= m3DHisto->GetZaxis()->GetNbins()) {
gridPoints = 4;
}
else {
return m3DHisto->Interpolate(x,y,z);
}
//
// Find bin centers
//
double xc = m3DHisto->GetXaxis()->GetBinCenter(nx);
double yc = m3DHisto->GetYaxis()->GetBinCenter(ny);
double zc = m3DHisto->GetZaxis()->GetBinCenter(nz);
//
// Define the scale for coordinate transformations
// grid_spline expects x,y,z in bin units
//
double xscale = m3DHisto->GetXaxis()->GetBinWidth(1);
double yscale = m3DHisto->GetYaxis()->GetBinWidth(1);
double zscale = m3DHisto->GetZaxis()->GetBinWidth(1);
//
// Prepare grid spline alogorithm
//
double result, splineOrder;
if (gridPoints == 8) {
local_scal_3D<8> lf;
for (int i=-3;i<=4;i++) {
for(int j=-3;j<=4;j++) {
lf(-3,i,j) = m3DHisto->GetBinContent(nx-3, ny+i, nz+j);
lf(-2,i,j) = m3DHisto->GetBinContent(nx-2, ny+i, nz+j);
lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
lf( 3,i,j) = m3DHisto->GetBinContent(nx+3, ny+i, nz+j);
lf( 4,i,j) = m3DHisto->GetBinContent(nx+4, ny+i, nz+j);
}
}
splineOrder = 7;
result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
}
else if (gridPoints == 6) {
local_scal_3D<6> lf;
for (int i=-2;i<=3;i++) {
for(int j=-2;j<=3;j++) {
lf(-2,i,j) = m3DHisto->GetBinContent(nx-2, ny+i, nz+j);
lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
lf( 3,i,j) = m3DHisto->GetBinContent(nx+3, ny+i, nz+j);
}
}
splineOrder = 5;
result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
}
else if (gridPoints == 4) {
local_scal_3D<4> lf;
for (int i=-1;i<=2;i++) {
for(int j=-1;j<=2;j++) {
lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
}
}
splineOrder = 3;
result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
}
else {
cout << "Table::InterpolateGridSpline(): Error, illegal number of grid points." << endl;
result = 0;
}
return result;
}
double Table::modInterpolation(double x, double y, double z) const{ // Q2, W2, t
//
// After testing many points: The best results is had for content in log and rest in lin.
// However, if any of the bins that participate in the interpolation are zero
// the interpolation has to be linear.
//
// Make sure the point is within the histogram:
if(m3DHisto->GetXaxis()->GetBinCenter(1) > x ||
m3DHisto->GetXaxis()->GetBinCenter(m3DHisto->GetXaxis()->GetNbins()) < x ||
m3DHisto->GetYaxis()->GetBinCenter(1) > y ||
m3DHisto->GetYaxis()->GetBinCenter(m3DHisto->GetYaxis()->GetNbins()) < y ||
m3DHisto->GetZaxis()->GetBinCenter(1) > z ||
m3DHisto->GetZaxis()->GetBinCenter(m3DHisto->GetZaxis()->GetNbins()) < z ){
cout<<"Table::myInterpolation Error: point lies outside the limits of the table!"<GetXaxis()->FindBin(x);
if( x < m3DHisto->GetXaxis()->GetBinCenter(lbx) ) lbx--;
int ubx=lbx+1;
int lby = m3DHisto->GetYaxis()->FindBin(y);
if( y < m3DHisto->GetYaxis()->GetBinCenter(lby) ) lby--;
int uby=lby+1;
int lbz = m3DHisto->GetZaxis()->FindBin(z);
if( z < m3DHisto->GetZaxis()->GetBinCenter(lbz) ) lbz--;
int ubz=lbz+1;
//The corresponding bin-centers:
double ux=m3DHisto->GetXaxis()->GetBinCenter(ubx);
double lx=m3DHisto->GetXaxis()->GetBinCenter(lbx);
double uy=m3DHisto->GetYaxis()->GetBinCenter(uby);
double ly=m3DHisto->GetYaxis()->GetBinCenter(lby);
double uz=m3DHisto->GetZaxis()->GetBinCenter(ubz);
double lz=m3DHisto->GetZaxis()->GetBinCenter(lbz);
ux = isLogQ2() ? exp(ux) : ux;
lx = isLogQ2() ? exp(lx) : lx;
uy = isLogW2() ? exp(uy) : uy;
ly = isLogW2() ? exp(ly) : ly;
uz = isLogT() ? exp(uz) : uz;
lz = isLogT() ? exp(lz) : lz;
x = isLogQ2() ? exp(x) : x;
y = isLogW2() ? exp(y) : y;
z = isLogT() ? exp(z) : z;
//Find the distance to the point:
double xd = (x-lx)/(ux-lx);
double yd = (y-ly)/(uy-ly);
double zd = (z-lz)/(uz-lz);
// Make a vector containing all the values of the 8 surrounding bins
double v[8] = { m3DHisto->GetBinContent(lbx, lby, lbz), m3DHisto->GetBinContent(lbx, lby, ubz),
m3DHisto->GetBinContent(lbx, uby, lbz), m3DHisto->GetBinContent(lbx, uby, ubz),
m3DHisto->GetBinContent(ubx, lby, lbz), m3DHisto->GetBinContent(ubx, lby, ubz),
m3DHisto->GetBinContent(ubx, uby, lbz), m3DHisto->GetBinContent(ubx, uby, ubz) };
bool logC=true;
for(int i=0; i<8; i++){
if(isLogContent() and v[i] <= log(numeric_limits::min()*2) )
logC=false;
if(!isLogContent() and v[i] == 0)
logC=false;
}
// Make interpolation in log(content) except for when at least one of the points are zero.
if(isLogContent() and !logC){
for(int i=0; i<8; i++){
if (v[i] <= log(numeric_limits::min()*2))
v[i] = 0;
else
v[i]=exp(v[i]);
}
}
if(!isLogContent() and logC){
for(int i=0; i<8; i++)
v[i]=log(v[i]);
}
// First make four 1D interpolations in the z-direction
double i1 = v[0] * (1 - zd) + v[1] * zd;
double i2 = v[2] * (1 - zd) + v[3] * zd;
double j1 = v[4] * (1 - zd) + v[5] * zd;
double j2 = v[6] * (1 - zd) + v[7] * zd;
// Secondly, make two 1D interpolations in the y-direction using the values obtained.
double w1 = i1 * (1 - yd) + i2 * yd;
double w2 = j1 * (1 - yd) + j2 * yd;
// Finaly, make 1D interpolation in the x-direction
double result= w1 * (1-xd) + w2 * xd;
//Reverse exp/log compensation to get real result back:
if(isLogContent() and !logC){
if(result == 0)
result=log(numeric_limits::min());
else
result=log(result);
}
if(!isLogContent() and logC)
result=exp(result);
return result;
}
//
// UPC version
//
double Table::modInterpolation(double x, double y) const{ // x, t
// Make sure the point is within the histogram:
if(m2DHisto->GetXaxis()->GetBinCenter(1) > x ||
m2DHisto->GetXaxis()->GetBinCenter(m2DHisto->GetXaxis()->GetNbins()) < x ||
m2DHisto->GetYaxis()->GetBinCenter(1) > y ||
m2DHisto->GetYaxis()->GetBinCenter(m2DHisto->GetYaxis()->GetNbins()) < y) {
cout<<"Table::myInterpolation Error: point lies outside the limits of the table!"<GetXaxis()->FindBin(x);
int bin_y = m2DHisto->GetYaxis()->FindBin(y);
// Find quadrant of the bin we are in
int quadrant = 0;
double dx = m2DHisto->GetXaxis()->GetBinUpEdge(bin_x)-x;
double dy = m2DHisto->GetYaxis()->GetBinUpEdge(bin_y)-y;
if (dx<=m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy<=m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 1; // upper right
if (dx>m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy<=m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 2; // upper left
if (dx>m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy>m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 3; // lower left
if (dx<=m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy>m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 4; // lower right
double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
switch(quadrant) {
case 1:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x+1);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y+1);
break;
case 2:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x-1);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y+1);
break;
case 3:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x-1);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y-1);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
break;
case 4:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y-1);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x+1);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
break;
}
// Now find the bins we interpolate with
int bin_x1 = m2DHisto->GetXaxis()->FindBin(x1);
if (bin_x1<1) bin_x1=1;
int bin_x2 = m2DHisto->GetXaxis()->FindBin(x2);
if (bin_x2>m2DHisto->GetXaxis()->GetNbins()) bin_x2=m2DHisto->GetXaxis()->GetNbins();
int bin_y1 = m2DHisto->GetYaxis()->FindBin(y1);
if (bin_y1<1) bin_y1=1;
int bin_y2 = m2DHisto->GetYaxis()->FindBin(y2);
if (bin_y2>m2DHisto->GetYaxis()->GetNbins()) bin_y2=m2DHisto->GetYaxis()->GetNbins();
// Get content
int bin_q22 = m2DHisto->GetBin(bin_x2,bin_y2);
int bin_q12 = m2DHisto->GetBin(bin_x1,bin_y2);
int bin_q11 = m2DHisto->GetBin(bin_x1,bin_y1);
int bin_q21 = m2DHisto->GetBin(bin_x2,bin_y1);
double q11 = m2DHisto->GetBinContent(bin_q11);
double q12 = m2DHisto->GetBinContent(bin_q12);
double q21 = m2DHisto->GetBinContent(bin_q21);
double q22 = m2DHisto->GetBinContent(bin_q22);
//
// As explained in the 3D version we interpolate on linear x, y
// but on log content, except when 0's are involved.
//
bool logC = true;
if ((isLogContent() && q11 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q11 <= 0)) logC = false;
if ((isLogContent() && q12 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q12 <= 0)) logC = false;
if ((isLogContent() && q21 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q21 <= 0)) logC = false;
if ((isLogContent() && q22 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q22 <= 0)) logC = false;
if(isLogContent() && !logC) { // 0's present, use linear content
if (q11 <= log(numeric_limits::min()*2))
q11 = 0;
else
q11=exp(q11);
if (q12 <= log(numeric_limits::min()*2))
q12 = 0;
else
q12=exp(q12);
if (q21 <= log(numeric_limits::min()*2))
q21 = 0;
else
q21=exp(q21);
if (q22 <= log(numeric_limits::min()*2))
q22 = 0;
else
q22=exp(q22);
}
if (!isLogContent() && logC) {
q11 = log(q11);
q12 = log(q12);
q21 = log(q21);
q22 = log(q22);
}
x = isLogX() ? exp(x) : x;
x1 = isLogX() ? exp(x1) : x1;
x2 = isLogX() ? exp(x2) : x2;
y = isLogT() ? exp(y) : y;
y1 = isLogT() ? exp(y1) : y1;
y2 = isLogT() ? exp(y2) : y2;
// Interpolation
double d = 1.0*(x2-x1)*(y2-y1);
double result = 1.0*q11/d*(x2-x)*(y2-y)+1.0*q21/d*(x-x1)*(y2-y)+1.0*q12/d*(x2-x)*(y-y1)+1.0*q22/d*(x-x1)*(y-y1);
// Reverse exp/log compensation to get real result back:
if (isLogContent() && !logC){
if (result == 0)
result=log(numeric_limits::min());
else
result=log(result);
}
if (!isLogContent() && logC)
result=exp(result);
return result;
}
void Table::setAutobackup(const char* prefix, int freq)
{
mBackupPrefix = string(prefix);
mBackupFrequence = freq;
}
void Table::backup(int backupBin)
{
ostringstream backupFilenameStream;
backupFilenameStream << mBackupPrefix.c_str() << "_backup." << static_cast(getpid())
<< '.' << mID << '.' << backupBin << ".root";
string filename = backupFilenameStream.str();
TFile file(filename.c_str(),"RECREATE");
m3DHisto->Write();
file.Close();
if (filename != mLastBackupFilename)
remove(mLastBackupFilename.c_str());
mLastBackupFilename = filename;
time_t now = time(0);
cout << "Table::backup(): autobackup performed, file = '" << mLastBackupFilename.c_str() << "', time = " << ctime(&now);
}
int Table::globalBin(int binx, int biny, int binz) const
{
int nbinx = m3DHisto->GetXaxis()->GetNbins();
int nbiny = m3DHisto->GetYaxis()->GetNbins();
return (binz-1)*nbinx*nbiny+(biny-1)*nbinx+(binx-1);
}
//
// UPC version
//
int Table::globalBin(int binx, int biny) const
{
int nbinx = m3DHisto->GetXaxis()->GetNbins();
return (biny-1)*nbinx+(binx-1);
}
Index: trunk/src/Settings.cpp
===================================================================
--- trunk/src/Settings.cpp (revision 392)
+++ trunk/src/Settings.cpp (revision 393)
@@ -1,580 +1,584 @@
//==============================================================================
// Settings.cpp
//
// 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 "Settings.h"
#include
#include
#include
#include
#include
#include "TParticlePDG.h"
#include
#define PR(x) cout << #x << " = " << (x) << endl;
TRandom3 Settings::mRandomGenerator;
Settings::Settings()
{
mPDG = new TDatabasePDG;
mPDG->ReadPDGTable();
//
// Setup name table (map) for nuclei
//
mPeriodicTable[1] = "H";
mPeriodicTable[2] = "He";
mPeriodicTable[3] = "Li";
mPeriodicTable[4] = "Be";
mPeriodicTable[5] = "B";
mPeriodicTable[6] = "C";
mPeriodicTable[7] = "N";
mPeriodicTable[8] = "O";
mPeriodicTable[9] = "F";
mPeriodicTable[10] = "Ne";
mPeriodicTable[11] = "Na";
mPeriodicTable[12] = "Mg";
mPeriodicTable[13] = "Al";
mPeriodicTable[14] = "Si";
mPeriodicTable[15] = "P";
mPeriodicTable[16] = "S";
mPeriodicTable[17] = "Cl";
mPeriodicTable[18] = "Ar";
mPeriodicTable[19] = "K";
mPeriodicTable[20] = "Ca";
mPeriodicTable[21] = "Sc";
mPeriodicTable[22] = "Ti";
mPeriodicTable[23] = "V";
mPeriodicTable[24] = "Cr";
mPeriodicTable[25] = "Mn";
mPeriodicTable[26] = "Fe";
mPeriodicTable[27] = "Co";
mPeriodicTable[28] = "Ni";
mPeriodicTable[29] = "Cu";
mPeriodicTable[30] = "Zn";
mPeriodicTable[31] = "Ga";
mPeriodicTable[32] = "Ge";
mPeriodicTable[33] = "As";
mPeriodicTable[34] = "Se";
mPeriodicTable[35] = "Br";
mPeriodicTable[36] = "Kr";
mPeriodicTable[37] = "Rb";
mPeriodicTable[38] = "Sr";
mPeriodicTable[39] = "Y";
mPeriodicTable[40] = "Zr";
mPeriodicTable[41] = "Nb";
mPeriodicTable[42] = "Mo";
mPeriodicTable[43] = "Tc";
mPeriodicTable[44] = "Ru";
mPeriodicTable[45] = "Rh";
mPeriodicTable[46] = "Pd";
mPeriodicTable[47] = "Ag";
mPeriodicTable[48] = "Cd";
mPeriodicTable[49] = "In";
mPeriodicTable[50] = "Sn";
mPeriodicTable[51] = "Sb";
mPeriodicTable[52] = "Te";
mPeriodicTable[53] = "I";
mPeriodicTable[54] = "Xe";
mPeriodicTable[55] = "Cs";
mPeriodicTable[56] = "Ba";
mPeriodicTable[57] = "La";
mPeriodicTable[58] = "Ce";
mPeriodicTable[59] = "Pr";
mPeriodicTable[60] = "Nd";
mPeriodicTable[61] = "Pm";
mPeriodicTable[62] = "Sm";
mPeriodicTable[63] = "Eu";
mPeriodicTable[64] = "Gd";
mPeriodicTable[65] = "Tb";
mPeriodicTable[66] = "Dy";
mPeriodicTable[67] = "Ho";
mPeriodicTable[68] = "Er";
mPeriodicTable[69] = "Tm";
mPeriodicTable[70] = "Yb";
mPeriodicTable[71] = "Lu";
mPeriodicTable[72] = "Hf";
mPeriodicTable[73] = "Ta";
mPeriodicTable[74] = "W";
mPeriodicTable[75] = "Re";
mPeriodicTable[76] = "Os";
mPeriodicTable[77] = "Ir";
mPeriodicTable[78] = "Pt";
mPeriodicTable[79] = "Au";
mPeriodicTable[80] = "Hg";
mPeriodicTable[81] = "Tl";
mPeriodicTable[82] = "Pb";
mPeriodicTable[83] = "Bi";
mPeriodicTable[84] = "Po";
mPeriodicTable[85] = "At";
mPeriodicTable[86] = "Rn";
mPeriodicTable[87] = "Fr";
mPeriodicTable[88] = "Ra";
mPeriodicTable[89] = "Ac";
mPeriodicTable[90] = "Th";
mPeriodicTable[91] = "Pa";
mPeriodicTable[92] = "U";
mPeriodicTable[93] = "Np";
mPeriodicTable[94] = "Pu";
mPeriodicTable[95] = "Am";
mPeriodicTable[96] = "Cm";
mPeriodicTable[97] = "Bk";
mPeriodicTable[251] = "Cf";
mPeriodicTable[252] = "Es";
mPeriodicTable[100] = "Fm";
mPeriodicTable[258] = "Md";
mPeriodicTable[102] = "No";
mPeriodicTable[103] = "Lr";
mPeriodicTable[261] = "Rf";
mPeriodicTable[105] = "Db";
mPeriodicTable[106] = "Sg";
mPeriodicTable[107] = "Bh";
mPeriodicTable[108] = "Hs";
mPeriodicTable[109] = "Mt";
//
// Register parameters (ptr, name, default)
//
registerParameter(&mUserInt, "userInt", 0);
registerParameter(&mUserDouble, "userDouble", 0.);
registerParameter(&mUserString, "userString", string(""));
registerParameter(&mA, "A", static_cast(1));
registerParameter(&mVectorMesonId, "vectorMesonId", 443); // J/psi
registerParameter(&mDipoleModelName, "dipoleModel", string("bSat"));
registerParameter(&mDipoleModelParameterSetName, "dipoleModelParameterSet", string("KMW"));
registerParameter(&mTableSetTypeName, "tableSetType", string("total_and_coherent"));
registerParameter(&mQ2min, "Q2min", 1000.); // no limits if max <= min
registerParameter(&mQ2max, "Q2max", 0.);
registerParameter(&mWmin, "Wmin", 1000.);
registerParameter(&mWmax, "Wmax", 0.);
registerParameter(&mXpomMin, "xpomMin", 1e-8);
registerParameter(&mXpomMax, "xpomMax", 1.);
registerParameter(&mVerbose, "verbose", false);
registerParameter(&mVerboseLevel, "verboseLevel", 0);
registerParameter(&mRootfile, "rootfile", string("sartre.root"));
registerParameter(&mSeed, "seed", static_cast(time(0)));
registerParameter(&mUPC, "UPC", false);
registerParameter(&mUPCA, "UPCA", static_cast(1));
}
Settings::~Settings()
{
for (unsigned int k=0; kSetSeed(mSeed); // needed for TH1::GetRandom()
}
bool Settings::readSettingsFromFile(const char *file)
{
if (!file) return false; // nothing to do
mRuncard = string(file);
mLines.clear();
//
// Open file
//
ifstream ifs(file);
if (!ifs) {
cout << "Settings::readSettingsFromFile(): error, cannot open file '"
<< file << "'." << endl;
return false;
}
//
// Read file into vector of strings, skip comments and empty lines
//
while (ifs.good() && !ifs.eof()) {
string line;
getline (ifs, line);
if (ifs.eof() && line.empty()) break;
// empty line
if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) continue;
// if first character is not a letter/digit, then taken to be a comment.
int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
if (!isalnum(line[firstChar])) continue;
mLines.push_back(line);
}
ifs.close(); // done with I/O
//
// Process vector of strings one at a time and use
// it to set registered variables.
//
for (unsigned int i=0; i> name;
// find value string
string valueString;
splitLine >> valueString;
if (!splitLine) {
cout << "Settings::readSettingsFromFile(): error, value of variable '"
<< name.c_str() << "' not recognized." << endl;
}
istringstream modeData(valueString);
//
// Loop over registered variables and see which fits.
// Not particular elegant but does the job and saves
// a lot of programming in derived classes.
//
bool isRegistered = false;
for (unsigned int k=0; k>)) { // test
SettingsParameter> *var = dynamic_cast>*> (mRegisteredParameters[k]);
if (var->name == name) {
var->address->push_back(atof(valueString.c_str())); // first value
while (splitLine.good() && !splitLine.eof()) { // get remaining
string nextValue;
splitLine >> nextValue;
var->address->push_back(atof(nextValue.c_str()));
if (splitLine.eof()) break;
}
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
if (var->name == name) {
*(var->address) = valueString;
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
if (var->name == name) {
isRegistered = true;
if (valueString == string("true") ||
valueString == string("True") ||
valueString == string("TRUE") ||
valueString == string("on") ||
valueString == string("On") ||
valueString == string("ON") ||
valueString == string("Yes") ||
valueString == string("yes") ||
valueString == string("YES") ||
valueString == string("T") ||
valueString == string("t") ||
valueString == string("1") ) {
*(var->address) = true;
}
else {
*(var->address) = false;
}
}
}
}
if (!isRegistered) {
cout << "Settings::readSettingsFromFile(): error, parameter identifier '" <<
name.c_str() << "' not recognized." << endl;
}
}
//
// Consolidate input (after burner)
//
consolidateCommonSettings();
consolidateSettings(); // overloaded
return true;
}
bool Settings::list(ostream& os)
{
const int fieldWidth = 28;
os << "\nRun Settings:" << endl;
for (unsigned int k=0; k)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << var->address->c_str() << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter)) {
SettingsParameter *var = dynamic_cast*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << (*(var->address) ? "true" : "false") << endl;
}
}
os << endl;
return true;
}
TParticlePDG* Settings::lookupPDG(int id) const
{
if (mPDG)
return mPDG->GetParticle(id);
else
return 0;
}
string Settings::particleName(int pdgID)
{
string name("unknown");
if (abs(pdgID) < 1000000000) { // particle
if (mPDG) {
TParticlePDG *part = lookupPDG(pdgID);
if (part) name = part->GetName();
}
if (abs(pdgID) == 990) name = "pomeron";
}
else { // nucleus in 10LZZZAAAI PDG format
int id = pdgID;
// int iso = id%10;
id /= 10;
int A = id%1000;
id /= 1000;
int Z = id%1000;
stringstream namestream;
namestream << mPeriodicTable[Z] << "(" << A << ")";
name = namestream.str();
}
return name;
}
void Settings::setVerbose(bool val) {
mVerbose = val;
if (mVerbose && mVerboseLevel == 0) mVerboseLevel = 1;
if (!mVerbose && mVerboseLevel != 0) mVerboseLevel = 0;
}
bool Settings::verbose() const {return mVerbose;}
void Settings::setVerboseLevel(int val) {
mVerboseLevel = val;
if (!mVerbose && mVerboseLevel != 0) mVerbose = true;
if (mVerbose && mVerboseLevel == 0) mVerbose = false;
}
int Settings::verboseLevel() const {return mVerboseLevel;}
void Settings::setQ2min(double val) { mQ2min = val;}
double Settings::Q2min() const {return mQ2min;}
double Settings::Qmin() const {return sqrt(mQ2min);}
void Settings::setQ2max(double val) { mQ2max = val;}
double Settings::Q2max() const {return mQ2max;}
double Settings::Qmax() const {return sqrt(mQ2max);}
void Settings::setW2min(double val) { mWmin = sqrt(val);}
void Settings::setWmin(double val) { mWmin = val;}
double Settings::Wmin() const {return mWmin;}
double Settings::W2min() const {return mWmin*mWmin;}
void Settings::setW2max(double val) { mWmax = sqrt(val);}
void Settings::setWmax(double val) { mWmax = val;}
double Settings::Wmax() const {return mWmax;}
double Settings::W2max() const {return mWmax*mWmax;}
void Settings::setXpomMin(double val) {mXpomMin = val;}
void Settings::setXpomMax(double val) {mXpomMax = val;}
double Settings::xpomMin() const {return mXpomMin;}
double Settings::xpomMax() const {return mXpomMax;}
int Settings::vectorMesonId() const {return mVectorMesonId;}
void Settings::setVectorMesonId(int val) {mVectorMesonId = val;}
string Settings::dipoleModelName() const {return mDipoleModelName;}
DipoleModelType Settings::dipoleModelType() const {return mDipoleModelType;}
void Settings::setDipoleModelType(DipoleModelType val)
{
mDipoleModelType = val;
if (mDipoleModelType == bSat)
mDipoleModelName = string("bSat");
else if (mDipoleModelType == bNonSat)
mDipoleModelName = string("bNonSat");
else if (mDipoleModelType == bCGC)
mDipoleModelName = string("bCGC");
}
unsigned int Settings::A() const {return mA;}
void Settings::setA(unsigned int val) {mA = val;}
void Settings::setRootfile(const char* val){ mRootfile = val; }
string Settings::rootfile() const { return mRootfile; }
string Settings::dipoleModelParameterSetName() const {return mDipoleModelParameterSetName;}
DipoleModelParameterSet Settings::dipoleModelParameterSet() const {return mDipoleModelParameterSet;}
void Settings::setDipoleModelParameterSet(DipoleModelParameterSet val)
{
mDipoleModelParameterSet = val;
if (mDipoleModelParameterSet == KMW)
mDipoleModelParameterSetName = string("KMW");
else if (mDipoleModelParameterSet == HMPZ)
mDipoleModelParameterSetName = string("HMPZ");
+ else if (mDipoleModelParameterSet == BSTT)
+ mDipoleModelParameterSetName = string("BSTT");
else if (mDipoleModelParameterSet == CUSTOM)
mDipoleModelParameterSetName = string("CUSTOM");
}
string Settings::tableSetTypeName() const {return mTableSetTypeName;}
TableSetType Settings::tableSetType() const {return mTableSetType;}
void Settings::setTableSetType(TableSetType val)
{
mTableSetType = val;
if (mTableSetType == total_and_coherent)
mTableSetTypeName = string("total_and_coherent");
else if (mTableSetType == coherent_and_incoherent)
mTableSetTypeName = string("coherent_and_incoherent");
}
void Settings::consolidateCommonSettings()
{
//
// Check if verbose levels and flags are consistent
// The verbose flag superseeds the verboseLevel.
//
if (mVerbose && mVerboseLevel == 0) mVerboseLevel = 1;
if (mVerboseLevel != 0 && !mVerbose) mVerboseLevel = 0;
//
// Set random generator seed
//
mRandomGenerator.SetSeed(mSeed);
gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()
//
// Dipole Model
//
if (mDipoleModelName == string("bSat"))
mDipoleModelType = bSat;
else if (mDipoleModelName == string("bNonSat"))
mDipoleModelType = bNonSat;
else if (mDipoleModelName == string("bCGC"))
mDipoleModelType = bCGC;
else {
cout << "Settings::consolidateCommonSettings(): Error, dipole model '"
<< mDipoleModelName.c_str() << "' is not defined." << endl;
exit(1);
}
//
// Dipole Model Parameter Set
//
if (mDipoleModelParameterSetName == string("KMW"))
mDipoleModelParameterSet = KMW;
else if (mDipoleModelParameterSetName == string("HMPZ"))
+ mDipoleModelParameterSet = BSTT;
+ else if (mDipoleModelParameterSetName == string("BSTT"))
mDipoleModelParameterSet = HMPZ;
else if (mDipoleModelParameterSetName == string("CUSTOM"))
mDipoleModelParameterSet = CUSTOM;
else {
cout << "Settings::consolidateCommonSettings(): Error, dipole model parameter set'"
<< mDipoleModelParameterSetName.c_str() << "' is not defined." << endl;
exit(1);
}
//
// Table Set Type
//
if (mTableSetTypeName == string("total_and_coherent"))
mTableSetType = total_and_coherent;
else if (mTableSetTypeName == string("coherent_and_incoherent"))
mTableSetType = coherent_and_incoherent;
else {
cout << "Settings::consolidateCommonSettings(): Error, table set type '"
<< mTableSetTypeName.c_str() << "' is not defined." << endl;
exit(1);
}
}
void Settings::setUPC(bool val){ mUPC = val; }
bool Settings::UPC() const { return mUPC; }
void Settings::setUPCA(unsigned int val){ mUPCA = val; }
unsigned int Settings::UPCA() const { return mUPCA; }
Index: trunk/src/DipoleModel.h
===================================================================
--- trunk/src/DipoleModel.h (revision 392)
+++ trunk/src/DipoleModel.h (revision 393)
@@ -1,104 +1,99 @@
//==============================================================================
// DipoleModel.h
//
// 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 double dsigmadb2ep(double, double, double); // ep
+ virtual double coherentDsigmadb2(double, double, double) {return 0;}; // eA
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);
+ void createConfiguration(int);
double dsigmadb2(double, double, double, double);
double bDependence(double, double);
double dsigmadb2ep(double, double, double);
-
+ double coherentDsigmadb2(double, double, double);
+
protected:
TH2F* mBDependence;
private:
double dsigmadb2epForIntegration(double*, double*);
- TH1F* mSigma_ep_LookupTable;
- double coherentDsigmadb2(double, double, double);
-
+ TH1F* mSigma_ep_LookupTable;
};
-class DipoleModel_bNonSat : public DipoleModel_bSat{
+class DipoleModel_bNonSat : public DipoleModel_bSat {
public:
DipoleModel_bNonSat();
- ~DipoleModel_bNonSat();
+ ~DipoleModel_bNonSat();
-private:
- double dsigmadb2(double, double, double, double);
+ 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);
+ void createConfiguration(int);
double dsigmadb2(double, double, double, double);
double dsigmadb2ep(double, double, double);
double bDependence(double);
};
#endif
Index: trunk/src/Enumerations.h
===================================================================
--- trunk/src/Enumerations.h (revision 392)
+++ trunk/src/Enumerations.h (revision 393)
@@ -1,33 +1,33 @@
//==============================================================================
// Enumerations.h
//
// 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 DipoleModelParameterSet {KMW, HMPZ, BSTT, CUSTOM};
enum TableSetType {total_and_coherent, coherent_and_incoherent};
#endif
Index: trunk/src/DipoleModelParameters.cpp
===================================================================
--- trunk/src/DipoleModelParameters.cpp (revision 392)
+++ trunk/src/DipoleModelParameters.cpp (revision 393)
@@ -1,474 +1,504 @@
//==============================================================================
// DipoleModelParameters.cpp
//
// 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 == BSTT) {
+ // Similar as HMPZ but with modified dipole sizes, http://arxiv.org/abs/0807.0325v1
+ cout << "DipoleModelParameters::setup_bSat(): BSTT parameters to be defined in future. Doesn't work yet." << endl;
+ exit(2);
+ mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.0; // u,d,s quarks
+ mQuarkMass[3] = 0; // c quark
+
+ mBG = 0;
+ mMu02 = 0; // Gev^2
+ mLambdaG = 0;
+ mAg = 0;
+ mC = 0;
+ mRMax = 0;
+ }
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 == BSTT) {
+ // Similar as HMPZ but with modified dipole sizes, http://arxiv.org/abs/0807.0325v1
+ cout << "DipoleModelParameters::setup_bNonSat(): BSTT parameters to be defined in future. Doesn't work yet." << endl;
+ exit(2);
+ mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.0; // u,d,s quarks
+ mQuarkMass[3] = 0; // c quark
+
+ mBG = 0;
+ mMu02 = 0; // Gev^2
+ mLambdaG = 0;
+ mAg = 0;
+ mC = 0;
+ mRMax = 0;
}
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;
+ mRMax = 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.
+ // and also do not distinguish 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) {
+ else if (mDipoleModelParameterSet == HMPZ || mDipoleModelParameterSet == BSTT) {
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;
+ if (mDipoleModelParameterSet == BSTT)
+ os << setw(fieldWidth) << "rMax: " << mRMax << 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/DipoleModel.cpp
===================================================================
--- trunk/src/DipoleModel.cpp (revision 392)
+++ trunk/src/DipoleModel.cpp (revision 393)
@@ -1,313 +1,337 @@
//==============================================================================
// DipoleModel.cpp
//
// 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;
+ 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;
+ 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 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));
+ 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;
+
+ }
+ return *this;
}
-
+
DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp)
{
- mBDependence = new TH2F(*(dp.mBDependence));
+ 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;
+ 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);
+ if (mParameters->dipoleModelParameterSet() == BSTT) {
+ double rm = mParameters->rMax();
+ r = rm*sqrt(log(1+r*r/(rm*rm)));
+ }
+ 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 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 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
+ if (mParameters->dipoleModelParameterSet() == BSTT) {
+ double rm = mParameters->rMax();
+ r = rm*sqrt(log(1+r*r/(rm*rm)));
+ }
+ 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 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;
+double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe) {
+ if (mParameters->dipoleModelParameterSet() == BSTT) {
+ double rm = mParameters->rMax();
+ r = rm*sqrt(log(1+r*r/(rm*rm)));
+ }
+ 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 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);
+ 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.
+ // 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
+{
+ if (mParameters->dipoleModelParameterSet() == BSTT) {
+ double rm = mParameters->rMax();
+ r = rm*sqrt(log(1+r*r/(rm*rm)));
+ }
+ 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;
+
+ return omega;
}
-
+
double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
- double bDep=bDependence(b, phi);
+ if (mParameters->dipoleModelParameterSet() == BSTT) {
+ double rm = mParameters->rMax();
+ r = rm*sqrt(log(1+r*r/(rm*rm)));
+ }
+ 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;
+
+ 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();
+ if (mParameters->dipoleModelParameterSet() == BSTT) {
+ double rm = mParameters->rMax();
+ r = rm*sqrt(log(1+r*r/(rm*rm)));
+ }
+ 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;
+ 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;
+ 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 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);
+ return exp(-0.5*b*b/Bcgc/gammas/hbarc2);
}
-
+