No OneTemporary
Actions

Size

84 KB

Subscribers

None

View Options

	Index: trunk/src/SartreInclusiveDiffraction.cpp
	===================================================================
	--- trunk/src/SartreInclusiveDiffraction.cpp (revision 550)
	+++ trunk/src/SartreInclusiveDiffraction.cpp (revision 551)
	@@ -1,988 +1,987 @@
	//==============================================================================
	// SartreInclusiveDiffraction.cpp
	//
	// Copyright (C) 2024 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: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
	// $Author: ullrich $
	//==============================================================================
	//
	// Note:
	// When not using runcards, the user must first create an instance of
	// SartreInclusiveDiffraction and then get the settings via one of:
	// SartreInclusiveDiffraction::runSettings()
	// EventGeneratorSettings::instance()
	// Once init() is called settings cannot be changed any more.
	//==============================================================================
	#include "Version.h"
	#include "Kinematics.h"
	#include "Constants.h"
	#include "SartreInclusiveDiffraction.h"
	#include "Math/Functor.h"
	#include "ModeFinderFunctor.h"
	#include "Math/BrentMinimizer1D.h"
	#include "Math/IntegratorMultiDim.h"
	#include "Math/GSLMinimizer.h"
	#include "TUnuranMultiContDist.h"
	#include <string>
	#include <limits>
	#include <cmath>
	#include <iomanip>
	#include "TH2D.h"
	#include "TFile.h"
	#include "TF3.h"
	#include "InclusiveDiffractiveCrossSectionsFromTables.h"
	#include "InclusiveDiffractiveCrossSections.h"
	#include "DglapEvolution.h"

	using namespace std;

	#define PR(x) cout << #x << " = " << (x) << endl;

	SartreInclusiveDiffraction::SartreInclusiveDiffraction()
	{
	mSettings = EventGeneratorSettings::instance();
	mRandom = mSettings->randomGenerator();
	mIsInitialized = false;
	mCurrentEvent = 0;
	mNucleus = 0;
	mUpcNucleus= 0;
	// mSettings = 0;
	mPDF_Functor = 0;
	mPDF = 0;
	mEventCounter = 0;
	mTriesCounter = 0;
	mTotalCrossSection = 0;
	mTotalCrossSection_qqg = 0;
	mTotalCrossSection_qq = 0;
	mCrossSection = 0;
	mUnuran = 0;
	mEvents = 0;
	mTries = 0;
	mS = 0;
	mA = 0;
	}

	SartreInclusiveDiffraction::~SartreInclusiveDiffraction()
	{
	delete mNucleus;
	delete mUpcNucleus;
	delete mPDF_Functor;
	delete mPDF;
	delete mCrossSection;
	delete mUnuran;
	delete mCurrentEvent;

	}

	bool SartreInclusiveDiffraction::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;
	mTotalCrossSection_qqg = 0;
	mTotalCrossSection_qq = 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 inclusive diffraction \|" << endl;
	cout << "\| in ep and eA collisions based on the dipole model. \|" << endl;
	cout << "\| \|" << endl;
	cout << "\| Copyright (C) 2010-2024 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
	mSettings->setInclusiveDiffractionMode(true);

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

	mUseCrossSectionTables=mSettings->useInclusiveCrossSectionTables();

	TableGeneratorSettings* tgsettings=TableGeneratorSettings::instance();
	tgsettings->setDipoleModelType(mSettings->dipoleModelType());
	tgsettings->setDipoleModelParameterSet(mSettings->dipoleModelParameterSet());
	// Set up DGLAP evolution for integrals:
	if(!mUseCrossSectionTables){
	DglapEvolution &dglap = DglapEvolution::instance();
	dglap.generateLookupTable(1000, 1000); dglap.useLookupTable(true);
	}
	//
	// 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);

	string upcNucleusName;
	cout << "Sartre is running in inclusive diffraction mode" << endl;
	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();
	mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
	cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
	cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
	cout << "Process is ";
	if (mA > 1)
	cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + X"<< endl;
	else
	cout << "e + p -> e' + p' + X"<< 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);
	}

	//
	// 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.
	// For UPC that's xpom.
	// 2. the range of the table(s)
	// 3. the kinematically allowed range
	//
	// Of course (3) is more complex than a simple cube/square.
	// 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/square
	// envelope that satifies (1)-(3).
	// Note, that they are correlated which makes the order
	// in which we do things a bit tricky.
	//


	//
	// Step 2:
	// Kinematic limits might overrule boundaries from step 1
	//
	// let:
	// 0: beta
	// 1: Q2
	// 2: W2
	// 3: z
	//
	// 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;
	if(!mUseCrossSectionTables){
	mCrossSection = new InclusiveDiffractiveCrossSectionsIntegrals;
	}
	else{
	mCrossSection = new InclusiveDiffractionCrossSectionsFromTables;
	}
	//
	// Here we need to put in Q2 > mu02.
	// The calculation goes haywire for Q2 becoming too small
	//
	double mu02 = mCrossSection->dipoleModel()->getParameters()->mu02();
	// double mf = Settings::quarkMass(0);
	double mf=tt_quarkMass[0];

	//Now for the diffractive mass:
	double kineMX2min = 4mfmf;//max(mu02, 4mfmf);
	// This gives the limits for W2 and for inelasticity y:
	double kineW2min = kineMX2min+protonMass2;
	double kineYmin = Kinematics::ymin(mS, kineMX2min);
	// ...and y gives the limit for Q2, unless it's smaller than the cut-off mu02:
	double kineQ2min = Kinematics::Q2min(kineYmin);
	kineQ2min = max(kineQ2min, mu02);
	//The maximum momenta are given by s:
	double kineQ2max = Kinematics::Q2max(mS);
	double kineW2max = mS-kineQ2min-protonMass2;

	//
	// Now for all the fractions, making sure they make sense
	//
	double kineBetamax = kineQ2max/(kineQ2max+kineMX2min);
	double kineBetamin = Kinematics::betamin(kineQ2min, mS, sqrt(kineW2min), 0);
	double kineMX2max = (1.-kineBetamin)/kineBetamin*kineQ2max;
	kineMX2max=min(kineMX2max, kineW2max-protonMass2);
	double kineZmin=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	double kineZmax=1-kineZmin;
	mLowerLimit[0] = kineBetamin;
	mUpperLimit[0] = kineBetamax;
	mLowerLimit[1] = kineQ2min;
	mUpperLimit[1] = kineQ2max;
	mLowerLimit[2] = kineW2min;
	mUpperLimit[2] = kineW2max;
	mLowerLimit[3] = kineZmin;
	mUpperLimit[3] = kineZmax;

	// Step 3:
	// Deal with user provided limits.
	// User settings are ignored (switched off) if min >= max.
	// Only allow user limits for Q2, W2, and beta,
	// not z.
	//
	bool W2orQ2changed=false;
	if (mSettings->Wmin() < mSettings->Wmax()) { // 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();
	W2orQ2changed=true;
	}

	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();
	W2orQ2changed=true;
	}
	}
	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();
	W2orQ2changed=true;
	}

	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();
	W2orQ2changed=true;
	}
	}
	if (W2orQ2changed){//mLowerLimit beta, Q2, W2, z
	//kineBetamax = kineQ2max/(kineQ2max+kineMX2min);
	mUpperLimit[0] = mUpperLimit[1]/(mUpperLimit[1]+kineMX2min);

	//kineBetamin = Kinematics::betamin(kineQ2min, mS, sqrt(kineW2min), 0);
	mLowerLimit[0] = Kinematics::betamin(mLowerLimit[1], mS, sqrt(mLowerLimit[2]), 0);

	// kineMX2max = (1.-kineBetamin)/kineBetamin*kineQ2max;
	// kineMX2max=min(kineMX2max, kineW2max-protonMass2);
	kineMX2max = (1.-mLowerLimit[0])/mLowerLimit[0]*mUpperLimit[1];
	kineMX2max=min(kineMX2max, mUpperLimit[2]-protonMass2);

	// kineZmin=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	// kineZmax=1-kineZmin;
	mLowerLimit[3]=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	mUpperLimit[3]=1-mLowerLimit[3];
	}

	if (mSettings->betamin() < mSettings->betamax()) { // beta

	if (mSettings->betamin() < mLowerLimit[0]) {
	cout << "Warning, requested lower limit of beta (" << mSettings->betamin() << ") "
	<< "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[0] << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mLowerLimit[0] = mSettings->betamin();
	}

	if (mSettings->betamax() > mUpperLimit[0]) {
	cout << "Warning, requested upper limit of Q2 (" << mSettings->betamax() << ") "
	<< "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[0] << "). ";
	cout << "Limit has no effect." << endl;
	}
	else {
	mUpperLimit[0] = mSettings->betamax();
	}
	}

	//
	// Check if any phase space is left
	//
	if (mLowerLimit[0] >= mUpperLimit[0]) {
	cout << "Invalid range in beta: beta=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
	exit(1);
	}
	if (mLowerLimit[1] >= mUpperLimit[1]) {
	cout << "Invalid range in Q2: Q2=[";
	cout << 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);
	}
	if (mLowerLimit[3] >= mUpperLimit[3]) {
	cout << "Invalid range in z: z=[" << sqrt(mLowerLimit[3]) << ", " << sqrt(mUpperLimit[3]) << "]." << endl;
	exit(1);
	}

	//
	// Print-out limits (all verbose levels)
	//
	if (true){//}(mSettings->verbose()) {
	cout << "Sartre was thrown into the world, restricted by kinematics and user inputs alike:" << endl;
	cout << setw(10) << " beta=[" << 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;
	cout << setw(10) << " z=[" << mLowerLimit[3] << ", " << mUpperLimit[3] << "]" << endl;
	cout << setw(10) << " MX2=[" << kineMX2min << ", " << kineMX2max << "]" << endl;
	cout << setw(10) << " y=[" << kineYmin << ", 1" << "]" << endl;

	}
	//
	// 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 beta is not obvious. The mode should be at z=0.5.
	// The approach here is to use the BrentMinimizer1D that
	// performs first a scan a then a Brent fit.
	//
	mCrossSection->setCheckKinematics(false);
	double theMode[4];

	//
	// Find the mode.
	// Assume that the mode is at W2=W2max, Q2=Q2min, z=0.5
	// Start with beta=0.5 => MX2=Q2
	// It is kinemtically easier to use MX2 instead of beta here.
	//
	if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
	theMode[1] = mLowerLimit[1]; // Q2
	theMode[2] = mUpperLimit[2]; // W2
	theMode[3] = 0.5; // z

	double MX2min=mLowerLimit[1]*(1-mUpperLimit[0])/mUpperLimit[0];
	MX2min=max(kineMX2min, MX2min);
	double MX2max=mLowerLimit[1]*(1-mLowerLimit[0])/mLowerLimit[0];
	MX2max=min(kineMX2max, MX2max);
	InclusiveDiffractionModeFinderFunctor modeFunctor(mCrossSection, theMode[1], theMode[2], theMode[3], MX2min, MX2max);

	ROOT::Math::BrentMinimizer1D minimizer;
	minimizer.SetFunction(modeFunctor, MX2min, MX2max);
	minimizer.SetNpx(static_cast<int>(mUpperLimit[0]-mLowerLimit[0]));
	ok = minimizer.Minimize(100000, 0, 1.e-8);
	if (! ok) {
	cout << "Error, failed to find mode of pdf." << endl;
	exit(1);
	}
	theMode[0] = minimizer.XMinimum(); // Mx2

	double MX2=theMode[0];
	theMode[0]=theMode[1]/(theMode[1]+MX2); //Mx2->beta
	double crossSectionAtMode = (*mCrossSection)(MX2, theMode[1], theMode[2], theMode[3]);
	if (mSettings->verbose()) {
	cout << "\tlocation: beta=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]) << ", z=" << theMode[3];
	cout << "; value: " << crossSectionAtMode << endl;
	}
	mCrossSection->setCheckKinematics(true);

	// 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, &InclusiveDiffractiveCrossSections::unuranPDF, 4);
	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[4];
	for (int i=0; i<100; i++) {
	mUnuran->SampleMulti(xrandom);
	// cout<<"QQ: beta="<<xrandom[0];
	// cout<<" Q2="<<exp(xrandom[1]);
	// cout<<" W="<<sqrt(xrandom[2]);
	// cout<<" z="<<xrandom[3]<<endl;
	}

	//
	// Set up QQG unu.ran:
	//
	//The mode is simpler for QQG, as it is at minimum beta
	//The limits are the same, except for z which is a different quantity in QQG.
	mLowerLimit_qqg[0]=mLowerLimit[0]; //beta
	mLowerLimit_qqg[1]=mLowerLimit[1]; //log(Q2)
	mLowerLimit_qqg[2]=mLowerLimit[2]; //W2
	//z>beta(1+4mf*mf/Q2);
	mLowerLimit_qqg[3]=mLowerLimit[0] * (1+4mfmf/exp(mUpperLimit[1]))+1e-10;

	mUpperLimit_qqg[0]=mUpperLimit[0]; //beta
	mUpperLimit_qqg[1]=mUpperLimit[1]; //log(Q2)
	mUpperLimit_qqg[2]=mUpperLimit[2]; //W2
	mUpperLimit_qqg[3]=mUpperLimit[3]; //z

	double theMode_qqg[4];
	theMode_qqg[0]=mLowerLimit_qqg[0]; //beta
	theMode_qqg[1]=mLowerLimit_qqg[1]; //log(Q2)
	theMode_qqg[2]=mUpperLimit_qqg[2]; //W2
	theMode_qqg[3]=theMode_qqg[0](1+4mf*mf/exp(theMode_qqg[1]))+1e-10; //z

	if (mPDF_Functor_qqg) delete mPDF_Functor_qqg;
	if (mPDF_qqg) delete mPDF_qqg;

	mPDF_Functor_qqg = new ROOT::Math::Functor(mCrossSection, &InclusiveDiffractiveCrossSections::unuranPDF_qqg, 4);
	mPDF_qqg = new TUnuranMultiContDist(*mPDF_Functor_qqg, true); // last arg = pdf in log or not
	//Rescale the domain in z (and fix it in the pdf)
	mPDF_qqg->SetDomain(mLowerLimit_qqg, mUpperLimit_qqg);
	mPDF_qqg->SetMode(theMode_qqg);

	mCrossSection->setCheckKinematics(false);
	if (mUnuran_qqg) delete mUnuran_qqg;
	mUnuran_qqg = new TUnuran;
	mUnuran_qqg->Init(*mPDF_qqg, "method=hitro");
	mUnuran_qqg->SetSeed(mSettings->seed());
	// mCrossSection->setCheckKinematics(true);
	//
	// Burn in QQG generator
	//
	for (int i=0; i<100; i++) {
	mUnuran_qqg->SampleMulti(xrandom);
	// cout<<"QQG: beta="<<xrandom[0];
	// cout<<" Q2="<<exp(xrandom[1]);
	// cout<<" W="<<sqrt(xrandom[2]);
	// cout<<" z="<<xrandom[3]<<endl;
	}
	mEventCounter = 0;
	mTriesCounter = 0;
	mIsInitialized = true;
	cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;


	//
	// Calculate total cross section in given
	// kinematic limits. This will be used to
	// decide which type of event to generate,
	// q-qbar or q-qbar-g
	//
	cout<<"Calculating cross sections over the given kinematic range..."<<endl;
	cout<<"(This may take while, use your radical freedom and do other things.)"<<endl;
	mCrossSection->setFockState(QQ);
	mTotalCrossSection=0;
	mTotalCrossSection_qq=7.4543e+08; //totalCrossSection();
	cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
	mCrossSection->setFockState(QQG);
	mTotalCrossSection=0;
	mTotalCrossSection_qqg=3.79917e+08; //totalCrossSection();
	cout<<"The QQG cross section is: "<<mTotalCrossSection_qqg<<" nb. (2/2)"<<endl;
	cout<<"Total Cross Section: "<<mTotalCrossSection_qqg+mTotalCrossSection_qq<<" nb"<<endl;
	mTotalCrossSection=mTotalCrossSection_qqg+mTotalCrossSection_qq;
	return true;
	}


	double SartreInclusiveDiffraction::cross_section_wrapper(double* xx, double* par){
	double arg[4]={xx[0], xx[1], xx[2], par[0]};
	double result=mCrossSection->unuranPDF(arg);
	return -result;
	}

	bool SartreInclusiveDiffraction::init(const string& str) // overloaded version of init()
	{
	if (str.empty())
	return init();
	else
	return init(str.c_str());
	}

	vector<pair<double,double> > SartreInclusiveDiffraction::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 or xpom
	if (!mSettings->UPC())
	array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W
	return array;
	}

	Event* SartreInclusiveDiffraction::generateEvent()
	{
	if (!mIsInitialized) {
	cout << "SartreInclusiveDiffraction::generateEvent(): Error, Sartre is not initialized yet." << endl;
	cout << " Call init() before trying to generate events." << endl;
	return 0;
	}
	//
	// Decide which Fock-state to generate:
	//
	bool isQQ=false;
	if(mRandom->Uniform(mTotalCrossSection) <= mTotalCrossSection_qq)
	isQQ=true;
	//
	// Generate one event
	//
	double xrandom[4];
	while (true) {
	mTriesCounter++;
	delete mCurrentEvent;
	mCurrentEvent = new Event;
	//
	// xrandom[i]
	// i=0: beta
	// i=1: Q2
	// i=2: W2
	// i=3: z
	//
	if(isQQ)
	mUnuran->SampleMulti(xrandom);
	else
	mUnuran_qqg->SampleMulti(xrandom);
	xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2
	double MX2=xrandom[1]*(1-xrandom[0])/xrandom[0];
	bool isValidEvent;
	isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2], xrandom[3], sqrt(MX2), true, (mSettings->verboseLevel() > 1));
	if (!isValidEvent) {
	if (mSettings->verboseLevel() > 2)
	cout << "SartreInclusiveDiffraction::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->beta = xrandom[0]; // beta
	mCurrentEvent->Q2 = xrandom[1]; // Q2
	mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
	mCurrentEvent->xpom=mCurrentEvent->x/xrandom[0]; //xpom=x/beta
	mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
	mCurrentEvent->s = mS; // s
	mCurrentEvent->W = sqrt(xrandom[2]);
	mCurrentEvent->z = xrandom[3];
	mCurrentEvent->MX = sqrt(MX2);
	if(isQQ)
	mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
	else
	mCurrentEvent->polarization = transverse;
	// mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
	mCurrentEvent->quarkSpecies = mCrossSection->quarkSpeciesOfLastCall();
	if(isQQ)
	mCurrentEvent->crossSectionRatioLT=mCrossSection->crossSectionRatioLTOfLastCall();
	else
	mCurrentEvent->crossSectionRatioLT=0;

	//
	// For generating final state azimuthal decorelation
	// we need to calculate Omega(b, r).
	// For this we need to generate b and r values.
	//
	// Start with b (written by Bhakta Naik and TT):
	// From the thickness function (a leading twist approximation)
	//
	double bmin = 0;
	double bmax = 3.5 * mCrossSection->dipoleModel()->nucleus()->radius();

	TF1* fbDistribution = new TF1("fbDistribution",this ,&SartreInclusiveDiffraction::bDistribution, bmin, bmax, 0);

	double bGenerated = fbDistribution->GetRandom()*hbarc; // fm
	//
	// Now generate r
	//
	double mf = mCrossSection->quarkMassOfLastCall();
	double rmin=1e-3;
	double rmax=3.5 * mCrossSection->dipoleModel()->nucleus()->radius();
	double rGenerated=0;
	if(isQQ)
	{
	TF1* frDistribution_qq;
	if(mCurrentEvent->polarization == transverse )
	{
	frDistribution_qq = new TF1("frDistribution_qq", this, &SartreInclusiveDiffraction::rDistribution_qqT, rmin, rmax, 6);
	}
	else
	{
	frDistribution_qq = new TF1("frDistribution_qq", this, &SartreInclusiveDiffraction::rDistribution_qqL, rmin, rmax, 6);
	}
	frDistribution_qq->SetParameters(mCurrentEvent->beta, mCurrentEvent->xpom, mCurrentEvent->z, mCurrentEvent->Q2, bGenerated, mf);
	rGenerated = frDistribution_qq->GetRandom();

	delete frDistribution_qq;
	}
	else{ //QQG
	// Create a TF2 object.
	// The range for 'r' is from 0 to max_r,
	// and for 'kappa' is from 0 to Q2.
	// The parameter array for TF2 is defined as:
	// {Q2, xpom, z, b_rand_value}
	//
	TF2* frDistribution_qqg = new TF2("frDistribution_qqg", this,&SartreInclusiveDiffraction::rDistribution_qqg, rmin, rmax, 0, mCurrentEvent->Q2, 4);

	frDistribution_qqg->SetParameters(mCurrentEvent->xpom, mCurrentEvent->z, bGenerated, sqrt(mCurrentEvent->Q2));

	double kGenerated;
	frDistribution_qqg->GetRandom2(rGenerated, kGenerated);
	delete frDistribution_qqg;
	}
	//Calculate Omega:
	double dsigmadb2=0;
	if (mA == 1) {
	dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(rGenerated, bGenerated, mCurrentEvent->xpom);
	}
	else {
	dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(rGenerated, bGenerated, mCurrentEvent->xpom);
	}
	mCurrentEvent->Omega=-2.*log(1-dsigmadb2/2.);
	// Calculate the saturation scale following
	// H. Kowalski, T. Lappi and R. Venugopalan, Phys. Rev. Lett. 100 (2008) 022303 [arXiv:0705.3047 [hep-ph]]
	double rS;
	for(rS=0.1; rS<20; rS+=0.001){
	double dsigmadb2=0;
	if (mA == 1) {
	dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(rS, bGenerated, mCurrentEvent->xpom);
	}
	else {
	dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(rS, bGenerated, mCurrentEvent->xpom);
	}
	if(dsigmadb2 >= 0.44239843) break;
	}
	mCurrentEvent->QS=1./rS*hbarc; //GeV
	- double QS=1./rS*hbarc;

	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;
	if (isQQ) {
	ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQ);
	}
	else {
	ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQG);
	}
	if (!ok) {
	if (mSettings->verboseLevel() > 1) cout << "SartreInclusiveDiffraction::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);
	}
	}
	return mCurrentEvent;
	}

	double SartreInclusiveDiffraction::totalCrossSection()
	{
	if (mTotalCrossSection == 0) {
	//
	// Limits of integration in beta, Q2, W2, z
	//
	double xmin[4];
	double xmax[4];
	copy(mLowerLimit, mLowerLimit+4, xmin);
	copy(mUpperLimit, mUpperLimit+4, xmax);

	//
	// At this point mLowerLimit[1] and mUpperLimit[1]
	// are in log(Q2)
	//
	xmin[1] = exp(xmin[1]); // log Q2 limit
	xmax[1] = exp(xmax[1]); // log Q2 limit

	mTotalCrossSection = calculateTotalCrossSection(xmin, xmax);
	}
	return mTotalCrossSection;
	}

	double SartreInclusiveDiffraction::totalCrossSection(double lower[4], double upper[4]) // beta, Q2, W, z
	{
	lower[2] = lower[2]; upper[2] = upper[2]; // W -> W2
	double result = calculateTotalCrossSection(lower, upper);
	return result;
	}

	double SartreInclusiveDiffraction::integrationVolume(){
	//
	// Limits of integration in beta, Q2, W2, z
	//
	double xmin[4];
	double xmax[4];
	copy(mLowerLimit, mLowerLimit+4, xmin);
	copy(mUpperLimit, mUpperLimit+4, xmax);
	//
	// At this point mLowerLimit[1] and mUpperLimit[1]
	// are in log(Q2)
	//
	xmin[1] = exp(xmin[1]); // log Q2 limit
	xmax[1] = exp(xmax[1]); // log Q2 limit

	double volume = (xmax[0]-xmin[0]) * (xmax[1]-xmin[1]) * (xmax[2]-xmin[2]) * (xmax[3]-xmin[3]);
	return volume; //GeV4
	}

	EventGeneratorSettings* SartreInclusiveDiffraction::runSettings()
	{
	return EventGeneratorSettings::instance();
	}

	double SartreInclusiveDiffraction::crossSectionsClassWrapper(const double* var){
	return (*mCrossSection)(var);
	}

	double SartreInclusiveDiffraction::calculateTotalCrossSection(double lower[4], double upper[4])
	{
	double result = 0;

	if (!mIsInitialized) {
	cout << "SartreInclusiveDiffraction::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;
	// unsigned int numberofcalls = 1000000;
	unsigned int numberofcalls = 0;//1e5;
	ROOT::Math::Functor wfL(this, &SartreInclusiveDiffraction::crossSectionsClassWrapper, 4);
	ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE, 0, precision, numberofcalls);

	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 << "SartreInclusiveDiffraction::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;
	}

	double SartreInclusiveDiffraction::bDistribution(double* var, double*)
	{
	double b = *var;

	double result=0;
	if(mA==1){
	double BG = mCrossSection->dipoleModel()->getParameters()->BG();
	double arg = bb/(2BG);
	arg /= hbarc2;
	result = 1/(2M_PIBG) * exp(-arg);
	}
	else{
	result = mNucleus->T(b);
	}
	return result;
	}
	// Three functions for randomly choosing a value of r:
	double SartreInclusiveDiffraction::rDistribution_qqT(double* var, double* par)
	{
	double result = mCrossSection->uiAmplitude_1(var, par);
	if(isnan(result)) return 0;
	return result * result;
	}

	double SartreInclusiveDiffraction::rDistribution_qqL(double* var, double* par)
	{
	double result = mCrossSection->uiAmplitude_0(var, par);

	if(isnan(result)) return 0;
	return result * result;
	}
	double SartreInclusiveDiffraction::rDistribution_qqg(double* var, double* par)
	{
	double r = var[0]; //fm
	double k = var[1]; //GeV
	if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
	double xpom=par[0];
	double z=par[1];

	double b=par[2];
	double Q2 = par[3];

	double para[4] = {xpom, z, k , b};
	double result = mCrossSection->uiAmplitude_qqg(&r, para);
	return pow(k, 4) * log(Q2 / (k * k)) * result * result ;
	}


	const FrangibleNucleus* SartreInclusiveDiffraction::nucleus() const {return mNucleus;}

	void SartreInclusiveDiffraction::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 SartreInclusiveDiffraction::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;
	}
	Index: trunk/src/InclusiveFinalStateGenerator.cpp
	===================================================================
	--- trunk/src/InclusiveFinalStateGenerator.cpp (revision 550)
	+++ trunk/src/InclusiveFinalStateGenerator.cpp (revision 551)
	@@ -1,1117 +1,815 @@
	//==============================================================================
	// InclusiveFinalStateGenerator.cpp
	//
	// Copyright (C) 2024 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: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
	// $Author: ullrich $
	//==============================================================================
	#include "InclusiveFinalStateGenerator.h"
	#include "EventGeneratorSettings.h"
	#include "Kinematics.h"
	#include "Event.h"
	#include "Math/BrentRootFinder.h"
	#include "Math/GSLRootFinder.h"
	#include "Math/RootFinderAlgorithms.h"
	#include "Math/IFunction.h"
	#include "Math/SpecFuncMathCore.h"
	#include <cmath>
	#include <algorithm>
	#include <limits>
	#include <iomanip>
	#include "Math/WrappedTF1.h"
	#include "Math/GaussIntegrator.h"
	#include "TF1.h"
	#include "Enumerations.h"
	#include "Constants.h"

	#define PR(x) cout << #x << " = " << (x) << endl;

	//-------------------------------------------------------------------------------
	//
	// Implementation of ExclusiveFinalStateGenerator
	//
	//-------------------------------------------------------------------------------

	InclusiveFinalStateGenerator::InclusiveFinalStateGenerator()
	{
	//
	// Setup Pythia8 which is need to hadronize the dipole (and higher Fock states)
	// We use a dirty hack here to suppress PYTHIA output by setting the cout failbit
	//
	cout.setstate(std::ios::failbit);

	mPythia = new Pythia8::Pythia("", false); // 2nd arg implies no banner
	mPythiaEvent = &(mPythia->event);

	mPdt = &(mPythia->particleData);
	mPythia->readString("ProcessLevel:all = off"); // needed
	mPythia->readString("Check:event = off");

	bool ok = mPythia->init();
	cout.clear();

	if (!ok) {
	cout << "Error initializing Pythia8. Stop here." << endl;
	exit(1);
	}
	else {
	cout << "Pythia8 initialized (" << PYTHIA_VERSION_INTEGER << ")" << endl;
	}
	}

	InclusiveFinalStateGenerator::~InclusiveFinalStateGenerator()
	{
	// mPythia->stat();
	delete mPythia;
	}

	-double InclusiveFinalStateGenerator::uiGamma_N(double* var, double* par) const
	-{
	- double t=var[0];
	- double s=par[0];
	- double result = pow(t, s-1)*exp(-t);
	- return result;
	-}
	-
	-double InclusiveFinalStateGenerator::gamma_N(unsigned int i, double val)
	-{
	- TF1 fGamma_N("fGamma_N", this, &InclusiveFinalStateGenerator::uiGamma_N, -10, 10, 1);
	- fGamma_N.SetParameter(0, i);
	- ROOT::Math::WrappedTF1 wfGamma_N(fGamma_N);
	- ROOT::Math::GaussIntegrator giGamma_N;
	- giGamma_N.SetFunction(wfGamma_N);
	- giGamma_N.SetAbsTolerance(0.);
	- giGamma_N.SetRelTolerance(1e-5);
	-
	- int sign=1;
	- double low=0, up=val;
	- if(val<0){
	- low=val;
	- up=0;
	- sign=-1;
	- }
	- double Tfact=1;
	- for(int j=i-1; j>0; j--){
	- Tfact*=j;
	- }
	- return sign*giGamma_N.Integral(low, up)/Tfact;
	-}
	-
	-bool InclusiveFinalStateGenerator::decayPseudoParticle(TLorentzVector In, TLorentzVector& Out1, TLorentzVector& Out2, double m_1, double m_2, double pt=0){
	+bool InclusiveFinalStateGenerator::decayPseudoParticle(TLorentzVector In, TLorentzVector& Out1, TLorentzVector& Out2, double m_1, double m_2, double z, bool isPplus=true){
	//
	- // Calculate the momenta in the In particle's restframe
	+ // Calculates the outgoing 4-momenta Out1 and Out2
	+ // with masses m_1 and m_2, when particle In
	+ // with four-momentum P and mass M decays.
	+ // Particle Out1 has P+ momentum fraction z of In+.
	+ // (P+ from proton side, P- from electron)
	+ // For P- fractions exchange z->(1-z)
	+ //
	+ // Formula:
	+ // pt2=(\vec p_{\perp 1}-z\vec P_\perp)^2 =
	+ // z(1-z)M^2-(1-z)m_1^2-zm_2^2
	+ // px=zPx+pt*cos(phi)
	+ // py=zPy+pt*sin(phi)
	+ // phi is a random angle
	//
	EventGeneratorSettings *settings = EventGeneratorSettings::instance();
	TRandom3 *rndm = settings->randomGenerator();
	-
	- TVector3 boost=TVector3(In.Px()/In.E(), In.Py()/In.E(), In.Pz()/In.E());
	+// if(!isPplus) z=1-z;
	double M2=In.M2();
	+ double Px=In.Px();
	+ double Py=In.Py();
	+
	double m2_1=m_1*m_1;
	double m2_2=m_2*m_2;
	- //In the rest-frame:
	- double p2=2m2_1m2_2 + 2(m2_1+m2_2)M2 - m2_1m2_1 - m2_2m2_2 - M2*M2;
	- p2/=-4*M2;
	- if(p2<0) return false;
	- double p=sqrt(p2);
	- double E_1, E_2, px, py, pz;
	- if(pt==0){ //no user preference
	- //Generate a direction for decay
	- double phi=rndm->Uniform(2*M_PI);
	- double cosTheta=rndm->Uniform(2)-1;
	- double theta=acos(cosTheta);
	-
	- px=pcos(phi)sin(theta);
	- py=psin(phi)sin(theta);
	- pz=p*cosTheta;
	-
	- E_1=sqrt(m2_1+p2);
	- E_2=sqrt(m2_2+p2);
	+
	+ //We want E, px, py, pz:
	+ double E, px, py, pz;
	+ double pt2;
	+ if(isPplus)
	+ pt2=z(1-z)M2-(1-z)m2_1-zm2_2;
	+ else
	+ pt2=z(1-z)M2-zm2_1-(1-z)m2_2;
	+ if(pt2<0) return false;
	+ double phi=rndm->Uniform(2*M_PI);
	+ px=zPx+sqrt(pt2)cos(phi);
	+ py=zPy+sqrt(pt2)sin(phi);
	+ double p1_plus, p1_minus;
	+ if(isPplus){
	+ p1_plus=z*(In.E()+In.Pz());
	+ p1_minus=(pxpx+pypy+m2_1)/p1_plus;
	}
	- else{ //user defined pt
	- double phi=rndm->Uniform(2*M_PI);
	- px=pt*cos(phi);
	- py=pt*sin(phi);
	- pz=sqrt(p2-pxpx-pypy);
	-
	- E_1=sqrt(m2_1+p2);
	- E_2=sqrt(m2_2+p2);
	+ else{
	+ p1_minus=z*(In.E()-In.Pz());
	+ p1_plus=(pxpx+pypy+m2_1)/p1_minus;
	}
	- Out1=TLorentzVector(px, py, pz, E_1);
	- Out2=TLorentzVector(-px, -py, -pz, E_2);
	-
	- Out1.Boost(boost);
	- Out2.Boost(boost);
	+ E=(p1_plus+p1_minus)/2;
	+ pz=(p1_plus-p1_minus)/2;
	+ Out1=TLorentzVector(px, py, pz, E);
	+ Out2=In-Out1;

	return true;
	}

	bool InclusiveFinalStateGenerator::generate(int A, Event *event, FockState fockstate)
	{
	//
	// Get generator settings and the random generator
	//
	EventGeneratorSettings *settings = EventGeneratorSettings::instance();
	TRandom3 *rndm = settings->randomGenerator();

	//
	// The beam particles must be present in the event list
	//
	int ePos = -1;
	int hPos = -1;
	bool parentsOK = true;
	if (event->particles.size() == 2) {
	if (abs(event->particles[0].pdgId) == 11) {
	ePos = 0;
	hPos = 1;
	}
	else if (abs(event->particles[1].pdgId) == 11) {
	ePos = 1;
	hPos = 0;
	}
	else
	parentsOK = false;
	}
	else
	parentsOK = false;

	if (!parentsOK) {
	cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
	return false;
	}
	//
	// Store arguments locally
	// (Some could also be obtained from the event structure)
	//
	mA = A;
	double xpom=event->xpom;
	double beta=event->beta;
	mT = Kinematics::tmax(xpom);//event->t;
	if (mT > 0) mT = -mT; // ensure t<0
	mQ2 = event->Q2;
	mY = event->y;
	mElectronBeam = event->particles[ePos].p;
	mHadronBeam = event->particles[hPos].p;
	mMX = event->MX;
	mS = Kinematics::s(mElectronBeam, mHadronBeam);
	double MX=event->MX;
	double MX2=MX*MX;

	//
	// Constants
	//
	double const twopi = 2*M_PI;
	double const hMass2 = mHadronBeam.M2();
	mMY2 = hMass2;

	//
	// Re-engineer scattered electron
	//
	// e'=(E', pt', pz') -> 3 unknowns
	//
	// Three equations:
	// 1: meme=E'E'-pt'pt'-pz'pz'
	// 2: Q2=-(e-e')^2=-2meme + 2(EE'-pz*pz')
	// 3: W2=(P+e-e')^2=mp2+2me2+2(EpE-Pzpz)-2(EpE'-Pzpz')-2(EE'-pzpz')
	//

	double Ee=mElectronBeam.E();
	double Pe=mElectronBeam.Pz();
	double Ep=mHadronBeam.E();
	double Pp=mHadronBeam.Pz();
	double W=event->W;
	double W2=W*W;
	double Q2=event->Q2;
	// Take masses from the beams in case they are not actually electrons or protons
	double me2=mElectronBeam.M2();
	double mp2=mHadronBeam.M2();
	//
	// What we want for each particle:
	//
	double E, pz, pt, px, py, phi;

	//
	// Equations 2 and 3 yield:
	//
	E = Pe(W2-mp2-2EeEp) + (Pp+Pe)Q2 + 2PePePp + 2me2*Pp;
	E /= 2(EePp-Ep*Pe);
	pz = Ee(W2-mp2) + (Ep+Ee)Q2 + 2EePePp + 2Epme2 - 2EeEeEp;
	pz /= 2(EePp-Ep*Pe);
	//
	// Equation 1:
	//
	pt = sqrt(EE-pzpz-me2);
	phi = rndm->Uniform(twopi);
	TLorentzVector theScatteredElectron(ptsin(phi), ptcos(phi), pz, E);
	//
	// Re-engineer virtual photon
	//
	// gamma=E-E'
	E=mElectronBeam.E()-theScatteredElectron.E();
	pz=mElectronBeam.Pz()-theScatteredElectron.Pz();
	px=mElectronBeam.Px()-theScatteredElectron.Px();
	py=mElectronBeam.Py()-theScatteredElectron.Py();
	TLorentzVector theVirtualPhoton = TLorentzVector(px, py, pz, E);

	//
	// Re-engineer scattered proton/dissociated proton
	//
	double Pplus=mHadronBeam.E()+mHadronBeam.Pz();
	double PplusPrime=(1-xpom)*Pplus;
	double PminusPrime=(mMY2-mT)/PplusPrime;
	E=(PplusPrime+PminusPrime)/2;
	pz=(PplusPrime-PminusPrime)/2;
	double pt_squared = EE-pzpz-mMY2;
	if (pt_squared < 0) {
	if (settings->verbose() && settings->verboseLevel() > 2) {
	cout<< "InclusiveFinalStateGenerator::generate(): No phase-space for scattered proton pt." << endl;
	cout<< " Abort event generation." << endl;
	}
	return false;
	}
	pt = sqrt(pt_squared);
	px = pt*cos(phi);
	py = pt*sin(phi);
	TLorentzVector theScatteredProton(px, py, pz, E);
	//
	// Use scattered proton to set up four momentum of pomeron:
	//
	TLorentzVector thePomeron=mHadronBeam-theScatteredProton;


	//
	// Finally the diffractive final state,
	// treat at first as a pseudo particle
	//
	TLorentzVector theXparticle((mHadronBeam + mElectronBeam) - (theScatteredElectron + theScatteredProton));

	//
	// The quark and anti quark
	//
	TLorentzVector theQuark;
	TLorentzVector theAntiQuark;
	TLorentzVector theGluon; //for QQG
	double z=event->z;
	double pxq=0., pyq=0., pxqbar=0., pyqbar=0.;
	double mf = Settings::quarkMass(event->quarkSpecies);
	if(fockstate==QQ){
	//
	- // Decay the X-particle using the correct pt
	+ // Decay the X-particle
	//
	- double pt2=z(1-z)Q2-mf*mf;
	- if (4(mfmf+pt2) > MX2) {
	+ if(!decayPseudoParticle(theXparticle, theQuark, theAntiQuark, mf, mf, z, false)){
	if (settings->verboseLevel() > 2) {
	- cout << "InclusiveFinalStateGenerator::generate(): 4(mfmf+pt2) > MX2, for quark flavour "<<event->quarkSpecies << endl;
	- cout << " Abort event generation" << endl;
	+ cout<<"InclusiveFinalStateGenerator::generate: Failed to decay X state, QQ."<<endl;
	}
	return false;
	}
	- if(!decayPseudoParticle(theXparticle, theQuark, theAntiQuark, mf, mf, sqrt(pt2))){
	- if (settings->verboseLevel() > 2) {
	- cout << "InclusiveFinalStateGenerator::generate(): No phase-space for quark pt" << endl;
	- cout << " Abort event generation" << endl;
	- }
	- return false;
	- }
	- if(z<0.5 && theQuark.Pz()<theAntiQuark.Pz()){
	- TLorentzVector tmp=theQuark;
	- theQuark=theAntiQuark;
	- theAntiQuark=tmp;
	- }
	}
	else if(fockstate==QQG){
	- //The knowns:
	double Mqq2=(z/beta-1)*Q2;
	if(Mqq2<4mfmf){
	if (settings->verboseLevel() > 2) {
	cout << "InclusiveFinalStateGenerator::generate(): 4mfmf > Mqq2" << endl;
	cout << " Abort event generation" << endl;
	}
	return false;
	}
	//
	// Given by definitions of z and beta from 2206.13161v1 fig. 12.
	- // Here pomeron comes from + direction, so we reverse the axis. (xi=beta/z)
	- // They make the assumption that the photon does not
	- // carry any +momentum, therefore all is momentum fractions of
	- // the pomeron +momentum. However, for us, the photon does carry
	- // a little bit of +momentum. Therefore we make momentum fractions
	- // of theXparticle instead.
	- //
	- double PXPlus=theXparticle.E()+theXparticle.Pz();
	- double PgPlus=(1-z)*PXPlus; //gluon
	- double PqPlus=beta*PXPlus; //quark
	- double PqbarPlus=(z-beta)*PXPlus; //antiquark
	-
	+ // Here pomeron comes from + direction, so we reverse the axis.
	+ // (xi=beta/z)
	//
	- // First solve for qqbar as a pseudo particle.
	- // Notation: 1:Photon, 2:Pomeron, 3:qqbar, 4:gluon
	- // The unknowns pt3, pt4, P3Minus, P4Minus
	- // Knowns: P1, P2, P3Plus, P4Plus
	- // Once we know the angle between X and g all else solves.
	- // However, finding this angle is tricky.
	- // We use a numerical bisection method.
	- //
	- double phi_totg;
	- bool found=false;
	- // First find the brackets with a rough scan:
	- double valold=evaluate(0., theXparticle, z, Mqq2, MX2, PgPlus, 0.);
	- double valnew=0;
	- double phi_low, phi_high;
	- double stepsize=2*M_PI/100.;
	- for(double phi=0; phi<2*M_PI; phi+=stepsize){
	- valnew=evaluate(phi, theXparticle, z, Mqq2, MX2, PgPlus, 0.);
	- if(valold*valnew==0){
	- phi_totg=phi;
	- cout<<"Found it!"<<endl;
	- found=true;
	- }
	- if(valold*valnew<0){
	- phi_high=phi;
	- phi_low=phi-stepsize;
	- break;
	- }
	- valold=valnew;
	- }
	- if(!found){
	- if(!bisection(phi_high, phi_low, theXparticle, z, Mqq2, MX2, PgPlus, 1e-5, phi_totg, 0.)){
	- // If the bisection fails we do a scan for best value.
	- // It could mean that there is no valid solution...
	- valold=fabs(evaluate(0., theXparticle, z, Mqq2, MX2, PgPlus, 0.));
	- valnew=0;
	- double eps=1e-4;
	- double stepsize=2*M_PI/100.;
	- while(true){
	- phi+=stepsize;
	- if(phi>=2*M_PI){
	- cout<<"no value found"<<endl;
	- return false;
	- }
	- valnew=fabs(evaluate(phi, theXparticle, z, Mqq2, MX2, PgPlus, 0.));
	- if(valnew>valold){
	- stepsize=-stepsize/10.;
	- }
	- if(fabs(stepsize)<eps){
	- phi_totg=phi;
	- break;
	- }
	- valold=valnew;
	- }
	- }
	- }
	- double pt=theXparticle.Pt();
	- double b=2(1-z)cos(phi_totg)*pt;
	- double c=(1-z)(1-z)ptpt+(1-z)Mqq2-z(1-z)MX2;
	- double ptg=-b/2(1-sqrt(1-4c/b/b));
	- double phi_tot=acos(theXparticle.Px()/pt);
	- double PgMinus=ptg*ptg/PgPlus;
	-
	- //Fill the gluon 4-vector:
	- E=(PgPlus+PgMinus)/2;
	- pz=(PgPlus-PgMinus)/2;
	- px=ptg*cos(phi_tot+phi_totg);
	- py=ptg*sin(phi_tot+phi_totg);
	- theGluon=TLorentzVector(px, py, pz, E);
	-
	- //Decay the pseudo particle into qqbar,
	- //Call the pseudo gluon Y
	- //TODO: Check that we get the correct P+ momenta, may have to use the same method again as above...
	- //My original method is not working. For now, therefore,
	- //just decay the pseudo gluon into the qq-bar without
	- //caring about the correct momentum fractions:
	- TLorentzVector thePseudoGluon = theXparticle - theGluon;
	- decayPseudoParticle(thePseudoGluon, theQuark, theAntiQuark, mf, mf);
	-
	-// TLorentzVector thePseudoGluon = theXparticle - theGluon;
	-// double xi=beta/z;
	-// // Angle between q and Y:
	-// double phi_qY;
	-// found=false;
	-// // First find the brackets with a rough scan:
	-// valold=evaluate(0., thePseudoGluon, xi, mfmf, Mqq2, PqPlus, mfmf);
	-// valnew=0;
	-// stepsize=2*M_PI/100.;
	-// for(double phi=0; phi<2*M_PI; phi+=stepsize){
	-// valnew=evaluate(phi, thePseudoGluon, xi, mfmf, Mqq2, PqPlus, mfmf);
	-// if(valold*valnew==0){
	-// phi_qY=phi;
	-// cout<<"Found it!"<<endl;
	-// found=true;
	-// }
	-// if(valold*valnew<0){
	-// phi_high=phi;
	-// phi_low=phi-stepsize;
	-// break;
	-// }
	-// valold=valnew;
	-// }
	-// if(!found){
	-// if(!bisection(phi_high, phi_low, theXparticle, z, Mqq2, MX2, PqPlus, 1e-5, phi_qY, mf*mf)){
	-// cout<<"bisection 2 failed"<<endl;
	-// // If the bisection fails we do a scan for best value.
	-// // It could mean that there is no valid solution...
	-// valold=fabs(evaluate(0., thePseudoGluon, xi, mfmf, Mqq2, PqPlus, mfmf));
	-// valnew=0;
	-// double eps=1e-4;
	-// double stepsize=2*M_PI/10.;
	-// phi=0;
	-// double valmin=valold;
	-// double phimin=0;
	-// double upperBracket=2*M_PI;
	-// while(true){
	-// phi+=stepsize;
	-// valnew=fabs(evaluate(phi, thePseudoGluon, xi, mfmf, Mqq2, PqPlus, mfmf));
	-// if(valnew<valmin){
	-// valmin=valnew;
	-// phimin=phi;
	-// }
	-// if(phi>upperBracket){
	-// phi=phimin-2*stepsize;
	-// upperBracket=phimin+stepsize;
	-// stepsize/=10.;
	-// valmin=valnew;
	-// }
	-// PR(valnew);
	-// PR(phi);
	-// if(fabs(stepsize)<eps){
	-// phi_qY=phi;
	-// break;
	-// }
	-// valold=valnew;
	-// }
	-// PR(valnew/mf/mf);
	-// }
	-// }
	-// pt=thePseudoGluon.Pt();
	-// b=2(1-xi)cos(phi_qY)*pt;
	-// c=(1-xi)(1-xi)ptpt+mfmf-xi(1-xi)Mqq2;
	-// double ptq=-b/2(1-sqrt(1-4c/b/b));
	-// phi_tot=acos(thePseudoGluon.Px()/pt);
	-// double PqMinus=ptq*ptq/PqPlus;
	-//
	-// //Fill the quark 4-vector:
	-// E=(PqPlus+PqMinus)/2;
	-// pz=(PqPlus-PqMinus)/2;
	-// px=ptq*cos(phi_tot+phi_qY);
	-// py=ptq*sin(phi_tot+phi_qY);
	-// theQuark=TLorentzVector(px, py, pz, E);
	-// theAntiQuark=thePseudoGluon-theQuark;
	-
	-// double k2=(thePomeron-theGluon).M2();
	-// double pt2=xi(1-xi)(-k2)-mf*mf;
	-// if(!decayPseudoParticle(thePseudoGluon, theQuark, theAntiQuark, mf, mf, sqrt(pt2))){
	-// if (settings->verboseLevel() > 2) {
	-// cout << "InclusiveFinalStateGenerator::generate(): No phase-space for quark pt" << endl;
	-// cout << " Abort event generation" << endl;
	-// }
	-// return false;
	-// }
	-/*
	-
	- phi = rndm->Uniform(twopi);
	- pxq=sqrt(pt2)*sin(phi);
	- pyq=sqrt(pt2)*cos(phi);
	- pxqbar=-pxq;
	- pyqbar=-pyq;
	- if (4(mfmf+pt2) > Mqq2) {
	- if (settings->verboseLevel() > 1) {
	- cout << "InclusiveFinalStateGenerator::generate(): 4(mfmf+pt2) > Mqq2" << endl;
	- cout << " Abort event generation" << endl;
	+ TLorentzVector thePseudoGluon;
	+ if(!decayPseudoParticle(theXparticle, theGluon, thePseudoGluon, 0., sqrt(Mqq2), 1-z)){
	+ if (settings->verboseLevel() > 2) {
	+ cout<<"InclusiveFinalStateGenerator::generate: Failed to decay X state."<<endl;
	}
	return false;
	}

	//
	- // Calculate the quark momenta in the X-particle's restframe
	+ // Decay the pseudo particle into qqbar
	//
	- TVector3 boost=TVector3(thePseudoGluon.Px()/thePseudoGluon.E(), thePseudoGluon.Py()/thePseudoGluon.E(), thePseudoGluon.Pz()/thePseudoGluon.E());
	- int sign=1.;
	- if(xi<0.5) sign=-1.;
	- theQuark=TLorentzVector(pxq, pyq, signsqrt(Mqq2/4-mfmf-pt2), sqrt(Mqq2)/2);
	- theAntiQuark=TLorentzVector(pxqbar, pyqbar, -signsqrt(Mqq2/4-mfmf-pt2), sqrt(Mqq2)/2);
	- theQuark.Boost(boost);
	- theAntiQuark.Boost(boost);
	-// PR(PpomPlus);
	- PR(PgPlus);
	- PR((theGluon.E()+theGluon.Pz()));
	- PR(PqPlus);
	- PR(theQuark.E()+theQuark.Pz());
	- PR(PqbarPlus);
	- PR(theAntiQuark.E()+theAntiQuark.Pz());
	- */
	- }
	+ double xi=beta/z;
	+ if(!decayPseudoParticle(thePseudoGluon, theQuark, theAntiQuark, mf, mf, xi)){
	+ if (settings->verboseLevel() > 2) {
	+ cout<<"InclusiveFinalStateGenerator::generate: Failed to decay pseudo-gluon."<<endl;
	+ }
	+ return false;
	+ }
	+ } //QQG
	else{
	cout<<"InclusiveFinalStateGenerator::generate: unknown fockstate: "<<fockstate<<" ending program."<<endl;
	exit(0);
	}
	// Generate decorrelation between the qqbar system.
	// This we do by deconvoluting the pomeron as 2T gluons
	// where T is the twist that we are at. These gluons
	// are assumed to have momenta distributed as a Gaussian
	// with width Qs around the momentum of the pomeron
	// in transverse space.
	//
	// We will take a random walk while keeping both particles
	// on shell, as well as (P1+P1)^2=MX2
	//
	// 1. First find which twist we are operating at,
	//
	int Twist=0;
	double Omega=event->Omega;
	double eps=1e-3;
	double sum=0;
	for(int n=1; n<100; n++){
	int nfactorial=1;
	for(int i=1; i<=n; i++){
	nfactorial*=i;
	}
	sum+=pow(-1, n-1)2/nfactorialpow(Omega/2, n);
	double trueVal=2*(1-exp(-Omega/2));
	if(abs(sum-trueVal)<eps){
	Twist=n;
	break;
	}
	}
	- /*
	- Twist=0;
	- for(int i=0; i<100; i++){
	- double gammaN=gamma_N(i+1, -Omega/2);
	- if(fabs(2exp(-Omega/2)gammaN)<eps){
	- Twist=i+1;
	- break;
	- }
	- }
	- PR(Twist);
	- exit(1);
	- */
	//
	// 2. Now generate 2*Twist number of gluons:
	//
	double Qs=event->QS;
	vector<double> pom_px, pom_py, pom_pz;
	double pom_px_tot=thePomeron.Px();
	double pom_py_tot=thePomeron.Py();
	double pom_pz_tot=thePomeron.Pz();
	for(int i=0; i<2*Twist-1; i++){ //-1 since the last gluon conserves total momentum
	//
	// Generate momenta from a Gaussian with width Qs
	//
	double p_x=rndm->Gaus(0., Qs);
	double p_y=rndm->Gaus(0., Qs);
	double p_z=rndm->Gaus(0., Qs);
	//
	// Distribute the individual gluon momenta around
	// the pomeron momentum (the pomeron is the combined
	// particle from all the exchange gluons
	//
	pom_px.push_back(thePomeron.Px()/(2*Twist)+p_x);
	pom_py.push_back(thePomeron.Py()/(2*Twist)+p_y);
	pom_pz.push_back(thePomeron.Pz()/(2*Twist)+p_z);
	pom_px_tot+=p_x;
	pom_py_tot+=p_y;
	pom_pz_tot+=p_z;
	}
	//conserve the total pomeron momentum:
	pom_px.push_back(thePomeron.Px()-pom_px_tot);
	pom_py.push_back(thePomeron.Py()-pom_py_tot);
	pom_pz.push_back(thePomeron.Pz()-pom_pz_tot);
	//
	// 3. Distribute the recoils to:
	// the quark and antiquark (QQ) or
	// the quark, antiquarkgluon and gluon (QQG).
	// Start with QQ
	//
	double Eq=0, Eqbar=0, pzq=0, pzqbar=0;
	if(fockstate==QQ){
	Eq=theQuark.E();
	Eqbar=theAntiQuark.E();
	pzq=theQuark.Pz();
	pzqbar=theAntiQuark.Pz();
	pxq=theQuark.Px();
	pyq=theQuark.Py();
	pxqbar=theAntiQuark.Px();
	pyqbar=theAntiQuark.Py();
	for(unsigned int i=0; i<pom_px.size(); i++){
	//
	//First calculate the relative velocity between q, qbar and gluon
	//
	double pomP=sqrt(pom_px[i]pom_px[i] + pom_py[i]pom_py[i] + pom_pz[i]*pom_pz[i]);
	double qP=sqrt(pxqpxq + pyqpyq + pzq*pzq);
	double qbarP=sqrt(pxqbarpxqbar + pyqbarpyqbar + pzqbar*pzqbar);
	double deltaVq=sqrt(1+qPqP/Eq/Eq - 2(pom_px[i]pxq + pom_py[i]pyq + pom_pz[i]*pzq) /Eq/pomP);
	double deltaVqbar=sqrt(1+qbarPqbarP/Eqbar/Eqbar - 2(pom_px[i]pxqbar + pom_py[i]pyqbar + pom_pz[i]*pzqbar) / Eqbar/pomP);
	//The momentum kick is inversely proportional to relative velocity
	double fraction=0.5;
	if(deltaVq+deltaVqbar!=0)
	fraction=deltaVqbar/(deltaVq+deltaVqbar);
	//Update momenta
	pxq+=fraction*pom_px[i];
	pyq+=fraction*pom_py[i];
	pzq+=fraction*pom_pz[i];
	pxqbar+=(1-fraction)*pom_px[i];
	pyqbar+=(1-fraction)*pom_py[i];
	pzqbar+=(1-fraction)*pom_pz[i];

	//Scale the transverse momenta in order to preserve MX2,
	//First find the scaling factor lambda
	double lower_bracket=0, upper_bracket=100;
	double lambda=computeScalingFactor(pxq, pyq, pzq, pxqbar, pyqbar, pzqbar, mf, mf, sqrt(MX2), lower_bracket, upper_bracket);
	if(lambda<=lower_bracket+1e-5 \|\| lambda>=upper_bracket-1e-5) {
	//failed to find the root
	if (settings->verboseLevel() > 2) {
	cout<<"InclusiveFinalStateGenerator::generate: Failed to scatter the quarks"<<endl;
	}
	return false;
	}
	if(lambda==(lower_bracket+upper_bracket)/2)
	PR(lambda);
	pxq*=lambda;
	pyq*=lambda;
	pzq*=lambda;
	pxqbar*=lambda;
	pyqbar*=lambda;
	pzqbar*=lambda;
	//Put on shell:
	Eq=sqrt(pxqpxq+pyqpyq+pzqpzq+mfmf);
	Eqbar=sqrt(pxqbarpxqbar+pyqbarpyqbar+pzqbarpzqbar+mfmf);
	}
	//
	//Update the 4-vectors:
	//
	//the Quark:
	theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
	//the Anti-Quark
	theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
	}
	else if(fockstate==QQG){
	Eq=theQuark.E();
	pzq=theQuark.Pz();
	pxq=theQuark.Px();
	pyq=theQuark.Py();
	Eqbar=theAntiQuark.E();
	pzqbar=theAntiQuark.Pz();
	pxqbar=theAntiQuark.Px();
	pyqbar=theAntiQuark.Py();
	double Eg=theGluon.E();
	double pzg=theGluon.Pz();
	double pxg=theGluon.Px();
	double pyg=theGluon.Py();

	for(unsigned int i=0; i<pom_px.size(); i++){
	//
	//First calculate the relative velocity between q, qbar and gluon
	//
	double pomP=sqrt(pom_px[i]pom_px[i] + pom_py[i]pom_py[i] + pom_pz[i]*pom_pz[i]);
	double qP=sqrt(pxqpxq + pyqpyq + pzq*pzq);
	double qbarP=sqrt(pxqbarpxqbar + pyqbarpyqbar + pzqbar*pzqbar);
	double gP=sqrt(pxgpxg + pygpyg + pzg*pzg);
	double deltaVg=sqrt(1+gPgP/Eg/Eg - 2(pom_px[i]pxg + pom_py[i]pyg + pom_pz[i]*pzg) /Eg/pomP);
	//Figure out which dipole the gluon is interacting with,
	// g-q or q-qbar:
	double p2_gq=(pxq+pxg)(pxq+pxg) + (pyq+pyg)(pyq+pyg) + (pzq+pzg)*(pzq+pzg);
	double p2_gqbar=(pxqbar+pxg)(pxqbar+pxg) + (pyqbar+pyg)(pyqbar+pyg) + (pzqbar+pzg)*(pzqbar+pzg);
	//Interaction probability propto dipole size:
	bool isGQbar = p2_gq/(p2_gq+p2_gqbar) > rndm->Uniform(1);
	if(isGQbar){
	//Interacting with the g-qbar dipole:
	double deltaVqbar=sqrt(1+qbarPqbarP/Eqbar/Eqbar - 2(pom_px[i]pxqbar + pom_py[i]pyqbar + pom_pz[i]*pzqbar) / Eqbar/pomP);
	//The momentum kick is inversely proportional to relative velocity
	double fraction=0.5;
	if(deltaVg+deltaVqbar!=0)
	fraction=deltaVg/(deltaVg+deltaVqbar);
	//Update momenta
	pxg+=fraction*pom_px[i];
	pyg+=fraction*pom_py[i];
	pzg+=fraction*pom_pz[i];
	pxqbar+=(1-fraction)*pom_px[i];
	pyqbar+=(1-fraction)*pom_py[i];
	pzqbar+=(1-fraction)*pom_pz[i];
	}
	else{
	//Interacting with the g-q dipole:
	double deltaVq=sqrt(1+qPqP/Eq/Eq - 2(pom_px[i]pxq + pom_py[i]pyq + pom_pz[i]*pzq) /Eq/pomP);
	//The momentum kick is inversely proportional to relative velocity
	double fraction=0.5;
	if(deltaVg+deltaVq!=0)
	fraction=deltaVq/(deltaVq+deltaVg);
	//Update momenta
	pxg+=fraction*pom_px[i];
	pyg+=fraction*pom_py[i];
	pzg+=fraction*pom_pz[i];
	pxq+=(1-fraction)*pom_px[i];
	pyq+=(1-fraction)*pom_py[i];
	pzq+=(1-fraction)*pom_pz[i];
	}

	//Scale the transverse momenta in order to preserve MX2 and Mqq2,
	//First find the scaling factor lambda
	// start by scaling the quarks:
	double Mqq2=(z/beta-1)*Q2;
	double lower_bracket=0, upper_bracket=100;
	double lambda=computeScalingFactor(pxq, pyq, pzq, pxqbar, pyqbar, pzqbar, mf, mf, sqrt(Mqq2), lower_bracket, upper_bracket);
	if(lambda<=lower_bracket+1e-5 \|\| lambda>=upper_bracket-1e-5) {
	if (settings->verboseLevel() > 2) {
	cout<<"InclusiveFinalStateGenerator::generate: Failed to scatter the gluons"<<endl;
	}
	return false;
	}
	pxq*=lambda;
	pyq*=lambda;
	pzq*=lambda;
	pxqbar*=lambda;
	pyqbar*=lambda;
	pzqbar*=lambda;
	//Put on shell:
	Eq=sqrt(pxqpxq+pyqpyq+pzqpzq+mfmf);
	Eqbar=sqrt(pxqbarpxqbar+pyqbarpyqbar+pzqbarpzqbar+mfmf);

	//
	// Create the pseudo gluon based on the current quarks
	//
	theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
	theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
	TLorentzVector thePseudoGluon=theQuark+theAntiQuark;

	double pxpg=thePseudoGluon.Px();
	double pypg=thePseudoGluon.Py();
	double pzpg=thePseudoGluon.Pz();
	double Mpg=thePseudoGluon.M();
	// Calculate momentum fractions for later:
	double pxFraction=pxq/pxpg;
	double pyFraction=pyq/pypg;
	double pzFraction=pzq/pzpg;

	lower_bracket=0; upper_bracket=100;
	lambda=computeScalingFactor(pxpg, pypg, pzpg, pxg, pyg, pzg, Mpg, 0, sqrt(MX2), lower_bracket, upper_bracket);
	if(lambda<=lower_bracket+1e-5 \|\| lambda>=upper_bracket-1e-5) {
	if (settings->verboseLevel() > 2) {
	cout<<"InclusiveFinalStateGenerator::generate: Failed to scatter the quarks"<<endl;
	}
	return false;
	}
	if(lambda==(lower_bracket+upper_bracket)/2)
	pxpg*=lambda;
	pypg*=lambda;
	pzpg*=lambda;

	pxq=pxFraction*pxpg;
	pyq=pyFraction*pypg;
	pzq=pzFraction*pzpg;
	Eq=sqrt(pxqpxq+pyqpyq+pzqpzq+mfmf);

	pxqbar=(1-pxFraction)*pxpg;
	pyqbar=(1-pyFraction)*pypg;
	pzqbar=(1-pzFraction)*pzpg;
	Eqbar=sqrt(pxqbarpxqbar+pyqbarpyqbar+pzqbarpzqbar+mfmf);

	pxg*=lambda;
	pyg*=lambda;
	pzg*=lambda;
	Eg=sqrt(pxgpxg+pygpyg+pzg*pzg);
	}
	//
	//Update the 4-vectors:
	//
	//the Quark:
	theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
	//the Anti-Quark
	theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
	//the Gluon:
	theGluon=TLorentzVector(pxg, pyg, pzg, Eg);
	}
	else{
	cout<<"InclusiveFinalStateGenerator::generate: unknown fockstate: "<<fockstate<<" ending program."<<endl;
	exit(9);
	}
	//
	// Decorrelations ended
	//

	//
	// Add particles to event record
	//
	double ifock;
	if(fockstate==QQ)
	ifock=0;
	else
	ifock=1;
	event->particles.resize(2+5+ifock);
	unsigned int eOut = 2;
	unsigned int gamma = 3;
	unsigned int hOut = 4;
	unsigned int quark = 5;
	unsigned int antiquark = 6;
	unsigned int gluon = 7;

	// Global indices
	event->particles[eOut].index = eOut;
	event->particles[gamma].index = gamma;
	event->particles[hOut].index = hOut;
	event->particles[quark].index = quark;
	event->particles[antiquark].index = antiquark;
	if(fockstate==QQG)
	event->particles[antiquark].index = gluon;

	// 4-vectors
	event->particles[eOut].p = theScatteredElectron;
	event->particles[hOut].p = theScatteredProton;
	event->particles[gamma].p = theVirtualPhoton;
	event->particles[quark].p = theQuark;
	event->particles[antiquark].p = theAntiQuark;
	if(fockstate==QQG)
	event->particles[gluon].p = theGluon;

	// PDG Ids
	event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
	event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
	event->particles[gamma].pdgId = 22;
	event->particles[quark].pdgId = event->quarkSpecies+1;
	event->particles[antiquark].pdgId = -event->particles[quark].pdgId;
	if(fockstate==QQG)
	event->particles[gluon].pdgId = 21;

	// status
	//
	// HepMC conventions (February 2009).
	// 0 : an empty entry, with no meaningful information
	// 1 : a final-state particle, i.e. a particle that is not decayed further by
	// the generator (may also include unstable particles that are to be decayed later);
	// 2 : a decayed hadron or tau or mu lepton
	// 3 : a documentation entry (not used in PYTHIA);
	// 4 : an incoming beam particle;
	// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
	// fulfill the criteria of status code 2

	event->particles[ePos].status = 4;
	event->particles[hPos].status = 4;
	event->particles[eOut].status = 1;
	event->particles[hOut].status = mIsIncoherent ? 2 : 1;
	event->particles[gamma].status = 2;
	event->particles[quark].status = 2;
	event->particles[antiquark].status = 2;
	if(fockstate==QQG)
	event->particles[gluon].status = 2;

	// parents (ignore dipole)
	event->particles[eOut].parents.push_back(ePos);
	event->particles[gamma].parents.push_back(ePos);
	event->particles[hOut].parents.push_back(hPos);
	event->particles[quark].parents.push_back(gamma);
	event->particles[antiquark].parents.push_back(gamma);
	if(fockstate==QQG){
	event->particles[gluon].parents.push_back(quark);
	event->particles[gluon].parents.push_back(antiquark);
	}

	// daughters (again ignore dipole)
	event->particles[ePos].daughters.push_back(eOut);
	event->particles[ePos].daughters.push_back(gamma);
	event->particles[gamma].daughters.push_back(quark);
	event->particles[gamma].daughters.push_back(antiquark);
	event->particles[hPos].daughters.push_back(hOut);
	if(fockstate==QQG){
	event->particles[quark].daughters.push_back(gluon);
	event->particles[antiquark].daughters.push_back(gluon);
	}


	//fill event structure
	double y=mHadronBeamtheVirtualPhoton/(mHadronBeammElectronBeam);
	W2=(theVirtualPhoton+mHadronBeam).M2();
	mQ2=-theVirtualPhoton.M2();
	MX2=(theQuark+theAntiQuark).M2();
	double Delta=thePomeron.Pt();

	event->t=-Delta*Delta;
	event->Q2=mQ2;
	event->W=sqrt(W2);
	event->y=y;
	event->MX=sqrt(MX2);

	//
	// Now the afterburner 'Pythia8' part for hadronization
	//

	mPythiaEvent->reset();

	// Pythia starts with u = 1, Sartre with u = 0
	int quarkID = event->quarkSpecies + 1;
	int status = 23;
	int color = 101; // 101 102 103

	if(fockstate==QQ){
	mPythiaEvent->append( quarkID, status, color, 0,
	theQuark.Px(), theQuark.Py(), theQuark.Pz(),
	theQuark.E(), theQuark.M() );
	mPythiaEvent->append(-quarkID, status, 0, color,
	theAntiQuark.Px(), theAntiQuark.Py(), theAntiQuark.Pz(),
	theAntiQuark.E(), theAntiQuark.M());
	}
	else if(fockstate==QQG){
	int acolor=102;
	mPythiaEvent->append( quarkID, status, color, 0,
	theQuark.Px(), theQuark.Py(), theQuark.Pz(),
	theQuark.E(), theQuark.M() );
	mPythiaEvent->append(-quarkID, status, 0, acolor,
	theAntiQuark.Px(), theAntiQuark.Py(), theAntiQuark.Pz(),
	theAntiQuark.E(), theAntiQuark.M());
	mPythiaEvent->append(21, status, acolor, color,
	theGluon.Px(), theGluon.Py(), theGluon.Pz(),
	theGluon.E(), theGluon.M());
	}
	else{
	cout<<"Fockstate is not recognised, ending program."<<endl;
	exit(1);
	}
	//
	// Generate event. Abort on failure.
	//
	cout.setstate(std::ios::failbit); // suppress all Pythia print-out
	bool ok = mPythia->next();
	cout.clear();

	if (!ok \|\| mPythiaEvent->size() == 0) {
	if (settings->verbose() && settings->verboseLevel() > 2) {
	cout << "InclusiveFinalStateGenerator::generate(): Pythia8 event generation aborted prematurely, owing to error!" << endl;
	}
	return false;
	}
	// mPythiaEvent->list(); //Print out the pythia event to the screen
	//
	// Loop over pythia particles and fill our particle list
	//
	int offset = 3+ifock; //pythia particles start at i=3(QQ), or i=4(QQG) in pythia list

	event->numberOfPythiaParticlesStored = 0;
	for (int i = offset; i < mPythiaEvent->size(); i++) {
	Particle pypart;
	pypart.index = event->particles.size();
	pypart.pdgId = mPythiaEvent->at(i).id();
	if (mPythiaEvent->at(i).status() < 0)
	pypart.status = 2;
	else
	pypart.status = 1;

	pypart.p = TLorentzVector(mPythiaEvent->at(i).px(),
	mPythiaEvent->at(i).py(),
	mPythiaEvent->at(i).pz(),
	mPythiaEvent->at(i).e());
	if (mPythiaEvent->at(i).mother1()) pypart.parents.push_back(mPythiaEvent->at(i).mother1()+offset+1);
	if (mPythiaEvent->at(i).mother2()) pypart.parents.push_back(mPythiaEvent->at(i).mother2()+offset+1);
	if (mPythiaEvent->at(i).daughter1()) pypart.daughters.push_back(mPythiaEvent->at(i).daughter1()+offset+1);
	if (mPythiaEvent->at(i).daughter2()) pypart.daughters.push_back(mPythiaEvent->at(i).daughter2()+offset+1);

	event->particles.push_back(pypart);
	event->numberOfPythiaParticlesStored++;
	}

	return true;
	}

	// Function to calculate the scaling factor lambda numerically (from ChatGPT)
	double InclusiveFinalStateGenerator::computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m_1, double m_2, double M, double lambda_min=0, double lambda_max=100) {
	// Define the function to compute the total energy after scaling
	auto totalEnergy = [&](double lambda) {
	double E1_prime = std::sqrt(lambdalambda(px1px1 + py1py1 + pz1pz1) + m_1m_1);
	double E2_prime = std::sqrt(lambdalambda(px2px2 + py2py2 + pz2pz2) + m_2m_2);
	return E1_prime + E2_prime;
	};

	// Function to compute the invariant mass condition for a given lambda
	auto invariantMassCondition = [&](double lambda) {
	double px_sum = lambda * (px1 + px2);
	double py_sum = lambda * (py1 + py2);
	double pz_sum = lambda * (pz1 + pz2);
	double energy_sum = totalEnergy(lambda);

	return energy_sumenergy_sum - px_sumpx_sum - py_sumpy_sum - pz_sumpz_sum - M*M;
	};

	// Solve for lambda using a numerical method (bisection method for simplicity)
	double tolerance = 1e-6;
	while (lambda_max - lambda_min > tolerance) {
	double lambda_mid = 0.5 * (lambda_min + lambda_max);
	if (invariantMassCondition(lambda_mid) > 0)
	lambda_max = lambda_mid;
	else
	lambda_min = lambda_mid;
	}

	return 0.5 * (lambda_min + lambda_max);
	}
	-
	-// Function that computes f(phi_totg) = (pseudo gluon mass)^2 - Mqq2
	-double InclusiveFinalStateGenerator::evaluate(double phi_totg, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus, double mg2=0) {
	- double pt = theXparticle.Pt();
	- double b = 2 * (1 - z) * cos(phi_totg) * pt;
	- double c = (1 - z) * (1 - z) * pt * pt + (1 - z) * Mqq2 - z * (1 - z) * MX2 + z * mg2;
	- double discriminant = 1 - 4 * c / (b * b);
	-
	- if (discriminant < 0) {
	- return 1e9; // Large value if not physical
	- }
	-
	- double ptg = -b / 2 * (1 - sqrt(discriminant));
	- double phi_tot = acos(theXparticle.Px() / pt);
	-
	- double PgMinus = ptg * ptg / PgPlus;
	- double E = (PgPlus + PgMinus) / 2;
	- double pz = (PgPlus - PgMinus) / 2;
	- double px = ptg * cos(phi_tot + phi_totg);
	- double py = ptg * sin(phi_tot + phi_totg);
	-
	- TLorentzVector theGluon(px, py, pz, E);
	- TLorentzVector thePseudoGluon = theXparticle - theGluon;
	-
	- return thePseudoGluon.M2() - Mqq2;
	-}
	-
	-bool InclusiveFinalStateGenerator::bisection(double phi_min, double phi_max, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus,
	- double tol, double& phi_root, double mg2=0) {
	-
	- double f_min = evaluate(phi_min, theXparticle, z, Mqq2, MX2, PgPlus, mg2);
	- double f_max = evaluate(phi_max, theXparticle, z, Mqq2, MX2, PgPlus, mg2);
	-
	- if (f_min * f_max > 0) {
	- // No sign change, cannot use bisection
	- return false;
	- }
	-
	- int max_iter = 100;
	- for (int iter = 0; iter < max_iter; ++iter) {
	- double phi_mid = 0.5 * (phi_min + phi_max);
	- double f_mid = evaluate(phi_mid, theXparticle, z, Mqq2, MX2, PgPlus);
	-
	- if (std::abs(f_mid) < tol) {
	- phi_root = phi_mid;
	- return true;
	- }
	-
	- if (f_min * f_mid < 0) {
	- phi_max = phi_mid;
	- f_max = f_mid;
	- } else {
	- phi_min = phi_mid;
	- f_min = f_mid;
	- }
	- }
	-
	- // If we get here, didn't converge
	- return false;
	-}
	-
	-
	Index: trunk/src/InclusiveFinalStateGenerator.h
	===================================================================
	--- trunk/src/InclusiveFinalStateGenerator.h (revision 550)
	+++ trunk/src/InclusiveFinalStateGenerator.h (revision 551)
	@@ -1,53 +1,47 @@
	//==============================================================================
	// ExclusiveFinalStateGenerator.h
	//
	// Copyright (C) 2024 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: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
	// $Author: ullrich $
	//==============================================================================
	#ifndef InclusiveFinalStateGenerator_h
	#define InclusiveFinalStateGenerator_h
	#include "Pythia8/Pythia.h"
	#include "FinalStateGenerator.h"
	#include "Constants.h"
	#include "Enumerations.h"

	class InclusiveFinalStateGenerator : public FinalStateGenerator {
	public:
	InclusiveFinalStateGenerator();
	~InclusiveFinalStateGenerator();

	bool generate(int A, Event *event, FockState fockstate);
	-
	- double uiGamma_N(double* var, double* par) const;
	- double gamma_N(unsigned int i, double val);
	-
	+
	private:
	Pythia8::Pythia *mPythia;
	Pythia8::Event *mPythiaEvent;
	Pythia8::ParticleData *mPdt;

	double computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m_1, double m_2, double M, double, double);

	- bool decayPseudoParticle(TLorentzVector, TLorentzVector&, TLorentzVector&, double, double, double);
	+ bool decayPseudoParticle(TLorentzVector, TLorentzVector&, TLorentzVector&, double, double, double, bool);

	- double evaluate(double phi_totg, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus, double mf2);
	- bool bisection(double phi_min, double phi_max, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus,
	- double tol, double& phi_root, double mf2);
	};
	#endif

File Metadata

Mime Type: application/octet-stream
Expires: Wed, May 14, 10:52 AM (2 d)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111280
Default Alt Text: (84 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions