Page MenuHomeHEPForge
Files F8310250

No OneTemporary
Actions

View Options

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 <http://www.gnu.org/licenses/>.
#
# 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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update: Sat Aug 13 02:18:37 2011
//==============================================================================
//
// Main program to create lambda lookup tables.
// [Developer only]
//==============================================================================
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include "Amplitudes.h"
#include "TH1D.h"
#include "TFile.h"
#include "Constants.h"
#include "Table.h"
#include "Settings.h"
#include "TableGeneratorSettings.h"
#include "Enumerations.h"
#include "Kinematics.h"
#include <limits>
#include <cmath>
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
int main(int argc, char *argv[]){
//To explicitly delete all new histograms myself:
TH1::AddDirectory(false);
//
// Check arguments
//
char* runcard;
if (argc != 4) {
cout<<"Usage: lambdaGeneratorMain runcard startBin endBin"<<endl;
return 2;
}
else{
runcard = argv[1];
settings->setStartBin(atoi(argv[2]));
settings->setEndBin(atoi(argv[3]));
}
settings->readSettingsFromFile(runcard);
settings->setA(1);
settings->setModesToCalculate(1);
settings->setNumberOfConfigurations(1);
settings->consolidateSettings();
int nBinQ2 = settings->Q2bins();
int nBinW2 = settings->W2bins();
int nBinT= settings->tbins();
double Q2min= settings->Q2min();
double Q2max= settings->Q2max();
double Wmin= settings->Wmin();
double Wmax= settings->Wmax();
double W2min=Wmin*Wmin;
double W2max=Wmax*Wmax;
double tmin= settings->tmin();
double tmax= settings->tmax();
unsigned int massA=settings->A();
int vmPDG = settings->vectorMesonId();
DipoleModelType model = settings->dipoleModelType();
DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
int startingBin=settings->startBin();
int endingBin=settings->endBin();
int modes=settings->modesToCalculate();
Table tableT;
Table tableL;
bool logQ2=true, logW2=false, logT=false, logC=false;
//set filenames for the tables
string rootfile=settings->rootfile();
ostringstream filenameT, filenameL, filenameT2, filenameL2;
filenameT.str("");
filenameT<<rootfile<<"_bin"<<startingBin<<"_T.root";
filenameL.str("");
filenameL<<rootfile<<"_bin"<<startingBin<<"_L.root";
(void) tableT.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, logC, // all bools
lambda_A, transverse,
massA, vmPDG, model, pset,
filenameT.str().c_str());
(void) tableL.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, logC, // all bools
lambda_A, longitudinal,
massA, vmPDG, model, pset,
filenameL.str().c_str());
//Create and initialize the amplitudes calculator:
Amplitudes amps;
//Generate the the nucleon configurations:
amps.generateConfigurations();
// Print out settings:
cout<<"Tables will be generated for:"<<endl;
cout<<"Nucleus mass A="<<massA<<endl;
cout<<"Mode to calculate: "<<modes<<endl;
cout<<"Vector Meson Id: "<<vmPDG<<endl;
cout<<"Bins: "<<startingBin<<"-"<<endingBin<<endl;
cout<<"Q2 range: ["<<Q2min<<", "<<Q2max<<"], "<<nBinQ2<<" bins."<<endl;
cout<<"W2 range: ["<<W2min<<", "<<W2max<<"], "<<nBinW2<<" bins."<<endl;
cout<<" t range: ["<<tmin<<", "<<tmax<<"], "<<nBinT<<" bins."<<endl;
cout<<"Dipole Model: "<<model<<endl;
cout<<endl;
cout<<"Filling tables..."<<endl;
time_t tableFillStart = time(0);
for(int i=startingBin; i<endingBin; i++) {
double Q2, W2, t;
tableT.binCenter(i, Q2, W2, t);
double kinematicPoint[3]={t, Q2, W2};
//calculate derivative of amplitude in point:
double hplus, hminus;
hplus=hminus=(W2max-W2min)/(4*1e4); //This value comes from tests of large and small W2
hminus=min(hminus, W2-W2min);
hplus=min(hplus, W2max-W2);
hminus -= numeric_limits<float>::epsilon();
hplus -= numeric_limits<float>::epsilon();
//calculate contents and fill tables:
double lambda[2]={0, 0};
kinematicPoint[2]= W2+hplus;
amps.calculate(kinematicPoint);
double ampPlus[2]={0, 0};
ampPlus[0]=amps.amplitudeT();
ampPlus[1]=amps.amplitudeL();
kinematicPoint[2]=W2-hminus;
amps.calculate(kinematicPoint);
double ampMinus[2]={0, 0};
ampMinus[0]=amps.amplitudeT();
ampMinus[1]=amps.amplitudeL();
for(int j=0; j<2; j++){
PR(ampMinus[j]);
PR(ampPlus[j]);
// Don't calculate unless numerically viable:
if(ampPlus[j]==0 || ampMinus[j]==0){
lambda[j]=0;
}
else{
//Calculate derivate d(logA/dW2)=log(A+/A-)/(h+-h-):
double derivate=log(abs(ampPlus[j]/ampMinus[j]))/(hplus+hminus);
// Finally calculate lambda:
double jacobian = (W2-protonMass2+Q2);
lambda[j] = jacobian*derivate;
}
}
tableT.fill(i, lambda[0]);
tableL.fill(i, lambda[1]);
if(i%static_cast<int>((endingBin-startingBin)/10+1)==0)
cout<<"Processed "<<100*(i-startingBin)/(endingBin-startingBin)<<" %"<<endl;
if(settings->verbose()){
cout << "bin = " << setw(8) << left << i;
cout << "Q2 = " << setw(10) << left << fixed << setprecision(3) << Q2;
cout << "W2 = " << setw(11) << left << fixed << W2;
cout << "t = " << setw(14) << left << scientific << t;
cout << "lambdaT = " << setw(12) << left << scientific << lambda[0]<<endl;
cout << "bin = " << setw(8) << left << i;
cout << "Q2 = " << setw(10) << left << fixed << setprecision(3) << Q2;
cout << "W2 = " << setw(11) << left << fixed << W2;
cout << "t = " << setw(14) << left << scientific << t;
cout << "lambdaL = " << setw(12) << left << scientific << lambda[1]<<endl;
}
}
cout<<"CPU Time/Entry: " << double(time(0)-tableFillStart)/endingBin<<" s/entry."<<endl;
cout<<"Total time: "<<double(time(0)-tableFillStart)/60./60.<<" h"<<endl;
tableT.write();
tableL.write();
return 0;
}
Index: trunk/examples/sartreTest.cpp
===================================================================
--- trunk/examples/sartreTest.cpp (revision 370)
+++ trunk/examples/sartreTest.cpp (revision 371)
@@ -1,536 +1,542 @@
//==============================================================================
// sartreTest.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 <http://www.gnu.org/licenses/>.
//
// 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 <iostream>
#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 <fstream>
#include <cmath>
#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; i<n1; i++) {
tableA.binCenter(i, xpom, t);
//cout << i << '\t' << xpom << '\t' << t << endl;
tableA.fill(i, testfun(xpom, t));
tableA2.fill(i, testfun(xpom, t));
}
cout << "Checking interpolation" << endl;
TRandom2 rndm;
for (int i=0; i<100; i++) {
xpom = rndm.Uniform(log(1e-5), log(0.2));
xpom = exp(xpom);
t = rndm.Uniform(-0.051, 0);
double cinter = tableA.get(xpom, t);
cout << "xp = " << xpom << ", t =" << t << ": real=" << testfun(xpom, t)
<< ", interpolated: " << cinter
<< " (" << 100*(cinter-testfun(xpom, t))/testfun(xpom, t) << "%)" << endl;
}
cout << "Writing UPC tables" << endl;
tableA.write();
tableA2.write(); // filename no needed - was passed along in create()
#endif
//
// Test creating and writing tables
//
/*
unsigned int Table::create(int nbinsQ2, double Q2min, double Q2max,
int nbinsW2, double W2min, double W2max,
int nbinsT, double tmin, double tmax,
bool logQ2, bool logW2, bool logt, bool logC,
AmplitudeMoment mom, GammaPolarization pol,
unsigned int A, int vm,
DipoleModelType model, DipoleModelParameterSet pset,
const char* filename, unsigned char priority)
*/
#if defined (TEST_WRITETABLE)
Table table1, table2, table3;
int n = table1.create(10, 1, 100, 10, 1, 80, 17, 0, 0.051,
true, true, false, false,
mean_A2, transverse, 197, 411, bSat, KMW); // no filename, no priority (implies 0)
n = table2.create(10, 1, 100, 10, 1, 80, 10, 0, 2,
false, false, false, false,
mean_A, longitudinal, 197, 411, bCGC, KMW, "table2.root"); // filename but no priority (implies 0)
n = table3.create(10, 1, 100, 10, 1, 80, 10, 0, 2,
false, false, false, false,
lambda_A, transverse, 197, 411, bCGC, KMW, 0, 2); // no filename but priority 2
PR(n);
table1.setAutobackup("table1", 5);
table2.setAutobackup("table2", 10);
table2.setAutobackup("table3", 10);
double Q2, W2, t;
for (int i=0; i<n; i++) {
table1.binCenter(i, Q2, W2, t);
// cout << i << '\t' << Q2 << '\t' << W2 << '\t' << t << endl;
table1.fill(i, i+1);
table2.fill(i, i+1);
table3.fill(i, i+1);
}
table1.write("table1.root");
table2.write(); // filename no needed - was passed along in create()
table3.write("table3.root");
#endif
//
// Test final state generator
//
#if defined(TEST_FINALSTATE)
ExclusiveFinalStateGenerator generator;
Event event;
event.eventNumber = 10345;
event.t = -0.2;
event.y = 0.6;
event.Q2 = 1;
event.W = 80;
event.x = Kinematics::x(event.Q2, event.W*event.W);
event.xpom = 0.003;
event.polarisation = 'T';
event.particles.resize(2);
event.particles[0].index = 0;
event.particles[1].index = 1;
event.particles[0].pdgId = 11;
event.particles[1].pdgId = -2212;
event.particles[0].status = 4;
event.particles[1].status = 4;
event.particles[0].p = settings->eBeam();
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<InterTestPoint> 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<InterTestPoint> 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<InterTestPoint> 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<InterTestPoint> 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; i<vector_L.size(); i++) {
double res = coll.get(vector_L[i].Q2, vector_L[i].W2, vector_L[i].t, longitudinal, mean_A);
histoL.Fill((vector_L[i].value - res)/vector_L[i].value, 1.);
}
TH1D histoL2("histoL2", "L2", 100, -0.1, 0.1);
for (unsigned int i=0; i<vector_L2.size(); i++) {
double res = coll.get(vector_L2[i].Q2, vector_L2[i].W2, vector_L2[i].t, longitudinal, mean_A2);
histoL2.Fill((vector_L2[i].value - res)/vector_L2[i].value, 1.);
}
TH1D histoT("histoT", "T", 100, -0.1, 0.1);
for (unsigned int i=0; i<vector_T.size(); i++) {
double res = coll.get(vector_T[i].Q2, vector_T[i].W2, vector_T[i].t, transverse, mean_A);
histoT.Fill((vector_T[i].value - res)/vector_T[i].value, 1.);
}
TH1D histoT2("histoT2", "T2", 100, -0.1, 0.1);
for (unsigned int i=0; i<vector_T2.size(); i++) {
double res = coll.get(vector_T2[i].Q2, vector_T2[i].W2, vector_T2[i].t, transverse, mean_A2);
histoT2.Fill((vector_T2[i].value - res)/vector_T2[i].value, 1.);
}
cout << "histos written to file 'interTest.root'" << endl;
hfile->Write();
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<InterTestSlope> 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<InterTestSlope> 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; i<vector_L.size(); i++) {
double sartreEstimate = xSection.logDerivateOfAmplitude(vector_L[i].t, vector_L[i].Q2, vector_L[i].W2, longitudinal);
double trueValue = vector_L[i].slope;
//cout << "true=" << trueValue << ", sartre=" << sartreEstimate << ", rel-diff=" << (trueValue-sartreEstimate)/trueValue << endl;
histoL.Fill((trueValue-sartreEstimate)/trueValue, 1.);
}
TH1D histoT("histoT", "T", 100, -0.1, 0.1);
for (unsigned int i=0; i<vector_T.size(); i++) {
double sartreEstimate = xSection.logDerivateOfAmplitude(vector_T[i].t, vector_T[i].Q2, vector_T[i].W2, transverse);
double trueValue = vector_T[i].slope;
//cout << "true=" << trueValue << ", sartre=" << sartreEstimate << ", rel-diff=" << (trueValue-sartreEstimate)/trueValue << endl;
histoT.Fill((trueValue-sartreEstimate)/trueValue, 1.);
}
cout << "histos written to file 'slopeTest.root'" << endl;
hfile->Write();
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<Nucleon> 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; k<nucleons.size(); k++) {
double radius = nucleons[k].position().Mag();
double weight = 1./(radius*radius*binsize)/numberOfNuclei;
h1.Fill(radius, weight);
h2xy.Fill(nucleons[k].position().X(), nucleons[k].position().Y(),1);
h2xz.Fill(nucleons[k].position().X(), nucleons[k].position().Z(),1);
h2yz.Fill(nucleons[k].position().Y(), nucleons[k].position().Z(),1);
if (i==1) {
h2xy_1.Fill(nucleons[k].position().X(), nucleons[k].position().Y(),1);
h2xz_1.Fill(nucleons[k].position().X(), nucleons[k].position().Z(),1);
h2yz_1.Fill(nucleons[k].position().Y(), nucleons[k].position().Z(),1);
}
if (i==2) {
h2xy_2.Fill(nucleons[k].position().X(), nucleons[k].position().Y(),1);
h2xz_2.Fill(nucleons[k].position().X(), nucleons[k].position().Z(),1);
h2yz_2.Fill(nucleons[k].position().Y(), nucleons[k].position().Z(),1);
}
if (i==3) {
h2xy_3.Fill(nucleons[k].position().X(), nucleons[k].position().Y(),1);
h2xz_3.Fill(nucleons[k].position().X(), nucleons[k].position().Z(),1);
h2yz_3.Fill(nucleons[k].position().Y(), nucleons[k].position().Z(),1);
}
}
}
cout << "Processed " << numberOfNuclei << " nuclei." << endl;
hfile->Write();
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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Main program to create impact parameter lookup table.
// [Developer only]
//==============================================================================
#include <iostream>
#include <cmath>
#include <sstream>
#include <cstdlib>
#include <vector>
#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="<<A<<" and for "<<numConf<<" configurations."<<endl;
createLookUpTableOnFile(A, numConf);
return 0;
}
double overlapFunctionT(double b)
{
// Gaussian shape for proton
// b in GeV
double BG = gParameters->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<TH2F*> 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; iNuclei<numberOfConfigurations; iNuclei++) {
histoName.str( "" );
histoName << "Configuration_" << iNuclei;
hConfigurations[iNuclei] = new TH2F(histoName.str().c_str(), "Sum_i^A T_p(b-b_i)",
numBbins, 0., bRange, numAngleBins, 0., angleRange);
//Generate a configuration:
while(!myNucleus.generate()){}
for(int ib=1; ib <= numBbins; ib++){
double b=hConfigurations[iNuclei]->GetXaxis()->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<<double(iNuclei)/numberOfConfigurations*100<<"% done."<<endl;
}
} //iNuclei
cout<<"100% done."<<endl<<endl;
cout<<"Writing configurations to file...";
//Open the file and write to it:
TFile *lufile = 0;
ostringstream filename;
filename.str("");
filename << "bSat_bDependence_A" << A <<".root";
lufile = new TFile(filename.str().c_str(), "RECREATE");
for(unsigned int i=0; i<numberOfConfigurations; i++){
hConfigurations[i]->Write();
delete hConfigurations[i];
hConfigurations[i] = 0;
}
lufile->Close();
cout<<" done, bye!"<<endl;
}
Index: trunk/examples/tableGeneratorMain.cpp
===================================================================
--- trunk/examples/tableGeneratorMain.cpp (revision 370)
+++ trunk/examples/tableGeneratorMain.cpp (revision 371)
@@ -1,475 +1,475 @@
//==============================================================================
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Main program to create amplitude lookup tables.
// [Developer only]
//
// Usage:
// tableGeneratorMain runcard startBin endBin
//
// Bins run from 0 to nbin-1.
// Loop fill all bins from startBIn to endBin (including endBin).
// If endBin exceeds the size of the table it set to the nbin-1.
//
//==============================================================================
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <cmath>
#include <iomanip>
#include "Amplitudes.h"
#include "TROOT.h"
#include "TH1D.h"
#include "TFile.h"
#include "Nucleus.h"
#include "Constants.h"
#include "Table.h"
#include "TableGeneratorSettings.h"
#include "Enumerations.h"
#include "DglapEvolution.h"
#include "Version.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
int main(int argc, char *argv[]) {
//
// Print header
//
time_t theStartTime = time(0);
string ctstr(ctime(&theStartTime));
ctstr.erase(ctstr.size()-1, 1);
cout << "/========================================================================\\" << endl;
cout << "| |" << endl;
cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
cout << "| |" << endl;
cout << "| An event generator for exclusive diffractive vector meson production |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
cout << "| Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich |" << endl;
cout << "| |" << endl;
cout << "| This program is free software: you can redistribute it and/or modify |" << endl;
cout << "| it under the terms of the GNU General Public License as published by |" << endl;
cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl;
cout << "| any later version. |" << endl;
cout << "| |" << endl;
cout << "| Code compiled on " << setw(12) << left << __DATE__;
cout << setw(41) << left << __TIME__ << right << '|' << endl;
cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
cout << "| |" << endl;
cout << "| -------------------------------------------------------------------- |" << endl;
cout << "| |" << endl;
cout << "| Sartre Table Generator (Experts only) |" << endl;
cout << "| |" << endl;
cout << "\\========================================================================/" << endl;
TH1::AddDirectory(false); // to explicitly delete all new histograms by hand
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
//
// Check arguments
//
char* runcard;
if (argc != 4) {
cout << "Usage: tableGeneratorMain runcard startBin endBin" << endl;
return 2;
}
else {
runcard = argv[1];
settings->setStartBin(atoi(argv[2]));
settings->setEndBin(atoi(argv[3]));
}
cout << "Reading settings from runcard." << endl;
settings->readSettingsFromFile(runcard);
settings->consolidateSettings();
int nBinQ2 = settings->Q2bins();
int nBinW2 = settings->W2bins();
int nBinT = settings->tbins();
double Q2min = settings->Q2min();
double Q2max = settings->Q2max();
double Wmin = settings->Wmin();
double Wmax = settings->Wmax();
double W2min = Wmin*Wmin;
double W2max = Wmax*Wmax;
double tmin = settings->tmin();
double tmax = settings->tmax();
unsigned int massA = settings->A();
int vmPDG = settings->vectorMesonId();
DipoleModelType model = settings->dipoleModelType();
DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
int startingBin = settings->startBin();
int endingBin = settings->endBin();
int modes = settings->modesToCalculate();
unsigned char priority = settings->priority();
settings->list();
//
// Check if lambda tables can be calculated
//
bool createLambdaTables = true;
if (massA == 1 && modes == 1 && settings->numberOfConfigurations() == 1) {
cout << "\nLambda tables will be generated." << endl;
}
else {
cout << "\nLambda tables will not be generated. Requires A = 1, mode = 1, and 1 configuration only." << endl;
createLambdaTables = false;
}
//
// Check bins
//
// Table's bin indices run from 0 ... nbins-1
//
int maxbins = nBinQ2*nBinW2*nBinT;
if (endingBin >= maxbins) {
cout << "Warning, given end bin (" << endingBin << ") exceeds the size of the table." << endl;
cout << " set to maximum value (" << maxbins-1 << ") now." << endl;
endingBin = maxbins-1;
}
//
// Define all tables. Depending on tableset type some
// will not be written but we define them all anyway.
//
Table tableT;
Table tableL;
Table tableT2;
Table tableL2;
Table tableVarT;
Table tableVarL;
Table tableLambdaT;
Table tableLambdaL;
bool logQ2=true, logW2=true, logT=false, logC=true;
//
// Set filenames for the tables.
// Be as desciptive as possible. Helps when mass producing tables.
//
// We create all tables and decide later what
// gets written and what not.
//
string rootfile=settings->rootfile();
rootfile += "_" + settings->dipoleModelName();
rootfile += "_" + settings->dipoleModelParameterSetName();
ostringstream filenameT, filenameL, filenameT2,
filenameL2, filenameVarT, filenameVarL,
filenameLambdaT, filenameLambdaL;
filenameT.str("");
filenameT << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin <<"_T.root";
filenameL.str("");
filenameL << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_L.root";
filenameT2.str("");
filenameT2 << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_T2.root";
filenameL2.str("");
filenameL2 << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_L2.root";
filenameVarT.str("");
filenameVarT << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_VarT.root";
filenameVarL.str("");
filenameVarL << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_VarL.root";
filenameLambdaT.str("");
filenameLambdaT << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_LambdaT.root";
filenameLambdaL.str("");
filenameLambdaL << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_LambdaL.root";
(void) tableT.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, logC, // all bools
mean_A, transverse,
massA, vmPDG, model, pset,
filenameT.str().c_str(),
priority);
(void) tableL.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, logC, // all bools
mean_A, longitudinal,
massA, vmPDG, model, pset,
filenameL.str().c_str(),
priority);
(void) tableT2.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, logC, // all bools
mean_A2, transverse,
massA, vmPDG, model, pset,
filenameT2.str().c_str(),
priority);
(void) tableL2.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, logC, // all bools
mean_A2, longitudinal,
massA, vmPDG, model, pset,
filenameL2.str().c_str(),
priority);
bool contentVar = modes == 1 ? false : logC; // if mode = 1 they will be 0 so lin is better
(void) tableVarT.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, contentVar, // all bools
variance_A, transverse,
massA, vmPDG, model, pset,
filenameVarT.str().c_str(),
priority);
(void) tableVarL.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, contentVar, // all bools
variance_A, longitudinal,
massA, vmPDG, model, pset,
filenameVarL.str().c_str(),
priority);
double contentLambda = false; // content lin
(void) tableLambdaT.create(nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinT, tmin, tmax,
logQ2, logW2, logT, contentLambda, // all bools
lambda_A, transverse,
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="<<massA<<endl;
cout << "\tModes to calculate: " << modes << endl;
cout << "\tVector Meson Id: " << vmPDG << endl;
cout << "\tBins: " << startingBin << "-" << endingBin << endl;
cout << "\tQ2 range: [" << Q2min << ", " << Q2max << "], " << nBinQ2 << " bins." << endl;
cout << "\tW2 range: [" << W2min << ", " << W2max << "], " << nBinW2 << " bins." << endl;
cout << "\t t range: [" << tmin << ", " << tmax << "], " << nBinT << " bins." << endl;
cout << "\tDipole model: " << settings->dipoleModelName() << endl;
cout << "\tDipole model parameter set: " << settings->dipoleModelParameterSetName() << endl;
cout << "\tTable set mode: " << settings->tableSetTypeName() << endl;
cout << endl;
//
// Calculate contents and fill tables
//
// Note that we fill all tables. What is written
// at the end is another story.
//
cout << "Start filling tables" << endl;
time_t tableFillStart = time(0);
int nShow = (endingBin - startingBin)/100;
for (int i=startingBin; i<=endingBin; i++) {
if (nShow && (i%nShow == 0 || i == startingBin || i == endingBin))
cout << "processing bin " << i << endl;
double Q2, W2, t;
tableT.binCenter(i, Q2, W2, t);
double kinematicPoint[3]={t, Q2, W2};
amps.calculate(kinematicPoint);
double aT = 0;
double aL = 0;
double aT2 = 0;
double aL2 = 0;
double aVarT = 0;
double aVarL = 0;
//
// For mode = 1 only coherent is
// calculated and incoherent is therefore
// assumed to be 0.
// Hence total == coherent
//
// Recall:
// coherent = A
// total = A2
// incoherent = variance = A2 - A*A
//
aL = amps.amplitudeL();
aT = amps.amplitudeT();
tableT.fill(i, aT);
tableL.fill(i, aL);
if (modes != 1) {
aT2 = amps.amplitudeT2();
aL2 = amps.amplitudeL2();
}
else {
aT2 = aT*aT;
aL2 = aL*aL;
}
tableT2.fill(i, aT2);
tableL2.fill(i, aL2);
if (modes != 1) {
aVarT = aT2-aT*aT;
aVarL = aL2-aL*aL;
}
tableVarT.fill(i, aVarT);
tableVarL.fill(i, aVarL);
//
// Calculate lambda and fill table
//
if (createLambdaTables) {
double hplus, hminus;
hplus=hminus=(W2max-W2min)/(4*1e4); //This value comes from tests of large and small W2
hminus=min(hminus, W2-W2min);
hplus=min(hplus, W2max-W2);
hminus -= numeric_limits<float>::epsilon();
hplus -= numeric_limits<float>::epsilon();
double lambda[2]={0, 0};
kinematicPoint[2]=W2+hplus;
amps.calculate(kinematicPoint);
double ampPlus[2]={0, 0};
ampPlus[0]=amps.amplitudeT();
ampPlus[1]=amps.amplitudeL();
kinematicPoint[2]=W2-hminus;
amps.calculate(kinematicPoint);
double ampMinus[2]={0, 0};
ampMinus[0]=amps.amplitudeT();
ampMinus[1]=amps.amplitudeL();
for (int j=0; j<2; j++) {
// Don't calculate unless numerically viable:
if (ampPlus[j]==0 || ampMinus[j]==0) {
lambda[j]=0;
}
else {
//Calculate derivate d(logA/dW2)=log(A+/A-)/(h+-h-):
double derivate=log(abs(ampPlus[j]/ampMinus[j]))/(hplus+hminus);
// Finally calculate lambda:
double jacobian = (W2-protonMass2+Q2);
lambda[j] = jacobian*derivate;
}
}
tableLambdaT.fill(i, lambda[0]);
tableLambdaL.fill(i, lambda[1]);
}
}
time_t tableFillEnd = time(0);
//
// Report CPU time/cell
//
cout << endl;
cout << "CPU time/bin: "
<< static_cast<double>(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
<< " s" << endl;
cout << "Total time: " << static_cast<double>(tableFillEnd-tableFillStart)/60./60. << " h" << endl << endl;
//
// We write out all tables.
//
// Whoever runs the production can then decide
// later what to keep and what to delete.
// If the desired run mode is total_and_coherent or
// coherent_and_incoherent, if this is just to improve
// a coherent table in some phase space, or if one wants
// to maintain redundancy, all these factor might affect
// your choice.
//
tableT.write();
tableL.write();
tableT2.write();
tableL2.write();
tableVarT.write();
tableVarL.write();
if (createLambdaTables) {
tableLambdaT.write();
tableLambdaL.write();
}
cout << "All tables written" << endl;
cout << "All done. Bye." << endl;
return 0;
}
Index: trunk/examples/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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Example main program. Use to get started.
//
//==============================================================================
#include <iostream>
#include <cmath>
#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; i<myEvent->particles.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 <http://www.gnu.org/licenses/>.
//
// 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 <A> analytical, <A2> averaged over configurations
# 1 <A> only, analytical
# 2 <A2>, and <A> 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 <http://www.gnu.org/licenses/>.
//
// 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 <A> analytical, <A2> averaged over configurations
# 1 <A> only, analytical
# 2 <A2>, and <A> 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 <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
//
// Main program to create amplitude lookup tables.
// [Developer only]
//
// Usage:
// tableGeneratorMain runcard startBin endBin
//
// Bins run from 0 to nbin-1.
// Loop fill all bins from startBIn to endBin (including endBin).
// If endBin exceeds the size of the table it set to the nbin-1.
//
//==============================================================================
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <cmath>
#include <iomanip>
#include "Amplitudes.h"
#include "TROOT.h"
#include "TH1D.h"
#include "TFile.h"
#include "Nucleus.h"
#include "Constants.h"
#include "Table.h"
#include "TableGeneratorSettings.h"
#include "Enumerations.h"
#include "DglapEvolution.h"
#include "Version.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
int main(int argc, char *argv[]) {
//
// Print header
//
time_t theStartTime = time(0);
string ctstr(ctime(&theStartTime));
ctstr.erase(ctstr.size()-1, 1);
cout << "/========================================================================\\" << endl;
cout << "| |" << endl;
cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
cout << "| |" << endl;
cout << "| An event generator for exclusive diffractive vector meson production |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
cout << "| Copyright (C) 2010-2018 Tobias Toll and Thomas Ullrich |" << endl;
cout << "| |" << endl;
cout << "| This program is free software: you can redistribute it and/or modify |" << endl;
cout << "| it under the terms of the GNU General Public License as published by |" << endl;
cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl;
cout << "| any later version. |" << endl;
cout << "| |" << endl;
cout << "| Code compiled on " << setw(12) << left << __DATE__;
cout << setw(41) << left << __TIME__ << right << '|' << endl;
cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
cout << "| |" << endl;
cout << "| -------------------------------------------------------------------- |" << endl;
cout << "| |" << endl;
cout << "| Sartre Table Generator (Experts only) |" << endl;
cout << "| |" << endl;
cout << "\\========================================================================/" << endl;
TH1::AddDirectory(false); // to explicitly delete all new histograms by hand
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
//
// Check arguments
//
char* runcard;
if (argc != 4) {
cout << "Usage: tableGeneratorMain runcard startBin endBin" << endl;
return 2;
}
else {
runcard = argv[1];
settings->setStartBin(atoi(argv[2]));
settings->setEndBin(atoi(argv[3]));
}
cout << "Reading settings from runcard." << endl;
settings->readSettingsFromFile(runcard);
settings->consolidateSettings();
int nBinT = settings->tbins();
int nBinX = settings->xbins();
double tmin = settings->tmin();
double tmax = settings->tmax();
double xmin = settings->xmin();
double xmax = settings->xmax();
unsigned int massA = settings->A();
int vmPDG = settings->vectorMesonId();
DipoleModelType model = settings->dipoleModelType();
DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
int startingBin = settings->startBin();
int endingBin = settings->endBin();
int modes = settings->modesToCalculate();
unsigned char priority = settings->priority();
settings->list();
//
// Check if lambda tables can be calculated
//
bool createLambdaTables = true;
if (massA == 1 && modes != 1 && settings->numberOfConfigurations() == 1) {
cout << "\nLambda tables will be generated." << endl;
}
else {
cout << "\nLambda tables will not be generated. Requires A = 1, mode != 1, and 1 configuration only." << endl;
createLambdaTables = false;
}
//
// Check bins
//
// Table's bin indices run from 0 ... nbins-1
//
int maxbins = nBinX*nBinT;
if (endingBin >= maxbins) {
cout << "Warning, given end bin (" << endingBin << ") exceeds the size of the table." << endl;
cout << " set to maximum value (" << maxbins-1 << ") now." << endl;
endingBin = maxbins-1;
}
//
// Define all tables. Depending on tableset type some
// will not be written but we define them all anyway.
//
Table tableT;
Table tableT2;
Table tableVarT;
Table tableLambdaRealT;
Table tableLambdaSkewT;
bool logX=true, logT=false, logC=true;
//
// Set filenames for the tables.
// Be as desciptive as possible. Helps when mass producing tables.
//
// We create all tables and decide later what
// gets written and what not.
//
string rootfile=settings->rootfile();
rootfile += "_"; rootfile += settings->dipoleModelName();
rootfile += "_"; rootfile += settings->dipoleModelParameterSetName();
rootfile += "_A"; rootfile += to_string(settings->A());
rootfile += "_UPC";
ostringstream filenameT, filenameT2, filenameVarT, filenameLambdaTReal, filenameLambdaTSkew;
filenameT.str("");
filenameT << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin <<"_T.root";
filenameT2.str("");
filenameT2 << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_T2.root";
filenameVarT.str("");
filenameVarT << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_VarT.root";
filenameLambdaTReal.str("");
filenameLambdaTReal << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_LambdaTReal.root";
filenameLambdaTSkew.str("");
filenameLambdaTSkew << rootfile << "_" << settings->vectorMesonId() << "_bin"
<< startingBin << "-" << endingBin << "_LambdaTSkew.root";
(void) tableT.create(nBinX, xmin, xmax,
nBinT, tmin, tmax,
logX, logT, logC, // all bools
mean_A, transverse,
massA, vmPDG, model, pset,
filenameT.str().c_str(),
priority);
(void) tableT2.create(nBinX, xmin, xmax,
nBinT, tmin, tmax,
logX, logT, logC, // all bools
mean_A2, transverse,
massA, vmPDG, model, pset,
filenameT2.str().c_str(),
priority);
bool contentVar = modes == 1 ? false : logC; // if mode = 1 they will be 0 so lin is better
(void) tableVarT.create(nBinX, xmin, xmax,
nBinT, tmin, tmax,
logX, logT, contentVar, // all bools
variance_A, transverse,
massA, vmPDG, model, pset,
filenameVarT.str().c_str(),
priority);
double contentLambda = false; // content lin
if(createLambdaTables){
(void) tableLambdaRealT.create(nBinX, xmin, xmax,
nBinT, tmin, tmax,
logX, logT, contentLambda, // all bools
lambda_A, transverse,
massA, vmPDG, model, pset,
filenameLambdaTReal.str().c_str(),
priority);
(void) tableLambdaSkewT.create(nBinX, xmin, xmax,
nBinT, tmin, tmax,
logX, logT, contentLambda, // all bools
lambda_skew, transverse,
massA, vmPDG, model, pset,
filenameLambdaTSkew.str().c_str(),
priority);
}
cout << "\nAll tables created:" << endl;
tableT.list();
tableT2.list();
tableVarT.list();
if(createLambdaTables){
tableLambdaRealT.list();
tableLambdaSkewT.list();
}
cout << "\nTables have " << maxbins << " bins each.\n" << endl;
if (settings->useBackupFile()) {
int ibin = settings->startingBinFromBackup();
tableT.setAutobackup("tableT", ibin);
tableT2.setAutobackup("tableT2", ibin);
tableVarT.setAutobackup("tableVarT", ibin);
tableLambdaRealT.setAutobackup("tableLambdaRealT", ibin);
tableLambdaSkewT.setAutobackup("tableLambdaSkewT", ibin);
cout << "Automatic backup of tables is enabled." << endl;
}
else
cout << "Automatic backup of tables is off.\n" << endl;
//
// DGLAP Evolution can be speed up by using lookup tables
//
DglapEvolution &dglap = DglapEvolution::instance();
// dglap.generateLookupTable(100, 100);
dglap.generateLookupTable(1000, 1000);
dglap.useLookupTable(true);
//
// Create and initialize the amplitudes calculator
//
Amplitudes amps;
//
// Generate the the nucleon configurations
//
amps.generateConfigurations();
//
// Print out settings
//
cout << endl;
cout << "Tables will be generated for:" << endl;
cout << "\tNucleus mass A="<<massA<<endl;
cout << "\tModes to calculate: " << modes << endl;
cout << "\tVector Meson Id: " << vmPDG << endl;
cout << "\tBins: " << startingBin << "-" << endingBin << endl;
cout << "\txpom range: [" << xmin << ", " << xmax << "], " << nBinX << " bins." << endl;
cout << "\t t range: [" << tmin << ", " << tmax << "], " << nBinT << " bins." << endl;
cout << "\tDipole model: " << settings->dipoleModelName() << endl;
cout << "\tDipole model parameter set: " << settings->dipoleModelParameterSetName() << endl;
cout << "\tTable set mode: " << settings->tableSetTypeName() << endl;
cout << endl;
//
// Calculate contents and fill tables
//
// Note that we fill all tables. What is written
// at the end is another story.
//
cout << "Start filling tables" << endl;
time_t tableFillStart = time(0);
int nShow = (endingBin - startingBin)/100;
if(nShow==0) nShow=1;
for (int i=startingBin; i<=endingBin; i++) {
if (i%nShow == 0 || i == startingBin || i == endingBin)
cout << "processing bin " << i << endl;
double xpom, t;
tableT.binCenter(i, xpom, t);
double kinematicPoint[2]={t, xpom};
amps.calculate(kinematicPoint);
double aT = 0;
double aT2 = 0;
double aVarT = 0;
//
// For mode = 1 only coherent is
// calculated and incoherent is therefore
// assumed to be 0.
// Hence total == coherent
//
// coherent = A
// total = A2
// incoherent = variance = A2 - A*A
//
aT = amps.amplitudeT();
tableT.fill(i, aT);
if (modes != 1) {
aT2 = amps.amplitudeT2();
}
else {
aT2 = aT*aT;
}
tableT2.fill(i, aT2);
aVarT = aT2-aT*aT;
tableVarT.fill(i, aVarT);
//
// Calculate lambda and fill table lambda=dlog(A)/dlog(1/x)
//
if (createLambdaTables) {
double hplus, hminus;
hplus=hminus=xpom*1e-4; //This value needs to be tested
double lambdaReal=0;
double lambdaSkew=0;
kinematicPoint[1]=xpom+hplus;
amps.calculate(kinematicPoint);
double ampPlus=amps.amplitudeT();
double ampSkewPlus=amps.amplitudeTForSkewednessCorrection();
kinematicPoint[1]=xpom-hminus;
amps.calculate(kinematicPoint);
double ampMinus=amps.amplitudeT();
double ampSkewMinus=amps.amplitudeTForSkewednessCorrection();
// Don't calculate unless numerically viable:
if (ampPlus==0 || ampMinus==0) {
lambdaReal=0;
}
else {
//Calculate derivate d(logA/dxpom)=log(A+/A-)/(h+-h-):
double derivate=log(abs(ampPlus/ampMinus))/(hplus+hminus);
// Finally calculate lambdas:
double jacobian = -xpom;
lambdaReal = jacobian*derivate;
}
// Don't calculate unless numerically viable:
if (ampSkewPlus==0 || ampSkewMinus==0) {
lambdaSkew=0;
}
else {
//Calculate derivate d(logA/dxpom)=log(A+/A-)/(h+-h-):
double derivate=log(abs(ampSkewPlus/ampSkewMinus))/(hplus+hminus);
// Finally calculate lambdas:
double jacobian = -xpom;
lambdaSkew = jacobian*derivate;
}
tableLambdaRealT.fill(i, lambdaReal);
tableLambdaSkewT.fill(i, lambdaSkew);
}
}
time_t tableFillEnd = time(0);
//
// Report CPU time/cell
//
cout << endl;
cout << "CPU time/bin: "
<< static_cast<double>(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
<< " s" << endl;
cout << "Total time: " << static_cast<double>(tableFillEnd-tableFillStart)/60./60. << " h" << endl << endl;
//
// We write out all tables.
//
// Whoever runs the production can then decide
// later what to keep and what to delete.
// If the desired run mode is total_and_coherent or
// coherent_and_incoherent, if this is just to improve
// a coherent table in some phase space, or if one wants
// to maintain redundancy, all these factor might affect
// your choice.
//
tableT.write();
tableT2.write();
tableVarT.write();
if (createLambdaTables) {
tableLambdaRealT.write();
tableLambdaSkewT.write();
}
cout << "All tables written" << endl;
cout << "All done. Bye." << endl;
return 0;
}
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 <http://www.gnu.org/licenses/>.
//
// 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 <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include <iostream>
#include <cmath>
#include <vector>
#include "TVector3.h"
#include <algorithm>
#include <unistd.h>
#include <cstdlib>
#include "Nucleus.h"
#include "TableGeneratorNucleus.h"
#include <glut.h>
#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; i<myNucleus->configuration.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; i<str.length(); i++) {
if(!isdigit(str[i])) {
cout << "Error, argument must be a digit." << endl;
return 2;
}
}
unsigned int A = atoi(argv[optind]);
myNucleus = new TableGeneratorNucleus(A);
while(!myNucleus->generate()){};
cout<<"List of Commands:"<<endl;
cout<<"\trotate view: left/right up/down arrow keys"<<endl;
cout<<"\tzoom in/out: a/z"<<endl;
cout<<"\tregenerate nucleus: n"<<endl;
cout<<"\tquit: q/esc"<<endl;
// init GLUT and create window
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
glutInitWindowPosition(100,100);
glutInitWindowSize(320,320);
glutCreateWindow("Sartre Nucleus");
glClearColor(.0f, .0f, .0f, 0.0f);
//define colours:
//light
const GLfloat white[]= { 1.0, 1.0, 1.0, 1.0 };
const GLfloat magenta[]= { 1.0, 0.0, 1.0, 1.0 };
const GLfloat gold[] = { 1.0, 1.0, 0.0, 1.0 };
const GLfloat red[] = { 1.0, 0.0, 0.0, 1.0 };
const GLfloat green[] = { 0.0, 1.0, 0.0, 1.0 };
const GLfloat blue[] = { 0.0, 0.0, 1.0, 1.0 };
// const GLfloat color[];
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
switch(A){
case 208:
case 16:
glLightfv(GL_LIGHT0, GL_DIFFUSE, blue);
glLightfv(GL_LIGHT0, GL_SPECULAR, blue);
break;
case 197:
glLightfv(GL_LIGHT0, GL_DIFFUSE, gold);
glLightfv(GL_LIGHT0, GL_SPECULAR, gold);
break;
case 63:
glLightfv(GL_LIGHT0, GL_DIFFUSE, red);
glLightfv(GL_LIGHT0, GL_SPECULAR, red);
break;
case 27:
case 40:
glLightfv(GL_LIGHT0, GL_DIFFUSE, white);
glLightfv(GL_LIGHT0, GL_SPECULAR, white);
break;
case 110:
glLightfv(GL_LIGHT0, GL_DIFFUSE, magenta);
glLightfv(GL_LIGHT0, GL_SPECULAR, magenta);
break;
case 1:
glLightfv(GL_LIGHT0, GL_DIFFUSE, green);
glLightfv(GL_LIGHT0, GL_SPECULAR, green);
break;
}
if (hasWhiteBackground)
glClearColor(1.0f, 1.0f, 1.0f, 0.0f);
// register callbacks
glutDisplayFunc(renderScene);
glutReshapeFunc(changeSize);
glutIdleFunc(idle);
glutSpecialFunc(pressKey);
// here are the new entries
glutIgnoreKeyRepeat(1);
glutSpecialUpFunc(releaseKey);
// OpenGL init
glEnable(GL_DEPTH_TEST);
// enter GLUT event processing cycle
glutMainLoop();
return 0;
}
Index: trunk/examples/tableGeneratorUPCRuncard.txt
===================================================================
--- trunk/examples/tableGeneratorUPCRuncard.txt (revision 370)
+++ trunk/examples/tableGeneratorUPCRuncard.txt (revision 371)
@@ -1,115 +1,115 @@
//==============================================================================
// 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 <http://www.gnu.org/licenses/>.
//
// 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 <A> analytical, <A2> averaged over configurations
# 1 <A> only, analytical
# 2 <A2>, and <A> 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 <http://www.gnu.org/licenses/>.
//
// 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 <http://www.gnu.org/licenses/>.
//
// 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 <iostream>
#include <cmath>
#include <vector>
#include <string>
#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<HeraDataPoint>&);
void loadPhiData(vector<HeraDataPoint>&);
void loadDVCSData(vector<HeraDataPoint>&, vector<HeraDataPoint>&);
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<HeraDataPoint> zeusDataPhi_table1;
loadPhiData(zeusDataPhi_table1);
vector<TGraphErrors*> graphVec;
vector<TGraphErrors*> sgraphVec;
vector<string> sgraphName;
vector<string> 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; i<zeusDataPhi_table1.size()-1 ;i++) {
// for (unsigned int i=0; i<7+6+6+6 ;i++) {
// fill graphs
pointCounter = i;
if(pointCounter<7+6+6+6+9)
graphCounter = 4;
if(pointCounter<7+6+6+6)
graphCounter = 3;
if(pointCounter<7+6+6)
graphCounter = 2;
if(pointCounter<7+6)
graphCounter = 1;
if(pointCounter<7)
graphCounter = 0;
if(graphCounter<4)
graphVec[graphCounter]->SetPoint(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<pair<double, double> > 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 << "+/-" <<zeusDataPhi_table1[i].xsection_sys << endl;
cout << "\t\tSartre/Zeus = " << CS/zeusDataPhi_table1[i].xsection << endl;
}
cout << "-------------------------------------------------" << endl;
cout << "Comparison HERA-Sartre:" << endl;
cout << "chi2 = " << chi2_table1 << endl;
cout << "ndof = " << ndof_table1 << endl;
cout << "chi2/ndof = " << chi2_table1/ndof_table1 << endl;
cout << "-------------------------------------------------" << endl;
cout<<endl;
cout << "-------------------------------------------------" << endl;
cout << "Comparison HERA-Sartre:" << endl;
cout << "chi2 = " << chi2_table3 << endl;
cout << "ndof = " << ndof_table3 << endl;
cout << "chi2/ndof = " << chi2_table3/ndof_table3 << endl;
cout << "-------------------------------------------------" << endl;
//
// ROOT file
//
TFile *hfile = 0;
hfile = new TFile(rootfile.c_str(),"RECREATE");
cout << "ROOT file is '" << rootfile.c_str() << "'." << endl;
for (unsigned int i=0; i<graphVec.size(); i++) {
graphVec[i]->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<HeraDataPoint> h1DataDVCS_table4;
vector<HeraDataPoint> h1DataDVCS_table1;
loadDVCSData(h1DataDVCS_table1, h1DataDVCS_table4);
vector<TGraphErrors*> graphVec;
vector<TGraphErrors*> sgraphVec;
vector<string> sgraphName;
vector<string> 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; i<h1DataDVCS_table4.size() ;i++) {
// fill graphs
graphCounter = i/4;
pointCounter = i-graphCounter*4;
graphVec[graphCounter]->SetPoint(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<pair<double, double> > 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 << "+/-" <<h1DataDVCS_table4[i].xsection_sys << endl;
cout << "\t\tSartre/H1 = " << CS/h1DataDVCS_table4[i].xsection << endl;
}
//
// Table 1
//
double sigma_a = 0;
double sigma_b = 0;
double h1_sigma_a = 0;
double h1_sigma_b = 0;
double diff_table1, chi2_table1 = 0;
int ndof_table1 = 0;
for (unsigned int i=0; i<h1DataDVCS_table1.size() ;i++) {
// fill graphs
double error = sqrt(h1DataDVCS_table1[i].xsection_stat*h1DataDVCS_table1[i].xsection_stat +
h1DataDVCS_table1[i].xsection_sys*h1DataDVCS_table1[i].xsection_sys);
if (i < 4) {
graph1a->SetPoint(i, fabs(h1DataDVCS_table1[i].Q2), h1DataDVCS_table1[i].xsection);
graph1a->SetPointError(i , 0., error);
}
else {
graph1b->SetPoint(i-4, fabs(h1DataDVCS_table1[i].W), h1DataDVCS_table1[i].xsection);
graph1b->SetPointError(i-4 , 0., error);
}
// set up sartre
settings->setVerbose(true);
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_table1[i].Q2-.1);
settings->setQ2max(h1DataDVCS_table1[i].Q2+.1);
settings->setWmin(h1DataDVCS_table1[i].W-.1);
settings->setWmax(h1DataDVCS_table1[i].W+.1);
settings->setMaxLambdaUsedInCorrections(maxLam);
if(settings->Wmax() > 142 || settings->Wmin()>142)
continue;
bool ok = sartre.init();
if (!ok) {
cout << "Initialization of sartre failed. Going to next data point." << endl;
continue;
}
//
// Sartre cross-section
//
vector<pair<double, double> > 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 << "+/-" <<h1DataDVCS_table1[i].xsection_sys << endl;
cout << "\t\tSartre/H1 = " << CS/h1DataDVCS_table1[i].xsection << endl;
}
cout << "final:\n";
PR(sigma_a);
PR(sigma_b);
PR(h1_sigma_a);
PR(h1_sigma_b);
cout << "-------------------------------------------------" << endl;
cout << "Comparison HERA-Sartre (table 1):" << endl;
cout << "chi2 = " << chi2_table1 << endl;
cout << "ndof = " << ndof_table1 << endl;
cout << "chi2/ndof = " << chi2_table1/ndof_table1 << endl;
cout << endl;
cout << "Comparison HERA-Sartre (table 4):" << endl;
cout << "chi2 = " << chi2_table4 << endl;
cout << "ndof = " << ndof_table4 << endl;
cout << "chi2/ndof = " << chi2_table4/ndof_table4 << endl;
cout << "-------------------------------------------------" << endl;
//
// ROOT file
//
TFile *hfile = 0;
hfile = new TFile(rootfile.c_str(),"RECREATE");
cout << "ROOT file is '" << rootfile.c_str() << "'." << endl;
// Table 4
for (unsigned int i=0; i<graphVec.size(); i++) {
graphVec[i]->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<HeraDataPoint> 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; i<zeusDataJpsi.size() ;i++) {
settings->setVerbose(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 << "+" <<zeusDataJpsi[i].xsection_sys_up << "-" << zeusDataJpsi[i].xsection_sys_lo << endl;
cout << "\t\tRatio: " << 1000*sartre.totalCrossSection()/zeusDataJpsi[i].xsection << endl;
TGraphAsymmErrors *graphS, *graphH;
TGraphAsymmErrors *graphSS, *graphHH;
graphS = graphSS = graphH = graphHH = 0;
double offset = 0;
if (i <= 2) {
graphS = graph1S;
graphH = graph1H;
}
else if (i >= 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<HeraDataPoint>& 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<HeraDataPoint> 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<HeraDataPoint>& vec1, vector<HeraDataPoint>& 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<HeraDataPoint> h1DataDVCS_table1(4+5);
vector<HeraDataPoint> 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<HeraDataPoint>& vec1)
{
vector<HeraDataPoint> 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;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:52 PM (9 h, 31 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023705
Default Alt Text
(189 KB)

Event Timeline