No OneTemporary
Actions

Size

112 KB

Subscribers

None

View Options

	Index: trunk/src/TableGeneratorSettings.cpp
	===================================================================
	--- trunk/src/TableGeneratorSettings.cpp (revision 538)
	+++ trunk/src/TableGeneratorSettings.cpp (revision 539)
	@@ -1,242 +1,225 @@
	//==============================================================================
	// TableGeneratorSettings.cpp
	//
	// Copyright (C) 2010-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$
	// $Author$
	//==============================================================================
	#include "Settings.h"
	#include "TableGeneratorSettings.h"
	#include "Constants.h"
	#include <cmath>
	#include <ctime>
	#include <vector>
	#include <stdlib.h>

	using namespace std;

	TableGeneratorSettings* TableGeneratorSettings::mInstance = 0; // initialize static

	TableGeneratorSettings* TableGeneratorSettings::instance()
	{
	if (mInstance == 0)
	mInstance = new TableGeneratorSettings;
	return mInstance;
	}

	TableGeneratorSettings::TableGeneratorSettings()
	{
	//
	// Register all the parameters that can be defined
	// via a runcard.
	// Arguments for registerParameter():
	// 1. pointer to data memeber
	// 2. text string to be used in the runcard
	// 3. default parameter set
	//
	registerParameter(&mBSatLookupPath, "bSatLookupPath", string("./"));

	registerParameter(&mTmin, "tmin", -2.);
	registerParameter(&mTmax, "tmax", 0.);

	registerParameter(&mXmin, "xmin", 1e-9); //UPC
	registerParameter(&mXmax, "xmax", 3e-2); //UPC

	- registerParameter(&mBetamin, "betamin", 1e-9); //inc.diff.
	- registerParameter(&mBetamax, "betamax", 0.9); //inc.diff.
	- registerParameter(&mZmin, "zmin", 1e-9); //inc.diff.
	- registerParameter(&mZmax, "zmax", 0.9); //inc.diff.
	-
	registerParameter(&mQ2bins, "Q2bins", static_cast<unsigned int>(1));
	registerParameter(&mW2bins, "W2bins", static_cast<unsigned int>(1));
	registerParameter(&mTbins, "tbins", static_cast<unsigned int>(1));
	registerParameter(&mXbins, "xbins", static_cast<unsigned int>(1)); //UPC
	registerParameter(&mBetabins, "betabins", static_cast<unsigned int>(1)); //inc.diff.
	registerParameter(&mZbins, "zbins", static_cast<unsigned int>(1)); //inc.diff.

	registerParameter(&mNumberOfConfigurations, "numberOfConfigurations", static_cast<unsigned int>(1000));
	vector<double> vec;
	registerParameter(&mDipoleModelCustomParameters, "dipoleModelCustomParameters", vec);
	registerParameter(&mUseBackupFile, "useBackupFile", false);
	registerParameter(&mStartingBinFromBackup, "startingBinFromBackup", 0);

	registerParameter(&mStartBin, "startBin", -1);
	registerParameter(&mEndBin, "endBin", -1);

	registerParameter(&mModesToCalculate, "modesToCalculate", 0);

	registerParameter(&mPriority, "priority", 0);

	registerParameter(&mHasSubstructure, "hasSubstructure", false);

	registerParameter(&mFractionOfBinsToFill, "fractionOfBinsToFill", 1.);

	}

	void TableGeneratorSettings::consolidateSettings() // called after runcard is read
	{
	//
	// Kinematic limits
	//
	if (mQ2min>=mQ2max && !mUPC) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, Q2min >= Q2max. Stopping" << endl;
	exit(1);
	}
	if (mWmin>=mWmax && !mUPC) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, Wmin >= Wmax. Stopping" << endl;
	exit(1);
	}
	if (mTmin>=mTmax) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, tmin >= tmax. Stopping" << endl;
	exit(1);
	}
	if (mTmin>0. \|\| mTmax >0.) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, t must be negative, please change t-limits. Stopping" << endl;
	exit(1);
	}
	if (mXmin>=mXmax && mUPC) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, xmin >= xmax. Stopping" << endl;
	exit(1);
	}
	if ((mXmin>=1 \|\| mXmin<=0 \|\| mXmax>=1 \|\| mXmax<=0) && mUPC) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, xmin or xmax out of range. Stopping" << endl;
	exit(1);
	}
	if (mXmax>0.01 && mUPC){
	cout << "TableGeneratorSettings::consolidateSettings(): Warning, xmax>1e-2, model may be unreliable." << endl;
	}
	if (mA==1 && !mHasSubstructure) mNumberOfConfigurations = 1;

	if (!mUseBackupFile) mStartingBinFromBackup = 0;

	if (!mUPC)
	mXbins=1;
	else {
	mQ2bins=1;
	mW2bins=1;
	}
	if (mStartBin >= 0 && mEndBin >= 0 && mStartBin>mEndBin) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, endBin < startBin : " << mEndBin << " <= " << mStartBin << "! Stopping." << endl;
	exit(1);
	}
	if ( mStartBin < 0 ) mStartBin=0;
	if ( mStartBin >= signed(mQ2binsmW2binsmTbins*mXbins) ) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, starting bin >= table! Stopping." << endl;
	exit(1);
	}
	if ( mEndBin > signed(mQ2binsmW2binsmTbins*mXbins) \|\| mEndBin < 0) {
	cout << "TableGeneratorSettings::consolidateSettings(): endBin is set to table size=" << mQ2binsmW2binsmTbins << endl;
	mEndBin=mQ2binsmW2binsmTbinsmXbinsmBetabins*mZbins;
	}
	if ( mModesToCalculate < 0 \|\| mModesToCalculate > 2 ) {
	cout << "TableGeneratorSettings::consolidateSettings(): Error, modesToCalculate can only take values 0, 1, or 2; not "
	<< mModesToCalculate << endl;
	exit(1);
	}

	//
	// Make sure the W range is allowed
	//
	double VMMass=lookupPDG(mVectorMesonId)->Mass();
	double W2min=VMMass*VMMass+protonMass2;
	double Wmin=sqrt(W2min);
	if (mWmin<Wmin){
	mWmin=Wmin;
	cout << "TableGeneratorSettings::consolidateSettings(): Warning, Wmin is smaller than allowed value." << endl;
	cout << " It has been changed to Wmin=" << mWmin << endl;
	}
	}

	//
	// Access functions
	//

	void TableGeneratorSettings::setTmax(double val){ mTmax=val; }
	double TableGeneratorSettings::tmax() const { return mTmax; }

	void TableGeneratorSettings::setTmin(double val){ mTmin=val; }
	double TableGeneratorSettings::tmin() const { return mTmin; }

	void TableGeneratorSettings::setXmax(double val){ mXmax=val; }
	double TableGeneratorSettings::xmax() const { return mXmax; }

	void TableGeneratorSettings::setXmin(double val){ mXmin=val; }
	double TableGeneratorSettings::xmin() const { return mXmin; }

	-void TableGeneratorSettings::setZmax(double val){ mZmax=val; }
	-double TableGeneratorSettings::zmax() const { return mZmax; }
	-
	-void TableGeneratorSettings::setZmin(double val){ mZmin=val; }
	-double TableGeneratorSettings::zmin() const { return mZmin; }
	-
	-void TableGeneratorSettings::setBetamax(double val){ mBetamax=val; }
	-double TableGeneratorSettings::betamax() const { return mBetamax; }
	-
	-void TableGeneratorSettings::setBetamin(double val){ mBetamin=val; }
	-double TableGeneratorSettings::betamin() const { return mBetamin; }
	-
	void TableGeneratorSettings::setQ2bins(unsigned int val){ mQ2bins=val; }
	unsigned int TableGeneratorSettings::Q2bins() const { return mQ2bins; }

	void TableGeneratorSettings::setW2bins(unsigned int val){ mW2bins=val; }
	unsigned int TableGeneratorSettings::W2bins() const { return mW2bins; }

	void TableGeneratorSettings::setTbins(unsigned int val){ mTbins=val; }
	unsigned int TableGeneratorSettings::tbins() const { return mTbins; }

	void TableGeneratorSettings::setXbins(unsigned int val){ mXbins=val; }
	unsigned int TableGeneratorSettings::xbins() const { return mXbins; }

	void TableGeneratorSettings::setBetabins(unsigned int val){ mBetabins=val; }
	unsigned int TableGeneratorSettings::betabins() const { return mBetabins; }

	void TableGeneratorSettings::setZbins(unsigned int val){ mZbins=val; }
	unsigned int TableGeneratorSettings::zbins() const { return mZbins; }

	void TableGeneratorSettings::setBSatLookupPath(string val){ mBSatLookupPath = val; }
	string TableGeneratorSettings::bSatLookupPath() const { return mBSatLookupPath; }

	void TableGeneratorSettings::setNumberOfConfigurations(unsigned int val){ mNumberOfConfigurations=val; }
	unsigned int TableGeneratorSettings::numberOfConfigurations() const { return mNumberOfConfigurations; }

	vector<double> TableGeneratorSettings::dipoleModelCustomParameters() const {return mDipoleModelCustomParameters;}

	void TableGeneratorSettings::setUseBackupFile(bool val){ mUseBackupFile = val; }
	bool TableGeneratorSettings::useBackupFile() const { return mUseBackupFile; }

	void TableGeneratorSettings::setStartingBinFromBackup(int val){ mStartingBinFromBackup = val; }
	int TableGeneratorSettings::startingBinFromBackup() const { return mStartingBinFromBackup; }

	void TableGeneratorSettings::setStartBin(int val){ mStartBin = val; }
	int TableGeneratorSettings::startBin() const { return mStartBin; }

	void TableGeneratorSettings::setEndBin(int val){ mEndBin = val; }
	int TableGeneratorSettings::endBin() const{ return mEndBin; }

	void TableGeneratorSettings::setModesToCalculate(int val){ mModesToCalculate = val; }
	int TableGeneratorSettings::modesToCalculate() const { return mModesToCalculate; }

	void TableGeneratorSettings::setPriority(unsigned char val){ mPriority = val; }
	unsigned char TableGeneratorSettings::priority() const { return mPriority; }

	void TableGeneratorSettings::setHasSubstructure(bool val){ mHasSubstructure = val; }
	bool TableGeneratorSettings::hasSubstructure() const { return mHasSubstructure; }

	void TableGeneratorSettings::setFractionOfBinsToFill(double val){ mFractionOfBinsToFill = val; }
	double TableGeneratorSettings::fractionOfBinsToFill() const { return mFractionOfBinsToFill; }

	Index: trunk/src/InclusiveDiffractiveCrossSections.h
	===================================================================
	--- trunk/src/InclusiveDiffractiveCrossSections.h (revision 538)
	+++ trunk/src/InclusiveDiffractiveCrossSections.h (revision 539)
	@@ -1,170 +1,145 @@
	//==============================================================================
	// InclusiveDiffractiveCrossSection.h
	//
	// 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 $
	//==============================================================================
	#ifndef InclusiveDiffractiveCrossSections_h
	#define InclusiveDiffractiveCrossSections_h
	#include "AlphaStrong.h"
	#include "TH1.h"
	#include "TMath.h"
	#include "TF1.h"
	#include "Math/WrappedTF1.h"
	#include "Math/GaussIntegrator.h"
	#include "DipoleModel.h"
	#include "PhotonFlux.h"
	#include "Enumerations.h"
	#include "Constants.h"
	#include "Table.h"
	#include "THn.h"

	class DipoleModel;
	class DipoleModelParameters;

	class InclusiveDiffractiveCrossSections{
	public:
	InclusiveDiffractiveCrossSections();
	virtual ~InclusiveDiffractiveCrossSections();
	- virtual double operator()(double beta, double Q2, double W2, double z);
	- virtual double operator()(const double* array);
	-
	- virtual DipoleModel* dipoleModel();
	- virtual double unuranPDF(const double*);
	- virtual double unuranPDF_qqg(const double*);
	- virtual void setCheckKinematics(bool);
	- virtual void setFockState(FockState);
	-
	- virtual GammaPolarization polarizationOfLastCall() ;
	- virtual double crossSectionOfLastCall();
	- virtual unsigned int quarkSpeciesOfLastCall();
	- virtual double quarkMassOfLastCall();
	- virtual double crossSectionRatioLTOfLastCall() const;
	-
	- virtual void setQuarkIndex(unsigned int);
	-
	- virtual double dsigmadbetadz_T(double, double, double, double) ;
	- virtual double dsigmadbetadz_L(double, double, double, double) ;
	-
	- virtual double dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z);
	- virtual double dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z);
	-
	- virtual double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;
	-
	- virtual double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
	-
	- virtual double* dsigdbetadQ2dWdz_T_total();
	- virtual double* dsigdbetadQ2dWdz_L_total();
	- virtual double* dsigdbetadQ2dWdz_T_qqg();
	-
	-};
	-
	-class InclusiveDiffractiveCrossSectionsIntegrals: public InclusiveDiffractiveCrossSections{
	-
	-public:
	- InclusiveDiffractiveCrossSectionsIntegrals();
	- ~InclusiveDiffractiveCrossSectionsIntegrals();
	-
	double operator()(double beta, double Q2, double W2, double z);
	double operator()(const double* array);
	-
	+
	+ DipoleModel* dipoleModel();
	double unuranPDF(const double*);
	double unuranPDF_qqg(const double*);
	+ void setCheckKinematics(bool);
	+ void setFockState(FockState);

	GammaPolarization polarizationOfLastCall() ;
	double crossSectionOfLastCall();
	-
	unsigned int quarkSpeciesOfLastCall();
	double quarkMassOfLastCall();
	double crossSectionRatioLTOfLastCall() const;
	- void setCheckKinematics(bool);
	- void setFockState(FockState);
	- FockState getFockState();
	-
	- double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;

	- double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
	+ void setQuarkIndex(unsigned int);
	+
	+ virtual double dsigmadbetadz_T(double, double, double, double) ;
	+ virtual double dsigmadbetadz_L(double, double, double, double) ;
	+ virtual double dsigmadbetadz_QQG(double, double, double, double) ;
	+

	-// InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections&);
	-// InclusiveDiffractiveCrossSections& operator=(const InclusiveDiffractiveCrossSections&);
	double dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z);
	double dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z);
	-
	- DipoleModel* dipoleModel();

	- double dsigmadbetadz_T(double, double, double, double) ;
	- double dsigmadbetadz_L(double, double, double, double) ;
	- double dsigmadbetadz_QQG(double, double, double, double) ;
	+ double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;

	- double* dsigdbetadQ2dWdz_T_total(){return mDsigdbetadQ2dWdz_L_total;}
	- double* dsigdbetadQ2dWdz_L_total(){return mDsigdbetadQ2dWdz_T_total;}
	- double* dsigdbetadQ2dWdz_T_qqg(){ return mDsigdbetadQ2dWdz_T_qqg;};
	+ double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);

	- void setQuarkIndex(unsigned int);
	+ FockState getFockState();
	+
	+ virtual double* dsigdbetadQ2dWdz_T_total();
	+ virtual double* dsigdbetadQ2dWdz_L_total();
	+ virtual double* dsigdbetadQ2dWdz_T_qqg();

	protected:
	double dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z);
	double dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol);

	-private:
	+
	+//private:
	TRandom3* mRandom;
	GammaPolarization mPolarization;
	PhotonFlux mPhotonFlux;
	EventGeneratorSettings* mSettings;
	unsigned int mQuarkSpecies;
	double mTotalCS, mTotalCS_qqg;
	FockState mFockState;

	double mS, mQ2;
	-
	+
	bool mCheckKinematics;

	+ unsigned int mA;
	+ unsigned int mQuarkIndex;
	+ double mCrossSectionRatioLT;
	+
	+ double mDsigdbetadQ2dWdz_L_total[5];
	+ double mDsigdbetadQ2dWdz_T_total[5];
	+ double mDsigdbetadQ2dWdz_T_qqg[5];
	+
	+ DipoleModel* mDipoleModel;
	+ DipoleModelParameterSet mDipoleModelParameterSet;
	+
	+};
	+
	+class InclusiveDiffractiveCrossSectionsIntegrals: public InclusiveDiffractiveCrossSections{
	+
	+public:
	+ InclusiveDiffractiveCrossSectionsIntegrals();
	+ ~InclusiveDiffractiveCrossSectionsIntegrals();
	+
	+ double dsigmadbetadz_T(double, double, double, double) ;
	+ double dsigmadbetadz_L(double, double, double, double) ;
	+ double dsigmadbetadz_QQG(double, double, double, double) ;
	+
	+protected:
	+
	+private:
	+
	double Phi_0(double, double, double, double, double);
	double Phi_1(double, double, double, double, double);
	double Phi_qqg(double, double, double);
	double uiPhi_0(double, double);
	double uiPhi_1(double, double);
	double Amplitude_0(double);
	double Amplitude_1(double);
	double uiAmplitude_0(double, double);
	double uiAmplitude_1(double, double);


	double uiPhi_qqg(const double*);
	double uiAmplitude_qqg(double, double);
	-
	- unsigned int mA;
	- unsigned int mQuarkIndex;
	- double mCrossSectionRatioLT;
	-
	- double mDsigdbetadQ2dWdz_L_total[5];
	- double mDsigdbetadQ2dWdz_T_total[5];
	- double mDsigdbetadQ2dWdz_T_qqg[5];

	- DipoleModel* mDipoleModel;
	- DipoleModelParameterSet mDipoleModelParameterSet;
	-
	TF1 Amp0, Amp1, *Ampqqg;
	ROOT::Math::WrappedTF1 WFAmp0, WFAmp1, *WFAmpqqg;
	ROOT::Math::GaussIntegrator GIAmp0, GIAmp1, GIAmpqqg;
	-
	-
	+
	};

	#endif
	Index: trunk/src/Settings.cpp
	===================================================================
	--- trunk/src/Settings.cpp (revision 538)
	+++ trunk/src/Settings.cpp (revision 539)
	@@ -1,618 +1,620 @@
	//==============================================================================
	// 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);
	}

	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::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 538)
	+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 539)
	@@ -1,579 +1,271 @@
	#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(){

	- mSettings = EventGeneratorSettings::instance();
	- mRandom = mSettings->randomGenerator();
	- mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
	- mPhotonFlux.setS(mS);
	- mCheckKinematics = true;
	- mCrossSectionRatioLT = 0;
	-
	- mQ2=0;
	-
	- mA=mSettings->A();
	-
	- 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);
	- }
	-
	- //
	- // Set up the lookup tables:
	- //
	- TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
	- TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
	- TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
	- TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
	-
	- TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
	- TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
	- TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
	- TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
	-
	- TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
	- TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
	- TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
	- TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
	+ 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");

	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);

	- // if (!hist4D) {
	- // std::cerr << "Error: Histogram 'hist4D' not found!" << std::endl;
	- // file->Close();
	- // return -1;
	- // }
	}
	InclusiveDiffractionCrossSectionsFromTables::~InclusiveDiffractionCrossSectionsFromTables(){
	}

	-void InclusiveDiffractionCrossSectionsFromTables::setCheckKinematics(bool val) {mCheckKinematics = val;}
	-
	-void InclusiveDiffractionCrossSectionsFromTables::setFockState(FockState val){
	- mFockState=val;
	-}
	-
	-FockState InclusiveDiffractionCrossSectionsFromTables::getFockState(){
	- return mFockState;
	-}
	-
	-double InclusiveDiffractionCrossSectionsFromTables::operator()(double MX2, double Q2, double W2, double z){
	- return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
	-}
	-
	-double InclusiveDiffractionCrossSectionsFromTables::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<<"InclusiveDiffractiveCrossSections::operator(): Fockstate is invalid, stopping"<<endl;
	- exit(0);
	- }
	- return result;
	-}
	-
	-double InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::crossSectionOfLastCall(){
	- return mTotalCS;
	-}
	-
	-GammaPolarization InclusiveDiffractionCrossSectionsFromTables::polarizationOfLastCall(){
	- return mPolarization;
	-}
	-
	-unsigned int InclusiveDiffractionCrossSectionsFromTables::quarkSpeciesOfLastCall(){
	- return mQuarkSpecies;
	-}
	-
	-double InclusiveDiffractionCrossSectionsFromTables::quarkMassOfLastCall(){
	- // return Settings::quarkMass(mQuarkSpecies); //#TT tmp
	- return tt_quarkMass[mQuarkSpecies];
	-}
	-
	-void InclusiveDiffractionCrossSectionsFromTables::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
	-
	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);
	+ 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;
	}

	-DipoleModel* InclusiveDiffractionCrossSectionsFromTables::dipoleModel(){
	- return mDipoleModel;
	-}
	-
	-double InclusiveDiffractionCrossSectionsFromTables::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
	-
	//===================================================
	// QQG Calculations
	//===================================================

	-double InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::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 << "InclusiveDiffractiveCrossSections::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 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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 << "InclusiveDiffractiveCrossSections::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 InclusiveDiffractionCrossSectionsFromTables::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
	-}
	-
	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){
	- THnF* hist4D=0;// = dynamic_cast<THnF*>(file->Get("hist4D"));
	+// double scaleFactor=(0.9-0.1)/10;
	+ beta = beta*0.8/0.9+0.1;
	+ THnF* hist4D=0;
	if(pol==transverse and fock==QQ){
	hist4D=tableQQ_T.at(quark);
	}
	else if(pol==longitudinal and fock==QQ){
	hist4D=tableQQ_L.at(quark);
	}
	else if(fock==QQG){
	hist4D=tableQQG.at(quark);
	}
	else{
	- cout<<"There is nothing to interpolate, ending"<<endl;
	+ 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);
	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<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(gridBeta.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;

	// 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;
	+// 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;
	-
	-
	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);
	count++;
	}
	}
	}
	}
	+// PR(count);
	// std::cout<<"done"<<std::endl;

	// 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);
	}

	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);
	}
	return gridtype;
	}

	Index: trunk/src/TableGeneratorSettings.h
	===================================================================
	--- trunk/src/TableGeneratorSettings.h (revision 538)
	+++ trunk/src/TableGeneratorSettings.h (revision 539)
	@@ -1,154 +1,139 @@
	//==============================================================================
	// TableGeneratorSettings.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: Tobias Toll
	// Last update:
	// $Date$
	// $Author$
	//==============================================================================
	//
	// Singleton class
	//
	//==============================================================================
	#ifndef TableGeneratorSettings_h
	#define TableGeneratorSettings_h
	#include "Settings.h"
	#include "Enumerations.h"

	using namespace std;

	class TableGeneratorSettings : public Settings {
	public:
	static TableGeneratorSettings* instance();

	void setTmin(double);
	double tmin() const;

	void setTmax(double);
	double tmax() const;

	void setXmin(double);
	double xmin() const;

	void setXmax(double);
	double xmax() const;

	- void setBetamin(double);
	- double betamin() const;
	-
	- void setBetamax(double);
	- double betamax() const;
	-
	- void setZmin(double);
	- double zmin() const;
	-
	- void setZmax(double);
	- double zmax() const;
	-
	void setQ2bins(unsigned int);
	unsigned int Q2bins() const;

	void setW2bins(unsigned int);
	unsigned int W2bins() const;

	void setTbins(unsigned int);
	unsigned int tbins() const;

	void setXbins(unsigned int);
	unsigned int xbins() const;

	void setBetabins(unsigned int);
	unsigned int betabins() const;

	void setZbins(unsigned int);
	unsigned int zbins() const;

	string bSatLookupPath() const;
	void setBSatLookupPath(string);

	unsigned int numberOfConfigurations() const;
	void setNumberOfConfigurations(unsigned int);

	vector<double> dipoleModelCustomParameters() const;

	bool useBackupFile() const;
	void setUseBackupFile(bool);

	int startingBinFromBackup() const;
	void setStartingBinFromBackup(int);

	int startBin() const;
	void setStartBin(int);

	int endBin() const;
	void setEndBin(int);

	int modesToCalculate() const;
	void setModesToCalculate(int);

	unsigned char priority() const;
	void setPriority(unsigned char);

	bool hasSubstructure() const;
	void setHasSubstructure(bool);

	double fractionOfBinsToFill() const;
	void setFractionOfBinsToFill(double);

	void consolidateSettings();

	private:
	TableGeneratorSettings();

	private:
	static TableGeneratorSettings* mInstance;

	private:
	unsigned int mNumberOfConfigurations;

	vector<double> mDipoleModelCustomParameters; // developer

	bool mUseBackupFile;

	int mStartingBinFromBackup;
	int mStartBin, mEndBin;
	int mModesToCalculate;

	double mTmin;
	double mTmax;
	double mXmin;
	double mXmax;
	-
	- double mBetamin, mBetamax;
	- double mZmin, mZmax;
	-
	+
	unsigned int mQ2bins;
	unsigned int mW2bins;
	unsigned int mTbins;
	unsigned int mXbins;
	unsigned int mBetabins;
	unsigned int mZbins;

	string mBSatLookupPath;

	int mPriority;

	bool mHasSubstructure;
	double mFractionOfBinsToFill;
	};

	#endif
	Index: trunk/src/InclusiveDiffractiveCrossSections.cpp
	===================================================================
	--- trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 538)
	+++ trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 539)
	@@ -1,719 +1,729 @@
	//==============================================================================
	// 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(){}
	-
	-InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){}
	-
	-double InclusiveDiffractiveCrossSections::operator()(double, double, double, double){return 0;}
	-double InclusiveDiffractiveCrossSections::operator()(const double*){ return 0;}
	-DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel(){return 0;}
	-double InclusiveDiffractiveCrossSections::unuranPDF(const double*){return 0;}
	-double InclusiveDiffractiveCrossSections::unuranPDF_qqg(const double*){return 0;}
	-void InclusiveDiffractiveCrossSections::setCheckKinematics(bool){}
	-void InclusiveDiffractiveCrossSections::setFockState(FockState){}
	-GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){return transverse;}
	-double InclusiveDiffractiveCrossSections::crossSectionOfLastCall(){return 0;}
	-unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){return 0;}
	-double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){return 0;}
	-double InclusiveDiffractiveCrossSections::crossSectionRatioLTOfLastCall() const{return 0;}
	-void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int){}
	-double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double, double, double, double){ return 0; }
	-double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double, double, double, double){ return 0; }
	-double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(double, double, double, double){return 0;}
	-double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(double, double, double, double){return 0;}
	-
	-double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double, double, double, double, GammaPolarization){return 0;}
	-
	-double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_qqg(double, double, double, double){return 0;}
	-
	-double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_total(){return 0;}
	-double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_L_total(){return 0;}
	-double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_qqg(){return 0;}
	-
	-
	-InclusiveDiffractiveCrossSectionsIntegrals::InclusiveDiffractiveCrossSectionsIntegrals()
	-{
	+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) {
	+ 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();
	- //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);
	-
	- //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;
	-}
	+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;};

	-void InclusiveDiffractiveCrossSectionsIntegrals::setCheckKinematics(bool val) {mCheckKinematics = val;}
	+//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 InclusiveDiffractiveCrossSectionsIntegrals::setFockState(FockState val){
	+
	+void InclusiveDiffractiveCrossSections::setCheckKinematics(bool val) {mCheckKinematics = val;}
	+
	+void InclusiveDiffractiveCrossSections::setFockState(FockState val){
	mFockState=val;
	}

	-FockState InclusiveDiffractiveCrossSectionsIntegrals::getFockState(){
	+FockState InclusiveDiffractiveCrossSections::getFockState(){
	return mFockState;
	}
	-double InclusiveDiffractiveCrossSectionsIntegrals::operator()(double MX2, double Q2, double W2, double z){
	+
	+double InclusiveDiffractiveCrossSections::operator()(double MX2, double Q2, double W2, double z){
	return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
	}

	-double InclusiveDiffractiveCrossSectionsIntegrals::operator()(const double* array){
	+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 InclusiveDiffractiveCrossSectionsIntegrals::unuranPDF(const double* array){
	+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 InclusiveDiffractiveCrossSectionsIntegrals::dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z){
	+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 InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z){
	+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 xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
	- if(mA>1)
	- dipoleModel()->createSigma_ep_LookupTable(xpom);
	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 InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
	+double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
	double result = 0;
	- for(int i=0; i<5; i++){
	+ 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 InclusiveDiffractiveCrossSectionsIntegrals::crossSectionOfLastCall(){
	+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 InclusiveDiffractiveCrossSectionsIntegrals::polarizationOfLastCall(){
	+GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){
	return mPolarization;
	}

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

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

	-void InclusiveDiffractiveCrossSectionsIntegrals::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
	+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);
	+
	+ //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;
	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;
	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;
	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
	+ else{
	dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
	+ }
	return rbesselK0besselJ0*dsigmadb2; //fm
	}

	double InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_1(double var, double par){
	double r=*var; //fm
	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
	}

	-DipoleModel* InclusiveDiffractiveCrossSectionsIntegrals::dipoleModel(){
	- return mDipoleModel;
	-}
	-
	-double InclusiveDiffractiveCrossSectionsIntegrals::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
	-
	-
	//===================================================
	// QQG Calculations
	//===================================================
	-double InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z) {
	- double result = 0;
	- for(int i=0; i<5; 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
	-}

	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.14;
	+ 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]={5., Q2}; //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;
	}
	Index: trunk/src/InclusiveDiffractiveCrossSectionsFromTables.h
	===================================================================
	--- trunk/src/InclusiveDiffractiveCrossSectionsFromTables.h (revision 538)
	+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.h (revision 539)
	@@ -1,95 +1,75 @@
	#ifndef InclusiveDiffractiveCrossSectionsFromTables_h
	#define InclusiveDiffractiveCrossSectionsFromTables_h
	#include "AlphaStrong.h"
	#include "TH1.h"
	#include "TMath.h"
	#include "TF1.h"
	#include "Math/WrappedTF1.h"
	#include "Math/GaussIntegrator.h"
	#include "DipoleModel.h"
	#include "PhotonFlux.h"
	#include "Enumerations.h"
	#include "Constants.h"
	#include "Table.h"
	#include "THn.h"
	#include <vector>
	#include "InclusiveDiffractiveCrossSections.h"

	using namespace std;

	class InclusiveDiffractionCrossSectionsFromTables : public InclusiveDiffractiveCrossSections {

	public:
	InclusiveDiffractionCrossSectionsFromTables();
	~InclusiveDiffractionCrossSectionsFromTables();

	double operator()(double beta, double Q2, double W2, double z);
	double operator()(const double* array);

	double unuranPDF(const double*);
	double unuranPDF_qqg(const double*);

	GammaPolarization polarizationOfLastCall() ;
	double crossSectionOfLastCall();

	unsigned int quarkSpeciesOfLastCall();
	double quarkMassOfLastCall();
	double crossSectionRatioLTOfLastCall() const;
	void setCheckKinematics(bool);
	void setFockState(FockState);
	FockState getFockState();

	double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;

	double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);

	double dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z);
	double dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z);

	double dsigmadbetadz_T(double, double, double, double) ;
	double dsigmadbetadz_L(double, double, double, double) ;
	double dsigmadbetadz_QQG(double, double, double, double) ;

	double mDsigdbetadQ2dWdz_L_total[5];
	double mDsigdbetadQ2dWdz_T_total[5];
	double mDsigdbetadQ2dWdz_T_qqg[5];

	void setQuarkIndex(unsigned int);

	double interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark);

	vector<double> mgrid(THnF* hist,int axis, int len);

	vector<THnF*> tableQQ_T;
	vector<THnF*> tableQQ_L;
	vector<THnF*> tableQQG;

	DipoleModel* dipoleModel();

	protected:
	double dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z);
	double dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol);

	private:
	- TRandom3* mRandom;
	- GammaPolarization mPolarization;
	- PhotonFlux mPhotonFlux;
	- EventGeneratorSettings* mSettings;
	- unsigned int mQuarkSpecies;
	- double mTotalCS, mTotalCS_qqg;
	- FockState mFockState;
	-
	- double mS, mQ2;
	-
	- bool mCheckKinematics;
	-
	- unsigned int mA;
	- unsigned int mQuarkIndex;
	- double mCrossSectionRatioLT;
	-
	- DipoleModel* mDipoleModel;
	- DipoleModelParameterSet mDipoleModelParameterSet;
	-
	-
	};
	#endif

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 3:09 PM (1 d, 51 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4014372
Default Alt Text: (112 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions