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