Page MenuHomeHEPForge
Files F8309114

sartreTest.cpp
No OneTemporary
Actions

sartreTest.cpp
View Options

//==============================================================================
// sartreTest.cpp
//
// Copyright (C) 2010-2013 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: 2017-12-20 17:08:30 +0000 (Wed, 20 Dec 2017) $
// $Author: ullrich $
//==============================================================================
//
// Developer test program. Checks individual components.
// Should not be in the distribution.
//
//==============================================================================
#include "EventGeneratorSettings.h"
#include "Table.h"
#include <iostream>
#include "TParticlePDG.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 "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
#if defined(TEST_INTERPOLATE)
struct InterTestPoint {
double t, Q2, W2, value;
};
#endif
#if defined(TEST_SLOPE)
struct InterTestSlope {
double t, Q2, W2, slope;
};
#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
//
// 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)
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;
}

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 2:22 PM (9 h, 19 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023097
Default Alt Text
sartreTest.cpp (16 KB)

Event Timeline