Index: trunk/examples/CMakeLists.txt
===================================================================
--- trunk/examples/CMakeLists.txt (revision 370)
+++ trunk/examples/CMakeLists.txt (revision 371)
@@ -1,162 +1,162 @@
#===============================================================================
# CMakeLists.txt (examples)
#
-# Copyright (C) 2010-2016 Tobias Toll and Thomas Ullrich
+# 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 .
#
# Author: Thomas Ullrich
# Last update:
# $Date$
# $Author$
#===============================================================================
cmake_minimum_required (VERSION 3.1)
project(SARTRE_EXAMPLES)
#
# Compiler options
#
set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -pedantic -Wno-long-long")
set(CMAKE_CXX_FLAGS_DEBUG "-g")
set(CMAKE_CXX_FLAGS_RELEASE "-O")
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
#
# Set the location of the installed sartre package
#
set(SARTRE_DIR "$ENV{SARTRE_DIR}")
if (IS_DIRECTORY "${SARTRE_DIR}" AND
IS_DIRECTORY "${SARTRE_DIR}/lib" AND
IS_DIRECTORY "${SARTRE_DIR}/include" AND
IS_DIRECTORY "${SARTRE_DIR}/gemini" AND
IS_DIRECTORY "${SARTRE_DIR}/tables")
message (STATUS "Using Sartre installed in ${SARTRE_DIR}")
else()
message(FATAL_ERROR "Environment variable SARTRE_DIR does not point to a valid Sartre installation.")
endif()
message(STATUS "System = ${CMAKE_SYSTEM}")
#
# User can select if Sartre is to be compiled
# in multithreaded mode or not.
# Default is off.
# cmake -DMULTITHREADED=ON .
#
option(MULTITHREADED "Switch on/off multithreaded mode" OFF)
if (MULTITHREADED)
add_definitions(-DSARTRE_IN_MULTITHREADED_MODE)
message (STATUS "MULTITHREADED mode is ON")
endif (MULTITHREADED)
#
# See if optional libraries (Glut & OpenGL) are present.
# If not skip the examples that require them.
#
# GLUT
#find_package(GLUT)
#if (GLUT_FOUND)
# message("blah")
# include_directories(${GLUT_INCLUDE_DIR})
# set(LIBS ${LIBS} ${GLUT_LIBRARIES})
#endif()
# OpenGL
find_package(OpenGL)
if (OPENGL_FOUND)
include_directories(${OPENGL_INCLUDE_DIR})
set(LIBS ${LIBS} ${OPENGL_LIBRARIES})
endif()
#
# Built all programs in example directory
#
add_executable(createBSatBDependenceTable createBSatBDependenceTable.cpp)
add_executable(heraCompare heraCompare.cpp)
add_executable(sartreMain sartreMain.cpp)
add_executable(sartreTest sartreTest.cpp)
add_executable(tableGeneratorMain tableGeneratorMain.cpp)
add_executable(tableGeneratorUPCMain tableGeneratorUPCMain.cpp)
add_executable(lambdaGeneratorMain lambdaGeneratorMain.cpp)
add_executable(convergenceTest convergenceTest.cpp)
#if (OPENGL_FOUND AND GLUT_FOUND)
# add_executable(nucleusAnimation nucleusAnimation.cpp)
#endif()
#
# Make sure the files to search for external libraries are available
#
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/modules")
#
# Compiler options
#
set(CMAKE_C_FLAGS_DEBUG " -g -W")
set(CMAKE_CXX_FLAGS_DEBUG " -g -W -Wall -Wextra -pedantic -Wno-long-long -std=c++11")
set(CMAKE_C_FLAGS_RELEASE " -O -W")
set(CMAKE_CXX_FLAGS_RELEASE " -O -W -Wall -Wextra -pedantic -Wno-long-long -std=c++11")
set(CMAKE_C_FLAGS " -O -W")
set(CMAKE_CXX_FLAGS " -O -std=c++11 -W -Wall -Wextra -pedantic -Wno-long-long")
#
# Include files locations
#
include_directories("${SARTRE_DIR}/include")
#
# Find external required packages
# (see also FindGSL.cmake and FindROOT.cmke in cmake/modules)
#
# GSL
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIR})
set(LIBS ${LIBS} ${GSL_LIBRARIES})
# ROOT
find_package(ROOT REQUIRED)
include_directories(${ROOT_INCLUDE_DIR})
set(LIBS ${LIBS} ${ROOT_LIBRARIES} "-lEG -lMathMore -lUnuran" )
#BOOST
if (MULTITHREADED)
set(Boost_USE_MULTITHREADED ON)
find_package(Boost 1.39 COMPONENTS thread REQUIRED)
if(Boost_FOUND)
include_directories(${Boost_INCLUDE_DIR})
set(LIBS ${LIBS} ${Boost_LIBRARIES} -lboost_system)
endif(Boost_FOUND)
endif (MULTITHREADED)
#
# Libraries to link with
#
#link_directories(${SARTRE_DIR}/lib)
file (GLOB ALLSARTRE ${SARTRE_DIR}/lib/*)
set(LIBS ${ALLSARTRE} ${ALLSARTRE} ${LIBS} ${ALLSARTRE})
target_link_libraries(createBSatBDependenceTable ${LIBS})
target_link_libraries(heraCompare ${LIBS})
target_link_libraries(sartreMain ${LIBS})
target_link_libraries(sartreTest ${LIBS})
target_link_libraries(tableGeneratorMain ${LIBS})
target_link_libraries(tableGeneratorUPCMain ${LIBS})
target_link_libraries(lambdaGeneratorMain ${LIBS})
target_link_libraries(convergenceTest ${LIBS})
#if (OPENGL_FOUND AND GLUT_FOUND)
# target_link_libraries(nucleusAnimation ${LIBS})
#endif()
Index: trunk/examples/lambdaGeneratorMain.cpp
===================================================================
--- trunk/examples/lambdaGeneratorMain.cpp (revision 370)
+++ trunk/examples/lambdaGeneratorMain.cpp (revision 371)
@@ -1,192 +1,192 @@
//==============================================================================
// lambdaGeneratorMain.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Thomas Ullrich
// Last update: Sat Aug 13 02:18:37 2011
//==============================================================================
//
// Main program to create lambda lookup tables.
// [Developer only]
//==============================================================================
#include
#include
#include
#include
#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
#include
#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"<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<::epsilon();
hplus -= numeric_limits::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((endingBin-startingBin)/10+1)==0)
cout<<"Processed "<<100*(i-startingBin)/(endingBin-startingBin)<<" %"<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]<.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Developer test program. Checks individual components.
// Should not be in the distribution.
//
//==============================================================================
#include "EventGeneratorSettings.h"
#include "Table.h"
#include
#include "TParticlePDG.h"
#include "TRandom2.h"
#include "Event.h"
#include "ExclusiveFinalStateGenerator.h"
#include "Kinematics.h"
#include "FrangibleNucleus.h"
#include "Enumerations.h"
#include "TableCollection.h"
#include "Sartre.h"
#include
#include
#include "TFile.h"
#include "TH3D.h"
#include "TH2D.h"
#include "Nucleus.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
//#define TEST_SETTINGS
//#define TEST_WRITETABLE
//#define TEST_READTABLE
//#define TEST_BREAKUP
//#define TEST_FINALSTATE
//#define TEST_TABLECOLLECTION
//#define TEST_GENERATOR
//#define TEST_INTERPOLATE
//#define TEST_SLOPE
//#define TEST_CROSSSECTION
-//#define TEST_GLAUBER
-#define TEST_UPC_TABLE
+#define TEST_GLAUBER
+//#define TEST_UPC_TABLE
#if defined(TEST_INTERPOLATE)
struct InterTestPoint {
double t, Q2, W2, value;
};
#endif
#if defined(TEST_SLOPE)
struct InterTestSlope {
double t, Q2, W2, slope;
};
#endif
#if defined(TEST_UPC_TABLE)
double testfun(double x, double t)
{
return fabs(log(x)*exp(t)*(exp(t)+x)+sin(t*300));
}
#endif
int main()
{
#if defined(TEST_SETTINGS)
//
// Test Settings
//
EventGeneratorSettings* settings = EventGeneratorSettings::instance();
settings->readSettingsFromFile("sartreRuncard.txt");
settings->setVerbose(true);
settings->setVerboseLevel(3);
settings->list();
//
// Test PDG lookup
//
TParticlePDG *part = settings->lookupPDG(113);
PR(part->Mass());
PR(part->ParticleClass());
PR(part->GetName());
#endif
//
// Test reading tables
//
#if defined(TEST_READTABLE)
Table table1, table2;
table1.read("table1.root");
table2.read("table2.root");
table3.read("table3.root");
table1.list(cout, false);
table2.list(cout, false);
table3.list(cout, false);
#endif
#if defined (TEST_UPC_TABLE)
cout << "Creating UPC tables" << endl;
Table tableA;
Table tableA2;
int n1 = tableA.create(50, 1e-5, 0.2, // x
100, -0.051, 0, // t
true, false, true, // bool logx, bool logt, bool logContent,
mean_A, transverse, 197, 411, bSat, HMPZ, "tableA.root");
int n2 = tableA2.create(1000, 1e-5, 0.2,
1000, -0.051, 0,
true, false, true,
mean_A2, transverse, 197, 411, bSat, HMPZ, "tableA2.root");
cout << "Filling UPC tables" << endl;
double xpom, t;
for (int i=0; ieBeam();
event.particles[1].p = settings->hBeam();
generator.generate(443, 0, -event.t, event.y, event.Q2, false, &event);
cout << endl;
event.list();
cout << endl;
#endif
#if defined(TEST_BREAKUP)
FrangibleNucleus kern(208, true);
kern.normalizationOfT();
TLorentzVector someVec(0.5, 0.3, 99., sqrt(0.5*0.5+0.3*0.3+99.*99.)+0.01);
kern.breakup(someVec);
kern.listBreakupProducts();
kern.breakup(someVec);
kern.listBreakupProducts();
kern.breakup(someVec);
kern.listBreakupProducts();
#endif
#if defined(TEST_TABLECOLLECTION)
TableCollection coll;
coll.init(197, bSat, 443);
PR(coll.minQ2());
PR(coll.maxQ2());
PR(coll.minW2());
PR(coll.maxW2());
PR(coll.minW());
PR(coll.maxW());
PR(coll.minT());
PR(coll.maxT());
#endif
#if defined(TEST_GENERATOR)
Sartre sartre;
sartre.init("sartreRuncard.txt");
Event *myEvent = sartre.generateEvent();
myEvent->list();
PR(sartre.totalCrossSection());
#endif
#if defined(TEST_INTERPOLATE)
TableCollection coll;
coll.init(1, bSat, 22);
vector vector_L;
ifstream ifs("../testing/dvcs/randomPoints_L.txt");
double t, Q2, W2, val;
while (ifs.good() && !ifs.eof()) {
ifs >> t >> Q2 >> W2 >> val;
if (ifs.eof()) break;
InterTestPoint point;
point.t = t;
point.Q2 = Q2;
point.W2 = W2;
point.value = val;
vector_L.push_back(point);
}
vector vector_L2;
ifstream ifs2("../testing/dvcs/randomPoints_L2.txt");
while (ifs2.good() && !ifs2.eof()) {
ifs2 >> t >> Q2 >> W2 >> val;
if (ifs2.eof()) break;
InterTestPoint point;
point.t = t;
point.Q2 = Q2;
point.W2 = W2;
point.value = val;
vector_L2.push_back(point);
}
vector vector_T;
ifstream ifs3("../testing/dvcs/randomPoints_T.txt");
while (ifs3.good() && !ifs3.eof()) {
ifs3 >> t >> Q2 >> W2 >> val;
if (ifs3.eof()) break;
InterTestPoint point;
point.t = t;
point.Q2 = Q2;
point.W2 = W2;
point.value = val;
vector_T.push_back(point);
}
vector vector_T2;
ifstream ifs4("../testing/dvcs/randomPoints_T2.txt");
while (ifs4.good() && !ifs4.eof()) {
ifs4 >> t >> Q2 >> W2 >> val;
if (ifs4.eof()) break;
InterTestPoint point;
point.t = t;
point.Q2 = Q2;
point.W2 = W2;
point.value = val;
vector_T2.push_back(point);
}
PR(vector_T.size());
PR(vector_T2.size());
PR(vector_L.size());
PR(vector_L2.size());
TFile *hfile = new TFile("interTest.root","RECREATE");
TH1D histoL("histoL", "L", 100, -0.1, 0.1);
for (unsigned int i=0; iWrite();
hfile->Close();
#endif
#if defined(TEST_SLOPE)
//
// For this test need to remove jacobian in CrossSection::logDerivateOfAmplitude()
// and make the method public
//
TableCollection coll;
coll.init(1, bSat, 443);
CrossSection xSection(&coll);
vector vector_L;
ifstream ifs("randomPoints_dAdW2L.txt");
double t, Q2, W2, val;
while (ifs.good() && !ifs.eof()) {
ifs >> t >> Q2 >> W2 >> val;
if (ifs.eof()) break;
InterTestSlope point;
point.t = t;
point.Q2 = Q2;
point.W2 = W2;
point.slope = val;
vector_L.push_back(point);
}
vector vector_T;
ifstream ifs2("randomPoints_dAdW2T.txt");
while (ifs2.good() && !ifs2.eof()) {
ifs2 >> t >> Q2 >> W2 >> val;
if (ifs2.eof()) break;
InterTestSlope point;
point.t = t;
point.Q2 = Q2;
point.W2 = W2;
point.slope = val;
vector_T.push_back(point);
}
PR(vector_T.size());
PR(vector_L.size());
TFile *hfile = new TFile("slopeTest.root","RECREATE");
TH1D histoL("histoL", "L", 100, -0.1, 0.1);
for (unsigned int i=0; iWrite();
hfile->Close();
#endif
#if defined(TEST_CROSSSECTION)
Sartre sartre;
EventGeneratorSettings* settings = sartre.runSettings();
bool useCorrections = true;
settings->setVerbose(true);
settings->setVerboseLevel(1);
settings->setNumberOfEvents(0);
settings->setTimesToShow(0);
settings->setQ2min(1);
settings->setQ2max(2);
settings->setWmin(64);
settings->setWmax(65);
settings->setVectorMesonId(333);
settings->setElectronBeamEnergy(20);
settings->setHadronBeamEnergy(100);
settings->setA(1);
settings->setDipoleModelType(bNonSat);
settings->setCorrectForRealAmplitude(true);
settings->setCorrectSkewedness(true);
settings->setEnableNuclearBreakup(false);
settings->setMaxLambdaUsedInCorrections(0.2);
// settings->list();
bool ok = sartre.init();
if (!ok) {
cout << "Initialization of sartre failed." << endl;
return 0;
}
// Normalize Sartre cross-section
double sartreCS = sartre.totalCrossSection();
cout << "sartreCS = " << sartreCS << endl;
-
+
#endif
#if defined(TEST_GLAUBER)
+ Nucleus deuteron(2);
+ Nucleus lead(208);
+ lead.normalizationOfT();
+ deuteron.normalizationOfT();
+ /*
TFile *hfile = new TFile("glauberTest.root","RECREATE");
TH1D h1("h1","Density", 300, 0, 15.);
TH2D h2xy("h2xy","xy", 60, -15, 15., 60, -15., 15.);
TH2D h2xz("h2xz","xz", 60, -15, 15., 60, -15., 15.);
TH2D h2yz("h2yz","yz", 60, -15, 15., 60, -15., 15.);
TH2D h2xy_1("h2xy_1","xy", 60, -15, 15., 60, -15., 15.);
TH2D h2xz_1("h2xz_1","xz", 60, -15, 15., 60, -15., 15.);
TH2D h2yz_1("h2yz_1","yz", 60, -15, 15., 60, -15., 15.);
TH2D h2xy_2("h2xy_2","xy", 60, -15, 15., 60, -15., 15.);
TH2D h2xz_2("h2xz_2","xz", 60, -15, 15., 60, -15., 15.);
TH2D h2yz_2("h2yz_2","yz", 60, -15, 15., 60, -15., 15.);
TH2D h2xy_3("h2xy_3","xy", 60, -15, 15., 60, -15., 15.);
TH2D h2xz_3("h2xz_3","xz", 60, -15, 15., 60, -15., 15.);
TH2D h2yz_3("h2yz_3","yz", 60, -15, 15., 60, -15., 15.);
double binsize = h1.GetBinWidth(1);
NewNucleus nucleus(208);
vector nucleons;
const unsigned int numberOfNuclei = 100000;
for (unsigned int i=0; i < numberOfNuclei; i++) {
if (i%10000 == 0) cout << i << endl;
nucleus.generate();
nucleons = nucleus.configuration();
for (unsigned int k=0; kWrite();
hfile->Close();
cout << "Histograms written to file 'glauberTest.root'" << endl;
+ */
#endif
return 0;
}
Index: trunk/examples/createBSatBDependenceTable.cpp
===================================================================
--- trunk/examples/createBSatBDependenceTable.cpp (revision 370)
+++ trunk/examples/createBSatBDependenceTable.cpp (revision 371)
@@ -1,132 +1,132 @@
//==============================================================================
// createBSatBDependenceTable.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Main program to create impact parameter lookup table.
// [Developer only]
//==============================================================================
#include
#include
#include
#include
#include
#include "TableGeneratorNucleus.h"
#include "TableGeneratorSettings.h"
#include "TFile.h"
#include "TH2F.h"
#include "TVector3.h"
#include "Constants.h"
#include "DipoleModelParameters.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
double overlapFunctionT(double);
void createLookUpTableOnFile(unsigned int, unsigned int);
DipoleModelParameters *gParameters;
int main(int argc, const char* argv[])
{
if (argc != 3) {
cout << "Usage: " << argv[0] << " A numberOfConfigurations" << endl;
return 2;
}
gParameters = new DipoleModelParameters(bSat, KMW);
//To explicitly delete all new histograms myself:
TH1::AddDirectory(false);
unsigned int A=atoi(argv[1]);
unsigned int numConf=atoi(argv[2]);
cout<<"Create table for A="<BG(); // GeV^-2
double arg = (b*b/hbarc2) / (2*BG);
return 1/(2*M_PI*BG) * exp(-arg);
}
void createLookUpTableOnFile(unsigned int A, unsigned int numberOfConfigurations)
{
//
// A function to create a lookup table of nucleon configurations and save them on a file
// This function is supposed to be independent from the rest of the code
// and be called in the Main program with the nucleus number as an argument.
//
TableGeneratorNucleus myNucleus(A);
// myNucleus.init(A);
TableGeneratorSettings::instance()->setSeed(42); // sets also seed for ROOT
vector hConfigurations(numberOfConfigurations);
int numBbins=1e3, numAngleBins=1e3;
double bRange=upperIntegrationLimit*myNucleus.radius()*1.1, angleRange=2*M_PI;
ostringstream histoName;
for(unsigned int iNuclei=0; iNucleiGetXaxis()->GetBinCenter(ib);
for(int iphi=1; iphi <= numAngleBins; iphi++) {
double phi=hConfigurations[iNuclei]->GetYaxis()->GetBinCenter(iphi);
TVector3 bvector=TVector3(b*cos(phi), b*sin(phi), 0.);
double sumOfT=0.;
for(unsigned int iA=0; iA < A; iA++){
sumOfT+=overlapFunctionT((bvector-myNucleus.configuration.at(iA).position()).Perp());
} //for iA
hConfigurations[iNuclei]->SetBinContent(ib, iphi, sumOfT);
}//for iphi
} //for ib
if(iNuclei % (numberOfConfigurations/10) == 0){
cout<Write();
delete hConfigurations[i];
hConfigurations[i] = 0;
}
lufile->Close();
cout<<" done, bye!"<.
//
// 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
#include
#include
#include
#include
#include
#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,
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,
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="<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::epsilon();
hplus -= numeric_limits::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(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
<< " s" << endl;
cout << "Total time: " << static_cast(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/sartreMain.cpp
===================================================================
--- trunk/examples/sartreMain.cpp (revision 370)
+++ trunk/examples/sartreMain.cpp (revision 371)
@@ -1,291 +1,291 @@
//==============================================================================
// sartreMain.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Example main program. Use to get started.
//
//==============================================================================
#include
#include
#include "TTree.h"
#include "TFile.h"
#include "Sartre.h"
-#include "DipoleModelParameters.h" // TMP
+#include "DipoleModelParameters.h"
#include "TGenPhaseSpace.h"
#include "Settings.h"
#include "TH1D.h"
#include "TROOT.h"
#include "TLorentzVector.h"
#include "TClonesArray.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
struct rootSartreEvent {
double t;
double Q2;
double x;
double s;
double y;
double W;
double xpom;
int iEvent;
int pol; // 0=transverse or 1=longitudinal
int dmode; // 0=coherent, 1=Incoherent
};
rootSartreEvent myRootSartreEvent;
void randomlyReverseBeams(Event* ); // for UPC mode
int main(int argc, char *argv[])
{
//
// Check command line arguments
//
if (! (argc == 2 || argc == 3) ) {
cout << "Usage: " << argv[0] << " runcard [rootfile]" << endl;
return 2;
}
string runcard = argv[1];
//
// Create the generator and initialize it.
// Once initialized you cannot (should not) change
// the settings w/o re-initialing sartre.
//
Sartre sartre;
bool ok = sartre.init(runcard);
if (!ok) {
cerr << "Initialization of sartre failed." << endl;
return 1;
}
EventGeneratorSettings* settings = sartre.runSettings();
//
// ROOT file
// Use the one from command line unless not given
// in which case the one in the runcard is used.
//
string rootfile;
if (argc == 3) {
rootfile = argv[2];
settings->setRootfile(argv[2]);
}
else
rootfile = settings->rootfile();
//
// Print out all Sartre settings
//
settings->list();
TFile *hfile = 0;
if (rootfile.size()) {
hfile = new TFile(rootfile.c_str(),"RECREATE");
cout << "ROOT file is '" << rootfile.c_str() << "'." << endl;
}
//
// Setup ROOT tree
//
TLorentzVector *eIn = new TLorentzVector;
TLorentzVector *pIn = new TLorentzVector;
TLorentzVector *vm = new TLorentzVector;
TLorentzVector *eOut = new TLorentzVector;
TLorentzVector *pOut = new TLorentzVector;
TLorentzVector *gamma = new TLorentzVector;
TLorentzVector *vmDaughter1 = new TLorentzVector;
TLorentzVector *vmDaughter2 = new TLorentzVector;
TClonesArray protons("TLorentzVector");
TClonesArray neutrons("TLorentzVector");
TClonesArray remnants("TLorentzVector");
TTree tree("tree","sartre");
tree.Branch("event", &myRootSartreEvent.t,
"t/D:Q2/D:x/D:s/D:y/D:W/D:xpom/D:iEvent/I:pol/I:dmode/I");
tree.Branch("eIn", "TLorentzVector", &eIn, 32000, 0);
tree.Branch("pIn", "TLorentzVector", &pIn, 32000, 0);
tree.Branch("vm", "TLorentzVector", &vm, 32000, 0);
tree.Branch("eOut", "TLorentzVector", &eOut, 32000, 0);
tree.Branch("pOut", "TLorentzVector", &pOut, 32000, 0);
tree.Branch("gamma","TLorentzVector", &gamma, 32000, 0);
tree.Branch("vmDaughter1", "TLorentzVector", &vmDaughter1, 32000, 0);
tree.Branch("vmDaughter2", "TLorentzVector", &vmDaughter2, 32000, 0);
if(settings->enableNuclearBreakup()){
tree.Branch("protons", &protons);
tree.Branch("neutrons", &neutrons);
tree.Branch("nuclearRemnants", &remnants);
}
//
// Prepare event generation
//
TGenPhaseSpace decay; // for VM decays
int daughterID = settings->userInt();
double daughterMasses[2] = {0, 0};
bool doPerformDecay = false;
if (daughterID && settings->vectorMesonId() != 22) {
doPerformDecay = true;
daughterMasses[0] = settings->lookupPDG(daughterID)->Mass();
daughterMasses[1] = settings->lookupPDG(-daughterID)->Mass();
cout << "Will decay vector meson: ";
cout << settings->lookupPDG(settings->vectorMesonId())->GetName();
cout << " -> ";
cout << settings->lookupPDG(daughterID)->GetName();
cout << " ";
cout << settings->lookupPDG(-daughterID)->GetName();
cout << endl;
}
int nPrint;
if (settings->timesToShow())
nPrint = settings->numberOfEvents()/settings->timesToShow();
else
nPrint = 0;
unsigned long maxEvents = settings->numberOfEvents();
cout << "Generating " << maxEvents << " events." << endl << endl;
//
// Event generation
//
for (unsigned long iEvent = 0; iEvent < maxEvents; iEvent++) {
//
// Generate one event
//
Event *event = sartre.generateEvent();
if (nPrint && (iEvent+1)%nPrint == 0 && iEvent != 0) {
cout << "processed " << iEvent+1 << " events" << endl;
}
//
// If Sartre is run in UPC mode, half of the events needs to be
// rotated around and axis perpendicular to z:
//
if(settings->UPC() and settings->A()==settings->UPCA()){
randomlyReverseBeams(event);
}
//
// Print out (here only for the first few events)
//
if (iEvent < 10) event->list();
//
// Fill ROOT tree
//
myRootSartreEvent.iEvent = event->eventNumber;
myRootSartreEvent.t = event->t;
myRootSartreEvent.Q2 = event->Q2;
myRootSartreEvent.x = event->x;
myRootSartreEvent.y = event->y;
myRootSartreEvent.s = event->s;
myRootSartreEvent.W = event->W;
myRootSartreEvent.xpom = event->xpom;
myRootSartreEvent.pol = event->polarization == transverse ? 0 : 1;
myRootSartreEvent.dmode = event->diffractiveMode == coherent ? 0 : 1;
eIn = &event->particles[0].p;
pIn = &event->particles[1].p;
eOut = &event->particles[2].p;
pOut = &event->particles[6].p;
vm = &event->particles[4].p;
gamma = &event->particles[3].p;
//If the event is incoherent, and nuclear breakup is enabled, fill the remnants to the tree
if(settings->enableNuclearBreakup() and event->diffractiveMode == incoherent){
protons.Clear();
neutrons.Clear();
remnants.Clear();
for(unsigned int iParticle=7; iParticle < event->particles.size(); iParticle++){
if(event->particles[iParticle].status == 1) { // Final-state particle
const Particle& particle = event->particles[iParticle];
switch (abs(particle.pdgId)) {
case 2212: // (Anti-)proton
new(protons[protons.GetEntries()]) TLorentzVector(particle.p);
break;
case 2112: // (Anti-)neutron
new(neutrons[neutrons.GetEntries()]) TLorentzVector(particle.p);
break;
default: // Any other remnant
new(remnants[remnants.GetEntries()]) TLorentzVector(particle.p);
break;
} // switch
} // if
} // for
} // if
//
// Decay the vector meson and fill the decay roducts in the tree:
//
if(doPerformDecay) {
if( decay.SetDecay(*vm, 2, daughterMasses) ){
double weight = decay.Generate(); // weight is always 1 here
if ( (weight-1) > FLT_EPSILON) {
cout << "sartreMain: Warning weight != 1, weight = " << weight << endl;
}
vmDaughter1 = decay.GetDecay(0);
vmDaughter2 = decay.GetDecay(1);
}
else {
cout << "sartreMain: Warning: Kinematics of Vector Meson does not allow decay!" << endl;
}
}
tree.Fill();
}
cout << "All events processed\n" << endl;
//
// That's it, finish up
//
double totalCS=sartre.totalCrossSection();
TH1D* histoForCSandNumberOfEvents = new TH1D("histoForCSandNumberOfEvents", "Cross-section and Number of Events", 2, 0., 1.);
histoForCSandNumberOfEvents->SetBinContent(1, totalCS);
histoForCSandNumberOfEvents->SetBinContent(2, maxEvents);
double runTime = sartre.runTime();
hfile->Write();
cout << rootfile.c_str() << " written." << endl;
cout << "Total cross-section: " << totalCS << " nb" << endl;
sartre.listStatus();
cout << "CPU Time/event: " << 1000*runTime/maxEvents << " msec/evt" << endl;
return 0;
}
// UPC only
void randomlyReverseBeams(Event* myEvent) {
TRandom3 *random = EventGeneratorSettings::randomGenerator();
if(random->Uniform(1) > 0.5){
for(unsigned int i=0; iparticles.size(); i++)
myEvent->particles.at(i).p.RotateX(M_PI);
}
}
Index: trunk/examples/convergenceTestRuncard.txt
===================================================================
--- trunk/examples/convergenceTestRuncard.txt (revision 370)
+++ trunk/examples/convergenceTestRuncard.txt (revision 371)
@@ -1,136 +1,136 @@
//==============================================================================
// tableGeneratorRuncard.txt
//
-// Copyright (C) 2010-2016 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Tobias Toll
// Last update:
// $Date: 2017-03-01 01:36:16 +0530 (Wed, 01 Mar 2017) $
// $Author: ullrich $
//==============================================================================
//
// Example Runcard for Sartre Table Generator.
//
// Comments start with a # or #
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Nucleus number
#
A = 2
#
# Which moments of the amplitude to calculate:
# 0 analytical, averaged over configurations
# 1 only, analytical
# 2 , and both averaged over configurations
#
modesToCalculate = 2
#
# Rootfile This is a prefix table type (T and L, or T2 and L2) and ".root" will be appended
#
numberOfConfigurations = 100
rootfile = test_100.root
#
# Model parameters
#
vectorMesonId = 443
dipoleModel = bSat
dipoleModelParameterSet = KMW
#
# Table Set Type (experts only)
#
tableSetType = total_and_coherent
#tableSetType = coherent_and_incoherent
#
# Developer only
# Only takes affect when dipoleModelParameterSet = CUSTOM
# All parameters must be given.
#
# For bSat and bNonSat (8 parameters):
# m_u m_d m_s m_c BG mu02 lambdaG Ag
# For bCGC (10 parameters):
# m_u m_d m_s m_c dkappa dN0 dx0 dlambda dgammas dBcgc
#
#dipoleModelCustomParameters = 0.14 0.14 0.14 1.27 4. 1.51 0.058 2.308
# Path for the b-dependence lookup table for bSat. Not used with bCGC.
bSatLookupPath = /home/tollto/sartreSVN/sartre/bSatLookupTables/
#
# Kinematic range
#
# EIC limits:
# Say we aim s = 20000 (30000 for pp)
# That implies:
# Q2 = 0.0001 - 20000 (or with x<1e-2, Q2<200)
# t = -3 - 0
# W2 = 1 - 20000 (W = 1-141)
#
Q2min = 0.01
Q2max = 20.
Wmin = 10
Wmax = 380
tmin = -.5
tmax = -.0
#H1 DVCS:
#Q2min = 6.5
#Q2max = 80.
#Wmin = 30
#Wmax = 140
#tmin = -1.
#tmax = -.0
#
# Number of bins in the tables:
#
#
# ZEUS test, bin density for
# 40x40x40: 64000/(99.85*190*1)=3.37 bins/GeV^5
#
# Or in different directions:
# Q2: 40/100=0.4 bins/GeV^2
# W: 40/190=0.21 bins/GeV
# t: 40/1= 40 bins/GeV^2
#
Q2bins = 20
W2bins = 20
tbins = 50
#
# If a run fails, use the backup and continue:
# If in doubt use the tableInspector to make sure the rest
# of the runcard is as it should be.
#
useBackupFile = false
startingBinFromBackup=2
#
# If the table you are making should take presendence
# over previous tables, make it with a higher priority
#
priority = 0
Index: trunk/examples/tableGeneratorRuncard.txt
===================================================================
--- trunk/examples/tableGeneratorRuncard.txt (revision 370)
+++ trunk/examples/tableGeneratorRuncard.txt (revision 371)
@@ -1,135 +1,135 @@
//==============================================================================
// tableGeneratorRuncard.txt
//
-// Copyright (C) 2010-2016 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Example Runcard for Sartre Table Generator.
//
// Comments start with a # or #
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Nucleus number
#
A = 197
#
# Which moments of the amplitude to calculate:
# 0 analytical, averaged over configurations
# 1 only, analytical
# 2 , and both averaged over configurations
#
modesToCalculate = 2
#
# Rootfile This is a prefix table type (T and L, or T2 and L2) and ".root" will be appended
#
rootfile = bintests
#
# Model parameters
#
vectorMesonId = 443
dipoleModel = bSat
dipoleModelParameterSet = KMW
#
# Table Set Type (experts only)
#
tableSetType = total_and_coherent
#tableSetType = coherent_and_incoherent
#
# Developer only
# Only takes affect when dipoleModelParameterSet = CUSTOM
# All parameters must be given.
#
# For bSat and bNonSat (8 parameters):
# m_u m_d m_s m_c BG mu02 lambdaG Ag
# For bCGC (10 parameters):
# m_u m_d m_s m_c dkappa dN0 dx0 dlambda dgammas dBcgc
#
#dipoleModelCustomParameters = 0.14 0.14 0.14 1.27 4. 1.51 0.058 2.308
# Path for the b-dependence lookup table for bSat. Not used with bCGC.
bSatLookupPath = ./
numberOfConfigurations = 100
#
# Kinematic range
#
# EIC limits:
# Say we aim s = 20000 (30000 for pp)
# That implies:
# Q2 = 0.0001 - 20000 (or with x<1e-2, Q2<200)
# t = -3 - 0
# W2 = 1 - 20000 (W = 1-141)
#
Q2min = 0.01
Q2max = 20.
Wmin = 10
Wmax = 380
tmin = -.5
tmax = -.0
#H1 DVCS:
#Q2min = 6.5
#Q2max = 80.
#Wmin = 30
#Wmax = 140
#tmin = -1.
#tmax = -.0
#
# Number of bins in the tables:
#
#
# ZEUS test, bin density for
# 40x40x40: 64000/(99.85*190*1)=3.37 bins/GeV^5
#
# Or in different directions:
# Q2: 40/100=0.4 bins/GeV^2
# W: 40/190=0.21 bins/GeV
# t: 40/1= 40 bins/GeV^2
#
Q2bins = 20
W2bins = 20
tbins = 50
#
# If a run fails, use the backup and continue:
# If in doubt use the tableInspector to make sure the rest
# of the runcard is as it should be.
#
useBackupFile = false
startingBinFromBackup=2
#
# If the table you are making should take presendence
# over previous tables, make it with a higher priority
#
priority = 0
Index: trunk/examples/tableGeneratorUPCMain.cpp
===================================================================
--- trunk/examples/tableGeneratorUPCMain.cpp (revision 370)
+++ trunk/examples/tableGeneratorUPCMain.cpp (revision 371)
@@ -1,417 +1,417 @@
//==============================================================================
// tableGeneratorMain.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// 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 .
//
// 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
#include
#include
#include
#include
#include
#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,
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="<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(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
<< " s" << endl;
cout << "Total time: " << static_cast(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;
}
Index: trunk/examples/sartreRuncard.txt
===================================================================
--- trunk/examples/sartreRuncard.txt (revision 370)
+++ trunk/examples/sartreRuncard.txt (revision 371)
@@ -1,124 +1,140 @@
//==============================================================================
// sartreRuncard.txt
//
-// Copyright (C) 2010-2016 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Example Runcard for Sartre Event Generator.
//
// Comments start with a # or //
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Define beams
#
eBeamEnergy = 20
hBeamEnergy = 100
-A = 197
+A = 1
#
# UPC settings, to run in UPC mode set UPC=true and UPCA into the photon emitting species
#
-#UPC=false
+UPC=false
#UPCA=197
#xpomMin = 1e-5
#xpomMax = 1
#
# Number of events and printout frequency
#
numberOfEvents = 10000
timesToShow = 20
#
# Set verbosity
#
verbose = true
verboseLevel = 2
-
#
# Rootfile
#
rootfile = example.root
#
# Dipole model
#
+# dipoleModel: can be bSat or bNonSat
+# dipoleModelParameterSet:
+# KMW: older set but most complete
+# HMPZ: new parameters with still a limited set
+# of tables. For eA use KMW.
+#
dipoleModel = bSat
-dipoleModelParameterSet = KMW
+dipoleModelParameterSet = HMPZ
#
-# Model parameters
-#
-# vectorMesonID: 22, 113, 333, 443
+# Vector mesons to create.
+# Note that not all vector mesons
+# have tables for each dipole model and
+# available A.
+#
+# vectorMesonID:
+# 22 DVCS
+# 113 rho
+# 333 phi
+# 443 J/psi
+# 553 Upsilon
#
vectorMesonId = 443
#
# Table Set Type (experts only)
#
tableSetType = total_and_coherent
#tableSetType = coherent_and_incoherent
#
-# User variable used for vector meson decays
+# User variable
+# If you use the provided sartreMain.cpp this is
+# used for vector meson decays
# PDG: pi+ = 211, K+ = 321, e- = 11, mu- = 13
+# with BR=100%
#
userInt = 11
#
# Kinematic range min > max means no limits (given by table range)
#
Q2min = 1000000
Q2max = 0
Wmin = 1000000
Wmax = 0
#
# Corrections
#
correctForRealAmplitude = true
correctSkewedness = true
-# maxLambdaUsedInCorrections = 0.65
+maxLambdaUsedInCorrections = 0.4
#
# Misc
#
enableNuclearBreakup = false
maxNuclearExcitationEnergy = 3.0
#
# Random generator seed (if not given current time is used)
#
#seed = 2011987
#
# User parameters
#
# userDouble = 0.
# userString = "Hello World!"
# userInt = 0
#
# Expert flags
#
# applyPhotonFlux = true
Index: trunk/examples/nucleusAnimation.cpp
===================================================================
--- trunk/examples/nucleusAnimation.cpp (revision 370)
+++ trunk/examples/nucleusAnimation.cpp (revision 371)
@@ -1,300 +1,300 @@
//==============================================================================
// nucleusAnimation.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include
#include
#include
#include "TVector3.h"
#include
#include
#include
#include "Nucleus.h"
#include "TableGeneratorNucleus.h"
#include
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
// actual vector representing the camera's direction
float lx=0.0f,lz=-1.0f;
// XZ position of the camera
float x=0.0f, y=40.0f, z=5.0f;
// the key states. These variables will be zero
// when no key is being presses
float deltalr=0.0f, deltaud=0.0f;
bool newNucleus=0;
float deltazoom=0.0f;
TableGeneratorNucleus* myNucleus;// = new TableGeneratorNucleus;
void changeSize(int w, int h) {
// Prevent a divide by zero, when window is too short
// (you cant make a window of zero width).
if (h == 0)
h = 1;
float ratio = w * 1.0 / h;
// Use the Projection Matrix
glMatrixMode(GL_PROJECTION);
// Reset Matrix
glLoadIdentity();
// Set the viewport to be the entire window
glViewport(0, 0, w, h);
// Set the correct perspective.
gluPerspective(45.0f, ratio, 0.1f, 100.0f);
// Get Back to the Modelview
glMatrixMode(GL_MODELVIEW);
}
void drawNucleon() {
// glColor3f(1.0f, .0f, .0f);
glTranslatef(0.0f ,0.75f, 0.0f);
glutSolidSphere(0.7f,20,20);
}
void computeDir(){
double sud=sin(deltaud), cud=cos(deltaud),
slr=sin(deltalr), clr=cos(deltalr);
double zoom=1.+deltazoom;
// theta=lr, phi=ud, psi=0;
x= (clr*x + sud*slr*y + cud*slr*z)*zoom;
y= ( cud*y + (-sud)*z)*zoom;
z=(-slr*x + sud*clr*y + cud*clr*z)*zoom;
}
// Position the light at the origin.
const GLfloat light_pos[] = { 0.0, 50.0, 0.0, 1.0 };
// Attenuation factors for light.
GLfloat const_att = 1.0;
void renderScene(void) {
if (deltalr or deltaud or deltazoom)
computeDir();
if( newNucleus )
while(!myNucleus->generate()){};
// Clear Color and Depth Buffers
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
// Reset transformations
glLoadIdentity();
// Set the camera
gluLookAt(x, y, z, // camera position
0., 0., 0., //camera focus
0.0f, 1., 0.0f); //Vector pointing up
for(unsigned int i=0; iconfiguration.size(); i++){
TVector3 pos=myNucleus->configuration.at(i).position();
glPushMatrix();
glLightf(GL_LIGHT0, GL_CONSTANT_ATTENUATION, const_att);
glLightfv(GL_LIGHT0, GL_POSITION, light_pos);
glTranslatef(pos.X(),pos.Y(),pos.Z());
drawNucleon();
glPopMatrix();
}
glutSwapBuffers();
}
void idle(void){
glutPostRedisplay();
}
void pressKey(int key, int xx, int yy) {
xx=yy;
switch (key) {
case GLUT_KEY_LEFT : deltalr = -0.01; break;
case GLUT_KEY_RIGHT : deltalr = 0.01; break;
case GLUT_KEY_UP : deltaud = 0.01; break;
case GLUT_KEY_DOWN : deltaud = -0.01; break;
case 27 :
case 81 :
case 113: exit(0); break;
case 110:
case 78: newNucleus=true; break;
case 122:
case 90: deltazoom = 0.01; break;
case 65:
case 97: deltazoom = -0.01; break;
}
}
void releaseKey(int key, int x, int y) {
x=y;
switch (key) {
case GLUT_KEY_LEFT :
case GLUT_KEY_RIGHT : deltalr = 0.0; break;
case GLUT_KEY_UP :
case GLUT_KEY_DOWN : deltaud = 0; break;
case 110:
case 78: newNucleus=false; break;
case 122:
case 90:
case 65:
case 97: deltazoom = 0; break;
}
}
void exitFunction(){
delete myNucleus;
}
void usage(const char* prog)
{
cout << "Usage: " << prog << "[-w] A" << endl;
cout << " Valid A are: 208, 197, 110, 63, 40, 27, 16, 1" << endl;
cout << " -w White instead of black background" << endl;
}
int main(int argc, char **argv) {
atexit(exitFunction);
if (argc < 2) {
usage(argv[0]);
return 2;
}
bool hasWhiteBackground = false;
int ch;
while ((ch = getopt(argc, argv, "w")) != -1) {
switch (ch) {
case 'w':
hasWhiteBackground = true;
break;
default:
usage(argv[0]);
return 2;
break;
}
}
if (optind == argc) {
usage(argv[0]);
return 2;
}
string str(argv[optind]);
for (unsigned int i=0; igenerate()){};
cout<<"List of Commands:"<.
//
// Author: Tobias Toll
// Last update:
// $Date: 2017-02-28 15:06:16 -0500 (Tue, 28 Feb 2017) $
// $Author: ullrich $
//==============================================================================
//
// Example Runcard for Sartre Table Generator.
//
// Comments start with a # or #
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Are tables for ultra-peripheral collisions
#
#
UPC=true
# Nucleus number
#
A = 1
#
# Which moments of the amplitude to calculate:
# 0 analytical, averaged over configurations
# 1 only, analytical
# 2 , and both averaged over configurations
#
modesToCalculate = 0
#
# Rootfile This is a prefix table type (T and L, or T2 and L2) and ".root" will be appended
#
rootfile = test
#
# Model parameters
#
vectorMesonId = 443
dipoleModel = bSat
dipoleModelParameterSet = HMPZ
#
# Table Set Type (experts only)
#
#tableSetType = total_and_coherent
tableSetType = coherent_and_incoherent
#
# Developer only
# Only takes affect when dipoleModelParameterSet = CUSTOM
# All parameters must be given.
#
# For bSat and bNonSat (8 parameters):
# m_u m_d m_s m_c BG mu02 lambdaG Ag
# For bCGC (10 parameters):
# m_u m_d m_s m_c dkappa dN0 dx0 dlambda dgammas dBcgc
#
#dipoleModelCustomParameters = 0.14 0.14 0.14 1.27 4. 1.51 0.058 2.308
# Path for the b-dependence lookup table for bSat. Not used with bCGC.
bSatLookupPath = /home/tollto/SNU/sartre-inclusive165/bSatLookupTables/GaussBg4/
numberOfConfigurations = 1
#
# Kinematic range
#
xmin = 1e-8
xmax = 2e-2
tmin = -.5
tmax = -.0
#
# Number of bins in the tables:
#
#
xbins = 5
tbins = 5
#
# If a run fails, use the backup and continue:
# If in doubt use the tableInspector to make sure the rest
# of the runcard is as it should be.
#
useBackupFile = false
startingBinFromBackup=2
#
# If the table you are making should take presendence
# over previous tables, make it with a higher priority
#
priority = 8
Index: trunk/examples/sartreUPCRuncard.txt
===================================================================
--- trunk/examples/sartreUPCRuncard.txt (revision 370)
+++ trunk/examples/sartreUPCRuncard.txt (revision 371)
@@ -1,124 +1,124 @@
//==============================================================================
// sartreRuncard.txt
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Thomas Ullrich
// Last update:
// $Date: 2018-06-20 16:34:51 -0400 (Wed, 20 Jun 2018) $
// $Author: ullrich $
//==============================================================================
//
// Example Runcard for Sartre Event Generator running in UPC mode
//
// Comments start with a # or //
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Define beams
#
# LHC UPC PbPb mode sqrt(s)=2.76 TeV
eBeamEnergy = 1380
hBeamEnergy = 1380
A = 208
# Examples:
# LHC UPC pPb mode
#eBeamEnergy = 1575
#hBeamEnergy = 4000
#A = 1
#
# UPC settings, to run in UPC mode set UPC=true and UPCA into the photon emitting species
#
UPC=true
UPCA=208
xpomMin = 1 // min > max means no user limits
xpomMax = 0
#
# Number of events and printout frequency
#
numberOfEvents = 10000
timesToShow = 20
#
# Set verbosity
#
verbose = true
verboseLevel = 2
#
# Rootfile
#
rootfile = example.root
#
# Dipole model
#
dipoleModel = bSat
dipoleModelParameterSet = HMPZ
#
# Model parameters
#
# vectorMesonID: 22, 113, 333, 443
#
vectorMesonId = 443
#
# Table Set Type (experts only)
#
tableSetType = coherent_and_incoherent
#
# User variable used for vector meson decays
# PDG: pi+ = 211, K+ = 321, e- = 11, mu- = 13
#
userInt = 11
#
# Corrections
#
correctForRealAmplitude = true
correctSkewedness = true
maxLambdaUsedInCorrections = 0.1
#
# Misc
#
enableNuclearBreakup = false
maxNuclearExcitationEnergy = 3.0
#
# Random generator seed (if not given current time is used)
#
seed = 2011987
#
# User parameters
#
# userDouble = 0.
# userString = "Hello World!"
# userInt = 0
#
# Expert flags
#
# applyPhotonFlux = true
Index: trunk/examples/heraCompare.cpp
===================================================================
--- trunk/examples/heraCompare.cpp (revision 370)
+++ trunk/examples/heraCompare.cpp (revision 371)
@@ -1,1702 +1,1702 @@
//==============================================================================
// heraCompare.cpp
//
-// Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich
+// 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 .
//
// Author: Thomas Ullrich and Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Main program to compare Sartre with HERA data (J/psi and DVCS)
// Usage:
// heraCompare rootfile jpsi|phi|dvcs
//
// Output are graphs with HERA and Sartre data to compare.
//==============================================================================
#include
#include
#include
#include
#include "Sartre.h"
#include "TFile.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "Constants.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
double betaMin=0.7;
double xpomMax=1e-2;
double maxLam=0.35;
DipoleModelType theDipoleModelType = bNonSat;
DipoleModelParameterSet theDipoleModelParameterSet = HMPZ;
TableSetType theTableSetType = total_and_coherent; // or coherent_and_incoherent;
struct HeraDataPoint {
double Q2;
double W;
double t;
double tmin;
double tmax;
double Q2min;
double Q2max;
double Wmin;
double Wmax;
double xsection;
double xsection_stat;
double xsection_sys;
double xsection_sys_up;
double xsection_sys_lo;
HeraDataPoint();
};
HeraDataPoint::HeraDataPoint()
{
t=-1;
tmin=-1; tmax=-1;
Q2=-1;
Q2min=-1; Q2max=-1;
W=-1;
Wmin=-1; Wmax=-1;
xsection=-1;
xsection_stat=-1;
xsection_sys=-1;
xsection_sys_up=-1; xsection_sys_lo=-1;
}
void loadJpsiData(vector&);
void loadPhiData(vector&);
void loadDVCSData(vector&, vector&);
void compareJpsi(string&);
void compareDVCS(string&);
void comparePhi(string&);
double chi2_global = 0;
int ndof_global = 0;
int main(int argc, char** argv)
{
//
// Check command line arguments
//
if (argc != 3) {
cout << "Usage: " << argv[0] << " rootfile jpsi|phi|dvcs|all" << endl;
return 2;
}
//
// Run the requested comparison(s)
//
string rootfile = argv[1];
string what = argv[2];
if (what.find("jpsi") != string::npos)
compareJpsi(rootfile);
else if (what.find("dvcs") != string::npos)
compareDVCS(rootfile);
else if (what.find("phi") != string::npos)
comparePhi(rootfile);
else if (what.find("all") != string::npos){
compareDVCS(rootfile);
comparePhi(rootfile);
compareJpsi(rootfile);
cout << "-------------------------------------------------" << endl;
cout << "Global Comparison HERA-Sartre:" << endl;
cout << "chi2 = " << chi2_global << endl;
cout << "ndof = " << ndof_global << endl;
cout << "chi2/ndof = " << chi2_global/ndof_global << endl;
cout << "-------------------------------------------------" << endl;
}
else {
cout << "Error: '" << what.c_str() << "' no such probe to compare." << endl;
return 2;
}
return 0;
}
void comparePhi(string& rootfile)
{
//
// PHI HERA/ZEUS data
//
double heraProtonEnergy = 920;
double heraElectronEnergy = 27.6;
vector zeusDataPhi_table1;
loadPhiData(zeusDataPhi_table1);
vector graphVec;
vector sgraphVec;
vector sgraphName;
vector graphName;
// Table 1
TGraphErrors *graph1a = new TGraphErrors(7);
graph1a->SetTitle("Q2=2.4");
graphVec.push_back(graph1a);
graphName.push_back("graph1a");
TGraphErrors *graph1b = new TGraphErrors(6);
graph1b->SetTitle("Q2=3.8");
graphVec.push_back(graph1b);
graphName.push_back("graph1b");
TGraphErrors *graph1c = new TGraphErrors(6);
graph1c->SetTitle("Q2=6.5");
graphVec.push_back(graph1c);
graphName.push_back("graph1c");
TGraphErrors *graph1d = new TGraphErrors(6);
graph1d->SetTitle("Q2=13.0");
graphVec.push_back(graph1d);
graphName.push_back("graph1d");
TGraphErrors *sgraph1a = new TGraphErrors(7);
sgraph1a->SetTitle("Q2=2.4");
sgraphVec.push_back(sgraph1a);
sgraphName.push_back("sgraph1a");
TGraphErrors *sgraph1b = new TGraphErrors(6);
sgraph1b->SetTitle("Q2=3.8");
sgraphVec.push_back(sgraph1b);
sgraphName.push_back("sgraph1b");
TGraphErrors *sgraph1c = new TGraphErrors(6);
sgraph1c->SetTitle("Q2=6.5");
sgraphVec.push_back(sgraph1c);
sgraphName.push_back("sgraph1c");
TGraphErrors *sgraph1d = new TGraphErrors(6);
sgraph1d->SetTitle("Q2=13.0");
sgraphVec.push_back(sgraph1d);
sgraphName.push_back("sgraph1d");
// Table 3
TGraphErrors *graph3 = new TGraphErrors(9);
graph3->SetTitle("Q2");
graphVec.push_back(graph3);
graphName.push_back("graph3");
TGraphErrors *sgraph3 = new TGraphErrors(9);
sgraph3->SetTitle("Q2");
sgraphVec.push_back(sgraph3);
sgraphName.push_back("sgraph3");
//
// Loop over data points and for each initialize Sartre and
// calculate the cross-section. We are using no runcard but
// set thigs for each data point
//
Sartre sartre;
EventGeneratorSettings* settings = sartre.runSettings();
//
// Table 1
//
double diff_table1, chi2_table1 = 0;
double diff_table3, chi2_table3 = 0;
int ndof_table1 = 0;
int ndof_table3 = 0;
int graphCounter = 0;
int pointCounter = 0;
for (unsigned int i=0; iSetPoint(pointCounter, zeusDataPhi_table1[i].W, zeusDataPhi_table1[i].xsection);
if(graphCounter==4)
graphVec[graphCounter]->SetPoint(pointCounter, zeusDataPhi_table1[i].Q2, zeusDataPhi_table1[i].xsection);
double error = sqrt(zeusDataPhi_table1[i].xsection_stat*zeusDataPhi_table1[i].xsection_stat + zeusDataPhi_table1[i].xsection_sys*zeusDataPhi_table1[i].xsection_sys);
graphVec[graphCounter]->SetPointError(pointCounter, 0., error);
// set up sartre
settings->setVerbose(true);
settings->setVerboseLevel(2);
settings->setNumberOfEvents(0);
settings->setTimesToShow(0);
settings->setVectorMesonId(333);
settings->setApplyPhotonFlux(false);
settings->setElectronBeamEnergy(heraElectronEnergy);
settings->setHadronBeamEnergy(heraProtonEnergy);
settings->setDipoleModelType(theDipoleModelType);
settings->setDipoleModelParameterSet(theDipoleModelParameterSet);
settings->setA(1);
settings->setCorrectForRealAmplitude(true);
settings->setCorrectSkewedness(true);
settings->setMaxLambdaUsedInCorrections(maxLam);
settings->setEnableNuclearBreakup(false);
settings->setQ2min(zeusDataPhi_table1[i].Q2min); double Q2min=settings->Q2min();
settings->setQ2max(zeusDataPhi_table1[i].Q2max); double Q2max=settings->Q2max();
settings->setWmin(zeusDataPhi_table1[i].Wmin); double Wmin=settings->Wmin();
settings->setWmax(zeusDataPhi_table1[i].Wmax); double Wmax=settings->Wmax();
if(Wmax > 142 || Wmin > 142) continue;
//calculate in the t=0 limits, mp=0:
TParticlePDG *vectorMesonPDG = settings->lookupPDG(settings->vectorMesonId());
double vmMass = vectorMesonPDG->Mass();
double beta=Q2min/(vmMass*vmMass + Q2min);
double xpom=(vmMass*vmMass + Q2max) / (Wmin*Wmin + Q2max - protonMass2);
if(xpom > xpomMax) continue;
if(beta < betaMin) continue;
bool ok = sartre.init();
if (!ok) {
cout << "Initialization of sartre failed. Going to next data point." << endl;
continue;
}
settings->list();
//
// Sartre cross-section
//
vector > limits = sartre.kinematicLimits(); // t, Q2, W
double lower[3], upper[3];
lower[0] = -zeusDataPhi_table1[i].tmax; // t
upper[0] = -zeusDataPhi_table1[i].tmin;
lower[1] = limits[1].first; // Q2
upper[1] = limits[1].second;
lower[2] = limits[2].first; // W
upper[2] = limits[2].second;
double CS = sartre.totalCrossSection(lower, upper);
CS /= upper[0]-lower[0];
CS /= limits[1].second - limits[1].first;
CS /= limits[2].second*limits[2].second - limits[2].first*limits[2].first;
if(graphCounter<4){
sgraphVec[graphCounter]->SetPoint(pointCounter, zeusDataPhi_table1[i].W, CS);
sgraphVec[graphCounter]->SetPointError(pointCounter, 0., 0.);
}
if(graphCounter==4){
sgraphVec[graphCounter]->SetPoint(pointCounter, zeusDataPhi_table1[i].Q2, CS);
sgraphVec[graphCounter]->SetPointError(pointCounter, 0., 0.);
}
//
// chi2 calculation
//
if(graphCounter<4){
diff_table1 = CS - zeusDataPhi_table1[i].xsection;
chi2_table1 += (diff_table1*diff_table1) / (error*error);
ndof_table1++;
chi2_global+=(diff_table1*diff_table1) / (error*error);
ndof_global++;
}
if(graphCounter==4){
diff_table3 = CS - zeusDataPhi_table1[i].xsection;
chi2_table3 += (diff_table3*diff_table3) / (error*error);
ndof_table3++;
chi2_global+=(diff_table3*diff_table3) / (error*error);
ndof_global++;
}
//
// Print out
//
cout << endl;
cout << "dsig/dt at |t|=" << fabs(zeusDataPhi_table1[i].t)
<< ", Q2=" << zeusDataPhi_table1[i].Q2
<< ", W=" << zeusDataPhi_table1[i].W << endl;
cout << "\t\tSartre: " << CS << endl; // in nb/GeV2
cout << "\t\tZeus: " << zeusDataPhi_table1[i].xsection << "+/-" << zeusDataPhi_table1[i].xsection_stat << "+/-" <Write(graphName[i].c_str());
sgraphVec[i]->Write(sgraphName[i].c_str());
}
// Table 1
// graph1a->Write("graph1a");
// sgraph1a->Write("sgraph1a");
// for(int i=0; i<3; i++){
// graphVec.at(i)->Write("graph1a");
// sgraphVec.at(i)->Write("sgraph1a");
// }
hfile->Close();
}
void compareDVCS(string& rootfile)
{
//
// DVCS HERA/H1 data
//
double heraProtonEnergy = 920;
double heraElectronEnergy = 27.6;
vector h1DataDVCS_table4;
vector h1DataDVCS_table1;
loadDVCSData(h1DataDVCS_table1, h1DataDVCS_table4);
vector graphVec;
vector sgraphVec;
vector sgraphName;
vector graphName;
// Table 4
TGraphErrors *graph4a = new TGraphErrors(4); graph4a->SetTitle("Q2=8, W=40"); graphVec.push_back(graph4a); graphName.push_back("graph4a");
TGraphErrors *graph4b = new TGraphErrors(4); graph4b->SetTitle("Q2=8, W=70"); graphVec.push_back(graph4b); graphName.push_back("graph4b");
TGraphErrors *graph4c = new TGraphErrors(4); graph4c->SetTitle("Q2=8, W=100"); graphVec.push_back(graph4c); graphName.push_back("graph4c");
TGraphErrors *graph4d = new TGraphErrors(4); graph4d->SetTitle("Q2=20, W=40"); graphVec.push_back(graph4d); graphName.push_back("graph4d");
TGraphErrors *graph4e = new TGraphErrors(4); graph4e->SetTitle("Q2=20, W=70"); graphVec.push_back(graph4e); graphName.push_back("graph4e");
TGraphErrors *graph4f = new TGraphErrors(4); graph4f->SetTitle("Q2=20, W=100"); graphVec.push_back(graph4f); graphName.push_back("graph4f");
TGraphErrors *sgraph4a = new TGraphErrors(4); sgraph4a->SetTitle("Q2=8, W=40"); sgraphVec.push_back(sgraph4a); sgraphName.push_back("sgraph4a");
TGraphErrors *sgraph4b = new TGraphErrors(4); sgraph4b->SetTitle("Q2=8, W=70"); sgraphVec.push_back(sgraph4b); sgraphName.push_back("sgraph4b");
TGraphErrors *sgraph4c = new TGraphErrors(4); sgraph4c->SetTitle("Q2=8, W=100"); sgraphVec.push_back(sgraph4c); sgraphName.push_back("sgraph4c");
TGraphErrors *sgraph4d = new TGraphErrors(4); sgraph4d->SetTitle("Q2=20, W=40"); sgraphVec.push_back(sgraph4d); sgraphName.push_back("sgraph4d");
TGraphErrors *sgraph4e = new TGraphErrors(4); sgraph4e->SetTitle("Q2=20, W=70"); sgraphVec.push_back(sgraph4e); sgraphName.push_back("sgraph4e");
TGraphErrors *sgraph4f = new TGraphErrors(4); sgraph4f->SetTitle("Q2=20, W=100"); sgraphVec.push_back(sgraph4f); sgraphName.push_back("sgraph4f");
// Table 1
TGraphErrors *graph1a = new TGraphErrors(4); graph1a->SetTitle("30 < W < 140, |t|<1");
TGraphErrors *sgraph1a = new TGraphErrors(4); sgraph1a->SetTitle("30 < W < 140, |t|<1");
TGraphErrors *graph1b = new TGraphErrors(5); graph1b->SetTitle("6.5 < Q2 < 80, |t|<1");
TGraphErrors *sgraph1b = new TGraphErrors(5); sgraph1b->SetTitle("6.5 < Q2 < 80, |t|<1");
//
// Loop over data points and for each initialize Sartre and
// calculate the cross-section. We are using no runcard but
// set thigs for each data point
//
Sartre sartre;
EventGeneratorSettings* settings = sartre.runSettings();
//
// Table 4
//
double diff_table4, chi2_table4 = 0;
int ndof_table4 = 0;
int graphCounter, pointCounter;
for (unsigned int i=0; iSetPoint(pointCounter, fabs(h1DataDVCS_table4[i].t), h1DataDVCS_table4[i].xsection);
double error = sqrt(h1DataDVCS_table4[i].xsection_stat*h1DataDVCS_table4[i].xsection_stat + h1DataDVCS_table4[i].xsection_sys*h1DataDVCS_table4[i].xsection_sys);
graphVec[graphCounter]->SetPointError(pointCounter, 0., error);
// set up sartre
settings->setVerbose(false);
settings->setVerboseLevel(0);
settings->setNumberOfEvents(0);
settings->setTimesToShow(0);
settings->setVectorMesonId(22);
settings->setApplyPhotonFlux(false);
settings->setElectronBeamEnergy(heraElectronEnergy);
settings->setHadronBeamEnergy(heraProtonEnergy);
settings->setDipoleModelType(theDipoleModelType);
settings->setDipoleModelParameterSet(theDipoleModelParameterSet);
settings->setA(1);
settings->setCorrectForRealAmplitude(false);
settings->setCorrectSkewedness(false);
settings->setEnableNuclearBreakup(false);
settings->setQ2min(h1DataDVCS_table4[i].Q2-.1);
settings->setQ2max(h1DataDVCS_table4[i].Q2+.1);
settings->setWmin(h1DataDVCS_table4[i].W-.1);
settings->setWmax(h1DataDVCS_table4[i].W+.1);
if (settings->Wmax() > 142 || settings->Wmin()>142) continue;
settings->setMaxLambdaUsedInCorrections(maxLam);
bool ok = sartre.init();
if (!ok) {
cout << "Initialization of sartre failed. Going to next data point." << endl;
continue;
}
//
// Sartre cross-section
//
vector > limits = sartre.kinematicLimits(); // t, Q2, W
double lower[3], upper[3];
lower[0] = -h1DataDVCS_table4[i].t-0.05; // t
upper[0] = -h1DataDVCS_table4[i].t+0.05;
lower[1] = limits[1].first; // Q2
upper[1] = limits[1].second;
lower[2] = limits[2].first; // W
upper[2] = limits[2].second;
double CS = sartre.totalCrossSection(lower, upper);
CS /= upper[0]-lower[0];
CS /= limits[1].second - limits[1].first;
CS /= limits[2].second*limits[2].second - limits[2].first*limits[2].first;
sgraphVec[graphCounter]->SetPoint(pointCounter, fabs(h1DataDVCS_table4[i].t), CS);
sgraphVec[graphCounter]->SetPointError(pointCounter, 0., 0.);
//
// chi2 calculation
//
diff_table4 = CS - h1DataDVCS_table4[i].xsection;
chi2_table4 += (diff_table4*diff_table4) / (error*error);
ndof_table4++;
//
// Print out
//
cout << endl;
cout << "dsig/dt at |t|=" << fabs(h1DataDVCS_table4[i].t)
<< ", Q2=" << h1DataDVCS_table4[i].Q2
<< ", W=" << h1DataDVCS_table4[i].W << endl;
cout << "\t\tSartre: " << CS << endl; // in nb/GeV2
cout << "\t\tH1: " << h1DataDVCS_table4[i].xsection << "+/-" << h1DataDVCS_table4[i].xsection_stat << "+/-" < > limits = sartre.kinematicLimits(); // t, Q2, W
double lower[3], upper[3];
lower[0] = -1; // t
upper[0] = limits[0].second;
lower[1] = limits[1].first; // Q2
upper[1] = limits[1].second;
lower[2] = limits[2].first; // W
upper[2] = limits[2].second;
double CS = sartre.totalCrossSection(lower, upper);
// H1 cross-section sum for crosschecks
double h1cs = h1DataDVCS_table1[i].xsection;
h1cs *= h1DataDVCS_table1[i].Q2max - h1DataDVCS_table1[i].Q2min; // Q2
h1cs *= h1DataDVCS_table1[i].Wmax - h1DataDVCS_table1[i].Wmin; // W
if (i < 4) {
sigma_a += CS;
h1_sigma_a += h1cs;
}
else {
sigma_b += CS;
h1_sigma_b += h1cs;
}
limits = sartre.kinematicLimits();
CS /= limits[1].second - limits[1].first; // Q2
CS /= limits[2].second*limits[2].second - limits[2].first*limits[2].first; // W2
// fill graphs
if (i < 4) {
sgraph1a->SetPoint(i, fabs(h1DataDVCS_table1[i].Q2), CS);
sgraph1a->SetPointError(i , 0., 0.);
}
else {
sgraph1b->SetPoint(i-4, fabs(h1DataDVCS_table1[i].W), CS);
sgraph1b->SetPointError(i-4 , 0., 0.);
}
//
// chi2 calculation
//
diff_table1 = CS - h1DataDVCS_table1[i].xsection;
chi2_table1 += (diff_table1*diff_table1) / (error*error);
chi2_global += (diff_table1*diff_table1) / (error*error);
ndof_table1++;
ndof_global++;
//
// Print out
//
cout << endl;
cout << "\t\tSartre: " << CS << endl;
cout << "\t\tH1: " << h1DataDVCS_table1[i].xsection << "+/-" << h1DataDVCS_table1[i].xsection_stat << "+/-" <Write(graphName[i].c_str());
sgraphVec[i]->Write(sgraphName[i].c_str());
}
// Table 1
graph1a->Write("graph1a");
graph1b->Write("graph1b");
sgraph1a->Write("sgraph1a");
sgraph1b->Write("sgraph1b");
hfile->Close();
}
void compareJpsi(string& rootfile)
{
//
// J/Psi HERA/ZEUS data
//
double heraProtonEnergy = 920;
double heraElectronEnergy = 27.5;
vector zeusDataJpsi;
loadJpsiData(zeusDataJpsi);
double chi2=0;
double ndof=0;
// HERA data - table 1
TGraphAsymmErrors *graph1H = new TGraphAsymmErrors(3); graph1H->SetTitle("0.15 < Q2 < 0.8");
TGraphAsymmErrors *graph2H = new TGraphAsymmErrors(6); graph2H->SetTitle("2 < Q2 < 5");
TGraphAsymmErrors *graph3H = new TGraphAsymmErrors(6); graph3H->SetTitle("5 < Q2 < 10");
TGraphAsymmErrors *graph4H = new TGraphAsymmErrors(6); graph4H->SetTitle("10 < Q2 < 100");
// Sartre results - table 1
TGraphAsymmErrors *graph1S = new TGraphAsymmErrors(3); graph1S->SetTitle("0.15 < Q2 < 0.8");
TGraphAsymmErrors *graph2S = new TGraphAsymmErrors(6); graph2S->SetTitle("2 < Q2 < 5");
TGraphAsymmErrors *graph3S = new TGraphAsymmErrors(6); graph3S->SetTitle("5 < Q2 < 10");
TGraphAsymmErrors *graph4S = new TGraphAsymmErrors(6); graph4S->SetTitle("10 < Q2 < 100");
// HERA data - table 2
TGraphAsymmErrors *graph5H = new TGraphAsymmErrors(8); graph5H->SetTitle("30 < W < 220");
TGraphAsymmErrors *graph6H = new TGraphAsymmErrors(7); graph6H->SetTitle("45 < W < 160");
// Sartre results - table 2
TGraphAsymmErrors *graph5S = new TGraphAsymmErrors(8); graph5S->SetTitle("30 < W < 220");
TGraphAsymmErrors *graph6S = new TGraphAsymmErrors(7); graph6S->SetTitle("45 < W < 160");
//
// Loop over data points and for each initialize Sartre and
// calculate the cross-section. We are using no runcard but
// set things for each data point
//
Sartre sartre;
EventGeneratorSettings* settings = sartre.runSettings();
double diff;
for (unsigned int i=0; isetVerbose(true);
settings->setVerboseLevel(2);
settings->setNumberOfEvents(0);
settings->setTimesToShow(0);
settings->setQ2min(zeusDataJpsi[i].Q2min);
settings->setQ2max(zeusDataJpsi[i].Q2max);
settings->setWmin(zeusDataJpsi[i].Wmin);
settings->setWmax(zeusDataJpsi[i].Wmax);
settings->setVectorMesonId(443);
settings->setApplyPhotonFlux(true);
settings->setElectronBeamEnergy(heraElectronEnergy);
settings->setHadronBeamEnergy(heraProtonEnergy);
settings->setDipoleModelType(theDipoleModelType);
settings->setDipoleModelParameterSet(theDipoleModelParameterSet);
settings->setA(1);
settings->setCorrectForRealAmplitude(true);
settings->setCorrectSkewedness(true);
settings->setEnableNuclearBreakup(false);
settings->setMaxLambdaUsedInCorrections(maxLam);
settings->setTableSetType(theTableSetType);
bool ok = sartre.init();
if (!ok) {
cout << "Initialization of sartre failed. Going to next data point." << endl;
continue;
}
settings->list();
cout << "Total cross-sections for Q2=[" << zeusDataJpsi[i].Q2min << ", " << zeusDataJpsi[i].Q2max << "], ";
cout << "W=[" << zeusDataJpsi[i].Wmin << ", " << zeusDataJpsi[i].Wmax<< "], |t| < 1" << endl;
cout << "\t\tSartre: " << sartre.totalCrossSection()*1000 << endl; // in pb
cout << "\t\tZEUS: " << zeusDataJpsi[i].xsection << "+/-" << zeusDataJpsi[i].xsection_stat << "+" <= 3 && i <= 8) {
graphS = graph2S;
graphH = graph2H;
offset = 3;
}
else if (i >= 9 && i <= 14) {
graphS = graph3S;
graphH = graph3H;
offset = 9;
}
else if (i >= 15 && i <= 20) {
graphS = graph4S;
graphH = graph4H;
offset = 15;
}
else if (i >= 21 && i <= 28) {
graphSS = graph5S;
graphHH = graph5H;
offset = 21;
}
else if (i >= 29 && i <= 35) {
graphSS = graph6S;
graphHH = graph6H;
offset = 29;
}
double meanW, meanQ2, exl, exu, eyl, eyu;
eyl = sqrt(zeusDataJpsi[i].xsection_stat*zeusDataJpsi[i].xsection_stat + zeusDataJpsi[i].xsection_sys_lo*zeusDataJpsi[i].xsection_sys_lo);
eyu = sqrt(zeusDataJpsi[i].xsection_stat*zeusDataJpsi[i].xsection_stat + zeusDataJpsi[i].xsection_sys_up*zeusDataJpsi[i].xsection_sys_up);
if (graphS && graphH) {
meanW = (zeusDataJpsi[i].Wmax+zeusDataJpsi[i].Wmin)/2;
graphS->SetPoint(i-offset, meanW, 1000*sartre.totalCrossSection());
graphS->SetPointError(i-offset, 0, 0, 0, 0);
exl = meanW-zeusDataJpsi[i].Wmin;
exu = zeusDataJpsi[i].Wmax-meanW;
graphH->SetPoint(i-offset, meanW, zeusDataJpsi[i].xsection);
graphH->SetPointError(i-offset, exl, exu, eyl, eyu);
}
else if (graphSS && graphHH) {
meanQ2 = (zeusDataJpsi[i].Q2max+zeusDataJpsi[i].Q2min)/2;
graphSS->SetPoint(i-offset, meanQ2, 1000*sartre.totalCrossSection());
graphSS->SetPointError(i-offset, 0, 0, 0, 0);
exl = meanQ2-zeusDataJpsi[i].Q2min;
exu = zeusDataJpsi[i].Q2max-meanQ2;
graphHH->SetPoint(i-offset, meanQ2, zeusDataJpsi[i].xsection);
graphHH->SetPointError(i-offset, exl, exu, eyl, eyu);
}
diff = 1000*sartre.totalCrossSection() - zeusDataJpsi[i].xsection;
if ( zeusDataJpsi[i].xsection > 1000*sartre.totalCrossSection() ) {
chi2 += (diff*diff) / (eyl*eyl);
chi2_global+= (diff*diff) / (eyl*eyl);
}
else {
chi2 += (diff*diff) / (eyu*eyu);
chi2_global += (diff*diff) / (eyu*eyu);
}
ndof++;
}
cout << "-------------------------------------------------" << endl;
cout << "Comparison HERA-Sartre:" << endl;
cout << "chi2 = " << chi2 << endl;
cout << "ndof = " << ndof << endl;
cout << "chi2/ndof = " << chi2/ndof << endl;
cout << "-------------------------------------------------" << endl;
//
// ROOT file
//
TFile *hfile = 0;
hfile = new TFile(rootfile.c_str(),"RECREATE");
cout << "ROOT file is '" << rootfile.c_str() << "'." << endl;
graph1H->Write("graph1H");
graph1S->Write("graph1S");
graph2H->Write("graph2H");
graph2S->Write("graph2S");
graph3H->Write("graph3H");
graph3S->Write("graph3S");
graph4H->Write("graph4H");
graph4S->Write("graph4S");
graph5H->Write("graph5H");
graph5S->Write("graph5S");
graph6H->Write("graph6H");
graph6S->Write("graph6S");
hfile->Close();
return;
}
void loadJpsiData(vector& vec)
{
//
// J/Psi HERA/ZEUS data
//
// Data for e p -> e' p' J/psi
// Journal reference: Nucl.Phys.B695:3-37,2004
// arXiv:hep-ex/0404008v1
//
// All data are for |t|<1 GeV2
// Cross section in pb
//
vector zeusDataJpsi(21+15);
unsigned int i = 0;
// table 1 (21 entries)
zeusDataJpsi[i].Q2min=0.15; zeusDataJpsi[i].Q2max=0.8; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=65; zeusDataJpsi[i].xsection=217; zeusDataJpsi[i].xsection_stat=53; zeusDataJpsi[i].xsection_sys_up=12; zeusDataJpsi[i].xsection_sys_lo=19; i++;
zeusDataJpsi[i].Q2min=0.15; zeusDataJpsi[i].Q2max=0.8; zeusDataJpsi[i].Wmin=65; zeusDataJpsi[i].Wmax=105; zeusDataJpsi[i].xsection=257; zeusDataJpsi[i].xsection_stat=46; zeusDataJpsi[i].xsection_sys_up=18; zeusDataJpsi[i].xsection_sys_lo=17; i++;
zeusDataJpsi[i].Q2min=0.15; zeusDataJpsi[i].Q2max=0.8; zeusDataJpsi[i].Wmin=105; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=498; zeusDataJpsi[i].xsection_stat=89; zeusDataJpsi[i].xsection_sys_up=37; zeusDataJpsi[i].xsection_sys_lo=38; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=45; zeusDataJpsi[i].xsection=41.5; zeusDataJpsi[i].xsection_stat=8.4; zeusDataJpsi[i].xsection_sys_up=5.6; zeusDataJpsi[i].xsection_sys_lo=6.6; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=70; zeusDataJpsi[i].xsection=48.8; zeusDataJpsi[i].xsection_stat=5.2; zeusDataJpsi[i].xsection_sys_up=3.1; zeusDataJpsi[i].xsection_sys_lo=3.9; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=70; zeusDataJpsi[i].Wmax=90; zeusDataJpsi[i].xsection=36.4; zeusDataJpsi[i].xsection_stat=4.1; zeusDataJpsi[i].xsection_sys_up=10.5; zeusDataJpsi[i].xsection_sys_lo=3.0; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=90; zeusDataJpsi[i].Wmax=112; zeusDataJpsi[i].xsection=35.4; zeusDataJpsi[i].xsection_stat=4.0; zeusDataJpsi[i].xsection_sys_up=3.0; zeusDataJpsi[i].xsection_sys_lo=4.5; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=112; zeusDataJpsi[i].Wmax=145; zeusDataJpsi[i].xsection=44.7; zeusDataJpsi[i].xsection_stat=5.0; zeusDataJpsi[i].xsection_sys_up=9.0; zeusDataJpsi[i].xsection_sys_lo=4.3; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=145; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=76.5; zeusDataJpsi[i].xsection_stat=10.3; zeusDataJpsi[i].xsection_sys_up=11.5; zeusDataJpsi[i].xsection_sys_lo=5.1; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=50; zeusDataJpsi[i].xsection=19.6; zeusDataJpsi[i].xsection_stat=4.1; zeusDataJpsi[i].xsection_sys_up=3.9; zeusDataJpsi[i].xsection_sys_lo=1.9; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=50; zeusDataJpsi[i].Wmax=74; zeusDataJpsi[i].xsection=19.3; zeusDataJpsi[i].xsection_stat=2.2; zeusDataJpsi[i].xsection_sys_up=2.9; zeusDataJpsi[i].xsection_sys_lo=1.3; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=74; zeusDataJpsi[i].Wmax=96; zeusDataJpsi[i].xsection=15.6; zeusDataJpsi[i].xsection_stat=1.8; zeusDataJpsi[i].xsection_sys_up=1.6; zeusDataJpsi[i].xsection_sys_lo=1.4; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=96; zeusDataJpsi[i].Wmax=120; zeusDataJpsi[i].xsection=13.5; zeusDataJpsi[i].xsection_stat=1.7; zeusDataJpsi[i].xsection_sys_up=1.1; zeusDataJpsi[i].xsection_sys_lo=0.7; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=120; zeusDataJpsi[i].Wmax=150; zeusDataJpsi[i].xsection=14.9; zeusDataJpsi[i].xsection_stat=1.9; zeusDataJpsi[i].xsection_sys_up=1.1; zeusDataJpsi[i].xsection_sys_lo=1.3; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=150; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=27.9; zeusDataJpsi[i].xsection_stat=4.1; zeusDataJpsi[i].xsection_sys_up=4.5; zeusDataJpsi[i].xsection_sys_lo=1.4; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=55; zeusDataJpsi[i].xsection=10.9; zeusDataJpsi[i].xsection_stat=3.1; zeusDataJpsi[i].xsection_sys_up=0.8; zeusDataJpsi[i].xsection_sys_lo=1.0; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=55; zeusDataJpsi[i].Wmax=78; zeusDataJpsi[i].xsection=8.4; zeusDataJpsi[i].xsection_stat=1.2; zeusDataJpsi[i].xsection_sys_up=1.4; zeusDataJpsi[i].xsection_sys_lo=0.4; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=78; zeusDataJpsi[i].Wmax=100; zeusDataJpsi[i].xsection=8.6; zeusDataJpsi[i].xsection_stat=1.1; zeusDataJpsi[i].xsection_sys_up=0.9; zeusDataJpsi[i].xsection_sys_lo=1.4; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=100; zeusDataJpsi[i].Wmax=124; zeusDataJpsi[i].xsection=8.4; zeusDataJpsi[i].xsection_stat=1.1; zeusDataJpsi[i].xsection_sys_up=0.4; zeusDataJpsi[i].xsection_sys_lo=1.2; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=124; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=10.8; zeusDataJpsi[i].xsection_stat=1.4; zeusDataJpsi[i].xsection_sys_up=2.1; zeusDataJpsi[i].xsection_sys_lo=0.8; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=160; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=25.1; zeusDataJpsi[i].xsection_stat=3.8; zeusDataJpsi[i].xsection_sys_up=1.7; zeusDataJpsi[i].xsection_sys_lo=1.2; i++;
// table 2 (15 entries)
zeusDataJpsi[i].Q2min=0.15; zeusDataJpsi[i].Q2max=0.8; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=954; zeusDataJpsi[i].xsection_stat=108; zeusDataJpsi[i].xsection_sys_up=63; zeusDataJpsi[i].xsection_sys_lo=74; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=3.2; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=150; zeusDataJpsi[i].xsection_stat=14; zeusDataJpsi[i].xsection_sys_up=53; zeusDataJpsi[i].xsection_sys_lo=8; i++;
zeusDataJpsi[i].Q2min=3.2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=132; zeusDataJpsi[i].xsection_stat=12; zeusDataJpsi[i].xsection_sys_up=8; zeusDataJpsi[i].xsection_sys_lo=17; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=7; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=59.9; zeusDataJpsi[i].xsection_stat=6.1; zeusDataJpsi[i].xsection_sys_up=5.5; zeusDataJpsi[i].xsection_sys_lo=3.6; i++;
zeusDataJpsi[i].Q2min=7; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=42.6; zeusDataJpsi[i].xsection_stat=4.3; zeusDataJpsi[i].xsection_sys_up=4.7; zeusDataJpsi[i].xsection_sys_lo=5.0; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=15; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=36.7; zeusDataJpsi[i].xsection_stat=4.0; zeusDataJpsi[i].xsection_sys_up=1.4; zeusDataJpsi[i].xsection_sys_lo=2.9; i++;
zeusDataJpsi[i].Q2min=15; zeusDataJpsi[i].Q2max=40; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=29.3; zeusDataJpsi[i].xsection_stat=3.7; zeusDataJpsi[i].xsection_sys_up=2.0; zeusDataJpsi[i].xsection_sys_lo=4.7; i++;
zeusDataJpsi[i].Q2min=40; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=30; zeusDataJpsi[i].Wmax=220; zeusDataJpsi[i].xsection=4.5; zeusDataJpsi[i].xsection_stat=1.5; zeusDataJpsi[i].xsection_sys_up=0.5; zeusDataJpsi[i].xsection_sys_lo=1.1; i++;
zeusDataJpsi[i].Q2min=2; zeusDataJpsi[i].Q2max=3.2; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=96; zeusDataJpsi[i].xsection_stat=11; zeusDataJpsi[i].xsection_sys_up=5; zeusDataJpsi[i].xsection_sys_lo=14; i++;
zeusDataJpsi[i].Q2min=3.2; zeusDataJpsi[i].Q2max=5; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=91.6; zeusDataJpsi[i].xsection_stat=8.9; zeusDataJpsi[i].xsection_sys_up=12.2; zeusDataJpsi[i].xsection_sys_lo=6.6; i++;
zeusDataJpsi[i].Q2min=5; zeusDataJpsi[i].Q2max=7; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=48.7; zeusDataJpsi[i].xsection_stat=5.2; zeusDataJpsi[i].xsection_sys_up=1.2; zeusDataJpsi[i].xsection_sys_lo=2.5; i++;
zeusDataJpsi[i].Q2min=7; zeusDataJpsi[i].Q2max=10; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=32.5; zeusDataJpsi[i].xsection_stat=3.4; zeusDataJpsi[i].xsection_sys_up=4.4; zeusDataJpsi[i].xsection_sys_lo=2.4; i++;
zeusDataJpsi[i].Q2min=10; zeusDataJpsi[i].Q2max=15; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=24.1; zeusDataJpsi[i].xsection_stat=2.8; zeusDataJpsi[i].xsection_sys_up=1.2; zeusDataJpsi[i].xsection_sys_lo=1.6; i++;
zeusDataJpsi[i].Q2min=15; zeusDataJpsi[i].Q2max=40; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=18.4; zeusDataJpsi[i].xsection_stat=2.4; zeusDataJpsi[i].xsection_sys_up=0.9; zeusDataJpsi[i].xsection_sys_lo=1.4; i++;
zeusDataJpsi[i].Q2min=40; zeusDataJpsi[i].Q2max=100; zeusDataJpsi[i].Wmin=45; zeusDataJpsi[i].Wmax=160; zeusDataJpsi[i].xsection=2.2; zeusDataJpsi[i].xsection_stat=0.9; zeusDataJpsi[i].xsection_sys_up=0.4; zeusDataJpsi[i].xsection_sys_lo=0.6;
vec = zeusDataJpsi;
}
void loadDVCSData(vector& vec1, vector& vec4)
{
//
// DVCS HERA/ZEUS data
//
// Origin:
// Journal reference: Phys.Lett.B659:796-806,2008
// Report number: DESY-07-142
// Cite as: arXiv:0709.4114v1 [hep-ex]
vector h1DataDVCS_table1(4+5);
vector h1DataDVCS_table4(12+12);
unsigned int i = 0;
// Table 1a (4 points versus Q2)
h1DataDVCS_table1[i].Wmin=30; h1DataDVCS_table1[i].Wmax=140; h1DataDVCS_table1[i].W=82; // W is ref only
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=6.5; h1DataDVCS_table1[i].Q2max=11; h1DataDVCS_table1[i].Q2=8.75; // Q2 is ref value only
h1DataDVCS_table1[i].xsection=3.59; h1DataDVCS_table1[i].xsection_stat=.21; h1DataDVCS_table1[i].xsection_sys=.41;
i++;
h1DataDVCS_table1[i].Wmin=30; h1DataDVCS_table1[i].Wmax=140;; h1DataDVCS_table1[i].W=82; // W is ref only
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=11; h1DataDVCS_table1[i].Q2max=20; h1DataDVCS_table1[i].Q2=15.5;
h1DataDVCS_table1[i].xsection=1.38; h1DataDVCS_table1[i].xsection_stat=.1; h1DataDVCS_table1[i].xsection_sys=.21;
i++;
h1DataDVCS_table1[i].Wmin=30; h1DataDVCS_table1[i].Wmax=140;; h1DataDVCS_table1[i].W=82; // W is ref only
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=20; h1DataDVCS_table1[i].Q2max=30; h1DataDVCS_table1[i].Q2=25;
h1DataDVCS_table1[i].xsection=.58; h1DataDVCS_table1[i].xsection_stat=.09; h1DataDVCS_table1[i].xsection_sys=.09;
i++;
h1DataDVCS_table1[i].Wmin=30; h1DataDVCS_table1[i].Wmax=140;; h1DataDVCS_table1[i].W=82; // W is ref only
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=30; h1DataDVCS_table1[i].Q2max=80; h1DataDVCS_table1[i].Q2=55;
h1DataDVCS_table1[i].xsection=.13; h1DataDVCS_table1[i].xsection_stat=.03; h1DataDVCS_table1[i].xsection_sys=.04;
i++;
// Table 1b (5 points versus W)
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=6.5; h1DataDVCS_table1[i].Q2max=80; h1DataDVCS_table1[i].Q2=8; // Q2 is ref only
h1DataDVCS_table1[i].Wmin=30; h1DataDVCS_table1[i].Wmax=60; h1DataDVCS_table1[i].W=45;
h1DataDVCS_table1[i].xsection=2.91; h1DataDVCS_table1[i].xsection_stat=.2; h1DataDVCS_table1[i].xsection_sys=.25;
i++;
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=6.5; h1DataDVCS_table1[i].Q2max=80; h1DataDVCS_table1[i].Q2=8; // Q2 is ref only
h1DataDVCS_table1[i].Wmin=60; h1DataDVCS_table1[i].Wmax=80; h1DataDVCS_table1[i].W=70;// W is ref value only
h1DataDVCS_table1[i].xsection=3.96; h1DataDVCS_table1[i].xsection_stat=.32; h1DataDVCS_table1[i].xsection_sys=.37;
i++;
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=6.5; h1DataDVCS_table1[i].Q2max=80; h1DataDVCS_table1[i].Q2=8; // Q2 is ref only
h1DataDVCS_table1[i].Wmin=80; h1DataDVCS_table1[i].Wmax=100; h1DataDVCS_table1[i].W=90;
h1DataDVCS_table1[i].xsection= 4.78; h1DataDVCS_table1[i].xsection_stat=.41; h1DataDVCS_table1[i].xsection_sys=.57;
i++;
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=6.5; h1DataDVCS_table1[i].Q2max=80; h1DataDVCS_table1[i].Q2=8; // Q2 is ref only
h1DataDVCS_table1[i].Wmin=100; h1DataDVCS_table1[i].Wmax=120; h1DataDVCS_table1[i].W=110;
h1DataDVCS_table1[i].xsection= 5.55; h1DataDVCS_table1[i].xsection_stat=.57; h1DataDVCS_table1[i].xsection_sys=.88;
i++;
h1DataDVCS_table1[i].tmin=0; h1DataDVCS_table1[i].tmax=1;
h1DataDVCS_table1[i].Q2min=6.5; h1DataDVCS_table1[i].Q2max=80; h1DataDVCS_table1[i].Q2=8;
h1DataDVCS_table1[i].Wmin=120; h1DataDVCS_table1[i].Wmax=140; h1DataDVCS_table1[i].W=130;
h1DataDVCS_table1[i].xsection= 6.56; h1DataDVCS_table1[i].xsection_stat=1.17; h1DataDVCS_table1[i].xsection_sys=1.77;
i++;
vec1 = h1DataDVCS_table1;
/*
// Table 3a
// PROBLEM: H1 gives W = 82 w/o specifying the bin range
// and the same for Q=10. The call this reference value
// which is pretty meaningless. The bins here are a guess
// and are apparently not working.
// We leave those out for the comparison.
//
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=11; h1DataDVCS[i].Q2=8.;
h1DataDVCS[i].tmin=0; h1DataDVCS[i].tmax=0.2; h1DataDVCS[i].t=0.1;
h1DataDVCS[i].xsection=13.1; h1DataDVCS[i].xsection_stat=1.10; h1DataDVCS[i].xsection_sys=1.85;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=11; h1DataDVCS[i].Q2=8.;
h1DataDVCS[i].tmin=0.2; h1DataDVCS[i].tmax=0.4; h1DataDVCS[i].t=0.3;
h1DataDVCS[i].xsection=4.69; h1DataDVCS[i].xsection_stat=0.45; h1DataDVCS[i].xsection_sys=0.55;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=11; h1DataDVCS[i].Q2=8.;
h1DataDVCS[i].tmin=0.4; h1DataDVCS[i].tmax=0.6; h1DataDVCS[i].t=0.5;
h1DataDVCS[i].xsection=1.37; h1DataDVCS[i].xsection_stat=0.21; h1DataDVCS[i].xsection_sys=0.23;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=11; h1DataDVCS[i].Q2=8.;
h1DataDVCS[i].tmin=0.6; h1DataDVCS[i].tmax=1; h1DataDVCS[i].t=0.8;
h1DataDVCS[i].xsection=0.19; h1DataDVCS[i].xsection_stat=0.04; h1DataDVCS[i].xsection_sys=0.06;
i++;
//-
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=11; h1DataDVCS[i].Q2max=20; h1DataDVCS[i].Q2=15.5;
h1DataDVCS[i].tmin=0; h1DataDVCS[i].tmax=0.2; h1DataDVCS[i].t=0.1;
h1DataDVCS[i].xsection=4.37; h1DataDVCS[i].xsection_stat=0.47; h1DataDVCS[i].xsection_sys=0.86;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=11; h1DataDVCS[i].Q2max=20; h1DataDVCS[i].Q2=15.5;
h1DataDVCS[i].tmin=0.2; h1DataDVCS[i].tmax=0.4; h1DataDVCS[i].t=0.3;
h1DataDVCS[i].xsection=1.02; h1DataDVCS[i].xsection_stat=0.16; h1DataDVCS[i].xsection_sys=0.18;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=11; h1DataDVCS[i].Q2max=20; h1DataDVCS[i].Q2=15.5;
h1DataDVCS[i].tmin=0.4; h1DataDVCS[i].tmax=0.6; h1DataDVCS[i].t=0.5;
h1DataDVCS[i].xsection=0.49; h1DataDVCS[i].xsection_stat=0.08; h1DataDVCS[i].xsection_sys=0.08;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=11; h1DataDVCS[i].Q2max=20; h1DataDVCS[i].Q2=15.5;
h1DataDVCS[i].tmin=0.6; h1DataDVCS[i].tmax=1; h1DataDVCS[i].t=0.8;
h1DataDVCS[i].xsection=0.12; h1DataDVCS[i].xsection_stat=0.02; h1DataDVCS[i].xsection_sys=0.02;
i++;
//-
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=20; h1DataDVCS[i].Q2max=80; h1DataDVCS[i].Q2=25;
h1DataDVCS[i].tmin=0; h1DataDVCS[i].tmax=0.2; h1DataDVCS[i].t=0.1;
h1DataDVCS[i].xsection=1.41; h1DataDVCS[i].xsection_stat=0.40; h1DataDVCS[i].xsection_sys=0.43;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=20; h1DataDVCS[i].Q2max=80; h1DataDVCS[i].Q2=25;
h1DataDVCS[i].tmin=0.2; h1DataDVCS[i].tmax=0.4; h1DataDVCS[i].t=0.3;
h1DataDVCS[i].xsection=0.71; h1DataDVCS[i].xsection_stat=0.16; h1DataDVCS[i].xsection_sys=0.08;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=20; h1DataDVCS[i].Q2max=80; h1DataDVCS[i].Q2=25;
h1DataDVCS[i].tmin=0.4; h1DataDVCS[i].tmax=0.6; h1DataDVCS[i].t=0.5;
h1DataDVCS[i].xsection=0.28; h1DataDVCS[i].xsection_stat=0.07; h1DataDVCS[i].xsection_sys=0.04;
i++;
h1DataDVCS[i].W=82; h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=140;
h1DataDVCS[i].Q2min=20; h1DataDVCS[i].Q2max=80; h1DataDVCS[i].Q2=25;
h1DataDVCS[i].tmin=0.6; h1DataDVCS[i].tmax=1; h1DataDVCS[i].t=0.8;
h1DataDVCS[i].xsection=0.04; h1DataDVCS[i].xsection_stat=0.01; h1DataDVCS[i].xsection_sys=0.02;
i++;
// Table 3b
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=50; h1DataDVCS[i].W=40;
h1DataDVCS[i].tmin=0; h1DataDVCS[i].tmax=0.2; h1DataDVCS[i].t=0.1;
h1DataDVCS[i].xsection=4.99; h1DataDVCS[i].xsection_stat=0.66; h1DataDVCS[i].xsection_sys=0.54;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=50; h1DataDVCS[i].W=40;
h1DataDVCS[i].tmin=0.2; h1DataDVCS[i].tmax=0.4; h1DataDVCS[i].t=0.3;
h1DataDVCS[i].xsection=1.45; h1DataDVCS[i].xsection_stat=0.29; h1DataDVCS[i].xsection_sys=0.18;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=50; h1DataDVCS[i].W=40;
h1DataDVCS[i].tmin=0.4; h1DataDVCS[i].tmax=0.6; h1DataDVCS[i].t=0.5;
h1DataDVCS[i].xsection=0.49; h1DataDVCS[i].xsection_stat=0.14; h1DataDVCS[i].xsection_sys=0.08;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=30; h1DataDVCS[i].Wmax=50; h1DataDVCS[i].W=40;
h1DataDVCS[i].tmin=0.6; h1DataDVCS[i].tmax=1; h1DataDVCS[i].t=0.8;
h1DataDVCS[i].xsection=0.12; h1DataDVCS[i].xsection_stat=0.03; h1DataDVCS[i].xsection_sys=0.03;
i++;
//-
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=50; h1DataDVCS[i].Wmax=85; h1DataDVCS[i].W=70;
h1DataDVCS[i].tmin=0; h1DataDVCS[i].tmax=0.2; h1DataDVCS[i].t=0.1;
h1DataDVCS[i].xsection=7.78; h1DataDVCS[i].xsection_stat=0.69; h1DataDVCS[i].xsection_sys=0.87;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=50; h1DataDVCS[i].Wmax=85; h1DataDVCS[i].W=70;
h1DataDVCS[i].tmin=0.2; h1DataDVCS[i].tmax=0.4; h1DataDVCS[i].t=0.3;
h1DataDVCS[i].xsection=2.74; h1DataDVCS[i].xsection_stat=0.31; h1DataDVCS[i].xsection_sys=0.3;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=50; h1DataDVCS[i].Wmax=85; h1DataDVCS[i].W=70;
h1DataDVCS[i].tmin=0.4; h1DataDVCS[i].tmax=0.6; h1DataDVCS[i].t=0.5;
h1DataDVCS[i].xsection=0.81; h1DataDVCS[i].xsection_stat=0.14; h1DataDVCS[i].xsection_sys=0.11;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=50; h1DataDVCS[i].Wmax=85; h1DataDVCS[i].W=70;
h1DataDVCS[i].tmin=0.6; h1DataDVCS[i].tmax=1; h1DataDVCS[i].t=0.8;
h1DataDVCS[i].xsection=0.19; h1DataDVCS[i].xsection_stat=0.03; h1DataDVCS[i].xsection_sys=0.03;
i++;
//-
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=85; h1DataDVCS[i].Wmax=140; h1DataDVCS[i].W=100;
h1DataDVCS[i].tmin=0; h1DataDVCS[i].tmax=0.2; h1DataDVCS[i].t=0.1;
h1DataDVCS[i].xsection=10.9; h1DataDVCS[i].xsection_stat=1.14; h1DataDVCS[i].xsection_sys=2.36;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=85; h1DataDVCS[i].Wmax=140; h1DataDVCS[i].W=100;
h1DataDVCS[i].tmin=0.2; h1DataDVCS[i].tmax=0.4; h1DataDVCS[i].t=0.3;
h1DataDVCS[i].xsection=3.47; h1DataDVCS[i].xsection_stat=0.42; h1DataDVCS[i].xsection_sys=0.53;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=85; h1DataDVCS[i].Wmax=140; h1DataDVCS[i].W=100;
h1DataDVCS[i].tmin=0.4; h1DataDVCS[i].tmax=0.6; h1DataDVCS[i].t=0.5;
h1DataDVCS[i].xsection=1.49; h1DataDVCS[i].xsection_stat=0.21; h1DataDVCS[i].xsection_sys=0.24;
i++;
h1DataDVCS[i].Q2=10; h1DataDVCS[i].Q2min=6.5; h1DataDVCS[i].Q2max=80;
h1DataDVCS[i].Wmin=85; h1DataDVCS[i].Wmax=140; h1DataDVCS[i].W=100;
h1DataDVCS[i].tmin=0.6; h1DataDVCS[i].tmax=1; h1DataDVCS[i].t=0.8;
h1DataDVCS[i].xsection=0.19; h1DataDVCS[i].xsection_stat=0.04; h1DataDVCS[i].xsection_sys=0.06;
i++;
*/
// Table 4a
i=0;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0; h1DataDVCS_table4[i].tmax=0.2; h1DataDVCS_table4[i].t=0.1;
h1DataDVCS_table4[i].xsection=8.10; h1DataDVCS_table4[i].xsection_stat=1.22; h1DataDVCS_table4[i].xsection_sys=0.82;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0.2; h1DataDVCS_table4[i].tmax=0.4; h1DataDVCS_table4[i].t=0.3;
h1DataDVCS_table4[i].xsection=2.30; h1DataDVCS_table4[i].xsection_stat=0.54; h1DataDVCS_table4[i].xsection_sys=0.28;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0.4; h1DataDVCS_table4[i].tmax=0.6; h1DataDVCS_table4[i].t=0.5;
h1DataDVCS_table4[i].xsection=0.45; h1DataDVCS_table4[i].xsection_stat=0.22; h1DataDVCS_table4[i].xsection_sys=0.10;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0.6; h1DataDVCS_table4[i].tmax=1; h1DataDVCS_table4[i].t=0.8;
h1DataDVCS_table4[i].xsection=0.16; h1DataDVCS_table4[i].xsection_stat=0.06; h1DataDVCS_table4[i].xsection_sys=0.03;
i++;
//-
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0; h1DataDVCS_table4[i].tmax=0.2; h1DataDVCS_table4[i].t=0.1;
h1DataDVCS_table4[i].xsection=10.0; h1DataDVCS_table4[i].xsection_stat=1.30; h1DataDVCS_table4[i].xsection_sys=1.27;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0.2; h1DataDVCS_table4[i].tmax=0.4; h1DataDVCS_table4[i].t=0.3;
h1DataDVCS_table4[i].xsection=4.35; h1DataDVCS_table4[i].xsection_stat=0.63; h1DataDVCS_table4[i].xsection_sys=0.46;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0.4; h1DataDVCS_table4[i].tmax=0.6; h1DataDVCS_table4[i].t=0.5;
h1DataDVCS_table4[i].xsection=1.08; h1DataDVCS_table4[i].xsection_stat=0.27; h1DataDVCS_table4[i].xsection_sys=0.17;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0.6; h1DataDVCS_table4[i].tmax=1; h1DataDVCS_table4[i].t=0.8;
h1DataDVCS_table4[i].xsection=0.13; h1DataDVCS_table4[i].xsection_stat=0.06; h1DataDVCS_table4[i].xsection_sys=0.04;
i++;
//-
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0; h1DataDVCS_table4[i].tmax=0.2; h1DataDVCS_table4[i].t=0.1;
h1DataDVCS_table4[i].xsection=16.0; h1DataDVCS_table4[i].xsection_stat=2.11; h1DataDVCS_table4[i].xsection_sys=2.74;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0.2; h1DataDVCS_table4[i].tmax=0.4; h1DataDVCS_table4[i].t=0.3;
h1DataDVCS_table4[i].xsection=5.45; h1DataDVCS_table4[i].xsection_stat=0.80; h1DataDVCS_table4[i].xsection_sys=0.73;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0.4; h1DataDVCS_table4[i].tmax=0.6; h1DataDVCS_table4[i].t=0.5;
h1DataDVCS_table4[i].xsection=1.96; h1DataDVCS_table4[i].xsection_stat=0.41; h1DataDVCS_table4[i].xsection_sys=0.35;
i++;
h1DataDVCS_table4[i].Q2min=6.5; h1DataDVCS_table4[i].Q2max=11; h1DataDVCS_table4[i].Q2=8.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0.6; h1DataDVCS_table4[i].tmax=1; h1DataDVCS_table4[i].t=0.8;
h1DataDVCS_table4[i].xsection=0.21; h1DataDVCS_table4[i].xsection_stat=0.09; h1DataDVCS_table4[i].xsection_sys=0.08;
i++;
// Table 4b
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0; h1DataDVCS_table4[i].tmax=0.2; h1DataDVCS_table4[i].t=0.1;
h1DataDVCS_table4[i].xsection=1.06; h1DataDVCS_table4[i].xsection_stat=0.28; h1DataDVCS_table4[i].xsection_sys=0.28;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0.2; h1DataDVCS_table4[i].tmax=0.4; h1DataDVCS_table4[i].t=0.3;
h1DataDVCS_table4[i].xsection=0.33; h1DataDVCS_table4[i].xsection_stat=0.07; h1DataDVCS_table4[i].xsection_sys=0.07;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0.4; h1DataDVCS_table4[i].tmax=0.6; h1DataDVCS_table4[i].t=0.5;
h1DataDVCS_table4[i].xsection=0.22; h1DataDVCS_table4[i].xsection_stat=0.06; h1DataDVCS_table4[i].xsection_sys=0.06;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=30; h1DataDVCS_table4[i].Wmax=50; h1DataDVCS_table4[i].W=40;
h1DataDVCS_table4[i].tmin=0.6; h1DataDVCS_table4[i].tmax=1; h1DataDVCS_table4[i].t=0.8;
h1DataDVCS_table4[i].xsection=0.04; h1DataDVCS_table4[i].xsection_stat=0.01; h1DataDVCS_table4[i].xsection_sys=0.01;
i++;
//-
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0; h1DataDVCS_table4[i].tmax=0.2; h1DataDVCS_table4[i].t=0.1;
h1DataDVCS_table4[i].xsection=2.38; h1DataDVCS_table4[i].xsection_stat=0.29; h1DataDVCS_table4[i].xsection_sys=0.26;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0.2; h1DataDVCS_table4[i].tmax=0.4; h1DataDVCS_table4[i].t=0.3;
h1DataDVCS_table4[i].xsection=0.67; h1DataDVCS_table4[i].xsection_stat=0.12; h1DataDVCS_table4[i].xsection_sys=0.07;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0.4; h1DataDVCS_table4[i].tmax=0.6; h1DataDVCS_table4[i].t=0.5;
h1DataDVCS_table4[i].xsection=0.24; h1DataDVCS_table4[i].xsection_stat=0.05; h1DataDVCS_table4[i].xsection_sys=0.03;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=50; h1DataDVCS_table4[i].Wmax=85; h1DataDVCS_table4[i].W=70;
h1DataDVCS_table4[i].tmin=0.6; h1DataDVCS_table4[i].tmax=1; h1DataDVCS_table4[i].t=0.8;
h1DataDVCS_table4[i].xsection=0.07; h1DataDVCS_table4[i].xsection_stat=0.01; h1DataDVCS_table4[i].xsection_sys=0.02;
i++;
//-
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0; h1DataDVCS_table4[i].tmax=0.2; h1DataDVCS_table4[i].t=0.1;
h1DataDVCS_table4[i].xsection=2.98; h1DataDVCS_table4[i].xsection_stat=0.49; h1DataDVCS_table4[i].xsection_sys=0.85;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0.2; h1DataDVCS_table4[i].tmax=0.4; h1DataDVCS_table4[i].t=0.3;
h1DataDVCS_table4[i].xsection=0.89; h1DataDVCS_table4[i].xsection_stat=0.17; h1DataDVCS_table4[i].xsection_sys=0.17;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0.4; h1DataDVCS_table4[i].tmax=0.6; h1DataDVCS_table4[i].t=0.5;
h1DataDVCS_table4[i].xsection=0.44; h1DataDVCS_table4[i].xsection_stat=0.08; h1DataDVCS_table4[i].xsection_sys=0.08;
i++;
h1DataDVCS_table4[i].Q2min=11; h1DataDVCS_table4[i].Q2max=80; h1DataDVCS_table4[i].Q2=20.;
h1DataDVCS_table4[i].Wmin=85; h1DataDVCS_table4[i].Wmax=140; h1DataDVCS_table4[i].W=100;
h1DataDVCS_table4[i].tmin=0.6; h1DataDVCS_table4[i].tmax=1; h1DataDVCS_table4[i].t=0.8;
h1DataDVCS_table4[i].xsection=0.06; h1DataDVCS_table4[i].xsection_stat=0.02; h1DataDVCS_table4[i].xsection_sys=0.02;
i++;
vec4 = h1DataDVCS_table4;
}
void loadPhiData(vector& vec1)
{
vector ZEUSDataPhi_table1(7+6+6+6+9);
unsigned int i = 0;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=45;
ZEUSDataPhi_table1[i].W=40;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=76.4;
ZEUSDataPhi_table1[i].xsection_stat=6.5;
ZEUSDataPhi_table1[i].xsection_sys=5.9;
i++;
ZEUSDataPhi_table1[i].Wmin=45;
ZEUSDataPhi_table1[i].Wmax=55;
ZEUSDataPhi_table1[i].W=50;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=101.2;
ZEUSDataPhi_table1[i].xsection_stat=7.7;
ZEUSDataPhi_table1[i].xsection_sys=8.1;
i++;
ZEUSDataPhi_table1[i].Wmin=55;
ZEUSDataPhi_table1[i].Wmax=65;
ZEUSDataPhi_table1[i].W=60;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=101.9;
ZEUSDataPhi_table1[i].xsection_stat=8.3;
ZEUSDataPhi_table1[i].xsection_sys=7.4;
i++;
ZEUSDataPhi_table1[i].Wmin=65;
ZEUSDataPhi_table1[i].Wmax=75;
ZEUSDataPhi_table1[i].W=70;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=112.8;
ZEUSDataPhi_table1[i].xsection_stat=9.4;
ZEUSDataPhi_table1[i].xsection_sys=7.6;
i++;
ZEUSDataPhi_table1[i].Wmin=75;
ZEUSDataPhi_table1[i].Wmax=85;
ZEUSDataPhi_table1[i].W=80;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=107.;
ZEUSDataPhi_table1[i].xsection_stat=11.;
ZEUSDataPhi_table1[i].xsection_sys=10.;
i++;
ZEUSDataPhi_table1[i].Wmin=85;
ZEUSDataPhi_table1[i].Wmax=95;
ZEUSDataPhi_table1[i].W=90;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=122;
ZEUSDataPhi_table1[i].xsection_stat=11.;
ZEUSDataPhi_table1[i].xsection_sys=10.;
i++;
ZEUSDataPhi_table1[i].Wmin=95;
ZEUSDataPhi_table1[i].Wmax=105;
ZEUSDataPhi_table1[i].W=100;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=110.;
ZEUSDataPhi_table1[i].xsection_stat=11.;
ZEUSDataPhi_table1[i].xsection_sys=8.;
i++;
//***** Q2 = 3.8 ******************
ZEUSDataPhi_table1[i].Wmin=40;
ZEUSDataPhi_table1[i].Wmax=50;
ZEUSDataPhi_table1[i].W=45;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=5;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=47.;
ZEUSDataPhi_table1[i].xsection_stat=4.1;
ZEUSDataPhi_table1[i].xsection_sys=2.5;
i++;
ZEUSDataPhi_table1[i].Wmin=50;
ZEUSDataPhi_table1[i].Wmax=60;
ZEUSDataPhi_table1[i].W=55;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=5;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=44.3;
ZEUSDataPhi_table1[i].xsection_stat=4.5;
ZEUSDataPhi_table1[i].xsection_sys=3.6;
i++;
ZEUSDataPhi_table1[i].Wmin=60;
ZEUSDataPhi_table1[i].Wmax=70;
ZEUSDataPhi_table1[i].W=65;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=5;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=56.7;
ZEUSDataPhi_table1[i].xsection_stat=5.1;
ZEUSDataPhi_table1[i].xsection_sys=4.;
i++;
ZEUSDataPhi_table1[i].Wmin=70;
ZEUSDataPhi_table1[i].Wmax=85;
ZEUSDataPhi_table1[i].W=77.5;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=5;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=62.3;
ZEUSDataPhi_table1[i].xsection_stat=5.;
ZEUSDataPhi_table1[i].xsection_sys=4.3;
i++;
ZEUSDataPhi_table1[i].Wmin=85;
ZEUSDataPhi_table1[i].Wmax=100;
ZEUSDataPhi_table1[i].W=92.5;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=5;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=57.4;
ZEUSDataPhi_table1[i].xsection_stat=5.4;
ZEUSDataPhi_table1[i].xsection_sys=4.3;
i++;
ZEUSDataPhi_table1[i].Wmin=100;
ZEUSDataPhi_table1[i].Wmax=115;
ZEUSDataPhi_table1[i].W=107.5;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=5;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=59.0;
ZEUSDataPhi_table1[i].xsection_stat=6.2;
ZEUSDataPhi_table1[i].xsection_sys=3.;
i++;
//*******Q2=6.5******************
ZEUSDataPhi_table1[i].Wmin=45;
ZEUSDataPhi_table1[i].Wmax=55;
ZEUSDataPhi_table1[i].W=50;
ZEUSDataPhi_table1[i].Q2min=5;
ZEUSDataPhi_table1[i].Q2max=9;
ZEUSDataPhi_table1[i].Q2=6.5;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=16.4;
ZEUSDataPhi_table1[i].xsection_stat=1.8;
ZEUSDataPhi_table1[i].xsection_sys=1.1;
i++;
ZEUSDataPhi_table1[i].Wmin=55;
ZEUSDataPhi_table1[i].Wmax=70;
ZEUSDataPhi_table1[i].W=62.5;
ZEUSDataPhi_table1[i].Q2min=5;
ZEUSDataPhi_table1[i].Q2max=9;
ZEUSDataPhi_table1[i].Q2=6.5;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=19.1;
ZEUSDataPhi_table1[i].xsection_stat=1.7;
ZEUSDataPhi_table1[i].xsection_sys=1.5;
i++;
ZEUSDataPhi_table1[i].Wmin=70;
ZEUSDataPhi_table1[i].Wmax=85;
ZEUSDataPhi_table1[i].W=77.5;
ZEUSDataPhi_table1[i].Q2min=5;
ZEUSDataPhi_table1[i].Q2max=9;
ZEUSDataPhi_table1[i].Q2=6.5;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=19.6;
ZEUSDataPhi_table1[i].xsection_stat=1.9;
ZEUSDataPhi_table1[i].xsection_sys=1.6;
i++;
ZEUSDataPhi_table1[i].Wmin=85;
ZEUSDataPhi_table1[i].Wmax=100;
ZEUSDataPhi_table1[i].W=92.5;
ZEUSDataPhi_table1[i].Q2min=5;
ZEUSDataPhi_table1[i].Q2max=9;
ZEUSDataPhi_table1[i].Q2=6.5;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=21.6;
ZEUSDataPhi_table1[i].xsection_stat=2.3;
ZEUSDataPhi_table1[i].xsection_sys=1.1;
i++;
ZEUSDataPhi_table1[i].Wmin=100;
ZEUSDataPhi_table1[i].Wmax=115;
ZEUSDataPhi_table1[i].W=107.5;
ZEUSDataPhi_table1[i].Q2min=5;
ZEUSDataPhi_table1[i].Q2max=9;
ZEUSDataPhi_table1[i].Q2=6.5;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=23.1;
ZEUSDataPhi_table1[i].xsection_stat=2.5;
ZEUSDataPhi_table1[i].xsection_sys=1.3;
i++;
ZEUSDataPhi_table1[i].Wmin=115;
ZEUSDataPhi_table1[i].Wmax=135;
ZEUSDataPhi_table1[i].W=125;
ZEUSDataPhi_table1[i].Q2min=5;
ZEUSDataPhi_table1[i].Q2max=9;
ZEUSDataPhi_table1[i].Q2=6.5;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=25.3;
ZEUSDataPhi_table1[i].xsection_stat=3.5;
ZEUSDataPhi_table1[i].xsection_sys=3.9;
i++;
//*******Q2=13.0******************
ZEUSDataPhi_table1[i].Wmin=50;
ZEUSDataPhi_table1[i].Wmax=60;
ZEUSDataPhi_table1[i].W=55;
ZEUSDataPhi_table1[i].Q2min=9;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=13.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=5.05;
ZEUSDataPhi_table1[i].xsection_stat=.73;
ZEUSDataPhi_table1[i].xsection_sys=.37;
i++;
ZEUSDataPhi_table1[i].Wmin=60;
ZEUSDataPhi_table1[i].Wmax=75;
ZEUSDataPhi_table1[i].W=67.5;
ZEUSDataPhi_table1[i].Q2min=9;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=13.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=4.96;
ZEUSDataPhi_table1[i].xsection_stat=.59;
ZEUSDataPhi_table1[i].xsection_sys=.27;
i++;
ZEUSDataPhi_table1[i].Wmin=75;
ZEUSDataPhi_table1[i].Wmax=90;
ZEUSDataPhi_table1[i].W=82.5;
ZEUSDataPhi_table1[i].Q2min=9;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=13.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=6.12;
ZEUSDataPhi_table1[i].xsection_stat=.81;
ZEUSDataPhi_table1[i].xsection_sys=.4;
i++;
ZEUSDataPhi_table1[i].Wmin=90;
ZEUSDataPhi_table1[i].Wmax=105;
ZEUSDataPhi_table1[i].W=97.5;
ZEUSDataPhi_table1[i].Q2min=9;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=13.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=6.82;
ZEUSDataPhi_table1[i].xsection_stat=.85;
ZEUSDataPhi_table1[i].xsection_sys=.4;
i++;
ZEUSDataPhi_table1[i].Wmin=105;
ZEUSDataPhi_table1[i].Wmax=125;
ZEUSDataPhi_table1[i].W=115;
ZEUSDataPhi_table1[i].Q2min=9;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=13.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=6.33;
ZEUSDataPhi_table1[i].xsection_stat=.78;
ZEUSDataPhi_table1[i].xsection_sys=.5;
i++;
ZEUSDataPhi_table1[i].Wmin=125;
ZEUSDataPhi_table1[i].Wmax=145;
ZEUSDataPhi_table1[i].W=135;
ZEUSDataPhi_table1[i].Q2min=9;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=13.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=6.65;
ZEUSDataPhi_table1[i].xsection_stat=.97;
ZEUSDataPhi_table1[i].xsection_sys=.52;
i++;
//*****Q2 cross-section**********
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=2;
ZEUSDataPhi_table1[i].Q2max=3;
ZEUSDataPhi_table1[i].Q2=2.4;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=105.5;
ZEUSDataPhi_table1[i].xsection_stat=3.4;
ZEUSDataPhi_table1[i].xsection_sys=6.;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=3;
ZEUSDataPhi_table1[i].Q2max=4.5;
ZEUSDataPhi_table1[i].Q2=3.6;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=57.6;
ZEUSDataPhi_table1[i].xsection_stat=2.4;
ZEUSDataPhi_table1[i].xsection_sys=3.5;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=4.5;
ZEUSDataPhi_table1[i].Q2max=6;
ZEUSDataPhi_table1[i].Q2=5.2;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=31.1;
ZEUSDataPhi_table1[i].xsection_stat=1.8;
ZEUSDataPhi_table1[i].xsection_sys=1.8;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=6;
ZEUSDataPhi_table1[i].Q2max=8;
ZEUSDataPhi_table1[i].Q2=6.9;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=17.9;
ZEUSDataPhi_table1[i].xsection_stat=1.1;
ZEUSDataPhi_table1[i].xsection_sys=1.;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=8;
ZEUSDataPhi_table1[i].Q2max=11;
ZEUSDataPhi_table1[i].Q2=9.2;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=11.06;
ZEUSDataPhi_table1[i].xsection_stat=0.73;
ZEUSDataPhi_table1[i].xsection_sys=0.56;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=11;
ZEUSDataPhi_table1[i].Q2max=15;
ZEUSDataPhi_table1[i].Q2=12.6;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=6.42;
ZEUSDataPhi_table1[i].xsection_stat=0.52;
ZEUSDataPhi_table1[i].xsection_sys=0.24;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=15;
ZEUSDataPhi_table1[i].Q2max=20;
ZEUSDataPhi_table1[i].Q2=17.1;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=2.5;
ZEUSDataPhi_table1[i].xsection_stat=0.37;
ZEUSDataPhi_table1[i].xsection_sys=0.22;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=20;
ZEUSDataPhi_table1[i].Q2max=30;
ZEUSDataPhi_table1[i].Q2=24.;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=0.98;
ZEUSDataPhi_table1[i].xsection_stat=0.19;
ZEUSDataPhi_table1[i].xsection_sys=0.05;
i++;
ZEUSDataPhi_table1[i].Wmin=35;
ZEUSDataPhi_table1[i].Wmax=142;
ZEUSDataPhi_table1[i].W=75;
ZEUSDataPhi_table1[i].Q2min=30;
ZEUSDataPhi_table1[i].Q2max=70;
ZEUSDataPhi_table1[i].Q2=38.8;
ZEUSDataPhi_table1[i].tmin=0;
ZEUSDataPhi_table1[i].tmax=1;
ZEUSDataPhi_table1[i].xsection=0.37;
ZEUSDataPhi_table1[i].xsection_stat=0.13;
ZEUSDataPhi_table1[i].xsection_sys=0.04;
i++;
vec1 = ZEUSDataPhi_table1;
}