Index: trunk/src/TableGeneratorSettings.h
===================================================================
--- trunk/src/TableGeneratorSettings.h	(revision 478)
+++ trunk/src/TableGeneratorSettings.h	(revision 479)
@@ -1,122 +1,127 @@
 //==============================================================================
 //  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 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;  
 
     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);
+
     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;         
              
     unsigned int mQ2bins;         
     unsigned int mW2bins;         
     unsigned int mTbins;         
     unsigned int mXbins;         
          
     string mBSatLookupPath;
  
     int mPriority;
-};         
+
+    bool mHasSubstructure;
+};
          
 #endif         
Index: trunk/src/Amplitudes.cpp
===================================================================
--- trunk/src/Amplitudes.cpp	(revision 478)
+++ trunk/src/Amplitudes.cpp	(revision 479)
@@ -1,334 +1,348 @@
 //==============================================================================
 //  Amplitudes.cpp
 //
 //  Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich 
 //
 //  This file is part of Sartre. 
 //
 //  This program is free software: you can redistribute it and/or modify 
 //  it under the terms of the GNU General Public License as published by 
 //  the Free Software Foundation.   
 //  This program is distributed in the hope that it will be useful, 
 //  but without any warranty; without even the implied warranty of 
 //  merchantability or fitness for a particular purpose. See the 
 //  GNU General Public License for more details. 
 //  You should have received a copy of the GNU General Public License
 //  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 //
 //  Author: Tobias Toll
 //  Last update: 
 //  $Date$
 //  $Author$
 //==============================================================================
 //#define SARTRE_IN_MULTITHREADED_MODE 1  
 #include <iostream>  
 #include <cstdio>  
 #include "Amplitudes.h"  
 #include "Constants.h"  
 #include "TableGeneratorSettings.h"  
 #include "DglapEvolution.h"  
 #include "Enumerations.h"  
 #include "Kinematics.h"  
 #include "Integrals.h"  
 #include "DipoleModel.h"  
 #if defined(SARTRE_IN_MULTITHREADED_MODE)
 #include <boost/thread.hpp>  
 #endif   
 #define PR(x) cout << #x << " = " << (x) << endl;   
   
 using namespace std;  
   
 Amplitudes::Amplitudes()
 {
     mAmplitudeT = 0;
     mAmplitudeL = 0;
     mAmplitudeT2 = 0;
     mAmplitudeL2 = 0;
     mNumberOfConfigurations = 0;
     mTheModes = 0;
     mA = 0;
     mErrorT = 0;
     mErrorL = 0;
     mErrorT2 = 0;
     mErrorL2 = 0;
 
     mAmplitudeTForSkewednessCorrection = 0;
     mAmplitudeLForSkewednessCorrection = 0;
     
     TableGeneratorSettings* settings = TableGeneratorSettings::instance();
     mNumberOfConfigurations = settings->numberOfConfigurations();
     mVerbose = settings->verbose();
     //
     // Create a vector containing instances of the Integrals class
     // and initialize them:
     //
     for (int i=0; i<=mNumberOfConfigurations; i++) {
         mIntegrals.push_back(new IntegralsExclusive);
     }
     
     mA = settings->A();
     mUPC = settings->UPC();
 
     isBNonSat = false;
     if (settings->dipoleModelType() == bNonSat)
       isBNonSat = true;
     //
     // Get the modes to calculate:
     // 0: <A> analytically <A2> averaged over configurations
     // 1: only <A> analytically
     // 2: Both <A> and <A2> averaged over configurations
     //
     mTheModes = settings->modesToCalculate();
 }
 
 Amplitudes& Amplitudes::operator=(const Amplitudes& amp)
 {  
     if (this != &amp) {
         for (unsigned int i=0; i<mIntegrals.size(); i++)
             delete mIntegrals[i];
         mIntegrals.clear();
         
         mAmplitudeT = amp.mAmplitudeT;
         mAmplitudeL = amp.mAmplitudeL;
         mAmplitudeT2 = amp.mAmplitudeT2;
         mAmplitudeL2 = amp.mAmplitudeL2;
         mNumberOfConfigurations = amp.mNumberOfConfigurations;
         mTheModes = amp.mTheModes;
         mA = amp.mA;
         mErrorT = amp.mErrorT;
         mErrorL = amp.mErrorL;
         mErrorT2 = amp.mErrorT2;
         mErrorL2 = amp.mErrorL2;
 
 	mAmplitudeTForSkewednessCorrection = amp.mAmplitudeTForSkewednessCorrection;
 	mAmplitudeLForSkewednessCorrection = amp.mAmplitudeLForSkewednessCorrection;
 
         for (unsigned int i=0; i<amp.mIntegrals.size(); i++) {    // deep copy
             mIntegrals.push_back(new IntegralsExclusive(*(amp.mIntegrals[i])));
         }
     }
     return *this;
 }  
 
 Amplitudes::Amplitudes(const Amplitudes& amp)   
 {   
     mAmplitudeT = amp.mAmplitudeT;
     mAmplitudeL = amp.mAmplitudeL;
     mAmplitudeT2 = amp.mAmplitudeT2;
     mAmplitudeL2 = amp.mAmplitudeL2;
     mErrorT = amp.mErrorT;
     mErrorL = amp.mErrorL;
     mErrorT2 = amp.mErrorT2;
     mErrorL2 = amp.mErrorL2;
     mNumberOfConfigurations = amp.mNumberOfConfigurations;
     mTheModes = amp.mTheModes;
     mA = amp.mA;
     mAmplitudeTForSkewednessCorrection = amp.mAmplitudeTForSkewednessCorrection;
     mAmplitudeLForSkewednessCorrection = amp.mAmplitudeLForSkewednessCorrection;
 
     for (unsigned int i=0; i<amp.mIntegrals.size(); i++) {    // deep copy
         mIntegrals.push_back(new IntegralsExclusive(*(amp.mIntegrals[i])));
     }
 }  
 
 Amplitudes::~Amplitudes()  
 {  
     for (unsigned int i=0; i<mIntegrals.size(); i++)
         delete mIntegrals[i];
 }  
 
 void Amplitudes::generateConfigurations()  
 {  
     for (int i=0; i<mNumberOfConfigurations; i++)
         mIntegrals[i]->dipoleModel()->createConfiguration(i);
 }  
 
 
 //void Amplitudes::calculate(double t, double Q2, double W2)
 void Amplitudes::calculate(double* kinematicPoint)
 {
     double t=0, Q2=0, W2=0, xpom=0;
     if (!mUPC){
         t  = kinematicPoint[0];
         Q2 = kinematicPoint[1];
         W2 = kinematicPoint[2];
     }
     else{
         t    = kinematicPoint[0];
         xpom = kinematicPoint[1];
     }
 #if defined(SARTRE_IN_MULTITHREADED_MODE) // multithreaded version
-    if (mA == 1) {
+    if (mA == 1 && mNumberOfConfigurations == 1) {
         cout << "Amplitudes::calculate(): Multithreaded mode (SARTRE_IN_MULTITHREADED_MODE)" << endl;
-        cout << "                         is not supported for ep (A=1). Stopping." << endl;
+        cout << "                         is not supported for ep (A=1)"<<endl;
+        cout << "                         without nucleonsubstructure(mNumberOfConfigurations=1)."<<endl;
+        cout << "                         Stopping." << endl;
         exit(1);
     }
     
     //
     //   Create a vector containing the threads:
     //
     std::vector<boost::thread*> vThreads;
     vThreads.clear();
     
     //
     //   Create the thread group:
     //
     boost::thread_group gThreads;
     if (mTheModes==0 || mTheModes == 2){
         //Start loop over configurations, each calculated on a separate thread:
         for (int i=0; i<mNumberOfConfigurations; i++){
             if (!mUPC)
                 vThreads.push_back(new boost::thread(boost::ref(*mIntegrals.at(i)),
                                                      t, Q2, W2));
             else
                 vThreads.push_back(new boost::thread(boost::ref(*mIntegrals.at(i)),
                                                      t, xpom));
             gThreads.add_thread(vThreads.at(i));
         }
     }
     
     
     //
     //   Calculate coherent cross-section according to eq.(47) in KT arXiv:hep-ph/0304189v3,
     //   this is done in the main thread in parallel with the other threads
     //   and only in eA:
     //
     if (mA>1 && (mTheModes==1 || mTheModes == 0)) {
         if (!mUPC)
             mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2);
         else
             mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom);
     }
