Index: trunk/examples/lambdaGeneratorMain.cpp
===================================================================
--- trunk/examples/lambdaGeneratorMain.cpp	(revision 375)
+++ trunk/examples/lambdaGeneratorMain.cpp	(revision 376)
@@ -1,192 +1,192 @@
 //==============================================================================
 //  lambdaGeneratorMain.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. See <http://www.gnu.org/licenses/>. 
 //
 //  Author: Thomas Ullrich
 //  Last update: Sat Aug 13 02:18:37 2011
 //==============================================================================
 // 
 //  Main program to create lambda lookup tables.
 //  [Developer only]
 //============================================================================== 
 #include <sstream>
 #include <cstdlib>
 #include <iostream>
 #include <iomanip>
 #include "Amplitudes.h"
 #include "TH1D.h"
 #include "TFile.h"
 #include "Constants.h"
 #include "Table.h"
 #include "Settings.h"
 #include "TableGeneratorSettings.h"
 #include "Enumerations.h"
 #include "Kinematics.h"
 #include <limits> 
 #include <cmath>   
 
 #define PR(x) cout << #x << " = " << (x) << endl;
 
 using namespace std;
 
 TableGeneratorSettings* settings = TableGeneratorSettings::instance();
 
 int main(int argc, char *argv[]){
     //To explicitly delete all new histograms myself:
     TH1::AddDirectory(false);
     //
     //  Check arguments
     //
     char* runcard;
     if (argc != 4) {
         cout<<"Usage: lambdaGeneratorMain runcard startBin endBin"<<endl;
         return 2;
     }
     else{
         runcard = argv[1];
         settings->setStartBin(atoi(argv[2]));
         settings->setEndBin(atoi(argv[3]));
     }
     
     settings->readSettingsFromFile(runcard);
     settings->setA(1);
     settings->setModesToCalculate(1);
     settings->setNumberOfConfigurations(1);
     settings->consolidateSettings();
     
     int nBinQ2 =       settings->Q2bins();
     int nBinW2 =       settings->W2bins();
     int nBinT=         settings->tbins();
     double Q2min=      settings->Q2min();
     double Q2max=      settings->Q2max();
     double Wmin=       settings->Wmin();
     double Wmax=       settings->Wmax();
     double W2min=Wmin*Wmin;
     double W2max=Wmax*Wmax;
     double tmin=       settings->tmin();
     double tmax=       settings->tmax();
     unsigned int massA=settings->A();
     int vmPDG =        settings->vectorMesonId();
     DipoleModelType model = settings->dipoleModelType();
     DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
     int startingBin=settings->startBin();
     int endingBin=settings->endBin();
     int modes=settings->modesToCalculate();
     
     Table tableT;
     Table tableL;
     bool logQ2=true, logW2=false, logT=false, logC=false;
     
     //set filenames for the tables
     string rootfile=settings->rootfile();
     ostringstream filenameT, filenameL, filenameT2, filenameL2;
     filenameT.str("");
     filenameT<<rootfile<<"_bin"<<startingBin<<"_T.root";
     filenameL.str("");
     filenameL<<rootfile<<"_bin"<<startingBin<<"_L.root";
     (void) tableT.create(nBinQ2, Q2min, Q2max,
                          nBinW2, W2min, W2max,
                          nBinT,  tmin,  tmax,
                          logQ2, logW2, logT, logC,       // all bools
-                         lambda_A, transverse,
+                         lambda_real, transverse,
                          massA, vmPDG, model, pset,
                          filenameT.str().c_str());
     (void) tableL.create(nBinQ2, Q2min, Q2max,
                          nBinW2, W2min, W2max,
                          nBinT,  tmin,  tmax,
                          logQ2, logW2, logT, logC,      // all bools
-                         lambda_A, longitudinal,
+                         lambda_real, longitudinal,
                          massA, vmPDG, model, pset,
                          filenameL.str().c_str());
     //Create and initialize the amplitudes calculator:
     Amplitudes amps;
     
     //Generate the the nucleon configurations:
     amps.generateConfigurations();
     
     // Print out settings:
     cout<<"Tables will be generated for:"<<endl;
     cout<<"Nucleus mass A="<<massA<<endl;
     cout<<"Mode to calculate: "<<modes<<endl;
     cout<<"Vector Meson Id: "<<vmPDG<<endl;
     cout<<"Bins: "<<startingBin<<"-"<<endingBin<<endl;
     cout<<"Q2 range: ["<<Q2min<<", "<<Q2max<<"], "<<nBinQ2<<" bins."<<endl;
     cout<<"W2 range: ["<<W2min<<", "<<W2max<<"], "<<nBinW2<<" bins."<<endl;
     cout<<" t range: ["<<tmin<<", "<<tmax<<"], "<<nBinT<<" bins."<<endl;
     cout<<"Dipole Model: "<<model<<endl;
     cout<<endl;
     cout<<"Filling tables..."<<endl;
     time_t tableFillStart = time(0);
     for(int i=startingBin; i<endingBin; i++) {
         double Q2, W2, t;
         tableT.binCenter(i, Q2, W2, t);
 	double kinematicPoint[3]={t, Q2, W2};				
 	
         //calculate derivative of amplitude in point:
         double hplus, hminus;
         hplus=hminus=(W2max-W2min)/(4*1e4); //This value comes from tests of large and small W2
         hminus=min(hminus, W2-W2min);
         hplus=min(hplus, W2max-W2);
         hminus -= numeric_limits<float>::epsilon();
         hplus -= numeric_limits<float>::epsilon();
         //calculate contents and fill tables:
         double lambda[2]={0, 0};
 	kinematicPoint[2]= W2+hplus;
         amps.calculate(kinematicPoint);
         double ampPlus[2]={0, 0};
         ampPlus[0]=amps.amplitudeT();
         ampPlus[1]=amps.amplitudeL();
         
 	
 	kinematicPoint[2]=W2-hminus;
         amps.calculate(kinematicPoint);
         double ampMinus[2]={0, 0};
         ampMinus[0]=amps.amplitudeT();
         ampMinus[1]=amps.amplitudeL();
         for(int j=0; j<2; j++){
             PR(ampMinus[j]);
             PR(ampPlus[j]);
             // Don't calculate unless numerically viable:
             if(ampPlus[j]==0 || ampMinus[j]==0){
                 lambda[j]=0;
             }
             else{
                 //Calculate derivate d(logA/dW2)=log(A+/A-)/(h+-h-):
                 double derivate=log(abs(ampPlus[j]/ampMinus[j]))/(hplus+hminus);
                 //  Finally calculate lambda:
                 double jacobian = (W2-protonMass2+Q2);
                 lambda[j] = jacobian*derivate;
             }
         }
         tableT.fill(i, lambda[0]);
         tableL.fill(i, lambda[1]);
         if(i%static_cast<int>((endingBin-startingBin)/10+1)==0)
             cout<<"Processed "<<100*(i-startingBin)/(endingBin-startingBin)<<" %"<<endl;
         
         if(settings->verbose()){
             cout << "bin = "     << setw(8)  << left << i;
             cout << "Q2 = "      << setw(10) << left << fixed << setprecision(3) << Q2;
             cout << "W2 = "      << setw(11) << left << fixed << W2;
             cout << "t = "       << setw(14) << left << scientific << t;
             cout << "lambdaT = " << setw(12) << left << scientific << lambda[0]<<endl;
             cout << "bin = "     << setw(8)  << left << i;  
             cout << "Q2 = "      << setw(10) << left << fixed << setprecision(3) << Q2;  
             cout << "W2 = "      << setw(11) << left << fixed << W2;  
             cout << "t = "       << setw(14) << left << scientific << t;
             cout << "lambdaL = " << setw(12) << left << scientific << lambda[1]<<endl;
         }
     }
     cout<<"CPU Time/Entry: " << double(time(0)-tableFillStart)/endingBin<<" s/entry."<<endl;
     cout<<"Total time: "<<double(time(0)-tableFillStart)/60./60.<<" h"<<endl;
     
     tableT.write();
     tableL.write();
     return 0;
 }
