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()) {delete tokens; 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     delete tokens;
  99   }
 100   infile.close();
 101
 102
 103 }
 104
 105 AliHEPDataParser::~AliHEPDataParser(){
 106   // dtor
 107   if(fHistStat) delete fHistStat;
 108   if(fHistSyst) delete fHistSyst;
 109   if(fGraphStat) delete fGraphStat;
 110   if(fGraphSyst) delete fGraphSyst;
 111   if(fHEPDataFileLines) delete fHEPDataFileLines;
 112 }
 113
 114 void AliHEPDataParser::SaveHEPDataFile(const char * hepfileName, Bool_t trueUseGraphFalesUseHisto) {
 115   // Fills fHEPDataFileLines and saves its content to a file
 116   if(!fHEPDataFileLines)   fHEPDataFileLines = new TObjArray;
 117   if(trueUseGraphFalesUseHisto) {
 118     if(!fGraphStat) {
 119       AliError("Graph not set");
 120       return;
 121     }
 122     Bool_t asym = kFALSE; // check if this has asymmetric errors
 123     if       (!strcmp(fGraphStat->ClassName(), "TGraphErrors")) asym = kFALSE;
 124     else if  (!strcmp(fGraphStat->ClassName(), "TGraphAsymmErrors")) asym = kTRUE;
 125     else     {AliError("Unsupported graph type"); return;}
 126     Int_t npoint = fGraphStat->GetN();
 127     if(asym) AliInfo("Assymmetric errors");
 128     for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
 129       if(ipoint == 0) {
 130         if(fGraphSyst) {
 131           if (asym)
 132             fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +stat -stat +syst -syst", fValueName.Data())));
 133           else
 134             fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +-stat +-syst", fValueName.Data())));
 135         }
 136         else {
 137           if(asym)
 138             fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +stat -stat", fValueName.Data())));
 139           else
 140             fHEPDataFileLines->Add(new TObjString(Form("BinCenter %s +-stat", fValueName.Data())));
 141         }
 142       }
 143       // Skip empty bins
 144       if(!fGraphStat->GetY()[ipoint]) continue;
 145       TObjString * line = new TObjString;
 146       if(fGraphSyst) {
 147         if (asym)
 148           line->String().Form("%f %f +%f -%f +%f -%f",
 149                               fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
 150                               ((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint],
 151                               ((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint],
 152                               ((TGraphAsymmErrors*)fGraphSyst)->GetEYhigh()[ipoint],
 153                               ((TGraphAsymmErrors*)fGraphSyst)->GetEYlow()[ipoint]);
 154         else
 155           line->String().Form("%f %f +-%f +-%f",
 156                               fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
 157                               ((TGraphErrors*)fGraphStat)->GetEY()[ipoint],
 158                               ((TGraphErrors*)fGraphSyst)->GetEY()[ipoint]);
 159
 160         fHEPDataFileLines->Add(line);
 161       } else {
 162         if (asym)
 163           line->String().Form("%f %f +%f -%f",
 164                               fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
 165                               ((TGraphAsymmErrors*)fGraphStat)->GetEYhigh()[ipoint], ((TGraphAsymmErrors*)fGraphStat)->GetEYlow()[ipoint]);
 166         else {
 167           line->String().Form("%f %f +-%f",
 168                               fGraphStat->GetX()[ipoint], fGraphStat->GetY()[ipoint],
 169                               ((TGraphErrors*)fGraphStat)->GetEY()[ipoint]);
 170         }
 171
 172         fHEPDataFileLines->Add(line);
 173       }
 174     }
 175   }
 176   else {
 177     if(!fHistStat) {
 178       AliError("Hist not set");
 179       return;
 180     }
 181     Int_t nbin = fHistStat->GetNbinsX();
 182
 183     for(Int_t ibin = 1; ibin <= nbin; ibin++){
 184       if(ibin == 1) {
 185         if(fHistSyst)
 186           fHEPDataFileLines->Add(new TObjString(Form("BinLow BinHigh %s +-stat +-syst", fValueName.Data())));
 187         else
 188           fHEPDataFileLines->Add(new TObjString(Form("BinLow BinHigh %s +-stat", fValueName.Data())));
 189       }
 190       // Skip empty bins
 191       if(!fHistStat->GetBinContent(ibin)) continue;
 192       TObjString * line = new TObjString;
 193       if(fHistSyst) {
 194         line->String().Form("%f %f %f +-%f +-%f",
 195                             fHistStat->GetBinLowEdge(ibin),
 196                             fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin),
 197                             fHistStat->GetBinContent(ibin), fHistStat->GetBinError(ibin), fHistSyst->GetBinError(ibin));
 198         fHEPDataFileLines->Add(line);
 199       } else {
 200         line->String().Form("%f %f %f +-%f",
 201                             fHistStat->GetBinLowEdge(ibin),
 202                             fHistStat->GetBinLowEdge(ibin)+fHistStat->GetBinWidth(ibin),
 203                             fHistStat->GetBinContent(ibin), fHistStat->GetBinError(ibin));
 204         fHEPDataFileLines->Add(line);
 205       }
 206       //      delete line;
 207     }
 208   }
 209
 210   TIterator * lineIter = fHEPDataFileLines->MakeIterator();
 211   TObjString * obj = 0;
 212   ofstream outfile;
 213   outfile.open (hepfileName);
 214   cout << "Saving HEP File " << hepfileName << endl;
 215
 216   while ((obj = (TObjString*) lineIter->Next())) {
 217     cout << obj->String().Data() << endl;
 218     outfile << obj->String().Data() << endl;
 219   }
 220   outfile.close();
 221 }
 222
 223