Index: trunk/src/Amplitudes.cpp
===================================================================
--- trunk/src/Amplitudes.cpp	(revision 383)
+++ trunk/src/Amplitudes.cpp	(revision 384)
@@ -1,334 +1,334 @@
 //==============================================================================
 //  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) {
         cout << "Amplitudes::calculate(): Multithreaded mode (SARTRE_IN_MULTITHREADED_MODE)" << endl;
         cout << "                         is not supported for ep (A=1). 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 (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) {
         //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);
     }
 #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 and probimT > probimL and probimT > probreL ?
                     cout<< "                         The probability for this is "<<probimT<<endl :
                     probreT > probimT and probreT > probimL and probreT > probreL ?
                     cout<< "                         The probability for this is "<<probreT<<endl :
                     probimL > probimT and probimL > probreT and 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;
+        mAmplitudeL=coherentL/mNumberOfConfigurations;
         mErrorT=errCoherentT/mNumberOfConfigurations;
         mErrorL=errCoherentL/mNumberOfConfigurations;
     }
-    if(mA==1 and mTheModes != 1 and mNumberOfConfigurations == 1){
-      if(isBNonSat){
-	mAmplitudeTForSkewednessCorrection = mAmplitudeT;
-	mAmplitudeLForSkewednessCorrection = mAmplitudeL;
-      }
-      else{
-	mAmplitudeTForSkewednessCorrection = mIntegrals.at(0)->integralTForSkewedness();
-	mAmplitudeLForSkewednessCorrection = mIntegrals.at(0)->integralLForSkewedness();
-      }
+    if (mA==1 and mTheModes != 1 and mNumberOfConfigurations == 1){
+        if (isBNonSat){
+            mAmplitudeTForSkewednessCorrection = mAmplitudeT;
+            mAmplitudeLForSkewednessCorrection = mAmplitudeL;
+        }
+        else {
+            mAmplitudeTForSkewednessCorrection = mIntegrals.at(0)->integralTForSkewedness();
+            mAmplitudeLForSkewednessCorrection = mIntegrals.at(0)->integralLForSkewedness();
+        }
     }
 }
Index: trunk/src/TableGeneratorSettings.cpp
===================================================================
--- trunk/src/TableGeneratorSettings.cpp	(revision 383)
+++ trunk/src/TableGeneratorSettings.cpp	(revision 384)
@@ -1,210 +1,210 @@
 //==============================================================================
 //  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);
     
 }    
 
 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) and 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 (!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, modeToCalculate can only take values 0, 1, or 2; not "
+        cout << "TableGeneratorSettings::consolidateSettings(): Error, modesToCalculate can only take values 0, 1, or 2; not "
         << mModesToCalculate << endl;
         exit(1);
     }
-    if ( (mModesToCalculate == 1) && (mA == 1) ){
-        cout << "TableGeneratorSettings::consolidateSettings(): Error, modeToCalculate and A cannot both be 1 simultaneously. This setting will produce NO tables." << endl;
+    if ( (mModesToCalculate == 1) && (mA == 1) ) {
+        cout << "TableGeneratorSettings::consolidateSettings(): Error, modesToCalculate and A cannot both be 1 simultaneously. This setting will produce NO tables." << 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; }
Index: trunk/src/Integrals.cpp
===================================================================
--- trunk/src/Integrals.cpp	(revision 383)
+++ trunk/src/Integrals.cpp	(revision 384)
@@ -1,745 +1,745 @@
 //==============================================================================
 //  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 PRf(x) printf(#x); printf("=%g \n", (x));
 #define PRi(x) printf(#x); printf("=%d \n", (x));
 #define PR(x) cout << #x << " = " << (x) << endl;
 
 using namespace std;  
 
 
 Integrals::Integrals()   
 {   
     mIsInitialized=false;
     mRelativePrecisionOfIntegration = 0;
     mWaveOverlap = 0;
     mDipoleModel = 0;
     mDipoleModelForSkewednessCorrection = 0;
     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 and settings->modesToCalculate()!=1 and settings->numberOfConfigurations()==1 and model==bSat) {
-      mCalculateSkewedness=true;
-      mDipoleModelForSkewednessCorrection = new DipoleModel_bNonSat;
+    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();
     //make sure the configurations have been generated:
     if (!mDipoleModel->configurationExists() and A!=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();
     }
     else {
         fillZeroes();
     }
 }
 
 void Integrals::operator() (double t, double xpom) //UPC
 {
     unsigned int A=dipoleModel()->nucleus()->A();
     //make sure the configurations have been generated:
     if (!mDipoleModel->configurationExists() and A!=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();
     }
     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();
 }  
 
 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=1e4, 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 and mVerbose)
         printf("IntegralsExclusive::calculate(): Warning: Integration TIm did not reach desired precision! Error code=%d \n", failTIm);
 
     //
     // 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 and mVerbose)
         printf("IntegralsExclusive::calculate(): Warning: Integration LIm did not reach desired precision! Error code=%d \n", failLIm);
     }
     Cuhre(4, 1, integrandWrapperTRe, this, nvec,
 	  epsrel, epsabs, flags,
 	  mineval, maxeval, key, statefile, 0, &nregionsTRe, &nevalTRe, &failTRe, &valTRe, &errTRe, &probTRe);
     if(failTRe!=0 and mVerbose)
       printf("IntegralsExclusive::calculate(): Warning: Integration TRe did not reach desired precision! Error code=%d \n", failTRe);
         
     //
     // 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 and mVerbose)
 	printf("IntegralsExclusive::calculate(): Warning: Integration LRe did not reach desired precision! Error code=%d \n", failLRe);
     }
     
     //
     //   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 and mVerbose)
         printf("IntegralsExclusive::calculateEp(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
     /*
     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 and mVerbose)
   	    printf("IntegralsExclusive::calculateEp(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
     }
     //
     //   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 and mVerbose)
         printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
     //
     // 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 and mVerbose)
 	  printf("IntegralsExclusive::calculateSkweedness(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
     }
     //
     //   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-6, 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 and mVerbose)
         printf("IntegralsExclusive::calculateCoherent(): Warning: Integration T did not reach desired precision! Error code=%d \n", failT);
     //
     // 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 and mVerbose)
             printf("IntegralsExclusive::calculateCoherent(): Warning: Integration L did not reach desired precision! Error code=%d \n", failL);
     }
     //
     //   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((1-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((1-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((1-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((1-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((1-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((1-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((1-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((1-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((1-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((1-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;
 }