Fix in the last caall to CleanOwnPrimaryVertex
[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"
25#include <fstream>
26#include <iostream>
27using namespace std;
28
29ClassImp(AliHEPDataParser)
30
31AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName("")
32{
33 // default ctor
34
35}
36
37AliHEPDataParser::AliHEPDataParser(TH1 * hStat, TH1 * hSyst): TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName("y")
38{
39 //ctor
40 fHistStat = hStat;
41 fHistSyst = hSyst;
42 fHEPDataFileLines = new TObjArray;
43
44}
45AliHEPDataParser::AliHEPDataParser(TGraph * grStat, TGraph * grSyst): TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName("")
46{
47 // ctor
48 fGraphStat = grStat;
49 fGraphSyst = grSyst;
50 fHEPDataFileLines = new TObjArray;
51}
52
53AliHEPDataParser::AliHEPDataParser(const char * hepfileName): TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName("y")
54{
55 //ctor
56 // Put results in graphs
57 fGraphSyst = new TGraphAsymmErrors();
58 fGraphStat = new TGraphAsymmErrors();
59 ifstream infile;
60 infile.open(hepfileName);
61 TString line;
62 Int_t ipoints = 0;
63 while (line.ReadLine(infile)) {
64 TObjArray * tokens = line.Tokenize(" ");
65 if(tokens->GetEntries() < 1) {
66 delete tokens;
67 AliError("not enough columns");
68 return;
69 }
70 // TODO: Assumes the format
71 // binmin binmax y +-stat +-syst. Try to make it smarter...
72 TObjString * binMin = (TObjString*) tokens->At(0);
73 TObjString * binMax = (TObjString*) tokens->At(1);
74 TObjString * value = (TObjString*) tokens->At(2);
75 TObjString * stat = (TObjString*) tokens->At(3);
76 TObjString * syst = (TObjString*) tokens->At(4);
77 stat->String().ReplaceAll("+-","");
78 if(syst) syst->String().ReplaceAll("+-","");
79 if (!binMin->String().Atof()) continue; // skip headers
80 Float_t binCenter = (binMax->String().Atof() + binMin->String().Atof())/2;
81 Float_t binWidth = (binMax->String().Atof() - binMin->String().Atof())/2;
82 cout << line.Data() << endl;//<< " " << binMin->String().Atof() <<" " << binCenter << " " << binWidth << endl;
83
84
85 fGraphStat->SetPoint(ipoints, binCenter, value->String().Atof());
86 fGraphSyst->SetPoint(ipoints, binCenter, value->String().Atof());
87 ((TGraphAsymmErrors*)fGraphStat)->SetPointError(ipoints,
88 binWidth,
89 binWidth,
90 stat->String().Atof(),
91 stat->String().Atof());
92 if(syst) ((TGraphAsymmErrors*)fGraphSyst)->SetPointError(ipoints,
93 binWidth,
94 binWidth,
95 syst->String().Atof(),
96 syst->String().Atof());
97 ipoints++;
98 }
99 infile.close();
100
101
102}
103
104AliHEPDataParser::~AliHEPDataParser(){
105 // dtor
106 if(fHistStat) delete fHistStat;
107 if(fHistSyst) delete fHistSyst;
108 if(fGraphStat) delete fGraphStat;
109 if(fGraphSyst) delete fGraphSyst;
110 if(fHEPDataFileLines) delete fHEPDataFileLines;
111}
112
113void AliHEPDataParser::SaveHEPDataFile(const char * hepfileName, Bool_t trueUseGraphFalesUseHisto) {
114 // Fills fHEPDataFileLines and saves its content to a file
115 if(!fHEPDataFileLines) fHEPDataFileLines = new TObjArray;
116 if(trueUseGraphFalesUseHisto) {
117 if(!fGraphStat) {
118 AliError("Graph not set");
119 return;
120 }
121 Bool_t asym = kFALSE; // check if this has asymmetric errors
122 if (!strcmp(fGraphStat->ClassName(), "TGraphErrors")) asym = kFALSE;
123 else if (!strcmp(fGraphStat->ClassName(), "TGraphAsymmErrors")) asym = kTRUE;
124 else {AliError("Unsupported graph type"); return;}
125 Int_t npoint = fGraphStat->GetN();
126 if(asym) AliInfo("Assymmetric errors");
127 for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
128 if(ipoint == 0) {
129 if(fGraphSyst) {
130 if (asym)
131 fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +stat -stat +syst -syst", fValueName.Data())));
132 else
133 fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +-stat +-syst", fValueName.Data())));
134 }
135 else {
136 if(asym)
137 fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +stat -stat", fValueName.Data())));
138 else
139 fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +-stat", fValueName.Data())));
140 }
141 }
142 // Skip empty bins
143 if(!fGraphStat->GetY()[ipoint]) continue;
144 TObjString * line = new TObjString;
145 if(fGraphSyst) {
146 if (asym)
147 line->String().Form("%f %f +%f -%f +%f -%f",
148 fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
149 ((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint],
150 ((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint],
151 ((TGraphAsymmErrors*)fGraphSyst)->GetEYhigh()[ipoint],
152 ((TGraphAsymmErrors*)fGraphSyst)->GetEYlow()[ipoint]);
153 else
154 line->String().Form("%f %f +-%f +-%f",
155 fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
156 ((TGraphErrors*)fGraphStat)->GetEY()[ipoint],
157 ((TGraphErrors*)fGraphSyst)->GetEY()[ipoint]);
158
159 fHEPDataFileLines->Add(line);
160 } else {
161 if (asym)
162 line->String().Form("%f %f +%f -%f",
163 fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
164 ((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], ((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint]);
165 else {
166 line->String().Form("%f %f +-%f",
167 fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
168 ((TGraphErrors*)fGraphStat)->GetEY()[ipoint]);
169 }
170
171 fHEPDataFileLines->Add(line);
172 }
173 }
174 }
175 else {
176 if(!fHistStat) {
177 AliError("Hist not set");
178 return;
179 }
180 Int_t nbin = fHistStat->GetNbinsX();
181
182 for(Int_t ibin = 1; ibin <= nbin; ibin++){
183 if(ibin == 1) {
184 if(fHistSyst)
185 fHEPDataFileLines->Add(new TObjString(Form("BinLow BinHigh %s +-stat +-syst", fValueName.Data())));
186 else
187 fHEPDataFileLines->Add(new TObjString(Form("BinLow BinHigh %s +-stat", fValueName.Data())));
188 }
189 // Skip empty bins
190 if(!fHistStat->GetBinContent(ibin)) continue;
191 TObjString * line = new TObjString;
192 if(fHistSyst) {
193 line->String().Form("%f %f %f +-%f +-%f",
194 fHistStat->GetBinLowEdge(ibin),
195 fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin),
196 fHistStat->GetBinContent(ibin), fHistStat->GetBinError(ibin), fHistSyst->GetBinError(ibin));
197 fHEPDataFileLines->Add(line);
198 } else {
199 line->String().Form("%f %f %f +-%f",
200 fHistStat->GetBinLowEdge(ibin),
201 fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin),
202 fHistStat->GetBinContent(ibin), fHistStat->GetBinError(ibin));
203 fHEPDataFileLines->Add(line);
204 }
205 // delete line;
206 }
207 }
208
209 TIterator * lineIter = fHEPDataFileLines->MakeIterator();
210 TObjString * obj = 0;
211 ofstream outfile;
212 outfile.open (hepfileName);
213 cout << "Saving HEP File " << hepfileName << endl;
214
215 while ((obj = (TObjString*) lineIter->Next())) {
216 cout << obj->String().Data() << endl;
217 outfile << obj->String().Data() << endl;
218 }
219 outfile.close();
220}
221
222