+    if (mA == 1 && (mTheModes==1 || mTheModes == 0)) {
+        if (!mUPC)
+            mIntegrals.at(mNumberOfConfigurations)->coherentIntegralsEp(t, Q2, W2);
+        else
+            mIntegrals.at(mNumberOfConfigurations)->coherentIntegralsEp(t, xpom);
+    }
     if (mTheModes==0 || mTheModes == 2) {
         //   Wait for all threads to finish before continuing main thread:
         gThreads.join_all();
         //   Clean up threads
         vThreads.clear();
     }
 #else // unforked version
-    if ((mTheModes==0 || mTheModes == 2) || mA==1) {
+    if ((mTheModes==0 || mTheModes == 2)) {
         //Start loop over configurations:
         for (int i=0; i<mNumberOfConfigurations; i++) {
             if (!mUPC)
                 mIntegrals.at(i)->operator()(t, Q2, W2);
             else
                 mIntegrals.at(i)->operator()(t, xpom);
         }
     }
     
     //
     //  Calculate coherent cross-section according to eq.(47) in KT arXiv:hep-ph/0304189v3,
     //  (only in eA)
     //
     if (mA>1 && (mTheModes==1 || mTheModes == 0)){
         if (!mUPC)
             mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, Q2, W2);
         else
             mIntegrals.at(mNumberOfConfigurations)->coherentIntegrals(t, xpom);
     }
+    if (mA == 1 && (mTheModes==1 || mTheModes == 0)) {
+        if (!mUPC)
+            mIntegrals.at(mNumberOfConfigurations)->coherentIntegralsEp(t, Q2, W2);
+        else
+            mIntegrals.at(mNumberOfConfigurations)->coherentIntegralsEp(t, xpom);
+    }
 #endif
     
     //
     //   Calculate the resulting <A2>:
     //
     double coherentT = 0, coherentL = 0;
     double errCoherentT = 0, errCoherentL = 0;
     if ((mTheModes==0 || mTheModes == 2) || mA==1) {
         double totalT = 0;
         double totalL = 0;
         double err2TotalT = 0, err2TotalL = 0;
         double probabilityCutOff = 1e-6;
         for (int i=0; i<mNumberOfConfigurations; i++) {
             double valimT = mIntegrals.at(i)->integralImT();
             double valreT = mIntegrals.at(i)->integralReT();
             double valimL = mIntegrals.at(i)->integralImL();
             double valreL = mIntegrals.at(i)->integralReL();
             
             double errimT = mIntegrals.at(i)->errorImT();
             double errimL = mIntegrals.at(i)->errorImL();
             double errreT = mIntegrals.at(i)->errorReT();
             double errreL = mIntegrals.at(i)->errorReL();
             
             double probimT = mIntegrals.at(i)->probImT();
             double probimL = mIntegrals.at(i)->probImL();
             double probreT = mIntegrals.at(i)->probReT();
             double probreL = mIntegrals.at(i)->probReL();
             
             if (probimT > probabilityCutOff || probimL > probabilityCutOff ||
                 probreT > probabilityCutOff || probreL > probabilityCutOff){
                 if (mVerbose) {
                     cout<< "Amplitudes::calculate(): Warning, Integrals may not have reached desired precision" <<endl;
                     //Print out the largest probability:
                     probimT > probreT && probimT > probimL && probimT > probreL ?
                     cout<< "                         The probability for this is "<<probimT<<endl :
                     probreT > probimT && probreT > probimL && probreT > probreL ?
                     cout<< "                         The probability for this is "<<probreT<<endl :
                     probimL > probimT && probimL > probreT && probimL > probreL ?
                     cout<< "                         The probability for this is "<<probimL<<endl :
                     cout<< "                         The probability for this is "<<probreL<<endl;
                 }
             }
             
             //
             //  Calculate the averages
             //
             totalT += (valimT*valimT + valreT*valreT);
             totalL += (valimL*valimL + valreL*valreL);
             coherentT += valimT;
             coherentL += valimL;
             
             //
             //  ...and their errors:
             //
             // err2Total = |dtotal/dval|^2*err^2
             // dtotal/dval = 2*val
             //
             err2TotalT += (4*valimT*valimT*errimT*errimT
                            + 4*valreT*valreT*errreT*errreT);
             err2TotalL += (4*valimL*valimL*errimL*errimL
                            + 4*valreL*valreL*errreL*errreL);
             
             errCoherentT += errimT;
             errCoherentL += errimL;
         }    //for
         
         //
         //   Store the results of the second moment of the amplitudes:
         //
         mAmplitudeT2 = totalT/mNumberOfConfigurations;
         mAmplitudeL2 = totalL/mNumberOfConfigurations;
         //...and it's error:
         mErrorT2 = sqrt(err2TotalT)/mNumberOfConfigurations;
         mErrorL2 = sqrt(err2TotalL)/mNumberOfConfigurations;
     }//if (theModes)
     
     //
     //   Store the results and error of the first moment of the amplitudes:
     //
     if (mA>1 && (mTheModes==1 || mTheModes == 0)) {
         double coherentKTT = mIntegrals.at(mNumberOfConfigurations)->integralImT();
         double coherentKTL = mIntegrals.at(mNumberOfConfigurations)->integralImL();
         double errCoherentKTT = mIntegrals.at(mNumberOfConfigurations)->errorImT();
         double errCoherentKTL = mIntegrals.at(mNumberOfConfigurations)->errorImL();
         mAmplitudeT = coherentKTT;
         mAmplitudeL = coherentKTL;
         mErrorT = errCoherentKTT;
         mErrorL = errCoherentKTL;
     }
     else {
         mAmplitudeT = coherentT/mNumberOfConfigurations;
         mAmplitudeL = coherentL/mNumberOfConfigurations;
         mErrorT = errCoherentT/mNumberOfConfigurations;
         mErrorL = errCoherentL/mNumberOfConfigurations;
     }
     if (mA == 1 && mTheModes != 1 && mNumberOfConfigurations == 1){
         if (isBNonSat){
             mAmplitudeTForSkewednessCorrection = mAmplitudeT;
             mAmplitudeLForSkewednessCorrection = mAmplitudeL;
         }
         else {
             mAmplitudeTForSkewednessCorrection = mIntegrals.at(0)->integralTForSkewedness();
             mAmplitudeLForSkewednessCorrection = mIntegrals.at(0)->integralLForSkewedness();
         }
     }
 }
