[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 "TMath.h"
#include <fstream>
#include <iostream>
using namespace std;

ClassImp(AliHEPDataParser)

AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0),  fHistSyst(0),  fGraphStat(0),  fGraphSyst(0),  fHEPDataFileLines(0), fValueName(""), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
{
  // default ctor

}

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)
{
  //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(""), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
{
  // 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"), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
{
  //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(" \t");
    if(!tokens) continue;
    if(! ((TObjString*) tokens->At(0))->String().Atof()) {
      //The first column is not a number: those are the headers: skip!
      delete tokens;
      continue;
    }
    if(tokens->GetEntries() != 8) {
      // this should never happen!
      delete tokens;
      AliError(Form("Wrong number of columns %d! Assuming [x xlow xhigh y dystat+ dystat- dysyst+ dysyst-]", tokens->GetEntries()));
      cout << line.Data() << endl;
      return;      
      //continue;
    }
    // FIXME: Assumes the format
    // x xlow xhigh y dystat+ dystat- dysyst+ dysyst-
    TObjString * xbin   = (TObjString*) tokens->At(0);
    TObjString * xlow   = (TObjString*) tokens->At(1);
    TObjString * xhigh  = (TObjString*) tokens->At(2);
    TObjString * value  = (TObjString*) tokens->At(3);
    TObjString * statp  = (TObjString*) tokens->At(4);
    TObjString * statm  = (TObjString*) tokens->At(5);
    TObjString * systp  = (TObjString*) tokens->At(6);
    TObjString * systm  = (TObjString*) tokens->At(7);
    statm->String().ReplaceAll("+-","");
    statp->String().ReplaceAll("+-","");
    if(systp) systp->String().ReplaceAll("+-","");
    if(systm) systm->String().ReplaceAll("+-","");
    //    if (!binMin->String().Atof()) {delete tokens; continue;} // skip headers
    Float_t binCenter = xbin->String().Atof();
    Float_t binWidth  =  (xlow->String().Atof() - xhigh->String().Atof())/2;


    fGraphStat->SetPoint(ipoints, binCenter, value->String().Atof());
    fGraphSyst->SetPoint(ipoints, binCenter, value->String().Atof());
    ((TGraphAsymmErrors*)fGraphStat)->SetPointError(ipoints, 
						    binWidth,
						    binWidth,
						    statp->String().Atof(),
						    statm->String().Atof());
    if(systp && systm) ((TGraphAsymmErrors*)fGraphSyst)->SetPointError(ipoints, 
							      binWidth,
							      binWidth,
							      systp->String().Atof(),
							      systm->String().Atof());
    ipoints++;
    delete tokens;
  }
  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;
  // Write headers if relevant
  if(fTitle.Length())         fHEPDataFileLines->Add(new TObjString(fTitle));
  if(fReaction.Length())      fHEPDataFileLines->Add(new TObjString(fReaction));
  if(fEnergy.Length())        fHEPDataFileLines->Add(new TObjString(fEnergy));
  if(fRapidityRange.Length()) fHEPDataFileLines->Add(new TObjString(fRapidityRange));
  if(!fValueName.Length() || !fXaxisName.Length()) AliFatal("At least x and y titles should be given!");
  fHEPDataFileLines->Add(new TObjString(Form("x: %s", fXaxisName.Data())));
  fHEPDataFileLines->Add(new TObjString(Form("y: %s", fValueName.Data())));


  if(trueUseGraphFalesUseHisto) {
    AliWarning("Graph saving not thoroughly tested!!");
    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) {
	  fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat\t+syst\t-syst"));
	}
	else {
	  fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat"));
	}
      }
      // Skip empty bins
      if(!fGraphStat->GetY()[ipoint]) continue;
      TObjString * line = new TObjString;    
      if(fGraphSyst) {
	if (asym)
	  line->String().Form("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s",
			      //	line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s", 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint],                                                                   10).Data(), 			    
			      GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetX()[ipoint]-((TGraphAsymmErrors*)fGraphStat)->GetEXlow()[ipoint],  fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetX()[ipoint]+((TGraphAsymmErrors*)fGraphStat)->GetEXhigh()[ipoint], fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint],                            fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint] , fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphSyst)->GetEYhigh()[ipoint], fPrecision), 10).Data(),
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphSyst)->GetEYlow()[ipoint] , fPrecision), 10).Data());
	else 
	  line->String().Form("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s",
			      //	line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s", 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint],                                                                   10).Data(), 			    
			      GetFixWidthCol(fGraphStat->GetX()[ipoint]-fGraphStat->GetEX()[ipoint],                                    10).Data(), 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint]+fGraphStat->GetEX()[ipoint],                                10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint],                            fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data(),
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphSyst)->GetEY()[ipoint], fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphSyst)->GetEY()[ipoint], fPrecision), 10).Data());
	  // line->String().Form("%f %f +-%f +-%f", 
	  // 		      fGraphStat->GetX()[ipoint], RoundToSignificantFigures(fGraphStat->GetY()[ipoint],fPrecision),
	  // 		      RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint],fPrecision), 
	  // 		      RoundToSignificantFigures(((TGraphErrors*)fGraphSyst)->GetEY()[ipoint],fPrecision));

	fHEPDataFileLines->Add(line);
      } else {
	if (asym)
	  line->String().Form("%s\t%s\t%s\t%s\t%s\t%s",
			      //	line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s", 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint],                                                                   10).Data(), 			    
			      GetFixWidthCol(fGraphStat->GetX()[ipoint]-fGraphStat->GetEXlow()[ipoint],                                    10).Data(), 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint]+fGraphStat->GetEXhigh()[ipoint],                                   10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint],                            fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint] , fPrecision), 10).Data()); 
	else 
	  line->String().Form("%s\t%s\t%s\t%s\t%s\t%s",
			      //	line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s", 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint],                                                                   10).Data(), 			    
			      GetFixWidthCol(fGraphStat->GetX()[ipoint]-fGraphStat->GetEX()[ipoint],                                    10).Data(), 
			      GetFixWidthCol(fGraphStat->GetX()[ipoint]+fGraphStat->GetEX()[ipoint],                                10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(fGraphStat->GetY()[ipoint],                            fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data(), 
			      GetFixWidthCol(RoundToSignificantFigures(((TGraphErrors*)fGraphStat)->GetEY()[ipoint], fPrecision), 10).Data());

	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("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-\t\tdysyst+\t\tdysyst-")); 
	else 
	  fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-"));       	
      }
      // Skip empty bins
      if(!fHistStat->GetBinContent(ibin)) continue;
      TObjString * line = new TObjString;      
      if(fHistSyst) {
	
	line->String().Form("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s",
			    //	line->String().Form("%10s %10s %10s %10s %10s %10s %10s %10s", 
			    GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin)/2,         12).Data(), 			    
			    GetFixWidthCol(fHistStat->GetBinLowEdge(ibin),                                        12).Data(), 
			    GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin),           12).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinContent(ibin),fPrecision),  12).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin),  fPrecision),  12).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin),  fPrecision),  12).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistSyst->GetBinError(ibin),  fPrecision),  12).Data(),
			    GetFixWidthCol(RoundToSignificantFigures(fHistSyst->GetBinError(ibin),  fPrecision),  12).Data());

	fHEPDataFileLines->Add(line);
      } else {
	line->String().Form("%s\t%s\t%s\t%s\t%s\t%s", 
			    GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin)/2,         10).Data(), 			    
			    GetFixWidthCol(fHistStat->GetBinLowEdge(ibin),                                        10).Data(), 
			    GetFixWidthCol(fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin),           10).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinContent(ibin),fPrecision),  10).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin),  fPrecision),  10).Data(), 
			    GetFixWidthCol(RoundToSignificantFigures(fHistStat->GetBinError(ibin),  fPrecision),  10).Data()); 
	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();
}


