Page MenuHomeHEPForge
Files F8723628

No OneTemporary
Actions

View Options

Index: trunk/src/Enumerations.h
===================================================================
--- trunk/src/Enumerations.h (revision 532)
+++ trunk/src/Enumerations.h (revision 533)
@@ -1,33 +1,35 @@
//==============================================================================
// Enumerations.h
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#ifndef Enumerations_h
#define Enumerations_h
enum DipoleModelType {bSat, bNonSat, bCGC};
enum GammaPolarization {transverse, longitudinal};
enum AmplitudeMoment {mean_A, mean_A2, variance_A, lambda_real, lambda_skew};
enum DiffractiveMode {coherent, incoherent};
enum DipoleModelParameterSet {KMW, HMPZ, STU, CUSTOM};
enum TableSetType {total_and_coherent, coherent_and_incoherent};
-#endif
+enum FockState {QQ, QQG, ALL};
+
+#endif
Index: trunk/src/InclusiveDiffractiveCrossSections.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 532)
+++ trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 533)
@@ -1,669 +1,687 @@
//==============================================================================
// InclusiveDiffractiveCrossSection.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 "Math/SpecFunc.h"
#include "TFile.h"
#include <iostream>
#include <cmath>
#include <algorithm>
+#include <algorithm>
#include "Constants.h"
#include "Nucleus.h"
#include "DipoleModel.h"
#include "AlphaStrong.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/Functor.h"
#include "TMath.h"
#include "WaveOverlap.h"
#include "Kinematics.h"
//#include "TableGeneratorSettings.h"
#include "Settings.h"
#include "Enumerations.h"
#include "TF1.h"
#include "TF1.h"
#include "TH1F.h"
#include "cuba.h"
#include "InclusiveDiffractiveCrossSections.h"
#include "DglapEvolution.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include "DipoleModelParameters.h"
#define PR(x) cout << #x << " = " << (x) << endl;
InclusiveDiffractiveCrossSections::InclusiveDiffractiveCrossSections()
{
mSettings = EventGeneratorSettings::instance();
mRandom = mSettings->randomGenerator();
mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
mPhotonFlux.setS(mS);
mCheckKinematics = true;
mCrossSectionRatioLT = 0;
mQ2=0;
// TableGeneratorSettings* settings = TableGeneratorSettings::instance();
if(mDipoleModel) delete mDipoleModel;
+ mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
DipoleModelType model=mSettings->dipoleModelType();
if (model==bSat) {
mDipoleModel = new DipoleModel_bSat(mSettings);
}
else if(model==bNonSat){
mDipoleModel = new DipoleModel_bNonSat(mSettings);
}
else {
cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
exit(1);
}
mA=mSettings->A();
//Set up the integrals
//Amplitude0:
if(Amp0) delete Amp0;
Amp0=new TF1("Amp0", this, &InclusiveDiffractiveCrossSections::uiAmplitude_0, 0, 30., 5);
if(WFAmp0) delete WFAmp0;
WFAmp0=new ROOT::Math::WrappedTF1(*Amp0);
GIAmp0.SetFunction(*WFAmp0);
- GIAmp0.SetRelTolerance(1e-1);
+ GIAmp0.SetRelTolerance(1e-2);
//Amplitude1:
if(Amp1) delete Amp1;
Amp1=new TF1("Amp1", this, &InclusiveDiffractiveCrossSections::uiAmplitude_1, 0, 30. , 5);
if(WFAmp1) delete WFAmp1;
WFAmp1=new ROOT::Math::WrappedTF1(*Amp1);
GIAmp1.SetFunction(*WFAmp1);
- GIAmp1.SetRelTolerance(1e-1);
+ GIAmp1.SetRelTolerance(1e-2);
//AmplitudeQQG:
if(Ampqqg) delete Ampqqg;
Ampqqg=new TF1("Ampqqg", this, &InclusiveDiffractiveCrossSections::uiAmplitude_qqg, 0, 30., 4);
if(WFAmpqqg) delete WFAmpqqg;
WFAmpqqg=new ROOT::Math::WrappedTF1(*Ampqqg);
GIAmpqqg.SetFunction(*WFAmpqqg);
- GIAmpqqg.SetRelTolerance(1e-1);
+ GIAmpqqg.SetRelTolerance(1e-2);
}
InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){
if(WFAmp0) delete WFAmp0;
if(Amp0) delete Amp0;
if(WFAmp1) delete WFAmp1;
if(Amp1) delete Amp1;
if(Ampqqg) delete Ampqqg;
}
void InclusiveDiffractiveCrossSections::setCheckKinematics(bool val) {mCheckKinematics = val;}
void InclusiveDiffractiveCrossSections::setFockState(FockState val){
mFockState=val;
}
FockState InclusiveDiffractiveCrossSections::getFockState(){
return mFockState;
}
double InclusiveDiffractiveCrossSections::operator()(double MX2, double Q2, double W2, double z){
return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
}
double InclusiveDiffractiveCrossSections::operator()(const double* array){
double result=0;
- if(mFockState==QQ or mFockState==ALL)
+ if(mFockState==QQ or mFockState==ALL){
result+=dsigdbetadQ2dW2dz_total_checked(array[0], array[1], array[2], array[3]);
- else if(mFockState==QQG or mFockState==ALL)
+ }
+ else if(mFockState==QQG or
+ mFockState==ALL){
result+=dsigdbetadQ2dW2dz_total_qqg_checked(array[0], array[1], array[2], array[3]);
+ }
else{
cout<<"InclusiveDiffractiveCrossSections::operator(): Fockstate is invalid, stopping"<<endl;
exit(0);
}
return result;
}
double InclusiveDiffractiveCrossSections::unuranPDF(const double* array){
//
// array is: beta, log(Q2), W2, z
//
double result = 0;
double Q2 = exp(array[1]);
result = dsigdbetadQ2dW2dz_total_checked(array[0], Q2, array[2], array[3]);
result *= Q2; // Jacobian
return log(result);
}
double InclusiveDiffractiveCrossSections::dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z){
double beta=Q2/(Q2+MX2);
double jacobian=beta*beta/Q2;
double result= dsigdbetadQ2dW2dz_total_checked(beta, Q2, W2, z);
return jacobian*result;
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z){
double result = 0;
//
// Check if in kinematically allowed region
double MX2=Kinematics::MX2(Q2, beta, 0);
if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
<< " is outside of kinematically allowed region. Return 0." << endl;
return 0;
}
double xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double csT = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, transverse);
double csL = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, longitudinal);
result = csT + csL;
mCrossSectionRatioLT=csL/csT;
//
// Polarization
//
if (mRandom->Uniform(result) <= csT)
mPolarization = transverse;
else
mPolarization = longitudinal;
//
// Quark Species
//
//
if(mPolarization == transverse){
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(csT-cs_rejected);
double csTi=mDsigdbetadQ2dWdz_T_total[i];
if(R <= csTi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_T_total[i];
}
if(!isChosen) mQuarkSpecies=4;
}
else{
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(csL-cs_rejected);
double csLi=mDsigdbetadQ2dWdz_L_total[i];
if(R <= csLi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_L_total[i];
}
if(!isChosen) mQuarkSpecies=4;
}
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): " << result;
cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z;
cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal");
cout << ')' << endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = NaN at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = inf at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (result < 0) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, negative cross-section at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
result = 0;
}
mTotalCS=result;
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
if(pol==transverse){
result = dsigmadbetadz_T(beta, Q2, W2, z);
}
else if(pol==longitudinal){
result = dsigmadbetadz_L(beta, Q2, W2, z);
}
return result; //nb
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
for(int i=0; i<5; i++){
setQuarkIndex(i);
double tmpresult = dsigdbetadz_total(beta, Q2, W2, z, pol); //nb
if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,pol); //nb/GeV4
if(pol == transverse)
mDsigdbetadQ2dWdz_T_total[i]=tmpresult;
if(pol == longitudinal)
mDsigdbetadQ2dWdz_L_total[i]=tmpresult;
result+=tmpresult;
}
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::crossSectionOfLastCall(){
return mTotalCS;
}
GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){
return mPolarization;
}
unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){
return mQuarkSpecies;
}
double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){
- return Settings::quarkMass(mQuarkSpecies);
+// return Settings::quarkMass(mQuarkSpecies); //#TT tmp
+ tt_quarkMass[mQuarkSpecies];
}
void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {
double MX2 = Q2*(1-beta)/beta;
double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
- double mf = Settings::quarkMass(mQuarkIndex);
+// double mf = Settings::quarkMass(mQuarkIndex);
+ double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
double sqrtArg = 1.-4.*mf*mf/MX2;
if (sqrtArg<0){
if (mSettings->verboseLevel() > 2)
cout<<"There is no phase-space for z, return 0."<<endl;
return 0;
}
double z0 = (1.-sqrt(sqrtArg))/2.;
if (z<z0 or z>1-z0){
if (mSettings->verboseLevel() > 2)
cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
return 0;
}
double pt2 = z*(1-z)*Q2-mf*mf;
if (pt2<0){
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
return 0;
}
if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark production, return 0."<<endl;
return 0;
}
double epsilon2=z*(1-z)*Q2+mf*mf;
double ef=quarkCharge[mQuarkIndex];
double Phi0=Phi_0(beta, xpom, z, Q2, mf);
double Phi1=Phi_1(beta, xpom, z, Q2, mf);
double prefactor=Nc*Q2*alpha_em/(4.*M_PI*beta*beta)*ef*ef*z*(1-z); //GeV2
prefactor/=2.; //expanded z-range
double term1=epsilon2*(z*z+(1-z)*(1-z))*Phi1; //GeV2*fm6
double term2=mf*mf*Phi0; //GeV2*fm6
double result=prefactor*(term1+term2); //GeV4*fm6
result/=hbarc2*hbarc2; //fm2
result *= 1e7; //nb
return result;//nb
}
double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double beta, double Q2, double W2, double z){
double MX2 = Q2*(1-beta)/beta;
double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
- double mf = Settings::quarkMass(mQuarkIndex);
+// double mf = Settings::quarkMass(mQuarkIndex);
+ double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
double pt2 = z*(1-z)*Q2-mf*mf;
if (pt2<0){
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
return 0;
}
double sqrtArg=1.-4.*mf*mf/MX2;
if (sqrtArg<0){
if (mSettings->verboseLevel() > 2)
cout<<"There is no phase-space for z, return 0."<<endl;
return 0;
}
double z0=(1.-sqrt(sqrtArg))/2.;
if (z<z0 or z>1-z0){
if (mSettings->verboseLevel() > 2)
cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
return 0;
}
if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark production, return 0."<<endl;
return 0;
}
double ef=quarkCharge[mQuarkIndex];
double Phi0=Phi_0(beta, xpom, z, Q2, mf); //fm6
double term=Nc*Q2*Q2*alpha_em/(M_PI*beta*beta); //GeV4
term*=ef*ef*z*z*z*(1-z)*(1-z)*(1-z);
term*=Phi0; //fm6*GeV4
term/=2.; //expanded z-range
double result=term;
result/=hbarc2*hbarc2; //fm2
result *= 1e7; //nb
return result; //nb
}
double InclusiveDiffractiveCrossSections::Phi_0(double beta, double xpom, double z, double Q2, double mf){
TF1* Phi0=new TF1("Phi0", this, &InclusiveDiffractiveCrossSections::uiPhi_0, 0, 30., 0);
ROOT::Math::WrappedTF1* WFPhi0=new ROOT::Math::WrappedTF1(*Phi0);
ROOT::Math::GaussIntegrator GIPhi0;
GIPhi0.SetFunction(*WFPhi0);
- GIPhi0.SetRelTolerance(1e-1);
+ GIPhi0.SetRelTolerance(1e-2);
Amp0->SetParameter(0, beta);
Amp0->SetParameter(1, xpom);
Amp0->SetParameter(2, z);
Amp0->SetParameter(3, Q2);
Amp0->SetParameter(5, mf);
double blow=0;
double result=GIPhi0.IntegralUp(blow); //fm6
delete Phi0;
delete WFPhi0;
return result;
}
double InclusiveDiffractiveCrossSections::Phi_1(double beta, double xpom, double z, double Q2, double mf){
TF1* Phi1=new TF1("Phi1", this, &InclusiveDiffractiveCrossSections::uiPhi_1, 0, 30., 0);
ROOT::Math::WrappedTF1* WFPhi1=new ROOT::Math::WrappedTF1(*Phi1);
ROOT::Math::GaussIntegrator GIPhi1;
GIPhi1.SetFunction(*WFPhi1);
- GIPhi1.SetRelTolerance(1e-1);
+ GIPhi1.SetRelTolerance(1e-2);
Amp1->SetParameter(0, beta);
Amp1->SetParameter(1, xpom);
Amp1->SetParameter(2, z);
Amp1->SetParameter(3, Q2);
Amp1->SetParameter(5, mf);
double blow=0;
double result=GIPhi1.IntegralUp(blow); //fm6
delete Phi1;
delete WFPhi1;
return result;
}
double InclusiveDiffractiveCrossSections::uiPhi_0(double *var, double *par){
double b=*var;
double amplitude=Amplitude_0(b); //fm2
double result=2*M_PI*b*amplitude*amplitude; //fm5
return result;
}
double InclusiveDiffractiveCrossSections::uiPhi_1(double *var, double *par){
double b=*var; //fm
double amplitude=Amplitude_1(b); //fm2
double result=2*M_PI*b*amplitude*amplitude; //fm5
return result;
}
double InclusiveDiffractiveCrossSections::Amplitude_0(double b){
//Gauss Integrator defined in header file and set in constructor
Amp0->SetParameter(4, b);
double rlow=1e-3;
double result=GIAmp0.IntegralUp(rlow); //fm2
return result;
}
double InclusiveDiffractiveCrossSections::Amplitude_1(double b){
//Gauss Integrator defined in header file and set in constructor
Amp1->SetParameter(4, b);
double rlow=1e-3;
double result=GIAmp1.IntegralUp(rlow); //fm2
return result;
}
double InclusiveDiffractiveCrossSections::uiAmplitude_0(double *var, double *par){
double r=*var;
double beta=par[0];
double xpom=par[1];
double z=par[2];
double Q2=par[3];
double b=par[4];
double mf=par[5];
double MX2=Q2*(1-beta)/beta;
double kappa2=z*(1-z)*MX2-mf*mf;
double kappa=sqrt(kappa2);
double epsilon2=z*(1-z)*Q2+mf*mf;
double epsilon=sqrt(epsilon2); //GeV
double besselK0=TMath::BesselK0(epsilon*r/hbarc);
double besselJ0=TMath::BesselJ0(kappa*r/hbarc);
double dsigmadb2=0;
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
return r*besselK0*besselJ0*dsigmadb2; //fm
}
double InclusiveDiffractiveCrossSections::uiAmplitude_1(double *var, double *par){
double r=*var; //fm
double beta=par[0];
double xpom=par[1];
double z=par[2];
double Q2=par[3]; //GeV2
double b=par[4];
double mf=par[5];
double MX2=Q2*(1-beta)/beta; //GeV2
double kappa2=z*(1-z)*MX2-mf*mf; //GeV2
double kappa=sqrt(kappa2); //GeV
double epsilon2=z*(1-z)*Q2+mf*mf;
double epsilon=sqrt(epsilon2); //GeV
double besselK1=TMath::BesselK1(epsilon*r/hbarc);
double besselJ1=TMath::BesselJ1(kappa*r/hbarc);
double dsigmadb2=0; //GeV0
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
return r*besselK1*besselJ1*dsigmadb2; //fm
}
DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel(){
return mDipoleModel;
}
double InclusiveDiffractiveCrossSections::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
//===================================================
// QQG Calculations
//===================================================
double InclusiveDiffractiveCrossSections::unuranPDF_qqg(const double* array){
//
// array is: beta, log(Q2), W2, z
//
double result = 0;
double Q2 = exp(array[1]);
- result = dsigdbetadQ2dW2dz_total_checked(array[0], Q2, array[2], array[3]);
- result *= Q2; // Jacobian
+
+ double z=array[3];
+ result = dsigdbetadQ2dW2dz_total_qqg_checked(array[0], Q2, array[2], z);
+ result *= Q2; // Jacobian log(Q2)->Q2
return log(result);
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z){
//
// Check if in kinematically allowed region
//
double MX2=Kinematics::MX2(Q2, beta, 0);
+ if(mCheckKinematics && z<beta) {
+ if (mSettings->verboseLevel() > 2)
+ cout<<"z<beta, beta="<<beta<<", z="<<z<<endl;
+ return 0; // For QQG beta < z < 1
+ }
if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
- cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
<< " is outside of kinematically allowed region. Return 0." << endl;
return 0;
}
double xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double result = dsigdbetadQ2dW2dz_qqg(beta, Q2, W2, z);
-
//
// Quark Species
//
//
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(result-cs_rejected);
- double csTi=mDsigdbetadQ2dWdz_T_total[i];
+ double csTi=mDsigdbetadQ2dWdz_T_qqg[i];
if(R <= csTi){
mQuarkSpecies=i;
isChosen=true;
break;
}
- cs_rejected+=mDsigdbetadQ2dWdz_T_total[i];
+ cs_rejected+=mDsigdbetadQ2dWdz_T_qqg[i];
}
if(!isChosen) mQuarkSpecies=4;
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): " << result;
cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z<< endl;
}
//
// Check validity of return value
//
if (std::isnan(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = NaN at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (std::isinf(result)) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = inf at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
result = 0;
}
if (result < 0) {
cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, negative cross-section at" << endl;
cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
result = 0;
}
mTotalCS_qqg=result;
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z) {
double result = 0;
for(int i=0; i<5; i++){
setQuarkIndex(i);
+// double mf = Settings::quarkMass(mQuarkIndex);
+ double mf = tt_quarkMass[i]; //#TT tmp
+ double Mqq2=(z/beta-1)*Q2;
+ if(Mqq2<4*mf*mf) continue; //Not enough phase-space for the q and qbar
double tmpresult = dsigmadbetadz_QQG(beta, Q2, W2, z); //nb
if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,transverse); //nb/GeV4
mDsigdbetadQ2dWdz_T_qqg[i]=tmpresult;
result+=tmpresult;
}
return result; //nb/GeV4
}
double InclusiveDiffractiveCrossSections::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
-
double MX2 = Q2*(1-beta)/beta;
double xpom = Kinematics::xpomeron(0., Q2, W2, sqrt(MX2));
double Phi=Phi_qqg(xpom, ztilde, Q2);
double alpha_S=0.14;
double ef=quarkCharge[mQuarkIndex];
double betafactor=(1-beta/ztilde)*
(1-beta/ztilde)+beta/ztilde*beta/ztilde;
double result=alpha_S*alpha_em/(2*M_PI*M_PI*Q2)*betafactor*ef*ef*Phi; //GeV-2
result*=hbarc2; //fm2
result *= 1e7; //nb
return result; //nb
}
double InclusiveDiffractiveCrossSections::Phi_qqg(double xpom, double z, double Q2){
Ampqqg->SetParameter(0, xpom);
Ampqqg->SetParameter(1, z);
mQ2=Q2;
ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiPhi_qqg, 2);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig.SetFunction(wf);
ig.SetAbsTolerance(0.);
- ig.SetRelTolerance(1e-1);
- double lo[2]={0.,0.}; //b, k2
+ ig.SetRelTolerance(1e-2);
+ double lo[2]={0.,1e-4}; //b, k2
double hi[2]={5., Q2}; //b, k2
double result = ig.Integral(lo, hi)/hbarc; //GeV0
return result; //GeV0
}
double InclusiveDiffractiveCrossSections::uiPhi_qqg(const double* var) {
double b=var[0];
double k2=var[1];
double Q2=mQ2;
Ampqqg->SetParameter(2, sqrt(k2));
Ampqqg->SetParameter(3, b);
double rlow=1e-3;
double amp_qqg=GIAmpqqg.IntegralUp(rlow); //fm2
amp_qqg/=hbarc2; //GeV-2
double result = (2*M_PI)*b/hbarc*k2*k2*log(Q2/k2)*amp_qqg*amp_qqg; //GeV-1
return result; //GeV-1
}
double InclusiveDiffractiveCrossSections::uiAmplitude_qqg(double *var, double *par){
double r=*var; //fm
+ 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 k=par[2]; //GeV
double b=par[3];
double besselK2=ROOT::Math::cyl_bessel_k(2, sqrt(z)*k*r/hbarc);
double besselJ2=ROOT::Math::cyl_bessel_j(2, sqrt(1-z)*k*r/hbarc);
double dsigmadb2=0; //GeV0
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
double term=1-0.5*dsigmadb2;
double dsigmadb2tilde=2*(1-term*term);
double result=r*besselK2*besselJ2*dsigmadb2tilde; //fm
return result;
}
Index: trunk/src/Constants.h
===================================================================
--- trunk/src/Constants.h (revision 532)
+++ trunk/src/Constants.h (revision 533)
@@ -1,49 +1,52 @@
//==============================================================================
// Constants.h
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date$
// $Author$
//==============================================================================
#ifndef Constants_h
#define Constants_h
//
// General constants
//
const double electronMass = 0.510998902E-3; // GeV
const double electronMass2 = electronMass*electronMass; // GeV^2
const double protonMass = 0.9382700; // GeV
const double protonMass2 = protonMass*protonMass; // GeV^2
const double alpha_em = 1/137.036;
const double hbarc = 0.197327; // GeV*fm
-const double hbarc2 = hbarc*hbarc; // GeV^2*fm^2
+//#TT a temporary hack:
+const double tt_quarkMass[6]={0.004008411318047, 0.004008411318047, 0.004008411318047, 1.428034604225, 4.2, 175.};//#TT tmp
+
+const double hbarc2 = hbarc*hbarc; // GeV^2*fm^2
//
// Constants used in dipole model
//
const double Nc = 3.;
const double quarkCharge[6] = {2./3., -1./3., -1./3., 2./3., -1./3., 2./3.}; // u, d, s, c, b, t
//
// Constants used in integartion and other numerical operations
//
const double upperIntegrationLimit = 2.5; // factor: nuclear radius -> integration limits (b, r)
#endif
Index: trunk/src/SartreInclusiveDiffraction.cpp
===================================================================
--- trunk/src/SartreInclusiveDiffraction.cpp (revision 532)
+++ trunk/src/SartreInclusiveDiffraction.cpp (revision 533)
@@ -1,764 +1,852 @@
//==============================================================================
// 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 "Math/Functor.h"
#include "TUnuranMultiContDist.h"
#include <string>
#include <limits>
#include <cmath>
#include <iomanip>
#include "TH2D.h"
#include "TFile.h"
#include "TF3.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;
+// 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;
//
// 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;
mCrossSection = new InclusiveDiffractiveCrossSections;
-
//
- // Here we need to put inQ2 > mu02.
+ // 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 = Settings::quarkMass(0);
+ double mf=tt_quarkMass[0];
//Now for the diffractive mass:
double kineMX2min = 4*mf*mf;//max(mu02, 4*mf*mf);
// 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-4*mf*mf/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 and W2,
- // not beta (for now) and z.
+ // 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-4*mf*mf/kineMX2max))/2.;
// kineZmax=1-kineZmin;
mLowerLimit[3]=(1-sqrt(1-4*mf*mf/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);
}
//
// 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 kinemtaically easier to use MX2 instead of beta here.
+ // 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<1; 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+4*mf*mf/Q2);
+ mLowerLimit_qqg[3]=mLowerLimit[0] * (1+4*mf*mf/exp(mUpperLimit[1]));
+
+ 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+4*mf*mf/exp(theMode_qqg[1])); //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<1; 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=totalCrossSection();
+ cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
+ mCrossSection->setFockState(QQG);
+ mTotalCrossSection=0;
+ mTotalCrossSection_qqg=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;
}
-
- double xrandom[4];
-
+ //
+ // 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
//
- mUnuran->SampleMulti(xrandom);
+ 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);
- mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
+ if(isQQ)
+ mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
+ else
+ mCurrentEvent->polarization = transverse;
// mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
mCurrentEvent->quarkSpecies = mCrossSection->quarkSpeciesOfLastCall();
- mCurrentEvent->crossSectionRatioLT=mCrossSection->crossSectionRatioLTOfLastCall();
+ if(isQQ)
+ mCurrentEvent->crossSectionRatioLT=mCrossSection->crossSectionRatioLTOfLastCall();
+ else
+ mCurrentEvent->crossSectionRatioLT=0;
+
+ //
+ // The following needs to be modified for QQG
+ //
double dsigmadb2=0;
double z=xrandom[3];
double eps2=z*(1-z)*(mCurrentEvent->Q2+MX2);
double r_average=1./sqrt(eps2);
if (mA == 1) {
dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(r_average*hbarc, 0, mCurrentEvent->xpom);
}
else {
dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(r_average*hbarc, 0, mCurrentEvent->xpom);
- if(mCurrentEvent->xpom<1e-2){
- cout<<"==============="<<endl;
- PR(dsigmadb2);
- PR(mCrossSection->dipoleModel()->dsigmadb2ep(r_average*hbarc, 0, mCurrentEvent->xpom));
- PR(mCurrentEvent->xpom);
- PR(r_average*hbarc);
- cout<<"==============="<<endl;
- }
}
mCurrentEvent->Omega=-2.*log(1-dsigmadb2/2.);
+
+
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 (mSettings->UPC()) {
- // also fills some of the undefined event variables
- ok = mFinalStateGenerator.generate(mA, mCurrentEvent);
+ if (isQQ) {
+ ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQ);
}
else {
- ok = mFinalStateGenerator.generate(mA, mCurrentEvent);
+ 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) or log(xpom).
+ // are in log(Q2)
//
- xmin[1] = exp(xmin[1]); // log Q2 limit -> Q2 limit or equivalent xpom for UPC
- xmax[1] = exp(xmax[1]); // log Q2 limit -> Q2 limit
-
+ 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 = 1e-1;//5e-4;
- unsigned int numberofcalls = 10000;//1000000;
- ROOT::Math::Functor wfL((*mCrossSection), 4);
- ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE,
- 0, precision, numberofcalls);
-
+ unsigned int numberofcalls = 1;//1000000;
+ 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;
}
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 532)
+++ trunk/src/InclusiveFinalStateGenerator.cpp (revision 533)
@@ -1,419 +1,773 @@
//==============================================================================
// 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::generate(int A, Event *event)
+bool InclusiveFinalStateGenerator::generate(int A, Event *event, FockState fockstate)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
+// cout<<"#TT Generating final state"<<endl;
//
// 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: 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();
double pt_squared = E*E-pz*pz-mMY2;
if (pt_squared < 0) {
- if (settings->verbose() && settings->verboseLevel() > 2) {
+ if (settings->verbose() && settings->verboseLevel() > 1) {
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
//
- double mf = Settings::quarkMass(event->quarkSpecies);
+ TLorentzVector theQuark;
+ TLorentzVector theAntiQuark;
+ TLorentzVector theGluon; //for QQG
double z=event->z;
- double pt2=z*(1-z)*Q2-mf*mf;
- phi = rndm->Uniform(twopi);
- double pxq=sqrt(pt2)*sin(phi);
- double pyq=sqrt(pt2)*cos(phi);
- double pxqbar=-pxq;
- double pyqbar=-pyq;
- if (4*(mf*mf+pt2) > MX2) {
- if (settings->verboseLevel() > 2) {
- cout << "InclusiveFinalStateGenerator::generate(): 4*(mf*mf+pt2) > MX2" << endl;
- cout << " Abort event generation" << endl;
+ double pxq=0., pyq=0., pxqbar=0., pyqbar=0.;
+// double mf = Settings::quarkMass(event->quarkSpecies);
+ double mf = tt_quarkMass[event->quarkSpecies]; //#TT tmp
+ if(fockstate==QQ){
+ double pt2=z*(1-z)*Q2-mf*mf;
+ phi = rndm->Uniform(twopi);
+ pxq=sqrt(pt2)*sin(phi)/2.;
+ pyq=sqrt(pt2)*cos(phi)/2.;
+ pxqbar=-pxq;
+ pyqbar=-pyq;
+ if (4*(mf*mf+pt2) > MX2) {
+ if (settings->verboseLevel() > 2) {
+ cout << "InclusiveFinalStateGenerator::generate(): 4*(mf*mf+pt2) > MX2" << endl;
+ cout << " Abort event generation" << endl;
+ }
+ return false;
+ }
+
+ //
+ // Calculate the quark momenta in the X-particle's restframe
+ //
+ TVector3 boost=TVector3(theXparticle.Px()/theXparticle.E(), theXparticle.Py()/theXparticle.E(), theXparticle.Pz()/theXparticle.E());
+ int sign=1.;
+ if(z<0.5) sign=-1.;
+ theQuark=TLorentzVector(pxq, pyq, sign*sqrt(MX2/4-mf*mf-pt2), sqrt(MX2)/2);
+ theAntiQuark=TLorentzVector(pxqbar, pyqbar, -sign*sqrt(MX2/4-mf*mf-pt2), sqrt(MX2)/2);
+ theQuark.Boost(boost);
+ theAntiQuark.Boost(boost);
+ }
+ else if(fockstate==QQG){
+ //The knowns:
+ double Mqq2=(z/beta-1)*Q2;
+ if(Mqq2<4*mf*mf){
+ if (settings->verboseLevel() > 1) {
+ cout << "InclusiveFinalStateGenerator::generate(): 4*mf*mf > 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.
+ //
+ double PpomPlus=thePomeron.E()+thePomeron.Pz();
+ double PgPlus=(1-z)*PpomPlus; //gluon
+ double PqPlus=beta*PpomPlus; //quark
+ double PqbarPlus=(z-beta)*PpomPlus; //antiquark
+
+ //
+ // 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
+ //
+ // First solve for pt4=ptg...
+ double pttot=(theVirtualPhoton+thePomeron).Pt();
+ double A = (PqPlus+PqbarPlus)/PgPlus +PgPlus/(PqPlus+PqbarPlus)+2;
+ double B = -2*pttot*PgPlus/(PqPlus+PqbarPlus)- 2*pttot;
+ double C = Mqq2+PgPlus/(PqPlus+PqbarPlus)*(Mqq2+pttot*pttot)-MX2;
+ //Generate angle in the kinematically allowed region:
+ double phi_totg=0;
+ if(A*C<=0) //There is always a solution
+ phi_totg = rndm->Uniform(twopi); //Angle between transverse total-g
+ else{
+ double cosphi_min=sqrt(4*A*C/B/B);
+ if(cosphi_min>1){ //no solution
+ if (settings->verboseLevel() > 1) {
+ cout << "InclusiveFinalStateGenerator::generate(): Not enough phasespace for angle X-g" << endl;
+ cout << " Abort event generation" << endl;
+ }
+ return false;
+ }
+ double phi_max=acos(cosphi_min);
+ phi_totg = rndm->Uniform(4*phi_max);
+ //find the quadrant:
+ if(phi_totg>4*phi_max)
+ phi_totg+=3*M_PI/2;
+ else if(phi_totg>3*phi_max)
+ phi_totg+=M_PI;
+ else
+ phi_totg+=M_PI/2;
+ }
+ B*=cos(phi_totg);
+ double ptg = B/2/A+sqrt(B*B/4/A/A-C/A);
+ // ...and use this to calculate the rest:
+ double pt32 = pttot*pttot+ptg*ptg-2*pttot*ptg*cos(phi_totg);
+ double pt3 = sqrt(pt32);
+ double PgMinus = ptg*ptg/PgPlus;
+ double PqqMinus = (Mqq2+pt32)/(PqPlus+PqbarPlus);
+ double cosphi34 = (Mqq2+PgMinus*(PqPlus+PqbarPlus)+PgPlus*PqqMinus-MX2)/ (2*pt3*ptg);
+ if(cosphi34*cosphi34>1){
+ if (settings->verboseLevel() > 1) {
+ cout << "InclusiveFinalStateGenerator::generate(): Not enough phasespace for angle between qqbar and gluon" << endl;
+ cout << " Abort event generation" << endl;
+ }
+ return false;
+ }
+ double phi_3g=acos(cosphi34); //angle between gluon and qqbar
+ double phi_tot=acos((theVirtualPhoton+thePomeron).Px() /(theVirtualPhoton+thePomeron).Pt());
+
+ //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,
+ //First create the pseudoparticle 4-vector
+ E=(PqPlus+PqbarPlus+PqqMinus)/2;
+ pz=(PqPlus+PqbarPlus-PqqMinus)/2;
+ px=pt3*cos(phi_tot-phi_totg+phi_3g);
+ py=pt3*sin(phi_tot-phi_totg+phi_3g);
+ TLorentzVector thePseudoGluon(px, py, pz, E);
+
+ double xi=beta/z;
+ double k2=(thePomeron-theGluon).M2();
+ double pt2=xi*(1-xi)*(-k2)-mf*mf;
+ phi = rndm->Uniform(twopi);
+ pxq=sqrt(pt2)*sin(phi)/2;
+ pyq=sqrt(pt2)*cos(phi)/2;
+ pxqbar=-pxq;
+ pyqbar=-pyq;
+ if (4*(mf*mf+pt2) > Mqq2) {
+ if (settings->verboseLevel() > 1) {
+ cout << "InclusiveFinalStateGenerator::generate(): 4*(mf*mf+pt2) > Mqq2" << endl;
+ cout << " Abort event generation" << endl;
+ }
+ return false;
+ }
+
+ //
+ // Calculate the quark momenta in the X-particle's restframe
+ //
+ 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, sign*sqrt(Mqq2/4-mf*mf-pt2), sqrt(Mqq2)/2);
+ theAntiQuark=TLorentzVector(pxqbar, pyqbar, -sign*sqrt(Mqq2/4-mf*mf-pt2), sqrt(Mqq2)/2);
+ theQuark.Boost(boost);
+ theAntiQuark.Boost(boost);
+ }
+ else{
+ cout<<"InclusiveFinalStateGenerator::generate: unknown fockstate: "<<fockstate<<" ending program."<<endl;
+ exit(0);
+ }
+ /*
+ // TODO: Turn this into a function
+ // 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;
+ 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;
}
- return false;
}
-
//
- // Calculate tha quark momenta in the X-particle's restframe
+ // 2. Now generate 2*Twist number of gluons:
+ // TODO: Merge this loop with the interaction loop below
+ 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++){
+ //
+ // 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()+p_x);
+ pom_py.push_back(thePomeron.Py()+p_y);
+ pom_px_tot+=p_x;
+ pom_py_tot+=p_y;
+ }
+ //conserve the total pomeron momentum:
+ pom_px.push_back(thePomeron.Px()-pom_px_tot);
+ pom_py.push_back(thePomeron.Py()-pom_py_tot);
+ //
+ // 3. Distribute the recoils to:
+ // the quark and antiquark (QQ) or gluon and pseudogluon (QQG):
+ // Call these q and qbar (for parton)
+ //
+ 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();
+ }
+ else if(fockstate==QQG){
+ Eq=theQuark.E()+theAntiQuark.E();
+ Eqbar=theGluon.E();
+ pzq=theQuark.Pz()+theAntiQuark.Pz();
+ pzqbar=theGluon.Pz();
+ pxq=theQuark.Px()+theAntiQuark.Px();
+ pyq=theQuark.Py()+theAntiQuark.Py();
+ pxqbar=theGluon.Px();
+ pyqbar=theGluon.Py();
+ }
+ else{
+ cout<<"InclusiveFinalStateGenerator::generate: unknown fockstate: "<<fockstate<<" ending program."<<endl;
+ exit(9);
+ }
+ for(unsigned int i=0; i<pom_px.size(); i++){
+ //TODO: Merge this loop to where we generate the gluons
+ //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]);
+ double qP=sqrt(pxq*pxq + pyq*pyq);
+ double qbarP=sqrt(pxqbar*pxqbar + pyqbar*pyqbar);
+ double deltaVq=sqrt(1+qP*qP/Eq/Eq - 2*(pom_px[i]*pxq+ pom_py[i]*pyq) /Eq/pomP);
+ double deltaVqbar=sqrt(1+qbarP*qbarP/Eqbar/Eqbar - 2*(pom_px[i]*pxqbar+pom_py[i]*pyqbar) / 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];
+
+ //Scale the transverse momenta in order to preserve MX2,
+ //First find the scaling factor lambda
+ //
+ // From ChatGPT
+ //
+ double mass=0;
+ if(fockstate==QQ)
+ mass=mf;
+ else
+ mass=sqrt((z/beta-1)*Q2);
+ PR(z);
+ PR(beta);
+ PR(mass);
+ PR(mf);
+ PR(fockstate);
+ double lower_bracket=0, upper_bracket=100;
+ double lambda=computeScalingFactor(pxq, pyq, pzq, pxqbar, pyqbar, pzqbar, mass, sqrt(MX2), lower_bracket, upper_bracket);
+ if(lambda=(lower_bracket+upper_bracket)
+ PR(lambda);
+ pxq*=lambda;
+ pyq*=lambda;
+ pxqbar*=lambda;
+ pyqbar*=lambda;
+ //Put on shell:
+ Eq=sqrt(pxq*pxq+pyq*pyq+pzq*pzq+mf*mf);
+ Eqbar=sqrt(pxqbar*pxqbar+pyqbar*pyqbar+pzqbar*pzqbar+mf*mf);
+ PR(MX2);
+ PR(2*mf*mf+2*Eq*Eqbar-2*pxq*pxqbar-2*pyq*pyqbar-2*pzq*pzqbar);
+ exit(1);
+ }
//
- TVector3 boost=TVector3(theXparticle.Px()/theXparticle.E(), theXparticle.Py()/theXparticle.E(), theXparticle.Pz()/theXparticle.E());
- int sign=1.;
- if(z<0.5) sign=-1.;
- TLorentzVector theQuark(pxq, pyq, sign*sqrt(MX2/4-mf*mf-pt2), sqrt(MX2)/2);
- TLorentzVector theAntiQuark(pxqbar, pyqbar, -sign*sqrt(MX2/4-mf*mf-pt2), sqrt(MX2)/2);
- theQuark.Boost(boost);
- theAntiQuark.Boost(boost);
-
+ //Update the 4-vectors:
+ //
+ if(fockstate==QQ){
+ //the Quark:
+ theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
+ //the Anti-Quark
+ theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
+ }
+ else{
+ //the Quark/antiquark:
+ //#TT TODO: These have to be split up so that Mqq2 is preserved
+ theQuark=TLorentzVector(pxq/2, pyq/2, pzq/2, Eq/2);
+ theAntiQuark=TLorentzVector(pxq/2, pyq/2, pzq/2, Eq/2);
+
+ //the gluon:
+ theGluon=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
+ }
+ //
+ // Decorrelations ended
+ //
+ */
//
// Add particles to event record
//
- event->particles.resize(2+5);
+ 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=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
//
mPythiaEvent->reset();
- int quarkID = event->quarkSpecies + 1; // Pythia starts with u = 1, Sartre with u = 0
+ // Pythia starts with u = 1, Sartre with u = 0
+ int quarkID = event->quarkSpecies + 1;
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());
-
+ 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; // generated stuff starts at i=3 in pythia 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
+double InclusiveFinalStateGenerator::computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m, double M) {
+ // Define the function to compute the total energy after scaling
+ auto totalEnergy = [&](double lambda) {
+ double E1_prime = std::sqrt(lambda*lambda*(px1*px1 + py1*py1) + pz1*pz1 + m*m);
+ double E2_prime = std::sqrt(lambda*lambda*(px2*px2 + py2*py2) + pz2*pz2 + m*m);
+ 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 = pz1 + pz2;
+ double energy_sum = totalEnergy(lambda);
+
+ return energy_sum*energy_sum - px_sum*px_sum - py_sum*py_sum - pz_sum*pz_sum - M*M;
+ };
+
+ // Solve for lambda using a numerical method (bisection method for simplicity)
+ double lambda_min = 0., lambda_max = 100.0;
+ 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);
+}
Index: trunk/src/DipoleModel.cpp
===================================================================
--- trunk/src/DipoleModel.cpp (revision 532)
+++ trunk/src/DipoleModel.cpp (revision 533)
@@ -1,369 +1,366 @@
//==============================================================================
// DipoleModel.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date$
// $Author$
//==============================================================================
#include <fstream>
#include <iostream>
#include <sstream>
#include <cmath>
#include "DipoleModel.h"
#include "TableGeneratorSettings.h"
#include "DglapEvolution.h"
#include "Constants.h"
#include "TFile.h"
#include "TVector3.h"
#include "TMath.h"
#include "TH2F.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
DipoleModel::DipoleModel()
{
mConfigurationExists = false;
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = settings->A();
mNucleus.init(A);
mIsInitialized = true;
mParameters = nullptr;
}
DipoleModel::~DipoleModel()
{
delete mParameters;
}
const TableGeneratorNucleus* DipoleModel::nucleus() const { return &mNucleus; }
bool DipoleModel::configurationExists() const { return mConfigurationExists; }
double DipoleModel::bDependence(double) { return 0; }
double DipoleModel::bDependence(double, double) { return 0; }
double DipoleModel::dsigmadb2ep(double, double, double) { return 0;}
//***********bSat:*****************************************************
DipoleModel_bSat::DipoleModel_bSat()
{
mBDependence = 0;
mSigma_ep_LookupTable = 0;
//
// Set the parameters. Note that we enforce here the bSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bSat, settings->dipoleModelParameterSet());
-
-}
+}
DipoleModel_bSat::DipoleModel_bSat(Settings* settings)
{
//
// Set the parameters. Note that we enforce here the bSat model
// independent of what the settings say.
//
mParameters = new DipoleModelParameters(settings);
}
DipoleModel_bSat::~DipoleModel_bSat()
{
delete mSigma_ep_LookupTable;
delete mBDependence;
delete mParameters;
}
DipoleModel_bSat& DipoleModel_bSat::operator=(const DipoleModel_bSat& dp)
{
if (this != &dp) {
delete mBDependence;
delete mSigma_ep_LookupTable;
DipoleModel::operator=(dp);
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
mSigma_ep_LookupTable = new TH1F(*(dp.mSigma_ep_LookupTable));
mSigma_ep_LookupTable->SetDirectory(0);
}
return *this;
}
DipoleModel_bSat::DipoleModel_bSat(const DipoleModel_bSat& dp) : DipoleModel(dp)
{
if (mBDependence) delete mBDependence;
mBDependence = new TH2F(*(dp.mBDependence));
mBDependence->SetDirectory(0);
}
void DipoleModel_bSat::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bSat::createConfiguration(): DipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
unsigned int A = mNucleus.A();
string path=settings->bSatLookupPath();
ostringstream filename;
filename.str("");
filename << path << "/bSat_bDependence_A" << A <<".root";
ifstream ifs(filename.str().c_str());
if (!ifs) {
cout << "DipoleModel_bSat::createConfiguration(): File does not exist: " << filename.str().c_str() << endl;
cout << "Stopping." << endl;
exit(1);
}
TFile* lufile= new TFile(filename.str().c_str());
ostringstream histoName;
histoName.str( "" );
histoName << "Configuration_" << iConfiguration;
lufile->GetObject( histoName.str().c_str(), mBDependence );
mBDependence->SetDirectory(0);
lufile->Close();
mConfigurationExists=true;
}
double DipoleModel_bSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
double bDep=bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
return result;
}
double DipoleModel_bSat::bDependence(double b, double phi)
{
double result = mBDependence->Interpolate(b, phi);
return result;
}
double DipoleModel_bSat::dsigmadb2ep(double r, double b, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
- if(muQ2>1e7 or muQ2<1) PR(muQ2);
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
double result = 2.*(1. - exp(-omega/2));
-
return result;
}
double DipoleModel_bSat::coherentDsigmadb2(double r, double b, double xprobe) {
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
double sigmap = mSigma_ep_LookupTable->Interpolate(r);
int A=nucleus()->A();
double TA=nucleus()->T(b)/A;
double result = 2 * ( 1 - pow(1 - TA/2.*sigmap, A) );
return result;
}
void DipoleModel_bSat::createSigma_ep_LookupTable(double xprobe)
{
double rbRange=3.*nucleus()->radius();
TF1* dsigmaForIntegration = new TF1("dsigmaForIntegration", this, &DipoleModel_bSat::dsigmadb2epForIntegration, 0., rbRange, 2);
if(mSigma_ep_LookupTable) delete mSigma_ep_LookupTable;
mSigma_ep_LookupTable = new TH1F("", "", 1000, 0, rbRange);
ROOT::Math::WrappedTF1* WFlookup=new ROOT::Math::WrappedTF1(*dsigmaForIntegration);
ROOT::Math::GaussIntegrator GIlookup;
GIlookup.SetFunction(*WFlookup);
GIlookup.SetRelTolerance(1e-8);
GIlookup.SetAbsTolerance(0.);
// dsigmaForIntegration->SetNpx(1000);
for (int iR=1; iR<=1000; iR++) {
double r=mSigma_ep_LookupTable->GetBinCenter(iR);
dsigmaForIntegration->SetParameter(0, r);
dsigmaForIntegration->SetParameter(1, xprobe);
// double result=dsigmaForIntegration->Integral(0, rbRange);
double result = GIlookup.IntegralUp(0.); //GeV-2
mSigma_ep_LookupTable->SetBinContent(iR, result);
}
delete dsigmaForIntegration; //GeV-2
delete WFlookup;
}
double DipoleModel_bSat::dsigmadb2epForIntegration(double *x, double* par)
{
double r=par[0]; //fm
double xprobe=par[1];
double b = *x; //fm
return 2*M_PI*b/hbarc2*dsigmadb2ep(r, b, xprobe); //GeV-2fm-1
}
//***********bNonSat:*************************************************
DipoleModel_bNonSat::DipoleModel_bNonSat()
{
//
// Set the parameters. Note that we need bNonSat to calculate
// the skewedness correction for bSat.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters( settings->dipoleModelType(), settings->dipoleModelParameterSet());
}
DipoleModel_bNonSat::DipoleModel_bNonSat(Settings* settings)
{
//
// Set the parameters. Note that we need bNonSat to calculate
// the skewedness correction for bSat.
//
mParameters = new DipoleModelParameters(settings->dipoleModelType(), settings->dipoleModelParameterSet());
}
DipoleModel_bNonSat::~DipoleModel_bNonSat(){
delete mParameters;
}
double DipoleModel_bNonSat::dsigmadb2ep(double r, double b, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
const double BG = mParameters->BG(); // GeV^-2
double arg = b*b/(2*BG);
arg /= hbarc2;
double bDep= 1/(2*M_PI*BG) * exp(-arg);
double Mu02 = mParameters->mu02(); // GeV^2
double muQ2 = mParameters->C()/(r*r/hbarc2) + Mu02;
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::dsigmadb2(double r, double b, double phi, double xprobe)
{
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
double bDep=bDependence(b, phi);
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double omega = ((M_PI*M_PI)/Nc)*(r*r/hbarc2)*asxg*bDep;
return omega;
}
double DipoleModel_bNonSat::coherentDsigmadb2(double r, double b, double xprobe){
if (mParameters->dipoleModelParameterSet() == STU) {
double rm = mParameters->rMax();
r = rm*sqrt(log(1+r*r/(rm*rm)));
}
int A=nucleus()->A();
double TA=nucleus()->T(b)/A;
double muQ2 = mParameters->C()/(r*r/hbarc2) + mParameters->mu02();
double asxg = DglapEvolution::instance().alphaSxG(xprobe, muQ2);
double result=A*TA*M_PI*M_PI/Nc*r*r/hbarc2*asxg;
return result;
}
//***********bCGC:*****************************************************
DipoleModel_bCGC::DipoleModel_bCGC()
{
//
// Set the parameters. Note that we enforce here the bNonSat model
// independent of what the settings say.
//
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
mParameters = new DipoleModelParameters(bCGC, settings->dipoleModelParameterSet());
}
void DipoleModel_bCGC::createConfiguration(int iConfiguration)
{
if (!mIsInitialized) {
cout << "DipoleModel_bCGC::createConfigurationDipoleModel class has not been initialized! Stopping." << endl;
exit(1);
}
mNucleus.generate();
mConfigurationExists=true;
}
double DipoleModel_bCGC::dsigmadb2(double r, double b, double phi, double x)
{
double result=1;
for (unsigned int iA=0; iA<mNucleus.A(); iA++) {
double absdeltab=( TVector3(b*cos(phi), b*sin(phi), 0.)
-mNucleus.configuration.at(iA).position() ).Perp();
result*=(1.-0.5*dsigmadb2ep(r, absdeltab, x));
}
return 2.*(1.-result);
}
double DipoleModel_bCGC::dsigmadb2ep(double r, double b, double xprobe)
{
double Y = log(1/xprobe);
double kappa = mParameters->kappa();
double N0 = mParameters->N0();
double x0 = mParameters->x0();
double lambda = mParameters->lambda();
double gammas = mParameters->gammas();
double A = -N0*N0*gammas*gammas/((1-N0)*(1-N0)*log(1-N0));
double B = 0.5*pow(1-N0,-(1-N0)/(N0*gammas));
double Qs = pow(x0/xprobe,lambda/2)*sqrt(DipoleModel_bCGC::bDependence(b));
double rQs = r*Qs/hbarc;
double result=0;
if (rQs <= 2)
result = 2*N0*pow(0.5*rQs, 2*(gammas+(1/(kappa*lambda*Y))*log(2/rQs)));
else
result = 2*(1 - exp(-A*log(B*rQs)*log(B*rQs)));
return result;
}
double DipoleModel_bCGC::bDependence(double b)
{
double gammas = mParameters->gammas();
double Bcgc = mParameters->Bcgc();
return exp(-0.5*b*b/Bcgc/gammas/hbarc2);
}
Index: trunk/src/SartreInclusiveDiffraction.h
===================================================================
--- trunk/src/SartreInclusiveDiffraction.h (revision 532)
+++ trunk/src/SartreInclusiveDiffraction.h (revision 533)
@@ -1,110 +1,120 @@
//==============================================================================
// SartreInclusiveDiffraction.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 SartreInclusiveDiffraction_h
#define SartreInclusiveDiffraction_h
#include "Event.h"
#include "EventGeneratorSettings.h"
#include "InclusiveFinalStateGenerator.h"
#include "InclusiveDiffractiveCrossSections.h"
#include "FrangibleNucleus.h"
#include "Enumerations.h"
#include "TLorentzVector.h"
#include "Math/Functor.h"
#include "TUnuran.h"
#include <ctime>
#include <iostream>
#include <vector>
using namespace std;
class TUnuranMultiContDist;
class SartreInclusiveDiffraction {
public:
SartreInclusiveDiffraction();
virtual ~SartreInclusiveDiffraction();
virtual bool init(const char* = 0);
virtual bool init(const string&);
virtual Event* generateEvent();
-
+
virtual double totalCrossSection(); // in kinematic limits used for generation
virtual double totalCrossSection(double lower[4], double upper[4]); // beta, Q2, W, z
virtual double integrationVolume();
EventGeneratorSettings* runSettings();
const FrangibleNucleus* nucleus() const;
void listStatus(ostream& os=cout) const;
time_t runTime() const;
vector<pair<double,double> > kinematicLimits(); // t, Q2, W
private:
virtual double calculateTotalCrossSection(double lower[4], double upper[4]); //beta, Q2, W, z
private:
+ TRandom3* mRandom;
SartreInclusiveDiffraction(const SartreInclusiveDiffraction&);
SartreInclusiveDiffraction operator=(const SartreInclusiveDiffraction&);
bool mIsInitialized;
time_t mStartTime;
unsigned long mEvents;
unsigned long mTries;
double mTotalCrossSection;
+ double mTotalCrossSection_qq, mTotalCrossSection_qqg;
+
TLorentzVector mElectronBeam;
TLorentzVector mHadronBeam;
double mS;
unsigned int mA;
int mVmID;
DipoleModelType mDipoleModelType;
DipoleModelParameterSet mDipoleModelParameterSet;
Event *mCurrentEvent;
FrangibleNucleus *mNucleus;
FrangibleNucleus *mUpcNucleus;
InclusiveDiffractiveCrossSections *mCrossSection;
EventGeneratorSettings *mSettings;
double mLowerLimit[4]; // beta, Q2, W2, z
double mUpperLimit[4]; // beta, Q2, W2, z
-
+ double mLowerLimit_qqg[4]; // beta, Q2, W2, z
+ double mUpperLimit_qqg[4]; // beta, Q2, W2, z
+
ROOT::Math::Functor *mPDF_Functor;
TUnuran *mUnuran;
TUnuranMultiContDist *mPDF;
-
+
+ ROOT::Math::Functor *mPDF_Functor_qqg;
+ TUnuran *mUnuran_qqg;
+ TUnuranMultiContDist *mPDF_qqg;
+
unsigned long mEventCounter;
unsigned long mTriesCounter;
InclusiveFinalStateGenerator mFinalStateGenerator;
double cross_section_wrapper(double*, double*);
+ double crossSectionsClassWrapper(const double*);
};
ostream& operator<<(ostream& os, const TLorentzVector&);
#endif
Index: trunk/src/InclusiveFinalStateGenerator.h
===================================================================
--- trunk/src/InclusiveFinalStateGenerator.h (revision 532)
+++ trunk/src/InclusiveFinalStateGenerator.h (revision 533)
@@ -1,44 +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);
-
+
+ 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, double M);
};
-#endif
+#endif
Index: trunk/src/InclusiveDiffractiveCrossSections.h
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.h (revision 532)
+++ trunk/src/InclusiveDiffractiveCrossSections.h (revision 533)
@@ -1,120 +1,126 @@
//==============================================================================
// InclusiveDiffractiveCrossSection.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 InclusiveDiffractiveCrossSections_h
#define InclusiveDiffractiveCrossSections_h
#include "AlphaStrong.h"
#include "TH1.h"
#include "TMath.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include "DipoleModel.h"
#include "PhotonFlux.h"
#include "Enumerations.h"
#include "Constants.h"
class DipoleModel;
class DipoleModelParameters;
class InclusiveDiffractiveCrossSections {
public:
InclusiveDiffractiveCrossSections();
~InclusiveDiffractiveCrossSections();
double operator()(double beta, double Q2, double W2, double z);
double operator()(const double* array);
double unuranPDF(const double*);
-
+ double unuranPDF_qqg(const double*);
+
GammaPolarization polarizationOfLastCall() ;
double crossSectionOfLastCall();
unsigned int quarkSpeciesOfLastCall();
double quarkMassOfLastCall();
double crossSectionRatioLTOfLastCall() const;
void setCheckKinematics(bool);
-
+ void setFockState(FockState);
+ FockState getFockState();
+
double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;
double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
// InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections&);
// InclusiveDiffractiveCrossSections& operator=(const InclusiveDiffractiveCrossSections&);
double dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z);
+ double dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z);
DipoleModel* dipoleModel();
double dsigmadbetadz_T(double, double, double, double) ;
double dsigmadbetadz_L(double, double, double, double) ;
double dsigmadbetadz_QQG(double, double, double, double) ;
double mDsigdbetadQ2dWdz_L_total[5];
double mDsigdbetadQ2dWdz_T_total[5];
double mDsigdbetadQ2dWdz_T_qqg[5];
void setQuarkIndex(unsigned int);
-
+
protected:
double dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z);
double dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol);
private:
TRandom3* mRandom;
GammaPolarization mPolarization;
PhotonFlux mPhotonFlux;
EventGeneratorSettings* mSettings;
unsigned int mQuarkSpecies;
- double mTotalCS;
-
+ double mTotalCS, mTotalCS_qqg;
+ FockState mFockState;
+
double mS, mQ2;
bool mCheckKinematics;
double Phi_0(double, double, double, double, double);
double Phi_1(double, double, double, double, double);
double Phi_qqg(double, double, double);
double uiPhi_0(double*, double*);
double uiPhi_1(double*, double*);
double Amplitude_0(double);
double Amplitude_1(double);
double uiAmplitude_0(double*, double*);
double uiAmplitude_1(double*, double*);
double uiPhi_qqg(const double*);
double uiAmplitude_qqg(double*, double*);
unsigned int mA;
unsigned int mQuarkIndex;
double mCrossSectionRatioLT;
DipoleModel* mDipoleModel;
+ DipoleModelParameterSet mDipoleModelParameterSet;
TF1 *Amp0, *Amp1, *Ampqqg;
ROOT::Math::WrappedTF1 *WFAmp0, *WFAmp1, *WFAmpqqg;
ROOT::Math::GaussIntegrator GIAmp0, GIAmp1, GIAmpqqg;
};
#endif

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 8:59 PM (23 h, 10 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242361
Default Alt Text
(117 KB)

Event Timeline