Index: trunk/src/CrossSection.h
===================================================================
--- trunk/src/CrossSection.h (revision 436)
+++ trunk/src/CrossSection.h (revision 437)
@@ -1,104 +1,106 @@
//==============================================================================
// CrossSection.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$
//==============================================================================
//
// Functor class.
//
// operator() returns d^3sig/(dt dQ2 dW2) in nb/GeV^6.
//
//===============================================================================
#ifndef CrossSection_h
#define CrossSection_h
#include "Enumerations.h"
#include "PhotonFlux.h"
class TRandom3;
class EventGeneratorSettings;
class TableCollection;
class CrossSection {
public:
CrossSection(TableCollection* = 0, TableCollection* = 0);
~CrossSection();
double operator()(double t, double Q2, double W2);
double operator()(double t, double xpom); // UPC version
double operator()(const double*); // array of t, Q2, W2 or t, xpom for UPC
double unuranPDF(const double*); // for UNU.RAN using log(Q2) (or log(xpom) for UPC)
// and returning log of cross-section
void setTableCollection(TableCollection*);
void setProtonTableCollection(TableCollection*);
GammaPolarization polarizationOfLastCall() const;
DiffractiveMode diffractiveModeOfLastCall() const;
-
+ double crossSectionRatioLTOfLastCall() const;
+
void setCheckKinematics(bool);
double dsigdtdQ2dW2_total(double t, double Q2, double W2, GammaPolarization) const;
protected:
friend class Sartre; // mostly for debugging and QA
double dsigdtdQ2dW2_total_checked(double t, double Q2, double W2);
double dsigdtdxp_total_checked(double t, double xpom);
double dsigdt_total(double t, double Q2, double W2, GammaPolarization) const; // modified
double dsigdt_coherent(double t, double Q2, double W2, GammaPolarization) const;
double dsigdt_incoherent(double t, double Q2, double W2, GammaPolarization) const; // new
double dsigdtdQ2dW2_coherent(double t, double Q2, double W2, GammaPolarization) const;
double dsigdt_total(double t, double xpom) const; // UPC version
double dsigdt_coherent(double t, double xpom) const; // UPC version
double dsigdt_incoherent(double t, double xpom) const; // UPC version
double dsigdtdxp_total(double t, double xpom) const; // UPC
double dsigdtdxp_coherent(double t, double xpom) const; // UPC
double logDerivateOfAmplitude(double t, double Q2, double W2, GammaPolarization) const;
double logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization) const;
double realAmplitudeCorrection(double t, double Q2, double W2, GammaPolarization pol) const;
double skewednessCorrection(double t, double Q2, double W2, GammaPolarization pol) const;
double logDerivateOfAmplitude(double t, double xpom) const; // UPC version
double logDerivateOfGluonDensity(double t, double xpom) const; // UPC version
double realAmplitudeCorrection(double t, double xpom) const; // UPC version
double skewednessCorrection(double t, double xpom) const; // UPC version
double skewednessCorrection(double lambda) const;
double realAmplitudeCorrection(double lambda) const;
double UPCPhotonFlux(double t, double xpom) const;
private:
TRandom3* mRandom;
GammaPolarization mPolarization;
DiffractiveMode mDiffractiveMode;
PhotonFlux mPhotonFlux;
TableCollection* mTableCollection;
TableCollection* mProtonTableCollection;
EventGeneratorSettings* mSettings;
double mS;
- double mVmMass;
-
+ double mVmMass;
+ double mCrossSectionRatioLT;
+
bool mCheckKinematics;
};
#endif
Index: trunk/src/CrossSection.cpp
===================================================================
--- trunk/src/CrossSection.cpp (revision 436)
+++ trunk/src/CrossSection.cpp (revision 437)
@@ -1,921 +1,926 @@
//==============================================================================
// CrossSection.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "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;
+ mCheckKinematics = true;
+ mCrossSectionRatioLT = 0;
}
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;
+ result = csT + csL;
+ mCrossSectionRatioLT = csL/csT;
//
// 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);
result *= realAmplitudeCorrection(lambda);
}
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 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_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);
}
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);
}
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;}
+DiffractiveMode CrossSection::diffractiveModeOfLastCall() const {return mDiffractiveMode;}
+
+double CrossSection::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
+
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_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 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.
//
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 (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;
cout << " Set to 0." << endl;
return 0;
}
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;
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());
}
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_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
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);
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;
else
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.
//
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;
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;
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);
}
if (mSettings->verboseLevel() > 3) {
cout << "CrossSection::logDerivateOfGluonDensity(): ";
if (lambdaFromTable)
cout << "Info, lambda taken from table." << endl;
else
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.
//
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;
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;
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 436)
+++ trunk/src/Sartre.cpp (revision 437)
@@ -1,1106 +1,1108 @@
//==============================================================================
// 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;
exit(1);
}
if (mTableCollection) delete mTableCollection;
mTableCollection = new TableCollection;
ok = mTableCollection->init(mA, mDipoleModelType, mDipoleModelParameterSet, mVmID);
if (!ok) {
cout << "Error, could not initialize lookup tables for requested process." << endl;
return false;
}
//
// Load in the p tables for the lambda lookup tables (or to calculate lambda if not available)
//
if (mSettings->correctForRealAmplitude() || mSettings->correctSkewedness()) {
if (mA == 1) {
mProtonTableCollection = mTableCollection;
}
else {
if (mProtonTableCollection) delete mProtonTableCollection;
mProtonTableCollection = new TableCollection;
ok = mProtonTableCollection->init(1, mDipoleModelType, mDipoleModelParameterSet, mVmID);
if (!ok) {
cout << "Error: could not initialize proton lookup tables for requested process." << endl;
cout << " These tables are needed for skewedness and real amplitude corrections." << endl;
return false;
}
}
}
else
mProtonTableCollection = 0;
//
// Kinematic limits and generator range
//
// There are 3 ranges we have to deal with
// 1. the kinematic range requested by the user
// if given.
// The user can only control Q2 and W but not t.
// For UPC that's xpom.
// 2. the range of the table(s)
// 3. the kinematically allowed range
//
// Of course (3) is more complex than a simple cube/square.
// However, we deal with the detailed shape of the kinematic
// range later using Kinematics::valid() when we generate the
// individual events.
// For setting up UNU.RAN we have to get the cubic/square
// envelope that satifies (1)-(3).
// Note, that they are correlated which makes the order
// in which we do things a bit tricky.
//
//
// Step 1:
// Set the limits to that of the table(s).
// Note, the indices 0-2 refer to t, Q2, and W2.
// For UPC we have only t and xpom.
//
if (mSettings->UPC()) {
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minX(); mUpperLimit[1] = mTableCollection->maxX();
mLowerLimit[2] = mUpperLimit[2] = 0;
}
else {
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minQ2(); mUpperLimit[1] = mTableCollection->maxQ2();
mLowerLimit[2] = mTableCollection->minW2(); mUpperLimit[2] = mTableCollection->maxW2();
}
//
// Step 2:
// Kinematic limits might overrule boundaries from step 1
//
if (mSettings->UPC()) {
double kineXpomMin = Kinematics::xpomMin(vectorMesonPDG->Mass(), mUpperLimit[0], mHadronBeam, mElectronBeam);
double kineXpomMax = 1;
mLowerLimit[1] = max(kineXpomMin, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineXpomMax);
double kineTmax = Kinematics::tmax(kineXpomMin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
}
else {
// double kineYmax = Kinematics::ymax(mS, vectorMesonPDG->Mass());
double kineYmin = Kinematics::ymin(mS, vectorMesonPDG->Mass());
double kineQ2min = Kinematics::Q2min(kineYmin);
double kineQ2max = Kinematics::Q2max(mS);
double kineW2min = Kinematics::W2min(vectorMesonPDG->Mass());
double kineW2max = Kinematics::W2max(mS);
kineQ2min = max(kineQ2min, mLowerLimit[1]);
kineQ2max = min(kineQ2max, mUpperLimit[1]);
kineW2min = max(kineW2min, mLowerLimit[2]);
kineW2max = min(kineW2max, mUpperLimit[2]);
double kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // first arg (t) set to 0 (recursive)
double kineTmax = Kinematics::tmax(kineXPmin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
mLowerLimit[1] = max(kineQ2min, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineQ2max);
mLowerLimit[2] = max(kineW2min, mLowerLimit[2]); mUpperLimit[2] = min(mUpperLimit[2], kineW2max);
}
//
// Step 3:
// Deal with user provided limits.
// User settings are ignored (switched off) if min >= max.
//
if (mSettings->UPC()) {
if (mSettings->xpomMin() < mSettings->xpomMax()) {
if (mSettings->xpomMin() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of xpomeron (" << mSettings->xpomMin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->xpomMin();
}
if (mSettings->xpomMax() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of xpomeron (" << mSettings->xpomMax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->xpomMax();
}
}
}
else {
if (mSettings->W2min() < mSettings->W2max()) { // W2 first
if (mSettings->W2min() < mLowerLimit[2]) {
cout << "Warning, requested lower limit of W (" << mSettings->Wmin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << sqrt(mLowerLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[2] = mSettings->W2min();
}
if (mSettings->W2max() > mUpperLimit[2]) {
cout << "Warning, requested upper limit of W (" << mSettings->Wmax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << sqrt(mUpperLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[2] = mSettings->W2max();
}
}
if (mSettings->Q2min() < mSettings->Q2max()) { // Q2
if (mSettings->Q2min() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of Q2 (" << mSettings->Q2min() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->Q2min();
// ????? kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // new Q2min changes tmax
// ????? mUpperLimit[0] = min(mUpperLimit[0], Kinematics::tmax(kineXPmin));
}
if (mSettings->Q2max() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of Q2 (" << mSettings->Q2max() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->Q2max();
}
}
}
//
// Check if any phase space is left
//
if (mLowerLimit[0] >= mUpperLimit[0]) {
cout << "Invalid range in t: t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
exit(1);
}
if (mLowerLimit[1] >= mUpperLimit[1]) {
if (mSettings->UPC())
cout << "Invalid range in xpomeron: xpomeron=[";
else
cout << "Invalid range in Q2: Q2=[";
cout << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl;
exit(1);
}
if (!mSettings->UPC() && mLowerLimit[2] >= mUpperLimit[2]) {
cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl;
exit(1);
}
//
// Print-out limits (all verbose levels)
//
if (mSettings->verbose()) {
cout << "Kinematic limits used for event generation:" << endl;
if (mSettings->UPC()) {
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "xpom=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
}
else {
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl;
}
}
//
// 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) {
if (mSettings->UPC()) {
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 {
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;
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;
}
}
// for EIC
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();
+ mCurrentEvent->crossSectionRatioLT = mCrossSection->crossSectionRatioLTOfLastCall();
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();
+ mCurrentEvent->crossSectionRatioLT = mCrossSection->crossSectionRatioLTOfLastCall();
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;
}