Page MenuHomeHEPForge
Files F10872023

No OneTemporary
Actions

View Options

Index: trunk/src/InclusiveDiffractiveCrossSections.h
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.h (revision 549)
+++ trunk/src/InclusiveDiffractiveCrossSections.h (revision 550)
@@ -1,145 +1,147 @@
//==============================================================================
// 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"
#include "Table.h"
#include "THn.h"
class DipoleModel;
class DipoleModelParameters;
class InclusiveDiffractiveCrossSections{
public:
InclusiveDiffractiveCrossSections();
virtual ~InclusiveDiffractiveCrossSections();
double operator()(double beta, double Q2, double W2, double z);
double operator()(const double* array);
DipoleModel* dipoleModel();
double unuranPDF(const double*);
double unuranPDF_qqg(const double*);
void setCheckKinematics(bool);
void setFockState(FockState);
GammaPolarization polarizationOfLastCall() ;
double crossSectionOfLastCall();
unsigned int quarkSpeciesOfLastCall();
double quarkMassOfLastCall();
double crossSectionRatioLTOfLastCall() const;
void setQuarkIndex(unsigned int);
virtual double dsigmadbetadz_T(double, double, double, double) ;
virtual double dsigmadbetadz_L(double, double, double, double) ;
virtual double dsigmadbetadz_QQG(double, double, double, double) ;
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);
double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;
double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
FockState getFockState();
virtual double* dsigdbetadQ2dWdz_T_total();
virtual double* dsigdbetadQ2dWdz_L_total();
virtual double* dsigdbetadQ2dWdz_T_qqg();
+ double uiAmplitude_1(double*, double*);
+ double uiAmplitude_0(double*, double*);
+ double uiAmplitude_qqg(double*, double*);
+
+
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, mTotalCS_qqg;
FockState mFockState;
double mS, mQ2;
bool mCheckKinematics;
unsigned int mA;
unsigned int mQuarkIndex;
double mCrossSectionRatioLT;
double mDsigdbetadQ2dWdz_L_total[5];
double mDsigdbetadQ2dWdz_T_total[5];
double mDsigdbetadQ2dWdz_T_qqg[5];
DipoleModel* mDipoleModel;
DipoleModelParameterSet mDipoleModelParameterSet;
};
class InclusiveDiffractiveCrossSectionsIntegrals: public InclusiveDiffractiveCrossSections{
public:
InclusiveDiffractiveCrossSectionsIntegrals();
~InclusiveDiffractiveCrossSectionsIntegrals();
double dsigmadbetadz_T(double, double, double, double) ;
double dsigmadbetadz_L(double, double, double, double) ;
double dsigmadbetadz_QQG(double, double, double, double) ;
protected:
private:
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*);
TF1 *Amp0, *Amp1, *Ampqqg;
ROOT::Math::WrappedTF1 *WFAmp0, *WFAmp1, *WFAmpqqg;
ROOT::Math::GaussIntegrator GIAmp0, GIAmp1, GIAmpqqg;
};
#endif
Index: trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 549)
+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 550)
@@ -1,306 +1,328 @@
#include "InclusiveDiffractiveCrossSectionsFromTables.h"
#include <TFile.h>
#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 "TableGeneratorSettings.h"
#include "Settings.h"
#include "Enumerations.h"
#include "DglapEvolution.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include "Kinematics.h"
#include "linterp.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
InclusiveDiffractionCrossSectionsFromTables::InclusiveDiffractionCrossSectionsFromTables(){
string path="/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/";
/*
TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
*/
//17x17x17 A1:
TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q0_bin0-83521_TQQ.root", "READ");
TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q1_bin0-83521_TQQ.root", "READ");
TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q2_bin0-83521_TQQ.root", "READ");
TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q3_bin0-83521_tQQ.root", "READ");
TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q0_bin0-83521_LQQ.root", "READ");
TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q1_bin0-83521_LQQ.root", "READ");
TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q2_bin0-83521_LQQ.root", "READ");
TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q3_bin0-83521_LQQ.root", "READ");
TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q0_bin0-83521_TQQG.root", "READ");
TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q1_bin0-83521_TQQG.root", "READ");
TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q2_bin0-83521_TQQG.root", "READ");
TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q3_bin0-83521_TQQG.root", "READ");
/*
// 10x10x10x10:
TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
*/
tableQQ_T.clear();
tableQQ_L.clear();
tableQQG.clear();
// Retrieve the histograms
THnF* histQQ_T0 = dynamic_cast<THnF*>(fileQQ_T0->Get("table"));
THnF* histQQ_T1 = dynamic_cast<THnF*>(fileQQ_T1->Get("table"));
THnF* histQQ_T2 = dynamic_cast<THnF*>(fileQQ_T2->Get("table"));
THnF* histQQ_T3 = dynamic_cast<THnF*>(fileQQ_T3->Get("table"));
THnF* histQQ_L0 = dynamic_cast<THnF*>(fileQQ_L0->Get("table"));
THnF* histQQ_L1 = dynamic_cast<THnF*>(fileQQ_L1->Get("table"));
THnF* histQQ_L2 = dynamic_cast<THnF*>(fileQQ_L2->Get("table"));
THnF* histQQ_L3 = dynamic_cast<THnF*>(fileQQ_L3->Get("table"));
THnF* histQQG_0 = dynamic_cast<THnF*>(fileQQG_0->Get("table"));
THnF* histQQG_1 = dynamic_cast<THnF*>(fileQQG_1->Get("table"));
THnF* histQQG_2 = dynamic_cast<THnF*>(fileQQG_2->Get("table"));
THnF* histQQG_3 = dynamic_cast<THnF*>(fileQQG_3->Get("table"));
tableQQ_T.push_back(histQQ_T0);
tableQQ_T.push_back(histQQ_T1);
tableQQ_T.push_back(histQQ_T2);
tableQQ_T.push_back(histQQ_T3);
tableQQ_L.push_back(histQQ_L0);
tableQQ_L.push_back(histQQ_L1);
tableQQ_L.push_back(histQQ_L2);
tableQQ_L.push_back(histQQ_L3);
tableQQG.push_back(histQQG_0);
tableQQG.push_back(histQQG_1);
tableQQG.push_back(histQQG_2);
tableQQG.push_back(histQQG_3);
}
InclusiveDiffractionCrossSectionsFromTables::~InclusiveDiffractionCrossSectionsFromTables(){
}
double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {
double MX2 = Q2*(1-beta)/beta;
double mf = Settings::quarkMass(mQuarkIndex);
- double sqrtArg = 1.-4.*mf*mf/MX2;
+ double pt2 = z*(1-z)*Q2-mf*mf;
+ double sqrtArg = 1.-4.*(mf*mf+pt2)/MX2;
if (sqrtArg<0){
- if (mSettings->verboseLevel() > 2)
- cout<<"There is no phase-space for z, return 0."<<endl;
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T: ";
+ cout<<"MX2<4(mf2_pt2), return 0."<<endl;
+ PR(MX2);
+ PR(4*mf*mf);
+ }
return 0;
}
- double z0 = (1.-sqrt(sqrtArg))/2.;
+ double z0 = (1.-sqrt(sqrtArg))/2.; //the limits of z depend on beta
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;
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T: ";
+ 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)
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T: ";
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
+ PR(z);
+ }
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;
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T: ";
+ cout<<"Not enough phase-space for quark production, z*(1-z)*MX2-mf*mf<0, return 0."<<endl;
+ }
return 0;
}
double result = interpolate(beta, Q2, W2, z, transverse, QQ, mQuarkIndex); //nb
return result;
}
double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L(double beta, double Q2, double W2, double z){
double MX2 = Q2*(1-beta)/beta;
double mf = Settings::quarkMass(mQuarkIndex);
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;
+ double sqrtArg = 1.-4.*(mf*mf+pt2)/MX2;
if (sqrtArg<0){
- if (mSettings->verboseLevel() > 2)
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L: ";
cout<<"There is no phase-space for z, return 0."<<endl;
+ PR(MX2);
+ PR(4*mf*mf);
+ }
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;
+ if (z<z0 or z>1-z0){ //the limits of z depend on beta
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L: ";
+ cout<<"z="<<z<<" which is smaller than z0="<<z0<<", or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
+ }
+ return 0;
+ }
+ if (pt2<0){
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L: ";
+ cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
+ PR(z);
+ }
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;
+ if (mSettings->verboseLevel() > 3 && mQuarkIndex==0){
+ cout<<"InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L: ";
+ cout<<"Not enough phase-space for quark production, z*(1-z)*MX2-mf*mf<0, return 0."<<endl;
+ }
return 0;
}
double result = interpolate(beta, Q2, W2, z, longitudinal, QQ, mQuarkIndex);
return result;
}
//===================================================
// QQG Calculations
//===================================================
double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
double result = interpolate(beta, Q2, W2, ztilde, transverse, QQG, mQuarkIndex);
return result;
}
double InclusiveDiffractionCrossSectionsFromTables::interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark){
THnF* hist4D=0;
//#TT TODO: There seem to be a strange swap happening
// between T QQ and L QQ...
if(pol==longitudinal and fock==QQ){
hist4D=tableQQ_T.at(quark);
}
else if(pol==transverse and fock==QQ){
hist4D=tableQQ_L.at(quark);
}
else if(fock==QQG){
hist4D=tableQQG.at(quark);
}
else{
cout<<"InclusiveDiffractionCrossSectionsFromTables::interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark): There is nothing to interpolate, ending"<<endl;
exit(2);
}
int nDims = hist4D->GetNdimensions();
vector <double>MinRange(nDims); //beta, Q2,W2, z
vector <double>MaxRange(nDims); //beta, Q2,W2, z
vector <int>Bins(nDims); //beta, Q2,W2, z
for (int dim = 0; dim < 4; dim++) {
TAxis* axis = hist4D->GetAxis(dim); // Get the axis for this dimension
if (!axis) {
std::cerr << "Error: Cannot retrieve axis for dimension " << dim << std::endl;
continue;
}
// Get the number of bins and range
int nBins = axis->GetNbins();
double minRange = axis->GetXmin();
double maxRange = axis->GetXmax();
Bins[dim]= nBins;
MinRange[dim]= minRange;
MaxRange[dim]= maxRange;
}
// construct the grid in each dimension.
// note that we will pass in a sequence of iterators pointing to the beginning of each grid
std::vector<double> gridBeta = mgrid(hist4D, 0, Bins[0]);
std::vector<double> gridQ2 = mgrid(hist4D, 1, Bins[1]) ;
std::vector<double> gridW2 = mgrid(hist4D, 2, Bins[2]);
std::vector<double> gridZ = mgrid(hist4D, 3, Bins[3]) ;
std::vector< std::vector<double>::iterator > grid_iter_list;
grid_iter_list.push_back(gridBeta.begin());
grid_iter_list.push_back(gridQ2.begin());
grid_iter_list.push_back(gridW2.begin());
grid_iter_list.push_back(gridZ.begin());
// the size of the grid in each dimension
array<int,4> grid_sizes;
for(unsigned long i=0; i<grid_sizes.size(); i++){
grid_sizes[i] = Bins[i];
}
// total number of elements
int num_elements = grid_sizes[0] * grid_sizes[1]* grid_sizes[2]* grid_sizes[3];
// fill in the values of f(x) at the gridpoints.
// we will pass in a contiguous sequence, values are assumed to be laid out C-style
std::vector<double> f_values(num_elements);
int count = 0;
int countFail=0;
for (int ibeta=1; ibeta<=grid_sizes[0]; ibeta++) {
for (int iQ2=1; iQ2<=grid_sizes[1]; iQ2++) {
for (int iW2=1; iW2<=grid_sizes[2]; iW2++) {
for (int iZ=1; iZ<=grid_sizes[3]; iZ++) {
int binIndex[4] = {ibeta-1, iQ2-1, iW2-1, iZ-1};
double val=hist4D->GetBinContent(binIndex);
// There seems to be some spikes
// This has to be adapted to lin/log content
// if(exp(val)>1e20){
// val=0;
// }
if(isnan(val) or isinf(val)){
countFail++;
val=0;
}
if(val!=0)
f_values[count] = val;//exp(val);
else
f_values[count] = 0.;
count++;
}
}
}
}
//
// Construct the interpolator.
// The last two arguments are pointers to the underlying data
//
InterpMultilinear <4, double> interp_ML(grid_iter_list.begin(), grid_sizes.begin(), f_values.data(), f_values.data() + num_elements);
// interpolate:
array<double,4> args = {beta, Q2, W2, z};
double result=interp_ML.interp(args.begin());
// return result; //For linear content
if(result==0)
return 0.;
else
return exp(result); //For log content
}
vector<double> InclusiveDiffractionCrossSectionsFromTables::mgrid(THnF* hist,int axis, int len){
vector<double> gridtype(len);
for(int i = 1; i <=len; i++){
double val = hist->GetAxis(axis)->GetBinCenter(i);
gridtype[i-1] = val;
}
return gridtype;
}
Index: trunk/src/Event.h
===================================================================
--- trunk/src/Event.h (revision 549)
+++ trunk/src/Event.h (revision 550)
@@ -1,85 +1,86 @@
//==============================================================================
// Event.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 Event_h
#define Event_h
#include <vector>
#include <string>
#include <iostream>
#include "TLorentzVector.h"
#include "Enumerations.h"
using namespace std;
class Particle {
public:
int index; // starts at 0 equals index in particle vector
int pdgId; // particle ID according to PDG scheme
int status; // 1 = not decayed/final, 2 = decayed
TLorentzVector p;
vector<int> parents;
vector<int> daughters;
};
class Event {
public:
unsigned long eventNumber;
//
// Event kinematics
//
double Q2;
double W;
double t;
double x;
double s;
double y;
double xpom;
double beta;
double z; // new for inclusive <
double MX; // new for inclusive
double Omega; // new for inclusive
unsigned int quarkSpecies; // new for inclusive
-
- //
+ double QS; //saturation scale
+
+ //
// Event traits
//
GammaPolarization polarization;
DiffractiveMode diffractiveMode;
double crossSectionRatioLT;
//
// List of particles in event.
// First two are always beam particles.
//
vector<Particle> particles;
int numberOfPythiaParticlesStored;
public:
void list(ostream& = cout) const;
};
#endif
Index: trunk/src/InclusiveFinalStateGenerator.h
===================================================================
--- trunk/src/InclusiveFinalStateGenerator.h (revision 549)
+++ trunk/src/InclusiveFinalStateGenerator.h (revision 550)
@@ -1,47 +1,53 @@
//==============================================================================
// ExclusiveFinalStateGenerator.h
//
// Copyright (C) 2024 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
// $Author: ullrich $
//==============================================================================
#ifndef InclusiveFinalStateGenerator_h
#define InclusiveFinalStateGenerator_h
#include "Pythia8/Pythia.h"
#include "FinalStateGenerator.h"
#include "Constants.h"
#include "Enumerations.h"
class InclusiveFinalStateGenerator : public FinalStateGenerator {
public:
InclusiveFinalStateGenerator();
~InclusiveFinalStateGenerator();
bool generate(int A, Event *event, FockState fockstate);
double uiGamma_N(double* var, double* par) const;
double gamma_N(unsigned int i, double val);
private:
Pythia8::Pythia *mPythia;
Pythia8::Event *mPythiaEvent;
Pythia8::ParticleData *mPdt;
- double computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m, double M);
+ double computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m_1, double m_2, double M, double, double);
+
+ bool decayPseudoParticle(TLorentzVector, TLorentzVector&, TLorentzVector&, double, double, double);
+
+ double evaluate(double phi_totg, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus, double mf2);
+ bool bisection(double phi_min, double phi_max, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus,
+ double tol, double& phi_root, double mf2);
};
#endif
Index: trunk/src/InclusiveFinalStateGenerator.cpp
===================================================================
--- trunk/src/InclusiveFinalStateGenerator.cpp (revision 549)
+++ trunk/src/InclusiveFinalStateGenerator.cpp (revision 550)
@@ -1,773 +1,1117 @@
//==============================================================================
// InclusiveFinalStateGenerator.cpp
//
// Copyright (C) 2024 Tobias Toll and Thomas Ullrich
//
-// This file is part of Sartre.
+// 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.
+// 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/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 <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()
+#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()
+
+InclusiveFinalStateGenerator::~InclusiveFinalStateGenerator()
{
// mPythia->stat();
delete mPythia;
}
-
-double InclusiveFinalStateGenerator::uiGamma_N(double* var, double* par) const
+
+double InclusiveFinalStateGenerator::uiGamma_N(double* var, double* par) const
{
double t=var[0];
double s=par[0];
double result = pow(t, s-1)*exp(-t);
return result;
}
double InclusiveFinalStateGenerator::gamma_N(unsigned int i, double val)
{
TF1 fGamma_N("fGamma_N", this, &InclusiveFinalStateGenerator::uiGamma_N, -10, 10, 1);
fGamma_N.SetParameter(0, i);
ROOT::Math::WrappedTF1 wfGamma_N(fGamma_N);
ROOT::Math::GaussIntegrator giGamma_N;
giGamma_N.SetFunction(wfGamma_N);
giGamma_N.SetAbsTolerance(0.);
giGamma_N.SetRelTolerance(1e-5);
int sign=1;
double low=0, up=val;
if(val<0){
low=val;
up=0;
sign=-1;
}
double Tfact=1;
for(int j=i-1; j>0; j--){
Tfact*=j;
}
return sign*giGamma_N.Integral(low, up)/Tfact;
}
+bool InclusiveFinalStateGenerator::decayPseudoParticle(TLorentzVector In, TLorentzVector& Out1, TLorentzVector& Out2, double m_1, double m_2, double pt=0){
+ //
+ // Calculate the momenta in the In particle's restframe
+ //
+ EventGeneratorSettings *settings = EventGeneratorSettings::instance();
+ TRandom3 *rndm = settings->randomGenerator();
+
+ TVector3 boost=TVector3(In.Px()/In.E(), In.Py()/In.E(), In.Pz()/In.E());
+ double M2=In.M2();
+ double m2_1=m_1*m_1;
+ double m2_2=m_2*m_2;
+ //In the rest-frame:
+ double p2=2*m2_1*m2_2 + 2*(m2_1+m2_2)*M2 - m2_1*m2_1 - m2_2*m2_2 - M2*M2;
+ p2/=-4*M2;
+ if(p2<0) return false;
+ double p=sqrt(p2);
+ double E_1, E_2, px, py, pz;
+ if(pt==0){ //no user preference
+ //Generate a direction for decay
+ double phi=rndm->Uniform(2*M_PI);
+ double cosTheta=rndm->Uniform(2)-1;
+ double theta=acos(cosTheta);
+
+ px=p*cos(phi)*sin(theta);
+ py=p*sin(phi)*sin(theta);
+ pz=p*cosTheta;
+
+ E_1=sqrt(m2_1+p2);
+ E_2=sqrt(m2_2+p2);
+ }
+ else{ //user defined pt
+ double phi=rndm->Uniform(2*M_PI);
+ px=pt*cos(phi);
+ py=pt*sin(phi);
+ pz=sqrt(p2-px*px-py*py);
+
+ E_1=sqrt(m2_1+p2);
+ E_2=sqrt(m2_2+p2);
+ }
+ Out1=TLorentzVector(px, py, pz, E_1);
+ Out2=TLorentzVector(-px, -py, -pz, E_2);
+
+ Out1.Boost(boost);
+ Out2.Boost(boost);
+
+ return true;
+}
+
bool InclusiveFinalStateGenerator::generate(int A, Event *event, FockState fockstate)
{
//
// Get generator settings and the random generator
//
EventGeneratorSettings *settings = EventGeneratorSettings::instance();
TRandom3 *rndm = settings->randomGenerator();
-// 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 Pplus=mHadronBeam.E()+mHadronBeam.Pz();
+ double PplusPrime=(1-xpom)*Pplus;
+ double PminusPrime=(mMY2-mT)/PplusPrime;
+ E=(PplusPrime+PminusPrime)/2;
+ pz=(PplusPrime-PminusPrime)/2;
double pt_squared = E*E-pz*pz-mMY2;
if (pt_squared < 0) {
- if (settings->verbose() && settings->verboseLevel() > 1) {
+ if (settings->verbose() && settings->verboseLevel() > 2) {
cout<< "InclusiveFinalStateGenerator::generate(): No phase-space for scattered proton pt." << endl;
cout<< " Abort event generation." << endl;
}
return false;
}
pt = sqrt(pt_squared);
px = pt*cos(phi);
py = pt*sin(phi);
TLorentzVector theScatteredProton(px, py, pz, E);
//
// Use scattered proton to set up four momentum of pomeron:
//
TLorentzVector thePomeron=mHadronBeam-theScatteredProton;
-
-
+
+
//
// Finally the diffractive final state,
// treat at first as a pseudo particle
//
TLorentzVector theXparticle((mHadronBeam + mElectronBeam) - (theScatteredElectron + theScatteredProton));
//
// The quark and anti quark
//
TLorentzVector theQuark;
TLorentzVector theAntiQuark;
TLorentzVector theGluon; //for QQG
double z=event->z;
double pxq=0., pyq=0., pxqbar=0., pyqbar=0.;
double mf = Settings::quarkMass(event->quarkSpecies);
- PR(mf);
if(fockstate==QQ){
+ //
+ // Decay the X-particle using the correct pt
+ //
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 << "InclusiveFinalStateGenerator::generate(): 4*(mf*mf+pt2) > MX2, for quark flavour "<<event->quarkSpecies << 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);
+ if(!decayPseudoParticle(theXparticle, theQuark, theAntiQuark, mf, mf, sqrt(pt2))){
+ if (settings->verboseLevel() > 2) {
+ cout << "InclusiveFinalStateGenerator::generate(): No phase-space for quark pt" << endl;
+ cout << " Abort event generation" << endl;
+ }
+ return false;
+ }
+ if(z<0.5 && theQuark.Pz()<theAntiQuark.Pz()){
+ TLorentzVector tmp=theQuark;
+ theQuark=theAntiQuark;
+ theAntiQuark=tmp;
+ }
}
else if(fockstate==QQG){
//The knowns:
double Mqq2=(z/beta-1)*Q2;
if(Mqq2<4*mf*mf){
- if (settings->verboseLevel() > 1) {
+ if (settings->verboseLevel() > 2) {
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.
+ // Here pomeron comes from + direction, so we reverse the axis. (xi=beta/z)
+ // They make the assumption that the photon does not
+ // carry any +momentum, therefore all is momentum fractions of
+ // the pomeron +momentum. However, for us, the photon does carry
+ // a little bit of +momentum. Therefore we make momentum fractions
+ // of theXparticle instead.
//
- double PpomPlus=thePomeron.E()+thePomeron.Pz();
- double PgPlus=(1-z)*PpomPlus; //gluon
- double PqPlus=beta*PpomPlus; //quark
- double PqbarPlus=(z-beta)*PpomPlus; //antiquark
+ double PXPlus=theXparticle.E()+theXparticle.Pz();
+ double PgPlus=(1-z)*PXPlus; //gluon
+ double PqPlus=beta*PXPlus; //quark
+ double PqbarPlus=(z-beta)*PXPlus; //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
+ // Once we know the angle between X and g all else solves.
+ // However, finding this angle is tricky.
+ // We use a numerical bisection method.
//
- // 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_totg;
+ bool found=false;
+ // First find the brackets with a rough scan:
+ double valold=evaluate(0., theXparticle, z, Mqq2, MX2, PgPlus, 0.);
+ double valnew=0;
+ double phi_low, phi_high;
+ double stepsize=2*M_PI/100.;
+ for(double phi=0; phi<2*M_PI; phi+=stepsize){
+ valnew=evaluate(phi, theXparticle, z, Mqq2, MX2, PgPlus, 0.);
+ if(valold*valnew==0){
+ phi_totg=phi;
+ cout<<"Found it!"<<endl;
+ found=true;
}
- 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;
+ if(valold*valnew<0){
+ phi_high=phi;
+ phi_low=phi-stepsize;
+ break;
+ }
+ valold=valnew;
+ }
+ if(!found){
+ if(!bisection(phi_high, phi_low, theXparticle, z, Mqq2, MX2, PgPlus, 1e-5, phi_totg, 0.)){
+ // If the bisection fails we do a scan for best value.
+ // It could mean that there is no valid solution...
+ valold=fabs(evaluate(0., theXparticle, z, Mqq2, MX2, PgPlus, 0.));
+ valnew=0;
+ double eps=1e-4;
+ double stepsize=2*M_PI/100.;
+ while(true){
+ phi+=stepsize;
+ if(phi>=2*M_PI){
+ cout<<"no value found"<<endl;
+ return false;
+ }
+ valnew=fabs(evaluate(phi, theXparticle, z, Mqq2, MX2, PgPlus, 0.));
+ if(valnew>valold){
+ stepsize=-stepsize/10.;
+ }
+ if(fabs(stepsize)<eps){
+ phi_totg=phi;
+ break;
+ }
+ valold=valnew;
+ }
}
- return false;
}
- double phi_3g=acos(cosphi34); //angle between gluon and qqbar
- double phi_tot=acos((theVirtualPhoton+thePomeron).Px() /(theVirtualPhoton+thePomeron).Pt());
+ double pt=theXparticle.Pt();
+ double b=2*(1-z)*cos(phi_totg)*pt;
+ double c=(1-z)*(1-z)*pt*pt+(1-z)*Mqq2-z*(1-z)*MX2;
+ double ptg=-b/2*(1-sqrt(1-4*c/b/b));
+ double phi_tot=acos(theXparticle.Px()/pt);
+ double PgMinus=ptg*ptg/PgPlus;
//Fill the gluon 4-vector:
E=(PgPlus+PgMinus)/2;
pz=(PgPlus-PgMinus)/2;
- px=ptg*cos(phi_tot-phi_totg);
- py=ptg*sin(phi_tot-phi_totg);
+ 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;
+ //Call the pseudo gluon Y
+ //TODO: Check that we get the correct P+ momenta, may have to use the same method again as above...
+ //My original method is not working. For now, therefore,
+ //just decay the pseudo gluon into the qq-bar without
+ //caring about the correct momentum fractions:
+ TLorentzVector thePseudoGluon = theXparticle - theGluon;
+ decayPseudoParticle(thePseudoGluon, theQuark, theAntiQuark, mf, mf);
+
+// TLorentzVector thePseudoGluon = theXparticle - theGluon;
+// double xi=beta/z;
+// // Angle between q and Y:
+// double phi_qY;
+// found=false;
+// // First find the brackets with a rough scan:
+// valold=evaluate(0., thePseudoGluon, xi, mf*mf, Mqq2, PqPlus, mf*mf);
+// valnew=0;
+// stepsize=2*M_PI/100.;
+// for(double phi=0; phi<2*M_PI; phi+=stepsize){
+// valnew=evaluate(phi, thePseudoGluon, xi, mf*mf, Mqq2, PqPlus, mf*mf);
+// if(valold*valnew==0){
+// phi_qY=phi;
+// cout<<"Found it!"<<endl;
+// found=true;
+// }
+// if(valold*valnew<0){
+// phi_high=phi;
+// phi_low=phi-stepsize;
+// break;
+// }
+// valold=valnew;
+// }
+// if(!found){
+// if(!bisection(phi_high, phi_low, theXparticle, z, Mqq2, MX2, PqPlus, 1e-5, phi_qY, mf*mf)){
+// cout<<"bisection 2 failed"<<endl;
+// // If the bisection fails we do a scan for best value.
+// // It could mean that there is no valid solution...
+// valold=fabs(evaluate(0., thePseudoGluon, xi, mf*mf, Mqq2, PqPlus, mf*mf));
+// valnew=0;
+// double eps=1e-4;
+// double stepsize=2*M_PI/10.;
+// phi=0;
+// double valmin=valold;
+// double phimin=0;
+// double upperBracket=2*M_PI;
+// while(true){
+// phi+=stepsize;
+// valnew=fabs(evaluate(phi, thePseudoGluon, xi, mf*mf, Mqq2, PqPlus, mf*mf));
+// if(valnew<valmin){
+// valmin=valnew;
+// phimin=phi;
+// }
+// if(phi>upperBracket){
+// phi=phimin-2*stepsize;
+// upperBracket=phimin+stepsize;
+// stepsize/=10.;
+// valmin=valnew;
+// }
+// PR(valnew);
+// PR(phi);
+// if(fabs(stepsize)<eps){
+// phi_qY=phi;
+// break;
+// }
+// valold=valnew;
+// }
+// PR(valnew/mf/mf);
+// }
+// }
+// pt=thePseudoGluon.Pt();
+// b=2*(1-xi)*cos(phi_qY)*pt;
+// c=(1-xi)*(1-xi)*pt*pt+mf*mf-xi*(1-xi)*Mqq2;
+// double ptq=-b/2*(1-sqrt(1-4*c/b/b));
+// phi_tot=acos(thePseudoGluon.Px()/pt);
+// double PqMinus=ptq*ptq/PqPlus;
+//
+// //Fill the quark 4-vector:
+// E=(PqPlus+PqMinus)/2;
+// pz=(PqPlus-PqMinus)/2;
+// px=ptq*cos(phi_tot+phi_qY);
+// py=ptq*sin(phi_tot+phi_qY);
+// theQuark=TLorentzVector(px, py, pz, E);
+// theAntiQuark=thePseudoGluon-theQuark;
+
+// double k2=(thePomeron-theGluon).M2();
+// double pt2=xi*(1-xi)*(-k2)-mf*mf;
+// if(!decayPseudoParticle(thePseudoGluon, theQuark, theAntiQuark, mf, mf, sqrt(pt2))){
+// if (settings->verboseLevel() > 2) {
+// cout << "InclusiveFinalStateGenerator::generate(): No phase-space for quark pt" << endl;
+// cout << " Abort event generation" << endl;
+// }
+// return false;
+// }
+/*
+
phi = rndm->Uniform(twopi);
- pxq=sqrt(pt2)*sin(phi)/2;
- pyq=sqrt(pt2)*cos(phi)/2;
+ pxq=sqrt(pt2)*sin(phi);
+ pyq=sqrt(pt2)*cos(phi);
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);
+// PR(PpomPlus);
+ PR(PgPlus);
+ PR((theGluon.E()+theGluon.Pz()));
+ PR(PqPlus);
+ PR(theQuark.E()+theQuark.Pz());
+ PR(PqbarPlus);
+ PR(theAntiQuark.E()+theAntiQuark.Pz());
+ */
}
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
+ // 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;
+ double sum=0;
+ for(int n=1; n<100; n++){
+ int nfactorial=1;
+ for(int i=1; i<=n; i++){
+ nfactorial*=i;
+ }
+ sum+=pow(-1, n-1)*2/nfactorial*pow(Omega/2, n);
+ double trueVal=2*(1-exp(-Omega/2));
+ if(abs(sum-trueVal)<eps){
+ Twist=n;
break;
}
}
+ /*
+ Twist=0;
+ for(int i=0; i<100; i++){
+ double gammaN=gamma_N(i+1, -Omega/2);
+ if(fabs(2*exp(-Omega/2)*gammaN)<eps){
+ Twist=i+1;
+ break;
+ }
+ }
+ PR(Twist);
+ exit(1);
+ */
//
// 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 Qs=event->QS;
+ vector<double> pom_px, pom_py, pom_pz;
double pom_px_tot=thePomeron.Px();
double pom_py_tot=thePomeron.Py();
- // double pom_pz_tot=thePomeron.Pz();
- for(int i=0; i<2*Twist-1; i++){
+ double pom_pz_tot=thePomeron.Pz();
+ for(int i=0; i<2*Twist-1; i++){ //-1 since the last gluon conserves total momentum
//
// Generate momenta from a Gaussian with width Qs
//
double p_x=rndm->Gaus(0., Qs);
double p_y=rndm->Gaus(0., Qs);
- // double p_z=rndm->Gaus(0., Qs);
+ 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.push_back(thePomeron.Px()/(2*Twist)+p_x);
+ pom_py.push_back(thePomeron.Py()/(2*Twist)+p_y);
+ pom_pz.push_back(thePomeron.Pz()/(2*Twist)+p_z);
pom_px_tot+=p_x;
pom_py_tot+=p_y;
+ pom_pz_tot+=p_z;
}
//conserve the total pomeron momentum:
pom_px.push_back(thePomeron.Px()-pom_px_tot);
pom_py.push_back(thePomeron.Py()-pom_py_tot);
+ pom_pz.push_back(thePomeron.Pz()-pom_pz_tot);
//
// 3. Distribute the recoils to:
- // the quark and antiquark (QQ) or gluon and pseudogluon (QQG):
- // Call these q and qbar (for parton)
+ // the quark and antiquark (QQ) or
+ // the quark, antiquarkgluon and gluon (QQG).
+ // Start with QQ
//
double Eq=0, Eqbar=0, pzq=0, pzqbar=0;
if(fockstate==QQ){
Eq=theQuark.E();
Eqbar=theAntiQuark.E();
pzq=theQuark.Pz();
pzqbar=theAntiQuark.Pz();
pxq=theQuark.Px();
pyq=theQuark.Py();
pxqbar=theAntiQuark.Px();
pyqbar=theAntiQuark.Py();
+ for(unsigned int i=0; i<pom_px.size(); i++){
+ //
+ //First calculate the relative velocity between q, qbar and gluon
+ //
+ double pomP=sqrt(pom_px[i]*pom_px[i] + pom_py[i]*pom_py[i] + pom_pz[i]*pom_pz[i]);
+ double qP=sqrt(pxq*pxq + pyq*pyq + pzq*pzq);
+ double qbarP=sqrt(pxqbar*pxqbar + pyqbar*pyqbar + pzqbar*pzqbar);
+ double deltaVq=sqrt(1+qP*qP/Eq/Eq - 2*(pom_px[i]*pxq + pom_py[i]*pyq + pom_pz[i]*pzq) /Eq/pomP);
+ double deltaVqbar=sqrt(1+qbarP*qbarP/Eqbar/Eqbar - 2*(pom_px[i]*pxqbar + pom_py[i]*pyqbar + pom_pz[i]*pzqbar) / Eqbar/pomP);
+ //The momentum kick is inversely proportional to relative velocity
+ double fraction=0.5;
+ if(deltaVq+deltaVqbar!=0)
+ fraction=deltaVqbar/(deltaVq+deltaVqbar);
+ //Update momenta
+ pxq+=fraction*pom_px[i];
+ pyq+=fraction*pom_py[i];
+ pzq+=fraction*pom_pz[i];
+ pxqbar+=(1-fraction)*pom_px[i];
+ pyqbar+=(1-fraction)*pom_py[i];
+ pzqbar+=(1-fraction)*pom_pz[i];
+
+ //Scale the transverse momenta in order to preserve MX2,
+ //First find the scaling factor lambda
+ double lower_bracket=0, upper_bracket=100;
+ double lambda=computeScalingFactor(pxq, pyq, pzq, pxqbar, pyqbar, pzqbar, mf, mf, sqrt(MX2), lower_bracket, upper_bracket);
+ if(lambda<=lower_bracket+1e-5 || lambda>=upper_bracket-1e-5) {
+ //failed to find the root
+ if (settings->verboseLevel() > 2) {
+ cout<<"InclusiveFinalStateGenerator::generate: Failed to scatter the quarks"<<endl;
+ }
+ return false;
+ }
+ if(lambda==(lower_bracket+upper_bracket)/2)
+ PR(lambda);
+ pxq*=lambda;
+ pyq*=lambda;
+ pzq*=lambda;
+ pxqbar*=lambda;
+ pyqbar*=lambda;
+ pzqbar*=lambda;
+ //Put on shell:
+ Eq=sqrt(pxq*pxq+pyq*pyq+pzq*pzq+mf*mf);
+ Eqbar=sqrt(pxqbar*pxqbar+pyqbar*pyqbar+pzqbar*pzqbar+mf*mf);
+ }
+ //
+ //Update the 4-vectors:
+ //
+ //the Quark:
+ theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
+ //the Anti-Quark
+ theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
}
else if(fockstate==QQG){
- Eq=theQuark.E()+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);
- }
- //
- //Update the 4-vectors:
- //
- if(fockstate==QQ){
+ Eq=theQuark.E();
+ pzq=theQuark.Pz();
+ pxq=theQuark.Px();
+ pyq=theQuark.Py();
+ Eqbar=theAntiQuark.E();
+ pzqbar=theAntiQuark.Pz();
+ pxqbar=theAntiQuark.Px();
+ pyqbar=theAntiQuark.Py();
+ double Eg=theGluon.E();
+ double pzg=theGluon.Pz();
+ double pxg=theGluon.Px();
+ double pyg=theGluon.Py();
+
+ for(unsigned int i=0; i<pom_px.size(); i++){
+ //
+ //First calculate the relative velocity between q, qbar and gluon
+ //
+ double pomP=sqrt(pom_px[i]*pom_px[i] + pom_py[i]*pom_py[i] + pom_pz[i]*pom_pz[i]);
+ double qP=sqrt(pxq*pxq + pyq*pyq + pzq*pzq);
+ double qbarP=sqrt(pxqbar*pxqbar + pyqbar*pyqbar + pzqbar*pzqbar);
+ double gP=sqrt(pxg*pxg + pyg*pyg + pzg*pzg);
+ double deltaVg=sqrt(1+gP*gP/Eg/Eg - 2*(pom_px[i]*pxg + pom_py[i]*pyg + pom_pz[i]*pzg) /Eg/pomP);
+ //Figure out which dipole the gluon is interacting with,
+ // g-q or q-qbar:
+ double p2_gq=(pxq+pxg)*(pxq+pxg) + (pyq+pyg)*(pyq+pyg) + (pzq+pzg)*(pzq+pzg);
+ double p2_gqbar=(pxqbar+pxg)*(pxqbar+pxg) + (pyqbar+pyg)*(pyqbar+pyg) + (pzqbar+pzg)*(pzqbar+pzg);
+ //Interaction probability propto dipole size:
+ bool isGQbar = p2_gq/(p2_gq+p2_gqbar) > rndm->Uniform(1);
+ if(isGQbar){
+ //Interacting with the g-qbar dipole:
+ double deltaVqbar=sqrt(1+qbarP*qbarP/Eqbar/Eqbar - 2*(pom_px[i]*pxqbar + pom_py[i]*pyqbar + pom_pz[i]*pzqbar) / Eqbar/pomP);
+ //The momentum kick is inversely proportional to relative velocity
+ double fraction=0.5;
+ if(deltaVg+deltaVqbar!=0)
+ fraction=deltaVg/(deltaVg+deltaVqbar);
+ //Update momenta
+ pxg+=fraction*pom_px[i];
+ pyg+=fraction*pom_py[i];
+ pzg+=fraction*pom_pz[i];
+ pxqbar+=(1-fraction)*pom_px[i];
+ pyqbar+=(1-fraction)*pom_py[i];
+ pzqbar+=(1-fraction)*pom_pz[i];
+ }
+ else{
+ //Interacting with the g-q dipole:
+ double deltaVq=sqrt(1+qP*qP/Eq/Eq - 2*(pom_px[i]*pxq + pom_py[i]*pyq + pom_pz[i]*pzq) /Eq/pomP);
+ //The momentum kick is inversely proportional to relative velocity
+ double fraction=0.5;
+ if(deltaVg+deltaVq!=0)
+ fraction=deltaVq/(deltaVq+deltaVg);
+ //Update momenta
+ pxg+=fraction*pom_px[i];
+ pyg+=fraction*pom_py[i];
+ pzg+=fraction*pom_pz[i];
+ pxq+=(1-fraction)*pom_px[i];
+ pyq+=(1-fraction)*pom_py[i];
+ pzq+=(1-fraction)*pom_pz[i];
+ }
+
+ //Scale the transverse momenta in order to preserve MX2 and Mqq2,
+ //First find the scaling factor lambda
+ // start by scaling the quarks:
+ double Mqq2=(z/beta-1)*Q2;
+ double lower_bracket=0, upper_bracket=100;
+ double lambda=computeScalingFactor(pxq, pyq, pzq, pxqbar, pyqbar, pzqbar, mf, mf, sqrt(Mqq2), lower_bracket, upper_bracket);
+ if(lambda<=lower_bracket+1e-5 || lambda>=upper_bracket-1e-5) {
+ if (settings->verboseLevel() > 2) {
+ cout<<"InclusiveFinalStateGenerator::generate: Failed to scatter the gluons"<<endl;
+ }
+ return false;
+ }
+ pxq*=lambda;
+ pyq*=lambda;
+ pzq*=lambda;
+ pxqbar*=lambda;
+ pyqbar*=lambda;
+ pzqbar*=lambda;
+ //Put on shell:
+ Eq=sqrt(pxq*pxq+pyq*pyq+pzq*pzq+mf*mf);
+ Eqbar=sqrt(pxqbar*pxqbar+pyqbar*pyqbar+pzqbar*pzqbar+mf*mf);
+
+ //
+ // Create the pseudo gluon based on the current quarks
+ //
+ theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
+ theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
+ TLorentzVector thePseudoGluon=theQuark+theAntiQuark;
+
+ double pxpg=thePseudoGluon.Px();
+ double pypg=thePseudoGluon.Py();
+ double pzpg=thePseudoGluon.Pz();
+ double Mpg=thePseudoGluon.M();
+ // Calculate momentum fractions for later:
+ double pxFraction=pxq/pxpg;
+ double pyFraction=pyq/pypg;
+ double pzFraction=pzq/pzpg;
+
+ lower_bracket=0; upper_bracket=100;
+ lambda=computeScalingFactor(pxpg, pypg, pzpg, pxg, pyg, pzg, Mpg, 0, sqrt(MX2), lower_bracket, upper_bracket);
+ if(lambda<=lower_bracket+1e-5 || lambda>=upper_bracket-1e-5) {
+ if (settings->verboseLevel() > 2) {
+ cout<<"InclusiveFinalStateGenerator::generate: Failed to scatter the quarks"<<endl;
+ }
+ return false;
+ }
+ if(lambda==(lower_bracket+upper_bracket)/2)
+ pxpg*=lambda;
+ pypg*=lambda;
+ pzpg*=lambda;
+
+ pxq=pxFraction*pxpg;
+ pyq=pyFraction*pypg;
+ pzq=pzFraction*pzpg;
+ Eq=sqrt(pxq*pxq+pyq*pyq+pzq*pzq+mf*mf);
+
+ pxqbar=(1-pxFraction)*pxpg;
+ pyqbar=(1-pyFraction)*pypg;
+ pzqbar=(1-pzFraction)*pzpg;
+ Eqbar=sqrt(pxqbar*pxqbar+pyqbar*pyqbar+pzqbar*pzqbar+mf*mf);
+
+ pxg*=lambda;
+ pyg*=lambda;
+ pzg*=lambda;
+ Eg=sqrt(pxg*pxg+pyg*pyg+pzg*pzg);
+ }
+ //
+ //Update the 4-vectors:
+ //
//the Quark:
theQuark=TLorentzVector(pxq, pyq, pzq, Eq);
//the Anti-Quark
theAntiQuark=TLorentzVector(pxqbar, pyqbar, pzqbar, Eqbar);
+ //the Gluon:
+ theGluon=TLorentzVector(pxg, pyg, pzg, Eg);
}
else{
- //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);
+ cout<<"InclusiveFinalStateGenerator::generate: unknown fockstate: "<<fockstate<<" ending program."<<endl;
+ exit(9);
}
//
// Decorrelations ended
//
- */
+
//
// Add particles to event record
//
double ifock;
if(fockstate==QQ)
ifock=0;
else
ifock=1;
event->particles.resize(2+5+ifock);
unsigned int eOut = 2;
unsigned int gamma = 3;
unsigned int hOut = 4;
unsigned int quark = 5;
unsigned int antiquark = 6;
unsigned int gluon = 7;
// Global indices
event->particles[eOut].index = eOut;
event->particles[gamma].index = gamma;
event->particles[hOut].index = hOut;
event->particles[quark].index = quark;
event->particles[antiquark].index = antiquark;
if(fockstate==QQG)
event->particles[antiquark].index = gluon;
-
+
// 4-vectors
event->particles[eOut].p = theScatteredElectron;
event->particles[hOut].p = theScatteredProton;
event->particles[gamma].p = theVirtualPhoton;
event->particles[quark].p = theQuark;
event->particles[antiquark].p = theAntiQuark;
if(fockstate==QQG)
event->particles[gluon].p = theGluon;
-
+
// PDG Ids
event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
event->particles[gamma].pdgId = 22;
event->particles[quark].pdgId = event->quarkSpecies+1;
event->particles[antiquark].pdgId = -event->particles[quark].pdgId;
if(fockstate==QQG)
event->particles[gluon].pdgId = 21;
-
+
// status
//
// HepMC conventions (February 2009).
// 0 : an empty entry, with no meaningful information
// 1 : a final-state particle, i.e. a particle that is not decayed further by
// the generator (may also include unstable particles that are to be decayed later);
// 2 : a decayed hadron or tau or mu lepton
// 3 : a documentation entry (not used in PYTHIA);
// 4 : an incoming beam particle;
// 11 - 200 : an intermediate (decayed/branched/...) particle that does not
// fulfill the criteria of status code 2
-
+
event->particles[ePos].status = 4;
event->particles[hPos].status = 4;
event->particles[eOut].status = 1;
event->particles[hOut].status = mIsIncoherent ? 2 : 1;
event->particles[gamma].status = 2;
event->particles[quark].status = 2;
event->particles[antiquark].status = 2;
if(fockstate==QQG)
event->particles[gluon].status = 2;
-
+
// parents (ignore dipole)
event->particles[eOut].parents.push_back(ePos);
event->particles[gamma].parents.push_back(ePos);
event->particles[hOut].parents.push_back(hPos);
event->particles[quark].parents.push_back(gamma);
event->particles[antiquark].parents.push_back(gamma);
if(fockstate==QQG){
event->particles[gluon].parents.push_back(quark);
event->particles[gluon].parents.push_back(antiquark);
}
// daughters (again ignore dipole)
event->particles[ePos].daughters.push_back(eOut);
event->particles[ePos].daughters.push_back(gamma);
event->particles[gamma].daughters.push_back(quark);
event->particles[gamma].daughters.push_back(antiquark);
event->particles[hPos].daughters.push_back(hOut);
if(fockstate==QQG){
event->particles[quark].daughters.push_back(gluon);
event->particles[antiquark].daughters.push_back(gluon);
}
-
-
+
+
//fill event structure
double y=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();
// Pythia starts with u = 1, Sartre with u = 0
int quarkID = event->quarkSpecies + 1;
int status = 23;
int color = 101; // 101 102 103
-
+
if(fockstate==QQ){
mPythiaEvent->append( quarkID, status, color, 0,
theQuark.Px(), theQuark.Py(), theQuark.Pz(),
theQuark.E(), theQuark.M() );
mPythiaEvent->append(-quarkID, status, 0, color,
theAntiQuark.Px(), theAntiQuark.Py(), theAntiQuark.Pz(),
theAntiQuark.E(), theAntiQuark.M());
}
else if(fockstate==QQG){
int acolor=102;
mPythiaEvent->append( quarkID, status, color, 0,
theQuark.Px(), theQuark.Py(), theQuark.Pz(),
theQuark.E(), theQuark.M() );
mPythiaEvent->append(-quarkID, status, 0, acolor,
theAntiQuark.Px(), theAntiQuark.Py(), theAntiQuark.Pz(),
theAntiQuark.E(), theAntiQuark.M());
mPythiaEvent->append(21, status, acolor, color,
theGluon.Px(), theGluon.Py(), theGluon.Pz(),
theGluon.E(), theGluon.M());
}
else{
cout<<"Fockstate is not recognised, ending program."<<endl;
exit(1);
}
//
// Generate event. Abort on failure.
//
cout.setstate(std::ios::failbit); // suppress all Pythia print-out
bool ok = mPythia->next();
cout.clear();
if (!ok || mPythiaEvent->size() == 0) {
if (settings->verbose() && settings->verboseLevel() > 2) {
cout << "InclusiveFinalStateGenerator::generate(): Pythia8 event generation aborted prematurely, owing to error!" << endl;
}
return false;
}
// mPythiaEvent->list(); //Print out the pythia event to the screen
//
// Loop over pythia particles and fill our particle list
//
int offset = 3+ifock; //pythia particles start at i=3(QQ), or i=4(QQG) in pythia list
event->numberOfPythiaParticlesStored = 0;
for (int i = offset; i < mPythiaEvent->size(); i++) {
Particle pypart;
pypart.index = event->particles.size();
pypart.pdgId = mPythiaEvent->at(i).id();
if (mPythiaEvent->at(i).status() < 0)
pypart.status = 2;
else
pypart.status = 1;
-
+
pypart.p = TLorentzVector(mPythiaEvent->at(i).px(),
mPythiaEvent->at(i).py(),
mPythiaEvent->at(i).pz(),
mPythiaEvent->at(i).e());
if (mPythiaEvent->at(i).mother1()) pypart.parents.push_back(mPythiaEvent->at(i).mother1()+offset+1);
if (mPythiaEvent->at(i).mother2()) pypart.parents.push_back(mPythiaEvent->at(i).mother2()+offset+1);
if (mPythiaEvent->at(i).daughter1()) pypart.daughters.push_back(mPythiaEvent->at(i).daughter1()+offset+1);
if (mPythiaEvent->at(i).daughter2()) pypart.daughters.push_back(mPythiaEvent->at(i).daughter2()+offset+1);
event->particles.push_back(pypart);
event->numberOfPythiaParticlesStored++;
}
return true;
}
-// Function to calculate the scaling factor lambda numerically
-double InclusiveFinalStateGenerator::computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m, double M) {
+// Function to calculate the scaling factor lambda numerically (from ChatGPT)
+double InclusiveFinalStateGenerator::computeScalingFactor(double px1, double py1, double pz1, double px2, double py2, double pz2, double m_1, double m_2, double M, double lambda_min=0, double lambda_max=100) {
// Define the function to compute the total energy after scaling
auto totalEnergy = [&](double lambda) {
- double E1_prime = std::sqrt(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);
+ double E1_prime = std::sqrt(lambda*lambda*(px1*px1 + py1*py1 + pz1*pz1) + m_1*m_1);
+ double E2_prime = std::sqrt(lambda*lambda*(px2*px2 + py2*py2 + pz2*pz2) + m_2*m_2);
return E1_prime + E2_prime;
};
-
+
// Function to compute the invariant mass condition for a given lambda
auto invariantMassCondition = [&](double lambda) {
double px_sum = lambda * (px1 + px2);
double py_sum = lambda * (py1 + py2);
- double pz_sum = pz1 + pz2;
+ double pz_sum = lambda * (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);
}
+
+// Function that computes f(phi_totg) = (pseudo gluon mass)^2 - Mqq2
+double InclusiveFinalStateGenerator::evaluate(double phi_totg, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus, double mg2=0) {
+ double pt = theXparticle.Pt();
+ double b = 2 * (1 - z) * cos(phi_totg) * pt;
+ double c = (1 - z) * (1 - z) * pt * pt + (1 - z) * Mqq2 - z * (1 - z) * MX2 + z * mg2;
+ double discriminant = 1 - 4 * c / (b * b);
+
+ if (discriminant < 0) {
+ return 1e9; // Large value if not physical
+ }
+
+ double ptg = -b / 2 * (1 - sqrt(discriminant));
+ double phi_tot = acos(theXparticle.Px() / pt);
+
+ double PgMinus = ptg * ptg / PgPlus;
+ double E = (PgPlus + PgMinus) / 2;
+ double pz = (PgPlus - PgMinus) / 2;
+ double px = ptg * cos(phi_tot + phi_totg);
+ double py = ptg * sin(phi_tot + phi_totg);
+
+ TLorentzVector theGluon(px, py, pz, E);
+ TLorentzVector thePseudoGluon = theXparticle - theGluon;
+
+ return thePseudoGluon.M2() - Mqq2;
+}
+
+bool InclusiveFinalStateGenerator::bisection(double phi_min, double phi_max, const TLorentzVector& theXparticle, double z, double Mqq2, double MX2, double PgPlus,
+ double tol, double& phi_root, double mg2=0) {
+
+ double f_min = evaluate(phi_min, theXparticle, z, Mqq2, MX2, PgPlus, mg2);
+ double f_max = evaluate(phi_max, theXparticle, z, Mqq2, MX2, PgPlus, mg2);
+
+ if (f_min * f_max > 0) {
+ // No sign change, cannot use bisection
+ return false;
+ }
+
+ int max_iter = 100;
+ for (int iter = 0; iter < max_iter; ++iter) {
+ double phi_mid = 0.5 * (phi_min + phi_max);
+ double f_mid = evaluate(phi_mid, theXparticle, z, Mqq2, MX2, PgPlus);
+
+ if (std::abs(f_mid) < tol) {
+ phi_root = phi_mid;
+ return true;
+ }
+
+ if (f_min * f_mid < 0) {
+ phi_max = phi_mid;
+ f_max = f_mid;
+ } else {
+ phi_min = phi_mid;
+ f_min = f_mid;
+ }
+ }
+
+ // If we get here, didn't converge
+ return false;
+}
+
+
Index: trunk/src/SartreInclusiveDiffraction.h
===================================================================
--- trunk/src/SartreInclusiveDiffraction.h (revision 549)
+++ trunk/src/SartreInclusiveDiffraction.h (revision 550)
@@ -1,124 +1,129 @@
//==============================================================================
// 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 "InclusiveDiffractiveCrossSectionsFromTables.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;
bool mUseCrossSectionTables;
InclusiveDiffractiveCrossSections *mCrossSection;
InclusiveDiffractionCrossSectionsFromTables *mCrossSectionFromTables;
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*);
+ double bDistribution(double*, double*);
+ double rDistribution_qqT(double*, double*);
+ double rDistribution_qqL(double*, double*);
+ double rDistribution_qqg(double*, double*);
+
};
ostream& operator<<(ostream& os, const TLorentzVector&);
#endif
Index: trunk/src/InclusiveDiffractiveCrossSections.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 549)
+++ trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 550)
@@ -1,722 +1,722 @@
//==============================================================================
// 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();
}
InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){}
double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_total(){return mDsigdbetadQ2dWdz_L_total;}
double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_L_total(){return mDsigdbetadQ2dWdz_T_total;}
double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_qqg(){ return mDsigdbetadQ2dWdz_T_qqg;};
//These three functions are unique for the sub-classes:
double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double, double, double, double){ return 0; }
double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double, double, double, double){ return 0; }
double InclusiveDiffractiveCrossSections::dsigmadbetadz_QQG(double, double, double, double){return 0;}
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){
result+=dsigdbetadQ2dW2dz_total_checked(array[0], array[1], array[2], array[3]);
}
else if(mFockState==QQG or
mFockState==ALL){
result+=dsigdbetadQ2dW2dz_total_qqg_checked(array[0], array[1], array[2], array[3]);
}
else{
cout<<"InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
for(int i=0; i<4; 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::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::crossSectionOfLastCall(){
return mTotalCS;
}
GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){
return mPolarization;
}
unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){
return mQuarkSpecies;
}
double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){
return Settings::quarkMass(mQuarkSpecies);
}
void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
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]);
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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 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_qqg[i];
if(R <= csTi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_T_qqg[i];
}
if(!isChosen) mQuarkSpecies=4;
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
cout << "InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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 << "InclusiveDiffractiveCrossSectionsIntegrals::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<4; i++){
setQuarkIndex(i);
double mf = Settings::quarkMass(mQuarkIndex);
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 (with flux) or nb (without)
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// InclusiveDiffractiveCrossSectionsIntegrals:
// Inheriting class from InclusiveDiffractiveCrossSections
//
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
InclusiveDiffractiveCrossSectionsIntegrals::InclusiveDiffractiveCrossSectionsIntegrals()
{
//Set up the integrals
//Amplitude0:
if(Amp0) delete Amp0;
- Amp0=new TF1("Amp0", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_0, 0, 30., 5);
+ 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-4);
GIAmp0.SetRelTolerance(1e-4);
//Amplitude1:
if(Amp1) delete Amp1;
- Amp1=new TF1("Amp1", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_1, 0, 30. , 5);
+ 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-4);
//AmplitudeQQG:
if(Ampqqg) delete Ampqqg;
- Ampqqg=new TF1("Ampqqg", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_qqg, 0, 30., 4);
+ 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-4);
}
InclusiveDiffractiveCrossSectionsIntegrals::~InclusiveDiffractiveCrossSectionsIntegrals(){
if(WFAmp0) delete WFAmp0;
if(Amp0) delete Amp0;
if(WFAmp1) delete WFAmp1;
if(Amp1) delete Amp1;
if(Ampqqg) delete Ampqqg;
}
double InclusiveDiffractiveCrossSectionsIntegrals::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 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. mf="<<mf<<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 InclusiveDiffractiveCrossSectionsIntegrals::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 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. mf="<<mf<<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 InclusiveDiffractiveCrossSectionsIntegrals::Phi_0(double beta, double xpom, double z, double Q2, double mf){
TF1* Phi0=new TF1("Phi0", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_0, 0, 30., 0);
ROOT::Math::WrappedTF1* WFPhi0=new ROOT::Math::WrappedTF1(*Phi0);
ROOT::Math::GaussIntegrator GIPhi0;
GIPhi0.SetFunction(*WFPhi0);
GIPhi0.SetRelTolerance(1e-4);
Amp0->SetParameter(0, beta);
Amp0->SetParameter(1, xpom);
Amp0->SetParameter(2, z);
Amp0->SetParameter(3, Q2);
Amp0->SetParameter(5, mf);
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double blow=0;
double result=GIPhi0.IntegralUp(blow); //fm6
delete Phi0;
delete WFPhi0;
return result;
}
double InclusiveDiffractiveCrossSectionsIntegrals::Phi_1(double beta, double xpom, double z, double Q2, double mf){
TF1* Phi1=new TF1("Phi1", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_1, 0, 30., 0);
ROOT::Math::WrappedTF1* WFPhi1=new ROOT::Math::WrappedTF1(*Phi1);
ROOT::Math::GaussIntegrator GIPhi1;
GIPhi1.SetFunction(*WFPhi1);
GIPhi1.SetRelTolerance(1e-4);
Amp1->SetParameter(0, beta);
Amp1->SetParameter(1, xpom);
Amp1->SetParameter(2, z);
Amp1->SetParameter(3, Q2);
Amp1->SetParameter(5, mf);
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double blow=0;
double result=GIPhi1.IntegralUp(blow); //fm6
delete Phi1;
delete WFPhi1;
return result;
}
double InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_0(double *var, double*){
double b=*var;
double amplitude=Amplitude_0(b); //fm2
double result=2*M_PI*b*amplitude*amplitude; //fm5
return result;
}
double InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_1(double *var, double*){
double b=*var; //fm
double amplitude=Amplitude_1(b); //fm2
double result=2*M_PI*b*amplitude*amplitude; //fm5
return result;
}
double InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_0(double *var, double *par){
+double InclusiveDiffractiveCrossSections::uiAmplitude_0(double *var, double *par){
double r=*var;
if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
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 InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_1(double *var, double *par){
+double InclusiveDiffractiveCrossSections::uiAmplitude_1(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 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
}
//===================================================
// QQG Calculations
//===================================================
double InclusiveDiffractiveCrossSectionsIntegrals::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.15;
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 InclusiveDiffractiveCrossSectionsIntegrals::Phi_qqg(double xpom, double z, double Q2){
Ampqqg->SetParameter(0, xpom);
Ampqqg->SetParameter(1, z);
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
mQ2=Q2;
ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSectionsIntegrals::uiPhi_qqg, 2);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig.SetFunction(wf);
ig.SetAbsTolerance(0.);
ig.SetRelTolerance(1e-4);
double bRange=3. * dipoleModel()->nucleus()->radius();
double lo[2]={0.,1e-4}; //b, k2
double hi[2]={bRange, Q2}; //b, k2
double result = ig.Integral(lo, hi)/hbarc; //GeV0
return result; //GeV0
}
double InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_qqg(double *var, double *par){
+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/SartreInclusiveDiffraction.cpp
===================================================================
--- trunk/src/SartreInclusiveDiffraction.cpp (revision 549)
+++ trunk/src/SartreInclusiveDiffraction.cpp (revision 550)
@@ -1,875 +1,988 @@
//==============================================================================
// SartreInclusiveDiffraction.cpp
//
// Copyright (C) 2024 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date: 2024-06-03 11:27:05 -0400 (Mon, 03 Jun 2024) $
// $Author: ullrich $
//==============================================================================
//
// Note:
// When not using runcards, the user must first create an instance of
// SartreInclusiveDiffraction and then get the settings via one of:
// SartreInclusiveDiffraction::runSettings()
// EventGeneratorSettings::instance()
// Once init() is called settings cannot be changed any more.
//==============================================================================
#include "Version.h"
#include "Kinematics.h"
#include "Constants.h"
#include "SartreInclusiveDiffraction.h"
#include "Math/Functor.h"
#include "ModeFinderFunctor.h"
#include "Math/BrentMinimizer1D.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/GSLMinimizer.h"
#include "TUnuranMultiContDist.h"
#include <string>
#include <limits>
#include <cmath>
#include <iomanip>
#include "TH2D.h"
#include "TFile.h"
#include "TF3.h"
#include "InclusiveDiffractiveCrossSectionsFromTables.h"
#include "InclusiveDiffractiveCrossSections.h"
#include "DglapEvolution.h"
using namespace std;
#define PR(x) cout << #x << " = " << (x) << endl;
SartreInclusiveDiffraction::SartreInclusiveDiffraction()
{
mSettings = EventGeneratorSettings::instance();
mRandom = mSettings->randomGenerator();
mIsInitialized = false;
mCurrentEvent = 0;
mNucleus = 0;
mUpcNucleus= 0;
// mSettings = 0;
mPDF_Functor = 0;
mPDF = 0;
mEventCounter = 0;
mTriesCounter = 0;
mTotalCrossSection = 0;
mTotalCrossSection_qqg = 0;
mTotalCrossSection_qq = 0;
mCrossSection = 0;
mUnuran = 0;
mEvents = 0;
mTries = 0;
mS = 0;
mA = 0;
}
SartreInclusiveDiffraction::~SartreInclusiveDiffraction()
{
delete mNucleus;
delete mUpcNucleus;
delete mPDF_Functor;
delete mPDF;
delete mCrossSection;
delete mUnuran;
delete mCurrentEvent;
}
bool SartreInclusiveDiffraction::init(const char* runcard)
{
mStartTime = time(0);
bool ok;
//
// Reset member variables.
// Note that one instance of Sartre should be able to get
// initialized multiple times.
//
mEvents = 0;
mTries = 0;
mTotalCrossSection = 0;
mTotalCrossSection_qqg = 0;
mTotalCrossSection_qq = 0;
//
// Print header
//
string ctstr(ctime(&mStartTime));
ctstr.erase(ctstr.size()-1, 1);
cout << "/========================================================================\\" << endl;
cout << "| |" << endl;
cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
cout << "| |" << endl;
cout << "| An event generator for inclusive diffraction |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
cout << "| Copyright (C) 2010-2024 Tobias Toll and Thomas Ullrich |" << endl;
cout << "| |" << endl;
cout << "| This program is free software: you can redistribute it and/or modify |" << endl;
cout << "| it under the terms of the GNU General Public License as published by |" << endl;
cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl;
cout << "| any later version. |" << endl;
cout << "| |" << endl;
cout << "| Code compiled on " << setw(12) << left << __DATE__;
cout << setw(41) << left << __TIME__ << right << '|' << endl;
cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
cout << "\\========================================================================/" << endl;
mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton
mSettings->setInclusiveDiffractionMode(true);
//
// Read runcard if available
//
if (runcard) {
if (!mSettings->readSettingsFromFile(runcard)) {
cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl;
exit(1);
}
else
cout << "Runcard is '" << runcard << "'." << endl;
}
else
cout << "No runcard provided." << endl;
mUseCrossSectionTables=mSettings->useInclusiveCrossSectionTables();
TableGeneratorSettings* tgsettings=TableGeneratorSettings::instance();
tgsettings->setDipoleModelType(mSettings->dipoleModelType());
tgsettings->setDipoleModelParameterSet(mSettings->dipoleModelParameterSet());
// Set up DGLAP evolution for integrals:
if(!mUseCrossSectionTables){
DglapEvolution &dglap = DglapEvolution::instance();
dglap.generateLookupTable(1000, 1000); dglap.useLookupTable(true);
}
//
// Set beam particles and center of mass energy
//
mElectronBeam = mSettings->eBeam();
mHadronBeam = mSettings->hBeam();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mA = mSettings->A();
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (allowBreakup) {
if (!getenv("SARTRE_DIR")) {
cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n"
"to locate tables needed for the generation if nuclear breakups." << endl;
exit(1);
}
}
if (mNucleus) delete mNucleus;
mNucleus = new FrangibleNucleus(mA, allowBreakup);
string upcNucleusName;
cout << "Sartre is running in inclusive diffraction mode" << endl;
cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron beam: " << mHadronBeam << endl;
cout << "Electron beam: " << mElectronBeam << endl;
//
// Get details about the processes and models
//
mDipoleModelType = mSettings->dipoleModelType();
mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
cout << "Process is ";
if (mA > 1)
cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + X"<< endl;
else
cout << "e + p -> e' + p' + X"<< endl;
//
// Print-out seed for reference
//
cout << "Random generator seed: " << mSettings->seed() << endl;
//
// Load in the tables containing the amplitude moments
//
if (!getenv("SARTRE_DIR")) {
cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
exit(1);
}
//
// Kinematic limits and generator range
//
// There are 3 ranges we have to deal with
// 1. the kinematic range requested by the user
// if given.
// The user can only control Q2 and W but not t.
// For UPC that's xpom.
// 2. the range of the table(s)
// 3. the kinematically allowed range
//
// Of course (3) is more complex than a simple cube/square.
// However, we deal with the detailed shape of the kinematic
// range later using Kinematics::valid() when we generate the
// individual events.
// For setting up UNU.RAN we have to get the cubic/square
// envelope that satifies (1)-(3).
// Note, that they are correlated which makes the order
// in which we do things a bit tricky.
//
//
// Step 2:
// Kinematic limits might overrule boundaries from step 1
//
// let:
// 0: beta
// 1: Q2
// 2: W2
// 3: z
//
// Setup cross-section functor
// It is this functor that is used by all other functors,
// functions, and wrappers when dealing with cross-sections.
//
if (mCrossSection) delete mCrossSection;
if(!mUseCrossSectionTables){
mCrossSection = new InclusiveDiffractiveCrossSectionsIntegrals;
}
else{
mCrossSection = new InclusiveDiffractionCrossSectionsFromTables;
}
//
// Here we need to put in Q2 > mu02.
// The calculation goes haywire for Q2 becoming too small
//
double mu02 = mCrossSection->dipoleModel()->getParameters()->mu02();
// double mf = Settings::quarkMass(0);
double mf=tt_quarkMass[0];
//Now for the diffractive mass:
double kineMX2min = 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, 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);
}
-
+ if (mLowerLimit[3] >= mUpperLimit[3]) {
+ cout << "Invalid range in z: z=[" << sqrt(mLowerLimit[3]) << ", " << sqrt(mUpperLimit[3]) << "]." << endl;
+ exit(1);
+ }
+
//
// Print-out limits (all verbose levels)
//
if (true){//}(mSettings->verbose()) {
cout << "Sartre was thrown into the world, restricted by kinematics and user inputs alike:" << endl;
cout << setw(10) << " beta=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl;
cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl;
cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl;
cout << setw(10) << " z=[" << mLowerLimit[3] << ", " << mUpperLimit[3] << "]" << endl;
cout << setw(10) << " MX2=[" << kineMX2min << ", " << kineMX2max << "]" << endl;
cout << setw(10) << " y=[" << kineYmin << ", 1" << "]" << endl;
}
//
// UNU.RAN needs the domain (boundaries) and the mode.
// The domain is already defined, here we find the mode, which is tricky.
// The max. cross-section is clearly at the domain boundary in Q2=Q2min.
// The position in W2 and beta is not obvious. The mode should be at z=0.5.
// The approach here is to use the BrentMinimizer1D that
// performs first a scan a then a Brent fit.
//
mCrossSection->setCheckKinematics(false);
double theMode[4];
//
// Find the mode.
// Assume that the mode is at W2=W2max, Q2=Q2min, z=0.5
// Start with beta=0.5 => MX2=Q2
// It is kinemtically easier to use MX2 instead of beta here.
//
if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
theMode[1] = mLowerLimit[1]; // Q2
theMode[2] = mUpperLimit[2]; // W2
theMode[3] = 0.5; // z
double MX2min=mLowerLimit[1]*(1-mUpperLimit[0])/mUpperLimit[0];
MX2min=max(kineMX2min, MX2min);
double MX2max=mLowerLimit[1]*(1-mLowerLimit[0])/mLowerLimit[0];
MX2max=min(kineMX2max, MX2max);
InclusiveDiffractionModeFinderFunctor modeFunctor(mCrossSection, theMode[1], theMode[2], theMode[3], MX2min, MX2max);
ROOT::Math::BrentMinimizer1D minimizer;
minimizer.SetFunction(modeFunctor, MX2min, MX2max);
minimizer.SetNpx(static_cast<int>(mUpperLimit[0]-mLowerLimit[0]));
ok = minimizer.Minimize(100000, 0, 1.e-8);
if (! ok) {
cout << "Error, failed to find mode of pdf." << endl;
exit(1);
}
theMode[0] = minimizer.XMinimum(); // Mx2
double MX2=theMode[0];
theMode[0]=theMode[1]/(theMode[1]+MX2); //Mx2->beta
double crossSectionAtMode = (*mCrossSection)(MX2, theMode[1], theMode[2], theMode[3]);
if (mSettings->verbose()) {
cout << "\tlocation: beta=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]) << ", z=" << theMode[3];
cout << "; value: " << crossSectionAtMode << endl;
}
mCrossSection->setCheckKinematics(true);
// domain and mode for Q2 -> log(Q2)
mLowerLimit[1] = log(mLowerLimit[1]);
mUpperLimit[1] = log(mUpperLimit[1]);
theMode[1] = log(theMode[1]);
if (mPDF_Functor) delete mPDF_Functor;
if (mPDF) delete mPDF;
mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &InclusiveDiffractiveCrossSections::unuranPDF, 4);
mPDF = new TUnuranMultiContDist(*mPDF_Functor, true); // last arg = pdf in log or not
mPDF->SetDomain(mLowerLimit, mUpperLimit);
mPDF->SetMode(theMode);
if (mUnuran) delete mUnuran;
mUnuran = new TUnuran;
mCrossSection->setCheckKinematics(false); // avoid numeric glitch in Init()
mUnuran->Init(*mPDF, "method=hitro");
mCrossSection->setCheckKinematics(true);
mUnuran->SetSeed(mSettings->seed());
//
// Burn in generator
//
double xrandom[4];
for (int i=0; i<100; i++) {
mUnuran->SampleMulti(xrandom);
// cout<<"QQ: beta="<<xrandom[0];
// cout<<" Q2="<<exp(xrandom[1]);
// cout<<" W="<<sqrt(xrandom[2]);
// cout<<" z="<<xrandom[3]<<endl;
}
//
// Set up QQG unu.ran:
//
//The mode is simpler for QQG, as it is at minimum beta
//The limits are the same, except for z which is a different quantity in QQG.
mLowerLimit_qqg[0]=mLowerLimit[0]; //beta
mLowerLimit_qqg[1]=mLowerLimit[1]; //log(Q2)
mLowerLimit_qqg[2]=mLowerLimit[2]; //W2
//z>beta*(1+4*mf*mf/Q2);
mLowerLimit_qqg[3]=mLowerLimit[0] * (1+4*mf*mf/exp(mUpperLimit[1]))+1e-10;
mUpperLimit_qqg[0]=mUpperLimit[0]; //beta
mUpperLimit_qqg[1]=mUpperLimit[1]; //log(Q2)
mUpperLimit_qqg[2]=mUpperLimit[2]; //W2
mUpperLimit_qqg[3]=mUpperLimit[3]; //z
double theMode_qqg[4];
theMode_qqg[0]=mLowerLimit_qqg[0]; //beta
theMode_qqg[1]=mLowerLimit_qqg[1]; //log(Q2)
theMode_qqg[2]=mUpperLimit_qqg[2]; //W2
theMode_qqg[3]=theMode_qqg[0]*(1+4*mf*mf/exp(theMode_qqg[1]))+1e-10; //z
if (mPDF_Functor_qqg) delete mPDF_Functor_qqg;
if (mPDF_qqg) delete mPDF_qqg;
mPDF_Functor_qqg = new ROOT::Math::Functor(mCrossSection, &InclusiveDiffractiveCrossSections::unuranPDF_qqg, 4);
mPDF_qqg = new TUnuranMultiContDist(*mPDF_Functor_qqg, true); // last arg = pdf in log or not
//Rescale the domain in z (and fix it in the pdf)
mPDF_qqg->SetDomain(mLowerLimit_qqg, mUpperLimit_qqg);
mPDF_qqg->SetMode(theMode_qqg);
mCrossSection->setCheckKinematics(false);
if (mUnuran_qqg) delete mUnuran_qqg;
mUnuran_qqg = new TUnuran;
mUnuran_qqg->Init(*mPDF_qqg, "method=hitro");
mUnuran_qqg->SetSeed(mSettings->seed());
// mCrossSection->setCheckKinematics(true);
//
// Burn in QQG generator
//
for (int i=0; i<100; i++) {
mUnuran_qqg->SampleMulti(xrandom);
// cout<<"QQG: beta="<<xrandom[0];
// cout<<" Q2="<<exp(xrandom[1]);
// cout<<" W="<<sqrt(xrandom[2]);
// cout<<" z="<<xrandom[3]<<endl;
}
mEventCounter = 0;
mTriesCounter = 0;
mIsInitialized = true;
cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;
//
// Calculate total cross section in given
// kinematic limits. This will be used to
// decide which type of event to generate,
// q-qbar or q-qbar-g
//
cout<<"Calculating cross sections over the given kinematic range..."<<endl;
cout<<"(This may take while, use your radical freedom and do other things.)"<<endl;
mCrossSection->setFockState(QQ);
mTotalCrossSection=0;
mTotalCrossSection_qq=7.4543e+08; //totalCrossSection();
cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
mCrossSection->setFockState(QQG);
mTotalCrossSection=0;
mTotalCrossSection_qqg=3.79917e+08; //totalCrossSection();
cout<<"The QQG cross section is: "<<mTotalCrossSection_qqg<<" nb. (2/2)"<<endl;
cout<<"Total Cross Section: "<<mTotalCrossSection_qqg+mTotalCrossSection_qq<<" nb"<<endl;
mTotalCrossSection=mTotalCrossSection_qqg+mTotalCrossSection_qq;
return true;
}
double SartreInclusiveDiffraction::cross_section_wrapper(double* xx, double* par){
double arg[4]={xx[0], xx[1], xx[2], par[0]};
double result=mCrossSection->unuranPDF(arg);
return -result;
}
bool SartreInclusiveDiffraction::init(const string& str) // overloaded version of init()
{
if (str.empty())
return init();
else
return init(str.c_str());
}
vector<pair<double,double> > SartreInclusiveDiffraction::kinematicLimits()
{
vector<pair<double,double> > array;
array.push_back(make_pair(mLowerLimit[0], mUpperLimit[0])); // t
array.push_back(make_pair(exp(mLowerLimit[1]), exp(mUpperLimit[1]))); // Q2 or xpom
if (!mSettings->UPC())
array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W
return array;
}
Event* SartreInclusiveDiffraction::generateEvent()
{
if (!mIsInitialized) {
cout << "SartreInclusiveDiffraction::generateEvent(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return 0;
}
//
// Decide which Fock-state to generate:
//
bool isQQ=false;
if(mRandom->Uniform(mTotalCrossSection) <= mTotalCrossSection_qq)
isQQ=true;
//
// Generate one event
//
double xrandom[4];
while (true) {
mTriesCounter++;
delete mCurrentEvent;
mCurrentEvent = new Event;
//
// xrandom[i]
// i=0: beta
// i=1: Q2
// i=2: W2
// i=3: z
//
if(isQQ)
mUnuran->SampleMulti(xrandom);
else
mUnuran_qqg->SampleMulti(xrandom);
xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2
double MX2=xrandom[1]*(1-xrandom[0])/xrandom[0];
bool isValidEvent;
isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2], xrandom[3], sqrt(MX2), true, (mSettings->verboseLevel() > 1));
if (!isValidEvent) {
if (mSettings->verboseLevel() > 2)
cout << "SartreInclusiveDiffraction::generateEvent(): event rejected, not within kinematic limits" << endl;
continue;
}
//
// Fill beam particles in Event structure
// Kinematics for eA is reported as 'per nucleon'
//
mCurrentEvent->eventNumber = mEventCounter;
mCurrentEvent->beta = xrandom[0]; // beta
mCurrentEvent->Q2 = xrandom[1]; // Q2
mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
mCurrentEvent->xpom=mCurrentEvent->x/xrandom[0]; //xpom=x/beta
mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
mCurrentEvent->s = mS; // s
mCurrentEvent->W = sqrt(xrandom[2]);
mCurrentEvent->z = xrandom[3];
mCurrentEvent->MX = sqrt(MX2);
if(isQQ)
mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
else
mCurrentEvent->polarization = transverse;
// mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
mCurrentEvent->quarkSpecies = mCrossSection->quarkSpeciesOfLastCall();
if(isQQ)
mCurrentEvent->crossSectionRatioLT=mCrossSection->crossSectionRatioLTOfLastCall();
else
mCurrentEvent->crossSectionRatioLT=0;
//
- // The following needs to be modified for QQG
- //
+ // For generating final state azimuthal decorelation
+ // we need to calculate Omega(b, r).
+ // For this we need to generate b and r values.
+ //
+ // Start with b (written by Bhakta Naik and TT):
+ // From the thickness function (a leading twist approximation)
+ //
+ double bmin = 0;
+ double bmax = 3.5 * mCrossSection->dipoleModel()->nucleus()->radius();
+
+ TF1* fbDistribution = new TF1("fbDistribution",this ,&SartreInclusiveDiffraction::bDistribution, bmin, bmax, 0);
+
+ double bGenerated = fbDistribution->GetRandom()*hbarc; // fm
+ //
+ // Now generate r
+ //
+ double mf = mCrossSection->quarkMassOfLastCall();
+ double rmin=1e-3;
+ double rmax=3.5 * mCrossSection->dipoleModel()->nucleus()->radius();
+ double rGenerated=0;
+ if(isQQ)
+ {
+ TF1* frDistribution_qq;
+ if(mCurrentEvent->polarization == transverse )
+ {
+ frDistribution_qq = new TF1("frDistribution_qq", this, &SartreInclusiveDiffraction::rDistribution_qqT, rmin, rmax, 6);
+ }
+ else
+ {
+ frDistribution_qq = new TF1("frDistribution_qq", this, &SartreInclusiveDiffraction::rDistribution_qqL, rmin, rmax, 6);
+ }
+ frDistribution_qq->SetParameters(mCurrentEvent->beta, mCurrentEvent->xpom, mCurrentEvent->z, mCurrentEvent->Q2, bGenerated, mf);
+ rGenerated = frDistribution_qq->GetRandom();
+
+ delete frDistribution_qq;
+ }
+ else{ //QQG
+ // Create a TF2 object.
+ // The range for 'r' is from 0 to max_r,
+ // and for 'kappa' is from 0 to Q2.
+ // The parameter array for TF2 is defined as:
+ // {Q2, xpom, z, b_rand_value}
+ //
+ TF2* frDistribution_qqg = new TF2("frDistribution_qqg", this,&SartreInclusiveDiffraction::rDistribution_qqg, rmin, rmax, 0, mCurrentEvent->Q2, 4);
+
+ frDistribution_qqg->SetParameters(mCurrentEvent->xpom, mCurrentEvent->z, bGenerated, sqrt(mCurrentEvent->Q2));
+
+ double kGenerated;
+ frDistribution_qqg->GetRandom2(rGenerated, kGenerated);
+ delete frDistribution_qqg;
+ }
+ //Calculate Omega:
double dsigmadb2=0;
- 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);
+ dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(rGenerated, bGenerated, mCurrentEvent->xpom);
}
else {
- dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(r_average*hbarc, 0, mCurrentEvent->xpom);
+ dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(rGenerated, bGenerated, mCurrentEvent->xpom);
}
mCurrentEvent->Omega=-2.*log(1-dsigmadb2/2.);
-
+ // Calculate the saturation scale following
+ // H. Kowalski, T. Lappi and R. Venugopalan, Phys. Rev. Lett. 100 (2008) 022303 [arXiv:0705.3047 [hep-ph]]
+ double rS;
+ for(rS=0.1; rS<20; rS+=0.001){
+ double dsigmadb2=0;
+ if (mA == 1) {
+ dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(rS, bGenerated, mCurrentEvent->xpom);
+ }
+ else {
+ dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(rS, bGenerated, mCurrentEvent->xpom);
+ }
+ if(dsigmadb2 >= 0.44239843) break;
+ }
+ mCurrentEvent->QS=1./rS*hbarc; //GeV
+ double QS=1./rS*hbarc;
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = 11; // e-
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
//
// Generate the final state particles
//
bool ok;
if (isQQ) {
ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQ);
}
else {
ok = mFinalStateGenerator.generate(mA, mCurrentEvent, QQG);
}
if (!ok) {
if (mSettings->verboseLevel() > 1) cout << "SartreInclusiveDiffraction::generateEvent(): failed to generate final state" << endl;
continue;
}
break;
}
mEventCounter++;
//
// Nuclear breakup
//
// If the event is incoherent the final state generator does produce a
// 'virtual' proton with m > m_p which is used in Nucleus to calculate
// the excitation energy and the boost.
//
int indexOfScatteredHadron = 6;
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (mNucleus) mNucleus->resetBreakup(); // clear previous event in any case
if (allowBreakup && mCurrentEvent->diffractiveMode == incoherent && mNucleus) {
int nFragments = mNucleus->breakup(mCurrentEvent->particles[indexOfScatteredHadron].p);
//
// Merge the list of products into the event list.
// We loose some information here. The user can always go back to
// the nucleus and check the decay products for more details.
// In the original list the energy is per nuclei, here we transform it
// to per nucleon to stay consistent with Sartre conventions.
//
const vector<BreakupProduct>& products = mNucleus->breakupProducts();
for (int i=0; i<nFragments; i++) {
Particle fragment;
fragment.index = mCurrentEvent->particles.size();
fragment.pdgId = products[i].pdgId;
fragment.status = 1;
fragment.p = products[i].p*(1/static_cast<double>(products[i].A));
fragment.parents.push_back(indexOfScatteredHadron);
mCurrentEvent->particles.push_back(fragment);
}
}
return mCurrentEvent;
}
double SartreInclusiveDiffraction::totalCrossSection()
{
if (mTotalCrossSection == 0) {
//
// Limits of integration in beta, Q2, W2, z
//
double xmin[4];
double xmax[4];
copy(mLowerLimit, mLowerLimit+4, xmin);
copy(mUpperLimit, mUpperLimit+4, xmax);
//
// At this point mLowerLimit[1] and mUpperLimit[1]
// are in log(Q2)
//
xmin[1] = exp(xmin[1]); // log Q2 limit
xmax[1] = exp(xmax[1]); // log Q2 limit
mTotalCrossSection = calculateTotalCrossSection(xmin, xmax);
}
return mTotalCrossSection;
}
double SartreInclusiveDiffraction::totalCrossSection(double lower[4], double upper[4]) // beta, Q2, W, z
{
lower[2] *= lower[2]; upper[2] *= upper[2]; // W -> W2
double result = calculateTotalCrossSection(lower, upper);
return result;
}
double SartreInclusiveDiffraction::integrationVolume(){
//
// Limits of integration in beta, Q2, W2, z
//
double xmin[4];
double xmax[4];
copy(mLowerLimit, mLowerLimit+4, xmin);
copy(mUpperLimit, mUpperLimit+4, xmax);
//
// At this point mLowerLimit[1] and mUpperLimit[1]
// are in log(Q2)
//
xmin[1] = exp(xmin[1]); // log Q2 limit
xmax[1] = exp(xmax[1]); // log Q2 limit
double volume = (xmax[0]-xmin[0]) * (xmax[1]-xmin[1]) * (xmax[2]-xmin[2]) * (xmax[3]-xmin[3]);
return volume; //GeV4
}
EventGeneratorSettings* SartreInclusiveDiffraction::runSettings()
{
return EventGeneratorSettings::instance();
}
double SartreInclusiveDiffraction::crossSectionsClassWrapper(const double* var){
return (*mCrossSection)(var);
}
double SartreInclusiveDiffraction::calculateTotalCrossSection(double lower[4], double upper[4])
{
double result = 0;
if (!mIsInitialized) {
cout << "SartreInclusiveDiffraction::calculateTotalCrossSection(): Error, Sartre is not initialized yet." << endl;
cout << " Call init() before trying to generate events." << endl;
return result;
}
//
// Calculate integral using adaptive numerical method
//
// Options: ADAPTIVE, kVEGAS, kPLAIN, kMISER
// no abs tolerance given -> relative only
const double precision = 5e-4;
// unsigned int numberofcalls = 1000000;
unsigned int numberofcalls = 0;//1e5;
ROOT::Math::Functor wfL(this, &SartreInclusiveDiffraction::crossSectionsClassWrapper, 4);
ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE, 0, precision, numberofcalls);
ig.SetFunction(wfL);
result = ig.Integral(lower, upper);
//
// If it fails we switch to a MC integration which is usually more robust
// although not as accurate. This should happen very rarely if at all.
//
if (result <= numeric_limits<float>::epsilon()) {
cout << "SartreInclusiveDiffraction::calculateTotalCrossSection(): warning, adaptive integration failed - switching to VEGAS method." << endl;
ROOT::Math::IntegratorMultiDim igAlt(ROOT::Math::IntegrationMultiDim::kVEGAS);
igAlt.SetFunction(wfL);
igAlt.SetRelTolerance(precision);
igAlt.SetAbsTolerance(0);
result = igAlt.Integral(lower, upper);
}
return result;
}
+double SartreInclusiveDiffraction::bDistribution(double* var, double*)
+{
+ double b = *var;
+
+ double result=0;
+ if(mA==1){
+ double BG = mCrossSection->dipoleModel()->getParameters()->BG();
+ double arg = b*b/(2*BG);
+ arg /= hbarc2;
+ result = 1/(2*M_PI*BG) * exp(-arg);
+ }
+ else{
+ result = mNucleus->T(b);
+ }
+ return result;
+}
+// Three functions for randomly choosing a value of r:
+double SartreInclusiveDiffraction::rDistribution_qqT(double* var, double* par)
+{
+ double result = mCrossSection->uiAmplitude_1(var, par);
+ if(isnan(result)) return 0;
+ return result * result;
+}
+
+double SartreInclusiveDiffraction::rDistribution_qqL(double* var, double* par)
+{
+ double result = mCrossSection->uiAmplitude_0(var, par);
+
+ if(isnan(result)) return 0;
+ return result * result;
+}
+double SartreInclusiveDiffraction::rDistribution_qqg(double* var, double* par)
+{
+ double r = var[0]; //fm
+ double k = var[1]; //GeV
+ if(r>100.) return 0; //IntegralUp sometimes tries too large values that break the dipole model
+ double xpom=par[0];
+ double z=par[1];
+
+ double b=par[2];
+ double Q2 = par[3];
+
+ double para[4] = {xpom, z, k , b};
+ double result = mCrossSection->uiAmplitude_qqg(&r, para);
+ return pow(k, 4) * log(Q2 / (k * k)) * result * result ;
+}
+
+
const FrangibleNucleus* SartreInclusiveDiffraction::nucleus() const {return mNucleus;}
void SartreInclusiveDiffraction::listStatus(ostream& os) const
{
os << "Event summary: " << mEventCounter<< " events generated, " << mTriesCounter << " tried" << endl;
time_t delta = runTime();
os << "Total time used: " << delta/60 << " min " << delta - 60*(delta/60) << " sec" << endl;
}
time_t SartreInclusiveDiffraction::runTime() const
{
time_t now = time(0);
return now-mStartTime;
}
//==============================================================================
//
// Utility functions and operators (helpers)
//
//==============================================================================
ostream& operator<<(ostream& os, const TLorentzVector& v)
{
os << v.Px() << '\t' << v.Py() << '\t' << v.Pz() << '\t' << v.E() << '\t';
double m2 = v*v;
if (m2 < 0)
os << '(' << -sqrt(-m2) << ')';
else
os << '(' << sqrt(m2) << ')';
return os;
}

File Metadata

Mime Type
application/octet-stream
Expires
Fri, May 2, 4:31 PM (2 d)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4978729
Default Alt Text
(149 KB)

Event Timeline