Index: trunk/src/tableVarianceMaker.cpp
===================================================================
--- trunk/src/tableVarianceMaker.cpp (revision 422)
+++ trunk/src/tableVarianceMaker.cpp (revision 423)
@@ -1,330 +1,336 @@
//==============================================================================
// tableVarianceMaker.cpp
//
// Copyright (C) 2016-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
// $Date$
// $Author: ullrich $
//==============================================================================
//
// Utility program to create variance_A tables from existing mean_A2 and
// mean_A tables.
//
// Several checks are performed to ensure the integrity of the created
// table but in any case the new tables should be checked with
// tableInspector.
//
// Usage:
// tableVarianceMaker [-n | -o output] file_mean_A2 file_mean_A
// -n Do not create the new tables but perform all checks
// -o file Output file where new table gets stored
//
// If neither -n nor -o output is given the content of the new file is
// written to the screen (stdout).
//
// Note it is a good idea to always check with option -n first before
// starting to create new tables.
//
//==============================================================================
#include "TFile.h"
#include "TH3F.h"
#include "Table.h"
#include <iostream>
#include <vector>
#include <sstream>
#include <set>
#include <unistd.h>
#include <limits>
#include <cmath>
#include <cstdlib>
#include <bitset>
#define PR(x) cout << #x << " = " << (x) << endl;
using namespace std;
void usage(const char* prog)
{
cout << "Usage: " << prog << " [-n | -o output] file_mean_A2 file_mean_A" << endl;
cout << "\t\t\t" << "-n Do not create the new tables but perform all checks" << endl;
cout << "\t\t\t" << "-o file Output file where new table gets stored" << endl;
}
TH3F* getHistoFromFile(string&);
int main(int argc, char **argv)
{
//
// Handle command line arguments
//
if (argc < 3) {
usage(argv[0]);
return 2;
}
bool noOutput = false;
bool printOnly = false;
string outputFile;
int ch;
while ((ch = getopt(argc, argv, "no:")) != -1) {
switch (ch) {
case 'o':
outputFile = string(optarg);
break;
case 'n':
noOutput = true;
break;
case '?':
default:
usage(argv[0]);
return 2;
break;
}
}
if (optind == argc) {
usage(argv[0]);
return 2;
}
if (noOutput && !outputFile.empty()) {
usage(argv[0]);
return 2;
}
printOnly = (!noOutput && outputFile.empty());
vector<string> inputFiles;
for (int index = optind; index < argc; index++)
inputFiles.push_back(string(argv[index]));
if (inputFiles.size() != 2) {
usage(argv[0]);
return 2;
}
//
// First we perform all possible checks to ensure the tables
// are identical other than their content. Using the Table
// class makes things easier.
//
Table t1, t2;
bool ok = t1.read(inputFiles[0].c_str());
if (!ok) {
cout << "Error: failed opening input file '" << inputFiles[0].c_str() << "'." << endl;
return 1;
}
ok = t2.read(inputFiles[1].c_str());
if (!ok) {
cout << "Error: failed opening input file '" << inputFiles[1].c_str() << "'." << endl;
return 1;
}
ok = (t1.isMeanA() && t2.isMeanA2()) || (t1.isMeanA2() && t2.isMeanA()); // xor
if (!ok) {
cout << "Error: one input file must contain an amplitude mode of mean_A data, the other mean_A2 data." << endl;
return 1;
}
if (t1.minQ2() != t2.minQ2()) {
cout << "Error: lower Q2 bound is different for the two input tables." << endl;
return 1;
}
if (t1.maxQ2() != t2.maxQ2()) {
cout << "Error: upper Q2 bound is different for the two input tables." << endl;
return 1;
}
if (t1.minW2() != t2.minW2()) {
cout << "Error: lower W2 bound is different for the two input tables." << endl;
return 1;
}
if (t1.maxW2() != t2.maxW2()) {
cout << "Error: upper W2 bound is different for the two input tables." << endl;
return 1;
}
if (t1.minT() != t2.minT()) {
cout << "Error: lower t bound is different for the two input tables." << endl;
return 1;
}
if (t1.maxT() != t2.maxT()) {
cout << "Error: upper t bound is different for the two input tables." << endl;
return 1;
}
if (t1.binWidthQ2() != t2.binWidthQ2()) {
cout << "Error: Bin width of Q2 axis is different for the two input tables." << endl;
return 1;
}
if (t1.binWidthW2() != t2.binWidthW2()) {
cout << "Error: Bin width of Q2 axis is different for the two input tables." << endl;
return 1;
}
if (t1.binWidthT() != t2.binWidthT()) {
cout << "Error: Bin width of Q2 axis is different for the two input tables." << endl;
return 1;
}
if (t1.polarization() != t2.polarization()) {
cout << "Error: The two input files hold amplitudes for different polarizations." << endl;
return 1;
}
if (t1.A() != t2.A()) {
cout << "Error: The two input files are for different nuclei A." << endl;
return 1;
}
if (t1.vectorMesonId() != t2.vectorMesonId()) {
cout << "Error: The two input files are for different vector mesons." << endl;
return 1;
}
if (t1.dipoleModelType() != t2.dipoleModelType()) {
cout << "Error: The two input files are for different dipole models." << endl;
return 1;
}
if (t1.dipoleModelParameterSet() != t2.dipoleModelParameterSet()) {
cout << "Error: The two input files are were generated with different dipole model parameter sets." << endl;
return 1;
}
if (t1.priority() != t2.priority()) {
cout << "Error: The two input files are were generated with different priorities." << endl;
return 1;
}
//
// If we get here the two histos should just work fine. Note that we
// did not check if the axis are lin/log since if the upper/lower bound
// and the bin width are the same we are good.
// The *only* thing that matters from here on is if the content is in
// lin or log.
//
cout << "OK - The two input files are consistent and can be used to create a variance_A table." << endl;
if (noOutput) return 0;
//
// Order input files: first = mean_A2, second = mean_A
//
pair<string, string> rootfiles(inputFiles[0], inputFiles[1]);
- if (rootfiles.first != t1.isMeanA2()) swap(rootfiles.first, rootfiles.second);
+ if (!t1.isMeanA2()) swap(rootfiles.first, rootfiles.second);
//
// Open the files and get the two tables to deal with.
//
//
pair<TH3F*, TH3F*> histograms;
histograms.first = getHistoFromFile(rootfiles.first);
histograms.second = getHistoFromFile(rootfiles.second);
//
// Create the variance histogram and set the ID.
// Also find our if the content is in log or lin.
//
TH3F* vhisto = reinterpret_cast<TH3F*>(histograms.first->Clone());
vhisto->SetDirectory(0);
uint64_t theID = atoll(vhisto->GetTitle());
// set bit 0 and 33 to 0, set bit 45
bitset<64> bits(theID);
bits[0] = 0;
bits[33] = 0;
bits[45] = 1;
theID = bits.to_ullong();
ostringstream titlestream;
titlestream << theID;
string title = titlestream.str();
vhisto->Reset("M");
vhisto->SetTitle(title.c_str());
bool varianceAndA2HaveLogContent = bits.test(32);
uint64_t meanA_ID = atoll(histograms.second->GetTitle());
bitset<64> meanA_Bits(meanA_ID);
bool meanAHasLogContent = meanA_Bits.test(32);
cout << "Created new table in memory." << endl;
//
// Fill the variance histo (all calculations happen here)
//
-
+ PR(t1.isMeanA());
+ PR(t1.isMeanA2());
+ PR(t2.isMeanA());
+ PR(t2.isMeanA2());
int nx = vhisto->GetNbinsX();
int ny = vhisto->GetNbinsY();
int nz = vhisto->GetNbinsZ();
for (int ix = 1; ix <= nx; ix++) {
for (int iy = 1; iy <= ny; iy++) {
for (int iz = 1; iz <= nz; iz++) {
// old table
double A2_Content = histograms.first->GetBinContent(ix, iy, iz);
double A_Content = histograms.second->GetBinContent(ix, iy, iz);
+
if (varianceAndA2HaveLogContent) A2_Content = exp(A2_Content);
if (meanAHasLogContent) A_Content = exp(A_Content);
double variance = A2_Content - A_Content*A_Content;
+ if (variance == 0 && varianceAndA2HaveLogContent) variance += numeric_limits<float>::min();
if (printOnly) {
cout << "bin = " << ix << ", " << iy << ", " << iz
<< "; mean_A2 = " << A2_Content
<< ", mean_A = " << A_Content
+ << ", mean_A^2 = " << A_Content*A_Content
<< ", variance_A = " << variance
<< " (" << (varianceAndA2HaveLogContent? log(variance) : variance) << ")" << endl;
}
if (varianceAndA2HaveLogContent) variance = log(variance);
vhisto->SetBinContent(ix, iy, iz, variance);
vhisto->SetBinError(ix, iy, iz, 0);
}
}
}
cout << "Variance table is filled. Content is stored in " << (varianceAndA2HaveLogContent ? "log" : "lin") << " scale." << endl;
//
// Store the table
//
if (!outputFile.empty()) {
cout << "Attempting to write table into new file '" << outputFile << "'." << endl;
TFile *hfile = new TFile(outputFile.c_str(),"NEW");
if (hfile && hfile->IsZombie()) return 1;
if (!hfile) {
cout << "Error writing variance table. Cannot open new file '" << outputFile.c_str() << "'." << endl;
}
else {
vhisto->Write();
hfile->Close();
cout << "Variance table written to file '" << outputFile.c_str() << "'." << endl;
delete vhisto;
delete histograms.first;
delete histograms.second;
}
}
cout << "All done. Bye." << endl;
return 0;
}
TH3F* getHistoFromFile(string& filename)
{
TFile *file = TFile::Open(filename.c_str(),"READ");
if (file == 0) {
cout << "Error: failed opening file '" << filename.c_str() << "'." << endl;
exit(1);
}
TH3F *histo = reinterpret_cast<TH3F*>(file->Get("table"));
if (histo == 0) {
cout << "Error: failed retrieving table from file '" << filename.c_str() << "'." << endl;
exit(1);
}
histo->SetDirectory(0);
file->Close();
return histo;
}