Index: trunk/src/TableCollection.cpp
===================================================================
--- trunk/src/TableCollection.cpp (revision 373)
+++ trunk/src/TableCollection.cpp (revision 374)
@@ -1,569 +1,569 @@
//==============================================================================
// TableCollection.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
-// Note that we do not take the lambda_A table into account when calculating
+// Note that we do not take the lambda_real table into account when calculating
// the range since there is a fall back solution to calculate lambda if the
// table is not present. See class CrossSection.
//
//==============================================================================
#include "TSystemDirectory.h"
#include "TSystem.h"
#include "TList.h"
#include "EventGeneratorSettings.h"
#include "TableCollection.h"
#include "Table.h"
#include
#include
#include
#include
#include
#define PR(x) cout << #x << " = " << (x) << endl;
TableCollection::TableCollection() {/* no op */}
TableCollection::TableCollection(int A, DipoleModelType typ, DipoleModelParameterSet set, int vmID)
{
init(A, typ, set, vmID);
}
TableCollection& TableCollection::operator=(const TableCollection& tc)
{
if (this != &tc) {
for (unsigned int i=0; iWorkingDirectory();
//
// Build directory path
//
stringstream pathstream;
pathstream << getenv("SARTRE_DIR") << "/tables/" << A << '/';
if (type == bSat)
pathstream << "bSat";
else if (type == bNonSat)
pathstream << "bNonSat";
else
pathstream << "bCGC";
pathstream << '/' << vmID;
string path = pathstream.str();
//
// Query list of all files in directory
// and create tables for each file ending
// in ".root", ignore others.
//
TSystemDirectory directory;
directory.SetDirectory(path.c_str());
TList *list = directory.GetListOfFiles();
if (!list) {
cout << "TableCollection::init(): Error, cannot find directory '" << path.c_str() << "' holding tables." << endl;
return false;
}
TIter next(list);
unsigned int numberOfTablesRead = 0;
bool upcMode = EventGeneratorSettings::instance()->UPC();
while (TSystemFile* file = dynamic_cast(next()) ) {
if (file->IsDirectory()) continue; // ignore directories
string name(file->GetName());
size_t pos = name.find(".root");
if (pos == string::npos || name.substr(pos) != string(".root")) continue; // ignore files not ending in .root
string fullpath = path + '/' + name;
Table *table = new Table;
if (table->read(fullpath.c_str()) && table->dipoleModelParameterSet() == set) {
if ( (!upcMode && !table->isUPC()) || (upcMode && table->isUPC())) mTables.push_back(table);
numberOfTablesRead++;
if (EventGeneratorSettings::instance()->verboseLevel() > 1)
cout << "Loaded table from file '" << fullpath.c_str() << "'." << endl;
}
}
//
// Cleanup
//
// Change ROOT directory back to directory we were before reading the tables.
// Otherwise reading tables interferes with user application.
//
list->Delete();
gSystem->ChangeDirectory(saveCWD.c_str());
if (!numberOfTablesRead) {
cout << "TableCollection::init(): Error, could not find any tables at '" << path.c_str() << "'" << endl;
cout << " that fit the requested parameters: A=" << A << ", type=" << type
<< ", set=" << set << " , vmID=" << vmID << endl;
return false;
}
return true;
}
bool TableCollection::tableExists(GammaPolarization pol, AmplitudeMoment mom) const
{
Table* currentTable;
for (unsigned int i=0; ipolarization() != pol) continue;
if (currentTable->moment() != mom) continue;
return true;
}
return false;
}
//
// UPC version
//
bool TableCollection::tableExists(AmplitudeMoment mom) const
{
Table* currentTable;
for (unsigned int i=0; imoment() != mom) continue;
return true;
}
return false;
}
bool TableCollection::available(double Q2, double W2, double t, GammaPolarization pol, AmplitudeMoment mom) const
{
//
// Check if table can provide this value
//
unsigned short nTables = 0;
Table* currentTable;
for (unsigned int i=0; ipolarization() != pol) continue;
if (currentTable->moment() != mom) continue;
if (t >= currentTable->minT() && t <= currentTable->maxT()) {
if (Q2 >= currentTable->minQ2() && Q2 <= currentTable->maxQ2()) {
if (W2 >= currentTable->minW2() && W2 <= currentTable->maxW2()) {
nTables++;
}
}
}
}
if (nTables)
return true;
else
return false;
}
bool TableCollection::available(double xpom, double t, AmplitudeMoment mom) const
{
//
// Check if table can provide this value
//
unsigned short nTables = 0;
Table* currentTable;
for (unsigned int i=0; imoment() != mom) continue;
if (t >= currentTable->minT() && t <= currentTable->maxT()) {
if (xpom >= currentTable->minX() && xpom <= currentTable->maxX()) {
nTables++;
}
}
}
if (nTables)
return true;
else
return false;
}
double TableCollection::get(double Q2, double W2, double t,
GammaPolarization pol, AmplitudeMoment mom) const
{
Table *table;
return get(Q2, W2, t, pol, mom, table);
}
//
// UPC version
//
double TableCollection::get(double xpom, double t, AmplitudeMoment mom) const
{
Table *table;
return get(xpom, t, mom, table);
}
double TableCollection::get(double Q2, double W2, double t,
GammaPolarization pol, AmplitudeMoment mom, Table *&table) const
{
static unsigned int errorCount = 0;
const unsigned int maxErrorCount = 10;
//
// First get the tables that contain the necessary info.
// Later this should be a bit refined, here we simply
// loop over all tables to collect the relevant one(s).
//
vector associatedTables;
Table* currentTable;
for (unsigned int i=0; ipolarization() != pol) continue;
if (currentTable->moment() != mom) continue;
if (t >= currentTable->minT() && t <= currentTable->maxT()) {
if (Q2 >= currentTable->minQ2() && Q2 <= currentTable->maxQ2()) {
if (W2 >= currentTable->minW2() && W2 <= currentTable->maxW2()) {
associatedTables.push_back(currentTable);
}
}
}
}
if (associatedTables.size() == 0) {
table = 0;
- if (mom != lambda_A && mom != lambda_skew) { // no warnings needed (can be calculated w/o tables)
+ if (mom != lambda_real && mom != lambda_skew) { // no warnings needed (can be calculated w/o tables)
string txt;
if (mom == mean_A)
txt = "mean_A";
else if (mom == mean_A2)
txt = "mean_A2";
else if (mom == variance_A)
txt = "variance_A";
else
txt = "unknown";
cout << "TableCollection::get(): Warning, could not find any table containing t=" << t
<< ", Q2=" << Q2 << ", W2=" << W2 << endl;
cout << " Tables searched were for moment = " << txt.c_str()
<< ", polarization = " << (pol == transverse ? 'T' : 'L') << endl;
errorCount++;
if (errorCount > maxErrorCount) {
cout << "TableCollection::get(): Error: Too many warnings (>"
<< maxErrorCount
<< "), possibly due to missing table(s)." << endl;
cout << " Stop execution now. Please check the installation of tables." << endl;
exit(1);
}
}
return 0;
}
//
// In case of overlap of tables the following
// policy applies:
// 1. Use the table with the highest priority.
// 2. If there's more than one high priority table
// we average their values (if > 0).
//
unsigned int maxPriority = 0;
for (unsigned int i=0; ipriority() > maxPriority)
maxPriority = associatedTables[i]->priority();
}
double result = 0;
int validCounter = 0;
table = 0;
for (unsigned int i=0; ipriority() == maxPriority) {
double value = associatedTables[i]->get(Q2, W2, t);
if (value > 0) {
validCounter++;
result += value;
table = associatedTables[i];
}
}
}
if (validCounter) result /= validCounter;
return result;
}
//
// UPC version
//
double TableCollection::get(double xpom, double t, AmplitudeMoment mom, Table *&table) const
{
static unsigned int errorCount = 0;
const unsigned int maxErrorCount = 10;
//
// First get the tables that contain the necessary info.
// Later this should be a bit refined, here we simply
// loop over all tables to collect the relevant one(s).
//
vector associatedTables;
Table* currentTable;
for (unsigned int i=0; imoment() != mom) continue;
if (t >= currentTable->minT() && t <= currentTable->maxT()) {
if (xpom >= currentTable->minX() && xpom <= currentTable->maxX()) {
associatedTables.push_back(currentTable);
}
}
}
if (associatedTables.size() == 0) {
table = 0;
- if ( mom != lambda_skew && mom != lambda_A ) { // no warnings needed (can be calculated w/o tables)
+ if ( mom != lambda_skew && mom != lambda_real ) { // no warnings needed (can be calculated w/o tables)
string txt;
if (mom == mean_A)
txt = "mean_A";
else if (mom == mean_A2)
txt = "mean_A2";
else if (mom == variance_A)
txt = "variance_A";
else
txt = "unknown";
cout << "TableCollection::get(): Warning, could not find any table containing t=" << t << ", xp=" << xpom
<< ". Moment = " << txt << ", " << mom << endl;
errorCount++;
if (errorCount > maxErrorCount) {
cout << "TableCollection::get(): Error: Too many warnings (>"
<< maxErrorCount
<< "), possibly due to missing table(s)." << endl;
cout << " Stop execution now. Please check the installation of tables." << endl;
exit(1);
}
}
return 0;
}
//
// In case of overlap of tables the following
// policy applies:
// 1. Use the table with the highest priority.
// 2. If there's more than one high priority table
// we average their values (if > 0).
//
unsigned int maxPriority = 0;
for (unsigned int i=0; ipriority() > maxPriority)
maxPriority = associatedTables[i]->priority();
}
double result = 0;
int validCounter = 0;
table = 0;
double value = 0;
for (unsigned int i=0; ipriority() == maxPriority) {
value = associatedTables[i]->get(xpom, t);
if (value > 0) {
validCounter++;
result += value;
table = associatedTables[i];
}
}
}
if (validCounter) result /= validCounter;
return result;
}
void TableCollection::list(ostream& os, bool opt) const
{
for (unsigned int i=0; ilist(os, opt);
}
double TableCollection::minQ2() const
{
if (EventGeneratorSettings::instance()->UPC()) return 0;
return minimumValue(0);
}
double TableCollection::maxQ2() const
{
if (EventGeneratorSettings::instance()->UPC()) return 0;
return maximumValue(0);
}
double TableCollection::minW2() const
{
if (EventGeneratorSettings::instance()->UPC()) return 0;
return minimumValue(1);
}
double TableCollection::maxW2() const
{
if (EventGeneratorSettings::instance()->UPC()) return 0;
return maximumValue(1);
}
double TableCollection::minW() const {return sqrt(minW2());}
double TableCollection::maxW() const {return sqrt(maxW2());}
double TableCollection::minT() const
{
if (EventGeneratorSettings::instance()->UPC())
return minimumValue(1);
else
return minimumValue(2);
}
double TableCollection::maxT() const
{
if (EventGeneratorSettings::instance()->UPC())
return maximumValue(1);
else
return maximumValue(2);
}
double TableCollection::minX() const
{
if (EventGeneratorSettings::instance()->UPC())
return minimumValue(0);
else
return 0;
}
double TableCollection::maxX() const
{
if (EventGeneratorSettings::instance()->UPC())
return maximumValue(0);
else
return 0;
}
//
// For regular tables: kind: Q2=0, W2=1, T=2
// For UPC tables: kind: xpom = 0, t=1
//
double TableCollection::minimumValue(unsigned int kind) const
{
double minPerTableType[4]; // L, L2, T, T2
fill(minPerTableType, minPerTableType+4, numeric_limits::max());
for (unsigned int i=0; iUPC()) {
switch (kind) {
case (0):
val = mTables[i]->minX();
break;
case (1):
val = mTables[i]->minT();
break;
}
}
else {
switch (kind) {
case (0):
val = mTables[i]->minQ2();
break;
case (1):
val = mTables[i]->minW2();
break;
default:
val = mTables[i]->minT();
break;
}
}
if (mTables[i]->isLongitudinal()) { // L or L2
if (mTables[i]->isMeanA())
minPerTableType[0] = min(minPerTableType[0],val); // L
else
minPerTableType[1] = min(minPerTableType[1],val); // L2
}
else { // T or T2
if (mTables[i]->isMeanA())
minPerTableType[2] = min(minPerTableType[2],val); // T
else
minPerTableType[3] = min(minPerTableType[3],val); // T2
}
}
int startElement = EventGeneratorSettings::instance()->UPC() ? 2 : 0;
double largestMin = *max_element(minPerTableType+startElement, minPerTableType+4);
return largestMin;
}
double TableCollection::maximumValue(unsigned int kind) const
{
double maxPerTableType[4]; // L, L2, T, T2
fill(maxPerTableType, maxPerTableType+4, -numeric_limits::max());
for (unsigned int i=0; iUPC()) {
switch (kind) {
case (0):
val = mTables[i]->maxX();
break;
case (1):
val = mTables[i]->maxT();
break;
}
}
else {
switch (kind) {
case (0):
val = mTables[i]->maxQ2();
break;
case (1):
val = mTables[i]->maxW2();
break;
default:
val = mTables[i]->maxT();
break;
}
}
if (mTables[i]->isLongitudinal()) { // L or L2
if (mTables[i]->isMeanA())
maxPerTableType[0] = max(maxPerTableType[0],val); // L
else
maxPerTableType[1] = max(maxPerTableType[1],val); // L2
}
else { // T or T2
if (mTables[i]->isMeanA())
maxPerTableType[2] = max(maxPerTableType[2],val); // T
else
maxPerTableType[3] = max(maxPerTableType[3],val); // T2
}
}
int startElement = EventGeneratorSettings::instance()->UPC() ? 2 : 0;
double smallestMax = *min_element(maxPerTableType+startElement, maxPerTableType+4);
return smallestMax;
}
Index: trunk/src/Sartre.h
===================================================================
--- trunk/src/Sartre.h (revision 373)
+++ trunk/src/Sartre.h (revision 374)
@@ -1,111 +1,113 @@
//==============================================================================
// Sartre.h
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#ifndef Sartre_h
#define Sartre_h
#include "Event.h"
#include "EventGeneratorSettings.h"
#include "ExclusiveFinalStateGenerator.h"
#include "CrossSection.h"
#include "TableCollection.h"
#include "FrangibleNucleus.h"
#include "Enumerations.h"
#include "TLorentzVector.h"
#include "Math/Functor.h"
#include "TUnuran.h"
#include
#include
#include
using namespace std;
class TUnuranMultiContDist;
class Sartre {
public:
Sartre();
virtual ~Sartre();
virtual bool init(const char* = 0);
virtual bool init(const string&);
virtual Event* generateEvent();
virtual double totalCrossSection(); // in kinematic limits used for generation
virtual double totalCrossSection(double lower[3], double upper[3]); // t, Q2, W
EventGeneratorSettings* runSettings();
const FrangibleNucleus* nucleus() const;
void listStatus(ostream& os=cout) const;
time_t runTime() const;
vector > kinematicLimits(); // t, Q2, W
private:
virtual double calculateTotalCrossSection(double lower[3], double upper[3]); // t, Q2, W2
-
+ virtual bool tableFitsKinematicRange(TableCollection*, AmplitudeMoment, GammaPolarization);
+ virtual bool tableFitsKinematicRange(TableCollection*, AmplitudeMoment); // UPC version
+
private:
Sartre(const Sartre&);
Sartre operator=(const Sartre&);
bool mIsInitialized;
time_t mStartTime;
unsigned long mEvents;
unsigned long mTries;
double mTotalCrossSection;
TLorentzVector mElectronBeam;
TLorentzVector mHadronBeam;
double mS;
unsigned int mA;
int mVmID;
DipoleModelType mDipoleModelType;
DipoleModelParameterSet mDipoleModelParameterSet;
Event *mCurrentEvent;
FrangibleNucleus *mNucleus;
FrangibleNucleus *mUpcNucleus;
TableCollection *mTableCollection;
TableCollection *mProtonTableCollection;
CrossSection *mCrossSection;
EventGeneratorSettings *mSettings;
double mLowerLimit[3]; // t, Q2, W2 (t, xpom for UPC)
double mUpperLimit[3]; // t, Q2, W2
ROOT::Math::Functor *mPDF_Functor;
TUnuran *mUnuran;
TUnuranMultiContDist *mPDF;
unsigned long mEventCounter;
unsigned long mTriesCounter;
ExclusiveFinalStateGenerator mFinalStateGenerator;
};
ostream& operator<<(ostream& os, const TLorentzVector&);
#endif
Index: trunk/src/Enumerations.h
===================================================================
--- trunk/src/Enumerations.h (revision 373)
+++ trunk/src/Enumerations.h (revision 374)
@@ -1,33 +1,33 @@
//==============================================================================
// Enumerations.h
//
// Copyright (C) 2010-2017 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#ifndef Enumerations_h
#define Enumerations_h
enum DipoleModelType {bSat, bNonSat, bCGC};
enum GammaPolarization {transverse, longitudinal};
-enum AmplitudeMoment {mean_A, mean_A2, lambda_A, variance_A, lambda_skew};
+enum AmplitudeMoment {mean_A, mean_A2, variance_A, lambda_real, lambda_skew};
enum DiffractiveMode {coherent, incoherent};
enum DipoleModelParameterSet {KMW, HMPZ, CUSTOM};
enum TableSetType {total_and_coherent, coherent_and_incoherent};
#endif
Index: trunk/src/CrossSection.cpp
===================================================================
--- trunk/src/CrossSection.cpp (revision 373)
+++ trunk/src/CrossSection.cpp (revision 374)
@@ -1,911 +1,921 @@
//==============================================================================
// CrossSection.cpp
//
// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "CrossSection.h"
#include "EventGeneratorSettings.h"
#include "Kinematics.h"
#include "TableCollection.h"
#include "Table.h"
#include "Constants.h"
#include "TH2F.h"
#include "TFile.h"
#include "DglapEvolution.h"
#include
#include
#include
#define PR(x) cout << #x << " = " << (x) << endl;
CrossSection::CrossSection(TableCollection* tc, TableCollection* ptc)
{
mSettings = EventGeneratorSettings::instance();
mRandom = mSettings->randomGenerator();
mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
mPhotonFlux.setS(mS);
mTableCollection = tc;
mProtonTableCollection = ptc;
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mSettings->vectorMesonId());
mVmMass = vectorMesonPDG->Mass();
mCheckKinematics = true;
}
CrossSection::~CrossSection() { /* no op */ }
void CrossSection::setTableCollection(TableCollection* tc) {mTableCollection = tc;}
void CrossSection::setProtonTableCollection(TableCollection* ptc) {mProtonTableCollection = ptc;}
void CrossSection::setCheckKinematics(bool val) {mCheckKinematics = val;}
double CrossSection::operator()(double t, double Q2, double W2)
{
return dsigdtdQ2dW2_total_checked(t, Q2, W2);
}
// UPC Version
double CrossSection::operator()(double t, double xpom)
{
return dsigdtdxp_total_checked(t, xpom);
}
double CrossSection::operator()(const double* array)
{
if (mSettings->UPC())
return dsigdtdxp_total_checked(array[0], array[1]);
else
return dsigdtdQ2dW2_total_checked(array[0], array[1], array[2]);
}
//
// PDF passed to UNURAN
// Array holds: t, log(Q2), W2 for e+p/A running
// t, log(xpom) for UPC
//
double CrossSection::unuranPDF(const double* array) // array is t, log(Q2), W2
{ // or t and log(xpom)
double result = 0;
if (mSettings->UPC()) {
double xpom = exp(array[1]);
result = dsigdtdxp_total_checked(array[0], xpom); // t, xpom
result *= xpom; // Jacobian
}
else {
double Q2 = exp(array[1]);
result = dsigdtdQ2dW2_total_checked(array[0], Q2, array[2]); // t, Q2, W2
result *= Q2; // Jacobian
}
return log(result);
}
double CrossSection::dsigdtdQ2dW2_total_checked(double t, double Q2, double W2)
{
double result = 0;
//
// Check if in kinematically allowed region
//
if (mCheckKinematics && !Kinematics::valid(mS, t, Q2, W2, mVmMass, false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): warning, t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2)
<< " is outside of kinematically allowed region. Return 0." << endl;
return result;
}
//
// Total cross-section dsig2/(dQ2 dW2 dt)
// This is the probability density function needed for UNU.RAN
//
double csT = dsigdtdQ2dW2_total(t, Q2, W2, transverse);
double csL = dsigdtdQ2dW2_total(t, Q2, W2, longitudinal);
result = csT + csL;
//
// Polarization
//
if (mRandom->Uniform(result) <= csT)
mPolarization = transverse;
else
mPolarization = longitudinal;
//
// Diffractive Mode
//
double sampleRange = (mPolarization == transverse ? csT : csL);
double sampleDivider = dsigdtdQ2dW2_coherent(t, Q2, W2, mPolarization);
if (mRandom->Uniform(sampleRange) <= sampleDivider)
mDiffractiveMode = coherent;
else
mDiffractiveMode = incoherent;
//
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): " << result;
cout << " at t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2);
cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal");
cout << " ," << (mDiffractiveMode == coherent ? "coherent" : "incoherent");
cout << ')' << endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, return value = NaN at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, return value = inf at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
if (result < 0) {
cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, negative cross-section at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
return result;
}
double CrossSection::dsigdtdxp_total_checked(double t, double xpom)
{
double result = 0;
//
// Check if in kinematically allowed region
//
if (mCheckKinematics && !Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
t, xpom, mVmMass,
(mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "CrossSection::dsigdtdxp_total_checked(): warning, t=" << t << ", xpom=" << xpom
<< " is outside of kinematically allowed region. Return 0." << endl;
return result;
}
//
// Total cross-section dsig2/(dt dxp)
// This is the probability density function needed for UNU.RAN
//
result = dsigdtdxp_total(t, xpom);
//
// Polarization
//
mPolarization = transverse; // always
//
// Diffractive Mode
//
double sampleRange = result;
double sampleDivider = dsigdtdxp_coherent(t, xpom);
if (mRandom->Uniform(sampleRange) <= sampleDivider)
mDiffractiveMode = coherent;
else
mDiffractiveMode = incoherent;
//
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) {
cout << "CrossSection::dsigdtdxp_total_checked(): " << result;
cout << " at t=" << t << ", xp=" << xpom;
cout << " (" << (mDiffractiveMode == coherent ? "coherent" : "incoherent");
cout << ')' << endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "CrossSection::dsigdtdxp_total_checked(): Error, return value = NaN at" << endl;
cout << " t=" << t << ", xp=" << xpom << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "CrossSection::dsigdtdxp_total_checked(): Error, return value = inf at" << endl;
cout << " t=" << t << ", xp=" << xpom << endl;
result = 0;
}
if (result < 0) {
cout << "CrossSection::dsigdtdxp_total_checked(): Error, negative cross-section at" << endl;
cout << " t=" << t << ", xp=" << xpom << endl;
result = 0;
}
return result;
}
double CrossSection::dsigdt_total(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = 0;
if (mSettings->tableSetType() == coherent_and_incoherent) {
result = dsigdt_coherent(t, Q2, W2, pol) + dsigdt_incoherent(t, Q2, W2, pol);
}
else if (mSettings->tableSetType() == total_and_coherent) {
result = mTableCollection->get(Q2, W2, t, pol, mean_A2); // units now are fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()){
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()){
lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
result *= skewednessCorrection(lambda);
}
}
return result;
}
//
// UPC Version
//
double CrossSection::dsigdt_total(double t, double xpom) const
{
double result = 0;
if (mSettings->tableSetType() == coherent_and_incoherent) {
result = dsigdt_coherent(t, xpom) + dsigdt_incoherent(t, xpom);
}
else if (mSettings->tableSetType() == total_and_coherent) {
result = mTableCollection->get(xpom, t, mean_A2); // units now are fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
- if (mSettings->correctForRealAmplitude()){
- lambda = logDerivateOfAmplitude(t, xpom);
+ if (mSettings->correctForRealAmplitude()) {
+ lambda = logDerivateOfAmplitude(t, xpom);
result *= realAmplitudeCorrection(lambda);
- }
- if (mSettings->correctSkewedness()){
- lambda = logDerivateOfGluonDensity(t, xpom);
+ }
+ if (mSettings->correctSkewedness()) {
+ lambda = logDerivateOfGluonDensity(t, xpom);
result *= skewednessCorrection(lambda);
- }
+ }
}
return result;
}
double CrossSection::dsigdt_coherent(double t, double Q2, double W2, GammaPolarization pol) const
{
double val = mTableCollection->get(Q2, W2, t, pol, mean_A); // fm^2
- double result = val*val/(16*M_PI); // units now are fm^4
- result /= hbarc2; // in fm^2/GeV^2
- result *= 1e7; // in nb/GeV^2
+ double result = val*val/(16*M_PI); // units now are fm^4
+ result /= hbarc2; // in fm^2/GeV^2
+ result *= 1e7; // in nb/GeV^2
double lambda = 0;
- if (mSettings->correctForRealAmplitude()){
- lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
- result *= realAmplitudeCorrection(lambda);
- }
- if (mSettings->correctSkewedness()){
- lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
- result *= skewednessCorrection(lambda);
+ if (mSettings->correctForRealAmplitude()) {
+ lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
+ result *= realAmplitudeCorrection(lambda);
+ }
+ if (mSettings->correctSkewedness()) {
+ lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
+ result *= skewednessCorrection(lambda);
}
return result;
}
//
// UPC Version
//
double CrossSection::dsigdt_coherent(double t, double xpom) const
{
double val = mTableCollection->get(xpom, t, mean_A); // fm^2
double result = val*val/(16*M_PI); // units now are fm^4
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
- if (mSettings->correctForRealAmplitude()){
- lambda = logDerivateOfAmplitude(t, xpom);
- result *= realAmplitudeCorrection(lambda);
- }
- if (mSettings->correctSkewedness()){
- lambda = logDerivateOfGluonDensity(t, xpom);
- result *= skewednessCorrection(lambda);
+ if (mSettings->correctForRealAmplitude()) {
+ lambda = logDerivateOfAmplitude(t, xpom);
+ result *= realAmplitudeCorrection(lambda);
+ }
+ if (mSettings->correctSkewedness()) {
+ lambda = logDerivateOfGluonDensity(t, xpom);
+ result *= skewednessCorrection(lambda);
}
return result;
}
double CrossSection::dsigdtdQ2dW2_total(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = dsigdt_total(t, Q2, W2, pol);
if (mSettings->applyPhotonFlux()) result *= mPhotonFlux(Q2,W2,pol);
return result;
}
//
// UPC version
//
double CrossSection::dsigdtdxp_total(double t, double xpom) const
{
double result = dsigdt_total(t, xpom);
if (mSettings->applyPhotonFlux()) {
result *= UPCPhotonFlux(t, xpom);
}
return result;
}
double CrossSection::dsigdt_incoherent(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = mTableCollection->get(Q2, W2, t, pol, variance_A); // fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
if (mSettings->correctForRealAmplitude()){
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
result *= realAmplitudeCorrection(lambda);
}
if (mSettings->correctSkewedness()){
lambda = logDerivateOfGluonDensity(t, Q2, W2, pol);
result *= skewednessCorrection(lambda);
}
return result;
}
//
// UPC Version
//
double CrossSection::dsigdt_incoherent(double t, double xpom) const
{
double result = mTableCollection->get(xpom, t, variance_A); // fm^4
result /= (16*M_PI);
result /= hbarc2; // in fm^2/GeV^2
result *= 1e7; // in nb/GeV^2
double lambda = 0;
- if (mSettings->correctForRealAmplitude()){
- lambda = logDerivateOfAmplitude(t, xpom);
- result *= realAmplitudeCorrection(lambda);
- }
- if (mSettings->correctSkewedness()){
- lambda = logDerivateOfGluonDensity(t, xpom);
- result *= skewednessCorrection(lambda);
+ if (mSettings->correctForRealAmplitude()) {
+ lambda = logDerivateOfAmplitude(t, xpom);
+ result *= realAmplitudeCorrection(lambda);
+ }
+ if (mSettings->correctSkewedness()) {
+ lambda = logDerivateOfGluonDensity(t, xpom);
+ result *= skewednessCorrection(lambda);
}
return result;
}
double CrossSection::dsigdtdQ2dW2_coherent(double t, double Q2, double W2, GammaPolarization pol) const
{
double result = dsigdt_coherent(t, Q2, W2, pol);
if (mSettings->applyPhotonFlux()) result *= mPhotonFlux(Q2,W2,pol);
return result;
}
//
// UPC Version
//
double CrossSection::dsigdtdxp_coherent(double t, double xpom) const
{
double result = dsigdt_coherent(t, xpom);
if (mSettings->applyPhotonFlux()) {
result *= UPCPhotonFlux(t, xpom);
}
return result;
}
double CrossSection::UPCPhotonFlux(double t, double xpom) const{
double Egamma = Kinematics::Egamma(xpom, t, mVmMass,
mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy());
double result = mPhotonFlux(Egamma);
//
// Jacobian = dEgamma/dN, such that dsig/dtdxp = dsigd/tdEgam*dEgam/dxp
//
double h=xpom*1e-3;
double Egamma_p = Kinematics::Egamma(xpom+h, t, mVmMass,
mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy());
double Egamma_m = Kinematics::Egamma(xpom-h, t, mVmMass,
mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy());
double jacobian = (Egamma_p-Egamma_m)/(2*h);
result *= fabs(jacobian);
return result;
}
GammaPolarization CrossSection::polarizationOfLastCall() const {return mPolarization;}
DiffractiveMode CrossSection::diffractiveModeOfLastCall() const {return mDiffractiveMode;}
double CrossSection::logDerivateOfAmplitude(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfAmplitude(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
// Note (TU): if the lambda value from a table is not valid and not > 0,
// get() returns lambda=0 and table=0. This enforces a renewed
// calculation.
//
- lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_A, table);
+ lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_real, table);
if (!table) { // no lambda value from correct table, use fallback solution
lambdaFromTable = false;
double value = mProtonTableCollection->get(Q2, W2, t, pol, mean_A, table); // use obtained table from here on
- /*
if (value <= 0) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << value << '.' << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", pol="
<< (pol == transverse ? 'T' : 'L') << endl;
return 0;
}
- */
+
if (!table) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid pointer to lookup table." << endl;
return 0;
}
//
// Note: the derivate taken from values in the table is a delicate issue.
// The standard interpolation method(s) used in Table::get() are
// at times not accurate and causes ripples on lambda(W2).
// However, interpolation methods or robust fitting is by far too
- // expensive in terms of time per point.
+ // expensive in terms of CPU time per point.
//
//
// Calculate the derivative using simple method.
//
double derivative;
double hplus, hminus;
double dW2 = table->binWidthW2();
hplus = hminus = 0.5*dW2; // half bin width is found to be the best choice after some testing
hplus = min(hplus, table->maxW2()-W2);
hminus = min(hminus, W2-table->minW2());
hminus -= numeric_limits::epsilon();
hplus -= numeric_limits::epsilon();
if (hminus < 0) hminus = 0;
if (hplus < 0) hplus = 0;
if (hplus == 0 && hminus == 0) {
cout << "CrossSection::logDerivateOfAmplitude(): Warning, cannot find derivative." << endl;
return 0;
}
double a = table->get(Q2, W2+hplus, t);
+ if (a <= 0) {
+ cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << a << '.' << endl;
+ cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2+hplus) << ", pol="
+ << (pol == transverse ? 'T' : 'L') << endl;
+ return 0;
+ }
double b = table->get(Q2, W2-hminus, t);
+ if (b <= 0) {
+ cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << b << '.' << endl;
+ cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2-hminus) << ", pol="
+ << (pol == transverse ? 'T' : 'L') << endl;
+ return 0;
+ }
derivative = (a-b)/(hplus+hminus);
//
// Finally calculate lambda
//
double jacobian = (W2-protonMass2+Q2)/value;
lambda = jacobian*derivative;
} // end fall back solution
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
if (lambdaFromTable)
cout << "Info, lambda taken from table." << endl;
else
cout << "Info, lambda derived numerically from proton amplitude table" << endl;
}
//
+ // Check lambda value.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
- // In all cases where something goes wrong we set lambda to
- // its most probably value (0.2)
- //
+ //
- const double lambdaInCaseOfError = 0.2;
-
-
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
- /*
- if (lambda < 0) {
- if (mSettings->verboseLevel() > 2) {
- cout << "CrossSection::logDerivateOfAmplitude(): ";
- cout << "Warning, lambda is negative (" << lambda << "). " ;
- cout << "Set to " << (-lambda) << "." << endl;
- }
- lambda = -lambda;
- }
- */
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
- lambda = lambdaInCaseOfError;
- }
- if (std::isnan(lambda)) {
- cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl;
- cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
- lambda = lambdaInCaseOfError;
- }
- return lambda;
-}
-
-double CrossSection::logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization pol) const
-{
- double lambda = 0;
- bool lambdaFromTable = true;
- Table *table = 0;
-
- if (!mProtonTableCollection) {
- cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl;
- cout << " Corrections not available. Should be off." << endl;
+ cout << " Set to 0." << endl;
return 0;
}
-
- //
- // If the lambda table is present we use the more accurate and numerically
- // stable table value. Otherwise we calculate it from the table(s).
- //
- lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_skew, table);
-
- if (!table) {
- //
- // no lambda value from correct table, use fallback solution
- // lambda_skew ~ lambda_real
- //
- lambda = logDerivateOfAmplitude(t, Q2, W2, pol) ; // This is an approximation in this case
- } // end fall back solution
-
- if (mSettings->verboseLevel() > 3) {
- cout << "CrossSection::logDerivateOfGluonDensity(): ";
- if (lambdaFromTable)
- cout << "Info, lambda taken from table." << endl;
- else
- cout << "Info, lambda taken from logDerivateOfAmplitude as an approximation" << endl;
- }
- //
- // At a lambda of ~0.6 both corrections have equal value
- // of around 2.9. This will yield excessive large (unphysical)
- // corrections. Large values are typically caused by fluctuations
- // and glitches in the tables and should be rare.
- // In all cases where something goes wrong we set lambda to
- // its most probably value (0.2)
- //
-
- const double lambdaInCaseOfError = 0.2;
-
-
- double maxLambda = mSettings->maxLambdaUsedInCorrections();
- if (fabs(lambda) > maxLambda) {
- if (mSettings->verboseLevel() > 2) {
- cout << "CrossSection::logDerivateOfGluonDensity(): ";
- cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
- cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl;
- cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
- }
- lambda = lambda > 0 ? maxLambda : -maxLambda;
- }
-
- if (std::isinf(lambda)) {
- cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl;
- cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
- lambda = lambdaInCaseOfError;
- }
if (std::isnan(lambda)) {
- cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl;
+ cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
- lambda = lambdaInCaseOfError;
+ cout << " Set to 0." << endl;
+ return 0;
}
-
return lambda;
}
//
// UPC only version
//
double CrossSection::logDerivateOfAmplitude(double t, double xpom) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfAmplitude(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// Usual numeric issues at boundaries.
// Subtracting an eps in log(xpom) does the trick.
//
if (xpom > mProtonTableCollection->maxX()) {
- xpom = exp(log(xpom)-numeric_limits::epsilon());
+ xpom = exp(log(xpom)-numeric_limits::epsilon());
}
if (xpom < mProtonTableCollection->minX()) {
xpom = exp(log(xpom)+numeric_limits::epsilon());
}
-
+
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
- lambda = mProtonTableCollection->get(xpom, t, lambda_A, table);
-
+ lambda = mProtonTableCollection->get(xpom, t, lambda_real, table);
+
if (!table) { // no lambda value from correct table, use fallback solution
lambdaFromTable = false;
- (void) mProtonTableCollection->get(xpom, t, mean_A, table); // use obtained table from here on
+ (void) mProtonTableCollection->get(xpom, t, mean_A, table); // use obtained table from here on
if (!table) {
cout << "CrossSection::logDerivateOfAmplitude(): got invalid pointer to lookup table." << endl;
return 0;
}
-
+
//
// Here's the tricky part (see comments in non-UPC version above)
//
double theLogxpom = log(xpom);
double dlogxpom = table->binWidthX(); // assuming table is in log x ???????? FIX later
double maxLogxpom = log(table->maxX());
double derivative;
double hplus, hminus;
hplus = hminus = 0.5*dlogxpom;
hplus = min(hplus, fabs(maxLogxpom-theLogxpom));
hminus = min(hminus, fabs(theLogxpom-log(table->minX())));
hminus -= numeric_limits::epsilon();
hplus -= numeric_limits::epsilon();
if (hminus < 0) hminus = 0;
if (hplus < 0) hplus = 0;
if (hplus == 0 && hminus == 0) {
cout << "CrossSection::logDerivateOfAmplitude(): Warning, cannot find derivative." << endl;
return 0;
}
double a = table->get(exp(theLogxpom+hplus), t);
+ if (a <= 0) {
+ cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << a << '.' << endl;
+ cout << " t=" << t << ", W=" << sqrt(exp(theLogxpom+hplus)) << endl;
+ return 0;
+ }
double b = table->get(exp(theLogxpom-hminus), t);
- derivative = log(a/b)/(hplus+hminus);
-
+ if (a <= 0) {
+ cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << b << '.' << endl;
+ cout << " t=" << t << ", W=" << sqrt(exp(theLogxpom-hminus)) << endl;
+ return 0;
+ }
+ derivative = log(a/b)/(hplus+hminus);
+
//
// Finally calculate lambda
// Directly dlog(A)/-dlog(xpom) here.
- //
+ //
lambda = -derivative;
}
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
if (lambdaFromTable)
- cout << "Info, lambda taken from table." << endl;
+ cout << "Info, lambda taken from table." << endl;
else
- cout << "Info, lambda derived numerically from proton amplitude table" << endl;
- cout << " t=" << t << ", xpom=" << xpom << endl;
+ cout << "Info, lambda derived numerically from proton amplitude table" << endl;
+ cout << " t=" << t << ", xpom=" << xpom << endl;
}
//
+ // Check lambda value.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
- // In all cases where something goes wrong we set lambda to
- // its most probably value (0.2)
//
- const double lambdaInCaseOfError = 0.2;
-
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "xpom=" << xpom << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
-
+
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda = infinity for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
- lambda = lambdaInCaseOfError;
+ cout << " Set to 0." << endl;
+ return 0;
}
+
if (std::isnan(lambda)) {
cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
- lambda = lambdaInCaseOfError;
+ cout << " Set to 0." << endl;
+ return 0;
}
-
+
return lambda;
}
+
+double CrossSection::logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization pol) const
+{
+ double lambda = 0;
+ bool lambdaFromTable = true;
+ Table *table = 0;
+
+ if (!mProtonTableCollection) {
+ cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl;
+ cout << " Corrections not available. Should be off." << endl;
+ return 0;
+ }
+
+ //
+ // If the lambda table is present we use the more accurate and numerically
+ // stable table value. Otherwise we calculate it from the table(s).
+ //
+ lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_skew, table);
+
+ if (!table) {
+ //
+ // no lambda value from correct table, use fallback solution
+ // lambda_skew ~ lambda_real
+ //
+ lambda = logDerivateOfAmplitude(t, Q2, W2, pol) ; // This is an approximation in this case
+ } // end fall back solution
+
+ if (mSettings->verboseLevel() > 3) {
+ cout << "CrossSection::logDerivateOfGluonDensity(): ";
+ if (lambdaFromTable)
+ cout << "Info, lambda taken from table." << endl;
+ else
+ cout << "Info, lambda taken from logDerivateOfAmplitude as an approximation" << endl;
+ }
+
+ //
+ // Checking lambda.
+ // At a lambda of ~0.6 both corrections have equal value
+ // of around 2.9. This will yield excessive large (unphysical)
+ // corrections. Large values are typically caused by fluctuations
+ // and glitches in the tables and should be rare.
+ //
+
+ double maxLambda = mSettings->maxLambdaUsedInCorrections();
+ if (fabs(lambda) > maxLambda) {
+ if (mSettings->verboseLevel() > 2) {
+ cout << "CrossSection::logDerivateOfGluonDensity(): ";
+ cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
+ cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl;
+ cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
+ }
+ lambda = lambda > 0 ? maxLambda : -maxLambda;
+ }
+
+ if (std::isinf(lambda)) {
+ cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl;
+ cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
+ cout << " Set to 0." << endl;
+ return 0;
+ }
+
+ if (std::isnan(lambda)) {
+ cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl;
+ cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
+ cout << " Set to 0." << endl;
+ return 0;
+ }
+
+ return lambda;
+}
+
//
// UPC version
//
double CrossSection::logDerivateOfGluonDensity(double t, double xpom) const
{
double lambda = 0;
bool lambdaFromTable = true;
Table *table = 0;
if (!mProtonTableCollection) {
cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl;
cout << " Corrections not available. Should be off." << endl;
return 0;
}
//
// Usual numeric issues at boundaries.
// Subtracting an eps in log(xpom) does the trick.
//
if (xpom > mProtonTableCollection->maxX()) {
xpom = exp(log(xpom)-numeric_limits::epsilon());
}
if (xpom < mProtonTableCollection->minX()) {
xpom = exp(log(xpom)+numeric_limits::epsilon());
}
-
+
//
// If the lambda table is present we use the more accurate and numerically
// stable table value. Otherwise we calculate it from the table(s).
//
lambda = mProtonTableCollection->get(xpom, t, lambda_skew, table);
if (!table) {
//
// no lambda value from correct table, use fallback solution
// lambda_skew ~ lambda_real
//
-
+
lambdaFromTable = false;
-
- lambda = logDerivateOfAmplitude(t, xpom);
+
+ lambda = logDerivateOfAmplitude(t, xpom);
}
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
if (lambdaFromTable)
- cout << "Info, lambda taken from table." << endl;
+ cout << "Info, lambda taken from table." << endl;
else
- cout << "Info, lambda taken from logDerivateOfAmplitude as approximation." << endl;
+ cout << "Info, lambda taken from logDerivateOfAmplitude as approximation." << endl;
}
//
+ // Check lambda value.
// At a lambda of ~0.6 both corrections have equal value
// of around 2.9. This will yield excessive large (unphysical)
// corrections. Large values are typically caused by fluctuations
// and glitches in the tables and should be rare.
- // In all cases where something goes wrong we set lambda to
- // its most probably value (0.2)
//
- const double lambdaInCaseOfError = 0.2;
-
double maxLambda = mSettings->maxLambdaUsedInCorrections();
if (fabs(lambda) > maxLambda) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
cout << "Warning, lambda is excessively large (" << lambda << ") at " ;
cout << "xpom=" << xpom << ", t=" << t << endl;
cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl;
}
lambda = lambda > 0 ? maxLambda : -maxLambda;
}
-
+
if (std::isinf(lambda)) {
cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = infinity for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
- lambda = lambdaInCaseOfError;
+ cout << " Set to 0." << endl;
+ return 0;
}
+
if (std::isnan(lambda)) {
cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for" << endl;
cout << " t=" << t << ", xpom=" << xpom << endl;
- lambda = lambdaInCaseOfError;
+ cout << " Set to 0." << endl;
+ return 0;
}
-
+
return lambda;
}
double CrossSection::realAmplitudeCorrection(double lambda) const
{
//
// Correction factor for real amplitude contribution
//
double beta = tan(lambda*M_PI/2.);
double correction = 1 + beta*beta;
return correction;
}
double CrossSection::realAmplitudeCorrection(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
double correction = realAmplitudeCorrection(lambda);
// correction *= exp(-10*Kinematics::x(Q2, W2)); // damped
return correction;
}
//
// UPC version
//
double CrossSection::realAmplitudeCorrection(double t, double xpom) const
{
double lambda = logDerivateOfAmplitude(t, xpom);
double correction = realAmplitudeCorrection(lambda);
return correction;
}
double CrossSection::skewednessCorrection(double lambda) const
{
//
// Skewedness correction
//
double R = pow(2.,2*lambda+3)*TMath::Gamma(lambda+2.5)/(sqrt(M_PI)*TMath::Gamma(lambda+4));
double correction = R*R;
return correction;
}
double CrossSection::skewednessCorrection(double t, double Q2, double W2, GammaPolarization pol) const
{
double lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
double correction = skewednessCorrection(lambda);
// correction *= exp(-10*Kinematics::x(Q2, W2)); // damped
return correction;
}
//
// UPC version
//
double CrossSection::skewednessCorrection(double t, double xpom) const
{
double lambda = logDerivateOfAmplitude(t, xpom);
double correction = skewednessCorrection(lambda);
return correction;
}
Index: trunk/src/Sartre.cpp
===================================================================
--- trunk/src/Sartre.cpp (revision 373)
+++ trunk/src/Sartre.cpp (revision 374)
@@ -1,1012 +1,1105 @@
//==============================================================================
// Sartre.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Note:
// When not using runcards, the user must first create an instance of Sartre
// then get the settings via one of:
// Sartre::runSettings()
// EventGeneratorSettings::instance()
// Once init() is called settings cannot be changed any more.
//==============================================================================
#include "Version.h"
#include "Kinematics.h"
#include "Sartre.h"
#include "ModeFinderFunctor.h"
#include "Math/BrentMinimizer1D.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/GSLMinimizer.h"
#include "Math/Functor.h"
#include "TUnuranMultiContDist.h"
#include
#include
#include
#include "TH2D.h"
#include "TFile.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
Sartre::Sartre()
{
mIsInitialized = false;
mCurrentEvent = 0;
mNucleus = 0;
mUpcNucleus= 0;
mSettings = 0;
mPDF_Functor = 0;
mPDF = 0;
mEventCounter = 0;
mTriesCounter = 0;
mTotalCrossSection = 0;
mCrossSection = 0;
mTableCollection = 0;
mProtonTableCollection = 0;
mUnuran = 0;
mEvents = 0;
mTries = 0;
mS = 0;
mA = 0;
}
Sartre::~Sartre()
{
delete mNucleus;
delete mUpcNucleus;
delete mPDF_Functor;
delete mPDF;
delete mCrossSection;
if (mTableCollection != mProtonTableCollection) {
delete mTableCollection;
delete mProtonTableCollection;
}
else
delete mTableCollection;
delete mUnuran;
}
bool Sartre::init(const char* runcard)
{
mStartTime = time(0);
bool ok;
//
// Reset member variables.
// Note that one instance of Sartre should be able to get
// initialized multiple times.
//
mEvents = 0;
mTries = 0;
mTotalCrossSection = 0;
//
// Print header
//
string ctstr(ctime(&mStartTime));
ctstr.erase(ctstr.size()-1, 1);
cout << "/========================================================================\\" << endl;
cout << "| |" << endl;
cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
cout << "| |" << endl;
cout << "| An event generator for exclusive diffractive vector meson production |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
cout << "| Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich |" << endl;
cout << "| |" << endl;
cout << "| This program is free software: you can redistribute it and/or modify |" << endl;
cout << "| it under the terms of the GNU General Public License as published by |" << endl;
cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl;
cout << "| any later version. |" << endl;
cout << "| |" << endl;
cout << "| Code compiled on " << setw(12) << left << __DATE__;
cout << setw(41) << left << __TIME__ << right << '|' << endl;
cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
cout << "\\========================================================================/" << endl;
mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton
//
// Read runcard if available
//
if (runcard) {
if (!mSettings->readSettingsFromFile(runcard)) {
cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl;
exit(1);
}
else
cout << "Runcard is '" << runcard << "'." << endl;
}
else
cout << "No runcard provided." << endl;
//
// Set beam particles and center of mass energy
// Note, that if we are in UPC mode the electron
// is actually treated as the hadron beam, i.e.
// the mass is the proton mass.
//
mElectronBeam = mSettings->eBeam();
mHadronBeam = mSettings->hBeam();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mA = mSettings->A();
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (allowBreakup) {
if (!getenv("SARTRE_DIR")) {
cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n"
"to locate tables needed for the generation if nuclear breakups." << endl;
exit(1);
}
}
if (mNucleus) delete mNucleus;
mNucleus = new FrangibleNucleus(mA, allowBreakup);
if (mSettings->UPC()) {
if (mUpcNucleus) delete mUpcNucleus;
mUpcNucleus = new FrangibleNucleus(mSettings->UPCA(), allowBreakup);
}
string upcNucleusName;
if (mSettings->UPC()) {
upcNucleusName = mUpcNucleus->name();
cout << "Sartre is running in UPC mode" << endl;
cout << "Hadron 1 beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron 1 beam: " << mHadronBeam << endl;
cout << "Hadron 2 beam species: " << upcNucleusName << " (" << mSettings->UPCA() << ")" << endl;
cout << "Hadron 2 beam: " << mElectronBeam << endl;
}
else {
cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron beam: " << mHadronBeam << endl;
cout << "Electron beam: " << mElectronBeam << endl;
}
//
// Get details about the processes and models
//
mDipoleModelType = mSettings->dipoleModelType();
mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
mVmID = mSettings->vectorMesonId();
cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
cout << "Process is ";
if (mSettings->UPC()) {
cout << mNucleus->name() << " + " << upcNucleusName << " -> "
<< mNucleus->name() << "' + " << upcNucleusName << "' + ";
}
else {
if (mA > 1)
cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + ";
else
cout << "e + p -> e' + p' + ";
}
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID);
cout << vectorMesonPDG->GetName() << endl;
//
// Print-out seed for reference
//
cout << "Random generator seed: " << mSettings->seed() << endl;
//
// Load in the tables containing the amplitude moments
//
if (!getenv("SARTRE_DIR")) {
- cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
+ cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
exit(1);
}
if (mTableCollection) delete mTableCollection;
mTableCollection = new TableCollection;
ok = mTableCollection->init(mA, mDipoleModelType, mDipoleModelParameterSet, mVmID);
if (!ok) {
- cout << "Error, could not initialize lookup tables for requested process." << endl;
+ cout << "Error, could not initialize lookup tables for requested process." << endl;
return false;
}
//
// Load in the p tables for the lambda lookup tables (or to calculate lambda if not available)
//
if (mSettings->correctForRealAmplitude() || mSettings->correctSkewedness()) {
if (mA == 1) {
mProtonTableCollection = mTableCollection;
}
else {
if (mProtonTableCollection) delete mProtonTableCollection;
mProtonTableCollection = new TableCollection;
ok = mProtonTableCollection->init(1, mDipoleModelType, mDipoleModelParameterSet, mVmID);
if (!ok) {
- cout << "Error: could not initialize proton lookup tables for requested process." << endl;
+ cout << "Error: could not initialize proton lookup tables for requested process." << endl;
cout << " These tables are needed for skewedness and real amplitude corrections." << endl;
return false;
}
}
}
else
mProtonTableCollection = 0;
//
// Kinematic limits and generator range
//
// There are 3 ranges we have to deal with
// 1. the kinematic range requested by the user
// if given.
// The user can only control Q2 and W but not t.
// For UPC that's xpom.
// 2. the range of the table(s)
// 3. the kinematically allowed range
//
// Of course (3) is more complex than a simple cube/square.
// However, we deal with the detailed shape of the kinematic
// range later using Kinematics::valid() when we generate the
// individual events.
// For setting up UNU.RAN we have to get the cubic/square
// envelope that satifies (1)-(3).
// Note, that they are correlated which makes the order
// in which we do things a bit tricky.
//
//
// Step 1:
// Set the limits to that of the table(s).
// Note, the indices 0-2 refer to t, Q2, and W2.
// For UPC we have only t and xpom.
//
if (mSettings->UPC()) {
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minX(); mUpperLimit[1] = mTableCollection->maxX();
mLowerLimit[2] = mUpperLimit[2] = 0;
}
else {
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minQ2(); mUpperLimit[1] = mTableCollection->maxQ2();
mLowerLimit[2] = mTableCollection->minW2(); mUpperLimit[2] = mTableCollection->maxW2();
}
//
// Step 2:
// Kinematic limits might overrule boundaries from step 1
//
if (mSettings->UPC()) {
double kineXpomMin = Kinematics::xpomMin(vectorMesonPDG->Mass(), mUpperLimit[0], mHadronBeam, mElectronBeam);
double kineXpomMax = 1;
mLowerLimit[1] = max(kineXpomMin, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineXpomMax);
double kineTmax = Kinematics::tmax(kineXpomMin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
}
else {
// double kineYmax = Kinematics::ymax(mS, vectorMesonPDG->Mass());
double kineYmin = Kinematics::ymin(mS, vectorMesonPDG->Mass());
double kineQ2min = Kinematics::Q2min(kineYmin);
double kineQ2max = Kinematics::Q2max(mS);
double kineW2min = Kinematics::W2min(vectorMesonPDG->Mass());
double kineW2max = Kinematics::W2max(mS);
kineQ2min = max(kineQ2min, mLowerLimit[1]);
kineQ2max = min(kineQ2max, mUpperLimit[1]);
kineW2min = max(kineW2min, mLowerLimit[2]);
kineW2max = min(kineW2max, mUpperLimit[2]);
double kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // first arg (t) set to 0 (recursive)
double kineTmax = Kinematics::tmax(kineXPmin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
mLowerLimit[1] = max(kineQ2min, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineQ2max);
mLowerLimit[2] = max(kineW2min, mLowerLimit[2]); mUpperLimit[2] = min(mUpperLimit[2], kineW2max);
}
//
// Step 3:
// Deal with user provided limits.
// User settings are ignored (switched off) if min >= max.
//
if (mSettings->UPC()) {
if (mSettings->xpomMin() < mSettings->xpomMax()) {
if (mSettings->xpomMin() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of xpomeron (" << mSettings->xpomMin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->xpomMin();
}
if (mSettings->xpomMax() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of xpomeron (" << mSettings->xpomMax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->xpomMax();
}
}
}
else {
if (mSettings->W2min() < mSettings->W2max()) { // W2 first
if (mSettings->W2min() < mLowerLimit[2]) {
cout << "Warning, requested lower limit of W (" << mSettings->Wmin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << sqrt(mLowerLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[2] = mSettings->W2min();
}
if (mSettings->W2max() > mUpperLimit[2]) {
cout << "Warning, requested upper limit of W (" << mSettings->Wmax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << sqrt(mUpperLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[2] = mSettings->W2max();
}
}
if (mSettings->Q2min() < mSettings->Q2max()) { // Q2
if (mSettings->Q2min() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of Q2 (" << mSettings->Q2min() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->Q2min();
// ????? kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // new Q2min changes tmax
// ????? mUpperLimit[0] = min(mUpperLimit[0], Kinematics::tmax(kineXPmin));
}
if (mSettings->Q2max() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of Q2 (" << mSettings->Q2max() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->Q2max();
}
}
}
//
// Check if any phase space is left
//
if (mLowerLimit[0] >= mUpperLimit[0]) {
cout << "Invalid range in t: t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
exit(1);
}
if (mLowerLimit[1] >= mUpperLimit[1]) {
if (mSettings->UPC())
cout << "Invalid range in xpomeron: xpomeron=[";
else
cout << "Invalid range in Q2: Q2=[";
cout << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl;
exit(1);
}
if (!mSettings->UPC() && mLowerLimit[2] >= mUpperLimit[2]) {
cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl;
exit(1);
}
//
- // Print-out limits (all verbose levels)
+ // Print-out limits (all verbose levels)
//
if (mSettings->verbose()) {
cout << "Kinematic limits used for event generation:" << endl;
if (mSettings->UPC()) {
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "xpom=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
}
else {
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl;
}
}
//
- // Check if the lambda table covers the kinematic range
- // or if we have to calculate lambda from the p mean_A table
- // directly.
- // This is just to inform the user, no action is taken.
- // Here we assume that transverse and longitudinal tables
- // always have the same range and vheck only one.
- //
- // If even the p mean_A table is not big enough we switch
- // corrections off and inform the user.
+ // Skewedness and real part of amplitude corrections rely
+ // on parameters that are stored in tables. Here we check
+ // if they exist and if the tables have the right size.
+ // This is mostly to inform the user what's going on
+ // for verbose level 2 and higher.
+ //
+ // Should they not exist or are too small they are calculated
+ // on the fly (see class CrossSection for details).
+ // If this happens we also check if the kinematic range
+ // of the proton amplitude tables is big enough to allow
+ // calculating lambda for the actual kinematic range.
+ // If not we stop execution.
//
if (mProtonTableCollection) {
- bool lambdaTableExists, tooSmallLambda, tooSmallMean;
-
if (mSettings->UPC()) {
- lambdaTableExists = mProtonTableCollection->tableExists(lambda_A);
-
- tooSmallLambda = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A);
-
- tooSmallMean = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A);
-
+ bool skewUPC_TableExists = false;
+ bool realUPC_TableExists = false;
+ bool skewUPC_TableFits = false;
+ bool realUPC_TableFits = false;
+ bool ampUPC_TableFits = false;
+
+ skewUPC_TableExists = mProtonTableCollection->tableExists(longitudinal, lambda_skew);
+ realUPC_TableExists = mProtonTableCollection->tableExists(transverse, lambda_real);
+ if (skewUPC_TableExists) skewUPC_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_skew);
+ if (realUPC_TableExists) realUPC_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_real);
+ ampUPC_TableFits = tableFitsKinematicRange(mProtonTableCollection, mean_A);
+
+ if (mSettings->verbose() && mSettings->verboseLevel() > 1) {
+ if (!realUPC_TableExists) {
+ cout << "Info: there is no table with lambda values to perform the real" << endl;
+ cout << " amplitude corrections. The missing lambda values will be" << endl;
+ cout << " calculated on the fly from the ep amplitude tables." << endl;
+ }
+ if (!skewUPC_TableExists) {
+ cout << "Info: there is no table with lambda values to perform the skewedness" << endl;
+ cout << " corrections. We will use the lambda values from the real amplitude" << endl;
+ cout << " corrections instead. This is a good approximation." << endl;
+ }
+ if (realUPC_TableExists && !realUPC_TableFits) {
+ cout << "Info: the table with the lambda values to perform the real part corrections" << endl;
+ cout << " does not have sufficient coverage of the current kinematic range." << endl;
+ cout << " The missing lambda values will be calculated on the fly from the ep" << endl;
+ cout << " amplitude tables." << endl;
+ }
+ if (skewUPC_TableExists && !skewUPC_TableFits) {
+ cout << "Info: the table with the lambda values to perform the skewedness corrections" << endl;
+ cout << " does not have sufficient coverage of the current kinematic range." << endl;
+ cout << " The missing lambda values will be replaced by the lambda values used" << endl;
+ cout << " for the real amplitude corrections. This is a good approximation." << endl;
+
+ }
+ }
+ if (!ampUPC_TableFits && (!realUPC_TableExists || (realUPC_TableExists && !realUPC_TableFits))) {
+ cout << "Error: due to missing correction tables or tables that do not cover the" << endl;
+ cout << " full kinematic range the needed lambda values have to be calculated" << endl;
+ cout << " on the fly. However, the ep mean amplitude table that is needed for" << endl;
+ cout << " this calculation does not cover the full kinematic range selected." << endl;
+ cout << " Corrections need to be switched off to run with the existing tables" << endl;
+ cout << " and the selected kinematic range." << endl;
+ exit(1);
+ }
+
}
else {
- lambdaTableExists = mProtonTableCollection->tableExists(transverse, lambda_A);
+ bool skewT_TableExists, skewL_TableExists;
+ bool realT_TableExists, realL_TableExists;
+ bool skewT_TableFits, skewL_TableFits;
+ bool realT_TableFits, realL_TableFits;
+ bool ampT_TableFits, ampL_TableFits;
- tooSmallLambda = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A);
-
- tooSmallMean = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) ||
- !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A);
- }
-
- if (!lambdaTableExists || tooSmallLambda) {
-
- if (tooSmallMean) {
-
- cout << "Warning: the kinematic coverage of the table containing the lambda values " << endl;
- cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
- cout << " then the required range (or does not exist). The proton amplitude " << endl;
- cout << " tables that could be used to calculate them on the fly are too" << endl;
- cout << " small as well. Corrections will be switched off." << endl;
-
- mSettings->setCorrectForRealAmplitude(false);
- mSettings->setCorrectSkewedness(false);
- delete mProtonTableCollection;
- mProtonTableCollection = 0;
- }
- else {
- if (lambdaTableExists) {
- cout << "Info: the kinematic coverage of the table containing the lambda values " << endl;
- cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
- cout << " then the required range. The missing lambda values will be " << endl;
- cout << " calculated using the proton amplitude tables directly." << endl;
- }
- else {
- cout << "Info: the table containing the lambda values needed for skewedness and/or real " << endl;
- cout << " amplitude corrections is not available. The missing lambda values will be" << endl;
- cout << " calculated using the referring proton amplitude tables directly." << endl;
+ skewT_TableExists = mProtonTableCollection->tableExists(transverse, lambda_skew);
+ skewL_TableExists = mProtonTableCollection->tableExists(longitudinal, lambda_skew);
+ realT_TableExists = mProtonTableCollection->tableExists(transverse, lambda_real);
+ realL_TableExists = mProtonTableCollection->tableExists(longitudinal, lambda_real);
+
+ if (skewT_TableExists) skewT_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_skew, transverse);
+ if (skewL_TableExists) skewL_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_skew, longitudinal);
+ if (realT_TableExists) realT_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_real, transverse);
+ if (realL_TableExists) realL_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_real, longitudinal);
+ ampT_TableFits = tableFitsKinematicRange(mProtonTableCollection, mean_A, transverse);
+ ampL_TableFits = tableFitsKinematicRange(mProtonTableCollection, mean_A, longitudinal);
+
+ if (mSettings->verbose() && mSettings->verboseLevel() > 1) {
+ if (!realT_TableExists) {
+ cout << "Info: there is no table with lambda values to perform the real" << endl;
+ cout << " amplitude corrections for tranverse polarized photons. The" << endl;
+ cout << " missing lambda values will be calculated on the fly from" << endl;
+ cout << " the referring ep amplitude tables." << endl;
+ }
+ if (!realL_TableExists) {
+ cout << "Info: there is no table with lambda values to perform the real" << endl;
+ cout << " amplitude corrections for longitudinal polarized photons. " << endl;
+ cout << " The missing lambda values will be calculated on the fly " << endl;
+ cout << " from the referring ep amplitude tables." << endl;
+ }
+ if (!skewT_TableExists) {
+ cout << "Info: there is no table with lambda values to perform the skewedness" << endl;
+ cout << " corrections for tranverse polarized photons. We will use the" << endl;
+ cout << " lambda values from the real amplitude corrections instead. This" << endl;
+ cout << " is a good approximation." << endl;
+ }
+ if (!skewL_TableExists) {
+ cout << "Info: there is no table with lambda values to perform the skewedness" << endl;
+ cout << " corrections for longitudinal polarized photons. We will use the" << endl;
+ cout << " lambda values from the real amplitude corrections instead. This" << endl;
+ cout << " is a good approximation." << endl;
+ }
+ if (realT_TableExists && !realT_TableFits) {
+ cout << "Info: the table with the lambda values to perform the real part corrections" << endl;
+ cout << " for tranverse polarized photons does not have sufficient coverage of" << endl;
+ cout << " the current kinematic range. The missing lambda values will be calculated" << endl;
+ cout << " on the fly from the referring ep amplitude tables." << endl;
+ }
+ if (realL_TableExists && !realL_TableFits) {
+ cout << "Info: the table with the lambda values to perform the real part corrections" << endl;
+ cout << " for longitudinal polarized photons does not have sufficient coverage of" << endl;
+ cout << " the current kinematic range. The missing lambda values will be calculated" << endl;
+ cout << " on the fly from the referring ep amplitude tables." << endl;
+ }
+ if (skewT_TableExists && !skewT_TableFits) {
+ cout << "Info: the table with the lambda values to perform the skewedness corrections" << endl;
+ cout << " for tranverse polarized photons does not have sufficient coverage for " << endl;
+ cout << " the current kinematic range. The missing lambda values will be replaced " << endl;
+ cout << " by the lambda values used for the real amplitude corrections. This is a" << endl;
+ cout << " good approximation." << endl;
+ }
+ if (skewL_TableExists && !skewL_TableFits) {
+ cout << "Info: the table with the lambda values to perform the skewedness corrections" << endl;
+ cout << " for longitudinal polarized photons does not have sufficient coverage for " << endl;
+ cout << " the current kinematic range. The missing lambda values will be replaced " << endl;
+ cout << " by the lambda values used for the real amplitude corrections. This is a" << endl;
+ cout << " good approximation." << endl;
}
}
+ if ( (!ampT_TableFits && (!realT_TableExists || (realT_TableExists && !realT_TableFits))) ||
+ (!ampL_TableFits && (!realL_TableExists || (realL_TableExists && !realL_TableFits)))) {
+ cout << "Error: due to missing correction tables or tables that do not cover the" << endl;
+ cout << " full kinematic range the needed lambda values have to be calculated" << endl;
+ cout << " on the fly. However, the ep mean amplitude table that is needed for" << endl;
+ cout << " this calculation does not cover the full kinematic range selected." << endl;
+ cout << " Corrections need to be switched off to run with the existing tables" << endl;
+ cout << " and the selected kinematic range." << endl;
+ exit(1);
+ }
}
}
-
+
//
// Setup cross-section functor
// It is this functor that is used by all other functors,
// functions, and wrappers when dealing with cross-sections.
//
if (mCrossSection) delete mCrossSection;
mCrossSection = new CrossSection;
mCrossSection->setTableCollection(mTableCollection);
if (mProtonTableCollection) mCrossSection->setProtonTableCollection(mProtonTableCollection);
//
// UNU.RAN needs the domain (boundaries) and the mode.
// The domain is already defined, here we find the mode, which is tricky.
// The max. cross-section is clearly at the domain boundary in Q2=Q2min.
// The position in W2 and t is not obvious. It sits along a line given
// by tmax(W2). The approach here is to use the BrentMinimizer1D that
// performs first a scan a then a Brent fit.
//
double theMode[3];
if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
if (mSettings->UPC()) {
//
// For now we keep the option to create a histogram of the
// generated cross-sections(t, log(xpom)) for QA purposes.
//
if (mSettings->verbose() && mSettings->verboseLevel() > 9) {
cout << "Creating 2D histogram of cross-section as fct of t and log(xpom)." << endl;
cout << "Be patient this might take a while. Histo stored in file 'landscape.root'" << endl;
TFile file("landscape.root", "RECREATE");
int nbins_t = 300;
int nbins_logx = 300;
TH2D histo("histo", "mode finder studies",
nbins_t, mLowerLimit[0], mUpperLimit[0], // t
nbins_logx, log(mLowerLimit[1]), log(mUpperLimit[1])); // log(xpom)
for (int ix = 1; ix <= nbins_logx; ix++) {
for (int it = 1; it <= nbins_t; it++) {
double logx = histo.GetYaxis()->GetBinCenter(ix);
double x = exp(logx);
double t = histo.GetXaxis()->GetBinCenter(it);
if (!Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
t, x, vectorMesonPDG->Mass(), false)) {
histo.SetBinContent(it, ix, -5.e9);
}
else {
histo.SetBinContent(it, ix, (*mCrossSection)(t, x));
}
}
}
histo.Write();
file.Close();
cout << "Histogram written to file. All done." << endl;
}
//
// Create functor
// Mode is at the max t available.
// Note that the functop works with log(xpom).
//
theMode[0] = mUpperLimit[0]; // t
theMode[1] = mUpperLimit[1]; // xpom for UPC
theMode[2] = 0; // not used for UPC
UPCModeFinderFunctor modeFunctor(mCrossSection, vectorMesonPDG->Mass(),
mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy());
ROOT::Math::BrentMinimizer1D minimizer;
minimizer.SetFunction(modeFunctor, log(mLowerLimit[1]), log(mUpperLimit[1]));
minimizer.SetNpx(100000);
ok = minimizer.Minimize(1000000, 1.e-8, 1.e-10);
if (! ok) {
cout << "Error, failed to find mode of pdf." << endl;
return false;
}
//
// Get the result
//
theMode[1] = exp(minimizer.XMinimum()); // xpom
theMode[0] = Kinematics::tmax(theMode[1]);
double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1]);
if (mSettings->verbose()) {
cout << "\tlocation: t=" << theMode[0] << ", xpom=" << theMode[1]
<< "; value: " << crossSectionAtMode << endl;
}
//
// Consistency check of mode
//
if (crossSectionAtMode <= 0) {
cout << "Error: cross-section at mode value is invalid." << endl;
return false;
}
if (!Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
theMode[0], theMode[1],
vectorMesonPDG->Mass(), true)) {
cout << "Error: mode of pdf is outside kinematic limits." << endl;
return false;
}
}
else {
theMode[0] = mUpperLimit[0]; // t
theMode[1] = mLowerLimit[1]; // Q2 (xpom for UPC)
theMode[2] = mLowerLimit[2]; // W2 (dummy for UPC)
ModeFinderFunctor modeFunctor(mCrossSection, theMode[1], vectorMesonPDG->Mass(), mLowerLimit[0], mUpperLimit[0]);
ROOT::Math::BrentMinimizer1D minimizer;
minimizer.SetFunction(modeFunctor, mLowerLimit[2], mUpperLimit[2]);
minimizer.SetNpx(static_cast(mUpperLimit[2]-mLowerLimit[2]));
ok = minimizer.Minimize(100000, 0, 1.e-8);
if (! ok) {
cout << "Error, failed to find mode of pdf." << endl;
exit(1);
}
theMode[2] = minimizer.XMinimum(); // W2
theMode[0] = Kinematics::tmax(0, theMode[1], theMode[2], vectorMesonPDG->Mass()); // first arg (t) must be 0 here
if (theMode[0] > mUpperLimit[0]) theMode[0] = mUpperLimit[0];
double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1], theMode[2]);
if (mSettings->verbose()) {
cout << "\tlocation: t=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]);
cout << "; value: " << crossSectionAtMode << endl;
}
}
//
// Initialize 2D (UPC) or 3D random generator
//
// Test show that UNU.RAN runs smoother in log(Q2)
// and log(cross-section). Functor CrossSection has
// a spezialized method for UNU.RAN, unuranPDF().
//
// In UPC mode the mPDF_Functor is using a different
// method and is only 2D since Q2=0
//
// domain and mode for Q2 -> log(Q2) or xpom -> log(xpom)
mLowerLimit[1] = log(mLowerLimit[1]);
mUpperLimit[1] = log(mUpperLimit[1]);
theMode[1] = log(theMode[1]);
if (mPDF_Functor) delete mPDF_Functor;
if (mPDF) delete mPDF;
if (mSettings->UPC())
mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 2);
else
mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 3);
mPDF = new TUnuranMultiContDist(*mPDF_Functor, true); // last arg = pdf in log or not
mPDF->SetDomain(mLowerLimit, mUpperLimit);
mPDF->SetMode(theMode);
if (mUnuran) delete mUnuran;
mUnuran = new TUnuran;
mCrossSection->setCheckKinematics(false); // avoid numeric glitch in Init()
mUnuran->Init(*mPDF, "method=hitro");
mCrossSection->setCheckKinematics(true);
mUnuran->SetSeed(mSettings->seed());
//
// Burn in generator
//
double xrandom[3];
for (int i=0; i<1000; i++) {
mUnuran->SampleMulti(xrandom);
}
mEventCounter = 0;
mTriesCounter = 0;
mIsInitialized = true;
cout << "Sartre is initialized." << endl << endl;
return true;
}
bool Sartre::init(const string& str) // overloaded version of init()
{
if (str.empty())
return init();
else
return init(str.c_str());
}
vector > Sartre::kinematicLimits()
{
vector > array;
array.push_back(make_pair(mLowerLimit[0], mUpperLimit[0])); // t
array.push_back(make_pair(exp(mLowerLimit[1]), exp(mUpperLimit[1]))); // Q2 or xpom
if (!mSettings->UPC())
array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W
return array;
}
Event* Sartre::generateEvent()
{
if (!mIsInitialized) {
cout << "Sartre::generateEvent(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return 0;
}
double xrandom[3];
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID);
double vmMass = vectorMesonPDG->Mass();
//
// Generate one event
//
while (true) {
mTriesCounter++;
delete mCurrentEvent;
mCurrentEvent = new Event;
//
// Get t, Q2, W2 from TUnuran and check for kinematics.
// Q2 is generated as log(Q2) so we transform it back first.
// This is the only place where Kinematics::valid() is called
// with the fully correct xpomeron calculation switched on.
//
mUnuran->SampleMulti(xrandom);
xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2 or log(xpom) -> xpom
bool isValidEvent;
if (mSettings->UPC()) {
isValidEvent = Kinematics::validUPC(mSettings->hadronBeamEnergy(),
mSettings->electronBeamEnergy(),
xrandom[0], xrandom[1],
vectorMesonPDG->Mass(), (mSettings->verboseLevel() > 1));
}
else {
isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2],
vmMass, true, (mSettings->verboseLevel() > 1));
}
if (!isValidEvent) {
if (mSettings->verboseLevel() > 2)
cout << "Sartre::generateEvent(): event rejected, not within kinematic limits" << endl;
continue;
}
//
// Fill beam particles in Event structure
// Kinematics for eA is reported as 'per nucleon'
//
if (mSettings->UPC()) {
//
// For UPC some of the event variables
// are filled in the final state generator.
// They are set to 0 here.
//
mCurrentEvent->eventNumber = mEventCounter;
mCurrentEvent->t = xrandom[0]; // t
mCurrentEvent->Q2 = 0;
mCurrentEvent->x = 0;
mCurrentEvent->y = 0;
mCurrentEvent->s = mS; // s
mCurrentEvent->W = 0;
mCurrentEvent->beta = 1;
mCurrentEvent->xpom = xrandom[1];
mCurrentEvent->polarization = GammaPolarization::transverse;
mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = mUpcNucleus->pdgID(); // misuse ebeam spot for this
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
}
else {
mCurrentEvent->eventNumber = mEventCounter;
mCurrentEvent->t = xrandom[0]; // t
mCurrentEvent->Q2 = xrandom[1]; // Q2
mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
mCurrentEvent->s = mS; // s
mCurrentEvent->W = sqrt(xrandom[2]);
mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = 11; // e-
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
}
//
// Generate the final state particles
//
bool ok;
if (mSettings->UPC()) {
// also fills some of the undefined event variables
ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->xpom,
(mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent);
}
else {
ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->y, mCurrentEvent->Q2,
(mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent);
}
if (!ok) {
if (mSettings->verboseLevel() > 1) cout << "Sartre::generateEvent(): failed to generate final state" << endl;
continue;
}
break;
}
mEventCounter++;
//
// Nuclear breakup
//
// If the event is incoherent the final state generator does produce a
// 'virtual' proton with m > m_p which is used in Nucleus to calculate
// the excitation energy and the boost.
//
int indexOfScatteredHadron = 6;
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (mNucleus) mNucleus->resetBreakup(); // clear previous event in any case
if (allowBreakup && mCurrentEvent->diffractiveMode == incoherent && mNucleus) {
int nFragments = mNucleus->breakup(mCurrentEvent->particles[indexOfScatteredHadron].p);
//
// Merge the list of products into the event list.
// We loose some information here. The user can always go back to
// the nucleus and check the decay products for more details.
// In the original list the energy is per nuclei, here we transform it
// to per nucleon to stay consistent with Sartre conventions.
//
const vector& products = mNucleus->breakupProducts();
for (int i=0; iparticles.size();
fragment.pdgId = products[i].pdgId;
fragment.status = 1;
fragment.p = products[i].p*(1/static_cast(products[i].A));
fragment.parents.push_back(indexOfScatteredHadron);
mCurrentEvent->particles.push_back(fragment);
}
}
//
// Complete event record
//
if (!mSettings->UPC()) {
mCurrentEvent->xpom = Kinematics::xpomeron(mCurrentEvent->t, mCurrentEvent->Q2, mCurrentEvent->W*mCurrentEvent->W, vmMass);
mCurrentEvent->beta = mCurrentEvent->x/mCurrentEvent->xpom;
}
return mCurrentEvent;
}
double Sartre::totalCrossSection()
{
if (mTotalCrossSection == 0) {
//
// Limits of integration in t, Q2, W2
// or in the case of UPC t, xpom, dummy
//
double xmin[3];
double xmax[3];
copy(mLowerLimit, mLowerLimit+3, xmin);
copy(mUpperLimit, mUpperLimit+3, xmax);
//
// At this point mLowerLimit[1] and mUpperLimit[1]
// are in log(Q2) or log(xpom).
//
xmin[1] = exp(xmin[1]); // log Q2 limit -> Q2 limit or equivalent xpom for UPC
xmax[1] = exp(xmax[1]); // log Q2 limit -> Q2 limit
mTotalCrossSection = calculateTotalCrossSection(xmin, xmax);
}
return mTotalCrossSection;
}
double Sartre::totalCrossSection(double lower[3], double upper[3]) // t, Q2, W (or t, xpom, dummy for UPC)
{
lower[2] *= lower[2]; upper[2] *= upper[2]; // W -> W2
double result = calculateTotalCrossSection(lower, upper);
return result;
}
EventGeneratorSettings* Sartre::runSettings()
{
return EventGeneratorSettings::instance();
}
double Sartre::calculateTotalCrossSection(double lower[3], double upper[3])
{
double result = 0;
if (!mIsInitialized) {
cout << "Sartre::calculateTotalCrossSection(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return result;
}
//
// Calculate integral using adaptive numerical method
//
// Options: ADAPTIVE, kVEGAS, kPLAIN, kMISER
// no abs tolerance given -> relative only
const double precision = 5e-4;
ROOT::Math::Functor wfL((*mCrossSection), mSettings->UPC() ? 2 : 3);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE,
0, precision, 1000000);
ig.SetFunction(wfL);
result = ig.Integral(lower, upper);
//
// If it fails we switch to a MC integration which is usually more robust
// although not as accurate. This should happen very rarely if at all.
//
if (result <= numeric_limits::epsilon()) {
cout << "Sartre::calculateTotalCrossSection(): warning, adaptive integration failed - switching to VEGAS method." << endl;
ROOT::Math::IntegratorMultiDim igAlt(ROOT::Math::IntegrationMultiDim::kVEGAS);
igAlt.SetFunction(wfL);
igAlt.SetRelTolerance(precision);
igAlt.SetAbsTolerance(0);
result = igAlt.Integral(lower, upper);
}
//
// UPC with symmetric beams requires the cross-section
// being multiplied by a factor of 2.
//
if (mSettings->UPC()) {
if (mSettings->A() == mSettings->UPCA()) result *= 2;
}
return result;
}
const FrangibleNucleus* Sartre::nucleus() const {return mNucleus;}
void Sartre::listStatus(ostream& os) const
{
os << "Event summary: " << mEventCounter<< " events generated, " << mTriesCounter << " tried" << endl;
time_t delta = runTime();
os << "Total time used: " << delta/60 << " min " << delta - 60*(delta/60) << " sec" << endl;
}
time_t Sartre::runTime() const
{
time_t now = time(0);
return now-mStartTime;
}
+
+bool Sartre::tableFitsKinematicRange(TableCollection* tc, AmplitudeMoment mom, GammaPolarization pol)
+{
+ if (!tc) return false;
+
+ bool fitsIn = tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom);
+
+ return fitsIn;
+ }
+
+// UPC version
+bool Sartre::tableFitsKinematicRange(TableCollection* tc, AmplitudeMoment mom)
+{
+ if (!tc) return false;
+
+ bool fitsIn = tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) &&
+ tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom) &&
+ tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) &&
+ tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) &&
+ tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom);
+
+ return fitsIn;
+}
+
//==============================================================================
//
// Utility functions and operators (helpers)
//
//==============================================================================
ostream& operator<<(ostream& os, const TLorentzVector& v)
{
os << v.Px() << '\t' << v.Py() << '\t' << v.Pz() << '\t' << v.E() << '\t';
double m2 = v*v;
if (m2 < 0)
os << '(' << -sqrt(-m2) << ')';
else
os << '(' << sqrt(m2) << ')';
return os;
}