Index: trunk/examples/tableGeneratorInclusiveDiffractionMain.cpp
===================================================================
--- trunk/examples/tableGeneratorInclusiveDiffractionMain.cpp (revision 535)
+++ trunk/examples/tableGeneratorInclusiveDiffractionMain.cpp (revision 536)
@@ -1,414 +1,414 @@
//==============================================================================
// 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;
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();
+// 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;
//
// 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,
+ 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,
+ 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,
+ 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;
//
// 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;
// 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++){
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;
}
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++){
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;
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]);
}
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;
}