Page MenuHomeHEPForge
Files F19244622

InclusiveFinalStateGenerator.cpp
No OneTemporary
Actions

Size
15 KB
Referenced Files
None
Subscribers
None

InclusiveFinalStateGenerator.cpp
View Options

//==============================================================================
// 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"
#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)
//
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");
if (!mPythia->init()) {
cout << "Error initializing Pythia8. Stop here." << endl;
exit(1);
}
else {
cout << "Pythia8 initialized (" << PYTHIA_VERSION_INTEGER << ")" << endl;
}
}
InclusiveFinalStateGenerator::~InclusiveFinalStateGenerator()
{
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::generate(int A, Event *event)
{
//
// 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;
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: me*me=E'*E'-pt'*pt'-pz'*pz'
// 2: Q2=-(e-e')^2=-2*me*me + 2*(E*E'-pz*pz')
// 3: W2=(P+e-e')^2=mp2+2*me2+2*(Ep*E-Pz*pz)-2*(Ep*E'-Pz*pz')-2*(E*E'-pz*pz')
//
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-2*Ee*Ep) + (Pp+Pe)*Q2 + 2*Pe*Pe*Pp + 2*me2*Pp;
E /= 2*(Ee*Pp-Ep*Pe);
pz = Ee*(W2-mp2) + (Ep+Ee)*Q2 + 2*Ee*Pe*Pp + 2*Ep*me2 - 2*Ee*Ee*Ep;
pz /= 2*(Ee*Pp-Ep*Pe);
//
// Equation 1:
//
pt = sqrt(E*E-pz*pz-me2);
phi = rndm->Uniform(twopi);
TLorentzVector theScatteredElectron(pt*sin(phi), pt*cos(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
//
E=(1-xpom)*mHadronBeam.Pz()*mHadronBeam.Pz()-mT/2.+0.5*hMass2+0.5*mMY2;
E=E/mHadronBeam.E();
pz = (1-xpom)*mHadronBeam.Pz();
pt = sqrt(E*E-pz*pz-mMY2);
if(pt!=pt){
if(settings->verboseLevel() > 2)
cout<<"No phase-space for scattered proton pt"<<endl;
return false;
}
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
//
double mf = quarkMass[event->quarkSpecies];
double z=event->z;
double pt2=z*(1-z)*Q2-mf*mf;
if(4*(mf*mf+pt2)/MX2>1) pt2=MX2/4.-mf*mf;
//generate angle
phi = rndm->Uniform(twopi);
double pxq=sqrt(pt2)*sin(phi);
double pyq=sqrt(pt2)*cos(phi);
double pxqbar=-pxq;
double pyqbar=-pyq;
double z0=0.5*(1-sqrt(1-4*(mf*mf+pt2)/MX2));
if(mElectronBeam.Pz()<0 and z<0.5) z0=1-z0;
if(mElectronBeam.Pz()>0 and z>0.5) z0=1-z0;
double qPplus=z0*(theXparticle.E()+theXparticle.Pz());
double qPminus=(1-z0)*(theXparticle.E()-theXparticle.Pz());
double qBarPplus=(1-z0)*(theXparticle.E()+theXparticle.Pz());
double qBarPminus=z0*(theXparticle.E()-theXparticle.Pz());
//
// Generate decorrelation between the qqbar system.
//
//
// First find which twist we are operating at
//
int Twist=0;
double Omega=event->Omega;
double eps=1e-3;
for(int i=0; i<100; i++){
double gammaN=gamma_N(i+1, -Omega/2);
if(fabs(2*exp(-Omega/2)*gammaN)<eps){
Twist=i+1;
break;
}
}
//
// Now generate 2*Twist number of gluons:
// Do they need a pz as well??
//
double eps2=z*(1-z)*(Q2+MX2);
double r_average=1./sqrt(eps2);
double Qs=sqrt(2)*Omega/r_average;
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++){
//#TT This part does not seem quite right yet.
double p_x=rndm->Gaus(0., Qs);
double p_y=rndm->Gaus(0., Qs);
double p_z=rndm->Gaus(0., Qs);
pom_px.push_back(thePomeron.Px()+p_x);
pom_py.push_back(thePomeron.Py()+p_y);
pom_pz.push_back(thePomeron.Pz()+p_z);
pom_px_tot+=px;
pom_py_tot+=py;
pom_pz_tot+=pz;
}
//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);
//Distribute the recoils to the quark and antiquark:
double Eq=(qPplus+qPminus)/2;
double Eqbar=(qBarPplus+qBarPminus)/2;
double pzq=(qPplus-qPminus)/2;
double pzqbar=(qBarPplus-qBarPminus)/2;
for(unsigned int i=0; i<pom_px.size(); i++){
//first calculate the relative velocity between q, qbar and pomeron
double pomP=sqrt(pom_px[i]*pom_px[i]+pom_py[i]*pom_py[i]+pom_pz[i]*pom_pz[i]);
double qP=sqrt(pxq*pxq+pyq*pyq+pzq*pzq);
double qbarP=sqrt(pxqbar*pxqbar+pyqbar*pyqbar+pzqbar*pzqbar);
double deltaVq=sqrt(1+qP*qP/Eq/Eq-2*(pom_px[i]*pxq+pom_py[i]*pyq+pom_pz[i]*pzq)/Eq/pomP);
double deltaVqbar=sqrt(1+qbarP*qbarP/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];
pxqbar+=(1-fraction)*pom_px[i];
pyqbar+=(1-fraction)*pom_py[i];
//Update the energies
Eq=sqrt(pxq*pxq+pyq*pyq+pzq*pzq+mf*mf);
Eqbar=sqrt(pxqbar*pxqbar+pyqbar*pyqbar+pzqbar*pzqbar+mf*mf);
}
//Fill the 4-vectors:
//the Quark:
TLorentzVector theQuark(pxq, pyq, pzq, Eq);
//the Anti-Quark
TLorentzVector theAntiQuark(pxqbar, pyqbar, pzqbar, Eqbar);
theXparticle=theQuark+theAntiQuark;
//
// Add particles to event record
//
event->particles.resize(2+7);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int X = 4;
unsigned int pomeron = 5;
unsigned int hOut = 6;
unsigned int quark = 7;
unsigned int antiquark = 8;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[X].index = X;
event->particles[pomeron].index = pomeron;
event->particles[hOut].index = hOut;
event->particles[quark].index = quark;
event->particles[antiquark].index = antiquark;
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[X].p = theXparticle;
event->particles[pomeron].p = thePomeron;
event->particles[quark].p = theQuark;
event->particles[antiquark].p = theAntiQuark;
// 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[X].pdgId = 90; //pseudoparticle
int iquark=event->quarkSpecies;
if(iquark<=1)
event->particles[pomeron].pdgId = 113;
if(iquark==2)
event->particles[pomeron].pdgId = 333;
if(iquark==3)
event->particles[pomeron].pdgId = 443;
if(iquark==4)
event->particles[pomeron].pdgId = 553;
event->particles[quark].pdgId = iquark+1;
event->particles[antiquark].pdgId = -event->particles[quark].pdgId;
// 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[X].status = 90;
event->particles[pomeron].status = 2;
event->particles[quark].status = 1;
event->particles[antiquark].status = 1;
// 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[hOut].parents.push_back(pomeron);
event->particles[pomeron].parents.push_back(gamma);
event->particles[pomeron].parents.push_back(gamma);
event->particles[X].parents.push_back(gamma);
event->particles[quark].parents.push_back(gamma);
event->particles[antiquark].parents.push_back(gamma);
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(X);
event->particles[hPos].daughters.push_back(pomeron);
event->particles[gamma].daughters.push_back(quark);
event->particles[gamma].daughters.push_back(antiquark);
event->particles[pomeron].daughters.push_back(hOut);
event->particles[hPos].daughters.push_back(hOut);
//fill event structure
double y=mHadronBeam*theVirtualPhoton/(mHadronBeam*mElectronBeam);
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
//
int quarkID = event->quarkSpecies + 1; // Pythia starts with u = 1, Sartre with u = 0
int status = 23;
int color = 101; // 101 102 103
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());
// Generate events. Quit if failure.
if (!mPythia->next()) {
cout << "InclusiveFinalStateGenerator::generate(): Pythia8 event generation aborted prematurely, owing to error!\n";
return false;
}
// Loop over pythia particles and fill our particle list
mPythiaEvent->list();
return true;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 4:44 AM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6564796
Default Alt Text
InclusiveFinalStateGenerator.cpp (15 KB)

Event Timeline

HEPForge.org · .