Page MenuHomeHEPForge
Files F8723629

No OneTemporary
Actions

View Options

Index: trunk/src/TableGeneratorSettings.h
===================================================================
--- trunk/src/TableGeneratorSettings.h (revision 533)
+++ trunk/src/TableGeneratorSettings.h (revision 534)
@@ -1,131 +1,154 @@
//==============================================================================
// TableGeneratorSettings.h
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// 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 setBetamin(double);
+ double betamin() const;
+
+ void setBetamax(double);
+ double betamax() const;
+
+ void setZmin(double);
+ double zmin() const;
+
+ void setZmax(double);
+ double zmax() 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);
+ void setBetabins(unsigned int);
+ unsigned int betabins() const;
+
+ void setZbins(unsigned int);
+ unsigned int zbins() const;
+
+ string bSatLookupPath() const;
+ void setBSatLookupPath(string);
unsigned int numberOfConfigurations() const;
void setNumberOfConfigurations(unsigned int);
vector<double> 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);
bool hasSubstructure() const;
void setHasSubstructure(bool);
double fractionOfBinsToFill() const;
void setFractionOfBinsToFill(double);
void consolidateSettings();
private:
TableGeneratorSettings();
private:
static TableGeneratorSettings* mInstance;
private:
unsigned int mNumberOfConfigurations;
vector<double> mDipoleModelCustomParameters; // developer
bool mUseBackupFile;
int mStartingBinFromBackup;
int mStartBin, mEndBin;
int mModesToCalculate;
double mTmin;
double mTmax;
double mXmin;
- double mXmax;
+ double mXmax;
+
+ double mBetamin, mBetamax;
+ double mZmin, mZmax;
unsigned int mQ2bins;
unsigned int mW2bins;
unsigned int mTbins;
- unsigned int mXbins;
+ unsigned int mXbins;
+ unsigned int mBetabins;
+ unsigned int mZbins;
string mBSatLookupPath;
int mPriority;
bool mHasSubstructure;
double mFractionOfBinsToFill;
};
#endif
Index: trunk/src/Settings.h
===================================================================
--- trunk/src/Settings.h (revision 533)
+++ trunk/src/Settings.h (revision 534)
@@ -1,210 +1,217 @@
//==============================================================================
// Settings.h
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#ifndef Settings_h
#define Settings_h
#include <iostream>
#include <string>
#include <vector>
#include <map>
#include "TDatabasePDG.h"
#include "TRandom3.h"
#include "Enumerations.h"
using namespace std;
class TParticlePDG;
class SettingsParameterBase {
public:
virtual ~SettingsParameterBase() {}
public:
string name;
};
template<typename T> class SettingsParameter : public SettingsParameterBase {
public:
T* address;
T defaultValue;
};
class Settings {
public:
Settings();
virtual ~Settings();
bool readSettingsFromFile(const char*);
virtual bool list(ostream& = cout);
TParticlePDG* lookupPDG(int) const;
string particleName(int pdgID);
static TRandom3* randomGenerator();
unsigned int seed() const;
void setSeed(unsigned int);
int userInt() const;
double userDouble() const;
string userString() const;
void setUserInt(int);
void setUserDouble(double);
void setUserString(const string&);
void setVerbose(bool);
bool verbose() const;
void setVerboseLevel(int);
int verboseLevel() const;
void setQ2min(double);
double Q2min() const;
double Qmin() const;
void setQ2max(double);
double Q2max() const;
double Qmax() const;
void setWmin(double);
void setW2min(double);
double Wmin() const;
double W2min() const;
void setWmax(double);
void setW2max(double);
double Wmax() const;
double W2max() const;
void setXpomMin(double); // UPC only
void setXpomMax(double);
double xpomMin() const;
double xpomMax() const;
void setbetamin(double); // Inclusive diffraction
double betamin() const;
void setbetamax(double);
double betamax() const;
+ void setZmin(double); // Inclusive diffraction
+ double zmin() const;
+ void setZmax(double);
+ double zmax() const;
+
int vectorMesonId() const;
void setVectorMesonId(int);
static double quarkMass(int);
static void setQuarkMass(int, double);
string dipoleModelName() const;
DipoleModelType dipoleModelType() const;
void setDipoleModelType(DipoleModelType);
string dipoleModelParameterSetName() const;
DipoleModelParameterSet dipoleModelParameterSet() const;
void setDipoleModelParameterSet(DipoleModelParameterSet);
string tableSetTypeName() const;
TableSetType tableSetType() const;
void setTableSetType(TableSetType);
unsigned int A() const;
void setA(unsigned int);
string rootfile() const;
void setRootfile(const char*);
void setInclusiveDiffractionMode(bool);
bool inclusiveDiffractionMode() const;
bool exclusiveDiffractionMode() const;
void setUPC(bool);
bool UPC() const;
void setUPCA(unsigned int);
unsigned int UPCA() const;
protected:
template<typename T> void registerParameter(T*, const char*, T);
virtual void consolidateSettings() = 0;
virtual void consolidateCommonSettings();
protected:
vector<SettingsParameterBase*> mRegisteredParameters;
static TRandom3 mRandomGenerator;
unsigned int mSeed;
static double mQuarkMass[6];
bool mVerbose;
int mVerboseLevel;
double mQ2min;
double mQ2max;
double mWmin;
double mWmax;
double mXpomMin;
double mXpomMax;
double mbetamin;
double mbetamax;
+ double mZmin;
+ double mZmax;
int mVectorMesonId;
string mDipoleModelName;
DipoleModelType mDipoleModelType;
DipoleModelParameterSet mDipoleModelParameterSet;
string mDipoleModelParameterSetName;
TableSetType mTableSetType;
string mTableSetTypeName;
unsigned int mA;
string mRootfile;
bool mUPC;
unsigned int mUPCA;
bool mInclusiveDiffractionMode; // if false exclusive
private:
string mRuncard;
vector<string> mLines;
TDatabasePDG *mPDG;
map<int, string> mPeriodicTable;
int mUserInt;
double mUserDouble;
string mUserString;
};
template<typename T> void Settings::registerParameter(T* var, const char* name, T def)
{
SettingsParameter<T>* sv = new SettingsParameter<T>;
sv->address = var;
sv->name = name;
sv->defaultValue = def;
*var = def;
mRegisteredParameters.push_back(sv);
}
inline TRandom3* Settings::randomGenerator() {return &mRandomGenerator;}
inline double Settings::quarkMass(int i) {return mQuarkMass[i];}
inline void Settings::setQuarkMass(int i, double val) {mQuarkMass[i] = val;}
#endif
Index: trunk/src/TableGeneratorSettings.cpp
===================================================================
--- trunk/src/TableGeneratorSettings.cpp (revision 533)
+++ trunk/src/TableGeneratorSettings.cpp (revision 534)
@@ -1,217 +1,242 @@
//==============================================================================
// TableGeneratorSettings.cpp
//
// Copyright (C) 2010-2024 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "Settings.h"
#include "TableGeneratorSettings.h"
#include "Constants.h"
#include <cmath>
#include <ctime>
#include <vector>
#include <stdlib.h>
using namespace std;
TableGeneratorSettings* TableGeneratorSettings::mInstance = 0; // initialize static
TableGeneratorSettings* TableGeneratorSettings::instance()
{
if (mInstance == 0)
mInstance = new TableGeneratorSettings;
return mInstance;
}
TableGeneratorSettings::TableGeneratorSettings()
{
//
// Register all the parameters that can be defined
// via a runcard.
// Arguments for registerParameter():
// 1. pointer to data memeber
// 2. text string to be used in the runcard
// 3. default parameter set
//
registerParameter(&mBSatLookupPath, "bSatLookupPath", string("./"));
registerParameter(&mTmin, "tmin", -2.);
registerParameter(&mTmax, "tmax", 0.);
registerParameter(&mXmin, "xmin", 1e-9); //UPC
registerParameter(&mXmax, "xmax", 3e-2); //UPC
-
+
+ registerParameter(&mBetamin, "betamin", 1e-9); //inc.diff.
+ registerParameter(&mBetamax, "betamax", 0.9); //inc.diff.
+ registerParameter(&mZmin, "zmin", 1e-9); //inc.diff.
+ registerParameter(&mZmax, "zmax", 0.9); //inc.diff.
+
registerParameter(&mQ2bins, "Q2bins", static_cast<unsigned int>(1));
registerParameter(&mW2bins, "W2bins", static_cast<unsigned int>(1));
registerParameter(&mTbins, "tbins", static_cast<unsigned int>(1));
registerParameter(&mXbins, "xbins", static_cast<unsigned int>(1)); //UPC
-
+ registerParameter(&mBetabins, "betabins", static_cast<unsigned int>(1)); //inc.diff.
+ registerParameter(&mZbins, "zbins", static_cast<unsigned int>(1)); //inc.diff.
+
registerParameter(&mNumberOfConfigurations, "numberOfConfigurations", static_cast<unsigned int>(1000));
vector<double> vec;
registerParameter(&mDipoleModelCustomParameters, "dipoleModelCustomParameters", vec);
registerParameter(&mUseBackupFile, "useBackupFile", false);
registerParameter(&mStartingBinFromBackup, "startingBinFromBackup", 0);
registerParameter(&mStartBin, "startBin", -1);
registerParameter(&mEndBin, "endBin", -1);
registerParameter(&mModesToCalculate, "modesToCalculate", 0);
registerParameter(&mPriority, "priority", 0);
registerParameter(&mHasSubstructure, "hasSubstructure", false);
registerParameter(&mFractionOfBinsToFill, "fractionOfBinsToFill", 1.);
}
void TableGeneratorSettings::consolidateSettings() // called after runcard is read
{
//
// Kinematic limits
//
if (mQ2min>=mQ2max && !mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, Q2min >= Q2max. Stopping" << endl;
exit(1);
}
if (mWmin>=mWmax && !mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, Wmin >= Wmax. Stopping" << endl;
exit(1);
}
if (mTmin>=mTmax) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, tmin >= tmax. Stopping" << endl;
exit(1);
}
if (mTmin>0. || mTmax >0.) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, t must be negative, please change t-limits. Stopping" << endl;
exit(1);
}
if (mXmin>=mXmax && mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, xmin >= xmax. Stopping" << endl;
exit(1);
}
if ((mXmin>=1 || mXmin<=0 || mXmax>=1 || mXmax<=0) && mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, xmin or xmax out of range. Stopping" << endl;
exit(1);
}
if (mXmax>0.01 && mUPC){
cout << "TableGeneratorSettings::consolidateSettings(): Warning, xmax>1e-2, model may be unreliable." << endl;
}
if (mA==1 && !mHasSubstructure) mNumberOfConfigurations = 1;
if (!mUseBackupFile) mStartingBinFromBackup = 0;
if (!mUPC)
mXbins=1;
else {
mQ2bins=1;
mW2bins=1;
}
if (mStartBin >= 0 && mEndBin >= 0 && mStartBin>mEndBin) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, endBin < startBin : " << mEndBin << " <= " << mStartBin << "! Stopping." << endl;
exit(1);
}
if ( mStartBin < 0 ) mStartBin=0;
if ( mStartBin >= signed(mQ2bins*mW2bins*mTbins*mXbins) ) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, starting bin >= table! Stopping." << endl;
exit(1);
}
if ( mEndBin > signed(mQ2bins*mW2bins*mTbins*mXbins) || mEndBin < 0) {
cout << "TableGeneratorSettings::consolidateSettings(): endBin is set to table size=" << mQ2bins*mW2bins*mTbins << endl;
- mEndBin=mQ2bins*mW2bins*mTbins*mXbins;
+ mEndBin=mQ2bins*mW2bins*mTbins*mXbins*mBetabins*mZbins;
}
if ( mModesToCalculate < 0 || mModesToCalculate > 2 ) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, modesToCalculate can only take values 0, 1, or 2; not "
<< mModesToCalculate << endl;
exit(1);
}
//
// Make sure the W range is allowed
//
double VMMass=lookupPDG(mVectorMesonId)->Mass();
double W2min=VMMass*VMMass+protonMass2;
double Wmin=sqrt(W2min);
if (mWmin<Wmin){
mWmin=Wmin;
cout << "TableGeneratorSettings::consolidateSettings(): Warning, Wmin is smaller than allowed value." << endl;
cout << " It has been changed to Wmin=" << mWmin << endl;
}
}
//
// Access functions
//
void TableGeneratorSettings::setTmax(double val){ mTmax=val; }
double TableGeneratorSettings::tmax() const { return mTmax; }
void TableGeneratorSettings::setTmin(double val){ mTmin=val; }
double TableGeneratorSettings::tmin() const { return mTmin; }
void TableGeneratorSettings::setXmax(double val){ mXmax=val; }
double TableGeneratorSettings::xmax() const { return mXmax; }
void TableGeneratorSettings::setXmin(double val){ mXmin=val; }
double TableGeneratorSettings::xmin() const { return mXmin; }
-
-void TableGeneratorSettings::setQ2bins(unsigned int val){ mQ2bins=val; }
+
+void TableGeneratorSettings::setZmax(double val){ mZmax=val; }
+double TableGeneratorSettings::zmax() const { return mZmax; }
+
+void TableGeneratorSettings::setZmin(double val){ mZmin=val; }
+double TableGeneratorSettings::zmin() const { return mZmin; }
+
+void TableGeneratorSettings::setBetamax(double val){ mBetamax=val; }
+double TableGeneratorSettings::betamax() const { return mBetamax; }
+
+void TableGeneratorSettings::setBetamin(double val){ mBetamin=val; }
+double TableGeneratorSettings::betamin() const { return mBetamin; }
+
+void TableGeneratorSettings::setQ2bins(unsigned int val){ mQ2bins=val; }
unsigned int TableGeneratorSettings::Q2bins() const { return mQ2bins; }
void TableGeneratorSettings::setW2bins(unsigned int val){ mW2bins=val; }
unsigned int TableGeneratorSettings::W2bins() const { return mW2bins; }
void TableGeneratorSettings::setTbins(unsigned int val){ mTbins=val; }
unsigned int TableGeneratorSettings::tbins() const { return mTbins; }
void TableGeneratorSettings::setXbins(unsigned int val){ mXbins=val; }
unsigned int TableGeneratorSettings::xbins() const { return mXbins; }
-
-void TableGeneratorSettings::setBSatLookupPath(string val){ mBSatLookupPath = val; }
+
+void TableGeneratorSettings::setBetabins(unsigned int val){ mBetabins=val; }
+unsigned int TableGeneratorSettings::betabins() const { return mBetabins; }
+
+void TableGeneratorSettings::setZbins(unsigned int val){ mZbins=val; }
+unsigned int TableGeneratorSettings::zbins() const { return mZbins; }
+
+void TableGeneratorSettings::setBSatLookupPath(string val){ mBSatLookupPath = val; }
string TableGeneratorSettings::bSatLookupPath() const { return mBSatLookupPath; }
void TableGeneratorSettings::setNumberOfConfigurations(unsigned int val){ mNumberOfConfigurations=val; }
unsigned int TableGeneratorSettings::numberOfConfigurations() const { return mNumberOfConfigurations; }
vector<double> TableGeneratorSettings::dipoleModelCustomParameters() const {return mDipoleModelCustomParameters;}
void TableGeneratorSettings::setUseBackupFile(bool val){ mUseBackupFile = val; }
bool TableGeneratorSettings::useBackupFile() const { return mUseBackupFile; }
void TableGeneratorSettings::setStartingBinFromBackup(int val){ mStartingBinFromBackup = val; }
int TableGeneratorSettings::startingBinFromBackup() const { return mStartingBinFromBackup; }
void TableGeneratorSettings::setStartBin(int val){ mStartBin = val; }
int TableGeneratorSettings::startBin() const { return mStartBin; }
void TableGeneratorSettings::setEndBin(int val){ mEndBin = val; }
int TableGeneratorSettings::endBin() const{ return mEndBin; }
void TableGeneratorSettings::setModesToCalculate(int val){ mModesToCalculate = val; }
int TableGeneratorSettings::modesToCalculate() const { return mModesToCalculate; }
void TableGeneratorSettings::setPriority(unsigned char val){ mPriority = val; }
unsigned char TableGeneratorSettings::priority() const { return mPriority; }
void TableGeneratorSettings::setHasSubstructure(bool val){ mHasSubstructure = val; }
bool TableGeneratorSettings::hasSubstructure() const { return mHasSubstructure; }
void TableGeneratorSettings::setFractionOfBinsToFill(double val){ mFractionOfBinsToFill = val; }
double TableGeneratorSettings::fractionOfBinsToFill() const { return mFractionOfBinsToFill; }
Index: trunk/src/Table.h
===================================================================
--- trunk/src/Table.h (revision 533)
+++ trunk/src/Table.h (revision 534)
@@ -1,214 +1,239 @@
//==============================================================================
// Table.h
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// To read/load an existing table:
//
// Table tbl;
// tbl.read("filename"); // read table into memory
// double content = tbl.get(Q2, W2, t); // return content for Q2, W2, t
//
// To create a table:
//
// Table tbl;
// int n = tbl.create(nbinQ2, Q2min, Q2max,
// nBinW2, W2min, W2max,
// nBinT, tmin, tmax,
// logQ2, logW2, logT, logC // all bools
// AmplitudeMoment m, GammaPolarization p,
// massA, vmPDG, model);
//
// for (int i=0; i<n; i++) { // filling table
// double Q2, W2, t, content;
// tbl.binCenter(i, Q2, W2, t);
// // calculate cross-section for Q2, W2, t...
// tbl.fill(i, content);
// }
// tbl.write("filename");
//
// Use Table::list() to query the definition and content of a table.
//
//
// Backup mechanism:
//
// void setAutobackup(const char* prefix, int freq);
//
// Will store the current table after each freq'th to fill().
// Filename is <prefix>_backup.<pid>.<id>.<bin>.root
// where
// <prefix> is a unique name passed to setAutobackup()
// <pid> is the process ID (see getpid())
// <bin> is the number of the last bin filled
// <id> is the table ID
//
// If a new backup file is created the old/previous one gets deleted.
// The last valid backup file is deleted after a successful call to write().
//==============================================================================
#ifndef Table_h
#define Table_h
#include <iostream>
#include <string>
#include "Enumerations.h"
#include <stdint.h>
using namespace std;
class TH2F;
class TH3F;
+class THn;
class Table {
public:
Table();
Table(const Table&);
~Table();
Table& operator=(const Table&);
//
// Reading and using a table
//
bool read(const char*);
double get(double Q2, double W2, double t) const;
double get(double xpom, double t) const; // UPC version
//
// Creating, filling, and storing a table
- //
- unsigned int create(int, double, double,
+ //
+ //Exclusive VM version:
+ unsigned int create(int, double, double,
int, double, double,
int, double, double,
bool logQ2, bool logW2, bool logt,
bool logContent,
AmplitudeMoment, GammaPolarization,
unsigned int A, int vm,
DipoleModelType model,
DipoleModelParameterSet pset,
const char* filename = 0,
unsigned char priority = 0,
bool hotspots = false);
-
- unsigned int create(int, double, double, // UPC version
+
+ //Inclusive Diffraction version:
+ unsigned int create(int, double, double, //beta
+ int, double, double, //Q2
+ int, double, double, //W2
+ int, double, double, //z
+ bool logbeta, bool logQ2, bool logW2, bool logz,
+ bool logContent,
+ AmplitudeMoment, GammaPolarization, FockState,
+ unsigned int A,
+ DipoleModelType model,
+ DipoleModelParameterSet pset,
+ const char* filename = 0,
+ unsigned char priority = 0,
+ bool hotspots = false);
+
+ //UPC version:
+ unsigned int create(int, double, double,
int, double, double,
bool logx, bool logt,
bool logContent,
AmplitudeMoment, GammaPolarization,
unsigned int A, int vm,
DipoleModelType model,
DipoleModelParameterSet pset,
const char* filename = 0,
unsigned char priority = 0,
bool hotspots = false);
+
void binCenter(int, double& Q2, double& W2, double& t) const;
void binCenter(int, double& xpom, double& t) const; // UPC version
void fill(int, double, double err = 0);
- void write(const char* filename = 0);
+ void fill(int, int, int, int, double, double err = 0); //Incl. Diff.
+ void write(const char* filename = 0);
//
// Backup mechanism
//
void setAutobackup(const char* prefix, int freq);
//
// Query functions
//
unsigned int A() const; // mass number
bool isTransverse() const;
bool isLongitudinal() const;
GammaPolarization polarization() const;
bool isMeanA2() const;
bool isMeanA() const;
bool isLambdaA() const;
bool isLambdaSkew() const;
bool isVarianceA() const;
AmplitudeMoment moment() const;
int vectorMesonId() const;
unsigned int dipoleModelType() const;
unsigned int dipoleModelParameterSet() const;
int numberOfEntries() const;
double minQ2() const;
double maxQ2() const;
double minW2() const;
double maxW2() const;
double minT() const;
double maxT() const;
double minX() const;
double maxX() const;
-
+ double minBeta() const;
+ double maxBeta() const;
+ double minZ() const;
+ double maxZ() const;
+
double binWidthQ2() const;
double binWidthW2() const;
double binWidthT() const;
double binWidthX() const;
unsigned int priority() const;
bool isUPC() const;
bool hasHotspots() const;
uint64_t id() const;
void list(ostream& = cout, bool = false, bool = false, int = 0, int = 0) const;
string filename() const;
int globalBin(int binx, int biny, int binz) const;
int globalBin(int binx, int biny) const; // UPC version
void binXYZ(int globalBin, int& binx, int& biny, int& binz) const;
void binXY(int globalBin, int& binx, int& biny) const; // UPC version
bool writeToBinaryFile(const string&, bool verbose=false) const;
private:
bool isLogQ2() const;
bool isLogW2() const;
bool isLogT() const;
bool isLogX() const;
bool isLogContent() const;
bool fexists(const char*) const;
double InterpolateGridSpline(double x, double y, double z) const;
double modInterpolation(double, double, double) const;
double modInterpolation(double, double) const; // UPC version
double centerOverHighestNeighborRatio(int, int, int = 0) const;
void backup(int);
private:
uint64_t mID;
TH3F *m3DHisto;
TH2F *m2DHisto; // UPC table
+ THn *m4DHisto; // Incl. Diff.
string mFilename;
int mFillCounter;
int mBackupFrequence;
string mBackupPrefix;
string mLastBackupFilename;
};
#endif
Index: trunk/src/Table.cpp
===================================================================
--- trunk/src/Table.cpp (revision 533)
+++ trunk/src/Table.cpp (revision 534)
@@ -1,1513 +1,1656 @@
//==============================================================================
// Table.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// 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 <A>, 1 for <A2> (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_<A> (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: 0 = no nucleon substructe fluctuations, 1 = MS-hs (hotspots)
// bit 50-63: not used
//
// log implies ln (not log10)
// 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 "THn.h"
#include <cstdio>
#include <cmath>
#include <ctime>
#include <algorithm>
#include <string>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <limits>
#include <unistd.h>
#include <sys/types.h>
#include "GridSpline.h"
#define PR(x) cout << #x << " = " << (x) << endl;
Table::Table()
{
mID = 0;
- m3DHisto = 0;
+ m4DHisto = 0;
+ m3DHisto = 0;
m2DHisto = 0;
mFillCounter = 0;
mBackupFrequence = 0;
}
Table& Table::operator=(const Table& tab)
{
if (this != &tab) {
- delete m3DHisto;
+ delete m4DHisto;
+ delete m3DHisto;
delete m2DHisto;
mID = tab.mID;
+ if (tab.m4DHisto) {
+ m4DHisto = new THnF();//*(tab.m4DHisto));
+// m4DHisto->SetDirectory(0);
+ }
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;
+ delete m4DHisto;
}
-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,
+//
+// Exclusive Diffraction
+//
+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,
+ int quarkSpecies,
const char* filename, unsigned char priority,
bool hotspots)
-{
+{
if (m3DHisto != 0) {
cout << "Table:create(): cannot create, table already exists." << endl;
return 0;
}
if (nbinsQ2 < 2) {
cout << "Table:create(): number of bins in Q2 must be at least 2." << endl;
return 0;
}
if (nbinsW2 < 2) {
cout << "Table:create(): number of bins in W2 must be at least 2." << endl;
return 0;
}
if (nbinsT < 2) {
cout << "Table:create(): number of bins in t must be at least 2." << endl;
return 0;
}
if (filename) mFilename = string(filename); // name of file where table is written to
//
// Create table ID first
//
mID = (vm << 16);
mID |= (A << 8);
mID |= (model << 5);
mID |= (static_cast<uint64_t>(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<uint64_t>(1) << 32);
if (mom == lambda_real) mID |= (static_cast<uint64_t>(1) << 33);
if (mom == lambda_skew) mID |= (static_cast<uint64_t>(1) << 48);
mID |= (static_cast<uint64_t>(priority) << 34);
if (mom == variance_A) mID |= (static_cast<uint64_t>(1) << 45);
if (hotspots) mID |= (static_cast<uint64_t>(1) << 49);
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 Inclusive Diffraction
+//
+unsigned int Table::create(int nbinsbeta, double betamin, double betamax,
+ int nbinsQ2, double Q2min, double Q2max,
+ int nbinsW2, double W2min, double W2max,
+ int nbinsz, double zmin, double zmax,
+ bool logbeta, bool logQ2, bool logW2, bool logz, bool logC,
+ AmplitudeMoment mom, GammaPolarization pol,
+ FockState fock,
+ unsigned int A,
+ DipoleModelType model, DipoleModelParameterSet pset,
+ const char* filename, unsigned char priority,
+ bool hotspots){
+ if (m4DHisto != 0) {
+ cout << "Table:create(): cannot create, table already exists." << endl;
+ return 0;
+ }
+ if (nbinsbeta < 2){
+ cout << "Table:create(): number of bins in beta must be at least 2." << endl;
+ return 0;
+ }
+ if (nbinsQ2 < 2) {
+ cout << "Table:create(): number of bins in Q2 must be at least 2." << endl;
+ return 0;
+ }
+ if (nbinsW2 < 2) {
+ cout << "Table:create(): number of bins in W2 must be at least 2." << endl;
+ return 0;
+ }
+ if (nbinsz < 2) {
+ cout << "Table:create(): number of bins in z must be at least 2." << endl;
+ return 0;
+ }
+
+ if (filename) mFilename = string(filename); // name of file where table is written to
+
+ //
+ // Create table ID first
+ // TODO: Figure this out for inclusive diffraction
+ // TODO: There is no VM, logt->logz, add logbeta,
+ // TODO: there are no corrections
+// mID = (vm << 16);
+ mID |= (A << 8);
+ mID |= (model << 5);
+ mID |= (static_cast<uint64_t>(pset) << 42);
+ if (mom == mean_A2) mID |= 1;
+ if (pol == transverse) mID |= (1 << 1);
+ if (logz) mID |= (1 << 2);
+ if (logW2) mID |= (1 << 3);
+ if (logQ2) mID |= (1 << 4);
+ if (logC) mID |= (static_cast<uint64_t>(1) << 32);
+ if (mom == lambda_real) mID |= (static_cast<uint64_t>(1) << 33);
+ if (mom == lambda_skew) mID |= (static_cast<uint64_t>(1) << 48);
+ mID |= (static_cast<uint64_t>(priority) << 34);
+ if (mom == variance_A) mID |= (static_cast<uint64_t>(1) << 45);
+ if (hotspots) mID |= (static_cast<uint64_t>(1) << 49);
+
+ 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.
+ //
+ if (logbeta) betamin = log(betamin);
+ if (logbeta) betamax = log(betamax);
+ if (logQ2) Q2min = log(Q2min);
+ if (logQ2) Q2max = log(Q2max);
+ if (logW2) W2min = log(W2min);
+ if (logW2) W2max = log(W2max);
+ if (logz) zmin = log(zmin);
+ if (logz) zmax = log(zmax);
+ double binWidthbeta = (betamax - betamin)/(nbinsbeta-1);
+ double binWidthQ2 = (Q2max - Q2min)/(nbinsQ2-1);
+ double binWidthW2 = (W2max - W2min)/(nbinsW2-1);
+ double binWidthZ = (zmax - zmin)/(nbinsz-1);
+
+ //
+ // Book 4D histo to hold table
+ //
+ const int nDims = 4; // Number of dimensions
+ const int nbins[nDims] = {nbinsbeta, nbinsQ2, nbinsW2, nbinsz}; // Number of bins for each axis
+ const double minRange[nDims] = {betamin-binWidthbeta/2, Q2min-binWidthQ2/2, W2min-binWidthW2/2, zmin-binWidthZ/2}; // Minimum range for each axis
+ const double maxRange[nDims] = {betamax+binWidthbeta/2, Q2max+binWidthQ2/2, W2max+binWidthW2/2, zmax+binWidthZ/2}; // Maximum range for each axis
+
+ // Create a 4D histogram
+ m4DHisto = new THnF("table", title.c_str(), nDims, nbins, minRange, maxRange);
+
+ return nbinsbeta*nbinsQ2*nbinsW2*nbinsz; // 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,
bool hotspots)
{
if (m2DHisto != 0) {
cout << "Table:create(): cannot create, table already exists." << endl;
return 0;
}
if (nbinsX < 2) {
cout << "Table:create(): number of bins in x must be at least 2." << endl;
return 0;
}
if (nbinsT < 2) {
cout << "Table:create(): number of bins in t must be at least 2." << endl;
return 0;
}
if (filename) mFilename = string(filename); // name of file where table is written to
//
// Create table ID first
//
mID = (vm << 16);
mID |= (A << 8);
mID |= (model << 5);
mID |= (static_cast<uint64_t>(pset) << 42);
if (mom == mean_A2) mID |= 1;
if (pol == transverse) mID |= (1 << 1);
if (logt) mID |= (1 << 2);
if (logC) mID |= (static_cast<uint64_t>(1) << 32);
if (mom == lambda_real) mID |= (static_cast<uint64_t>(1) << 33);
if (mom == lambda_skew) mID |= (static_cast<uint64_t>(1) << 48);
mID |= (static_cast<uint64_t>(priority) << 34);
if (mom == variance_A) mID |= (static_cast<uint64_t>(1) << 45);
mID |= (static_cast<uint64_t>(1) << 46); // flag UPC
if (logX) mID |= (static_cast<uint64_t>(1) << 47);
if (hotspots) mID |= (static_cast<uint64_t>(1) << 49);
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;
}
+ if (m4DHisto) {
+ return m4DHisto->GetNbins();
+ }
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)
+//Inclusive Diffraction version:
+void Table::fill(int betaBin, int Q2Bin, int W2Bin, int zBin, double val, double err){
+ if (isLogContent()) {
+ if (val == 0) {
+ // cout << "Table::fill(): warning, log scale requested but value = 0." << endl;
+ val = numeric_limits<float>::min();
+ }
+ val = log(fabs(val));
+ }
+ const int currentBins[4]={betaBin, Q2Bin, W2Bin, zBin};
+ int currentBin = m4DHisto->GetBin(currentBins);
+ m4DHisto->SetBinContent(currentBins, val);
+ if (err!=0.) {
+ m4DHisto->SetBinError2(currentBin, err);
+ }
+ mFillCounter++;
+ //
+ // Check if backup is due
+ //
+ if (mBackupFrequence)
+ if (mFillCounter%mBackupFrequence == 0) backup(currentBin);
+}
+
+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<float>::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();
+ else{
+ if(m3DHisto)
+ m3DHisto->Write();
+ if(m4DHisto)
+ m4DHisto->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 m4DHisto; m4DHisto = 0;
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<TH2F*>(ptr);
m2DHisto->SetDirectory(0);
}
else {
m3DHisto = reinterpret_cast<TH3F*>(ptr);
m3DHisto->SetDirectory(0);
}
if (m2DHisto == 0 && m3DHisto == 0) {
cout << "Table::read(): failed retrieving table from file '" << filename << "'." << endl;
return false;
}
file->Close();
mFilename = string(filename);
return true;
}
int Table::vectorMesonId() const {return ((mID >> 16) & 0xFFFF);}
unsigned int Table::dipoleModelType() const {return ((mID >> 5) & 0x7);}
unsigned int Table::dipoleModelParameterSet() const {return ((mID >> 42) & 0x7);}
unsigned int Table::A() const {return ((mID >> 8) & 0xFF);}
bool Table::isLongitudinal() const {return !isTransverse();}
bool Table::isTransverse() const {return (mID & (1 << 1));}
GammaPolarization Table::polarization() const {return (isTransverse() ? transverse : longitudinal);}
AmplitudeMoment Table::moment() const {
if (isLambdaA())
return lambda_real;
else if (isMeanA())
return mean_A;
else if (isVarianceA())
return variance_A;
else if (isLambdaSkew())
return lambda_skew;
else
return mean_A2;
}
bool Table::isMeanA() const {
// bits 0, 33, 45 and 48 are not set
return !(mID & 1) && !(mID & (static_cast<uint64_t>(1) << 33)) && !(mID & (static_cast<uint64_t>(1) << 45))
&& !(mID & (static_cast<uint64_t>(1) << 48));
}
bool Table::isMeanA2() const {
// bit 0 is set, bits 33, 45 and 48 are not set
return (mID & 1) && !(mID & (static_cast<uint64_t>(1) << 33)) && !(mID & (static_cast<uint64_t>(1) << 45))
&& !(mID & (static_cast<uint64_t>(1) << 48));
}
bool Table::isLambdaA() const {
// bit 33 is set, bits 0 and 45 are not set
return !(mID & 1) && (mID & (static_cast<uint64_t>(1) << 33)) && !(mID & (static_cast<uint64_t>(1) << 45))
&& !(mID & (static_cast<uint64_t>(1) << 48));
}
bool Table::isLambdaSkew() const {
// bit 48 is set, bits 0, 33, and 45 are not set
return !(mID & 1) && (mID & (static_cast<uint64_t>(1) << 48)) && !(mID & (static_cast<uint64_t>(1) << 33))
&& !(mID & (static_cast<uint64_t>(1) << 45));
}
bool Table::isVarianceA() const {
// bit 45 is set, bitss 0, 33, and 48 are not set
return !(mID & 1) && !(mID & (static_cast<uint64_t>(1) << 33)) && (mID & (static_cast<uint64_t>(1) << 45))
&& !(mID & (static_cast<uint64_t>(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<uint64_t>(1) << 47));}
bool Table::isLogContent() const {return (mID & (static_cast<uint64_t>(1) << 32));}
unsigned int Table::priority() const {return ((mID >> 34) & 0xFF);}
bool Table::isUPC() const {return (mID & (static_cast<uint64_t>(1) << 46));}
bool Table::hasHotspots() const {return (mID & (static_cast<uint64_t>(1) << 49));}
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<float>::epsilon());
x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::epsilon());
axis = m3DHisto->GetYaxis();
y = max(y, axis->GetBinCenter(1)+numeric_limits<float>::epsilon());
y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::epsilon());
axis = m3DHisto->GetZaxis();
z = max(z, axis->GetBinCenter(1)+numeric_limits<float>::epsilon());
z = min(z, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::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<float>::epsilon());
x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::epsilon());
axis = m2DHisto->GetYaxis();
y = max(y, axis->GetBinCenter(1)+numeric_limits<float>::epsilon());
y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits<float>::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; i<numberOfEntries(); i++) {
if (isUPC()) {
binXY(i, binx, biny);
binCenter(i, xpom, t);
rawContent = m2DHisto->GetBinContent(binx, biny);
if (isLogContent()) {
rawContent = exp(rawContent);
}
binaryFile.write(reinterpret_cast<const char *>(&i), sizeof(i));
binaryFile.write(reinterpret_cast<const char *>(&xpom), sizeof(xpom));
binaryFile.write(reinterpret_cast<const char *>(&t), sizeof(t));
binaryFile.write(reinterpret_cast<const char *>(&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<const char *>(&i), sizeof(i));
binaryFile.write(reinterpret_cast<const char *>(&Q2), sizeof(Q2));
binaryFile.write(reinterpret_cast<const char *>(&W2), sizeof(W2));
binaryFile.write(reinterpret_cast<const char *>(&t), sizeof(t));
binaryFile.write(reinterpret_cast<const char *>(&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;
+ if ((isUPC() && !m2DHisto) || (!isUPC() && !m3DHisto && !m4DHisto)) {
+ 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 << "<A>";
else if (isMeanA2())
os << "<A^2>";
else if (isLambdaA())
os << "lambda_<A>";
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() == STU)
os << "STU";
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;
//
// Fourth row: hotspots
//
os << setw(34) << right << "hotspots = " << (hasHotspots() ? "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 {
+ else if(m3DHisto) {
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;
}
-
+ else{
+ cout<<"Information is currently missing..."<<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; i<endBin; i++) {
if (isUPC()) {
binXY(i, binx, biny);
binCenter(i, xpom, t);
double value = get(xpom, t);
rawContent = m2DHisto->GetBinContent(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 << " ";
os << " peak = " << fixed << setprecision(3) << centerOverHighestNeighborRatio(binx, biny) << ')';
}
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 << " ";
os << "peak = " << fixed << setprecision(3) << centerOverHighestNeighborRatio(binx, biny, binz) << ')';
}
if ( (isLogContent() && rawContent <= log(numeric_limits<float>::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<float>::max();
int minBin, maxBin;
minBin = maxBin = 0;
double c = 0;
for (int i=0; i<numberOfEntries(); i++) {
if (isUPC()) {
binXY(i, binx, biny);
c = m2DHisto->GetBinContent(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::centerOverHighestNeighborRatio(int binx, int biny, int binz) const
{
double maxRatio = 0;
if (isUPC()) {
double centerValue = m2DHisto->GetBinContent(binx, biny);
int maxBinX = m2DHisto->GetXaxis()->GetNbins();
int maxBinY = m2DHisto->GetYaxis()->GetNbins();
for (int ix = -1; ix < 2; ix++) {
for (int iy = -1; iy < 2; iy++) {
int ixx = binx + ix;
int iyy = biny + iy;
if (ixx > 0 && ixx <= maxBinX && iyy > 0 && iyy <= maxBinY) {
double neighborValue = m2DHisto->GetBinContent(ixx, iyy);
double ratio = centerValue/neighborValue;
if (ratio > maxRatio) maxRatio = ratio;
}
}
}
}
else {
double centerValue = m3DHisto->GetBinContent(binx, biny, binz);
int maxBinX = m3DHisto->GetXaxis()->GetNbins();
int maxBinY = m3DHisto->GetYaxis()->GetNbins();
int maxBinZ = m3DHisto->GetZaxis()->GetNbins();
for (int ix = -1; ix < 2; ix++) {
for (int iy = -1; iy < 2; iy++) {
for (int iz = -1; iz < 2; iz++) {
int ixx = binx + ix;
int iyy = biny + iy;
int izz = binz + iz;
if (ixx > 0 && ixx <= maxBinX && iyy > 0 && iyy <= maxBinY && izz > 0 && izz <= maxBinZ) {
double neighborValue = m3DHisto->GetBinContent(ixx, iyy, izz);
double ratio = centerValue/neighborValue;
if (ratio > maxRatio) maxRatio = ratio;
}
}
}
}
}
return maxRatio;
}
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!"<<endl;
return 0;
}
// Find the bins surrounding the point, lower bins and upper bins
int lbx = m3DHisto->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() && v[i] <= log(numeric_limits<float>::min()*2) )
logC=false;
if (!isLogContent() && v[i] == 0)
logC=false;
}
// Make interpolation in log(content) except for when at least one of the points are zero.
if (isLogContent() && !logC){
for (int i=0; i<8; i++){
if (v[i] <= log(numeric_limits<float>::min()*2))
v[i] = 0;
else
v[i]=exp(v[i]);
}
}
if (!isLogContent() && 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() && !logC){
if (result == 0)
result=log(numeric_limits<float>::min());
else
result=log(result);
}
if (!isLogContent() && 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!"<<endl;
return 0;
}
// Bins that contain the requested values
int bin_x = m2DHisto->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<float>::min()*2)) ||
(!isLogContent() && q11 <= 0)) logC = false;
if ((isLogContent() && q12 <= log(numeric_limits<float>::min()*2)) ||
(!isLogContent() && q12 <= 0)) logC = false;
if ((isLogContent() && q21 <= log(numeric_limits<float>::min()*2)) ||
(!isLogContent() && q21 <= 0)) logC = false;
if ((isLogContent() && q22 <= log(numeric_limits<float>::min()*2)) ||
(!isLogContent() && q22 <= 0)) logC = false;
if (isLogContent() && !logC) { // 0's present, use linear content
if (q11 <= log(numeric_limits<float>::min()*2))
q11 = 0;
else
q11=exp(q11);
if (q12 <= log(numeric_limits<float>::min()*2))
q12 = 0;
else
q12=exp(q12);
if (q21 <= log(numeric_limits<float>::min()*2))
q21 = 0;
else
q21=exp(q21);
if (q22 <= log(numeric_limits<float>::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<float>::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<int>(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/InclusiveDiffractiveCrossSections.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 533)
+++ trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 534)
@@ -1,687 +1,688 @@
//==============================================================================
// InclusiveDiffractiveCrossSection.cpp
//
// Copyright (C) 2024 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
// $Author: ullrich $
//==============================================================================
#include "Math/SpecFunc.h"
#include "TFile.h"
#include <iostream>
#include <cmath>
#include <algorithm>
#include <algorithm>
#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 "Settings.h"
#include "Enumerations.h"
#include "TF1.h"
#include "TF1.h"
#include "TH1F.h"
#include "cuba.h"
#include "InclusiveDiffractiveCrossSections.h"
#include "DglapEvolution.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include "DipoleModelParameters.h"
#define PR(x) cout << #x << " = " << (x) << endl;
InclusiveDiffractiveCrossSections::InclusiveDiffractiveCrossSections()
{
mSettings = EventGeneratorSettings::instance();
mRandom = mSettings->randomGenerator();
mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
mPhotonFlux.setS(mS);
mCheckKinematics = true;
mCrossSectionRatioLT = 0;
mQ2=0;
- // TableGeneratorSettings* settings = TableGeneratorSettings::instance();
+// TableGeneratorSettings* settings = TableGeneratorSettings::instance();
if(mDipoleModel) delete mDipoleModel;
+
mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
DipoleModelType model=mSettings->dipoleModelType();
- if (model==bSat) {
+ if (model==bSat) {
mDipoleModel = new DipoleModel_bSat(mSettings);
}
else if(model==bNonSat){
mDipoleModel = new DipoleModel_bNonSat(mSettings);
}
else {
cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
exit(1);
}
mA=mSettings->A();
//Set up the integrals
//Amplitude0:
if(Amp0) delete Amp0;
Amp0=new TF1("Amp0", this, &InclusiveDiffractiveCrossSections::uiAmplitude_0, 0, 30., 5);
if(WFAmp0) delete WFAmp0;
WFAmp0=new ROOT::Math::WrappedTF1(*Amp0);
GIAmp0.SetFunction(*WFAmp0);
- GIAmp0.SetRelTolerance(1e-2);
-
+ GIAmp0.SetRelTolerance(1e-4);
+
//Amplitude1:
if(Amp1) delete Amp1;
Amp1=new TF1("Amp1", this, &InclusiveDiffractiveCrossSections::uiAmplitude_1, 0, 30. , 5);
if(WFAmp1) delete WFAmp1;
WFAmp1=new ROOT::Math::WrappedTF1(*Amp1);
GIAmp1.SetFunction(*WFAmp1);
- GIAmp1.SetRelTolerance(1e-2);
+ GIAmp1.SetRelTolerance(1e-4);
//AmplitudeQQG:
if(Ampqqg) delete Ampqqg;
Ampqqg=new TF1("Ampqqg", this, &InclusiveDiffractiveCrossSections::uiAmplitude_qqg, 0, 30., 4);
if(WFAmpqqg) delete WFAmpqqg;
WFAmpqqg=new ROOT::Math::WrappedTF1(*Ampqqg);
GIAmpqqg.SetFunction(*WFAmpqqg);
- GIAmpqqg.SetRelTolerance(1e-2);
+ GIAmpqqg.SetRelTolerance(1e-4);
}
InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){
if(WFAmp0) delete WFAmp0;
if(Amp0) delete Amp0;
if(WFAmp1) delete WFAmp1;
if(Amp1) delete Amp1;
if(Ampqqg) delete Ampqqg;
}
void InclusiveDiffractiveCrossSections::setCheckKinematics(bool val) {mCheckKinematics = val;}
void InclusiveDiffractiveCrossSections::setFockState(FockState val){
mFockState=val;
}
FockState InclusiveDiffractiveCrossSections::getFockState(){
return mFockState;
}
double InclusiveDiffractiveCrossSections::operator()(double MX2, double Q2, double W2, double z){
return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
}
double InclusiveDiffractiveCrossSections::operator()(const double* array){
double result=0;
if(mFockState==QQ or mFockState==ALL){
result+=dsigdbetadQ2dW2dz_total_checked(array[0], array[1], array[2], array[3]);
}
else if(mFockState==QQG or
mFockState==ALL){
result+=dsigdbetadQ2dW2dz_total_qqg_checked(array[0], array[1], array[2], array[3]);
}
else{
cout<<"InclusiveDiffractiveCrossSections::operator(): Fockstate is invalid, stopping"<<endl;
exit(0);
}
return result;
}
double InclusiveDiffractiveCrossSections::unuranPDF(const double* array){
//
// array is: beta, log(Q2), W2, z
//
double result = 0;
double Q2 = exp(array[1]);
result = dsigdbetadQ2dW2dz_total_checked(array[0], Q2, array[2], array[3]);
result *= Q2; // Jacobian
return log(result);
}
double InclusiveDiffractiveCrossSections::dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z){
double beta=Q2/(Q2+MX2);
double jacobian=beta*beta/Q2;
double result= dsigdbetadQ2dW2dz_total_checked(beta, Q2, W2, z);
return jacobian*result;
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z){
double result = 0;
//
// Check if in kinematically allowed region
double MX2=Kinematics::MX2(Q2, beta, 0);
if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
<< " is outside of kinematically allowed region. Return 0." << endl;
return 0;
}
double xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double csT = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, transverse);
double csL = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, longitudinal);
result = csT + csL;
mCrossSectionRatioLT=csL/csT;
//
// Polarization
//
if (mRandom->Uniform(result) <= csT)
mPolarization = transverse;
else
mPolarization = longitudinal;
//
// Quark Species
//
//
if(mPolarization == transverse){
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(csT-cs_rejected);
double csTi=mDsigdbetadQ2dWdz_T_total[i];
if(R <= csTi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_T_total[i];
}
if(!isChosen) mQuarkSpecies=4;
}
else{
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(csL-cs_rejected);
double csLi=mDsigdbetadQ2dWdz_L_total[i];
if(R <= csLi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_L_total[i];
}
if(!isChosen) mQuarkSpecies=4;
}
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): " << result;
cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z;
cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal");
cout << ')' << endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = NaN at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = inf at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (result < 0) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, negative cross-section at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
result = 0;
}
mTotalCS=result;
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
if(pol==transverse){
result = dsigmadbetadz_T(beta, Q2, W2, z);
}
else if(pol==longitudinal){
result = dsigmadbetadz_L(beta, Q2, W2, z);
}
return result; //nb
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
for(int i=0; i<5; i++){
setQuarkIndex(i);
double tmpresult = dsigdbetadz_total(beta, Q2, W2, z, pol); //nb
if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,pol); //nb/GeV4
if(pol == transverse)
mDsigdbetadQ2dWdz_T_total[i]=tmpresult;
if(pol == longitudinal)
mDsigdbetadQ2dWdz_L_total[i]=tmpresult;
result+=tmpresult;
}
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::crossSectionOfLastCall(){
return mTotalCS;
}
GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){
return mPolarization;
}
unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){
return mQuarkSpecies;
}
double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){
// return Settings::quarkMass(mQuarkSpecies); //#TT tmp
tt_quarkMass[mQuarkSpecies];
}
void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {
double MX2 = Q2*(1-beta)/beta;
double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
// double mf = Settings::quarkMass(mQuarkIndex);
double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
double sqrtArg = 1.-4.*mf*mf/MX2;
if (sqrtArg<0){
if (mSettings->verboseLevel() > 2)
cout<<"There is no phase-space for z, return 0."<<endl;
return 0;
}
double z0 = (1.-sqrt(sqrtArg))/2.;
if (z<z0 or z>1-z0){
if (mSettings->verboseLevel() > 2)
cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
return 0;
}
double pt2 = z*(1-z)*Q2-mf*mf;
if (pt2<0){
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
return 0;
}
if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark production, return 0."<<endl;
return 0;
}
double epsilon2=z*(1-z)*Q2+mf*mf;
double ef=quarkCharge[mQuarkIndex];
double Phi0=Phi_0(beta, xpom, z, Q2, mf);
double Phi1=Phi_1(beta, xpom, z, Q2, mf);
double prefactor=Nc*Q2*alpha_em/(4.*M_PI*beta*beta)*ef*ef*z*(1-z); //GeV2
prefactor/=2.; //expanded z-range
double term1=epsilon2*(z*z+(1-z)*(1-z))*Phi1; //GeV2*fm6
double term2=mf*mf*Phi0; //GeV2*fm6
double result=prefactor*(term1+term2); //GeV4*fm6
result/=hbarc2*hbarc2; //fm2
result *= 1e7; //nb
return result;//nb
}
double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double beta, double Q2, double W2, double z){
double MX2 = Q2*(1-beta)/beta;
double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
// double mf = Settings::quarkMass(mQuarkIndex);
double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
double pt2 = z*(1-z)*Q2-mf*mf;
if (pt2<0){
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
return 0;
}
double sqrtArg=1.-4.*mf*mf/MX2;
if (sqrtArg<0){
if (mSettings->verboseLevel() > 2)
cout<<"There is no phase-space for z, return 0."<<endl;
return 0;
}
double z0=(1.-sqrt(sqrtArg))/2.;
if (z<z0 or z>1-z0){
if (mSettings->verboseLevel() > 2)
cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
return 0;
}
if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark production, return 0."<<endl;
return 0;
}
double ef=quarkCharge[mQuarkIndex];
double Phi0=Phi_0(beta, xpom, z, Q2, mf); //fm6
double term=Nc*Q2*Q2*alpha_em/(M_PI*beta*beta); //GeV4
term*=ef*ef*z*z*z*(1-z)*(1-z)*(1-z);
term*=Phi0; //fm6*GeV4
term/=2.; //expanded z-range
double result=term;
result/=hbarc2*hbarc2; //fm2
result *= 1e7; //nb
return result; //nb
}
double InclusiveDiffractiveCrossSections::Phi_0(double beta, double xpom, double z, double Q2, double mf){
TF1* Phi0=new TF1("Phi0", this, &InclusiveDiffractiveCrossSections::uiPhi_0, 0, 30., 0);
ROOT::Math::WrappedTF1* WFPhi0=new ROOT::Math::WrappedTF1(*Phi0);
ROOT::Math::GaussIntegrator GIPhi0;
GIPhi0.SetFunction(*WFPhi0);
- GIPhi0.SetRelTolerance(1e-2);
+ GIPhi0.SetRelTolerance(1e-4);
Amp0->SetParameter(0, beta);
Amp0->SetParameter(1, xpom);
Amp0->SetParameter(2, z);
Amp0->SetParameter(3, Q2);
Amp0->SetParameter(5, mf);
double blow=0;
double result=GIPhi0.IntegralUp(blow); //fm6
delete Phi0;
delete WFPhi0;
return result;
}
double InclusiveDiffractiveCrossSections::Phi_1(double beta, double xpom, double z, double Q2, double mf){
TF1* Phi1=new TF1("Phi1", this, &InclusiveDiffractiveCrossSections::uiPhi_1, 0, 30., 0);
ROOT::Math::WrappedTF1* WFPhi1=new ROOT::Math::WrappedTF1(*Phi1);
ROOT::Math::GaussIntegrator GIPhi1;
GIPhi1.SetFunction(*WFPhi1);
- GIPhi1.SetRelTolerance(1e-2);
+ GIPhi1.SetRelTolerance(1e-4);
Amp1->SetParameter(0, beta);
Amp1->SetParameter(1, xpom);
Amp1->SetParameter(2, z);
Amp1->SetParameter(3, Q2);
Amp1->SetParameter(5, mf);
double blow=0;
double result=GIPhi1.IntegralUp(blow); //fm6
delete Phi1;
delete WFPhi1;
return result;
}
double InclusiveDiffractiveCrossSections::uiPhi_0(double *var, double *par){
double b=*var;
double amplitude=Amplitude_0(b); //fm2
double result=2*M_PI*b*amplitude*amplitude; //fm5
return result;
}
double InclusiveDiffractiveCrossSections::uiPhi_1(double *var, double *par){
double b=*var; //fm
double amplitude=Amplitude_1(b); //fm2
double result=2*M_PI*b*amplitude*amplitude; //fm5
return result;
}
double InclusiveDiffractiveCrossSections::Amplitude_0(double b){
//Gauss Integrator defined in header file and set in constructor
Amp0->SetParameter(4, b);
double rlow=1e-3;
double result=GIAmp0.IntegralUp(rlow); //fm2
return result;
}
double InclusiveDiffractiveCrossSections::Amplitude_1(double b){
//Gauss Integrator defined in header file and set in constructor
Amp1->SetParameter(4, b);
double rlow=1e-3;
double result=GIAmp1.IntegralUp(rlow); //fm2
return result;
}
double InclusiveDiffractiveCrossSections::uiAmplitude_0(double *var, double *par){
double r=*var;
double beta=par[0];
double xpom=par[1];
double z=par[2];
double Q2=par[3];
double b=par[4];
double mf=par[5];
double MX2=Q2*(1-beta)/beta;
double kappa2=z*(1-z)*MX2-mf*mf;
double kappa=sqrt(kappa2);
double epsilon2=z*(1-z)*Q2+mf*mf;
double epsilon=sqrt(epsilon2); //GeV
double besselK0=TMath::BesselK0(epsilon*r/hbarc);
double besselJ0=TMath::BesselJ0(kappa*r/hbarc);
double dsigmadb2=0;
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
return r*besselK0*besselJ0*dsigmadb2; //fm
}
double InclusiveDiffractiveCrossSections::uiAmplitude_1(double *var, double *par){
double r=*var; //fm
double beta=par[0];
double xpom=par[1];
double z=par[2];
double Q2=par[3]; //GeV2
double b=par[4];
double mf=par[5];
double MX2=Q2*(1-beta)/beta; //GeV2
double kappa2=z*(1-z)*MX2-mf*mf; //GeV2
double kappa=sqrt(kappa2); //GeV
double epsilon2=z*(1-z)*Q2+mf*mf;
double epsilon=sqrt(epsilon2); //GeV
double besselK1=TMath::BesselK1(epsilon*r/hbarc);
double besselJ1=TMath::BesselJ1(kappa*r/hbarc);
double dsigmadb2=0; //GeV0
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
return r*besselK1*besselJ1*dsigmadb2; //fm
}
DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel(){
return mDipoleModel;
}
double InclusiveDiffractiveCrossSections::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
//===================================================
// QQG Calculations
//===================================================
double InclusiveDiffractiveCrossSections::unuranPDF_qqg(const double* array){
//
// array is: beta, log(Q2), W2, z
//
double result = 0;
double Q2 = exp(array[1]);
double z=array[3];
result = dsigdbetadQ2dW2dz_total_qqg_checked(array[0], Q2, array[2], z);
result *= Q2; // Jacobian log(Q2)->Q2
return log(result);
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z){
//
// Check if in kinematically allowed region
//
double MX2=Kinematics::MX2(Q2, beta, 0);
if(mCheckKinematics && z<beta) {
if (mSettings->verboseLevel() > 2)
cout<<"z<beta, beta="<<beta<<", z="<<z<<endl;
return 0; // For QQG beta < z < 1
}
if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
<< " is outside of kinematically allowed region. Return 0." << endl;
return 0;
}
double xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double result = dsigdbetadQ2dW2dz_qqg(beta, Q2, W2, z);
//
// Quark Species
//
//
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(result-cs_rejected);
double csTi=mDsigdbetadQ2dWdz_T_qqg[i];
if(R <= csTi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_T_qqg[i];
}
if(!isChosen) mQuarkSpecies=4;
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): " << result;
cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z<< endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = NaN at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = inf at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (result < 0) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, negative cross-section at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
result = 0;
}
mTotalCS_qqg=result;
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z) {
double result = 0;
for(int i=0; i<5; i++){
setQuarkIndex(i);
// double mf = Settings::quarkMass(mQuarkIndex);
double mf = tt_quarkMass[i]; //#TT tmp
double Mqq2=(z/beta-1)*Q2;
if(Mqq2<4*mf*mf) continue; //Not enough phase-space for the q and qbar
double tmpresult = dsigmadbetadz_QQG(beta, Q2, W2, z); //nb
if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,transverse); //nb/GeV4
mDsigdbetadQ2dWdz_T_qqg[i]=tmpresult;
result+=tmpresult;
}
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
double MX2 = Q2*(1-beta)/beta;
double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
double Phi=Phi_qqg(xpom, ztilde, Q2);
double alpha_S=0.14;
double ef=quarkCharge[mQuarkIndex];
double betafactor=(1-beta/ztilde)*
(1-beta/ztilde)+beta/ztilde*beta/ztilde;
double result=alpha_S*alpha_em/(2*M_PI*M_PI*Q2)*betafactor*ef*ef*Phi; //GeV-2
result*=hbarc2; //fm2
result *= 1e7; //nb
return result; //nb
}
double InclusiveDiffractiveCrossSections::Phi_qqg(double xpom, double z, double Q2){
Ampqqg->SetParameter(0, xpom);
Ampqqg->SetParameter(1, z);
mQ2=Q2;
ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiPhi_qqg, 2);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig.SetFunction(wf);
ig.SetAbsTolerance(0.);
- ig.SetRelTolerance(1e-2);
+ ig.SetRelTolerance(1e-4);
double lo[2]={0.,1e-4}; //b, k2
double hi[2]={5., Q2}; //b, k2
double result = ig.Integral(lo, hi)/hbarc; //GeV0
return result; //GeV0
}
double InclusiveDiffractiveCrossSections::uiPhi_qqg(const double* var) {
double b=var[0];
double k2=var[1];
double Q2=mQ2;
Ampqqg->SetParameter(2, sqrt(k2));
Ampqqg->SetParameter(3, b);
double rlow=1e-3;
double amp_qqg=GIAmpqqg.IntegralUp(rlow); //fm2
amp_qqg/=hbarc2; //GeV-2
double result = (2*M_PI)*b/hbarc*k2*k2*log(Q2/k2)*amp_qqg*amp_qqg; //GeV-1
return result; //GeV-1
}
double InclusiveDiffractiveCrossSections::uiAmplitude_qqg(double *var, double *par){
double r=*var; //fm
if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
double xpom=par[0];
double z=par[1];
double k=par[2]; //GeV
double b=par[3];
double besselK2=ROOT::Math::cyl_bessel_k(2, sqrt(z)*k*r/hbarc);
double besselJ2=ROOT::Math::cyl_bessel_j(2, sqrt(1-z)*k*r/hbarc);
double dsigmadb2=0; //GeV0
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
double term=1-0.5*dsigmadb2;
double dsigmadb2tilde=2*(1-term*term);
double result=r*besselK2*besselJ2*dsigmadb2tilde; //fm
return result;
}
Index: trunk/src/Settings.cpp
===================================================================
--- trunk/src/Settings.cpp (revision 533)
+++ trunk/src/Settings.cpp (revision 534)
@@ -1,612 +1,618 @@
//==============================================================================
// Settings.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "Settings.h"
#include "Constants.h"
#include <typeinfo>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <ctype.h>
#include <cstdlib>
#include "TParticlePDG.h"
#include "TError.h"
#define PR(x) cout << #x << " = " << (x) << endl;
TRandom3 Settings::mRandomGenerator;
// initialize static member
double Settings::mQuarkMass[] = {0.0022, 0.0047, 0.096, 1.27, 4.18, 173.21}; // PDG values
Settings::Settings()
{
mPDG = new TDatabasePDG;
mPDG->ReadPDGTable();
//
// Setup name table (map) for nuclei
//
mPeriodicTable[1] = "H";
mPeriodicTable[2] = "He";
mPeriodicTable[3] = "Li";
mPeriodicTable[4] = "Be";
mPeriodicTable[5] = "B";
mPeriodicTable[6] = "C";
mPeriodicTable[7] = "N";
mPeriodicTable[8] = "O";
mPeriodicTable[9] = "F";
mPeriodicTable[10] = "Ne";
mPeriodicTable[11] = "Na";
mPeriodicTable[12] = "Mg";
mPeriodicTable[13] = "Al";
mPeriodicTable[14] = "Si";
mPeriodicTable[15] = "P";
mPeriodicTable[16] = "S";
mPeriodicTable[17] = "Cl";
mPeriodicTable[18] = "Ar";
mPeriodicTable[19] = "K";
mPeriodicTable[20] = "Ca";
mPeriodicTable[21] = "Sc";
mPeriodicTable[22] = "Ti";
mPeriodicTable[23] = "V";
mPeriodicTable[24] = "Cr";
mPeriodicTable[25] = "Mn";
mPeriodicTable[26] = "Fe";
mPeriodicTable[27] = "Co";
mPeriodicTable[28] = "Ni";
mPeriodicTable[29] = "Cu";
mPeriodicTable[30] = "Zn";
mPeriodicTable[31] = "Ga";
mPeriodicTable[32] = "Ge";
mPeriodicTable[33] = "As";
mPeriodicTable[34] = "Se";
mPeriodicTable[35] = "Br";
mPeriodicTable[36] = "Kr";
mPeriodicTable[37] = "Rb";
mPeriodicTable[38] = "Sr";
mPeriodicTable[39] = "Y";
mPeriodicTable[40] = "Zr";
mPeriodicTable[41] = "Nb";
mPeriodicTable[42] = "Mo";
mPeriodicTable[43] = "Tc";
mPeriodicTable[44] = "Ru";
mPeriodicTable[45] = "Rh";
mPeriodicTable[46] = "Pd";
mPeriodicTable[47] = "Ag";
mPeriodicTable[48] = "Cd";
mPeriodicTable[49] = "In";
mPeriodicTable[50] = "Sn";
mPeriodicTable[51] = "Sb";
mPeriodicTable[52] = "Te";
mPeriodicTable[53] = "I";
mPeriodicTable[54] = "Xe";
mPeriodicTable[55] = "Cs";
mPeriodicTable[56] = "Ba";
mPeriodicTable[57] = "La";
mPeriodicTable[58] = "Ce";
mPeriodicTable[59] = "Pr";
mPeriodicTable[60] = "Nd";
mPeriodicTable[61] = "Pm";
mPeriodicTable[62] = "Sm";
mPeriodicTable[63] = "Eu";
mPeriodicTable[64] = "Gd";
mPeriodicTable[65] = "Tb";
mPeriodicTable[66] = "Dy";
mPeriodicTable[67] = "Ho";
mPeriodicTable[68] = "Er";
mPeriodicTable[69] = "Tm";
mPeriodicTable[70] = "Yb";
mPeriodicTable[71] = "Lu";
mPeriodicTable[72] = "Hf";
mPeriodicTable[73] = "Ta";
mPeriodicTable[74] = "W";
mPeriodicTable[75] = "Re";
mPeriodicTable[76] = "Os";
mPeriodicTable[77] = "Ir";
mPeriodicTable[78] = "Pt";
mPeriodicTable[79] = "Au";
mPeriodicTable[80] = "Hg";
mPeriodicTable[81] = "Tl";
mPeriodicTable[82] = "Pb";
mPeriodicTable[83] = "Bi";
mPeriodicTable[84] = "Po";
mPeriodicTable[85] = "At";
mPeriodicTable[86] = "Rn";
mPeriodicTable[87] = "Fr";
mPeriodicTable[88] = "Ra";
mPeriodicTable[89] = "Ac";
mPeriodicTable[90] = "Th";
mPeriodicTable[91] = "Pa";
mPeriodicTable[92] = "U";
mPeriodicTable[93] = "Np";
mPeriodicTable[94] = "Pu";
mPeriodicTable[95] = "Am";
mPeriodicTable[96] = "Cm";
mPeriodicTable[97] = "Bk";
mPeriodicTable[251] = "Cf";
mPeriodicTable[252] = "Es";
mPeriodicTable[100] = "Fm";
mPeriodicTable[258] = "Md";
mPeriodicTable[102] = "No";
mPeriodicTable[103] = "Lr";
mPeriodicTable[261] = "Rf";
mPeriodicTable[105] = "Db";
mPeriodicTable[106] = "Sg";
mPeriodicTable[107] = "Bh";
mPeriodicTable[108] = "Hs";
mPeriodicTable[109] = "Mt";
//
// Register parameters (ptr, name, default)
//
registerParameter(&mUserInt, "userInt", 0);
registerParameter(&mUserDouble, "userDouble", 0.);
registerParameter(&mUserString, "userString", string(""));
registerParameter(&mA, "A", static_cast<unsigned int>(1));
registerParameter(&mVectorMesonId, "vectorMesonId", 443); // J/psi
registerParameter(&mDipoleModelName, "dipoleModel", string("bSat"));
registerParameter(&mDipoleModelParameterSetName, "dipoleModelParameterSet", string("KMW"));
registerParameter(&mTableSetTypeName, "tableSetType", string("total_and_coherent"));
registerParameter(&mQ2min, "Q2min", 1000.); // no limits if max <= min
registerParameter(&mQ2max, "Q2max", 0.);
registerParameter(&mWmin, "Wmin", 1000.);
registerParameter(&mWmax, "Wmax", 0.);
registerParameter(&mXpomMin, "xpomMin", 1e-8);
registerParameter(&mXpomMax, "xpomMax", 1.);
registerParameter(&mbetamin, "betamin", 0.);
registerParameter(&mbetamax, "betamax", 1.);
registerParameter(&mVerbose, "verbose", false);
registerParameter(&mVerboseLevel, "verboseLevel", 0);
registerParameter(&mRootfile, "rootfile", string("sartre.root"));
registerParameter(&mSeed, "seed", static_cast<unsigned int>(time(0)));
registerParameter(&mUPC, "UPC", false);
registerParameter(&mUPCA, "UPCA", static_cast<unsigned int>(1));
registerParameter(&mInclusiveDiffractionMode, "inclusiveDiffractionMode", false);
}
Settings::~Settings()
{
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
delete mRegisteredParameters[k];
}
mRegisteredParameters.clear();
delete mPDG;
}
int Settings::userInt() const {return mUserInt;}
double Settings::userDouble() const {return mUserDouble;}
string Settings::userString() const {return mUserString;}
void Settings::setUserInt(int val) {mUserInt = val;}
void Settings::setUserDouble(double val) {mUserDouble = val;}
void Settings::setUserString(const string& val) {mUserString = val;}
unsigned int Settings::seed() const {return mSeed;}
void Settings::setSeed(unsigned int val)
{
mSeed = val;
mRandomGenerator.SetSeed(mSeed);
gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()
}
bool Settings::readSettingsFromFile(const char *file)
{
if (!file) return false; // nothing to do
mRuncard = string(file);
mLines.clear();
//
// Open file
//
ifstream ifs(file);
if (!ifs) {
cout << "Settings::readSettingsFromFile(): error, cannot open file '"
<< file << "'." << endl;
return false;
}
//
// Read file into vector of strings, skip comments and empty lines
//
while (ifs.good() && !ifs.eof()) {
string line;
getline (ifs, line);
if (ifs.eof() && line.empty()) break;
// empty line
if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) continue;
// if first character is not a letter/digit, then taken to be a comment.
int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
if (!isalnum(line[firstChar])) continue;
mLines.push_back(line);
}
ifs.close(); // done with I/O
//
// Process vector of strings one at a time and use
// it to set registered variables.
//
for (unsigned int i=0; i<mLines.size(); i++) {
// replace '=' by blank to make parsing simpler.
while (mLines[i].find("=") != string::npos) {
int firstEqual = mLines[i].find_first_of("=");
mLines[i].replace(firstEqual, 1, " ");
}
// replace ':' by blank to make parsing simpler.
while (mLines[i].find(":") != string::npos) {
int firstEqual = mLines[i].find_first_of(":");
mLines[i].replace(firstEqual, 1, " ");
}
// get identifier string
istringstream splitLine(mLines[i]);
string name;
splitLine >> name;
// find value string
string valueString;
splitLine >> valueString;
if (!splitLine) {
cout << "Settings::readSettingsFromFile(): error, value of variable '"
<< name.c_str() << "' not recognized." << endl;
}
istringstream modeData(valueString);
//
// Loop over registered variables and see which fits.
// Not particular elegant but does the job and saves
// a lot of programming in derived classes.
//
bool isRegistered = false;
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<vector<double>>)) { // test
SettingsParameter<vector<double>> *var = dynamic_cast<SettingsParameter<vector<double>>*> (mRegisteredParameters[k]);
if (var->name == name) {
var->address->push_back(atof(valueString.c_str())); // first value
while (splitLine.good() && !splitLine.eof()) { // get remaining
string nextValue;
splitLine >> nextValue;
var->address->push_back(atof(nextValue.c_str()));
if (splitLine.eof()) break;
}
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<double>)) {
SettingsParameter<double> *var = dynamic_cast<SettingsParameter<double>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<int>)) {
SettingsParameter<int> *var = dynamic_cast<SettingsParameter<int>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned int>)) {
SettingsParameter<unsigned int> *var = dynamic_cast<SettingsParameter<unsigned int>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned long>)) {
SettingsParameter<unsigned long> *var = dynamic_cast<SettingsParameter<unsigned long>*> (mRegisteredParameters[k]);
if (var->name == name) {
modeData >> (*(var->address));
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<string>)) {
SettingsParameter<string> *var = dynamic_cast<SettingsParameter<string>*> (mRegisteredParameters[k]);
if (var->name == name) {
*(var->address) = valueString;
isRegistered = true;
}
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<bool>)) {
SettingsParameter<bool> *var = dynamic_cast<SettingsParameter<bool>*> (mRegisteredParameters[k]);
if (var->name == name) {
isRegistered = true;
if (valueString == string("true") ||
valueString == string("True") ||
valueString == string("TRUE") ||
valueString == string("on") ||
valueString == string("On") ||
valueString == string("ON") ||
valueString == string("Yes") ||
valueString == string("yes") ||
valueString == string("YES") ||
valueString == string("T") ||
valueString == string("t") ||
valueString == string("1") ) {
*(var->address) = true;
}
else {
*(var->address) = false;
}
}
}
}
if (!isRegistered) {
cout << "Settings::readSettingsFromFile(): error, parameter identifier '" <<
name.c_str() << "' not recognized." << endl;
}
}
//
// Consolidate input (after burner)
//
consolidateCommonSettings();
consolidateSettings(); // overloaded
return true;
}
bool Settings::list(ostream& os)
{
const int fieldWidth = 28;
os << "\nRun Settings:" << endl;
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<double>)) {
SettingsParameter<double> *var = dynamic_cast<SettingsParameter<double>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<int>)) {
SettingsParameter<int> *var = dynamic_cast<SettingsParameter<int>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned int>)) {
SettingsParameter<unsigned int> *var = dynamic_cast<SettingsParameter<unsigned int>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned long>)) {
SettingsParameter<unsigned long> *var = dynamic_cast<SettingsParameter<unsigned long>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<string>)) {
SettingsParameter<string> *var = dynamic_cast<SettingsParameter<string>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << var->address->c_str() << endl;
}
else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<bool>)) {
SettingsParameter<bool> *var = dynamic_cast<SettingsParameter<bool>*> (mRegisteredParameters[k]);
os << setw(fieldWidth) << var->name.c_str() << "\t" << (*(var->address) ? "true" : "false") << endl;
}
}
os << endl;
return true;
}
TParticlePDG* Settings::lookupPDG(int id) const
{
if (mPDG)
return mPDG->GetParticle(id);
else
return 0;
}
string Settings::particleName(int pdgID)
{
string name("unknown");
if (abs(pdgID) < 1000000000) { // particle
if (mPDG) {
TParticlePDG *part = lookupPDG(pdgID);
if (part) name = part->GetName();
}
if (abs(pdgID) == 990) name = "pomeron";
}
else { // nucleus in 10LZZZAAAI PDG format
int id = pdgID;
// int iso = id%10;
id /= 10;
int A = id%1000;
id /= 1000;
int Z = id%1000;
stringstream namestream;
namestream << mPeriodicTable[Z] << "(" << A << ")";
name = namestream.str();
}
return name;
}
void Settings::setVerbose(bool val) {
mVerbose = val;
if (mVerbose && mVerboseLevel == 0) setVerboseLevel(1);
if (!mVerbose && mVerboseLevel != 0) setVerboseLevel(0);
}
bool Settings::verbose() const {return mVerbose;}
void Settings::setVerboseLevel(int val) {
mVerboseLevel = val;
if (!mVerbose && mVerboseLevel != 0) mVerbose = true;
if (mVerbose && mVerboseLevel == 0) mVerbose = false;
//
// Unless verbose level is 5 or higher we suppress the many redundant
// ROOT messages from algorithms since most of their errors are handled
// internally anyway.
//
if (mVerboseLevel < 5) gErrorIgnoreLevel = 5000;
}
int Settings::verboseLevel() const {return mVerboseLevel;}
void Settings::setQ2min(double val) { mQ2min = val;}
double Settings::Q2min() const {return mQ2min;}
double Settings::Qmin() const {return sqrt(mQ2min);}
void Settings::setQ2max(double val) { mQ2max = val;}
double Settings::Q2max() const {return mQ2max;}
double Settings::Qmax() const {return sqrt(mQ2max);}
void Settings::setW2min(double val) { mWmin = sqrt(val);}
void Settings::setWmin(double val) { mWmin = val;}
double Settings::Wmin() const {return mWmin;}
double Settings::W2min() const {return mWmin*mWmin;}
void Settings::setW2max(double val) { mWmax = sqrt(val);}
void Settings::setWmax(double val) { mWmax = val;}
double Settings::Wmax() const {return mWmax;}
double Settings::W2max() const {return mWmax*mWmax;}
void Settings::setXpomMin(double val) {mXpomMin = val;}
void Settings::setXpomMax(double val) {mXpomMax = val;}
double Settings::xpomMin() const {return mXpomMin;}
double Settings::xpomMax() const {return mXpomMax;}
void Settings::setbetamin(double val) {mbetamin = val;}
void Settings::setbetamax(double val) {mbetamax = val;}
double Settings::betamin() const {return mbetamin;}
double Settings::betamax() const {return mbetamax;}
+void Settings::setZmin(double val) {mZmin = val;}
+void Settings::setZmax(double val) {mZmax = val;}
+double Settings::zmin() const {return mZmin;}
+double Settings::zmax() const {return mZmax;}
+
+
int Settings::vectorMesonId() const {return mVectorMesonId;}
void Settings::setVectorMesonId(int val) {mVectorMesonId = val;}
string Settings::dipoleModelName() const {return mDipoleModelName;}
DipoleModelType Settings::dipoleModelType() const {return mDipoleModelType;}
void Settings::setDipoleModelType(DipoleModelType val)
{
mDipoleModelType = val;
if (mDipoleModelType == bSat)
mDipoleModelName = string("bSat");
else if (mDipoleModelType == bNonSat)
mDipoleModelName = string("bNonSat");
else if (mDipoleModelType == bCGC)
mDipoleModelName = string("bCGC");
}
unsigned int Settings::A() const {return mA;}
void Settings::setA(unsigned int val) {mA = val;}
void Settings::setRootfile(const char* val){ mRootfile = val; }
string Settings::rootfile() const { return mRootfile; }
string Settings::dipoleModelParameterSetName() const {return mDipoleModelParameterSetName;}
DipoleModelParameterSet Settings::dipoleModelParameterSet() const {return mDipoleModelParameterSet;}
void Settings::setDipoleModelParameterSet(DipoleModelParameterSet val)
{
mDipoleModelParameterSet = val;
if (mDipoleModelParameterSet == KMW)
mDipoleModelParameterSetName = string("KMW");
else if (mDipoleModelParameterSet == HMPZ)
mDipoleModelParameterSetName = string("HMPZ");
else if (mDipoleModelParameterSet == STU)
mDipoleModelParameterSetName = string("STU");
else if (mDipoleModelParameterSet == CUSTOM)
mDipoleModelParameterSetName = string("CUSTOM");
}
string Settings::tableSetTypeName() const {return mTableSetTypeName;}
TableSetType Settings::tableSetType() const {return mTableSetType;}
void Settings::setTableSetType(TableSetType val)
{
mTableSetType = val;
if (mTableSetType == total_and_coherent)
mTableSetTypeName = string("total_and_coherent");
else if (mTableSetType == coherent_and_incoherent)
mTableSetTypeName = string("coherent_and_incoherent");
}
void Settings::consolidateCommonSettings()
{
//
// Check if verbose levels and flags are consistent
// The verbose flag superseeds the verboseLevel.
//
if (mVerbose && mVerboseLevel == 0) setVerboseLevel(1);
if (mVerboseLevel != 0 && !mVerbose) setVerboseLevel(0);
setVerboseLevel(verboseLevel()); // invoke method no matter what
//
// Set random generator seed
//
mRandomGenerator.SetSeed(mSeed);
gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()
//
// Dipole Model
//
if (mDipoleModelName == string("bSat"))
mDipoleModelType = bSat;
else if (mDipoleModelName == string("bNonSat"))
mDipoleModelType = bNonSat;
else if (mDipoleModelName == string("bCGC"))
mDipoleModelType = bCGC;
else {
cout << "Settings::consolidateCommonSettings(): Error, dipole model '"
<< mDipoleModelName.c_str() << "' is not defined." << endl;
exit(1);
}
//
// Dipole Model Parameter Set
//
if (mDipoleModelParameterSetName == string("KMW"))
mDipoleModelParameterSet = KMW;
else if (mDipoleModelParameterSetName == string("HMPZ"))
mDipoleModelParameterSet = HMPZ;
else if (mDipoleModelParameterSetName == string("STU"))
mDipoleModelParameterSet = STU;
else if (mDipoleModelParameterSetName == string("CUSTOM"))
mDipoleModelParameterSet = CUSTOM;
else {
cout << "Settings::consolidateCommonSettings(): Error, dipole model parameter set'"
<< mDipoleModelParameterSetName.c_str() << "' is not defined." << endl;
exit(1);
}
//
// Table Set Type
//
if (mTableSetTypeName == string("total_and_coherent"))
mTableSetType = total_and_coherent;
else if (mTableSetTypeName == string("coherent_and_incoherent"))
mTableSetType = coherent_and_incoherent;
else {
cout << "Settings::consolidateCommonSettings(): Error, table set type '"
<< mTableSetTypeName.c_str() << "' is not defined." << endl;
exit(1);
}
}
void Settings::setInclusiveDiffractionMode(bool val) { mInclusiveDiffractionMode = val; }
bool Settings::inclusiveDiffractionMode() const { return mInclusiveDiffractionMode; }
bool Settings::exclusiveDiffractionMode() const { return !mInclusiveDiffractionMode; }
void Settings::setUPC(bool val) { mUPC = val; }
bool Settings::UPC() const { return mUPC; }
void Settings::setUPCA(unsigned int val){ mUPCA = val; }
unsigned int Settings::UPCA() const { return mUPCA; }

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:00 PM (23 h, 6 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242362
Default Alt Text
(141 KB)

Event Timeline