Page MenuHomeHEPForge
Files F8723635

No OneTemporary
Actions

View Options

Index: trunk/src/InclusiveDiffractiveCrossSections.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 537)
+++ trunk/src/InclusiveDiffractiveCrossSections.cpp (revision 538)
@@ -1,688 +1,719 @@
//==============================================================================
// 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()
+InclusiveDiffractiveCrossSections::InclusiveDiffractiveCrossSections(){}
+
+InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){}
+
+double InclusiveDiffractiveCrossSections::operator()(double, double, double, double){return 0;}
+double InclusiveDiffractiveCrossSections::operator()(const double*){ return 0;}
+DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel(){return 0;}
+double InclusiveDiffractiveCrossSections::unuranPDF(const double*){return 0;}
+double InclusiveDiffractiveCrossSections::unuranPDF_qqg(const double*){return 0;}
+void InclusiveDiffractiveCrossSections::setCheckKinematics(bool){}
+void InclusiveDiffractiveCrossSections::setFockState(FockState){}
+GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){return transverse;}
+double InclusiveDiffractiveCrossSections::crossSectionOfLastCall(){return 0;}
+unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){return 0;}
+double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){return 0;}
+double InclusiveDiffractiveCrossSections::crossSectionRatioLTOfLastCall() const{return 0;}
+void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int){}
+double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double, double, double, double){ return 0; }
+double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double, double, double, double){ return 0; }
+double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(double, double, double, double){return 0;}
+double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(double, double, double, double){return 0;}
+
+double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double, double, double, double, GammaPolarization){return 0;}
+
+double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_qqg(double, double, double, double){return 0;}
+
+double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_total(){return 0;}
+double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_L_total(){return 0;}
+double* InclusiveDiffractiveCrossSections::dsigdbetadQ2dWdz_T_qqg(){return 0;}
+
+
+InclusiveDiffractiveCrossSectionsIntegrals::InclusiveDiffractiveCrossSectionsIntegrals()
{
mSettings = EventGeneratorSettings::instance();
mRandom = mSettings->randomGenerator();
mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
mPhotonFlux.setS(mS);
mCheckKinematics = true;
mCrossSectionRatioLT = 0;
mQ2=0;
// TableGeneratorSettings* settings = TableGeneratorSettings::instance();
if(mDipoleModel) delete mDipoleModel;
mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
DipoleModelType model=mSettings->dipoleModelType();
if (model==bSat) {
mDipoleModel = new DipoleModel_bSat(mSettings);
}
else if(model==bNonSat){
mDipoleModel = new DipoleModel_bNonSat(mSettings);
}
else {
cout << "Integrals::init(): Error, model not implemented: "<< model << endl;
exit(1);
}
mA=mSettings->A();
//Set up the integrals
//Amplitude0:
if(Amp0) delete Amp0;
- Amp0=new TF1("Amp0", this, &InclusiveDiffractiveCrossSections::uiAmplitude_0, 0, 30., 5);
+ Amp0=new TF1("Amp0", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_0, 0, 30., 5);
if(WFAmp0) delete WFAmp0;
WFAmp0=new ROOT::Math::WrappedTF1(*Amp0);
GIAmp0.SetFunction(*WFAmp0);
GIAmp0.SetRelTolerance(1e-4);
//Amplitude1:
if(Amp1) delete Amp1;
- Amp1=new TF1("Amp1", this, &InclusiveDiffractiveCrossSections::uiAmplitude_1, 0, 30. , 5);
+ Amp1=new TF1("Amp1", this, &InclusiveDiffractiveCrossSectionsIntegrals::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, &InclusiveDiffractiveCrossSections::uiAmplitude_qqg, 0, 30., 4);
+ Ampqqg=new TF1("Ampqqg", this, &InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_qqg, 0, 30., 4);
if(WFAmpqqg) delete WFAmpqqg;
WFAmpqqg=new ROOT::Math::WrappedTF1(*Ampqqg);
GIAmpqqg.SetFunction(*WFAmpqqg);
GIAmpqqg.SetRelTolerance(1e-4);
}
-InclusiveDiffractiveCrossSections::~InclusiveDiffractiveCrossSections(){
+InclusiveDiffractiveCrossSectionsIntegrals::~InclusiveDiffractiveCrossSectionsIntegrals(){
if(WFAmp0) delete WFAmp0;
if(Amp0) delete Amp0;
if(WFAmp1) delete WFAmp1;
if(Amp1) delete Amp1;
if(Ampqqg) delete Ampqqg;
}
-void InclusiveDiffractiveCrossSections::setCheckKinematics(bool val) {mCheckKinematics = val;}
+void InclusiveDiffractiveCrossSectionsIntegrals::setCheckKinematics(bool val) {mCheckKinematics = val;}
-void InclusiveDiffractiveCrossSections::setFockState(FockState val){
+void InclusiveDiffractiveCrossSectionsIntegrals::setFockState(FockState val){
mFockState=val;
}
-FockState InclusiveDiffractiveCrossSections::getFockState(){
+FockState InclusiveDiffractiveCrossSectionsIntegrals::getFockState(){
return mFockState;
}
-double InclusiveDiffractiveCrossSections::operator()(double MX2, double Q2, double W2, double z){
+double InclusiveDiffractiveCrossSectionsIntegrals::operator()(double MX2, double Q2, double W2, double z){
return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
}
-double InclusiveDiffractiveCrossSections::operator()(const double* array){
+double InclusiveDiffractiveCrossSectionsIntegrals::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<<"InclusiveDiffractiveCrossSections::operator(): Fockstate is invalid, stopping"<<endl;
+ cout<<"InclusiveDiffractiveCrossSectionsIntegrals::operator(): Fockstate is invalid, stopping"<<endl;
exit(0);
}
return result;
}
-double InclusiveDiffractiveCrossSections::unuranPDF(const double* array){
+double InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::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 InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z){
double result = 0;
//
// Check if in kinematically allowed region
double MX2=Kinematics::MX2(Q2, beta, 0);
if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
- cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
+ 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 xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double csT = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, transverse);
double csL = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, longitudinal);
result = csT + csL;
mCrossSectionRatioLT=csL/csT;
//
// Polarization
//
if (mRandom->Uniform(result) <= csT)
mPolarization = transverse;
else
mPolarization = longitudinal;
//
// Quark Species
//
//
if(mPolarization == transverse){
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(csT-cs_rejected);
double csTi=mDsigdbetadQ2dWdz_T_total[i];
if(R <= csTi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_T_total[i];
}
if(!isChosen) mQuarkSpecies=4;
}
else{
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(csL-cs_rejected);
double csLi=mDsigdbetadQ2dWdz_L_total[i];
if(R <= csLi){
mQuarkSpecies=i;
isChosen=true;
break;
}
cs_rejected+=mDsigdbetadQ2dWdz_L_total[i];
}
if(!isChosen) mQuarkSpecies=4;
}
// Print-out at high verbose levels
//
if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
- cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): " << result;
+ cout << "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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = NaN at" << endl;
+ 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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = inf at" << endl;
+ 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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, negative cross-section at" << endl;
+ 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::dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
+double InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
if(pol==transverse){
result = dsigmadbetadz_T(beta, Q2, W2, z);
}
else if(pol==longitudinal){
result = dsigmadbetadz_L(beta, Q2, W2, z);
}
return result; //nb
}
-double InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
+double InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization pol) {
double result = 0;
for(int i=0; i<5; i++){
setQuarkIndex(i);
double tmpresult = dsigdbetadz_total(beta, Q2, W2, z, pol); //nb
if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,pol); //nb/GeV4
if(pol == transverse)
mDsigdbetadQ2dWdz_T_total[i]=tmpresult;
if(pol == longitudinal)
mDsigdbetadQ2dWdz_L_total[i]=tmpresult;
result+=tmpresult;
}
return result; //nb/GeV4
}
-double InclusiveDiffractiveCrossSections::crossSectionOfLastCall(){
+double InclusiveDiffractiveCrossSectionsIntegrals::crossSectionOfLastCall(){
return mTotalCS;
}
-GammaPolarization InclusiveDiffractiveCrossSections::polarizationOfLastCall(){
+GammaPolarization InclusiveDiffractiveCrossSectionsIntegrals::polarizationOfLastCall(){
return mPolarization;
}
-unsigned int InclusiveDiffractiveCrossSections::quarkSpeciesOfLastCall(){
+unsigned int InclusiveDiffractiveCrossSectionsIntegrals::quarkSpeciesOfLastCall(){
return mQuarkSpecies;
}
-double InclusiveDiffractiveCrossSections::quarkMassOfLastCall(){
+double InclusiveDiffractiveCrossSectionsIntegrals::quarkMassOfLastCall(){
// return Settings::quarkMass(mQuarkSpecies); //#TT tmp
- tt_quarkMass[mQuarkSpecies];
+ return tt_quarkMass[mQuarkSpecies];
}
-void InclusiveDiffractiveCrossSections::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
+void InclusiveDiffractiveCrossSectionsIntegrals::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
-double InclusiveDiffractiveCrossSections::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {
+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 mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
double sqrtArg = 1.-4.*mf*mf/MX2;
if (sqrtArg<0){
if (mSettings->verboseLevel() > 2)
cout<<"There is no phase-space for z, return 0."<<endl;
return 0;
}
double z0 = (1.-sqrt(sqrtArg))/2.;
if (z<z0 or z>1-z0){
if (mSettings->verboseLevel() > 2)
cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
return 0;
}
double pt2 = z*(1-z)*Q2-mf*mf;
if (pt2<0){
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
return 0;
}
if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark production, return 0."<<endl;
return 0;
}
double epsilon2=z*(1-z)*Q2+mf*mf;
double ef=quarkCharge[mQuarkIndex];
double Phi0=Phi_0(beta, xpom, z, Q2, mf);
double Phi1=Phi_1(beta, xpom, z, Q2, mf);
double prefactor=Nc*Q2*alpha_em/(4.*M_PI*beta*beta)*ef*ef*z*(1-z); //GeV2
prefactor/=2.; //expanded z-range
double term1=epsilon2*(z*z+(1-z)*(1-z))*Phi1; //GeV2*fm6
double term2=mf*mf*Phi0; //GeV2*fm6
double result=prefactor*(term1+term2); //GeV4*fm6
result/=hbarc2*hbarc2; //fm2
result *= 1e7; //nb
return result;//nb
}
-double InclusiveDiffractiveCrossSections::dsigmadbetadz_L(double beta, double Q2, double W2, double z){
+double 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 mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
double pt2 = z*(1-z)*Q2-mf*mf;
if (pt2<0){
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
return 0;
}
double sqrtArg=1.-4.*mf*mf/MX2;
if (sqrtArg<0){
if (mSettings->verboseLevel() > 2)
cout<<"There is no phase-space for z, return 0."<<endl;
return 0;
}
double z0=(1.-sqrt(sqrtArg))/2.;
if (z<z0 or z>1-z0){
if (mSettings->verboseLevel() > 2)
cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
return 0;
}
if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
if (mSettings->verboseLevel() > 2)
cout<<"Not enough phase-space for quark production, return 0."<<endl;
return 0;
}
double ef=quarkCharge[mQuarkIndex];
double Phi0=Phi_0(beta, xpom, z, Q2, mf); //fm6
double term=Nc*Q2*Q2*alpha_em/(M_PI*beta*beta); //GeV4
term*=ef*ef*z*z*z*(1-z)*(1-z)*(1-z);
term*=Phi0; //fm6*GeV4
term/=2.; //expanded z-range
double result=term;
result/=hbarc2*hbarc2; //fm2
result *= 1e7; //nb
return result; //nb
}
-double InclusiveDiffractiveCrossSections::Phi_0(double beta, double xpom, double z, double Q2, double mf){
- TF1* Phi0=new TF1("Phi0", this, &InclusiveDiffractiveCrossSections::uiPhi_0, 0, 30., 0);
+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);
double blow=0;
double result=GIPhi0.IntegralUp(blow); //fm6
delete Phi0;
delete WFPhi0;
return result;
}
-double InclusiveDiffractiveCrossSections::Phi_1(double beta, double xpom, double z, double Q2, double mf){
- TF1* Phi1=new TF1("Phi1", this, &InclusiveDiffractiveCrossSections::uiPhi_1, 0, 30., 0);
+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);
double blow=0;
double result=GIPhi1.IntegralUp(blow); //fm6
delete Phi1;
delete WFPhi1;
return result;
}
-double InclusiveDiffractiveCrossSections::uiPhi_0(double *var, double *par){
+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 InclusiveDiffractiveCrossSections::uiPhi_1(double *var, double *par){
+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 InclusiveDiffractiveCrossSections::Amplitude_0(double b){
+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 InclusiveDiffractiveCrossSections::Amplitude_1(double b){
+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 InclusiveDiffractiveCrossSections::uiAmplitude_0(double *var, double *par){
+double InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_0(double *var, double *par){
double r=*var;
double beta=par[0];
double xpom=par[1];
double z=par[2];
double Q2=par[3];
double b=par[4];
double mf=par[5];
double MX2=Q2*(1-beta)/beta;
double kappa2=z*(1-z)*MX2-mf*mf;
double kappa=sqrt(kappa2);
double epsilon2=z*(1-z)*Q2+mf*mf;
double epsilon=sqrt(epsilon2); //GeV
double besselK0=TMath::BesselK0(epsilon*r/hbarc);
double besselJ0=TMath::BesselJ0(kappa*r/hbarc);
double dsigmadb2=0;
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
return r*besselK0*besselJ0*dsigmadb2; //fm
}
-double InclusiveDiffractiveCrossSections::uiAmplitude_1(double *var, double *par){
+double InclusiveDiffractiveCrossSectionsIntegrals::uiAmplitude_1(double *var, double *par){
double r=*var; //fm
double beta=par[0];
double xpom=par[1];
double z=par[2];
double Q2=par[3]; //GeV2
double b=par[4];
double mf=par[5];
double MX2=Q2*(1-beta)/beta; //GeV2
double kappa2=z*(1-z)*MX2-mf*mf; //GeV2
double kappa=sqrt(kappa2); //GeV
double epsilon2=z*(1-z)*Q2+mf*mf;
double epsilon=sqrt(epsilon2); //GeV
double besselK1=TMath::BesselK1(epsilon*r/hbarc);
double besselJ1=TMath::BesselJ1(kappa*r/hbarc);
double dsigmadb2=0; //GeV0
if(mA==1)
dsigmadb2=mDipoleModel->dsigmadb2ep(r, b, xpom);
else
dsigmadb2=mDipoleModel->coherentDsigmadb2(r, b, xpom);
return r*besselK1*besselJ1*dsigmadb2; //fm
}
-DipoleModel* InclusiveDiffractiveCrossSections::dipoleModel(){
+DipoleModel* InclusiveDiffractiveCrossSectionsIntegrals::dipoleModel(){
return mDipoleModel;
}
-double InclusiveDiffractiveCrossSections::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
+double InclusiveDiffractiveCrossSectionsIntegrals::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
//===================================================
// QQG Calculations
//===================================================
-double InclusiveDiffractiveCrossSections::unuranPDF_qqg(const double* array){
+double InclusiveDiffractiveCrossSectionsIntegrals::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){
+double InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z){
//
// Check if in kinematically allowed region
//
double MX2=Kinematics::MX2(Q2, beta, 0);
if(mCheckKinematics && z<beta) {
if (mSettings->verboseLevel() > 2)
cout<<"z<beta, beta="<<beta<<", z="<<z<<endl;
return 0; // For QQG beta < z < 1
}
if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
if (mSettings->verboseLevel() > 2)
- cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
+ 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 xpom=Kinematics::xpomeron(0, Q2, W2, sqrt(MX2));
if(mA>1)
dipoleModel()->createSigma_ep_LookupTable(xpom);
double result = dsigdbetadQ2dW2dz_qqg(beta, Q2, W2, z);
//
// Quark Species
//
//
bool isChosen=false;
double cs_rejected=0;
for(int i=0; i<4; i++){
double R=mRandom->Uniform(result-cs_rejected);
double csTi=mDsigdbetadQ2dWdz_T_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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): " << result;
+ 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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = NaN at" << endl;
+ 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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = inf at" << endl;
+ 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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, negative cross-section at" << endl;
+ 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 InclusiveDiffractiveCrossSectionsIntegrals::dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z) {
double result = 0;
for(int i=0; i<5; i++){
setQuarkIndex(i);
// double mf = Settings::quarkMass(mQuarkIndex);
double mf = tt_quarkMass[i]; //#TT tmp
double Mqq2=(z/beta-1)*Q2;
if(Mqq2<4*mf*mf) continue; //Not enough phase-space for the q and qbar
double tmpresult = dsigmadbetadz_QQG(beta, Q2, W2, z); //nb
if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,transverse); //nb/GeV4
mDsigdbetadQ2dWdz_T_qqg[i]=tmpresult;
result+=tmpresult;
}
return result; //nb/GeV4
}
-double InclusiveDiffractiveCrossSections::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
+double 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.14;
double ef=quarkCharge[mQuarkIndex];
double betafactor=(1-beta/ztilde)*
(1-beta/ztilde)+beta/ztilde*beta/ztilde;
double result=alpha_S*alpha_em/(2*M_PI*M_PI*Q2)*betafactor*ef*ef*Phi; //GeV-2
result*=hbarc2; //fm2
result *= 1e7; //nb
return result; //nb
}
-double InclusiveDiffractiveCrossSections::Phi_qqg(double xpom, double z, double Q2){
+double InclusiveDiffractiveCrossSectionsIntegrals::Phi_qqg(double xpom, double z, double Q2){
Ampqqg->SetParameter(0, xpom);
Ampqqg->SetParameter(1, z);
mQ2=Q2;
- ROOT::Math::Functor wf(this, &InclusiveDiffractiveCrossSections::uiPhi_qqg, 2);
+ ROOT::Math::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 lo[2]={0.,1e-4}; //b, k2
double hi[2]={5., Q2}; //b, k2
double result = ig.Integral(lo, hi)/hbarc; //GeV0
return result; //GeV0
}
-double InclusiveDiffractiveCrossSections::uiPhi_qqg(const double* var) {
+double 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 InclusiveDiffractiveCrossSections::uiAmplitude_qqg(double *var, double *par){
+double InclusiveDiffractiveCrossSectionsIntegrals::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/InclusiveDiffractiveCrossSectionsFromTables.h
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSectionsFromTables.h (revision 0)
+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.h (revision 538)
@@ -0,0 +1,95 @@
+#ifndef InclusiveDiffractiveCrossSectionsFromTables_h
+#define InclusiveDiffractiveCrossSectionsFromTables_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"
+#include <vector>
+#include "InclusiveDiffractiveCrossSections.h"
+
+using namespace std;
+
+class InclusiveDiffractionCrossSectionsFromTables : public InclusiveDiffractiveCrossSections {
+
+public:
+ InclusiveDiffractionCrossSectionsFromTables();
+ ~InclusiveDiffractionCrossSectionsFromTables();
+
+ double operator()(double beta, double Q2, double W2, double z);
+ double operator()(const double* array);
+
+ double unuranPDF(const double*);
+ double unuranPDF_qqg(const double*);
+
+ GammaPolarization polarizationOfLastCall() ;
+ double crossSectionOfLastCall();
+
+ unsigned int quarkSpeciesOfLastCall();
+ double quarkMassOfLastCall();
+ double crossSectionRatioLTOfLastCall() const;
+ void setCheckKinematics(bool);
+ void setFockState(FockState);
+ FockState getFockState();
+
+ double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;
+
+ double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
+
+ 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 dsigmadbetadz_T(double, double, double, double) ;
+ double dsigmadbetadz_L(double, double, double, double) ;
+ double dsigmadbetadz_QQG(double, double, double, double) ;
+
+ double mDsigdbetadQ2dWdz_L_total[5];
+ double mDsigdbetadQ2dWdz_T_total[5];
+ double mDsigdbetadQ2dWdz_T_qqg[5];
+
+ void setQuarkIndex(unsigned int);
+
+ double interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark);
+
+ vector<double> mgrid(THnF* hist,int axis, int len);
+
+ vector<THnF*> tableQQ_T;
+ vector<THnF*> tableQQ_L;
+ vector<THnF*> tableQQG;
+
+ DipoleModel* dipoleModel();
+
+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;
+
+ DipoleModel* mDipoleModel;
+ DipoleModelParameterSet mDipoleModelParameterSet;
+
+
+};
+#endif
Index: trunk/src/linterp.h
===================================================================
--- trunk/src/linterp.h (revision 0)
+++ trunk/src/linterp.h (revision 538)
@@ -0,0 +1,421 @@
+
+//
+// Copyright (c) 2012 Ronaldo Carpio
+//
+// Permission to use, copy, modify, distribute and sell this software
+// and its documentation for any purpose is hereby granted without fee,
+// provided that the above copyright notice appear in all copies and
+// that both that copyright notice and this permission notice appear
+// in supporting documentation. The authors make no representations
+// about the suitability of this software for any purpose.
+// It is provided "as is" without express or implied warranty.
+//
+
+/*
+This is a C++ header-only library for N-dimensional linear interpolation on a rectangular grid. Implements two methods:
+* Multilinear: Interpolate using the N-dimensional hypercube containing the point. Interpolation step is O(2^N)
+* Simplicial: Interpolate using the N-dimensional simplex containing the point. Interpolation step is O(N log N), but less accurate.
+Requires boost/multi_array library.
+*/
+
+#ifndef _linterp_h
+#define _linterp_h
+
+#include <assert.h>
+#include <math.h>
+#include <stdarg.h>
+#include <float.h>
+#include <cstdarg>
+#include <string>
+#include <vector>
+#include <array>
+#include <functional>
+#include <boost/multi_array.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/storage.hpp>
+
+using std::vector;
+using std::array;
+typedef unsigned int uint;
+typedef vector<int> iVec;
+typedef vector<double> dVec;
+
+
+// TODO:
+// - specify behavior past grid boundaries.
+// 1) clamp
+// 2) return a pre-determined value (e.g. NaN)
+
+// compile-time params:
+// 1) number of dimensions
+// 2) scalar type T
+// 3) copy data or not (default: false). The grids will always be copied
+// 4) ref count class (default: none)
+// 5) continuous or not
+
+// run-time constructor params:
+// 1) f
+// 2) grids
+// 3) behavior outside grid: default=clamp
+// 4) value to return outside grid, defaut=nan
+
+struct EmptyClass {};
+
+template <int N, class T, bool CopyData = true, bool Continuous = true, class ArrayRefCountT = EmptyClass, class GridRefCountT = EmptyClass>
+class NDInterpolator {
+public:
+ typedef T value_type;
+ typedef ArrayRefCountT array_ref_count_type;
+ typedef GridRefCountT grid_ref_count_type;
+
+ static const int m_N = N;
+ static const bool m_bCopyData = CopyData;
+ static const bool m_bContinuous = Continuous;
+
+ typedef boost::numeric::ublas::array_adaptor<T> grid_type;
+ typedef boost::const_multi_array_ref<T, N> array_type;
+ typedef std::unique_ptr<array_type> array_type_ptr;
+
+ array_type_ptr m_pF;
+ ArrayRefCountT m_ref_F; // reference count for m_pF
+ vector<T> m_F_copy; // if CopyData == true, this holds the copy of F
+
+ vector<grid_type> m_grid_list;
+ vector<GridRefCountT> m_grid_ref_list; // reference counts for grids
+ vector<vector<T> > m_grid_copy_list; // if CopyData == true, this holds the copies of the grids
+
+ // constructors assume that [f_begin, f_end) is a contiguous array in C-order
+ // non ref-counted constructor.
+ template <class IterT1, class IterT2, class IterT3>
+ NDInterpolator(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end) {
+ init(grids_begin, grids_len_begin, f_begin, f_end);
+ }
+
+ // ref-counted constructor
+ template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
+ NDInterpolator(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT grid_refs_begin) {
+ init_refcount(grids_begin, grids_len_begin, f_begin, f_end, refF, grid_refs_begin);
+ }
+
+ template <class IterT1, class IterT2, class IterT3>
+ void init(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end) {
+ set_grids(grids_begin, grids_len_begin, m_bCopyData);
+ set_f_array(f_begin, f_end, m_bCopyData);
+ }
+ template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
+ void init_refcount(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT grid_refs_begin) {
+ set_grids(grids_begin, grids_len_begin, m_bCopyData);
+ set_grids_refcount(grid_refs_begin, grid_refs_begin + N);
+ set_f_array(f_begin, f_end, m_bCopyData);
+ set_f_refcount(refF);
+ }
+
+ template <class IterT1, class IterT2>
+ void set_grids(IterT1 grids_begin, IterT2 grids_len_begin, bool bCopy) {
+ m_grid_list.clear();
+ m_grid_ref_list.clear();
+ m_grid_copy_list.clear();
+ for (int i=0; i<N; i++) {
+ int gridLength = grids_len_begin[i];
+ if (bCopy == false) {
+ T const *grid_ptr = &(*grids_begin[i]);
+ m_grid_list.push_back(grid_type(gridLength, (T*) grid_ptr )); // use the given pointer
+ } else {
+ m_grid_copy_list.push_back( vector<T>(grids_begin[i], grids_begin[i] + grids_len_begin[i]) ); // make our own copy of the grid
+ T *begin = &(m_grid_copy_list[i][0]);
+ m_grid_list.push_back(grid_type(gridLength, begin)); // use our copy
+ }
+ }
+ }
+ template <class IterT1, class RefCountIterT>
+ void set_grids_refcount(RefCountIterT refs_begin, RefCountIterT refs_end) {
+ assert(refs_end - refs_begin == N);
+ m_grid_ref_list.assign(refs_begin, refs_begin + N);
+ }
+
+ // assumes that [f_begin, f_end) is a contiguous array in C-order
+ template <class IterT>
+ void set_f_array(IterT f_begin, IterT f_end, bool bCopy) {
+ unsigned int nGridPoints = 1;
+ array<int,N> sizes;
+ for (unsigned int i=0; i<m_grid_list.size(); i++) {
+ sizes[i] = m_grid_list[i].size();
+ nGridPoints *= sizes[i];
+ }
+
+ int f_len = f_end - f_begin;
+ if ( (m_bContinuous && f_len != nGridPoints) || (!m_bContinuous && f_len != 2 * nGridPoints) ) {
+ throw std::invalid_argument("f has wrong size");
+ }
+ for (unsigned int i=0; i<m_grid_list.size(); i++) {
+ if (!m_bContinuous) { sizes[i] *= 2; }
+ }
+
+ m_F_copy.clear();
+ if (bCopy == false) {
+ m_pF.reset(new array_type(f_begin, sizes));
+ } else {
+ m_F_copy = vector<T>(f_begin, f_end);
+ m_pF.reset(new array_type(&m_F_copy[0], sizes));
+ }
+ }
+ void set_f_refcount(ArrayRefCountT &refF) {
+ m_ref_F = refF;
+ }
+
+ // -1 is before the first grid point
+ // N-1 (where grid.size() == N) is after the last grid point
+ int find_cell(int dim, T x) const {
+ grid_type const &grid(m_grid_list[dim]);
+ if (x < *(grid.begin())) return -1;
+ else if (x >= *(grid.end()-1)) return grid.size()-1;
+ else {
+ auto i_upper = std::upper_bound(grid.begin(), grid.end(), x);
+ return i_upper - grid.begin() - 1;
+ }
+ }
+
+ // return the value of f at the given cell and vertex
+ T get_f_val(array<int,N> const &cell_index, array<int,N> const &v_index) const {
+ array<int,N> f_index;
+
+ if (m_bContinuous) {
+ for (int i=0; i<N; i++) {
+ if (cell_index[i] < 0) {
+ f_index[i] = 0;
+ } else if (cell_index[i] >= m_grid_list[i].size()-1) {
+ f_index[i] = m_grid_list[i].size()-1;
+ } else {
+ f_index[i] = cell_index[i] + v_index[i];
+ }
+ }
+ } else {
+ for (int i=0; i<N; i++) {
+ if (cell_index[i] < 0) {
+ f_index[i] = 0;
+ } else if (cell_index[i] >= m_grid_list[i].size()-1) {
+ f_index[i] = (2*m_grid_list[i].size())-1;
+ } else {
+ f_index[i] = 1 + (2*cell_index[i]) + v_index[i];
+ }
+ }
+ }
+ return (*m_pF)(f_index);
+ }
+
+ T get_f_val(array<int,N> const &cell_index, int v) const {
+ array<int,N> v_index;
+ for (int dim=0; dim<N; dim++) {
+ v_index[dim] = (v >> (N-dim-1)) & 1; // test if the i-th bit is set
+ }
+ return get_f_val(cell_index, v_index);
+ }
+};
+
+template <int N, class T, bool CopyData = true, bool Continuous = true, class ArrayRefCountT = EmptyClass, class GridRefCountT = EmptyClass>
+class InterpSimplex : public NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> {
+public:
+ typedef NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> super;
+
+ template <class IterT1, class IterT2, class IterT3>
+ InterpSimplex(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end)
+ : super(grids_begin, grids_len_begin, f_begin, f_end)
+ {}
+ template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
+ InterpSimplex(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT ref_begins)
+ : super(grids_begin, grids_len_begin, f_begin, f_end, refF, ref_begins)
+ {}
+
+ template <class IterT>
+ T interp(IterT x_begin) const {
+ array<T,1> result;
+ array< array<T,1>, N > coord_iter;
+ for (int i=0; i<N; i++) {
+ coord_iter[i][0] = x_begin[i];
+ }
+ interp_vec(1, coord_iter.begin(), coord_iter.end(), result.begin());
+ return result[0];
+ }
+
+ template <class IterT1, class IterT2>
+ void interp_vec(int n, IterT1 coord_iter_begin, IterT1 coord_iter_end, IterT2 i_result) const {
+ assert(N == coord_iter_end - coord_iter_begin);
+
+ array<int,N> cell_index, v_index;
+ array<std::pair<T, int>,N> xipair;
+ int c;
+ T y, v0, v1;
+ //mexPrintf("%d\n", n);
+ for (int i=0; i<n; i++) { // for each point
+ for (int dim=0; dim<N; dim++) {
+ typename super::grid_type const &grid(super::m_grid_list[dim]);
+ c = this->find_cell(dim, coord_iter_begin[dim][i]);
+ //mexPrintf("%d\n", c);
+ if (c == -1) { // before first grid point
+ y = 1.0;
+ } else if (c == grid.size()-1) { // after last grid point
+ y = 0.0;
+ } else {
+ //mexPrintf("%f %f\n", grid[c], grid[c+1]);
+ y = (coord_iter_begin[dim][i] - grid[c]) / (grid[c + 1] - grid[c]);
+ if (y < 0.0) y=0.0;
+ else if (y > 1.0) y=1.0;
+ }
+ xipair[dim].first = y;
+ xipair[dim].second = dim;
+ cell_index[dim] = c;
+ }
+ // sort xi's and get the permutation
+ std::sort(xipair.begin(), xipair.end(), [](std::pair<T, int> const &a, std::pair<T, int> const &b) {
+ return (a.first < b.first);
+ });
+ // walk the vertices of the simplex determined by the permutation
+ for (int j=0; j<N; j++) {
+ v_index[j] = 1;
+ }
+ v0 = this->get_f_val(cell_index, v_index);
+ y = v0;
+ for (int j=0; j<N; j++) {
+ v_index[xipair[j].second]--;
+ v1 = this->get_f_val(cell_index, v_index);
+ y += (1.0 - xipair[j].first) * (v1-v0); // interpolate
+ v0 = v1;
+ }
+ *i_result++ = y;
+ }
+ }
+};
+
+
+template <int N, class T, bool CopyData = true, bool Continuous = true, class ArrayRefCountT = EmptyClass, class GridRefCountT = EmptyClass>
+class InterpMultilinear : public NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> {
+public:
+ typedef NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> super;
+
+ template <class IterT1, class IterT2, class IterT3>
+ InterpMultilinear(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end)
+ : super(grids_begin, grids_len_begin, f_begin, f_end)
+ {}
+ template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
+ InterpMultilinear(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT ref_begins)
+ : super(grids_begin, grids_len_begin, f_begin, f_end, refF, ref_begins)
+ {}
+
+ template <class IterT1, class IterT2>
+ static T linterp_nd_unitcube(IterT1 f_begin, IterT1 f_end, IterT2 xi_begin, IterT2 xi_end) {
+ int n = xi_end - xi_begin;
+ int f_len = f_end - f_begin;
+ assert(1 << n == f_len);
+ T sub_lower, sub_upper;
+ if (n == 1) {
+ sub_lower = f_begin[0];
+ sub_upper = f_begin[1];
+ } else {
+ sub_lower = linterp_nd_unitcube(f_begin, f_begin + (f_len/2), xi_begin + 1, xi_end);
+ sub_upper = linterp_nd_unitcube(f_begin + (f_len/2), f_end, xi_begin + 1, xi_end);
+ }
+ T result = sub_lower + (*xi_begin)*(sub_upper - sub_lower);
+ return result;
+ }
+
+ template <class IterT>
+ T interp(IterT x_begin) const {
+ array<T,1> result;
+ array< array<T,1>, N > coord_iter;
+ for (int i=0; i<N; i++) {
+ coord_iter[i][0] = x_begin[i];
+ }
+ interp_vec(1, coord_iter.begin(), coord_iter.end(), result.begin());
+ return result[0];
+ }
+
+ template <class IterT1, class IterT2>
+ void interp_vec(int n, IterT1 coord_iter_begin, IterT1 coord_iter_end, IterT2 i_result) const {
+ assert(N == coord_iter_end - coord_iter_begin);
+ array<int,N> index;
+ int c;
+ T y, xi;
+ vector<T> f(1 << N);
+ array<T,N> x;
+
+ for (int i=0; i<n; i++) { // loop over each point
+ for (int dim=0; dim<N; dim++) { // loop over each dimension
+ auto const &grid(super::m_grid_list[dim]);
+ xi = coord_iter_begin[dim][i];
+ c = this->find_cell(dim, coord_iter_begin[dim][i]);
+ if (c == -1) { // before first grid point
+ y = 1.0;
+ } else if (c == grid.size()-1) { // after last grid point
+ y = 0.0;
+ } else {
+ y = (coord_iter_begin[dim][i] - grid[c]) / (grid[c + 1] - grid[c]);
+ if (y < 0.0) y=0.0;
+ else if (y > 1.0) y=1.0;
+ }
+ index[dim] = c;
+ x[dim] = y;
+ }
+ // copy f values at vertices
+ for (int v=0; v < (1 << N); v++) { // loop over each vertex of hypercube
+ f[v] = this->get_f_val(index, v);
+ }
+ *i_result++ = linterp_nd_unitcube(f.begin(), f.end(), x.begin(), x.end());
+ }
+ }
+};
+
+
+typedef InterpSimplex<1,double> NDInterpolator_1_S;
+typedef InterpSimplex<2,double> NDInterpolator_2_S;
+typedef InterpSimplex<3,double> NDInterpolator_3_S;
+typedef InterpSimplex<4,double> NDInterpolator_4_S;
+typedef InterpSimplex<5,double> NDInterpolator_5_S;
+
+typedef InterpMultilinear<1,double> NDInterpolator_1_ML;
+typedef InterpMultilinear<2,double> NDInterpolator_2_ML;
+typedef InterpMultilinear<3,double> NDInterpolator_3_ML;
+typedef InterpMultilinear<4,double> NDInterpolator_4_ML;
+typedef InterpMultilinear<5,double> NDInterpolator_5_ML;
+
+
+// C interface
+extern "C" {
+ void linterp_simplex_1(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult);
+ void linterp_simplex_2(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult);
+ void linterp_simplex_3(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult);
+}
+
+void inline linterp_simplex_1(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult) {
+ const int N=1;
+ size_t total_size = 1;
+ for (int i=0; i<N; i++) {
+ total_size *= grid_len_begin[i];
+ }
+ InterpSimplex<N, double, false> interp_obj(grids_begin, grid_len_begin, pF, pF + total_size);
+ interp_obj.interp_vec(xi_len, xi_begin, xi_begin + N, pResult);
+}
+
+void inline linterp_simplex_2(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult) {
+ const int N=2;
+ size_t total_size = 1;
+ for (int i=0; i<N; i++) {
+ total_size *= grid_len_begin[i];
+ }
+ InterpSimplex<N, double, false> interp_obj(grids_begin, grid_len_begin, pF, pF + total_size);
+ interp_obj.interp_vec(xi_len, xi_begin, xi_begin + N, pResult);
+}
+
+void inline linterp_simplex_3(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult) {
+ const int N=3;
+ size_t total_size = 1;
+ for (int i=0; i<N; i++) {
+ total_size *= grid_len_begin[i];
+ }
+ InterpSimplex<N, double, false> interp_obj(grids_begin, grid_len_begin, pF, pF + total_size);
+ interp_obj.interp_vec(xi_len, xi_begin, xi_begin + N, pResult);
+}
+
+
+
+#endif //_linterp_h
Index: trunk/src/InclusiveDiffractiveCrossSections.h
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSections.h (revision 537)
+++ trunk/src/InclusiveDiffractiveCrossSections.h (revision 538)
@@ -1,126 +1,170 @@
//==============================================================================
// 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 {
-
+
+class InclusiveDiffractiveCrossSections{
public:
InclusiveDiffractiveCrossSections();
- ~InclusiveDiffractiveCrossSections();
+ virtual ~InclusiveDiffractiveCrossSections();
+ virtual double operator()(double beta, double Q2, double W2, double z);
+ virtual double operator()(const double* array);
+
+ virtual DipoleModel* dipoleModel();
+ virtual double unuranPDF(const double*);
+ virtual double unuranPDF_qqg(const double*);
+ virtual void setCheckKinematics(bool);
+ virtual void setFockState(FockState);
+
+ virtual GammaPolarization polarizationOfLastCall() ;
+ virtual double crossSectionOfLastCall();
+ virtual unsigned int quarkSpeciesOfLastCall();
+ virtual double quarkMassOfLastCall();
+ virtual double crossSectionRatioLTOfLastCall() const;
+
+ virtual void setQuarkIndex(unsigned int);
+
+ virtual double dsigmadbetadz_T(double, double, double, double) ;
+ virtual double dsigmadbetadz_L(double, double, double, double) ;
+
+ virtual double dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z);
+ virtual double dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z);
+
+ virtual double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;
+
+ virtual double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
+
+ virtual double* dsigdbetadQ2dWdz_T_total();
+ virtual double* dsigdbetadQ2dWdz_L_total();
+ virtual double* dsigdbetadQ2dWdz_T_qqg();
+
+};
+
+class InclusiveDiffractiveCrossSectionsIntegrals: public InclusiveDiffractiveCrossSections{
+
+public:
+ InclusiveDiffractiveCrossSectionsIntegrals();
+ ~InclusiveDiffractiveCrossSectionsIntegrals();
double operator()(double beta, double Q2, double W2, double z);
double operator()(const double* array);
double unuranPDF(const double*);
double unuranPDF_qqg(const double*);
GammaPolarization polarizationOfLastCall() ;
double crossSectionOfLastCall();
unsigned int quarkSpeciesOfLastCall();
double quarkMassOfLastCall();
double crossSectionRatioLTOfLastCall() const;
void setCheckKinematics(bool);
void setFockState(FockState);
FockState getFockState();
double dsigdbetadQ2dW2dz_total(double beta, double Q2, double W2, double z, GammaPolarization) ;
double dsigdbetadQ2dW2dz_qqg(double beta, double Q2, double W2, double z);
// InclusiveDiffractiveCrossSections(const InclusiveDiffractiveCrossSections&);
// InclusiveDiffractiveCrossSections& operator=(const InclusiveDiffractiveCrossSections&);
double dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z);
double dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z);
DipoleModel* dipoleModel();
double dsigmadbetadz_T(double, double, double, double) ;
double dsigmadbetadz_L(double, double, double, double) ;
double dsigmadbetadz_QQG(double, double, double, double) ;
- double mDsigdbetadQ2dWdz_L_total[5];
- double mDsigdbetadQ2dWdz_T_total[5];
- double mDsigdbetadQ2dWdz_T_qqg[5];
+ double* dsigdbetadQ2dWdz_T_total(){return mDsigdbetadQ2dWdz_L_total;}
+ double* dsigdbetadQ2dWdz_L_total(){return mDsigdbetadQ2dWdz_T_total;}
+ double* dsigdbetadQ2dWdz_T_qqg(){ return mDsigdbetadQ2dWdz_T_qqg;};
void setQuarkIndex(unsigned int);
protected:
double dsigdMX2dQ2dW2dz_total_checked(double MX2, double Q2, double W2, double z);
double dsigdbetadz_total(double beta, double Q2, double W2, double z, GammaPolarization pol);
private:
TRandom3* mRandom;
GammaPolarization mPolarization;
PhotonFlux mPhotonFlux;
EventGeneratorSettings* mSettings;
unsigned int mQuarkSpecies;
double mTotalCS, mTotalCS_qqg;
FockState mFockState;
double mS, mQ2;
bool mCheckKinematics;
double Phi_0(double, double, double, double, double);
double Phi_1(double, double, double, double, double);
double Phi_qqg(double, double, double);
double uiPhi_0(double*, double*);
double uiPhi_1(double*, double*);
double Amplitude_0(double);
double Amplitude_1(double);
double uiAmplitude_0(double*, double*);
double uiAmplitude_1(double*, double*);
double uiPhi_qqg(const double*);
double uiAmplitude_qqg(double*, double*);
unsigned int mA;
unsigned int mQuarkIndex;
double mCrossSectionRatioLT;
+ double mDsigdbetadQ2dWdz_L_total[5];
+ double mDsigdbetadQ2dWdz_T_total[5];
+ double mDsigdbetadQ2dWdz_T_qqg[5];
+
DipoleModel* mDipoleModel;
DipoleModelParameterSet mDipoleModelParameterSet;
TF1 *Amp0, *Amp1, *Ampqqg;
ROOT::Math::WrappedTF1 *WFAmp0, *WFAmp1, *WFAmpqqg;
ROOT::Math::GaussIntegrator GIAmp0, GIAmp1, GIAmpqqg;
};
#endif
Index: trunk/src/SartreInclusiveDiffraction.h
===================================================================
--- trunk/src/SartreInclusiveDiffraction.h (revision 537)
+++ trunk/src/SartreInclusiveDiffraction.h (revision 538)
@@ -1,120 +1,124 @@
//==============================================================================
// 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*);
};
ostream& operator<<(ostream& os, const TLorentzVector&);
#endif
Index: trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 0)
+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 538)
@@ -0,0 +1,579 @@
+#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(){
+
+ mSettings = EventGeneratorSettings::instance();
+ mRandom = mSettings->randomGenerator();
+ mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam());
+ mPhotonFlux.setS(mS);
+ mCheckKinematics = true;
+ mCrossSectionRatioLT = 0;
+
+ mQ2=0;
+
+ mA=mSettings->A();
+
+ 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);
+ }
+
+ //
+ // Set up the lookup tables:
+ //
+ TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
+
+ TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
+
+ TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables/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);
+
+ // if (!hist4D) {
+ // std::cerr << "Error: Histogram 'hist4D' not found!" << std::endl;
+ // file->Close();
+ // return -1;
+ // }
+}
+InclusiveDiffractionCrossSectionsFromTables::~InclusiveDiffractionCrossSectionsFromTables(){
+}
+
+void InclusiveDiffractionCrossSectionsFromTables::setCheckKinematics(bool val) {mCheckKinematics = val;}
+
+void InclusiveDiffractionCrossSectionsFromTables::setFockState(FockState val){
+ mFockState=val;
+}
+
+FockState InclusiveDiffractionCrossSectionsFromTables::getFockState(){
+ return mFockState;
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::operator()(double MX2, double Q2, double W2, double z){
+ return dsigdMX2dQ2dW2dz_total_checked(MX2, Q2, W2, z);
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::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<<"InclusiveDiffractiveCrossSections::operator(): Fockstate is invalid, stopping"<<endl;
+ exit(0);
+ }
+ return result;
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::dsigdbetadQ2dW2dz_total_checked(double beta, double Q2, double W2, double z){
+ double result = 0;
+
+ //
+ // Check if in kinematically allowed region
+ double MX2=Kinematics::MX2(Q2, beta, 0);
+ if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
+ if (mSettings->verboseLevel() > 2)
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): warning, beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z="<< z
+ << " is outside of kinematically allowed region. Return 0." << endl;
+ return 0;
+ }
+ double csT = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, transverse);
+ double csL = dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, longitudinal);
+ result = csT + csL;
+
+ mCrossSectionRatioLT=csL/csT;
+ //
+ // Polarization
+ //
+ if (mRandom->Uniform(result) <= csT)
+ mPolarization = transverse;
+ else
+ mPolarization = longitudinal;
+ //
+ // Quark Species
+ //
+ //
+ if(mPolarization == transverse){
+ bool isChosen=false;
+ double cs_rejected=0;
+ for(int i=0; i<4; i++){
+ double R=mRandom->Uniform(csT-cs_rejected);
+ double csTi=mDsigdbetadQ2dWdz_T_total[i];
+ if(R <= csTi){
+ mQuarkSpecies=i;
+ isChosen=true;
+ break;
+ }
+ cs_rejected+=mDsigdbetadQ2dWdz_T_total[i];
+ }
+ if(!isChosen) mQuarkSpecies=4;
+ }
+ else{
+ bool isChosen=false;
+ double cs_rejected=0;
+ for(int i=0; i<4; i++){
+ double R=mRandom->Uniform(csL-cs_rejected);
+ double csLi=mDsigdbetadQ2dWdz_L_total[i];
+ if(R <= csLi){
+ mQuarkSpecies=i;
+ isChosen=true;
+ break;
+ }
+ cs_rejected+=mDsigdbetadQ2dWdz_L_total[i];
+ }
+ if(!isChosen) mQuarkSpecies=4;
+ }
+ // Print-out at high verbose levels
+ //
+ if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-)
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): " << result;
+ cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z;
+ cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal");
+ cout << ')' << endl;
+ }
+ //
+ // Check validity of return value
+ //
+ if (std::isnan(result)) {
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = NaN at" << endl;
+ cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
+ result = 0;
+ }
+ if (std::isinf(result)) {
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, return value = inf at" << endl;
+ cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
+ result = 0;
+ }
+ if (result < 0) {
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_checked(): Error, negative cross-section at" << endl;
+ cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
+ result = 0;
+ }
+ mTotalCS=result;
+ return result; //nb/GeV4
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::crossSectionOfLastCall(){
+ return mTotalCS;
+}
+
+GammaPolarization InclusiveDiffractionCrossSectionsFromTables::polarizationOfLastCall(){
+ return mPolarization;
+}
+
+unsigned int InclusiveDiffractionCrossSectionsFromTables::quarkSpeciesOfLastCall(){
+ return mQuarkSpecies;
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::quarkMassOfLastCall(){
+ // return Settings::quarkMass(mQuarkSpecies); //#TT tmp
+ return tt_quarkMass[mQuarkSpecies];
+}
+
+void InclusiveDiffractionCrossSectionsFromTables::setQuarkIndex(unsigned int val){mQuarkIndex=val;}
+
+double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {
+
+ double MX2 = Q2*(1-beta)/beta;
+ // double mf = Settings::quarkMass(mQuarkIndex);
+ double mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
+ double sqrtArg = 1.-4.*mf*mf/MX2;
+ if (sqrtArg<0){
+ if (mSettings->verboseLevel() > 2)
+ cout<<"There is no phase-space for z, return 0."<<endl;
+ return 0;
+ }
+ double z0 = (1.-sqrt(sqrtArg))/2.;
+ if (z<z0 or z>1-z0){
+ if (mSettings->verboseLevel() > 2)
+ cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
+ return 0;
+ }
+ double pt2 = z*(1-z)*Q2-mf*mf;
+ if (pt2<0){
+ if (mSettings->verboseLevel() > 2)
+ cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
+ return 0;
+ }
+ if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
+ if (mSettings->verboseLevel() > 2)
+ cout<<"Not enough phase-space for quark production, return 0."<<endl;
+ return 0;
+ }
+ double result = interpolate(beta, Q2, W2, z, transverse, QQ, mQuarkIndex);
+ 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 mf = tt_quarkMass[mQuarkIndex]; //#TT tmp
+ double pt2 = z*(1-z)*Q2-mf*mf;
+ if (pt2<0){
+ if (mSettings->verboseLevel() > 2)
+ cout<<"Not enough phase-space for quark transverse momenta, return 0."<<endl;
+ return 0;
+ }
+ double sqrtArg=1.-4.*mf*mf/MX2;
+ if (sqrtArg<0){
+ if (mSettings->verboseLevel() > 2)
+ cout<<"There is no phase-space for z, return 0."<<endl;
+ return 0;
+ }
+ double z0=(1.-sqrt(sqrtArg))/2.;
+ if (z<z0 or z>1-z0){
+ if (mSettings->verboseLevel() > 2)
+ cout<<"z="<<z<<" which is smaller than z0="<<z0<<"or larger than (1-z0)="<<1-z0<<", return 0."<<endl;
+ return 0;
+ }
+ if(z*(1-z)*MX2-mf*mf<0){ //kappa2<0
+ if (mSettings->verboseLevel() > 2)
+ cout<<"Not enough phase-space for quark production, return 0."<<endl;
+ return 0;
+ }
+ double result = interpolate(beta, Q2, W2, z, longitudinal, QQ, mQuarkIndex);
+ return result;
+}
+
+DipoleModel* InclusiveDiffractionCrossSectionsFromTables::dipoleModel(){
+ return mDipoleModel;
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;}
+
+//===================================================
+// QQG Calculations
+//===================================================
+
+double InclusiveDiffractionCrossSectionsFromTables::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 InclusiveDiffractionCrossSectionsFromTables::dsigdbetadQ2dW2dz_total_qqg_checked(double beta, double Q2, double W2, double z){
+ //
+ // Check if in kinematically allowed region
+ //
+ double MX2=Kinematics::MX2(Q2, beta, 0);
+ if(mCheckKinematics && z<beta) {
+ if (mSettings->verboseLevel() > 2)
+ cout<<"z<beta, beta="<<beta<<", z="<<z<<endl;
+ return 0; // For QQG beta < z < 1
+ }
+ if (mCheckKinematics && !Kinematics::valid(mS, beta, Q2, W2, z, sqrt(MX2), false, (mSettings->verboseLevel() > 2) )) {
+ if (mSettings->verboseLevel() > 2)
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_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 << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): " << result;
+ cout << " at beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z<< endl;
+ }
+ //
+ // Check validity of return value
+ //
+ if (std::isnan(result)) {
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = NaN at" << endl;
+ cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
+ result = 0;
+ }
+ if (std::isinf(result)) {
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, return value = inf at" << endl;
+ cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << z << endl;
+ result = 0;
+ }
+ if (result < 0) {
+ cout << "InclusiveDiffractiveCrossSections::dsigdbetadQ2dW2dz_total_qqg_checked(): Error, negative cross-section at" << endl;
+ cout << " beta=" << beta << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", z=" << endl;
+ result = 0;
+ }
+ mTotalCS_qqg=result;
+ return result; //nb/GeV4
+}
+
+double InclusiveDiffractionCrossSectionsFromTables::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 mf = tt_quarkMass[i]; //#TT tmp
+ double Mqq2=(z/beta-1)*Q2;
+ if(Mqq2<4*mf*mf) {
+ continue; //Not enough phase-space for the q and qbar
+ }
+ double tmpresult = dsigmadbetadz_QQG(beta, Q2, W2, z); //nb
+ if (mSettings->applyPhotonFlux()) tmpresult *= mPhotonFlux(Q2,W2,transverse); //nb/GeV4
+ mDsigdbetadQ2dWdz_T_qqg[i]=tmpresult;
+ result+=tmpresult;
+ }
+ return result; //nb/GeV4
+}
+
+double 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;// = dynamic_cast<THnF*>(file->Get("hist4D"));
+ if(pol==transverse and fock==QQ){
+ hist4D=tableQQ_T.at(quark);
+ }
+ else if(pol==longitudinal and fock==QQ){
+ hist4D=tableQQ_L.at(quark);
+ }
+ else if(fock==QQG){
+ hist4D=tableQQG.at(quark);
+ }
+ else{
+ cout<<"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]) ;
+
+ // print_grid_vector(gridQ2, "gridQ2");
+ // print_grid_vector(gridW2, "gridW2");
+ // print_grid_vector(gridBeta, "gridBeta");
+ // print_grid_vector(gridZ, "gridZ");
+
+ std::vector< std::vector<double>::iterator > grid_iter_list;
+ grid_iter_list.push_back(gridQ2.begin());
+ grid_iter_list.push_back(gridW2.begin());
+ grid_iter_list.push_back(gridBeta.begin());
+ grid_iter_list.push_back(gridZ.begin());
+
+ // cout<<"grid\t"<<*(gridQ2.begin())<<endl;
+ // the size of the grid in each dimension
+ array<int,4> grid_sizes;
+ grid_sizes[0] = Bins[0];
+ grid_sizes[1] = Bins[1];
+ grid_sizes[2] = Bins[2];
+ grid_sizes[3] = Bins[3];
+
+
+ // for(unsigned long i=0;i<grid_sizes.size(); i++)
+ // std::cout<<"Grid size \t"<<grid_sizes[i]<<std::endl;
+
+ // total number of elements
+ int num_elements = grid_sizes[0] * grid_sizes[1]* grid_sizes[2]* grid_sizes[3];
+ // std::cout<<"num_elements = "<<num_elements<<std::endl;
+ // 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;
+
+
+ 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[] = {ibeta, iQ2, iW2, iZ};
+ f_values[count] = hist4D->GetBinContent(binIndex);
+ count++;
+ }
+ }
+ }
+ }
+ // std::cout<<"done"<<std::endl;
+
+ // 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);
+
+ // InterpMultilinear 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 exp(result);
+}
+
+vector<double> InclusiveDiffractionCrossSectionsFromTables::mgrid(THnF* hist,int axis, int len){
+ vector<double> gridtype(len);
+ for(int i = 0; i <len; i++){
+ double val = hist->GetAxis(axis)->GetBinCenter(i+1);
+ gridtype[i] = val;
+ }
+ return gridtype;
+}
+
Index: trunk/src/SartreInclusiveDiffraction.cpp
===================================================================
--- trunk/src/SartreInclusiveDiffraction.cpp (revision 537)
+++ trunk/src/SartreInclusiveDiffraction.cpp (revision 538)
@@ -1,852 +1,861 @@
//==============================================================================
// 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"
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;
+ mUseCrossSectionTables=true;
}
SartreInclusiveDiffraction::~SartreInclusiveDiffraction()
{
delete mNucleus;
delete mUpcNucleus;
delete mPDF_Functor;
delete mPDF;
delete mCrossSection;
delete mUnuran;
delete mCurrentEvent;
}
bool SartreInclusiveDiffraction::init(const char* runcard)
{
mStartTime = time(0);
bool ok;
//
// Reset member variables.
// Note that one instance of Sartre should be able to get
// initialized multiple times.
//
mEvents = 0;
mTries = 0;
mTotalCrossSection = 0;
mTotalCrossSection_qqg = 0;
mTotalCrossSection_qq = 0;
//
// Print header
//
string ctstr(ctime(&mStartTime));
ctstr.erase(ctstr.size()-1, 1);
cout << "/========================================================================\\" << endl;
cout << "| |" << endl;
cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl;
cout << "| |" << endl;
cout << "| An event generator for inclusive diffraction |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
cout << "| Copyright (C) 2010-2024 Tobias Toll and Thomas Ullrich |" << endl;
cout << "| |" << endl;
cout << "| This program is free software: you can redistribute it and/or modify |" << endl;
cout << "| it under the terms of the GNU General Public License as published by |" << endl;
cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl;
cout << "| any later version. |" << endl;
cout << "| |" << endl;
cout << "| Code compiled on " << setw(12) << left << __DATE__;
cout << setw(41) << left << __TIME__ << right << '|' << endl;
cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl;
cout << "\\========================================================================/" << endl;
mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton
mSettings->setInclusiveDiffractionMode(true);
//
// Read runcard if available
//
if (runcard) {
if (!mSettings->readSettingsFromFile(runcard)) {
cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl;
exit(1);
}
else
cout << "Runcard is '" << runcard << "'." << endl;
}
else
cout << "No runcard provided." << endl;
//
// Set beam particles and center of mass energy
//
mElectronBeam = mSettings->eBeam();
mHadronBeam = mSettings->hBeam();
mS = Kinematics::s(mElectronBeam, mHadronBeam);
mA = mSettings->A();
-
+
bool allowBreakup = mSettings->enableNuclearBreakup();
if (mA == 1) allowBreakup = false;
if (allowBreakup) {
if (!getenv("SARTRE_DIR")) {
cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n"
"to locate tables needed for the generation if nuclear breakups." << endl;
exit(1);
}
}
if (mNucleus) delete mNucleus;
mNucleus = new FrangibleNucleus(mA, allowBreakup);
string upcNucleusName;
cout << "Sartre is running in inclusive diffraction mode" << endl;
cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl;
cout << "Hadron beam: " << mHadronBeam << endl;
cout << "Electron beam: " << mElectronBeam << endl;
//
// Get details about the processes and models
//
mDipoleModelType = mSettings->dipoleModelType();
mDipoleModelParameterSet = mSettings->dipoleModelParameterSet();
cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
cout << "Process is ";
if (mA > 1)
cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + X"<< endl;
else
cout << "e + p -> e' + p' + X"<< endl;
//
// Print-out seed for reference
//
cout << "Random generator seed: " << mSettings->seed() << endl;
//
// Load in the tables containing the amplitude moments
//
if (!getenv("SARTRE_DIR")) {
cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
exit(1);
}
//
// Kinematic limits and generator range
//
// There are 3 ranges we have to deal with
// 1. the kinematic range requested by the user
// if given.
// The user can only control Q2 and W but not t.
// For UPC that's xpom.
// 2. the range of the table(s)
// 3. the kinematically allowed range
//
// Of course (3) is more complex than a simple cube/square.
// However, we deal with the detailed shape of the kinematic
// range later using Kinematics::valid() when we generate the
// individual events.
// For setting up UNU.RAN we have to get the cubic/square
// envelope that satifies (1)-(3).
// Note, that they are correlated which makes the order
// in which we do things a bit tricky.
//
//
// Step 2:
// Kinematic limits might overrule boundaries from step 1
//
// let:
// 0: beta
// 1: Q2
// 2: W2
// 3: z
//
// Setup cross-section functor
// It is this functor that is used by all other functors,
// functions, and wrappers when dealing with cross-sections.
//
if (mCrossSection) delete mCrossSection;
- mCrossSection = new InclusiveDiffractiveCrossSections;
+ if(!mUseCrossSectionTables){
+ mCrossSection = new InclusiveDiffractiveCrossSections;
+ }
+ 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);
}
//
// 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<1; i++) {
+ 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;
+// 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]));
+ 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])); //z
+ 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);
+// mCrossSection->setCheckKinematics(true);
//
// Burn in QQG generator
//
- for (int i=0; i<1; i++) {
+ 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;
+// cout<<"QQG: beta="<<xrandom[0];
+// cout<<" Q2="<<exp(xrandom[1]);
+// cout<<" W="<<sqrt(xrandom[2]);
+// cout<<" z="<<xrandom[3]<<endl;
}
- mEventCounter = 0;
- mTriesCounter = 0;
- mIsInitialized = true;
- cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;
//
// Calculate total cross section in given
// kinematic limits. This will be used to
// decide which type of event to generate,
// q-qbar or q-qbar-g
//
- cout<<"Calculating cross sections over the given kinematic range..."<<endl;
- cout<<"(This may take while, use your radical freedom and do other things.)"<<endl;
- mCrossSection->setFockState(QQ);
- mTotalCrossSection=0;
- mTotalCrossSection_qq=totalCrossSection();
- cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
- mCrossSection->setFockState(QQG);
- mTotalCrossSection=0;
- mTotalCrossSection_qqg=totalCrossSection();
- cout<<"The QQG cross section is: "<<mTotalCrossSection_qqg<<" nb. (2/2)"<<endl;
- cout<<"Total Cross Section: "<<mTotalCrossSection_qqg+mTotalCrossSection_qq<<" nb"<<endl;
- mTotalCrossSection=mTotalCrossSection_qqg+mTotalCrossSection_qq;
+// cout<<"Calculating cross sections over the given kinematic range..."<<endl;
+// cout<<"(This may take while, use your radical freedom and do other things.)"<<endl;
+// mCrossSection->setFockState(QQ);
+// mTotalCrossSection=0;
+// mTotalCrossSection_qq=totalCrossSection();
+// cout<<" The QQ cross section is: "<<mTotalCrossSection_qq<<" nb. (1/2)"<<endl;
+// mCrossSection->setFockState(QQG);
+// mTotalCrossSection=0;
+// mTotalCrossSection_qqg=totalCrossSection();
+// cout<<"The QQG cross section is: "<<mTotalCrossSection_qqg<<" nb. (2/2)"<<endl;
+// cout<<"Total Cross Section: "<<mTotalCrossSection_qqg+mTotalCrossSection_qq<<" nb"<<endl;
+// mTotalCrossSection=mTotalCrossSection_qqg+mTotalCrossSection_qq;
+
+ mEventCounter = 0;
+ mTriesCounter = 0;
+ mIsInitialized = true;
+ cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;
+
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
//
double dsigmadb2=0;
double z=xrandom[3];
double eps2=z*(1-z)*(mCurrentEvent->Q2+MX2);
double r_average=1./sqrt(eps2);
if (mA == 1) {
dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(r_average*hbarc, 0, mCurrentEvent->xpom);
}
else {
dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(r_average*hbarc, 0, mCurrentEvent->xpom);
}
mCurrentEvent->Omega=-2.*log(1-dsigmadb2/2.);
Particle eIn, hIn;
eIn.index = 0;
eIn.pdgId = 11; // e-
eIn.status = 1;
eIn.p = mElectronBeam;
hIn.index = 1;
hIn.pdgId = mNucleus->pdgID();
hIn.status = 1;
hIn.p = mHadronBeam;
mCurrentEvent->particles.push_back(eIn);
mCurrentEvent->particles.push_back(hIn);
//
// Generate the final state particles
//
bool ok;
if (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 = 1e-1;//5e-4;
- unsigned int numberofcalls = 1;//1000000;
+ 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;
}
const FrangibleNucleus* SartreInclusiveDiffraction::nucleus() const {return mNucleus;}
void SartreInclusiveDiffraction::listStatus(ostream& os) const
{
os << "Event summary: " << mEventCounter<< " events generated, " << mTriesCounter << " tried" << endl;
time_t delta = runTime();
os << "Total time used: " << delta/60 << " min " << delta - 60*(delta/60) << " sec" << endl;
}
time_t SartreInclusiveDiffraction::runTime() const
{
time_t now = time(0);
return now-mStartTime;
}
//==============================================================================
//
// Utility functions and operators (helpers)
//
//==============================================================================
ostream& operator<<(ostream& os, const TLorentzVector& v)
{
os << v.Px() << '\t' << v.Py() << '\t' << v.Pz() << '\t' << v.E() << '\t';
double m2 = v*v;
if (m2 < 0)
os << '(' << -sqrt(-m2) << ')';
else
os << '(' << sqrt(m2) << ')';
return os;
}
Index: trunk/src/CMakeLists.txt
===================================================================
--- trunk/src/CMakeLists.txt (revision 537)
+++ trunk/src/CMakeLists.txt (revision 538)
@@ -1,181 +1,182 @@
#===============================================================================
# CMakeLists.txt (src)
#
# Copyright (C) 2010-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: Thomas Ullrich
# Last update:
# $Date$
# $Author$
#===============================================================================
include(ExternalProject)
cmake_minimum_required (VERSION 3.5)
#
# Compiler flags for release and debug version
#
set(CMAKE_C_FLAGS "-W")
set(CMAKE_C_FLAGS_DEBUG "-g")
set(CMAKE_C_FLAGS_RELEASE "-O")
message("COMPILER = ${CMAKE_CXX_COMPILER_ID}")
set(CMAKE_CXX_FLAGS "-W -Wall -Wextra -pedantic -Wno-long-long")
if(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-potentially-evaluated-expression")
endif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
set(CMAKE_CXX_FLAGS_DEBUG "-g")
set(CMAKE_CXX_FLAGS_RELEASE "-O")
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
#
# Find external required packages
# (see also FindGSL.cmake and FindROOT.cmke in cmake/modules)
# Herer we need only the root library to create the table
# tools.
#
# GSL
find_package(GSL REQUIRED)
include_directories(${GSL_INCLUDE_DIR})
# ROOT
find_package(ROOT REQUIRED)
include_directories(${ROOT_INCLUDE_DIR})
set(LIBS ${LIBS} ${ROOT_LIBRARIES})
# strip leading and trailing blanks
string(REGEX REPLACE "\n$" "" LIBS "${LIBS}")
string(STRIP "${LIBS}" LIBS)
#BOOST
-if (MULTITHREADED)
+#if (MULTITHREADED)
set(Boost_USE_STATIC_LIBS ON)
set(Boost_USE_MULTITHREADED ON)
find_package(Boost 1.39 COMPONENTS thread REQUIRED)
if(Boost_FOUND)
include_directories(${Boost_INCLUDE_DIR})
endif(Boost_FOUND)
-endif (MULTITHREADED)
+#endif (MULTITHREADED)
# strip leading and trailing blanks
string(REGEX REPLACE "\n$" "" LIBS "${LIBS}")
string(STRIP "${LIBS}" LIBS)
# PYTHIA8
find_package(Pythia8 REQUIRED)
include_directories(${PYTHIA8_INCLUDE_DIR})
set(LIBS ${LIBS} ${PYTHIA8_LIBRARIES})
# strip leading and trailing blanks
string(REGEX REPLACE "\n$" "" LIBS "${LIBS}")
string(STRIP "${LIBS}" LIBS)
#
# Include files from sartre package
#
include_directories(${PROJECT_SOURCE_DIR}/src)
include_directories(${PROJECT_SOURCE_DIR}/gemini)
include_directories(${PROJECT_SOURCE_DIR}/cuba)
#
# Cuba is built using the autoconf shipped with it
#
ExternalProject_Add(
cuba
DOWNLOAD_COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/cuba ${PROJECT_BINARY_DIR}/cuba
SOURCE_DIR ${PROJECT_BINARY_DIR}/cuba
PREFIX ${PROJECT_BINARY_DIR}/cuba
CONFIGURE_COMMAND ${PROJECT_BINARY_DIR}/cuba/configure --prefix=${PROJECT_BINARY_DIR}/cuba/build
BUILD_COMMAND make lib
BUILD_IN_SOURCE 1
)
#
# Defines source files for sartre library
#
set(SARTRE_SRC "AlphaStrong.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Amplitudes.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "BreakupProduct.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "CrossSection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DglapEvolution.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModel.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "DipoleModelParameters.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "EicSmearFormatWriter.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Event.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "EventGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ExclusiveFinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "FrangibleNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractiveCrossSections.cpp")
+set(SARTRE_SRC ${SARTRE_SRC} "InclusiveDiffractiveCrossSectionsFromTables.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "InclusiveFinalStateGenerator.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Integrals.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Kinematics.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "ModeFinderFunctor.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleon.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Nucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "PhotonFlux.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Sartre.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "SartreInclusiveDiffraction.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Settings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "Table.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableCollection.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorNucleus.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TableGeneratorSettings.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "TwoBodyVectorMesonDecay.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "VectorMesonDecayMass.cpp")
set(SARTRE_SRC ${SARTRE_SRC} "WaveOverlap.cpp")
add_library(sartre ${SARTRE_SRC})
#
# Table tools (all stand alone programs)
#
add_executable(tableInspector tableInspector.cpp)
add_executable(tableMerger tableMerger.cpp)
add_executable(tableQuery tableQuery.cpp)
add_executable(tableDumper tableDumper.cpp)
add_executable(tableVarianceMaker tableVarianceMaker.cpp)
target_link_libraries(tableInspector sartre ${LIBS})
target_link_libraries(tableMerger sartre ${LIBS})
target_link_libraries(tableQuery sartre ${LIBS})
target_link_libraries(tableDumper sartre ${LIBS})
target_link_libraries(tableVarianceMaker sartre ${LIBS})
add_dependencies(tableInspector sartre)
add_dependencies(tableMerger sartre)
add_dependencies(tableQuery sartre)
add_dependencies(tableDumper sartre)
add_dependencies(tableVarianceMaker sartre)
#
# Install library and include files (make install) within
# the distribution tree. Top level CMakeLists.txt will install
# Sartre in final destination.
#
install(TARGETS sartre DESTINATION sartre/lib)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/lib/libcuba.a" DESTINATION sartre/lib)
FILE(GLOB AllIncludeFiles *.h)
install(FILES ${AllIncludeFiles} DESTINATION sartre/include)
install(FILES "${PROJECT_BINARY_DIR}/cuba/build/include/cuba.h" DESTINATION sartre/include)
install(TARGETS tableInspector DESTINATION sartre/bin)
install(TARGETS tableMerger DESTINATION sartre/bin)
install(TARGETS tableQuery DESTINATION sartre/bin)
install(TARGETS tableDumper DESTINATION sartre/bin)
install(TARGETS tableVarianceMaker DESTINATION sartre/bin)

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:01 PM (23 h, 9 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242367
Default Alt Text
(126 KB)

Event Timeline