Index: trunk/src/Integrals.h
===================================================================
--- trunk/src/Integrals.h	(revision 478)
+++ trunk/src/Integrals.h	(revision 479)
@@ -1,187 +1,190 @@
 //==============================================================================
 //  Integrals.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$
 //==============================================================================
 //  
 //  This class calculate the integrals over the   
 //  unintegrated amplitude for \gamma-A collisions  
 //  It calculates four 4dimensional integrals for a given (t, Q2, W2):  
 //  Imaginary part for T photon  
 //  Real part for T photon  
 //  Imaginary part for L photon  
 //  Real part for L photon  
 //  
 //  The results are accessed via:  
 //    double integralImT()  
 //    double integralReT()  
 //    double integralImL()  
 //    double integralReL()  
 //  
 //==============================================================================  
 #ifndef Integrals_h  
 #define Integrals_h  
   
 class WaveOverlap;  
 class DipoleModel;  
 class TH1F;  
   
 class Integrals {  
 public:  
     Integrals();  
     virtual ~Integrals();  
     Integrals(const Integrals&);   
     Integrals& operator=(const Integrals&);   
   
     void operator() (double, double, double);  
     void operator() (double, double); //#TT UPC only takes two arguments  
       
     double integralImT() const;  
     double integralImL() const;  
     double integralReT() const;  
     double integralReL() const;  
       
     double errorImT() const;  
     double errorImL() const;  
     double errorReT() const;  
     double errorReL() const;  
   
     double probImT() const;  
     double probImL() const;  
     double probReT() const;  
     double probReL() const;
 
     double integralTForSkewedness() const;
     double integralLForSkewedness() const;
     double errorTForSkewedness() const;
     double errorLForSkewedness() const;
     
     DipoleModel* dipoleModel() const;
     DipoleModel* dipoleModelForSkewednessCorrection() const;
 
 protected:  
     virtual void calculate() = 0;  
     virtual bool setKinematicPoint(double, double, double) = 0;  
     virtual bool setKinematicPoint(double, double) = 0;  
     virtual void calculateCoherent() = 0;  
     virtual void calculateSkewedness() = 0;
     virtual void calculateEp() = 0;
       
 protected:  
     //the wave-overlap:  
     WaveOverlap* mWaveOverlap;  
     //the dipole model:  
     DipoleModel* mDipoleModel;  
     DipoleModel* mDipoleModelForSkewednessCorrection;  
       
     bool         mIsInitialized;  
     double       mRelativePrecisionOfIntegration;  
     double       mMV;  
       
     //The result from each integral:  
     double mIntegralImT;  
     double mIntegralImL;  
     double mIntegralReT;  
     double mIntegralReL;  
     //The (absolute) error from each integral:  
     double mErrorImT;  
     double mErrorImL;  
     double mErrorReT;  
     double mErrorReL;  
     //The probability for each integral:  
     double mProbImT;  
     double mProbImL;  
     double mProbReT;  
     double mProbReL;  
 
     bool   mVerbose;
     void fillZeroes();  
 
     bool   mIsUPC;
     bool   mCalculateSkewedness;
     double mIntegralTForSkewedness;
     double mIntegralLForSkewedness;
     double mErrorTForSkewedness;
     double mErrorLForSkewedness;
     double mProbImTForSkewedness;
     double mProbImLForSkewedness;    
 };  
   
   
 class IntegralsExclusive : public Integrals {  
 public:  
     IntegralsExclusive();  
     IntegralsExclusive(const IntegralsExclusive&);  
     IntegralsExclusive& operator=(const IntegralsExclusive&);  
       
     double uiAmplitudeTIm(double, double, double, double, double, double, double);  
     double uiAmplitudeLIm(double, double, double, double, double, double, double);  
     double uiAmplitudeTRe(double, double, double, double, double, double, double);  
     double uiAmplitudeLRe(double, double, double, double, double, double, double);  
     double uiAmplitudeTep(double, double, double, double, double, double);  
     double uiAmplitudeLep(double, double, double, double, double, double);  
     double uiCoherentAmplitudeT(double, double, double, double, double);  
     double uiCoherentAmplitudeL(double, double, double, double, double);  
     void coherentIntegrals(double, double, double);  
-    void coherentIntegrals(double, double);  
+    void coherentIntegrals(double, double);
+    void coherentIntegralsEp(double, double, double);
+    void coherentIntegralsEp(double, double);
+
 
     double uiAmplitudeTForSkewedness(double, double, double, double, double, double);  
     double uiAmplitudeLForSkewedness(double, double, double, double, double, double);  
     
 public:      
     double kinematicPoint[4];  
       
 private:  
     void   calculate();
     void   calculateSkewedness();
     void   calculateCoherent();
     void   calculateEp();
     bool   setKinematicPoint(double, double, double);  
     bool   setKinematicPoint(double, double);  
 
 };  
   
 inline double Integrals::integralImT() const { return mIntegralImT; }  
 inline double Integrals::integralImL() const { return mIntegralImL; }  
 inline double Integrals::integralReT() const { return mIntegralReT; }  
 inline double Integrals::integralReL() const { return mIntegralReL; }  
   
 inline double Integrals::errorImT() const { return mErrorImT; }  
 inline double Integrals::errorImL() const { return mErrorImL; }  
 inline double Integrals::errorReT() const { return mErrorReT; }  
 inline double Integrals::errorReL() const { return mErrorReL; }  
   
 inline double Integrals::probImT() const {return mProbImT; }  
 inline double Integrals::probImL() const {return mProbImL; }  
 inline double Integrals::probReT() const {return mProbReT; }  
 inline double Integrals::probReL() const {return mProbReL; }  
 
 inline double Integrals::integralTForSkewedness() const {return mIntegralTForSkewedness;}
 inline double Integrals::integralLForSkewedness() const {return mIntegralLForSkewedness;}
 
 inline double Integrals::errorTForSkewedness() const {return mErrorTForSkewedness;}
 inline double Integrals::errorLForSkewedness() const {return mErrorLForSkewedness;}
 
 
 inline DipoleModel* Integrals::dipoleModel() const { return mDipoleModel; }  
 inline DipoleModel* Integrals::dipoleModelForSkewednessCorrection()
   const { return mDipoleModelForSkewednessCorrection; }
   
 #endif  
