Page MenuHomeHEPForge
Files F10881488

No OneTemporary
Actions

View Options

Index: trunk/src/ExclusiveFinalStateGenerator.h
===================================================================
--- trunk/src/ExclusiveFinalStateGenerator.h (revision 490)
+++ trunk/src/ExclusiveFinalStateGenerator.h (revision 491)
@@ -1,41 +1,43 @@
//==============================================================================
// ExclusiveFinalStateGenerator.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 ExclusiveFinalStateGenerator_h
#define ExclusiveFinalStateGenerator_h
#include "FinalStateGenerator.h"
class ExclusiveFinalStateGenerator : public FinalStateGenerator {
public:
ExclusiveFinalStateGenerator();
~ExclusiveFinalStateGenerator();
-
- bool generate(int id, double t, double y, double Q2,
- bool isIncoherent, int A, Event *event);
+
+ bool generate(int id, double t, double y, double Q2, bool isIncoherent, int A, Event *event);
//UPC version:
- bool generate(int id, double t, double xpom,
- bool isIncoherent, int A, Event *event);
+ bool generate(int id, double t, double xpom, bool isIncoherent, int A, Event *event);
+
double xpomMin(double massVM, double t, TLorentzVector hBeam, TLorentzVector eBeam);
-};
+
+ double uiGamma_N(double* var, double* par) const;
+ double gamma_N(unsigned int i, double val);
+};
#endif
Index: trunk/src/Settings.cpp
===================================================================
--- trunk/src/Settings.cpp (revision 490)
+++ trunk/src/Settings.cpp (revision 491)
@@ -1,592 +1,603 @@
//==============================================================================
// 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 <typeinfo>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <ctype.h>
+#include <cstdlib>
#include "TParticlePDG.h"
#include "TError.h"
-#include <cstdlib>
#define PR(x) cout << #x << " = " << (x) << endl;
TRandom3 Settings::mRandomGenerator;
Settings::Settings()
{
- mPDG = new TDatabasePDG;
- mPDG->ReadPDGTable();
+ 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));
}
Settings::~Settings()
{
for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
delete mRegisteredParameters[k];
}
- mRegisteredParameters.clear();
-}
+ 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(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;}
+
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::setUPC(bool val){ mUPC = val; }
bool Settings::UPC() const { return mUPC; }
void Settings::setUPCA(unsigned int val){ mUPCA = val; }
unsigned int Settings::UPCA() const { return mUPCA; }
Index: trunk/src/DipoleModel.h
===================================================================
--- trunk/src/DipoleModel.h (revision 490)
+++ trunk/src/DipoleModel.h (revision 491)
@@ -1,99 +1,102 @@
//==============================================================================
// DipoleModel.h
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// 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$
//==============================================================================
#ifndef DipoleModel_h
#define DipoleModel_h
#include "AlphaStrong.h"
#include "TableGeneratorNucleus.h"
#include "DipoleModelParameters.h"
class TH2F;
class TH1F;
class DipoleModel {
public:
DipoleModel();
virtual ~DipoleModel();
const TableGeneratorNucleus* nucleus() const;
bool configurationExists() const;
virtual void createConfiguration(int)=0;
virtual double dsigmadb2(double, double, double, double)=0;
virtual double bDependence(double);
virtual double bDependence(double, double);
virtual double dsigmadb2ep(double, double, double); // ep
virtual double coherentDsigmadb2(double, double, double) {return 0;}; // eA
virtual void createSigma_ep_LookupTable(double) {/* nothing */};
+ DipoleModelParameters* getParameters(){return mParameters;}
protected:
TableGeneratorNucleus mNucleus;
DipoleModelParameters *mParameters;
AlphaStrong mAs;
bool mConfigurationExists;
bool mIsInitialized;
};
class DipoleModel_bSat : public DipoleModel {
public:
- DipoleModel_bSat();
+ DipoleModel_bSat();
+ DipoleModel_bSat(Settings*);
DipoleModel_bSat(const DipoleModel_bSat&);
DipoleModel_bSat& operator=(const DipoleModel_bSat&);
~DipoleModel_bSat();
void createSigma_ep_LookupTable(double);
void createConfiguration(int);
double dsigmadb2(double, double, double, double);
double bDependence(double, double);
double dsigmadb2ep(double, double, double);
double coherentDsigmadb2(double, double, double);
protected:
TH2F* mBDependence;
private:
double dsigmadb2epForIntegration(double*, double*);
TH1F* mSigma_ep_LookupTable;
};
class DipoleModel_bNonSat : public DipoleModel_bSat {
public:
- DipoleModel_bNonSat();
+ DipoleModel_bNonSat();
+ DipoleModel_bNonSat(Settings*);
~DipoleModel_bNonSat();
double dsigmadb2(double, double, double, double);
double dsigmadb2ep(double, double, double);
double coherentDsigmadb2(double, double, double);
};
class DipoleModel_bCGC : public DipoleModel {
public:
DipoleModel_bCGC();
void createConfiguration(int);
double dsigmadb2(double, double, double, double);
double dsigmadb2ep(double, double, double);
double bDependence(double);
};
#endif
Index: trunk/src/ModeFinderFunctor.h
===================================================================
--- trunk/src/ModeFinderFunctor.h (revision 490)
+++ trunk/src/ModeFinderFunctor.h (revision 491)
@@ -1,70 +1,97 @@
//==============================================================================
// ModeFinderFunctor.h
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// 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: Thomas Ullrich
-// Last update:
+// Author: Tobias Toll
+// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Functor class.
//
// Used by BrentMinimizer1D algotithm to find maxima in cross-section
// in W2 for a given t and Q2.
//==============================================================================
#ifndef ModeFinderFunctor_h
#define ModeFinderFunctor_h
-#include "Math/IFunction.h"
+#include "Math/IFunction.h"
+#include "InclusiveDiffractiveCrossSections.h"
class CrossSection;
class ModeFinderFunctor : public ROOT::Math::IGenFunction {
public:
ModeFinderFunctor();
- ModeFinderFunctor(CrossSection*, double Q2, double vmMass, double tmin, double tmax);
+ ModeFinderFunctor(CrossSection*, double Q2, double vmMass, double tmin, double tmax);
double DoEval(double) const;
- ROOT::Math::IGenFunction* Clone() const;
- void setQ2(double);
+ void setQ2(double);
void setVmMass(double);
void setMaxT(double);
void setMinT(double);
-
-private:
- CrossSection *mCrossSection;
- double mQ2;
+
+ ROOT::Math::IGenFunction* Clone() const;
+
+private:
+ double mQ2;
double mVmMass;
double mMinT;
double mMaxT;
+ CrossSection *mCrossSection;
+};
+
+class InclusiveDiffractionModeFinderFunctor : public ROOT::Math::IGenFunction {
+public:
+ InclusiveDiffractionModeFinderFunctor();
+ InclusiveDiffractionModeFinderFunctor(InclusiveDiffractiveCrossSections* cs, double Q2, double W2, double z, double s, double MX2min, double MX2max, double Q2min, double Q2max, double W2min, double W2max, double ymin, double betamin, double betamax, double mu02);
+ InclusiveDiffractionModeFinderFunctor(InclusiveDiffractiveCrossSections* cs, double Q2, double W2, double z, double MX2min, double MX2max);
+ double DoEval(double) const;
+
+ void setQuarkMass(double);
+
+ ROOT::Math::IGenFunction* Clone() const;
+
+private:
+ InclusiveDiffractiveCrossSections *mIDCrossSection;
+ double mZ;
+ double mQ2;
+ double mMinMX2;
+ double mMaxMX2;
+ double mW2;
+ double mHBeamEnergy;
+ double mEBeamEnergy;
+ double mMq;
+ double mS;
+ double mQ2max, mQ2min, mW2max, mW2min, mYmin, mBetamin, mBetamax, mMu02, mMX2min, mMX2max;
};
class UPCModeFinderFunctor : public ROOT::Math::IGenFunction {
public:
UPCModeFinderFunctor();
UPCModeFinderFunctor(CrossSection*, double vmMass, double hBeamEnergy, double eBeamEnergy);
double DoEval(double) const; // xpom
ROOT::Math::IGenFunction* Clone() const;
private:
CrossSection *mCrossSection;
double mVmMass;
double mHBeamEnergy;
double mEBeamEnergy;
};
#endif
Index: trunk/src/Kinematics.h
===================================================================
--- trunk/src/Kinematics.h (revision 490)
+++ trunk/src/Kinematics.h (revision 491)
@@ -1,105 +1,110 @@
//==============================================================================
// Kinematics.h
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
-// This file is part of Sartre.
+// 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.
+// 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:
+// Last update:
// $Date$
// $Author$
//
//------------------------------------------------------------------------------
//
// Note that Kinematics is only meant for ep and eA diffactive events.
// It is not applicable for UPC kinematics. Use UpcKinematics instead.
//
//==============================================================================
-#ifndef Kinematics_h
-#define Kinematics_h
-#include "Constants.h"
-#include "Settings.h"
-#include "TLorentzVector.h"
-#include "Math/IFunction.h"
-
-class Kinematics {
-public:
- static double xprobe(double Q2, double W2, double vmM);
+#ifndef Kinematics_h
+#define Kinematics_h
+#include "Constants.h"
+#include "Settings.h"
+#include "TLorentzVector.h"
+#include "Math/IFunction.h"
+
+class Kinematics {
+public:
+ static double xprobe(double Q2, double W2, double vmM);
static TLorentzVector electronBeam(double eE, bool upc = false);
- static TLorentzVector hadronBeam(double eH);
- static double s(const TLorentzVector& e, const TLorentzVector& p);
+ static TLorentzVector hadronBeam(double eH);
+ static double s(const TLorentzVector& e, const TLorentzVector& p);
static double s(double eE, double hE, bool upc = false);
- static double x(double Q2, double W2);
- static double y(double Q2, double x, double s);
- static double ymin(double s, double vmM);
- static double ymax(double s, double vmM);
- static double W2(double Q2, double x);
- static double W2(double Q2, double xprobe, double vmM);
- static double W2min(double vmM);
- static double W2max(double s);
- static double Q2min(double y);
- static double Q2max(double s);
- static double xpomeron(double t, double Q2, double W2, double vmM);
- static double tmax(double xpom); // often referred to as tmin implying min |T|
- static double tmax(double t, double Q2, double W2, double vmM);
- static double tmin(double hE); // the smallest t value possible
-
+ static double x(double Q2, double W2);
+ static double y(double Q2, double x, double s);
+ static double ymin(double s, double vmM);
+ static double ymax(double s, double vmM);
+ static double W2(double Q2, double x);
+ static double W2(double Q2, double xprobe, double vmM);
+ static double W2min(double vmM);
+ static double W2max(double s);
+ static double W2max(double s, double Q2, double MX);
+ static double Q2min(double y);
+ static double Q2max(double s);
+ static double xpomeron(double t, double Q2, double W2, double vmM);
+ static double MX2(double Q2, double beta, double t);
+ static double tmax(double xpom); // often referred to as tmin implying min |T|
+ static double tmax(double t, double Q2, double W2, double vmM);
+ static double tmin(double hE); // the smallest t value possible
+ static double betamin(double Q2, double s, double MX, double t);
static double Egamma(double xpom, double t, double vmM, double hBeamEnergy, double eBeamEnergy, double MY2minusM2 = 0);
-
+
static bool validUPC(double hBeamEnergy, double eBeamEnergy, double t,
- double xpom, double vmMass, bool verbose = false); // UPC only
- static bool valid(double s, double t, double Q2, double W2, double vmMass,
- bool useTrueXp = false,
- bool verbose = false);
+ double xpom, double vmMass, bool verbose = false); // UPC only
+ static bool valid(double s, double t, double Q2, double W2, double vmMass,
+ bool useTrueXp = false,
+ bool verbose = false);
+ static bool valid(double s, double beta, double Q2, double W2, double z, double MC,
+ bool useTrueXp = false,
+ bool verbose = false);
static bool error();
-
+
static double xpomMin(double massVM, double t, TLorentzVector hBeam, TLorentzVector eBeam, double MY2minusM2 = 0); // UPC only
-
+
private:
static double xpomMin(double s, double vmM, double t); // UPC only
-
-private:
+
+private:
static bool mError;
-
+
static double mXpomMin;
static double mTold;
static bool mXpomMinIsEvaluated;
-};
+};
-//-------------------------------------------------------------------------------
-//
+//-------------------------------------------------------------------------------
+//
// Helper class needed to find root in
// Kinematics::validUPC()
//
-//-------------------------------------------------------------------------------
+//-------------------------------------------------------------------------------
class LowerXpomeronFormula : public ROOT::Math::IBaseFunctionOneDim
{
-public:
+public:
double DoEval(double) const;
ROOT::Math::IBaseFunctionOneDim* Clone() const;
void calculateValidRange(double&, double&);
-public:
+public:
double mT;
double mVmMass2;
double mXpomMin;
double mMY2;
TLorentzVector mElectronBeam;
TLorentzVector mHadronBeam;
double mMY2minusM2;
};
-#endif
+#endif
Index: trunk/src/FinalStateGenerator.cpp
===================================================================
--- trunk/src/FinalStateGenerator.cpp (revision 490)
+++ trunk/src/FinalStateGenerator.cpp (revision 491)
@@ -1,50 +1,51 @@
//==============================================================================
// FinalStateGenerator.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 "FinalStateGenerator.h"
#include "Event.h"
#include <cmath>
using namespace std;
FinalStateGenerator::FinalStateGenerator()
{
mT = 0;
mQ2 = 0;
mY = 0;
mS = 0;
mMY2 = 0;
- mMassVM = 0;
+ mMX = 0;
+ mMassVM = 0;
mA = 0;
mIsIncoherent = false;
}
FinalStateGenerator::~FinalStateGenerator() {/* no op */}
bool FinalStateGenerator::isValid(TLorentzVector & v) const
{
for (int i=0; i<4; i++) {
if (std::isnan(v[0])) return false;
if (std::isinf(v[0])) return false;
}
return true;
}
Index: trunk/src/Nucleus.h
===================================================================
--- trunk/src/Nucleus.h (revision 490)
+++ trunk/src/Nucleus.h (revision 491)
@@ -1,76 +1,76 @@
//==============================================================================
// Nucleus.h
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
-// This program is free software: you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation.
-// This program is distributed in the hope that it will be useful,
-// but without any warranty; without even the implied warranty of
-// merchantability or fitness for a particular purpose. See the
-// GNU General Public License for more details.
+// 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:
+// Last update:
// $Date$
// $Author$
//==============================================================================
-#ifndef Nucleus_h
-#define Nucleus_h
-#include <string>
+#ifndef Nucleus_h
+#define Nucleus_h
+#include <string>
#include <memory>
using namespace std;
-
-class TH1D;
-
-class Nucleus {
-public:
- Nucleus();
- Nucleus(unsigned int A);
- Nucleus(const Nucleus&);
- virtual ~Nucleus();
-
- Nucleus& operator=(const Nucleus&);
-
- virtual void init(unsigned int A);
-
- double T(double b) const; // b in fm, returns in GeV^2
- double TofProton(double b);
- unsigned int A() const;
- unsigned int Z() const;
- float spin() const; // in hbar
- double radius() const; // in fm
- string name() const;
- int pdgID() const; // id of this nucleus
- int pdgID(int Z, int A) const;
- double nuclearMass() const;
- void normalizationOfT(double eps = 1.e-8); // for checks only
+
+class TH1D;
+
+class Nucleus {
+public:
+ Nucleus();
+ Nucleus(unsigned int A);
+ Nucleus(const Nucleus&);
+ virtual ~Nucleus();
+
+ Nucleus& operator=(const Nucleus&);
+
+ virtual void init(unsigned int A);
+
+ double T(double b) const; // b in fm, returns in GeV^2
+ double TofProton(double b);
+ unsigned int A() const;
+ unsigned int Z() const;
+ float spin() const; // in hbar
+ double radius() const; // in fm
+ string name() const;
+ int pdgID() const; // id of this nucleus
+ int pdgID(int Z, int A) const;
+ double atomicMass() const;
+ void normalizationOfT(double eps = 1.e-8); // for checks only
double rho0() const; //in fm
-protected:
- double rho(double, double);
- double rhoForIntegration(double*, double*);
- double TForIntegration(double*, double*) const;
-
-protected:
- unsigned int mA;
- unsigned int mZ;
- float mSpin;
- double mMass; // nuclear mass in GeV
+protected:
+ double rho(double, double);
+ double rhoForIntegration(double*, double*);
+ double TForIntegration(double*, double*) const;
+
+protected:
+ unsigned int mA;
+ unsigned int mZ;
+ float mSpin;
+ double mMass; // atomic mass in GeV
double mRadius; // fm - Wood-Saxon
- double mSurfaceThickness; // fm - Wood-Saxon
+ double mSurfaceThickness; // fm - Wood-Saxon
double mRho0; // fm^-3 - Wood-Saxon
double mOmega; // Wood-Saxon
double mHulthenA; // for Hulthen distribution
- double mHulthenB;
+ double mHulthenB;
string mName;
unique_ptr<TH1D> mLookupTable; // T lookup table
};
-
-#endif
-
+
+#endif
+
Index: trunk/src/CMakeLists.txt
===================================================================
--- trunk/src/CMakeLists.txt (revision 490)
+++ trunk/src/CMakeLists.txt (revision 491)
@@ -1,163 +1,169 @@
#===============================================================================
# CMakeLists.txt (src)
#
-# Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
-#
+# 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: Thomas Ullrich
# Last update:
# $Date$
# $Author$
#===============================================================================
include(ExternalProject)
cmake_minimum_required (VERSION 3.1)
#
# Compiler flags for release and debug version
#
set(CMAKE_C_FLAGS "-W")
set(CMAKE_C_FLAGS_DEBUG "-g")
set(CMAKE_C_FLAGS_RELEASE "-O")
message("COMPILER = ${CMAKE_CXX_COMPILER_ID}")
-set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -Wpedantic -Wno-long-long")
-if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -pedantic -Wno-long-long")
+if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-potentially-evaluated-expression")
endif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS_DEBUG "-g")
set(CMAKE_CXX_FLAGS_RELEASE "-O")
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
#
# Find external required packages
# (see also FindGSL.cmake and FindROOT.cmke in cmake/modules)
# Herer we need only the root library to create the table
# tools.
#
# GSL
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIR})
# ROOT
find_package(ROOT REQUIRED)
include_directories(${ROOT_INCLUDE_DIR})
set(LIBS ${LIBS} ${ROOT_LIBRARIES})
+# PYTHIA8
+find_package(Pythia8 REQUIRED)
+include_directories(${PYTHIA8_INCLUDE_DIR})
+set(LIBS ${LIBS} ${PYTHIA8_LIBRARIES})
+
#BOOST
if (MULTITHREADED)
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED ON)
find_package(Boost 1.39 COMPONENTS thread REQUIRED)
- if (Boost_FOUND)
+ if(Boost_FOUND)
include_directories(${Boost_INCLUDE_DIR})
endif(Boost_FOUND)
endif (MULTITHREADED)
#
# Include files from sartre package
#
include_directories(${PROJECT_SOURCE_DIR}/src)
include_directories(${PROJECT_SOURCE_DIR}/gemini)
include_directories(${PROJECT_SOURCE_DIR}/cuba)
#
# Cuba is built using the autoconf shipped with it
#
ExternalProject_Add(
cuba
DOWNLOAD_COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/cuba ${PROJECT_BINARY_DIR}/cuba
SOURCE_DIR ${PROJECT_BINARY_DIR}/cuba
PREFIX ${PROJECT_BINARY_DIR}/cuba
CONFIGURE_COMMAND ${PROJECT_BINARY_DIR}/cuba/configure --prefix=${PROJECT_BINARY_DIR}/cuba/build
BUILD_COMMAND make lib
BUILD_IN_SOURCE 1
)
#
# Defines source files for sartre library
#
set(SARTRE_SRC "AlphaStrong.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Amplitudes.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "BreakupProduct.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "CrossSection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DglapEvolution.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModel.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModelParameters.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "EicSmearFormatWriter.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Event.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "EventGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ExclusiveFinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FrangibleNucleus.cpp")
-set(SARTRE_SRC ${SARTRE_SRC} "EventGeneratorSettings.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractionSettings.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractiveCrossSections.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "InclusiveFinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Integrals.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Kinematics.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ModeFinderFunctor.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleon.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "PhotonFlux.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Sartre.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "SartreInclusiveDiffraction.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Settings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Table.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableCollection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorSettings.cpp")
-set(SARTRE_SRC ${SARTRE_SRC} "WaveOverlap.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TwoBodyVectorMesonDecay.cpp")
-set(SARTRE_SRC ${SARTRE_SRC} "EicSmearFormatWriter.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "VectorMesonDecayMass.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "WaveOverlap.cpp")
+
add_library(sartre ${SARTRE_SRC})
#
# Table tools (all stand alone programs)
#
add_executable(tableInspector tableInspector.cpp)
add_executable(tableMerger tableMerger.cpp)
add_executable(tableQuery tableQuery.cpp)
add_executable(tableDumper tableDumper.cpp)
add_executable(tableVarianceMaker tableVarianceMaker.cpp)
-add_executable(tablePriority tablePriority.cpp)
target_link_libraries(tableInspector sartre ${LIBS})
target_link_libraries(tableMerger sartre ${LIBS})
target_link_libraries(tableQuery sartre ${LIBS})
target_link_libraries(tableDumper sartre ${LIBS})
target_link_libraries(tableVarianceMaker sartre ${LIBS})
-target_link_libraries(tablePriority sartre ${LIBS})
add_dependencies(tableInspector sartre)
add_dependencies(tableMerger sartre)
add_dependencies(tableQuery sartre)
add_dependencies(tableDumper sartre)
add_dependencies(tableVarianceMaker sartre)
-add_dependencies(tablePriority sartre)
#
# Install library and include files (make install) within
# the distribution tree. Top level CMakeLists.txt will install
# Sartre in final destination.
#
install(TARGETS sartre DESTINATION sartre/lib)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/lib/libcuba.a" DESTINATION sartre/lib)
FILE(GLOB AllIncludeFiles *.h)
install(FILES ${AllIncludeFiles} DESTINATION sartre/include)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/include/cuba.h" DESTINATION sartre/include)
install(TARGETS tableInspector DESTINATION sartre/bin)
install(TARGETS tableMerger DESTINATION sartre/bin)
install(TARGETS tableQuery DESTINATION sartre/bin)
install(TARGETS tableDumper DESTINATION sartre/bin)
install(TARGETS tableVarianceMaker DESTINATION sartre/bin)
-install(TARGETS tablePriority DESTINATION sartre/bin)
Index: trunk/src/FinalStateGenerator.h
===================================================================
--- trunk/src/FinalStateGenerator.h (revision 490)
+++ trunk/src/FinalStateGenerator.h (revision 491)
@@ -1,59 +1,54 @@
//==============================================================================
// FinalStateGenerator.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 FinalStateGenerator_h
#define FinalStateGenerator_h
#include "TLorentzVector.h"
class Event;
class FinalStateGenerator {
public:
FinalStateGenerator();
virtual ~FinalStateGenerator();
-
- virtual bool generate(int id, double t, double y, double Q2,
- bool isIncoherent, int A, Event *event) = 0;
-
- virtual bool generate(int id, double t, double xpom,
- bool isIncoherent, int A, Event *event) = 0;
-
- bool isValid(TLorentzVector &) const;
+
+ bool isValid(TLorentzVector &) const;
protected:
double mT;
double mQ2;
double mY;
double mS;
- double mXp; //UPC
- double mEgam; //UPC
+ double mXp; // UPC
+ double mEgam; // UPC
+ double mMX; // inclusive diffraction
double mMY2;
double mMassVM;
double mA;
bool mIsIncoherent;
TLorentzVector mElectronBeam;
TLorentzVector mHadronBeam;
};
#endif
Index: trunk/src/PhotonFlux.cpp
===================================================================
--- trunk/src/PhotonFlux.cpp (revision 490)
+++ trunk/src/PhotonFlux.cpp (revision 491)
@@ -1,273 +1,273 @@
//==============================================================================
// PhotonFlux.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, Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "PhotonFlux.h"
#include "Constants.h"
#include "Kinematics.h"
#include "EventGeneratorSettings.h"
#include "Nucleus.h"
#include <cmath>
#include "TF1.h"
#include "Math/Functor.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#define PR(x) cout << #x << " = " << (x) << endl;
PhotonFlux::PhotonFlux()
{
mS = 0;
mIsUPC = false;
mSettings = EventGeneratorSettings::instance();
if (mSettings->UPC()){
mNucleusUPC = new Nucleus(mSettings->UPCA());
mNucleus = new Nucleus(mSettings->A());
if (mSettings->UPCA()>1 && mSettings->A()>1) {
cout << "PhotonFlux::PhotonFlux(): Calculating TAA table ... ";
setupTAAlookupTable();
cout << "done." << endl;
}
mEBeamEnergy = mSettings->electronBeamEnergy(); // the beam that shines
mIsUPC = true;
}
mSIsSet = false;
}
PhotonFlux::PhotonFlux(double s)
{
mS = s;
mIsUPC = false;
mSettings = EventGeneratorSettings::instance();
if (mSettings->UPC()) {
mNucleusUPC = new Nucleus(mSettings->UPCA());
mNucleus = new Nucleus(mSettings->A());
if (mSettings->UPCA()>1 && mSettings->A()>1){
cout << "PhotonFlux::PhotonFlux(): Calculating TAA table ... ";
setupTAAlookupTable();
cout << "done" << endl;
}
mEBeamEnergy = mSettings->electronBeamEnergy(); //The beam that shines
mIsUPC = true;
cout<<"PhotonFlux::PhotonFlux(): Warning Sartre is not yet ready for UPC..." << endl;
}
mSIsSet = true;
}
void PhotonFlux::setS(double s)
{
mS = s;
mSIsSet = true;
}
double PhotonFlux::operator()(double Q2, double W2, GammaPolarization p) const
{
if (!mSIsSet) cout<<"Warning in PhotonFlux::operator() s is not set!"<<endl;
if (p == transverse)
return fluxTransverse(Q2, W2);
else
return fluxLongitudinal(Q2, W2);
}
double PhotonFlux::operator()(double Egamma) const
{
if (!mSIsSet) cout << "PhotonFlux::operator(): Warning, s is not set!" << endl;
return nuclearPhotonFlux(Egamma);
}
double PhotonFlux::fluxTransverse(double Q2, double W2) const
{
//
// Transverse photon flux
//
double y = Kinematics::y(Q2, Kinematics::x(Q2, W2), mS);
if (y<0 || y>1 || Kinematics::error()) return 0;
if (Q2 > Kinematics::Q2max(mS)) return 0;
double result = alpha_em/(2*M_PI*Q2*mS*y);
result *= (1 + (1-y)*(1-y)) - (2*(1-y)*Kinematics::Q2min(y)/Q2);
return result;
}
double PhotonFlux::fluxLongitudinal(double Q2, double W2) const
{
//
// Longitudinal photon flux
//
double y = Kinematics::y(Q2, Kinematics::x(Q2, W2), mS);
if (y<0 || y>1 || Kinematics::error()) return 0;
if (Q2 > Kinematics::Q2max(mS)) return 0;
double result = alpha_em/(2*M_PI*Q2*mS*y);
result *= 2*(1-y);
return result;
}
double PhotonFlux::nuclearPhotonFlux(double Egamma) const
{
if (Egamma < 0) {
if (mSettings->verbose() && mSettings->verboseLevel() > 4)
cout << "PhotonFlux::nuclearPhotonFlux(): Warning, negative Egamma as argument.\n"
<< " Return 0 for photon flux." << endl;
return 0;
}
double bmin = 0.;
double bmax = 1e10; //Could be arbitrary large (TF1 doesn't care)
TF1 fluxFunction("fluxFunction", this, &PhotonFlux::unintegratedNuclearPhotonFlux, bmin, bmax, 1);
fluxFunction.SetParameter(0, Egamma);
ROOT::Math::WrappedTF1 wrappedFluxFunction(fluxFunction);
ROOT::Math::GaussIntegrator fluxIntegrator;
fluxIntegrator.SetFunction(wrappedFluxFunction);
fluxIntegrator.SetAbsTolerance(0.);
fluxIntegrator.SetRelTolerance(1e-5);
return fluxIntegrator.IntegralUp(bmin)/hbarc; //dN/dEgamma GeV^-1
}
double PhotonFlux::unintegratedNuclearPhotonFlux(double* var, double* par) const
{
//
// https://arxiv.org/abs/1607.03838v1
//
double b = var[0]; //fm
double eGamma = par[0]; //GeV
int Apom = mNucleus->A();
int AUPC = mNucleusUPC->A();
double Z = mNucleusUPC->Z();
- double ebeamMass = mNucleusUPC->nuclearMass()/mNucleusUPC->A(); //Mass per nucleon //GeV
+ double ebeamMass = mNucleusUPC->atomicMass()/mNucleusUPC->A(); //Mass per nucleon //GeV
double beamLorentzGamma = mEBeamEnergy/ebeamMass;
double beamLorentzGamma2 = beamLorentzGamma*beamLorentzGamma;
double xi = eGamma*b/hbarc/beamLorentzGamma; //GeV0
double K0 = TMath::BesselK0(xi);
double K1 = TMath::BesselK1(xi);
double prefactor = Z*Z*alpha_em/(M_PI*M_PI)*eGamma/beamLorentzGamma2; //GeV1
double N = prefactor*(K1*K1+K0*K0/beamLorentzGamma2); //GeV1
//Spectrum is calculated under the assumption that there is no
//hadronic interaction between the beam particles.
//Therefore it has to be multiplied by the probability for this:
double Pnohad = 1;
if (Apom>1 && AUPC>1) { //nucleus-nucleus interaction
Pnohad = exp(-sigma_nn(mS)*mTAA_of_b->Interpolate(b)); //GeV0
}
else if (Apom == 1 && AUPC > 1) { //proton-nucleus
Pnohad=exp(-sigma_nn(mS)*mNucleusUPC->T(b));
}
else if (Apom > 1 && AUPC == 1) { //nucleus-proton
Pnohad = exp(-sigma_nn(mS)*mNucleus->T(b));
}
else { //proton-proton
double b02 = 19.8; //GeV^-2
double scamp = 1-exp(-b*b/hbarc2/2./b02);
Pnohad = scamp*scamp;
}
double result = 2*M_PI*b/hbarc*N*Pnohad; //GeV0
return result;
}
void PhotonFlux::setupTAAlookupTable()
{
double rbhigh = upperIntegrationLimit*(mNucleus->radius()+mNucleusUPC->radius());
mTAA_of_b = new TH1D("mTAA_of_b", "mTAA_of_b", 1000, 0, rbhigh);
for (unsigned int i=1; i<=1000; i++) {
double b = mTAA_of_b->GetBinCenter(i);
mTAA_of_b->SetBinContent(i, TAA(b));
}
}
double PhotonFlux::sigma_nn(double s) const {
//
// For parameters see
// http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20080014212.pdf table 2
//
// Following KNSGB: http://arxiv.org/abs/1607.03838v1 (eq.6)
//
double Z = 33.73; //mb
double B = 0.2838; //mb
double s0 = 1.;//GeV2
double Y1 = 13.67; //mb
double s1 = 1.; //GeV2
double eta1 = 0.412;
double Y2 = 7.77; //mb
double eta2 = 0.5626;
double result = Z+B*log(s/s0)*log(s/s0)+Y1*pow(s1/s, eta1)-Y2*pow(s1/s, eta2); //mb = 1e-3b = 1e-1 fm2
result *= 1e-1; //fm2 (1 barn = 1e2 fm2)
return result/hbarc2; //GeV^-2
}
double PhotonFlux::TAA(double b)
{
//
// Integral over dr^2=2*pi*r*dr
//
double rbhigh=upperIntegrationLimit*(mNucleus->radius()+mNucleusUPC->radius());
double lo[2] = {0., 0.}; //r, theta
double hi[2] = {rbhigh, 2*M_PI};
ROOT::Math::Functor wf(this, &PhotonFlux::TAAForIntegration, 2);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig.SetFunction(wf);
ig.SetAbsTolerance(0.);
ig.SetRelTolerance(1e-4);
mB=b;
double result = ig.Integral(lo, hi); //GeV^3*fm
result /= hbarc; //GeV^2
return result;
}
double PhotonFlux::TAAForIntegration(const double* var) const
{
double r = var[0]; //fm
double phi = var[1];
double b = mB; //fm
//Put A1 in the origin, and A2 on the y-axis at distance b:
double arg1 = r;
double arg2 = sqrt(r*r+b*b-2*b*r*sin(phi)); //fm
double T1 = mNucleusUPC->T(arg1); //GeV^2
double T2 = mNucleus->T(arg2); //GeV^2
return r/hbarc*T1*T2; //GeV^3
}
Index: trunk/src/Nucleus.cpp
===================================================================
--- trunk/src/Nucleus.cpp (revision 490)
+++ trunk/src/Nucleus.cpp (revision 491)
@@ -1,355 +1,353 @@
//==============================================================================
// Nucleus.cpp
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
-// This program is free software: you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation.
-// This program is distributed in the hope that it will be useful,
-// but without any warranty; without even the implied warranty of
-// merchantability or fitness for a particular purpose. See the
-// GNU General Public License for more details.
+// 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:
+// Last update:
// $Date$
// $Author$
//==============================================================================
-#include "Nucleus.h"
-#include "Constants.h"
-#include "TF1.h"
-#include "TH1.h"
+#include "Nucleus.h"
+#include "Constants.h"
+#include "TF1.h"
+#include "TH1.h"
#include "TH1D.h"
#include <cmath>
-#include <cstdlib>
-
-#define PR(x) cout << #x << " = " << (x) << endl;
+#include <cstdlib>
+#include "Math/WrappedTF1.h"
+#include "Math/GaussIntegrator.h"
+
+
+#define PR(x) cout << #x << " = " << (x) << endl;
Nucleus::Nucleus()
-{
- mA = mZ = 0;
- mRadius = 0;
- mSpin = 0;
- mSurfaceThickness = 0;
- mRho0 = 0;
- mOmega = 0;
+{
+ mA = mZ = 0;
+ mRadius = 0;
+ mSpin = 0;
+ mSurfaceThickness = 0;
+ mRho0 = 0;
+ mOmega = 0;
mMass = 0;
mHulthenA = 0;
mHulthenB = 0;
}
-
-Nucleus::Nucleus(unsigned int A)
+
+Nucleus::Nucleus(unsigned int A)
{
mA = mZ = 0;
mRadius = 0;
mSpin = 0;
mSurfaceThickness = 0;
mRho0 = 0;
mOmega = 0;
mMass = 0;
mHulthenA = 0;
mHulthenB = 0;
Nucleus::init(A);
-}
-
-Nucleus& Nucleus::operator=(const Nucleus& n)
-{
+}
+
+Nucleus& Nucleus::operator=(const Nucleus& n)
+{
if (this != &n) {
mA = n.mA;
- mZ = n.mZ;
- mRadius = n.mRadius;
- mSpin = n.mSpin;
- mSurfaceThickness = n.mSurfaceThickness;
- mRho0 = n.mRho0;
- mOmega = n.mOmega;
- mMass = n.mMass;
- mName = n.mName;
+ mZ = n.mZ;
+ mRadius = n.mRadius;
+ mSpin = n.mSpin;
+ mSurfaceThickness = n.mSurfaceThickness;
+ mRho0 = n.mRho0;
+ mOmega = n.mOmega;
+ mMass = n.mMass;
+ mName = n.mName;
mHulthenA = n.mHulthenA;
mHulthenB = n.mHulthenB;
mLookupTable = unique_ptr<TH1D>(new TH1D(*(n.mLookupTable)));
- mLookupTable->SetDirectory(0);
- }
- return *this;
-}
-
-Nucleus::Nucleus(const Nucleus& n)
-{
+ mLookupTable->SetDirectory(0);
+ }
+ return *this;
+}
+
+Nucleus::Nucleus(const Nucleus& n)
+{
mA = n.mA;
- mZ = n.mZ;
- mRadius = n.mRadius;
- mSpin = n.mSpin;
- mSurfaceThickness = n.mSurfaceThickness;
- mRho0 = n.mRho0;
- mOmega = n.mOmega;
- mMass = n.mMass;
- mName = n.mName;
+ mZ = n.mZ;
+ mRadius = n.mRadius;
+ mSpin = n.mSpin;
+ mSurfaceThickness = n.mSurfaceThickness;
+ mRho0 = n.mRho0;
+ mOmega = n.mOmega;
+ mMass = n.mMass;
+ mName = n.mName;
mHulthenA = n.mHulthenA;
mHulthenB = n.mHulthenB;
mLookupTable = unique_ptr<TH1D>(new TH1D(*(n.mLookupTable)));
mLookupTable->SetDirectory(0);
-}
-
+}
+
Nucleus::~Nucleus()
{
//
// This is a temorary fix of a problem in the code that is not understood.
// In the table generation code after all is calculated the ~TH1 of the
// look up table histogram causes a crash. This appears (?) to be related
// to issues in GSL and ROOT. Releasing the pointer causes a memory leak
// but this is not a big problem since it happens at the end of the program.
//
mLookupTable.release();
}
+
+void Nucleus::init(unsigned int A)
+{
+ mA = A;
-void Nucleus::init(unsigned int A)
-{
- mA = A;
-
- //
- // Set the parameters of the nucleus
- // Woods-Saxon parameter are from Ramona Vogt's paper on nuclear geometry (table 1).
- // The parametrization is a 3 parameter fermi model (or 2 parameter with omega=0).
- // One exception is the deuteron where a Hulthen function is used (nucl-ex/0701025v1).
+ //
+ // Set the parameters of the nucleus
+ // Woods-Saxon parameter are from Ramona Vogt's paper on nuclear geometry (table 1)
+ // One exception is the deuteron were a Hulthen function is used (nucl-ex/0701025v1).
// All version are defined such that Integral(d3r rho(r)) = A
- // Masses from https://wwwndc.jaea.go.jp/NuC/
//
- const double unifiedAtomicMass = 0.931494013;
- switch (mA) {
- case 208: // Pb
- mZ = 82;
- mRadius = 6.624;
- mSurfaceThickness = 0.549;
- mOmega = 0;
- mRho0 = 0.160;
- mName = "Pb";
- mMass = 207.976651189*unifiedAtomicMass-mZ*electronMass;
- mSpin = 0;
- break;
- case 197: // Au
- mZ = 79;
- mRadius = 6.38;
- mSurfaceThickness = 0.535;
- mOmega = 0;
- mRho0 = 0.1693;
- mName = "Au";
- mMass = 196.966568812*unifiedAtomicMass-mZ*electronMass;
- mSpin = 3./2.;
- break;
- case 110: // Cd
- mZ = 48;
- mRadius = 5.33;
- mSurfaceThickness = 0.535;
- mOmega = 0;
- mRho0 = 0.1577;
- mName = "Cd";
- mMass = 109.903004927*unifiedAtomicMass-mZ*electronMass;
- mSpin = 0;
+ const double unifiedAtomicMass = 0.931494013;
+ switch (mA) {
+ case 208: // Pb
+ mZ = 82;
+ mRadius = 6.624;
+ mSurfaceThickness = 0.549;
+ mOmega = 0;
+ mRho0 = 0.160;
+ mName = "Pb";
+ mMass = 207.2*unifiedAtomicMass;
+ mSpin = 1./2.;
+ break;
+ case 197: // Au
+ mZ = 79;
+ mRadius = 6.38;
+ mSurfaceThickness = 0.535;
+ mOmega = 0;
+ mRho0 = 0.1693;
+ mName = "Au";
+ mMass = 196.96655*unifiedAtomicMass;
+ mSpin = 3./2.;
+ break;
+ case 110: // Cd
+ mZ = 48;
+ mRadius = 5.33;
+ mSurfaceThickness = 0.535;
+ mOmega = 0;
+ mRho0 = 0.1577;
+ mName = "Cd";
+ mMass = 112.411*unifiedAtomicMass;
+ mSpin = 1./2.;
+ break;
+ case 63: // Cu
+ mZ = 29;
+ mRadius = 4.214;
+ mSurfaceThickness = 0.586;
+ mOmega = 0;
+ mRho0 = 0.1701;
+ mName = "Cu";
+ mMass = 63.546*unifiedAtomicMass;
+ mSpin = 3./2.;
+ break;
+ case 40: // Ca
+ mZ = 20;
+ mRadius = 3.766;
+ mSurfaceThickness = 0.586;
+ mOmega = -0.161;
+ mRho0 = 0.1699;
+ mName = "Ca";
+ mMass= 40.078*unifiedAtomicMass;
+ mSpin = 7./2.;
+ break;
+ case 27: // Al
+ mZ = 13;
+ mRadius = 3.07;
+ mSurfaceThickness = 0.519;
+ mOmega = 0;
+ mRho0 = 0.1739;
+ mName = "Al";
+ mMass = 26.981538*unifiedAtomicMass;
+ mSpin = 5./2.;
break;
- case 90: // Zr
- mZ = 40;
- mRadius = 4.9736093; // from R. Vogt 1.19*A^(1/3)-1.61*A^(-1/3)
- mSurfaceThickness = 0.54; // from R. Vogt
- mOmega = 0; // from R. Vogt
- mRho0 = 0.15643742;
- mName = "Zr";
- mMass = 89.904696939*unifiedAtomicMass-mZ*electronMass;
- mSpin = 0;
+ case 16: // O
+ mZ = 8;
+ mRadius = 2.608;
+ mSurfaceThickness = 0.513;
+ mOmega = -0.051;
+ mRho0 = 0.1654;
+ mName = "O";
+ mMass = 15.9994*unifiedAtomicMass;
+ mSpin = 5./2.;
break;
- case 63: // Cu
- mZ = 29;
- mRadius = 4.214;
- mSurfaceThickness = 0.586;
- mOmega = 0;
- mRho0 = 0.1701;
- mName = "Cu";
- mMass = 62.929597770*unifiedAtomicMass-mZ*electronMass;
- mSpin = 3./2.;
- break;
- case 40: // Ca
- mZ = 20;
- mRadius = 3.766;
- mSurfaceThickness = 0.586;
- mOmega = -0.161;
- mRho0 = 0.1699;
- mName = "Ca";
- mMass= 39.962590863*unifiedAtomicMass-mZ*electronMass;
- mSpin = 0;
- break;
- case 27: // Al
- mZ = 13;
- mRadius = 3.07;
- mSurfaceThickness = 0.519;
- mOmega = 0;
- mRho0 = 0.1739;
- mName = "Al";
- mMass = 26.981538578*unifiedAtomicMass-mZ*electronMass;
- mSpin = 5./2.;
- break;
- case 16: // O
- mZ = 8;
- mRadius = 2.608;
- mSurfaceThickness = 0.513;
- mOmega = -0.051;
- mRho0 = 0.1654;
- mName = "O";
- mMass = 15.99491461957*unifiedAtomicMass-mZ*electronMass;
- mSpin = 0;
- break;
case 2: // D
mZ = 1;
mHulthenA = 0.456;
mHulthenB = 2.36;
mRho0 = 0.39946312;
mName = "D";
- mMass = 2.01410177812*unifiedAtomicMass-mZ*electronMass;
+ mMass = 2.014102*unifiedAtomicMass;
mSpin = 1;
mRadius = 2.116; // not used for density function - taken from Jager et al.
break;
case 1:
- // When generating this nucleus, only one nucleon is put at (0, 0, 0)
- mZ = 1;
- mRadius = 1.;
- mMass = protonMass; // 1.00794*unifiedAtomicMass
- mName = "p";
- //The following 3 are dummies, needed for radial distribution etc.:
- mSurfaceThickness = 0.513;
- mOmega = -0.051;
- mRho0 = 0.1654;
- mSpin = 1./2.;
- break;
- default:
- cout << "Nucleus::init(): Error, cannot handle A=" << mA << ", no parameters defined for this mass number." << endl;
- exit(1);
- break;
- }
-
- //
- // Generate lookup table for overlap integral T
- // (b and z in fm)
- //
- double range = 3*mRadius;
+ // When generating this nucleus, only one nucleon is put at (0, 0, 0)
+ mZ = 1;
+ mRadius = 1.;
+ mMass = protonMass; // 1.00794*unifiedAtomicMass
+ mName = "p";
+ //The following 3 are dummies, needed for radial distribution etc.:
+ mSurfaceThickness = 0.513;
+ mOmega = -0.051;
+ mRho0 = 0.1654;
+ mSpin = 1./2.;
+ break;
+ default:
+ cout << "Nucleus::init(): Error, cannot handle A=" << mA << ", no parameters defined for this mass number." << endl;
+ exit(1);
+ break;
+ }
+
+ //
+ // Generate lookup table for overlap integral T
+ // (b and z in fm)
+ //
+ double range = 3*mRadius;
TF1* funcToIntegrate = new TF1("funcToIntegrate",this, &Nucleus::rhoForIntegration, -range, range, 1);
int nbins = 5000; // size of lookup table
mLookupTable = unique_ptr<TH1D>(new TH1D("mLookupTable","T lookup table", nbins, 0, range));
mLookupTable->SetDirectory(0);
+ ROOT::Math::WrappedTF1* WFlookup=new ROOT::Math::WrappedTF1(*funcToIntegrate);
+ ROOT::Math::GaussIntegrator GIlookup;
+ GIlookup.SetFunction(*WFlookup);
+ GIlookup.SetRelTolerance(1e-8);
+
for (int i=1; i<=nbins; i++) {
- double b = mLookupTable->GetBinCenter(i);
- funcToIntegrate->SetNpx(1000);
- funcToIntegrate->SetParameter(0, b);
- double res = funcToIntegrate->Integral(-range, range, 1e-8);
- mLookupTable->SetBinContent(i, res);
- }
+ double b = mLookupTable->GetBinCenter(i);
+ funcToIntegrate->SetNpx(1000);
+ funcToIntegrate->SetParameter(0, b);
+ // double res = funcToIntegrate->Integral(-range, range, 1e-8);
+ double res = GIlookup.Integral(-range, range);
+ mLookupTable->SetBinContent(i, res);
+ }
delete funcToIntegrate;
+ delete WFlookup;
}
-
-unsigned int Nucleus::A() const { return mA; }
-
-unsigned int Nucleus::Z() const { return mZ; }
-
-float Nucleus::spin() const { return mSpin; }
-
-double Nucleus::nuclearMass() const { return mMass; }
-
-string Nucleus::name() const { return mName; }
-
-double Nucleus::radius() const { return mRadius; }
-
-double Nucleus::rho(double b, double z) // returns value in units of fm^-3
-{
- // b and z in fm
+
+unsigned int Nucleus::A() const { return mA; }
+
+unsigned int Nucleus::Z() const { return mZ; }
+
+float Nucleus::spin() const { return mSpin; }
+
+double Nucleus::atomicMass() const { return mMass; }
+
+string Nucleus::name() const { return mName; }
+
+double Nucleus::radius() const { return mRadius; }
+
+double Nucleus::rho(double b, double z) // returns value in units of fm^-3
+{
+ // b and z in fm
double r = sqrt(b*b+z*z);
if (mA == 2) {
double term = (exp(-mHulthenA*r)/r) - (exp(-mHulthenB*r)/r); // Hulthen
return mRho0*term*term;
}
else {
return mRho0*(1+mOmega*(r/mRadius)*(r/mRadius))/(1+exp((r-mRadius)/mSurfaceThickness)); // 3pF
}
-}
-
-double Nucleus::rhoForIntegration(double *x, double* par)
-{
- // b and z in fm
- double b = par[0];
- double z = *x;
- return rho(b, z);
-}
-
-double Nucleus::T(double b) const
-{
- //
- // Returns overlap integral
- // b in fm, overlap function T in GeV^2
- // Returns 0 if b exceeds table range
- //
- if (mA == 0) {
- cout << "Nucleus::T(): Error, instance is not initialized - cannot calculate overlap integral." << endl;
- return 0;
- }
+}
+
+double Nucleus::rhoForIntegration(double *x, double* par)
+{
+ // b and z in fm
+ double b = par[0];
+ double z = *x;
+ return rho(b, z);
+}
+
+double Nucleus::T(double b) const
+{
+ //
+ // Returns overlap integral
+ // b in fm, overlap function T in GeV^2
+ // Returns 0 if b exceeds table range
+ //
+ if (mA == 0) {
+ cout << "Nucleus::T(): Error, instance is not initialized - cannot calculate overlap integral." << endl;
+ return 0;
+ }
int bin = mLookupTable->FindBin(b);
double res = mLookupTable->GetBinContent(bin);
return res*hbarc2;
-}
-
-double Nucleus::TForIntegration(double *x, double*) const
-{
- double b = x[0];
- return 2*b*M_PI*T(b)/hbarc2;
-}
-
+}
+
+double Nucleus::TForIntegration(double *x, double*) const
+{
+ double b = x[0];
+ return 2*b*M_PI*T(b)/hbarc2;
+}
+
void Nucleus::normalizationOfT(double eps)
-{
- //
- // This is for internal checks only.
- // Normalization of rho/T is such that fully integrated
- // function should yield A.
- //
- double range = 3*mRadius;
- TF1* func = new TF1("func",this, &Nucleus::TForIntegration, 0, range, 0);
+{
+ //
+ // This is for internal checks only.
+ // Normalization of rho/T is such that fully integrated
+ // function should yield A.
+ //
+ double range = 3*mRadius;
+ TF1* func = new TF1("func",this, &Nucleus::TForIntegration, 0, range, 0);
double res = func->Integral(0, range, eps);
cout << "Normalization of T: = " << res << endl;
delete func;
-}
-
+}
+
double Nucleus::rho0() const { return mRho0; }
-double Nucleus::TofProton(double b)
-{
- //
- // Gaussian shape for proton
- // Units:
- // b in fm
+double Nucleus::TofProton(double b)
+{
+ //
+ // Gaussian shape for proton
+ // Units:
+ // b in fm
//
- const double BG = 4; // GeV^-2
- double arg = b*b/(2*BG);
- arg /= hbarc2;
- return 1/(2*M_PI*BG) * exp(-arg);
-}
-
-int Nucleus::pdgID(int Z, int A) const
-{
- // PDG for nuclei 10LZZZAAAI
- // L = number of strange quark (for hypernuclei)
- // I = isomer level, with I = 0 corresponding to the
- // ground state and I > 0 to excitations,
- if (Z==1 && A==1)
- return 2212;
- else if (Z==0 && A==1)
- return 2112;
- else if (A > 1)
- return 1000000000 + 10000*Z + 10*A;
- else
- return 0;
-}
-
-int Nucleus::pdgID() const
-{
- return pdgID(mZ, mA);
-}
-
+ const double BG = 4; // GeV^-2
+ double arg = b*b/(2*BG);
+ arg /= hbarc2;
+ return 1/(2*M_PI*BG) * exp(-arg);
+}
+
+int Nucleus::pdgID(int Z, int A) const
+{
+ // PDG for nuclei 10LZZZAAAI
+ // L = number of strange quark (for hypernuclei)
+ // I = isomer level, with I = 0 corresponding to the
+ // ground state and I > 0 to excitations,
+ if (Z==1 && A==1)
+ return 2212;
+ else if (Z==0 && A==1)
+ return 2112;
+ else if (A > 1)
+ return 1000000000 + 10000*Z + 10*A;
+ else
+ return 0;
+}
+
+int Nucleus::pdgID() const
+{
+ return pdgID(mZ, mA);
+}
+
Index: trunk/src/FrangibleNucleus.cpp
===================================================================
--- trunk/src/FrangibleNucleus.cpp (revision 490)
+++ trunk/src/FrangibleNucleus.cpp (revision 491)
@@ -1,204 +1,204 @@
//==============================================================================
// FrangibleNucleus.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 "FrangibleNucleus.h"
#include "CNucleus.h"
#include "Constants.h"
#include "EventGeneratorSettings.h"
#include "TRandom3.h"
#define PR(x) cout << #x << " = " << (x) << endl;
FrangibleNucleus::FrangibleNucleus()
{
mGeminiNucleus = nullptr;
mExcitationEnergy = 0;
}
FrangibleNucleus& FrangibleNucleus::operator=(const FrangibleNucleus& gn)
{
if (this != &gn) {
delete mGeminiNucleus;
mProducts.clear();
Nucleus::operator=(gn);
mExcitationEnergy = gn.mExcitationEnergy;
copy(gn.mProducts.begin(), gn.mProducts.end(), mProducts.begin());
if (gn.mGeminiNucleus)
mGeminiNucleus = new CNucleus(gn.mZ, gn.mA);
else
mGeminiNucleus = nullptr;
}
return *this;
}
FrangibleNucleus::FrangibleNucleus(const FrangibleNucleus& gn) : Nucleus(gn)
{
mExcitationEnergy = gn.mExcitationEnergy;
copy(gn.mProducts.begin(), gn.mProducts.end(), mProducts.begin());
if (gn.mGeminiNucleus)
mGeminiNucleus = new CNucleus(gn.mZ, gn.mA);
else
- mGeminiNucleus = 0;
-}
+ mGeminiNucleus = nullptr;
+}
void FrangibleNucleus::init(unsigned int A)
{
Nucleus::init(A);
if (mGeminiNucleus) delete mGeminiNucleus;
mGeminiNucleus = nullptr;
}
void FrangibleNucleus::init(unsigned int A, bool enableBreakup)
{
Nucleus::init(A);
if (mGeminiNucleus) delete mGeminiNucleus;
//
// Init Gemini.
//
if (enableBreakup)
mGeminiNucleus = new CNucleus(mZ, mA);
else
mGeminiNucleus = nullptr;
mExcitationEnergy = 0;
}
FrangibleNucleus::FrangibleNucleus(unsigned int A, bool enableBreakup) : Nucleus(A)
{
//
// Init Gemini.
// Only if requested, not needed otherwise.
//
if (enableBreakup)
mGeminiNucleus = new CNucleus(mZ, mA);
else
mGeminiNucleus = nullptr;
mExcitationEnergy = 0;
}
FrangibleNucleus::~FrangibleNucleus()
{
delete mGeminiNucleus;
}
void FrangibleNucleus::resetBreakup()
{
mExcitationEnergy = 0;
mProducts.clear();
}
int FrangibleNucleus::breakup(const TLorentzVector& dissSystem)
{
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
//
// Estimate excitation energy
// Note that the dissSystem is given in units of 'per nucleon'.
// Hence we have to multiply the excitation energy with A.
//
mExcitationEnergy = (dissSystem.M()-protonMass)*mA;
double Ex = mExcitationEnergy;
double maxEx = settings->maxNuclearExcitationEnergy();
if (Ex > maxEx) {
if (settings->verboseLevel() > 1)
cout << "FrangibleNucleus::breakup(): Actual excitation energy (" << Ex
<< ") exceeded upper limit (" << maxEx << "). Reset to maximum value." << endl;
Ex = maxEx;
}
Ex *= 1000; // Gemini uses the total energy in MeV
//
// Setup excited nucleus
//
// Pass excitation energy and spin to Gemini
mGeminiNucleus->setCompoundNucleus(Ex,mSpin); //specify the excitation energy and spin
mGeminiNucleus->setVelocityCartesian(); // set initial velocity to zero (CMS)
// Set the direction of the spin vector (random)
TRandom3 *rndm = Settings::randomGenerator();
double phi = rndm->Uniform(2*M_PI);
double theta = acos(rndm->Uniform(-1, 1));
CAngle spin(theta, phi);
mGeminiNucleus->setSpinAxis(spin);
//
// Let the nucleus breakup
//
mGeminiNucleus->decay();
mProducts.clear();
if (mGeminiNucleus->abortEvent) {
cout << "Nucleus::breakup(): Error, decay aborted in Gemini++." << endl;
mGeminiNucleus->reset();
return 0;
}
int nfragments = mGeminiNucleus->getNumberOfProducts();
mProducts.resize(nfragments);
for (int i=0; i<nfragments; i++) {
CNucleus *product = mGeminiNucleus->getProducts(i); //set pointer to first stable product
mProducts[i].Z = product->iZ;
mProducts[i].A = product->iA;
mProducts[i].emissionTime = product->getTime();
mProducts[i].p = product->getLorentzVector();
mProducts[i].p.Boost(dissSystem.BoostVector());
mProducts[i].name = product->getName();
if (product->iZ == 1 && product->iA == 1)
mProducts[i].pdgId = 2212; // p
else if (product->iZ == 0 && product->iA == 1)
mProducts[i].pdgId = 2112; // n
else
mProducts[i].pdgId = pdgID(product->iZ, product->iA);
}
mGeminiNucleus->reset();
return nfragments;
}
const vector<BreakupProduct>& FrangibleNucleus::breakupProducts() const
{
return mProducts;
}
void FrangibleNucleus::listBreakupProducts(ostream& os) const
{
TLorentzVector sys;
cout << "Excitation energy = " << mExcitationEnergy << endl;
cout << "Number of fragments = " << mProducts.size() << endl;
for (unsigned int i=0; i<mProducts.size(); i++) {
os << mProducts[i] << endl;
sys += mProducts[i].p;
}
cout << "Invariant mass of dissociated system = " << sys.M() << endl;
cout << "Dissociated system 4-vector p=(" << sys.Px() << ", "
<< sys.Py() << ", " << sys.Pz() << ", " << sys.E() << ')' << endl;
cout << endl;
}
Index: trunk/src/ExclusiveFinalStateGenerator.cpp
===================================================================
--- trunk/src/ExclusiveFinalStateGenerator.cpp (revision 490)
+++ trunk/src/ExclusiveFinalStateGenerator.cpp (revision 491)
@@ -1,612 +1,617 @@
//==============================================================================
// ExclusiveFinalStateGenerator.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 "ExclusiveFinalStateGenerator.h"
-#include "EventGeneratorSettings.h"
-#include "Kinematics.h"
-#include "Event.h"
-#include "VectorMesonDecayMass.h"
-#include "Math/BrentRootFinder.h"
-#include "Math/GSLRootFinder.h"
-#include "Math/RootFinderAlgorithms.h"
-#include "Math/IFunction.h"
-#include <cmath>
-#include <algorithm>
-#include <limits>
-#include <iomanip>
-
+ #include <algorithm>
+ #include <limits>
+ #include <iomanip>
+ #include <cmath>
+ #include "ExclusiveFinalStateGenerator.h"
+ #include "EventGeneratorSettings.h"
+ #include "Kinematics.h"
+ #include "Event.h"
+ #include "VectorMesonDecayMass.h"
+ #include "Math/BrentRootFinder.h"
+ #include "Math/GSLRootFinder.h"
+ #include "Math/RootFinderAlgorithms.h"
+ #include "Math/IFunction.h"
+ #include "Math/SpecFuncMathCore.h"
+ #include "Math/WrappedTF1.h"
+ #include "Math/GaussIntegrator.h"
+ #include "TF1.h"
+
#define PR(x) cout << #x << " = " << (x) << endl;
//-------------------------------------------------------------------------------
//
// Helper class needed to find root in
// ExclusiveFinalStateGenerator::generate()
//
//-------------------------------------------------------------------------------
class ScatteredProtonEnergyFormula : public ROOT::Math::IBaseFunctionOneDim
{
public:
double DoEval(double) const;
ROOT::Math::IBaseFunctionOneDim* Clone() const;
void calculateValidRange(double&, double&);
public:
double mT;
double mVmMass2;
double mMY2;
double mPhi; // azimuthal angle for scattered proton
TLorentzVector mProtonIn;
TLorentzVector mVirtualPhoton;
};
ROOT::Math::IBaseFunctionOneDim* ScatteredProtonEnergyFormula::Clone() const
{
return new ScatteredProtonEnergyFormula();
}
double ScatteredProtonEnergyFormula::DoEval(double Ep) const
{
double m2 = mProtonIn.M2();
double pzp = (mT - m2 - mMY2 + 2*mProtonIn.E() * Ep)/(2*mProtonIn.Pz());
double term = Ep*Ep-pzp*pzp-mMY2;
if (term < 0) return 99999; // out of kinematically allowed range
double ptp = sqrt(term);
TLorentzVector p_out(ptp*cos(mPhi), ptp*sin(mPhi), pzp, Ep);
double f = (mVirtualPhoton + mProtonIn - p_out)*(mVirtualPhoton + mProtonIn - p_out) - mVmMass2;
return f;
}
void ScatteredProtonEnergyFormula::calculateValidRange(double& lower, double& upper)
{
double m2 = mProtonIn.M2();
double term1 = mT-m2-mMY2;
double termA = mProtonIn.E()*term1;
double termB = sqrt(mProtonIn.Pz()*mProtonIn.Pz()*(term1*term1-4*m2*mMY2));
double termC = -2*m2;
lower = (termA+termB)/termC;
upper = (termA-termB)/termC;
if (lower > upper) swap(lower, upper);
lower += numeric_limits<float>::epsilon();
upper -= numeric_limits<float>::epsilon();
}
//-------------------------------------------------------------------------------
//
// Implementation of ExclusiveFinalStateGenerator
//
//-------------------------------------------------------------------------------
ExclusiveFinalStateGenerator::ExclusiveFinalStateGenerator() {/* no op */}
ExclusiveFinalStateGenerator::~ExclusiveFinalStateGenerator() {/* no op */}
bool ExclusiveFinalStateGenerator::generate(int id, double t, double y, double Q2,
bool isIncoherent, int A, Event *event)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
//
// The beam particles must be present in the event list
//
int ePos = -1;
int hPos = -1;
bool parentsOK = true;
if (event->particles.size() == 2) {
if (abs(event->particles[0].pdgId) == 11) {
ePos = 0;
hPos = 1;
}
else if (abs(event->particles[1].pdgId) == 11) {
ePos = 1;
hPos = 0;
}
else
parentsOK = false;
}
else
parentsOK = false;
if (!parentsOK) {
cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
return false;
}
//
// Store arguments locally
// (Some could also be obtained from the event structure)
//
mA = A;
mT = t;
if (mT > 0) mT = -mT; // ensure t<0
mQ2 = Q2;
mY = y;
mIsIncoherent = isIncoherent;
mElectronBeam = event->particles[ePos].p;
mHadronBeam = event->particles[hPos].p;
mMassVM = VectorMesonDecayMass::mass(id);
mS = Kinematics::s(mElectronBeam, mHadronBeam);
//
// Constants
//
double const twopi = 2*M_PI;
double const hMass2 = mHadronBeam.M2();
//
// Incoherent diffarction
//
// Generate hadron dissociation mass according to
// dN/dM2 ~ 1/M2. Lower bound is of course the hadron
// mass and upper bound is some arbitrary value (for now).
//
// Note that we calculate and quote eA kinematics always in
// units of 'per nucleon'. Our model of incoherence is that the
// difference of the diffractive mass of one (1) proton out
// of the nucleus gives the final excitation energy E*.
// Hence we have to calculate E* and divide it by A to keep
// the kinematic consistent.
//
if (mIsIncoherent && mA > 1) {
const double lower = hMass2;
const double upper = 9; // GeV2
mMY2 = lower*upper/(upper - rndm->Uniform()*(upper-lower));
double MY_per_nucleon = (sqrt(hMass2)*(mA-1) + sqrt(mMY2))/mA;
mMY2 = MY_per_nucleon*MY_per_nucleon;
if (mMY2 < hMass2) mMY2 = hMass2;
}
else {
mMY2 = hMass2;
}
//
// Re-engineer scattered electron
//
// e'=(E', pt', pz') -> 3 unknowns
//
// Three equations:
// 1: me*me=E'*E'-pt'*pt'-pz'*pz'
// 2: Q2=-(e-e')^2=-2*me*me + 2*(E*E'-pz*pz')
// 3: W2=(P+e-e')^2=mp2+2*me2+2*(Ep*E-Pz*pz)-2*(Ep*E'-Pz*pz')-2*(E*E'-pz*pz')
//
double Ee = mElectronBeam.E();
double Pe = mElectronBeam.Pz();
double Ep = mHadronBeam.E();
double Pp = mHadronBeam.Pz();
double W = event->W;
double W2 = W*W;
// Take masses from the beams in case they are not actually electrons or protons
double me2 = mElectronBeam.M2();
double mp2 = mHadronBeam.M2();
//
// What we want for each particle:
//
double E, pz, pt, px, py, phi;
//
// Equations 2 and 3 yield:
//
E = Pe*(W2-mp2-2*Ee*Ep) + (Pp+Pe)*Q2 + 2*Pe*Pe*Pp + 2*me2*Pp;
E /= 2*(Ee*Pp-Ep*Pe);
pz = Ee*(W2-mp2) + (Ep+Ee)*Q2 + 2*Ee*Pe*Pp + 2*Ep*me2 - 2*Ee*Ee*Ep;
pz /= 2*(Ee*Pp-Ep*Pe);
//
// Equation 1:
//
pt = sqrt(E*E-pz*pz-me2);
phi = rndm->Uniform(twopi);
TLorentzVector theScatteredElectron(pt*sin(phi), pt*cos(phi), pz, E);
//
// Re-engineer virtual photon
//
// gamma=E-E'
E = mElectronBeam.E()-theScatteredElectron.E();
pz = mElectronBeam.Pz()-theScatteredElectron.Pz();
px = mElectronBeam.Px()-theScatteredElectron.Px();
py = mElectronBeam.Py()-theScatteredElectron.Py();
TLorentzVector theVirtualPhoton = TLorentzVector(px, py, pz, E);
//
// Re-engineer scattered proton/dissociated proton
//
// No analytic solution. Need to run a root finder that does
// not need derivates but uses a bracketing algorithm (Brent).
// Correct brackets are crucial since ScatteredProtonEnergyFormula
// produces sqrt(-x) if outside the kinematically allowed range (it
// actually catches it and returns a large positive number, 0 doesn't
// work).
//
//
// Setup formula to solve root
//
phi = rndm->Uniform(twopi);
ScatteredProtonEnergyFormula formula;
formula.mT = mT;
formula.mVmMass2 = mMassVM*mMassVM;
formula.mPhi = phi;
formula.mProtonIn = mHadronBeam;
formula.mVirtualPhoton = theVirtualPhoton;
formula.mMY2 = mMY2;
//
// Find correct brackets to start with
//
double lower, upper;
formula.calculateValidRange(lower, upper);
if (upper > mHadronBeam.E() + theVirtualPhoton.E()) // limit excessive values
upper = mHadronBeam.E() + theVirtualPhoton.E(); // make it easier for Brent
//
// Run root finder
//
ROOT::Math::BrentRootFinder rootfinder;
rootfinder.SetFunction(formula, lower, upper);
rootfinder.Solve(10000, 0, 1.e-12);
E = rootfinder.Root();
if (/* rootfinder.Status() || */ fabs(formula(E)) > 1e-6) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, cannot find root. No final state defined." << endl;
- return false;
- }
+ return false;
+ }
//
// Outgoing proton (hadron) system
//
pz = (mT- hMass2 - mMY2 + 2*mHadronBeam.E()*E)/(2*mHadronBeam.Pz());
pt = sqrt(E*E-pz*pz-mMY2);
px = pt*cos(phi);
py = pt*sin(phi);
TLorentzVector theScatteredProton(px, py, pz, E);
//
// Finally the vector meson
//
TLorentzVector theVectorMeson((mHadronBeam + mElectronBeam) - (theScatteredElectron + theScatteredProton));
//
// Check for numerical glitches
//
if (!isValid(theScatteredElectron)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered electron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theScatteredProton)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered hadron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theVectorMeson)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, vector meson 4-vector is invalid." << endl;
return false;
}
//
// Add particles to event record
//
event->particles.resize(2+5);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int vm = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[vm].index = vm;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[vm].p = theVectorMeson;
event->particles[pomeron].p = theScatteredProton - mHadronBeam;
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[vm].pdgId = id;
event->particles[pomeron].pdgId = 990;
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[vm].status = 1;
event->particles[pomeron].status = 2;
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[vm].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(vm);
event->particles[gamma].daughters.push_back(pomeron);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
return true;
}
//
// UPC version
//
// Note: We call the beam particle that emits the photon "electron"
// even if its a proton/nucleus
//
bool ExclusiveFinalStateGenerator::generate(int id, double t, double xpom,
bool isIncoherent, int A, Event *event)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
//
// The beam particles must be present in the event list
//
int ePos = 0;
int hPos = 1;
bool parentsOK = true;
if (event->particles.size() != 2)
parentsOK = false;
if (!parentsOK) {
cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
return false;
}
//
// Store arguments locally
// (Some could also be obtained from the event structure)
//
mA = A;
mT = t;
if (mT > 0) mT = -mT; // ensure t<0
mIsIncoherent = isIncoherent;
mElectronBeam = event->particles[ePos].p;
mHadronBeam = event->particles[hPos].p;
+
//
// For the final state we pick a VM mass according to the
// specific line shape/width.
//
mMassVM = VectorMesonDecayMass::mass(id);
// mMassVM = settings->lookupPDG(id)->Mass(); // old fixed mass version
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mXp = xpom;
//
// Constants
//
double const twopi = 2*M_PI;
double const hMass2 = mHadronBeam.M2();
//
// Incoherent diffraction
//
// Generate hadron dissociation mass according to
// dN/dM2 ~ 1/M2. Lower bound is of course the hadron
// mass and upper bound is some arbitrary value (for now).
//
// Note that we calculate and quote eA kinematics always in
// units of 'per nucleon'. Our model of incoherence is that the
// difference of the diffractive mass of one (1) proton out
// of the nucleus gives the final excitation energy E*.
// Hence we have to calculate E* and divide it by A to keep
// the kinematic consistent.
//
if (mIsIncoherent && mA > 1) {
const double lower = hMass2;
const double upper = 9; // GeV2
mMY2 = lower*upper/(upper - rndm->Uniform()*(upper-lower));
double MY_per_nucleon = (sqrt(hMass2)*(mA-1) + sqrt(mMY2))/mA;
mMY2 = MY_per_nucleon*MY_per_nucleon;
if (mMY2 < hMass2) mMY2 = hMass2;
}
else {
mMY2 = hMass2;
}
//
// Check for smallest allowed xpom
//
double xpom_min = Kinematics::xpomMin(mMassVM, mT, mHadronBeam, mElectronBeam, mMY2-hMass2);
if (xpom < xpom_min){
if (settings->verboseLevel() > 2) cout<<"xpom = "<<xpom<<", which is smaller than xpom_min="<<xpom_min<<endl;
return false;
}
//
// We have three unknown: E_gamma, P_z,gamma, E_pomeron
// and three equations M^2(scattered electrons), M^2(scattered proton),
// and M^2(vector meson)
//
//Start by setting up the knowns: Incoming beams, xpom, and t:
double Ee = mElectronBeam.E();
double Pe = mElectronBeam.Pz();
double Pp = mHadronBeam.Pz();
double Ep = mHadronBeam.E();
//
// What we want for each particle:
//
double E, pz, pt, phi;
//
// Use scattered proton to set up four momentum of pomeron:
//
double sqrtarg= 1 - ( (2*xpom - xpom*xpom)*Pp*Pp + mT + hMass2 - mMY2) / (Ep*Ep);
E = Ep*(1-sqrt(sqrtarg)); //this is close to xpom*Ep
phi = rndm->Uniform(twopi);
pt = sqrt(-t);
pz = xpom*Pp;
TLorentzVector thePomeron = TLorentzVector(pt*sin(phi), pt*cos(phi), xpom*Pp, E);
//
// That give scattered proton:
//
TLorentzVector theScatteredProton = mHadronBeam-thePomeron;
//
// Fill Virtual Photon:
//
mEgam = Kinematics::Egamma(xpom, t, mMassVM, Ep, Ee, mMY2-hMass2);
E = mEgam;
pz = Pe*(1 - sqrt( 1 - ( 2*Ee*E - E*E )/(Pe*Pe) ) ); //this is close to Egamma
TLorentzVector theVirtualPhoton = TLorentzVector(0, 0, pz, E);
//
// Scattered Electron:
//
TLorentzVector theScatteredElectron = mElectronBeam - theVirtualPhoton ;
//
// Vector Meson:
//
TLorentzVector theVectorMeson = theVirtualPhoton + thePomeron;
//
// Check for numerical glitches
//
if (!isValid(theScatteredElectron)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered electron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theScatteredProton)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered hadron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theVectorMeson)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, vector meson 4-vector is invalid." << endl;
return false;
}
//
// Add particles to event record
//
event->particles.resize(2+5);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int vm = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[vm].index = vm;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[vm].p = theVectorMeson;
event->particles[pomeron].p = thePomeron;
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[vm].pdgId = id;
event->particles[pomeron].pdgId = 990;
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[vm].status = 1;
event->particles[pomeron].status = 2;
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[vm].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(vm);
event->particles[gamma].daughters.push_back(pomeron);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
//fill event structure
double y = mHadronBeam*theVirtualPhoton/(mHadronBeam*mElectronBeam);
double W2 = (theVirtualPhoton+mHadronBeam).M2();
mQ2 = -theVirtualPhoton.M2();
event->Q2 = mQ2;
event->W = sqrt(W2);
event->y = y;
return true;
}
Index: trunk/src/Version.h
===================================================================
--- trunk/src/Version.h (revision 490)
+++ trunk/src/Version.h (revision 491)
@@ -1,29 +1,29 @@
//==============================================================================
// Version.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 Version_h
#define Version_h
-#define VERSION "1.40"
-#define VERSION_NUMBER 1.40
+#define VERSION "2.00"
+#define VERSION_NUMBER 2.00
#endif
Index: trunk/src/DipoleModel.cpp
===================================================================
--- trunk/src/DipoleModel.cpp (revision 490)
+++ trunk/src/DipoleModel.cpp (revision 491)
@@ -1,337 +1,369 @@
//==============================================================================
// DipoleModel.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include <fstream>
#include <iostream>
#include <sstream>
#include <cmath>
#include "DipoleModel.h"
#include "TableGeneratorSettings.h"
#include "DglapEvolution.h"
#include "Constants.h"
#include "TFile.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH2F.h"
#include "TF1.h"
+#include "Math/WrappedTF1.h"
+#include "Math/GaussIntegrator.h"
+
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
DipoleModel::DipoleModel()
{
mConfigurationExists = false;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = settings->A();
mNucleus.init(A);
mIsInitialized = true;
mParameters = nullptr;
-}
+}
DipoleModel::~DipoleModel()
{
delete mParameters;
}
const TableGeneratorNucleus* DipoleModel::nucleus() const { return &mNucleus; }
bool DipoleModel::configurationExists() const { return mConfigurationExists; }
double DipoleModel::bDependence(double) { return 0; }
double DipoleModel::bDependence(double, double) { return 0; }
double DipoleModel::dsigmadb2ep(double, double, double) { return 0;}
//***********bSat:*****************************************************
DipoleModel_bSat::DipoleModel_bSat()
{
mBDependence = 0;
mSigma_ep_LookupTable = 0;
//
// Set the parameters. Note that we enforce here the bSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bSat, settings->dipoleModelParameterSet());
}
+DipoleModel_bSat::DipoleModel_bSat(Settings* settings)
+{
+ //
+ // Set the parameters. Note that we enforce here the bSat model
+ // independent of what the settings say.
+ //
+ mParameters = new DipoleModelParameters(settings);
+}
+
+
DipoleModel_bSat::~DipoleModel_bSat()
{
delete mSigma_ep_LookupTable;
delete mBDependence;
-}
+ delete mParameters;
+}
DipoleModel_bSat& DipoleModel_bSat::operator=(const DipoleModel_bSat& dp)
{
if (this != &dp) {
delete mBDependence;
delete mSigma_ep_LookupTable;
DipoleModel::operator=(dp);
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable));
mSigma_ep_LookupTable->SetDirectory(0);
}
return *this;
}
DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp)
{
+ if (mBDependence) delete mBDependence;
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
}
void DipoleModel_bSat::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = mNucleus.A();
- string path = settings->bSatLookupPath();
+ string path=settings->bSatLookupPath();
ostringstream filename;
filename.str("");
- filename << path;// << "/bSat_bDependence_A" << A <<".root";
+ filename << path << "/bSat_bDependence_A" << A <<".root";
ifstream ifs(filename.str().c_str());
if (!ifs) {
cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl;
cout << "Stopping." << endl;
exit(1);
}
TFile* lufile= new TFile(filename.str().c_str());
ostringstream histoName;
histoName.str( "" );
histoName << "Configuration_" << iConfiguration;
lufile->GetObject( histoName.str().c_str(), mBDependence );
mBDependence->SetDirectory(0);
lufile->Close();
- mConfigurationExists = true;
+ mConfigurationExists=true;
}
double DipoleModel_bSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- double bDep = bDependence(b, phi);
+ double bDep=bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
double DipoleModel_bSat::bDependence(double b, double phi)
{
double result = mBDependence->Interpolate(b, phi);
return result;
}
double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
+ if(muQ2>1e7 or muQ2<1) PR(muQ2);
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
-double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe) {
+double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe) {
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- xprobe*=1.;
double sigmap = mSigma_ep_LookupTable->Interpolate(r);
- int A = nucleus()->A();
- double TA = nucleus()->T(b)/A;
+ int A=nucleus()->A();
+ double TA=nucleus()->T(b)/A;
double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) );
return result;
}
void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe)
{
- double rbRange = 3.*nucleus()->radius();
+ double rbRange=3.*nucleus()->radius();
TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2);
+ if(mSigma_ep_LookupTable) delete mSigma_ep_LookupTable;
mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange);
- dsigmaForIntegration->SetNpx(1000);
+ ROOT::Math::WrappedTF1* WFlookup=new ROOT::Math::WrappedTF1(*dsigmaForIntegration);
+ ROOT::Math::GaussIntegrator GIlookup;
+ GIlookup.SetFunction(*WFlookup);
+ GIlookup.SetRelTolerance(1e-8);
+ GIlookup.SetAbsTolerance(0.);
+
+// dsigmaForIntegration->SetNpx(1000);
for (int iR=1; iR<=1000; iR++) {
- double r = mSigma_ep_LookupTable->GetBinCenter(iR);
+ double r=mSigma_ep_LookupTable->GetBinCenter(iR);
dsigmaForIntegration->SetParameter(0, r);
dsigmaForIntegration->SetParameter(1, xprobe);
- double result = dsigmaForIntegration->Integral(0, rbRange);
+// double result=dsigmaForIntegration->Integral(0, rbRange);
+ double result = GIlookup.IntegralUp(0.); //GeV-2
mSigma_ep_LookupTable->SetBinContent(iR, result);
}
- delete dsigmaForIntegration;
-}
+ delete dsigmaForIntegration; //GeV-2
+ delete WFlookup;
+}
double DipoleModel_bSat::dsigmadb2epForIntegration(double *x, double* par)
{
- double r = par[0];
- double xprobe = par[1];
- double b = *x;
- return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe);
-}
+ double r=par[0]; //fm
+ double xprobe=par[1];
+ double b = *x; //fm
+ return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe); //GeV-2fm-1
+}
//***********bNonSat:*************************************************
DipoleModel_bNonSat::DipoleModel_bNonSat()
{
//
// Set the parameters. Note that we need bNonSat to calculate
// the skewedness correction for bSat.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
+ mParameters = new DipoleModelParameters( settings->dipoleModelType(), settings->dipoleModelParameterSet());
+}
+DipoleModel_bNonSat::DipoleModel_bNonSat(Settings* settings)
+{
+ //
+ // Set the parameters. Note that we need bNonSat to calculate
+ // the skewedness correction for bSat.
+ //
mParameters = new DipoleModelParameters(settings->dipoleModelType(), settings->dipoleModelParameterSet());
}
-DipoleModel_bNonSat::~DipoleModel_bNonSat(){}
+DipoleModel_bNonSat::~DipoleModel_bNonSat(){
+ delete mParameters;
+}
double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- double bDep = bDependence(b, phi);
+ double bDep=bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- int A = nucleus()->A();
- double TA = nucleus()->T(b)/A;
+ int A=nucleus()->A();
+ double TA=nucleus()->T(b)/A;
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
- double result = A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
+ double result=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
return result;
}
//***********bCGC:*****************************************************
DipoleModel_bCGC::DipoleModel_bCGC()
{
//
// Set the parameters. Note that we enforce here the bNonSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bCGC, settings->dipoleModelParameterSet());
}
void DipoleModel_bCGC::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
- //iConfiguration is a dummy for bCGC. Do this to avoid compiler warnings:
- iConfiguration++;
mNucleus.generate();
- mConfigurationExists = true;
+ mConfigurationExists=true;
}
double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x)
{
- double result = 1;
+ double result=1;
for (unsigned int iA=0; iA<mNucleus.A(); iA++) {
double absdeltab=( TVector3(b*cos(phi), b*sin(phi), 0.)
-mNucleus.configuration.at(iA).position() ).Perp();
result*=(1.-0.5*dsigmadb2ep(r, absdeltab, x));
}
return 2.*(1.-result);
}
double DipoleModel_bCGC::dsigmadb2ep(double r, double b, double xprobe)
{
double Y = log(1/xprobe);
double kappa = mParameters->kappa();
double N0 = mParameters->N0();
double x0 = mParameters->x0();
double lambda = mParameters->lambda();
double gammas = mParameters->gammas();
double A = -N0*N0*gammas*gammas/((1-N0)*(1-N0)*log(1-N0));
double B = 0.5*pow(1-N0,-(1-N0)/(N0*gammas));
double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b));
double rQs = r*Qs/hbarc;
- double result = 0;
+ double result=0;
if (rQs <= 2)
result = 2*N0*pow(0.5*rQs, 2*(gammas+(1/(kappa*lambda*Y))*log(2/rQs)));
else
result = 2*(1 - exp(-A*log(B*rQs)*log(B*rQs)));
return result;
}
double DipoleModel_bCGC::bDependence(double b)
{
double gammas = mParameters->gammas();
double Bcgc = mParameters->Bcgc();
return exp(-0.5*b*b/Bcgc/gammas/hbarc2);
}
Index: trunk/src/Settings.h
===================================================================
--- trunk/src/Settings.h (revision 490)
+++ trunk/src/Settings.h (revision 491)
@@ -1,190 +1,197 @@
//==============================================================================
// 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;
+
int vectorMesonId() const;
void setVectorMesonId(int);
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 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;
bool mVerbose;
int mVerboseLevel;
double mQ2min;
double mQ2max;
double mWmin;
double mWmax;
double mXpomMin;
double mXpomMax;
+ double mbetamin;
+ double mbetamax;
int mVectorMesonId;
string mDipoleModelName;
DipoleModelType mDipoleModelType;
DipoleModelParameterSet mDipoleModelParameterSet;
string mDipoleModelParameterSetName;
TableSetType mTableSetType;
string mTableSetTypeName;
unsigned int mA;
string mRootfile;
bool mUPC;
unsigned int mUPCA;
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;}
#endif
Index: trunk/src/ModeFinderFunctor.cpp
===================================================================
--- trunk/src/ModeFinderFunctor.cpp (revision 490)
+++ trunk/src/ModeFinderFunctor.cpp (revision 491)
@@ -1,110 +1,183 @@
//==============================================================================
// ModeFinderFunctor.cpp
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// 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: Thomas Ullrich
-// Last update:
+// Author: Tobias Toll
+// Last update:
// $Date$
// $Author$
//==============================================================================
#include "ModeFinderFunctor.h"
-#include "CrossSection.h"
-#include "Kinematics.h"
+#include "CrossSection.h"
+#include "Kinematics.h"
#include <iostream>
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
ModeFinderFunctor::ModeFinderFunctor()
{
mCrossSection = 0;
mQ2 = mVmMass = mMinT = mMaxT = 0;
}
ModeFinderFunctor::ModeFinderFunctor(CrossSection* cs, double Q2, double vmMass, double tmin, double tmax)
{
mCrossSection = cs;
mQ2 = Q2;
mVmMass = vmMass;
mMinT = tmin;
mMaxT = tmax;
}
-
+
+InclusiveDiffractionModeFinderFunctor::InclusiveDiffractionModeFinderFunctor(InclusiveDiffractiveCrossSections* cs, double Q2, double W2, double z, double MX2min, double MX2max)
+{
+ mIDCrossSection = cs;
+ mQ2 = Q2;
+ mW2 = W2;
+ mZ = z;
+
+ mMX2max=MX2max;
+ mMX2min=MX2min;
+}
+
+InclusiveDiffractionModeFinderFunctor::InclusiveDiffractionModeFinderFunctor(InclusiveDiffractiveCrossSections* cs, double Q2, double W2, double z, double s, double MX2min, double MX2max, double Q2min, double Q2max, double W2min, double W2max, double ymin, double betamin, double betamax, double mu02){
+ mIDCrossSection = cs;
+ mQ2 = Q2;
+ mW2 = W2;
+ mZ = z;
+ mS = s;
+ mMinMX2 = MX2min;
+ mMaxMX2 = MX2max;
+
+ mQ2max=Q2max;
+ mQ2min=Q2min;
+ mW2max=W2max;
+ mW2min=W2min;
+ mYmin=ymin;
+ mBetamin=betamin;
+ mBetamax=betamax;
+ mMu02=mu02;
+ mMX2max=MX2max;
+ mMX2min=MX2min;
+}
+
+ROOT::Math::IGenFunction* InclusiveDiffractionModeFinderFunctor::Clone() const
+{
+ return new InclusiveDiffractionModeFinderFunctor(mIDCrossSection, mQ2, mW2, mZ, mS, mMinMX2, mMaxMX2, mQ2min, mQ2max, mW2min, mW2max, mYmin, mBetamin, mBetamax, mMu02);
+}
+
+double InclusiveDiffractionModeFinderFunctor::DoEval(double MX2) const {
+// if(mQ2<mQ2min or mQ2>mQ2max){
+// cout<<"Q2 out of range! Q2="<<mQ2<<", ["<<mQ2min<<", "<<mQ2max<<"]"<<endl;
+// return 0;
+// }
+ if(MX2>mMX2max or MX2<mMX2min){
+// if(cout<<"MX2 out of range! MX2="<<MX2<<", ["<<mMX2min<<", "<<mMX2max<<"]"<<endl;
+ return 0;
+ }
+// if(mW2>mW2max or mW2<mW2min){
+// cout<<"MX2 out of range! MX2="<<MX2<<", ["<<mMX2min<<", "<<mMX2max<<"]"<<endl;
+// return 0;
+// }
+// double x = Kinematics::x(mQ2, mW2);
+// double y = Kinematics::y(mQ2, x, mS);
+// double beta=mQ2/(mQ2+MX2);
+// double xpom=x/beta;
+// if(x<0 or x>1 or beta<0 or beta>1 or xpom<0 or xpom>1){
+// PR(x);
+// PR(y);
+// PR(beta);
+// PR(xpom);
+// return 0;
+// }
+ if (mIDCrossSection) {
+ double result = (*mIDCrossSection)(MX2, mQ2, mW2, mZ);
+ return -result;
+ }
+ else {
+ return 0;
+ }
+}
+
+void InclusiveDiffractionModeFinderFunctor::setQuarkMass(double val) {mMq=val;}
+
void ModeFinderFunctor::setVmMass(double val) {mVmMass = val;}
void ModeFinderFunctor::setQ2(double val) {mQ2 = val;}
void ModeFinderFunctor::setMinT(double val) {mMinT = val;}
void ModeFinderFunctor::setMaxT(double val) {mMaxT = val;}
double ModeFinderFunctor::DoEval(double W2) const
{
if (mCrossSection) {
double t = Kinematics::tmax(0, mQ2, W2, mVmMass); // first arg (t) set to 0 here
if (t > mMaxT) t = mMaxT; // don't exceed given (table) maximum (smallest |t|)
if (t < mMinT) return 0; // lower table limits (treat different than max limit)
double result = (*mCrossSection)(t, mQ2, W2); // t, Q2, W2
return -result; // minimum=maximum
}
else {
return 0;
}
}
ROOT::Math::IGenFunction* ModeFinderFunctor::Clone() const
{
return new ModeFinderFunctor(mCrossSection, mQ2, mVmMass, mMinT, mMaxT);
}
UPCModeFinderFunctor::UPCModeFinderFunctor()
{
mCrossSection = 0;
mVmMass = mHBeamEnergy = mEBeamEnergy = 0;
}
UPCModeFinderFunctor::UPCModeFinderFunctor(CrossSection* cs, double vmMass, double hEnergy, double eEnergy)
{
mCrossSection = cs;
mVmMass = vmMass;
mHBeamEnergy = hEnergy;
mEBeamEnergy = eEnergy;
}
double UPCModeFinderFunctor::DoEval(double val) const
{
if (mCrossSection) {
double xpom = exp(val); // x comes as log(x)
double t = Kinematics::tmax(xpom);
bool ok = Kinematics::validUPC(mHBeamEnergy, mEBeamEnergy, t, xpom, mVmMass, false);
if (ok) {
double result = (*mCrossSection)(t, xpom);
return -result; // minimum=maximum
}
else {
return 0;
}
}
else {
return 0;
}
}
ROOT::Math::IGenFunction* UPCModeFinderFunctor::Clone() const
{
return new UPCModeFinderFunctor(mCrossSection, mVmMass, mHBeamEnergy, mEBeamEnergy);
}
+
Index: trunk/src/TableGeneratorSettings.cpp
===================================================================
--- trunk/src/TableGeneratorSettings.cpp (revision 490)
+++ trunk/src/TableGeneratorSettings.cpp (revision 491)
@@ -1,217 +1,217 @@
//==============================================================================
// TableGeneratorSettings.cpp
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// 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(&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(&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;
}
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; }
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; }
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/Kinematics.cpp
===================================================================
--- trunk/src/Kinematics.cpp (revision 490)
+++ trunk/src/Kinematics.cpp (revision 491)
@@ -1,455 +1,531 @@
//==============================================================================
// Kinematics.cpp
//
-// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
-// This file is part of Sartre.
+// 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.
+// 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:
+// Last update:
// $Date$
// $Author$
//------------------------------------------------------------------------------
//
// Note that Kinematics is only meant for ep and eA scattering.
// It is not applicable for UPC kinematics. Use UpcKinematics instead.
//
//==============================================================================
-#include "Kinematics.h"
-#include <cmath>
-#include <iostream>
+#include "Kinematics.h"
+#include <cmath>
+#include <iostream>
#include "Math/RootFinder.h"
#include "Math/BrentRootFinder.h"
#include "Math/GSLRootFinder.h"
-#include "Math/RootFinderAlgorithms.h"
+#include "Math/RootFinderAlgorithms.h"
#include "Math/IFunction.h"
#include "Math/BrentMinimizer1D.h"
// TMP
#include "TH1D.h"
#include "TFile.h"
-#define PR(x) cout << #x << " = " << (x) << endl;
-
-using namespace std;
-
-bool Kinematics::mError = false;
-
-bool Kinematics::error() {return mError;}
+#define PR(x) cout << #x << " = " << (x) << endl;
+
+using namespace std;
+
+bool Kinematics::mError = false;
+
+bool Kinematics::error() {return mError;}
//UPC:
double Kinematics::mXpomMin=0;
double Kinematics::mTold=0;
bool Kinematics::mXpomMinIsEvaluated=false;
//
-double Kinematics::xprobe(double Q2, double W2, double vmMass)
-{
- mError = false;
- double val = ( Q2 + vmMass*vmMass ) / (W2 - protonMass2 + Q2);
- if (val > 1 || val < 0) mError = true;
- return val;
-}
+double Kinematics::xprobe(double Q2, double W2, double vmMass)
+{
+ mError = false;
+ double val = ( Q2 + vmMass*vmMass ) / (W2 - protonMass2 + Q2);
+ if (val > 1 || val < 0) mError = true;
+ return val;
+}
TLorentzVector Kinematics::electronBeam(double eE, bool upc)
-{
+{
mError = false;
if (upc)
return TLorentzVector(0, 0, -sqrt(eE*eE-protonMass2), eE);
else
return TLorentzVector(0, 0, -sqrt(eE*eE-electronMass2), eE);
-}
-
-TLorentzVector Kinematics::hadronBeam(double hE)
-{
- mError = false;
- return TLorentzVector(0, 0, sqrt(hE*hE-protonMass2), hE);
-}
-
-double Kinematics::s(const TLorentzVector& e, const TLorentzVector& p)
-{
- mError = false;
- return (e+p)*(e+p);
-}
-
-double Kinematics::s(double eE, double hE, bool upc)
-{
- mError = false;
- TLorentzVector e = electronBeam(eE, upc);
- TLorentzVector p = hadronBeam(hE);
- return (e+p)*(e+p);
-}
-
-double Kinematics::x(double Q2, double W2)
-{
- mError = false;
- double val = Q2/(W2-protonMass2+Q2);
- if (val > 1 || val < 0) mError = true;
+}
+
+TLorentzVector Kinematics::hadronBeam(double hE)
+{
+ mError = false;
+ return TLorentzVector(0, 0, sqrt(hE*hE-protonMass2), hE);
+}
+
+double Kinematics::s(const TLorentzVector& e, const TLorentzVector& p)
+{
+ mError = false;
+ return (e+p)*(e+p);
+}
+
+double Kinematics::s(double eE, double hE, bool upc)
+{
+ mError = false;
+ TLorentzVector e = electronBeam(eE, upc);
+ TLorentzVector p = hadronBeam(hE);
+ return (e+p)*(e+p);
+}
+
+double Kinematics::x(double Q2, double W2)
+{
+ mError = false;
+ double val = Q2/(W2-protonMass2+Q2);
+ if (val > 1 || val < 0) mError = true;
return val;
-}
-
-double Kinematics::y(double Q2, double x, double s)
-{
- mError = false;
- double val = Q2/(x*(s-electronMass2-protonMass2));
- if (val > 1 || val < 0) mError = true;
- return val;
-}
-
-double Kinematics::ymin(double s, double vmMass)
-{
- mError = false;
- double W2min = Kinematics::W2min(vmMass);
- double val = (s + W2min - sqrt((s-W2min)*(s-W2min)-4*electronMass2*W2min))/2/(s+electronMass2);
- if (val > 1 || val < 0) mError = true;
- return val;
-}
-
-double Kinematics::ymax(double s, double vmMass)
-{
- mError = false;
- double W2min = Kinematics::W2min(vmMass);
- double val = (s + W2min + sqrt((s-W2min)*(s-W2min)-4*electronMass2*W2min))/2/(s+electronMass2);
- if (val > 1 || val < 0) mError = true;
- return val;
-}
-
-double Kinematics::W2(double Q2, double x)
-{
- //
- // Returns W2 at a given Q2 and Bjorken x
- //
- mError = false;
- return protonMass2 + Q2*(1-x)/x;
-}
-
-double Kinematics::W2(double Q2, double xprobe, double vmMass)
-{
- //
- // Returns W2 at a given Q2 and xprobe
- // where xprobe is
- // xprobe = xBJ * (1 + MV^2/Q2)
- //
- mError = false;
- return protonMass2 + vmMass*vmMass/xprobe + Q2*(1-xprobe)/xprobe;
-}
-
-double Kinematics::W2min(double vmMass)
-{
- mError = false;
- // return vmMass*vmMass + protonMass2;
- return (vmMass + protonMass)*(vmMass + protonMass);
-}
-
-double Kinematics::W2max(double s)
-{
- mError = false;
- return s; // check !?
-}
-
-double Kinematics::Q2min(double y)
-{
- mError = false;
- return electronMass2*y*y/(1-y);
-}
-
-double Kinematics::Q2max(double s)
-{
- mError = false;
- return s-electronMass2-protonMass2;
-}
-
-double Kinematics::xpomeron(double t, double Q2, double W2, double vmM)
-{
- mError = false;
- double val = (vmM*vmM + Q2 - t)/(W2 + Q2 - protonMass2); // full expression
- if (val > 1 || val < 0) mError = true;
- return val;
-}
-
-double Kinematics::tmax(double xpom)
-{
- mError = false;
- return (- xpom*xpom * protonMass2 / (1-xpom)); // only if m_p (in) = m_p' (out) - to be fixed
-}
-
-double Kinematics::tmax(double t, double Q2, double W2, double vmM)
-{
- return tmax(xpomeron(t, Q2, W2, vmM));
-}
-
-double Kinematics::tmin(double eE)
-{
- mError = false;
- return -2*(eE*protonMass-protonMass2); // check ?!
-}
+}
+
+double Kinematics::y(double Q2, double x, double s)
+{
+ mError = false;
+ double val = Q2/(x*(s-electronMass2-protonMass2));
+ if (val > 1 || val < 0) mError = true;
+ if( val>1 ) val=0.9999; //#TT hack
+ return val;
+}
+
+double Kinematics::ymin(double s, double vmMass)
+{
+ mError = false;
+ double W2min = Kinematics::W2min(vmMass);
+ double val = (s + W2min - sqrt((s-W2min)*(s-W2min)-4*electronMass2*W2min))
+ /2/(s+electronMass2);
+ if (val > 1 || val < 0) mError = true;
+ return val;
+}
+
+double Kinematics::ymax(double s, double vmMass)
+{
+ mError = false;
+ double W2min = Kinematics::W2min(vmMass);
+ double val = (s + W2min + sqrt((s-W2min)*(s-W2min)-4*electronMass2*W2min))/2/(s+electronMass2);
+ if (val > 1 || val < 0) mError = true;
+ return val;
+}
+
+double Kinematics::W2(double Q2, double x)
+{
+ //
+ // Returns W2 at a given Q2 and Bjorken x
+ //
+ mError = false;
+ return protonMass2 + Q2*(1-x)/x;
+}
+
+double Kinematics::W2(double Q2, double xprobe, double vmMass)
+{
+ //
+ // Returns W2 at a given Q2 and xprobe
+ // where xprobe is
+ // xprobe = xBJ * (1 + MV^2/Q2)
+ //
+ mError = false;
+ return protonMass2 + vmMass*vmMass/xprobe + Q2*(1-xprobe)/xprobe;
+}
+
+double Kinematics::W2min(double vmMass)
+{
+ mError = false;
+ return vmMass*vmMass + protonMass2;
+}
+
+double Kinematics::W2max(double s)
+{
+ mError = false;
+ return s; // check !?
+}
+double Kinematics::W2max(double s, double Q2, double MX)
+{
+ double ymax=Kinematics::ymax(s, MX);
+ mError = false;
+ return ymax*s-Q2+protonMass2; // check !?
+}
+double Kinematics::betamin(double Q2, double s, double MX, double t){
+ double ymax=Kinematics::ymax(s, MX);
+ return Q2/(ymax*s-protonMass2-t);
+}
+
+double Kinematics::Q2min(double y)
+{
+ mError = false;
+ return electronMass2*y*y/(1-y);
+}
+
+double Kinematics::Q2max(double s)
+{
+ mError = false;
+ return s-electronMass2-protonMass2;
+}
+double Kinematics::xpomeron(double t, double Q2, double W2, double vmM)
+{
+ mError = false;
+ double val = (vmM*vmM + Q2 - t)/(W2 + Q2 - protonMass2); // full expression
+ if (val > 1 || val < 0) mError = true;
+ return val;
+}
+
+double Kinematics::tmax(double xpom)
+{
+ mError = false;
+ return (- xpom*xpom * protonMass2 / (1-xpom)); // only if m_p (in) = m_p' (out) - to be fixed
+}
+
+double Kinematics::tmax(double t, double Q2, double W2, double vmM)
+{
+ return tmax(xpomeron(t, Q2, W2, vmM));
+}
+
+double Kinematics::tmin(double eE)
+{
+ mError = false;
+ return -2*(eE*protonMass-protonMass2); // check ?!
+}
+double Kinematics::MX2(double Q2, double beta, double t){
+ return Q2*(1-beta)/beta+t;
+}
//
// UPC only
//
double Kinematics::xpomMin(double massVM, double t, TLorentzVector hBeam, TLorentzVector eBeam, double MY2minusM2){
- if (mXpomMinIsEvaluated && mTold == t) return mXpomMin;
+ if(mXpomMinIsEvaluated and mTold == t) return mXpomMin;
LowerXpomeronFormula formula;
mTold = t;
formula.mT = t;
formula.mVmMass2=massVM*massVM;
formula.mElectronBeam = eBeam;
formula.mHadronBeam = hBeam;
formula.mMY2minusM2 = MY2minusM2;
//
//Find starting brackets
//
double lower=1e-12;
double upper=1e-1;
//
// Find minimum and use as lower
// bracket to avoid two roots
//
ROOT::Math::BrentMinimizer1D bm;
bm.SetFunction(formula, log(lower), log(upper));
bm.Minimize(100000, 1e-8, 1e-10);
double theMin = bm.XMinimum();
-
+
//
// Run root finder
//
ROOT::Math::RootFinder rootfinder(ROOT::Math::RootFinder::kBRENT);
rootfinder.SetFunction(formula, theMin, log(upper));
rootfinder.Solve(100000000, 1e-14, 0);
-
+
//
// Now some fine adjustements
//
double bestGuess = rootfinder.Root();
int ntimes = 0;
while (formula.DoEval(bestGuess) < 0) {
bestGuess += 1e-14;
ntimes++;
if (ntimes> 10000) break;
}
double result = exp(bestGuess);
double s = (hBeam+eBeam).M2();
mXpomMin = max(result, xpomMin(s, massVM, t));
mXpomMinIsEvaluated = true;
return mXpomMin;
}
double Kinematics::xpomMin(double s, double vmM, double t) {
- //
- // This is the threshold value.
- // In some cases a more restrictive value may be needed
- //
- return (electronMass2 + protonMass2 + vmM*vmM - t)/s;
+ //
+ // This is the threshold value.
+ // In some cases a more restrictive value may be needed
+ //
+ return (electronMass2 + protonMass2 + vmM*vmM - t)/s;
}
double Kinematics::Egamma(double xpom, double t, double vmM, double hBeamEnergy, double eBeamEnergy, double MY2minusM2)
{
//
// We use that the vector meson = photon + pomeron to calculate Egamma
// (keeping the pomeron variables from above)
// We solve a quadratic equation to get the photon energy such that
// the vector meson has the correct mass.
//
double Ep=hBeamEnergy;
double Pp=-sqrt(Ep*Ep-protonMass2);
double Ee=eBeamEnergy;
double Pe=sqrt(Ee*Ee-protonMass2);
double sqrtarg= 1 - ( (2*xpom - xpom*xpom)*Pp*Pp + t - MY2minusM2) / (Ep*Ep);
double Epom=Ep*(1-sqrt(sqrtarg)); //this is close to xpom*Ep
double Ppom=xpom*Pp;
double AA = vmM*vmM - Epom*Epom - t + Ppom*Ppom + 2*Pe*(Pe+Ppom); //GeV^2
double aa = 4*(Pe+Ppom)*(Pe+Ppom) - 4*(Ee+Epom)*(Ee+Epom); //GeV^2
double bb = 4*AA*(Ee+Epom) - 8*(Pe+Ppom)*(Pe+Ppom)*Ee; //GeV^3
double cc = 4*(Pe+Ppom)*(Pe+Ppom)*Pe*Pe - AA*AA;
sqrtarg = bb*bb-4*aa*cc;
if (sqrtarg < 0) {
//This is used by validUPC as a sign of failure, needs to return something negative
return -1.;
}
return ( -bb + sqrt(sqrtarg) )/(2*aa); //GeV
-}
+}
//
// UPC only
//
bool Kinematics::validUPC(double hBeamEnergy, double eBeamEnergy, double t, double xpom, double vmMass, bool verbose)
{
//
// Here we perform all complete test of all kinematic variables
//
//
// xpom
//
double s = Kinematics::s(eBeamEnergy, hBeamEnergy, true);
double Egam = Egamma(xpom, t, vmMass, hBeamEnergy, eBeamEnergy);
double minxpom = xpomMin(s, vmMass, t);
if (xpom < minxpom) {
- if (verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
+ if(verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
<< " because xpom < xpomMin (" << minxpom << ")" << endl;
return false;
}
if (Egam < 0) {
- if (verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
+ if(verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
<< " because Egamma < 0 (" << Egam << ")" << endl;
return false;
}
if (xpom > 1) {
- if (verbose) cout << "Kinematics::validUPC(): reject xpom="<<xpom<< " because xpom > 1" << endl;
+ if(verbose) cout << "Kinematics::validUPC(): reject xpom="<<xpom<< " because xpom > 1" << endl;
return false;
}
//
// t
//
double t_max=tmax(xpom);
if (t > t_max) {
if (verbose) cout << "Kinematics::validUPC(): reject t=" << t << " because t > tmax (" << t_max << ")" << endl;
return false;
}
return true; //survived all tests
}
-bool Kinematics::valid(double s, double t, double Q2, double W2, double vmMass,
- bool useTrueXp, bool verbose)
-{
- //
- // Here we perform all complete test of all kinematic variables
- //
- double x = Kinematics::x(Q2, W2);
- double y = Kinematics::y(Q2, x, s);
-
- //
- // W
- //
- if (W2 < Kinematics::W2min(vmMass)) {
- if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W < Wmin" << endl;
- return false;
- }
- if (W2 > Kinematics::W2max(s)) {
- if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W > Wmax" << endl;
- return false;
- }
-
- //
- // y range
- //
- double ymin = Kinematics::ymin(s, vmMass);
- double ymax = Kinematics::ymax(s, vmMass);
- if (y < ymin || y > ymax) {
- if (verbose) cout << "Kinematics::valid(): reject y=" << y << " because y < ymin (" << ymin
- << ") || y > ymax (" << ymax << ")" << endl;
- return false;
- }
-
- //
- // x and y
- //
- if (x < 0 || x > 1) {
- if (verbose) cout << "Kinematics::valid(): reject x=" << x << " because x < 0 || x > 1" << endl;
- return false;
- }
-
- //
- // Q2
- //
- if (Q2 > Kinematics::Q2max(s)) {
- if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 > Q2max" << endl;
- return false;
- }
- if (Q2 < Kinematics::Q2min(y)) {
- if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 < Q2min" << endl;
- return false;
- }
-
- //
- // t
- //
- // On calculating xpomeron. xpomeron() provides the correct value
- // (not neglecting t or nucleon mass). However, during initialization
- // to avoid recursion we set t=0. Calling valid() with t != 0 would
- // reject these values. In this case useTrueXp is set to false.
- // For final checks if an event has the correct kinematics we have
- // to set useTrueXp = true.
- //
- double xpom;
- if (useTrueXp)
- xpom = Kinematics::xpomeron(t, Q2, W2, vmMass);
- else
- xpom = Kinematics::xpomeron(0, Q2, W2, vmMass);
- double tmax = Kinematics::tmax(xpom);
- if (t > tmax) {
- if (verbose) cout << "Kinematics::valid(): reject t=" << t << " because t > tmax (" << tmax << ")" << endl;
- return false;
- }
-
-
- //
- // xpom
- //
- if (xpom < 0 || xpom > 1) {
- if (verbose) cout << "Kinematics::valid(): reject xpom=" << xpom << " because xpom < 0 || xpom > 1" << endl;
- return false;
- }
- if (x > xpom) {
- if (verbose) cout << "Kinematics::valid(): reject since x=" << x << " > xpom = " << xpom << endl;
- return false;
- }
-
- return true; // survived all tests
-}
+bool Kinematics::valid(double s, double beta, double Q2, double W2, double z, double MX, bool useTrueXp, bool verbose){
+ double x = Kinematics::x(Q2, W2);
+ double y = Kinematics::y(Q2, x, s);
+ double xpom=x/beta;
+ if(beta<x){
+ if (verbose) cout << "Kinematics::valid(): reject beta=" << beta << " because beta < x=" << x<<endl;
+ return false;
+ }
+ if(beta>1.){
+ if (verbose) cout << "Kinematics::valid(): reject beta=" << beta << " because beta > 1" <<endl;
+ return false;
+ }
+ if(z<0 or z>1)
+ return false;
+ if (W2 < Kinematics::W2min(MX)) {
+ if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W < Wmin" << endl;
+ return false;
+ }
+ if (W2 > Kinematics::W2max(s, Q2, MX)) {
+ if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W > Wmax="<<sqrt(Kinematics::W2max(s, Q2, MX)) << endl;
+ return false;
+ }
+ double ymin = Kinematics::ymin(s, MX);
+ double ymax = Kinematics::ymax(s, MX);
+ if (y < ymin || y > ymax) {
+ if (verbose) cout << "Kinematics::valid(): reject y=" << y << " because y < ymin (" << ymin
+ << ") || y > ymax (" << ymax << ")" << endl;
+ return false;
+ }
+ //
+ // x and y
+ //
+ if (x < 0 || x > 1) {
+ if (verbose) cout << "Kinematics::valid(): reject x=" << x << " because x < 0 || x > 1" << endl;
+ return false;
+ }
+
+ //
+ // Q2
+ //
+ if (Q2 > Kinematics::Q2max(s)) {
+ if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 > Q2max (" << Kinematics::Q2max(s) << ")" << endl;
+ return false;
+ }
+ if (Q2 < Kinematics::Q2min(y)) {
+ if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 < Q2min (" <<Kinematics::Q2min(y)<<")" << endl;
+ return false;
+ }
+ //
+ // xpom
+ //
+ if (xpom < 0 || xpom > 1) {
+ if (verbose) cout << "Kinematics::valid(): reject xpom=" << xpom << " because xpom < 0 || xpom > 1" << endl;
+ return false;
+ }
+ if (x > xpom) {
+ if (verbose) cout << "Kinematics::valid(): reject since x=" << x << " > xpom = " << xpom << endl;
+ return false;
+ }
+
+ return true;
+}
+
+bool Kinematics::valid(double s, double t, double Q2, double W2, double vmMass,
+ bool useTrueXp, bool verbose)
+{
+ //
+ // Here we perform all complete test of all kinematic variables
+ //
+ double x = Kinematics::x(Q2, W2);
+ double y = Kinematics::y(Q2, x, s);
+
+ //
+ // W
+ //
+ if (W2 < Kinematics::W2min(vmMass)) {
+ if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W < Wmin" << endl;
+ return false;
+ }
+ if (W2 > Kinematics::W2max(s)) {
+ if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W > Wmax" << endl;
+ return false;
+ }
+
+ //
+ // y range
+ //
+ double ymin = Kinematics::ymin(s, vmMass);
+ double ymax = Kinematics::ymax(s, vmMass);
+ if (y < ymin || y > ymax) {
+ if (verbose) cout << "Kinematics::valid(): reject y=" << y << " because y < ymin (" << ymin
+ << ") || y > ymax (" << ymax << ")" << endl;
+ return false;
+ }
+
+ //
+ // x and y
+ //
+ if (x < 0 || x > 1) {
+ if (verbose) cout << "Kinematics::valid(): reject x=" << x << " because x < 0 || x > 1" << endl;
+ return false;
+ }
+
+ //
+ // Q2
+ //
+ if (Q2 > Kinematics::Q2max(s)) {
+ if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 > Q2max" << endl;
+ return false;
+ }
+ if (Q2 < Kinematics::Q2min(y)) {
+ if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 < Q2min" << endl;
+ return false;
+ }
+
+ //
+ // t
+ //
+ // On calculating xpomeron. xpomeron() provides the correct value
+ // (not neglecting t or nucleon mass). However, during initialization
+ // to avoid recursion we set t=0. Calling valid() with t != 0 would
+ // reject these values. In this case useTrueXp is set to false.
+ // For final checks if an event has the correct kinematics we have
+ // to set useTrueXp = true.
+ //
+ double xpom;
+ if (useTrueXp)
+ xpom = Kinematics::xpomeron(t, Q2, W2, vmMass);
+ else
+ xpom = Kinematics::xpomeron(0, Q2, W2, vmMass);
+ double tmax = Kinematics::tmax(xpom);
+ if (t > tmax) {
+ if (verbose) cout << "Kinematics::valid(): reject t=" << t << " because t > tmax (" << tmax << ")" << endl;
+ return false;
+ }
+
+
+ //
+ // xpom
+ //
+ if (xpom < 0 || xpom > 1) {
+ if (verbose) cout << "Kinematics::valid(): reject xpom=" << xpom << " because xpom < 0 || xpom > 1" << endl;
+ return false;
+ }
+ if (x > xpom) {
+ if (verbose) cout << "Kinematics::valid(): reject since x=" << x << " > xpom = " << xpom << endl;
+ return false;
+ }
+
+ return true; // survived all tests
+}
ROOT::Math::IBaseFunctionOneDim* LowerXpomeronFormula::Clone() const
{
return new LowerXpomeronFormula();
}
double LowerXpomeronFormula::DoEval(double logxpom) const
{
double xpom = exp(logxpom);
//Start by setting up the knowns: Incoming beams, xpom, and t:
double Ee=mElectronBeam.E();
double Pe=mElectronBeam.Pz();
double Pp=mHadronBeam.Pz();
double Ep=mHadronBeam.E();
double E, pz;
//
// Use scattered proton to set up four momentum of pomeron (so far assuming no break-up):
//
double sqrtarg= 1 - ( (2*xpom - xpom*xpom)*Pp*Pp + mT - mMY2minusM2) / (Ep*Ep);
E=Ep*(1-sqrt(sqrtarg)); //this is close to xpom*Ep
pz=xpom*Pp;
//
// Next we use that the vector meson = photon + pomeron to calculate Egamma
// (keeping the pomeron variables from above)
// We solve a quadratic equation to get the photon energy such that
// the vector meson has the correct mass.
//
double AA = mVmMass2 - E*E - mT + pz*pz + 2*Pe*(Pe+pz); //GeV^2
double aa = 4*(Pe+pz)*(Pe+pz) - 4*(Ee+E)*(Ee+E); //GeV^2
double bb = 4*AA*(Ee+E) - 8*(Pe+pz)*(Pe+pz)*Ee; //GeV^3
double cc = 4*(Pe+pz)*(Pe+pz)*Pe*Pe - AA*AA;
sqrtarg = bb*bb-4*aa*cc;
return sqrtarg;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:27 AM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983028
Default Alt Text
(176 KB)

Event Timeline