Page MenuHomeHEPForge

No OneTemporary

Index: trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp
===================================================================
--- trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 541)
+++ trunk/src/InclusiveDiffractiveCrossSectionsFromTables.cpp (revision 542)
@@ -1,276 +1,309 @@
#include "InclusiveDiffractiveCrossSectionsFromTables.h"
#include <TFile.h>
#include "Math/SpecFunc.h"
#include "TFile.h"
#include <iostream>
#include <cmath>
#include <algorithm>
#include <algorithm>
#include "Constants.h"
#include "Nucleus.h"
#include "DipoleModel.h"
#include "AlphaStrong.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/Functor.h"
#include "TMath.h"
#include "TableGeneratorSettings.h"
#include "Settings.h"
#include "Enumerations.h"
#include "DglapEvolution.h"
#include "Math/WrappedTF1.h"
#include "Math/GaussIntegrator.h"
#include "Kinematics.h"
#include "linterp.h"
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
InclusiveDiffractionCrossSectionsFromTables::InclusiveDiffractionCrossSectionsFromTables(){
+ string path="/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/";
+ /*
+ TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
+
+ TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
- TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
- TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
- TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
- TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
-
- TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
- TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
- TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
- TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_A208_dec15/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
+*/
+ //17x17x17 A1:
+ TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q0_bin0-83521_TQQ.root", "READ");
+ TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q1_bin0-83521_TQQ.root", "READ");
+ TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q2_bin0-83521_TQQ.root", "READ");
+ TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q3_bin0-83521_tQQ.root", "READ");
+
+ TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q0_bin0-83521_LQQ.root", "READ");
+ TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q1_bin0-83521_LQQ.root", "READ");
+ TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q2_bin0-83521_LQQ.root", "READ");
+ TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q3_bin0-83521_LQQ.root", "READ");
- TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
- TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
- TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
- TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_linear_10x10x10_A208/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q0_bin0-83521_TQQG.root", "READ");
+ TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q1_bin0-83521_TQQG.root", "READ");
+ TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q2_bin0-83521_TQQG.root", "READ");
+ TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_17x17x17x17_dec14/incdifftest_bSat_STU_q3_bin0-83521_TQQG.root", "READ");
+ /*
+ // 10x10x10x10:
+ TFile* fileQQ_T0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q0_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q1_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q2_bin0-10000_TQQ.root", "READ");
+ TFile* fileQQ_T3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q3_bin0-10000_TQQ.root", "READ");
+
+ TFile* fileQQ_L0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q0_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q1_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q2_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQ_L3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q3_bin0-10000_LQQ.root", "READ");
+ TFile* fileQQG_0 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q0_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_1 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q1_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_2 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q2_bin0-10000_TQQG.root", "READ");
+ TFile* fileQQG_3 = TFile::Open("/Users/tbu/sartre_SVN/sartre/examples/inclusiveTables_10x10x10_dec14/incdifftest_bSat_STU_q3_bin0-10000_TQQG.root", "READ");
+ */
tableQQ_T.clear();
tableQQ_L.clear();
tableQQG.clear();
// Retrieve the histograms
THnF* histQQ_T0 = dynamic_cast<THnF*>(fileQQ_T0->Get("table"));
THnF* histQQ_T1 = dynamic_cast<THnF*>(fileQQ_T1->Get("table"));
THnF* histQQ_T2 = dynamic_cast<THnF*>(fileQQ_T2->Get("table"));
THnF* histQQ_T3 = dynamic_cast<THnF*>(fileQQ_T3->Get("table"));
THnF* histQQ_L0 = dynamic_cast<THnF*>(fileQQ_L0->Get("table"));
THnF* histQQ_L1 = dynamic_cast<THnF*>(fileQQ_L1->Get("table"));
THnF* histQQ_L2 = dynamic_cast<THnF*>(fileQQ_L2->Get("table"));
THnF* histQQ_L3 = dynamic_cast<THnF*>(fileQQ_L3->Get("table"));
THnF* histQQG_0 = dynamic_cast<THnF*>(fileQQG_0->Get("table"));
THnF* histQQG_1 = dynamic_cast<THnF*>(fileQQG_1->Get("table"));
THnF* histQQG_2 = dynamic_cast<THnF*>(fileQQG_2->Get("table"));
THnF* histQQG_3 = dynamic_cast<THnF*>(fileQQG_3->Get("table"));
tableQQ_T.push_back(histQQ_T0);
tableQQ_T.push_back(histQQ_T1);
tableQQ_T.push_back(histQQ_T2);
tableQQ_T.push_back(histQQ_T3);
tableQQ_L.push_back(histQQ_L0);
tableQQ_L.push_back(histQQ_L1);
tableQQ_L.push_back(histQQ_L2);
tableQQ_L.push_back(histQQ_L3);
tableQQG.push_back(histQQG_0);
tableQQG.push_back(histQQG_1);
tableQQG.push_back(histQQG_2);
tableQQG.push_back(histQQG_3);
}
InclusiveDiffractionCrossSectionsFromTables::~InclusiveDiffractionCrossSectionsFromTables(){
}
double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_T(double beta, double Q2, double W2, double z) {
double MX2 = Q2*(1-beta)/beta;
// double mf = Settings::quarkMass(mQuarkIndex);
double 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); //nb
return result;
}
double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_L(double beta, double Q2, double W2, double z){
double MX2 = Q2*(1-beta)/beta;
// double mf = Settings::quarkMass(mQuarkIndex);
double 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;
}
//===================================================
// QQG Calculations
//===================================================
double InclusiveDiffractionCrossSectionsFromTables::dsigmadbetadz_QQG(double beta, double Q2, double W2, double ztilde) {
double result = interpolate(beta, Q2, W2, ztilde, transverse, QQG, mQuarkIndex);
return result;
}
double InclusiveDiffractionCrossSectionsFromTables::interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark){
THnF* hist4D=0;
//#TT TODO: There seem to be a strange swap happening
// between T QQ and L QQ...
if(pol==longitudinal and fock==QQ){
hist4D=tableQQ_T.at(quark);
}
else if(pol==transverse and fock==QQ){
hist4D=tableQQ_L.at(quark);
}
else if(fock==QQG){
hist4D=tableQQG.at(quark);
}
else{
cout<<"InclusiveDiffractionCrossSectionsFromTables::interpolate(double beta, double Q2, double W2, double z, GammaPolarization pol, FockState fock, int quark): There is nothing to interpolate, ending"<<endl;
exit(2);
}
int nDims = hist4D->GetNdimensions();
-// PR(nDims);
-
vector <double>MinRange(nDims); //beta, Q2,W2, z
vector <double>MaxRange(nDims); //beta, Q2,W2, z
vector <int>Bins(nDims); //beta, Q2,W2, z
for (int dim = 0; dim < 4; dim++) {
TAxis* axis = hist4D->GetAxis(dim); // Get the axis for this dimension
if (!axis) {
std::cerr << "Error: Cannot retrieve axis for dimension " << dim << std::endl;
continue;
}
// Get the number of bins and range
int nBins = axis->GetNbins();
double minRange = axis->GetXmin();
double maxRange = axis->GetXmax();
Bins[dim]= nBins;
MinRange[dim]= minRange;
MaxRange[dim]= maxRange;
}
// construct the grid in each dimension.
// note that we will pass in a sequence of iterators pointing to the beginning of each grid
std::vector<double> gridBeta = mgrid(hist4D, 0, Bins[0]);
std::vector<double> gridQ2 = mgrid(hist4D, 1, Bins[1]) ;
std::vector<double> gridW2 = mgrid(hist4D, 2, Bins[2]);
std::vector<double> gridZ = mgrid(hist4D, 3, Bins[3]) ;
std::vector< std::vector<double>::iterator > grid_iter_list;
grid_iter_list.push_back(gridBeta.begin());
grid_iter_list.push_back(gridQ2.begin());
grid_iter_list.push_back(gridW2.begin());
grid_iter_list.push_back(gridZ.begin());
// the size of the grid in each dimension
array<int,4> grid_sizes;
for(unsigned long i=0; i<grid_sizes.size(); i++){
grid_sizes[i] = Bins[i];
}
// total number of elements
int num_elements = grid_sizes[0] * grid_sizes[1]* grid_sizes[2]* grid_sizes[3];
// fill in the values of f(x) at the gridpoints.
// we will pass in a contiguous sequence, values are assumed to be laid out C-style
std::vector<double> f_values(num_elements);
int count = 0;
- int countExtreme=0;
+ int countFail=0;
for (int ibeta=1; ibeta<=grid_sizes[0]; ibeta++) {
for (int iQ2=1; iQ2<=grid_sizes[1]; iQ2++) {
for (int iW2=1; iW2<=grid_sizes[2]; iW2++) {
for (int iZ=1; iZ<=grid_sizes[3]; iZ++) {
int binIndex[4] = {ibeta-1, iQ2-1, iW2-1, iZ-1};
double val=hist4D->GetBinContent(binIndex);
- if(exp(val)>1e20){ //There seems to be some spikes
- countExtreme++;
+ // There seems to be some spikes
+ // This has to be adapted to lin/log content
+ // if(exp(val)>1e20){
+ // val=0;
+ // }
+ if(isnan(val) or isinf(val)){
+ countFail++;
val=0;
}
if(val!=0)
f_values[count] = val;//exp(val);
else
f_values[count] = 0.;
count++;
+ valOld=val;
}
}
}
}
-// if(countExtreme){
-// PR(countExtreme);
-//// exit(3);
-// }
//
// Construct the interpolator.
// The last two arguments are pointers to the underlying data
//
InterpMultilinear <4, double> interp_ML(grid_iter_list.begin(), grid_sizes.begin(), f_values.data(), f_values.data() + num_elements);
// interpolate:
array<double,4> args = {beta, Q2, W2, z};
double result=interp_ML.interp(args.begin());
-// return result;
+// return result; //For linear content
if(result==0)
return 0.;
else
- return exp(result);
+ return exp(result); //For log content
}
vector<double> InclusiveDiffractionCrossSectionsFromTables::mgrid(THnF* hist,int axis, int len){
vector<double> gridtype(len);
for(int i = 1; i <=len; i++){
double val = hist->GetAxis(axis)->GetBinCenter(i);
gridtype[i-1] = val;
}
return gridtype;
}

File Metadata

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

Event Timeline