Index: trunk/src/TableGeneratorSettings.cpp
===================================================================
--- trunk/src/TableGeneratorSettings.cpp	(revision 478)
+++ trunk/src/TableGeneratorSettings.cpp	(revision 479)
@@ -1,206 +1,213 @@
 //==============================================================================
 //  TableGeneratorSettings.cpp
 //
 //  Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich 
 //
 //  This file is part of Sartre.
 //
 //  This program is free software: you can redistribute it and/or modify 
 //  it under the terms of the GNU General Public License as published by 
 //  the Free Software Foundation.   
 //  This program is distributed in the hope that it will be useful, 
 //  but without any warranty; without even the implied warranty of 
 //  merchantability or fitness for a particular purpose. See the 
 //  GNU General Public License for more details. 
 //  You should have received a copy of the GNU General Public License
 //  along with this program.  If not, see <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(&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(&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);
+
+    
 }    
 
 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) mNumberOfConfigurations = 1;
+    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(mQ2bins*mW2bins*mTbins*mXbins) ) {
         cout << "TableGeneratorSettings::consolidateSettings(): Error, starting bin >= table! Stopping." << endl;
         exit(1);
     }
     if ( mEndBin > signed(mQ2bins*mW2bins*mTbins*mXbins) || mEndBin < 0) {
         cout << "TableGeneratorSettings::consolidateSettings(): endBin is set to table size=" << mQ2bins*mW2bins*mTbins << endl;
         mEndBin=mQ2bins*mW2bins*mTbins*mXbins;
     }
     if ( mModesToCalculate < 0 || mModesToCalculate > 2 ) {
         cout << "TableGeneratorSettings::consolidateSettings(): Error, modesToCalculate can only take values 0, 1, or 2; not "
         << mModesToCalculate << endl;
         exit(1);
     }
     
     //
     //  Make sure the W range is allowed
     //
     double VMMass=lookupPDG(mVectorMesonId)->Mass();
     double W2min=VMMass*VMMass+protonMass2;
     double Wmin=sqrt(W2min);
     if (mWmin<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::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::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; }
+
Index: trunk/src/Integrals.cpp
===================================================================
--- trunk/src/Integrals.cpp	(revision 478)
+++ trunk/src/Integrals.cpp	(revision 479)
@@ -1,753 +1,743 @@
 //==============================================================================
 //  Integrals.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: Tobias Toll
 //  Last update: 
 //  $Date$
 //  $Author$
 //==============================================================================
 #include <iostream>  
 #include <cmath>  
 #include <algorithm>  
 #include "Integrals.h"  
 #include "Constants.h"  
 #include "Nucleus.h"  
 #include "DipoleModel.h"  
 #include "AlphaStrong.h"  
 #include "Math/IntegratorMultiDim.h"  
 #include "Math/Functor.h"  
 #include "TMath.h"  
 #include "WaveOverlap.h"  
 #include "Kinematics.h"  
 #include "TableGeneratorSettings.h"  
 #include "Enumerations.h"  
 #include "IntegrandWrappers.h"  
 #include "TF1.h"  
 #include "TH1F.h"  
 #include "cuba.h"  
 
 #define PRs(x) cout << #x << " = " << scientific << (x) << endl;
 #define PR(x) cout << #x << " = " << (x) << endl;
 
 using namespace std;  
 
 
 Integrals::Integrals()   
 {   
     mIsInitialized = false;
     mRelativePrecisionOfIntegration = 0;
     mWaveOverlap = nullptr;
     mDipoleModel = nullptr;
     mDipoleModelForSkewednessCorrection = nullptr;
     mIntegralImT = 0;
     mIntegralImL = 0;
     mIntegralReT = 0;
     mIntegralReL = 0;
     mErrorImT = 0;
     mErrorImL = 0;
     mErrorReT = 0;
     mErrorReL = 0;
     mProbImT = 0;
     mProbImL = 0;
     mProbReT = 0;
     mProbReL = 0;
     
     mIntegralTForSkewedness = 0;
     mIntegralLForSkewedness = 0;
     mErrorTForSkewedness = 0;
     mErrorLForSkewedness = 0;
     
     TableGeneratorSettings* settings = TableGeneratorSettings::instance();
     
     mVerbose = settings->verbose();
     
     int VMId = settings->vectorMesonId();
     mMV = settings->lookupPDG(VMId)->Mass();
     mIsUPC = settings->UPC();
     
     if (VMId==113 || VMId==333 || VMId == 443 || VMId == 553) {
         mWaveOverlap = new WaveOverlapVM;
         mWaveOverlap->setProcess(VMId);
         mWaveOverlap->setWaveOverlapFunctionParameters(VMId);
         if (mVerbose) mWaveOverlap->testBoostedGaussianParameters(VMId);
     }
     else if (VMId==22) {
         mWaveOverlap = new WaveOverlapDVCS;
     }
     else {
         cout << "Integrals::init(): Error, no exclusive production implemented for: "<< VMId << endl;
         exit(1);
     }
     DipoleModelType model = settings->dipoleModelType();
     if (model == bSat) {
         mDipoleModel = new DipoleModel_bSat;
     }
     else if (model == bNonSat){
         mDipoleModel = new DipoleModel_bNonSat;
     }
     else if (model == bCGC) {
         mDipoleModel = new DipoleModel_bCGC;
     }
     else {
         cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
         exit(1);
     }
     mCalculateSkewedness = false;
     if (settings->A()==1 && settings->modesToCalculate()!=1 && settings->numberOfConfigurations()==1 && model==bSat) {
         mCalculateSkewedness = true;
         mDipoleModelForSkewednessCorrection = new DipoleModel_bNonSat;
     }
     mIsInitialized = true;
 }  
 
 Integrals::Integrals(const Integrals& integrals)  
 {  
     mIsInitialized = integrals.mIsInitialized;
     mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
     if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS))
         mWaveOverlap  = new WaveOverlapDVCS;
     else
         mWaveOverlap  = new WaveOverlapVM;
     if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat))
         mDipoleModel  = new DipoleModel_bSat;
     else if (typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
         mDipoleModel = new DipoleModel_bNonSat;
     else
         mDipoleModel  = new DipoleModel_bCGC;
     if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
         mDipoleModel  = new DipoleModel_bNonSat;
     mIntegralImT  = integrals.mIntegralImT;
     mIntegralImL  = integrals.mIntegralImL;
     mIntegralReT  = integrals.mIntegralReT;
     mIntegralReL  = integrals.mIntegralReL;
     mErrorImT = integrals.mErrorImT;
     mErrorImL = integrals.mErrorImL;
     mErrorReT = integrals.mErrorReT;
     mErrorReL = integrals.mErrorReL;
     mProbImT = integrals.mProbImT;
     mProbImL = integrals.mProbImL;
     mProbReT = integrals.mProbReT;
     mProbReL = integrals.mProbReL;
     mMV = integrals.mMV;
 }  
 
 Integrals& Integrals::operator=(const Integrals& integrals)  
 {  
     if (this != &integrals) {
         delete mWaveOverlap;
         delete mDipoleModel;
         delete mDipoleModelForSkewednessCorrection;
         
         if (typeid(*mWaveOverlap) == typeid(WaveOverlapDVCS))
             mWaveOverlap  = new WaveOverlapDVCS;
         else
             mWaveOverlap  = new WaveOverlapVM;
         if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat))
             mDipoleModel  = new DipoleModel_bSat;
         else if (typeid(*mDipoleModel) == typeid(DipoleModel_bNonSat))
             mDipoleModel  = new DipoleModel_bNonSat;
         else
             mDipoleModel  = new DipoleModel_bCGC;
 	if (typeid(*mDipoleModelForSkewednessCorrection) == typeid(DipoleModel_bNonSat))
 	    mDipoleModel  = new DipoleModel_bNonSat;
 	
         mIntegralImT  = integrals.mIntegralImT;
         mIntegralImL  = integrals.mIntegralImL;
         mIntegralReT  = integrals.mIntegralReT;
         mIntegralReL  = integrals.mIntegralReL;
         mIntegralImT  = integrals.mIntegralImT;
         mIntegralImL  = integrals.mIntegralImL;
         mIntegralReT  = integrals.mIntegralReT;
         mIntegralReL  = integrals.mIntegralReL;
         mErrorImT = integrals.mErrorImT;
         mErrorImL = integrals.mErrorImL;
         mErrorReT = integrals.mErrorReT;
         mErrorReL = integrals.mErrorReL;
         mProbImT = integrals.mProbImT;
         mProbImL = integrals.mProbImL;
         mProbReT = integrals.mProbReT;
         mProbReL = integrals.mProbReL;
         mMV = integrals.mMV;
         mIsInitialized = integrals.mIsInitialized;
         mRelativePrecisionOfIntegration = integrals.mRelativePrecisionOfIntegration;
     }
     return *this;
 }  
 
 Integrals::~Integrals()   
 {   
     delete mWaveOverlap;
     delete mDipoleModel;
     if (mDipoleModelForSkewednessCorrection)
       delete mDipoleModelForSkewednessCorrection;
 }  
 
 void Integrals::operator() (double t, double Q2, double W2)  
 {  
     unsigned int A = dipoleModel()->nucleus()->A();
+    TableGeneratorSettings* settings = TableGeneratorSettings::instance();
+    unsigned int numberOfConfigurations=settings->numberOfConfigurations();
     //make sure the configurations have been generated:
-    if (!mDipoleModel->configurationExists() && A!=1) {
+    if (!mDipoleModel->configurationExists() && (A!=1 || (A==1 && numberOfConfigurations>1)) ) {
         // do not use cout
         cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
         exit(1);
     }
     if (setKinematicPoint(t, Q2, W2)) {
-  	if (A==1){
-	    calculateEp();	  
-	    if (mCalculateSkewedness){
-	        calculateSkewedness();
-	    }
-	}
-	else
-	    calculate();
+        if (A==1 && numberOfConfigurations == 1){
+            calculateEp();
+            if (mCalculateSkewedness){
+                calculateSkewedness();
+            }
+        }
+        else
+            calculate();
     }
     else {
         fillZeroes();
     }
 }
 
 void Integrals::operator() (double t, double xpom) //UPC
 {
+    TableGeneratorSettings* settings = TableGeneratorSettings::instance();
+    unsigned int numberOfConfigurations=settings->numberOfConfigurations();
     unsigned int A = dipoleModel()->nucleus()->A();
     //make sure the configurations have been generated:
-    if (!mDipoleModel->configurationExists() && A!=1) {
+    if (!mDipoleModel->configurationExists() && (A!=1 || (A==1 && numberOfConfigurations>1)) ) {
         // do not use cout
         cout << "Integrals::init(): Error, configuration has not been generated. Stopping." << endl;
         exit(1);
     }
     if (setKinematicPoint(t, xpom)) {
-  	if (A==1){
-	    calculateEp();	  
-	    if (mCalculateSkewedness){
-	        calculateSkewedness();
-	    }
-	}
-	else
-  	    calculate();
+        if (A==1 && numberOfConfigurations == 1){
+            calculateEp();
+            if (mCalculateSkewedness){
+                calculateSkewedness();
+            }
+        }
+        else
+            calculate();
     }
     else {
         fillZeroes();
     }
 }
 
 void Integrals::fillZeroes(){  
     //Store the results
     mIntegralImT = 0;
     mIntegralReT = 0;
     mIntegralImL = 0;
     mIntegralReL = 0;
     
     //Store the errors:
     mErrorImT = 0;
     mErrorImL = 0;
     mErrorReT = 0;
     mErrorReL = 0;
     
     //Store the probabilities:
     mProbImT = 0;
     mProbImL = 0;
     mProbReT = 0;
     mProbReL = 0;
     
 }  
 
 //*********EXCLUSIVE VECTOR MESONS OR DVCS: ********************************  
 
 void IntegralsExclusive::coherentIntegrals(double t, double Q2, double W2)  
 {  
     if (setKinematicPoint(t, Q2, W2)){
         if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){
             //store present kinematic point:
             double xprobe = kinematicPoint[3];
             dipoleModel()->createSigma_ep_LookupTable(xprobe);
         }
         calculateCoherent();
     }
     else
         fillZeroes();
 }  
 
 void IntegralsExclusive::coherentIntegrals(double t, double xpom)  
 {  
     if (setKinematicPoint(t, xpom)){
         if (typeid(*mDipoleModel) == typeid(DipoleModel_bSat)){
             dipoleModel()->createSigma_ep_LookupTable(xpom);
         }
         calculateCoherent();
     }
     else
         fillZeroes();
 }  
 
+void IntegralsExclusive::coherentIntegralsEp(double t, double Q2, double W2)
+{
+    if (setKinematicPoint(t, Q2, W2))
+        calculateEp();
+    else
+        fillZeroes();
+}
+
+void IntegralsExclusive::coherentIntegralsEp(double t, double xpom)
+{
+    if (setKinematicPoint(t, xpom))
+        calculateEp();
+    else
+        fillZeroes();
+}
+
 IntegralsExclusive::IntegralsExclusive()   
 {  
 }  
 
 IntegralsExclusive& IntegralsExclusive::operator=(const IntegralsExclusive& cobj)  
 {    
     if (this != &cobj) {
         Integrals::operator=(cobj);
         copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
     }
     return *this;
 }  
 
 IntegralsExclusive::IntegralsExclusive(const IntegralsExclusive& cobj) : Integrals(cobj)  
 {  
     copy(cobj.kinematicPoint, cobj.kinematicPoint+4, kinematicPoint);
 }  
 
 
 bool IntegralsExclusive::setKinematicPoint(double t, double xpom) //UPC
 {
     bool result = true;
     kinematicPoint[0] = t;
     kinematicPoint[1] = 0; //Q2
     kinematicPoint[2] = 0; //W2 is not used
     kinematicPoint[3] = xpom;
     
     return result;
 }  
 
 bool IntegralsExclusive::setKinematicPoint(double t, double Q2, double W2)
 {
     bool result = true;
     kinematicPoint[0] = t;
     kinematicPoint[1] = Q2;
     kinematicPoint[2] = W2;
     double xprobe = Kinematics::xpomeron(t, Q2, W2, mMV);
     
     if (xprobe<0 || xprobe>1)
         result = false;
     kinematicPoint[3] = xprobe;
     return result;
 }  
 
 
 void IntegralsExclusive::calculate()  
 {  
     //
     // This function calls a wrapper from where the
     // integral is calculated with the Cuhre method.
     // Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
     //
     const double epsrel = 1.e-2, epsabs = 1e-12;
     const int flags=0, mineval=3e6, maxeval=1e9, key=0;
     int nregionsTIm, nevalTIm, failTIm;
     int nregionsTRe, nevalTRe, failTRe;
     int nregionsLIm, nevalLIm, failLIm;
     int nregionsLRe, nevalLRe, failLRe;
     double valTIm=0, errTIm=0, probTIm=0;
     double valLIm=0, errLIm=0, probLIm=0;
     double valTRe=0, errTRe=0, probTRe=0;
     double valLRe=0, errLRe=0, probLRe=0;
     
     const char* statefile=0;
     
     const int nvec=1;
 
     //    double probabilityCutOff=1e-6;
     
     //
     //   Do the integrations
     //
     Cuhre(4, 1, integrandWrapperTIm, this, nvec,
           epsrel, epsabs, flags,
           mineval, maxeval, key, statefile, 0, &nregionsTIm, &nevalTIm, &failTIm, &valTIm, &errTIm, &probTIm);
     if (failTIm != 0 && mVerbose) {
         cout << "IntegralsExclusive::calculate(): Warning: Integration TIm did not reach desired precision! Error code=" << failTIm << endl;
     }
 
     //
     // For UPC, calculate only transverse polarisation case
     //
     if (!mIsUPC) {
         Cuhre(4, 1, integrandWrapperLIm, this, nvec,
               epsrel, epsabs, flags,
               mineval, maxeval, key, statefile, 0, &nregionsLIm, &nevalLIm, &failLIm, &valLIm, &errLIm, &probLIm);
         if (failLIm!=0 && mVerbose) {
             cout << "IntegralsExclusive::calculate(): Warning: Integration LIm did not reach desired precision! Error code=" << failLIm << endl;
         }
     }
     Cuhre(4, 1, integrandWrapperTRe, this, nvec,
           epsrel, epsabs, flags,
           mineval, maxeval, key, statefile, 0, &nregionsTRe, &nevalTRe, &failTRe, &valTRe, &errTRe, &probTRe);
     if (failTRe!=0 && mVerbose) {
         cout << "IntegralsExclusive::calculate(): Warning: Integration TRe did not reach desired precision! Error code=" << failTRe << endl;
     }
     
     //
     // For UPC, calculate only transverse polarisation case
     //
     if (!mIsUPC) {
         Cuhre(4, 1, integrandWrapperLRe, this, nvec,
               epsrel, epsabs, flags,
               mineval, maxeval, key, statefile, 0, &nregionsLRe, &nevalLRe, &failLRe, &valLRe, &errLRe, &probLRe);
         if (failLRe!=0 && mVerbose) {
             cout << "IntegralsExclusive::calculate(): Warning: Integration LRe did not reach desired precision! Error code=" << failLRe << endl;
         }
     }
     
     //
     //   Store the results:
     //
     mIntegralImT = valTIm;
     mIntegralReT = valTRe;
     mIntegralImL = valLIm;
     mIntegralReL = valLRe;
     
     //
     //   Store the errors:
     //
     mErrorImT = errTIm;
     mErrorImL = errLIm;
     mErrorReT = errTRe;
     mErrorReL = errLRe;
     
     //
     //   Store the probabilities:
     //
     mProbImT = probTIm;
     mProbImL = probLIm;
     mProbReT = probTRe;
     mProbReL = probLRe;
 }
 
 void IntegralsExclusive::calculateEp()
 {
     //
     // This function calls a wrapper from where the
     // integral is calculated with the Cuhre method.
     // Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
     //
     const double epsrel = 1.e-4, epsabs = 1e-12;
     const int flags=0, maxeval=1e9, key=0;
     const int mineval=3e6;
     int nregionsT, nevalT, failT;
     int nregionsL, nevalL, failL;
     double valT=0, errT=0, probT=0;
     double valL=0, errL=0, probL=0;
     
     const char* statefile=0;
     
     const int nvec=1;
     //
     //   Do the integrations
     //
-    /*
-    bool bContinue=true;
-    bool isFirst=true;
-    while(bContinue){
-      double valTOld=valT;
-      Cuhre(4, 1, integrandWrapperTep, this, nvec,
-	    epsrel, epsabs, flags,
-	    mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
-      if (abs(valT-valTOld)/valT > epsrel){
-	mineval*=3;
-	if (isFirst)
-	  mineval*=30;
-	isFirst=false;
-	if (mineval>1e4)
-	  PR(mineval);
-      }
-      else
-	bContinue=false;
-    }
-    */
     Cuhre(4, 1, integrandWrapperTep, this, nvec,
           epsrel, epsabs, flags,
           mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
     if (failT!=0 && mVerbose) {
         cout << "IntegralsExclusive::calculateEp(): Warning: Integration T did not reach desired precision! Error code=" << failT << endl;
     }
-    /*
-    if (errT/valT > epsrel)
-      PR(errT/valT);
-    if (errT < epsabs)
-      PR(errT);
-    if (probT>0.5)
-      PR(probT);
-    PR(nevalT);
-    PR(nregionsT);
-    */
     //
     // For UPC, calculate only transverse polarisation case
     //
     if (!mIsUPC){
         Cuhre(4, 1, integrandWrapperLep, this, nvec,
 	      epsrel, epsabs, flags,
 	      mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
         if (failL!=0 && mVerbose) {
             cout << "IntegralsExclusive::calculateEp(): Warning: Integration L did not reach desired precision! Error code=" << failL << endl;
         }
     }
     //
     //   Store the results:
     //
     mIntegralImT=valT;
     mIntegralImL=valL;
     
     //
     //   Store the errors:
     //
     mErrorImT=errT;
     mErrorImL=errL;
     
     //
     //   Store the probabilities:
     //
     mProbImT=probT;
     mProbImL=probL;
 }
 
 void IntegralsExclusive::calculateSkewedness()
 {
     //
     // This function calls a wrapper from where the
     // integral is calculated with the Cuhre method.
     // Pass this Integrals object as the fourth (void*) argument of the Cuhre function.
     //
     const double epsrel=1.e-4, epsabs=1e-12;
     const int flags=0, maxeval=1e9, key=0;
     const int  mineval=3e6;
     int nregionsT, nevalT, failT;
     int nregionsL, nevalL, failL;
     double valT=0, errT=0, probT=0;
     double valL=0, errL=0, probL=0;
 
     const char* statefile=0;
     
     const int nvec=1;
 
     //
     //   Do the integrations
     //
     /*
     bool bContinue=true;
     bool isFirst=true;
     while(bContinue){
       double valTOld=valT;
       Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
 	    epsrel, epsabs, flags,
 	    mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
       if (abs(valT-valTOld)/valT > epsrel){
 	mineval*=3;
 	if (isFirst)
 	  mineval*=30;
 	isFirst=false;
 	if (mineval>1e4)
 	  PR(mineval);
       }
       else
 	bContinue=false;
     }
     */
     Cuhre(4, 1, integrandWrapperTForSkewedness, this, nvec,
 	  epsrel, epsabs, flags,
 	  mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
     
     if (failT!=0 && mVerbose) {
         cout << "IntegralsExclusive::calculateSkweedness(): Warning: Integration T did not reach desired precision! Error code" << failT << endl;
     }
     //
     // For UPC, calculate only transverse polarisation case
     //
     if (!mIsUPC) {
         Cuhre(4, 1, integrandWrapperLForSkewedness, this, nvec,
               epsrel, epsabs, flags,
               mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
         if (failL!=0 && mVerbose) {
             cout << "IntegralsExclusive::calculateSkweedness(): Warning: Integration L did not reach desired precision! Error code" << failL << endl;
         }
     }
     //
     //   Store the results:
     //
     mIntegralTForSkewedness=valT;
     mIntegralLForSkewedness=valL;
     //
     //   Store the errors:
     //
     mErrorTForSkewedness=errT;
     mErrorLForSkewedness=errL;
     //
     //   Store the probabilities:
     //
     mProbImTForSkewedness=probT;
     mProbImLForSkewedness=probL;
 }
 
 void IntegralsExclusive::calculateCoherent()  
 {
     //
     //  As calculate() but for the coherent case
     //
     const double epsrel=1e-5, epsabs=1e-12;
     const int flags=0, mineval=1e1, maxeval=1e9, key=0;
     int nregionsT, nevalT, failT;
     int nregionsL, nevalL, failL;
     double valT, errT, probT;
     double valL=0, errL=0, probL=0;
     const int nvec=1;
     const char* statefile=0;
 
     //
     //   Do the integrations
     //
     Cuhre(3, 1, integrandWrapperCoherentAmplitudeT, this, nvec,
           epsrel, epsabs, flags,
           mineval, maxeval, key, statefile, 0, &nregionsT, &nevalT, &failT, &valT, &errT, &probT);
     if (failT!=0 && mVerbose) {
         cout << "IntegralsExclusive::calculateCoherent(): Warning: Integration T did not reach desired precision! Error code=" <<  failT << endl;
     }
     //
     // For UPC, calculate only transverse polarisation case
     //
     if (!mIsUPC){
         Cuhre(3, 1, integrandWrapperCoherentAmplitudeL, this, nvec,
               epsrel, epsabs, flags,
               mineval, maxeval, key, statefile, 0, &nregionsL, &nevalL, &failL, &valL, &errL, &probL);
         if (failL!=0 && mVerbose) {
             cout << "IntegralsExclusive::calculateCoherent(): Warning: Integration L did not reach desired precision! Error=" << failL << endl;
         }
     }
     //
     //   Store the results
     //
     mIntegralImT=valT;
     mIntegralImL=valL;
     
     //
     //   Store the errors
     //
     mErrorImT=errT;
     mErrorImL=errL;
 }
 
 //  
 //   The following functions are the Integrands in the Amplitudes:
 //  
 double IntegralsExclusive::uiAmplitudeTIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)  
 {  
     double cosArg = (b/hbarc)*Delta*cos(phi);
     double waveOverlap = mWaveOverlap->T(z, Q2, r);
     double dsigdb2 = dipoleModel()->dsigmadb2(r  , b  , phi, xprobe);
     double BesselJ0 = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*
     BesselJ0*b*cos(cosArg)*dsigdb2;
     return result;
 }  
 
 double IntegralsExclusive::uiAmplitudeTRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)  
 {  
     double sinArg = b*Delta*cos(phi)/hbarc;
     double waveOverlap = mWaveOverlap->T(z, Q2, r);
     double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
     double BesselJ0 = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*
     BesselJ0*b*sin(sinArg)*dsigdb2;
     return result;
 }  
 double IntegralsExclusive::uiAmplitudeLIm(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->L(z, Q2, r);
     double cosArg = b*Delta*cos(phi)/hbarc;
     double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
     double BesselJ0 = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*
     BesselJ0*b*cos(cosArg)*dsigdb2;
     return result;
 }  
 
 double IntegralsExclusive::uiAmplitudeLRe(double b, double z, double r, double phi, double Q2, double xprobe, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->L(z, Q2, r);
     double sinArg = b*Delta*cos(phi)/hbarc;
     double dsigdb2 = dipoleModel()->dsigmadb2(r, b, phi, xprobe);
     double BesselJ0 = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*
     BesselJ0*b*sin(sinArg)*dsigdb2;
     return result;
 }  
 
 double IntegralsExclusive::uiCoherentAmplitudeT(double b, double z, double r, double Q2, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->T(z, Q2, r);
     double BesselJ0r = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
     double xprobe=kinematicPoint[3];
     double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe);
     double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean;
     return result;
 }  
 
 double IntegralsExclusive::uiCoherentAmplitudeL(double b, double z, double r, double Q2, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->L(z, Q2, r);
     double BesselJ0r = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
     double xprobe=kinematicPoint[3];
     double dsigmadb2Mean=dipoleModel()->coherentDsigmadb2(r, b, xprobe);
     double result = M_PI*r*b/hbarc2*waveOverlap*BesselJ0r*BesselJ0b*dsigmadb2Mean;
     return result;
 }  
 
 double IntegralsExclusive::uiAmplitudeTep(double b, double z, double r, double Q2, double xprobe, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->T(z, Q2, r);
     double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
     double BesselJ0r = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
     result*=2*M_PI; 
     return result;
 }
 double IntegralsExclusive::uiAmplitudeLep(double b, double z, double r, double Q2, double xprobe, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->L(z, Q2, r);
     double dsigdb2 = dipoleModel()->dsigmadb2ep(r , b, xprobe);
     double BesselJ0r = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
     result*=2*M_PI; 
     return result;
 }
 
 //
 // Only for calculating the lamdba for Skewedness Corrections:
 //
 double IntegralsExclusive::uiAmplitudeTForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->T(z, Q2, r);
     double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r , b, xprobe);
     double BesselJ0r = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
     result*=2*M_PI; 
     return result;
 }
 
 double IntegralsExclusive::uiAmplitudeLForSkewedness(double b, double z, double r, double Q2, double xprobe, double Delta)  
 {  
     double waveOverlap = mWaveOverlap->L(z, Q2, r);
     double dsigdb2 = dipoleModelForSkewednessCorrection()->dsigmadb2ep(r, b, xprobe);
     double BesselJ0r = TMath::BesselJ0((0.5-z)*r*Delta/hbarc);
     double BesselJ0b = TMath::BesselJ0(b*Delta/hbarc);
     double result=0.5*r/hbarc2*waveOverlap*BesselJ0r*b*BesselJ0b*dsigdb2;
     result*=2*M_PI; 
     return result;
 }