/[PAMELA software]/PamCut/CollectionActions/Histo2DActions/Histo2DAction/Histo2DAction.h
ViewVC logotype

Contents of /PamCut/CollectionActions/Histo2DActions/Histo2DAction/Histo2DAction.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (show annotations) (download)
Tue Dec 8 17:38:39 2009 UTC (15 years ago) by pam-fi
Branch: MAIN
Changes since 1.4: +15 -15 lines
File MIME type: text/plain
Accessor for ROOT histogram added. Creation of default axis added. Some bugs fixed.

1 /*
2 * Histo2DAction.h
3 *
4 * Created on: 01/set/2009
5 * Author: Nicola Mori
6 */
7
8 /*! @file Histo2DAction.h The Histo2DAction class declaration file. */
9
10 #ifndef HISTO2DACTION_H_
11 #define HISTO2DACTION_H_
12
13 #include "../../CollectionAction/CollectionAction.h"
14
15 #include <TH2I.h>
16 #include <TH2F.h>
17 #include <TH2D.h>
18
19 /*! @brief An abstract action that fills a 2-dimensional histogram.
20 *
21 * This abstract class provides all the common methods needed to produce 2-dimensional histograms,
22 * like reading the bins from a file, setting up the root histogram and so on. It handles histograms
23 * both as a ROOT histogram and as a matrix. The OnGood() method has to be implemented in children
24 * classes, with the specific physical variables to fill the histogram. Eg., a dE/dx vs. rigidity
25 * histogram class derived from Histo2DAction could have this OnGood() implementation:
26 *
27 * Fill(event->GetTrack(0)->GetTrkTrack()->GetRigidity(), event->GetTrack(0)->GetTrkTrack()->GetDEDX());
28 *
29 * Fill() automatically fills the ROOT and the matrix histogram.
30 * The template argument HistoType is the type of histogram (Int_t, Double_t,...). In current
31 * implementation, ROOT histograms are only supported for Int_t, Float_t and Double_t. For all the
32 * other types, only the the matrix histogram with text file output will be produced. If you want to
33 * implement some new types of ROOT histograms, please specialize the template definitions of _CreateHisto()
34 * and ~Histo2DAction() (see the specializations in the .cpp file).
35 *
36 * */
37 template<class HistoType>
38 class Histo2DAction: public CollectionAction {
39 public:
40
41 /*! @brief Constructor.
42 *
43 * The histogram binning is defined by the content of the bins vector.
44 *
45 * @param actionName The action's name.
46 * @param title The ROOT histogram title.
47 * @param outFileBase The file base name for the histogram output. If "", no output file will be produced.
48 * @param mode The mode of ROOT file creation (see documentation of TFile constructor
49 * in ROOT's reference guide).
50 * @param outRoot If true, an output ROOT file named outFileBase + ".root" will be produced.
51 * @param outText If true, an output text file named outFileBase + ".txt" will be produced. It will overwrite an
52 * eventually existing file with the same name.
53 */
54 Histo2DAction(const char *actionName, TString title = "", TString outFileBase = "", TString mode = "UPDATE",
55 bool outRoot = false, bool outText = false);
56
57 /*! @brief Destructor
58 *
59 * This has to be specialized for each kind of ROOT histogram.
60 *
61 * */
62 ~Histo2DAction();
63
64 /*! @brief Sets the X axis using a binning read from a a vector.
65 *
66 * @param label The axis' label.
67 * @param bins A vector containing the histogram binning (with sign, eventually; from lowest to highest).
68 */
69 void SetXAxis(TString label, vector<float> &bins);
70 /*! @brief Sets the X axis reading the binning from a file.
71 *
72 * @param label The axis' label.
73 * @param binsFile The file containing the bins (with sign, eventually; from lowest to highest).
74 */
75 void SetXAxis(TString label, TString binsFile);
76 /*! @brief Sets the X axis specifying the binning parameters.
77 *
78 * @param label The axis' label.
79 * @param nBins The number of bins.
80 * @param min The lower limit of the axis.
81 * @param max The upper limit of the axis.
82 * @param logBinning If true, bins will be logarithmically spaced.
83 */
84 void SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning = false);
85
86 /*! @brief Sets the Y axis using a binning read from a a vector.
87 *
88 * @param label The axis' label.
89 * @param bins A vector containing the histogram binning (with sign, eventually; from lowest to highest).
90 */
91 void SetYAxis(TString label, vector<float> &bins);
92 /*! @brief Sets the Y axis reading the binning from a file.
93 *
94 * @param label The axis' label.
95 * @param binsFile The file containing the bins (with sign, eventually; from lowest to highest).
96 */
97 void SetYAxis(TString label, TString binsFile);
98
99 /*! @brief Sets the Y axis specifying the binning parameters.
100 *
101 * @param label The axis' label.
102 * @param nBins The number of bins.
103 * @param min The lower limit of the axis.
104 * @param max The upper limit of the axis.
105 * @param logBinning If true, bins will be logarithmically spaced.
106 */
107 void SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning = false);
108
109 /*! @brief Sets the ROOT histogram's title.
110 *
111 * @param title The histogram title as it will appear on the histogram itself.
112 */
113 void SetTitle(TString &title) {
114 _title = title;
115 }
116 /*! @brief Sets up the histogram
117 *
118 * This routine effectively prepares the histogram, after the desired parameters has been set by #SetXAxis() and #SetYAxis().
119 *
120 * @param events Pointer to PamLevel2 events (unused).
121 */
122 void Setup(PamLevel2 *events) {
123 CollectionAction::Setup(events);
124 _InitHistos();
125
126 }
127
128 /*! @ brief Post-analysis tasks.
129 *
130 * This method writes the output on a ROOT file if the proper option was set in the constructor.
131 */
132 void Finalize();
133
134 /*! @brief Returns the histogram.
135 *
136 * Element [i][j] is the value of the i-th bin on the Y axis and on the j-th bin on the X axis. This sort
137 * of index inversion is due to logically maintain the Y values on the vertical axis of the matrix, which
138 * is indexed by the row number (which conventionally is the first index of a matrix element).
139 *
140 * @return A reference to a matrix containing the values of the bins of the histogram.
141 */
142 SimpleMatrix<HistoType> &GetHisto() {
143 return _histo;
144 }
145
146 /*! @brief Returns a pointer to the ROOT histogram.
147 *
148 * @return A pointer to the root histogram
149 */
150 TH2 *GetRootHisto() {
151 return _rootHisto;
152 }
153
154 /*! Fills the ROOT and the vector histogram.
155 *
156 * @param xValue The value of the X coordinate associated to the event.
157 * @param yValue The value of the Y coordinate associated to the event.
158 * @param weight The weight which will be applied to the event.
159 */
160 void Fill(double xValue, double yValue, double weight = 1.);
161
162 /*! @brief Gets the X overflow histogram.
163 *
164 * The returned vector contains the histogram filled with the X overflow events, eg.,
165 * those events whose X variable falls above the upper limit of the X axis. This histogram
166 * is binned like the Y axis.
167 *
168 * @return The X overflow histogram.
169 */
170 vector<HistoType>& GetXOverflow() {
171 return _xOverflow;
172 }
173
174 /*! @brief Gets the X underflow histogram.
175 *
176 * The returned vector contains the histogram filled with the X underflow events, eg.,
177 * those events whose X variable falls below the lower limit of the X axis. This histogram
178 * is binned like the Y axis.
179 *
180 * @return The X underflow histogram.
181 */
182 vector<HistoType>& GetXUnderflow() {
183 return _xUnderflow;
184 }
185
186 /*! @brief Gets the Y overflow histogram.
187 *
188 * The returned vector contains the histogram filled with the Y overflow events, eg.,
189 * those events whose Y variable falls above the upper limit of the Y axis. This histogram
190 * is binned like the X axis.
191 *
192 * @return The Y overflow histogram.
193 */
194 vector<HistoType>& GetYOverflow() {
195 return _yOverflow;
196 }
197
198 /*! @brief Gets the Y underflow histogram.
199 *
200 * The returned vector contains the histogram filled with the Y underflow events, eg.,
201 * those events whose Y variable falls below the lower limit of the Y axis. This histogram
202 * is binned like the X axis.
203 *
204 * @return The Y underflow histogram.
205 */
206 vector<HistoType>& GetYUnderflow() {
207 return _yUnderflow;
208 }
209
210 /*! @brief Gets the X-Y overflow histogram.
211 *
212 * This single-bin histogram is built with events which fall above the lower limits of both the
213 * X and Y axis.
214 *
215 * @return The X-Y overflow histogram.
216 */
217 HistoType GetXOverYOverflow() {
218 return _xOverYOverflow;
219 }
220
221 /*! @brief Gets the X-Y underflow histogram.
222 *
223 * This single-bin histogram is built with events fall below the lower limits of both the X and Y axis.
224 *
225 * @return The X-Y overflow histogram.
226 */
227 HistoType GetXUnderYUnderflow() {
228 return _xUnderYUnderflow;
229 }
230 /*! @brief Gets the (X over - Y under)flow histogram.
231 *
232 * This single-bin histogram is built with events fall above the lower limit of the Xaxis and below
233 * that of the Y axis.
234 *
235 * @return The (X over -Y under)flow histogram.
236 */
237 HistoType GetXOverYUnderflow() {
238 return _xOverYUnderflow;
239 }
240
241 /*! @brief Gets the (X under - Y over)flow histogram.
242 *
243 * This single-bin histogram is built with events fall below the lower limit of the X axis and above
244 * that of the Y axis.
245 *
246 * @return The (X under -Y over)flow histogram.
247 */
248 HistoType GetXUnderYOverflow() {
249 return _xUnderYOverflow;
250 }
251
252 protected:
253
254 /*! @brief The vector containing the limits of the X bins(from lower to higher). */
255 std::vector<float> _xBins;
256 /*! @brief The vector containing the limits of the Y bins(from lower to higher). */
257 std::vector<float> _yBins;
258 /*! @brief A matrix containing the value of the histogram for each X-Y bin. */
259 SimpleMatrix<HistoType> _histo;
260 /*! @brief The ROOT histogram. */
261 TH2 *_rootHisto;
262
263 private:
264
265 vector<HistoType> _xUnderflow, _xOverflow, _yUnderflow, _yOverflow;
266 HistoType _xUnderYUnderflow, _xOverYOverflow, _xUnderYOverflow, _xOverYUnderflow;
267 TString _outFileBase;
268 TString _mode;
269 TString _title, _xLabel, _yLabel;
270 bool _outRoot;
271 bool _outText;
272
273 void _CreateHisto();
274 void _InitHistos();
275 };
276
277 template<class HistoType>
278 void Histo2DAction<HistoType>::_CreateHisto() {
279
280 // No ROOT histogram for generic type; see template specializations in .cpp file.
281 _rootHisto = NULL;
282 }
283
284 // Specializations for _CreateHistos(). See Histo2DAction.cpp
285 template<>
286 void Histo2DAction<Int_t>::_CreateHisto();
287
288 template<>
289 void Histo2DAction<Float_t>::_CreateHisto();
290
291 template<>
292 void Histo2DAction<Double_t>::_CreateHisto();
293
294 template<class HistoType>
295 void Histo2DAction<HistoType>::_InitHistos() {
296
297 _CreateHisto();
298 if (_xBins.size() < 2) // SetXAxis not called by the main program, or wrongly filled (only 1 bin limit)
299 SetXAxis("Default X", 10, 0., 1.);
300 if (_yBins.size() < 2) // SetYAxis not called by the main program, or wrongly filled (only 1 bin limit)
301 SetYAxis("Default Y", 10, 0., 1.);
302
303 if (_rootHisto) {
304 Double_t *auxXArray = new Double_t[_xBins.size()];
305
306 for (unsigned int i = 0; i < _xBins.size(); i++) {
307 auxXArray[i] = _xBins[i];
308 }
309
310 Double_t *auxYArray = new Double_t[_yBins.size()];
311
312 for (unsigned int i = 0; i < _yBins.size(); i++) {
313 auxYArray[i] = _yBins[i];
314 }
315
316 _rootHisto->SetBins(_xBins.size() - 1, auxXArray, _yBins.size() - 1, auxYArray);
317 _rootHisto->SetName(GetName());
318 _rootHisto->SetTitle(_title);
319 _rootHisto->SetXTitle(_xLabel);
320 _rootHisto->SetYTitle(_yLabel);
321
322 delete[] auxXArray;
323 delete[] auxYArray;
324 }
325
326 /* The row index (first) corresponds to the position on the vertical (Y) axis. */
327 _histo.Resize(_yBins.size() - 1, _xBins.size() - 1);
328
329 _xUnderflow.resize(_yBins.size());
330 _xOverflow.resize(_yBins.size());
331 _yUnderflow.resize(_xBins.size());
332 _yOverflow.resize(_xBins.size());
333
334 }
335
336 template<class HistoType>
337 Histo2DAction<HistoType>::~Histo2DAction() {
338
339 delete _rootHisto;
340 _rootHisto = NULL;
341 }
342
343 template<class HistoType>
344 Histo2DAction<HistoType>::Histo2DAction(const char *actionName, TString title, TString outFileBase, TString mode,
345 bool outRoot, bool outText) :
346 CollectionAction(actionName), _xBins(0), _yBins(0), _histo(0, 0), _rootHisto(NULL), _xUnderflow(0), _xOverflow(0),
347 _yUnderflow(0), _yOverflow(0), _outFileBase(outFileBase), _mode(mode), _title(title), _xLabel(""), _yLabel(""),
348 _outRoot(outRoot), _outText(outText) {
349
350 }
351
352 template<class HistoType>
353 void Histo2DAction<HistoType>::SetXAxis(TString label, vector<float> &bins) {
354
355 _xBins = bins;
356 _xLabel = label;
357
358 }
359
360 template<class HistoType>
361 void Histo2DAction<HistoType>::SetYAxis(TString label, vector<float> &bins) {
362
363 _yBins = bins;
364 _yLabel = label;
365
366 }
367
368 template<class HistoType>
369 void Histo2DAction<HistoType>::SetXAxis(TString label, TString binsFile) {
370
371 // Reading the bins from file
372 ifstream binListFile;
373 binListFile.open(binsFile);
374
375 TString auxString;
376 _xBins.resize(0);
377 while (!binListFile.eof()) {
378 binListFile >> auxString;
379 if (auxString != "") {
380 _xBins.push_back(auxString.Atof());
381 }
382 }
383 binListFile.close();
384
385 _xLabel = label;
386
387 }
388
389 template<class HistoType>
390 void Histo2DAction<HistoType>::SetYAxis(TString label, TString binsFile) {
391
392 // Reading the bins from file
393 ifstream binListFile;
394 binListFile.open(binsFile);
395
396 TString auxString;
397 _yBins.resize(0);
398 while (!binListFile.eof()) {
399 binListFile >> auxString;
400 if (auxString != "") {
401 _yBins.push_back(auxString.Atof());
402 }
403 }
404 binListFile.close();
405
406 _yLabel = label;
407
408 }
409
410 template<class HistoType>
411 void Histo2DAction<HistoType>::SetXAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {
412
413 _xBins.resize(nBins + 1);
414
415 if (!logBinning || (logBinning && min <= 0.)) {
416
417 #ifdef DEBUGPAMCUT
418 if (logbinning && rigMin <= 0.)
419 cout << "Warning: logarithmic binning was chosen for X axis but min <= 0. Using linear binning."
420 #endif
421
422 float step = (max - min) / nBins;
423 for (unsigned int i = 0; i < nBins + 1; i++) {
424 _xBins[i] = min + i * step;
425 }
426
427 }
428 else {
429
430 double maxExp = log10(max / min);
431 for (unsigned int i = 0; i < nBins + 1; i++) {
432 _xBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
433 }
434
435 }
436
437 _xLabel = label;
438 }
439
440 template<class HistoType>
441 void Histo2DAction<HistoType>::SetYAxis(TString label, unsigned int nBins, float min, float max, bool logBinning) {
442
443 _yBins.resize(nBins + 1);
444
445 if (!logBinning || (logBinning && min <= 0.)) {
446
447 #ifdef DEBUGPAMCUT
448 if (logbinning && rigMin <= 0.)
449 cout << "Warning: logarithmic binning was chosen for Y axis but min <= 0. Using linear binning."
450 #endif
451
452 float step = (max - min) / nBins;
453 for (unsigned int i = 0; i < nBins + 1; i++) {
454 _yBins[i] = min + i * step;
455 }
456
457 }
458 else {
459
460 double maxExp = log10(max / min);
461 for (unsigned int i = 0; i < nBins + 1; i++) {
462 _yBins[i] = min * pow(10., (double) i / ((double) nBins) * maxExp);
463 }
464
465 }
466
467 _yLabel = label;
468 }
469
470 template<class HistoType>
471 void Histo2DAction<HistoType>::Fill(double xValue, double yValue, double weight) {
472
473 _rootHisto->Fill(xValue, yValue, weight);
474
475 int xBin = 0, yBin = 0;
476
477 bool UOflow = false;
478
479 if (xValue < _xBins.front()) {
480 xBin = -1;
481 UOflow = true;
482 }
483
484 if (xValue >= _xBins.back()) {
485 xBin = _xBins.size();
486 UOflow = true;
487 }
488
489 if (yValue < _yBins.front()) {
490 yBin = -1;
491 UOflow = true;
492 }
493
494 if (yValue >= _yBins.back()) {
495 yBin = _yBins.size();
496 UOflow = true;
497 }
498
499 if (xBin == 0) {
500
501 while (xValue >= _xBins[xBin]) {
502 xBin++;
503 }
504 xBin--;
505 }
506
507 if (yBin == 0) {
508
509 while (yValue >= _yBins[yBin]) {
510 yBin++;
511 }
512 yBin--;
513 }
514
515 if (!UOflow) {
516 // No under or over flow.
517 _histo[yBin][xBin] += (HistoType) weight;
518 return;
519 }
520 else {
521 // An under or over flow has happened.
522 if (xBin == -1) {
523 if (yBin == -1) {
524 _xUnderYUnderflow += (HistoType) weight;
525 return;
526 }
527 else {
528 if (yBin == (int) _yBins.size()) {
529 _xUnderYOverflow += (HistoType) weight;
530 return;
531 }
532 else {
533 _xUnderflow[yBin] += (HistoType) weight;
534 return;
535 }
536 }
537 }
538
539 if (xBin == (int) _xBins.size()) {
540 if (yBin == -1) {
541 _xOverYUnderflow += (HistoType) weight;
542 return;
543 }
544 else {
545 if (yBin == (int) _yBins.size()) {
546 _xOverYOverflow += (HistoType) weight;
547 return;
548 }
549 else {
550 _xOverflow[yBin] += (HistoType) weight;
551 return;
552 }
553 }
554 }
555
556 if (yBin == -1) {
557 _yUnderflow[xBin] += (HistoType) weight;
558 return;
559 }
560
561 if (yBin == (int) _yBins.size()) {
562 _yOverflow[xBin] += (HistoType) weight;
563 return;
564 }
565
566 }
567 }
568
569 template<class HistoType>
570 void Histo2DAction<HistoType>::Finalize() {
571
572 if (_outFileBase != "") {
573 // Write the ROOT file
574 if (_rootHisto && _outRoot) {
575 TFile outRootFile((_outFileBase + ".root"), _mode);
576 outRootFile.cd();
577 _rootHisto->Write();
578 outRootFile.Close();
579 }
580
581 //Write the text file
582 if (_outText) {
583 ofstream outTextFile((_outFileBase + ".txt").Data(), ios_base::out);
584 for (unsigned int i = 0; i < _histo.GetNRows(); i++) {
585 for (unsigned int j = 0; j < _histo.GetNCols(); j++) {
586 outTextFile << _histo[i][j] << " ";
587 }
588 outTextFile << "\n";
589 }
590 outTextFile << endl;
591 outTextFile.close();
592 }
593 }
594
595 }
596 #endif /* HISTO2DACTION_H_ */

  ViewVC Help
Powered by ViewVC 1.1.23