Double_t AliHEPDataParser::RoundToSignificantFigures(double num, int n) {
  // Rounds num to n significant digits.
  // Recipe from :http://stackoverflow.com/questions/202302/rounding-to-an-arbitrary-number-of-significant-digits
  // Basically the log is used to determine the number of leading 0s, than convert to an integer by multipliing by the expo, 
  // round the integer and shift back.
  if(num == 0) {
    return 0;
  }

  Double_t d = TMath::Ceil(TMath::Log10(num < 0 ? -num: num));
  Int_t power = n - (int) d;

  Double_t magnitude = TMath::Power(10, power);
  Long_t shifted = TMath::Nint(num*magnitude);
  return shifted/magnitude;

}

TString AliHEPDataParser::GetFixWidthCol(Double_t number, Int_t width) {

  // Formats a column to fixed width
  TString col;
  char format[100];
  snprintf(format,100,"%%%d#g", fPrecision);
  col.Form(format, number);
  if(col.Length()>width) AliError("larger than width, cannot align!");

  if(col.Contains("e"))
    while (col.Length() < width) col.Append(" ");
  else
    while (col.Length() < width) col.Append("0");
  
  return col;
}
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"
d5e9ea5e	25	#include "TMath.h"
a5ae512a	26	#include <fstream>
	27	#include <iostream>
	28	using namespace std;
	29
	30	ClassImp(AliHEPDataParser)
	31
