Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8308865
Table.cpp
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
26 KB
Subscribers
None
Table.cpp
View Options
//==============================================================================
// Table.cpp
//
// Copyright (C) 2010-2013 Tobias Toll and Thomas Ullrich
//
// This file is part of Sartre version: 1.1
//
// 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: 2013-05-29 21:26:19 +0100 (Wed, 29 May 2013) $
// $Author: thomas.ullrich@bnl.gov $
//==============================================================================
//
// Table data is stored in a TH3F.
//
// All information is stored in the histogram title.
//
// Table ID is encoded as follows:
//
// histogram title is a number that is to be interpreted
// as an uint64_t with the bits set as follows:
//
// bit 0: content type: 0 for <A>, 1 for <A2>
// bit 1: polarization: T=1, L=0
// bit 2: t encoding: 0 for |t|, 1 for log(|t|)
// bit 3: W2 encoding: 0 for lin, 1 for log
// bit 4: Q2 encoding: 0 for lin, 1 for log
// bit 5-7: dipole model type (Enumerations.h)
// bit 8-15: mass number A
// bit 16-31: vector meson ID (PDG)
// bit 32: content encoding: 0 in lin, 1 in log
// bit 33: content type is lambda_<A> (bit 0 == 0 in this case)
// bit 34-63: not used
//
// log implies ln
// bit 0 is the LSB
//
// Internal histogram: x = Q2, y = W2, z=t (or the logs)
//
// Actual value is taken at the center of the bin.
// min/max functions (e.g. minQ2()) return the value of the first/last
// bin center.
//==============================================================================
#include
"Table.h"
#include
"TFile.h"
#include
"TH3F.h"
#include
<cstdio>
#include
<cmath>
#include
<ctime>
#include
<algorithm>
#include
<string>
#include
<sstream>
#include
<fstream>
#include
<iomanip>
#include
<limits>
#include
<unistd.h>
#include
<sys/types.h>
#include
"GridSpline.h"
#define PR(x) cout << #x << " = " << (x) << endl;
Table
::
Table
()
{
mID
=
0
;
m3DHisto
=
0
;
mFillCounter
=
0
;
mBackupFrequence
=
0
;
}
Table
&
Table
::
operator
=
(
const
Table
&
tab
)
{
if
(
this
!=
&
tab
)
{
delete
m3DHisto
;
mID
=
tab
.
mID
;
m3DHisto
=
new
TH3F
(
*
(
tab
.
m3DHisto
));
m3DHisto
->
SetDirectory
(
0
);
mFilename
=
tab
.
mID
;;
mFillCounter
=
tab
.
mID
;;
mBackupFrequence
=
tab
.
mID
;;
mBackupPrefix
=
tab
.
mID
;;
mLastBackupFilename
=
tab
.
mID
;;
}
return
*
this
;
}
Table
::
Table
(
const
Table
&
tab
)
{
mID
=
tab
.
mID
;
m3DHisto
=
new
TH3F
(
*
(
tab
.
m3DHisto
));
m3DHisto
->
SetDirectory
(
0
);
mFilename
=
tab
.
mID
;;
mFillCounter
=
tab
.
mID
;;
mBackupFrequence
=
tab
.
mID
;;
mBackupPrefix
=
tab
.
mID
;;
mLastBackupFilename
=
tab
.
mID
;;
}
Table
::~
Table
()
{
delete
m3DHisto
;
}
unsigned
int
Table
::
create
(
int
nbinsQ2
,
double
Q2min
,
double
Q2max
,
int
nbinsW2
,
double
W2min
,
double
W2max
,
int
nbinsT
,
double
tmin
,
double
tmax
,
bool
logQ2
,
bool
logW2
,
bool
logt
,
bool
logC
,
AmplitudeMoment
mom
,
GammaPolarization
pol
,
unsigned
int
A
,
int
vm
,
DipoleModelType
model
,
const
char
*
filename
)
{
if
(
m3DHisto
!=
0
)
{
cout
<<
"Table:create(): cannot create, table already exists."
<<
endl
;
return
0
;
}
if
(
filename
)
mFilename
=
string
(
filename
);
// name of file where table is written to
//
// Create table ID first
//
mID
=
(
vm
<<
16
);
mID
|=
(
A
<<
8
);
mID
|=
(
model
<<
5
);
if
(
mom
==
mean_A2
)
mID
|=
1
;
if
(
pol
==
transverse
)
mID
|=
(
1
<<
1
);
if
(
logt
)
mID
|=
(
1
<<
2
);
if
(
logW2
)
mID
|=
(
1
<<
3
);
if
(
logQ2
)
mID
|=
(
1
<<
4
);
if
(
logC
)
mID
|=
(
static_cast
<
uint64_t
>
(
1
)
<<
32
);
if
(
mom
==
lambda_A
)
mID
|=
(
static_cast
<
uint64_t
>
(
1
)
<<
33
);
ostringstream
titlestream
;
titlestream
<<
mID
;
string
title
=
titlestream
.
str
();
//
// Binning
//
// Interpolate() only works up to the bin center
// of the first and last bin. To guarantee that
// the full range is accessible we shift the min
// and max of each axis.
//
tmin
=
fabs
(
tmin
);
tmax
=
fabs
(
tmax
);
if
(
tmin
>
tmax
)
swap
(
tmin
,
tmax
);
if
(
logQ2
)
Q2min
=
log
(
Q2min
);
if
(
logQ2
)
Q2max
=
log
(
Q2max
);
if
(
logW2
)
W2min
=
log
(
W2min
);
if
(
logW2
)
W2max
=
log
(
W2max
);
if
(
logt
)
tmin
=
log
(
tmin
);
if
(
logt
)
tmax
=
log
(
tmax
);
double
binWidthQ2
=
(
Q2max
-
Q2min
)
/
(
nbinsQ2
-
1
);
double
binWidthW2
=
(
W2max
-
W2min
)
/
(
nbinsW2
-
1
);
double
binWidthT
=
(
tmax
-
tmin
)
/
(
nbinsT
-
1
);
//
// Book 3D histo to hold table
//
m3DHisto
=
new
TH3F
(
"table"
,
title
.
c_str
(),
nbinsQ2
,
Q2min
-
binWidthQ2
/
2.
,
Q2max
+
binWidthQ2
/
2.
,
nbinsW2
,
W2min
-
binWidthW2
/
2.
,
W2max
+
binWidthW2
/
2.
,
nbinsT
,
tmin
-
binWidthT
/
2.
,
tmax
+
binWidthT
/
2.
);
m3DHisto
->
SetDirectory
(
0
);
return
nbinsQ2
*
nbinsW2
*
nbinsT
;
// return total number of bins
}
int
Table
::
numberOfEntries
()
const
{
if
(
m3DHisto
)
{
int
nx
=
m3DHisto
->
GetXaxis
()
->
GetNbins
();
int
ny
=
m3DHisto
->
GetYaxis
()
->
GetNbins
();
int
nz
=
m3DHisto
->
GetZaxis
()
->
GetNbins
();
return
nx
*
ny
*
nz
;
}
else
return
0
;
}
void
Table
::
binXYZ
(
int
globalBin
,
int
&
binx
,
int
&
biny
,
int
&
binz
)
const
{
//
// Find correct bins for each axis.
// Don't use ROOT global bins here,
// they are a mess since they cross
// over underflow and overflow bins.
//
// All bins returned count starting
// at 1. The global bin starts at
// 0 as usual in C++.
//
int
nx
=
m3DHisto
->
GetXaxis
()
->
GetNbins
();
int
ny
=
m3DHisto
->
GetYaxis
()
->
GetNbins
();
binx
=
globalBin
%
nx
;
biny
=
((
globalBin
-
binx
)
/
nx
)
%
ny
;
binz
=
((
globalBin
-
binx
)
/
nx
-
biny
)
/
ny
;
binx
++
;
biny
++
;
binz
++
;
}
void
Table
::
binCenter
(
int
i
,
double
&
Q2
,
double
&
W2
,
double
&
t
)
const
{
int
binx
,
biny
,
binz
;
binXYZ
(
i
,
binx
,
biny
,
binz
);
// cout << i << '\t' << binx << '\t' << biny << '\t' << binz << endl;
double
x
=
m3DHisto
->
GetXaxis
()
->
GetBinCenter
(
binx
);
double
y
=
m3DHisto
->
GetYaxis
()
->
GetBinCenter
(
biny
);
double
z
=
m3DHisto
->
GetZaxis
()
->
GetBinCenter
(
binz
);
Q2
=
isLogQ2
()
?
exp
(
x
)
:
x
;
W2
=
isLogW2
()
?
exp
(
y
)
:
y
;
t
=
isLogT
()
?
-
exp
(
z
)
:
-
z
;
if
(
t
>
0
)
t
=
0
;
// avoid numeric glitch
}
void
Table
::
fill
(
int
i
,
double
val
,
double
err
)
{
int
binx
,
biny
,
binz
;
binXYZ
(
i
,
binx
,
biny
,
binz
);
if
(
isLogContent
())
{
if
(
val
==
0
)
{
cout
<<
"Table::fill(): warning, log scale requested but value = 0."
<<
endl
;
val
=
numeric_limits
<
float
>::
min
();
}
val
=
log
(
fabs
(
val
));
}
m3DHisto
->
SetBinContent
(
binx
,
biny
,
binz
,
val
);
if
(
err
!=
0.
)
m3DHisto
->
SetBinError
(
binx
,
biny
,
binz
,
err
);
mFillCounter
++
;
//
// Check if backup is due
//
if
(
mBackupFrequence
)
if
(
mFillCounter
%
mBackupFrequence
==
0
)
backup
(
i
);
}
bool
Table
::
fexists
(
const
char
*
filename
)
const
{
ifstream
ifs
(
filename
);
return
ifs
;
}
void
Table
::
write
(
const
char
*
filename
)
{
//
// 'filename' is optional argument. Null value implies
// that one was already passed via create(). Check
// that this is the case. If one is given, it is used
// no matter if it already was defined or not.
//
if
(
filename
)
mFilename
=
string
(
filename
);
else
{
if
(
mFilename
.
empty
())
{
// should be defined but is not
cout
<<
"Table::write(): Warning, no filename supplied. Will use 'sartre_table.root'."
<<
endl
;
mFilename
=
string
(
"sartre_table.root"
);
}
}
//
// Check if file exist. If so alter
// the filename. Creating tables is
// a time consuming business so we do
// not want to lose precious data if
// something goes wrong here and we do
// want to prevent accidents.
//
while
(
fexists
(
mFilename
.
c_str
()))
{
mFilename
+=
string
(
".save"
);
}
TFile
hfile
(
mFilename
.
c_str
(),
"RECREATE"
);
m3DHisto
->
Write
();
hfile
.
Close
();
cout
<<
"Table::write(): table stored in file '"
<<
mFilename
.
c_str
()
<<
"'."
<<
endl
;
if
(
mBackupFrequence
&&
!
mLastBackupFilename
.
empty
())
remove
(
mLastBackupFilename
.
c_str
());
}
bool
Table
::
read
(
const
char
*
filename
)
{
//
// Read histogram into memory
//
TFile
*
file
=
TFile
::
Open
(
filename
,
"READ"
);
if
(
file
==
0
)
{
cout
<<
"Table::read(): failed opening file '"
<<
filename
<<
"'."
<<
endl
;
return
false
;
}
m3DHisto
=
reinterpret_cast
<
TH3F
*>
(
file
->
Get
(
"table"
));
if
(
m3DHisto
==
0
)
{
cout
<<
"Table::read(): failed retrieving table from file '"
<<
filename
<<
"'."
<<
endl
;
return
false
;
}
m3DHisto
->
SetDirectory
(
0
);
file
->
Close
();
//
// Get table ID
//
mID
=
atoll
(
m3DHisto
->
GetTitle
());
mFilename
=
string
(
filename
);
return
true
;
}
int
Table
::
vectorMesonId
()
const
{
return
((
mID
>>
16
)
&
0xFFFF
);}
unsigned
int
Table
::
dipoleModelType
()
const
{
return
((
mID
>>
5
)
&
0x7
);}
unsigned
int
Table
::
A
()
const
{
return
((
mID
>>
8
)
&
0xFF
);}
bool
Table
::
isLongitudinal
()
const
{
return
!
isTransverse
();}
bool
Table
::
isTransverse
()
const
{
return
(
mID
&
(
1
<<
1
));}
GammaPolarization
Table
::
polarization
()
const
{
return
(
isTransverse
()
?
transverse
:
longitudinal
);}
AmplitudeMoment
Table
::
moment
()
const
{
if
(
isLambdaA
())
return
lambda_A
;
else
if
(
isMeanA
())
return
mean_A
;
else
return
mean_A2
;
}
bool
Table
::
isMeanA
()
const
{
return
!
(
mID
&
1
)
&&
!
(
mID
&
(
static_cast
<
uint64_t
>
(
1
)
<<
33
));}
bool
Table
::
isMeanA2
()
const
{
return
(
mID
&
1
)
&&
!
(
mID
&
(
static_cast
<
uint64_t
>
(
1
)
<<
33
));}
bool
Table
::
isLogQ2
()
const
{
return
(
mID
&
(
1
<<
4
));}
bool
Table
::
isLogW2
()
const
{
return
(
mID
&
(
1
<<
3
));}
bool
Table
::
isLogT
()
const
{
return
(
mID
&
(
1
<<
2
));}
bool
Table
::
isLogContent
()
const
{
return
(
mID
&
(
static_cast
<
uint64_t
>
(
1
)
<<
32
));}
bool
Table
::
isLambdaA
()
const
{
return
!
(
mID
&
1
)
&&
(
mID
&
(
static_cast
<
uint64_t
>
(
1
)
<<
33
));}
uint64_t
Table
::
id
()
const
{
return
mID
;}
double
Table
::
binWidthW2
()
const
{
return
m3DHisto
->
GetYaxis
()
->
GetBinWidth
(
1
);
}
double
Table
::
get
(
double
Q2
,
double
W2
,
double
t
)
const
{
if
(
m3DHisto
==
0
)
return
0
;
//
// Transform variables to how they are stored in the table
//
double
x
=
isLogQ2
()
?
log
(
Q2
)
:
Q2
;
double
y
=
isLogW2
()
?
log
(
W2
)
:
W2
;
t
=
fabs
(
t
);
double
z
=
isLogT
()
?
log
(
t
)
:
t
;
//
// Tiny numerical glitches will cause TH3F::Interpolate() to fail
// since it requires that the variables lie between the first
// and last bin center, excluding the centers.
// Here we enforce that this is the case. The downside of this
// is that *all* values of Q2, W2, t will be forced to lie within.
// Hence the user has to make sure that the values are within
// the boundaries (see minQ2(), maxQ2(), minW2() etc.) before
// calling get().
// In this case the corrections below are tiny and are only
// applied to avoid minor glitches that do not affect the results.
//
TAxis
*
axis
=
m3DHisto
->
GetXaxis
();
x
=
max
(
x
,
axis
->
GetBinCenter
(
1
)
+
numeric_limits
<
float
>::
epsilon
());
x
=
min
(
x
,
axis
->
GetBinCenter
(
axis
->
GetNbins
())
-
numeric_limits
<
float
>::
epsilon
());
axis
=
m3DHisto
->
GetYaxis
();
y
=
max
(
y
,
axis
->
GetBinCenter
(
1
)
+
numeric_limits
<
float
>::
epsilon
());
y
=
min
(
y
,
axis
->
GetBinCenter
(
axis
->
GetNbins
())
-
numeric_limits
<
float
>::
epsilon
());
axis
=
m3DHisto
->
GetZaxis
();
z
=
max
(
z
,
axis
->
GetBinCenter
(
1
)
+
numeric_limits
<
float
>::
epsilon
());
z
=
min
(
z
,
axis
->
GetBinCenter
(
axis
->
GetNbins
())
-
numeric_limits
<
float
>::
epsilon
());
// double result = InterpolateGridSpline(x, y, z); // tmp uncommented until 0's in tables are cleared
double
result
=
m3DHisto
->
Interpolate
(
x
,
y
,
z
);
if
(
result
==
0
&&
isLogContent
())
{
cout
<<
"Table::get(): warning, 0 is not a valid table content when working in log scale."
<<
endl
;
}
if
(
isLogContent
())
result
=
exp
(
result
);
return
result
;
}
double
Table
::
minQ2
()
const
{
if
(
m3DHisto
==
0
)
return
0
;
TAxis
*
axis
=
m3DHisto
->
GetXaxis
();
double
val
=
axis
->
GetBinCenter
(
1
);
return
isLogQ2
()
?
exp
(
val
)
:
val
;
}
double
Table
::
maxQ2
()
const
{
if
(
m3DHisto
==
0
)
return
0
;
TAxis
*
axis
=
m3DHisto
->
GetXaxis
();
double
val
=
axis
->
GetBinCenter
(
axis
->
GetNbins
());
return
isLogQ2
()
?
exp
(
val
)
:
val
;
}
double
Table
::
minW2
()
const
{
if
(
m3DHisto
==
0
)
return
0
;
TAxis
*
axis
=
m3DHisto
->
GetYaxis
();
double
val
=
axis
->
GetBinCenter
(
1
);
return
isLogW2
()
?
exp
(
val
)
:
val
;
}
double
Table
::
maxW2
()
const
{
if
(
m3DHisto
==
0
)
return
0
;
TAxis
*
axis
=
m3DHisto
->
GetYaxis
();
double
val
=
axis
->
GetBinCenter
(
axis
->
GetNbins
());
return
isLogW2
()
?
exp
(
val
)
:
val
;
}
double
Table
::
minT
()
const
{
if
(
m3DHisto
==
0
)
return
0
;
TAxis
*
axis
=
m3DHisto
->
GetZaxis
();
double
val
=
axis
->
GetBinCenter
(
axis
->
GetNbins
());
// t always as |t| in table
return
isLogT
()
?
-
exp
(
val
)
:
-
val
;
}
double
Table
::
maxT
()
const
{
if
(
m3DHisto
==
0
)
return
0
;
TAxis
*
axis
=
m3DHisto
->
GetZaxis
();
double
val
=
axis
->
GetBinCenter
(
1
);
// t always as |t| in table
double
result
=
isLogT
()
?
-
exp
(
val
)
:
-
val
;
if
(
result
>
0
)
{
// cout << "Table::maxT(): warning, t > 0, t (" << result << ") set to 0 now." << endl;
result
=
0
;
}
return
result
;
}
string
Table
::
filename
()
const
{
return
mFilename
;}
void
Table
::
list
(
ostream
&
os
,
bool
printContent
,
bool
printStatistics
)
const
{
ios
::
fmtflags
fmt
=
os
.
flags
();
// store current I/O flags
if
(
!
m3DHisto
)
{
os
<<
"Table::list(): table is undefined."
<<
endl
;
return
;
}
os
<<
setw
(
11
)
<<
"Table ID = "
<<
setw
(
14
)
<<
left
<<
mID
;
os
<<
setw
(
6
)
<<
right
<<
"A = "
<<
setw
(
4
)
<<
left
<<
A
();
os
<<
setw
(
18
)
<<
right
<<
"content = "
<<
(
isLogContent
()
?
"log"
:
""
);
if
(
isMeanA
())
os
<<
"<A>"
<<
endl
;
else
if
(
isMeanA2
())
os
<<
"<A^2>"
<<
endl
;
else
if
(
isLambdaA
())
os
<<
"lambda_<A>"
<<
endl
;
else
os
<<
"unknown"
<<
endl
;
os
<<
setw
(
6
+
11
+
14
)
<<
right
<<
"vmId = "
<<
setw
(
4
)
<<
left
<<
vectorMesonId
();
os
<<
setw
(
18
)
<<
right
<<
"polarization = "
<<
(
isTransverse
()
?
'T'
:
'L'
)
<<
endl
;
os
<<
setw
(
6
+
11
+
14
)
<<
right
<<
"model = "
;
if
(
dipoleModelType
()
==
bSat
)
os
<<
"bSat"
<<
endl
;
else
if
(
dipoleModelType
()
==
bNonSat
)
os
<<
"bNonSat"
<<
endl
;
else
os
<<
"bCGC"
<<
endl
;
os
<<
setw
(
6
+
11
+
14
+
1
)
<<
right
<<
"Q2 = ["
<<
minQ2
()
<<
", "
<<
maxQ2
()
<<
"; bins = "
<<
m3DHisto
->
GetXaxis
()
->
GetNbins
()
<<
(
isLogQ2
()
?
"; log]"
:
"; lin]"
)
<<
endl
;
os
<<
setw
(
6
+
11
+
14
+
1
)
<<
right
<<
"W2 = ["
<<
minW2
()
<<
", "
<<
maxW2
()
<<
"; bins = "
<<
m3DHisto
->
GetYaxis
()
->
GetNbins
()
<<
(
isLogW2
()
?
"; log]"
:
"; lin]"
)
<<
endl
;
os
<<
setw
(
6
+
11
+
14
+
1
)
<<
right
<<
"t = ["
<<
minT
()
<<
", "
<<
maxT
()
<<
"; bins = "
<<
m3DHisto
->
GetZaxis
()
->
GetNbins
()
<<
(
isLogT
()
?
"; log]"
:
"; lin]"
)
<<
endl
;
os
<<
setw
(
6
+
11
+
14
)
<<
right
<<
"file = "
<<
mFilename
.
c_str
()
<<
endl
;
os
<<
endl
;
double
Q2
,
W2
,
t
;
int
binx
,
biny
,
binz
;
if
(
printContent
)
{
for
(
int
i
=
0
;
i
<
numberOfEntries
();
i
++
)
{
binXYZ
(
i
,
binx
,
biny
,
binz
);
binCenter
(
i
,
Q2
,
W2
,
t
);
double
value
=
get
(
Q2
,
W2
,
t
);
double
rawContent
=
m3DHisto
->
GetBinContent
(
binx
,
biny
,
binz
);
double
error
=
m3DHisto
->
GetBinError
(
binx
,
biny
,
binz
);
os
<<
"bin = "
<<
setw
(
8
)
<<
left
<<
i
;
os
<<
"Q2 = "
<<
setw
(
10
)
<<
left
<<
fixed
<<
setprecision
(
3
)
<<
Q2
;
os
<<
"W2 = "
<<
setw
(
11
)
<<
left
<<
fixed
<<
W2
;
os
<<
"t = "
<<
setw
(
14
)
<<
left
<<
scientific
<<
t
;
os
<<
"value = "
<<
setw
(
12
)
<<
left
<<
scientific
<<
value
;
os
<<
"(binx = "
<<
setw
(
5
)
<<
left
<<
binx
;
os
<<
"biny = "
<<
setw
(
5
)
<<
left
<<
biny
;
os
<<
"binz = "
<<
setw
(
5
)
<<
left
<<
binz
;
os
<<
"content = "
<<
setw
(
12
)
<<
left
<<
scientific
<<
rawContent
;
os
<<
"error = "
<<
left
<<
error
<<
')'
;
if
(
(
isLogContent
()
&&
rawContent
<=
log
(
numeric_limits
<
float
>::
min
()
*
2
))
||
(
rawContent
==
0
))
cout
<<
" Z"
;
cout
<<
endl
;
}
os
<<
endl
;
}
if
(
printStatistics
)
{
int
nEmpty
=
0
;
int
nInvalid
=
0
;
int
nNegative
=
0
;
double
sum
=
0
;
double
maxContent
=
0
;
double
minContent
=
numeric_limits
<
float
>::
max
();
int
minBin
,
maxBin
;
minBin
=
maxBin
=
0
;
for
(
int
i
=
0
;
i
<
numberOfEntries
();
i
++
)
{
binXYZ
(
i
,
binx
,
biny
,
binz
);
double
c
=
m3DHisto
->
GetBinContent
(
binx
,
biny
,
binz
);
if
(
isLogContent
())
{
if
(
c
<=
log
(
numeric_limits
<
float
>::
min
()
*
2
))
nEmpty
++
;
}
else
{
if
(
c
==
0
)
nEmpty
++
;
}
if
(
!
isLogContent
()
&&
c
<
0
)
nNegative
++
;
if
(
isnan
(
c
)
||
isinf
(
c
))
nInvalid
++
;
if
(
isLogContent
())
c
=
exp
(
c
);
if
(
c
>
maxContent
)
{
maxContent
=
c
;
maxBin
=
i
;}
if
(
c
>
0
&&
c
<
minContent
)
{
minContent
=
c
;
minBin
=
i
;}
if
(
c
>
0
)
sum
+=
c
;
}
os
<<
setw
(
31
)
<<
right
<<
"total number of cells = "
<<
numberOfEntries
()
<<
endl
;
os
<<
setw
(
31
)
<<
right
<<
"cells with negative content = "
<<
nNegative
<<
endl
;
os
<<
setw
(
31
)
<<
right
<<
"cells with no (0) content = "
<<
nEmpty
<<
endl
;
os
<<
setw
(
31
)
<<
right
<<
"cells with invalid content = "
<<
nInvalid
<<
endl
;
binCenter
(
minBin
,
Q2
,
W2
,
t
);
os
<<
setw
(
31
)
<<
right
<<
"minimum content = "
<<
minContent
<<
" at Q2 = "
<<
Q2
<<
", W2 = "
<<
W2
<<
", t = "
<<
t
<<
endl
;
binCenter
(
maxBin
,
Q2
,
W2
,
t
);
os
<<
setw
(
31
)
<<
right
<<
"maximum content = "
<<
maxContent
<<
" at Q2 = "
<<
Q2
<<
", W2 = "
<<
W2
<<
", t = "
<<
t
<<
endl
;
os
<<
setw
(
31
)
<<
right
<<
"sum = "
<<
sum
<<
endl
;
os
<<
endl
;
}
os
.
flags
(
fmt
);
// restore I/O flags
}
double
Table
::
InterpolateGridSpline
(
double
x
,
double
y
,
double
z
)
const
// Q2, W2, t
{
//
// The algorithm was taken from Cristian Constantin Lalescu.
// See http://arxiv.org/abs/0905.3564 for details.
//
// Grid points: 4, 6, or 8
// Spline order: 3,5 for order 4
// 3,5,7,9 for order 6
// 3,5,7,9,11,13 for order 8
//
//
// Find cell that refer to the position
//
int
nx
=
m3DHisto
->
GetXaxis
()
->
FindBin
(
x
);
int
ny
=
m3DHisto
->
GetYaxis
()
->
FindBin
(
y
);
int
nz
=
m3DHisto
->
GetZaxis
()
->
FindBin
(
z
);
//
// Define the # of grid points depending on how far we are away
// from the edge. If at the edge we fall back to linear interpolation
// as provided by TH3.
// There must be a smarter way doing this than the code below
//
int
gridPoints
=
0
;
if
(
nx
-
3
>=
1
&&
nx
+
4
<=
m3DHisto
->
GetXaxis
()
->
GetNbins
()
&&
ny
-
3
>=
1
&&
ny
+
4
<=
m3DHisto
->
GetYaxis
()
->
GetNbins
()
&&
nz
-
3
>=
1
&&
nz
+
4
<=
m3DHisto
->
GetZaxis
()
->
GetNbins
())
{
gridPoints
=
8
;
}
else
if
(
nx
-
2
>=
1
&&
nx
+
3
<=
m3DHisto
->
GetXaxis
()
->
GetNbins
()
&&
ny
-
2
>=
1
&&
ny
+
3
<=
m3DHisto
->
GetYaxis
()
->
GetNbins
()
&&
nz
-
2
>=
1
&&
nz
+
3
<=
m3DHisto
->
GetZaxis
()
->
GetNbins
())
{
gridPoints
=
6
;
}
else
if
(
nx
-
1
>=
1
&&
nx
+
2
<=
m3DHisto
->
GetXaxis
()
->
GetNbins
()
&&
ny
-
1
>=
1
&&
ny
+
2
<=
m3DHisto
->
GetYaxis
()
->
GetNbins
()
&&
nz
-
1
>=
1
&&
nz
+
2
<=
m3DHisto
->
GetZaxis
()
->
GetNbins
())
{
gridPoints
=
4
;
}
else
{
return
m3DHisto
->
Interpolate
(
x
,
y
,
z
);
}
//
// Find bin centers
//
double
xc
=
m3DHisto
->
GetXaxis
()
->
GetBinCenter
(
nx
);
double
yc
=
m3DHisto
->
GetYaxis
()
->
GetBinCenter
(
ny
);
double
zc
=
m3DHisto
->
GetZaxis
()
->
GetBinCenter
(
nz
);
//
// Define the scale for coordinate transformations
// grid_spline expects x,y,z in bin units
//
double
xscale
=
m3DHisto
->
GetXaxis
()
->
GetBinWidth
(
1
);
double
yscale
=
m3DHisto
->
GetYaxis
()
->
GetBinWidth
(
1
);
double
zscale
=
m3DHisto
->
GetZaxis
()
->
GetBinWidth
(
1
);
//
// Prepare grid spline alogorithm
//
double
result
,
splineOrder
;
if
(
gridPoints
==
8
)
{
local_scal_3D
<
8
>
lf
;
for
(
int
i
=-
3
;
i
<=
4
;
i
++
)
{
for
(
int
j
=-
3
;
j
<=
4
;
j
++
)
{
lf
(
-
3
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
-
3
,
ny
+
i
,
nz
+
j
);
lf
(
-
2
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
-
2
,
ny
+
i
,
nz
+
j
);
lf
(
-
1
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
-
1
,
ny
+
i
,
nz
+
j
);
lf
(
0
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
,
ny
+
i
,
nz
+
j
);
lf
(
1
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
1
,
ny
+
i
,
nz
+
j
);
lf
(
2
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
2
,
ny
+
i
,
nz
+
j
);
lf
(
3
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
3
,
ny
+
i
,
nz
+
j
);
lf
(
4
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
4
,
ny
+
i
,
nz
+
j
);
}
}
splineOrder
=
7
;
result
=
grid_spline
(
splineOrder
,
lf
,
(
x
-
xc
)
/
xscale
,(
y
-
yc
)
/
yscale
,(
z
-
zc
)
/
zscale
);
}
else
if
(
gridPoints
==
6
)
{
local_scal_3D
<
6
>
lf
;
for
(
int
i
=-
2
;
i
<=
3
;
i
++
)
{
for
(
int
j
=-
2
;
j
<=
3
;
j
++
)
{
lf
(
-
2
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
-
2
,
ny
+
i
,
nz
+
j
);
lf
(
-
1
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
-
1
,
ny
+
i
,
nz
+
j
);
lf
(
0
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
,
ny
+
i
,
nz
+
j
);
lf
(
1
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
1
,
ny
+
i
,
nz
+
j
);
lf
(
2
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
2
,
ny
+
i
,
nz
+
j
);
lf
(
3
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
3
,
ny
+
i
,
nz
+
j
);
}
}
splineOrder
=
5
;
result
=
grid_spline
(
splineOrder
,
lf
,
(
x
-
xc
)
/
xscale
,(
y
-
yc
)
/
yscale
,(
z
-
zc
)
/
zscale
);
}
else
if
(
gridPoints
==
4
)
{
local_scal_3D
<
4
>
lf
;
for
(
int
i
=-
1
;
i
<=
2
;
i
++
)
{
for
(
int
j
=-
1
;
j
<=
2
;
j
++
)
{
lf
(
-
1
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
-
1
,
ny
+
i
,
nz
+
j
);
lf
(
0
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
,
ny
+
i
,
nz
+
j
);
lf
(
1
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
1
,
ny
+
i
,
nz
+
j
);
lf
(
2
,
i
,
j
)
=
m3DHisto
->
GetBinContent
(
nx
+
2
,
ny
+
i
,
nz
+
j
);
}
}
splineOrder
=
3
;
result
=
grid_spline
(
splineOrder
,
lf
,
(
x
-
xc
)
/
xscale
,(
y
-
yc
)
/
yscale
,(
z
-
zc
)
/
zscale
);
}
else
{
cout
<<
"Table::InterpolateGridSpline(): Error, illegal number of grid points."
<<
endl
;
result
=
0
;
}
return
result
;
}
void
Table
::
setAutobackup
(
const
char
*
prefix
,
int
freq
)
{
mBackupPrefix
=
string
(
prefix
);
mBackupFrequence
=
freq
;
}
void
Table
::
backup
(
int
backupBin
)
{
ostringstream
backupFilenameStream
;
backupFilenameStream
<<
mBackupPrefix
.
c_str
()
<<
"_backup."
<<
static_cast
<
int
>
(
getpid
())
<<
'.'
<<
mID
<<
'.'
<<
backupBin
<<
".root"
;
string
filename
=
backupFilenameStream
.
str
();
TFile
file
(
filename
.
c_str
(),
"RECREATE"
);
m3DHisto
->
Write
();
file
.
Close
();
if
(
filename
!=
mLastBackupFilename
)
remove
(
mLastBackupFilename
.
c_str
());
mLastBackupFilename
=
filename
;
time_t
now
=
time
(
0
);
cout
<<
"Table::backup(): autobackup performed, file = '"
<<
mLastBackupFilename
.
c_str
()
<<
"', time = "
<<
ctime
(
&
now
);
}
int
Table
::
globalBin
(
int
binx
,
int
biny
,
int
binz
)
const
{
int
nbinx
=
m3DHisto
->
GetXaxis
()
->
GetNbins
();
int
nbiny
=
m3DHisto
->
GetYaxis
()
->
GetNbins
();
return
(
binz
-
1
)
*
nbinx
*
nbiny
+
(
biny
-
1
)
*
nbinx
+
(
binx
-
1
);
}
File Metadata
Details
Attached
Mime Type
text/x-c
Expires
Sat, Dec 21, 1:25 PM (19 h, 35 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4015342
Default Alt Text
Table.cpp (26 KB)
Attached To
rSARTRESVN sartresvn
Event Timeline
Log In to Comment