File indexing completed on 2024-04-06 11:56:33
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063 #include <fstream>
0064 #include <vector>
0065 #include <utility> // for std::pair
0066 #include <iostream>
0067
0068 #include <TROOT.h>
0069 #include <TFile.h>
0070 #include <TDirectory.h>
0071 #include <TError.h>
0072 #include <TH1.h>
0073 #include <TGraph.h>
0074 #include <TGraphErrors.h>
0075 #include <TCanvas.h>
0076 #include <TMath.h>
0077 #include <TPaveText.h>
0078
0079
0080 void readPedeHists(Option_t *option = "", const char *txtFile = "millepede.his");
0081
0082
0083
0084
0085
0086 class ReadPedeHists
0087 {
0088 public:
0089 explicit ReadPedeHists(const char *txtFile);
0090 ~ReadPedeHists() {}
0091 void Draw();
0092 void Write(const char *rootFileName);
0093 void Print(const char *printFileName);
0094
0095 private:
0096 template<class T>
0097 bool readBasic(std::ifstream &aStream, T &outValue);
0098 template<class T>
0099 bool readBasicVector(std::ifstream &aStream, std::vector<T> &outVector);
0100 bool proceedTo(std::ifstream &aStream, const TString &pattern);
0101 void readNumVersTypeTitle(std::ifstream &file, Int_t &num, Int_t &version, Int_t &type,
0102 TString &title);
0103 TH1 *readNextHist(std::ifstream &file);
0104 std::pair<TGraph *,Option_t*> readNextGraph(std::ifstream &file);
0105 bool readNext(std::ifstream &file, TH1 *&hist, std::pair<TGraph*,Option_t*> &graphOpt);
0106 void read(std::ifstream &stream);
0107
0108 std::ifstream theStream;
0109 std::vector<TH1*> theHists;
0110 std::vector<std::pair<TGraph*, Option_t*> > theGraphOpts;
0111 std::vector<TCanvas*> theCanvases;
0112 };
0113
0114
0115 ReadPedeHists::ReadPedeHists(const char *txtFile) :
0116 theStream(txtFile, ios::in)
0117 {
0118 if (!theStream.is_open()) {
0119 ::Error("ReadPedeHists::ReadPedeHists", "file %s could not be opened", txtFile);
0120 } else {
0121 this->read(theStream);
0122 }
0123 }
0124
0125
0126
0127 template<class T>
0128 bool ReadPedeHists::readBasic(std::ifstream &aStream, T &outValue)
0129 {
0130 while (true) {
0131 const int aChar = aStream.get();
0132 if (!aStream.good()) return false;
0133
0134 switch(aChar) {
0135 case ' ':
0136 case '\t':
0137 case '\n':
0138 if (aStream.eof()) return false;
0139 continue;
0140 default:
0141 aStream.unget();
0142 aStream >> outValue;
0143 if (aStream.fail()) {
0144 aStream.clear();
0145 return false;
0146 } else {
0147 return true;
0148 }
0149 }
0150 }
0151
0152 ::Error("ReadPedeHists::readBasic", "Should never come here!");
0153 return false;
0154 }
0155
0156
0157 template<class T>
0158 bool ReadPedeHists::readBasicVector(std::ifstream &aStream, std::vector<T> &outVector)
0159 {
0160
0161 for (unsigned int i = 0; i < outVector.size(); ++i) {
0162 if (!readBasic(aStream, outVector[i])) return false;
0163 }
0164
0165 return true;
0166 }
0167
0168
0169 bool ReadPedeHists::proceedTo(std::ifstream &aStream, const TString &pattern)
0170 {
0171 if (pattern.IsNull()) return true;
0172 const char *method = "ReadPedeHists::proceedTo";
0173
0174 TString line;
0175 do {
0176 line.ReadLine(aStream);
0177 if (line.Contains(pattern)) {
0178 line.ReplaceAll(pattern, "");
0179 line.ReplaceAll(" ", "");
0180 if (!line.IsNull()) {
0181 ::Warning(method, "line contains also '%s'", line.Data());
0182 }
0183 return true;
0184 } else {
0185 ::Warning(method, "skipping line '%s'", line.Data());
0186 }
0187 } while (!aStream.eof());
0188
0189 ::Error(method, "pattern '%s' not found", pattern.Data());
0190 return false;
0191 }
0192
0193
0194 void ReadPedeHists::readNumVersTypeTitle(std::ifstream &file, Int_t &num,
0195 Int_t &version, Int_t &type, TString &title)
0196 {
0197 std::string key;
0198
0199 const char *method = "ReadPedeHists::readNumVersTypeTitle";
0200 if (!readBasic(file, num)) {
0201 ::Error(method, "failed reading hist number");
0202 }
0203
0204 if (!readBasic(file, key) || key != "version") {
0205 ::Error(method, "expect key 'version', got '%s'", key.c_str());
0206 }
0207 if (!readBasic(file, version)) {
0208 ::Error(method, "failed reading version");
0209 }
0210
0211 if (!readBasic(file, key) || key != "type") {
0212 ::Error(method, "expect key 'type', got '%s'", key.c_str());
0213 }
0214 if (!readBasic(file, type)) ::Error(method, "failed reading type");
0215
0216 title.ReadLine(file);
0217 Ssiz_t len = title.Length();
0218 while (len != kNPOS && len > 0 && title[--len] == ' ') {}
0219 title.Resize(len+1);
0220 title += Form(" (version %d)", version);
0221 }
0222
0223
0224 TH1 *ReadPedeHists::readNextHist(std::ifstream &file)
0225 {
0226
0227
0228
0229 Int_t num = -1;
0230 Int_t version = -1;
0231 Int_t type = -1;
0232 TString title;
0233
0234 const char *method = "ReadPedeHists::readNextHist";
0235
0236 readNumVersTypeTitle(file, num, version, type, title);
0237 if (num == -1 || version == -1 || type == -1) {
0238 ::Error(method, "Problems reading hist number, version or type, so skip it.");
0239 proceedTo(file, "end of histogram");
0240 return 0;
0241 }
0242
0243
0244
0245
0246
0247 std::vector<Float_t> nBinsUpLow(3, -1.);
0248 std::vector<Int_t> underInOver(3, -1);
0249 std::vector<Float_t> binContent;
0250 Float_t min = 0., max = 0., mean = 0., sigma = 0.;
0251
0252 std::string key;
0253 while (readBasic(file, key)) {
0254 if (key == "bins,") {
0255
0256 if (!readBasic(file, key) || key != "limits") {
0257 ::Error(method, "expect key 'limits', got (%s)", key.c_str());
0258 } else if (!readBasicVector(file, nBinsUpLow)) {
0259 ::Error(method, "failed reading nBins, xLow, xUp (%f %f %f)",
0260 nBinsUpLow[0], nBinsUpLow[1], nBinsUpLow[2]);
0261 } else {
0262 binContent.resize(static_cast<unsigned int>(nBinsUpLow[0]));
0263 }
0264 } else if (key == "out-low") {
0265
0266 if (!readBasic(file, key) || key != "inside"
0267 || !readBasic(file, key) || key != "out-high") {
0268 ::Error(method, "expected keys 'inside' and 'out-high', got (%s)", key.c_str());
0269 } else if (!readBasicVector(file, underInOver) || underInOver[0] == -1
0270 || underInOver[1] == -1 || underInOver[2] == -1) {
0271 ::Error(method, "failed reading under-, 'in-' and overflow (%d %d %d)",
0272 underInOver[0], underInOver[1], underInOver[2]);
0273 }
0274 } else if (key == "bincontent") {
0275
0276 if (nBinsUpLow[0] == -1.) {
0277 ::Error(method, "n(bins) (key 'bins') not yet set, bin content cannot be read");
0278 } else if (!readBasicVector(file, binContent)) {
0279 ::Error(method, "failed reading bincontent ");
0280 }
0281 } else if (key == "minmax") {
0282
0283 if (!readBasic(file, min) || !readBasic(file, max)) {
0284 ::Error(method, "failed reading min or max (%f %f)", min, max);
0285 }
0286 } else if (key == "meansigma") {
0287
0288 if (!readBasic(file, mean) || !readBasic(file, sigma)) {
0289 ::Error(method, "failed reading mean or sigma (%f %f)", mean, sigma);
0290 }
0291 } else if (key == "end") {
0292
0293 proceedTo(file, "of histogram");
0294 break;
0295 } else {
0296 ::Error(method, "unknown key '%s', try next word", key.c_str());
0297 }
0298 }
0299
0300
0301 if (nBinsUpLow[1] == nBinsUpLow[2]) {
0302 nBinsUpLow[2] = nBinsUpLow[1] + 1.;
0303 ::Error(method, "Hist %d (version %d): same upper and lower edge (%f), set upper %f.",
0304 num, version, nBinsUpLow[1], nBinsUpLow[2]);
0305 }
0306 TH1 *h = new TH1F(Form("hist%d_version%d", num, version), title,
0307 binContent.size(), nBinsUpLow[1], nBinsUpLow[2]);
0308 h->SetBinContent(0, underInOver[0]);
0309 for (UInt_t iBin = 1; iBin <= binContent.size(); ++iBin) {
0310 h->SetBinContent(iBin, binContent[iBin - 1]);
0311 }
0312 h->SetBinContent(binContent.size() + 1, underInOver[2]);
0313 h->SetEntries(underInOver[0] + underInOver[1] + underInOver[2]);
0314
0315 if (type == 2) {
0316
0317 h->SetXTitle("log_{10}");
0318 } else if (type != 1) {
0319 ::Warning(method, "Unknown histogram type %d.", type);
0320 }
0321
0322 if (mean || sigma) {
0323 Double_t stats[11] = {0.};
0324 h->GetStats(stats);
0325 stats[0] = stats[1] = h->GetEntries();
0326 stats[2] = mean * stats[0];
0327 stats[3] = (sigma * sigma + mean * mean) * stats[0];
0328 h->PutStats(stats);
0329 }
0330 if (min || max) {
0331 TPaveText *text = new TPaveText(.175, .675, .45, .875, "NDC");
0332 text->AddText(Form("min = %g", min));
0333 text->AddText(Form("max = %g", max));
0334 text->SetTextAlign(12);
0335 text->SetBorderSize(1);
0336 h->GetListOfFunctions()->Add(text);
0337 }
0338
0339 return h;
0340 }
0341
0342
0343 std::pair<TGraph*, Option_t*> ReadPedeHists::readNextGraph(std::ifstream &file)
0344 {
0345
0346
0347
0348 TGraph *graph = 0;
0349 Option_t *drawOpt = 0;
0350
0351
0352 Int_t num = -1;
0353 Int_t version = -1;
0354 Int_t type = -1;
0355 TString title;
0356
0357 const char *method = "ReadPedeHists::readNextGraph";
0358
0359 readNumVersTypeTitle(file, num, version, type, title);
0360 if (num == -1 || version == -1 || type == -1) {
0361 ::Error(method, "Problems reading graph number, version or type, so skip it.");
0362 proceedTo(file, "end of xy-data");
0363 return std::make_pair(graph, drawOpt);
0364 }
0365
0366
0367
0368
0369
0370 if (type < 1 || type > 5) {
0371 ::Error(method, "Unknown xy-data type %d, so skip graph.", type);
0372 proceedTo(file, "end of xy-data");
0373 }
0374
0375
0376 UInt_t numPoints = 0;
0377 std::vector<Float_t> content;
0378
0379 std::string key;
0380 while (readBasic(file, key)) {
0381 if (key == "stored") {
0382 if (!readBasic(file, key) || key != "not-stored") {
0383 ::Error(method, "expected key 'not-stored', got '%s'", key.c_str());
0384 } else if (!readBasic(file, numPoints)) {
0385 ::Error(method, "failed reading number of points (%d)", numPoints);
0386 }
0387 } else if (key == "x-y") {
0388 if (type < 1 || type > 3) {
0389 ::Error(method, "expect key x-y-dx-dy for type %d, found x-y", type);
0390 }
0391 content.resize(numPoints * 2);
0392 if (!readBasicVector(file, content) || !numPoints) {
0393 ::Error(method, "failed reading x-y content%s",
0394 (!numPoints ? " since n(points) (key 'stored') not yet set" : ""));
0395 }
0396 } else if (key == "x-y-dx-dy") {
0397 if (type < 4 || type > 5) {
0398 ::Error(method, "expect key x-y for type %d, found x-y-dx-dy", type);
0399 }
0400 content.resize(numPoints * 4);
0401 if (!readBasicVector(file, content) || !numPoints) {
0402 ::Error(method, "failed reading x-y-dx-dy content%s",
0403 (!numPoints ? " since n(points) (key 'stored') not yet set" : ""));
0404 }
0405 } else if (key == "end") {
0406 proceedTo(file, "of xy-data");
0407 break;
0408 } else {
0409 break;
0410 ::Error(method, "unknown key '%s', try next word", key.c_str());
0411 }
0412 }
0413
0414
0415 if (type == 4 || type == 5) {
0416 TGraphErrors *graphE = new TGraphErrors(numPoints);
0417 for (unsigned int i = 0; i < content.size(); i += 4) {
0418 graphE->SetPoint (i/4, content[i] , content[i+1]);
0419 graphE->SetPointError(i/4, content[i+2], content[i+3]);
0420 }
0421 drawOpt = "AP";
0422 graph = graphE;
0423 } else if (type >= 1 && type <= 3) {
0424 graph = new TGraph(numPoints);
0425 for (unsigned int i = 0; i < content.size(); i += 2) {
0426 graph->SetPoint(i/2, content[i], content[i+1]);
0427 }
0428 if (type == 1) {
0429 drawOpt = "AP";
0430 } else if (type == 2) {
0431 drawOpt = "AL";
0432 } else if (type == 3) {
0433 drawOpt = "ALP";
0434 }
0435 if (TString(drawOpt).Contains("P")) graph->SetMarkerStyle(20);
0436 }
0437
0438 if (graph) graph->SetNameTitle(Form("graph%d_version%d", num, version), title);
0439 return std::make_pair(graph, drawOpt);
0440 }
0441
0442
0443 bool ReadPedeHists::readNext(std::ifstream &file, TH1 *&hist,
0444 std::pair<TGraph*,Option_t*> &graphOpt)
0445 {
0446 hist = 0;
0447 graphOpt.first = 0;
0448 graphOpt.second = 0;
0449
0450 TString type;
0451 while (true) {
0452 if (file.eof()) break;
0453 file >> type;
0454 if (file.fail() || (type != "Histogram" && type != "XY-Data")) {
0455 TString line;
0456 line.ReadLine(file);
0457 if (line != "" && line.Length() != line.CountChar(' ')) {
0458 ::Error("ReadPedeHists::readNext",
0459 "Expect 'Histogram' or 'XY-Data', but failed, line is '%s'",
0460 line.Data());
0461 if (proceedTo(file, "end of")) line.ReadLine(file);
0462 }
0463 }
0464
0465 if (type == "Histogram") hist = readNextHist(file);
0466 if (type == "XY-Data") graphOpt = readNextGraph(file);
0467 if (hist || graphOpt.first) break;
0468 }
0469
0470 return (hist || graphOpt.first);
0471 }
0472
0473
0474 void ReadPedeHists::read(std::ifstream &file)
0475 {
0476 theHists.clear();
0477 theGraphOpts.clear();
0478
0479 TH1 *hist = 0;
0480 std::pair<TGraph*, Option_t*> graphOpt(0,0);
0481 while (readNext(file, hist, graphOpt)) {
0482 if (hist) theHists.push_back(hist);
0483 if (graphOpt.first) theGraphOpts.push_back(graphOpt);
0484 }
0485 }
0486
0487
0488 void ReadPedeHists::Draw()
0489 {
0490 theCanvases.clear();
0491
0492 const Int_t nHistX = 3;
0493 const Int_t nPixelX = 700;
0494 const Int_t nHistY = 2;
0495 const Int_t nPixelY = 500;
0496 Int_t last = nHistX * nHistY;
0497 unsigned int iH = 0;
0498
0499 while (iH < theHists.size()) {
0500 if (last >= nHistX * nHistY) {
0501 unsigned int canCorner = theCanvases.size() * 20;
0502 theCanvases.push_back(new TCanvas(Form("hists%d", iH), "",
0503 canCorner, canCorner, nPixelX, nPixelY));
0504 theCanvases.back()->Divide(nHistX, nHistY);
0505 last = 0;
0506 }
0507 theCanvases.back()->cd(++last);
0508 theHists[iH]->Draw();
0509 ++iH;
0510 }
0511
0512 last = nHistX * nHistY;
0513 iH = 0;
0514 while (iH < theGraphOpts.size()) {
0515 if (last >= nHistX * nHistY) {
0516 unsigned int canCorner = theCanvases.size() * 20;
0517 theCanvases.push_back(new TCanvas(Form("graphs%d", iH), "",
0518 canCorner, canCorner, nPixelX, nPixelY));
0519 theCanvases.back()->Divide(nHistX, nHistY);
0520 last = 0;
0521 }
0522 theCanvases.back()->cd(++last);
0523 theGraphOpts[iH].first->Draw(theGraphOpts[iH].second);
0524 ++iH;
0525 }
0526 }
0527
0528
0529 void ReadPedeHists::Print(const char *printFileName)
0530 {
0531 std::vector<TCanvas*>::iterator iC = theCanvases.begin(), iE = theCanvases.end();
0532 if (iC == iE) return;
0533
0534 theCanvases.front()->Print(Form("%s[", printFileName));
0535 while(iC != iE) {
0536 (*iC)->Print(printFileName);
0537 ++iC;
0538 }
0539 theCanvases.front()->Print(Form("%s]", printFileName));
0540
0541 }
0542
0543
0544 void ReadPedeHists::Write(const char *rootFileName)
0545 {
0546 if (theHists.empty() && theGraphOpts.empty()) return;
0547
0548 ::Info("ReadPedeHists::Write", "(Re-)Creating ROOT file %s.", rootFileName);
0549
0550 TDirectory *oldDir = gDirectory;
0551 TFile *rootFile = TFile::Open(rootFileName, "RECREATE");
0552
0553 for (std::vector<TH1*>::iterator iH = theHists.begin(), iE = theHists.end();
0554 iH != iE; ++iH) {
0555 (*iH)->Write();
0556 }
0557
0558 for (std::vector<std::pair<TGraph*,Option_t*> >::iterator iG = theGraphOpts.begin(),
0559 iE = theGraphOpts.end(); iG != iE; ++iG) {
0560 (*iG).first->Write();
0561 }
0562
0563 delete rootFile;
0564 oldDir->cd();
0565 }
0566
0567
0568
0569
0570 void readPedeHists(Option_t *option, const char *txtFile)
0571 {
0572 ReadPedeHists reader(txtFile);
0573 TString opt(option);
0574 opt.ToLower();
0575
0576 const bool oldBatch = gROOT->IsBatch();
0577 if (opt.Contains("nodraw")) {
0578 opt.ReplaceAll("nodraw", "");
0579 gROOT->SetBatch(true);
0580 }
0581
0582 reader.Draw();
0583
0584 if (opt.Contains("print")) {
0585 opt.ReplaceAll("print", "");
0586 reader.Print(TString(Form("%s.pdf", txtFile)));
0587 }
0588
0589 if (opt.Contains("write")) {
0590 opt.ReplaceAll("write", "");
0591 reader.Write(TString(Form("%s.root", txtFile)));
0592 }
0593
0594 gROOT->SetBatch(oldBatch);
0595 opt.ReplaceAll(" ", "");
0596 if (!opt.IsNull()) {
0597 ::Warning("readPedeHists", "Unknown option '%s', know 'nodraw', 'print' and 'write'.",
0598 opt.Data());
0599 }
0600 }