Page MenuHomeHEPForge
Files F7879379

No OneTemporary
Actions

View Options

Index: trunk/src/Sartre.cpp
===================================================================
--- trunk/src/Sartre.cpp (revision 168)
+++ trunk/src/Sartre.cpp (revision 169)
@@ -1,746 +1,746 @@
//==============================================================================
// Sartre.cpp
//
-// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
+// Copyright (C) 2010-2015 Tobias Toll and Thomas Ullrich
//
-// This file is part of Sartre version: 1.1
+// This file is part of Sartre version: 1.2
//
// 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$
//==============================================================================
//
// Note:
// When not using runcards, the user must first create an instance of Sartre
// then get the settings via one of:
// Sartre::runSettings()
// EventGeneratorSettings::instance()
// Once init() is called settings cannot be changed any more.
//==============================================================================
#include "Version.h"
#include "Kinematics.h"
#include "Sartre.h"
#include "ModeFinderFunctor.h"
#include "Math/BrentMinimizer1D.h"
#include "Math/IntegratorMultiDim.h"
#include "TUnuranMultiContDist.h"
#include <limits>
#include <cmath>
#include <iomanip>
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
Sartre::Sartre()
{
mIsInitialized = false;
mCurrentEvent = 0;
mNucleus = 0;
mSettings = 0;
mPDF_Functor = 0;
mPDF = 0;
mEventCounter = 0;
mTriesCounter = 0;
mTotalCrossSection = 0;
mCrossSection = 0;
mTableCollection = 0;
mProtonTableCollection = 0;
mUnuran = 0;
mEvents = 0;
mTries = 0;
mS = 0;
mA = 0;
}
Sartre::~Sartre()
{
delete mNucleus;
delete mPDF_Functor;
delete mPDF;
delete mCrossSection;
if (mTableCollection != mProtonTableCollection) {
delete mTableCollection;
delete mProtonTableCollection;
}
else
delete mTableCollection;
delete mUnuran;
}
bool Sartre::init(const char* runcard)
{
mStartTime = time(0);
bool ok;
//
// Reset member variables.
// Note that one instance of Sartre should be able to get
// initialized multiple times.
//
mEvents = 0;
mTries = 0;
mTotalCrossSection = 0;
//
// Print header
//
string ctstr(ctime(&mStartTime));
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-2013 Tobias Toll and Thomas Ullrich |" << endl;
+ cout << "| Copyright (C) 2010-2015 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;
mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton
//
// Read runcard if available
//
if (runcard) {
if (!mSettings->readSettingsFromFile(runcard)) {
cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl;
exit(1);
}
else
cout << "Runcard is '" << runcard << "'." << endl;
}
else
cout << "No runcard provided." << endl;
//
// Set beam particles and center of mass energy
//
mElectronBeam = mSettings->eBeam();
mHadronBeam = mSettings->hBeam();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mA = mSettings->A();
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (allowBreakup) {
if (!getenv("SARTRE_DIR")) {
cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n"
"to locate tables needed for the generation if nuclear breakups." << endl;
exit(1);
}
}
if (mNucleus) delete mNucleus;
mNucleus = new FrangibleNucleus(mA, allowBreakup);
cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron beam: " << mHadronBeam << endl;
cout << "Electron beam: " << mElectronBeam << endl;
//
// Get details about the processes and models
//
mDipoleModelType = mSettings->dipoleModelType();
mVmID = mSettings->vectorMesonId();
cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
cout << "Process is ";
if (mA > 1)
cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + ";
else
cout << "e + p -> e' + p' + ";
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID);
cout << vectorMesonPDG->GetName() << endl;
//
// Print-out seed for reference
//
cout << "Random generator seed: " << mSettings->seed() << endl;
//
// Load in the tables containing the amplitude moments
//
if (!getenv("SARTRE_DIR")) {
cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
exit(1);
}
if (mTableCollection) delete mTableCollection;
mTableCollection = new TableCollection;
ok = mTableCollection->init(mA, mDipoleModelType, mVmID);
if (!ok) {
cout << "Error, could not initialize lookup tables for requested process." << endl;
return false;
}
//
// Load in the p tables for the lambda lookup tables (or to calculate lambda if not available)
//
if (mSettings->correctForRealAmplitude() || mSettings->correctSkewedness()) {
if (mA == 1) {
mProtonTableCollection = mTableCollection;
}
else {
if (mProtonTableCollection) delete mProtonTableCollection;
mProtonTableCollection = new TableCollection;
ok = mProtonTableCollection->init(1, mDipoleModelType, mVmID);
if (!ok) {
cout << "Error: could not initialize proton lookup tables for requested process." << endl;
cout << " These tables are needed for skewedness and real amplitude corrections." << endl;
return false;
}
}
}
else
mProtonTableCollection = 0;
//
// Kinematic limits and generator range
//
// There are 3 ranges we have to deal with
// 1. the kinematic range requested by the user
// if given.
// The user can only control Q2 and W but not t.
// 2. the range of the table(s)
// 3. the kinematically allowed range
//
// Of course (3) is more complex than a simple cube.
// However, we deal with the detailed shape of the kinematic
// range later using Kinematics::valid() when we generate the
// individual events.
// For setting up UNU.RAN we have to get the cubic envelope
// that satifies (1)-(3).
// Note, that they are correlated which makes the order
// in which we do things a bit tricky.
//
//
// Step 1:
// Set the limits to that of the table(s).
// Note, the indices 0-2 refer to t, Q2, and W2.
//
mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT();
mLowerLimit[1] = mTableCollection->minQ2(); mUpperLimit[1] = mTableCollection->maxQ2();
mLowerLimit[2] = mTableCollection->minW2(); mUpperLimit[2] = mTableCollection->maxW2();
//
// Step 2:
// Kinematic limits might overrule boundaries from step 1
//
// double kineYmax = Kinematics::ymax(mS, vectorMesonPDG->Mass());
double kineYmin = Kinematics::ymin(mS, vectorMesonPDG->Mass());
double kineQ2min = Kinematics::Q2min(kineYmin);
double kineQ2max = Kinematics::Q2max(mS);
double kineW2min = Kinematics::W2min(vectorMesonPDG->Mass());
double kineW2max = Kinematics::W2max(mS);
kineQ2min = max(kineQ2min, mLowerLimit[1]);
kineQ2max = min(kineQ2max, mUpperLimit[1]);
kineW2min = max(kineW2min, mLowerLimit[2]);
kineW2max = min(kineW2max, mUpperLimit[2]);
double kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // first arg (t) set to 0 (recursive)
double kineTmax = Kinematics::tmax(kineXPmin);
double kineTmin = Kinematics::tmin(mHadronBeam.E());
mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax);
mLowerLimit[1] = max(kineQ2min, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineQ2max);
mLowerLimit[2] = max(kineW2min, mLowerLimit[2]); mUpperLimit[2] = min(mUpperLimit[2], kineW2max);
//
// Step 3:
// Deal with user provided limits.
// User settings are ignored (switched off) if min >= max.
//
if (mSettings->W2min() < mSettings->W2max()) { // W2 first
if (mSettings->W2min() < mLowerLimit[2]) {
cout << "Warning, requested lower limit of W (" << mSettings->Wmin() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << sqrt(mLowerLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[2] = mSettings->W2min();
}
if (mSettings->W2max() > mUpperLimit[2]) {
cout << "Warning, requested upper limit of W (" << mSettings->Wmax() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << sqrt(mUpperLimit[2]) << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[2] = mSettings->W2max();
}
}
if (mSettings->Q2min() < mSettings->Q2max()) { // Q2
if (mSettings->Q2min() < mLowerLimit[1]) {
cout << "Warning, requested lower limit of Q2 (" << mSettings->Q2min() << ") "
<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mLowerLimit[1] = mSettings->Q2min();
kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // new Q2min changes tmax
mUpperLimit[0] = min(mUpperLimit[0], Kinematics::tmax(kineXPmin));
}
if (mSettings->Q2max() > mUpperLimit[1]) {
cout << "Warning, requested upper limit of Q2 (" << mSettings->Q2max() << ") "
<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). ";
cout << "Limit has no effect." << endl;
}
else {
mUpperLimit[1] = mSettings->Q2max();
}
}
//
// Check if any phase space is left
//
if (mLowerLimit[0] >= mUpperLimit[0]) {
cout << "Invalid range in t: t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
exit(1);
}
if (mLowerLimit[1] >= mUpperLimit[1]) {
cout << "Invalid range in Q2: Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl;
exit(1);
}
if (mLowerLimit[2] >= mUpperLimit[2]) {
cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl;
exit(1);
}
//
// Print-out limits (all verbose levels)
//
if (mSettings->verbose()) {
cout << "Kinematic limits used for event generation:" << endl;
cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl;
}
//
// Check if the lambda table covers the kinematic range
// or if we have to calculate lambda from the p mean_A table
// directly.
// This is just to inform the user, no action is taken.
// Here we assume that transverse and longitudinal tables
// always have the same range and vheck only one.
//
// If even the p mean_A table is not big enough we switch
// corrections off and inform the user.
//
if (mProtonTableCollection) {
bool lambdaTableExists = mProtonTableCollection->tableExists(transverse, lambda_A);
if (!lambdaTableExists ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, lambda_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, lambda_A)) {
if (!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mLowerLimit[1]+numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mLowerLimit[0]+numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mLowerLimit[2]+numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, mean_A) ||
!mProtonTableCollection->available(mUpperLimit[1]-numeric_limits<float>::epsilon(), mUpperLimit[2]-numeric_limits<float>::epsilon(), mUpperLimit[0]-numeric_limits<float>::epsilon(), transverse, mean_A)) {
cout << "Warning: the kinematic coverage of the table containing the lambda values " << endl;
cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
cout << " then the required range (or does not exist). The proton amplitude " << endl;
cout << " tables that could be used to calculate them on the fly are too" << endl;
cout << " small as well. Corrections will be switched off." << endl;
mSettings->setCorrectForRealAmplitude(false);
mSettings->setCorrectSkewedness(false);
delete mProtonTableCollection;
mProtonTableCollection = 0;
}
else {
if (lambdaTableExists) {
cout << "Warning: the kinematic coverage of the table containing the lambda values " << endl;
cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl;
cout << " then the required range. The missing lambda values will be " << endl;
cout << " calculated using the proton amplitude tables directly." << endl;
}
else {
cout << "Warning: the table containing the lambda values needed for skewedness and/or real " << endl;
cout << " amplitude corrections is not available. The missing lambda values will be" << endl;
cout << " calculated using the referring proton amplitude tables directly." << endl;
}
}
}
}
//
// Setup cross-section functor
// It is this functor that is used by all other functors,
// functions, and wrappers when dealing with cross-sections.
//
if (mCrossSection) delete mCrossSection;
mCrossSection = new CrossSection;
mCrossSection->setTableCollection(mTableCollection);
if (mProtonTableCollection) mCrossSection->setProtonTableCollection(mProtonTableCollection);
//
// UNU.RAN needs the domain (boundaries) and the mode.
// The domain is already defined, here we find the mode, which is tricky.
// The max. cross-section is clearly at the domain boundary in Q2=Q2min.
// The position in W2 and t is not obvious. It sits along a line given
// by tmax(W2). The approach here is to use the BrentMinimizer1D that
// performs first a scan a then a Brent fit.
//
double theMode[3];
theMode[0] = mUpperLimit[0]; // t
theMode[1] = mLowerLimit[1]; // Q2
theMode[2] = mLowerLimit[2]; // W2
ModeFinderFunctor modeFunctor(mCrossSection, theMode[1], vectorMesonPDG->Mass(), mLowerLimit[0], mUpperLimit[0]);
if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
ROOT::Math::BrentMinimizer1D minimizer;
minimizer.SetFunction(modeFunctor, mLowerLimit[2], mUpperLimit[2]);
minimizer.SetNpx(static_cast<int>(mUpperLimit[2]-mLowerLimit[2]));
ok = minimizer.Minimize(10000, 0, 1.e-6);
if (! ok) {
cout << "Error, failed to find mode of pdf." << endl;
exit(1);
}
theMode[2] = minimizer.XMinimum(); // W2
theMode[0] = Kinematics::tmax(0, theMode[1], theMode[2], vectorMesonPDG->Mass()); // first arg (t) must be 0 here
if (theMode[0] > mUpperLimit[0]) theMode[0] = mUpperLimit[0];
double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1], theMode[2]);
if (mSettings->verbose()) {
cout << "\tlocation: t=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]);
cout << "; value: " << crossSectionAtMode << endl;
}
//
// Consistency check of mode
//
if (crossSectionAtMode <= 0) {
cout << "Error: cross-section at mode value is invalid." << endl;
return false;
}
if (! Kinematics::valid(mS, theMode[0], theMode[1], theMode[2], vectorMesonPDG->Mass()) ) {
cout << "Error: mode of pdf is outside kinematic limits." << endl;
return false;
}
//
// Initialize 3D random generator
//
// Test show that UNU.RAN runs smoother in log(Q2)
// and log(cross-section). Functor CrossSection has
// a spezialized method for UNU.RAN, unuranPDF().
//
// domain and mode for Q2 -> log(Q2)
mLowerLimit[1] = log(mLowerLimit[1]);
mUpperLimit[1] = log(mUpperLimit[1]);
theMode[1] = log(theMode[1]);
if (mPDF_Functor) delete mPDF_Functor;
if (mPDF) delete mPDF;
mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 3);
mPDF = new TUnuranMultiContDist(*mPDF_Functor, true); // last arg = pdf in log or not
mPDF->SetDomain(mLowerLimit, mUpperLimit);
mPDF->SetMode(theMode);
if (mUnuran) delete mUnuran;
mUnuran = new TUnuran;
mCrossSection->setCheckKinematics(false); // avoid numeric glitch in Init()
mUnuran->Init(*mPDF, "method=hitro");
mCrossSection->setCheckKinematics(true);
mUnuran->SetSeed(mSettings->seed());
//
// Burn in generator
//
double xrandom[3];
for (int i=0; i<1000; i++) mUnuran->SampleMulti(xrandom);
mEventCounter = 0;
mTriesCounter = 0;
mIsInitialized = true;
cout << "Sartre is initialized." << endl << endl;
return true;
}
bool Sartre::init(const string& str) // overloaded version of init()
{
if (str.empty())
return init();
else
return init(str.c_str());
}
vector<pair<double,double> > Sartre::kinematicLimits()
{
vector<pair<double,double> > array;
array.push_back(make_pair(mLowerLimit[0], mUpperLimit[0])); // t
array.push_back(make_pair(exp(mLowerLimit[1]), exp(mUpperLimit[1]))); // Q2
array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W
return array;
}
Event* Sartre::generateEvent()
{
if (!mIsInitialized) {
cout << "Sartre::generateEvent(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return 0;
}
double xrandom[3];
TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID);
double vmMass = vectorMesonPDG->Mass();
//
// Generate one event
//
while (true) {
mTriesCounter++;
delete mCurrentEvent;
mCurrentEvent = new Event;
//
// Get t, Q2, W2 from TUnuran and check for kinematics.
// Q2 is generated as log(Q2) so we transform it back first.
// This is the only place where Kinematics::valid() is called
// with the fully correct xpomeron calculation switched on.
//
mUnuran->SampleMulti(xrandom);
xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2
if (!Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2], vmMass, true, (mSettings->verboseLevel() > 1))) {
if (mSettings->verboseLevel() > 2) cout << "Sartre::generateEvent(): event rejected, not within kinematic limits" << endl;
continue;
}
//
// Fill beam particles in Event structure
// Kinematics for eA is reported as 'per nucleon'
//
mCurrentEvent->eventNumber = mEventCounter;
mCurrentEvent->t = xrandom[0]; // t
mCurrentEvent->Q2 = xrandom[1]; // Q2
mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
mCurrentEvent->s = mS; // s
mCurrentEvent->W = sqrt(xrandom[2]);
mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = 11; // e-
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
//
// Generate the final state particles
//
bool ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->y, mCurrentEvent->Q2,
(mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent);
if (!ok) {
if (mSettings->verboseLevel() > 1) cout << "Sartre::generateEvent(): failed to generate final state" << endl;
continue;
}
break;
}
mEventCounter++;
//
// Nuclear breakup
//
// If the event is incoherent the final state generator does produce a
// 'virtual' proton with m > m_p which is used in Nucleus to calculate
// the excitation energy and the boost.
//
int indexOfScatteredHadron = 6;
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (mNucleus) mNucleus->resetBreakup(); // clear previous event in any case
if (allowBreakup && mCurrentEvent->diffractiveMode == incoherent && mNucleus) {
int nFragments = mNucleus->breakup(mCurrentEvent->particles[indexOfScatteredHadron].p);
//
// Merge the list of products into the event list.
// We loose some information here. The user can always go back to
// the nucleus and check the decay products for more details.
// In the original list the energy is per nuclei, here we transform it
// to per nucleon to stay consistent with Sartre conventions.
//
const vector<BreakupProduct>& products = mNucleus->breakupProducts();
for (int i=0; i<nFragments; i++) {
Particle fragment;
fragment.index = mCurrentEvent->particles.size();
fragment.pdgId = products[i].pdgId;
fragment.status = 1;
fragment.p = products[i].p*(1/static_cast<double>(products[i].A));
fragment.parents.push_back(indexOfScatteredHadron);
mCurrentEvent->particles.push_back(fragment);
}
}
//
// Complete event record
//
mCurrentEvent->xpom = Kinematics::xpomeron(mCurrentEvent->t, mCurrentEvent->Q2, mCurrentEvent->W*mCurrentEvent->W, vmMass);
mCurrentEvent->beta = mCurrentEvent->x/mCurrentEvent->xpom;
return mCurrentEvent;
}
double Sartre::totalCrossSection()
{
if (mTotalCrossSection == 0) {
//
// Limits of integration in t, Q2, W2
//
double xmin[3];
double xmax[3];
copy(mLowerLimit, mLowerLimit+3, xmin);
copy(mUpperLimit, mUpperLimit+3, xmax);
xmin[1] = exp(xmin[1]); // log Q2 limit -> Q2 limit
xmax[1] = exp(xmax[1]); // log Q2 limit -> Q2 limit
mTotalCrossSection = calculateTotalCrossSection(xmin, xmax);
}
return mTotalCrossSection;
}
double Sartre::totalCrossSection(double lower[3], double upper[3]) // t, Q2, W
{
lower[2] *= lower[2]; upper[2] *= upper[2]; // W -> W2
double result = calculateTotalCrossSection(lower, upper);
return result;
}
EventGeneratorSettings* Sartre::runSettings()
{
return EventGeneratorSettings::instance();
}
double Sartre::calculateTotalCrossSection(double lower[3], double upper[3])
{
double result = 0;
if (!mIsInitialized) {
cout << "Sartre::calculateTotalCrossSection(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return result;
}
//
// Calculate integral using adaptive numerical method
//
// Options: ADAPTIVE, kVEGAS, kPLAIN, kMISER
// no abs tolerance given -> relative only
const double precision = 5e-4;
ROOT::Math::Functor wfL((*mCrossSection), 3);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE,
0, precision, 1000000);
ig.SetFunction(wfL);
result = ig.Integral(lower, upper);
//
// If it fails we switch to a MC integration which is usually more robust
// although not as accurate. This should happen very rarely if at all.
//
if (result <= numeric_limits<float>::epsilon()) {
cout << "Sartre::calculateTotalCrossSection(): warning, adaptive integration failed - switching to VEGAS method." << endl;
ROOT::Math::IntegratorMultiDim igAlt(ROOT::Math::IntegrationMultiDim::kVEGAS);
igAlt.SetFunction(wfL);
igAlt.SetRelTolerance(precision);
igAlt.SetAbsTolerance(0);
result = igAlt.Integral(lower, upper);
}
return result;
}
const FrangibleNucleus* Sartre::nucleus() const {return mNucleus;}
void Sartre::listStatus(ostream& os) const
{
os << "Event summary: " << mEventCounter<< " events generated, " << mTriesCounter << " tried" << endl;
time_t delta = runTime();
os << "Total time used: " << delta/60 << " min " << delta - 60*(delta/60) << " sec" << endl;
}
time_t Sartre::runTime() const
{
time_t now = time(0);
return now-mStartTime;
}
//==============================================================================
//
// Utility functions and operators (helpers)
//
//==============================================================================
ostream& operator<<(ostream& os, const TLorentzVector& v)
{
os << v.Px() << '\t' << v.Py() << '\t' << v.Pz() << '\t' << v.E() << '\t';
double m2 = v*v;
if (m2 < 0)
os << '(' << -sqrt(-m2) << ')';
else
os << '(' << sqrt(m2) << ')';
return os;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:09 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805976
Default Alt Text
(34 KB)

Event Timeline