Millepede-II V04-17-03
readPedeHists.C
Go to the documentation of this file.
1
87#include <fstream>
88#include <vector>
89#include <utility> // for std::pair
90#include <iostream>
91
92#include <TROOT.h>
93#include <TFile.h>
94#include <TDirectory.h>
95#include <TError.h>
96#include <TH1.h>
97#include <TGraph.h>
98#include <TGraphErrors.h>
99#include <TCanvas.h>
100#include <TMath.h>
101#include <TPaveText.h>
102
103//__________________________________________________________________________
104void readPedeHists(Option_t *option = "", const char *txtFile = "millepede.his");//"nodraw" skips drawing, "print" produces ps file, "write" writes object in ROOT file
105
106//__________________________________________________________________________
107//__________________________________________________________________________
108//__________________________________________________________________________
109
111{
112 public:
113 explicit ReadPedeHists(const char *txtFile);
114 ~ReadPedeHists() {} // non virtual: do not inherit from this class
115 void Draw(); // draw hists and graphs read from file
116 void Write(const char *rootFileName); // write hists/graphs to file
117 void Print(const char *printFileName); // if drawn, print into file, e.g. ps
118
119 private:
120 template<class T>
121 bool readBasic(std::ifstream &aStream, T &outValue);
122 template<class T>
123 bool readBasicVector(std::ifstream &aStream, std::vector<T> &outVector);
124 bool proceedTo(std::ifstream &aStream, const TString &pattern);
125 void readNumVersTypeTitle(std::ifstream &file, Int_t &num, Int_t &version, Int_t &type,
126 TString &title);
127 TH1 *readNextHist(std::ifstream &file);
128 std::pair<TGraph *,Option_t*> readNextGraph(std::ifstream &file);
129 bool readNext(std::ifstream &file, TH1 *&hist, std::pair<TGraph*,Option_t*> &graphOpt);
130 void read(std::ifstream &stream);
131
132 std::ifstream theStream;
133 std::vector<TH1*> theHists;
134 std::vector<std::pair<TGraph*, Option_t*> > theGraphOpts;
135 std::vector<TCanvas*> theCanvases;
136};
137
138//__________________________________________________________________________
139ReadPedeHists::ReadPedeHists(const char *txtFile) :
140 theStream(txtFile, ios::in)
141{
142 if (!theStream.is_open()) {
143 ::Error("ReadPedeHists::ReadPedeHists", "file %s could not be opened", txtFile);
144 } else {
145 this->read(theStream);
146 }
147}
148
149
150//__________________________________________________________________________
151template<class T>
152bool ReadPedeHists::readBasic(std::ifstream &aStream, T &outValue)
153{
154 while (true) {
155 const int aChar = aStream.get();
156 if (!aStream.good()) return false;
157
158 switch(aChar) {
159 case ' ':
160 case '\t':
161 case '\n':
162 if (aStream.eof()) return false;
163 continue; // to next character
164 default:
165 aStream.unget();
166 aStream >> outValue;
167 if (aStream.fail()) {// not correct type 'T' (!aStream.good() is true also in case of EOF)
168 aStream.clear();
169 return false;
170 } else {
171 return true;
172 }
173 } // switch
174 } // while
175
176 ::Error("ReadPedeHists::readBasic", "Should never come here!");
177 return false;
178}
179
180//__________________________________________________________________________
181template<class T>
182bool ReadPedeHists::readBasicVector(std::ifstream &aStream, std::vector<T> &outVector)
183{
184 // vector must have desired size
185 for (unsigned int i = 0; i < outVector.size(); ++i) {
186 if (!readBasic(aStream, outVector[i])) return false;
187 }
188
189 return true;
190}
191
192//__________________________________________________________________________
193bool ReadPedeHists::proceedTo(std::ifstream &aStream, const TString &pattern)
194{
195 if (pattern.IsNull()) return true;
196 const char *method = "ReadPedeHists::proceedTo";
197
198 TString line;
199 do {
200 line.ReadLine(aStream);
201 if (line.Contains(pattern)) {
202 line.ReplaceAll(pattern, "");
203 line.ReplaceAll(" ", "");
204 if (!line.IsNull()) {
205 ::Warning(method, "line contains also '%s'", line.Data());
206 }
207 return true;
208 } else {
209 ::Warning(method, "skipping line '%s'", line.Data());
210 }
211 } while (!aStream.eof());
212
213 ::Error(method, "pattern '%s' not found", pattern.Data());
214 return false; // did not find pattern
215}
216
217//__________________________________________________________________________
218void ReadPedeHists::readNumVersTypeTitle(std::ifstream &file, Int_t &num,
219 Int_t &version, Int_t &type, TString &title)
220{
221 std::string key; // key word
222
223 const char *method = "ReadPedeHists::readNumVersTypeTitle";
224 if (!readBasic(file, num)) {
225 ::Error(method, "failed reading hist number");
226 }
227
228 if (!readBasic(file, key) || key != "version") {
229 ::Error(method, "expect key 'version', got '%s'", key.c_str());
230 }
231 if (!readBasic(file, version)) {
232 ::Error(method, "failed reading version");
233 }
234
235 if (!readBasic(file, key) || key != "type") {
236 ::Error(method, "expect key 'type', got '%s'", key.c_str());
237 }
238 if (!readBasic(file, type)) ::Error(method, "failed reading type");
239
240 title.ReadLine(file); // Title is a full line without key after the type!
241 Ssiz_t len = title.Length();
242 while (len != kNPOS && len > 0 && title[--len] == ' ') {} // empty loop
243 title.Resize(len+1); // remove trailing space
244 title += Form(" (version %d)", version);
245}
246
247//__________________________________________________________________________
248TH1 *ReadPedeHists::readNextHist(std::ifstream &file)
249{
250 // Key 'Histogram' assumed to be already read!
251
252 // Until histogram title we have a fixed order to read in these numbers:
253 Int_t num = -1; // hist number
254 Int_t version = -1; // version (is it iteration?)
255 Int_t type = -1; // type, e.g. x-axis in log scale
256 TString title; // skip spaces??
257
258 const char *method = "ReadPedeHists::readNextHist"; // for errors/warnings
259
260 readNumVersTypeTitle(file, num, version, type, title);
261 if (num == -1 || version == -1 || type == -1) {
262 ::Error(method, "Problems reading hist number, version or type, so skip it.");
263 proceedTo(file, "end of histogram");
264 return 0;
265 }
266 // type 1: normal 1D histogram
267 // 2: 1D histogram with bins in log_10
268
269 // For the remaining information we accept any order, but there will be problems
270 // in case number of bins (key 'bins,') comes after 'bincontent'...
271 std::vector<Float_t> nBinsUpLow(3, -1.); // nBins (int...), lower and upper edge
272 std::vector<Int_t> underInOver(3, -1); // underflow : between lower/upper : overflow
273 std::vector<Float_t> binContent; // do not yet know length
274 Float_t min = 0., max = 0., mean = 0., sigma = 0.; // min/max of x-axis, mean/sigma of distrib.
275
276 std::string key; // key word
277 while (readBasic(file, key)) {
278 if (key == "bins,") {
279 // read nBins with borders
280 if (!readBasic(file, key) || key != "limits") {
281 ::Error(method, "expect key 'limits', got (%s)", key.c_str());
282 } else if (!readBasicVector(file, nBinsUpLow)) {
283 ::Error(method, "failed reading nBins, xLow, xUp (%f %f %f)",
284 nBinsUpLow[0], nBinsUpLow[1], nBinsUpLow[2]);
285 } else {
286 binContent.resize(static_cast<unsigned int>(nBinsUpLow[0]));
287 }
288 } else if (key == "out-low") {
289 // read under-/overflow with what is 'in between'
290 if (!readBasic(file, key) || key != "inside"
291 || !readBasic(file, key) || key != "out-high") {
292 ::Error(method, "expected keys 'inside' and 'out-high', got (%s)", key.c_str());
293 } else if (!readBasicVector(file, underInOver) || underInOver[0] == -1
294 || underInOver[1] == -1 || underInOver[2] == -1) {
295 ::Error(method, "failed reading under-, 'in-' and overflow (%d %d %d)",
296 underInOver[0], underInOver[1], underInOver[2]);
297 }
298 } else if (key == "bincontent") {
299 // read bin content - problem if lenght not yet set!
300 if (nBinsUpLow[0] == -1.) {
301 ::Error(method, "n(bins) (key 'bins') not yet set, bin content cannot be read");
302 } else if (!readBasicVector(file, binContent)) {
303 ::Error(method, "failed reading bincontent ");
304 }
305 } else if (key == "minmax") {
306 // read minimal and maximal x-value
307 if (!readBasic(file, min) || !readBasic(file, max)) {
308 ::Error(method, "failed reading min or max (%f %f)", min, max);
309 }
310 } else if (key == "meansigma") {
311 // read mean and sigma as calculated in pede
312 if (!readBasic(file, mean) || !readBasic(file, sigma)) {
313 ::Error(method, "failed reading mean or sigma (%f %f)", mean, sigma);
314 }
315 } else if (key == "end") {
316 // reached end - hopefully all has been read...
317 proceedTo(file, "of histogram");
318 break; // ...the while reading the next key
319 } else {
320 ::Error(method, "unknown key '%s', try next word", key.c_str());
321 }
322 }
323
324 // now create histogram
325 if (nBinsUpLow[1] == nBinsUpLow[2]) { // causes ROOT drawing errors
326 nBinsUpLow[2] = nBinsUpLow[1] + 1.;
327 ::Error(method, "Hist %d (version %d): same upper and lower edge (%f), set upper %f.",
328 num, version, nBinsUpLow[1], nBinsUpLow[2]);
329 }
330 TH1 *h = new TH1F(Form("hist%d_version%d", num, version), title,
331 binContent.size(), nBinsUpLow[1], nBinsUpLow[2]);
332 h->SetBinContent(0, underInOver[0]);
333 for (UInt_t iBin = 1; iBin <= binContent.size(); ++iBin) {
334 h->SetBinContent(iBin, binContent[iBin - 1]);
335 }
336 h->SetBinContent(binContent.size() + 1, underInOver[2]);
337 h->SetEntries(underInOver[0] + underInOver[1] + underInOver[2]);
338
339 if (type == 2) {
340 // could do more fancy stuff for nicer display...
341 h->SetXTitle("log_{10}");
342 } else if (type != 1) {
343 ::Warning(method, "Unknown histogram type %d.", type);
344 }
345
346 if (mean || sigma) { // overwrite ROOT's approximations from bin contents
347 Double_t stats[11] = {0.}; // no way to get this '11' from TH1... :-(
348 h->GetStats(stats);
349 stats[0] = stats[1] = h->GetEntries();// sum w and w^2
350 stats[2] = mean * stats[0]; // sum wx
351 stats[3] = (sigma * sigma + mean * mean) * stats[0]; // sum wx^2
352 h->PutStats(stats);
353 }
354 if (min || max) {
355 TPaveText *text = new TPaveText(.175, .675, .45, .875, "NDC");
356 text->AddText(Form("min = %g", min));
357 text->AddText(Form("max = %g", max));
358 text->SetTextAlign(12);
359 text->SetBorderSize(1);
360 h->GetListOfFunctions()->Add(text);// 'hack' to get it drawn with the hist
361 }
362
363 return h;
364}
365
366//__________________________________________________________________________
367std::pair<TGraph*, Option_t*> ReadPedeHists::readNextGraph(std::ifstream &file)
368{
369 // graph and drawing option...
370 // Key 'XY-Data' assumed to be already read!
371
372 TGraph *graph = 0;
373 Option_t *drawOpt = 0; // fine to use simple pointer since assigned only hardcoded strings
374
375 // Until graph title we have a fixed order to read in these numbers:
376 Int_t num = -1; // graph number
377 Int_t version = -1; // version (is it iteration?)
378 Int_t type = -1; // cf. below
379 TString title;
380
381 const char *method = "ReadPedeHists::readNextGraph"; // for errors/warnings
382
383 readNumVersTypeTitle(file, num, version, type, title);
384 if (num == -1 || version == -1 || type == -1) {
385 ::Error(method, "Problems reading graph number, version or type, so skip it.");
386 proceedTo(file, "end of xy-data");
387 return std::make_pair(graph, drawOpt);
388 }
389 // graph types: 1 dots (x,y)
390 // 2 polyline
391 // 3 dots and polyline
392 // 4 symbols with (x,y) and dx, dy
393 // 5 same as 5
394 if (type < 1 || type > 5) {
395 ::Error(method, "Unknown xy-data type %d, so skip graph.", type);
396 proceedTo(file, "end of xy-data");
397 }
398
399 // now read number of points and content
400 UInt_t numPoints = 0;
401 std::vector<Float_t> content; // do not yet know length (need two/four values per point!)
402
403 std::string key;
404 while (readBasic(file, key)) {
405 if (key == "stored") {
406 if (!readBasic(file, key) || key != "not-stored") {
407 ::Error(method, "expected key 'not-stored', got '%s'", key.c_str());
408 } else if (!readBasic(file, numPoints)) {
409 ::Error(method, "failed reading number of points (%d)", numPoints);
410 }
411 } else if (key == "x-y") {
412 if (type < 1 || type > 3) {
413 ::Error(method, "expect key x-y-dx-dy for type %d, found x-y", type);
414 }
415 content.resize(numPoints * 2);
416 if (!readBasicVector(file, content) || !numPoints) {
417 ::Error(method, "failed reading x-y content%s",
418 (!numPoints ? " since n(points) (key 'stored') not yet set" : ""));
419 }
420 } else if (key == "x-y-dx-dy") {
421 if (type < 4 || type > 5) {
422 ::Error(method, "expect key x-y for type %d, found x-y-dx-dy", type);
423 }
424 content.resize(numPoints * 4);
425 if (!readBasicVector(file, content) || !numPoints) {
426 ::Error(method, "failed reading x-y-dx-dy content%s",
427 (!numPoints ? " since n(points) (key 'stored') not yet set" : ""));
428 }
429 } else if (key == "end") {
430 proceedTo(file, "of xy-data");
431 break;
432 } else {
433 break;
434 ::Error(method, "unknown key '%s', try next word", key.c_str());
435 }
436 }
437
438 // now create TGraph(Error) and fill drawOpt
439 if (type == 4 || type == 5) {
440 TGraphErrors *graphE = new TGraphErrors(numPoints);
441 for (unsigned int i = 0; i < content.size(); i += 4) {
442 graphE->SetPoint (i/4, content[i] , content[i+1]);
443 graphE->SetPointError(i/4, content[i+2], content[i+3]);
444 }
445 drawOpt = "AP";
446 graph = graphE;
447 } else if (type >= 1 && type <= 3) {
448 graph = new TGraph(numPoints);
449 for (unsigned int i = 0; i < content.size(); i += 2) {
450 graph->SetPoint(i/2, content[i], content[i+1]);
451 }
452 if (type == 1) {
453 drawOpt = "AP";
454 } else if (type == 2) {
455 drawOpt = "AL";
456 } else if (type == 3) {
457 drawOpt = "ALP";
458 }
459 if (TString(drawOpt).Contains("P")) graph->SetMarkerStyle(20); //
460 } // 'else' not needed, tested above
461
462 if (graph) graph->SetNameTitle(Form("graph%d_version%d", num, version), title);
463 return std::make_pair(graph, drawOpt);
464}
465
466//__________________________________________________________________________
467bool ReadPedeHists::readNext(std::ifstream &file, TH1 *&hist,
468 std::pair<TGraph*,Option_t*> &graphOpt)
469{
470 hist = 0;
471 graphOpt.first = 0;
472 graphOpt.second = 0;
473
474 TString type;
475 while (true) {
476 if (file.eof()) break;
477 file >> type;
478 if (file.fail() || (type != "Histogram" && type != "XY-Data")) {
479 TString line;
480 line.ReadLine(file);
481 if (line != "" && line.Length() != line.CountChar(' ')) { // not empty
482 ::Error("ReadPedeHists::readNext",
483 "Expect 'Histogram' or 'XY-Data', but failed, line is '%s'",
484 line.Data());
485 if (proceedTo(file, "end of")) line.ReadLine(file); // just skip rest of line...
486 }
487 }
488
489 if (type == "Histogram") hist = readNextHist(file);
490 if (type == "XY-Data") graphOpt = readNextGraph(file);
491 if (hist || graphOpt.first) break;
492 }
493
494 return (hist || graphOpt.first);
495}
496
497//__________________________________________________________________________
498void ReadPedeHists::read(std::ifstream &file)
499{
500 theHists.clear();
501 theGraphOpts.clear();
502
503 TH1 *hist = 0;
504 std::pair<TGraph*, Option_t*> graphOpt(0,0); // graph and its drawing option
505 while (readNext(file, hist, graphOpt)) {
506 if (hist) theHists.push_back(hist);
507 if (graphOpt.first) theGraphOpts.push_back(graphOpt);
508 }
509}
510
511//__________________________________________________________________________
513{
514 theCanvases.clear(); // memory leak?
515
516 const Int_t nHistX = 3;
517 const Int_t nPixelX = 700;
518 const Int_t nHistY = 2;
519 const Int_t nPixelY = 500;
520 Int_t last = nHistX * nHistY;
521 unsigned int iH = 0;
522
523 while (iH < theHists.size()) {
524 if (last >= nHistX * nHistY) {
525 unsigned int canCorner = theCanvases.size() * 20;
526 theCanvases.push_back(new TCanvas(Form("hists%d", iH), "",
527 canCorner, canCorner, nPixelX, nPixelY));
528 theCanvases.back()->Divide(nHistX, nHistY);
529 last = 0;
530 }
531 theCanvases.back()->cd(++last);
532 theHists[iH]->Draw();
533 ++iH;
534 }
535
536 last = nHistX * nHistY;
537 iH = 0;
538 while (iH < theGraphOpts.size()) {
539 if (last >= nHistX * nHistY) {
540 unsigned int canCorner = theCanvases.size() * 20;
541 theCanvases.push_back(new TCanvas(Form("graphs%d", iH), "",
542 canCorner, canCorner, nPixelX, nPixelY));
543 theCanvases.back()->Divide(nHistX, nHistY);
544 last = 0;
545 }
546 theCanvases.back()->cd(++last);
547 theGraphOpts[iH].first->Draw(theGraphOpts[iH].second);
548 ++iH;
549 }
550}
551
552//__________________________________________________________________________
553void ReadPedeHists::Print(const char *printFileName)
554{
555 std::vector<TCanvas*>::iterator iC = theCanvases.begin(), iE = theCanvases.end();
556 if (iC == iE) return; // empty...
557
558 theCanvases.front()->Print(Form("%s[", printFileName)); // just open ps
559 while(iC != iE) {
560 (*iC)->Print(printFileName);
561 ++iC;
562 }
563 theCanvases.front()->Print(Form("%s]", printFileName)); // just close ps
564
565}
566
567//__________________________________________________________________________
568void ReadPedeHists::Write(const char *rootFileName)
569{
570 if (theHists.empty() && theGraphOpts.empty()) return;
571
572 ::Info("ReadPedeHists::Write", "(Re-)Creating ROOT file %s.", rootFileName);
573
574 TDirectory *oldDir = gDirectory;
575 TFile *rootFile = TFile::Open(rootFileName, "RECREATE");
576
577 for (std::vector<TH1*>::iterator iH = theHists.begin(), iE = theHists.end();
578 iH != iE; ++iH) {
579 (*iH)->Write();
580 }
581
582 for (std::vector<std::pair<TGraph*,Option_t*> >::iterator iG = theGraphOpts.begin(),
583 iE = theGraphOpts.end(); iG != iE; ++iG) {
584 (*iG).first->Write();
585 }
586
587 delete rootFile;
588 oldDir->cd();
589}
590
591//__________________________________________________________________________
592//__________________________________________________________________________
593//__________________________________________________________________________
594void readPedeHists(Option_t *option, const char *txtFile)
595{
596 ReadPedeHists reader(txtFile);
597 TString opt(option);
598 opt.ToLower();
599
600 const bool oldBatch = gROOT->IsBatch();
601 if (opt.Contains("nodraw")) {
602 opt.ReplaceAll("nodraw", "");
603 gROOT->SetBatch(true);
604 }
605
606 reader.Draw();
607
608 if (opt.Contains("print")) {
609 opt.ReplaceAll("print", "");
610 reader.Print(TString(Form("%s.ps", txtFile)));
611 }
612
613 if (opt.Contains("write")) {
614 opt.ReplaceAll("write", "");
615 reader.Write(TString(Form("%s.root", txtFile)));
616 }
617
618 gROOT->SetBatch(oldBatch);
619 opt.ReplaceAll(" ", "");
620 if (!opt.IsNull()) {
621 ::Warning("readPedeHists", "Unknown option '%s', know 'nodraw', 'print' and 'write'.",
622 opt.Data());
623 }
624}
void readNumVersTypeTitle(std::ifstream &file, Int_t &num, Int_t &version, Int_t &type, TString &title)
std::vector< TCanvas * > theCanvases
bool readBasicVector(std::ifstream &aStream, std::vector< T > &outVector)
bool readNext(std::ifstream &file, TH1 *&hist, std::pair< TGraph *, Option_t * > &graphOpt)
std::vector< TH1 * > theHists
TH1 * readNextHist(std::ifstream &file)
void read(std::ifstream &stream)
std::vector< std::pair< TGraph *, Option_t * > > theGraphOpts
void Write(const char *rootFileName)
ReadPedeHists(const char *txtFile)
bool proceedTo(std::ifstream &aStream, const TString &pattern)
void Print(const char *printFileName)
std::ifstream theStream
std::pair< TGraph *, Option_t * > readNextGraph(std::ifstream &file)
bool readBasic(std::ifstream &aStream, T &outValue)
real(mps), dimension(nplan) sigma
measurement sigma (hit)
Definition: mptest1.f90:65
void readPedeHists(Option_t *option="", const char *txtFile="millepede.his")