Resolving the symbols in each library
[u/mrichter/AliRoot.git] / ANALYSIS / ANALYSISalice / AliHEPDataParser.cxx
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 ClassImp(AliHEPDataParser)
31
32 AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName(""), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
33 {
34   // default ctor
35
36 }
37
38 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 {
40   //ctor
41   fHistStat = hStat;
42   fHistSyst = hSyst;
43   fHEPDataFileLines = new TObjArray;
44
45 }
46 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 {
48   // ctor
49   fGraphStat = grStat;
50   fGraphSyst = grSyst;
51   fHEPDataFileLines = new TObjArray;
52 }
53
54 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 {
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)) {
65     TObjArray * tokens = line.Tokenize(" \t");
66     if(!tokens) continue;
67     if(! ((TObjString*) tokens->At(0))->String().Atof()) {
68       //The first column is not a number: those are the headers: skip!
69       delete tokens;
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;
77       return;      
78       //continue;
79     }
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;
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,
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());
111     ipoints++;
112     delete tokens;
113   }
114   infile.close();
115     
116
117 }
118
119 AliHEPDataParser::~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   
128 void AliHEPDataParser::SaveHEPDataFile(const char * hepfileName, Bool_t trueUseGraphFalesUseHisto) {
129
130   // Fills fHEPDataFileLines and saves its content to a file
131   if(!fHEPDataFileLines)   fHEPDataFileLines = new TObjArray;
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
142   if(trueUseGraphFalesUseHisto) {
143     AliWarning("Graph saving not thoroughly tested!!");
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) {
157           fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat\t+syst\t-syst"));
158         }
159         else {
160           fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat"));
161         }
162       }
163       // Skip empty bins
164       if(!fGraphStat->GetY()[ipoint]) continue;
165       TObjString * line = new TObjString;    
166       if(fGraphSyst) {
167         if (asym)
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());
178         else 
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));
193
194         fHEPDataFileLines->Add(line);
195       } else {
196         if (asym)
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());
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)
229           fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-\t\tdysyst+\t\tdysyst-")); 
230         else 
231           fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-"));              
232       }
233       // Skip empty bins
234       if(!fHistStat->GetBinContent(ibin)) continue;
235       TObjString * line = new TObjString;      
236       if(fHistSyst) {
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
249         fHEPDataFileLines->Add(line);
250       } else {
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()); 
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
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   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
296 TString 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 }