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.*mnLP*hbarc; //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)
         // ~2pi*sum(T'(bi)*bi*Delta b)
         //
         norm+=2.*M_PI*b/hbarc*LP_of_b*binWidth/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/(4*M_PI*mnLP*mnLP))*0.5/(mnLP*hbarc)*TMath::BesselK1(0.5/(mnLP*hbarc)));
     // PR((1/(4*M_PI*mnLP*mnLP))*1/(mnLP*hbarc)*TMath::BesselK1(1/(mnLP*hbarc)));
     // PR((1/(4*M_PI*mnLP*mnLP))*2/(mnLP*hbarc)*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+=b*b*b*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((z*z+b*b))/(hbarc*r));
     return (1/(8*M_PI*r*r*r))*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)
             // ~2pi*sum(T'(bi)*bi*Delta b)
             //
             norm+=2.*M_PI*b/hbarc*CE_of_b*binWidthb/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)/(4*M_PI*TMath::Gamma(3/mEP))));
     
     // PR((1/pow(2*M_PI*mBg, 1.5)));
     
     // // double u=muCE;
     // // double v=mEP;
     // // double z=0;
     // // double b=0.5;
     // // double arg=u*(pow((sqrt(z*z+b*b)/hbarc),v));
     // // PR(((pow(u,3/v)*v)/(4*M_PI*TMath::Gamma(3/v)))*exp(-arg));
     // // b=1;
     // // arg=u*(pow((sqrt(z*z+b*b)/hbarc),v));
     // // PR(((pow(u,3/v)*v)/(4*M_PI*TMath::Gamma(3/v)))*exp(-arg));
     // // b=2;
     // // arg=u*(pow((sqrt(z*z+b*b)/hbarc),v));
     // // PR(((pow(u,3/v)*v)/(4*M_PI*TMath::Gamma(3/v)))*exp(-arg));
     
     // PR((1/pow(2*M_PI*mBg, 1.5)) * exp(-(0.5*0.5)/(hbarc2*2*mBg)));
     // PR((1/pow(2*M_PI*mBg, 1.5)) * exp(-(1*1)/(hbarc2*2*mBg)));
     // PR((1/pow(2*M_PI*mBg, 1.5)) * exp(-(2*2)/(hbarc2*2*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/(4*M_PI*mnLP*mnLP))*0.5/(mnLP*hbarc)*TMath::BesselK1(0.5/(mnLP*hbarc)));
     // // PR((1/(4*M_PI*mnLP*mnLP))*1/(mnLP*hbarc)*TMath::BesselK1(1/(mnLP*hbarc)));
     // // PR((1/(4*M_PI*mnLP*mnLP))*2/(mnLP*hbarc)*TMath::BesselK1(2/(mnLP*hbarc)));
     
     // // PR(1/(2*M_PI*mBg) * exp(-0*0/((hbarc2)*(2*mBg))));
     // PR((1/pow(2*M_PI*mBg, 1)) * exp(-(0.5*0.5)/(hbarc2*2*mBg)));
     // PR((1/pow(2*M_PI*mBg, 1)) * exp(-(1*1)/(hbarc2*2*mBg)));
     // PR((1/pow(2*M_PI*mBg, 1)) * exp(-(2*2)/(hbarc2*2*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+=b*b*b*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(-exponentparameter*r*r);
     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(z*z+b*b)/hbarc),v));
     double arg=u*(pow((z*z+b*b)/hbarc2, 0.5*v));
     return pow(u,3./v)*v / (4.*M_PI*TMath::Gamma(3./v)) * exp(-arg); //GeV^
 }
 
 void DipoleModel::createWoodsSaxonThickness(){
     double zbhigh=3.*params_R0WS*hbarc; //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)
         // ~2pi*sum(T'(bi)*bi*Delta b)
         //
         norm+=2.*M_PI*b/hbarc*WS_of_b*binWidth/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+=b*b*b*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((z*z+b*b))/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_PI*M_PI)/Nc)*(r*r/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 = b*b/(2*BG);
      arg /= hbarc2;
      double bDep= 1/(2*M_PI*BG) * 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_PI*M_PI)/Nc)*(r*r/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 = omegaFactor*ddx_xg*exp(-omegaFactor*xprobe*Gofx/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/2pi*1/R^3*sqrt(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_Rhs*params_Rhs-b*b/hbarc2)/(M_PI*pow(params_Rhs, 3));
             break;
         case Laplace:
             result = (1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*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 = b*b/(2*BG); //fm2*GeV^2
             arg /= hbarc2; //GeV^0
             result= 1/(2*M_PI*BG) * 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./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg);
         double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*TMath::BesselK1(b/(mnLP*hbarc));
         double T=(1-w)*T_gauss+w*T_laplace;
         //        T=0.5*(T+T_old);
         //        T_old=T;
         norm+=2*M_PI*b*T/hbarc2*step;
     }
     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./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg);
             double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*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+weight*T_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+=2*M_PI*b*bDep*db/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=mRmax*sqrt(log(1+r*r/(mRmax*mRmax)));
