No OneTemporary
Actions

Size

38 KB

Subscribers

None

View Options

	Index: trunk/examples/nucleusAnimation.cpp
	===================================================================
	--- trunk/examples/nucleusAnimation.cpp (revision 245)
	+++ trunk/examples/nucleusAnimation.cpp (revision 246)
	@@ -1,300 +1,300 @@
	//==============================================================================
	// nucleusAnimation.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: 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 <GL/glut.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= (clrx + sudslry + cudslrz)zoom;
	y= ( cudy + (-sud)z)*zoom;
	z=(-slrx + sudclry + cudclrz)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/sartreTest.cpp
	===================================================================
	--- trunk/examples/sartreTest.cpp (revision 245)
	+++ trunk/examples/sartreTest.cpp (revision 246)
	@@ -1,471 +1,481 @@
	//==============================================================================
	// 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$
	// $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 "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_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
	- //
	-#if defined (TEST_WRITETABLE)
	- Table table1, table2, table3;
	- int n = table1.create(10, 1, 100, 10, 1, 80, 10, 0.00001, 2,
	- true, true, false, false,
	+ //
	+ /*
	+ 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.00001, 2,
	+ 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.00001, 2,
	+ 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.50.5+0.30.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./(radiusradiusbinsize)/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/sartreMain.cpp
	===================================================================
	--- trunk/examples/sartreMain.cpp (revision 245)
	+++ trunk/examples/sartreMain.cpp (revision 246)
	@@ -1,289 +1,289 @@
	//==============================================================================
	// sartreMain.cpp
	//
	// Copyright (C) 2010-2015 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 "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;

	// 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);
	}
	}


	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;
	}
	Index: trunk/examples/sartreRuncard.txt
	===================================================================
	--- trunk/examples/sartreRuncard.txt (revision 245)
	+++ trunk/examples/sartreRuncard.txt (revision 246)
	@@ -1,122 +1,122 @@
	//==============================================================================
	// sartreRuncard.txt
	//
	// Copyright (C) 2010-2016 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 = 10
	+eBeamEnergy = 20
	hBeamEnergy = 100
	A = 197

	#
	# UPC settings, to run in UPC mode set UPC=true and UPCA into the photon emitting species
	#
	UPC=false
	UPCA=197

	#
	# Number of events and printout frequency
	#
	-numberOfEvents = 100000
	+numberOfEvents = 10000
	timesToShow = 20

	#
	# Set verbosity
	#
	verbose = true
	-verboseLevel = 1
	+verboseLevel = 2

	#
	# Rootfile
	#
	rootfile = example.root

	#
	# Dipole model
	#
	dipoleModel = bSat
	dipoleModelParameterSet = KMW

	#
	# Model parameters
	#
	# vectorMesonID: 22, 113, 333, 443
	#
	vectorMesonId = 443

	#
	# Table Set Type (experts only)
	#
	tableSetType = total_and_coherent
	#tableSetType = coherent_and_incoherent

	#
	# User variable used for vector meson decays
	# PDG: pi+ = 211, K+ = 321, e- = 11, mu- = 13
	#
	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

	#
	# Misc
	#
	-enableNuclearBreakup = true
	+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

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 1:59 PM (2 h, 24 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486419
Default Alt Text: (38 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions