|
Home · Overview · Users Guide · Reference Guide |
The Sartre distribution contains a set of example programs in the sartre/examples directory. The probable most important one is sartreMain, a short but rather complete example program for event generation of exclusive vector mesons and DVCS in ep or eA collisions. It consist of two files, the main program in sartreMain.cpp and a runcard sartreRuncard.txt. The program sets up Sartre, generates events, and stores them in a self-defined ROOT tree that is written to file for later analysis. Here we go in detail through this example, show what it does and how the output looks like. It is always easier to start from a running example and modify it for your purposes than to start from scratch. Note that the actual program contained in the distribution might look slightly different.
TTree.h and TFile.h files.
For Sartre all we need to include is Sartre.h.
The structure rootSartreEvent is just an example on how
to store the event generator output into a ROOT tree.
#include <iostream>
#include <cmath>
#include "TTree.h"
#include "TFile.h"
#include "Sartre.h"
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;
The example program sartreMain takes 2 arguments. The first is the file name of the runcard and is mandatory the second is the rootfile which is optional.
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];
To initialize the generator all that is needed is to create an instance of class Sartre and initialize it using the parameter defined in the runcard. We also obtain a pointer to EventGeneratorSettings so we can make use of the parameters in the runcard such as file names, number of events, and many more. It's always good to list all setup parameters (settings->list()) for checks and record keeping.
//
// 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();
settings->list();
Now we are in principle ready to generate events. However we want to store these events for later analysis in a ROOT tree and save it to file. This is happening next.
In the ROOT tree we keep the general event parameters in the rootSartreEvent structure and the 4-momenta of in- and outgoing electron and proton/nucleus, as well as the vector meson and the virtual photon in different branches. The tree is stored in a ROOT file which name is either coming from the runcard or (if present) from the command line. Reading in this file in ROOT allows one to create directly plots, apply cuts, or just filter the events for further processing.
//
// ROOT file
//
string rootfile;
if (argc == 3)
rootfile = argv[2];
else
rootfile = settings->rootfile();
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;
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);
Now we generate the events. The number of events (EventGeneratorSettings::numberOfEvents()) and the frequency with which to print out the status (EventGeneratorSettings::timesToShow()) we get from the runcard. These are user parameters provided for convenience and not directly used by Sartre. The user can safely ignore them and hardwire the numbers instead.
Events are generated via Sartre::generateEvent(), which returns the event record. If you want to filter certain events that would be the place to do it. After that all what is left is to store the obtained information into the ROOT tree.
//
// Event generation
//
int nPrint;
if (settings->timesToShow())
nPrint = settings->numberOfEvents()/settings->timesToShow();
else
nPrint = 0;
unsigned long maxEvents = settings->numberOfEvents();
cout << "Generating " << maxEvents << " events." << endl << endl;
for (unsigned long iEvent = 0; iEvent < maxEvents; iEvent++) {
//
// Generate one event
//
Event *event = sartre.generateEvent();
if (nPrint && iEvent%nPrint == 0 && iEvent != 0) {
cout << "processed " << iEvent << " events" << endl;
}
//
// Print out (here only for the first few events)
//
if (iEvent < 4) 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;
tree.Fill();
}
Now we are essentially done. One important thing left to do is to calculate the total cross-section in the kinematic range the events were generated. Without knowledge of that number you will not be able to absolutely normalize the generated output. This is done using Sartre::TotalCrossSection(). This might actually take a while depending on the kinematic range chosen, especially in eA because of the complex structure along t. Note that the returned cross-section is in "nb". Just for bookkeeping purposes we also print the time used per events.
//
// That's it, finish up
//
double runTime = sartre.runTime();
hfile->Write();
cout << rootfile.c_str() << " written." << endl;
cout << "Total cross-section: " << sartre.totalCrossSection() << " nb" << endl;
sartre.listStatus();
cout << "CPU Time/event: " << 1000*runTime/maxEvents << " msec/evt" << endl;
return 0;
}
Here we assume you copied the sartre/examples directory from the installation directory
into you private area.
Make sure that the environment variable SARTRE_DIR is
defined and points to the installation directory. For example:
# example using bash
export SARTRE_DIR=/usr/local/sartre
Then within the example directory:
cmake .
make sartreMain
Your example generator is now ready for use.
Below is the runcard we are using for a first test run (examples/sartreRuncard.txt).
There is a complete description of the runcard on the runcard reference page. We
skip here all comment lines, that is those lines starting with a "#" or
"//" and focus only on the definition of the parameters.
eBeamEnergy = 20 hBeamEnergy = 100 A = 197We will generate only 1000 events and print the status (how many generated so far) 20 times.
numberOfEvents = 10000 timesToShow = 20We switch the verbosity on but keep it at the default level (minimal amount of additional messages).
verbose = true verboseLevel = 1Events are written into a root tree stored in file
example.root.
rootfile = example.rootWe generate J/ψ mesons using the bSat model (saturation).
vectorMesonId = 443 dipoleModel = bSatHere we set the kinematic limits. Note that when setting Q2min > Q2max the maximum available kinematic range is automatically selected. Same for Wmin and Wmax.
Q2min = 1 Q2max = 10 Wmin = 10 Wmax = 60We switch all corrections on to get the most reliable cross-section.
correctForRealAmplitude = true correctSkewedness = trueWe enable the nuclear breakup of the Au ion for incoherent events but limit the maximum excitation energy to 0.5 GeV. This affects the result very little but reduces the time Gemini (the used evaporation model) takes to generate the breakup products. When producing a large number of events, the breakup should probably be switched off.
enableNuclearBreakup = true maxNuclearExcitationEnergy = 0.5
Launch the program sarteMain with the name of the runcard as argument:
sartreMain sartreRuncard.txt
This will produce the following printout
/========================================================================\
| |
| Sartre, Version 1.10 |
| |
| An event generator for exclusive diffractive vector meson production |
| in ep and eA collisions based on the dipole model. |
| |
| Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich |
| |
| 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, either version 3 of the License, or |
| any later version. |
| |
| Code compiled on Dec 14 2012 12:36:03 |
| Run started at Fri Dec 14 12:38:09 2012 |
\========================================================================/
Runcard is 'sartreRuncard.txt'.
Hadron beam species: Au (197)
Hadron beam: 0 0 99.9956 100 (0.93827)
Electron beam: 0 0 -20 20 (0.000510999)
Dipole model: bSat
Process is e + Au -> e' + Au' + J/psi
Random generator seed: 1355506689
Kinematic limits used for event generation:
t=[-0.5, -1.2669e-06]
Q2=[1, 10]
W=[10, 60]
Finding mode of pdf:
location: t=-0.000693239, Q2=1, W=19.5609; value: 47.6573
Sartre is initialized.
Run Settings:
userInt 0
userDouble 0
userString
eBeamEnergy 20
hBeamEnergy 100
A 197
numberOfEvents 1000
timesToShow 20
vectorMesonId 443
dipoleModel bSat
Q2min 1
Q2max 10
Wmin 10
Wmax 60
correctForRealAmplitude true
correctSkewedness true
enableNuclearBreakup true
maxLambdaUsedInCorrections 0.65
maxNuclearExcitationEnergy 0.5
applyPhotonFlux true
verbose true
verboseLevel 1
rootfile example.root
seed 1355506689
ROOT file is 'example.root'.
Generating 1000 events.
evt = 0 Q2 = 5.475 x = 1.242e-02
W = 20.883 y = 0.055
t = -0.004 xpom = 3.420e-02
pol = L diff = coherent
# id name status parents daughters px py pz E m
0 11 e- 4 - - 2 3 0.000 0.000 -20.000 20.000 5.110e-04
1 1000791970 Au(197) 4 - - 6 - 0.000 0.000 99.996 100.000 0.938
2 11 e- 1 0 - - - -1.938 1.191 -18.830 18.967 5.110e-04
3 22 gamma 2 0 - 4 5 1.938 -1.191 -1.170 1.033 -2.340
4 443 J/psi 1 3 - - - 1.990 -1.212 2.306 4.510 3.097
5 990 pomeron 2 3 3 6 - -0.053 0.020 -3.477 -3.476 -0.066
6 1000791970 Au(197) 1 1 5 - - -0.053 0.020 96.519 96.524 0.938
evt = 1 Q2 = 2.150 x = 4.916e-03
W = 20.883 y = 0.055
t = -0.004 xpom = 2.685e-02
pol = T diff = coherent
# id name status parents daughters px py pz E m
0 11 e- 4 - - 2 3 0.000 0.000 -20.000 20.000 5.110e-04
1 1000791970 Au(197) 4 - - 6 - 0.000 0.000 99.996 100.000 0.938
2 11 e- 1 0 - - - 1.399 -0.276 -18.880 18.933 5.110e-04
3 22 gamma 2 0 - 4 5 -1.399 0.276 -1.120 1.067 -1.466
4 443 J/psi 1 3 - - - -1.417 0.218 1.569 3.756 3.097
5 990 pomeron 2 3 3 6 - 0.018 0.057 -2.690 -2.690 -0.066
6 1000791970 Au(197) 1 1 5 - - 0.018 0.057 97.306 97.310 0.938
evt = 2 Q2 = 2.150 x = 2.961e-03
W = 26.926 y = 0.091
t = -0.004 xpom = 1.617e-02
pol = T diff = coherent
# id name status parents daughters px py pz E m
0 11 e- 4 - - 2 3 0.000 0.000 -20.000 20.000 5.110e-04
1 1000791970 Au(197) 4 - - 6 - 0.000 0.000 99.996 100.000 0.938
2 11 e- 1 0 - - - -1.238 -0.650 -18.157 18.211 5.110e-04
3 22 gamma 2 0 - 4 5 1.238 0.650 -1.843 1.789 -1.466
4 443 J/psi 1 3 - - - 1.301 0.662 -0.202 3.429 3.097
5 990 pomeron 2 3 3 6 - -0.063 -0.012 -1.641 -1.641 -0.066
6 1000791970 Au(197) 1 1 5 - - -0.063 -0.012 98.355 98.359 0.938
evt = 3 Q2 = 1.683 x = 2.731e-03
W = 24.806 y = 0.077
t = -0.003 xpom = 1.830e-02
pol = T diff = coherent
# id name status parents daughters px py pz E m
0 11 e- 4 - - 2 3 0.000 0.000 -20.000 20.000 5.110e-04
1 1000791970 Au(197) 4 - - 6 - 0.000 0.000 99.996 100.000 0.938
2 11 e- 1 0 - - - -1.217 0.270 -18.439 18.481 5.110e-04
3 22 gamma 2 0 - 4 5 1.217 -0.270 -1.561 1.519 -1.297
4 443 J/psi 1 3 - - - 1.200 -0.219 0.258 3.339 3.097
5 990 pomeron 2 3 3 6 - 0.016 -0.051 -1.819 -1.819 -0.057
6 1000791970 Au(197) 1 1 5 - - 0.016 -0.051 98.176 98.181 0.938
processed 50 events
processed 100 events
processed 150 events
processed 200 events
processed 250 events
processed 300 events
processed 350 events
processed 400 events
processed 450 events
processed 500 events
processed 550 events
processed 600 events
processed 650 events
processed 700 events
processed 750 events
processed 800 events
processed 850 events
processed 900 events
processed 950 events
example.root written.
Total cross-section: 98.8 nb
Event summary: 1000 events generated, 1000 tried
Total time used: 1 min 0 sec
CPU Time/event: 10 msec/evt
It is easiest to start with this example and modify it for your purposes.
Last Update: May 29, 2013 |