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_ */ |