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