Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / ANALYSIS / AliHEPDataParser.cxx
CommitLineData
a5ae512a 1//-------------------------------------------------------------------------
2// Implementation of Class AliHEPDataParser
3//
4// This class is used to save the content of hisograms/graphs in the
5// HEP data format and viceversa. The HEP data format is not strictly
6// defined and there are many variants, the class only support a few
7// of them. More will be added as needed. The input can be a set of
8// 2 TH1, TGraphAsymmErrors or TGraphErrors (one for the stat and one
9// for the syst error). If the second one is a null pointer, only the
10// stat error is printed. The class can also import hepdata ascii
11// file (very preliminary)
12//
13// Author: Michele Floris, CERN
14//-------------------------------------------------------------------------
15
16
17#include "AliHEPDataParser.h"
18#include "AliLog.h"
19#include "TGraphAsymmErrors.h"
20#include "TGraph.h"
21#include "TGraphErrors.h"
22#include "TH1.h"
23#include "TObjArray.h"
24#include "TObjString.h"
d5e9ea5e 25#include "TMath.h"
a5ae512a 26#include <fstream>
27#include <iostream>
28using namespace std;
29
30ClassImp(AliHEPDataParser)
31
d5e9ea5e 32AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName(""), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
a5ae512a 33{
34 // default ctor
35
36}
37
d5e9ea5e 38AliHEPDataParser::AliHEPDataParser(TH1 * hStat, TH1 * hSyst): TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName("y"), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
a5ae512a 39{
40 //ctor
41 fHistStat = hStat;
42 fHistSyst = hSyst;
43 fHEPDataFileLines = new TObjArray;
44
45}
d5e9ea5e 46AliHEPDataParser::AliHEPDataParser(TGraph * grStat, TGraph * grSyst): TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName(""), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
a5ae512a 47{
48 // ctor
49 fGraphStat = grStat;
50 fGraphSyst = grSyst;
51 fHEPDataFileLines = new TObjArray;
52}
53
4bc92dcb 54AliHEPDataParser::AliHEPDataParser(const char * hepfileName): TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName("y"), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
a5ae512a 55{
56 //ctor
57 // Put results in graphs
58 fGraphSyst = new TGraphAsymmErrors();
59 fGraphStat = new TGraphAsymmErrors();
60 ifstream infile;
61 infile.open(hepfileName);
62 TString line;
63 Int_t ipoints = 0;
64 while (line.ReadLine(infile)) {
d5e9ea5e 65 TObjArray * tokens = line.Tokenize(" \t");
4bc92dcb 66 if(!tokens) continue;
d5e9ea5e 67 if(! ((TObjString*) tokens->At(0))->String().Atof()) {
68 //The first column is not a number: those are the headers: skip!
a5ae512a 69 delete tokens;
d5e9ea5e 70 continue;
71 }
72 if(tokens->GetEntries() != 8) {
73 // this should never happen!
74 delete tokens;
75 AliError(Form("Wrong number of columns %d! Assuming [x xlow xhigh y dystat+ dystat- dysyst+ dysyst-]", tokens->GetEntries()));
76 cout << line.Data() << endl;
a5ae512a 77 return;
d5e9ea5e 78 //continue;
a5ae512a 79 }
d5e9ea5e 80 // FIXME: Assumes the format
81 // x xlow xhigh y dystat+ dystat- dysyst+ dysyst-
82 TObjString * xbin = (TObjString*) tokens->At(0);
83 TObjString * xlow = (TObjString*) tokens->At(1);
84 TObjString * xhigh = (TObjString*) tokens->At(2);
85 TObjString * value = (TObjString*) tokens->At(3);
86 TObjString * statp = (TObjString*) tokens->At(4);
87 TObjString * statm = (TObjString*) tokens->At(5);
88 TObjString * systp = (TObjString*) tokens->At(6);
89 TObjString * systm = (TObjString*) tokens->At(7);
90 statm->String().ReplaceAll("+-","");
91 statp->String().ReplaceAll("+-","");
92 if(systp) systp->String().ReplaceAll("+-","");
93 if(systm) systm->String().ReplaceAll("+-","");
94 // if (!binMin->String().Atof()) {delete tokens; continue;} // skip headers
95 Float_t binCenter = xbin->String().Atof();
96 Float_t binWidth = (xlow->String().Atof() - xhigh->String().Atof())/2;
a5ae512a 97
98
99 fGraphStat->SetPoint(ipoints, binCenter, value->String().Atof());
100 fGraphSyst->SetPoint(ipoints, binCenter, value->String().Atof());
101 ((TGraphAsymmErrors*)fGraphStat)->SetPointError(ipoints,
102 binWidth,
103 binWidth,
d5e9ea5e 104 statp->String().Atof(),
105 statm->String().Atof());
106 if(systp && systm) ((TGraphAsymmErrors*)fGraphSyst)->SetPointError(ipoints,
107 binWidth,
108 binWidth,
109 systp->String().Atof(),
110 systm->String().Atof());
a5ae512a 111 ipoints++;
09d5920f 112 delete tokens;
a5ae512a 113 }
114 infile.close();
115
116
117}
118
119AliHEPDataParser::~AliHEPDataParser(){
120 // dtor
121 if(fHistStat) delete fHistStat;
122 if(fHistSyst) delete fHistSyst;
123 if(fGraphStat) delete fGraphStat;
124 if(fGraphSyst) delete fGraphSyst;
125 if(fHEPDataFileLines) delete fHEPDataFileLines;
126}
127
128void AliHEPDataParser::SaveHEPDataFile(const char * hepfileName, Bool_t trueUseGraphFalesUseHisto) {
d5e9ea5e 129
a5ae512a 130 // Fills fHEPDataFileLines and saves its content to a file
131 if(!fHEPDataFileLines) fHEPDataFileLines = new TObjArray;
d5e9ea5e 132 // Write headers if relevant
133 if(fTitle.Length()) fHEPDataFileLines->Add(new TObjString(fTitle));
134 if(fReaction.Length()) fHEPDataFileLines->Add(new TObjString(fReaction));
135 if(fEnergy.Length()) fHEPDataFileLines->Add(new TObjString(fEnergy));
136 if(fRapidityRange.Length()) fHEPDataFileLines->Add(new TObjString(fRapidityRange));
137 if(!fValueName.Length() || !fXaxisName.Length()) AliFatal("At least x and y titles should be given!");
138 fHEPDataFileLines->Add(new TObjString(Form("x: %s", fXaxisName.Data())));
139 fHEPDataFileLines->Add(new TObjString(Form("y: %s", fValueName.Data())));
140
141
a5ae512a 142 if(trueUseGraphFalesUseHisto) {
d5e9ea5e 143 AliWarning("Graph saving not thoroughly tested!!");
a5ae512a 144 if(!fGraphStat) {
145 AliError("Graph not set");
146 return;
147 }
148 Bool_t asym = kFALSE; // check if this has asymmetric errors
149 if (!strcmp(fGraphStat->ClassName(), "TGraphErrors")) asym = kFALSE;
150 else if (!strcmp(fGraphStat->ClassName(), "TGraphAsymmErrors")) asym = kTRUE;
151 else {AliError("Unsupported graph type"); return;}
152 Int_t npoint = fGraphStat->GetN();
153 if(asym) AliInfo("Assymmetric errors");
154 for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
155 if(ipoint == 0) {
156 if(fGraphSyst) {
d5e9ea5e 157 fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat\t+syst\t-syst"));
a5ae512a 158 }
159 else {
d5e9ea5e 160 fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat"));
a5ae512a 161 }
162 }
163 // Skip empty bins
164 if(!fGraphStat->GetY()[ipoint]) continue;
d5e9ea5e 165 TObjString * line = new TObjString;
a5ae512a 166 if(fGraphSyst) {
167 if (asym)
d5e9ea5e 168 line->String().Form("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s",
169 // line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s",
170 GetFixWidthCol(fGraphStat->GetX()[ipoint], 10).Data(),
171 GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetX()[ipoint]-((TGraphAsymmErrors*)fGraphStat)->GetEXlow()[ipoint], fPrecision), 10).Data(),
172 GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetX()[ipoint]+((TGraphAsymmErrors*)fGraphStat)->GetEXhigh()[ipoint], fPrecision), 10).Data(),
173 GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint], fPrecision), 10).Data(),
174 GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], fPrecision), 10).Data(),
175 GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint] , fPrecision), 10).Data(),
176 GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphSyst)->GetEYhigh()[ipoint], fPrecision), 10).Data(),
177 GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphSyst)->GetEYlow()[ipoint] , fPrecision), 10).Data());
a5ae512a 178 else
d5e9ea5e 179 line->String().Form("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s",
180 // line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s",
181 GetFixWidthCol(fGraphStat->GetX()[ipoint], 10).Data(),
182 GetFixWidthCol(fGraphStat->GetX()[ipoint]-fGraphStat->GetEX()[ipoint], 10).Data(),
183 GetFixWidthCol(fGraphStat->GetX()[ipoint]+fGraphStat->GetEX()[ipoint], 10).Data(),
184 GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint], fPrecision), 10).Data(),
185 GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data(),
186 GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data(),
187 GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphSyst)->GetEY()[ipoint], fPrecision), 10).Data(),
188 GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphSyst)->GetEY()[ipoint], fPrecision), 10).Data());
189 // line->String().Form("%f %f +-%f +-%f",
190 // fGraphStat->GetX()[ipoint], RoundToSignificantFigures(fGraphStat->GetY()[ipoint],fPrecision),
191 // RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint],fPrecision),
192 // RoundToSignificantFigures(((TGraphErrors*)fGraphSyst)->GetEY()[ipoint],fPrecision));
a5ae512a 193
194 fHEPDataFileLines->Add(line);
195 } else {
196 if (asym)
d5e9ea5e 197 line->String().Form("%s\t%s\t%s\t%s\t%s\t%s",
198 // line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s",
199 GetFixWidthCol(fGraphStat->GetX()[ipoint], 10).Data(),
200 GetFixWidthCol(fGraphStat->GetX()[ipoint]-fGraphStat->GetEXlow()[ipoint], 10).Data(),
201 GetFixWidthCol(fGraphStat->GetX()[ipoint]+fGraphStat->GetEXhigh()[ipoint], 10).Data(),
202 GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint], fPrecision), 10).Data(),
203 GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], fPrecision), 10).Data(),
204 GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint] , fPrecision), 10).Data());
205 else
206 line->String().Form("%s\t%s\t%s\t%s\t%s\t%s",
207 // line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s",
208 GetFixWidthCol(fGraphStat->GetX()[ipoint], 10).Data(),
209 GetFixWidthCol(fGraphStat->GetX()[ipoint]-fGraphStat->GetEX()[ipoint], 10).Data(),
210 GetFixWidthCol(fGraphStat->GetX()[ipoint]+fGraphStat->GetEX()[ipoint], 10).Data(),
211 GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint], fPrecision), 10).Data(),
212 GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data(),
213 GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data());
a5ae512a 214
215 fHEPDataFileLines->Add(line);
216 }
217 }
218 }
219 else {
220 if(!fHistStat) {
221 AliError("Hist not set");
222 return;
223 }
224 Int_t nbin = fHistStat->GetNbinsX();
225
226 for(Int_t ibin = 1; ibin <= nbin; ibin++){
227 if(ibin == 1) {
228 if(fHistSyst)
d5e9ea5e 229 fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-\t\tdysyst+\t\tdysyst-"));
a5ae512a 230 else
d5e9ea5e 231 fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-"));
a5ae512a 232 }
233 // Skip empty bins
234 if(!fHistStat->GetBinContent(ibin)) continue;
235 TObjString * line = new TObjString;
236 if(fHistSyst) {
d5e9ea5e 237
238 line->String().Form("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s",
239 // line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s",
240 GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin)/2, 12).Data(),
241 GetFixWidthCol(fHistStat->GetBinLowEdge(ibin), 12).Data(),
242 GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin), 12).Data(),
243 GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinContent(ibin),fPrecision), 12).Data(),
244 GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin), fPrecision), 12).Data(),
245 GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin), fPrecision), 12).Data(),
246 GetFixWidthCol(RoundToSignificantFigures(fHistSyst->GetBinError(ibin), fPrecision), 12).Data(),
247 GetFixWidthCol(RoundToSignificantFigures(fHistSyst->GetBinError(ibin), fPrecision), 12).Data());
248
a5ae512a 249 fHEPDataFileLines->Add(line);
250 } else {
d5e9ea5e 251 line->String().Form("%s\t%s\t%s\t%s\t%s\t%s",
252 GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin)/2, 10).Data(),
253 GetFixWidthCol(fHistStat->GetBinLowEdge(ibin), 10).Data(),
254 GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin), 10).Data(),
255 GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinContent(ibin),fPrecision), 10).Data(),
256 GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin), fPrecision), 10).Data(),
257 GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin), fPrecision), 10).Data());
a5ae512a 258 fHEPDataFileLines->Add(line);
259 }
260 // delete line;
261 }
262 }
263
264 TIterator * lineIter = fHEPDataFileLines->MakeIterator();
265 TObjString * obj = 0;
266 ofstream outfile;
267 outfile.open (hepfileName);
268 cout << "Saving HEP File " << hepfileName << endl;
269
270 while ((obj = (TObjString*) lineIter->Next())) {
271 cout << obj->String().Data() << endl;
272 outfile << obj->String().Data() << endl;
273 }
274 outfile.close();
275}
276
277
d5e9ea5e 278Double_t AliHEPDataParser::RoundToSignificantFigures(double num, int n) {
279 // Rounds num to n significant digits.
280 // Recipe from :http://stackoverflow.com/questions/202302/rounding-to-an-arbitrary-number-of-significant-digits
281 // Basically the log is used to determine the number of leading 0s, than convert to an integer by multipliing by the expo,
282 // round the integer and shift back.
283 if(num == 0) {
284 return 0;
285 }
286
287 Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num));
288 Int_t power = n - (int) d;
289
290 Double_t magnitude = TMath::Power(10, power);
291 Long_t shifted = TMath::Nint(num*magnitude);
292 return shifted/magnitude;
293
294}
295
296TString AliHEPDataParser::GetFixWidthCol(Double_t number, Int_t width) {
297
298 // Formats a column to fixed width
299 TString col;
300 char format[100];
301 snprintf(format,100,"%%%d#g", fPrecision);
302 col.Form(format, number);
303 if(col.Length()>width) AliError("larger than width, cannot align!");
304
305 if(col.Contains("e"))
306 while (col.Length() < width) col.Append(" ");
307 else
308 while (col.Length() < width) col.Append("0");
309
310 return col;
311}