Index: trunk/src/TableGeneratorSettings.h
===================================================================
--- trunk/src/TableGeneratorSettings.h (revision 369)
+++ trunk/src/TableGeneratorSettings.h (revision 370)
@@ -1,122 +1,122 @@
//==============================================================================
// TableGeneratorSettings.h
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Singleton class
//
//==============================================================================
#ifndef TableGeneratorSettings_h
#define TableGeneratorSettings_h
#include "Settings.h"
#include "Enumerations.h"
using namespace std;
class TableGeneratorSettings : public Settings {
public:
static TableGeneratorSettings* instance();
void setTmin(double);
double tmin() const;
void setTmax(double);
double tmax() const;
void setXmin(double);
double xmin() const;
void setXmax(double);
double xmax() const;
void setQ2bins(unsigned int);
unsigned int Q2bins() const;
void setW2bins(unsigned int);
unsigned int W2bins() const;
void setTbins(unsigned int);
unsigned int tbins() const;
void setXbins(unsigned int);
unsigned int xbins() const;
string bSatLookupPath() const;
void setBSatLookupPath(string);
unsigned int numberOfConfigurations() const;
void setNumberOfConfigurations(unsigned int);
vector dipoleModelCustomParameters() const;
bool useBackupFile() const;
void setUseBackupFile(bool);
int startingBinFromBackup() const;
void setStartingBinFromBackup(int);
int startBin() const;
void setStartBin(int);
int endBin() const;
void setEndBin(int);
int modesToCalculate() const;
void setModesToCalculate(int);
unsigned char priority() const;
void setPriority(unsigned char);
void consolidateSettings();
private:
TableGeneratorSettings();
private:
static TableGeneratorSettings* mInstance;
private:
unsigned int mNumberOfConfigurations;
vector mDipoleModelCustomParameters; // developer
bool mUseBackupFile;
int mStartingBinFromBackup;
int mStartBin, mEndBin;
int mModesToCalculate;
double mTmin;
double mTmax;
double mXmin;
double mXmax;
unsigned int mQ2bins;
unsigned int mW2bins;
unsigned int mTbins;
unsigned int mXbins;
string mBSatLookupPath;
int mPriority;
};
#endif
Index: trunk/src/Table.cpp
===================================================================
--- trunk/src/Table.cpp (revision 369)
+++ trunk/src/Table.cpp (revision 370)
@@ -1,1431 +1,1431 @@
//==============================================================================
// Table.cpp
//
-// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// 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 (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_A) 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 (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_A) 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_A;
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() == 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/FrangibleNucleus.cpp
===================================================================
--- trunk/src/FrangibleNucleus.cpp (revision 369)
+++ trunk/src/FrangibleNucleus.cpp (revision 370)
@@ -1,204 +1,204 @@
//==============================================================================
// FrangibleNucleus.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "FrangibleNucleus.h"
#include "CNucleus.h"
#include "Constants.h"
#include "EventGeneratorSettings.h"
#include "TRandom3.h"
#define PR(x) cout << #x << " = " << (x) << endl;
FrangibleNucleus::FrangibleNucleus()
{
mGeminiNucleus = 0;
mExcitationEnergy = 0;
}
FrangibleNucleus& FrangibleNucleus::operator=(const FrangibleNucleus& gn)
{
if (this != &gn) {
delete mGeminiNucleus;
mProducts.clear();
Nucleus::operator=(gn);
mExcitationEnergy = gn.mExcitationEnergy;
copy(gn.mProducts.begin(), gn.mProducts.end(), mProducts.begin());
if (gn.mGeminiNucleus)
mGeminiNucleus = new CNucleus(gn.mZ, gn.mA);
else
mGeminiNucleus = 0;
}
return *this;
}
FrangibleNucleus::FrangibleNucleus(const FrangibleNucleus& gn) : Nucleus(gn)
{
mExcitationEnergy = gn.mExcitationEnergy;
copy(gn.mProducts.begin(), gn.mProducts.end(), mProducts.begin());
if (gn.mGeminiNucleus)
mGeminiNucleus = new CNucleus(gn.mZ, gn.mA);
else
mGeminiNucleus = 0;
}
void FrangibleNucleus::init(unsigned int A)
{
Nucleus::init(A);
if (mGeminiNucleus) delete mGeminiNucleus;
mGeminiNucleus = 0;
}
void FrangibleNucleus::init(unsigned int A, bool enableBreakup)
{
Nucleus::init(A);
if (mGeminiNucleus) delete mGeminiNucleus;
//
// Init Gemini.
//
if (enableBreakup)
mGeminiNucleus = new CNucleus(mZ, mA);
else
mGeminiNucleus = 0;
mExcitationEnergy = 0;
}
FrangibleNucleus::FrangibleNucleus(unsigned int A, bool enableBreakup) : Nucleus(A)
{
//
// Init Gemini.
// Only if requested, not needed otherwise.
//
if (enableBreakup)
mGeminiNucleus = new CNucleus(mZ, mA);
else
mGeminiNucleus = 0;
mExcitationEnergy = 0;
}
FrangibleNucleus::~FrangibleNucleus()
{
delete mGeminiNucleus;
}
void FrangibleNucleus::resetBreakup()
{
mExcitationEnergy = 0;
mProducts.clear();
}
int FrangibleNucleus::breakup(const TLorentzVector& dissSystem)
{
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
//
// Estimate excitation energy
// Note that the dissSystem is given in units of 'per nucleon'.
// Hence we have to multiply the excitation energy with A.
//
mExcitationEnergy = (dissSystem.M()-protonMass)*mA;
double Ex = mExcitationEnergy;
double maxEx = settings->maxNuclearExcitationEnergy();
if (Ex > maxEx) {
if (settings->verboseLevel() > 1)
cout << "FrangibleNucleus::breakup(): Actual excitation energy (" << Ex
<< ") exceeded upper limit (" << maxEx << "). Reset to maximum value." << endl;
Ex = maxEx;
}
Ex *= 1000; // Gemini uses the total energy in MeV
//
// Setup excited nucleus
//
// Pass excitation energy and spin to Gemini
mGeminiNucleus->setCompoundNucleus(Ex,mSpin); //specify the excitation energy and spin
mGeminiNucleus->setVelocityCartesian(); // set initial velocity to zero (CMS)
// Set the direction of the spin vector (random)
TRandom3 *rndm = Settings::randomGenerator();
double phi = rndm->Uniform(2*M_PI);
double theta = acos(rndm->Uniform(-1, 1));
CAngle spin(theta, phi);
mGeminiNucleus->setSpinAxis(spin);
//
// Let the nucleus breakup
//
mGeminiNucleus->decay();
mProducts.clear();
if (mGeminiNucleus->abortEvent) {
cout << "Nucleus::breakup(): Error, decay aborted in Gemini++." << endl;
mGeminiNucleus->reset();
return 0;
}
int nfragments = mGeminiNucleus->getNumberOfProducts();
mProducts.resize(nfragments);
for(int i=0; igetProducts(i); //set pointer to first stable product
mProducts[i].Z = product->iZ;
mProducts[i].A = product->iA;
mProducts[i].emissionTime = product->getTime();
mProducts[i].p = product->getLorentzVector();
mProducts[i].p.Boost(dissSystem.BoostVector());
mProducts[i].name = product->getName();
if (product->iZ == 1 && product->iA == 1)
mProducts[i].pdgId = 2212; // p
else if (product->iZ == 0 && product->iA == 1)
mProducts[i].pdgId = 2112; // n
else
mProducts[i].pdgId = pdgID(product->iZ, product->iA);
}
mGeminiNucleus->reset();
return nfragments;
}
const vector& FrangibleNucleus::breakupProducts() const
{
return mProducts;
}
void FrangibleNucleus::listBreakupProducts(ostream& os) const
{
TLorentzVector sys;
cout << "Excitation energy = " << mExcitationEnergy << endl;
cout << "Number of fragments = " << mProducts.size() << endl;
for (unsigned int i=0; i.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include
#include
#include
#include "Integrals.h"
#include "Constants.h"
#include "Nucleus.h"
#include "DipoleModel.h"
#include "AlphaStrong.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/Functor.h"
#include "TMath.h"
#include "WaveOverlap.h"
#include "Kinematics.h"
#include "TableGeneratorSettings.h"
#include "Enumerations.h"
#include "IntegrandWrappers.h"
#include "TF1.h"
#include "TH1F.h"
#include "cuba.h"
#define PRf(x) printf(#x); printf("=%g \n", (x));
#define PRi(x) printf(#x); printf("=%d \n", (x));
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
Integrals::Integrals()
{
mIsInitialized=false;
mRelativePrecisionOfIntegration = 0;
mWaveOverlap = 0;
mDipoleModel = 0;
mDipoleModelForSkewednessCorrection = 0;
mIntegralImT = 0;
mIntegralImL = 0;
mIntegralReT = 0;
mIntegralReL = 0;
mErrorImT = 0;
mErrorImL = 0;
mErrorReT = 0;
mErrorReL = 0;
mProbImT = 0;
mProbImL = 0;
mProbReT = 0;
mProbReL = 0;
mIntegralTForSkewedness = 0;
mIntegralLForSkewedness = 0;
mErrorTForSkewedness = 0;
mErrorLForSkewedness = 0;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mVerbose=settings->verbose();
int VMId=settings->vectorMesonId();
mMV = settings->lookupPDG(VMId)->Mass();
mIsUPC = settings->UPC();
if (VMId==113 || VMId==333 || VMId == 443 || VMId == 553) {
mWaveOverlap = new WaveOverlapVM;
mWaveOverlap->setProcess(VMId);
mWaveOverlap->setWaveOverlapFunctionParameters(VMId);
if(mVerbose)
mWaveOverlap->testBoostedGaussianParameters(VMId);
}
else if (VMId==22) {
mWaveOverlap = new WaveOverlapDVCS;
}
else {
cout << "Integrals::init(): Error, no exclusive production implemented for: "<< VMId << endl;
exit(1);
}
DipoleModelType model=settings->dipoleModelType();
if (model==bSat) {
mDipoleModel = new DipoleModel_bSat;
}
else if(model==bNonSat){
mDipoleModel = new DipoleModel_bNonSat;
}
else if (model==bCGC) {
mDipoleModel = new DipoleModel_bCGC;
}
else {
cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
exit(1);
}
mCalculateSkewedness=false;
if(settings->A()==1 and settings->modesToCalculate()!=1 and settings->numberOfConfigurations()==1 and model==bSat) {
mCalculateSkewedness=true;
mDipoleModelForSkewednessCorrection = new DipoleModel_bNonSat;
}
mIsInitialized=true;
}
Integrals::Integrals(const Integrals& integrals)
{
mIsInitialized = integrals.mIsInitialized;
mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS))
mWaveOverlap = new WaveOverlapDVCS;
else
mWaveOverlap = new WaveOverlapVM;
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat))
mDipoleModel = new DipoleModel_bSat;
else if(typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
else
mDipoleModel = new DipoleModel_bCGC;
if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mErrorImT = integrals.mErrorImT;
mErrorImL = integrals.mErrorImL;
mErrorReT = integrals.mErrorReT;
mErrorReL = integrals.mErrorReL;
mProbImT = integrals.mProbImT;
mProbImL = integrals.mProbImL;
mProbReT = integrals.mProbReT;
mProbReL = integrals.mProbReL;
mMV = integrals.mMV;
}
Integrals& Integrals::operator=(const Integrals& integrals)
{
if (this != &integrals) {
delete mWaveOverlap;
delete mDipoleModel;
delete mDipoleModelForSkewednessCorrection;
if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS))
mWaveOverlap = new WaveOverlapDVCS;
else
mWaveOverlap = new WaveOverlapVM;
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat))
mDipoleModel = new DipoleModel_bSat;
else if (typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
else
mDipoleModel = new DipoleModel_bCGC;
if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mErrorImT = integrals.mErrorImT;
mErrorImL = integrals.mErrorImL;
mErrorReT = integrals.mErrorReT;
mErrorReL = integrals.mErrorReL;
mProbImT = integrals.mProbImT;
mProbImL = integrals.mProbImL;
mProbReT = integrals.mProbReT;
mProbReL = integrals.mProbReL;
mMV = integrals.mMV;
mIsInitialized = integrals.mIsInitialized;
mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
}
return *this;
}
Integrals::~Integrals()
{
delete mWaveOverlap;
delete mDipoleModel;
if(mDipoleModelForSkewednessCorrection)
delete mDipoleModelForSkewednessCorrection;
}
void Integrals::operator() (double t, double Q2, double W2)
{
unsigned int A=dipoleModel()->nucleus()->A();
//make sure the configurations have been generated:
if (!mDipoleModel->configurationExists() and A!=1) {
// do not use cout
cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
exit(1);
}
if (setKinematicPoint(t, Q2, W2)) {
if(A==1){
calculateEp();
if(mCalculateSkewedness){
calculateSkewedness();
}
}
else
calculate();
}
else {
fillZeroes();
}
}
void Integrals::operator() (double t, double xpom) //UPC
{
unsigned int A=dipoleModel()->nucleus()->A();
//make sure the configurations have been generated:
if (!mDipoleModel->configurationExists() and A!=1) {
// do not use cout
cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
exit(1);
}
if (setKinematicPoint(t, xpom)) {
if(A==1){
calculateEp();
if(mCalculateSkewedness){
calculateSkewedness();
}
}
else
calculate();
}
else {
fillZeroes();
}
}
void Integrals::fillZeroes(){
//Store the results
mIntegralImT=0;
mIntegralReT=0;
mIntegralImL=0;
mIntegralReL=0;
//Store the errors:
mErrorImT=0;
mErrorImL=0;
mErrorReT=0;
mErrorReL=0;
//Store the probabilities:
mProbImT=0;
mProbImL=0;
mProbReT=0;
mProbReL=0;
}
//*********EXCLUSIVE VECTOR MESONS OR DVCS: ********************************
void IntegralsExclusive::coherentIntegrals(double t, double Q2, double W2)
{
if(setKinematicPoint(t, Q2, W2)){
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){
//store present kinematic point:
double xprobe=kinematicPoint[3];
dipoleModel()->createSigma_ep_LookupTable(xprobe);
}
calculateCoherent();
}
else
fillZeroes();
}
void IntegralsExclusive::coherentIntegrals(double t, double xpom)
{
if(setKinematicPoint(t, xpom)){
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){
dipoleModel()->createSigma_ep_LookupTable(xpom);
}
calculateCoherent();
}
else
fillZeroes();
}
IntegralsExclusive::IntegralsExclusive()
{
}
IntegralsExclusive& IntegralsExclusive::operator=(const IntegralsExclusive& cobj)
{
if (this != &cobj) {
Integrals::operator=(cobj);
copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
}
return *this;
}
IntegralsExclusive::IntegralsExclusive(const IntegralsExclusive& cobj) : Integrals(cobj)
{
copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
}
bool IntegralsExclusive::setKinematicPoint(double t, double xpom) //UPC
{
bool result = true;
kinematicPoint[0]=t;
kinematicPoint[1]=0; //Q2
kinematicPoint[2]=0; //W2 is not used
kinematicPoint[3]=xpom;
return result;
}
bool IntegralsExclusive::setKinematicPoint(double t, double Q2, double W2)
{
bool result = true;
kinematicPoint[0]=t;
kinematicPoint[1]=Q2;
kinematicPoint[2]=W2;
double xprobe=Kinematics::xpomeron(t, Q2, W2, mMV);
if (xprobe<0 || xprobe>1)
result = false;
kinematicPoint[3]=xprobe;
return result;
}
void IntegralsExclusive::calculate()
{
//
// This function calls a wrapper from where the
// integral is calculated with the Cuhre method.
// Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
//
const double epsrel=1.e-2, epsabs=1e-12;
const int flags=0, mineval=1e4, maxeval=1e9, key=0;
int nregionsTIm, nevalTIm, failTIm;
int nregionsTRe, nevalTRe, failTRe;
int nregionsLIm, nevalLIm, failLIm;
int nregionsLRe, nevalLRe, failLRe;
double valTIm=0, errTIm=0, probTIm=0;
double valLIm=0, errLIm=0, probLIm=0;
double valTRe=0, errTRe=0, probTRe=0;
double valLRe=0, errLRe=0, probLRe=0;
const char* statefile=0;
const int nvec=1;
// double probabilityCutOff=1e-6;
//
// Do the integrations
//
Cuhre(4, 1, integrandWrapperTIm, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsTIm, &nevalTIm, &failTIm, &valTIm, &errTIm, &probTIm);
if(failTIm!=0 and mVerbose)
printf("IntegralsExclusive::calculate(): Warning: Integration TIm did not reach desired precision! Error code=%d \n", failTIm);
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
Cuhre(4, 1, integrandWrapperLIm, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsLIm, &nevalLIm, &failLIm, &valLIm, &errLIm, &probLIm);
if(failLIm!=0 and mVerbose)
printf("IntegralsExclusive::calculate(): Warning: Integration LIm did not reach desired precision! Error code=%d \n", failLIm);
}
Cuhre(4, 1, integrandWrapperTRe, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsTRe, &nevalTRe, &failTRe, &valTRe, &errTRe, &probTRe);
if(failTRe!=0 and mVerbose)
printf("IntegralsExclusive::calculate(): Warning: Integration TRe did not reach desired precision! Error code=%d \n", failTRe);
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
Cuhre(4, 1, integrandWrapperLRe, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsLRe, &nevalLRe, &failLRe, &valLRe, &errLRe, &probLRe);
if(failLRe!=0 and mVerbose)
printf("IntegralsExclusive::calculate(): Warning: Integration LRe did not reach desired precision! Error code=%d \n", failLRe);
}
//
// Store the results:
//
mIntegralImT=valTIm;
mIntegralReT=valTRe;
mIntegralImL=valLIm;
mIntegralReL=valLRe;
//
// Store the errors:
//
mErrorImT=errTIm;
mErrorImL=errLIm;
mErrorReT=errTRe;
mErrorReL=errLRe;
//
// Store the probabilities:
//
mProbImT=probTIm;
mProbImL=probLIm;
mProbReT=probTRe;
mProbReL=probLRe;
}
void IntegralsExclusive::calculateEp()
{
//
// This function calls a wrapper from where the
// integral is calculated with the Cuhre method.
// Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
//
const double epsrel=1.e-4, epsabs=1e-12;
const int flags=0, maxeval=1e9, key=0;
const int mineval=3e6;
int nregionsT, nevalT, failT;
int nregionsL, nevalL, failL;
double valT=0, errT=0, probT=0;
double valL=0, errL=0, probL=0;
const char* statefile=0;
const int nvec=1;
//
// Do the integrations
//
/*
bool bContinue=true;
bool isFirst=true;
while(bContinue){
double valTOld=valT;
Cuhre(4, 1, integrandWrapperTep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(abs(valT-valTOld)/valT > epsrel){
mineval*=3;
if(isFirst)
mineval*=30;
isFirst=false;
if(mineval>1e4)
PR(mineval);
}
else
bContinue=false;
}
*/
Cuhre(4, 1, integrandWrapperTep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(failT!=0 and mVerbose)
printf("IntegralsExclusive::calculateEp(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
/*
if(errT/valT > epsrel)
PR(errT/valT);
if(errT < epsabs)
PR(errT);
if(probT>0.5)
PR(probT);
PR(nevalT);
PR(nregionsT);
*/
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
Cuhre(4, 1, integrandWrapperLep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
if(failL!=0 and mVerbose)
printf("IntegralsExclusive::calculateEp(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
}
//
// Store the results:
//
mIntegralImT=valT;
mIntegralImL=valL;
//
// Store the errors:
//
mErrorImT=errT;
mErrorImL=errL;
//
// Store the probabilities:
//
mProbImT=probT;
mProbImL=probL;
}
void IntegralsExclusive::calculateSkewedness()
{
//
// This function calls a wrapper from where the
// integral is calculated with the Cuhre method.
// Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
//
const double epsrel=1.e-4, epsabs=1e-12;
const int flags=0, maxeval=1e9, key=0;
const int mineval=3e6;
int nregionsT, nevalT, failT;
int nregionsL, nevalL, failL;
double valT=0, errT=0, probT=0;
double valL=0, errL=0, probL=0;
const char* statefile=0;
const int nvec=1;
//
// Do the integrations
//
/*
bool bContinue=true;
bool isFirst=true;
while(bContinue){
double valTOld=valT;
Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(abs(valT-valTOld)/valT > epsrel){
mineval*=3;
if(isFirst)
mineval*=30;
isFirst=false;
if(mineval>1e4)
PR(mineval);
}
else
bContinue=false;
}
*/
Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(failT!=0 and mVerbose)
printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
//
// For UPC, calculate only transverse polarisation case
//
if(!mIsUPC){
Cuhre(4, 1, integrandWrapperLForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
if(failL!=0 and mVerbose)
printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
}
//
// Store the results:
//
mIntegralTForSkewedness=valT;
mIntegralLForSkewedness=valL;
//
// Store the errors:
//
mErrorTForSkewedness=errT;
mErrorLForSkewedness=errL;
//
// Store the probabilities:
//
mProbImTForSkewedness=probT;
mProbImLForSkewedness=probL;
}
void IntegralsExclusive::calculateCoherent()
{
//
// As calculate() but for the coherent case
//
const double epsrel=1e-6, epsabs=1e-12;
const int flags=0, mineval=1e1, maxeval=1e9, key=0;
int nregionsT, nevalT, failT;
int nregionsL, nevalL, failL;
double valT, errT, probT;
double valL, errL, probL;
const int nvec=1;
const char* statefile=0;
//
// Do the integrations
//
Cuhre(3, 1, integrandWrapperCoherentAmplitudeT, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
if(failT!=0 and mVerbose)
printf("IntegralsExclusive::calculateCoherent(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
Cuhre(3, 1, integrandWrapperCoherentAmplitudeL, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
if(failL!=0 and mVerbose)
printf("IntegralsExclusive::calculateCoherent(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
//
// Store the results
//
mIntegralImT=valT;
mIntegralImL=valL;
//
// Store the errors
//
mErrorImT=errT;
mErrorImL=errL;
}
//
// The following functions are the Integrands in the Amplitudes:
//
double IntegralsExclusive::uiAmplitudeTIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double cosArg = (b/hbarc)*Delta*cos(phi);
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2(r , b , phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*cos(cosArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiAmplitudeTRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double sinArg = b*Delta*cos(phi)/hbarc;
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*sin(sinArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiAmplitudeLIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double cosArg = b*Delta*cos(phi)/hbarc;
double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*cos(cosArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiAmplitudeLRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double sinArg = b*Delta*cos(phi)/hbarc;
double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*sin(sinArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiCoherentAmplitudeT(double b, double z, double r, double Q2, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double xprobe=kinematicPoint[3];
double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe);
double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean;
return result;
}
double IntegralsExclusive::uiCoherentAmplitudeL(double b, double z, double r, double Q2, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double xprobe=kinematicPoint[3];
double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe);
double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean;
return result;
}
double IntegralsExclusive::uiAmplitudeTep(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
double IntegralsExclusive::uiAmplitudeLep(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
//
// Only for calculating the lamdba for Skewedness Corrections:
//
double IntegralsExclusive::uiAmplitudeTForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
double IntegralsExclusive::uiAmplitudeLForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r, b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
Index: trunk/src/tableMerger.cpp
===================================================================
--- trunk/src/tableMerger.cpp (revision 369)
+++ trunk/src/tableMerger.cpp (revision 370)
@@ -1,558 +1,558 @@
//==============================================================================
// tableMerger.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Utility program to merge tables.
//
// Several checks are performed to ensure the integrity of the merged
// tables but in any case the merged table(s) should be examined with
// tableInspector.
//
// Usage:
// tableMerger [-n] [-0] [-p prefix] file(s) ...
// -n Do not create the new tables but perform all checks
// -p prefix Set prefix for output files
// -0 Do NOT copy bins with 0 content (ignore them)
// Use when merging partially filled tables
// -m List missing bins (if any)
//
// The output files are named _merged.root
// where ID is the unique identifier of the table.
//
// Note it is a good idea to always check with option -n on first.
//
//==============================================================================
#include "Table.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include "TFile.h"
#include "TH2F.h"
#include "TH3F.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
void usage(const char* prog)
{
cout << "Usage: " << prog << " [-n] [-0] [-p prefix] [-m] file(s) ..." << endl;
cout << " -n Do not create the new tables but perform all checks" << endl;
cout << " -p prefix Set prefix for output files" << endl;
cout << " -0 Do NOT copy bins with 0 content (ignore them)" << endl;
cout << " Use when merging partially filled tables" << endl;
cout << " -m List missing bins (if any)" << endl;
}
class Flags {
public:
bool listMissingBins;
bool noOutput;
bool ignoreZero;
string prefix;
};
bool isUPC(uint64_t id) {return (id & (static_cast(1) << 46));}
void createHistogramLists(vector& files, vector& all3DHistograms,
vector& all2DHistograms)
{
for (unsigned i=0; iGet("table");
if (ptr == 0) {
cout << "Error: failed retrieving table from file '" << files[i].c_str() << "'." << endl;
exit(1);
}
uint64_t id = atoll(ptr->GetTitle());
if (isUPC(id)) {
TH2F *histo = reinterpret_cast(ptr);
histo->SetDirectory(0);
all2DHistograms.push_back(histo);
}
else {
TH3F *histo = reinterpret_cast(ptr);
histo->SetDirectory(0);
all3DHistograms.push_back(histo);
}
file->Close();
}
}
int merge3DHistos(Flags &flags, vector& allHistograms)
{
set idSet;
for (unsigned i=0; iGetTitle());
idSet.insert(theID);
}
//
// Store histos according to their ID
//
vector vectorID(idSet.size());
copy(idSet.begin(), idSet.end(), vectorID.begin());
vector > vectorHistos(vectorID.size());
for (unsigned int i=0; iGetTitle());
if (theID == vectorID[i])
vectorHistos[i].push_back(allHistograms[j]);
}
}
cout << "Will merge: " << endl;
for (unsigned int i=0; i " << vectorHistos[i].size() << " tables" << endl;
//
// Merge
//
for (unsigned int i=0; iGetXaxis()->GetBinWidth(1);
double wy = refHisto->GetYaxis()->GetBinWidth(1);
double wz = refHisto->GetZaxis()->GetBinWidth(1);
bool notSameBinWidth = false;
for (unsigned int j=1; jGetXaxis()->GetBinWidth(1) - wx) > numeric_limits::epsilon() ||
fabs(vectorHistos[i][j]->GetYaxis()->GetBinWidth(1) - wy) > numeric_limits::epsilon() ||
fabs(vectorHistos[i][j]->GetZaxis()->GetBinWidth(1) - wz) > numeric_limits::epsilon()) {
cout << "Error: tables for id = " << vectorID[i] << " have different bin width. Cannot merge." << endl;
notSameBinWidth = true;
break;
}
}
if (notSameBinWidth) continue;
//
// Range of new histogram
//
double lowerEdgeX = refHisto->GetXaxis()->GetXmin();
double upperEdgeX = refHisto->GetXaxis()->GetXmax();
double lowerEdgeY = refHisto->GetYaxis()->GetXmin();
double upperEdgeY = refHisto->GetYaxis()->GetXmax();
double lowerEdgeZ = refHisto->GetZaxis()->GetXmin();
double upperEdgeZ = refHisto->GetZaxis()->GetXmax();
for (unsigned int j=1; jGetXaxis()->GetXmin(), lowerEdgeX);
upperEdgeX = max(vectorHistos[i][j]->GetXaxis()->GetXmax(), upperEdgeX);
lowerEdgeY = min(vectorHistos[i][j]->GetYaxis()->GetXmin(), lowerEdgeY);
upperEdgeY = max(vectorHistos[i][j]->GetYaxis()->GetXmax(), upperEdgeY);
lowerEdgeZ = min(vectorHistos[i][j]->GetZaxis()->GetXmin(), lowerEdgeZ);
upperEdgeZ = max(vectorHistos[i][j]->GetZaxis()->GetXmax(), upperEdgeZ);
}
//
// New binning
//
int nbinX = static_cast((upperEdgeX-lowerEdgeX)/wx + 0.5);
int nbinY = static_cast((upperEdgeY-lowerEdgeY)/wy + 0.5);
int nbinZ = static_cast((upperEdgeZ-lowerEdgeZ)/wz + 0.5);
//
// Create new histogram
//
TH3F *mergedHisto = new TH3F("table", refHisto->GetTitle(),
nbinX, lowerEdgeX, upperEdgeX,
nbinY, lowerEdgeY, upperEdgeY,
nbinZ, lowerEdgeZ, upperEdgeZ);
//
// Fill new histo
// Here we also make sure that bin centers between old and new match.
//
double x, y, z, xx, yy, zz, cellContent;
for (unsigned int k=0; kGetNbinsX();
int ny = vectorHistos[i][k]->GetNbinsY();
int nz = vectorHistos[i][k]->GetNbinsZ();
for (int ix = 1; ix <= nx; ix++) {
for (int iy = 1; iy <= ny; iy++) {
for (int iz = 1; iz <= nz; iz++) {
// old table
cellContent = vectorHistos[i][k]->GetBinContent(ix, iy, iz);
if (flags.ignoreZero && cellContent == 0) continue; // ignore 0 content
x = vectorHistos[i][k]->GetXaxis()->GetBinCenter(ix);
y = vectorHistos[i][k]->GetYaxis()->GetBinCenter(iy);
z = vectorHistos[i][k]->GetZaxis()->GetBinCenter(iz);
// new table
int globalBin = mergedHisto->FindBin(x, y, z);
if (globalBin == -1) {
cout << "Error while copying data from old to new table. Cannot find referring bin." << endl;
return 1;
}
int kx, ky, kz;
mergedHisto->GetBinXYZ(globalBin, kx, ky, kz);
xx = mergedHisto->GetXaxis()->GetBinCenter(kx);
yy = mergedHisto->GetYaxis()->GetBinCenter(ky);
zz = mergedHisto->GetZaxis()->GetBinCenter(kz);
if (fabs(xx-x) > numeric_limits::epsilon() ||
fabs(yy-y) > numeric_limits::epsilon() ||
fabs(zz-z) > numeric_limits::epsilon()) {
cout << "Error while copying data from old to new table. Bin center mismatch." << endl;
return 1;
}
mergedHisto->SetBinContent(kx, ky, kz, cellContent);
}
}
}
}
//
// At this point the merged histo is filled.
// Last check is to look for "holes" in the new table.
//
int nEmpty = 0;
int globalBin = 0;
vector missingBins;
for (int iz = 1; iz <= nbinZ; iz++)
for (int iy = 1; iy <= nbinY; iy++)
for (int ix = 1; ix <= nbinX; ix++) {
if (mergedHisto->GetBinContent(ix, iy, iz) == 0) {
nEmpty++;
if (flags.listMissingBins) missingBins.push_back(globalBin);
}
globalBin++;
}
if (nEmpty) {
cout << "Warning: the new merged table has " << nEmpty << " cells that are empty. This might indicate that\n";
cout << " there are not enough tables available to fully defines the full cubic range.\n";
cout << " Merged table (" << mergedHisto->GetTitle() << ") is not usable in Sartre." << endl;
if (flags.listMissingBins) {
cout << "List of missing/empty bins:" << endl;
for (unsigned int k=0; kGetTitle() << "_merged.root";
string filename = flags.prefix + filenameStream.str();
TFile *hfile = new TFile(filename.c_str(),"RECREATE");
if (!hfile) {
cout << "Error writing merged table. Cannot open file '" << filename.c_str() << "'." << endl;
}
else {
mergedHisto->Write();
hfile->Close();
cout << "Merged table written to file '" << filename.c_str() << "'." << endl;
delete mergedHisto;
mergedHisto = 0;
}
}
}
cout << "All 3D tables processed." << endl;
return 0;
}
int merge2DHistos(Flags &flags, vector& allHistograms)
{
set idSet;
for (unsigned i=0; iGetTitle());
idSet.insert(theID);
}
//
// Store histos according to their ID
//
vector vectorID(idSet.size());
copy(idSet.begin(), idSet.end(), vectorID.begin());
vector > vectorHistos(vectorID.size());
for (unsigned int i=0; iGetTitle());
if (theID == vectorID[i])
vectorHistos[i].push_back(allHistograms[j]);
}
}
cout << "Will merge: " << endl;
for (unsigned int i=0; i " << vectorHistos[i].size() << " tables" << endl;
cout << "---------------------------------------------------------------------" << endl;
//
// Merge
//
for (unsigned int i=0; iGetXaxis()->GetBinWidth(1);
double wy = refHisto->GetYaxis()->GetBinWidth(1);
bool notSameBinWidth = false;
for (unsigned int j=1; jGetXaxis()->GetBinWidth(1) - wx) > numeric_limits::epsilon() ||
fabs(vectorHistos[i][j]->GetYaxis()->GetBinWidth(1) - wy) > numeric_limits::epsilon() ) {
cout << "Error: tables for id = " << vectorID[i] << " have different bin width. Cannot merge." << endl;
notSameBinWidth = true;
break;
}
}
if (notSameBinWidth) continue;
//
// Range of new histogram
//
double lowerEdgeX = refHisto->GetXaxis()->GetXmin();
double upperEdgeX = refHisto->GetXaxis()->GetXmax();
double lowerEdgeY = refHisto->GetYaxis()->GetXmin();
double upperEdgeY = refHisto->GetYaxis()->GetXmax();
for (unsigned int j=1; jGetXaxis()->GetXmin(), lowerEdgeX);
upperEdgeX = max(vectorHistos[i][j]->GetXaxis()->GetXmax(), upperEdgeX);
lowerEdgeY = min(vectorHistos[i][j]->GetYaxis()->GetXmin(), lowerEdgeY);
upperEdgeY = max(vectorHistos[i][j]->GetYaxis()->GetXmax(), upperEdgeY);
}
//
// New binning
//
int nbinX = static_cast((upperEdgeX-lowerEdgeX)/wx + 0.5);
int nbinY = static_cast((upperEdgeY-lowerEdgeY)/wy + 0.5);
//
// Create new histogram
//
TH2F *mergedHisto = new TH2F("table", refHisto->GetTitle(),
nbinX, lowerEdgeX, upperEdgeX,
nbinY, lowerEdgeY, upperEdgeY);
//
// Fill new histo
// Here we also make sure that bin centers between old and new match.
//
double x, y, xx, yy, cellContent;
for (unsigned int k=0; kGetNbinsX();
int ny = vectorHistos[i][k]->GetNbinsY();
for (int ix = 1; ix <= nx; ix++) {
for (int iy = 1; iy <= ny; iy++) {
// old table
cellContent = vectorHistos[i][k]->GetBinContent(ix, iy);
if (flags.ignoreZero && cellContent == 0) continue; // ignore 0 content
x = vectorHistos[i][k]->GetXaxis()->GetBinCenter(ix);
y = vectorHistos[i][k]->GetYaxis()->GetBinCenter(iy);
// new table
int globalBin = mergedHisto->FindBin(x, y);
if (globalBin == -1) {
cout << "Error while copying data from old to new table. Cannot find referring bin." << endl;
return 1;
}
int kx, ky, kdummy;
mergedHisto->GetBinXYZ(globalBin, kx, ky, kdummy);
xx = mergedHisto->GetXaxis()->GetBinCenter(kx);
yy = mergedHisto->GetYaxis()->GetBinCenter(ky);
if (fabs(xx-x) > numeric_limits::epsilon() ||
fabs(yy-y) > numeric_limits::epsilon() ) {
cout << "Error while copying data from old to new table. Bin center mismatch." << endl;
return 1;
}
mergedHisto->SetBinContent(kx, ky, cellContent);
}
}
}
//
// At this point the merged histo is filled.
// Last check is to look for "holes" in the new table.
//
int nEmpty = 0;
int globalBin = 0;
vector missingBins;
for (int iy = 1; iy <= nbinY; iy++) {
for (int ix = 1; ix <= nbinX; ix++) {
if (mergedHisto->GetBinContent(ix, iy) == 0) {
nEmpty++;
if (flags.listMissingBins) missingBins.push_back(globalBin);
}
globalBin++;
}
}
if (nEmpty) {
cout << "Warning: the new merged table has " << nEmpty << " cells that are empty. This might indicate that\n";
cout << " there are not enough tables available to fully defines the full square range.\n";
cout << " Merged table (" << mergedHisto->GetTitle() << ") is not usable in Sartre." << endl;
if (flags.listMissingBins) {
cout << "List of missing/empty bins:" << endl;
for (unsigned int k=0; kGetTitle() << "_merged.root";
string filename = flags.prefix + filenameStream.str();
TFile *hfile = new TFile(filename.c_str(),"RECREATE");
if (!hfile) {
cout << "Error writing merged table. Cannot open file '" << filename.c_str() << "'." << endl;
}
else {
mergedHisto->Write();
hfile->Close();
cout << "Merged table written to file '" << filename.c_str() << "'." << endl;
delete mergedHisto;
mergedHisto = 0;
}
}
}
cout << "All 2D tables processed." << endl;
return 0;
}
int main(int argc, char **argv)
{
Flags flags;
//
// Handle command line arguments
//
if (argc == 1) {
usage(argv[0]);
return 2;
}
flags.listMissingBins = false;
flags.noOutput = false;
flags.ignoreZero = false;
int ch;
while ((ch = getopt(argc, argv, "m0np:")) != -1) {
switch (ch) {
case 'm':
flags.listMissingBins = true;
break;
case '0':
flags.ignoreZero = true;
break;
case 'p':
flags.prefix = string(optarg);
break;
case 'n':
flags.noOutput = true;
break;
case '?':
default:
usage(argv[0]);
return 2;
break;
}
}
if (optind == argc) {
usage(argv[0]);
return 2;
}
vector allFiles;
for (int index = optind; index < argc; index++)
allFiles.push_back(string(argv[index]));
//
// Open the files and get all tables/histograms to be merged
//
vector all3DHistograms;
vector all2DHistograms;
createHistogramLists(allFiles, all3DHistograms, all2DHistograms);
if (all3DHistograms.size() && all3DHistograms.size()) {
cout << "The given files contain 2D and 3D tables." << endl;
cout << "Will do the 3D tables first then the 2D tables." << endl;
}
if (all3DHistograms.size()) {
cout << "Starting process for 3D histograms (standard tables)" << endl;
int rc = merge3DHistos(flags, all3DHistograms);
if (rc != 0) {
cout << "Errors during merging of 3D histos." << endl;
return rc;
}
}
if (all2DHistograms.size()) {
cout << "Starting process for 2D histograms (UPC tables)" << endl;
int rc = merge2DHistos(flags, all2DHistograms);
if (rc != 0) {
cout << "Errors during merging of 2D histos." << endl;
return rc;
}
}
return 0;
}
Index: trunk/src/WaveOverlap.cpp
===================================================================
--- trunk/src/WaveOverlap.cpp (revision 369)
+++ trunk/src/WaveOverlap.cpp (revision 370)
@@ -1,274 +1,274 @@
//==============================================================================
// WaveOverlap.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "WaveOverlap.h"
#include "Constants.h"
#include "TableGeneratorSettings.h"
#include "DipoleModelParameters.h"
#include "TMath.h"
#include
#include
#include "TF1.h"
#include "Math/Functor.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
WaveOverlap::WaveOverlap()
{
mParameters = new DipoleModelParameters(TableGeneratorSettings::instance());
}
WaveOverlap::~WaveOverlap() {/* no op*/}
//
// VECTOR MESONS
//
WaveOverlapVM::WaveOverlapVM()
{
mNT = mRT2 = 0;
mMf = 0;
mMf2 = 0;
mEf = 0;
mMV = 0;
mNL = mRL2 = 0;
mBoostedGaussianMf = 0;
}
double WaveOverlapVM::uiNormT(const double* var){
double z=var[0];
double r=var[1];
double phi=transverseWaveFunction(r, z); //GeV0
double dphidr=dDrTransverseWaveFunction(r, z); //GeV1
return 3.*r/hbarc/(z*z*(1-z)*(1-z))
*(mMf2*phi*phi + (z*z + (1-z)*(1-z))*dphidr*dphidr); //GeV1
}
double WaveOverlapVM::uiNormL(const double* var){
double z=var[0];
double r=var[1];
double phi=longitudinalWaveFunction(r, z); //GeV0
double d2phidr2=laplaceRLongitudinalWaveFunction(r, z); //GeV2
double term=mMV*phi + (mMf2*phi - d2phidr2)/(mMV*z*(1-z)); //GeV1
return 3.*r/hbarc*term*term; //GeV1
}
double WaveOverlapVM::uiDecayWidth(double* var, double*)
{
double z=*var;
double phi=longitudinalWaveFunction(0., z); //GeV0
double d2phidr2=laplaceRLongitudinalWaveFunction(0., z); //GeV0
double result = mEf*3./M_PI*(mMV*phi+(mMf2*phi-d2phidr2)/(mMV*z*(1-z))); //GeV1
return result;
}
void WaveOverlapVM::testBoostedGaussianParameters(int id)
{
//
// This function calculates the resulting normalisation
// and decay width resulting from the boosted gaussian parameters
// and compare with the actual values. This can be used to
// test or modify the boosted gaussian parameters.
//
// KMW paper hep-ph/0606272 Eqs. (24)-(28)
//
// Start with decay width:
TF1 fDW("fDW", this, &WaveOverlapVM::uiDecayWidth, 0., 1., 0.);
ROOT::Math::WrappedTF1 wfDW(fDW);
ROOT::Math::GaussIntegrator giDW;
giDW.SetFunction(wfDW);
giDW.SetAbsTolerance(0.);
giDW.SetRelTolerance(1e-5);
double f_VL=giDW.Integral(0., 1.); //GeV
PR(f_VL);
cout<<"The e+e- decay width is: "<ee)="<<4*M_PI*alpha_em*alpha_em*f_VL*f_VL/3/mMV*1e6
<<" keV"<ee)=1.340 +- 0.018 keV"<ee)=5.55 +- 0.14 +- 0.02 keV"<ee)=1.27 +- 0.04 keV"<ee)=7.04 +- 0.06 keV"<boostedGaussianR2(val);
mRT2 = mRL2;
mNT = mParameters->boostedGaussianNT(val);
mNL = mParameters->boostedGaussianNL(val);
mBoostedGaussianMf = mParameters->boostedGaussianQuarkMass(val);
mBoostedGaussianMf2 = mBoostedGaussianMf*mBoostedGaussianMf;
}
void WaveOverlapVM::setProcess(int val)
{
switch (val) {
case 113:
mMf = mParameters->quarkMass(1);
mMV = 0.776;
mEf = 1./sqrt(2.);
break;
case 333:
mMf = mParameters->quarkMass(2);
mMV = 1.019;
mEf = 1./3.;
break;
case 443:
mMf = mParameters->quarkMass(3);
mMV = 3.096916;
mEf = 2./3.;
break;
case 553:
mMf = mParameters->quarkMass(4);
mMV = 9.46;
mEf = -1./3.;
break;
default:
cerr << "WaveOverlap::setProcess(): error no such type: " << val << endl;
break;
}
mMf2 = mMf*mMf;
}
//
// DVCS
//
double WaveOverlapDVCS::T(double z, double Q2, double r)
{
// KMW paper hep-ph/0606272 Eq. 17
double term0, term1, term2;
double result=0;
for (int iFlav=0; iFlav<4; iFlav++) {
double mf=mParameters->quarkMass(iFlav);
double ef=quarkCharge[iFlav];
double eps2 = z*(1-z)*Q2 + mf*mf;
double eps = sqrt(eps2);
term0=2.*Nc/M_PI*alpha_em*ef*ef;
term1=( z*z+(1-z)*(1-z) )*
eps*TMath::BesselK1(eps*r/hbarc)*mf*TMath::BesselK1(mf*r/hbarc);
term2=mf*mf*TMath::BesselK0(eps*r/hbarc)*TMath::BesselK0(mf*r/hbarc);
result += term0*(term1+term2);
}
return result;
}
double WaveOverlapDVCS::L(double, double, double) {return 0;}
Index: trunk/src/Kinematics.h
===================================================================
--- trunk/src/Kinematics.h (revision 369)
+++ trunk/src/Kinematics.h (revision 370)
@@ -1,105 +1,105 @@
//==============================================================================
// Kinematics.h
//
-// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//
//------------------------------------------------------------------------------
//
// Note that Kinematics is only meant for ep and eA diffactive events.
// It is not applicable for UPC kinematics. Use UpcKinematics instead.
//
//==============================================================================
#ifndef Kinematics_h
#define Kinematics_h
#include "Constants.h"
#include "Settings.h"
#include "TLorentzVector.h"
#include "Math/IFunction.h"
class Kinematics {
public:
static double xprobe(double Q2, double W2, double vmM);
static TLorentzVector electronBeam(double eE, bool upc = false);
static TLorentzVector hadronBeam(double eH);
static double s(const TLorentzVector& e, const TLorentzVector& p);
static double s(double eE, double hE, bool upc = false);
static double x(double Q2, double W2);
static double y(double Q2, double x, double s);
static double ymin(double s, double vmM);
static double ymax(double s, double vmM);
static double W2(double Q2, double x);
static double W2(double Q2, double xprobe, double vmM);
static double W2min(double vmM);
static double W2max(double s);
static double Q2min(double y);
static double Q2max(double s);
static double xpomeron(double t, double Q2, double W2, double vmM);
static double tmax(double xpom); // often referred to as tmin implying min |T|
static double tmax(double t, double Q2, double W2, double vmM);
static double tmin(double hE); // the smallest t value possible
static double Egamma(double xpom, double t, double vmM, double hBeamEnergy, double eBeamEnergy, double MY2minusM2 = 0);
static bool validUPC(double hBeamEnergy, double eBeamEnergy, double t,
double xpom, double vmMass, bool verbose = false); // UPC only
static bool valid(double s, double t, double Q2, double W2, double vmMass,
bool useTrueXp = false,
bool verbose = false);
static bool error();
static double xpomMin(double massVM, double t, TLorentzVector hBeam, TLorentzVector eBeam, double MY2minusM2 = 0); // UPC only
private:
static double xpomMin(double s, double vmM, double t); // UPC only
private:
static bool mError;
static double mXpomMin;
static double mTold;
static bool mXpomMinIsEvaluated;
};
//-------------------------------------------------------------------------------
//
// Helper class needed to find root in
// Kinematics::validUPC()
//
//-------------------------------------------------------------------------------
class LowerXpomeronFormula : public ROOT::Math::IBaseFunctionOneDim
{
public:
double DoEval(double) const;
ROOT::Math::IBaseFunctionOneDim* Clone() const;
void calculateValidRange(double&, double&);
public:
double mT;
double mVmMass2;
double mXpomMin;
double mMY2;
TLorentzVector mElectronBeam;
TLorentzVector mHadronBeam;
double mMY2minusM2;
};
#endif
Index: trunk/src/Sartre.cpp
===================================================================
--- trunk/src/Sartre.cpp (revision 369)
+++ trunk/src/Sartre.cpp (revision 370)
@@ -1,1012 +1,1012 @@
//==============================================================================
// Sartre.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Note:
// When not using runcards, the user must first create an instance of Sartre
// then get the settings via one of:
// Sartre::runSettings()
// EventGeneratorSettings::instance()
// Once init() is called settings cannot be changed any more.
//==============================================================================
#include "Version.h"
#include "Kinematics.h"
#include "Sartre.h"
#include "ModeFinderFunctor.h"
#include "Math/BrentMinimizer1D.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/GSLMinimizer.h"
#include "Math/Functor.h"
#include "TUnuranMultiContDist.h"
#include
#include
#include
#include "TH2D.h"
#include "TFile.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
Sartre::Sartre()
{
mIsInitialized = false;
mCurrentEvent = 0;
mNucleus = 0;
mUpcNucleus= 0;
mSettings = 0;
mPDF_Functor = 0;
mPDF = 0;
mEventCounter = 0;
mTriesCounter = 0;
mTotalCrossSection = 0;
mCrossSection = 0;
mTableCollection = 0;
mProtonTableCollection = 0;
mUnuran = 0;
mEvents = 0;
mTries = 0;
mS = 0;
mA = 0;
}
Sartre::~Sartre()
{
delete mNucleus;
delete mUpcNucleus;
delete mPDF_Functor;
delete mPDF;
delete mCrossSection;
if (mTableCollection != mProtonTableCollection) {
delete mTableCollection;
delete mProtonTableCollection;
}
else
delete mTableCollection;
delete mUnuran;
}
bool Sartre::init(const char* runcard)
{
mStartTime = time(0);
bool ok;
//
// Reset member variables.
// Note that one instance of Sartre should be able to get
// initialized multiple times.
//
mEvents = 0;
mTries = 0;
mTotalCrossSection = 0;
//
// Print header
//
string ctstr(ctime(&mStartTime));
ctstr.erase(ctstr.size()-1, 1);
cout << "/========================================================================\\" << endl;
cout << "| |" << endl;
cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
cout << "| |" << endl;
cout << "| An event generator for exclusive diffractive vector meson production |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
- cout << "| Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich |" << endl;
+ cout << "| Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich |" << endl;
cout << "| |" << endl;
cout << "| This program is free software: you can redistribute it and/or modify |" << endl;
cout << "| it under the terms of the GNU General Public License as published by |" << endl;
cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl;
cout << "| any later version. |" << endl;
cout << "| |" << endl;
cout << "| Code compiled on " << setw(12) << left << __DATE__;
cout << setw(41) << left << __TIME__ << right << '|' << endl;
cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
cout << "\\========================================================================/" << endl;
mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton
//
// Read runcard if available
//
if (runcard) {
if (!mSettings->readSettingsFromFile(runcard)) {
cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl;
exit(1);
}
else
cout << "Runcard is '" << runcard << "'." << endl;
}
else
cout << "No runcard provided." << endl;
//
// Set beam particles and center of mass energy
// Note, that if we are in UPC mode the electron
// is actually treated as the hadron beam, i.e.
// the mass is the proton mass.
//
mElectronBeam = mSettings->eBeam();
mHadronBeam = mSettings->hBeam();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mA = mSettings->A();
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (allowBreakup) {
if (!getenv("SARTRE_DIR")) {
cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n"
"to locate tables needed for the generation if nuclear breakups." << endl;
exit(1);
}
}
if (mNucleus) delete mNucleus;
mNucleus = new FrangibleNucleus(mA, allowBreakup);
if (mSettings->UPC()) {
if (mUpcNucleus) delete mUpcNucleus;
mUpcNucleus = new FrangibleNucleus(mSettings->UPCA(), allowBreakup);
}
string upcNucleusName;
if (mSettings->UPC()) {
upcNucleusName = mUpcNucleus->name();
cout << "Sartre is running in UPC mode" << endl;
cout << "Hadron 1 beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron 1 beam: " << mHadronBeam << endl;
cout << "Hadron 2 beam species: " << upcNucleusName << " (" << mSettings->UPCA() << ")" << endl;
cout << "Hadron 2 beam: " << mElectronBeam << endl;
}
else {
cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron beam: " << mHadronBeam << endl;
cout << "Electron beam: " << mElectronBeam << endl;
}
//
// Get details about the processes and models
//
mDipoleModelType = mSettings->dipoleModelType();
mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
mVmID = mSettings->vectorMesonId();
cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
cout << "Process is ";
if (mSettings->UPC()) {
cout << mNucleus->name() << " + " << upcNucleusName << " -> "
<< mNucleus->name() << "' + " << upcNucleusName << "' + ";
}
else {
if (mA > 1)
cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + ";
else
cout << "e + p -> e' + p' + ";
}
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID);
cout << vectorMesonPDG->GetName() << endl;
//
// Print-out seed for reference
//
cout << "Random generator seed: " << mSettings->seed() << endl;
//
// Load in the tables containing the amplitude moments
//
if (!getenv("SARTRE_DIR")) {
cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
exit(1);
}
if (mTableCollection) delete mTableCollection;
mTableCollection = new TableCollection;
ok = mTableCollection->init(mA, mDipoleModelType, mDipoleModelParameterSet, mVmID);
if (!ok) {
cout << "Error, could not initialize lookup tables for requested process." << endl;
return false;
}
//
// Load in the p tables for the lambda lookup tables (or to calculate lambda if not available)
//
if (mSettings->correctForRealAmplitude() || mSettings->correctSkewedness()) {
if (mA == 1) {
mProtonTableCollection = mTableCollection;
}
else {
if (mProtonTableCollection) delete mProtonTableCollection;
mProtonTableCollection = new TableCollection;
ok = mProtonTableCollection->init(1, mDipoleModelType, mDipoleModelParameterSet, mVmID);
if (!ok) {
cout << "Error: could not initialize proton lookup tables for requested process." << endl;
cout << " These tables are needed for skewedness and real amplitude corrections." << endl;
return false;
}
}
}
else
mProtonTableCollection = 0;
//
// Kinematic limits and generator range
//
// There are 3 ranges we have to deal with
// 1. the kinematic range requested by the user
// if given.
// The user can only control Q2 and W but not t.
// For UPC that's xpom.
// 2. the range of the table(s)
// 3. the kinematically allowed range
//
// Of course (3) is more complex than a simple cube/square.
// However, we deal with the detailed shape of the kinematic
// range later using Kinematics::valid() when we generate the
// individual events.
// For setting up UNU.RAN we have to get the cubic/square
// envelope that satifies (1)-(3).
// Note, that they are correlated which makes the order
// in which we do things a bit tricky.
//
//
// Step 1:
// Set the limits to that of the table(s).
// Note, the indices 0-2 refer to t, Q2, and W2.
// For UPC we have only t and xpom.
//
if (mSettings->UPC()) {
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minX(); mUpperLimit[1] = mTableCollection->maxX();
mLowerLimit[2] = mUpperLimit[2] = 0;
}
else {
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minQ2(); mUpperLimit[1] = mTableCollection->maxQ2();
mLowerLimit[2] = mTableCollection->minW2(); mUpperLimit[2] = mTableCollection->maxW2();
}
//
// Step 2:
// Kinematic limits might overrule boundaries from step 1
//
if (mSettings->UPC()) {
double kineXpomMin = Kinematics::xpomMin(vectorMesonPDG->Mass(), mUpperLimit[0], mHadronBeam, mElectronBeam);
double kineXpomMax = 1;
mLowerLimit[1] = max(kineXpomMin, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineXpomMax);
double kineTmax = Kinematics::tmax(kineXpomMin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
}
else {
// double kineYmax = Kinematics::ymax(mS, vectorMesonPDG->Mass());
double kineYmin = Kinematics::ymin(mS, vectorMesonPDG->Mass());
double kineQ2min = Kinematics::Q2min(kineYmin);
double kineQ2max = Kinematics::Q2max(mS);
double kineW2min = Kinematics::W2min(vectorMesonPDG->Mass());
double kineW2max = Kinematics::W2max(mS);
kineQ2min = max(kineQ2min, mLowerLimit[1]);
kineQ2max = min(kineQ2max, mUpperLimit[1]);
kineW2min = max(kineW2min, mLowerLimit[2]);
kineW2max = min(kineW2max, mUpperLimit[2]);
double kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // first arg (t) set to 0 (recursive)
double kineTmax = Kinematics::tmax(kineXPmin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
mLowerLimit[1] = max(kineQ2min, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineQ2max);
mLowerLimit[2] = max(kineW2min, mLowerLimit[2]); mUpperLimit[2] = min(mUpperLimit[2], kineW2max);
}
//
// Step 3:
// Deal with user provided limits.
// User settings are ignored (switched off) if min >= max.
//
if (mSettings->UPC()) {
if (mSettings->xpomMin() < mSettings->xpomMax()) {
if (mSettings->xpomMin() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of xpomeron (" << mSettings->xpomMin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->xpomMin();
}
if (mSettings->xpomMax() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of xpomeron (" << mSettings->xpomMax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->xpomMax();
}
}
}
else {
if (mSettings->W2min() < mSettings->W2max()) { // W2 first
if (mSettings->W2min() < mLowerLimit[2]) {
cout << "Warning, requested lower limit of W (" << mSettings->Wmin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << sqrt(mLowerLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[2] = mSettings->W2min();
}
if (mSettings->W2max() > mUpperLimit[2]) {
cout << "Warning, requested upper limit of W (" << mSettings->Wmax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << sqrt(mUpperLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[2] = mSettings->W2max();
}
}
if (mSettings->Q2min() < mSettings->Q2max()) { // Q2
if (mSettings->Q2min() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of Q2 (" << mSettings->Q2min() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->Q2min();
// ????? kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // new Q2min changes tmax
// ????? mUpperLimit[0] = min(mUpperLimit[0], Kinematics::tmax(kineXPmin));
}
if (mSettings->Q2max() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of Q2 (" << mSettings->Q2max() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->Q2max();
}
}
}
//
// Check if any phase space is left
//
if (mLowerLimit[0] >= mUpperLimit[0]) {
cout << "Invalid range in t: t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
exit(1);
}
if (mLowerLimit[1] >= mUpperLimit[1]) {
if (mSettings->UPC())
cout << "Invalid range in xpomeron: xpomeron=[";
else
cout << "Invalid range in Q2: Q2=[";
cout << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl;
exit(1);
}
if (!mSettings->UPC() && mLowerLimit[2] >= mUpperLimit[2]) {
cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl;
exit(1);
}
//
// Print-out limits (all verbose levels)
//
if (mSettings->verbose()) {
cout << "Kinematic limits used for event generation:" << endl;
if (mSettings->UPC()) {
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "xpom=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
}
else {
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl;
}
}
//
// Check if the lambda table covers the kinematic range
// or if we have to calculate lambda from the p mean_A table
// directly.
// This is just to inform the user, no action is taken.
// Here we assume that transverse and longitudinal tables
// always have the same range and vheck only one.
//
// If even the p mean_A table is not big enough we switch
// corrections off and inform the user.
//
if (mProtonTableCollection) {
bool lambdaTableExists, tooSmallLambda, tooSmallMean;
if (mSettings->UPC()) {
lambdaTableExists = mProtonTableCollection->tableExists(lambda_A);
tooSmallLambda = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A);
tooSmallMean = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A);
}
else {
lambdaTableExists = mProtonTableCollection->tableExists(transverse, lambda_A);
tooSmallLambda = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A);
tooSmallMean = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A);
}
if (!lambdaTableExists || tooSmallLambda) {
if (tooSmallMean) {
cout << "Warning: the kinematic coverage of the table containing the lambda values " << endl;
cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
cout << " then the required range (or does not exist). The proton amplitude " << endl;
cout << " tables that could be used to calculate them on the fly are too" << endl;
cout << " small as well. Corrections will be switched off." << endl;
mSettings->setCorrectForRealAmplitude(false);
mSettings->setCorrectSkewedness(false);
delete mProtonTableCollection;
mProtonTableCollection = 0;
}
else {
if (lambdaTableExists) {
- cout << "Warning: the kinematic coverage of the table containing the lambda values " << endl;
- cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
- cout << " then the required range. The missing lambda values will be " << endl;
- cout << " calculated using the proton amplitude tables directly." << endl;
+ cout << "Info: the kinematic coverage of the table containing the lambda values " << endl;
+ cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
+ cout << " then the required range. The missing lambda values will be " << endl;
+ cout << " calculated using the proton amplitude tables directly." << endl;
}
else {
- cout << "Warning: the table containing the lambda values needed for skewedness and/or real " << endl;
- cout << " amplitude corrections is not available. The missing lambda values will be" << endl;
- cout << " calculated using the referring proton amplitude tables directly." << endl;
+ cout << "Info: the table containing the lambda values needed for skewedness and/or real " << endl;
+ cout << " amplitude corrections is not available. The missing lambda values will be" << endl;
+ cout << " calculated using the referring proton amplitude tables directly." << endl;
}
}
}
}
//
// Setup cross-section functor
// It is this functor that is used by all other functors,
// functions, and wrappers when dealing with cross-sections.
//
if (mCrossSection) delete mCrossSection;
mCrossSection = new CrossSection;
mCrossSection->setTableCollection(mTableCollection);
if (mProtonTableCollection) mCrossSection->setProtonTableCollection(mProtonTableCollection);
//
// UNU.RAN needs the domain (boundaries) and the mode.
// The domain is already defined, here we find the mode, which is tricky.
// The max. cross-section is clearly at the domain boundary in Q2=Q2min.
// The position in W2 and t is not obvious. It sits along a line given
// by tmax(W2). The approach here is to use the BrentMinimizer1D that
// performs first a scan a then a Brent fit.
//
double theMode[3];
if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
if (mSettings->UPC()) {
//
// For now we keep the option to create a histogram of the
// generated cross-sections(t, log(xpom)) for QA purposes.
//
if (mSettings->verbose() && mSettings->verboseLevel() > 9) {
cout << "Creating 2D histogram of cross-section as fct of t and log(xpom)." << endl;
cout << "Be patient this might take a while. Histo stored in file 'landscape.root'" << endl;
TFile file("landscape.root", "RECREATE");
int nbins_t = 300;
int nbins_logx = 300;
TH2D histo("histo", "mode finder studies",
nbins_t, mLowerLimit[0], mUpperLimit[0], // t
nbins_logx, log(mLowerLimit[1]), log(mUpperLimit[1])); // log(xpom)
for (int ix = 1; ix <= nbins_logx; ix++) {
for (int it = 1; it <= nbins_t; it++) {
double logx = histo.GetYaxis()->GetBinCenter(ix);
double x = exp(logx);
double t = histo.GetXaxis()->GetBinCenter(it);
if (!Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
t, x, vectorMesonPDG->Mass(), false)) {
histo.SetBinContent(it, ix, -5.e9);
}
else {
histo.SetBinContent(it, ix, (*mCrossSection)(t, x));
}
}
}
histo.Write();
file.Close();
cout << "Histogram written to file. All done." << endl;
}
//
// Create functor
// Mode is at the max t available.
// Note that the functop works with log(xpom).
//
theMode[0] = mUpperLimit[0]; // t
theMode[1] = mUpperLimit[1]; // xpom for UPC
theMode[2] = 0; // not used for UPC
UPCModeFinderFunctor modeFunctor(mCrossSection, vectorMesonPDG->Mass(),
mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy());
ROOT::Math::BrentMinimizer1D minimizer;
minimizer.SetFunction(modeFunctor, log(mLowerLimit[1]), log(mUpperLimit[1]));
minimizer.SetNpx(100000);
ok = minimizer.Minimize(1000000, 1.e-8, 1.e-10);
if (! ok) {
cout << "Error, failed to find mode of pdf." << endl;
return false;
}
//
// Get the result
//
theMode[1] = exp(minimizer.XMinimum()); // xpom
theMode[0] = Kinematics::tmax(theMode[1]);
double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1]);
if (mSettings->verbose()) {
cout << "\tlocation: t=" << theMode[0] << ", xpom=" << theMode[1]
<< "; value: " << crossSectionAtMode << endl;
}
//
// Consistency check of mode
//
if (crossSectionAtMode <= 0) {
cout << "Error: cross-section at mode value is invalid." << endl;
return false;
}
if (!Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
theMode[0], theMode[1],
vectorMesonPDG->Mass(), true)) {
cout << "Error: mode of pdf is outside kinematic limits." << endl;
return false;
}
}
else {
theMode[0] = mUpperLimit[0]; // t
theMode[1] = mLowerLimit[1]; // Q2 (xpom for UPC)
theMode[2] = mLowerLimit[2]; // W2 (dummy for UPC)
ModeFinderFunctor modeFunctor(mCrossSection, theMode[1], vectorMesonPDG->Mass(), mLowerLimit[0], mUpperLimit[0]);
ROOT::Math::BrentMinimizer1D minimizer;
minimizer.SetFunction(modeFunctor, mLowerLimit[2], mUpperLimit[2]);
minimizer.SetNpx(static_cast(mUpperLimit[2]-mLowerLimit[2]));
ok = minimizer.Minimize(100000, 0, 1.e-8);
if (! ok) {
cout << "Error, failed to find mode of pdf." << endl;
exit(1);
}
theMode[2] = minimizer.XMinimum(); // W2
theMode[0] = Kinematics::tmax(0, theMode[1], theMode[2], vectorMesonPDG->Mass()); // first arg (t) must be 0 here
if (theMode[0] > mUpperLimit[0]) theMode[0] = mUpperLimit[0];
double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1], theMode[2]);
if (mSettings->verbose()) {
cout << "\tlocation: t=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]);
cout << "; value: " << crossSectionAtMode << endl;
}
}
//
// Initialize 2D (UPC) or 3D random generator
//
// Test show that UNU.RAN runs smoother in log(Q2)
// and log(cross-section). Functor CrossSection has
// a spezialized method for UNU.RAN, unuranPDF().
//
// In UPC mode the mPDF_Functor is using a different
// method and is only 2D since Q2=0
//
// domain and mode for Q2 -> log(Q2) or xpom -> log(xpom)
mLowerLimit[1] = log(mLowerLimit[1]);
mUpperLimit[1] = log(mUpperLimit[1]);
theMode[1] = log(theMode[1]);
if (mPDF_Functor) delete mPDF_Functor;
if (mPDF) delete mPDF;
if (mSettings->UPC())
mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 2);
else
mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 3);
mPDF = new TUnuranMultiContDist(*mPDF_Functor, true); // last arg = pdf in log or not
mPDF->SetDomain(mLowerLimit, mUpperLimit);
mPDF->SetMode(theMode);
if (mUnuran) delete mUnuran;
mUnuran = new TUnuran;
mCrossSection->setCheckKinematics(false); // avoid numeric glitch in Init()
mUnuran->Init(*mPDF, "method=hitro");
mCrossSection->setCheckKinematics(true);
mUnuran->SetSeed(mSettings->seed());
//
// Burn in generator
//
double xrandom[3];
for (int i=0; i<1000; i++) {
mUnuran->SampleMulti(xrandom);
}
mEventCounter = 0;
mTriesCounter = 0;
mIsInitialized = true;
cout << "Sartre is initialized." << endl << endl;
return true;
}
bool Sartre::init(const string& str) // overloaded version of init()
{
if (str.empty())
return init();
else
return init(str.c_str());
}
vector > Sartre::kinematicLimits()
{
vector > array;
array.push_back(make_pair(mLowerLimit[0], mUpperLimit[0])); // t
array.push_back(make_pair(exp(mLowerLimit[1]), exp(mUpperLimit[1]))); // Q2 or xpom
if (!mSettings->UPC())
array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W
return array;
}
Event* Sartre::generateEvent()
{
if (!mIsInitialized) {
cout << "Sartre::generateEvent(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return 0;
}
double xrandom[3];
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID);
double vmMass = vectorMesonPDG->Mass();
//
// Generate one event
//
while (true) {
mTriesCounter++;
delete mCurrentEvent;
mCurrentEvent = new Event;
//
// Get t, Q2, W2 from TUnuran and check for kinematics.
// Q2 is generated as log(Q2) so we transform it back first.
// This is the only place where Kinematics::valid() is called
// with the fully correct xpomeron calculation switched on.
//
mUnuran->SampleMulti(xrandom);
xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2 or log(xpom) -> xpom
bool isValidEvent;
if (mSettings->UPC()) {
isValidEvent = Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
xrandom[0], xrandom[1],
vectorMesonPDG->Mass(), (mSettings->verboseLevel() > 1));
}
else {
isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2],
vmMass, true, (mSettings->verboseLevel() > 1));
}
if (!isValidEvent) {
if (mSettings->verboseLevel() > 2)
cout << "Sartre::generateEvent(): event rejected, not within kinematic limits" << endl;
continue;
}
//
// Fill beam particles in Event structure
// Kinematics for eA is reported as 'per nucleon'
//
if (mSettings->UPC()) {
//
// For UPC some of the event variables
// are filled in the final state generator.
// They are set to 0 here.
//
mCurrentEvent->eventNumber = mEventCounter;
mCurrentEvent->t = xrandom[0]; // t
mCurrentEvent->Q2 = 0;
mCurrentEvent->x = 0;
mCurrentEvent->y = 0;
mCurrentEvent->s = mS; // s
mCurrentEvent->W = 0;
mCurrentEvent->beta = 1;
mCurrentEvent->xpom = xrandom[1];
mCurrentEvent->polarization = GammaPolarization::transverse;
mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = mUpcNucleus->pdgID(); // misuse ebeam spot for this
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
}
else {
mCurrentEvent->eventNumber = mEventCounter;
mCurrentEvent->t = xrandom[0]; // t
mCurrentEvent->Q2 = xrandom[1]; // Q2
mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
mCurrentEvent->s = mS; // s
mCurrentEvent->W = sqrt(xrandom[2]);
mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = 11; // e-
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
}
//
// Generate the final state particles
//
bool ok;
if (mSettings->UPC()) {
// also fills some of the undefined event variables
ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->xpom,
(mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent);
}
else {
ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->y, mCurrentEvent->Q2,
(mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent);
}
if (!ok) {
if (mSettings->verboseLevel() > 1) cout << "Sartre::generateEvent(): failed to generate final state" << endl;
continue;
}
break;
}
mEventCounter++;
//
// Nuclear breakup
//
// If the event is incoherent the final state generator does produce a
// 'virtual' proton with m > m_p which is used in Nucleus to calculate
// the excitation energy and the boost.
//
int indexOfScatteredHadron = 6;
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (mNucleus) mNucleus->resetBreakup(); // clear previous event in any case
if (allowBreakup && mCurrentEvent->diffractiveMode == incoherent && mNucleus) {
int nFragments = mNucleus->breakup(mCurrentEvent->particles[indexOfScatteredHadron].p);
//
// Merge the list of products into the event list.
// We loose some information here. The user can always go back to
// the nucleus and check the decay products for more details.
// In the original list the energy is per nuclei, here we transform it
// to per nucleon to stay consistent with Sartre conventions.
//
const vector& products = mNucleus->breakupProducts();
for (int i=0; iparticles.size();
fragment.pdgId = products[i].pdgId;
fragment.status = 1;
fragment.p = products[i].p*(1/static_cast(products[i].A));
fragment.parents.push_back(indexOfScatteredHadron);
mCurrentEvent->particles.push_back(fragment);
}
}
//
// Complete event record
//
if (!mSettings->UPC()) {
mCurrentEvent->xpom = Kinematics::xpomeron(mCurrentEvent->t, mCurrentEvent->Q2, mCurrentEvent->W*mCurrentEvent->W, vmMass);
mCurrentEvent->beta = mCurrentEvent->x/mCurrentEvent->xpom;
}
return mCurrentEvent;
}
double Sartre::totalCrossSection()
{
if (mTotalCrossSection == 0) {
//
// Limits of integration in t, Q2, W2
// or in the case of UPC t, xpom, dummy
//
double xmin[3];
double xmax[3];
copy(mLowerLimit, mLowerLimit+3, xmin);
copy(mUpperLimit, mUpperLimit+3, xmax);
//
// At this point mLowerLimit[1] and mUpperLimit[1]
// are in log(Q2) or log(xpom).
//
xmin[1] = exp(xmin[1]); // log Q2 limit -> Q2 limit or equivalent xpom for UPC
xmax[1] = exp(xmax[1]); // log Q2 limit -> Q2 limit
mTotalCrossSection = calculateTotalCrossSection(xmin, xmax);
}
return mTotalCrossSection;
}
double Sartre::totalCrossSection(double lower[3], double upper[3]) // t, Q2, W (or t, xpom, dummy for UPC)
{
lower[2] *= lower[2]; upper[2] *= upper[2]; // W -> W2
double result = calculateTotalCrossSection(lower, upper);
return result;
}
EventGeneratorSettings* Sartre::runSettings()
{
return EventGeneratorSettings::instance();
}
double Sartre::calculateTotalCrossSection(double lower[3], double upper[3])
{
double result = 0;
if (!mIsInitialized) {
cout << "Sartre::calculateTotalCrossSection(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return result;
}
//
// Calculate integral using adaptive numerical method
//
// Options: ADAPTIVE, kVEGAS, kPLAIN, kMISER
// no abs tolerance given -> relative only
const double precision = 5e-4;
ROOT::Math::Functor wfL((*mCrossSection), mSettings->UPC() ? 2 : 3);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE,
0, precision, 1000000);
ig.SetFunction(wfL);
result = ig.Integral(lower, upper);
//
// If it fails we switch to a MC integration which is usually more robust
// although not as accurate. This should happen very rarely if at all.
//
if (result <= numeric_limits