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 "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 }