[u/mrichter/AliRoot.git] / ANALYSIS / AliHEPDataParser.cxx

//-------------------------------------------------------------------------
//           Implementation of Class AliHEPDataParser
//
//  This class is used to save the content of hisograms/graphs in the
//  HEP data format and viceversa. The HEP data format is not strictly
//  defined and there are many variants, the class only support a few
//  of them. More will be added as needed.  The input can be a set of
//  2 TH1, TGraphAsymmErrors or TGraphErrors (one for the stat and one
//  for the syst error). If the second one is a null pointer, only the
//  stat error is printed. The class can also import hepdata ascii
//  file (very preliminary)
// 
//  Author: Michele Floris, CERN
//-------------------------------------------------------------------------


#include "AliHEPDataParser.h"
#include "AliLog.h"
#include "TGraphAsymmErrors.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TH1.h"
#include "TObjArray.h"
#include "TObjString.h"
#include <fstream>
#include <iostream>
using namespace std;

ClassImp(AliHEPDataParser)

AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName("")
{
  // default ctor

}

AliHEPDataParser::AliHEPDataParser(TH1 * hStat, TH1 * hSyst): TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName("y")
{
  //ctor
  fHistStat = hStat;
  fHistSyst = hSyst;
  fHEPDataFileLines = new TObjArray;

}
AliHEPDataParser::AliHEPDataParser(TGraph * grStat, TGraph * grSyst): TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName("")
{
  // ctor
  fGraphStat = grStat;
  fGraphSyst = grSyst;
  fHEPDataFileLines = new TObjArray;
}

AliHEPDataParser::AliHEPDataParser(const char * hepfileName): TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName("y")
{
  //ctor
  // Put results in graphs
  fGraphSyst = new TGraphAsymmErrors();
  fGraphStat = new TGraphAsymmErrors();
  ifstream infile;
  infile.open(hepfileName);
  TString line;
  Int_t ipoints = 0;
  while (line.ReadLine(infile)) {
    TObjArray * tokens = line.Tokenize(" ");
    if(tokens->GetEntries() < 1) {
      delete tokens;
      AliError("not enough columns");
      return;      
    }
    // TODO: Assumes the format
    // binmin binmax y +-stat +-syst. Try to make it smarter...
    TObjString * binMin = (TObjString*) tokens->At(0);
    TObjString * binMax = (TObjString*) tokens->At(1);
    TObjString * value  = (TObjString*) tokens->At(2);
    TObjString * stat   = (TObjString*) tokens->At(3);
    TObjString * syst   = (TObjString*) tokens->At(4);
    stat->String().ReplaceAll("+-","");
    if(syst) syst->String().ReplaceAll("+-","");
    if (!binMin->String().Atof()) continue; // skip headers
    Float_t binCenter = (binMax->String().Atof() + binMin->String().Atof())/2;
    Float_t binWidth  =  (binMax->String().Atof() - binMin->String().Atof())/2;
    cout << line.Data() << endl;//<< " " << binMin->String().Atof() <<" " <<  binCenter << " " << binWidth << endl;


    fGraphStat->SetPoint(ipoints, binCenter, value->String().Atof());
    fGraphSyst->SetPoint(ipoints, binCenter, value->String().Atof());
    ((TGraphAsymmErrors*)fGraphStat)->SetPointError(ipoints, 
						    binWidth,
						    binWidth,
						    stat->String().Atof(),
						    stat->String().Atof());
    if(syst) ((TGraphAsymmErrors*)fGraphSyst)->SetPointError(ipoints, 
							     binWidth,
							     binWidth,
							     syst->String().Atof(),
							     syst->String().Atof());
    ipoints++;
  }
  infile.close();
    

}

AliHEPDataParser::~AliHEPDataParser(){
  // dtor
  if(fHistStat) delete fHistStat;
  if(fHistSyst) delete fHistSyst;
  if(fGraphStat) delete fGraphStat;
  if(fGraphSyst) delete fGraphSyst;
  if(fHEPDataFileLines) delete fHEPDataFileLines;
}
  
