| 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_h1 4 #ifndef tools_histo_h1
5 #define tools_histo_h1 5 #define tools_histo_h1
6 6
7 #include "b1" 7 #include "b1"
8 8
9 namespace tools { 9 namespace tools {
10 namespace histo { //have that for h1 ? 10 namespace histo { //have that for h1 ?
11 11
12 //TC is for a coordinate. 12 //TC is for a coordinate.
13 //TO is for an offset used to identify a bin. 13 //TO is for an offset used to identify a bin.
14 //TN is for a number of entries. 14 //TN is for a number of entries.
15 //TW is for a weight. 15 //TW is for a weight.
16 //TH is for a height. Should be the same as TW 16 //TH is for a height. Should be the same as TW.
17 17
18 template <class TC,class TO,class TN,class TW, 18 template <class TC,class TO,class TN,class TW,class TH>
19 class h1 : public b1<TC,TO,TN,TW,TH> { 19 class h1 : public b1<TC,TO,TN,TW,TH> {
20 typedef b1<TC,TO,TN,TW,TH> parent; 20 typedef b1<TC,TO,TN,TW,TH> parent;
21 public: 21 public:
22 typedef histo_data<TC,TO,TN,TW> hd_t; 22 typedef histo_data<TC,TO,TN,TW> hd_t;
23 typedef typename parent::bn_t bn_t; 23 typedef typename parent::bn_t bn_t;
24 typedef typename parent::axis_t axis_t; 24 typedef typename parent::axis_t axis_t;
25 protected: 25 protected:
26 virtual TH get_bin_height(TO a_offset) const 26 virtual TH get_bin_height(TO a_offset) const { //TH should be the same as TW
27 return parent::m_bin_Sw[a_offset]; 27 return parent::m_bin_Sw[a_offset];
28 } 28 }
29 public: 29 public:
30 virtual TH bin_error(int aI) const { //TH sh 30 virtual TH bin_error(int aI) const { //TH should be the same as TW
31 TO offset; 31 TO offset;
32 if(!parent::_find_offset(aI,offset)) retur 32 if(!parent::_find_offset(aI,offset)) return 0;
33 return ::sqrt(parent::m_bin_Sw2[offset]); 33 return ::sqrt(parent::m_bin_Sw2[offset]);
34 } 34 }
35 35
36 public: 36 public:
37 bool multiply(TW a_factor){return parent::ba 37 bool multiply(TW a_factor){return parent::base_multiply(a_factor);}
38 bool scale(TW a_factor) {return multiply(a_f 38 bool scale(TW a_factor) {return multiply(a_factor);}
39 39
40 void copy_from_data(const hd_t& a_from) {par 40 void copy_from_data(const hd_t& a_from) {parent::base_from_data(a_from);}
41 hd_t get_histo_data() const {return *this;} 41 hd_t get_histo_data() const {return *this;} //deprecated. Keep it for g4tools.
42 42
43 bool reset() { 43 bool reset() {
44 parent::base_reset(); 44 parent::base_reset();
45 return true; 45 return true;
46 } 46 }
47 47
48 bool fill(TC aX,TW aWeight = 1) { 48 bool fill(TC aX,TW aWeight = 1) {
49 if(parent::m_dimension!=1) return false; 49 if(parent::m_dimension!=1) return false;
50 50
51 bn_t ibin; 51 bn_t ibin;
52 if(!parent::m_axes[0].coord_to_absolute_in 52 if(!parent::m_axes[0].coord_to_absolute_index(aX,ibin)) return false;
53 53
54 TO offset = ibin; 54 TO offset = ibin;
55 55
56 parent::m_bin_entries[offset]++; 56 parent::m_bin_entries[offset]++;
57 parent::m_bin_Sw[offset] += aWeight; 57 parent::m_bin_Sw[offset] += aWeight;
58 parent::m_bin_Sw2[offset] += aWeight * aWe 58 parent::m_bin_Sw2[offset] += aWeight * aWeight;
59 59
60 TC xw = aX * aWeight; 60 TC xw = aX * aWeight;
61 TC x2w = aX * xw; 61 TC x2w = aX * xw;
62 parent::m_bin_Sxw[offset][0] += xw; 62 parent::m_bin_Sxw[offset][0] += xw;
63 parent::m_bin_Sx2w[offset][0] += x2w; 63 parent::m_bin_Sx2w[offset][0] += x2w;
64 64
65 bool inRange = true; 65 bool inRange = true;
66 if(ibin==0) inRange = false; 66 if(ibin==0) inRange = false;
67 else if(ibin==(parent::m_axes[0].m_number_ 67 else if(ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
68 68
69 parent::m_all_entries++; 69 parent::m_all_entries++;
70 if(inRange) { 70 if(inRange) {
71 // fast getters : 71 // fast getters :
72 parent::m_in_range_entries++; 72 parent::m_in_range_entries++;
73 parent::m_in_range_Sw += aWeight; 73 parent::m_in_range_Sw += aWeight;
74 parent::m_in_range_Sw2 += aWeight*aWeigh 74 parent::m_in_range_Sw2 += aWeight*aWeight;
75 75
76 parent::m_in_range_Sxw[0] += xw; 76 parent::m_in_range_Sxw[0] += xw;
77 parent::m_in_range_Sx2w[0] += x2w; 77 parent::m_in_range_Sx2w[0] += x2w;
78 } 78 }
79 79
80 return true; 80 return true;
81 } 81 }
82 82
83 bool set_bin_content(bn_t a_ibin,TN a_entrie 83 bool set_bin_content(bn_t a_ibin,TN a_entries,TW a_Sw,TW a_Sw2,TC a_Sxw,TC a_Sx2w) {
84 if(parent::m_dimension!=1) return false; 84 if(parent::m_dimension!=1) return false;
85 if(a_ibin>(parent::m_axes[0].m_number_of_b 85 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) return false;
86 86
87 bool inRange = true; 87 bool inRange = true;
88 if(a_ibin==0) inRange = false; 88 if(a_ibin==0) inRange = false;
89 else if(a_ibin==(parent::m_axes[0].m_numbe 89 else if(a_ibin==(parent::m_axes[0].m_number_of_bins+1)) inRange = false;
90 90
91 TO offset = a_ibin; 91 TO offset = a_ibin;
92 92
93 parent::m_all_entries -= parent::m_bin_ent 93 parent::m_all_entries -= parent::m_bin_entries[offset];
94 if(inRange) { 94 if(inRange) {
95 parent::m_in_range_entries -= parent::m_ 95 parent::m_in_range_entries -= parent::m_bin_entries[offset];
96 parent::m_in_range_Sw -= parent::m_bin_S 96 parent::m_in_range_Sw -= parent::m_bin_Sw[offset];
97 parent::m_in_range_Sw2 -= parent::m_bin_ 97 parent::m_in_range_Sw2 -= parent::m_bin_Sw2[offset];
98 parent::m_in_range_Sxw[0] -= parent::m_b 98 parent::m_in_range_Sxw[0] -= parent::m_bin_Sxw[offset][0];
99 parent::m_in_range_Sx2w[0] -= parent::m_ 99 parent::m_in_range_Sx2w[0] -= parent::m_bin_Sx2w[offset][0];
100 } 100 }
101 101
102 parent::m_bin_entries[offset] = a_entries; 102 parent::m_bin_entries[offset] = a_entries;
103 parent::m_bin_Sw[offset] = a_Sw; 103 parent::m_bin_Sw[offset] = a_Sw;
104 parent::m_bin_Sw2[offset] = a_Sw2; 104 parent::m_bin_Sw2[offset] = a_Sw2;
105 105
106 parent::m_bin_Sxw[offset][0] = a_Sxw; 106 parent::m_bin_Sxw[offset][0] = a_Sxw;
107 parent::m_bin_Sx2w[offset][0] = a_Sx2w; 107 parent::m_bin_Sx2w[offset][0] = a_Sx2w;
108 108
109 parent::m_all_entries += a_entries; 109 parent::m_all_entries += a_entries;
110 if(inRange) { 110 if(inRange) {
111 parent::m_in_range_entries += a_entries; 111 parent::m_in_range_entries += a_entries;
112 parent::m_in_range_Sw += a_Sw; 112 parent::m_in_range_Sw += a_Sw;
113 parent::m_in_range_Sw2 += a_Sw2; 113 parent::m_in_range_Sw2 += a_Sw2;
114 114
115 parent::m_in_range_Sxw[0] += a_Sxw; 115 parent::m_in_range_Sxw[0] += a_Sxw;
116 parent::m_in_range_Sx2w[0] += a_Sx2w; 116 parent::m_in_range_Sx2w[0] += a_Sx2w;
117 } 117 }
118 118
119 return true; 119 return true;
120 } 120 }
121 121
122 bool get_bin_content(bn_t a_ibin,TN& a_entri 122 bool get_bin_content(bn_t a_ibin,TN& a_entries,TW& a_Sw,TW& a_Sw2,TC& a_Sxw,TC& a_Sx2w) {
123 if(parent::m_dimension!=1) { 123 if(parent::m_dimension!=1) {
124 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 124 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 0;a_Sx2w = 0;
125 return false; 125 return false;
126 } 126 }
127 if(a_ibin>(parent::m_axes[0].m_number_of_b 127 if(a_ibin>(parent::m_axes[0].m_number_of_bins+1)) {
128 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 128 a_entries = 0;a_Sw = 0;a_Sw2 = 0;a_Sxw = 0;a_Sx2w = 0;
129 return false; 129 return false;
130 } 130 }
131 131
132 TO offset = a_ibin; 132 TO offset = a_ibin;
133 133
134 a_entries = parent::m_bin_entries[offset]; 134 a_entries = parent::m_bin_entries[offset];
135 a_Sw = parent::m_bin_Sw[offset]; 135 a_Sw = parent::m_bin_Sw[offset];
136 a_Sw2 = parent::m_bin_Sw2[offset]; 136 a_Sw2 = parent::m_bin_Sw2[offset];
137 137
138 a_Sxw = parent::m_bin_Sxw[offset][0]; 138 a_Sxw = parent::m_bin_Sxw[offset][0];
139 a_Sx2w = parent::m_bin_Sx2w[offset][0]; 139 a_Sx2w = parent::m_bin_Sx2w[offset][0];
140 140
141 return true; 141 return true;
142 } 142 }
143 143
144 bool add(const h1& a_histo){ 144 bool add(const h1& a_histo){
145 parent::base_add(a_histo); 145 parent::base_add(a_histo);
146 return true; 146 return true;
147 } 147 }
148 bool subtract(const h1& a_histo){ 148 bool subtract(const h1& a_histo){
149 parent::base_subtract(a_histo); 149 parent::base_subtract(a_histo);
150 return true; 150 return true;
151 } 151 }
152 152
153 bool multiply(const h1& a_histo) { 153 bool multiply(const h1& a_histo) {
154 return parent::base_multiply(a_histo); 154 return parent::base_multiply(a_histo);
155 } 155 }
156 156
157 bool divide(const h1& a_histo) { 157 bool divide(const h1& a_histo) {
158 return parent::base_divide(a_histo); 158 return parent::base_divide(a_histo);
159 } 159 }
160 160
161 bool gather_bins(unsigned int a_factor) { // 161 bool gather_bins(unsigned int a_factor) { //for exa 2,3.
162 if(!a_factor) return false; 162 if(!a_factor) return false;
163 163
164 // actual bin number must be a multiple of 164 // actual bin number must be a multiple of a_factor.
165 165
166 const axis_t& _axis = parent::axis(); 166 const axis_t& _axis = parent::axis();
167 167
168 bn_t n = _axis.bins(); 168 bn_t n = _axis.bins();
169 if(!n) return false; 169 if(!n) return false;
170 170
171 bn_t new_n = n/a_factor; 171 bn_t new_n = n/a_factor;
172 if(a_factor*new_n!=n) return false; 172 if(a_factor*new_n!=n) return false;
173 173
174 h1* new_h = 0; 174 h1* new_h = 0;
175 if(_axis.is_fixed_binning()) { 175 if(_axis.is_fixed_binning()) {
176 new_h = new h1(parent::m_title,new_n,_ax 176 new_h = new h1(parent::m_title,new_n,_axis.lower_edge(),_axis.upper_edge());
177 } else { 177 } else {
178 const std::vector<TC>& _edges = _axis.ed 178 const std::vector<TC>& _edges = _axis.edges();
179 std::vector<TC> new_edges(new_n+1); 179 std::vector<TC> new_edges(new_n+1);
180 for(bn_t ibin=0;ibin<new_n;ibin++) { 180 for(bn_t ibin=0;ibin<new_n;ibin++) {
181 new_edges[ibin] = _edges[ibin*a_factor 181 new_edges[ibin] = _edges[ibin*a_factor];
182 } 182 }
183 new_edges[new_n] = _edges[n]; //upper ed 183 new_edges[new_n] = _edges[n]; //upper edge.
184 new_h = new h1(parent::m_title,new_edges 184 new_h = new h1(parent::m_title,new_edges);
185 } 185 }
186 if(!new_h) return false; 186 if(!new_h) return false;
187 187
188 TO offset,new_offset,offac; 188 TO offset,new_offset,offac;
189 for(bn_t ibin=0;ibin<new_n;ibin++) { 189 for(bn_t ibin=0;ibin<new_n;ibin++) {
190 new_offset = ibin+1; 190 new_offset = ibin+1;
191 offset = a_factor*ibin+1; 191 offset = a_factor*ibin+1;
192 for(unsigned int ifac=0;ifac<a_factor;if 192 for(unsigned int ifac=0;ifac<a_factor;ifac++) {
193 offac = offset+ifac; 193 offac = offset+ifac;
194 new_h->m_bin_entries[new_offset] += pa 194 new_h->m_bin_entries[new_offset] += parent::m_bin_entries[offac];
195 new_h->m_bin_Sw[new_offset] += parent: 195 new_h->m_bin_Sw[new_offset] += parent::m_bin_Sw[offac];
196 new_h->m_bin_Sw2[new_offset] += parent 196 new_h->m_bin_Sw2[new_offset] += parent::m_bin_Sw2[offac];
197 new_h->m_bin_Sxw[new_offset][0] += par 197 new_h->m_bin_Sxw[new_offset][0] += parent::m_bin_Sxw[offac][0];
198 new_h->m_bin_Sx2w[new_offset][0] += pa 198 new_h->m_bin_Sx2w[new_offset][0] += parent::m_bin_Sx2w[offac][0];
199 } 199 }
200 } 200 }
201 201
202 //underflow : 202 //underflow :
203 new_offset = 0; 203 new_offset = 0;
204 offac = 0; 204 offac = 0;
205 new_h->m_bin_entries[new_offset] = parent: 205 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac];
206 new_h->m_bin_Sw[new_offset] = parent::m_bi 206 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac];
207 new_h->m_bin_Sw2[new_offset] = parent::m_b 207 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac];
208 new_h->m_bin_Sxw[new_offset][0] = parent:: 208 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0];
209 new_h->m_bin_Sx2w[new_offset][0] = parent: 209 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0];
210 210
211 //overflow : 211 //overflow :
212 new_offset = new_n+1; 212 new_offset = new_n+1;
213 offac = n+1; 213 offac = n+1;
214 new_h->m_bin_entries[new_offset] = parent: 214 new_h->m_bin_entries[new_offset] = parent::m_bin_entries[offac];
215 new_h->m_bin_Sw[new_offset] = parent::m_bi 215 new_h->m_bin_Sw[new_offset] = parent::m_bin_Sw[offac];
216 new_h->m_bin_Sw2[new_offset] = parent::m_b 216 new_h->m_bin_Sw2[new_offset] = parent::m_bin_Sw2[offac];
217 new_h->m_bin_Sxw[new_offset][0] = parent:: 217 new_h->m_bin_Sxw[new_offset][0] = parent::m_bin_Sxw[offac][0];
218 new_h->m_bin_Sx2w[new_offset][0] = parent: 218 new_h->m_bin_Sx2w[new_offset][0] = parent::m_bin_Sx2w[offac][0];
219 219
220 *this = *new_h; 220 *this = *new_h;
221 return true; 221 return true;
222 } 222 }
223 223
224 bool equals_TH(const h1& a_from,const TW& a_ 224 bool equals_TH(const h1& a_from,const TW& a_prec,TW(*a_fabs)(TW)) const {
225 if(!parent::equals_TH(a_from,a_prec,a_fabs 225 if(!parent::equals_TH(a_from,a_prec,a_fabs,true)) return false;
226 return true; 226 return true;
227 } 227 }
228 228
229 void not_a_profile() const {} 229 void not_a_profile() const {}
230 230
231 public: //CERN-ROOT API (for MEMPHYS sim). 231 public: //CERN-ROOT API (for MEMPHYS sim).
232 bool Fill(TC aX,TW aWeight = 1) {return fill 232 bool Fill(TC aX,TW aWeight = 1) {return fill(aX,aWeight);}
233 233
234 public: 234 public:
235 h1(const std::string& a_title,bn_t aXnumber, 235 h1(const std::string& a_title,bn_t aXnumber,TC aXmin,TC aXmax)
236 :parent(a_title,aXnumber,aXmin,aXmax){} 236 :parent(a_title,aXnumber,aXmin,aXmax){}
237 237
238 h1(const std::string& a_title,const std::vec 238 h1(const std::string& a_title,const std::vector<TC>& a_edges)
239 :parent(a_title,a_edges){} 239 :parent(a_title,a_edges){}
240 240
241 virtual ~h1(){} 241 virtual ~h1(){}
242 public: 242 public:
243 h1(const h1& a_from):parent(a_from){} 243 h1(const h1& a_from):parent(a_from){}
244 h1& operator=(const h1& a_from){ 244 h1& operator=(const h1& a_from){
245 if(&a_from==this) return *this; 245 if(&a_from==this) return *this;
246 parent::operator=(a_from); 246 parent::operator=(a_from);
247 return *this; 247 return *this;
248 } 248 }
249 }; 249 };
250 250
251 }} 251 }}
252 252
253 #endif 253 #endif
254 254
255 255
256 256
257 257