No OneTemporary
Actions

Size

107 KB

Subscribers

None

View Options

	Index: trunk/src/SartreInclusiveDiffraction.cpp
	===================================================================
	--- trunk/src/SartreInclusiveDiffraction.cpp (revision 539)
	+++ trunk/src/SartreInclusiveDiffraction.cpp (revision 540)
	@@ -1,861 +1,865 @@
	//==============================================================================
	// SartreInclusiveDiffraction.cpp
	//
	// Copyright (C) 2024 Tobias Toll and Thomas Ullrich
	//
	// This file is part of Sartre.
	//
	// This program is free software: you can redistribute it and/or modify
	// it under the terms of the GNU General Public License as published by
	// the Free Software Foundation.
	// This program is distributed in the hope that it will be useful,
	// but without any warranty; without even the implied warranty of
	// merchantability or fitness for a particular purpose. See the
	// GNU General Public License for more details.
	// You should have received a copy of the GNU General Public License
	// along with this program. If not, see <http://www.gnu.org/licenses/>.
	//
	// Author: Tobias Toll
	// Last update:
	// $Date: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
	// $Author: ullrich $
	//==============================================================================
	//
	// Note:
	// When not using runcards, the user must first create an instance of
	// SartreInclusiveDiffraction and then get the settings via one of:
	// SartreInclusiveDiffraction::runSettings()
	// EventGeneratorSettings::instance()
	// Once init() is called settings cannot be changed any more.
	//==============================================================================
	#include "Version.h"
	#include "Kinematics.h"
	#include "Constants.h"
	#include "SartreInclusiveDiffraction.h"
	#include "Math/Functor.h"
	#include "ModeFinderFunctor.h"
	#include "Math/BrentMinimizer1D.h"
	#include "Math/IntegratorMultiDim.h"
	#include "Math/GSLMinimizer.h"
	#include "TUnuranMultiContDist.h"
	#include <string>
	#include <limits>
	#include <cmath>
	#include <iomanip>
	#include "TH2D.h"
	#include "TFile.h"
	#include "TF3.h"
	#include "InclusiveDiffractiveCrossSectionsFromTables.h"

	using namespace std;

	#define PR(x) cout << #x << " = " << (x) << endl;

	SartreInclusiveDiffraction::SartreInclusiveDiffraction()
	{
	mSettings = EventGeneratorSettings::instance();
	mRandom = mSettings->randomGenerator();
	mIsInitialized = false;
	mCurrentEvent = 0;
	mNucleus = 0;
	mUpcNucleus= 0;
	// mSettings = 0;
	mPDF_Functor = 0;
	mPDF = 0;
	mEventCounter = 0;
	mTriesCounter = 0;
	mTotalCrossSection = 0;
	mTotalCrossSection_qqg = 0;
	mTotalCrossSection_qq = 0;
	mCrossSection = 0;
	mUnuran = 0;
	mEvents = 0;
	mTries = 0;
	mS = 0;
	mA = 0;
	- mUseCrossSectionTables=true;
	}

	SartreInclusiveDiffraction::~SartreInclusiveDiffraction()
	{
	delete mNucleus;
	delete mUpcNucleus;
	delete mPDF_Functor;
	delete mPDF;
	delete mCrossSection;
	delete mUnuran;
	delete mCurrentEvent;

	}

	bool SartreInclusiveDiffraction::init(const char* runcard)
	{
	mStartTime = time(0);
	bool ok;

	//
	// Reset member variables.
	// Note that one instance of Sartre should be able to get
	// initialized multiple times.
	//
	mEvents = 0;
	mTries = 0;
	mTotalCrossSection = 0;
	mTotalCrossSection_qqg = 0;
	mTotalCrossSection_qq = 0;

	//
	// Print header
	//
	string ctstr(ctime(&mStartTime));
	ctstr.erase(ctstr.size()-1, 1);
	cout << "/========================================================================\\" << endl;
	cout << "\| \|" << endl;
	cout << "\| Sartre, Version " << setw(54) << left << VERSION << right << '\|' << endl;
	cout << "\| \|" << endl;
	cout << "\| An event generator for inclusive diffraction \|" << endl;
	cout << "\| in ep and eA collisions based on the dipole model. \|" << endl;
	cout << "\| \|" << endl;
	cout << "\| Copyright (C) 2010-2024 Tobias Toll and Thomas Ullrich \|" << endl;
	cout << "\| \|" << endl;
	cout << "\| This program is free software: you can redistribute it and/or modify \|" << endl;
	cout << "\| it under the terms of the GNU General Public License as published by \|" << endl;
	cout << "\| the Free Software Foundation, either version 3 of the License, or \|" << endl;
	cout << "\| any later version. \|" << endl;
	cout << "\| \|" << endl;
	cout << "\| Code compiled on " << setw(12) << left << __DATE__;
	cout << setw(41) << left << __TIME__ << right << '\|' << endl;
	cout << "\| Run started at " << setw(55) << left << ctstr.c_str() << right << '\|' << endl;
	cout << "\\========================================================================/" << endl;

	mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton
	mSettings->setInclusiveDiffractionMode(true);

	//
	// Read runcard if available
	//
	if (runcard) {
	if (!mSettings->readSettingsFromFile(runcard)) {
	cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl;
	exit(1);
	}
	else
	cout << "Runcard is '" << runcard << "'." << endl;
	}
	else
	cout << "No runcard provided." << endl;

	+ mUseCrossSectionTables=mSettings->useInclusiveCrossSectionTables();
	+
	//
	// Set beam particles and center of mass energy
	//
	mElectronBeam = mSettings->eBeam();
	mHadronBeam = mSettings->hBeam();
	mS = Kinematics::s(mElectronBeam, mHadronBeam);
	mA = mSettings->A();

	bool allowBreakup = mSettings->enableNuclearBreakup();
	if (mA == 1) allowBreakup = false;

	if (allowBreakup) {
	if (!getenv("SARTRE_DIR")) {
	cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n"
	"to locate tables needed for the generation if nuclear breakups." << endl;
	exit(1);
	}
	}
	if (mNucleus) delete mNucleus;
	mNucleus = new FrangibleNucleus(mA, allowBreakup);

	string upcNucleusName;
	cout << "Sartre is running in inclusive diffraction mode" << endl;
	cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
	cout << "Hadron beam: " << mHadronBeam << endl;
	cout << "Electron beam: " << mElectronBeam << endl;

	//
	// Get details about the processes and models
	//
	mDipoleModelType = mSettings->dipoleModelType();
	mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
	cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
	cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
	cout << "Process is ";
	if (mA > 1)
	cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + X"<< endl;
	else
	cout << "e + p -> e' + p' + X"<< endl;

	//
	// Print-out seed for reference
	//
	cout << "Random generator seed: " << mSettings->seed() << endl;

	//
	// Load in the tables containing the amplitude moments
	//
	if (!getenv("SARTRE_DIR")) {
	cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
	exit(1);
	}

	//
	// Kinematic limits and generator range
	//
	// There are 3 ranges we have to deal with
	// 1. the kinematic range requested by the user
	// if given.
	// The user can only control Q2 and W but not t.
	// For UPC that's xpom.
	// 2. the range of the table(s)
	// 3. the kinematically allowed range
	//
	// Of course (3) is more complex than a simple cube/square.
	// However, we deal with the detailed shape of the kinematic
	// range later using Kinematics::valid() when we generate the
	// individual events.
	// For setting up UNU.RAN we have to get the cubic/square
	// envelope that satifies (1)-(3).
	// Note, that they are correlated which makes the order
	// in which we do things a bit tricky.
	//


	//
	// Step 2:
	// Kinematic limits might overrule boundaries from step 1
	//
	// let:
	// 0: beta
	// 1: Q2
	// 2: W2
	// 3: z
	//
	// Setup cross-section functor
	// It is this functor that is used by all other functors,
	// functions, and wrappers when dealing with cross-sections.
	//
	+
	if (mCrossSection) delete mCrossSection;
	if(!mUseCrossSectionTables){
	mCrossSection = new InclusiveDiffractiveCrossSections;
	}
	else{
	mCrossSection = new InclusiveDiffractionCrossSectionsFromTables;
	}
	//
	// Here we need to put in Q2 > mu02.
	// The calculation goes haywire for Q2 becoming too small
	//
	double mu02 = mCrossSection->dipoleModel()->getParameters()->mu02();
	// double mf = Settings::quarkMass(0);
	double mf=tt_quarkMass[0];

	//Now for the diffractive mass:
	double kineMX2min = 4mfmf;//max(mu02, 4mfmf);
	// This gives the limits for W2 and for inelasticity y:
	double kineW2min = kineMX2min+protonMass2;
	double kineYmin = Kinematics::ymin(mS, kineMX2min);
	// ...and y gives the limit for Q2, unless it's smaller than the cut-off mu02:
	double kineQ2min = Kinematics::Q2min(kineYmin);
	kineQ2min = max(kineQ2min, mu02);
	//The maximum momenta are given by s:
	double kineQ2max = Kinematics::Q2max(mS);
	double kineW2max = mS-kineQ2min-protonMass2;

	//
	// Now for all the fractions, making sure they make sense
	//
	double kineBetamax = kineQ2max/(kineQ2max+kineMX2min);
	double kineBetamin = Kinematics::betamin(kineQ2min, mS, sqrt(kineW2min), 0);
	double kineMX2max = (1.-kineBetamin)/kineBetamin*kineQ2max;
	kineMX2max=min(kineMX2max, kineW2max-protonMass2);
	double kineZmin=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	double kineZmax=1-kineZmin;
	mLowerLimit[0] = kineBetamin;
	mUpperLimit[0] = kineBetamax;
	mLowerLimit[1] = kineQ2min;
	mUpperLimit[1] = kineQ2max;
	mLowerLimit[2] = kineW2min;
	mUpperLimit[2] = kineW2max;
	mLowerLimit[3] = kineZmin;
	mUpperLimit[3] = kineZmax;

	// Step 3:
	// Deal with user provided limits.
	// User settings are ignored (switched off) if min >= max.
	// Only allow user limits for Q2, W2, and beta,
	// not z.
	//
	bool W2orQ2changed=false;
	if (mSettings->Wmin() < mSettings->Wmax()) { // W2 first

	if (mSettings->W2min() < mLowerLimit[2]) {
	cout << "Warning, requested lower limit of W (" << mSettings->Wmin() << ") "
	<< "is smaller than limit given by lookup tables and/or kinematic range (" << sqrt(mLowerLimit[2]) << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mLowerLimit[2] = mSettings->W2min();
	W2orQ2changed=true;
	}

	if (mSettings->W2max() > mUpperLimit[2]) {
	cout << "Warning, requested upper limit of W (" << mSettings->Wmax() << ") "
	<< "exceeds limit given by lookup tables and/or kinematic range (" << sqrt(mUpperLimit[2]) << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mUpperLimit[2] = mSettings->W2max();
	W2orQ2changed=true;
	}
	}
	if (mSettings->Q2min() < mSettings->Q2max()) { // Q2

	if (mSettings->Q2min() < mLowerLimit[1]) {
	cout << "Warning, requested lower limit of Q2 (" << mSettings->Q2min() << ") "
	<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mLowerLimit[1] = mSettings->Q2min();
	W2orQ2changed=true;
	}

	if (mSettings->Q2max() > mUpperLimit[1]) {
	cout << "Warning, requested upper limit of Q2 (" << mSettings->Q2max() << ") "
	<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mUpperLimit[1] = mSettings->Q2max();
	W2orQ2changed=true;
	}
	}
	if (W2orQ2changed){//mLowerLimit beta, Q2, W2, z
	//kineBetamax = kineQ2max/(kineQ2max+kineMX2min);
	mUpperLimit[0] = mUpperLimit[1]/(mUpperLimit[1]+kineMX2min);

	//kineBetamin = Kinematics::betamin(kineQ2min, mS, sqrt(kineW2min), 0);
	mLowerLimit[0] = Kinematics::betamin(mLowerLimit[1], mS, sqrt(mLowerLimit[2]), 0);

	// kineMX2max = (1.-kineBetamin)/kineBetamin*kineQ2max;
	// kineMX2max=min(kineMX2max, kineW2max-protonMass2);
	kineMX2max = (1.-mLowerLimit[0])/mLowerLimit[0]*mUpperLimit[1];
	kineMX2max=min(kineMX2max, mUpperLimit[2]-protonMass2);

	// kineZmin=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	// kineZmax=1-kineZmin;
	mLowerLimit[3]=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	mUpperLimit[3]=1-mLowerLimit[3];
	}

	if (mSettings->betamin() < mSettings->betamax()) { // beta

	if (mSettings->betamin() < mLowerLimit[0]) {
	cout << "Warning, requested lower limit of beta (" << mSettings->betamin() << ") "
	<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[0] << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mLowerLimit[0] = mSettings->betamin();
	}

	if (mSettings->betamax() > mUpperLimit[0]) {
	cout << "Warning, requested upper limit of Q2 (" << mSettings->betamax() << ") "
	<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[0] << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mUpperLimit[0] = mSettings->betamax();
	}
	}

	//
	// Check if any phase space is left
	//
	if (mLowerLimit[0] >= mUpperLimit[0]) {
	cout << "Invalid range in beta: beta=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
	exit(1);
	}
	if (mLowerLimit[1] >= mUpperLimit[1]) {
	cout << "Invalid range in Q2: Q2=[";
	cout << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl;
	exit(1);
	}
	if (mLowerLimit[2] >= mUpperLimit[2]) {
	cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl;
	exit(1);
	}

	//
	// Print-out limits (all verbose levels)
	//
	if (true){//}(mSettings->verbose()) {
	cout << "Sartre was thrown into the world, restricted by kinematics and user inputs alike:" << endl;
	cout << setw(10) << " beta=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
	cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
	cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl;
	cout << setw(10) << " z=[" << mLowerLimit[3] << ", " << mUpperLimit[3] << "]" << endl;
	cout << setw(10) << " MX2=[" << kineMX2min << ", " << kineMX2max << "]" << endl;
	cout << setw(10) << " y=[" << kineYmin << ", 1" << "]" << endl;

	}
	//
	// UNU.RAN needs the domain (boundaries) and the mode.
	// The domain is already defined, here we find the mode, which is tricky.
	// The max. cross-section is clearly at the domain boundary in Q2=Q2min.
	// The position in W2 and beta is not obvious. The mode should be at z=0.5.
	// The approach here is to use the BrentMinimizer1D that
	// performs first a scan a then a Brent fit.
	//
	mCrossSection->setCheckKinematics(false);
	double theMode[4];

	//
	// Find the mode.
	// Assume that the mode is at W2=W2max, Q2=Q2min, z=0.5
	// Start with beta=0.5 => MX2=Q2
	// It is kinemtically easier to use MX2 instead of beta here.
	//
	if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
	theMode[1] = mLowerLimit[1]; // Q2
	theMode[2] = mUpperLimit[2]; // W2
	theMode[3] = 0.5; // z

	double MX2min=mLowerLimit[1]*(1-mUpperLimit[0])/mUpperLimit[0];
	MX2min=max(kineMX2min, MX2min);
	double MX2max=mLowerLimit[1]*(1-mLowerLimit[0])/mLowerLimit[0];
	MX2max=min(kineMX2max, MX2max);
	InclusiveDiffractionModeFinderFunctor modeFunctor(mCrossSection, theMode[1], theMode[2], theMode[3], MX2min, MX2max);

	ROOT::Math::BrentMinimizer1D minimizer;
	minimizer.SetFunction(modeFunctor, MX2min, MX2max);
	minimizer.SetNpx(static_cast<int>(mUpperLimit[0]-mLowerLimit[0]));
	ok = minimizer.Minimize(100000, 0, 1.e-8);
	if (! ok) {
	cout << "Error, failed to find mode of pdf." << endl;
	exit(1);
	}
	theMode[0] = minimizer.XMinimum(); // Mx2
	+
	double MX2=theMode[0];
	theMode[0]=theMode[1]/(theMode[1]+MX2); //Mx2->beta
	double crossSectionAtMode = (*mCrossSection)(MX2, theMode[1], theMode[2], theMode[3]);
	if (mSettings->verbose()) {
	cout << "\tlocation: beta=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]) << ", z=" << theMode[3];
	cout << "; value: " << crossSectionAtMode << endl;
	}
	mCrossSection->setCheckKinematics(true);

	// domain and mode for Q2 -> log(Q2)
	mLowerLimit[1] = log(mLowerLimit[1]);
	mUpperLimit[1] = log(mUpperLimit[1]);
	theMode[1] = log(theMode[1]);


	if (mPDF_Functor) delete mPDF_Functor;
	if (mPDF) delete mPDF;

	mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &InclusiveDiffractiveCrossSections::unuranPDF, 4);
	mPDF = new TUnuranMultiContDist(*mPDF_Functor, true); // last arg = pdf in log or not
	mPDF->SetDomain(mLowerLimit, mUpperLimit);
	mPDF->SetMode(theMode);

	if (mUnuran) delete mUnuran;
	mUnuran = new TUnuran;

	mCrossSection->setCheckKinematics(false); // avoid numeric glitch in Init()
	mUnuran->Init(*mPDF, "method=hitro");
	mCrossSection->setCheckKinematics(true);
	mUnuran->SetSeed(mSettings->seed());

	//
	// Burn in generator
	//
	+
	double xrandom[4];
	for (int i=0; i<100; i++) {
	mUnuran->SampleMulti(xrandom);
	// cout<<"QQ: beta="<<xrandom[0];
	// cout<<" Q2="<<exp(xrandom[1]);
	// cout<<" W="<<sqrt(xrandom[2]);
	// cout<<" z="<<xrandom[3]<<endl;
	}
	+
	//
	// Set up QQG unu.ran:
	//
	//The mode is simpler for QQG, as it is at minimum beta
	//The limits are the same, except for z which is a different quantity in QQG.
	mLowerLimit_qqg[0]=mLowerLimit[0]; //beta
	mLowerLimit_qqg[1]=mLowerLimit[1]; //log(Q2)
	mLowerLimit_qqg[2]=mLowerLimit[2]; //W2
	//z>beta(1+4mf*mf/Q2);
	mLowerLimit_qqg[3]=mLowerLimit[0] * (1+4mfmf/exp(mUpperLimit[1]))+1e-10;

	mUpperLimit_qqg[0]=mUpperLimit[0]; //beta
	mUpperLimit_qqg[1]=mUpperLimit[1]; //log(Q2)
	mUpperLimit_qqg[2]=mUpperLimit[2]; //W2
	mUpperLimit_qqg[3]=mUpperLimit[3]; //z
	-
	+
	double theMode_qqg[4];
	theMode_qqg[0]=mLowerLimit_qqg[0]; //beta
	theMode_qqg[1]=mLowerLimit_qqg[1]; //log(Q2)
	theMode_qqg[2]=mUpperLimit_qqg[2]; //W2
	theMode_qqg[3]=theMode_qqg[0](1+4mf*mf/exp(theMode_qqg[1]))+1e-10; //z

	if (mPDF_Functor_qqg) delete mPDF_Functor_qqg;
	if (mPDF_qqg) delete mPDF_qqg;

	mPDF_Functor_qqg = new ROOT::Math::Functor(mCrossSection, &InclusiveDiffractiveCrossSections::unuranPDF_qqg, 4);
	mPDF_qqg = new TUnuranMultiContDist(*mPDF_Functor_qqg, true); // last arg = pdf in log or not
	//Rescale the domain in z (and fix it in the pdf)
	mPDF_qqg->SetDomain(mLowerLimit_qqg, mUpperLimit_qqg);
	mPDF_qqg->SetMode(theMode_qqg);

	mCrossSection->setCheckKinematics(false);
	if (mUnuran_qqg) delete mUnuran_qqg;
	mUnuran_qqg = new TUnuran;
	mUnuran_qqg->Init(*mPDF_qqg, "method=hitro");
	mUnuran_qqg->SetSeed(mSettings->seed());
	// mCrossSection->setCheckKinematics(true);
	//
	// Burn in QQG generator
	//
	for (int i=0; i<100; i++) {
	mUnuran_qqg->SampleMulti(xrandom);
	// cout<<"QQG: beta="<<xrandom[0];
	// cout<<" Q2="<<exp(xrandom[1]);
	// cout<<" W="<<sqrt(xrandom[2]);
	// cout<<" z="<<xrandom[3]<<endl;
	}
	+ mEventCounter = 0;
	+ mTriesCounter = 0;
	+ mIsInitialized = true;
	+ cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;
	+

	//
	// Calculate total cross section in given
	// kinematic limits. This will be used to
	// decide which type of event to generate,
	// q-qbar or q-qbar-g
	//
	-// cout<<"Calculating cross sections over the given kinematic range..."<<endl;
	-// cout<<"(This may take while, use your radical freedom and do other things.)"<<endl;
	-// mCrossSection->setFockState(QQ);
	-// mTotalCrossSection=0;
	-// mTotalCrossSection_qq=totalCrossSection();
	-// cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
	-// mCrossSection->setFockState(QQG);
	-// mTotalCrossSection=0;
	-// mTotalCrossSection_qqg=totalCrossSection();
	-// cout<<"The QQG cross section is: "<<mTotalCrossSection_qqg<<" nb. (2/2)"<<endl;
	-// cout<<"Total Cross Section: "<<mTotalCrossSection_qqg+mTotalCrossSection_qq<<" nb"<<endl;
	-// mTotalCrossSection=mTotalCrossSection_qqg+mTotalCrossSection_qq;
	-
	- mEventCounter = 0;
	- mTriesCounter = 0;
	- mIsInitialized = true;
	- cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;
	-
	+ cout<<"Calculating cross sections over the given kinematic range..."<<endl;
	+ cout<<"(This may take while, use your radical freedom and do other things.)"<<endl;
	+ mCrossSection->setFockState(QQ);
	+ mTotalCrossSection=0;
	+ mTotalCrossSection_qq=totalCrossSection();
	+ cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
	+ mCrossSection->setFockState(QQG);
	+ mTotalCrossSection=0;
	+ mTotalCrossSection_qqg=totalCrossSection();
	+ cout<<"The QQG cross section is: "<<mTotalCrossSection_qqg<<" nb. (2/2)"<<endl;
	+ cout<<"Total Cross Section: "<<mTotalCrossSection_qqg+mTotalCrossSection_qq<<" nb"<<endl;
	+ mTotalCrossSection=mTotalCrossSection_qqg+mTotalCrossSection_qq;
	return true;
	}


	double SartreInclusiveDiffraction::cross_section_wrapper(double* xx, double* par){
	double arg[4]={xx[0], xx[1], xx[2], par[0]};
	double result=mCrossSection->unuranPDF(arg);
	return -result;
	}

	bool SartreInclusiveDiffraction::init(const string& str) // overloaded version of init()
	{
	if (str.empty())
	return init();
	else
	return init(str.c_str());
	}

	vector<pair<double,double> > SartreInclusiveDiffraction::kinematicLimits()
	{
	vector<pair<double,double> > array;
	array.push_back(make_pair(mLowerLimit[0], mUpperLimit[0])); // t
	array.push_back(make_pair(exp(mLowerLimit[1]), exp(mUpperLimit[1]))); // Q2 or xpom
	if (!mSettings->UPC())
	array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W
	return array;
	}

	Event* SartreInclusiveDiffraction::generateEvent()
	{
	if (!mIsInitialized) {
	cout << "SartreInclusiveDiffraction::generateEvent(): Error, Sartre is not initialized yet." << endl;
	cout << " Call init() before trying to generate events." << endl;
	return 0;
	}
	//
	// Decide which Fock-state to generate:
	//
	bool isQQ=false;
	if(mRandom->Uniform(mTotalCrossSection) <= mTotalCrossSection_qq)
	isQQ=true;
	//
	// Generate one event
	//
	double xrandom[4];
	while (true) {
	mTriesCounter++;
	delete mCurrentEvent;
	mCurrentEvent = new Event;
	//
	// xrandom[i]
	// i=0: beta
	// i=1: Q2
	// i=2: W2
	// i=3: z
	//
	if(isQQ)
	mUnuran->SampleMulti(xrandom);
	else
	mUnuran_qqg->SampleMulti(xrandom);
	xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2
	double MX2=xrandom[1]*(1-xrandom[0])/xrandom[0];
	bool isValidEvent;
	isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2], xrandom[3], sqrt(MX2), true, (mSettings->verboseLevel() > 1));
	if (!isValidEvent) {
	if (mSettings->verboseLevel() > 2)
	cout << "SartreInclusiveDiffraction::generateEvent(): event rejected, not within kinematic limits" << endl;
	continue;
	}
	//
	// Fill beam particles in Event structure
	// Kinematics for eA is reported as 'per nucleon'
	//
	mCurrentEvent->eventNumber = mEventCounter;
	mCurrentEvent->beta = xrandom[0]; // beta
	mCurrentEvent->Q2 = xrandom[1]; // Q2
	mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
	mCurrentEvent->xpom=mCurrentEvent->x/xrandom[0]; //xpom=x/beta
	mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
	mCurrentEvent->s = mS; // s
	mCurrentEvent->W = sqrt(xrandom[2]);
	mCurrentEvent->z = xrandom[3];
	mCurrentEvent->MX = sqrt(MX2);
	if(isQQ)
	mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
	else
	mCurrentEvent->polarization = transverse;
	// mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
	mCurrentEvent->quarkSpecies = mCrossSection->quarkSpeciesOfLastCall();
	if(isQQ)
	mCurrentEvent->crossSectionRatioLT=mCrossSection->crossSectionRatioLTOfLastCall();
	else
	mCurrentEvent->crossSectionRatioLT=0;

	//
	// The following needs to be modified for QQG
	//
	double dsigmadb2=0;
	double z=xrandom[3];
	double eps2=z(1-z)(mCurrentEvent->Q2+MX2);
	double r_average=1./sqrt(eps2);
	if (mA == 1) {
	dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(r_average*hbarc, 0, mCurrentEvent->xpom);
	}
	else {
	dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(r_average*hbarc, 0, mCurrentEvent->xpom);
	}
	mCurrentEvent->Omega=-2.*log(1-dsigmadb2/2.);


	Particle eIn, hIn;
	eIn.index = 0;
	eIn.pdgId = 11; // e-
	eIn.status = 1;
	eIn.p = mElectronBeam;
	hIn.index = 1;
	hIn.pdgId = mNucleus->pdgID();
	hIn.status = 1;
	hIn.p = mHadronBeam;
	mCurrentEvent->particles.push_back(eIn);
	mCurrentEvent->particles.push_back(hIn);

	//
	// Generate the final state particles
	//
	bool ok;
	if (isQQ) {
	ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQ);
	}
	else {
	ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQG);
	}
	if (!ok) {
	if (mSettings->verboseLevel() > 1) cout << "SartreInclusiveDiffraction::generateEvent(): failed to generate final state" << endl;
	continue;
	}
	break;
	}

	mEventCounter++;

	//
	// Nuclear breakup
	//
	// If the event is incoherent the final state generator does produce a
	// 'virtual' proton with m > m_p which is used in Nucleus to calculate
	// the excitation energy and the boost.
	//
	int indexOfScatteredHadron = 6;
	bool allowBreakup = mSettings->enableNuclearBreakup();
	if (mA == 1) allowBreakup = false;

	if (mNucleus) mNucleus->resetBreakup(); // clear previous event in any case

	if (allowBreakup && mCurrentEvent->diffractiveMode == incoherent && mNucleus) {

	int nFragments = mNucleus->breakup(mCurrentEvent->particles[indexOfScatteredHadron].p);

	//
	// Merge the list of products into the event list.
	// We loose some information here. The user can always go back to
	// the nucleus and check the decay products for more details.
	// In the original list the energy is per nuclei, here we transform it
	// to per nucleon to stay consistent with Sartre conventions.
	//
	const vector<BreakupProduct>& products = mNucleus->breakupProducts();
	for (int i=0; i<nFragments; i++) {
	Particle fragment;
	fragment.index = mCurrentEvent->particles.size();
	fragment.pdgId = products[i].pdgId;
	fragment.status = 1;
	fragment.p = products[i].p*(1/static_cast<double>(products[i].A));
	fragment.parents.push_back(indexOfScatteredHadron);
	mCurrentEvent->particles.push_back(fragment);
	}
	}
	return mCurrentEvent;
	}

	double SartreInclusiveDiffraction::totalCrossSection()
	{
	if (mTotalCrossSection == 0) {
	//
	// Limits of integration in beta, Q2, W2, z
	//
	double xmin[4];
	double xmax[4];
	copy(mLowerLimit, mLowerLimit+4, xmin);
	copy(mUpperLimit, mUpperLimit+4, xmax);

	//
	// At this point mLowerLimit[1] and mUpperLimit[1]
	// are in log(Q2)
	//
	xmin[1] = exp(xmin[1]); // log Q2 limit
	xmax[1] = exp(xmax[1]); // log Q2 limit

	mTotalCrossSection = calculateTotalCrossSection(xmin, xmax);
	}
	return mTotalCrossSection;
	}

	double SartreInclusiveDiffraction::totalCrossSection(double lower[4], double upper[4]) // beta, Q2, W, z
	{
	lower[2] = lower[2]; upper[2] = upper[2]; // W -> W2
	double result = calculateTotalCrossSection(lower, upper);
	return result;
	}

	double SartreInclusiveDiffraction::integrationVolume(){
	//
	// Limits of integration in beta, Q2, W2, z
	//
	double xmin[4];
	double xmax[4];
	copy(mLowerLimit, mLowerLimit+4, xmin);
	copy(mUpperLimit, mUpperLimit+4, xmax);
	//
	// At this point mLowerLimit[1] and mUpperLimit[1]
	// are in log(Q2)
	//
	xmin[1] = exp(xmin[1]); // log Q2 limit
	xmax[1] = exp(xmax[1]); // log Q2 limit

	double volume = (xmax[0]-xmin[0]) * (xmax[1]-xmin[1]) * (xmax[2]-xmin[2]) * (xmax[3]-xmin[3]);
	return volume; //GeV4
	}

	EventGeneratorSettings* SartreInclusiveDiffraction::runSettings()
	{
	return EventGeneratorSettings::instance();
	}

	double SartreInclusiveDiffraction::crossSectionsClassWrapper(const double* var){
	return (*mCrossSection)(var);
	}

	double SartreInclusiveDiffraction::calculateTotalCrossSection(double lower[4], double upper[4])
	{
	double result = 0;

	if (!mIsInitialized) {
	cout << "SartreInclusiveDiffraction::calculateTotalCrossSection(): Error, Sartre is not initialized yet." << endl;
	cout << " Call init() before trying to generate events." << endl;
	return result;
	}
	//
	// Calculate integral using adaptive numerical method
	//
	// Options: ADAPTIVE, kVEGAS, kPLAIN, kMISER
	// no abs tolerance given -> relative only
	const double precision = 5e-4;
	// unsigned int numberofcalls = 1000000;
	unsigned int numberofcalls = 0;//1e5;
	ROOT::Math::Functor wfL(this, &SartreInclusiveDiffraction::crossSectionsClassWrapper, 4);
	ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE, 0, precision, numberofcalls);

	ig.SetFunction(wfL);
	result = ig.Integral(lower, upper);

	//
	// If it fails we switch to a MC integration which is usually more robust
	// although not as accurate. This should happen very rarely if at all.
	//
	if (result <= numeric_limits<float>::epsilon()) {
	cout << "SartreInclusiveDiffraction::calculateTotalCrossSection(): warning, adaptive integration failed - switching to VEGAS method." << endl;
	ROOT::Math::IntegratorMultiDim igAlt(ROOT::Math::IntegrationMultiDim::kVEGAS);
	igAlt.SetFunction(wfL);
	igAlt.SetRelTolerance(precision);
	igAlt.SetAbsTolerance(0);
	result = igAlt.Integral(lower, upper);
	}

	return result;
	}

	const FrangibleNucleus* SartreInclusiveDiffraction::nucleus() const {return mNucleus;}

	void SartreInclusiveDiffraction::listStatus(ostream& os) const
	{
	os << "Event summary: " << mEventCounter<< " events generated, " << mTriesCounter << " tried" << endl;
	time_t delta = runTime();
	os << "Total time used: " << delta/60 << " min " << delta - 60*(delta/60) << " sec" << endl;

	}

	time_t SartreInclusiveDiffraction::runTime() const
	{
	time_t now = time(0);
	return now-mStartTime;
	}

	//==============================================================================
	//
	// Utility functions and operators (helpers)
	//
	//==============================================================================

	ostream& operator<<(ostream& os, const TLorentzVector& v)
	{
	os << v.Px() << '\t' << v.Py() << '\t' << v.Pz() << '\t' << v.E() << '\t';
	double m2 = v*v;
	if (m2 < 0)
	os << '(' << -sqrt(-m2) << ')';
	else
	os << '(' << sqrt(m2) << ')';

	return os;
	}
	Index: trunk/src/Settings.h
	===================================================================
	--- trunk/src/Settings.h (revision 539)
	+++ trunk/src/Settings.h (revision 540)
	@@ -1,217 +1,223 @@
	//==============================================================================
	// Settings.h
	//
	// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
	//
	// This file is part of Sartre.
	//
	// This program is free software: you can redistribute it and/or modify
	// it under the terms of the GNU General Public License as published by
	// the Free Software Foundation.
	// This program is distributed in the hope that it will be useful,
	// but without any warranty; without even the implied warranty of
	// merchantability or fitness for a particular purpose. See the
	// GNU General Public License for more details.
	// You should have received a copy of the GNU General Public License
	// along with this program. If not, see <http://www.gnu.org/licenses/>.
	//
	// Author: Thomas Ullrich
	// Last update:
	// $Date$
	// $Author$
	//==============================================================================
	#ifndef Settings_h
	#define Settings_h
	#include <iostream>
	#include <string>
	#include <vector>
	#include <map>
	#include "TDatabasePDG.h"
	#include "TRandom3.h"
	#include "Enumerations.h"

	using namespace std;

	class TParticlePDG;

	class SettingsParameterBase {
	public:
	virtual ~SettingsParameterBase() {}

	public:
	string name;
	};

	template<typename T> class SettingsParameter : public SettingsParameterBase {
	public:
	T* address;
	T defaultValue;
	};

	class Settings {
	public:
	Settings();
	virtual ~Settings();

	bool readSettingsFromFile(const char*);
	virtual bool list(ostream& = cout);

	TParticlePDG* lookupPDG(int) const;
	string particleName(int pdgID);
	static TRandom3* randomGenerator();

	unsigned int seed() const;
	void setSeed(unsigned int);

	int userInt() const;
	double userDouble() const;
	string userString() const;

	void setUserInt(int);
	void setUserDouble(double);
	void setUserString(const string&);

	void setVerbose(bool);
	bool verbose() const;

	void setVerboseLevel(int);
	int verboseLevel() const;

	void setQ2min(double);
	double Q2min() const;
	double Qmin() const;

	void setQ2max(double);
	double Q2max() const;
	double Qmax() const;

	void setWmin(double);
	void setW2min(double);
	double Wmin() const;
	double W2min() const;

	void setWmax(double);
	void setW2max(double);
	double Wmax() const;
	double W2max() const;

	void setXpomMin(double); // UPC only
	void setXpomMax(double);
	double xpomMin() const;
	double xpomMax() const;

	void setbetamin(double); // Inclusive diffraction
	double betamin() const;
	void setbetamax(double);
	double betamax() const;

	void setZmin(double); // Inclusive diffraction
	double zmin() const;
	void setZmax(double);
	double zmax() const;

	int vectorMesonId() const;
	void setVectorMesonId(int);

	static double quarkMass(int);
	static void setQuarkMass(int, double);

	string dipoleModelName() const;

	DipoleModelType dipoleModelType() const;
	void setDipoleModelType(DipoleModelType);

	string dipoleModelParameterSetName() const;
	DipoleModelParameterSet dipoleModelParameterSet() const;
	void setDipoleModelParameterSet(DipoleModelParameterSet);

	string tableSetTypeName() const;
	TableSetType tableSetType() const;
	void setTableSetType(TableSetType);

	unsigned int A() const;
	void setA(unsigned int);

	string rootfile() const;
	void setRootfile(const char*);

	void setInclusiveDiffractionMode(bool);
	bool inclusiveDiffractionMode() const;
	bool exclusiveDiffractionMode() const;

	void setUPC(bool);
	bool UPC() const;

	+ void setUseInclusiveCrossSectionTables(bool);
	+ bool useInclusiveCrossSectionTables() const;
	+
	+
	void setUPCA(unsigned int);
	unsigned int UPCA() const;

	protected:
	template<typename T> void registerParameter(T, const char, T);
	virtual void consolidateSettings() = 0;
	virtual void consolidateCommonSettings();

	protected:
	vector<SettingsParameterBase*> mRegisteredParameters;
	static TRandom3 mRandomGenerator;
	unsigned int mSeed;

	static double mQuarkMass[6];

	bool mVerbose;
	int mVerboseLevel;

	double mQ2min;
	double mQ2max;
	double mWmin;
	double mWmax;
	double mXpomMin;
	double mXpomMax;
	double mbetamin;
	double mbetamax;
	double mZmin;
	double mZmax;

	int mVectorMesonId;

	string mDipoleModelName;
	DipoleModelType mDipoleModelType;
	DipoleModelParameterSet mDipoleModelParameterSet;
	string mDipoleModelParameterSetName;

	TableSetType mTableSetType;
	string mTableSetTypeName;

	unsigned int mA;

	string mRootfile;

	bool mUPC;
	unsigned int mUPCA;

	+
	bool mInclusiveDiffractionMode; // if false exclusive
	+ bool mUseInclusiveCrossSectionTables;

	private:
	string mRuncard;
	vector<string> mLines;
	TDatabasePDG *mPDG;
	map<int, string> mPeriodicTable;
	int mUserInt;
	double mUserDouble;
	string mUserString;
	};

	template<typename T> void Settings::registerParameter(T* var, const char* name, T def)
	{
	SettingsParameter<T>* sv = new SettingsParameter<T>;
	sv->address = var;
	sv->name = name;
	sv->defaultValue = def;
	*var = def;
	mRegisteredParameters.push_back(sv);
	}

	inline TRandom3* Settings::randomGenerator() {return &mRandomGenerator;}
	inline double Settings::quarkMass(int i) {return mQuarkMass[i];}
	inline void Settings::setQuarkMass(int i, double val) {mQuarkMass[i] = val;}

	#endif
	Index: trunk/src/Settings.cpp
	===================================================================
	--- trunk/src/Settings.cpp (revision 539)
	+++ trunk/src/Settings.cpp (revision 540)
	@@ -1,620 +1,625 @@
	//==============================================================================
	// Settings.cpp
	//
	// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
	//
	// This file is part of Sartre.
	//
	// This program is free software: you can redistribute it and/or modify
	// it under the terms of the GNU General Public License as published by
	// the Free Software Foundation.
	// This program is distributed in the hope that it will be useful,
	// but without any warranty; without even the implied warranty of
	// merchantability or fitness for a particular purpose. See the
	// GNU General Public License for more details.
	// You should have received a copy of the GNU General Public License
	// along with this program. If not, see <http://www.gnu.org/licenses/>.
	//
	// Author: Thomas Ullrich
	// Last update:
	// $Date$
	// $Author$
	//==============================================================================
	#include "Settings.h"
	#include "Constants.h"
	#include <typeinfo>
	#include <fstream>
	#include <sstream>
	#include <iomanip>
	#include <ctype.h>
	#include <cstdlib>
	#include "TParticlePDG.h"
	#include "TError.h"

	#define PR(x) cout << #x << " = " << (x) << endl;

	TRandom3 Settings::mRandomGenerator;

	// initialize static member
	double Settings::mQuarkMass[] = {0.0022, 0.0047, 0.096, 1.27, 4.18, 173.21}; // PDG values

	Settings::Settings()
	{
	mPDG = new TDatabasePDG;
	mPDG->ReadPDGTable();

	//
	// Setup name table (map) for nuclei
	//
	mPeriodicTable[1] = "H";
	mPeriodicTable[2] = "He";
	mPeriodicTable[3] = "Li";
	mPeriodicTable[4] = "Be";
	mPeriodicTable[5] = "B";
	mPeriodicTable[6] = "C";
	mPeriodicTable[7] = "N";
	mPeriodicTable[8] = "O";
	mPeriodicTable[9] = "F";
	mPeriodicTable[10] = "Ne";
	mPeriodicTable[11] = "Na";
	mPeriodicTable[12] = "Mg";
	mPeriodicTable[13] = "Al";
	mPeriodicTable[14] = "Si";
	mPeriodicTable[15] = "P";
	mPeriodicTable[16] = "S";
	mPeriodicTable[17] = "Cl";
	mPeriodicTable[18] = "Ar";
	mPeriodicTable[19] = "K";
	mPeriodicTable[20] = "Ca";
	mPeriodicTable[21] = "Sc";
	mPeriodicTable[22] = "Ti";
	mPeriodicTable[23] = "V";
	mPeriodicTable[24] = "Cr";
	mPeriodicTable[25] = "Mn";
	mPeriodicTable[26] = "Fe";
	mPeriodicTable[27] = "Co";
	mPeriodicTable[28] = "Ni";
	mPeriodicTable[29] = "Cu";
	mPeriodicTable[30] = "Zn";
	mPeriodicTable[31] = "Ga";
	mPeriodicTable[32] = "Ge";
	mPeriodicTable[33] = "As";
	mPeriodicTable[34] = "Se";
	mPeriodicTable[35] = "Br";
	mPeriodicTable[36] = "Kr";
	mPeriodicTable[37] = "Rb";
	mPeriodicTable[38] = "Sr";
	mPeriodicTable[39] = "Y";
	mPeriodicTable[40] = "Zr";
	mPeriodicTable[41] = "Nb";
	mPeriodicTable[42] = "Mo";
	mPeriodicTable[43] = "Tc";
	mPeriodicTable[44] = "Ru";
	mPeriodicTable[45] = "Rh";
	mPeriodicTable[46] = "Pd";
	mPeriodicTable[47] = "Ag";
	mPeriodicTable[48] = "Cd";
	mPeriodicTable[49] = "In";
	mPeriodicTable[50] = "Sn";
	mPeriodicTable[51] = "Sb";
	mPeriodicTable[52] = "Te";
	mPeriodicTable[53] = "I";
	mPeriodicTable[54] = "Xe";
	mPeriodicTable[55] = "Cs";
	mPeriodicTable[56] = "Ba";
	mPeriodicTable[57] = "La";
	mPeriodicTable[58] = "Ce";
	mPeriodicTable[59] = "Pr";
	mPeriodicTable[60] = "Nd";
	mPeriodicTable[61] = "Pm";
	mPeriodicTable[62] = "Sm";
	mPeriodicTable[63] = "Eu";
	mPeriodicTable[64] = "Gd";
	mPeriodicTable[65] = "Tb";
	mPeriodicTable[66] = "Dy";
	mPeriodicTable[67] = "Ho";
	mPeriodicTable[68] = "Er";
	mPeriodicTable[69] = "Tm";
	mPeriodicTable[70] = "Yb";
	mPeriodicTable[71] = "Lu";
	mPeriodicTable[72] = "Hf";
	mPeriodicTable[73] = "Ta";
	mPeriodicTable[74] = "W";
	mPeriodicTable[75] = "Re";
	mPeriodicTable[76] = "Os";
	mPeriodicTable[77] = "Ir";
	mPeriodicTable[78] = "Pt";
	mPeriodicTable[79] = "Au";
	mPeriodicTable[80] = "Hg";
	mPeriodicTable[81] = "Tl";
	mPeriodicTable[82] = "Pb";
	mPeriodicTable[83] = "Bi";
	mPeriodicTable[84] = "Po";
	mPeriodicTable[85] = "At";
	mPeriodicTable[86] = "Rn";
	mPeriodicTable[87] = "Fr";
	mPeriodicTable[88] = "Ra";
	mPeriodicTable[89] = "Ac";
	mPeriodicTable[90] = "Th";
	mPeriodicTable[91] = "Pa";
	mPeriodicTable[92] = "U";
	mPeriodicTable[93] = "Np";
	mPeriodicTable[94] = "Pu";
	mPeriodicTable[95] = "Am";
	mPeriodicTable[96] = "Cm";
	mPeriodicTable[97] = "Bk";
	mPeriodicTable[251] = "Cf";
	mPeriodicTable[252] = "Es";
	mPeriodicTable[100] = "Fm";
	mPeriodicTable[258] = "Md";
	mPeriodicTable[102] = "No";
	mPeriodicTable[103] = "Lr";
	mPeriodicTable[261] = "Rf";
	mPeriodicTable[105] = "Db";
	mPeriodicTable[106] = "Sg";
	mPeriodicTable[107] = "Bh";
	mPeriodicTable[108] = "Hs";
	mPeriodicTable[109] = "Mt";

	//
	// Register parameters (ptr, name, default)
	//
	registerParameter(&mUserInt, "userInt", 0);
	registerParameter(&mUserDouble, "userDouble", 0.);
	registerParameter(&mUserString, "userString", string(""));
	registerParameter(&mA, "A", static_cast<unsigned int>(1));
	registerParameter(&mVectorMesonId, "vectorMesonId", 443); // J/psi
	registerParameter(&mDipoleModelName, "dipoleModel", string("bSat"));
	registerParameter(&mDipoleModelParameterSetName, "dipoleModelParameterSet", string("KMW"));
	registerParameter(&mTableSetTypeName, "tableSetType", string("total_and_coherent"));
	registerParameter(&mQ2min, "Q2min", 1000.); // no limits if max <= min
	registerParameter(&mQ2max, "Q2max", 0.);
	registerParameter(&mWmin, "Wmin", 1000.);
	registerParameter(&mWmax, "Wmax", 0.);
	registerParameter(&mXpomMin, "xpomMin", 1e-8);
	registerParameter(&mXpomMax, "xpomMax", 1.);
	registerParameter(&mbetamin, "betamin", 0.);
	registerParameter(&mbetamax, "betamax", 1.);
	registerParameter(&mZmin, "zmin", 0.);
	registerParameter(&mZmax, "zmax", 1.);

	registerParameter(&mVerbose, "verbose", false);
	registerParameter(&mVerboseLevel, "verboseLevel", 0);
	registerParameter(&mRootfile, "rootfile", string("sartre.root"));
	registerParameter(&mSeed, "seed", static_cast<unsigned int>(time(0)));
	registerParameter(&mUPC, "UPC", false);
	registerParameter(&mUPCA, "UPCA", static_cast<unsigned int>(1));
	registerParameter(&mInclusiveDiffractionMode, "inclusiveDiffractionMode", false);
	+ registerParameter(&mUseInclusiveCrossSectionTables, "useInclusiveCrossSectionTables", false);
	+
	}

	Settings::~Settings()
	{
	for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
	delete mRegisteredParameters[k];
	}
	mRegisteredParameters.clear();
	delete mPDG;
	}

	int Settings::userInt() const {return mUserInt;}
	double Settings::userDouble() const {return mUserDouble;}
	string Settings::userString() const {return mUserString;}

	void Settings::setUserInt(int val) {mUserInt = val;}
	void Settings::setUserDouble(double val) {mUserDouble = val;}
	void Settings::setUserString(const string& val) {mUserString = val;}

	unsigned int Settings::seed() const {return mSeed;}

	void Settings::setSeed(unsigned int val)
	{
	mSeed = val;
	mRandomGenerator.SetSeed(mSeed);
	gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()
	}

	bool Settings::readSettingsFromFile(const char *file)
	{
	if (!file) return false; // nothing to do
	mRuncard = string(file);
	mLines.clear();

	//
	// Open file
	//
	ifstream ifs(file);
	if (!ifs) {
	cout << "Settings::readSettingsFromFile(): error, cannot open file '"
	<< file << "'." << endl;
	return false;
	}

	//
	// Read file into vector of strings, skip comments and empty lines
	//
	while (ifs.good() && !ifs.eof()) {
	string line;
	getline (ifs, line);
	if (ifs.eof() && line.empty()) break;
	// empty line
	if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) continue;
	// if first character is not a letter/digit, then taken to be a comment.
	int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
	if (!isalnum(line[firstChar])) continue;

	mLines.push_back(line);
	}
	ifs.close(); // done with I/O

	//
	// Process vector of strings one at a time and use
	// it to set registered variables.
	//
	for (unsigned int i=0; i<mLines.size(); i++) {

	// replace '=' by blank to make parsing simpler.
	while (mLines[i].find("=") != string::npos) {
	int firstEqual = mLines[i].find_first_of("=");
	mLines[i].replace(firstEqual, 1, " ");
	}

	// replace ':' by blank to make parsing simpler.
	while (mLines[i].find(":") != string::npos) {
	int firstEqual = mLines[i].find_first_of(":");
	mLines[i].replace(firstEqual, 1, " ");
	}

	// get identifier string
	istringstream splitLine(mLines[i]);
	string name;
	splitLine >> name;

	// find value string
	string valueString;
	splitLine >> valueString;
	if (!splitLine) {
	cout << "Settings::readSettingsFromFile(): error, value of variable '"
	<< name.c_str() << "' not recognized." << endl;
	}
	istringstream modeData(valueString);

	//
	// Loop over registered variables and see which fits.
	// Not particular elegant but does the job and saves
	// a lot of programming in derived classes.
	//
	bool isRegistered = false;
	for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
	if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<vector<double>>)) { // test
	SettingsParameter<vector<double>> var = dynamic_cast<SettingsParameter<vector<double>>> (mRegisteredParameters[k]);
	if (var->name == name) {
	var->address->push_back(atof(valueString.c_str())); // first value
	while (splitLine.good() && !splitLine.eof()) { // get remaining
	string nextValue;
	splitLine >> nextValue;
	var->address->push_back(atof(nextValue.c_str()));
	if (splitLine.eof()) break;
	}
	isRegistered = true;
	}
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<double>)) {
	SettingsParameter<double> var = dynamic_cast<SettingsParameter<double>> (mRegisteredParameters[k]);
	if (var->name == name) {
	modeData >> (*(var->address));
	isRegistered = true;
	}
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<int>)) {
	SettingsParameter<int> var = dynamic_cast<SettingsParameter<int>> (mRegisteredParameters[k]);
	if (var->name == name) {
	modeData >> (*(var->address));
	isRegistered = true;
	}
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned int>)) {
	SettingsParameter<unsigned int> var = dynamic_cast<SettingsParameter<unsigned int>> (mRegisteredParameters[k]);
	if (var->name == name) {
	modeData >> (*(var->address));
	isRegistered = true;
	}
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned long>)) {
	SettingsParameter<unsigned long> var = dynamic_cast<SettingsParameter<unsigned long>> (mRegisteredParameters[k]);
	if (var->name == name) {
	modeData >> (*(var->address));
	isRegistered = true;
	}
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<string>)) {
	SettingsParameter<string> var = dynamic_cast<SettingsParameter<string>> (mRegisteredParameters[k]);
	if (var->name == name) {
	*(var->address) = valueString;
	isRegistered = true;
	}
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<bool>)) {
	SettingsParameter<bool> var = dynamic_cast<SettingsParameter<bool>> (mRegisteredParameters[k]);
	if (var->name == name) {
	isRegistered = true;
	if (valueString == string("true") \|\|
	valueString == string("True") \|\|
	valueString == string("TRUE") \|\|
	valueString == string("on") \|\|
	valueString == string("On") \|\|
	valueString == string("ON") \|\|
	valueString == string("Yes") \|\|
	valueString == string("yes") \|\|
	valueString == string("YES") \|\|
	valueString == string("T") \|\|
	valueString == string("t") \|\|
	valueString == string("1") ) {
	*(var->address) = true;
	}
	else {
	*(var->address) = false;
	}
	}
	}
	}
	if (!isRegistered) {
	cout << "Settings::readSettingsFromFile(): error, parameter identifier '" <<
	name.c_str() << "' not recognized." << endl;
	}
	}

	//
	// Consolidate input (after burner)
	//
	consolidateCommonSettings();
	consolidateSettings(); // overloaded

	return true;
	}

	bool Settings::list(ostream& os)
	{
	const int fieldWidth = 28;
	os << "\nRun Settings:" << endl;
	for (unsigned int k=0; k<mRegisteredParameters.size(); k++) {
	if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<double>)) {
	SettingsParameter<double> var = dynamic_cast<SettingsParameter<double>> (mRegisteredParameters[k]);
	os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<int>)) {
	SettingsParameter<int> var = dynamic_cast<SettingsParameter<int>> (mRegisteredParameters[k]);
	os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned int>)) {
	SettingsParameter<unsigned int> var = dynamic_cast<SettingsParameter<unsigned int>> (mRegisteredParameters[k]);
	os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<unsigned long>)) {
	SettingsParameter<unsigned long> var = dynamic_cast<SettingsParameter<unsigned long>> (mRegisteredParameters[k]);
	os << setw(fieldWidth) << var->name.c_str() << "\t" << *(var->address) << endl;
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<string>)) {
	SettingsParameter<string> var = dynamic_cast<SettingsParameter<string>> (mRegisteredParameters[k]);
	os << setw(fieldWidth) << var->name.c_str() << "\t" << var->address->c_str() << endl;
	}
	else if (typeid(*mRegisteredParameters[k]) == typeid(SettingsParameter<bool>)) {
	SettingsParameter<bool> var = dynamic_cast<SettingsParameter<bool>> (mRegisteredParameters[k]);
	os << setw(fieldWidth) << var->name.c_str() << "\t" << (*(var->address) ? "true" : "false") << endl;
	}
	}
	os << endl;
	return true;
	}

	TParticlePDG* Settings::lookupPDG(int id) const
	{
	if (mPDG)
	return mPDG->GetParticle(id);
	else
	return 0;
	}

	string Settings::particleName(int pdgID)
	{
	string name("unknown");

	if (abs(pdgID) < 1000000000) { // particle
	if (mPDG) {
	TParticlePDG *part = lookupPDG(pdgID);
	if (part) name = part->GetName();
	}
	if (abs(pdgID) == 990) name = "pomeron";
	}
	else { // nucleus in 10LZZZAAAI PDG format
	int id = pdgID;
	// int iso = id%10;
	id /= 10;
	int A = id%1000;
	id /= 1000;
	int Z = id%1000;
	stringstream namestream;
	namestream << mPeriodicTable[Z] << "(" << A << ")";
	name = namestream.str();
	}
	return name;
	}

	void Settings::setVerbose(bool val) {
	mVerbose = val;
	if (mVerbose && mVerboseLevel == 0) setVerboseLevel(1);
	if (!mVerbose && mVerboseLevel != 0) setVerboseLevel(0);
	}

	bool Settings::verbose() const {return mVerbose;}

	void Settings::setVerboseLevel(int val) {
	mVerboseLevel = val;
	if (!mVerbose && mVerboseLevel != 0) mVerbose = true;
	if (mVerbose && mVerboseLevel == 0) mVerbose = false;
	//
	// Unless verbose level is 5 or higher we suppress the many redundant
	// ROOT messages from algorithms since most of their errors are handled
	// internally anyway.
	//
	if (mVerboseLevel < 5) gErrorIgnoreLevel = 5000;
	}

	int Settings::verboseLevel() const {return mVerboseLevel;}

	void Settings::setQ2min(double val) { mQ2min = val;}
	double Settings::Q2min() const {return mQ2min;}
	double Settings::Qmin() const {return sqrt(mQ2min);}

	void Settings::setQ2max(double val) { mQ2max = val;}
	double Settings::Q2max() const {return mQ2max;}
	double Settings::Qmax() const {return sqrt(mQ2max);}

	void Settings::setW2min(double val) { mWmin = sqrt(val);}
	void Settings::setWmin(double val) { mWmin = val;}
	double Settings::Wmin() const {return mWmin;}
	double Settings::W2min() const {return mWmin*mWmin;}

	void Settings::setW2max(double val) { mWmax = sqrt(val);}
	void Settings::setWmax(double val) { mWmax = val;}
	double Settings::Wmax() const {return mWmax;}
	double Settings::W2max() const {return mWmax*mWmax;}

	void Settings::setXpomMin(double val) {mXpomMin = val;}
	void Settings::setXpomMax(double val) {mXpomMax = val;}
	double Settings::xpomMin() const {return mXpomMin;}
	double Settings::xpomMax() const {return mXpomMax;}

	void Settings::setbetamin(double val) {mbetamin = val;}
	void Settings::setbetamax(double val) {mbetamax = val;}
	double Settings::betamin() const {return mbetamin;}
	double Settings::betamax() const {return mbetamax;}

	void Settings::setZmin(double val) {mZmin = val;}
	void Settings::setZmax(double val) {mZmax = val;}
	double Settings::zmin() const {return mZmin;}
	double Settings::zmax() const {return mZmax;}


	int Settings::vectorMesonId() const {return mVectorMesonId;}
	void Settings::setVectorMesonId(int val) {mVectorMesonId = val;}

	string Settings::dipoleModelName() const {return mDipoleModelName;}
	DipoleModelType Settings::dipoleModelType() const {return mDipoleModelType;}
	void Settings::setDipoleModelType(DipoleModelType val)
	{
	mDipoleModelType = val;
	if (mDipoleModelType == bSat)
	mDipoleModelName = string("bSat");
	else if (mDipoleModelType == bNonSat)
	mDipoleModelName = string("bNonSat");
	else if (mDipoleModelType == bCGC)
	mDipoleModelName = string("bCGC");
	}

	unsigned int Settings::A() const {return mA;}
	void Settings::setA(unsigned int val) {mA = val;}

	void Settings::setRootfile(const char* val){ mRootfile = val; }
	string Settings::rootfile() const { return mRootfile; }

	string Settings::dipoleModelParameterSetName() const {return mDipoleModelParameterSetName;}

	DipoleModelParameterSet Settings::dipoleModelParameterSet() const {return mDipoleModelParameterSet;}

	void Settings::setDipoleModelParameterSet(DipoleModelParameterSet val)
	{
	mDipoleModelParameterSet = val;
	if (mDipoleModelParameterSet == KMW)
	mDipoleModelParameterSetName = string("KMW");
	else if (mDipoleModelParameterSet == HMPZ)
	mDipoleModelParameterSetName = string("HMPZ");
	else if (mDipoleModelParameterSet == STU)
	mDipoleModelParameterSetName = string("STU");
	else if (mDipoleModelParameterSet == CUSTOM)
	mDipoleModelParameterSetName = string("CUSTOM");
	}

	string Settings::tableSetTypeName() const {return mTableSetTypeName;}

	TableSetType Settings::tableSetType() const {return mTableSetType;}

	void Settings::setTableSetType(TableSetType val)
	{
	mTableSetType = val;
	if (mTableSetType == total_and_coherent)
	mTableSetTypeName = string("total_and_coherent");
	else if (mTableSetType == coherent_and_incoherent)
	mTableSetTypeName = string("coherent_and_incoherent");
	}

	void Settings::consolidateCommonSettings()
	{
	//
	// Check if verbose levels and flags are consistent
	// The verbose flag superseeds the verboseLevel.
	//
	if (mVerbose && mVerboseLevel == 0) setVerboseLevel(1);
	if (mVerboseLevel != 0 && !mVerbose) setVerboseLevel(0);
	setVerboseLevel(verboseLevel()); // invoke method no matter what

	//
	// Set random generator seed
	//
	mRandomGenerator.SetSeed(mSeed);
	gRandom->SetSeed(mSeed); // needed for TH1::GetRandom()

	//
	// Dipole Model
	//
	if (mDipoleModelName == string("bSat"))
	mDipoleModelType = bSat;
	else if (mDipoleModelName == string("bNonSat"))
	mDipoleModelType = bNonSat;
	else if (mDipoleModelName == string("bCGC"))
	mDipoleModelType = bCGC;
	else {
	cout << "Settings::consolidateCommonSettings(): Error, dipole model '"
	<< mDipoleModelName.c_str() << "' is not defined." << endl;
	exit(1);
	}

	//
	// Dipole Model Parameter Set
	//
	if (mDipoleModelParameterSetName == string("KMW"))
	mDipoleModelParameterSet = KMW;
	else if (mDipoleModelParameterSetName == string("HMPZ"))
	mDipoleModelParameterSet = HMPZ;
	else if (mDipoleModelParameterSetName == string("STU"))
	mDipoleModelParameterSet = STU;
	else if (mDipoleModelParameterSetName == string("CUSTOM"))
	mDipoleModelParameterSet = CUSTOM;
	else {
	cout << "Settings::consolidateCommonSettings(): Error, dipole model parameter set'"
	<< mDipoleModelParameterSetName.c_str() << "' is not defined." << endl;
	exit(1);
	}

	//
	// Table Set Type
	//
	if (mTableSetTypeName == string("total_and_coherent"))
	mTableSetType = total_and_coherent;
	else if (mTableSetTypeName == string("coherent_and_incoherent"))
	mTableSetType = coherent_and_incoherent;
	else {
	cout << "Settings::consolidateCommonSettings(): Error, table set type '"
	<< mTableSetTypeName.c_str() << "' is not defined." << endl;
	exit(1);
	}
	}

	void Settings::setInclusiveDiffractionMode(bool val) { mInclusiveDiffractionMode = val; }
	bool Settings::inclusiveDiffractionMode() const { return mInclusiveDiffractionMode; }
	bool Settings::exclusiveDiffractionMode() const { return !mInclusiveDiffractionMode; }

	+void Settings::setUseInclusiveCrossSectionTables(bool val) { mUseInclusiveCrossSectionTables = val; }
	+bool Settings::useInclusiveCrossSectionTables() const { return mUseInclusiveCrossSectionTables; }
	+
	void Settings::setUPC(bool val) { mUPC = val; }
	bool Settings::UPC() const { return mUPC; }

	void Settings::setUPCA(unsigned int val){ mUPCA = val; }
	unsigned int Settings::UPCA() const { return mUPCA; }
	Index: trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp
	===================================================================
	--- trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 539)
	+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 540)
	@@ -1,271 +1,276 @@
	#include "InclusiveDiffractiveCrossSectionsFromTables.h"
	#include <TFile.h>
	#include "Math/SpecFunc.h"
	#include "TFile.h"
	#include <iostream>
	#include <cmath>
	#include <algorithm>
	#include <algorithm>
	#include "Constants.h"
	#include "Nucleus.h"
	#include "DipoleModel.h"
	#include "AlphaStrong.h"
	#include "Math/IntegratorMultiDim.h"
	#include "Math/Functor.h"
	#include "TMath.h"
	#include "TableGeneratorSettings.h"
	#include "Settings.h"
	#include "Enumerations.h"
	#include "DglapEvolution.h"
	#include "Math/WrappedTF1.h"
	#include "Math/GaussIntegrator.h"
	#include "Kinematics.h"
	#include "linterp.h"

	#define PR(x) cout << #x << " = " << (x) << endl;

	using namespace std;

	InclusiveDiffractionCrossSectionsFromTables::InclusiveDiffractionCrossSectionsFromTables(){

	- TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
	- TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
	- TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
	- TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
	-
	- TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
	- TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
	- TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
	- TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
	-
	- TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
	- TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
	- TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
	- TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
	+ TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
	+ TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
	+ TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
	+ TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
	+
	+ TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
	+ TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
	+ TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
	+ TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
	+
	+ TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
	+ TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
	+ TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
	+ TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");

	+
	tableQQ_T.clear();
	tableQQ_L.clear();
	tableQQG.clear();
	-
	+
	// Retrieve the histograms
	THnF* histQQ_T0 = dynamic_cast<THnF*>(fileQQ_T0->Get("table"));
	THnF* histQQ_T1 = dynamic_cast<THnF*>(fileQQ_T1->Get("table"));
	THnF* histQQ_T2 = dynamic_cast<THnF*>(fileQQ_T2->Get("table"));
	THnF* histQQ_T3 = dynamic_cast<THnF*>(fileQQ_T3->Get("table"));
	-
	+
	THnF* histQQ_L0 = dynamic_cast<THnF*>(fileQQ_L0->Get("table"));
	THnF* histQQ_L1 = dynamic_cast<THnF*>(fileQQ_L1->Get("table"));
	THnF* histQQ_L2 = dynamic_cast<THnF*>(fileQQ_L2->Get("table"));
	THnF* histQQ_L3 = dynamic_cast<THnF*>(fileQQ_L3->Get("table"));

	THnF* histQQG_0 = dynamic_cast<THnF*>(fileQQG_0->Get("table"));
	THnF* histQQG_1 = dynamic_cast<THnF*>(fileQQG_1->Get("table"));
	THnF* histQQG_2 = dynamic_cast<THnF*>(fileQQG_2->Get("table"));
	THnF* histQQG_3 = dynamic_cast<THnF*>(fileQQG_3->Get("table"));
	-
	+
	tableQQ_T.push_back(histQQ_T0);
	tableQQ_T.push_back(histQQ_T1);
	tableQQ_T.push_back(histQQ_T2);
	tableQQ_T.push_back(histQQ_T3);

	tableQQ_L.push_back(histQQ_L0);
	tableQQ_L.push_back(histQQ_L1);
	tableQQ_L.push_back(histQQ_L2);
	tableQQ_L.push_back(histQQ_L3);

	tableQQG.push_back(histQQG_0);
	tableQQG.push_back(histQQG_1);
	tableQQG.push_back(histQQG_2);
	tableQQG.push_back(histQQG_3);
	-
	+
	}
	InclusiveDiffractionCrossSectionsFromTables::~InclusiveDiffractionCrossSectionsFromTables(){
	}

	double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {

	double MX2 = Q2*(1-beta)/beta;
	// double mf = Settings::quarkMass(mQuarkIndex);
	double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
	double sqrtArg = 1.-4.mfmf/MX2;
	if (sqrtArg<0){
	if (mSettings->verboseLevel() > 2)
	cout<<"There is no phase-space for z, return 0."<<endl;
	return 0;
	}
	double z0 = (1.-sqrt(sqrtArg))/2.;
	if (z<z0 or z>1-z0){
	if (mSettings->verboseLevel() > 2)
	cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
	return 0;
	}
	double pt2 = z(1-z)Q2-mf*mf;
	if (pt2<0){
	if (mSettings->verboseLevel() > 2)
	cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
	return 0;
	}
	if(z(1-z)MX2-mf*mf<0){ //kappa2<0
	if (mSettings->verboseLevel() > 2)
	cout<<"Not enough phase-space for quark production, return 0."<<endl;
	return 0;
	}
	double result = interpolate(beta, Q2, W2, z, transverse, QQ, mQuarkIndex); //nb

	return result;
	}

	double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L(double beta, double Q2, double W2, double z){

	double MX2 = Q2*(1-beta)/beta;
	// double mf = Settings::quarkMass(mQuarkIndex);
	double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
	double pt2 = z(1-z)Q2-mf*mf;
	if (pt2<0){
	if (mSettings->verboseLevel() > 2)
	cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
	return 0;
	}
	double sqrtArg=1.-4.mfmf/MX2;
	if (sqrtArg<0){
	if (mSettings->verboseLevel() > 2)
	cout<<"There is no phase-space for z, return 0."<<endl;
	return 0;
	}
	double z0=(1.-sqrt(sqrtArg))/2.;
	if (z<z0 or z>1-z0){
	if (mSettings->verboseLevel() > 2)
	cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
	return 0;
	}
	if(z(1-z)MX2-mf*mf<0){ //kappa2<0
	if (mSettings->verboseLevel() > 2)
	cout<<"Not enough phase-space for quark production, return 0."<<endl;
	return 0;
	}
	double result = interpolate(beta, Q2, W2, z, longitudinal, QQ, mQuarkIndex);
	return result;
	}

	//===================================================
	// QQG Calculations
	//===================================================

	double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
	double result = interpolate(beta, Q2, W2, ztilde, transverse, QQG, mQuarkIndex);
	return result;
	}

	double InclusiveDiffractionCrossSectionsFromTables::interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark){
	-// double scaleFactor=(0.9-0.1)/10;
	- beta = beta*0.8/0.9+0.1;
	THnF* hist4D=0;
	- if(pol==transverse and fock==QQ){
	+
	+ //#TT TODO: There seem to be a strange swap happening
	+ // between T QQ and L QQ...
	+ if(pol==longitudinal and fock==QQ){
	hist4D=tableQQ_T.at(quark);
	}
	- else if(pol==longitudinal and fock==QQ){
	+ else if(pol==transverse and fock==QQ){
	hist4D=tableQQ_L.at(quark);
	}
	else if(fock==QQG){
	hist4D=tableQQG.at(quark);
	}
	else{
	cout<<"InclusiveDiffractionCrossSectionsFromTables::interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark): There is nothing to interpolate, ending"<<endl;
	exit(2);
	}
	-
	+
	int nDims = hist4D->GetNdimensions();
	-
	-// PR(nDims);
	+
	+// PR(nDims);
	+
	vector <double>MinRange(nDims); //beta, Q2,W2, z
	vector <double>MaxRange(nDims); //beta, Q2,W2, z
	vector <int>Bins(nDims); //beta, Q2,W2, z
	for (int dim = 0; dim < 4; dim++) {
	TAxis* axis = hist4D->GetAxis(dim); // Get the axis for this dimension
	if (!axis) {
	std::cerr << "Error: Cannot retrieve axis for dimension " << dim << std::endl;
	continue;
	}
	// Get the number of bins and range
	int nBins = axis->GetNbins();
	double minRange = axis->GetXmin();
	double maxRange = axis->GetXmax();
	-// PR(minRange);
	-// PR(maxRange);
	Bins[dim]= nBins;
	-// PR(nBins);
	MinRange[dim]= minRange;
	MaxRange[dim]= maxRange;
	}
	// construct the grid in each dimension.
	// note that we will pass in a sequence of iterators pointing to the beginning of each grid
	std::vector<double> gridBeta = mgrid(hist4D, 0, Bins[0]);
	std::vector<double> gridQ2 = mgrid(hist4D, 1, Bins[1]) ;
	std::vector<double> gridW2 = mgrid(hist4D, 2, Bins[2]);
	std::vector<double> gridZ = mgrid(hist4D, 3, Bins[3]) ;
	-
	-// print_grid_vector(gridQ2, "gridQ2");
	-// print_grid_vector(gridW2, "gridW2");
	-// print_grid_vector(gridBeta, "gridBeta");
	-// print_grid_vector(gridZ, "gridZ");
	-
	+
	std::vector< std::vector<double>::iterator > grid_iter_list;
	grid_iter_list.push_back(gridBeta.begin());
	grid_iter_list.push_back(gridQ2.begin());
	grid_iter_list.push_back(gridW2.begin());
	grid_iter_list.push_back(gridZ.begin());
	-
	- // cout<<"grid\t"<<*(gridQ2.begin())<<endl;
	+
	// the size of the grid in each dimension
	array<int,4> grid_sizes;
	- grid_sizes[0] = Bins[0];
	- grid_sizes[1] = Bins[1];
	- grid_sizes[2] = Bins[2];
	- grid_sizes[3] = Bins[3];
	-
	-
	- // for(unsigned long i=0;i<grid_sizes.size(); i++)
	- // std::cout<<"Grid size \t"<<grid_sizes[i]<<std::endl;
	+ for(unsigned long i=0; i<grid_sizes.size(); i++){
	+ grid_sizes[i] = Bins[i];
	+ }

	// total number of elements
	int num_elements = grid_sizes[0] * grid_sizes[1]* grid_sizes[2]* grid_sizes[3];
	-// std::cout<<"num_elements = "<<num_elements<<std::endl;
	+
	// fill in the values of f(x) at the gridpoints.
	// we will pass in a contiguous sequence, values are assumed to be laid out C-style
	std::vector<double> f_values(num_elements);
	int count = 0;
	+ int countExtreme=0;
	for (int ibeta=1; ibeta<=grid_sizes[0]; ibeta++) {
	for (int iQ2=1; iQ2<=grid_sizes[1]; iQ2++) {
	for (int iW2=1; iW2<=grid_sizes[2]; iW2++) {
	for (int iZ=1; iZ<=grid_sizes[3]; iZ++) {
	- int binIndex[] = {ibeta, iQ2, iW2, iZ};
	- f_values[count] = hist4D->GetBinContent(binIndex);
	+ int binIndex[4] = {ibeta-1, iQ2-1, iW2-1, iZ-1};
	+ double val=hist4D->GetBinContent(binIndex);
	+ if(exp(val)>1e20){ //There seems to be some spikes
	+ countExtreme++;
	+ val=0;
	+ }
	+ if(val!=0)
	+ f_values[count] = val;//exp(val);
	+ else
	+ f_values[count] = 0.;
	count++;
	}
	}
	}
	+
	}
	-// PR(count);
	- // std::cout<<"done"<<std::endl;
	-
	- // construct the interpolator. the last two arguments are pointers to the underlying data
	+// if(countExtreme){
	+// PR(countExtreme);
	+//// exit(3);
	+// }
	+ //
	+ // Construct the interpolator.
	+ // The last two arguments are pointers to the underlying data
	+ //
	InterpMultilinear <4, double> interp_ML(grid_iter_list.begin(), grid_sizes.begin(), f_values.data(), f_values.data() + num_elements);

	- // InterpMultilinear interp_ML(grid_iter_list.begin(), grid_sizes.begin(), f_values.data(), f_values.data() + num_elements);
	-
	// interpolate:
	array<double,4> args = {beta, Q2, W2, z};
	double result=interp_ML.interp(args.begin());
	- return exp(result);
	+
	+// return result;
	+ if(result==0)
	+ return 0.;
	+ else
	+ return exp(result);
	}

	vector<double> InclusiveDiffractionCrossSectionsFromTables::mgrid(THnF* hist,int axis, int len){
	vector<double> gridtype(len);
	- for(int i = 0; i <len; i++){
	- double val = hist->GetAxis(axis)->GetBinCenter(i+1);
	- gridtype[i] = val;
	-// PR(val);
	+ for(int i = 1; i <=len; i++){
	+ double val = hist->GetAxis(axis)->GetBinCenter(i);
	+ gridtype[i-1] = val;
	}
	return gridtype;
	}

	Index: trunk/src/InclusiveDiffractiveCrossSections.cpp
	===================================================================
	--- trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 539)
	+++ trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 540)
	@@ -1,729 +1,728 @@
	//==============================================================================
	// InclusiveDiffractiveCrossSection.cpp
	//
	// Copyright (C) 2024 Tobias Toll and Thomas Ullrich
	//
	// This file is part of Sartre.
	//
	// This program is free software: you can redistribute it and/or modify
	// it under the terms of the GNU General Public License as published by
	// the Free Software Foundation.
	// This program is distributed in the hope that it will be useful,
	// but without any warranty; without even the implied warranty of
	// merchantability or fitness for a particular purpose. See the
	// GNU General Public License for more details.
	// You should have received a copy of the GNU General Public License
	// along with this program. If not, see <http://www.gnu.org/licenses/>.
	//
	// Author: Tobias Toll
	// Last update:
	// $Date: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
	// $Author: ullrich $
	//==============================================================================
	#include "Math/SpecFunc.h"
	#include "TFile.h"
	#include <iostream>
	#include <cmath>
	#include <algorithm>
	#include <algorithm>
	#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 "Settings.h"
	#include "Enumerations.h"
	#include "TF1.h"
	#include "TF1.h"
	#include "TH1F.h"
	#include "cuba.h"
	#include "InclusiveDiffractiveCrossSections.h"
	#include "DglapEvolution.h"
	#include "Math/WrappedTF1.h"
	#include "Math/GaussIntegrator.h"
	#include "DipoleModelParameters.h"

	#define PR(x) cout << #x << " = " << (x) << endl;

	InclusiveDiffractiveCrossSections::InclusiveDiffractiveCrossSections(){
	mSettings = EventGeneratorSettings::instance();
	mRandom = mSettings->randomGenerator();
	mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
	mPhotonFlux.setS(mS);
	mCheckKinematics = true;
	mCrossSectionRatioLT = 0;

	mQ2=0;

	// TableGeneratorSettings* settings = TableGeneratorSettings::instance();
	if(mDipoleModel) delete mDipoleModel;

	mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
	DipoleModelType model=mSettings->dipoleModelType();
	if (model==bSat) {
	mDipoleModel = new DipoleModel_bSat(mSettings);
	}
	else if(model==bNonSat){
	mDipoleModel = new DipoleModel_bNonSat(mSettings);
	}
	else {
	cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
	exit(1);
	}
	mA=mSettings->A();

	}

	InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){}

	double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_total(){return mDsigdbetadQ2dWdz_L_total;}
	double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_L_total(){return mDsigdbetadQ2dWdz_T_total;}
	double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_qqg(){ return mDsigdbetadQ2dWdz_T_qqg;};

	//These three functions are unique for the sub-classes:
	double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double, double, double, double){ return 0; }
	double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double, double, double, double){ return 0; }
	double InclusiveDiffractiveCrossSections::dsigmadbetadz_QQG(double, double, double, double){return 0;}


	void InclusiveDiffractiveCrossSections::setCheckKinematics(bool val) {mCheckKinematics = val;}

	void InclusiveDiffractiveCrossSections::setFockState(FockState val){
	mFockState=val;
	}

	FockState InclusiveDiffractiveCrossSections::getFockState(){
	return mFockState;
	}

	double InclusiveDiffractiveCrossSections::operator()(double MX2, double Q2, double W2, double z){
	return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
	}

	double InclusiveDiffractiveCrossSections::operator()(const double* array){
	double result=0;
	if(mFockState==QQ or mFockState==ALL){
	result+=dsigdbetadQ2dW2dz_total_checked(array[0], array[1], array[2], array[3]);
	}
	else if(mFockState==QQG or
	mFockState==ALL){
	result+=dsigdbetadQ2dW2dz_total_qqg_checked(array[0], array[1], array[2], array[3]);
	}
	else{
	cout<<"InclusiveDiffractiveCrossSectionsIntegrals::operator(): Fockstate is invalid, stopping"<<endl;
	exit(0);
	}
	return result;
	}

	double InclusiveDiffractiveCrossSections::unuranPDF(const double* array){
	//
	// array is: beta, log(Q2), W2, z
	//
	double result = 0;
	double Q2 = exp(array[1]);
	result = dsigdbetadQ2dW2dz_total_checked(array[0], Q2, array[2], array[3]);
	result *= Q2; // Jacobian
	return log(result);
	}

	double InclusiveDiffractiveCrossSections::dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z){
	double beta=Q2/(Q2+MX2);
	double jacobian=beta*beta/Q2;
	double result= dsigdbetadQ2dW2dz_total_checked(beta, Q2, W2, z);
	return jacobian*result;
	}

	double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z){
	double result = 0;
	//
	// Check if in kinematically allowed region
	double MX2=Kinematics::MX2(Q2, beta, 0);
	if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
	if (mSettings->verboseLevel() > 2)
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
	<< " is outside of kinematically allowed region. Return 0." << endl;
	return 0;
	}
	double csT = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, transverse);
	double csL = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, longitudinal);
	result = csT + csL;


	mCrossSectionRatioLT=csL/csT;
	//
	// Polarization
	//
	if (mRandom->Uniform(result) <= csT)
	mPolarization = transverse;
	else
	mPolarization = longitudinal;
	//
	// Quark Species
	//
	//
	if(mPolarization == transverse){
	bool isChosen=false;
	double cs_rejected=0;
	for(int i=0; i<4; i++){
	double R=mRandom->Uniform(csT-cs_rejected);
	double csTi=mDsigdbetadQ2dWdz_T_total[i];
	if(R <= csTi){
	mQuarkSpecies=i;
	isChosen=true;
	break;
	}
	cs_rejected+=mDsigdbetadQ2dWdz_T_total[i];
	}
	if(!isChosen) mQuarkSpecies=4;
	}
	else{
	bool isChosen=false;
	double cs_rejected=0;
	for(int i=0; i<4; i++){
	double R=mRandom->Uniform(csL-cs_rejected);
	double csLi=mDsigdbetadQ2dWdz_L_total[i];
	if(R <= csLi){
	mQuarkSpecies=i;
	isChosen=true;
	break;
	}
	cs_rejected+=mDsigdbetadQ2dWdz_L_total[i];
	}
	if(!isChosen) mQuarkSpecies=4;
	}
	// Print-out at high verbose levels
	//
	if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(): " << result;
	cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z;
	cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal");
	cout << ')' << endl;
	}
	//
	// Check validity of return value
	//
	if (std::isnan(result)) {
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(): Error, return value = NaN at" << endl;
	cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
	result = 0;
	}
	if (std::isinf(result)) {
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(): Error, return value = inf at" << endl;
	cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
	result = 0;
	}
	if (result < 0) {
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(): Error, negative cross-section at" << endl;
	cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
	result = 0;
	}
	mTotalCS=result;
	return result; //nb/GeV4
	}

	double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
	double result = 0;
	for(int i=0; i<4; i++){
	setQuarkIndex(i);
	double tmpresult = dsigdbetadz_total(beta, Q2, W2, z, pol); //nb
	if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,pol); //nb/GeV4
	if(pol == transverse)
	mDsigdbetadQ2dWdz_T_total[i]=tmpresult;
	if(pol == longitudinal)
	mDsigdbetadQ2dWdz_L_total[i]=tmpresult;
	result+=tmpresult;
	}
	return result; //nb/GeV4
	}

	double InclusiveDiffractiveCrossSections::dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
	double result = 0;
	if(pol==transverse){
	result = dsigmadbetadz_T(beta, Q2, W2, z);
	}
	else if(pol==longitudinal){
	result = dsigmadbetadz_L(beta, Q2, W2, z);
	}
	return result; //nb
	}


	double InclusiveDiffractiveCrossSections::crossSectionOfLastCall(){
	return mTotalCS;
	}

	GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){
	return mPolarization;
	}

	unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){
	return mQuarkSpecies;
	}

	double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){
	// return Settings::quarkMass(mQuarkSpecies); //#TT tmp
	return tt_quarkMass[mQuarkSpecies];
	}

	void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int val){mQuarkIndex=val;}

	DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel(){
	return mDipoleModel;
	}

	double InclusiveDiffractiveCrossSections::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}


	//===================================================
	// QQG Calculations
	//===================================================
	double InclusiveDiffractiveCrossSections::unuranPDF_qqg(const double* array){
	//
	// array is: beta, log(Q2), W2, z
	//
	double result = 0;
	double Q2 = exp(array[1]);

	double z=array[3];
	result = dsigdbetadQ2dW2dz_total_qqg_checked(array[0], Q2, array[2], z);
	result *= Q2; // Jacobian log(Q2)->Q2
	return log(result);
	}

	double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z){
	//
	// Check if in kinematically allowed region
	//
	double MX2=Kinematics::MX2(Q2, beta, 0);
	if(mCheckKinematics && z<beta) {
	if (mSettings->verboseLevel() > 2)
	cout<<"z<beta, beta="<<beta<<", z="<<z<<endl;
	return 0; // For QQG beta < z < 1
	}
	if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
	if (mSettings->verboseLevel() > 2)
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_qqg_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
	<< " is outside of kinematically allowed region. Return 0." << endl;
	return 0;
	}
	- double xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
	- if(mA>1)
	- dipoleModel()->createSigma_ep_LookupTable(xpom);
	double result = dsigdbetadQ2dW2dz_qqg(beta, Q2, W2, z);
	//
	// Quark Species
	//
	//
	bool isChosen=false;
	double cs_rejected=0;
	for(int i=0; i<4; i++){
	double R=mRandom->Uniform(result-cs_rejected);
	double csTi=mDsigdbetadQ2dWdz_T_qqg[i];
	if(R <= csTi){
	mQuarkSpecies=i;
	isChosen=true;
	break;
	}
	cs_rejected+=mDsigdbetadQ2dWdz_T_qqg[i];
	}
	if(!isChosen) mQuarkSpecies=4;
	// Print-out at high verbose levels
	//
	if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_qqg_checked(): " << result;
	cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z<< endl;
	}
	//
	// Check validity of return value
	//
	if (std::isnan(result)) {
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = NaN at" << endl;
	cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
	result = 0;
	}
	if (std::isinf(result)) {
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = inf at" << endl;
	cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
	result = 0;
	}
	if (result < 0) {
	cout << "InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, negative cross-section at" << endl;
	cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
	result = 0;
	}
	mTotalCS_qqg=result;
	return result; //nb/GeV4
	}

	double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z) {
	double result = 0;
	for(int i=0; i<4; i++){
	setQuarkIndex(i);
	// double mf = Settings::quarkMass(mQuarkIndex);
	double mf = tt_quarkMass[i]; //#TT tmp
	double Mqq2=(z/beta-1)*Q2;
	if(Mqq2<4mfmf) continue; //Not enough phase-space for the q and qbar
	double tmpresult = dsigmadbetadz_QQG(beta, Q2, W2, z); //nb
	if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,transverse); //nb/GeV4
	mDsigdbetadQ2dWdz_T_qqg[i]=tmpresult;
	result+=tmpresult;
	}
	return result; //nb/GeV4 (with flux) or nb (without)
	}

	////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	//
	// InclusiveDiffractiveCrossSectionsIntegrals:
	// Inheriting class from InclusiveDiffractiveCrossSections
	//
	////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	InclusiveDiffractiveCrossSectionsIntegrals::InclusiveDiffractiveCrossSectionsIntegrals()
	{
	//Set up the integrals
	//Amplitude0:
	if(Amp0) delete Amp0;
	Amp0=new TF1("Amp0", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_0, 0, 30., 5);

	if(WFAmp0) delete WFAmp0;
	WFAmp0=new ROOT::Math::WrappedTF1(*Amp0);

	GIAmp0.SetFunction(*WFAmp0);
	GIAmp0.SetRelTolerance(1e-4);
	+ GIAmp0.SetRelTolerance(1e-4);

	//Amplitude1:
	if(Amp1) delete Amp1;
	Amp1=new TF1("Amp1", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_1, 0, 30. , 5);

	if(WFAmp1) delete WFAmp1;
	WFAmp1=new ROOT::Math::WrappedTF1(*Amp1);

	GIAmp1.SetFunction(*WFAmp1);
	GIAmp1.SetRelTolerance(1e-4);

	//AmplitudeQQG:
	if(Ampqqg) delete Ampqqg;
	Ampqqg=new TF1("Ampqqg", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_qqg, 0, 30., 4);

	if(WFAmpqqg) delete WFAmpqqg;
	WFAmpqqg=new ROOT::Math::WrappedTF1(*Ampqqg);

	GIAmpqqg.SetFunction(*WFAmpqqg);
	GIAmpqqg.SetRelTolerance(1e-4);

	}

	InclusiveDiffractiveCrossSectionsIntegrals::~InclusiveDiffractiveCrossSectionsIntegrals(){
	if(WFAmp0) delete WFAmp0;
	if(Amp0) delete Amp0;
	if(WFAmp1) delete WFAmp1;
	if(Amp1) delete Amp1;
	if(Ampqqg) delete Ampqqg;
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {

	double MX2 = Q2*(1-beta)/beta;
	double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
	// double mf = Settings::quarkMass(mQuarkIndex);
	double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
	double sqrtArg = 1.-4.mfmf/MX2;
	if (sqrtArg<0){
	if (mSettings->verboseLevel() > 2)
	cout<<"There is no phase-space for z, return 0."<<endl;
	return 0;
	}
	double z0 = (1.-sqrt(sqrtArg))/2.;
	if (z<z0 or z>1-z0){
	if (mSettings->verboseLevel() > 2)
	cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
	return 0;
	}
	double pt2 = z(1-z)Q2-mf*mf;
	if (pt2<0){
	if (mSettings->verboseLevel() > 2)
	- cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
	+ cout<<"Not enough phase-space for quark transverse momenta, return 0. mf="<<mf<<endl;
	return 0;
	}
	if(z(1-z)MX2-mf*mf<0){ //kappa2<0
	if (mSettings->verboseLevel() > 2)
	cout<<"Not enough phase-space for quark production, return 0."<<endl;
	return 0;
	}
	double epsilon2=z(1-z)Q2+mf*mf;
	double ef=quarkCharge[mQuarkIndex];

	double Phi0=Phi_0(beta, xpom, z, Q2, mf);
	double Phi1=Phi_1(beta, xpom, z, Q2, mf);

	double prefactor=NcQ2alpha_em/(4.M_PIbetabeta)efefz*(1-z); //GeV2
	prefactor/=2.; //expanded z-range
	double term1=epsilon2(zz+(1-z)(1-z))Phi1; //GeV2*fm6
	double term2=mfmfPhi0; //GeV2*fm6

	double result=prefactor(term1+term2); //GeV4fm6
	result/=hbarc2*hbarc2; //fm2
	result *= 1e7; //nb
	return result;//nb
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::dsigmadbetadz_L(double beta, double Q2, double W2, double z){

	double MX2 = Q2*(1-beta)/beta;
	double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
	// double mf = Settings::quarkMass(mQuarkIndex);
	double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
	double pt2 = z(1-z)Q2-mf*mf;
	if (pt2<0){
	if (mSettings->verboseLevel() > 2)
	- cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
	+ cout<<"Not enough phase-space for quark transverse momenta, return 0. mf="<<mf<<endl;
	return 0;
	}
	double sqrtArg=1.-4.mfmf/MX2;
	if (sqrtArg<0){
	if (mSettings->verboseLevel() > 2)
	cout<<"There is no phase-space for z, return 0."<<endl;
	return 0;
	}
	double z0=(1.-sqrt(sqrtArg))/2.;
	if (z<z0 or z>1-z0){
	if (mSettings->verboseLevel() > 2)
	cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
	return 0;
	}
	if(z(1-z)MX2-mf*mf<0){ //kappa2<0
	if (mSettings->verboseLevel() > 2)
	cout<<"Not enough phase-space for quark production, return 0."<<endl;
	return 0;
	}
	double ef=quarkCharge[mQuarkIndex];
	double Phi0=Phi_0(beta, xpom, z, Q2, mf); //fm6

	double term=NcQ2Q2alpha_em/(M_PIbeta*beta); //GeV4
	term=efefzzz(1-z)(1-z)(1-z);
	term=Phi0; //fm6GeV4
	term/=2.; //expanded z-range
	double result=term;
	result/=hbarc2*hbarc2; //fm2
	result *= 1e7; //nb
	return result; //nb
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::Phi_0(double beta, double xpom, double z, double Q2, double mf){
	TF1* Phi0=new TF1("Phi0", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_0, 0, 30., 0);
	ROOT::Math::WrappedTF1* WFPhi0=new ROOT::Math::WrappedTF1(*Phi0);
	ROOT::Math::GaussIntegrator GIPhi0;
	GIPhi0.SetFunction(*WFPhi0);
	GIPhi0.SetRelTolerance(1e-4);

	Amp0->SetParameter(0, beta);
	Amp0->SetParameter(1, xpom);
	Amp0->SetParameter(2, z);
	Amp0->SetParameter(3, Q2);
	Amp0->SetParameter(5, mf);

	if(mA>1)
	dipoleModel()->createSigma_ep_LookupTable(xpom);

	double blow=0;
	double result=GIPhi0.IntegralUp(blow); //fm6
	delete Phi0;
	delete WFPhi0;
	return result;
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::Phi_1(double beta, double xpom, double z, double Q2, double mf){
	TF1* Phi1=new TF1("Phi1", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_1, 0, 30., 0);
	ROOT::Math::WrappedTF1* WFPhi1=new ROOT::Math::WrappedTF1(*Phi1);
	ROOT::Math::GaussIntegrator GIPhi1;
	GIPhi1.SetFunction(*WFPhi1);
	GIPhi1.SetRelTolerance(1e-4);

	Amp1->SetParameter(0, beta);
	Amp1->SetParameter(1, xpom);
	Amp1->SetParameter(2, z);
	Amp1->SetParameter(3, Q2);
	Amp1->SetParameter(5, mf);

	if(mA>1)
	dipoleModel()->createSigma_ep_LookupTable(xpom);

	double blow=0;
	double result=GIPhi1.IntegralUp(blow); //fm6
	delete Phi1;
	delete WFPhi1;
	return result;
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_0(double var, double){
	double b=*var;

	double amplitude=Amplitude_0(b); //fm2

	double result=2M_PIbamplitudeamplitude; //fm5
	return result;
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_1(double var, double){
	double b=*var; //fm

	double amplitude=Amplitude_1(b); //fm2
	double result=2M_PIbamplitudeamplitude; //fm5
	return result;
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::Amplitude_0(double b){
	//Gauss Integrator defined in header file and set in constructor

	Amp0->SetParameter(4, b);

	double rlow=1e-3;
	double result=GIAmp0.IntegralUp(rlow); //fm2

	return result;
	}


	double InclusiveDiffractiveCrossSectionsIntegrals::Amplitude_1(double b){
	//Gauss Integrator defined in header file and set in constructor
	Amp1->SetParameter(4, b);

	double rlow=1e-3;
	double result=GIAmp1.IntegralUp(rlow); //fm2

	return result;
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_0(double var, double par){
	double r=*var;
	+ if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
	double beta=par[0];
	double xpom=par[1];
	double z=par[2];
	double Q2=par[3];
	double b=par[4];
	double mf=par[5];
	double MX2=Q2*(1-beta)/beta;
	double kappa2=z(1-z)MX2-mf*mf;
	double kappa=sqrt(kappa2);
	double epsilon2=z(1-z)Q2+mf*mf;
	double epsilon=sqrt(epsilon2); //GeV
	double besselK0=TMath::BesselK0(epsilon*r/hbarc);
	double besselJ0=TMath::BesselJ0(kappa*r/hbarc);
	double dsigmadb2=0;
	if(mA==1)
	dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
	else{
	dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
	}
	return rbesselK0besselJ0*dsigmadb2; //fm
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_1(double var, double par){
	double r=*var; //fm
	+ if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
	double beta=par[0];
	double xpom=par[1];
	double z=par[2];
	double Q2=par[3]; //GeV2
	double b=par[4];
	double mf=par[5];
	double MX2=Q2*(1-beta)/beta; //GeV2
	double kappa2=z(1-z)MX2-mf*mf; //GeV2
	double kappa=sqrt(kappa2); //GeV
	double epsilon2=z(1-z)Q2+mf*mf;
	double epsilon=sqrt(epsilon2); //GeV
	double besselK1=TMath::BesselK1(epsilon*r/hbarc);
	double besselJ1=TMath::BesselJ1(kappa*r/hbarc);
	double dsigmadb2=0; //GeV0
	if(mA==1)
	dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
	else
	dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
	return rbesselK1besselJ1*dsigmadb2; //fm
	}

	//===================================================
	// QQG Calculations
	//===================================================


	double InclusiveDiffractiveCrossSectionsIntegrals::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
	double MX2 = Q2*(1-beta)/beta;
	double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));

	double Phi=Phi_qqg(xpom, ztilde, Q2);

	double alpha_S=0.15;
	double ef=quarkCharge[mQuarkIndex];

	double betafactor=(1-beta/ztilde)*
	(1-beta/ztilde)+beta/ztilde*beta/ztilde;
	double result=alpha_Salpha_em/(2M_PIM_PIQ2)betafactorefefPhi; //GeV-2
	result*=hbarc2; //fm2
	result *= 1e7; //nb
	return result; //nb
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::Phi_qqg(double xpom, double z, double Q2){

	Ampqqg->SetParameter(0, xpom);
	Ampqqg->SetParameter(1, z);

	if(mA>1)
	dipoleModel()->createSigma_ep_LookupTable(xpom);

	mQ2=Q2;

	ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_qqg, 2);
	ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
	ig.SetFunction(wf);
	ig.SetAbsTolerance(0.);
	ig.SetRelTolerance(1e-4);
	double bRange=3. * dipoleModel()->nucleus()->radius();
	-
	double lo[2]={0.,1e-4}; //b, k2
	double hi[2]={bRange, Q2}; //b, k2
	double result = ig.Integral(lo, hi)/hbarc; //GeV0
	return result; //GeV0
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_qqg(const double* var) {
	double b=var[0];
	double k2=var[1];
	double Q2=mQ2;

	Ampqqg->SetParameter(2, sqrt(k2));
	Ampqqg->SetParameter(3, b);

	double rlow=1e-3;
	double amp_qqg=GIAmpqqg.IntegralUp(rlow); //fm2
	amp_qqg/=hbarc2; //GeV-2

	double result = (2M_PI)b/hbarck2k2log(Q2/k2)amp_qqg*amp_qqg; //GeV-1
	return result; //GeV-1
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_qqg(double var, double par){
	double r=*var; //fm
	if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
	double xpom=par[0];
	double z=par[1];
	double k=par[2]; //GeV
	double b=par[3];

	double besselK2=ROOT::Math::cyl_bessel_k(2, sqrt(z)kr/hbarc);
	double besselJ2=ROOT::Math::cyl_bessel_j(2, sqrt(1-z)kr/hbarc);
	double dsigmadb2=0; //GeV0
	if(mA==1)
	dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
	else
	dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
	double term=1-0.5*dsigmadb2;
	double dsigmadb2tilde=2(1-termterm);
	double result=rbesselK2besselJ2*dsigmadb2tilde; //fm
	return result;
	}

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:01 PM (23 h, 1 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4242370
Default Alt Text: (107 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions