coverity
[u/mrichter/AliRoot.git] / ANALYSIS / 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 <fstream>
26 #include <iostream>
27 using namespace std;
28
29 ClassImp(AliHEPDataParser)
30
31 AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName("")
32 {
33   // default ctor
34
35 }
36
37 AliHEPDataParser::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 }
45 AliHEPDataParser::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
53 AliHEPDataParser::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
104 AliHEPDataParser::~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   
113 void 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