+    double T_gauss=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg);
+    double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*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+weight*T_laplace;
+        double omega = ((M_PI*M_PI)/Nc)*(r*r/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_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg;
+    bDep/=mThicknessNormTable->Interpolate(log10(X));
+    dsigdb2 = 2.*(1. - exp(-X*bDep/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=mRmax*sqrt(log(1+r*r/(mRmax*mRmax)));
-//    mBg=4;
-//    mnLP=1.8282;
-//    double T_gauss=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg);
-//    double T_laplace=(1/(4*M_PI*mnLP*mnLP))*b/(mnLP*hbarc)*TMath::BesselK1(b/(mnLP*hbarc));
+//    double bDep=1./(2*M_PI*mBg)*exp(-b*b/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+weight*T_laplace;
-//        double omega = ((M_PI*M_PI)/Nc)*(r*r/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_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg;
-//    //   bDep/=mThicknessNormTable->Interpolate(log10(X));
-//    dsigdb2 = 2.*(1. - exp(-X*bDep/2));
-//
-//    return bDep;
+//    double omega = ((M_PI*M_PI)/Nc)*(r*r/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=mRmax*sqrt(log(1+r*r/(mRmax*mRmax)));
-    double bDep=1./(2*M_PI*mBg)*exp(-b*b/hbarc2/2./mBg);
-    double muQ2 = mC/(r*r/hbarc2) + mMu_02; //GeV2
-    double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
-    double omega = ((M_PI*M_PI)/Nc)*(r*r/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 2*M_PI*b/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=mRmax*sqrt(log(1+r*r/(mRmax*mRmax)));
     //const double BG = bNonSat_BG; // GeV^-2
     //double arg = b*b/(2*BG);
     //arg /= hbarc2;
     // double bDep= 1/(2*M_PI*BG) * 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_PI*M_PI)/Nc)*(r*r/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_PI*M_PI)/Nc)*(r*r/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=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
     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(b*cos(phi), b*sin(phi), 0.)
 //                          -mNucleus.configuration.at(iA).position() ).Perp();
 //        result*=(1.-0.5*dsigmadb2ep(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 = -N0*N0*gammas*gammas/((1-N0)*(1-N0)*log(1-N0));
 //    double B = 0.5*pow(1-N0,-(1-N0)/(N0*gammas));
 //    double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b));
 //    double rQs = r*Qs/hbarc;
 //    double result=0;
 //    if (rQs <= 2)
 //        result = 2*N0*pow(0.5*rQs, 2*(gammas+(1/(kappa*lambda*Y))*log(2/rQs)));
 //    else
 //        result = 2*(1 - exp(-A*log(B*rQs)*log(B*rQs)));
 //    return result;
 //}
 //
 //double DipoleModel_bCGC::bDependence(double b)
 //{
 //    double gammas = bCGC_gammas, Bcgc = bCGC_Bcgc;
 //    return exp(-0.5*b*b/Bcgc/gammas/hbarc2);
 //}
 //