Page MenuHomeHEPForge
Files F8308926

No OneTemporary
Actions

View Options

Index: trunk/examples/tableGeneratorInclusiveDiffractionMain.cpp
===================================================================
--- trunk/examples/tableGeneratorInclusiveDiffractionMain.cpp (revision 540)
+++ trunk/examples/tableGeneratorInclusiveDiffractionMain.cpp (revision 541)
@@ -1,414 +1,419 @@
//==============================================================================
// tableGeneratorMain.cpp
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date: 2023-07-26 11:48:08 +0530 (Wed, 26 Jul 2023) $
// $Author: ttoll $
//==============================================================================
//
// Main program to create amplitude lookup tables.
// [Developer only]
//
// Usage:
// tableGeneratorMain runcard startBin endBin
//
// Bins run from 0 to nbin-1.
// Loop fill all bins from startBIn to endBin (including endBin).
// If endBin exceeds the size of the table it set to the nbin-1.
//
//==============================================================================
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <cmath>
#include <iomanip>
#include "Amplitudes.h"
#include "TROOT.h"
#include "TH1D.h"
#include "TFile.h"
#include "Nucleus.h"
#include "Constants.h"
#include "Table.h"
#include "TableGeneratorSettings.h"
#include "Enumerations.h"
#include "DglapEvolution.h"
#include "Version.h"
#include "TRandom3.h"
#include "InclusiveDiffractiveCrossSections.h"
#include "DipoleModel.h"
//void myFunction(){
// DipoleModel myDM;
// ///intialise myDM
//
// double b=0;
// double xprobe=1e-3;
//
// //solve for r:
// 2*myDM.dsigmadb2ep(r, b, xprobe)=1-exp(-0.5);
//
//}
//
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
+double testFunction(double beta, double Q2, double W2, double z){
+ return beta+Q2+W2+z;
+}
int main(int argc, char *argv[]) {
//
// Print header
//
time_t theStartTime = time(0);
string ctstr(ctime(&theStartTime));
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 exclusive diffractive vector meson production |" << endl;
cout << "| in ep and eA collisions based on the dipole model. |" << endl;
cout << "| |" << endl;
cout << "| Copyright (C) 2010-2018 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;
cout << "| -------------------------------------------------------------------- |" << endl;
cout << "| |" << endl;
cout << "| Sartre Table Generator (Experts only) |" << endl;
cout << "| |" << endl;
cout << "\\========================================================================/" << endl;
TH1::AddDirectory(false); // to explicitly delete all new histograms by hand
TableGeneratorSettings* settings = TableGeneratorSettings::instance();
//
// Check arguments
//
char* runcard;
if (argc != 4) {
cout << "Usage: tableGeneratorMain runcard startBin endBin" << endl;
return 2;
}
else {
runcard = argv[1];
settings->setStartBin(atoi(argv[2]));
settings->setEndBin(atoi(argv[3]));
}
cout << "Reading settings from runcard." << endl;
settings->readSettingsFromFile(runcard);
settings->consolidateSettings();
//dipoleModelParameterSet,
// TODO: eBeam, hBeam,
//dipoleModelType, A, verboseLevel, applyPhotonFlux,
//
// Transfer some of the settings to EventGeneratorSettings:
//
EventGeneratorSettings* egsettings = EventGeneratorSettings::instance();
egsettings->setDipoleModelParameterSet(settings->dipoleModelParameterSet());
egsettings->setDipoleModelType(settings->dipoleModelType());
egsettings->setA(settings->A());
egsettings->setVerboseLevel(settings->verboseLevel());
egsettings->setApplyPhotonFlux(false);
int nBinbeta = settings->betabins();
int nBinQ2 = settings->Q2bins();
int nBinW2 = settings->W2bins();
int nBinz = settings->zbins();
double betamin = settings->betamin();
double betamax = settings->betamax();
double Q2min = settings->Q2min();
double Q2max = settings->Q2max();
double Wmin = settings->Wmin();
double Wmax = settings->Wmax();
double W2min = Wmin*Wmin;
double W2max = Wmax*Wmax;
double zmin = settings->zmin();
double zmax = settings->zmax();
unsigned int massA = settings->A();
DipoleModelType model = settings->dipoleModelType();
DipoleModelParameterSet pset = settings->dipoleModelParameterSet();
int startingBin = settings->startBin();
int endingBin = settings->endBin();
// int modes = settings->modesToCalculate();
unsigned char priority = settings->priority();
bool hotspots = settings->hasSubstructure();
double fractionOfBinsToFill=settings->fractionOfBinsToFill();
settings->list();
//
// Check bins
//
// Table's bin indices run from 0 ... nbins-1
//
int maxbins = nBinbeta*nBinQ2*nBinW2*nBinz;
if (endingBin > maxbins) {
cout << "Warning, given end bin (" << endingBin << ") exceeds the size of the table." << endl;
cout << " set to maximum value (" << maxbins << ") now." << endl;
endingBin = maxbins;
}
//
// Define all tables. Depending on tableset type some
// will not be written but we define them all anyway.
//
Table tableTQQ[4];
Table tableLQQ[4];
Table tableTQQG[4];
- bool logbeta=false, logQ2=true, logW2=true, logz=false, logC=true;
+ bool logbeta=false, logQ2=false, logW2=false, logz=false, logC=false;
//
// Set filenames for the tables.
// Be as desciptive as possible. Helps when mass producing tables.
//
// We create all tables and decide later what
// gets written and what not.
//
string rootfile=settings->rootfile();
rootfile += "_" + settings->dipoleModelName();
rootfile += "_" + settings->dipoleModelParameterSetName();
ostringstream filenameTQQ[4], filenameLQQ[4], filenameTQQG[4];
for(int i=0; i<4; i++){
filenameTQQ[i].str("");
filenameTQQ[i] << rootfile << "_q" << i << "_bin"
<< startingBin << "-" << endingBin <<"_TQQ.root";
filenameLQQ[i].str("");
filenameLQQ[i] << rootfile << "_q" << i << "_bin"
<< startingBin << "-" << endingBin << "_LQQ.root";
filenameTQQG[i].str("");
filenameTQQG[i] << rootfile << "_q" << i << "_bin"
<< startingBin << "-" << endingBin << "_TQQG.root";
(void) tableTQQ[i].create(nBinbeta, betamin, betamax,
nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinz, zmin, zmax,
logbeta, logQ2, logW2, logz, logC, // all bools
mean_A, transverse, QQ,
massA, model, pset, i,
filenameTQQ[i].str().c_str(),
priority, hotspots);
(void) tableLQQ[i].create(nBinbeta, betamin, betamax,
nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinz, zmin, zmax,
logbeta, logQ2, logW2, logz, logC, // all bools
mean_A, longitudinal, QQ,
massA, model, pset, i,
filenameLQQ[i].str().c_str(),
priority, hotspots);
(void) tableTQQG[i].create(nBinbeta, betamin, betamax,
nBinQ2, Q2min, Q2max,
nBinW2, W2min, W2max,
nBinz, zmin, zmax,
logbeta, logQ2, logW2, logz, logC, // all bools
mean_A, transverse, QQG,
massA, model, pset, i,
filenameTQQG[i].str().c_str(),
priority, hotspots);
}
cout << "\nAll 3x4 tables created, here for up quarks:" << endl;
tableTQQ[0].list();
tableLQQ[0].list();
tableTQQG[0].list();
cout << "\nTables have " << maxbins << " bins each.\n" << endl;
// if (settings->useBackupFile()) {
// int ibin = settings->startingBinFromBackup();
// tableTQQ.setAutobackup("tableTQQ", ibin);
// tableLQQ.setAutobackup("tableLQQ", ibin);
// tableTQQG.setAutobackup("tableTQQG", ibin);
// cout << "Automatic backup of tables is enabled." << endl;
// }
// else
// cout << "Automatic backup of tables is off.\n" << endl;
//
// Setup random generator for creating partial tables.
//
TRandom3 *random = TableGeneratorSettings::randomGenerator();
random->SetSeed();
//
// DGLAP Evolution can be speed up by using lookup tables
//
DglapEvolution &dglap = DglapEvolution::instance();
dglap.generateLookupTable(1000, 1000);
dglap.useLookupTable(true);
//
// Print out settings
//
cout << endl;
cout << "Tables will be generated for:" << endl;
cout << "\tNucleus mass A="<<massA<<endl;
cout << "\tBins: " << startingBin << "-" << endingBin << endl;
cout << "\tbeta range: [" << betamin << ", " << betamax << "], " << nBinbeta << " bins." << endl;
cout << "\tQ2 range: [" << Q2min << ", " << Q2max << "], " << nBinQ2 << " bins." << endl;
cout << "\tW2 range: [" << W2min << ", " << W2max << "], " << nBinW2 << " bins." << endl;
cout << "\t z range: [" << zmin << ", " << zmax << "], " << nBinz << " bins." << endl;
cout << "\tDipole model: " << settings->dipoleModelName() << endl;
cout << "\tDipole model parameter set: " << settings->dipoleModelParameterSetName() << endl;
cout << "\tTable set mode: " << settings->tableSetTypeName() << endl;
cout << endl;
//
// Create and initialize the Inclusive Diffaction Cross Sections
//
- InclusiveDiffractiveCrossSections* idcs= new InclusiveDiffractiveCrossSections;
+ InclusiveDiffractiveCrossSections* idcs= new InclusiveDiffractiveCrossSectionsIntegrals;
//
// Calculate contents and fill tables
//
// Note that we fill all tables. What is written
// at the end is another story.
//
int nShow = (endingBin - startingBin)/100;
if(nShow==0) nShow=1;
- double binWidthBeta=(betamax-betamin)/nBinbeta;
- double binWidthQ2=(Q2max-Q2min)/nBinQ2;
- double binWidthW2=(W2max-W2min)/nBinW2;
- double binWidthZ=(zmax-zmin)/nBinz;
+ double binWidthBeta=(betamax-betamin)/(nBinbeta-1);
+ double binWidthQ2=(Q2max-Q2min)/(nBinQ2-1);
+ double binWidthW2=(W2max-W2min)/(nBinW2-1);
+ double binWidthZ=(zmax-zmin)/(nBinz-1);
// Crude way to find the start and finish bins:
// Loop over all the bins
int iBetaStart=0, iQ2Start=0, iW2Start=0, iZStart=0;
int iBetaEnd=0, iQ2End=0, iW2End=0, iZEnd=0;
bool endbin=false, startbin=false;
- for(int ibeta=1; ibeta<=nBinbeta; ibeta++){
- for(int iQ2=1; iQ2<=nBinQ2; iQ2++){
- for(int iW2=1; iW2<=nBinW2; iW2++){
- for(int iz=1; iz<=nBinz; iz++){
+ for(int ibeta=1; ibeta<nBinbeta; ibeta++){
+ for(int iQ2=1; iQ2<nBinQ2; iQ2++){
+ for(int iW2=1; iW2<nBinW2; iW2++){
+ for(int iz=1; iz<nBinz; iz++){
int currentBin=ibeta*iQ2*iW2*iz;
if(currentBin==startingBin){
iBetaStart=ibeta;
iQ2Start=iQ2;
iW2Start=iW2;
iZStart=ibeta;
startbin=true;
}
if(currentBin==endingBin){
iBetaEnd=ibeta;
iQ2End=iQ2;
iW2End=iW2;
iZEnd=ibeta;
endbin=true;
}
}
}
}
}
if(!endbin){
iBetaEnd=nBinbeta;
iQ2End=nBinQ2;
iW2End=nBinW2;
iZEnd=nBinz;
}
if(!startbin){
- iBetaStart=1;
- iQ2Start=1;
- iW2Start=1;
- iZStart=1;
+ iBetaStart=0;
+ iQ2Start=0;
+ iW2Start=0;
+ iZStart=0;
}
cout << "Start filling tables:" << endl;
time_t tableFillStart = time(0);
//Loop over all the bins
int ibin=0;
- for(int ibeta=iBetaStart; ibeta<=iBetaEnd; ibeta++){
- double beta=betamin+(ibeta-0.5)*binWidthBeta;
- for(int iQ2=iQ2Start; iQ2<=iQ2End; iQ2++){
- double Q2=Q2min+(iQ2-0.5)*binWidthQ2;
- for(int iW2=iW2Start; iW2<=iW2End; iW2++){
- double W2=W2min+(iW2-0.5)*binWidthW2;
- for(int iz=iZStart; iz<=iZEnd; iz++){
+ for(int ibeta=iBetaStart; ibeta<iBetaEnd; ibeta++){
+ double beta=betamin+ibeta*binWidthBeta;
+ for(int iQ2=iQ2Start; iQ2<iQ2End; iQ2++){
+ double Q2=Q2min+iQ2*binWidthQ2;
+ for(int iW2=iW2Start; iW2<iW2End; iW2++){
+ double W2=W2min+iW2*binWidthW2;
+ for(int iz=iZStart; iz<iZEnd; iz++){
+ double z=zmin+iz*binWidthZ;
if (ibin%nShow == 0 || ibin == startingBin || ibin == endingBin)
cout << "processing bin " << ibin << endl;
- double z=zmin+(iz-0.5)*binWidthZ;
-
- if(fractionOfBinsToFill <= random->Uniform()) //Choose fractionOfBinsToFill of the bins;
+ //Choose fractionOfBinsToFill of the bins:
+ if(fractionOfBinsToFill <= random->Uniform())
continue;
idcs->dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, transverse);
idcs->dsigdbetadQ2dW2dz_total(beta, Q2, W2, z, longitudinal);
idcs->dsigdbetadQ2dW2dz_qqg(beta, Q2, W2, z);
//Loop over quark species
for(int iq=0; iq<4; iq++){
- tableTQQ[iq].fill(ibeta, iQ2, iW2, iz, idcs->mDsigdbetadQ2dWdz_T_total[iq]);
- tableLQQ[iq].fill(ibeta, iQ2, iW2, iz, idcs->mDsigdbetadQ2dWdz_L_total[iq]);
- tableTQQG[iq].fill(ibeta, iQ2, iW2, iz, idcs->mDsigdbetadQ2dWdz_T_qqg[iq]);
+ double valTQQ = idcs->dsigdbetadQ2dWdz_T_total()[iq];//nb
+ double valLQQ = idcs->dsigdbetadQ2dWdz_L_total()[iq];//nb
+ double valQQG = idcs->dsigdbetadQ2dWdz_T_qqg()[iq];//nb
+ tableTQQ[iq].fill(ibeta, iQ2, iW2, iz, valTQQ);
+ tableLQQ[iq].fill(ibeta, iQ2, iW2, iz, valLQQ);
+ tableTQQG[iq].fill(ibeta, iQ2, iW2, iz, valQQG);
}
ibin++;
}
}
}
}
-
time_t tableFillEnd = time(0);
//
// Report CPU time/cell
//
cout << endl;
cout << "CPU time/bin: "
<< static_cast<double>(tableFillEnd-tableFillStart)/(endingBin-startingBin+1)
<< " s" << endl;
cout << "Total time: " << static_cast<double>(tableFillEnd-tableFillStart)/60./60. << " h" << endl << endl;
//
// We write out all tables.
//
// Whoever runs the production can then decide
// later what to keep and what to delete.
// If the desired run mode is total_and_coherent or
// coherent_and_incoherent, if this is just to improve
// a coherent table in some phase space, or if one wants
// to maintain redundancy, all these factor might affect
// your choice.
//
for(int i=0; i<4; i++){
tableTQQ[i].write();
tableLQQ[i].write();
tableTQQG[i].write();
}
cout << "All tables written" << endl;
cout << "All done. Bye." << endl;
return 0;
}
Index: trunk/examples/tableGeneratorInclusiveDiffractionRuncard.txt
===================================================================
--- trunk/examples/tableGeneratorInclusiveDiffractionRuncard.txt (revision 540)
+++ trunk/examples/tableGeneratorInclusiveDiffractionRuncard.txt (revision 541)
@@ -1,134 +1,134 @@
//==============================================================================
// tableGeneratorRuncard.txt
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Tobias Toll
// Last update:
// $Date: 2019-03-09 00:43:19 +0530 (Sat, 09 Mar 2019) $
// $Author: ullrich $
//==============================================================================
//
// Example Runcard for Sartre Table Generator.
//
// Comments start with a # or #
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Nucleus number
#
-A = 1
+A = 208
#
# Which moments of the amplitude to calculate:
# 0 <A> analytical, <A2> averaged over configurations
# 1 <A> only, analytical
# 2 <A2>, and <A> both averaged over configurations
# Dummy for inc. diff.
modesToCalculate = 2
#
# Rootfile This is a prefix table type (T and L, or T2 and L2) and ".root" will be appended
#
rootfile = incdifftest
#
# Model parameters
#
#vectorMesonId = 443
dipoleModel = bSat
dipoleModelParameterSet = STU
#
# Table Set Type (experts only)
#
tableSetType = total_and_coherent
#tableSetType = coherent_and_incoherent
#
# Developer only
# Only takes affect when dipoleModelParameterSet = CUSTOM
# All parameters must be given.
#
# For bSat and bNonSat (8 parameters):
# m_u m_d m_s m_c BG mu02 lambdaG Ag
# For bCGC (10 parameters):
# m_u m_d m_s m_c dkappa dN0 dx0 dlambda dgammas dBcgc
#
#dipoleModelCustomParameters = 0.14 0.14 0.14 1.27 4. 1.51 0.058 2.308
# Path for the b-dependence lookup table for bSat. Not used with bCGC.
bSatLookupPath = ./
numberOfConfigurations = 1
#Whether or not there is subnucleon structure (hotspot model)
hasSubstructure = false
#
# Kinematic range
#
# EIC limits:
# Say we aim s = 20000 (30000 for pp)
# That implies:
# Q2 = 0.0001 - 20000 (or with x<1e-2, Q2<200)
# t = -3 - 0
# W2 = 1 - 20000 (W = 1-141)
#
-betamin = 0
-betamax = 1
-Q2min = 1.
-Q2max = 20.
-Wmin = 2.
-Wmax = 380.
-zmin = 0.
-zmax = 1.
+betamin = 0.1
+betamax = 0.9
+Q2min = 2.
+Q2max = 201.
+Wmin = 20.
+Wmax = 240.
+zmin = 0.05
+zmax = .95
#
# Number of bins in the tables:
#
#
# ZEUS test, bin density for
# 40x40x40: 64000/(99.85*190*1)=3.37 bins/GeV^5
#
# Or in different directions:
# Q2: 40/100=0.4 bins/GeV^2
# W: 40/190=0.21 bins/GeV
# t: 40/1= 40 bins/GeV^2
#
fractionOfBinsToFill = 1.
-betabins = 10
+betabins = 10
Q2bins = 10
W2bins = 10
zbins = 10
#
# If a run fails, use the backup and continue:
# If in doubt use the tableInspector to make sure the rest
# of the runcard is as it should be.
#
useBackupFile = false
startingBinFromBackup=2
#
# If the table you are making should take presendence
# over previous tables, make it with a higher priority
#
priority = 0
Index: trunk/examples/sartreInclusiveDiffractionRuncard.txt
===================================================================
--- trunk/examples/sartreInclusiveDiffractionRuncard.txt (revision 540)
+++ trunk/examples/sartreInclusiveDiffractionRuncard.txt (revision 541)
@@ -1,108 +1,118 @@
//==============================================================================
// sartreRuncard.txt
//
// Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation.
// This program is distributed in the hope that it will be useful,
// but without any warranty; without even the implied warranty of
// merchantability or fitness for a particular purpose. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// Author: Thomas Ullrich
// Last update:
// $Date: 2019-06-10 23:53:14 +0530 (Mon, 10 Jun 2019) $
// $Author: ullrich $
//==============================================================================
//
// Example Runcard for Sartre Event Generator.
//
// Comments start with a # or //
//
// Name and value are separated by a "=": name = value
// (alternatively ":" can be used as well)
//==============================================================================
#
# Define beams
#
eBeamEnergy = 30
hBeamEnergy = 920
A = 1
#
# Number of events and printout frequency
#
-numberOfEvents = 10
-timesToShow = 10
+numberOfEvents = 100000
+timesToShow = 100
#
# Set verbosity
#
verbose = false
verboseLevel = 2
#
# Rootfile
#
-rootfile = example.root
+rootfile = example_tables.root
#
# Dipole model
#
# dipoleModel: can be bSat or bNonSat
# dipoleModelParameterSet:
# KMW: older set but most complete
# HMPZ: new parameters with still limited set
# of tables. For eA use KMW.
# STU: new parameters with still limited set
# (uses large dipole damping)
#
dipoleModel = bSat
dipoleModelParameterSet = STU
#
# User variable
# If you use the provided sartreMain.cpp this is
# used for vector meson decays
# PDG: pi+ = 211, K+ = 321, e- = 11, mu- = 13
# with BR=100%
#
userInt = 11
#
# Kinematic range min > max means no limits (given by table range)
#
-Q2min = 10
-Q2max = -11
-Wmin = 10.
-Wmax = -10.049876
-betamin = 0.4
-betamax = 1
-
+#Q2min = 10
+#Q2max = -11
+#Wmin = 10.
+#Wmax = -10.049876
+#betamin = 0.1
+#betamax = 1
+
+betamin = 0.1
+betamax = 0.9
+Q2min = 2.
+Q2max = 201.
+Wmin = 20.
+Wmax = 240.
+zmin = 0.
+zmax = 1.
+
+useInclusiveCrossSectionTables = true
#
# Misc
#
enableNuclearBreakup = false
maxNuclearExcitationEnergy = 3.0
#
# Random generator seed (if not given current time is used)
#
-#seed = 2011987
+seed = 0 #2011987
#
# User parameters
#
# userDouble = 0.
# userString = "Hello World!"
# userInt = 0
#
# Expert flags
#
-# applyPhotonFlux = true
+applyPhotonFlux = true

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:38 PM (19 h, 31 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4014886
Default Alt Text
(25 KB)

Event Timeline