| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms. 2 // See the file tools.license for terms.
3 3
4 #ifndef tools_histo_h2 4 #ifndef tools_histo_h2
5 #define tools_histo_h2 5 #define tools_histo_h2
6 6
7 #include "b2" 7 #include "b2"
8 8
9 namespace tools { 9 namespace tools {
10 namespace histo { 10 namespace histo {
11 11
12 template <class TC,class TO,class TN,class TW, 12 template <class TC,class TO,class TN,class TW,class TH>
13 class h2 : public b2<TC,TO,TN,TW,TH> { 13 class h2 : public b2<TC,TO,TN,TW,TH> {
14 typedef b2<TC,TO,TN,TW,TH> parent; 14 typedef b2<TC,TO,TN,TW,TH> parent;
15 public: 15 public:
16 typedef histo_data<TC,TO,TN,TW> hd_t; 16 typedef histo_data<TC,TO,TN,TW> hd_t;
17 typedef typename b2<TC,TO,TN,TW,TH>::bn_t bn 17 typedef typename b2<TC,TO,TN,TW,TH>::bn_t bn_t;
18 protected: 18 protected:
19 virtual TH get_bin_height(TO a_offset) const 19 virtual TH get_bin_height(TO a_offset) const { //TH should be the same as TW
20 return parent::m_bin_Sw[a_offset]; 20 return parent::m_bin_Sw[a_offset];
21 } 21 }
22 public: 22 public:
23 23
24 virtual TH bin_error(int aI,int aJ) const { 24 virtual TH bin_error(int aI,int aJ) const {
25 TO offset; 25 TO offset;
26 if(!parent::_find_offset(aI,aJ,offset)) re 26 if(!parent::_find_offset(aI,aJ,offset)) return 0;
27 return ::sqrt(parent::m_bin_Sw2[offset]); 27 return ::sqrt(parent::m_bin_Sw2[offset]);
28 } 28 }
29 29
30 public: 30 public:
31 bool multiply(TW a_factor){return parent::ba 31 bool multiply(TW a_factor){return parent::base_multiply(a_factor);}
32 bool scale(TW a_factor) {return multiply(a_f 32 bool scale(TW a_factor) {return multiply(a_factor);}
33 33
34 void copy_from_data(const hd_t& a_from) {par 34 void copy_from_data(const hd_t& a_from) {parent::base_from_data(a_from);}
35 hd_t get_histo_data() const {return *this;} 35 hd_t get_histo_data() const {return *this;} //deprecated. Keep it for g4tools.
36 36
37 bool reset() { 37 bool reset() {
38 parent::base_reset(); 38 parent::base_reset();
39 return true; 39 return true;
40 } 40 }
41 41
42 bool fill(TC aX,TC aY,TW aWeight = 1) { 42 bool fill(TC aX,TC aY,TW aWeight = 1) {
43 if(parent::m_dimension!=2) return false; 43 if(parent::m_dimension!=2) return false;
44 44
45 bn_t ibin,jbin; 45 bn_t ibin,jbin;
46 if(!parent::m_axes[0].coord_to_absolute_in 46 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
47 if(!parent::m_axes[1].coord_to_absolute_in 47 if(!parent::m_axes[1].coord_to_absolute_index(aY,jbin)) return false;
48 TO offset = ibin + jbin * parent::m_axes[1 48 TO offset = ibin + jbin * parent::m_axes[1].m_offset;
49 49
50 parent::m_bin_entries[offset]++; 50 parent::m_bin_entries[offset]++;
51 parent::m_bin_Sw[offset] += aWeight; 51 parent::m_bin_Sw[offset] += aWeight;
52 parent::m_bin_Sw2[offset] += aWeight * aWe 52 parent::m_bin_Sw2[offset] += aWeight * aWeight;
53 53
54 TC xw = aX * aWeight; 54 TC xw = aX * aWeight;
55 TC x2w = aX * xw; 55 TC x2w = aX * xw;
56 parent::m_bin_Sxw[offset][0] += xw; 56 parent::m_bin_Sxw[offset][0] += xw;
57 parent::m_bin_Sx2w[offset][0] += x2w; 57 parent::m_bin_Sx2w[offset][0] += x2w;
58 58
59 TC yw = aY * aWeight; 59 TC yw = aY * aWeight;
60 TC y2w = aY * yw; 60 TC y2w = aY * yw;
61 parent::m_bin_Sxw[offset][1] += yw; 61 parent::m_bin_Sxw[offset][1] += yw;
62 parent::m_bin_Sx2w[offset][1] += y2w; 62 parent::m_bin_Sx2w[offset][1] += y2w;
63 63
64 bool inRange = true; 64 bool inRange = true;
65 if(ibin==0) inRange = false; 65 if(ibin==0) inRange = false;
66 else if(ibin==(parent::m_axes[0].m_number_ 66 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
67 67
68 if(jbin==0) inRange = false; 68 if(jbin==0) inRange = false;
69 else if(jbin==(parent::m_axes[1].m_number_ 69 else if(jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
70 70
71 parent::m_all_entries++; 71 parent::m_all_entries++;
72 if(inRange) { 72 if(inRange) {
73 parent::m_in_range_plane_Sxyw[0] += aX * 73 parent::m_in_range_plane_Sxyw[0] += aX * aY * aWeight;
74 74
75 // fast getters : 75 // fast getters :
76 parent::m_in_range_entries++; 76 parent::m_in_range_entries++;
77 parent::m_in_range_Sw += aWeight; 77 parent::m_in_range_Sw += aWeight;
78 parent::m_in_range_Sw2 += aWeight*aWeigh 78 parent::m_in_range_Sw2 += aWeight*aWeight;
79 79
80 parent::m_in_range_Sxw[0] += xw; 80 parent::m_in_range_Sxw[0] += xw;
81 parent::m_in_range_Sx2w[0] += x2w; 81 parent::m_in_range_Sx2w[0] += x2w;
82 82
83 parent::m_in_range_Sxw[1] += yw; 83 parent::m_in_range_Sxw[1] += yw;
84 parent::m_in_range_Sx2w[1] += y2w; 84 parent::m_in_range_Sx2w[1] += y2w;
85 } 85 }
86 86
87 return true; 87 return true;
88 } 88 }
89 89
90 bool set_bin_content(bn_t a_ibin,bn_t a_jbin 90 bool set_bin_content(bn_t a_ibin,bn_t a_jbin,
91 TN a_entries,TW a_Sw,TW 91 TN a_entries,TW a_Sw,TW a_Sw2,
92 TC a_Sxw,TC a_Sx2w,TC a 92 TC a_Sxw,TC a_Sx2w,TC a_Syw,TC a_Sy2w) {
93 if(parent::m_dimension!=2) return false; 93 if(parent::m_dimension!=2) return false;
94 if(a_ibin>(parent::m_axes[0].m_number_of_b 94 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) return false;
95 if(a_jbin>(parent::m_axes[1].m_number_of_b 95 if(a_jbin>(parent::m_axes[1].m_number_of_bins+1)) return false;
96 96
97 bool inRange = true; 97 bool inRange = true;
98 if(a_ibin==0) inRange = false; 98 if(a_ibin==0) inRange = false;
99 else if(a_ibin==(parent::m_axes[0].m_numbe 99 else if(a_ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
100 if(a_jbin==0) inRange = false; 100 if(a_jbin==0) inRange = false;
101 else if(a_jbin==(parent::m_axes[1].m_numbe 101 else if(a_jbin==(parent::m_axes[1].m_number_of_bins+1)) inRange = false;
102 102
103 TO offset = a_ibin + a_jbin * parent::m_ax 103 TO offset = a_ibin + a_jbin * parent::m_axes[1].m_offset;
104 104
105 parent::m_all_entries -= parent::m_bin_ent 105 parent::m_all_entries -= parent::m_bin_entries[offset];
106 if(inRange) { 106 if(inRange) {
107 parent::m_in_range_entries -= parent::m_ 107 parent::m_in_range_entries -= parent::m_bin_entries[offset];
108 parent::m_in_range_Sw -= parent::m_bin_S 108 parent::m_in_range_Sw -= parent::m_bin_Sw[offset];
109 parent::m_in_range_Sw2 -= parent::m_bin_ 109 parent::m_in_range_Sw2 -= parent::m_bin_Sw2[offset];
110 parent::m_in_range_Sxw[0] -= parent::m_b 110 parent::m_in_range_Sxw[0] -= parent::m_bin_Sxw[offset][0];
111 parent::m_in_range_Sx2w[0] -= parent::m_ 111 parent::m_in_range_Sx2w[0] -= parent::m_bin_Sx2w[offset][0];
112 parent::m_in_range_Sxw[1] -= parent::m_b 112 parent::m_in_range_Sxw[1] -= parent::m_bin_Sxw[offset][1];
113 parent::m_in_range_Sx2w[1] -= parent::m_ 113 parent::m_in_range_Sx2w[1] -= parent::m_bin_Sx2w[offset][1];
114 } 114 }
115 115
116 parent::m_bin_entries[offset] = a_entries; 116 parent::m_bin_entries[offset] = a_entries;
117 parent::m_bin_Sw[offset] = a_Sw; 117 parent::m_bin_Sw[offset] = a_Sw;
118 parent::m_bin_Sw2[offset] = a_Sw2; 118 parent::m_bin_Sw2[offset] = a_Sw2;
119 119
120 parent::m_bin_Sxw[offset][0] = a_Sxw; 120 parent::m_bin_Sxw[offset][0] = a_Sxw;
121 parent::m_bin_Sx2w[offset][0] = a_Sx2w; 121 parent::m_bin_Sx2w[offset][0] = a_Sx2w;
122 parent::m_bin_Sxw[offset][1] = a_Syw; 122 parent::m_bin_Sxw[offset][1] = a_Syw;
123 parent::m_bin_Sx2w[offset][1] = a_Sy2w; 123 parent::m_bin_Sx2w[offset][1] = a_Sy2w;
124 124
125 parent::m_all_entries += a_entries; 125 parent::m_all_entries += a_entries;
126 if(inRange) { 126 if(inRange) {
127 //parent::m_in_range_plane_Sxyw[0] ??? i 127 //parent::m_in_range_plane_Sxyw[0] ??? ill-defined.
128 128
129 parent::m_in_range_entries += a_entries; 129 parent::m_in_range_entries += a_entries;
130 parent::m_in_range_Sw += a_Sw; 130 parent::m_in_range_Sw += a_Sw;
131 parent::m_in_range_Sw2 += a_Sw2; 131 parent::m_in_range_Sw2 += a_Sw2;
132 132
133 parent::m_in_range_Sxw[0] += a_Sxw; 133 parent::m_in_range_Sxw[0] += a_Sxw;
134 parent::m_in_range_Sx2w[0] += a_Sx2w; 134 parent::m_in_range_Sx2w[0] += a_Sx2w;
135 parent::m_in_range_Sxw[1] += a_Syw; 135 parent::m_in_range_Sxw[1] += a_Syw;
136 parent::m_in_range_Sx2w[1] += a_Sy2w; 136 parent::m_in_range_Sx2w[1] += a_Sy2w;
137 } 137 }
138 138
139 return true; 139 return true;
140 } 140 }
141 141
142 bool get_bin_content(bn_t a_ibin,bn_t a_jbin 142 bool get_bin_content(bn_t a_ibin,bn_t a_jbin,
143 TN& a_entries,TW& a_Sw, 143 TN& a_entries,TW& a_Sw,TW& a_Sw2,
144 TC& a_Sxw,TC& a_Sx2w, 144 TC& a_Sxw,TC& a_Sx2w,
145 TC& a_Syw,TC& a_Sy2w) { 145 TC& a_Syw,TC& a_Sy2w) {
146 if(parent::m_dimension!=2) { 146 if(parent::m_dimension!=2) {
147 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 147 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
148 a_Sxw = 0;a_Sx2w = 0; 148 a_Sxw = 0;a_Sx2w = 0;
149 a_Syw = 0;a_Sy2w = 0; 149 a_Syw = 0;a_Sy2w = 0;
150 return false; 150 return false;
151 } 151 }
152 if(a_ibin>(parent::m_axes[0].m_number_of_b 152 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) {
153 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 153 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
154 a_Sxw = 0;a_Sx2w = 0; 154 a_Sxw = 0;a_Sx2w = 0;
155 a_Syw = 0;a_Sy2w = 0; 155 a_Syw = 0;a_Sy2w = 0;
156 return false; 156 return false;
157 } 157 }
158 if(a_jbin>(parent::m_axes[1].m_number_of_b 158 if(a_jbin>(parent::m_axes[1].m_number_of_bins+1)) {
159 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 159 a_entries = 0;a_Sw = 0;a_Sw2 = 0;
160 a_Sxw = 0;a_Sx2w = 0; 160 a_Sxw = 0;a_Sx2w = 0;
161 a_Syw = 0;a_Sy2w = 0; 161 a_Syw = 0;a_Sy2w = 0;
162 return false; 162 return false;
163 } 163 }
164 164
165 TO offset = a_ibin + a_jbin * parent::m_ax 165 TO offset = a_ibin + a_jbin * parent::m_axes[1].m_offset;
166 166
167 a_entries = parent::m_bin_entries[offset]; 167 a_entries = parent::m_bin_entries[offset];
168 a_Sw = parent::m_bin_Sw[offset]; 168 a_Sw = parent::m_bin_Sw[offset];
169 a_Sw2 = parent::m_bin_Sw2[offset]; 169 a_Sw2 = parent::m_bin_Sw2[offset];
170 170
171 a_Sxw = parent::m_bin_Sxw[offset][0]; 171 a_Sxw = parent::m_bin_Sxw[offset][0];
172 a_Sx2w = parent::m_bin_Sx2w[offset][0]; 172 a_Sx2w = parent::m_bin_Sx2w[offset][0];
173 a_Syw = parent::m_bin_Sxw[offset][1]; 173 a_Syw = parent::m_bin_Sxw[offset][1];
174 a_Sy2w = parent::m_bin_Sx2w[offset][1]; 174 a_Sy2w = parent::m_bin_Sx2w[offset][1];
175 175
176 return true; 176 return true;
177 } 177 }
178 178
179 bool add(const h2& a_histo){ 179 bool add(const h2& a_histo){
180 parent::base_add(a_histo); 180 parent::base_add(a_histo);
181 return true; 181 return true;
182 } 182 }
183 bool subtract(const h2& a_histo){ 183 bool subtract(const h2& a_histo){
184 parent::base_subtract(a_histo); 184 parent::base_subtract(a_histo);
185 return true; 185 return true;
186 } 186 }
187 187
188 bool multiply(const h2& a_histo) { 188 bool multiply(const h2& a_histo) {
189 return parent::base_multiply(a_histo); 189 return parent::base_multiply(a_histo);
190 } 190 }
191 191
192 bool divide(const h2& a_histo) { 192 bool divide(const h2& a_histo) {
193 return parent::base_divide(a_histo); 193 return parent::base_divide(a_histo);
194 } 194 }
195 195
196 bool equals_TH(const h2& a_from,const TW& a_ 196 bool equals_TH(const h2& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
197 if(!parent::equals_TH(a_from,a_prec,a_fabs 197 if(!parent::equals_TH(a_from,a_prec,a_fabs,true)) return false;
198 return true; 198 return true;
199 } 199 }
200 200
201 void not_a_profile() const {} 201 void not_a_profile() const {}
202 202
203 public: //CERN-ROOT API 203 public: //CERN-ROOT API
204 bool Fill(TC aX,TC aY,TW aWeight = 1) {retur 204 bool Fill(TC aX,TC aY,TW aWeight = 1) {return fill(aX,aY,aWeight);}
205 205
206 public: 206 public:
207 h2(const std::string& a_title,bn_t aXnumber, 207 h2(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax,bn_t aYnumber,TC aYmin,TC aYmax)
208 :parent(a_title,aXnumber,aXmin,aXmax,aYnumbe 208 :parent(a_title,aXnumber,aXmin,aXmax,aYnumber,aYmin,aYmax)
209 {} 209 {}
210 h2(const std::string& a_title,const std::vec 210 h2(const std::string& a_title,const std::vector<TC>& a_edges_x,const std::vector<TC>& a_edges_y)
211 :parent(a_title,a_edges_x,a_edges_y) 211 :parent(a_title,a_edges_x,a_edges_y)
212 {} 212 {}
213 213
214 virtual ~h2(){} 214 virtual ~h2(){}
215 public: 215 public:
216 h2(const h2& a_from): parent(a_from){} 216 h2(const h2& a_from): parent(a_from){}
217 h2& operator=(const h2& a_from){ 217 h2& operator=(const h2& a_from){
218 parent::operator=(a_from); 218 parent::operator=(a_from);
219 return *this; 219 return *this;
220 } 220 }
221 }; 221 };
222 222
223 }} 223 }}
224 224
225 #endif 225 #endif
226 226