d5e9ea5e	32	AliHEPDataParser::AliHEPDataParser() : TObject(), fHistStat(0), fHistSyst(0), fGraphStat(0), fGraphSyst(0), fHEPDataFileLines(0), fValueName(""), fXaxisName(""), fTitle(""), fReaction(""), fEnergy(""), fRapidityRange(""), fPrecision(5)
a5ae512a	33	{
	34	// default ctor
	35
	36	}
	37
d5e9ea5e	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)
a5ae512a	39	{
	40	//ctor
	41	fHistStat = hStat;
	42	fHistSyst = hSyst;
	43	fHEPDataFileLines = new TObjArray;
	44
	45	}
d5e9ea5e	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)
a5ae512a	47	{
	48	// ctor
	49	fGraphStat = grStat;
	50	fGraphSyst = grSyst;
	51	fHEPDataFileLines = new TObjArray;
	52	}
	53
4bc92dcb	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)
a5ae512a	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)) {
d5e9ea5e	65	TObjArray * tokens = line.Tokenize(" \t");
4bc92dcb	66	if(!tokens) continue;
d5e9ea5e	67	if(! ((TObjString*) tokens->At(0))->String().Atof()) {
d5e9ea5e	68	//The first column is not a number: those are the headers: skip!
a5ae512a	69	delete tokens;
d5e9ea5e	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;
a5ae512a	77	return;
d5e9ea5e	78	//continue;
a5ae512a	79	}
d5e9ea5e	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;
a5ae512a	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,
d5e9ea5e	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());
a5ae512a	111	ipoints++;
09d5920f	112	delete tokens;
a5ae512a	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) {
d5e9ea5e	129
a5ae512a	130	// Fills fHEPDataFileLines and saves its content to a file
a5ae512a	131	if(!fHEPDataFileLines) fHEPDataFileLines = new TObjArray;
d5e9ea5e	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
a5ae512a	142	if(trueUseGraphFalesUseHisto) {
d5e9ea5e	143	AliWarning("Graph saving not thoroughly tested!!");
a5ae512a	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) {
d5e9ea5e	157	fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat\t+syst\t-syst"));
a5ae512a	158	}
a5ae512a	159	else {
d5e9ea5e	160	fHEPDataFileLines->Add(new TObjString("x\txlow\txhigh\t+stat\t-stat"));
a5ae512a	161	}
	162	}
	163	// Skip empty bins
	164	if(!fGraphStat->GetY()[ipoint]) continue;
d5e9ea5e	165	TObjString * line = new TObjString;
a5ae512a	166	if(fGraphSyst) {
a5ae512a	167	if (asym)
d5e9ea5e	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());
a5ae512a	178	else
d5e9ea5e	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));
a5ae512a	193
	194	fHEPDataFileLines->Add(line);
	195	} else {
	196	if (asym)
d5e9ea5e	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());
a5ae512a	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)
d5e9ea5e	229	fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-\t\tdysyst+\t\tdysyst-"));
a5ae512a	230	else
d5e9ea5e	231	fHEPDataFileLines->Add(new TObjString("x\t\txlow\t\txhigh\t\ty\t\tdystat+\t\tdystat-"));
a5ae512a	232	}
	233	// Skip empty bins
	234	if(!fHistStat->GetBinContent(ibin)) continue;
	235	TObjString * line = new TObjString;
	236	if(fHistSyst) {
d5e9ea5e	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
a5ae512a	249	fHEPDataFileLines->Add(line);
a5ae512a	250	} else {
d5e9ea5e	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());
a5ae512a	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
d5e9ea5e	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	}