Index: trunk/examples/tableGeneratorMain.cpp
===================================================================
--- trunk/examples/tableGeneratorMain.cpp	(revision 375)
+++ trunk/examples/tableGeneratorMain.cpp	(revision 376)
@@ -1,475 +1,475 @@
 //==============================================================================
 //  tableGeneratorMain.cpp
 //
 //  Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
 //
 //  This file is part of Sartre.
 //
 //  This program is free software: you can redistribute it and/or modify 
 //  it under the terms of the GNU General Public License as published by 
 //  the Free Software Foundation.   
 //  This program is distributed in the hope that it will be useful, 
 //  but without any warranty; without even the implied warranty of 
 //  merchantability or fitness for a particular purpose. See the 
 //  GNU General Public License for more details. 
 //  You should have received a copy of the GNU General Public License
 //  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 //
 //  Author: Thomas Ullrich
 //  Last update: 
 //  $Date$
 //  $Author$
 //==============================================================================
 //   
 //  Main program to create amplitude lookup tables.   
 //  [Developer only]
 //
 //  Usage:
 //     tableGeneratorMain runcard startBin endBin
 //
 //  Bins run from 0 to nbin-1.
 //  Loop fill all bins from startBIn to endBin (including endBin).
 //  If endBin exceeds the size of the table it set to the nbin-1.
 //
 //==============================================================================   
 #include <sstream>  
 #include <cstdlib>  
 #include <iostream>
 #include <limits>
 #include <cmath>
 #include <iomanip>
 #include "Amplitudes.h"  
 #include "TROOT.h"
 #include "TH1D.h"  
 #include "TFile.h"  
 #include "Nucleus.h"
 #include "Constants.h"
 #include "Table.h"  
 #include "TableGeneratorSettings.h"  
 #include "Enumerations.h"  
 #include "DglapEvolution.h"
 #include "Version.h"
 
 #define PR(x) cout << #x << " = " << (x) << endl;  
 
 using namespace std;  
 
 int main(int argc, char *argv[]) {
     
     //
     //  Print header
     //
     time_t theStartTime = time(0);
     string ctstr(ctime(&theStartTime));
     ctstr.erase(ctstr.size()-1, 1);
     cout << "/========================================================================\\" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  An event generator for exclusive diffractive vector meson production  |" << endl;
     cout << "|  in ep and eA collisions based on the dipole model.                    |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich                |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  This program is free software: you can redistribute it and/or modify  |" << endl;
     cout << "|  it under the terms of the GNU General Public License as published by  |" << endl;
     cout << "|  the Free Software Foundation, either version 3 of the License, or     |" << endl;
     cout << "|  any later version.                                                    |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Code compiled on " << setw(12) << left << __DATE__;
     cout << setw(41) << left << __TIME__ << right << '|' << endl;
     cout << "|  Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  --------------------------------------------------------------------  |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Sartre Table Generator (Experts only)                                 |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "\\========================================================================/" << endl;
 
     TH1::AddDirectory(false);  // to explicitly delete all new histograms by hand
     TableGeneratorSettings* settings = TableGeneratorSettings::instance();  
     
     //
     //  Check arguments  
     //  
     char* runcard;  
     if (argc != 4) {  
         cout << "Usage: tableGeneratorMain runcard startBin endBin" << endl;  
         return 2;
     }  
     else {
         runcard = argv[1];  
         settings->setStartBin(atoi(argv[2]));  
         settings->setEndBin(atoi(argv[3]));  
     }  
     
     cout << "Reading settings from runcard." << endl;
     settings->readSettingsFromFile(runcard);  
     settings->consolidateSettings();  
     
     int nBinQ2 = settings->Q2bins();
     int nBinW2 = settings->W2bins();
     int nBinT  = settings->tbins();
     double Q2min =  settings->Q2min();
     double Q2max =  settings->Q2max();
     double Wmin =   settings->Wmin();
     double Wmax =   settings->Wmax();
     double W2min =  Wmin*Wmin;
     double W2max =  Wmax*Wmax;
     double tmin =   settings->tmin();
     double tmax =   settings->tmax();
     unsigned int massA = settings->A();
     int vmPDG = settings->vectorMesonId();
     DipoleModelType model = settings->dipoleModelType();
     DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
     int startingBin = settings->startBin();
     int endingBin = settings->endBin();
     int modes = settings->modesToCalculate();
     unsigned char priority = settings->priority();
    
     settings->list();
     
     //
     //   Check if lambda tables can be calculated
     //
     bool createLambdaTables = true;
     if (massA == 1 && modes == 1 && settings->numberOfConfigurations() == 1) {
         cout << "\nLambda tables will be generated." << endl;
     }
     else {
         cout << "\nLambda tables will not be generated. Requires A = 1, mode = 1, and 1 configuration only." << endl;
         createLambdaTables = false;
     }
     
     
     //
     //  Check bins
     //
     //  Table's bin indices run from 0 ... nbins-1
     //
     int maxbins = nBinQ2*nBinW2*nBinT;
     if (endingBin >= maxbins) {
         cout << "Warning, given end bin (" << endingBin << ") exceeds the size of the table." << endl;
         cout << "         set to maximum value (" << maxbins-1 << ") now." << endl;
         endingBin = maxbins-1;
     }
     
     //
     //   Define all tables. Depending on tableset type some
     //   will not be written but we define them all anyway.
     //
     Table tableT;  
     Table tableL;  
     Table tableT2;  
     Table tableL2;  
     Table tableVarT;
     Table tableVarL;
     Table tableLambdaT;
     Table tableLambdaL;
     bool logQ2=true, logW2=true, logT=false, logC=true;
     
     //
     //  Set filenames for the tables.
     //  Be as desciptive as possible. Helps when mass producing tables.
     //
     //  We create all tables and decide later what
     //  gets written and what not.
     //
     string rootfile=settings->rootfile();
     rootfile += "_" + settings->dipoleModelName();
     rootfile += "_" + settings->dipoleModelParameterSetName();
     
     ostringstream filenameT, filenameL, filenameT2,
                   filenameL2, filenameVarT, filenameVarL,
                   filenameLambdaT, filenameLambdaL;
     filenameT.str("");  
     filenameT << rootfile << "_" << settings->vectorMesonId() << "_bin"
               << startingBin << "-" << endingBin <<"_T.root";
     filenameL.str("");  
     filenameL << rootfile << "_" << settings->vectorMesonId() << "_bin"
               << startingBin << "-" << endingBin << "_L.root";
     filenameT2.str("");  
     filenameT2 << rootfile << "_" << settings->vectorMesonId() << "_bin"
                << startingBin << "-" << endingBin << "_T2.root";
     filenameL2.str("");  
     filenameL2 << rootfile << "_" << settings->vectorMesonId() << "_bin"
                << startingBin << "-" << endingBin << "_L2.root";
     filenameVarT.str("");
     filenameVarT << rootfile << "_" << settings->vectorMesonId() << "_bin"
                  << startingBin << "-" << endingBin << "_VarT.root";
     filenameVarL.str("");
     filenameVarL << rootfile << "_" << settings->vectorMesonId() << "_bin"
                  << startingBin << "-" << endingBin << "_VarL.root";
     filenameLambdaT.str("");
     filenameLambdaT << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin << "_LambdaT.root";
     filenameLambdaL.str("");
     filenameLambdaL << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin << "_LambdaL.root";
 
     (void) tableT.create(nBinQ2, Q2min, Q2max,
                          nBinW2, W2min, W2max,  
                          nBinT,  tmin,  tmax,  
                          logQ2, logW2, logT, logC,       // all bools  
                          mean_A, transverse,  
                          massA, vmPDG, model, pset,
                          filenameT.str().c_str(),
 			             priority);
     (void) tableL.create(nBinQ2, Q2min, Q2max,  
                          nBinW2, W2min, W2max,  
                          nBinT,  tmin,  tmax,  
                          logQ2, logW2, logT, logC,      // all bools  
                          mean_A, longitudinal,  
                          massA, vmPDG, model, pset,
                          filenameL.str().c_str(),
 			             priority);
     (void) tableT2.create(nBinQ2, Q2min, Q2max,  
                           nBinW2, W2min, W2max,  
                           nBinT,  tmin,  tmax,  
                           logQ2, logW2, logT, logC,      // all bools  
                           mean_A2, transverse,  
                           massA, vmPDG, model, pset,
                           filenameT2.str().c_str(),
 			              priority);
     (void) tableL2.create(nBinQ2, Q2min, Q2max,  
                           nBinW2, W2min, W2max,  
                           nBinT,  tmin,  tmax,  
                           logQ2, logW2, logT, logC,      // all bools  
                           mean_A2, longitudinal,  
                           massA, vmPDG, model, pset,  
                           filenameL2.str().c_str(),
 			              priority);
     bool contentVar = modes == 1 ? false : logC;  // if mode = 1 they will be 0 so lin is better
     (void) tableVarT.create(nBinQ2, Q2min, Q2max,
                             nBinW2, W2min, W2max,
                             nBinT,  tmin,  tmax,
                             logQ2, logW2, logT, contentVar,      // all bools
                             variance_A, transverse,
                             massA, vmPDG, model, pset,
                             filenameVarT.str().c_str(),
                             priority);
     (void) tableVarL.create(nBinQ2, Q2min, Q2max,
                             nBinW2, W2min, W2max,
                             nBinT,  tmin,  tmax,
                             logQ2, logW2, logT, contentVar,      // all bools
                             variance_A, longitudinal,
                             massA, vmPDG, model, pset,
                             filenameVarL.str().c_str(),
                             priority);
     double contentLambda = false; // content lin
     (void) tableLambdaT.create(nBinQ2, Q2min, Q2max,
                             nBinW2, W2min, W2max,
                             nBinT,  tmin,  tmax,
                             logQ2, logW2, logT, contentLambda,      // all bools
-                            lambda_A, transverse,
+                            lambda_real, transverse,
                             massA, vmPDG, model, pset,
                             filenameLambdaT.str().c_str(),
                             priority);
     (void) tableLambdaL.create(nBinQ2, Q2min, Q2max,
                             nBinW2, W2min, W2max,
                             nBinT,  tmin,  tmax,
                             logQ2, logW2, logT, contentLambda,      // all bools
-                            lambda_A, longitudinal,
+                            lambda_real, longitudinal,
                             massA, vmPDG, model, pset,
                             filenameLambdaL.str().c_str(),
                             priority);
 
     cout << "\nAll tables created:" << endl;
     tableT.list();
     tableL.list();
     tableT2.list();
     tableL2.list();
     tableVarT.list();
     tableVarL.list();
     tableLambdaT.list();
     tableLambdaL.list();
     cout << "\nTables have " << maxbins << " bins each.\n" << endl;
     
     if (settings->useBackupFile()) {
         int ibin = settings->startingBinFromBackup();
         tableT.setAutobackup("tableT", ibin);
         tableL.setAutobackup("tableL", ibin);
         tableT2.setAutobackup("tableT2", ibin);
         tableL2.setAutobackup("tableL2", ibin);
         tableVarT.setAutobackup("tableVarT", ibin);
         tableVarL.setAutobackup("tableVarL", ibin);
         tableLambdaT.setAutobackup("tableLambdaT", ibin);
         tableLambdaL.setAutobackup("tableLambdaL", ibin);
         cout << "Automatic backup of tables is enabled." << endl;
     }
     else
         cout << "Automatic backup of tables is off.\n" << endl;
 
     //
     //   DGLAP Evolution can be speed up by using lookup tables
     //
     DglapEvolution &dglap = DglapEvolution::instance();
     dglap.generateLookupTable(1000, 1000);
     dglap.useLookupTable(true);
 
     //
     //   Create and initialize the amplitudes calculator
     //
     Amplitudes amps;
  
     //
     //   Generate the the nucleon configurations
     //
     amps.generateConfigurations();
     
     //
     //   Print out settings
     //
     cout << endl;
     cout << "Tables will be generated for:" << endl;  
     cout << "\tNucleus mass A="<<massA<<endl;
     cout << "\tModes to calculate: " << modes << endl;
     cout << "\tVector Meson Id: " << vmPDG << endl;
     cout << "\tBins: " << startingBin << "-" << endingBin << endl;
     cout << "\tQ2 range: [" << Q2min << ", " << Q2max << "], " << nBinQ2 << " bins." << endl;
     cout << "\tW2 range: [" << W2min << ", " << W2max << "], " << nBinW2 << " bins." << endl;
     cout << "\t t range: [" << tmin << ", " << tmax << "], " << nBinT << " bins." << endl;
     cout << "\tDipole model: " << settings->dipoleModelName() << endl;
     cout << "\tDipole model parameter set: " << settings->dipoleModelParameterSetName() << endl;
     cout << "\tTable set mode: " << settings->tableSetTypeName() << endl;
     cout << endl;
     
     //
     //   Calculate contents and fill tables
     //
     //   Note that we fill all tables. What is written
     //   at the end is another story.
     //
     
     cout << "Start filling tables" << endl;
     time_t tableFillStart = time(0);
     
     int nShow = (endingBin - startingBin)/100;
     
     for (int i=startingBin; i<=endingBin; i++) {
         
         if (nShow && (i%nShow == 0 || i == startingBin || i == endingBin))
             cout << "processing bin " << i << endl;
         double Q2, W2, t;
         tableT.binCenter(i, Q2, W2, t);
 	double kinematicPoint[3]={t, Q2, W2};
         amps.calculate(kinematicPoint);
         
         double aT = 0;
         double aL = 0;
         double aT2 = 0;
         double aL2 = 0;
         double aVarT = 0;
         double aVarL = 0;
         
         //
         //  For mode = 1 only coherent is
         //  calculated and incoherent is therefore
         //  assumed to be 0.
         //  Hence total == coherent
         //
         //  Recall:
         //  coherent = A
         //  total = A2
         //  incoherent = variance = A2 - A*A
         //
         aL = amps.amplitudeL();
         aT = amps.amplitudeT();
         tableT.fill(i, aT);
         tableL.fill(i, aL);
         
         if (modes != 1) {
             aT2 = amps.amplitudeT2();
             aL2 = amps.amplitudeL2();
         }
         else {
             aT2 = aT*aT;
             aL2 = aL*aL;
         }
         tableT2.fill(i, aT2);
         tableL2.fill(i, aL2);
         
         if (modes != 1) {
             aVarT = aT2-aT*aT;
             aVarL = aL2-aL*aL;
         }
         tableVarT.fill(i, aVarT);
         tableVarL.fill(i, aVarL);
         
         //
         //  Calculate lambda and fill table
         //
         if (createLambdaTables) {
             double hplus, hminus;
             hplus=hminus=(W2max-W2min)/(4*1e4); //This value comes from tests of large and small W2
             hminus=min(hminus, W2-W2min);
             hplus=min(hplus, W2max-W2);
             hminus -= numeric_limits<float>::epsilon();
             hplus -= numeric_limits<float>::epsilon();
             double lambda[2]={0, 0};
 	    kinematicPoint[2]=W2+hplus;
             amps.calculate(kinematicPoint);
             double ampPlus[2]={0, 0};
             ampPlus[0]=amps.amplitudeT();
             ampPlus[1]=amps.amplitudeL();
 
 	    kinematicPoint[2]=W2-hminus;
             amps.calculate(kinematicPoint);
             double ampMinus[2]={0, 0};
             ampMinus[0]=amps.amplitudeT();
             ampMinus[1]=amps.amplitudeL();
             for (int j=0; j<2; j++) {
                 // Don't calculate unless numerically viable:
                 if (ampPlus[j]==0 || ampMinus[j]==0) {
                      lambda[j]=0;
                 }
                 else {
                     //Calculate derivate d(logA/dW2)=log(A+/A-)/(h+-h-):
                     double derivate=log(abs(ampPlus[j]/ampMinus[j]))/(hplus+hminus);
                     //  Finally calculate lambda:
                     double jacobian = (W2-protonMass2+Q2);
                     lambda[j] = jacobian*derivate;
                 }
             }
             tableLambdaT.fill(i, lambda[0]);
             tableLambdaL.fill(i, lambda[1]);
         }
     }
     time_t tableFillEnd = time(0);
     
     //
     //   Report CPU time/cell
     //
     cout << endl;
     cout << "CPU time/bin: "
          <<  static_cast<double>(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
          << " s" << endl;
     cout << "Total time: " << static_cast<double>(tableFillEnd-tableFillStart)/60./60. << " h" << endl << endl;
     
     //
     //  We write out all tables.
     //
     //  Whoever runs the production can then decide
     //  later what to keep and what to delete.
     //  If the desired run mode is total_and_coherent or
     //  coherent_and_incoherent, if this is just to improve
     //  a coherent table in some phase space, or if one wants
     //  to maintain redundancy, all these factor might affect
     //  your choice.
     //
     tableT.write();
     tableL.write();
     tableT2.write();
     tableL2.write();
     tableVarT.write();
     tableVarL.write();
     if (createLambdaTables) {
         tableLambdaT.write();
         tableLambdaL.write();
     }
     cout << "All tables written" << endl;
     
     cout << "All done. Bye." << endl;
     
     return 0;
 }  
Index: trunk/examples/tableGeneratorUPCMain.cpp
===================================================================
--- trunk/examples/tableGeneratorUPCMain.cpp	(revision 375)
+++ trunk/examples/tableGeneratorUPCMain.cpp	(revision 376)
@@ -1,417 +1,417 @@
 //==============================================================================
 //  tableGeneratorMain.cpp
 //
 //  Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
 //
 //  This file is part of Sartre.
 //
 //  This program is free software: you can redistribute it and/or modify 
 //  it under the terms of the GNU General Public License as published by 
 //  the Free Software Foundation.   
 //  This program is distributed in the hope that it will be useful, 
 //  but without any warranty; without even the implied warranty of 
 //  merchantability or fitness for a particular purpose. See the 
 //  GNU General Public License for more details. 
 //  You should have received a copy of the GNU General Public License
 //  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 //
 //  Author: Thomas Ullrich
 //  Last update: 
 //  $Date$
 //  $Author$
 //==============================================================================
 //   
 //  Main program to create amplitude lookup tables.   
 //  [Developer only]
 //
 //  Usage:
 //     tableGeneratorMain runcard startBin endBin
 //
 //  Bins run from 0 to nbin-1.
 //  Loop fill all bins from startBIn to endBin (including endBin).
 //  If endBin exceeds the size of the table it set to the nbin-1.
 //
 //==============================================================================   
 #include <sstream>  
 #include <cstdlib>  
 #include <iostream>
 #include <limits>
 #include <cmath>
 #include <iomanip>
 #include "Amplitudes.h"  
 #include "TROOT.h"
 #include "TH1D.h"  
 #include "TFile.h"  
 #include "Nucleus.h"
 #include "Constants.h"
 #include "Table.h"  
 #include "TableGeneratorSettings.h"  
 #include "Enumerations.h"  
 #include "DglapEvolution.h"
 #include "Version.h"
 
 #define PR(x) cout << #x << " = " << (x) << endl;  
 
 using namespace std;  
 
 int main(int argc, char *argv[]) {
     
     //
     //  Print header
     //
     time_t theStartTime = time(0);
     string ctstr(ctime(&theStartTime));
     ctstr.erase(ctstr.size()-1, 1);
     cout << "/========================================================================\\" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  An event generator for exclusive diffractive vector meson production  |" << endl;
     cout << "|  in ep and eA collisions based on the dipole model.                    |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich                |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  This program is free software: you can redistribute it and/or modify  |" << endl;
     cout << "|  it under the terms of the GNU General Public License as published by  |" << endl;
     cout << "|  the Free Software Foundation, either version 3 of the License, or     |" << endl;
     cout << "|  any later version.                                                    |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Code compiled on " << setw(12) << left << __DATE__;
     cout << setw(41) << left << __TIME__ << right << '|' << endl;
     cout << "|  Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  --------------------------------------------------------------------  |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "|  Sartre Table Generator (Experts only)                                 |" << endl;
     cout << "|                                                                        |" << endl;
     cout << "\\========================================================================/" << endl;
     
     TH1::AddDirectory(false);  // to explicitly delete all new histograms by hand
     TableGeneratorSettings* settings = TableGeneratorSettings::instance();
     
     //
     //  Check arguments
     //
     char* runcard;
     if (argc != 4) {
         cout << "Usage: tableGeneratorMain runcard startBin endBin" << endl;
         return 2;
     }
     else {
         runcard = argv[1];
         settings->setStartBin(atoi(argv[2]));
         settings->setEndBin(atoi(argv[3]));
     }
     
     cout << "Reading settings from runcard." << endl;
     settings->readSettingsFromFile(runcard);
     settings->consolidateSettings();
     
     int nBinT  = settings->tbins();
     int nBinX  = settings->xbins();
     double tmin =   settings->tmin();
     double tmax =   settings->tmax();
     double xmin =   settings->xmin();
     double xmax =   settings->xmax();
     unsigned int massA = settings->A();
     int vmPDG = settings->vectorMesonId();
     DipoleModelType model = settings->dipoleModelType();
     DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
     int startingBin = settings->startBin();
     int endingBin = settings->endBin();
     int modes = settings->modesToCalculate();
     unsigned char priority = settings->priority();
     
     settings->list();
     
     //
     //   Check if lambda tables can be calculated
     //
     bool createLambdaTables = true;
     if (massA == 1 && modes != 1 && settings->numberOfConfigurations() == 1) {
         cout << "\nLambda tables will be generated." << endl;
     }
     else {
         cout << "\nLambda tables will not be generated. Requires A = 1, mode != 1, and 1 configuration only." << endl;
         createLambdaTables = false;
     }
     
     
     //
     //  Check bins
     //
     //  Table's bin indices run from 0 ... nbins-1
     //
     int maxbins = nBinX*nBinT;
     if (endingBin >= maxbins) {
         cout << "Warning, given end bin (" << endingBin << ") exceeds the size of the table." << endl;
         cout << "         set to maximum value (" << maxbins-1 << ") now." << endl;
         endingBin = maxbins-1;
     }
     
     //
     //   Define all tables. Depending on tableset type some
     //   will not be written but we define them all anyway.
     //
     Table tableT;
     Table tableT2;
     Table tableVarT;
     Table tableLambdaRealT;
     Table tableLambdaSkewT;
     bool  logX=true, logT=false, logC=true;
     
     //
     //  Set filenames for the tables.
     //  Be as desciptive as possible. Helps when mass producing tables.
     //
     //  We create all tables and decide later what
     //  gets written and what not.
     //
     string rootfile=settings->rootfile();
     rootfile += "_"; rootfile += settings->dipoleModelName();
     rootfile += "_"; rootfile += settings->dipoleModelParameterSetName();
     rootfile += "_A"; rootfile += to_string(settings->A());
     rootfile += "_UPC";
     
     ostringstream filenameT, filenameT2, filenameVarT, filenameLambdaTReal, filenameLambdaTSkew;
     filenameT.str("");
     filenameT << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin <<"_T.root";
     filenameT2.str("");
     filenameT2 << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin << "_T2.root";
     filenameVarT.str("");
     filenameVarT << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin << "_VarT.root";
     filenameLambdaTReal.str("");
     filenameLambdaTReal << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin << "_LambdaTReal.root";
     filenameLambdaTSkew.str("");
     filenameLambdaTSkew << rootfile << "_" << settings->vectorMesonId() << "_bin"
     << startingBin << "-" << endingBin << "_LambdaTSkew.root";
     
     (void) tableT.create(nBinX, xmin, xmax,
                          nBinT,  tmin,  tmax,
                          logX, logT, logC,       // all bools
                          mean_A, transverse,
                          massA, vmPDG, model, pset,
                          filenameT.str().c_str(),
                          priority);
     (void) tableT2.create(nBinX, xmin, xmax,
                           nBinT,  tmin,  tmax,
                           logX, logT, logC,      // all bools
                           mean_A2, transverse,
                           massA, vmPDG, model, pset,
                           filenameT2.str().c_str(),
                           priority);
     bool contentVar = modes == 1 ? false : logC;  // if mode = 1 they will be 0 so lin is better
     (void) tableVarT.create(nBinX, xmin, xmax,
                             nBinT,  tmin,  tmax,
                             logX, logT, contentVar,      // all bools
                             variance_A, transverse,
                             massA, vmPDG, model, pset,
                             filenameVarT.str().c_str(),
                             priority);
     double contentLambda = false; // content lin
     if(createLambdaTables){
         (void) tableLambdaRealT.create(nBinX, xmin, xmax,
 				       nBinT,  tmin,  tmax,
 				       logX, logT, contentLambda,      // all bools
-				       lambda_A, transverse,
+				       lambda_real, transverse,
 				       massA, vmPDG, model, pset,
 				       filenameLambdaTReal.str().c_str(),
 				       priority);
 	(void) tableLambdaSkewT.create(nBinX, xmin, xmax,
 				       nBinT,  tmin,  tmax,
 				       logX, logT, contentLambda,      // all bools
 				       lambda_skew, transverse,
 				       massA, vmPDG, model, pset,
 				       filenameLambdaTSkew.str().c_str(),
 				       priority);
     }    
     cout << "\nAll tables created:" << endl;
     tableT.list();
     tableT2.list();
     tableVarT.list();
     if(createLambdaTables){
       tableLambdaRealT.list();
       tableLambdaSkewT.list();
     }
     cout << "\nTables have " << maxbins << " bins each.\n" << endl;
     
     if (settings->useBackupFile()) {
         int ibin = settings->startingBinFromBackup();
         tableT.setAutobackup("tableT", ibin);
         tableT2.setAutobackup("tableT2", ibin);
         tableVarT.setAutobackup("tableVarT", ibin);
         tableLambdaRealT.setAutobackup("tableLambdaRealT", ibin);
         tableLambdaSkewT.setAutobackup("tableLambdaSkewT", ibin);
         cout << "Automatic backup of tables is enabled." << endl;
     }
     else
         cout << "Automatic backup of tables is off.\n" << endl;
     
     //
     //   DGLAP Evolution can be speed up by using lookup tables
     //
     DglapEvolution &dglap = DglapEvolution::instance();
     //    dglap.generateLookupTable(100, 100);
     dglap.generateLookupTable(1000, 1000);
     dglap.useLookupTable(true);
     
     //
     //   Create and initialize the amplitudes calculator
     //
     Amplitudes amps;
     
     //
     //   Generate the the nucleon configurations
     //
     amps.generateConfigurations();
     
     //
     //   Print out settings
     //
     cout << endl;
     cout << "Tables will be generated for:" << endl;
     cout << "\tNucleus mass A="<<massA<<endl;
     cout << "\tModes to calculate: " << modes << endl;
     cout << "\tVector Meson Id: " << vmPDG << endl;
     cout << "\tBins: " << startingBin << "-" << endingBin << endl;
     cout << "\txpom range: [" << xmin << ", " << xmax << "], " << nBinX << " bins." << endl;
     cout << "\t t range: [" << tmin << ", " << tmax << "], " << nBinT << " bins." << endl;
     cout << "\tDipole model: " << settings->dipoleModelName() << endl;
     cout << "\tDipole model parameter set: " << settings->dipoleModelParameterSetName() << endl;
     cout << "\tTable set mode: " << settings->tableSetTypeName() << endl;
     cout << endl;
     
     //
     //   Calculate contents and fill tables
     //
     //   Note that we fill all tables. What is written
     //   at the end is another story.
     //
     
     cout << "Start filling tables" << endl;
     time_t tableFillStart = time(0);
     
     int nShow = (endingBin - startingBin)/100;
     if(nShow==0) nShow=1;
     
     for (int i=startingBin; i<=endingBin; i++) {
         if (i%nShow == 0 || i == startingBin || i == endingBin)
             cout << "processing bin " << i << endl;
         double xpom, t;
         
         tableT.binCenter(i, xpom, t);
         double kinematicPoint[2]={t, xpom};
         amps.calculate(kinematicPoint);
         
         double aT = 0;
         double aT2 = 0;
         double aVarT = 0;
         
         //
         //  For mode = 1 only coherent is
         //  calculated and incoherent is therefore
         //  assumed to be 0.
         //  Hence total == coherent
         //
         //  coherent = A
         //  total = A2
         //  incoherent = variance = A2 - A*A
         //
         aT = amps.amplitudeT();
         tableT.fill(i, aT);
         
         if (modes != 1) {
             aT2 = amps.amplitudeT2();
         }
         else {
             aT2 = aT*aT;
         }
         tableT2.fill(i, aT2);
         
         aVarT = aT2-aT*aT;
         
         tableVarT.fill(i, aVarT);
 
         //
         //  Calculate lambda and fill table lambda=dlog(A)/dlog(1/x)
         //
         if (createLambdaTables) {
             double hplus, hminus;
             hplus=hminus=xpom*1e-4; //This value needs to be tested
             double lambdaReal=0;
 	    double lambdaSkew=0;
             kinematicPoint[1]=xpom+hplus;
             amps.calculate(kinematicPoint);
             double ampPlus=amps.amplitudeT();
 	    double ampSkewPlus=amps.amplitudeTForSkewednessCorrection();
 	    
             kinematicPoint[1]=xpom-hminus;
             amps.calculate(kinematicPoint);
             double ampMinus=amps.amplitudeT();
 	    double ampSkewMinus=amps.amplitudeTForSkewednessCorrection();
 
             // Don't calculate unless numerically viable:
             if (ampPlus==0 || ampMinus==0) {
                 lambdaReal=0;
             }
             else {
                 //Calculate derivate d(logA/dxpom)=log(A+/A-)/(h+-h-):
                 double derivate=log(abs(ampPlus/ampMinus))/(hplus+hminus);
                 //  Finally calculate lambdas:
                 double jacobian = -xpom;
                 lambdaReal = jacobian*derivate;
             }
             // Don't calculate unless numerically viable:
             if (ampSkewPlus==0 || ampSkewMinus==0) {
                 lambdaSkew=0;
             }
             else {
                 //Calculate derivate d(logA/dxpom)=log(A+/A-)/(h+-h-):
                 double derivate=log(abs(ampSkewPlus/ampSkewMinus))/(hplus+hminus);
                 //  Finally calculate lambdas:
                 double jacobian = -xpom;
                 lambdaSkew = jacobian*derivate;
             }
             tableLambdaRealT.fill(i, lambdaReal);
 	    tableLambdaSkewT.fill(i, lambdaSkew);
         }
     }
     time_t tableFillEnd = time(0);
     
     //
     //   Report CPU time/cell
     //
     cout << endl;
     cout << "CPU time/bin: "
          <<  static_cast<double>(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
          << " s" << endl;
     cout << "Total time: " << static_cast<double>(tableFillEnd-tableFillStart)/60./60. << " h" << endl << endl;
     
     //
     //  We write out all tables.
     //
     //  Whoever runs the production can then decide
     //  later what to keep and what to delete.
     //  If the desired run mode is total_and_coherent or
     //  coherent_and_incoherent, if this is just to improve
     //  a coherent table in some phase space, or if one wants
     //  to maintain redundancy, all these factor might affect
     //  your choice.
     //
     tableT.write();
     tableT2.write();
     tableVarT.write();
     if (createLambdaTables) {
         tableLambdaRealT.write();
         tableLambdaSkewT.write();
     }
     cout << "All tables written" << endl;
     
     cout << "All done. Bye." << endl;
     
     return 0;
 }