Index: trunk/src/Amplitudes.cpp
===================================================================
--- trunk/src/Amplitudes.cpp (revision 432)
+++ trunk/src/Amplitudes.cpp (revision 433)
@@ -1,334 +1,334 @@
//==============================================================================
// Amplitudes.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: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
//#define SARTRE_IN_MULTITHREADED_MODE 1
#include
#include
#include "Amplitudes.h"
#include "Constants.h"
#include "TableGeneratorSettings.h"
#include "DglapEvolution.h"
#include "Enumerations.h"
#include "Kinematics.h"
#include "Integrals.h"
#include "DipoleModel.h"
#if defined(SARTRE_IN_MULTITHREADED_MODE)
#include
#endif
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
Amplitudes::Amplitudes()
{
mAmplitudeT = 0;
mAmplitudeL = 0;
mAmplitudeT2 = 0;
mAmplitudeL2 = 0;
mNumberOfConfigurations = 0;
mTheModes = 0;
mA = 0;
mErrorT = 0;
mErrorL = 0;
mErrorT2 = 0;
mErrorL2 = 0;
mAmplitudeTForSkewednessCorrection = 0;
mAmplitudeLForSkewednessCorrection = 0;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mNumberOfConfigurations = settings->numberOfConfigurations();
mVerbose = settings->verbose();
//
// Create a vector containing instances of the Integrals class
// and initialize them:
//
for (int i=0; i<=mNumberOfConfigurations; i++) {
mIntegrals.push_back(new IntegralsExclusive);
}
- mA=settings->A();
- mUPC=settings->UPC();
+ mA = settings->A();
+ mUPC = settings->UPC();
isBNonSat = false;
- if(settings->dipoleModelType() == bNonSat)
+ if (settings->dipoleModelType() == bNonSat)
isBNonSat = true;
//
// Get the modes to calculate:
// 0: analytically averaged over configurations
// 1: only analytically
// 2: Both and averaged over configurations
//
- mTheModes=settings->modesToCalculate();
+ mTheModes = settings->modesToCalculate();
}
Amplitudes& Amplitudes::operator=(const Amplitudes& amp)
{
if (this != &) {
for (unsigned int i=0; idipoleModel()->createConfiguration(i);
}
//void Amplitudes::calculate(double t, double Q2, double W2)
void Amplitudes::calculate(double* kinematicPoint)
{
double t=0, Q2=0, W2=0, xpom=0;
- if(!mUPC){
- t =kinematicPoint[0];
- Q2=kinematicPoint[1];
- W2=kinematicPoint[2];
+ if (!mUPC){
+ t = kinematicPoint[0];
+ Q2 = kinematicPoint[1];
+ W2 = kinematicPoint[2];
}
else{
- t =kinematicPoint[0];
- xpom=kinematicPoint[1];
+ t = kinematicPoint[0];
+ xpom = kinematicPoint[1];
}
#if defined(SARTRE_IN_MULTITHREADED_MODE) // multithreaded version
if (mA == 1) {
cout << "Amplitudes::calculate(): Multithreaded mode (SARTRE_IN_MULTITHREADED_MODE)" << endl;
cout << " is not supported for ep (A=1). Stopping." << endl;
exit(1);
}
//
// Create a vector containing the threads:
//
std::vector vThreads;
vThreads.clear();
//
// Create the thread group:
//
boost::thread_group gThreads;
if (mTheModes==0 || mTheModes == 2){
//Start loop over configurations, each calculated on a separate thread:
for (int i=0; i1 && (mTheModes==1 || mTheModes == 0)) {
- if(!mUPC)
+ if (!mUPC)
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2);
else
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom);
}
if (mTheModes==0 || mTheModes == 2) {
// Wait for all threads to finish before continuing main thread:
gThreads.join_all();
// Clean up threads
vThreads.clear();
}
#else // unforked version
if ((mTheModes==0 || mTheModes == 2) || mA==1) {
//Start loop over configurations:
for (int i=0; ioperator()(t, Q2, W2);
else
mIntegrals.at(i)->operator()(t, xpom);
}
}
//
// Calculate coherent cross-section according to eq.(47) in KT arXiv:hep-ph/0304189v3,
// (only in eA)
//
if (mA>1 && (mTheModes==1 || mTheModes == 0)){
- if(!mUPC)
+ if (!mUPC)
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2);
else
mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom);
}
#endif
//
// Calculate the resulting :
//
- double coherentT=0, coherentL=0;
- double errCoherentT=0, errCoherentL=0;
+ double coherentT = 0, coherentL = 0;
+ double errCoherentT = 0, errCoherentL = 0;
if ((mTheModes==0 || mTheModes == 2) || mA==1) {
double totalT = 0;
double totalL = 0;
- double err2TotalT=0, err2TotalL=0;
+ double err2TotalT = 0, err2TotalL = 0;
double probabilityCutOff = 1e-6;
for (int i=0; iintegralImT();
- double valreT=mIntegrals.at(i)->integralReT();
- double valimL=mIntegrals.at(i)->integralImL();
- double valreL=mIntegrals.at(i)->integralReL();
+ double valimT = mIntegrals.at(i)->integralImT();
+ double valreT = mIntegrals.at(i)->integralReT();
+ double valimL = mIntegrals.at(i)->integralImL();
+ double valreL = mIntegrals.at(i)->integralReL();
- double errimT=mIntegrals.at(i)->errorImT();
- double errimL=mIntegrals.at(i)->errorImL();
- double errreT=mIntegrals.at(i)->errorReT();
- double errreL=mIntegrals.at(i)->errorReL();
+ double errimT = mIntegrals.at(i)->errorImT();
+ double errimL = mIntegrals.at(i)->errorImL();
+ double errreT = mIntegrals.at(i)->errorReT();
+ double errreL = mIntegrals.at(i)->errorReL();
- double probimT=mIntegrals.at(i)->probImT();
- double probimL=mIntegrals.at(i)->probImL();
- double probreT=mIntegrals.at(i)->probReT();
- double probreL=mIntegrals.at(i)->probReL();
+ double probimT = mIntegrals.at(i)->probImT();
+ double probimL = mIntegrals.at(i)->probImL();
+ double probreT = mIntegrals.at(i)->probReT();
+ double probreL = mIntegrals.at(i)->probReL();
if (probimT > probabilityCutOff || probimL > probabilityCutOff ||
probreT > probabilityCutOff || probreL > probabilityCutOff){
- if(mVerbose) {
+ if (mVerbose) {
cout<< "Amplitudes::calculate(): Warning, Integrals may not have reached desired precision" < probreT and probimT > probimL and probimT > probreL ?
+ probimT > probreT && probimT > probimL && probimT > probreL ?
cout<< " The probability for this is "< probimT and probreT > probimL and probreT > probreL ?
+ probreT > probimT && probreT > probimL && probreT > probreL ?
cout<< " The probability for this is "< probimT and probimL > probreT and probimL > probreL ?
+ probimL > probimT && probimL > probreT && probimL > probreL ?
cout<< " The probability for this is "<1 && (mTheModes==1 || mTheModes == 0)) {
- double coherentKTT=mIntegrals.at(mNumberOfConfigurations)->integralImT();
- double coherentKTL=mIntegrals.at(mNumberOfConfigurations)->integralImL();
- double errCoherentKTT=mIntegrals.at(mNumberOfConfigurations)->errorImT();
- double errCoherentKTL=mIntegrals.at(mNumberOfConfigurations)->errorImL();
- mAmplitudeT=coherentKTT;
- mAmplitudeL=coherentKTL;
- mErrorT=errCoherentKTT;
- mErrorL=errCoherentKTL;
+ double coherentKTT = mIntegrals.at(mNumberOfConfigurations)->integralImT();
+ double coherentKTL = mIntegrals.at(mNumberOfConfigurations)->integralImL();
+ double errCoherentKTT = mIntegrals.at(mNumberOfConfigurations)->errorImT();
+ double errCoherentKTL = mIntegrals.at(mNumberOfConfigurations)->errorImL();
+ mAmplitudeT = coherentKTT;
+ mAmplitudeL = coherentKTL;
+ mErrorT = errCoherentKTT;
+ mErrorL = errCoherentKTL;
}
else {
- mAmplitudeT=coherentT/mNumberOfConfigurations;
- mAmplitudeL=coherentL/mNumberOfConfigurations;
- mErrorT=errCoherentT/mNumberOfConfigurations;
- mErrorL=errCoherentL/mNumberOfConfigurations;
+ mAmplitudeT = coherentT/mNumberOfConfigurations;
+ mAmplitudeL = coherentL/mNumberOfConfigurations;
+ mErrorT = errCoherentT/mNumberOfConfigurations;
+ mErrorL = errCoherentL/mNumberOfConfigurations;
}
- if (mA==1 and mTheModes != 1 and mNumberOfConfigurations == 1){
+ if (mA == 1 && mTheModes != 1 && mNumberOfConfigurations == 1){
if (isBNonSat){
mAmplitudeTForSkewednessCorrection = mAmplitudeT;
mAmplitudeLForSkewednessCorrection = mAmplitudeL;
}
else {
mAmplitudeTForSkewednessCorrection = mIntegrals.at(0)->integralTForSkewedness();
mAmplitudeLForSkewednessCorrection = mIntegrals.at(0)->integralLForSkewedness();
}
}
}
Index: trunk/src/EventGeneratorSettings.cpp
===================================================================
--- trunk/src/EventGeneratorSettings.cpp (revision 432)
+++ trunk/src/EventGeneratorSettings.cpp (revision 433)
@@ -1,116 +1,116 @@
//==============================================================================
// EventGeneratorSettings.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 "EventGeneratorSettings.h"
#include "Constants.h"
#include
#include
EventGeneratorSettings* EventGeneratorSettings::mInstance = 0; // initialize static
EventGeneratorSettings* EventGeneratorSettings::instance()
{
if (mInstance == 0)
mInstance = new EventGeneratorSettings;
return mInstance;
}
EventGeneratorSettings::EventGeneratorSettings()
{
//
// Register all parameters.
// list() will print them in the order they are defined here.
//
registerParameter(&mElectronBeamEnergy, "eBeamEnergy", 10.);
registerParameter(&mHadronBeamEnergy, "hBeamEnergy", 100.);
registerParameter(&mNumberOfEvents, "numberOfEvents", static_cast(10000));
registerParameter(&mTimesToShow, "timesToShow", static_cast(100));
registerParameter(&mCorrectForRealAmplitude, "correctForRealAmplitude", true);
registerParameter(&mCorrectSkewedness, "correctSkewedness", true);
registerParameter(&mEnableNuclearBreakup, "enableNuclearBreakup", false);
registerParameter(&mMaxLambdaUsedInCorrections, "maxLambdaUsedInCorrections", 0.65);
registerParameter(&mMaxNuclearExcitationEnergy, "maxNuclearExcitationEnergy", 1.4);
registerParameter(&mApplyPhotonFlux, "applyPhotonFlux", true);
}
void EventGeneratorSettings::consolidateSettings() // called after runcard is read
{
//
// Disable nuclear breakup in UPC mode for now (July 2018)
//
if (mUPC && mEnableNuclearBreakup) {
cout << "EventGeneratorSettings::consolidateSettings(): Nuclear breakup is currently not possible in UPC mode.\n"
<< " Switched off now." << endl;
mEnableNuclearBreakup = false;
}
}
//
// Access functions
//
void EventGeneratorSettings::setNumberOfEvents(unsigned long val) { mNumberOfEvents = val;}
unsigned long EventGeneratorSettings::numberOfEvents() const {return mNumberOfEvents;}
void EventGeneratorSettings::setTimesToShow(unsigned int val) { mTimesToShow = val;}
unsigned int EventGeneratorSettings::timesToShow() const {return mTimesToShow;}
double EventGeneratorSettings::electronBeamEnergy() const {return mElectronBeamEnergy;}
void EventGeneratorSettings::setElectronBeamEnergy(double val) {mElectronBeamEnergy = val;}
double EventGeneratorSettings::hadronBeamEnergy() const {return mHadronBeamEnergy;}
void EventGeneratorSettings::setHadronBeamEnergy(double val) {mHadronBeamEnergy = val;}
TLorentzVector EventGeneratorSettings::eBeam() const
{
double mass2;
- if(mUPC)
+ if (mUPC)
mass2 = protonMass2; //#TT should be mass/nucleon in UPC...?
else
mass2 = electronMass2;
return TLorentzVector(0, 0, -sqrt(mElectronBeamEnergy*mElectronBeamEnergy-mass2), mElectronBeamEnergy);
}
TLorentzVector EventGeneratorSettings::hBeam() const
{
return TLorentzVector(0, 0, sqrt(mHadronBeamEnergy*mHadronBeamEnergy-protonMass2), mHadronBeamEnergy);
}
bool EventGeneratorSettings::correctForRealAmplitude() const {return mCorrectForRealAmplitude;}
void EventGeneratorSettings::setCorrectForRealAmplitude(bool val) {mCorrectForRealAmplitude = val;}
bool EventGeneratorSettings::correctSkewedness() const {return mCorrectSkewedness;}
void EventGeneratorSettings::setCorrectSkewedness(bool val) {mCorrectSkewedness = val;}
bool EventGeneratorSettings::enableNuclearBreakup() const {return mEnableNuclearBreakup;}
void EventGeneratorSettings::setEnableNuclearBreakup(bool val) {mEnableNuclearBreakup = val;}
double EventGeneratorSettings::maxNuclearExcitationEnergy() const {return mMaxNuclearExcitationEnergy;}
void EventGeneratorSettings::setMaxNuclearExcitationEnergy(double val) {mMaxNuclearExcitationEnergy = val;}
double EventGeneratorSettings::maxLambdaUsedInCorrections() const {return mMaxLambdaUsedInCorrections;}
void EventGeneratorSettings::setMaxLambdaUsedInCorrections(double val) {mMaxLambdaUsedInCorrections = val;}
bool EventGeneratorSettings::applyPhotonFlux() const {return mApplyPhotonFlux;}
void EventGeneratorSettings::setApplyPhotonFlux(bool val) {mApplyPhotonFlux = val;}
Index: trunk/src/PhotonFlux.h
===================================================================
--- trunk/src/PhotonFlux.h (revision 432)
+++ trunk/src/PhotonFlux.h (revision 433)
@@ -1,72 +1,72 @@
//==============================================================================
// PhotonFlux.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.
// Photon flux is given in: d2sig/(dQ2 dW2)
//==============================================================================
#ifndef PhotonFlux_h
#define PhotonFlux_h
#include "Enumerations.h"
#include "EventGeneratorSettings.h"
#include "TH1D.h"
class Nucleus;
class PhotonFlux {
public:
PhotonFlux();
PhotonFlux(double s);
void setS(double);
double operator()(double Q2, double W2, GammaPolarization p) const;
- //UPC:
+
+ // UPC only:
double operator()(double Egamma) const;
double nuclearPhotonFlux(double Egamma) const;
-
-private:
- double fluxTransverse(double Q2, double W2) const;
- double fluxLongitudinal(double Q2, double W2) const;
-
- double mS;
- bool mSIsSet;
-
- //UPC:
- double uiNuclearPhotonFlux(double*, double*) const;
- double sigma_nn(double) const;
+
+protected:
+ double fluxTransverse(double Q2, double W2) const;
+ double fluxLongitudinal(double Q2, double W2) const;
double TAA(double);
double TAAForIntegration(const double*) const;
- void calculateTAAlookupTable();
+ void setupTAAlookupTable();
+
+ // UPC only:
+ double unintegratedNuclearPhotonFlux(double*, double*) const;
+ double sigma_nn(double) const;
+private:
+ double mS;
+ bool mSIsSet;
bool mIsUPC;
double mEBeamEnergy;
- TH1D* mTAA_of_b;
double mB;
+ TH1D* mTAA_of_b;
Nucleus* mNucleus;
Nucleus* mNucleusUPC;
EventGeneratorSettings *mSettings;
-
};
#endif
Index: trunk/src/Nucleus.cpp
===================================================================
--- trunk/src/Nucleus.cpp (revision 432)
+++ trunk/src/Nucleus.cpp (revision 433)
@@ -1,343 +1,343 @@
//==============================================================================
// Nucleus.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 "Nucleus.h"
#include "Constants.h"
#include "TF1.h"
#include "TH1.h"
#include "TH1D.h"
#include
#include
#define PR(x) cout << #x << " = " << (x) << endl;
Nucleus::Nucleus()
{
mA = mZ = 0;
mRadius = 0;
mSpin = 0;
mSurfaceThickness = 0;
mRho0 = 0;
mOmega = 0;
mMass = 0;
mHulthenA = 0;
mHulthenB = 0;
}
Nucleus::Nucleus(unsigned int A)
{
mA = mZ = 0;
mRadius = 0;
mSpin = 0;
mSurfaceThickness = 0;
mRho0 = 0;
mOmega = 0;
mMass = 0;
mHulthenA = 0;
mHulthenB = 0;
Nucleus::init(A);
}
Nucleus& Nucleus::operator=(const Nucleus& n)
{
if (this != &n) {
mA = n.mA;
mZ = n.mZ;
mRadius = n.mRadius;
mSpin = n.mSpin;
mSurfaceThickness = n.mSurfaceThickness;
mRho0 = n.mRho0;
mOmega = n.mOmega;
mMass = n.mMass;
mName = n.mName;
mHulthenA = n.mHulthenA;
mHulthenB = n.mHulthenB;
mLookupTable = unique_ptr(new TH1D(*(n.mLookupTable)));
mLookupTable->SetDirectory(0);
}
return *this;
}
Nucleus::Nucleus(const Nucleus& n)
{
mA = n.mA;
mZ = n.mZ;
mRadius = n.mRadius;
mSpin = n.mSpin;
mSurfaceThickness = n.mSurfaceThickness;
mRho0 = n.mRho0;
mOmega = n.mOmega;
mMass = n.mMass;
mName = n.mName;
mHulthenA = n.mHulthenA;
mHulthenB = n.mHulthenB;
mLookupTable = unique_ptr(new TH1D(*(n.mLookupTable)));
mLookupTable->SetDirectory(0);
}
Nucleus::~Nucleus()
{
//
// This is a temorary fix of a problem in the code that is not understood.
// In the table generation code after all is calculated the ~TH1 of the
// look up table histogram causes a crash. This appears (?) to be related
// to issues in GSL and ROOT. Releasing the pointer causes a memory leak
// but this is not a big problem since it happens at the end of the program.
//
mLookupTable.release();
}
void Nucleus::init(unsigned int A)
{
mA = A;
//
// Set the parameters of the nucleus
// Woods-Saxon parameter are from Ramona Vogt's paper on nuclear geometry (table 1)
- // One exception is the deuteron were a Hulthen function is used (nucl-ex/0701025v1).
+ // One exception is the deuteron where a Hulthen function is used (nucl-ex/0701025v1).
// All version are defined such that Integral(d3r rho(r)) = A
//
const double unifiedAtomicMass = 0.931494013;
switch (mA) {
case 208: // Pb
mZ = 82;
mRadius = 6.624;
mSurfaceThickness = 0.549;
mOmega = 0;
mRho0 = 0.160;
mName = "Pb";
mMass = 207.2*unifiedAtomicMass;
mSpin = 1./2.;
break;
case 197: // Au
mZ = 79;
mRadius = 6.38;
mSurfaceThickness = 0.535;
mOmega = 0;
mRho0 = 0.1693;
mName = "Au";
mMass = 196.96655*unifiedAtomicMass;
mSpin = 3./2.;
break;
case 110: // Cd
mZ = 48;
mRadius = 5.33;
mSurfaceThickness = 0.535;
mOmega = 0;
mRho0 = 0.1577;
mName = "Cd";
mMass = 112.411*unifiedAtomicMass;
mSpin = 1./2.;
break;
case 63: // Cu
mZ = 29;
mRadius = 4.214;
mSurfaceThickness = 0.586;
mOmega = 0;
mRho0 = 0.1701;
mName = "Cu";
mMass = 63.546*unifiedAtomicMass;
mSpin = 3./2.;
break;
case 40: // Ca
mZ = 20;
mRadius = 3.766;
mSurfaceThickness = 0.586;
mOmega = -0.161;
mRho0 = 0.1699;
mName = "Ca";
mMass= 40.078*unifiedAtomicMass;
mSpin = 7./2.;
break;
case 27: // Al
mZ = 13;
mRadius = 3.07;
mSurfaceThickness = 0.519;
mOmega = 0;
mRho0 = 0.1739;
mName = "Al";
mMass = 26.981538*unifiedAtomicMass;
mSpin = 5./2.;
break;
case 16: // O
mZ = 8;
mRadius = 2.608;
mSurfaceThickness = 0.513;
mOmega = -0.051;
mRho0 = 0.1654;
mName = "O";
mMass = 15.9994*unifiedAtomicMass;
mSpin = 5./2.;
break;
case 2: // D
mZ = 1;
mHulthenA = 0.456;
mHulthenB = 2.36;
mRho0 = 0.39946312;
mName = "D";
mMass = 2.014102*unifiedAtomicMass;
mSpin = 1;
mRadius = 2.116; // not used for density function - taken from Jager et al.
break;
case 1:
// When generating this nucleus, only one nucleon is put at (0, 0, 0)
mZ = 1;
mRadius = 1.;
mMass = protonMass; // 1.00794*unifiedAtomicMass
mName = "p";
//The following 3 are dummies, needed for radial distribution etc.:
mSurfaceThickness = 0.513;
mOmega = -0.051;
mRho0 = 0.1654;
mSpin = 1./2.;
break;
default:
cout << "Nucleus::init(): Error, cannot handle A=" << mA << ", no parameters defined for this mass number." << endl;
exit(1);
break;
}
//
// Generate lookup table for overlap integral T
// (b and z in fm)
//
double range = 3*mRadius;
TF1* funcToIntegrate = new TF1("funcToIntegrate",this, &Nucleus::rhoForIntegration, -range, range, 1);
int nbins = 5000; // size of lookup table
mLookupTable = unique_ptr(new TH1D("mLookupTable","T lookup table", nbins, 0, range));
mLookupTable->SetDirectory(0);
for (int i=1; i<=nbins; i++) {
double b = mLookupTable->GetBinCenter(i);
funcToIntegrate->SetNpx(1000);
funcToIntegrate->SetParameter(0, b);
double res = funcToIntegrate->Integral(-range, range, 1e-8);
mLookupTable->SetBinContent(i, res);
}
delete funcToIntegrate;
}
unsigned int Nucleus::A() const { return mA; }
unsigned int Nucleus::Z() const { return mZ; }
float Nucleus::spin() const { return mSpin; }
double Nucleus::atomicMass() const { return mMass; }
string Nucleus::name() const { return mName; }
double Nucleus::radius() const { return mRadius; }
double Nucleus::rho(double b, double z) // returns value in units of fm^-3
{
// b and z in fm
double r = sqrt(b*b+z*z);
if (mA == 2) {
double term = (exp(-mHulthenA*r)/r) - (exp(-mHulthenB*r)/r); // Hulthen
return mRho0*term*term;
}
else {
return mRho0*(1+mOmega*(r/mRadius)*(r/mRadius))/(1+exp((r-mRadius)/mSurfaceThickness)); // 3pF
}
}
double Nucleus::rhoForIntegration(double *x, double* par)
{
// b and z in fm
double b = par[0];
double z = *x;
return rho(b, z);
}
double Nucleus::T(double b) const
{
//
// Returns overlap integral
// b in fm, overlap function T in GeV^2
// Returns 0 if b exceeds table range
//
if (mA == 0) {
cout << "Nucleus::T(): Error, instance is not initialized - cannot calculate overlap integral." << endl;
return 0;
}
int bin = mLookupTable->FindBin(b);
double res = mLookupTable->GetBinContent(bin);
return res*hbarc2;
}
double Nucleus::TForIntegration(double *x, double*) const
{
double b = x[0];
return 2*b*M_PI*T(b)/hbarc2;
}
void Nucleus::normalizationOfT(double eps)
{
//
// This is for internal checks only.
// Normalization of rho/T is such that fully integrated
// function should yield A.
//
double range = 3*mRadius;
TF1* func = new TF1("func",this, &Nucleus::TForIntegration, 0, range, 0);
double res = func->Integral(0, range, eps);
cout << "Normalization of T: = " << res << endl;
delete func;
}
double Nucleus::rho0() const { return mRho0; }
double Nucleus::TofProton(double b)
{
//
// Gaussian shape for proton
// Units:
// b in fm
//
const double BG = 4; // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
return 1/(2*M_PI*BG) * exp(-arg);
}
int Nucleus::pdgID(int Z, int A) const
{
// PDG for nuclei 10LZZZAAAI
// L = number of strange quark (for hypernuclei)
// I = isomer level, with I = 0 corresponding to the
// ground state and I > 0 to excitations,
if (Z==1 && A==1)
return 2212;
else if (Z==0 && A==1)
return 2112;
else if (A > 1)
return 1000000000 + 10000*Z + 10*A;
else
return 0;
}
int Nucleus::pdgID() const
{
return pdgID(mZ, mA);
}
Index: trunk/src/TableGeneratorSettings.cpp
===================================================================
--- trunk/src/TableGeneratorSettings.cpp (revision 432)
+++ trunk/src/TableGeneratorSettings.cpp (revision 433)
@@ -1,206 +1,206 @@
//==============================================================================
// TableGeneratorSettings.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: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "Settings.h"
#include "TableGeneratorSettings.h"
#include "Constants.h"
#include
#include
#include
#include
using namespace std;
TableGeneratorSettings* TableGeneratorSettings::mInstance = 0; // initialize static
TableGeneratorSettings* TableGeneratorSettings::instance()
{
if (mInstance == 0)
mInstance = new TableGeneratorSettings;
return mInstance;
}
TableGeneratorSettings::TableGeneratorSettings()
{
//
// Register all the parameters that can be defined
// via a runcard.
// Arguments for registerParameter():
// 1. pointer to data memeber
// 2. text string to be used in the runcard
// 3. default parameter set
//
registerParameter(&mBSatLookupPath, "bSatLookupPath", string("./"));
registerParameter(&mTmin, "tmin", -2.);
registerParameter(&mTmax, "tmax", 0.);
registerParameter(&mXmin, "xmin", 1e-9); //UPC
registerParameter(&mXmax, "xmax", 3e-2); //UPC
registerParameter(&mQ2bins, "Q2bins", static_cast(1));
registerParameter(&mW2bins, "W2bins", static_cast(1));
registerParameter(&mTbins, "tbins", static_cast(1));
registerParameter(&mXbins, "xbins", static_cast(1)); //UPC
registerParameter(&mNumberOfConfigurations, "numberOfConfigurations", static_cast(1000));
vector vec;
registerParameter(&mDipoleModelCustomParameters, "dipoleModelCustomParameters", vec);
registerParameter(&mUseBackupFile, "useBackupFile", false);
registerParameter(&mStartingBinFromBackup, "startingBinFromBackup", 0);
registerParameter(&mStartBin, "startBin", -1);
registerParameter(&mEndBin, "endBin", -1);
registerParameter(&mModesToCalculate, "modesToCalculate", 0);
registerParameter(&mPriority, "priority", 0);
}
void TableGeneratorSettings::consolidateSettings() // called after runcard is read
{
//
// Kinematic limits
//
if (mQ2min>=mQ2max && !mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, Q2min >= Q2max. Stopping" << endl;
exit(1);
}
if (mWmin>=mWmax && !mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, Wmin >= Wmax. Stopping" << endl;
exit(1);
}
if (mTmin>=mTmax) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, tmin >= tmax. Stopping" << endl;
exit(1);
}
if (mTmin>0. || mTmax >0.) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, t must be negative, please change t-limits. Stopping" << endl;
exit(1);
}
if (mXmin>=mXmax && mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, xmin >= xmax. Stopping" << endl;
exit(1);
}
- if ((mXmin>=1 || mXmin<=0 || mXmax>=1 || mXmax<=0) and mUPC) {
+ if ((mXmin>=1 || mXmin<=0 || mXmax>=1 || mXmax<=0) && mUPC) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, xmin or xmax out of range. Stopping" << endl;
exit(1);
}
if (mXmax>0.01 && mUPC){
cout << "TableGeneratorSettings::consolidateSettings(): Warning, xmax>1e-2, model may be unreliable." << endl;
}
if (mA==1) mNumberOfConfigurations = 1;
if (!mUseBackupFile) mStartingBinFromBackup = 0;
if (!mUPC)
mXbins=1;
else {
mQ2bins=1;
mW2bins=1;
}
if (mStartBin >= 0 && mEndBin >= 0 && mStartBin>mEndBin) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, endBin < startBin : " << mEndBin << " <= " << mStartBin << "! Stopping." << endl;
exit(1);
}
if ( mStartBin < 0 ) mStartBin=0;
if ( mStartBin >= signed(mQ2bins*mW2bins*mTbins*mXbins) ) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, starting bin >= table! Stopping." << endl;
exit(1);
}
if ( mEndBin > signed(mQ2bins*mW2bins*mTbins*mXbins) || mEndBin < 0) {
cout << "TableGeneratorSettings::consolidateSettings(): endBin is set to table size=" << mQ2bins*mW2bins*mTbins << endl;
mEndBin=mQ2bins*mW2bins*mTbins*mXbins;
}
if ( mModesToCalculate < 0 || mModesToCalculate > 2 ) {
cout << "TableGeneratorSettings::consolidateSettings(): Error, modesToCalculate can only take values 0, 1, or 2; not "
<< mModesToCalculate << endl;
exit(1);
}
//
// Make sure the W range is allowed
//
double VMMass=lookupPDG(mVectorMesonId)->Mass();
double W2min=VMMass*VMMass+protonMass2;
double Wmin=sqrt(W2min);
- if(mWmin
#include
#include "Math/RootFinder.h"
#include "Math/BrentRootFinder.h"
#include "Math/GSLRootFinder.h"
#include "Math/RootFinderAlgorithms.h"
#include "Math/IFunction.h"
#include "Math/BrentMinimizer1D.h"
// TMP
#include "TH1D.h"
#include "TFile.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
bool Kinematics::mError = false;
bool Kinematics::error() {return mError;}
//UPC:
double Kinematics::mXpomMin=0;
double Kinematics::mTold=0;
bool Kinematics::mXpomMinIsEvaluated=false;
//
double Kinematics::xprobe(double Q2, double W2, double vmMass)
{
mError = false;
double val = ( Q2 + vmMass*vmMass ) / (W2 - protonMass2 + Q2);
if (val > 1 || val < 0) mError = true;
return val;
}
TLorentzVector Kinematics::electronBeam(double eE, bool upc)
{
mError = false;
if (upc)
return TLorentzVector(0, 0, -sqrt(eE*eE-protonMass2), eE);
else
return TLorentzVector(0, 0, -sqrt(eE*eE-electronMass2), eE);
}
TLorentzVector Kinematics::hadronBeam(double hE)
{
mError = false;
return TLorentzVector(0, 0, sqrt(hE*hE-protonMass2), hE);
}
double Kinematics::s(const TLorentzVector& e, const TLorentzVector& p)
{
mError = false;
return (e+p)*(e+p);
}
double Kinematics::s(double eE, double hE, bool upc)
{
mError = false;
TLorentzVector e = electronBeam(eE, upc);
TLorentzVector p = hadronBeam(hE);
return (e+p)*(e+p);
}
double Kinematics::x(double Q2, double W2)
{
mError = false;
double val = Q2/(W2-protonMass2+Q2);
if (val > 1 || val < 0) mError = true;
return val;
}
double Kinematics::y(double Q2, double x, double s)
{
mError = false;
double val = Q2/(x*(s-electronMass2-protonMass2));
if (val > 1 || val < 0) mError = true;
return val;
}
double Kinematics::ymin(double s, double vmMass)
{
mError = false;
double W2min = Kinematics::W2min(vmMass);
double val = (s + W2min - sqrt((s-W2min)*(s-W2min)-4*electronMass2*W2min))/2/(s+electronMass2);
if (val > 1 || val < 0) mError = true;
return val;
}
double Kinematics::ymax(double s, double vmMass)
{
mError = false;
double W2min = Kinematics::W2min(vmMass);
double val = (s + W2min + sqrt((s-W2min)*(s-W2min)-4*electronMass2*W2min))/2/(s+electronMass2);
if (val > 1 || val < 0) mError = true;
return val;
}
double Kinematics::W2(double Q2, double x)
{
//
// Returns W2 at a given Q2 and Bjorken x
//
mError = false;
return protonMass2 + Q2*(1-x)/x;
}
double Kinematics::W2(double Q2, double xprobe, double vmMass)
{
//
// Returns W2 at a given Q2 and xprobe
// where xprobe is
// xprobe = xBJ * (1 + MV^2/Q2)
//
mError = false;
return protonMass2 + vmMass*vmMass/xprobe + Q2*(1-xprobe)/xprobe;
}
double Kinematics::W2min(double vmMass)
{
mError = false;
return vmMass*vmMass + protonMass2;
}
double Kinematics::W2max(double s)
{
mError = false;
return s; // check !?
}
double Kinematics::Q2min(double y)
{
mError = false;
return electronMass2*y*y/(1-y);
}
double Kinematics::Q2max(double s)
{
mError = false;
return s-electronMass2-protonMass2;
}
double Kinematics::xpomeron(double t, double Q2, double W2, double vmM)
{
mError = false;
double val = (vmM*vmM + Q2 - t)/(W2 + Q2 - protonMass2); // full expression
if (val > 1 || val < 0) mError = true;
return val;
}
double Kinematics::tmax(double xpom)
{
mError = false;
return (- xpom*xpom * protonMass2 / (1-xpom)); // only if m_p (in) = m_p' (out) - to be fixed
}
double Kinematics::tmax(double t, double Q2, double W2, double vmM)
{
return tmax(xpomeron(t, Q2, W2, vmM));
}
double Kinematics::tmin(double eE)
{
mError = false;
return -2*(eE*protonMass-protonMass2); // check ?!
}
//
// UPC only
//
double Kinematics::xpomMin(double massVM, double t, TLorentzVector hBeam, TLorentzVector eBeam, double MY2minusM2){
- if(mXpomMinIsEvaluated and mTold == t) return mXpomMin;
+ if (mXpomMinIsEvaluated && mTold == t) return mXpomMin;
LowerXpomeronFormula formula;
mTold = t;
formula.mT = t;
formula.mVmMass2=massVM*massVM;
formula.mElectronBeam = eBeam;
formula.mHadronBeam = hBeam;
formula.mMY2minusM2 = MY2minusM2;
//
//Find starting brackets
//
double lower=1e-12;
double upper=1e-1;
//
// Find minimum and use as lower
// bracket to avoid two roots
//
ROOT::Math::BrentMinimizer1D bm;
bm.SetFunction(formula, log(lower), log(upper));
bm.Minimize(100000, 1e-8, 1e-10);
double theMin = bm.XMinimum();
//
// Run root finder
//
ROOT::Math::RootFinder rootfinder(ROOT::Math::RootFinder::kBRENT);
rootfinder.SetFunction(formula, theMin, log(upper));
rootfinder.Solve(100000000, 1e-14, 0);
//
// Now some fine adjustements
//
double bestGuess = rootfinder.Root();
int ntimes = 0;
while (formula.DoEval(bestGuess) < 0) {
bestGuess += 1e-14;
ntimes++;
if (ntimes> 10000) break;
}
double result = exp(bestGuess);
double s = (hBeam+eBeam).M2();
mXpomMin = max(result, xpomMin(s, massVM, t));
mXpomMinIsEvaluated = true;
return mXpomMin;
}
double Kinematics::xpomMin(double s, double vmM, double t) {
//
// This is the threshold value.
// In some cases a more restrictive value may be needed
//
return (electronMass2 + protonMass2 + vmM*vmM - t)/s;
}
double Kinematics::Egamma(double xpom, double t, double vmM, double hBeamEnergy, double eBeamEnergy, double MY2minusM2)
{
//
// We use that the vector meson = photon + pomeron to calculate Egamma
// (keeping the pomeron variables from above)
// We solve a quadratic equation to get the photon energy such that
// the vector meson has the correct mass.
//
double Ep=hBeamEnergy;
double Pp=-sqrt(Ep*Ep-protonMass2);
double Ee=eBeamEnergy;
double Pe=sqrt(Ee*Ee-protonMass2);
double sqrtarg= 1 - ( (2*xpom - xpom*xpom)*Pp*Pp + t - MY2minusM2) / (Ep*Ep);
double Epom=Ep*(1-sqrt(sqrtarg)); //this is close to xpom*Ep
double Ppom=xpom*Pp;
double AA = vmM*vmM - Epom*Epom - t + Ppom*Ppom + 2*Pe*(Pe+Ppom); //GeV^2
double aa = 4*(Pe+Ppom)*(Pe+Ppom) - 4*(Ee+Epom)*(Ee+Epom); //GeV^2
double bb = 4*AA*(Ee+Epom) - 8*(Pe+Ppom)*(Pe+Ppom)*Ee; //GeV^3
double cc = 4*(Pe+Ppom)*(Pe+Ppom)*Pe*Pe - AA*AA;
sqrtarg = bb*bb-4*aa*cc;
if (sqrtarg < 0) {
//This is used by validUPC as a sign of failure, needs to return something negative
return -1.;
}
return ( -bb + sqrt(sqrtarg) )/(2*aa); //GeV
}
//
// UPC only
//
bool Kinematics::validUPC(double hBeamEnergy, double eBeamEnergy, double t, double xpom, double vmMass, bool verbose)
{
//
// Here we perform all complete test of all kinematic variables
//
//
// xpom
//
double s = Kinematics::s(eBeamEnergy, hBeamEnergy, true);
double Egam = Egamma(xpom, t, vmMass, hBeamEnergy, eBeamEnergy);
double minxpom = xpomMin(s, vmMass, t);
if (xpom < minxpom) {
- if(verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
+ if (verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
<< " because xpom < xpomMin (" << minxpom << ")" << endl;
return false;
}
if (Egam < 0) {
- if(verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
+ if (verbose) cout << "Kinematics::validUPC(): reject t=" << t << " and xpom=" << xpom
<< " because Egamma < 0 (" << Egam << ")" << endl;
return false;
}
if (xpom > 1) {
- if(verbose) cout << "Kinematics::validUPC(): reject xpom="< 1" << endl;
+ if (verbose) cout << "Kinematics::validUPC(): reject xpom="< 1" << endl;
return false;
}
//
// t
//
double t_max=tmax(xpom);
if (t > t_max) {
if (verbose) cout << "Kinematics::validUPC(): reject t=" << t << " because t > tmax (" << t_max << ")" << endl;
return false;
}
return true; //survived all tests
}
bool Kinematics::valid(double s, double t, double Q2, double W2, double vmMass,
bool useTrueXp, bool verbose)
{
//
// Here we perform all complete test of all kinematic variables
//
double x = Kinematics::x(Q2, W2);
double y = Kinematics::y(Q2, x, s);
//
// W
//
if (W2 < Kinematics::W2min(vmMass)) {
if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W < Wmin" << endl;
return false;
}
if (W2 > Kinematics::W2max(s)) {
if (verbose) cout << "Kinematics::valid(): reject W=" << sqrt(W2) << " because W > Wmax" << endl;
return false;
}
//
// y range
//
double ymin = Kinematics::ymin(s, vmMass);
double ymax = Kinematics::ymax(s, vmMass);
if (y < ymin || y > ymax) {
if (verbose) cout << "Kinematics::valid(): reject y=" << y << " because y < ymin (" << ymin
<< ") || y > ymax (" << ymax << ")" << endl;
return false;
}
//
// x and y
//
if (x < 0 || x > 1) {
if (verbose) cout << "Kinematics::valid(): reject x=" << x << " because x < 0 || x > 1" << endl;
return false;
}
//
// Q2
//
if (Q2 > Kinematics::Q2max(s)) {
if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 > Q2max" << endl;
return false;
}
if (Q2 < Kinematics::Q2min(y)) {
if (verbose) cout << "Kinematics::valid(): reject Q2=" << Q2 << " because Q2 < Q2min" << endl;
return false;
}
//
// t
//
// On calculating xpomeron. xpomeron() provides the correct value
// (not neglecting t or nucleon mass). However, during initialization
// to avoid recursion we set t=0. Calling valid() with t != 0 would
// reject these values. In this case useTrueXp is set to false.
// For final checks if an event has the correct kinematics we have
// to set useTrueXp = true.
//
double xpom;
if (useTrueXp)
xpom = Kinematics::xpomeron(t, Q2, W2, vmMass);
else
xpom = Kinematics::xpomeron(0, Q2, W2, vmMass);
double tmax = Kinematics::tmax(xpom);
if (t > tmax) {
if (verbose) cout << "Kinematics::valid(): reject t=" << t << " because t > tmax (" << tmax << ")" << endl;
return false;
}
//
// xpom
//
if (xpom < 0 || xpom > 1) {
if (verbose) cout << "Kinematics::valid(): reject xpom=" << xpom << " because xpom < 0 || xpom > 1" << endl;
return false;
}
if (x > xpom) {
if (verbose) cout << "Kinematics::valid(): reject since x=" << x << " > xpom = " << xpom << endl;
return false;
}
return true; // survived all tests
}
ROOT::Math::IBaseFunctionOneDim* LowerXpomeronFormula::Clone() const
{
return new LowerXpomeronFormula();
}
double LowerXpomeronFormula::DoEval(double logxpom) const
{
double xpom = exp(logxpom);
//Start by setting up the knowns: Incoming beams, xpom, and t:
double Ee=mElectronBeam.E();
double Pe=mElectronBeam.Pz();
double Pp=mHadronBeam.Pz();
double Ep=mHadronBeam.E();
double E, pz;
//
// Use scattered proton to set up four momentum of pomeron (so far assuming no break-up):
//
double sqrtarg= 1 - ( (2*xpom - xpom*xpom)*Pp*Pp + mT - mMY2minusM2) / (Ep*Ep);
E=Ep*(1-sqrt(sqrtarg)); //this is close to xpom*Ep
pz=xpom*Pp;
//
// Next we use that the vector meson = photon + pomeron to calculate Egamma
// (keeping the pomeron variables from above)
// We solve a quadratic equation to get the photon energy such that
// the vector meson has the correct mass.
//
double AA = mVmMass2 - E*E - mT + pz*pz + 2*Pe*(Pe+pz); //GeV^2
double aa = 4*(Pe+pz)*(Pe+pz) - 4*(Ee+E)*(Ee+E); //GeV^2
double bb = 4*AA*(Ee+E) - 8*(Pe+pz)*(Pe+pz)*Ee; //GeV^3
double cc = 4*(Pe+pz)*(Pe+pz)*Pe*Pe - AA*AA;
sqrtarg = bb*bb-4*aa*cc;
return sqrtarg;
}
Index: trunk/src/CMakeLists.txt
===================================================================
--- trunk/src/CMakeLists.txt (revision 432)
+++ trunk/src/CMakeLists.txt (revision 433)
@@ -1,156 +1,156 @@
#===============================================================================
# CMakeLists.txt (src)
#
# 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(ExternalProject)
cmake_minimum_required (VERSION 3.1)
#
# Compiler flags for release and debug version
#
set(CMAKE_C_FLAGS "-W")
set(CMAKE_C_FLAGS_DEBUG "-g")
set(CMAKE_C_FLAGS_RELEASE "-O")
message("COMPILER = ${CMAKE_CXX_COMPILER_ID}")
set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -pedantic -Wno-long-long")
-if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-potentially-evaluated-expression")
endif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS_DEBUG "-g")
set(CMAKE_CXX_FLAGS_RELEASE "-O")
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
#
# Find external required packages
# (see also FindGSL.cmake and FindROOT.cmke in cmake/modules)
# Herer we need only the root library to create the table
# tools.
#
# GSL
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIR})
# ROOT
find_package(ROOT REQUIRED)
include_directories(${ROOT_INCLUDE_DIR})
set(LIBS ${LIBS} ${ROOT_LIBRARIES})
#BOOST
if (MULTITHREADED)
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED ON)
find_package(Boost 1.39 COMPONENTS thread REQUIRED)
- if(Boost_FOUND)
+ if (Boost_FOUND)
include_directories(${Boost_INCLUDE_DIR})
endif(Boost_FOUND)
endif (MULTITHREADED)
#
# Include files from sartre package
#
include_directories(${PROJECT_SOURCE_DIR}/src)
include_directories(${PROJECT_SOURCE_DIR}/gemini)
include_directories(${PROJECT_SOURCE_DIR}/cuba)
#
# Cuba is built using the autoconf shipped with it
#
ExternalProject_Add(
cuba
DOWNLOAD_COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/cuba ${PROJECT_BINARY_DIR}/cuba
SOURCE_DIR ${PROJECT_BINARY_DIR}/cuba
PREFIX ${PROJECT_BINARY_DIR}/cuba
CONFIGURE_COMMAND ${PROJECT_BINARY_DIR}/cuba/configure --prefix=${PROJECT_BINARY_DIR}/cuba/build
BUILD_COMMAND make lib
BUILD_IN_SOURCE 1
)
#
# Defines source files for sartre library
#
set(SARTRE_SRC "AlphaStrong.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Amplitudes.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "BreakupProduct.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "CrossSection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DglapEvolution.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModel.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModelParameters.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Event.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ExclusiveFinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FrangibleNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "EventGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Integrals.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Kinematics.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ModeFinderFunctor.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleon.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "PhotonFlux.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Sartre.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Settings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Table.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableCollection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "WaveOverlap.cpp")
add_library(sartre ${SARTRE_SRC})
#
# Table tools (all stand alone programs)
#
add_executable(tableInspector tableInspector.cpp)
add_executable(tableMerger tableMerger.cpp)
add_executable(tableQuery tableQuery.cpp)
add_executable(tableDumper tableDumper.cpp)
add_executable(tableVarianceMaker tableVarianceMaker.cpp)
target_link_libraries(tableInspector sartre ${LIBS})
target_link_libraries(tableMerger sartre ${LIBS})
target_link_libraries(tableQuery sartre ${LIBS})
target_link_libraries(tableDumper sartre ${LIBS})
target_link_libraries(tableVarianceMaker sartre ${LIBS})
add_dependencies(tableInspector sartre)
add_dependencies(tableMerger sartre)
add_dependencies(tableQuery sartre)
add_dependencies(tableDumper sartre)
add_dependencies(tableVarianceMaker sartre)
#
# Install library and include files (make install) within
# the distribution tree. Top level CMakeLists.txt will install
# Sartre in final destination.
#
install(TARGETS sartre DESTINATION sartre/lib)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/lib/libcuba.a" DESTINATION sartre/lib)
FILE(GLOB AllIncludeFiles *.h)
install(FILES ${AllIncludeFiles} DESTINATION sartre/include)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/include/cuba.h" DESTINATION sartre/include)
install(TARGETS tableInspector DESTINATION sartre/bin)
install(TARGETS tableMerger DESTINATION sartre/bin)
install(TARGETS tableQuery DESTINATION sartre/bin)
install(TARGETS tableDumper DESTINATION sartre/bin)
install(TARGETS tableVarianceMaker DESTINATION sartre/bin)
Index: trunk/src/DipoleModel.cpp
===================================================================
--- trunk/src/DipoleModel.cpp (revision 432)
+++ trunk/src/DipoleModel.cpp (revision 433)
@@ -1,337 +1,337 @@
//==============================================================================
// DipoleModel.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include
#include
#include
#include
#include "DipoleModel.h"
#include "TableGeneratorSettings.h"
#include "DglapEvolution.h"
#include "Constants.h"
#include "TFile.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH2F.h"
#include "TF1.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
DipoleModel::DipoleModel()
{
- mConfigurationExists=false;
+ mConfigurationExists = false;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = settings->A();
mNucleus.init(A);
- mIsInitialized=true;
- mParameters = 0;
+ mIsInitialized = true;
+ mParameters = nullptr;
}
DipoleModel::~DipoleModel()
{
delete mParameters;
}
const TableGeneratorNucleus* DipoleModel::nucleus() const { return &mNucleus; }
bool DipoleModel::configurationExists() const { return mConfigurationExists; }
double DipoleModel::bDependence(double) { return 0; }
double DipoleModel::bDependence(double, double) { return 0; }
double DipoleModel::dsigmadb2ep(double, double, double) { return 0;}
//***********bSat:*****************************************************
DipoleModel_bSat::DipoleModel_bSat()
{
mBDependence = 0;
mSigma_ep_LookupTable = 0;
//
// Set the parameters. Note that we enforce here the bSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bSat, settings->dipoleModelParameterSet());
}
DipoleModel_bSat::~DipoleModel_bSat()
{
delete mSigma_ep_LookupTable;
delete mBDependence;
}
DipoleModel_bSat& DipoleModel_bSat::operator=(const DipoleModel_bSat& dp)
{
if (this != &dp) {
delete mBDependence;
delete mSigma_ep_LookupTable;
DipoleModel::operator=(dp);
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable));
mSigma_ep_LookupTable->SetDirectory(0);
}
return *this;
}
DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp)
{
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
}
void DipoleModel_bSat::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = mNucleus.A();
- string path=settings->bSatLookupPath();
+ string path = settings->bSatLookupPath();
ostringstream filename;
filename.str("");
filename << path << "/bSat_bDependence_A" << A <<".root";
ifstream ifs(filename.str().c_str());
if (!ifs) {
cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl;
cout << "Stopping." << endl;
exit(1);
}
TFile* lufile= new TFile(filename.str().c_str());
ostringstream histoName;
histoName.str( "" );
histoName << "Configuration_" << iConfiguration;
lufile->GetObject( histoName.str().c_str(), mBDependence );
mBDependence->SetDirectory(0);
lufile->Close();
- mConfigurationExists=true;
+ mConfigurationExists = true;
}
double DipoleModel_bSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- double bDep=bDependence(b, phi);
+ double bDep = bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
double DipoleModel_bSat::bDependence(double b, double phi)
{
double result = mBDependence->Interpolate(b, phi);
return result;
}
double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe) {
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
xprobe*=1.;
double sigmap = mSigma_ep_LookupTable->Interpolate(r);
- int A=nucleus()->A();
- double TA=nucleus()->T(b)/A;
+ int A = nucleus()->A();
+ double TA = nucleus()->T(b)/A;
double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) );
return result;
}
void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe)
{
- double rbRange=3.*nucleus()->radius();
+ double rbRange = 3.*nucleus()->radius();
TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2);
mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange);
dsigmaForIntegration->SetNpx(1000);
for (int iR=1; iR<=1000; iR++) {
- double r=mSigma_ep_LookupTable->GetBinCenter(iR);
+ double r = mSigma_ep_LookupTable->GetBinCenter(iR);
dsigmaForIntegration->SetParameter(0, r);
dsigmaForIntegration->SetParameter(1, xprobe);
- double result=dsigmaForIntegration->Integral(0, rbRange);
+ double result = dsigmaForIntegration->Integral(0, rbRange);
mSigma_ep_LookupTable->SetBinContent(iR, result);
}
delete dsigmaForIntegration;
}
double DipoleModel_bSat::dsigmadb2epForIntegration(double *x, double* par)
{
- double r=par[0];
- double xprobe=par[1];
+ double r = par[0];
+ double xprobe = par[1];
double b = *x;
return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe);
}
//***********bNonSat:*************************************************
DipoleModel_bNonSat::DipoleModel_bNonSat()
{
//
// Set the parameters. Note that we need bNonSat to calculate
// the skewedness correction for bSat.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(settings->dipoleModelType(), settings->dipoleModelParameterSet());
}
DipoleModel_bNonSat::~DipoleModel_bNonSat(){}
double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- double bDep=bDependence(b, phi);
+ double bDep = bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
- int A=nucleus()->A();
- double TA=nucleus()->T(b)/A;
+ int A = nucleus()->A();
+ double TA = nucleus()->T(b)/A;
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
- double result=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
+ double result = A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
return result;
}
//***********bCGC:*****************************************************
DipoleModel_bCGC::DipoleModel_bCGC()
{
//
// Set the parameters. Note that we enforce here the bNonSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bCGC, settings->dipoleModelParameterSet());
}
void DipoleModel_bCGC::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
//iConfiguration is a dummy for bCGC. Do this to avoid compiler warnings:
iConfiguration++;
mNucleus.generate();
- mConfigurationExists=true;
+ mConfigurationExists = true;
}
double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x)
{
- double result=1;
+ double result = 1;
for (unsigned int iA=0; iAkappa();
double N0 = mParameters->N0();
double x0 = mParameters->x0();
double lambda = mParameters->lambda();
double gammas = mParameters->gammas();
double A = -N0*N0*gammas*gammas/((1-N0)*(1-N0)*log(1-N0));
double B = 0.5*pow(1-N0,-(1-N0)/(N0*gammas));
double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b));
double rQs = r*Qs/hbarc;
- double result=0;
+ double result = 0;
if (rQs <= 2)
result = 2*N0*pow(0.5*rQs, 2*(gammas+(1/(kappa*lambda*Y))*log(2/rQs)));
else
result = 2*(1 - exp(-A*log(B*rQs)*log(B*rQs)));
return result;
}
double DipoleModel_bCGC::bDependence(double b)
{
double gammas = mParameters->gammas();
double Bcgc = mParameters->Bcgc();
return exp(-0.5*b*b/Bcgc/gammas/hbarc2);
}
Index: trunk/src/TableGeneratorNucleus.cpp
===================================================================
--- trunk/src/TableGeneratorNucleus.cpp (revision 432)
+++ trunk/src/TableGeneratorNucleus.cpp (revision 433)
@@ -1,197 +1,197 @@
//==============================================================================
// TableGeneratorNucleus.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 "TableGeneratorNucleus.h"
#include "TableGeneratorSettings.h"
#include "TH1D.h"
#include
#include
TableGeneratorNucleus::TableGeneratorNucleus()
{
mRadialDistributionHistogram = 0;
}
TableGeneratorNucleus::TableGeneratorNucleus(unsigned int A) : Nucleus(A)
{
//
// Generate a radial distribution histograms according to the Woods Saxon potential * Volume:
// dN/dr=4*pi*r^2*dN/dV; dN/dV=rho:
//
int numRadialBins = 10000;
mRadialDistributionHistogram = new TH1D("mRadialDistributionHistogram",
"Woods Saxon for Random Generator",
numRadialBins, 0., 3.*mRadius);
for (int i=1; i <= numRadialBins; i++) {
double radius=mRadialDistributionHistogram->GetBinCenter(i);
mRadialDistributionHistogram->SetBinContent(i, 4.*M_PI*radius*radius*rho(radius, 0.));
}
}
TableGeneratorNucleus& TableGeneratorNucleus::operator=(const TableGeneratorNucleus& n)
{
if (this != &n) {
delete mRadialDistributionHistogram;
configuration.clear();
Nucleus::operator=(n); // copy assign for base class
mRadialDistributionHistogram = new TH1D(*(n.mRadialDistributionHistogram));
mRadialDistributionHistogram->SetDirectory(0);
copy(n.configuration.begin(), n.configuration.end(), configuration.begin());
}
return *this;
}
TableGeneratorNucleus::TableGeneratorNucleus(const TableGeneratorNucleus& n) : Nucleus(n)
{
mRadialDistributionHistogram = new TH1D(*(n.mRadialDistributionHistogram));
mRadialDistributionHistogram->SetDirectory(0);
copy(n.configuration.begin(), n.configuration.end(), configuration.begin());
}
TableGeneratorNucleus::~TableGeneratorNucleus()
{
delete mRadialDistributionHistogram;
}
const TH1D* TableGeneratorNucleus::getRHisto() const {return mRadialDistributionHistogram;}
bool TableGeneratorNucleus::generate() {
//
// This function generate a Nucleus in the format of a vector of nucleons of size A
// It generates the position of each nucleon according to the Woods-Saxon potential,
// or in the case of deuterium, the Hulthen potential.
//
// Must clean up each event such that the new Nucleus is not just added to the last one.
//
TRandom3 *random = TableGeneratorSettings::randomGenerator();
configuration.clear();
if (mA==2) { //Deuterium
Nucleon nucleon1, nucleon2;
//
// Generate the distance between the nucleons
//
double distance=mRadialDistributionHistogram->GetRandom();
//
// Generate x, y, and z given radius=distance/2
// and set nucleons position.
//
double x, y, z;
random->Sphere(x, y, z, distance/2.);
nucleon1.setPosition(TVector3(x, y, z));
nucleon2.setPosition(TVector3(-x, -y, -z));
//
// Add nucleons to the configuration:
//
configuration.push_back(nucleon1);
configuration.push_back(nucleon2);
}
else{
TVector3 centerOfMass;
Nucleon tmpNucleon;
double *radiusArray = new double [mA];
const double dCore=0.7; // core size of each nucleon
const double dCore2 = dCore*dCore;
//
// Generate radii according to Woods-Saxon*Volume and
// sort them for a more linear check of the distance between nucleons
for (unsigned int iR=0; iRGetRandom();
}
sort(radiusArray, radiusArray+mA);
for (unsigned int iA=0; iA 10) {
+ if (count > 10) {
delete [] radiusArray;
return false;
}
//
// Generate x, y, and z given radius
// and set nucleons position.
//
double x, y, z;
random->Sphere(x, y, z, radius);
tmpNucleon.setPosition(TVector3(x, y, z));
//
// Find out how many previous nucleons are within dCore
// from this one in radius...
//
unsigned int iTooClose=0;
unsigned int ii=1;
while ( iA > ii && (radius-radiusArray[iA-ii]) < dCore) {
iTooClose++;
ii++;
}
loop = false;
//
// Check if any of those are overlapping in 3D.
// If so regenerate the angles.
//
- for(unsigned int j=1; j<=iTooClose; j++){
- if((configuration[iA-j].position()-tmpNucleon.position()).Mag2() < dCore2) {
+ for (unsigned int j=1; j<=iTooClose; j++){
+ if ((configuration[iA-j].position()-tmpNucleon.position()).Mag2() < dCore2) {
loop = true;
break;
}
}
}//while loop
//
// A nucleon has been generated.
// Add it to the configuration and
// to the center of mass vector
//
configuration.push_back( tmpNucleon );
centerOfMass.SetXYZ((centerOfMass.X()+tmpNucleon.position().X())/mA,
(centerOfMass.Y()+tmpNucleon.position().Y())/mA,
(centerOfMass.Z()+tmpNucleon.position().Z())/mA);
}// iA loop
//
// Adjust the origin to the center of mass
//
for (unsigned int iA=0; iA.
//
-// Author: Thomas Ullrich
+// Author: Thomas Ullrich, Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "PhotonFlux.h"
#include "Constants.h"
#include "Kinematics.h"
#include "EventGeneratorSettings.h"
#include "Nucleus.h"
#include
#include "TF1.h"
#include "Math/Functor.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#define PR(x) cout << #x << " = " << (x) << endl;
-PhotonFlux::PhotonFlux() {
+PhotonFlux::PhotonFlux()
+{
mS = 0;
- mIsUPC=false;
+ mIsUPC = false;
mSettings = EventGeneratorSettings::instance();
- if(mSettings->UPC()){
+
+ if (mSettings->UPC()){
mNucleusUPC = new Nucleus(mSettings->UPCA());
mNucleus = new Nucleus(mSettings->A());
- if(mSettings->UPCA()>1 and mSettings->A()>1){
- cout<<"Calculating TAA table..."<UPCA()>1 && mSettings->A()>1) {
+ cout << "PhotonFlux::PhotonFlux(): Calculating TAA table ... ";
+ setupTAAlookupTable();
+ cout << "done." << endl;
}
- mEBeamEnergy=mSettings->electronBeamEnergy(); //The beam that shines
- mIsUPC=true;
+ mEBeamEnergy = mSettings->electronBeamEnergy(); // the beam that shines
+ mIsUPC = true;
}
- mSIsSet=false;
+
+ mSIsSet = false;
}
-PhotonFlux::PhotonFlux(double s) {
+PhotonFlux::PhotonFlux(double s)
+{
mS = s;
- mIsUPC=false;
+ mIsUPC = false;
mSettings = EventGeneratorSettings::instance();
+
if (mSettings->UPC()) {
mNucleusUPC = new Nucleus(mSettings->UPCA());
mNucleus = new Nucleus(mSettings->A());
- if(mSettings->UPCA()>1 and mSettings->A()>1){
- cout<<"Calculating TAA table..."<UPCA()>1 && mSettings->A()>1){
+ cout << "PhotonFlux::PhotonFlux(): Calculating TAA table ... ";
+ setupTAAlookupTable();
+ cout << "done" << endl;
}
- mEBeamEnergy=mSettings->electronBeamEnergy(); //The beam that shines
- mIsUPC=true;
- cout<<"Warning in PhotonFlux::PhotonFlux() Sartre is not yet ready for UPC..."<electronBeamEnergy(); //The beam that shines
+ mIsUPC = true;
+ cout<<"PhotonFlux::PhotonFlux(): Warning Sartre is not yet ready for UPC..." << endl;
}
- mSIsSet=true;
+
+ mSIsSet = true;
}
-void PhotonFlux::setS(double s) {
+void PhotonFlux::setS(double s)
+{
mS = s;
- mSIsSet=true;
+ mSIsSet = true;
}
double PhotonFlux::operator()(double Q2, double W2, GammaPolarization p) const
{
- if(!mSIsSet) cout<<"Warning in PhotonFlux::operator() s is not set!"<1 || Kinematics::error()) return 0;
if (Q2 > Kinematics::Q2max(mS)) return 0;
double result = alpha_em/(2*M_PI*Q2*mS*y);
result *= (1 + (1-y)*(1-y)) - (2*(1-y)*Kinematics::Q2min(y)/Q2);
return result;
}
double PhotonFlux::fluxLongitudinal(double Q2, double W2) const
{
//
// Longitudinal photon flux
//
double y = Kinematics::y(Q2, Kinematics::x(Q2, W2), mS);
if (y<0 || y>1 || Kinematics::error()) return 0;
if (Q2 > Kinematics::Q2max(mS)) return 0;
double result = alpha_em/(2*M_PI*Q2*mS*y);
result *= 2*(1-y);
return result;
}
double PhotonFlux::nuclearPhotonFlux(double Egamma) const
{
if (Egamma < 0) {
if (mSettings->verbose() && mSettings->verboseLevel() > 4)
cout << "PhotonFlux::nuclearPhotonFlux(): Warning, negative Egamma as argument.\n"
- << " Return 0 for photon flux." << endl;
+ << " Return 0 for photon flux." << endl;
return 0;
}
- double bmin=0.;
- double bmax=1e10; //Could be arbitrary large (TF1 doesn't care)
- TF1 fFluxAA("fluxAA", this, &PhotonFlux::uiNuclearPhotonFlux, bmin, bmax, 1);
- fFluxAA.SetParameter(0, Egamma);
- ROOT::Math::WrappedTF1 wfFluxAA(fFluxAA);
- ROOT::Math::GaussIntegrator giFluxAA;
- giFluxAA.SetFunction(wfFluxAA);
- giFluxAA.SetAbsTolerance(0.);
- giFluxAA.SetRelTolerance(1e-5);
+ double bmin = 0.;
+ double bmax = 1e10; //Could be arbitrary large (TF1 doesn't care)
+ TF1 fluxFunction("fluxFunction", this, &PhotonFlux::unintegratedNuclearPhotonFlux, bmin, bmax, 1);
+ fluxFunction.SetParameter(0, Egamma);
+ ROOT::Math::WrappedTF1 wrappedFluxFunction(fluxFunction);
+ ROOT::Math::GaussIntegrator fluxIntegrator;
+ fluxIntegrator.SetFunction(wrappedFluxFunction);
+ fluxIntegrator.SetAbsTolerance(0.);
+ fluxIntegrator.SetRelTolerance(1e-5);
- return giFluxAA.IntegralUp(bmin)/hbarc; //dN/dEgamma GeV-1
+ return fluxIntegrator.IntegralUp(bmin)/hbarc; //dN/dEgamma GeV^-1
}
-double PhotonFlux::uiNuclearPhotonFlux(double* var, double* par) const {
+double PhotonFlux::unintegratedNuclearPhotonFlux(double* var, double* par) const
+{
//
// https://arxiv.org/abs/1607.03838v1
//
- double b=var[0]; //fm
- double eGamma=par[0]; //GeV
+ double b = var[0]; //fm
+ double eGamma = par[0]; //GeV
- int Apom=mNucleus->A();
- int AUPC=mNucleusUPC->A();
+ int Apom = mNucleus->A();
+ int AUPC = mNucleusUPC->A();
- double Z=mNucleusUPC->Z();
- double ebeamMass=mNucleusUPC->atomicMass()/mNucleusUPC->A(); //Mass per nucleon //GeV
- double beamLorentzGamma=mEBeamEnergy/ebeamMass;
- double beamLorentzGamma2=beamLorentzGamma*beamLorentzGamma;
+ double Z = mNucleusUPC->Z();
+ double ebeamMass = mNucleusUPC->atomicMass()/mNucleusUPC->A(); //Mass per nucleon //GeV
+ double beamLorentzGamma = mEBeamEnergy/ebeamMass;
+ double beamLorentzGamma2 = beamLorentzGamma*beamLorentzGamma;
- double xi=eGamma*b/hbarc/beamLorentzGamma; //GeV0
+ double xi = eGamma*b/hbarc/beamLorentzGamma; //GeV0
- double K0=TMath::BesselK0(xi);
- double K1=TMath::BesselK1(xi);
+ double K0 = TMath::BesselK0(xi);
+ double K1 = TMath::BesselK1(xi);
- double prefactor=Z*Z*alpha_em/(M_PI*M_PI)*eGamma/beamLorentzGamma2; //GeV1
- double N=prefactor*(K1*K1+K0*K0/beamLorentzGamma2); //GeV1
+ double prefactor = Z*Z*alpha_em/(M_PI*M_PI)*eGamma/beamLorentzGamma2; //GeV1
+ double N = prefactor*(K1*K1+K0*K0/beamLorentzGamma2); //GeV1
//Spectrum is calculated under the assumption that there is no
//hadronic interaction between the beam particles.
//Therefore it has to be multiplied by the probability for this:
- double Pnohad=1;
- if(Apom>1 and AUPC>1){ //nucleus-nucleus interaction
- Pnohad=exp(-sigma_nn(mS)*mTAA_of_b->Interpolate(b)); //GeV0
+ double Pnohad = 1;
+ if (Apom>1 && AUPC>1) { //nucleus-nucleus interaction
+ Pnohad = exp(-sigma_nn(mS)*mTAA_of_b->Interpolate(b)); //GeV0
}
- else if(Apom==1 and AUPC>1){ //proton-nucleus
+ else if (Apom == 1 && AUPC > 1) { //proton-nucleus
Pnohad=exp(-sigma_nn(mS)*mNucleusUPC->T(b));
}
- else if(Apom>1 and AUPC==1){ //nucleus-proton
- Pnohad=exp(-sigma_nn(mS)*mNucleus->T(b));
+ else if (Apom > 1 && AUPC == 1) { //nucleus-proton
+ Pnohad = exp(-sigma_nn(mS)*mNucleus->T(b));
}
- else{ //proton-proton
- double b02=19.8; //GeV-2
- double scamp=1-exp(-b*b/hbarc2/2./b02);
- Pnohad=scamp*scamp;
+ else { //proton-proton
+ double b02 = 19.8; //GeV^-2
+ double scamp = 1-exp(-b*b/hbarc2/2./b02);
+ Pnohad = scamp*scamp;
}
- double result=2*M_PI*b/hbarc*N*Pnohad; //GeV0
+ double result = 2*M_PI*b/hbarc*N*Pnohad; //GeV0
return result;
}
-void PhotonFlux::calculateTAAlookupTable(){
- double rbhigh=upperIntegrationLimit*(mNucleus->radius()+mNucleusUPC->radius());
- mTAA_of_b=new TH1D("mTAA_of_b", "mTAA_of_b", 1000, 0, rbhigh);
- for(unsigned int i=1; i<=1000; i++){
- double b=mTAA_of_b->GetBinCenter(i);
+void PhotonFlux::setupTAAlookupTable()
+{
+ double rbhigh = upperIntegrationLimit*(mNucleus->radius()+mNucleusUPC->radius());
+ mTAA_of_b = new TH1D("mTAA_of_b", "mTAA_of_b", 1000, 0, rbhigh);
+ for (unsigned int i=1; i<=1000; i++) {
+ double b = mTAA_of_b->GetBinCenter(i);
mTAA_of_b->SetBinContent(i, TAA(b));
}
}
-double PhotonFlux::sigma_nn(double s) const{
- //parameters http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20080014212.pdf table 2
+double PhotonFlux::sigma_nn(double s) const {
+ //
+ // For parameters see
+ // http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20080014212.pdf table 2
//
// Following KNSGB: http://arxiv.org/abs/1607.03838v1 (eq.6)
//
- double Z=33.73; //mb
- double B=0.2838; //mb
- double s0=1.;//GeV2
- double Y1=13.67; //mb
- double s1=1.; //GeV2
- double eta1=0.412;
- double Y2=7.77; //mb
- double eta2=0.5626;
+ double Z = 33.73; //mb
+ double B = 0.2838; //mb
+ double s0 = 1.;//GeV2
+ double Y1 = 13.67; //mb
+ double s1 = 1.; //GeV2
+ double eta1 = 0.412;
+ double Y2 = 7.77; //mb
+ double eta2 = 0.5626;
double result = Z+B*log(s/s0)*log(s/s0)+Y1*pow(s1/s, eta1)-Y2*pow(s1/s, eta2); //mb = 1e-3b = 1e-1 fm2
result *= 1e-1; //fm2 (1 barn = 1e2 fm2)
return result/hbarc2; //GeV^-2
}
-double PhotonFlux::TAA(double b){
-
-
- //Integral over dr^2=2*pi*r*dr for
+double PhotonFlux::TAA(double b)
+{
+ //
+ // Integral over dr^2=2*pi*r*dr
+ //
double rbhigh=upperIntegrationLimit*(mNucleus->radius()+mNucleusUPC->radius());
- double lo[2]={0., 0.}; //r, theta
- double hi[2]={rbhigh, 2*M_PI};
+ double lo[2] = {0., 0.}; //r, theta
+ double hi[2] = {rbhigh, 2*M_PI};
ROOT::Math::Functor wf(this, &PhotonFlux::TAAForIntegration, 2);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig.SetFunction(wf);
ig.SetAbsTolerance(0.);
ig.SetRelTolerance(1e-4);
mB=b;
- double result=ig.Integral(lo, hi); //GeV^3*fm
- result/=hbarc; //GeV^2
+ double result = ig.Integral(lo, hi); //GeV^3*fm
+ result /= hbarc; //GeV^2
return result;
}
-double PhotonFlux::TAAForIntegration(const double* var) const{
- double r=var[0]; //fm
- double phi=var[1];
- double b=mB; //fm
+double PhotonFlux::TAAForIntegration(const double* var) const
+{
+ double r = var[0]; //fm
+ double phi = var[1];
+ double b = mB; //fm
//Put A1 in the origin, and A2 on the y-axis at distance b:
- double arg1=r;
- double arg2=sqrt(r*r+b*b-2*b*r*sin(phi)); //fm
+ double arg1 = r;
+ double arg2 = sqrt(r*r+b*b-2*b*r*sin(phi)); //fm
double T1 = mNucleusUPC->T(arg1); //GeV^2
double T2 = mNucleus->T(arg2); //GeV^2
return r/hbarc*T1*T2; //GeV^3
}
Index: trunk/src/CrossSection.cpp
===================================================================
--- trunk/src/CrossSection.cpp (revision 432)
+++ trunk/src/CrossSection.cpp (revision 433)
@@ -1,921 +1,921 @@
//==============================================================================
// 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;
}
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);
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 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_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/Constants.h
===================================================================
--- trunk/src/Constants.h (revision 432)
+++ trunk/src/Constants.h (revision 433)
@@ -1,49 +1,49 @@
//==============================================================================
// Constants.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 Constants_h
#define Constants_h
//
// General constants
//
const double electronMass = 0.510998902E-3; // GeV
const double electronMass2 = electronMass*electronMass; // GeV^2
const double protonMass = 0.9382700; // GeV
const double protonMass2 = protonMass*protonMass; // GeV^2
const double alpha_em = 1/137.036;
const double hbarc = 0.197327; // GeV*fm
const double hbarc2 = hbarc*hbarc; // GeV^2*fm^2
//
// Constants used in dipole model
//
-const double Nc=3.;
+const double Nc = 3.;
const double quarkCharge[6] = {2./3., -1./3., -1./3., 2./3., -1./3., 2./3.}; // u, d, s, c, b, t
//
// Constants used in integartion and other numerical operations
//
-const double upperIntegrationLimit=2.5; // factor: nuclear radius -> integration limits (b, r)
+const double upperIntegrationLimit = 2.5; // factor: nuclear radius -> integration limits (b, r)
#endif
Index: trunk/src/Table.cpp
===================================================================
--- trunk/src/Table.cpp (revision 432)
+++ trunk/src/Table.cpp (revision 433)
@@ -1,1501 +1,1501 @@
//==============================================================================
// Table.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Table data is stored in a TH3F.
//
// All information is stored in the histogram title.
//
// Table ID is encoded as follows:
//
// histogram title is a number that is to be interpreted
// as an uint64_t with the bits set as follows:
//
// bit 0: content type: 0 for , 1 for (see also bits 33 and 45)
// bit 1: polarization: T=1, L=0
// bit 2: t encoding: 0 for |t|, 1 for log(|t|)
// bit 3: W2 encoding: 0 for lin, 1 for log
// bit 4: Q2 encoding: 0 for lin, 1 for log
// bit 5-7: dipole model type (Enumerations.h)
// bit 8-15: mass number A
// bit 16-31: vector meson ID (PDG)
// bit 32: content encoding: 0 in lin, 1 in log
// bit 33: content type is lambda_ (bits 0 and 45 are zero in this case)
// bit 34-41: priority
// bit 42-44: dipole model parameter set (Enumerations.h)
// bit 45: content type is variance_A (bit 33 and 45 are zero in this case)
// bit 46: upc only table (tables with one Q2 bin only)
// bit 47: xpomeron encoding: 0 for lin, 1 for log (UPC only)
// bit 48: content type is lambda_skewedness (bits 0, 33, and 45 are zero in this case).
// bit 49-63: not used
//
// log implies ln (not log10)
// bit 0 is the LSB
//
// Internal histogram: x = Q2, y = W2, z = t (or the logs)
// x = xpomeron, y = t (or the logs) // UPC only
//
// Actual value is taken at the center of the bin.
// min/max functions (e.g. minQ2()) return the value of the first/last
// bin center.
//==============================================================================
#include "Table.h"
#include "TFile.h"
#include "TH2F.h"
#include "TH3F.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "GridSpline.h"
#define PR(x) cout << #x << " = " << (x) << endl;
Table::Table()
{
mID = 0;
m3DHisto = 0;
m2DHisto = 0;
mFillCounter = 0;
mBackupFrequence = 0;
}
Table& Table::operator=(const Table& tab)
{
if (this != &tab) {
delete m3DHisto;
delete m2DHisto;
mID = tab.mID;
if (tab.m3DHisto) {
m3DHisto = new TH3F(*(tab.m3DHisto));
m3DHisto->SetDirectory(0);
}
if (tab.m2DHisto) {
m2DHisto = new TH2F(*(tab.m2DHisto));
m2DHisto->SetDirectory(0);
}
mFilename = tab.mID;
mFillCounter = tab.mID;
mBackupFrequence = tab.mID;
mBackupPrefix = tab.mID;
mLastBackupFilename = tab.mID;
}
return *this;
}
Table::Table(const Table& tab)
{
mID = tab.mID;
if (tab.m3DHisto) {
m3DHisto = new TH3F(*(tab.m3DHisto));
m3DHisto->SetDirectory(0);
}
if (tab.m2DHisto) {
m2DHisto = new TH2F(*(tab.m2DHisto));
m2DHisto->SetDirectory(0);
}
mFilename = tab.mID;
mFillCounter = tab.mID;
mBackupFrequence = tab.mID;
mBackupPrefix = tab.mID;
mLastBackupFilename = tab.mID;
}
Table::~Table()
{
delete m2DHisto;
delete m3DHisto;
}
unsigned int Table::create(int nbinsQ2, double Q2min, double Q2max,
int nbinsW2, double W2min, double W2max,
int nbinsT, double tmin, double tmax,
bool logQ2, bool logW2, bool logt, bool logC,
AmplitudeMoment mom, GammaPolarization pol,
unsigned int A, int vm,
DipoleModelType model, DipoleModelParameterSet pset,
const char* filename, unsigned char priority)
{
if (m3DHisto != 0) {
cout << "Table:create(): cannot create, table already exists." << endl;
return 0;
}
if (nbinsQ2 < 2) {
cout << "Table:create(): number of bins in Q2 must be at least 2." << endl;
return 0;
}
if (nbinsW2 < 2) {
cout << "Table:create(): number of bins in W2 must be at least 2." << endl;
return 0;
}
if (nbinsT < 2) {
cout << "Table:create(): number of bins in t must be at least 2." << endl;
return 0;
}
if (filename) mFilename = string(filename); // name of file where table is written to
//
// Create table ID first
//
mID = (vm << 16);
mID |= (A << 8);
mID |= (model << 5);
mID |= (static_cast(pset) << 42);
if (mom == mean_A2) mID |= 1;
if (pol == transverse) mID |= (1 << 1);
if (logt) mID |= (1 << 2);
if (logW2) mID |= (1 << 3);
if (logQ2) mID |= (1 << 4);
if (logC) mID |= (static_cast(1) << 32);
if (mom == lambda_real) mID |= (static_cast(1) << 33);
if (mom == lambda_skew) mID |= (static_cast(1) << 48);
mID |= (static_cast(priority) << 34);
if (mom == variance_A) mID |= (static_cast(1) << 45);
ostringstream titlestream;
titlestream << mID;
string title = titlestream.str();
//
// Binning
//
// Interpolate() only works up to the bin center
// of the first and last bin. To guarantee that
// the full range is accessible we shift the min
// and max of each axis.
//
tmin = fabs(tmin);
tmax = fabs(tmax);
if (tmin > tmax) swap(tmin, tmax);
if (logQ2) Q2min = log(Q2min);
if (logQ2) Q2max = log(Q2max);
if (logW2) W2min = log(W2min);
if (logW2) W2max = log(W2max);
if (logt) tmin = log(tmin);
if (logt) tmax = log(tmax);
double binWidthQ2 = (Q2max - Q2min)/(nbinsQ2-1);
double binWidthW2 = (W2max - W2min)/(nbinsW2-1);
double binWidthT = (tmax - tmin)/(nbinsT-1);
//
// Book 3D histo to hold table
//
m3DHisto = new TH3F("table", title.c_str(),
nbinsQ2, Q2min-binWidthQ2/2., Q2max+binWidthQ2/2.,
nbinsW2, W2min-binWidthW2/2., W2max+binWidthW2/2.,
nbinsT, tmin-binWidthT/2., tmax+binWidthT/2.);
m3DHisto->SetDirectory(0);
return nbinsQ2*nbinsW2*nbinsT; // return total number of bins
}
//
// Overloaded function for UPC only
//
unsigned int Table::create(int nbinsX, double Xmin, double Xmax,
int nbinsT, double tmin, double tmax,
bool logX, bool logt, bool logC,
AmplitudeMoment mom, GammaPolarization pol,
unsigned int A, int vm,
DipoleModelType model, DipoleModelParameterSet pset,
const char* filename, unsigned char priority)
{
if (m2DHisto != 0) {
cout << "Table:create(): cannot create, table already exists." << endl;
return 0;
}
if (nbinsX < 2) {
cout << "Table:create(): number of bins in x must be at least 2." << endl;
return 0;
}
if (nbinsT < 2) {
cout << "Table:create(): number of bins in t must be at least 2." << endl;
return 0;
}
if (filename) mFilename = string(filename); // name of file where table is written to
//
// Create table ID first
//
mID = (vm << 16);
mID |= (A << 8);
mID |= (model << 5);
mID |= (static_cast(pset) << 42);
if (mom == mean_A2) mID |= 1;
if (pol == transverse) mID |= (1 << 1);
if (logt) mID |= (1 << 2);
if (logC) mID |= (static_cast(1) << 32);
if (mom == lambda_real) mID |= (static_cast(1) << 33);
if (mom == lambda_skew) mID |= (static_cast(1) << 48);
mID |= (static_cast(priority) << 34);
if (mom == variance_A) mID |= (static_cast(1) << 45);
mID |= (static_cast(1) << 46); // flag UPC
if (logX) mID |= (static_cast(1) << 47);
ostringstream titlestream;
titlestream << mID;
string title = titlestream.str();
//
// Binning
//
// Interpolate() only works up to the bin center
// of the first and last bin. To guarantee that
// the full range is accessible we shift the min
// and max of each axis.
//
tmin = fabs(tmin);
tmax = fabs(tmax);
if (tmin > tmax) swap(tmin, tmax);
if (logX) Xmin = log(Xmin);
if (logX) Xmax = log(Xmax);
if (logt) tmin = log(tmin);
if (logt) tmax = log(tmax);
double binWidthX = (Xmax - Xmin)/(nbinsX-1);
double binWidthT = (tmax - tmin)/(nbinsT-1);
//
// Book 2D histo to hold table
//
m2DHisto = new TH2F("table", title.c_str(),
nbinsX, Xmin-binWidthX/2., Xmax+binWidthX/2.,
nbinsT, tmin-binWidthT/2., tmax+binWidthT/2.);
m2DHisto->SetDirectory(0);
return nbinsX*nbinsT; // return total number of bins
}
int Table::numberOfEntries() const
{
if (m3DHisto) {
int nx = m3DHisto->GetXaxis()->GetNbins();
int ny = m3DHisto->GetYaxis()->GetNbins();
int nz = m3DHisto->GetZaxis()->GetNbins();
return nx*ny*nz;
}
else if (m2DHisto) { // UPC case
int nx = m2DHisto->GetXaxis()->GetNbins();
int ny = m2DHisto->GetYaxis()->GetNbins();
return nx*ny;
}
else
return 0;
}
void Table::binXYZ(int globalBin, int& binx, int& biny, int& binz) const
{
//
// Find correct bins for each axis.
// Don't use ROOT global bins here,
// they are a mess since they cross
// over underflow and overflow bins.
//
// All bins returned count starting
// at 1. The global bin starts at
// 0 as usual in C++.
//
int nx = m3DHisto->GetXaxis()->GetNbins();
int ny = m3DHisto->GetYaxis()->GetNbins();
binx = globalBin%nx;
biny = ((globalBin-binx)/nx)%ny;
binz = ((globalBin-binx)/nx -biny)/ny;
binx++; biny++; binz++;
}
//
// UPC version of binXYZ()
//
void Table::binXY(int globalBin, int& binx, int& biny) const
{
//
// Find correct bins for each axis.
// Don't use ROOT global bins here,
// they are a mess since they cross
// over underflow and overflow bins.
//
// All bins returned count starting
// at 1. The global bin starts at
// 0 as usual in C++.
//
int nx = m2DHisto->GetXaxis()->GetNbins();
binx = globalBin%nx;
biny = (globalBin-binx)/nx;
binx++; biny++;
}
void Table::binCenter(int i, double& Q2, double& W2, double& t) const
{
int binx, biny, binz;
binXYZ(i, binx, biny, binz);
// cout << i << '\t' << binx << '\t' << biny << '\t' << binz << endl;
double x = m3DHisto->GetXaxis()->GetBinCenter(binx);
double y = m3DHisto->GetYaxis()->GetBinCenter(biny);
double z = m3DHisto->GetZaxis()->GetBinCenter(binz);
Q2 = isLogQ2() ? exp(x) : x;
W2 = isLogW2() ? exp(y) : y;
t = isLogT() ? -exp(z) : -z;
if (t > 0) t = 0; // avoid numeric glitch
}
//
// UPC overload
//
void Table::binCenter(int i, double& xpom, double& t) const
{
int binx, biny;
binXY(i, binx, biny);
// cout << i << '\t' << binx << '\t' << biny << endl;
double x = m2DHisto->GetXaxis()->GetBinCenter(binx);
double y = m2DHisto->GetYaxis()->GetBinCenter(biny);
xpom = isLogX() ? exp(x) : x;
t = isLogT() ? -exp(y) : -y;
if (t > 0) t = 0; // avoid numeric glitch
}
void Table::fill(int i, double val, double err)
{
int binx, biny, binz;
if (isUPC())
binXY(i, binx, biny);
else
binXYZ(i, binx, biny, binz);
if (isLogContent()) {
if (val == 0) {
// cout << "Table::fill(): warning, log scale requested but value = 0." << endl;
val = numeric_limits::min();
}
val = log(fabs(val));
}
if (isUPC())
m2DHisto->SetBinContent(binx, biny, val);
else
m3DHisto->SetBinContent(binx, biny, binz, val);
- if(err!=0.) {
+ if (err!=0.) {
if (isUPC())
m2DHisto->SetBinError(binx, biny, err);
else
m3DHisto->SetBinError(binx, biny, binz, err);
}
mFillCounter++;
//
// Check if backup is due
//
if (mBackupFrequence)
if (mFillCounter%mBackupFrequence == 0) backup(i);
}
bool Table::fexists(const char* filename) const
{
ifstream ifs(filename);
return !ifs.fail();
}
void Table::write(const char* filename)
{
//
// 'filename' is optional argument. Null value implies
// that one was already passed via create(). Check
// that this is the case. If one is given, it is used
// no matter if it already was defined or not.
//
if (filename)
mFilename = string(filename);
else {
if (mFilename.empty()) { // should be defined but is not
cout << "Table::write(): Warning, no filename supplied. Will use 'sartre_table.root'." << endl;
mFilename = string("sartre_table.root");
}
}
//
// Check if file exist. If so alter
// the filename. Creating tables is
// a time consuming business so we do
// not want to lose precious data if
// something goes wrong here and we do
// want to prevent accidents.
//
while (fexists(mFilename.c_str())) {
mFilename += string(".save");
}
TFile hfile(mFilename.c_str(),"RECREATE");
if (isUPC())
m2DHisto->Write();
else
m3DHisto->Write();
hfile.Close();
cout << "Table::write(): table stored in file '" << mFilename.c_str()<< "'." << endl;
if (mBackupFrequence && !mLastBackupFilename.empty())
remove(mLastBackupFilename.c_str());
}
bool Table::read(const char* filename)
{
//
// Read histogram into memory
//
TFile *file = TFile::Open(filename,"READ");
if (file == 0) {
cout << "Table::read(): failed opening file '" << filename << "'." << endl;
return false;
}
//
// Clean up
//
delete m3DHisto; m3DHisto = 0;
delete m2DHisto; m2DHisto = 0;
//
// Read plain object first and check title which
// tells us if it is an UPC table or a std one.
// Then cast into the proper object.
//
auto ptr = file->Get("table");
if (ptr == 0) {
cout << "Table::read(): failed retrieving table from file '" << filename << "'." << endl;
return false;
}
mID = atoll(ptr->GetTitle()); // for isUPC() to work
if (isUPC()) {
m2DHisto = reinterpret_cast(ptr);
m2DHisto->SetDirectory(0);
}
else {
m3DHisto = reinterpret_cast(ptr);
m3DHisto->SetDirectory(0);
}
if (m2DHisto == 0 && m3DHisto == 0) {
cout << "Table::read(): failed retrieving table from file '" << filename << "'." << endl;
return false;
}
file->Close();
mFilename = string(filename);
return true;
}
int Table::vectorMesonId() const {return ((mID >> 16) & 0xFFFF);}
unsigned int Table::dipoleModelType() const {return ((mID >> 5) & 0x7);}
unsigned int Table::dipoleModelParameterSet() const {return ((mID >> 42) & 0x7);}
unsigned int Table::A() const {return ((mID >> 8) & 0xFF);}
bool Table::isLongitudinal() const {return !isTransverse();}
bool Table::isTransverse() const {return (mID & (1 << 1));}
GammaPolarization Table::polarization() const {return (isTransverse() ? transverse : longitudinal);}
AmplitudeMoment Table::moment() const {
if (isLambdaA())
return lambda_real;
else if (isMeanA())
return mean_A;
else if (isVarianceA())
return variance_A;
else if (isLambdaSkew())
return lambda_skew;
else
return mean_A2;
}
bool Table::isMeanA() const {
// bits 0, 33, 45 and 48 are not set
return !(mID & 1) && !(mID & (static_cast(1) << 33)) && !(mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isMeanA2() const {
// bit 0 is set, bits 33, 45 and 48 are not set
return (mID & 1) && !(mID & (static_cast(1) << 33)) && !(mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isLambdaA() const {
// bit 33 is set, bits 0 and 45 are not set
return !(mID & 1) && (mID & (static_cast(1) << 33)) && !(mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isLambdaSkew() const {
// bit 48 is set, bits 0, 33, and 45 are not set
return !(mID & 1) && (mID & (static_cast(1) << 48)) && !(mID & (static_cast(1) << 33))
&& !(mID & (static_cast(1) << 45));
}
bool Table::isVarianceA() const {
// bit 45 is set, bitss 0, 33, and 48 are not set
return !(mID & 1) && !(mID & (static_cast(1) << 33)) && (mID & (static_cast(1) << 45))
&& !(mID & (static_cast(1) << 48));
}
bool Table::isLogQ2() const {return (mID & (1 << 4));}
bool Table::isLogW2() const {return (mID & (1 << 3));}
bool Table::isLogT() const {return (mID & (1 << 2));}
bool Table::isLogX() const {return (mID & (static_cast(1) << 47));}
bool Table::isLogContent() const {return (mID & (static_cast(1) << 32));}
unsigned int Table::priority() const {return ((mID >> 34) & 0xFF);}
bool Table::isUPC() const {return (mID & (static_cast(1) << 46));}
uint64_t Table::id() const {return mID;}
double Table::binWidthQ2() const
{
if (!m3DHisto) return 0;
return m3DHisto->GetXaxis()->GetBinWidth(1);
}
double Table::binWidthW2() const
{
if (!m3DHisto) return 0;
return m3DHisto->GetYaxis()->GetBinWidth(1);
}
double Table::binWidthT() const
{
if (isUPC()) {
return m2DHisto->GetYaxis()->GetBinWidth(1);
}
else
return m3DHisto->GetZaxis()->GetBinWidth(1);
}
double Table::binWidthX() const
{
if (!m2DHisto) return 0;
return m2DHisto->GetXaxis()->GetBinWidth(1);
}
double Table::get(double Q2, double W2, double t) const
{
if (m3DHisto == 0) return 0;
//
// Transform variables to how they are stored in the table
//
double x = isLogQ2() ? log(Q2) : Q2;
double y = isLogW2() ? log(W2) : W2;
t = fabs(t);
double z = isLogT() ? log(t) : t;
//
// Tiny numerical glitches will cause TH3F::Interpolate() to fail
// since it requires that the variables lie between the first
// and last bin center, excluding the centers.
// Here we enforce that this is the case. The downside of this
// is that *all* values of Q2, W2, t will be forced to lie within.
// Hence the user has to make sure that the values are within
// the boundaries (see minQ2(), maxQ2(), minW2() etc.) before
// calling get().
// In this case the corrections below are tiny and are only
// applied to avoid minor glitches that do not affect the results.
//
TAxis *axis = m3DHisto->GetXaxis();
x = max(x, axis->GetBinCenter(1)+numeric_limits::epsilon());
x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
axis = m3DHisto->GetYaxis();
y = max(y, axis->GetBinCenter(1)+numeric_limits::epsilon());
y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
axis = m3DHisto->GetZaxis();
z = max(z, axis->GetBinCenter(1)+numeric_limits::epsilon());
z = min(z, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
// double result = InterpolateGridSpline(x, y, z); // tmp uncommented until 0's in tables are cleared
// double result = m3DHisto->Interpolate(x, y, z);
double result = modInterpolation(x, y, z);
if (result == 0 && isLogContent()) {
cout << "Table::get(): warning, 0 is not a valid table content when working in log scale." << endl;
}
if (isLogContent()) result = exp(result);
return result;
}
//
// UPC version
//
double Table::get(double xpom, double t) const
{
if (m2DHisto == 0) return 0;
//
// Transform variables to how they are stored in the table
//
double x = isLogX() ? log(xpom) : xpom;
t = fabs(t);
double y = isLogT() ? log(t) : t;
//
// Tiny numerical glitches will cause TH2F::Interpolate() to fail
// since it requires that the variables lie between the first
// and last bin center, excluding the centers.
// Here we enforce that this is the case. The downside of this
// is that *all* values of x, t will be forced to lie within.
// Hence the user has to make sure that the values are within
// the boundaries (see minX(), maxX(), minT() etc.) before
// calling get().
// In this case the corrections below are tiny and are only
// applied to avoid minor glitches that do not affect the results.
//
TAxis *axis = m2DHisto->GetXaxis();
x = max(x, axis->GetBinCenter(1)+numeric_limits::epsilon());
x = min(x, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
axis = m2DHisto->GetYaxis();
y = max(y, axis->GetBinCenter(1)+numeric_limits::epsilon());
y = min(y, axis->GetBinCenter(axis->GetNbins())-numeric_limits::epsilon());
double result = modInterpolation(x, y); // xpom, t
if (result == 0 && isLogContent()) {
cout << "Table::get(): warning, 0 is not a valid table content when working in log scale." << endl;
}
if (isLogContent()) result = exp(result);
return result;
}
double Table::minQ2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetXaxis();
double val = axis->GetBinCenter(1);
return isLogQ2() ? exp(val) : val;
}
double Table::maxQ2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetXaxis();
double val = axis->GetBinCenter(axis->GetNbins());
return isLogQ2() ? exp(val) : val;
}
double Table::minW2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetYaxis();
double val = axis->GetBinCenter(1);
return isLogW2() ? exp(val) : val;
}
double Table::maxW2() const
{
if (m3DHisto == 0) return 0;
TAxis *axis = m3DHisto->GetYaxis();
double val = axis->GetBinCenter(axis->GetNbins());
return isLogW2() ? exp(val) : val;
}
double Table::minT() const
{
if (m2DHisto == 0 && m3DHisto == 0) return 0;
TAxis *axis = 0;
if (isUPC())
axis = m2DHisto->GetYaxis();
else
axis = m3DHisto->GetZaxis();
double val = axis->GetBinCenter(axis->GetNbins()); // t always as |t| in table
return isLogT() ? -exp(val) : -val;
}
double Table::maxT() const
{
if (m2DHisto == 0 && m3DHisto == 0) return 0;
TAxis *axis = 0;
if (isUPC())
axis = m2DHisto->GetYaxis();
else
axis = m3DHisto->GetZaxis();
double val = axis->GetBinCenter(1); // t always as |t| in table
double result = isLogT() ? -exp(val) : -val;
if (result > 0) {
// cout << "Table::maxT(): warning, t > 0, t (" << result << ") set to 0 now." << endl;
result = 0;
}
return result;
}
double Table::minX() const
{
if (m2DHisto == 0) return 0;
TAxis *axis = m2DHisto->GetXaxis();
double val = axis->GetBinCenter(1);
return isLogX() ? exp(val) : val;
}
double Table::maxX() const
{
if (m2DHisto == 0) return 0;
TAxis *axis = m2DHisto->GetXaxis();
double val = axis->GetBinCenter(axis->GetNbins());
return isLogX() ? exp(val) : val;
}
string Table::filename() const {return mFilename;}
bool Table::writeToBinaryFile(const string& filename, bool verbose) const
{
//
// Open binary file
//
ofstream binaryFile;
binaryFile.open (filename, ios::out | ios::binary);
if (binaryFile.fail()) {
cout << "Table::writeToBinaryFile: failed to open file '" << filename << "'." << endl;
return false;
}
//
// Loop over all entries
//
double Q2, W2, xpom, t;
int binx, biny, binz;
double rawContent;
for (int i=0; iGetBinContent(binx, biny);
if (isLogContent()) {
rawContent = exp(rawContent);
}
binaryFile.write(reinterpret_cast(&i), sizeof(i));
binaryFile.write(reinterpret_cast(&xpom), sizeof(xpom));
binaryFile.write(reinterpret_cast(&t), sizeof(t));
binaryFile.write(reinterpret_cast(&rawContent), sizeof(rawContent));
if (verbose)
cout << "i=" << i << "\txpom=" << xpom << "\tt=" << t << "\tvalue=" << rawContent << endl;
}
else {
binXYZ(i, binx, biny, binz);
binCenter(i, Q2, W2, t);
rawContent = m3DHisto->GetBinContent(binx, biny, binz);
if (isLogContent()) {
rawContent = exp(rawContent);
}
binaryFile.write(reinterpret_cast(&i), sizeof(i));
binaryFile.write(reinterpret_cast(&Q2), sizeof(Q2));
binaryFile.write(reinterpret_cast(&W2), sizeof(W2));
binaryFile.write(reinterpret_cast(&t), sizeof(t));
binaryFile.write(reinterpret_cast(&rawContent), sizeof(rawContent));
if (verbose)
cout << "i=" << i << "\tQ2=" << Q2 << "\tW2=" << W2 << "\tt=" << t << "\tvalue=" << rawContent << endl;
}
}
binaryFile.close();
return true;
}
void Table::list(ostream& os, bool printContent, bool printStatistics, int startBin, int endBin) const
{
//
// List table info as compact as possible
//
ios::fmtflags fmt = os.flags(); // store current I/O flags
if ((isUPC() && !m2DHisto) || (!isUPC() && !m3DHisto)) {
os << "Table::list(): table is undefined." << endl;
return;
}
//
// First row: table id, A, and content
//
os << setw(11) << "Table ID = " << setw(16) << left << mID;
os << setw(7) << right << "A = " << setw(4) << left << A();
os << setw(20) << right << "content = " << (isLogContent() ? "log(" : "");
if (isMeanA())
os << "";
else if (isMeanA2())
os << "";
else if (isLambdaA())
os << "lambda_";
else if (isVarianceA())
os << "variance_A";
else if (isLambdaSkew())
os << "lambda_skew";
else
os << "unknown";
os << (isLogContent() ? ")" : "") << endl;
//
// Second row: vm and polarization
//
os << setw(34) << right << "vmId = " << setw(4) << left << vectorMesonId();
os << setw(20) << right << "polarization = " << (isTransverse() ? 'T' : 'L') << endl;
//
// Third row: dipole model and parameter set
//
os << setw(34) << right << "model = " << setw(7) << left;
if (dipoleModelType() == bSat)
os << "bSat ";
else if (dipoleModelType() == bNonSat)
os << "bNonSat";
else
os << "bCGC ";
os << setw(17) << right << "parameter set = ";
if (dipoleModelParameterSet() == KMW)
os << "KMW ";
else if (dipoleModelParameterSet() == HMPZ)
os << "HMPZ";
else if (dipoleModelParameterSet() == STU)
os << "STU";
else if (dipoleModelParameterSet() == CUSTOM)
os << "CUSTOM";
else
os << "? ";
os << endl;
//
// Third row: priority and UPC
//
os << setw(34) << right << "priority = " << setw(4) << left << priority();
os << setw(20) << right << "UPC = " << (isUPC() ? "yes" : "no") << endl;
//
// Next three rows: Q2, W2, t bins, range and log/linear
// or for UPC: x, t
//
if (isUPC()) {
os << setw(35) << right << "xp = [" << minX() << ", " << maxX()
<< "; bins = " << m2DHisto->GetXaxis()->GetNbins() << (isLogX() ? "; log]" : "; lin]") << endl;
os << setw(35) << right << "t = [" << minT() << ", " << maxT()
<< "; bins = " << m2DHisto->GetYaxis()->GetNbins() << (isLogT() ? "; log]" : "; lin]") << endl;
}
else {
os << setw(35) << right << "Q2 = [" << minQ2() << ", " << maxQ2()
<< "; bins = " << m3DHisto->GetXaxis()->GetNbins() << (isLogQ2() ? "; log]" : "; lin]") << endl;
os << setw(35) << right << "W2 = [" << minW2() << ", " << maxW2()
<< "; bins = " << m3DHisto->GetYaxis()->GetNbins() << (isLogW2() ? "; log]" : "; lin]") << endl;
os << setw(35) << right << "t = [" << minT() << ", " << maxT()
<< "; bins = " << m3DHisto->GetZaxis()->GetNbins() << (isLogT() ? "; log]" : "; lin]") << endl;
}
//
// Filename at the end
//
os << setw(34) << right << "file = " << mFilename.c_str() << endl;
os << endl;
//
// Print content (optional)
//
double Q2, W2, t, xpom;
int binx, biny, binz;
if (printContent) {
int thePrecision = os.precision();
double rawContent = 0;
- if(endBin==0) endBin=numberOfEntries();
+ if (endBin==0) endBin=numberOfEntries();
for (int i=startBin; iGetBinContent(binx, biny);
double error = m2DHisto->GetBinError(binx, biny);
os << "bin = " << setw(8) << left << i;
os << "xp = " << setw(13) << left << fixed << setprecision(5) << xpom;
os << "t = " << setw(16) << left << scientific << t;
os << "value = " << setw(15) << left << scientific << value;
os << "(binx = " << setw(5) << left << binx;
os << "biny = " << setw(5) << left << biny;
os << "content = " << setw(14) << left << scientific << rawContent;
os << "error = " << left << error << " ";
os << " peak = " << fixed << setprecision(3) << centerOverHighestNeighborRatio(binx, biny) << ')';
}
else {
binXYZ(i, binx, biny, binz);
binCenter(i, Q2, W2, t);
double value = get(Q2, W2, t);
rawContent = m3DHisto->GetBinContent(binx, biny, binz);
double error = m3DHisto->GetBinError(binx, biny, binz);
os << "bin = " << setw(8) << left << i;
os << "Q2 = " << setw(13) << left << fixed << setprecision(5) << Q2;
os << "W2 = " << setw(12) << left << fixed << W2;
os << "t = " << setw(16) << left << scientific << t;
os << "value = " << setw(15) << left << scientific << value;
os << "(binx = " << setw(5) << left << binx;
os << "biny = " << setw(5) << left << biny;
os << "binz = " << setw(5) << left << binz;
os << "content = " << setw(14) << left << scientific << rawContent;
os << "error = " << left << error << " ";
os << "peak = " << fixed << setprecision(3) << centerOverHighestNeighborRatio(binx, biny, binz) << ')';
}
if ( (isLogContent() && rawContent <= log(numeric_limits::min()*2)) || (rawContent == 0)) cout << " Z";
cout << endl;
}
os << endl;
os.precision(thePrecision);
}
//
// Print statistics (optional)
//
if (printStatistics) {
int nEmpty = 0;
int nInvalid = 0;
int nNegative = 0;
double sum = 0;
double maxContent = 0;
double minContent = numeric_limits::max();
int minBin, maxBin;
minBin = maxBin = 0;
double c = 0;
for (int i=0; iGetBinContent(binx, biny);
}
else {
binXYZ(i, binx, biny, binz);
c = m3DHisto->GetBinContent(binx, biny, binz);
}
if (c == 0) nEmpty++;
if ( !isLogContent() && c < 0) nNegative++;
if (std::isnan(c) || std::isinf(c) || !std::isfinite(c)) nInvalid++;
if (isLogContent()) c = exp(c);
if (c > maxContent) {maxContent = c; maxBin = i;}
if (c >= 0 && c < minContent) {minContent = c; minBin = i;}
if (c > 0) sum += c;
}
os << setw(34) << right << "total number of cells = " << numberOfEntries() << endl;
os << setw(34) << right << "cells with negative content = " << nNegative << endl;
os << setw(34) << right << "cells with no (0) content = " << nEmpty << endl;
os << setw(34) << right << "cells with invalid content = " << nInvalid << endl;
if (isUPC()) {
binCenter(minBin, xpom, t);
os << setw(34) << right << "minimum content = " << minContent << " at xp = " << xpom << ", t = " << t << endl;
binCenter(maxBin, xpom, t);
os << setw(34) << right << "maximum content = " << maxContent << " at xp = " << xpom << ", t = " << t << endl;
}
else {
binCenter(minBin, Q2, W2, t);
os << setw(34) << right << "minimum content = " << minContent << " at Q2 = " << Q2 << ", W2 = " << W2 << ", t = " << t << endl;
binCenter(maxBin, Q2, W2, t);
os << setw(34) << right << "maximum content = " << maxContent << " at Q2 = " << Q2 << ", W2 = " << W2 << ", t = " << t << endl;
}
os << setw(34) << right << "sum = " << sum << endl;
os << endl;
}
os.flags(fmt); // restore I/O flags
}
double Table::centerOverHighestNeighborRatio(int binx, int biny, int binz) const
{
double maxRatio = 0;
if (isUPC()) {
double centerValue = m2DHisto->GetBinContent(binx, biny);
int maxBinX = m2DHisto->GetXaxis()->GetNbins();
int maxBinY = m2DHisto->GetYaxis()->GetNbins();
for (int ix = -1; ix < 2; ix++) {
for (int iy = -1; iy < 2; iy++) {
int ixx = binx + ix;
int iyy = biny + iy;
if (ixx > 0 && ixx <= maxBinX && iyy > 0 && iyy <= maxBinY) {
double neighborValue = m2DHisto->GetBinContent(ixx, iyy);
double ratio = centerValue/neighborValue;
if (ratio > maxRatio) maxRatio = ratio;
}
}
}
}
else {
double centerValue = m3DHisto->GetBinContent(binx, biny, binz);
int maxBinX = m3DHisto->GetXaxis()->GetNbins();
int maxBinY = m3DHisto->GetYaxis()->GetNbins();
int maxBinZ = m3DHisto->GetZaxis()->GetNbins();
for (int ix = -1; ix < 2; ix++) {
for (int iy = -1; iy < 2; iy++) {
for (int iz = -1; iz < 2; iz++) {
int ixx = binx + ix;
int iyy = biny + iy;
int izz = binz + iz;
if (ixx > 0 && ixx <= maxBinX && iyy > 0 && iyy <= maxBinY && izz > 0 && izz <= maxBinZ) {
double neighborValue = m3DHisto->GetBinContent(ixx, iyy, izz);
double ratio = centerValue/neighborValue;
if (ratio > maxRatio) maxRatio = ratio;
}
}
}
}
}
return maxRatio;
}
double Table::InterpolateGridSpline(double x, double y, double z) const // Q2, W2, t
{
//
// The algorithm was taken from Cristian Constantin Lalescu.
// See http://arxiv.org/abs/0905.3564 for details.
//
// Grid points: 4, 6, or 8
// Spline order: 3,5 for order 4
// 3,5,7,9 for order 6
// 3,5,7,9,11,13 for order 8
//
//
// Find cell that refer to the position
//
int nx = m3DHisto->GetXaxis()->FindBin(x);
int ny = m3DHisto->GetYaxis()->FindBin(y);
int nz = m3DHisto->GetZaxis()->FindBin(z);
//
// Define the # of grid points depending on how far we are away
// from the edge. If at the edge we fall back to linear interpolation
// as provided by TH3.
// There must be a smarter way doing this than the code below
//
int gridPoints = 0;
if (nx-3 >= 1 && nx+4 <= m3DHisto->GetXaxis()->GetNbins() &&
ny-3 >= 1 && ny+4 <= m3DHisto->GetYaxis()->GetNbins() &&
nz-3 >= 1 && nz+4 <= m3DHisto->GetZaxis()->GetNbins()) {
gridPoints = 8;
}
else if (nx-2 >= 1 && nx+3 <= m3DHisto->GetXaxis()->GetNbins() &&
ny-2 >= 1 && ny+3 <= m3DHisto->GetYaxis()->GetNbins() &&
nz-2 >= 1 && nz+3 <= m3DHisto->GetZaxis()->GetNbins()) {
gridPoints = 6;
}
else if (nx-1 >= 1 && nx+2 <= m3DHisto->GetXaxis()->GetNbins() &&
ny-1 >= 1 && ny+2 <= m3DHisto->GetYaxis()->GetNbins() &&
nz-1 >= 1 && nz+2 <= m3DHisto->GetZaxis()->GetNbins()) {
gridPoints = 4;
}
else {
return m3DHisto->Interpolate(x,y,z);
}
//
// Find bin centers
//
double xc = m3DHisto->GetXaxis()->GetBinCenter(nx);
double yc = m3DHisto->GetYaxis()->GetBinCenter(ny);
double zc = m3DHisto->GetZaxis()->GetBinCenter(nz);
//
// Define the scale for coordinate transformations
// grid_spline expects x,y,z in bin units
//
double xscale = m3DHisto->GetXaxis()->GetBinWidth(1);
double yscale = m3DHisto->GetYaxis()->GetBinWidth(1);
double zscale = m3DHisto->GetZaxis()->GetBinWidth(1);
//
// Prepare grid spline alogorithm
//
double result, splineOrder;
if (gridPoints == 8) {
local_scal_3D<8> lf;
for (int i=-3;i<=4;i++) {
- for(int j=-3;j<=4;j++) {
+ for (int j=-3;j<=4;j++) {
lf(-3,i,j) = m3DHisto->GetBinContent(nx-3, ny+i, nz+j);
lf(-2,i,j) = m3DHisto->GetBinContent(nx-2, ny+i, nz+j);
lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
lf( 3,i,j) = m3DHisto->GetBinContent(nx+3, ny+i, nz+j);
lf( 4,i,j) = m3DHisto->GetBinContent(nx+4, ny+i, nz+j);
}
}
splineOrder = 7;
result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
}
else if (gridPoints == 6) {
local_scal_3D<6> lf;
for (int i=-2;i<=3;i++) {
- for(int j=-2;j<=3;j++) {
+ for (int j=-2;j<=3;j++) {
lf(-2,i,j) = m3DHisto->GetBinContent(nx-2, ny+i, nz+j);
lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
lf( 3,i,j) = m3DHisto->GetBinContent(nx+3, ny+i, nz+j);
}
}
splineOrder = 5;
result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
}
else if (gridPoints == 4) {
local_scal_3D<4> lf;
for (int i=-1;i<=2;i++) {
- for(int j=-1;j<=2;j++) {
+ for (int j=-1;j<=2;j++) {
lf(-1,i,j) = m3DHisto->GetBinContent(nx-1, ny+i, nz+j);
lf( 0,i,j) = m3DHisto->GetBinContent(nx , ny+i, nz+j);
lf( 1,i,j) = m3DHisto->GetBinContent(nx+1, ny+i, nz+j);
lf( 2,i,j) = m3DHisto->GetBinContent(nx+2, ny+i, nz+j);
}
}
splineOrder = 3;
result = grid_spline(splineOrder,lf, (x-xc)/xscale,(y-yc)/yscale,(z-zc)/zscale);
}
else {
cout << "Table::InterpolateGridSpline(): Error, illegal number of grid points." << endl;
result = 0;
}
return result;
}
double Table::modInterpolation(double x, double y, double z) const{ // Q2, W2, t
//
// After testing many points: The best results is had for content in log and rest in lin.
// However, if any of the bins that participate in the interpolation are zero
// the interpolation has to be linear.
//
// Make sure the point is within the histogram:
- if(m3DHisto->GetXaxis()->GetBinCenter(1) > x ||
+ if (m3DHisto->GetXaxis()->GetBinCenter(1) > x ||
m3DHisto->GetXaxis()->GetBinCenter(m3DHisto->GetXaxis()->GetNbins()) < x ||
m3DHisto->GetYaxis()->GetBinCenter(1) > y ||
m3DHisto->GetYaxis()->GetBinCenter(m3DHisto->GetYaxis()->GetNbins()) < y ||
m3DHisto->GetZaxis()->GetBinCenter(1) > z ||
m3DHisto->GetZaxis()->GetBinCenter(m3DHisto->GetZaxis()->GetNbins()) < z ){
cout<<"Table::myInterpolation Error: point lies outside the limits of the table!"<GetXaxis()->FindBin(x);
- if( x < m3DHisto->GetXaxis()->GetBinCenter(lbx) ) lbx--;
+ if ( x < m3DHisto->GetXaxis()->GetBinCenter(lbx) ) lbx--;
int ubx=lbx+1;
int lby = m3DHisto->GetYaxis()->FindBin(y);
- if( y < m3DHisto->GetYaxis()->GetBinCenter(lby) ) lby--;
+ if ( y < m3DHisto->GetYaxis()->GetBinCenter(lby) ) lby--;
int uby=lby+1;
int lbz = m3DHisto->GetZaxis()->FindBin(z);
- if( z < m3DHisto->GetZaxis()->GetBinCenter(lbz) ) lbz--;
+ if ( z < m3DHisto->GetZaxis()->GetBinCenter(lbz) ) lbz--;
int ubz=lbz+1;
//The corresponding bin-centers:
double ux=m3DHisto->GetXaxis()->GetBinCenter(ubx);
double lx=m3DHisto->GetXaxis()->GetBinCenter(lbx);
double uy=m3DHisto->GetYaxis()->GetBinCenter(uby);
double ly=m3DHisto->GetYaxis()->GetBinCenter(lby);
double uz=m3DHisto->GetZaxis()->GetBinCenter(ubz);
double lz=m3DHisto->GetZaxis()->GetBinCenter(lbz);
ux = isLogQ2() ? exp(ux) : ux;
lx = isLogQ2() ? exp(lx) : lx;
uy = isLogW2() ? exp(uy) : uy;
ly = isLogW2() ? exp(ly) : ly;
uz = isLogT() ? exp(uz) : uz;
lz = isLogT() ? exp(lz) : lz;
x = isLogQ2() ? exp(x) : x;
y = isLogW2() ? exp(y) : y;
z = isLogT() ? exp(z) : z;
//Find the distance to the point:
double xd = (x-lx)/(ux-lx);
double yd = (y-ly)/(uy-ly);
double zd = (z-lz)/(uz-lz);
// Make a vector containing all the values of the 8 surrounding bins
double v[8] = { m3DHisto->GetBinContent(lbx, lby, lbz), m3DHisto->GetBinContent(lbx, lby, ubz),
m3DHisto->GetBinContent(lbx, uby, lbz), m3DHisto->GetBinContent(lbx, uby, ubz),
m3DHisto->GetBinContent(ubx, lby, lbz), m3DHisto->GetBinContent(ubx, lby, ubz),
m3DHisto->GetBinContent(ubx, uby, lbz), m3DHisto->GetBinContent(ubx, uby, ubz) };
bool logC=true;
- for(int i=0; i<8; i++){
- if(isLogContent() and v[i] <= log(numeric_limits::min()*2) )
+ for (int i=0; i<8; i++){
+ if (isLogContent() && v[i] <= log(numeric_limits::min()*2) )
logC=false;
- if(!isLogContent() and v[i] == 0)
+ if (!isLogContent() && v[i] == 0)
logC=false;
}
// Make interpolation in log(content) except for when at least one of the points are zero.
- if(isLogContent() and !logC){
- for(int i=0; i<8; i++){
+ if (isLogContent() && !logC){
+ for (int i=0; i<8; i++){
if (v[i] <= log(numeric_limits::min()*2))
v[i] = 0;
else
v[i]=exp(v[i]);
}
}
- if(!isLogContent() and logC){
- for(int i=0; i<8; i++)
+ if (!isLogContent() && logC){
+ for (int i=0; i<8; i++)
v[i]=log(v[i]);
}
// First make four 1D interpolations in the z-direction
double i1 = v[0] * (1 - zd) + v[1] * zd;
double i2 = v[2] * (1 - zd) + v[3] * zd;
double j1 = v[4] * (1 - zd) + v[5] * zd;
double j2 = v[6] * (1 - zd) + v[7] * zd;
// Secondly, make two 1D interpolations in the y-direction using the values obtained.
double w1 = i1 * (1 - yd) + i2 * yd;
double w2 = j1 * (1 - yd) + j2 * yd;
// Finaly, make 1D interpolation in the x-direction
double result= w1 * (1-xd) + w2 * xd;
//Reverse exp/log compensation to get real result back:
- if(isLogContent() and !logC){
- if(result == 0)
+ if (isLogContent() && !logC){
+ if (result == 0)
result=log(numeric_limits::min());
else
result=log(result);
}
- if(!isLogContent() and logC)
+ if (!isLogContent() && logC)
result=exp(result);
return result;
}
//
// UPC version
//
double Table::modInterpolation(double x, double y) const{ // x, t
// Make sure the point is within the histogram:
- if(m2DHisto->GetXaxis()->GetBinCenter(1) > x ||
+ if (m2DHisto->GetXaxis()->GetBinCenter(1) > x ||
m2DHisto->GetXaxis()->GetBinCenter(m2DHisto->GetXaxis()->GetNbins()) < x ||
m2DHisto->GetYaxis()->GetBinCenter(1) > y ||
m2DHisto->GetYaxis()->GetBinCenter(m2DHisto->GetYaxis()->GetNbins()) < y) {
cout<<"Table::myInterpolation Error: point lies outside the limits of the table!"<GetXaxis()->FindBin(x);
int bin_y = m2DHisto->GetYaxis()->FindBin(y);
// Find quadrant of the bin we are in
int quadrant = 0;
double dx = m2DHisto->GetXaxis()->GetBinUpEdge(bin_x)-x;
double dy = m2DHisto->GetYaxis()->GetBinUpEdge(bin_y)-y;
if (dx<=m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy<=m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 1; // upper right
if (dx>m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy<=m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 2; // upper left
if (dx>m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy>m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 3; // lower left
if (dx<=m2DHisto->GetXaxis()->GetBinWidth(bin_x)/2 && dy>m2DHisto->GetYaxis()->GetBinWidth(bin_y)/2)
quadrant = 4; // lower right
double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
switch(quadrant) {
case 1:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x+1);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y+1);
break;
case 2:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x-1);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y+1);
break;
case 3:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x-1);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y-1);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
break;
case 4:
x1 = m2DHisto->GetXaxis()->GetBinCenter(bin_x);
y1 = m2DHisto->GetYaxis()->GetBinCenter(bin_y-1);
x2 = m2DHisto->GetXaxis()->GetBinCenter(bin_x+1);
y2 = m2DHisto->GetYaxis()->GetBinCenter(bin_y);
break;
}
// Now find the bins we interpolate with
int bin_x1 = m2DHisto->GetXaxis()->FindBin(x1);
if (bin_x1<1) bin_x1=1;
int bin_x2 = m2DHisto->GetXaxis()->FindBin(x2);
if (bin_x2>m2DHisto->GetXaxis()->GetNbins()) bin_x2=m2DHisto->GetXaxis()->GetNbins();
int bin_y1 = m2DHisto->GetYaxis()->FindBin(y1);
if (bin_y1<1) bin_y1=1;
int bin_y2 = m2DHisto->GetYaxis()->FindBin(y2);
if (bin_y2>m2DHisto->GetYaxis()->GetNbins()) bin_y2=m2DHisto->GetYaxis()->GetNbins();
// Get content
int bin_q22 = m2DHisto->GetBin(bin_x2,bin_y2);
int bin_q12 = m2DHisto->GetBin(bin_x1,bin_y2);
int bin_q11 = m2DHisto->GetBin(bin_x1,bin_y1);
int bin_q21 = m2DHisto->GetBin(bin_x2,bin_y1);
double q11 = m2DHisto->GetBinContent(bin_q11);
double q12 = m2DHisto->GetBinContent(bin_q12);
double q21 = m2DHisto->GetBinContent(bin_q21);
double q22 = m2DHisto->GetBinContent(bin_q22);
//
// As explained in the 3D version we interpolate on linear x, y
// but on log content, except when 0's are involved.
//
bool logC = true;
if ((isLogContent() && q11 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q11 <= 0)) logC = false;
if ((isLogContent() && q12 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q12 <= 0)) logC = false;
if ((isLogContent() && q21 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q21 <= 0)) logC = false;
if ((isLogContent() && q22 <= log(numeric_limits::min()*2)) ||
(!isLogContent() && q22 <= 0)) logC = false;
- if(isLogContent() && !logC) { // 0's present, use linear content
+ if (isLogContent() && !logC) { // 0's present, use linear content
if (q11 <= log(numeric_limits::min()*2))
q11 = 0;
else
q11=exp(q11);
if (q12 <= log(numeric_limits::min()*2))
q12 = 0;
else
q12=exp(q12);
if (q21 <= log(numeric_limits::min()*2))
q21 = 0;
else
q21=exp(q21);
if (q22 <= log(numeric_limits::min()*2))
q22 = 0;
else
q22=exp(q22);
}
if (!isLogContent() && logC) {
q11 = log(q11);
q12 = log(q12);
q21 = log(q21);
q22 = log(q22);
}
x = isLogX() ? exp(x) : x;
x1 = isLogX() ? exp(x1) : x1;
x2 = isLogX() ? exp(x2) : x2;
y = isLogT() ? exp(y) : y;
y1 = isLogT() ? exp(y1) : y1;
y2 = isLogT() ? exp(y2) : y2;
// Interpolation
double d = 1.0*(x2-x1)*(y2-y1);
double result = 1.0*q11/d*(x2-x)*(y2-y)+1.0*q21/d*(x-x1)*(y2-y)+1.0*q12/d*(x2-x)*(y-y1)+1.0*q22/d*(x-x1)*(y-y1);
// Reverse exp/log compensation to get real result back:
if (isLogContent() && !logC){
if (result == 0)
result=log(numeric_limits::min());
else
result=log(result);
}
if (!isLogContent() && logC)
result=exp(result);
return result;
}
void Table::setAutobackup(const char* prefix, int freq)
{
mBackupPrefix = string(prefix);
mBackupFrequence = freq;
}
void Table::backup(int backupBin)
{
ostringstream backupFilenameStream;
backupFilenameStream << mBackupPrefix.c_str() << "_backup." << static_cast(getpid())
<< '.' << mID << '.' << backupBin << ".root";
string filename = backupFilenameStream.str();
TFile file(filename.c_str(),"RECREATE");
m3DHisto->Write();
file.Close();
if (filename != mLastBackupFilename)
remove(mLastBackupFilename.c_str());
mLastBackupFilename = filename;
time_t now = time(0);
cout << "Table::backup(): autobackup performed, file = '" << mLastBackupFilename.c_str() << "', time = " << ctime(&now);
}
int Table::globalBin(int binx, int biny, int binz) const
{
int nbinx = m3DHisto->GetXaxis()->GetNbins();
int nbiny = m3DHisto->GetYaxis()->GetNbins();
return (binz-1)*nbinx*nbiny+(biny-1)*nbinx+(binx-1);
}
//
// UPC version
//
int Table::globalBin(int binx, int biny) const
{
int nbinx = m3DHisto->GetXaxis()->GetNbins();
return (biny-1)*nbinx+(binx-1);
}
Index: trunk/src/FrangibleNucleus.cpp
===================================================================
--- trunk/src/FrangibleNucleus.cpp (revision 432)
+++ trunk/src/FrangibleNucleus.cpp (revision 433)
@@ -1,204 +1,204 @@
//==============================================================================
// FrangibleNucleus.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "FrangibleNucleus.h"
#include "CNucleus.h"
#include "Constants.h"
#include "EventGeneratorSettings.h"
#include "TRandom3.h"
#define PR(x) cout << #x << " = " << (x) << endl;
FrangibleNucleus::FrangibleNucleus()
{
- mGeminiNucleus = 0;
+ mGeminiNucleus = nullptr;
mExcitationEnergy = 0;
}
FrangibleNucleus& FrangibleNucleus::operator=(const FrangibleNucleus& gn)
{
if (this != &gn) {
delete mGeminiNucleus;
mProducts.clear();
Nucleus::operator=(gn);
mExcitationEnergy = gn.mExcitationEnergy;
copy(gn.mProducts.begin(), gn.mProducts.end(), mProducts.begin());
if (gn.mGeminiNucleus)
mGeminiNucleus = new CNucleus(gn.mZ, gn.mA);
else
- mGeminiNucleus = 0;
+ mGeminiNucleus = nullptr;
}
return *this;
}
FrangibleNucleus::FrangibleNucleus(const FrangibleNucleus& gn) : Nucleus(gn)
{
mExcitationEnergy = gn.mExcitationEnergy;
copy(gn.mProducts.begin(), gn.mProducts.end(), mProducts.begin());
if (gn.mGeminiNucleus)
mGeminiNucleus = new CNucleus(gn.mZ, gn.mA);
else
mGeminiNucleus = 0;
}
void FrangibleNucleus::init(unsigned int A)
{
Nucleus::init(A);
if (mGeminiNucleus) delete mGeminiNucleus;
- mGeminiNucleus = 0;
+ mGeminiNucleus = nullptr;
}
void FrangibleNucleus::init(unsigned int A, bool enableBreakup)
{
Nucleus::init(A);
if (mGeminiNucleus) delete mGeminiNucleus;
//
// Init Gemini.
//
if (enableBreakup)
mGeminiNucleus = new CNucleus(mZ, mA);
else
- mGeminiNucleus = 0;
+ mGeminiNucleus = nullptr;
mExcitationEnergy = 0;
}
FrangibleNucleus::FrangibleNucleus(unsigned int A, bool enableBreakup) : Nucleus(A)
{
//
// Init Gemini.
// Only if requested, not needed otherwise.
//
if (enableBreakup)
mGeminiNucleus = new CNucleus(mZ, mA);
else
- mGeminiNucleus = 0;
+ mGeminiNucleus = nullptr;
mExcitationEnergy = 0;
}
FrangibleNucleus::~FrangibleNucleus()
{
delete mGeminiNucleus;
}
void FrangibleNucleus::resetBreakup()
{
mExcitationEnergy = 0;
mProducts.clear();
}
int FrangibleNucleus::breakup(const TLorentzVector& dissSystem)
{
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
//
// Estimate excitation energy
// Note that the dissSystem is given in units of 'per nucleon'.
// Hence we have to multiply the excitation energy with A.
//
mExcitationEnergy = (dissSystem.M()-protonMass)*mA;
double Ex = mExcitationEnergy;
double maxEx = settings->maxNuclearExcitationEnergy();
if (Ex > maxEx) {
if (settings->verboseLevel() > 1)
cout << "FrangibleNucleus::breakup(): Actual excitation energy (" << Ex
<< ") exceeded upper limit (" << maxEx << "). Reset to maximum value." << endl;
Ex = maxEx;
}
Ex *= 1000; // Gemini uses the total energy in MeV
//
// Setup excited nucleus
//
// Pass excitation energy and spin to Gemini
mGeminiNucleus->setCompoundNucleus(Ex,mSpin); //specify the excitation energy and spin
mGeminiNucleus->setVelocityCartesian(); // set initial velocity to zero (CMS)
// Set the direction of the spin vector (random)
TRandom3 *rndm = Settings::randomGenerator();
double phi = rndm->Uniform(2*M_PI);
double theta = acos(rndm->Uniform(-1, 1));
CAngle spin(theta, phi);
mGeminiNucleus->setSpinAxis(spin);
//
// Let the nucleus breakup
//
mGeminiNucleus->decay();
mProducts.clear();
if (mGeminiNucleus->abortEvent) {
cout << "Nucleus::breakup(): Error, decay aborted in Gemini++." << endl;
mGeminiNucleus->reset();
return 0;
}
int nfragments = mGeminiNucleus->getNumberOfProducts();
mProducts.resize(nfragments);
- for(int i=0; igetProducts(i); //set pointer to first stable product
mProducts[i].Z = product->iZ;
mProducts[i].A = product->iA;
mProducts[i].emissionTime = product->getTime();
mProducts[i].p = product->getLorentzVector();
mProducts[i].p.Boost(dissSystem.BoostVector());
mProducts[i].name = product->getName();
if (product->iZ == 1 && product->iA == 1)
mProducts[i].pdgId = 2212; // p
else if (product->iZ == 0 && product->iA == 1)
mProducts[i].pdgId = 2112; // n
else
mProducts[i].pdgId = pdgID(product->iZ, product->iA);
}
mGeminiNucleus->reset();
return nfragments;
}
const vector& FrangibleNucleus::breakupProducts() const
{
return mProducts;
}
void FrangibleNucleus::listBreakupProducts(ostream& os) const
{
TLorentzVector sys;
cout << "Excitation energy = " << mExcitationEnergy << endl;
cout << "Number of fragments = " << mProducts.size() << endl;
for (unsigned int i=0; i.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include
#include
#include
#include "Integrals.h"
#include "Constants.h"
#include "Nucleus.h"
#include "DipoleModel.h"
#include "AlphaStrong.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/Functor.h"
#include "TMath.h"
#include "WaveOverlap.h"
#include "Kinematics.h"
#include "TableGeneratorSettings.h"
#include "Enumerations.h"
#include "IntegrandWrappers.h"
#include "TF1.h"
#include "TH1F.h"
#include "cuba.h"
-#define PRf(x) printf(#x); printf("=%g \n", (x));
-#define PRi(x) printf(#x); printf("=%d \n", (x));
+#define PRs(x) cout << #x << " = " << scientific << (x) << endl;
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
Integrals::Integrals()
{
- mIsInitialized=false;
+ mIsInitialized = false;
mRelativePrecisionOfIntegration = 0;
- mWaveOverlap = 0;
- mDipoleModel = 0;
- mDipoleModelForSkewednessCorrection = 0;
+ mWaveOverlap = nullptr;
+ mDipoleModel = nullptr;
+ mDipoleModelForSkewednessCorrection = nullptr;
mIntegralImT = 0;
mIntegralImL = 0;
mIntegralReT = 0;
mIntegralReL = 0;
mErrorImT = 0;
mErrorImL = 0;
mErrorReT = 0;
mErrorReL = 0;
mProbImT = 0;
mProbImL = 0;
mProbReT = 0;
mProbReL = 0;
-
+
mIntegralTForSkewedness = 0;
mIntegralLForSkewedness = 0;
mErrorTForSkewedness = 0;
mErrorLForSkewedness = 0;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
- mVerbose=settings->verbose();
-
- int VMId=settings->vectorMesonId();
+ mVerbose = settings->verbose();
+
+ int VMId = settings->vectorMesonId();
mMV = settings->lookupPDG(VMId)->Mass();
mIsUPC = settings->UPC();
if (VMId==113 || VMId==333 || VMId == 443 || VMId == 553) {
mWaveOverlap = new WaveOverlapVM;
mWaveOverlap->setProcess(VMId);
mWaveOverlap->setWaveOverlapFunctionParameters(VMId);
- if(mVerbose)
- mWaveOverlap->testBoostedGaussianParameters(VMId);
+ if (mVerbose) mWaveOverlap->testBoostedGaussianParameters(VMId);
}
else if (VMId==22) {
mWaveOverlap = new WaveOverlapDVCS;
}
else {
cout << "Integrals::init(): Error, no exclusive production implemented for: "<< VMId << endl;
exit(1);
}
- DipoleModelType model=settings->dipoleModelType();
- if (model==bSat) {
+ DipoleModelType model = settings->dipoleModelType();
+ if (model == bSat) {
mDipoleModel = new DipoleModel_bSat;
}
- else if(model==bNonSat){
+ else if (model == bNonSat){
mDipoleModel = new DipoleModel_bNonSat;
}
- else if (model==bCGC) {
+ else if (model == bCGC) {
mDipoleModel = new DipoleModel_bCGC;
}
else {
cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
exit(1);
}
- mCalculateSkewedness=false;
- if(settings->A()==1 && settings->modesToCalculate()!=1 && settings->numberOfConfigurations()==1 && model==bSat) {
- mCalculateSkewedness=true;
+ mCalculateSkewedness = false;
+ if (settings->A()==1 && settings->modesToCalculate()!=1 && settings->numberOfConfigurations()==1 && model==bSat) {
+ mCalculateSkewedness = true;
mDipoleModelForSkewednessCorrection = new DipoleModel_bNonSat;
}
- mIsInitialized=true;
+ mIsInitialized = true;
}
Integrals::Integrals(const Integrals& integrals)
{
mIsInitialized = integrals.mIsInitialized;
mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS))
mWaveOverlap = new WaveOverlapDVCS;
else
mWaveOverlap = new WaveOverlapVM;
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat))
mDipoleModel = new DipoleModel_bSat;
- else if(typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
+ else if (typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
else
mDipoleModel = new DipoleModel_bCGC;
if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mErrorImT = integrals.mErrorImT;
mErrorImL = integrals.mErrorImL;
mErrorReT = integrals.mErrorReT;
mErrorReL = integrals.mErrorReL;
mProbImT = integrals.mProbImT;
mProbImL = integrals.mProbImL;
mProbReT = integrals.mProbReT;
mProbReL = integrals.mProbReL;
mMV = integrals.mMV;
}
Integrals& Integrals::operator=(const Integrals& integrals)
{
if (this != &integrals) {
delete mWaveOverlap;
delete mDipoleModel;
delete mDipoleModelForSkewednessCorrection;
if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS))
mWaveOverlap = new WaveOverlapDVCS;
else
mWaveOverlap = new WaveOverlapVM;
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat))
mDipoleModel = new DipoleModel_bSat;
else if (typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
else
mDipoleModel = new DipoleModel_bCGC;
if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
mDipoleModel = new DipoleModel_bNonSat;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mIntegralImT = integrals.mIntegralImT;
mIntegralImL = integrals.mIntegralImL;
mIntegralReT = integrals.mIntegralReT;
mIntegralReL = integrals.mIntegralReL;
mErrorImT = integrals.mErrorImT;
mErrorImL = integrals.mErrorImL;
mErrorReT = integrals.mErrorReT;
mErrorReL = integrals.mErrorReL;
mProbImT = integrals.mProbImT;
mProbImL = integrals.mProbImL;
mProbReT = integrals.mProbReT;
mProbReL = integrals.mProbReL;
mMV = integrals.mMV;
mIsInitialized = integrals.mIsInitialized;
mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
}
return *this;
}
Integrals::~Integrals()
{
delete mWaveOverlap;
delete mDipoleModel;
- if(mDipoleModelForSkewednessCorrection)
+ if (mDipoleModelForSkewednessCorrection)
delete mDipoleModelForSkewednessCorrection;
}
void Integrals::operator() (double t, double Q2, double W2)
{
- unsigned int A=dipoleModel()->nucleus()->A();
+ unsigned int A = dipoleModel()->nucleus()->A();
//make sure the configurations have been generated:
- if (!mDipoleModel->configurationExists() and A!=1) {
+ if (!mDipoleModel->configurationExists() && A!=1) {
// do not use cout
cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
exit(1);
}
if (setKinematicPoint(t, Q2, W2)) {
- if(A==1){
+ if (A==1){
calculateEp();
- if(mCalculateSkewedness){
+ if (mCalculateSkewedness){
calculateSkewedness();
}
}
else
calculate();
}
else {
fillZeroes();
}
}
void Integrals::operator() (double t, double xpom) //UPC
{
- unsigned int A=dipoleModel()->nucleus()->A();
+ unsigned int A = dipoleModel()->nucleus()->A();
//make sure the configurations have been generated:
- if (!mDipoleModel->configurationExists() and A!=1) {
+ if (!mDipoleModel->configurationExists() && A!=1) {
// do not use cout
cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
exit(1);
}
if (setKinematicPoint(t, xpom)) {
- if(A==1){
+ if (A==1){
calculateEp();
- if(mCalculateSkewedness){
+ if (mCalculateSkewedness){
calculateSkewedness();
}
}
else
calculate();
}
else {
fillZeroes();
}
}
void Integrals::fillZeroes(){
//Store the results
- mIntegralImT=0;
- mIntegralReT=0;
- mIntegralImL=0;
- mIntegralReL=0;
+ mIntegralImT = 0;
+ mIntegralReT = 0;
+ mIntegralImL = 0;
+ mIntegralReL = 0;
//Store the errors:
- mErrorImT=0;
- mErrorImL=0;
- mErrorReT=0;
- mErrorReL=0;
+ mErrorImT = 0;
+ mErrorImL = 0;
+ mErrorReT = 0;
+ mErrorReL = 0;
//Store the probabilities:
- mProbImT=0;
- mProbImL=0;
- mProbReT=0;
- mProbReL=0;
+ mProbImT = 0;
+ mProbImL = 0;
+ mProbReT = 0;
+ mProbReL = 0;
}
//*********EXCLUSIVE VECTOR MESONS OR DVCS: ********************************
void IntegralsExclusive::coherentIntegrals(double t, double Q2, double W2)
{
- if(setKinematicPoint(t, Q2, W2)){
+ if (setKinematicPoint(t, Q2, W2)){
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){
//store present kinematic point:
- double xprobe=kinematicPoint[3];
+ double xprobe = kinematicPoint[3];
dipoleModel()->createSigma_ep_LookupTable(xprobe);
}
calculateCoherent();
}
else
fillZeroes();
}
void IntegralsExclusive::coherentIntegrals(double t, double xpom)
{
- if(setKinematicPoint(t, xpom)){
+ if (setKinematicPoint(t, xpom)){
if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){
dipoleModel()->createSigma_ep_LookupTable(xpom);
}
calculateCoherent();
}
else
fillZeroes();
}
IntegralsExclusive::IntegralsExclusive()
{
}
IntegralsExclusive& IntegralsExclusive::operator=(const IntegralsExclusive& cobj)
{
if (this != &cobj) {
Integrals::operator=(cobj);
copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
}
return *this;
}
IntegralsExclusive::IntegralsExclusive(const IntegralsExclusive& cobj) : Integrals(cobj)
{
copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
}
bool IntegralsExclusive::setKinematicPoint(double t, double xpom) //UPC
{
bool result = true;
- kinematicPoint[0]=t;
- kinematicPoint[1]=0; //Q2
- kinematicPoint[2]=0; //W2 is not used
- kinematicPoint[3]=xpom;
+ kinematicPoint[0] = t;
+ kinematicPoint[1] = 0; //Q2
+ kinematicPoint[2] = 0; //W2 is not used
+ kinematicPoint[3] = xpom;
return result;
}
bool IntegralsExclusive::setKinematicPoint(double t, double Q2, double W2)
{
bool result = true;
- kinematicPoint[0]=t;
- kinematicPoint[1]=Q2;
- kinematicPoint[2]=W2;
- double xprobe=Kinematics::xpomeron(t, Q2, W2, mMV);
+ kinematicPoint[0] = t;
+ kinematicPoint[1] = Q2;
+ kinematicPoint[2] = W2;
+ double xprobe = Kinematics::xpomeron(t, Q2, W2, mMV);
if (xprobe<0 || xprobe>1)
result = false;
- kinematicPoint[3]=xprobe;
+ kinematicPoint[3] = xprobe;
return result;
}
void IntegralsExclusive::calculate()
{
//
// This function calls a wrapper from where the
// integral is calculated with the Cuhre method.
// Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
//
- const double epsrel=1.e-2, epsabs=1e-12;
+ const double epsrel = 1.e-2, epsabs = 1e-12;
const int flags=0, mineval=3e6, maxeval=1e9, key=0;
int nregionsTIm, nevalTIm, failTIm;
int nregionsTRe, nevalTRe, failTRe;
int nregionsLIm, nevalLIm, failLIm;
int nregionsLRe, nevalLRe, failLRe;
double valTIm=0, errTIm=0, probTIm=0;
double valLIm=0, errLIm=0, probLIm=0;
double valTRe=0, errTRe=0, probTRe=0;
double valLRe=0, errLRe=0, probLRe=0;
const char* statefile=0;
const int nvec=1;
// double probabilityCutOff=1e-6;
//
// Do the integrations
//
Cuhre(4, 1, integrandWrapperTIm, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsTIm, &nevalTIm, &failTIm, &valTIm, &errTIm, &probTIm);
- if(failTIm!=0 and mVerbose)
- printf("IntegralsExclusive::calculate(): Warning: Integration TIm did not reach desired precision! Error code=%d \n", failTIm);
+ if (failTIm != 0 && mVerbose) {
+ cout << "IntegralsExclusive::calculate(): Warning: Integration TIm did not reach desired precision! Error code=" << failTIm << endl;
+ }
//
// For UPC, calculate only transverse polarisation case
//
- if(!mIsUPC){
- Cuhre(4, 1, integrandWrapperLIm, this, nvec,
- epsrel, epsabs, flags,
- mineval, maxeval, key, statefile, 0, &nregionsLIm, &nevalLIm, &failLIm, &valLIm, &errLIm, &probLIm);
- if(failLIm!=0 and mVerbose)
- printf("IntegralsExclusive::calculate(): Warning: Integration LIm did not reach desired precision! Error code=%d \n", failLIm);
+ if (!mIsUPC) {
+ Cuhre(4, 1, integrandWrapperLIm, this, nvec,
+ epsrel, epsabs, flags,
+ mineval, maxeval, key, statefile, 0, &nregionsLIm, &nevalLIm, &failLIm, &valLIm, &errLIm, &probLIm);
+ if (failLIm!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculate(): Warning: Integration LIm did not reach desired precision! Error code=" << failLIm << endl;
+ }
}
Cuhre(4, 1, integrandWrapperTRe, this, nvec,
- epsrel, epsabs, flags,
- mineval, maxeval, key, statefile, 0, &nregionsTRe, &nevalTRe, &failTRe, &valTRe, &errTRe, &probTRe);
- if(failTRe!=0 and mVerbose)
- printf("IntegralsExclusive::calculate(): Warning: Integration TRe did not reach desired precision! Error code=%d \n", failTRe);
-
+ epsrel, epsabs, flags,
+ mineval, maxeval, key, statefile, 0, &nregionsTRe, &nevalTRe, &failTRe, &valTRe, &errTRe, &probTRe);
+ if (failTRe!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculate(): Warning: Integration TRe did not reach desired precision! Error code=" << failTRe << endl;
+ }
+
//
// For UPC, calculate only transverse polarisation case
//
- if(!mIsUPC){
- Cuhre(4, 1, integrandWrapperLRe, this, nvec,
- epsrel, epsabs, flags,
- mineval, maxeval, key, statefile, 0, &nregionsLRe, &nevalLRe, &failLRe, &valLRe, &errLRe, &probLRe);
- if(failLRe!=0 and mVerbose)
- printf("IntegralsExclusive::calculate(): Warning: Integration LRe did not reach desired precision! Error code=%d \n", failLRe);
+ if (!mIsUPC) {
+ Cuhre(4, 1, integrandWrapperLRe, this, nvec,
+ epsrel, epsabs, flags,
+ mineval, maxeval, key, statefile, 0, &nregionsLRe, &nevalLRe, &failLRe, &valLRe, &errLRe, &probLRe);
+ if (failLRe!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculate(): Warning: Integration LRe did not reach desired precision! Error code=" << failLRe << endl;
+ }
}
//
// Store the results:
//
- mIntegralImT=valTIm;
- mIntegralReT=valTRe;
- mIntegralImL=valLIm;
- mIntegralReL=valLRe;
+ mIntegralImT = valTIm;
+ mIntegralReT = valTRe;
+ mIntegralImL = valLIm;
+ mIntegralReL = valLRe;
//
// Store the errors:
//
- mErrorImT=errTIm;
- mErrorImL=errLIm;
- mErrorReT=errTRe;
- mErrorReL=errLRe;
+ mErrorImT = errTIm;
+ mErrorImL = errLIm;
+ mErrorReT = errTRe;
+ mErrorReL = errLRe;
//
// Store the probabilities:
//
- mProbImT=probTIm;
- mProbImL=probLIm;
- mProbReT=probTRe;
- mProbReL=probLRe;
+ mProbImT = probTIm;
+ mProbImL = probLIm;
+ mProbReT = probTRe;
+ mProbReL = probLRe;
}
void IntegralsExclusive::calculateEp()
{
//
// This function calls a wrapper from where the
// integral is calculated with the Cuhre method.
// Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
//
- const double epsrel=1.e-4, epsabs=1e-12;
+ const double epsrel = 1.e-4, epsabs = 1e-12;
const int flags=0, maxeval=1e9, key=0;
const int mineval=3e6;
int nregionsT, nevalT, failT;
int nregionsL, nevalL, failL;
double valT=0, errT=0, probT=0;
double valL=0, errL=0, probL=0;
const char* statefile=0;
const int nvec=1;
//
// Do the integrations
//
/*
bool bContinue=true;
bool isFirst=true;
while(bContinue){
double valTOld=valT;
Cuhre(4, 1, integrandWrapperTep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
- if(abs(valT-valTOld)/valT > epsrel){
+ if (abs(valT-valTOld)/valT > epsrel){
mineval*=3;
- if(isFirst)
+ if (isFirst)
mineval*=30;
isFirst=false;
- if(mineval>1e4)
+ if (mineval>1e4)
PR(mineval);
}
else
bContinue=false;
}
*/
Cuhre(4, 1, integrandWrapperTep, this, nvec,
- epsrel, epsabs, flags,
- mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
- if(failT!=0 and mVerbose)
- printf("IntegralsExclusive::calculateEp(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
+ epsrel, epsabs, flags,
+ mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
+ if (failT!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculateEp(): Warning: Integration T did not reach desired precision! Error code=" << failT << endl;
+ }
/*
- if(errT/valT > epsrel)
+ if (errT/valT > epsrel)
PR(errT/valT);
- if(errT < epsabs)
+ if (errT < epsabs)
PR(errT);
- if(probT>0.5)
+ if (probT>0.5)
PR(probT);
PR(nevalT);
PR(nregionsT);
*/
//
// For UPC, calculate only transverse polarisation case
//
- if(!mIsUPC){
+ if (!mIsUPC){
Cuhre(4, 1, integrandWrapperLep, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
- if(failL!=0 and mVerbose)
- printf("IntegralsExclusive::calculateEp(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
+ if (failL!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculateEp(): Warning: Integration L did not reach desired precision! Error code=" << failL << endl;
+ }
}
//
// Store the results:
//
mIntegralImT=valT;
mIntegralImL=valL;
//
// Store the errors:
//
mErrorImT=errT;
mErrorImL=errL;
//
// Store the probabilities:
//
mProbImT=probT;
mProbImL=probL;
}
void IntegralsExclusive::calculateSkewedness()
{
//
// This function calls a wrapper from where the
// integral is calculated with the Cuhre method.
// Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
//
const double epsrel=1.e-4, epsabs=1e-12;
const int flags=0, maxeval=1e9, key=0;
const int mineval=3e6;
int nregionsT, nevalT, failT;
int nregionsL, nevalL, failL;
double valT=0, errT=0, probT=0;
double valL=0, errL=0, probL=0;
const char* statefile=0;
const int nvec=1;
//
// Do the integrations
//
/*
bool bContinue=true;
bool isFirst=true;
while(bContinue){
double valTOld=valT;
Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
- if(abs(valT-valTOld)/valT > epsrel){
+ if (abs(valT-valTOld)/valT > epsrel){
mineval*=3;
- if(isFirst)
+ if (isFirst)
mineval*=30;
isFirst=false;
- if(mineval>1e4)
+ if (mineval>1e4)
PR(mineval);
}
else
bContinue=false;
}
*/
Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
- if(failT!=0 and mVerbose)
- printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
+ if (failT!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculateSkweedness(): Warning: Integration T did not reach desired precision! Error code" << failT << endl;
+ }
//
// For UPC, calculate only transverse polarisation case
//
- if(!mIsUPC){
+ if (!mIsUPC) {
Cuhre(4, 1, integrandWrapperLForSkewedness, this, nvec,
- epsrel, epsabs, flags,
- mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
- if(failL!=0 and mVerbose)
- printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
+ epsrel, epsabs, flags,
+ mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
+ if (failL!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculateSkweedness(): Warning: Integration L did not reach desired precision! Error code" << failL << endl;
+ }
}
//
// Store the results:
//
mIntegralTForSkewedness=valT;
mIntegralLForSkewedness=valL;
//
// Store the errors:
//
mErrorTForSkewedness=errT;
mErrorLForSkewedness=errL;
//
// Store the probabilities:
//
mProbImTForSkewedness=probT;
mProbImLForSkewedness=probL;
}
void IntegralsExclusive::calculateCoherent()
{
//
// As calculate() but for the coherent case
//
const double epsrel=1e-5, epsabs=1e-12;
const int flags=0, mineval=1e1, maxeval=1e9, key=0;
int nregionsT, nevalT, failT;
int nregionsL, nevalL, failL;
double valT, errT, probT;
double valL=0, errL=0, probL=0;
const int nvec=1;
const char* statefile=0;
//
// Do the integrations
//
Cuhre(3, 1, integrandWrapperCoherentAmplitudeT, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
- if(failT!=0 and mVerbose)
- printf("IntegralsExclusive::calculateCoherent(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
+ if (failT!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculateCoherent(): Warning: Integration T did not reach desired precision! Error code=" << failT << endl;
+ }
//
// For UPC, calculate only transverse polarisation case
//
- if(!mIsUPC){
+ if (!mIsUPC){
Cuhre(3, 1, integrandWrapperCoherentAmplitudeL, this, nvec,
epsrel, epsabs, flags,
mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
- if(failL!=0 and mVerbose)
- printf("IntegralsExclusive::calculateCoherent(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
+ if (failL!=0 && mVerbose) {
+ cout << "IntegralsExclusive::calculateCoherent(): Warning: Integration L did not reach desired precision! Error=" << failL << endl;
+ }
}
//
// Store the results
//
mIntegralImT=valT;
mIntegralImL=valL;
//
// Store the errors
//
mErrorImT=errT;
mErrorImL=errL;
}
//
// The following functions are the Integrands in the Amplitudes:
//
double IntegralsExclusive::uiAmplitudeTIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double cosArg = (b/hbarc)*Delta*cos(phi);
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2(r , b , phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*cos(cosArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiAmplitudeTRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double sinArg = b*Delta*cos(phi)/hbarc;
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*sin(sinArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiAmplitudeLIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double cosArg = b*Delta*cos(phi)/hbarc;
double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*cos(cosArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiAmplitudeLRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double sinArg = b*Delta*cos(phi)/hbarc;
double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
double BesselJ0 = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*
BesselJ0*b*sin(sinArg)*dsigdb2;
return result;
}
double IntegralsExclusive::uiCoherentAmplitudeT(double b, double z, double r, double Q2, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double xprobe=kinematicPoint[3];
double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe);
double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean;
return result;
}
double IntegralsExclusive::uiCoherentAmplitudeL(double b, double z, double r, double Q2, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double xprobe=kinematicPoint[3];
double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe);
double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean;
return result;
}
double IntegralsExclusive::uiAmplitudeTep(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
double IntegralsExclusive::uiAmplitudeLep(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
//
// Only for calculating the lamdba for Skewedness Corrections:
//
double IntegralsExclusive::uiAmplitudeTForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->T(z, Q2, r);
double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r , b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
double IntegralsExclusive::uiAmplitudeLForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)
{
double waveOverlap = mWaveOverlap->L(z, Q2, r);
double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r, b, xprobe);
double BesselJ0r = TMath::BesselJ0((1-z)*r*Delta/hbarc);
double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
result*=2*M_PI;
return result;
}
Index: trunk/src/ExclusiveFinalStateGenerator.cpp
===================================================================
--- trunk/src/ExclusiveFinalStateGenerator.cpp (revision 432)
+++ trunk/src/ExclusiveFinalStateGenerator.cpp (revision 433)
@@ -1,607 +1,607 @@
//==============================================================================
// ExclusiveFinalStateGenerator.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "ExclusiveFinalStateGenerator.h"
#include "EventGeneratorSettings.h"
#include "Kinematics.h"
#include "Event.h"
#include "Math/BrentRootFinder.h"
#include "Math/GSLRootFinder.h"
#include "Math/RootFinderAlgorithms.h"
#include "Math/IFunction.h"
#include
#include
#include
#include
#define PR(x) cout << #x << " = " << (x) << endl;
//-------------------------------------------------------------------------------
//
// Helper class needed to find root in
// ExclusiveFinalStateGenerator::generate()
//
//-------------------------------------------------------------------------------
class ScatteredProtonEnergyFormula : public ROOT::Math::IBaseFunctionOneDim
{
public:
double DoEval(double) const;
ROOT::Math::IBaseFunctionOneDim* Clone() const;
void calculateValidRange(double&, double&);
public:
double mT;
double mVmMass2;
double mMY2;
double mPhi; // azimuthal angle for scattered proton
TLorentzVector mProtonIn;
TLorentzVector mVirtualPhoton;
};
ROOT::Math::IBaseFunctionOneDim* ScatteredProtonEnergyFormula::Clone() const
{
return new ScatteredProtonEnergyFormula();
}
double ScatteredProtonEnergyFormula::DoEval(double Ep) const
{
double m2 = mProtonIn.M2();
double pzp = (mT - m2 - mMY2 + 2*mProtonIn.E() * Ep)/(2*mProtonIn.Pz());
double term = Ep*Ep-pzp*pzp-mMY2;
if (term < 0) return 99999; // out of kinematically allowed range
double ptp = sqrt(term);
TLorentzVector p_out(ptp*cos(mPhi), ptp*sin(mPhi), pzp, Ep);
double f = (mVirtualPhoton + mProtonIn - p_out)*(mVirtualPhoton + mProtonIn - p_out) - mVmMass2;
return f;
}
void ScatteredProtonEnergyFormula::calculateValidRange(double& lower, double& upper)
{
double m2 = mProtonIn.M2();
double term1 = mT-m2-mMY2;
double termA = mProtonIn.E()*term1;
double termB = sqrt(mProtonIn.Pz()*mProtonIn.Pz()*(term1*term1-4*m2*mMY2));
double termC = -2*m2;
lower = (termA+termB)/termC;
upper = (termA-termB)/termC;
if (lower > upper) swap(lower, upper);
lower += numeric_limits::epsilon();
upper -= numeric_limits::epsilon();
}
//-------------------------------------------------------------------------------
//
// Implementation of ExclusiveFinalStateGenerator
//
//-------------------------------------------------------------------------------
ExclusiveFinalStateGenerator::ExclusiveFinalStateGenerator() {/* no op */}
ExclusiveFinalStateGenerator::~ExclusiveFinalStateGenerator() {/* no op */}
bool ExclusiveFinalStateGenerator::generate(int id, double t, double y, double Q2,
bool isIncoherent, int A, Event *event)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
//
// The beam particles must be present in the event list
//
int ePos = -1;
int hPos = -1;
bool parentsOK = true;
if (event->particles.size() == 2) {
if (abs(event->particles[0].pdgId) == 11) {
ePos = 0;
hPos = 1;
}
else if (abs(event->particles[1].pdgId) == 11) {
ePos = 1;
hPos = 0;
}
else
parentsOK = false;
}
else
parentsOK = false;
if (!parentsOK) {
cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
return false;
}
//
// Store arguments locally
// (Some could also be obtained from the event structure)
//
mA = A;
mT = t;
if (mT > 0) mT = -mT; // ensure t<0
mQ2 = Q2;
mY = y;
mIsIncoherent = isIncoherent;
mElectronBeam = event->particles[ePos].p;
mHadronBeam = event->particles[hPos].p;
mMassVM = settings->lookupPDG(id)->Mass();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
//
// Constants
//
double const twopi = 2*M_PI;
double const hMass2 = mHadronBeam.M2();
//
// Incoherent diffarction
//
// Generate hadron dissociation mass according to
// dN/dM2 ~ 1/M2. Lower bound is of course the hadron
// mass and upper bound is some arbitrary value (for now).
//
// Note that we calculate and quote eA kinematics always in
// units of 'per nucleon'. Our model of incoherence is that the
// difference of the diffractive mass of one (1) proton out
// of the nucleus gives the final excitation energy E*.
// Hence we have to calculate E* and divide it by A to keep
// the kinematic consistent.
//
if (mIsIncoherent && mA > 1) {
const double lower = hMass2;
const double upper = 9; // GeV2
mMY2 = lower*upper/(upper - rndm->Uniform()*(upper-lower));
double MY_per_nucleon = (sqrt(hMass2)*(mA-1) + sqrt(mMY2))/mA;
mMY2 = MY_per_nucleon*MY_per_nucleon;
if (mMY2 < hMass2) mMY2 = hMass2;
}
else {
mMY2 = hMass2;
}
//
// Re-engineer scattered electron
//
// e'=(E', pt', pz') -> 3 unknowns
//
// Three equations:
// 1: me*me=E'*E'-pt'*pt'-pz'*pz'
// 2: Q2=-(e-e')^2=-2*me*me + 2*(E*E'-pz*pz')
// 3: W2=(P+e-e')^2=mp2+2*me2+2*(Ep*E-Pz*pz)-2*(Ep*E'-Pz*pz')-2*(E*E'-pz*pz')
//
- double Ee=mElectronBeam.E();
- double Pe=mElectronBeam.Pz();
- double Ep=mHadronBeam.E();
- double Pp=mHadronBeam.Pz();
- double W=event->W;
- double W2=W*W;
+ double Ee = mElectronBeam.E();
+ double Pe = mElectronBeam.Pz();
+ double Ep = mHadronBeam.E();
+ double Pp = mHadronBeam.Pz();
+ double W = event->W;
+ double W2 = W*W;
// Take masses from the beams in case they are not actually electrons or protons
- double me2=mElectronBeam.M2();
- double mp2=mHadronBeam.M2();
+ double me2 = mElectronBeam.M2();
+ double mp2 = mHadronBeam.M2();
//
// What we want for each particle:
//
double E, pz, pt, px, py, phi;
//
// Equations 2 and 3 yield:
//
E = Pe*(W2-mp2-2*Ee*Ep) + (Pp+Pe)*Q2 + 2*Pe*Pe*Pp + 2*me2*Pp;
E /= 2*(Ee*Pp-Ep*Pe);
pz = Ee*(W2-mp2) + (Ep+Ee)*Q2 + 2*Ee*Pe*Pp + 2*Ep*me2 - 2*Ee*Ee*Ep;
pz /= 2*(Ee*Pp-Ep*Pe);
//
// Equation 1:
//
pt = sqrt(E*E-pz*pz-me2);
phi = rndm->Uniform(twopi);
TLorentzVector theScatteredElectron(pt*sin(phi), pt*cos(phi), pz, E);
//
// Re-engineer virtual photon
//
// gamma=E-E'
- E=mElectronBeam.E()-theScatteredElectron.E();
- pz=mElectronBeam.Pz()-theScatteredElectron.Pz();
- px=mElectronBeam.Px()-theScatteredElectron.Px();
- py=mElectronBeam.Py()-theScatteredElectron.Py();
+ E = mElectronBeam.E()-theScatteredElectron.E();
+ pz = mElectronBeam.Pz()-theScatteredElectron.Pz();
+ px = mElectronBeam.Px()-theScatteredElectron.Px();
+ py = mElectronBeam.Py()-theScatteredElectron.Py();
TLorentzVector theVirtualPhoton = TLorentzVector(px, py, pz, E);
//
// Re-engineer scattered proton/dissociated proton
//
// No analytic solution. Need to run a root finder that does
// not need derivates but uses a bracketing algorithm (Brent).
// Correct brackets are crucial since ScatteredProtonEnergyFormula
// produces sqrt(-x) if outside the kinematically allowed range (it
// actually catches it and returns a large positive number, 0 doesn't
// work).
//
//
// Setup formula to solve root
//
phi = rndm->Uniform(twopi);
ScatteredProtonEnergyFormula formula;
formula.mT = mT;
formula.mVmMass2 = mMassVM*mMassVM;
formula.mPhi = phi;
formula.mProtonIn = mHadronBeam;
formula.mVirtualPhoton = theVirtualPhoton;
formula.mMY2 = mMY2;
//
// Find correct brackets to start with
//
double lower, upper;
formula.calculateValidRange(lower, upper);
if (upper > mHadronBeam.E() + theVirtualPhoton.E()) // limit excessive values
upper = mHadronBeam.E() + theVirtualPhoton.E(); // make it easier for Brent
//
// Run root finder
//
ROOT::Math::BrentRootFinder rootfinder;
rootfinder.SetFunction(formula, lower, upper);
rootfinder.Solve(10000, 0, 1.e-12);
E = rootfinder.Root();
if (/* rootfinder.Status() || */ fabs(formula(E)) > 1e-6) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, cannot find root. No final state defined." << endl;
return false;
}
//
// Outgoing proton (hadron) system
//
pz = (mT- hMass2 - mMY2 + 2*mHadronBeam.E()*E)/(2*mHadronBeam.Pz());
pt = sqrt(E*E-pz*pz-mMY2);
px = pt*cos(phi);
py = pt*sin(phi);
TLorentzVector theScatteredProton(px, py, pz, E);
//
// Finally the vector meson
//
TLorentzVector theVectorMeson((mHadronBeam + mElectronBeam) - (theScatteredElectron + theScatteredProton));
//
// Check for numerical glitches
//
if (!isValid(theScatteredElectron)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered electron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theScatteredProton)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered hadron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theVectorMeson)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, vector meson 4-vector is invalid." << endl;
return false;
}
//
// Add particles to event record
//
event->particles.resize(2+5);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int vm = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[vm].index = vm;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[vm].p = theVectorMeson;
event->particles[pomeron].p = theScatteredProton - mHadronBeam;
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[vm].pdgId = id;
event->particles[pomeron].pdgId = 990;
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[vm].status = 1;
event->particles[pomeron].status = 2;
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[vm].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(vm);
event->particles[gamma].daughters.push_back(pomeron);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
return true;
}
//
// UPC version
//
// Note: We call the beam particle that emits the photon "electron"
// even if its a proton/nucleus
//
bool ExclusiveFinalStateGenerator::generate(int id, double t, double xpom,
bool isIncoherent, int A, Event *event)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
//
// The beam particles must be present in the event list
//
int ePos = 0;
int hPos = 1;
bool parentsOK = true;
if (event->particles.size() != 2)
parentsOK = false;
if (!parentsOK) {
cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
return false;
}
//
// Store arguments locally
// (Some could also be obtained from the event structure)
//
mA = A;
mT = t;
if (mT > 0) mT = -mT; // ensure t<0
mIsIncoherent = isIncoherent;
mElectronBeam = event->particles[ePos].p;
mHadronBeam = event->particles[hPos].p;
mMassVM = settings->lookupPDG(id)->Mass();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
- mXp=xpom;
+ mXp = xpom;
//
// Constants
//
double const twopi = 2*M_PI;
double const hMass2 = mHadronBeam.M2();
//
// Incoherent diffraction
//
// Generate hadron dissociation mass according to
// dN/dM2 ~ 1/M2. Lower bound is of course the hadron
// mass and upper bound is some arbitrary value (for now).
//
// Note that we calculate and quote eA kinematics always in
// units of 'per nucleon'. Our model of incoherence is that the
// difference of the diffractive mass of one (1) proton out
// of the nucleus gives the final excitation energy E*.
// Hence we have to calculate E* and divide it by A to keep
// the kinematic consistent.
//
if (mIsIncoherent && mA > 1) {
const double lower = hMass2;
const double upper = 9; // GeV2
mMY2 = lower*upper/(upper - rndm->Uniform()*(upper-lower));
double MY_per_nucleon = (sqrt(hMass2)*(mA-1) + sqrt(mMY2))/mA;
mMY2 = MY_per_nucleon*MY_per_nucleon;
if (mMY2 < hMass2) mMY2 = hMass2;
}
else {
mMY2 = hMass2;
}
//
// Check for smallest allowed xpom
//
- double xpom_min=Kinematics::xpomMin(mMassVM, mT, mHadronBeam, mElectronBeam, mMY2-hMass2);
- if(xpom < xpom_min){
+ double xpom_min = Kinematics::xpomMin(mMassVM, mT, mHadronBeam, mElectronBeam, mMY2-hMass2);
+ if (xpom < xpom_min){
if (settings->verboseLevel() > 2) cout<<"xpom = "<verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, scattered hadron 4-vector is invalid." << endl;
return false;
}
if (!isValid(theVectorMeson)) {
if (settings->verboseLevel() > 2) cout << "ExclusiveFinalStateGenerator::generate(): error, vector meson 4-vector is invalid." << endl;
return false;
}
//
// Add particles to event record
//
event->particles.resize(2+5);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int vm = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[vm].index = vm;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[vm].p = theVectorMeson;
event->particles[pomeron].p = thePomeron;
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[vm].pdgId = id;
event->particles[pomeron].pdgId = 990;
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[vm].status = 1;
event->particles[pomeron].status = 2;
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[vm].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(vm);
event->particles[gamma].daughters.push_back(pomeron);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
//fill event structure
- double y=mHadronBeam*theVirtualPhoton/(mHadronBeam*mElectronBeam);
- double W2=(theVirtualPhoton+mHadronBeam).M2();
- mQ2=-theVirtualPhoton.M2();
-
- event->Q2=mQ2;
- event->W=sqrt(W2);
- event->y=y;
+ double y = mHadronBeam*theVirtualPhoton/(mHadronBeam*mElectronBeam);
+ double W2 = (theVirtualPhoton+mHadronBeam).M2();
+ mQ2 = -theVirtualPhoton.M2();
+
+ event->Q2 = mQ2;
+ event->W = sqrt(W2);
+ event->y = y;
return true;
}
Index: trunk/src/WaveOverlap.cpp
===================================================================
--- trunk/src/WaveOverlap.cpp (revision 432)
+++ trunk/src/WaveOverlap.cpp (revision 433)
@@ -1,273 +1,272 @@
//==============================================================================
// WaveOverlap.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: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include "WaveOverlap.h"
#include "Constants.h"
#include "TableGeneratorSettings.h"
#include "DipoleModelParameters.h"
#include "TMath.h"
#include
#include
#include "TF1.h"
#include "Math/Functor.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
WaveOverlap::WaveOverlap()
{
mParameters = new DipoleModelParameters(TableGeneratorSettings::instance());
}
WaveOverlap::~WaveOverlap() {/* no op*/}
//
// VECTOR MESONS
//
WaveOverlapVM::WaveOverlapVM()
{
mNT = mRT2 = 0;
mMf = 0;
mMf2 = 0;
mEf = 0;
mMV = 0;
mNL = mRL2 = 0;
mBoostedGaussianMf = 0;
}
double WaveOverlapVM::uiNormT(const double* var) const{
double z=var[0];
double r=var[1];
double phi=transverseWaveFunction(r, z); //GeV0
double dphidr=dDrTransverseWaveFunction(r, z); //GeV1
return 3.*r/hbarc/(z*z*(1-z)*(1-z))
*(mMf2*phi*phi + (z*z + (1-z)*(1-z))*dphidr*dphidr); //GeV1
}
double WaveOverlapVM::uiNormL(const double* var) const{
- double z=var[0];
- double r=var[1];
- double phi=longitudinalWaveFunction(r, z); //GeV0
- double d2phidr2=laplaceRLongitudinalWaveFunction(r, z); //GeV2
- double term=mMV*phi + (mMf2*phi - d2phidr2)/(mMV*z*(1-z)); //GeV1
+ double z = var[0];
+ double r = var[1];
+ double phi = longitudinalWaveFunction(r, z); //GeV0
+ double d2phidr2 = laplaceRLongitudinalWaveFunction(r, z); //GeV2
+ double term = mMV*phi + (mMf2*phi - d2phidr2)/(mMV*z*(1-z)); //GeV1
return 3.*r/hbarc*term*term; //GeV1
}
double WaveOverlapVM::uiDecayWidth(double* var, double*) const
{
- double z=*var;
- double phi=longitudinalWaveFunction(0., z); //GeV0
- double d2phidr2=laplaceRLongitudinalWaveFunction(0., z); //GeV0
+ double z = *var;
+ double phi = longitudinalWaveFunction(0., z); //GeV0
+ double d2phidr2 = laplaceRLongitudinalWaveFunction(0., z); //GeV0
double result = mEf*3./M_PI*(mMV*phi+(mMf2*phi-d2phidr2)/(mMV*z*(1-z))); //GeV1
return result;
}
void WaveOverlapVM::testBoostedGaussianParameters(int id) const
{
//
// This function calculates the resulting normalisation
// and decay width resulting from the boosted gaussian parameters
// and compare with the actual values. This can be used to
// test or modify the boosted gaussian parameters.
//
// KMW paper hep-ph/0606272 Eqs. (24)-(28)
//
// Start with decay width:
TF1 fDW("fDW", this, &WaveOverlapVM::uiDecayWidth, 0., 1., 0.);
ROOT::Math::WrappedTF1 wfDW(fDW);
ROOT::Math::GaussIntegrator giDW;
giDW.SetFunction(wfDW);
giDW.SetAbsTolerance(0.);
giDW.SetRelTolerance(1e-5);
double f_VL=giDW.Integral(0., 1.); //GeV
PR(f_VL);
cout<<"The e+e- decay width is: "<ee)="<<4*M_PI*alpha_em*alpha_em*f_VL*f_VL/3/mMV*1e6
<<" keV"<ee)=1.340 +- 0.018 keV"<ee)=5.55 +- 0.14 +- 0.02 keV"<ee)=1.27 +- 0.04 keV"<ee)=7.04 +- 0.06 keV"<boostedGaussianR2(val);
mRT2 = mRL2;
mNT = mParameters->boostedGaussianNT(val);
mNL = mParameters->boostedGaussianNL(val);
mBoostedGaussianMf = mParameters->boostedGaussianQuarkMass(val);
mBoostedGaussianMf2 = mBoostedGaussianMf*mBoostedGaussianMf;
}
void WaveOverlapVM::setProcess(int val)
{
switch (val) {
case 113:
mMf = mParameters->quarkMass(1);
mMV = 0.776;
mEf = 1./sqrt(2.);
break;
case 333:
mMf = mParameters->quarkMass(2);
mMV = 1.019;
mEf = 1./3.;
break;
case 443:
mMf = mParameters->quarkMass(3);
mMV = 3.096916;
mEf = 2./3.;
break;
case 553:
mMf = mParameters->quarkMass(4);
mMV = 9.46;
mEf = 1./3.;
break;
default:
cerr << "WaveOverlap::setProcess(): error no such type: " << val << endl;
break;
}
mMf2 = mMf*mMf;
}
//
// DVCS
//
double WaveOverlapDVCS::T(double z, double Q2, double r) const
{
// KMW paper hep-ph/0606272 Eq. 17
double term0, term1, term2;
double result=0;
for (int iFlav=0; iFlav<4; iFlav++) {
double mf=mParameters->quarkMass(iFlav);
double ef=quarkCharge[iFlav];
double eps2 = z*(1-z)*Q2 + mf*mf;
double eps = sqrt(eps2);
- term0=2.*Nc/M_PI*alpha_em*ef*ef;
- term1=( z*z+(1-z)*(1-z) )*
- eps*TMath::BesselK1(eps*r/hbarc)*mf*TMath::BesselK1(mf*r/hbarc);
- term2=mf*mf*TMath::BesselK0(eps*r/hbarc)*TMath::BesselK0(mf*r/hbarc);
+ term0 = 2.*Nc/M_PI*alpha_em*ef*ef;
+ term1 = (z*z+(1-z)*(1-z))*eps*TMath::BesselK1(eps*r/hbarc)*mf*TMath::BesselK1(mf*r/hbarc);
+ term2 = mf*mf*TMath::BesselK0(eps*r/hbarc)*TMath::BesselK0(mf*r/hbarc);
result += term0*(term1+term2);
}
return result;
}
double WaveOverlapDVCS::L(double, double, double) const {return 0;}
Index: trunk/src/DglapEvolution.cpp
===================================================================
--- trunk/src/DglapEvolution.cpp (revision 432)
+++ trunk/src/DglapEvolution.cpp (revision 433)
@@ -1,524 +1,524 @@
//==============================================================================
// DglapEvolution.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
// $Date$
// $Author$
//==============================================================================
#include "DglapEvolution.h"
#include "TableGeneratorSettings.h"
#include "DipoleModelParameters.h"
#include
#include
#include
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
DglapEvolution* DglapEvolution::mInstance = 0;
DglapEvolution::DglapEvolution()
{
DipoleModelParameters parameters(TableGeneratorSettings::instance());
//
// For speed purposes the key parameters are held as data member.
// Here we assign the proper ones depending on the model and the
// parameters set choosen.
//
mAg = parameters.Ag();
mLambdaG = parameters.lambdaG();
mMu02 = parameters.mu02();
//
// Define alpha_s functor
//
mAlphaStrong = new AlphaStrong;
//
// Set alpha_s at c, b, t mass
//
mMC = mAlphaStrong->massCharm();
mMB = mAlphaStrong->massBottom();
mMT = mAlphaStrong->massTop();
const double FourPi = 4.*M_PI;
mALPC = mAlphaStrong->at(mMC*mMC)/FourPi;
mALPB = mAlphaStrong->at(mMB*mMB)/FourPi;
mALPT = mAlphaStrong->at(mMT*mMT)/FourPi;
mALPS = mAlphaStrong->at(mMu02)/FourPi;
//
// Initialization of support points in n-space and weights for the
// Gauss quadrature and of the anomalous dimensions for the RG
// evolution at these n-values.
//
//
// Weights and support points for nomalized 8 point gauss quadrature
//
double wz[9] = {0, 0.101228536290376,0.222381034453374,0.313706645877887,
0.362683783378362,0.362683783378362,0.313706645877887,
0.222381034453374,0.101228536290376};
double zs[9] = {0, -0.960289856497536,-0.796666477413627,-0.525532409916329,
-0.183434642495650,0.183434642495650,0.525532409916329,
0.796666477413627,0.960289856497536};
//
// Integration contour parameters
//
double down[18] = {0, 0., 0.5, 1., 2., 3., 4., 6., 8.,
1.e1, 1.2e1, 1.5e1, 1.8e1, 2.1e1, 2.4e1, 2.7e1, 3.e1, 3.3e1};
double up[18];
mC = 0.8;
double phi = M_PI * 3./4.;
double co = cos(phi);
double si = sin(phi);
mCC = complex(co, si);
for (int i=1; i <=16; i++) up[i] = down[i+1];
up[17] = 36.;
//
// Support points and weights for the gauss integration
// (the factor (up-down)/2 is included in the weights)
//
int k = 0;
double sum, diff, z;
for (int i=1; i <=17; i++) {
sum = up[i] + down[i];
diff = up[i] - down[i];
for (int j=1; j <=8; j++) {
k++;
z = 0.5 * (sum + diff * zs[j]);
mWN[k] = diff / 2.* wz[j];
mN[k] = complex(mC+co*z+1.,si*z);
}
}
anom();
//
// Defaults for lookup table
//
mTableMinX = 1.e-10;
mTableMaxX = 1;
mTableMinQ2 = 1.;
mTableMaxQ2 = 1e6;
mLookupTableIsFilled = false;
mUseLookupTable = false;
mNumberOfNodesInX = mNumberOfNodesInQ2 = 0;
}
DglapEvolution::~DglapEvolution()
{
if (mAlphaStrong) delete mAlphaStrong;
if (mLookupTableIsFilled) {
for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
delete [] mLookupTable[i];
delete [] mLookupTable;
}
}
double DglapEvolution::alphaSxG(double x, double Q2)
{
double val = xG(x, Q2);
return val*mAlphaStrong->at(Q2);
}
double DglapEvolution::xG(double x, double Q2)
{
double result = 0;
if (!mUseLookupTable) {
result = xG_Engine(x, Q2);
}
else if (mUseLookupTable && !mLookupTableIsFilled) {
cout << "DglapEvolution::xG(): Warning, use of lookup table requested\n"
<< " but table is not setup. First use \n"
<< " generateLookupTable() to fill table.\n"
<< " Will fall back to numeric calculation."
<< endl;
result = xG_Engine(x, Q2);
}
else
result = xG_Interpolator(x, Q2);
return result;
}
double DglapEvolution::xG_Engine(double x, double Q2)
{
//
// These are provided by AlphaStrong ensuring
// consistency between evolution and alpha_s.
//
if (mAlphaStrong->order() != 0) {
cout << "DglapEvolution::xG_Engine(): Fatal error, alpha_s is in order "
<< mAlphaStrong->order()
<< " but DglapEvolution only support order = 0. Stop here." << endl;
exit(1);
}
double alpq = mAlphaStrong->at(Q2)/(4*M_PI);
//
// Q2 and x dependent quantities
//
double ax = log(x);
double ex = exp(-mC * ax);
//
// integration length parameter for the mellin inversion
//
const int nmax = 136;
//
// Gluon density and output
//
complex fn[137];
reno(fn, alpq, nmax, mAg, mLambdaG);
long double fun = 0; // long double is needed
long double fz;
complex xnm,cex;
for (int i=1; i <= nmax; i++) {
xnm = (mC - mN[i]+1.) * ax;
cex = exp(xnm) / M_PI * mCC;
fz = imag(fn[i]*cex);
fun = fun + mWN[i] * fz;
}
double pa = fun * ex;
return pa;
}
void DglapEvolution::generateLookupTable(int nx, int nq2)
{
//
// Delete old lookup table (if it exists)
//
if (mLookupTableIsFilled) {
for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
delete [] mLookupTable[i];
delete [] mLookupTable;
}
mNumberOfNodesInX = nx;
mNumberOfNodesInQ2 = nq2;
//
// Create new table
//
mLookupTable = new double*[mNumberOfNodesInQ2];
- for(unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
+ for (unsigned int i = 0; i < mNumberOfNodesInQ2; ++i)
mLookupTable[i] = new double[mNumberOfNodesInX];
//
// Fill lookup table
//
int printEach = mNumberOfNodesInQ2*mNumberOfNodesInX/10;
int nCount = 0;
cout << "DglapEvolution: generating lookup table "; cout.flush();
for (unsigned int i = 0; i < mNumberOfNodesInQ2; i++) {
double Q2 = pow(mTableMaxQ2/mTableMinQ2,i*1./(mNumberOfNodesInQ2-1.))*mTableMinQ2;
for (unsigned int j = 0; j < mNumberOfNodesInX; j++) {
double x = pow(mTableMaxX/mTableMinX,j*1./(mNumberOfNodesInX-1.))*mTableMinX;
mLookupTable[i][j] = xG_Engine(x, Q2);
nCount++;
if (nCount%printEach == 0) cout << '.'; cout.flush();
}
}
cout << " done" << endl;
mLookupTableIsFilled = true;
}
double DglapEvolution::luovi(double f[4], double arg[4], double z)
{
double cof[4];
for (int i=0; i < 4; i++ ) cof[i]=f[i];
for (int i=0; i < 3 ; i++) {
for (int j=i; j < 3 ; j++) {
int jndex = 2 - j;
int index = jndex + 1 + i;
cof[index]= (cof[index]-cof[index-1])/(arg[index]-arg[jndex]);
}
}
- double sum=cof[3];
+ double sum = cof[3];
for (int i=0; i < 3; i++ ) {
int index = 2 - i;
sum = (z-arg[index])*sum + cof[index];
}
return sum;
}
double DglapEvolution::xG_Interpolator(double x, double Q2)
{
int q2steps = mNumberOfNodesInQ2-1;
int xsteps = mNumberOfNodesInX-1;
//
// Position in the Q2 grid
//
double realq = q2steps * log(Q2/mTableMinQ2)/log(mTableMaxQ2/mTableMinQ2);
int n_q2 = static_cast(realq);
if (n_q2 <= 0) {n_q2 = 1;}
if (n_q2 > q2steps-2) {n_q2 = n_q2-2;}
//
// Position in the x grid
//
double realx = xsteps * log(x/mTableMinX)/log(mTableMaxX/mTableMinX);
int n_x = static_cast(realx);
if (n_x <= 0) { n_x = 1;}
if (n_x > xsteps-2){ n_x = n_x-2;}
//
// Starting the interpolation
//
double arg[4];
for (int i=0; i < 4; i++) { arg[i] = n_x-1+i;}
double fu[4], fg[4];
for (int i=0; i < 4; i++) {
fu[0] = mLookupTable[n_q2-1+i][n_x-1];
fu[1] = mLookupTable[n_q2-1+i][n_x];
fu[2] = mLookupTable[n_q2-1+i][n_x+1];
fu[3] = mLookupTable[n_q2-1+i][n_x+2];
fg[i] = luovi(fu,arg,realx);
}
for (int i=0; i < 4; i++) { arg[i] = n_q2-1+i;}
return luovi(fg, arg, realq);
}
void DglapEvolution::reno(complex *fn, double alpq, int nmax, double ag, double lambdag)
{
//
// Mellin-n space Q**2 - evolution of the gluon at LO
//
// The moments are calculated on an array of moments, mN, suitable
// for a (non-adaptive) Gauss quadrature.
//
// Currently this takes the simplest possible fit form:
// xg = A_g x^(-lambdag) (1-x)^(5.6), following Amir&Raju
//
for (int k1 = 1; k1 <= nmax; k1++) {
//
// Input moments of the parton densities
// at the low scale
//
complex xn = mN[k1];
complex gln = ag * (1.0 / (xn + 5.0 - lambdag)
- 6.0 / (xn + 4.0 - lambdag)
+ 15.0 / (xn + 3.0 - lambdag)
- 20.0 / (xn + 2.0 -lambdag)
+ 15.0 / (xn + 1.0 - lambdag)
- 6.0 / (xn - lambdag)
+ 1.0 / (xn - lambdag - 1.0));
int f;
double xl, s, alp;
complex ep, gl;
if (alpq >= mALPC) { // evolution below the charm threshold
f = 3;
xl = mALPS / alpq;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
}
else if ((alpq < mALPC) && (alpq >= mALPB)) { // between thresholds
f = 3;
xl = mALPS / mALPC;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
f = 4;
xl = mALPC / alpq;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
}
else if (alpq < mALPB) { // above bottom threshold
f = 3;
xl = mALPS / mALPC;
s = log (xl);
ep = exp (- mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
f = 4;
alp = mALPB;
xl = mALPC / mALPB;
s = log (xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
f = 5;
xl = mALPB / alpq;
s = log(xl);
ep = exp(-mAP[k1][f]*s);
gl = gln;
gln = ep * gl;
}
fn[k1] = gln;
}
}
void DglapEvolution::anom() {
//
// Anomalous dimensions for leading order evolution of parton densities.
// The moments are calculated on an externally given array of mellin
// moments, mN, suitable for a (non-adaptive) quadrature.
//
// Present version: the number of active flavours in the factorization
// is fixed to ff=3, in the beta-function it varies between f=3 and
// f=5. The dimension of the moment array is 136.
//
double b0, b02;
complex ggi, ggf;
complex xn, xap;
complex gg;
for (int k1=1; k1 <= 136; k1++) {
xn = mN[k1];
anCalc(ggi, ggf, xn);
for (int k2=3; k2 <= 5; k2++) {
double f = k2;
// anomalous dimensions and related quantities in leading order
b0 = 11.- 2./3.* f;
b02 = 2.* b0;
gg = ggi + f * ggf;
xap = gg / b02;
mAP[k1][k2] = xap;
}
}
}
void DglapEvolution::anCalc(complex& ggi,
complex& ggf, complex& xn)
{
complex xn1 = xn + 1.;
complex xn2 = xn + 2.;
complex xnm = xn - 1.;
//
// Leading order
//
complex cpsi = psiFunction(xn1) + 0.577216;
ggi = -22.- 24./(xn * xnm) - 24./(xn1 * xn2) + 24.* cpsi;
ggf = 4./3.;
}
complex DglapEvolution::psiFunction(complex z)
{
//
// psi - function for complex argument
//
complex sub;
while (real(z) < 10) {
sub = sub - 1./z;
z += 1.;
}
complex rz = 1./z;
complex dz = rz * rz;
complex result = sub + log(z) - rz/2.- dz/2520. * ( 210.+ dz * (-21.+10.*dz ));
return result;
}
double DglapEvolution::logDerivative(double x, double Q2)
{
double alpq = mAlphaStrong->at(Q2)/(4*M_PI);
//
// Q2 and x dependent quantities
//
double ax = log(x);
double ex = exp(-mC * ax);
//
// integration length parameter for the mellin inversion
//
const int nmax = 136;
//
// Gluon density and output
//
complex fn[137];
reno(fn, alpq, nmax, mAg, mLambdaG);
long double fun = 0; // long double is needed
long double fz;
complex xnm,cex;
for (int i=1; i <= nmax; i++) {
xnm = (mC - mN[i]+1.) * ax;
cex = exp(xnm) / M_PI * mCC;
fz = imag(fn[i]*cex);
fun = fun + mWN[i] * fz;
}
double glue = fun * ex;
fun = 0;
for (int i=1; i <= nmax; i++) {
xnm = (mC - mN[i]+1.) * ax;
cex = exp(xnm) / M_PI * mCC;
fz = imag(fn[i]*cex*mN[i]);
fun = fun + mWN[i] * fz;
}
double pa = fun * ex;
double lambda = -(1-pa/glue);
return lambda;
}
void DglapEvolution::list(ostream& os)
{
os << "DglapEvolution parameters:" << endl;
os << " Ag = " << mAg << endl;
os << " lambdaG = " << mLambdaG << endl;
os << " mu02 = " << mMu02 << endl;
os << " m_c = " << mMC << " (from AlphaStrong)" << endl;
os << " m_b = " << mMB << " (from AlphaStrong)" << endl;
os << " m_t = " << mMT << " (from AlphaStrong)" << endl;
if (mUseLookupTable)
os << " Lookup table in use" << endl;
else
os << " Lookup table not used" << endl;
}