void AliHEPDataParser::SaveHEPDataFile(const char * hepfileName, Bool_t trueUseGraphFalesUseHisto) {
  // Fills fHEPDataFileLines and saves its content to a file
  if(!fHEPDataFileLines)   fHEPDataFileLines = new TObjArray;
  if(trueUseGraphFalesUseHisto) {
    if(!fGraphStat) {
      AliError("Graph not set");
      return;
    }
    Bool_t asym = kFALSE; // check if this has asymmetric errors
    if       (!strcmp(fGraphStat->ClassName(), "TGraphErrors")) asym = kFALSE;
    else if  (!strcmp(fGraphStat->ClassName(), "TGraphAsymmErrors")) asym = kTRUE;
    else     {AliError("Unsupported graph type"); return;}
    Int_t npoint = fGraphStat->GetN();
    if(asym) AliInfo("Assymmetric errors");
    for(Int_t ipoint = 0; ipoint < npoint; ipoint++){            
      if(ipoint == 0) {
	if(fGraphSyst) {
	  if (asym)
	    fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +stat -stat +syst -syst", fValueName.Data()))); 
	  else 
	    fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +-stat +-syst", fValueName.Data()))); 
	}
	else {
	  if(asym)
	    fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +stat -stat", fValueName.Data()))); 
	  else 
	    fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +-stat", fValueName.Data()))); 
	}
      }
      // Skip empty bins
      if(!fGraphStat->GetY()[ipoint]) continue;
      TObjString * line = new TObjString;      
      if(fGraphSyst) {
	if (asym)
	  line->String().Form("%f %f +%f -%f +%f -%f", 
			      fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
			      ((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], 
			      ((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint], 
			      ((TGraphAsymmErrors*)fGraphSyst)->GetEYhigh()[ipoint], 
			      ((TGraphAsymmErrors*)fGraphSyst)->GetEYlow()[ipoint]);
	else 
	  line->String().Form("%f %f +-%f +-%f", 
			      fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
			      ((TGraphErrors*)fGraphStat)->GetEY()[ipoint], 
			      ((TGraphErrors*)fGraphSyst)->GetEY()[ipoint]);

	fHEPDataFileLines->Add(line);
      } else {
	if (asym)
	  line->String().Form("%f %f +%f -%f", 
			      fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
			      ((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], ((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint]);
	else { 
	  line->String().Form("%f %f +-%f", 
			      fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
			      ((TGraphErrors*)fGraphStat)->GetEY()[ipoint]);
	}

	fHEPDataFileLines->Add(line);
      }
    }    
  }
  else {
    if(!fHistStat) {
      AliError("Hist not set");
      return;
    }
    Int_t nbin = fHistStat->GetNbinsX();
    
    for(Int_t ibin = 1; ibin <= nbin; ibin++){
      if(ibin == 1) {
	if(fHistSyst)
	  fHEPDataFileLines->Add(new TObjString(Form("BinLow BinHigh %s +-stat +-syst", fValueName.Data()))); 
	else 
	  fHEPDataFileLines->Add(new TObjString(Form("BinLow BinHigh %s +-stat", fValueName.Data())));       	
      }
      // Skip empty bins
      if(!fHistStat->GetBinContent(ibin)) continue;
      TObjString * line = new TObjString;      
      if(fHistSyst) {
	line->String().Form("%f %f %f +-%f +-%f", 
			    fHistStat->GetBinLowEdge(ibin), 
			    fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin), 
			    fHistStat->GetBinContent(ibin), fHistStat->GetBinError(ibin), fHistSyst->GetBinError(ibin));
	fHEPDataFileLines->Add(line);
      } else {
	line->String().Form("%f %f %f +-%f", 
			    fHistStat->GetBinLowEdge(ibin), 
			    fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin), 
			    fHistStat->GetBinContent(ibin), fHistStat->GetBinError(ibin));
	fHEPDataFileLines->Add(line);
      }
      //      delete line;      
    }           
  }

  TIterator * lineIter = fHEPDataFileLines->MakeIterator();
  TObjString * obj = 0;
  ofstream outfile;
  outfile.open (hepfileName);
  cout << "Saving HEP File " << hepfileName << endl;
  
  while ((obj = (TObjString*) lineIter->Next())) {    
    cout << obj->String().Data() << endl;    
    outfile << obj->String().Data() << endl;
  }
  outfile.close();
}
Commit	Line	Data
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>
	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