Page MenuHomeHEPForge
Files F8309030

CrossSection.cpp
No OneTemporary
Actions

CrossSection.cpp
View Options

//==============================================================================
// CrossSection.cpp
//
// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre version: 1.1
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date: 2013-05-29 21:26:19 +0100 (Wed, 29 May 2013) $
// $Author: thomas.ullrich@bnl.gov $
//==============================================================================
#include "CrossSection.h"
#include "EventGeneratorSettings.h"
#include "Kinematics.h"
#include "TableCollection.h"
#include "Table.h"
#include "Constants.h"
#include <cmath>
#include <limits>
#include "TH2F.h"
#include "TFile.h"
#include "Math/Interpolator.h"
#include "Math/InterpolationTypes.h"
//#define DO_NOT_APPLY_PHOTON_FLUX 1
#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(t, Q2, W2); }
double CrossSection::operator()(const double* array) { return dsigdtdQ2dW2_total(array[0], array[1], array[2]); }
double CrossSection::unuranPDF(const double* array)
{
double Q2 = exp(array[1]);
double result = dsigdtdQ2dW2_total(array[0], Q2, array[2]);
result *= Q2; // Jacobian
return log(result);
}
double CrossSection::dsigdtdQ2dW2_total(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(): 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() > 9) {
cout << "CrossSection::dsigdtdQ2dW2_total(): " << 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 (isnan(result)) {
cout << "CrossSection::dsigdtdQ2dW2_total(): error, return value = NaN at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
if (isinf(result)) {
cout << "CrossSection::dsigdtdQ2dW2_total(): error, return value = inf at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
if (result < 0) {
cout << "CrossSection::dsigdtdQ2dW2_total(): error, negative cross-section at" << endl;
cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl;
result = 0;
}
return result;
}
double CrossSection::dsigdt_total(double t, double Q2, double W2, GammaPolarization pol) const
{
double 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() || mSettings->correctSkewedness()) // calculate lambda only once
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
if (mSettings->correctForRealAmplitude())
result *= realAmplitudeCorrection(lambda);
if (mSettings->correctSkewedness())
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() || mSettings->correctSkewedness()) // calculate lambda only once
lambda = logDerivateOfAmplitude(t, Q2, W2, pol);
if (mSettings->correctForRealAmplitude())
result *= realAmplitudeCorrection(lambda);
if (mSettings->correctSkewedness())
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;
}
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;
}
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 <A> table(s).
//
lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_A, 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);
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 the table is a delicate issue.
// The standard interpolation method(s) used in Table::get() are
// not sufficiently accurate and cause enormous ripples on lambda(W2).
// Here we use a more sophisticated interpolation along W2 (1dim only).
//
//
// Fill vectors with all available central points (bin center)
// an set up Interpolator. Akima method requires at least 5 points.
//
vector<double> W2_values;
vector<double> table_values;
double dW2 = table->binWidthW2();
double W2val = table->minW2();
while (W2val <= table->maxW2()) {
if (fabs(W2val-W2) < 6*dW2) {
W2_values.push_back(W2val);
table_values.push_back(mProtonTableCollection->get(Q2, W2val, t, pol, mean_A));
}
W2val += dW2;
}
//
// Setup the interpolation engine.
// Available alogorithms are kLINEAR, kPOLYNOMIAL, kCSPLINE, kAKIMA.
// In extensive tests on tables for ep and eA it turns out akima is the
// best choice.
//
ROOT::Math::Interpolator ip(W2_values, table_values, ROOT::Math::Interpolation::kAKIMA);
//
// Calculate the derivative using simple method.
// Turns out to be better than Interpolator::Deriv().
//
double derivate;
double hplus, hminus;
hplus = hminus = 0.5*dW2; // half bin width is found to be the best choice after some testing
hplus = min(hplus, mProtonTableCollection->maxW2()-W2);
hminus = min(hminus, W2-mProtonTableCollection->minW2());
hminus -= numeric_limits<float>::epsilon();
hplus -= numeric_limits<float>::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 = ip.Eval(W2+hplus);
double b = ip.Eval(W2-hminus);
derivate=(a-b)/(hplus+hminus);
//
// Finally calculate lambda
//
double jacobian = (W2-protonMass2+Q2)/value;
lambda = jacobian*derivate;
} // 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;
}
//
// 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;
if (lambda < 0) {
if (mSettings->verboseLevel() > 2) {
cout << "CrossSection::logDerivateOfAmplitude(): ";
cout << "Warning, lambda is negative (" << lambda << "). " ;
cout << "Set to " << (-lambda) << "." << endl;
}
lambda = -lambda;
}
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 (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 (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::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;
}
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;
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 2:00 PM (15 h, 41 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4015104
Default Alt Text
CrossSection.cpp (15 KB)

Event Timeline