No OneTemporary
Actions

Size

45 KB

Subscribers

None

View Options

	Index: branches/fitting/src/DipoleModelParameters.cpp
	===================================================================
	--- branches/fitting/src/DipoleModelParameters.cpp (revision 385)
	+++ branches/fitting/src/DipoleModelParameters.cpp (revision 386)
	@@ -1,493 +1,493 @@
	//==============================================================================
	// DipoleModelParameters.cpp
	//
	// Copyright (C) 2016-2018 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
	// $Date$
	// $Author: ullrich $
	//==============================================================================
	#include "DipoleModelParameters.h"
	#include "Settings.h"
	#include "TableGeneratorSettings.h"
	#include <iomanip>

	using namespace std;

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


	DipoleModelParameters::DipoleModelParameters(Settings* settings)
	{
	mDipoleModelType = settings->dipoleModelType();
	mDipoleModelParameterSetName = settings->dipoleModelParameterSetName();
	mDipoleModelParameterSet = settings->dipoleModelParameterSet();
	setupParameters();
	}

	DipoleModelParameters::DipoleModelParameters(DipoleModelType mtype, DipoleModelParameterSet pset) :
	mDipoleModelType(mtype), mDipoleModelParameterSet(pset)
	{
	setupParameters();
	}

	void DipoleModelParameters::setDipoleModelType(DipoleModelType val)
	{
	mDipoleModelType = val;
	setupParameters();
	}

	void DipoleModelParameters::setDipoleModelParameterSet(DipoleModelParameterSet val)
	{
	mDipoleModelParameterSet = val;
	setupParameters();
	}

	DipoleModelType DipoleModelParameters::dipoleModelType() const {return mDipoleModelType;}

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

	void DipoleModelParameters::setup_bSat()
	{
	if (mDipoleModelParameterSet == KMW) {
	// KMW paper (arXiv:hep-ph/0606272), Table 3
	mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
	mQuarkMass[3] = 1.4;

	mBG = 4.;
	mMu02 = 1.17; // Gev^2
	mLambdaG = 0.02;
	mAg = 2.55;
	mC = 4;
	}
	else if (mDipoleModelParameterSet == HMPZ) {
	// Internal note by Heikki Mantysaari an Pia Zurita, arXiv pending
	mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.03; // u,d,s quarks
	mQuarkMass[3] = 1.3528;

	mBG = 4.;
	mMu02 = 1.1; // Gev^2
	mLambdaG = 0.08289;
	mAg = 2.1953;
	mC = 2.2894;
	}
	else if (mDipoleModelParameterSet == CUSTOM) {
	if (mCustomParameters.size() < 10) {
	cout << "DipoleModelParameters::setupParameters(): Error, require 8 custom parameters for bSAT when" << endl;
	cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	mQuarkMass[0] = mCustomParameters[0];
	mQuarkMass[1] = mCustomParameters[1];
	mQuarkMass[2] = mCustomParameters[2];
	mQuarkMass[3] = mCustomParameters[3];
	mQuarkMass[4] = mCustomParameters[4];

	mBG = mCustomParameters[5];
	mMu02 = mCustomParameters[6];
	mLambdaG = mCustomParameters[7];
	mAg = mCustomParameters[8];
	mC = mCustomParameters[9];
	if(mCustomParameters.size() > 10 and mCustomParameters.size() < 13){
	mShape = (ProtonShape) mCustomParameters[10];
	mV = mCustomParameters[11];
	mRshrink = mCustomParameters[12];
	}
	}
	else {
	cout << "DipoleModelParameters::setup_bSat(): Error, no known parameters for given dipole model" << endl;
	cout << " and requested parmeter set "
	<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	void DipoleModelParameters::setup_bNonSat()
	{
	if (mDipoleModelParameterSet == KMW) {
	// KT paper (arXiv:hep-ph/0304189v3), page 11
	mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
	mQuarkMass[3] = 1.4;

	mBG = 4.;
	mMu02 = 0.8;
	mLambdaG = -0.13;
	mAg = 3.5;
	mC = 4;
	}
	else if (mDipoleModelParameterSet == HMPZ) {
	// Internal note by Heikki Mantysaari an Pia Zurita, arXiv pending
	mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.1516; // u,d,s quarks
	mQuarkMass[3] = 1.3504;

	mBG = 4.;
	mMu02 = 1.1;
	mLambdaG = -0.006657;
	mAg = 3.0391;
	mC = 4.2974;

	}
	else if (mDipoleModelParameterSet == CUSTOM) {
	if (mCustomParameters.size() < 10) {
	cout << "DipoleModelParameters::setupParameters(): Error, require 8 custom parameters for bNonSAT when" << endl;
	cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	mQuarkMass[0] = mCustomParameters[0];
	mQuarkMass[1] = mCustomParameters[1];
	mQuarkMass[2] = mCustomParameters[2];
	mQuarkMass[3] = mCustomParameters[3];
	mQuarkMass[4] = mCustomParameters[4];

	mBG = mCustomParameters[5];
	mMu02 = mCustomParameters[6];
	mLambdaG = mCustomParameters[7];
	mAg = mCustomParameters[8];
	mC = mCustomParameters[9];
	if(mCustomParameters.size() > 10 and mCustomParameters.size() < 13){
	mShape = (ProtonShape)mCustomParameters[10];
	mV = mCustomParameters[11];
	}
	}
	else {
	cout << "DipoleModelParameters::setup_bNonSat(): Error, no known parameters for given dipole model" << endl;
	cout << " and requested parmeter set "
	<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	void DipoleModelParameters::setup_bCGC()
	{
	if (mDipoleModelParameterSet == KMW) {
	// WK paper (arXiv:0712.2670), Table II
	mQuarkMass[0] = mQuarkMass[1] = mQuarkMass[2] = 0.14; // u,d,s quarks
	mQuarkMass[3] = 1.4;

	mKappa = 9.9;
	mN0 = 0.558;
	mX0 = 1.84e-6;
	mLambda = 0.119;
	mGammas = 0.46;
	mBcgc = 7.5;
	}
	else if (mDipoleModelParameterSet == CUSTOM) {
	if (mCustomParameters.size() < 10) {
	cout << "DipoleModelParameters::setup_bCGC(): Error, require 10 custom parameters for bCGC when" << endl;
	cout << " dipole-model parameter set is 'CUSTOM'. Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	mQuarkMass[0] = mCustomParameters[0];
	mQuarkMass[1] = mCustomParameters[1];
	mQuarkMass[2] = mCustomParameters[2];
	mQuarkMass[3] = mCustomParameters[3];

	mKappa = mCustomParameters[4];
	mN0 = mCustomParameters[5];
	mX0 = mCustomParameters[6];
	mLambda = mCustomParameters[7];
	mGammas = mCustomParameters[8];
	mBcgc = mCustomParameters[9];
	}
	else {
	cout << "DipoleModelParameters::setup_bCGC(): Error, no known parameters for given dipole model"
	<< " and requested parmeter set "
	<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	void DipoleModelParameters::setupParameters()
	{
	TableGeneratorSettings *settings = TableGeneratorSettings::instance();
	if (mDipoleModelParameterSet == CUSTOM) mCustomParameters = settings->dipoleModelCustomParameters();

	//
	// Init
	//
	mKappa = 0;
	mN0 = 0;
	mX0 = 0;
	mLambda = 0;
	mGammas = 0;
	mBcgc = 0;

	mBG = 0;
	mMu02 = 0;
	mLambdaG = 0;
	mAg = 0;
	mC = 0;

	mV = 0;

	mBoostedGaussianR2_rho = 0;
	mBoostedGaussianNL_rho = 0;
	mBoostedGaussianNT_rho = 0;
	mBoostedGaussianQuarkMass_rho = 0;
	mBoostedGaussianR2_phi = 0;
	mBoostedGaussianNL_phi = 0;
	mBoostedGaussianNT_phi = 0;
	mBoostedGaussianQuarkMass_phi = 0;
	mBoostedGaussianR2_jpsi = 0;
	mBoostedGaussianNL_jpsi = 0;
	mBoostedGaussianNT_jpsi = 0;
	mBoostedGaussianQuarkMass_jpsi = 0;

	//
	// t masses (not used, just for completeness)
	//
	mQuarkMass[5] = 175.; // t quark consistent with HMPZ

	//
	// Parameters for boosted Gaussian wave function
	//
	setup_boostedGaussiansWaveFunction();

	//
	// Model parameters
	//
	if (mDipoleModelType == bSat) {
	setup_bSat();
	}
	else if (mDipoleModelType == bNonSat) {
	setup_bNonSat();
	}
	else if (mDipoleModelType == bCGC) {
	setup_bCGC();
	}
	else {
	cout << "DipoleModelParameters::setupParameters(): Error, no known parameters for given dipole model" << endl;
	cout << " and requested parmeter set "
	<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
	- exit(1);
	- }
	+// exit(1);
	+ }
	}

	double DipoleModelParameters::boostedGaussianR2(int vm)
	{
	if (vm == 113)
	return mBoostedGaussianR2_rho;
	else if (vm == 333)
	return mBoostedGaussianR2_phi;
	else if (vm == 443)
	return mBoostedGaussianR2_jpsi;
	else {
	cout << "DipoleModelParameters::boostedGaussianR2(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	double DipoleModelParameters::boostedGaussianNL(int vm)
	{
	if (vm == 113)
	return mBoostedGaussianNL_rho;
	else if (vm == 333)
	return mBoostedGaussianNL_phi;
	else if (vm == 443)
	return mBoostedGaussianNL_jpsi;
	else {
	cout << "DipoleModelParameters::boostedGaussianNL(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	double DipoleModelParameters::boostedGaussianNT(int vm)
	{
	if (vm == 113)
	return mBoostedGaussianNT_rho;
	else if (vm == 333)
	return mBoostedGaussianNT_phi;
	else if (vm == 443)
	return mBoostedGaussianNT_jpsi;
	else {
	cout << "DipoleModelParameters::boostedGaussianNT(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	void DipoleModelParameters::setup_boostedGaussiansWaveFunction()
	{
	if (mDipoleModelParameterSet == KMW) {
	- //
	+
	// KMW: bSat, bNonSat, and bCGC use the same parameters
	// and also does not distingiosh between T and L.
	//
	mBoostedGaussianR2_rho = 12.9;
	mBoostedGaussianNL_rho = 0.853;
	mBoostedGaussianNT_rho = 0.911;
	mBoostedGaussianQuarkMass_rho = 0.14;

	mBoostedGaussianR2_phi = 11.2;
	mBoostedGaussianNL_phi = 0.825;
	mBoostedGaussianNT_phi= 0.919;
	mBoostedGaussianQuarkMass_phi = 0.14;

	mBoostedGaussianR2_jpsi = 2.3;
	mBoostedGaussianNL_jpsi = 0.575;
	mBoostedGaussianNT_jpsi = 0.578;
	mBoostedGaussianQuarkMass_jpsi = 1.4;
	}
	else if (mDipoleModelParameterSet == HMPZ) {
	if (mDipoleModelType == bSat) {
	mBoostedGaussianR2_rho = 3.6376*3.6376;
	mBoostedGaussianNL_rho = 0.8926;
	mBoostedGaussianNT_rho = 0.9942;
	mBoostedGaussianQuarkMass_rho = 0.03;

	mBoostedGaussianR2_phi = 3.3922*3.3922;
	mBoostedGaussianNL_phi = 0.8400;
	mBoostedGaussianNT_phi = 0.9950;
	mBoostedGaussianQuarkMass_phi = 0.03;

	mBoostedGaussianR2_jpsi = 1.5070*1.5070;
	mBoostedGaussianNL_jpsi = 0.5860;
	mBoostedGaussianNT_jpsi = 0.5890;
	mBoostedGaussianQuarkMass_jpsi = 1.3528;
	}
	else if (mDipoleModelType == bNonSat) {
	mBoostedGaussianR2_rho = 3.5750*3.5750;
	mBoostedGaussianNL_rho = 0.8467;
	mBoostedGaussianNT_rho = 0.8978;
	mBoostedGaussianQuarkMass_rho = 0.1516;

	mBoostedGaussianR2_phi = 3.3530*3.3530;
	mBoostedGaussianNL_phi = 0.8196;
	mBoostedGaussianNT_phi = 0.9072;
	mBoostedGaussianQuarkMass_phi = 0.1516;

	mBoostedGaussianR2_jpsi = 1.5071*1.5071;
	mBoostedGaussianNL_jpsi = 0.5868;
	mBoostedGaussianNT_jpsi = 0.5899;
	mBoostedGaussianQuarkMass_jpsi = 1.3504;
	}
	else if (mDipoleModelType == bCGC) {
	cout << "DipoleModelParameters::setup_boostedGaussiansWaveFunction(): "
	"Error, no HMPZ wave function parameters for CGC model. Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	}
	else {
	cout << "DipoleModelParameters::setup_boostedGaussiansWaveFunction(): Error, no known parameters for given dipole model "
	<< " parmeter set "
	<< "(" << mDipoleModelType << "/" << mDipoleModelParameterSet << "). Stop." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	double DipoleModelParameters::boostedGaussianQuarkMass(int vm)
	{
	if (vm == 113)
	return mBoostedGaussianQuarkMass_rho;
	else if (vm == 333)
	return mBoostedGaussianQuarkMass_phi;
	else if (vm == 443)
	return mBoostedGaussianQuarkMass_jpsi;
	else {
	cout << "DipoleModelParameters::boostedGaussianQuarkMass(): Error, no boosted Gaussian parameter parameters for given vector meson (" << vm << ")." << endl;
	- exit(1);
	+// exit(1);
	}
	}

	bool DipoleModelParameters::list(ostream& os)
	{
	const int fieldWidth = 32;

	os << "\nDipole Model Parameters:" << endl;

	os << setw(fieldWidth) << "Set: " << mDipoleModelParameterSetName << endl;
	os << setw(fieldWidth) << "Quark masses: "
	<< "u=" << mQuarkMass[0]
	<< ", d=" << mQuarkMass[1]
	<< ", s=" << mQuarkMass[2]
	<< ", c=" << mQuarkMass[3]
	<< ", b=" << mQuarkMass[4]
	<< ", t=" << mQuarkMass[5] << endl;

	os << setw(fieldWidth) << "BG: " << mBG << endl;
	os << setw(fieldWidth) << "Mu02: " << mMu02 << endl;
	os << setw(fieldWidth) << "LambdaG: " << mLambdaG << endl;
	os << setw(fieldWidth) << "Ag: " << mAg << endl;
	os << setw(fieldWidth) << "C: " << mC << endl;

	os << setw(fieldWidth) << "Kappa: " << mKappa << endl;
	os << setw(fieldWidth) << "N0: " << mN0 << endl;
	os << setw(fieldWidth) << "X0: " << mX0 << endl;
	os << setw(fieldWidth) << "Lambda: " << mLambda << endl;
	os << setw(fieldWidth) << "Gammas: " << mGammas << endl;
	os << setw(fieldWidth) << "Bcgc: " << mBcgc << endl;

	os << setw(fieldWidth) << "BoostedGaussianR2_rho: " << mBoostedGaussianR2_rho << endl;
	os << setw(fieldWidth) << "BoostedGaussianNL_rho: " << mBoostedGaussianNL_rho << endl;
	os << setw(fieldWidth) << "BoostedGaussianNT_rho: " << mBoostedGaussianNT_rho << endl;
	os << setw(fieldWidth) << "BoostedGaussianQuarkMass_rho: " << mBoostedGaussianQuarkMass_rho << endl;
	os << setw(fieldWidth) << "BoostedGaussianR2_phi: " << mBoostedGaussianR2_phi << endl;
	os << setw(fieldWidth) << "BoostedGaussianNL_phi: " << mBoostedGaussianNL_phi << endl;
	os << setw(fieldWidth) << "BoostedGaussianNT_phi: " << mBoostedGaussianNT_phi << endl;
	os << setw(fieldWidth) << "BoostedGaussianQuarkMass_phi: " << mBoostedGaussianQuarkMass_phi << endl;
	os << setw(fieldWidth) << "BoostedGaussianR2_jpsi: " << mBoostedGaussianR2_jpsi << endl;
	os << setw(fieldWidth) << "BoostedGaussianNL_jpsi: " << mBoostedGaussianNL_jpsi << endl;
	os << setw(fieldWidth) << "BoostedGaussianNT_jpsi: " << mBoostedGaussianNT_jpsi << endl;
	os << setw(fieldWidth) << "BoostedGaussianQuarkMass_jpsi: " << mBoostedGaussianQuarkMass_jpsi << endl;

	os << endl;

	return true;
	}

	void DipoleModelParameters::setBG(double val) {
	mBG = val;
	mCustomParameters[5] = val;
	}

	void DipoleModelParameters::setMu02(double val) {
	mMu02 = val;
	mCustomParameters[6] = val;
	}

	void DipoleModelParameters::setLambdaG(double val) {
	mLambdaG = val;
	mCustomParameters[7] = val;
	}

	void DipoleModelParameters::setAg(double val) {
	mAg = val;
	mCustomParameters[8] = val;
	}

	void DipoleModelParameters::setC(double val) {
	mC = val;
	mCustomParameters[9] = val;
	}

	void DipoleModelParameters::setRshrink(double val) {
	mRshrink = val;
	mCustomParameters[12] = val;
	}

	void DipoleModelParameters::setV(double val) {
	mV = val;
	mCustomParameters[11] = val;
	}


	void DipoleModelParameters::setQuarkMass(int index, double val) {
	mQuarkMass[index] = val;
	mCustomParameters[index] = val;
	}





	Index: branches/fitting/src/DipoleModel.cpp
	===================================================================
	--- branches/fitting/src/DipoleModel.cpp (revision 385)
	+++ branches/fitting/src/DipoleModel.cpp (revision 386)
	@@ -1,751 +1,749 @@
	//==============================================================================
	// DipoleModel.cpp
	//
	// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
	//
	// This file is part of Sartre version: 1.1
	//
	// This program is free software: you can redistribute it and/or modify
	// it under the terms of the GNU General Public License as published by
	// the Free Software Foundation.
	// This program is distributed in the hope that it will be useful,
	// but without any warranty; without even the implied warranty of
	// merchantability or fitness for a particular purpose. See the
	// GNU General Public License for more details.
	// You should have received a copy of the GNU General Public License
	// along with this program. If not, see <http://www.gnu.org/licenses/>.
	//
	// Author: Tobias Toll
	// Last update:
	// $Date$
	// $Author$
	//==============================================================================
	#include <fstream>
	#include <iostream>
	#include <sstream>
	#include <cmath>
	#include "DipoleModel.h"
	#include "TableGeneratorSettings.h"
	#include "DglapEvolution.h"
	#include "Parameters_Minuit2.h"
	#include "Constants.h"
	#include "TFile.h"
	#include "TVector3.h"
	#include "TMath.h"
	#include "TH2F.h"
	#include "TF1.h"
	#include "Math/WrappedTF1.h"
	#include "Math/GaussIntegrator.h"
	#include "TObject.h"
	#include "Minuit2/MnUserParameters.h"
	//#include "LHAPDF/LHAPDF.h"


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

	using namespace std;

	DipoleModel::DipoleModel()
	{
	mConfigurationExists=false;
	TableGeneratorSettings* settings = TableGeneratorSettings::instance();
	unsigned int A = settings->A();
	mNucleus.init(A);
	mIsInitialized=true;
	//comparisons with previous version in ep yields
	//best agreement with as(MZ)=0.1196 at NLO:
	// mAs.init(0.1196, 2);
	// mAs.init(0.13939, 1); //Default Pythia 8 Too large
	// mAs.init(params_coupling, params_asOrder);
	// mAs.init(0.126, 1);//

	//
	// Initialise DGLAP evolution
	//
	cout << "DipoleModel(): Initializing the DGLAP engine:" << endl;
	//Upper scale in the DGLAP evoultion. Should be high enough to cover LHeC:
	// DglapEvolution::instance().setS(1e10); //this should do it...
	// DglapEvolution::instance().G(0.01, 4.);
	mUseSkewednessCorrection=false;
	}

	DipoleModel::~DipoleModel() { /* no op */}

	const TableGeneratorNucleus* DipoleModel::nucleus() const { return &mNucleus; }

	bool DipoleModel::configurationExists() const { return mConfigurationExists; }

	double DipoleModel::bDependence(double) { return 0; }

	double DipoleModel::bDependence(double, double) { return 0; }

	double DipoleModel::dsigmadb2ep(double, double, double) { return 0;}

	double DipoleModel::ddxdsigmadb2(double, double, double) { return 0;}

	//void DipoleModel::setParams(const std::vector<double>& par) {/this needs to return nothing/}

	//*********bSat:***************************************************

	void DipoleModel::setParams(const std::vector<double>& par)
	{

	mBg = par[0]; //par[0];//BS For MZ, BG is kept fixed
	mLmass = par[1];
	mCmass = par[2];
	mBmass = par[3];
	mC = par[4];
	mMu_02 = par[5]; //par[5];//BS For MZ, Mu02 is kept fixed
	mAg = par[6];
	mLamg = par[7];
	mRmax = par[8];
	mAs_qqg = par[9];
	mRhs = par[10];
	mnLP = par[11];
	mR0WS = par[12];
	mdWS = par[13];
	muCE = par[14];
	mEP = par[15];
	BV = par[16];
	RV = par[17];
	wV = par[18];

	qMass[0] = mLmass;
	qMass[1] = mLmass;
	qMass[2] = mLmass;
	qMass[3] = mCmass;
	qMass[4] = mBmass;
	qMass[5] = params_t_mass;
	//cout << "DipoleModel mAg = " + std::to_string(mAg) << endl;


	}

	DipoleModel_bSat::DipoleModel_bSat()
	{
	// mPDF=LHAPDF::mkPDF("HERAPDF20_LO_EIG", 0); //LHAPDF
	mBDependence = 0;
	mSigma_ep_LookupTable = 0;
	mThicknessNormTable = 0;
	}

	DipoleModel_bSat::~DipoleModel_bSat()
	{
	// delete mPDF;
	delete mSigma_ep_LookupTable;
	delete mThicknessNormTable;
	delete mBDependence;
	}

	DipoleModel_bSat& DipoleModel_bSat::operator=(const DipoleModel_bSat& dp)
	{
	if (this != &dp) {
	delete mBDependence;
	delete mSigma_ep_LookupTable;
	delete mThicknessNormTable;

	DipoleModel::operator=(dp);
	mBDependence = new TH2F(*(dp.mBDependence));
	mBDependence->SetDirectory(0);
	mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable));
	mThicknessNormTable = new TH1F(*(dp.mThicknessNormTable));
	mSigma_ep_LookupTable->SetDirectory(0);
	mThicknessNormTable->SetDirectory(0);

	}
	return *this;
	}

	DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp)
	{
	mBDependence = new TH2F(*(dp.mBDependence));
	mBDependence->SetDirectory(0);
	}

	void DipoleModel::createLaplaceThickness(){
	double zbhigh=30.mnLPhbarc; //fm
	int numBins=100000;
	PR(zbhigh);
	mLP_lu=new TH1D("LP_lu", "LP_lu", numBins, 0, zbhigh);
	double binWidth=zbhigh/numBins; //fm
	double norm=0;
	for(int i=1; i<=numBins; i++){
	double b=mLP_lu->GetBinCenter(i); //fm
	TF1 LP("LP", this, &DipoleModel_bSat::uiLaplace, -zbhigh, zbhigh, 3);
	LP.SetParameter(0, b);
	LP.SetParameter(1, mnLP);
	//WS.SetParameter(2, params_dWS);
	ROOT::Math::WrappedTF1 wf(LP);
	ROOT::Math::GaussIntegrator ig;
	ig.SetFunction(wf);
	ig.SetRelTolerance(1e-3);
	double LP_of_b=ig.Integral(-zbhigh, zbhigh)/hbarc; //GeV^2
	mLP_lu->SetBinContent(i, LP_of_b);
	//
	// Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db)
	// ~2pisum(T'(bi)bi*Delta b)
	//
	norm+=2.M_PIb/hbarcLP_of_bbinWidth/hbarc; //GeV^0
	}
	PR(norm);
	mLP_lu->Scale(1./norm); //normalise, GeV^2
	//
	// Calculate <b^2>=integral(b^2 T(b)d^2b)=integral(2pi b^3T(b)db)
	//
	// PR(mLP_lu->Interpolate(0.5));
	// PR(mLP_lu->Interpolate(1));
	// PR(mLP_lu->Interpolate(2));

	// PR((1/(4M_PImnLPmnLP))0.5/(mnLPhbarc)TMath::BesselK1(0.5/(mnLP*hbarc)));
	// PR((1/(4M_PImnLPmnLP))1/(mnLPhbarc)TMath::BesselK1(1/(mnLP*hbarc)));
	// PR((1/(4M_PImnLPmnLP))2/(mnLPhbarc)TMath::BesselK1(2/(mnLP*hbarc)));

	double b2_mean=0.;
	for(int i=1; i<=numBins; i++){
	double b=mLP_lu->GetBinCenter(i)/hbarc;
	double T_of_b=mLP_lu->GetBinContent(i);
	b2_mean+=bbb*T_of_b; //GeV-1
	}
	b2_mean=2.M_PI*binWidth/hbarc; //GeV^-2
	PR(sqrt(b2_mean*hbarc2));
	}//BS


	double DipoleModel::uiLaplace(double* x, double* par){
	//All in fm
	double z=*x; //fm
	double b=par[0]; //fm
	double r=par[1]; //Radius GeV^-1
	//double d=par[2]; //thickness GeV^-1

	double arg=(sqrt((zz+bb))/(hbarc*r));
	return (1/(8M_PIrrr))*exp(-arg); //GeV^3
	}

	void DipoleModel::createCompressedExponentThickness(){
	double zbhigh=30;//30.(1/muCE)pow(hbarc, 1); //fm
	double rhigh=sqrt(1/1.5)/hbarc;
	int numBinsb=10000;
	int numBinsr=1000;
	// int numBinsb=100;
	// int numBinsr=10;
	PR(zbhigh);
	PR(rhigh);
	mCE_lu=new TH2D("CE_lu", "CE_lu", numBinsb, 0, zbhigh, numBinsr, 0, rhigh); //r Upper limit can be changed
	//mCE_lu=new TH1D("CE_lu", "CE_lu", numBinsb, 0, zbhigh);//, numBinsr, 0, rhigh); //r Upper limit can be changed
	double binWidthb=zbhigh/numBinsb; //fm
	double binWidthr=rhigh/numBinsr;
	double norm=0;
	PR(mEP);
	PR(muCE);
	for(int i=1; i<=numBinsb; i++){
	for(int j=1; j<=numBinsr; j++){
	double b=mCE_lu->GetXaxis()->GetBinCenter(i); //fm
	double r=mCE_lu->GetYaxis()->GetBinCenter(j); //fm
	TF1 CE("CE", this, &DipoleModel::uiCompressedExponent, -zbhigh, zbhigh, 4);
	CE.SetParameter(0, b);
	CE.SetParameter(1, r);
	CE.SetParameter(2, muCE);
	CE.SetParameter(3, mEP);
	//if(r>0.8) PR(r);
	ROOT::Math::WrappedTF1 wf(CE);
	ROOT::Math::GaussIntegrator ig;
	ig.SetFunction(wf);
	ig.SetRelTolerance(1e-6);
	double CE_of_b=2.*ig.IntegralUp(0.)/hbarc;//(-zbhigh, zbhigh)/hbarc; //GeV^2
	//PR(CE_of_b);
	// PR(b);
	// PR(r);
	// PR(muCE);
	// PR(mEP);
	// PR(zbhigh);
	// PR(rhigh);
	// PR(CE_of_b);
	mCE_lu->SetBinContent(i, j, CE_of_b);
	//
	// Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db)
	// ~2pisum(T'(bi)bi*Delta b)
	//
	norm+=2.M_PIb/hbarcCE_of_bbinWidthb/hbarc; //GeV^0
	}
	}
	// PR(norm);
	// mCE_lu->Scale(1./norm); //normalise, GeV^2
	//
	// Calculate <b^2>=integral(b^2 T(b)d^2b)=integral(2pi b^3T(b)db)
	// //
	// double* z=0;
	// double parameters[3]={1, muCE, mEP};

	// PR(DipoleModel::uiCompressedExponent(z, parameters));
	// PR(mCE_lu->Interpolate(0.5, 0));//, 1.));
	// PR(mCE_lu->Interpolate(1, 0));//, 2.));
	// PR(mCE_lu->Interpolate(2, 0));//, 4.));

	// PR(((pow(muCE,3/mEP)mEP)/(4M_PI*TMath::Gamma(3/mEP))));

	// PR((1/pow(2M_PImBg, 1.5)));

	// // double u=muCE;
	// // double v=mEP;
	// // double z=0;
	// // double b=0.5;
	// // double arg=u(pow((sqrt(zz+b*b)/hbarc),v));
	// // PR(((pow(u,3/v)v)/(4M_PITMath::Gamma(3/v)))exp(-arg));
	// // b=1;
	// // arg=u(pow((sqrt(zz+b*b)/hbarc),v));
	// // PR(((pow(u,3/v)v)/(4M_PITMath::Gamma(3/v)))exp(-arg));
	// // b=2;
	// // arg=u(pow((sqrt(zz+b*b)/hbarc),v));
	// // PR(((pow(u,3/v)v)/(4M_PITMath::Gamma(3/v)))exp(-arg));

	// PR((1/pow(2M_PImBg, 1.5)) * exp(-(0.50.5)/(hbarc22*mBg)));
	// PR((1/pow(2M_PImBg, 1.5)) * exp(-(11)/(hbarc22*mBg)));
	// PR((1/pow(2M_PImBg, 1.5)) * exp(-(22)/(hbarc22*mBg)));

	// double zvalue=0;
	// double values1[4]={0.5, 0., muCE, mEP};
	// double values2[4]={1., 0., muCE, mEP};
	// double values3[4]={2., 0., muCE, mEP};
	// PR(uiCompressedExponent(&zvalue, values1));
	// PR(uiCompressedExponent(&zvalue, values2));
	// PR(uiCompressedExponent(&zvalue, values3));
	// // PR((1/(4M_PImnLPmnLP))0.5/(mnLPhbarc)TMath::BesselK1(0.5/(mnLP*hbarc)));
	// // PR((1/(4M_PImnLPmnLP))1/(mnLPhbarc)TMath::BesselK1(1/(mnLP*hbarc)));
	// // PR((1/(4M_PImnLPmnLP))2/(mnLPhbarc)TMath::BesselK1(2/(mnLP*hbarc)));

	// // PR(1/(2M_PImBg) * exp(-00/((hbarc2)(2*mBg))));
	// PR((1/pow(2M_PImBg, 1)) * exp(-(0.50.5)/(hbarc22*mBg)));
	// PR((1/pow(2M_PImBg, 1)) * exp(-(11)/(hbarc22*mBg)));
	// PR((1/pow(2M_PImBg, 1)) * exp(-(22)/(hbarc22*mBg)));
	// // /double b2_mean=0.;
	// // for(int i=1; i<=numBins; i++){
	// // double b=mCE_lu->GetBinCenter(i)/hbarc;
	// // double T_of_b=mCE_lu->GetBinContent(i);
	// // b2_mean+=bbb*T_of_b; //GeV-1
	// // }
	// // b2_mean=2.M_PI*binWidth/hbarc; //GeV^-2
	// // PR(sqrt(b2_mean*hbarc2));
	}//BS


	double DipoleModel::uiCompressedExponent(double* x, double* par){
	//All in fm
	double z=*x; //fm
	double b=par[0];
	double r=par[1]; //fm
	double u=par[2]; //Radius GeV^-1
	double exponentparameter=par[3];
	//PR(exponentparameter);
	//double v=2-exp(-exponentparameter*1e-2/(pow(r,6.)));
	//double v=1+exp(-exponentparameterrr);
	double v=2-exp(-exponentparameter*1e-4/(pow(r,4.)));
	//double v=1+exp(-exponentparameter*pow(r,4.));
	//double v=max(2.95673-exponentparameter*r, 1.);
	//double v=1+exp(-exponentparameter*pow(r,9));
	//double v=par[2];
	//double v=exponentparameter;
	// double arg=u(pow((sqrt(zz+b*b)/hbarc),v));
	double arg=u(pow((zz+bb)/hbarc2, 0.5v));
	return pow(u,3./v)v / (4.M_PITMath::Gamma(3./v)) exp(-arg); //GeV^
	}

	void DipoleModel::createWoodsSaxonThickness(){
	double zbhigh=3.params_R0WShbarc; //fm
	int numBins=100000;

	mWS_lu=new TH1D("WS_lu", "WS_lu", numBins, 0, zbhigh);
	double binWidth=zbhigh/numBins; //fm
	double norm=0;
	for(int i=1; i<=numBins; i++){
	double b=mWS_lu->GetBinCenter(i); //fm
	TF1 WS("WS", this, &DipoleModel_bSat::uiWoodsSaxon, -zbhigh, zbhigh, 3);
	WS.SetParameter(0, b);
	WS.SetParameter(1, params_R0WS);
	WS.SetParameter(2, params_dWS);
	ROOT::Math::WrappedTF1 wf(WS);
	ROOT::Math::GaussIntegrator ig;
	ig.SetFunction(wf);
	ig.SetRelTolerance(1e-2);
	double WS_of_b=ig.Integral(-zbhigh, zbhigh)/hbarc; //GeV^-1
	mWS_lu->SetBinContent(i, WS_of_b);
	//
	// Calculate the norm=integral(T'(b)d^2b=integral(2pibT'(b)db)
	// ~2pisum(T'(bi)bi*Delta b)
	//
	norm+=2.M_PIb/hbarcWS_of_bbinWidth/hbarc; //GeV^-3
	}
	PR(norm);
	mWS_lu->Scale(1./norm); //normalise, GeV^2
	//
	// Calculate <b^2>=integral(b^2 T(b)d^2b)=integral(2pi b^3T(b)db)
	//
	double b2_mean=0.;
	for(int i=1; i<=numBins; i++){
	double b=mWS_lu->GetBinCenter(i)/hbarc;
	double T_of_b=mWS_lu->GetBinContent(i);
	b2_mean+=bbb*T_of_b; //GeV-1
	}
	b2_mean=2.M_PI*binWidth/hbarc; //GeV^-2
	PR(sqrt(b2_mean));
	}


	double DipoleModel::uiWoodsSaxon(double* x, double* par){
	//All in fm
	double z=*x; //fm
	double b=par[0]; //fm
	double R0=par[1]; //Radius GeV^-1
	double d=par[2]; //thickness GeV^-1

	double arg=(sqrt((zz+bb))/hbarc-R0)/d;
	return 1./(1+exp(arg)); //GeV^0
	}

	void DipoleModel_bSat::createConfiguration(int iConfiguration)
	{
	if (!mIsInitialized) {
	cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl;
	exit(1);
	}
	TableGeneratorSettings* settings = TableGeneratorSettings::instance();
	unsigned int A = mNucleus.A();
	string path=settings->bSatLookupPath();
	ostringstream filename;
	filename.str("");
	filename << path << "/bSat_bDependence_A" << A <<".root";
	ifstream ifs(filename.str().c_str());
	if (!ifs) {
	cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl;
	cout << "Stopping." << endl;
	exit(1);
	}
	TFile* lufile= new TFile(filename.str().c_str());
	ostringstream histoName;
	histoName.str( "" );
	histoName << "Configuration_" << iConfiguration;
	lufile->GetObject( histoName.str().c_str(), mBDependence );
	mBDependence->SetDirectory(0);
	lufile->Close();
	mConfigurationExists=true;
	}


	double DipoleModel_bSat::dsigmadb2(double r, double b, double phi, double xprobe)
	{
	double bDep=bDependence(b, phi);
	double muQ2 = mC/(r*r/hbarc2) + mMu_02;
	double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;
	double result = 2.*(1. - exp(-omega/2));
	return result;
	}

	double DipoleModel_bSat::ddxdsigmadb2(double r, double b, double xprobe){
	/*#TT
	const double BG = mBg; // GeV^-2
	double arg = bb/(2BG);
	arg /= hbarc2;
	double bDep= 1/(2M_PIBG) * exp(-arg);
	double muQ2 = mC/(r*r/hbarc2) + mMu_02;
	double alphas=mAs.alphaS(muQ2);
	double Gofx=DglapEvolution::instance().G(xprobe, muQ2);
	// double Gofx=mPDF->xfxQ2(21, xprobe, muQ2)/xprobe;
	double omegaFactor = ((M_PIM_PI)/Nc)(rr/hbarc2)alphas*bDep;
	double h=xprobe*1e-3;
	double Gofxp=DglapEvolution::instance().G(xprobe+h, muQ2);
	double Gofxm=DglapEvolution::instance().G(xprobe-h, muQ2);
	// double Gofxp=mPDF->xfxQ2(21, xprobe+h, muQ2)/(xprobe+h);
	// double Gofxm=mPDF->xfxQ2(21, xprobe-h, muQ2)/(xprobe-h);
	double ddx_xg=((xprobe+h)Gofxp-(xprobe-h)Gofxm)/(2*h);
	double result = omegaFactorddx_xgexp(-omegaFactorxprobeGofx/2);
	return result;
	*/
	return 0.;
	}

	double DipoleModel_bSat::bDependence(double b, double phi)
	{
	double result = mBDependence->Interpolate(b, phi);
	return result;
	}


	double DipoleModel::protonThickness(double b, double r=0){
	double result=0;
	switch(params_Shape){
	case HardSphere:
	//
	// T(b)=3/2pi1/R^3sqrt(R^2-b^2)
	//
	// <b>=0, <b^2>=2/5*R^2
	// sigma(b)=0.63*R
	//
	if(b/hbarc<params_Rhs)
	result=3./2.sqrt(params_Rhsparams_Rhs-bb/hbarc2)/(M_PIpow(params_Rhs, 3));
	break;
	case Laplace:
	result = (1/(4M_PImnLPmnLP))b/(mnLPhbarc)TMath::BesselK1(b/(mnLP*hbarc)); //GeV^2
	break;
	case CompressedExponent:
	result= mCE_lu->Interpolate(b, r); //GeV^2
	break;
	case Gauss:
	default:
	//
	// T(b)=1/(2piBG)exp(-b2/(2BG))
	//
	// <b>=0, <b^2>=BG
	// sigma(b)=sqrt(BG)
	//
	const double BG = mBg; // GeV^-2
	double arg = bb/(2BG); //fm2*GeV^2
	arg /= hbarc2; //GeV^0
	result= 1/(2M_PIBG) * exp(-arg); //GeV^2
	break;
	}
	return result;
	}

	void DipoleModel_bSat::setMu02(double val){ mMu02=val; }

	void DipoleModel_bSat::useSkewednessCorrection(){ mUseSkewednessCorrection=true; }

	double DipoleModel_bSat::calculateNorm(double w){
	double norm=0;
	double step=1e-2;
	mBg=4;
	mnLP=1.8282;

	// double T_old=0;
	for(double b=step; b<3.5; b+=step){
	double T_gauss=1./(2M_PImBg)exp(-bb/hbarc2/2./mBg);
	double T_laplace=(1/(4M_PImnLPmnLP))b/(mnLPhbarc)TMath::BesselK1(b/(mnLP*hbarc));
	double T=(1-w)T_gauss+wT_laplace;
	// T=0.5*(T+T_old);
	// T_old=T;
	norm+=2M_PIbT/hbarc2step;
	}
	return norm;
	}

	void DipoleModel_bSat::createNormTable(){
	mBg=4; //#TT
	mnLP=1.8282; //#TT
	mThicknessNormTable = new TH1F("", "", 1000, -5., 5.);
	for(int i=1; i<=1000; i++){
	//Omega = X(r, b_ * T(b), this is a function of X(r, xprobe)
	double X=pow(10, mThicknessNormTable->GetBinCenter(i));
	double norm=0;
	double db=1e-3;
	for(double b=1e-4; b<4.; b+=db){
	//find weight:
	double T_gauss=1./(2M_PImBg)exp(-bb/hbarc2/2./mBg);
	double T_laplace=(1/(4M_PImnLPmnLP))b/(mnLPhbarc)TMath::BesselK1(b/(mnLP*hbarc));
	double weight=0.5;
	double dsigdb2=0;
	int count=0;
	int countMax=1e4;
	double bDep=0;
	while(true){
	double weight_old=weight;
	bDep= (1-weight)T_gauss+weightT_laplace;
	double omega = X*bDep;
	dsigdb2 = 2.*(1. - exp(-omega/2));
	weight=0.5*dsigdb2;
	if(abs(weight_old-weight)/weight < 1e-5 or weight==0 or count>countMax) break;
	}
	norm+=2M_PIbbDepdb/hbarc2;
	}
	mThicknessNormTable->SetBinContent(i, norm);
	}
	}

	//
	// Version which interpolates between different shapes:
	//
	+double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
	+{
	+ //
	+ //modify dipole size as http://arxiv.org/abs/0807.0325v1
	+ //
	+ r=mRmaxsqrt(log(1+rr/(mRmax*mRmax)));
	+ double T_gauss=1./(2M_PImBg)exp(-bb/hbarc2/2./mBg);
	+ double T_laplace=(1/(4M_PImnLPmnLP))b/(mnLPhbarc)TMath::BesselK1(b/(mnLP*hbarc));
	+ double muQ2 = mC/(r*r/hbarc2) + mMu_02; //GeV2
	+ double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	+
	+ //Use an iterative approach to find the thickness
	+ double weight=0.5;
	+ double dsigdb2=0;
	+ int count=0;
	+ int countMax=1e4;
	+ double bDep=0;
	+ while(true){
	+ double weight_old=weight;
	+ bDep= (1-weight)T_gauss+weightT_laplace;
	+ double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;
	+ dsigdb2 = 2.*(1. - exp(-omega/2));
	+ weight=0.5*dsigdb2;
	+ if(abs(weight_old-weight)/weight < 1e-5 or weight==0 or count>countMax) break;
	+ }
	+ //We need to normalise the thickness
	+ // double norm=calculateNorm(weight);
	+ double X=((M_PIM_PI)/Nc)(rr/hbarc2)asxg;
	+ bDep/=mThicknessNormTable->Interpolate(log10(X));
	+ dsigdb2 = 2.(1. - exp(-XbDep/2));
	+
	+// return bDep;
	+ return dsigdb2;
	+}
	+
	+//original version:
	//double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
	//{
	// //
	// //modify dipole size as http://arxiv.org/abs/0807.0325v1
	// //
	// //r=mRmaxsqrt(log(1+rr/(mRmax*mRmax)));
	-// mBg=4;
	-// mnLP=1.8282;
	-// double T_gauss=1./(2M_PImBg)exp(-bb/hbarc2/2./mBg);
	-// double T_laplace=(1/(4M_PImnLPmnLP))b/(mnLPhbarc)TMath::BesselK1(b/(mnLP*hbarc));
	+// double bDep=1./(2M_PImBg)exp(-bb/hbarc2/2./mBg);
	// double muQ2 = mC/(r*r/hbarc2) + mMu_02; //GeV2
	// double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	-//
	-// //Use an iterative approach to find the thickness
	-// double weight=0.5;
	-// double dsigdb2=0;
	-// int count=0;
	-// int countMax=1e4;
	-// double bDep=0;
	-// while(true){
	-// double weight_old=weight;
	-// bDep= (1-weight)T_gauss+weightT_laplace;
	-// double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;
	-// dsigdb2 = 2.*(1. - exp(-omega/2));
	-// weight=0.5*dsigdb2;
	-// if(abs(weight_old-weight)/weight < 1e-5 or weight==0 or count>countMax) break;
	-// }
	-// //We need to normalise the thickness
	-// // double norm=calculateNorm(weight);
	-// double X=((M_PIM_PI)/Nc)(rr/hbarc2)asxg;
	-// // bDep/=mThicknessNormTable->Interpolate(log10(X));
	-// dsigdb2 = 2.(1. - exp(-XbDep/2));
	-//
	-// return bDep;
	+// double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;
	+// double dsigdb2 = 2.*(1. - exp(-omega/2));
	// return dsigdb2;
	//}

	-//original version:
	-double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
	-{
	- //
	- //modify dipole size as http://arxiv.org/abs/0807.0325v1
	- //
	- //r=mRmaxsqrt(log(1+rr/(mRmax*mRmax)));
	- double bDep=1./(2M_PImBg)exp(-bb/hbarc2/2./mBg);
	- double muQ2 = mC/(r*r/hbarc2) + mMu_02; //GeV2
	- double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	- double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;
	- double dsigdb2 = 2.*(1. - exp(-omega/2));
	- return dsigdb2;
	-}
	-
	//Skewedness: x1+x2=xp
	// sigma(xp)^2 = \int_0^xp dx1\int_0^xp dx2 sigma(x1)sigma(x2)delta(xp-x1-x2)=
	// = \int_0^xp sigma(x1)sigma(xp-x1)


	double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe){
	xprobe*=1.;
	double sigmap = mSigma_ep_LookupTable->Interpolate(r);
	int A=nucleus()->A();
	double TA=nucleus()->T(b)/A;
	double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) );
	return result;
	}

	void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe)
	{
	double rbRange=3.*nucleus()->radius();
	TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2);
	mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange);
	dsigmaForIntegration->SetNpx(1000);
	for (int iR=1; iR<=1000; iR++) {
	double r=mSigma_ep_LookupTable->GetBinCenter(iR);
	dsigmaForIntegration->SetParameter(0, r);
	dsigmaForIntegration->SetParameter(1, xprobe);
	double result=dsigmaForIntegration->Integral(0, rbRange);
	mSigma_ep_LookupTable->SetBinContent(iR, result);
	}
	delete dsigmaForIntegration;
	}

	double DipoleModel_bSat::dsigmadb2epForIntegration(double x, double par)
	{
	double r=par[0];
	double xprobe=par[1];
	double b = *x;
	return 2M_PIb/hbarc2*dsigmadb2ep(r, b, xprobe);
	}


	//*********bNonSat:***********************************************
	DipoleModel_bNonSat::DipoleModel_bNonSat(){}

	DipoleModel_bNonSat::~DipoleModel_bNonSat(){}

	double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe)
	{
	//
	//modify dipole size as http://arxiv.org/abs/0807.0325v1
	//
	r=mRmaxsqrt(log(1+rr/(mRmax*mRmax)));
	//const double BG = bNonSat_BG; // GeV^-2
	//double arg = bb/(2BG);
	//arg /= hbarc2;
	// double bDep= 1/(2M_PIBG) * exp(-arg);
	double bDep= protonThickness(b, r);
	double Mu02 = mMu_02; // GeV^2
	double muQ2 = mC/(r*r/hbarc2) + Mu02;
	double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;

	return omega;
	}

	double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe)
	{
	double bDep=bDependence(b, phi);
	double muQ2 = params_C/(r*r/hbarc2) + mMu_02;
	double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	double omega = ((M_PIM_PI)/Nc)(rr/hbarc2)asxg*bDep;
	return omega;
	}

	double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){
	int A=nucleus()->A();
	double TA=nucleus()->T(b)/A;
	double muQ2 = params_C/(r*r/hbarc2) + mMu_02;
	double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
	double result=ATAM_PIM_PI/Ncrr/hbarc2asxg;
	return result;
	}

	//
	////*********bCGC:***************************************************
	//
	//void DipoleModel_bCGC::createConfiguration(int iConfiguration)
	//{
	// if (!mIsInitialized) {
	// cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl;
	// exit(1);
	// }
	// //iConfiguration is a dummy for bCGC. Do this to avoid compiler warnings:
	// iConfiguration++;
	// mNucleus.generate();
	// mConfigurationExists=true;
	//}
	//
	//double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x)
	//{
	// double result=1;
	// for (unsigned int iA=0; iA<mNucleus.A(); iA++) {
	// double absdeltab=( TVector3(bcos(phi), bsin(phi), 0.)
	// -mNucleus.configuration.at(iA).position() ).Perp();
	// result=(1.-0.5dsigmadb2ep(r, absdeltab, x));
	// }
	// return 2.*(1.-result);
	//}
	//
	//double DipoleModel_bCGC::dsigmadb2ep(double r, double b, double xprobe)
	//{
	// double Y = log(1/xprobe);
	// double kappa = bCGC_kappa;
	// double N0 = bCGC_N0;
	// double x0 = bCGC_x0;
	// double lambda = bCGC_lambda;
	// double gammas = bCGC_gammas;
	// double A = -N0N0gammasgammas/((1-N0)(1-N0)*log(1-N0));
	// double B = 0.5pow(1-N0,-(1-N0)/(N0gammas));
	// double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b));
	// double rQs = r*Qs/hbarc;
	// double result=0;
	// if (rQs <= 2)
	// result = 2N0pow(0.5rQs, 2(gammas+(1/(kappalambdaY))*log(2/rQs)));
	// else
	// result = 2(1 - exp(-Alog(BrQs)log(B*rQs)));
	// return result;
	//}
	//
	//double DipoleModel_bCGC::bDependence(double b)
	//{
	// double gammas = bCGC_gammas, Bcgc = bCGC_Bcgc;
	// return exp(-0.5bb/Bcgc/gammas/hbarc2);
	//}
	//

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:54 PM (1 d, 7 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805935
Default Alt Text: (45 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions