| 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_c1d 4 #ifndef tools_histo_c1d
5 #define tools_histo_c1d 5 #define tools_histo_c1d
6 6
7 #include "base_cloud" 7 #include "base_cloud"
8 8
9 #include "../mnmx" 9 #include "../mnmx"
10 10
11 #include "h1d" 11 #include "h1d"
12 12
13 namespace tools { 13 namespace tools {
14 namespace histo { 14 namespace histo {
15 15
16 class c1d : public base_cloud { 16 class c1d : public base_cloud {
17 public: 17 public:
18 static const std::string& s_class() { 18 static const std::string& s_class() {
19 static const std::string s_v("tools::histo 19 static const std::string s_v("tools::histo::c1d");
20 return s_v; 20 return s_v;
21 } 21 }
22 public: 22 public:
23 bool set_title(const std::string& a_title){ 23 bool set_title(const std::string& a_title){
24 m_title = a_title; 24 m_title = a_title;
25 if(m_histo) m_histo->set_title(a_title); 25 if(m_histo) m_histo->set_title(a_title);
26 return true; 26 return true;
27 } 27 }
28 28
29 unsigned int dimension() const {return 1;} 29 unsigned int dimension() const {return 1;}
30 bool reset() { 30 bool reset() {
31 clear(); 31 clear();
32 delete m_histo; 32 delete m_histo;
33 m_histo = 0; 33 m_histo = 0;
34 return true; 34 return true;
35 } 35 }
36 unsigned int entries() const { 36 unsigned int entries() const {
37 return m_histo ? m_histo->all_entries() : 37 return m_histo ? m_histo->all_entries() : (unsigned int)m_ws.size();
38 } 38 }
39 public: 39 public:
40 double sum_of_weights() const { 40 double sum_of_weights() const {
41 return (m_histo ? m_histo->sum_bin_heights 41 return (m_histo ? m_histo->sum_bin_heights() : m_Sw);
42 } 42 }
43 43
44 bool convert_to_histogram(){ 44 bool convert_to_histogram(){
45 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cn 45 if( (m_cnv_x_num<=0) || (m_cnv_x_max<=m_cnv_x_min) ) {
46 // Cloud min, max should be included in 46 // Cloud min, max should be included in the histo.
47 double dx = 0.01 * (upper_edge() - lower 47 double dx = 0.01 * (upper_edge() - lower_edge())/BINS();
48 return convert(BINS(),lower_edge(),upper 48 return convert(BINS(),lower_edge(),upper_edge() + dx);
49 } else { 49 } else {
50 return convert(m_cnv_x_num,m_cnv_x_min,m 50 return convert(m_cnv_x_num,m_cnv_x_min,m_cnv_x_max);
51 } 51 }
52 } 52 }
53 53
54 bool is_converted() const {return m_histo ? 54 bool is_converted() const {return m_histo ? true : false;}
55 bool scale(double a_scale) { 55 bool scale(double a_scale) {
56 if(m_histo) { 56 if(m_histo) {
57 return m_histo->scale(a_scale); 57 return m_histo->scale(a_scale);
58 } else { 58 } else {
59 size_t number = m_ws.size(); 59 size_t number = m_ws.size();
60 for(size_t index=0;index<number;index++) 60 for(size_t index=0;index<number;index++) m_ws[index] *= a_scale;
61 m_Sw *= a_scale; 61 m_Sw *= a_scale;
62 m_Sxw *= a_scale; 62 m_Sxw *= a_scale;
63 m_Sx2w *= a_scale; 63 m_Sx2w *= a_scale;
64 return true; 64 return true;
65 } 65 }
66 } 66 }
67 67
68 bool set_histogram(h1d* a_histo){ //we take 68 bool set_histogram(h1d* a_histo){ //we take ownership of a_histo.
69 reset(); 69 reset();
70 m_histo = a_histo; 70 m_histo = a_histo;
71 return true; 71 return true;
72 } 72 }
73 73
74 public: 74 public:
75 bool fill(double aX,double aW = 1){ 75 bool fill(double aX,double aW = 1){
76 if(!m_histo && (m_limit!=UNLIMITED()) && 76 if(!m_histo && (m_limit!=UNLIMITED()) &&
77 ((int)m_xs.size()>=m_limit)){ 77 ((int)m_xs.size()>=m_limit)){
78 convert_to_histogram(); 78 convert_to_histogram();
79 } 79 }
80 80
81 if(m_histo) { 81 if(m_histo) {
82 return m_histo->fill(aX,aW); 82 return m_histo->fill(aX,aW);
83 } else { 83 } else {
84 if(m_xs.size()) { 84 if(m_xs.size()) {
85 m_lower_x = mn<double>(aX,m_lower_x); 85 m_lower_x = mn<double>(aX,m_lower_x);
86 m_upper_x = mx<double>(aX,m_upper_x); 86 m_upper_x = mx<double>(aX,m_upper_x);
87 } else { 87 } else {
88 m_lower_x = aX; 88 m_lower_x = aX;
89 m_upper_x = aX; 89 m_upper_x = aX;
90 } 90 }
91 m_xs.push_back(aX); 91 m_xs.push_back(aX);
92 m_ws.push_back(aW); 92 m_ws.push_back(aW);
93 m_Sw += aW; 93 m_Sw += aW;
94 double xw = aX * aW; 94 double xw = aX * aW;
95 m_Sxw += xw; 95 m_Sxw += xw;
96 m_Sx2w += aX * xw; 96 m_Sx2w += aX * xw;
97 return true; 97 return true;
98 } 98 }
99 } 99 }
100 100
101 double lower_edge() const { 101 double lower_edge() const {
102 return (m_histo ? m_histo->axis().lower_ed 102 return (m_histo ? m_histo->axis().lower_edge() : m_lower_x);
103 } 103 }
104 double upper_edge() const { 104 double upper_edge() const {
105 return (m_histo ? m_histo->axis().upper_ed 105 return (m_histo ? m_histo->axis().upper_edge() : m_upper_x);
106 } 106 }
107 double value(unsigned int a_index) const {re 107 double value(unsigned int a_index) const {return (m_histo ?0:m_xs[a_index]);}
108 double weight(unsigned int a_index) const {r 108 double weight(unsigned int a_index) const {return (m_histo ?0:m_ws[a_index]);}
109 double mean() const { 109 double mean() const {
110 return (m_histo ? m_histo->mean() : (m_Sw? 110 return (m_histo ? m_histo->mean() : (m_Sw?m_Sxw/m_Sw:0));
111 } 111 }
112 double rms() const { 112 double rms() const {
113 double _rms = 0; //FIXME nan. 113 double _rms = 0; //FIXME nan.
114 if(m_histo) { 114 if(m_histo) {
115 _rms = m_histo->rms(); 115 _rms = m_histo->rms();
116 } else { 116 } else {
117 if(m_Sw==0) { 117 if(m_Sw==0) {
118 } else { 118 } else {
119 double _mean = m_Sxw / m_Sw; 119 double _mean = m_Sxw / m_Sw;
120 _rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) 120 _rms = ::sqrt(::fabs( (m_Sx2w / m_Sw) - _mean * _mean));
121 } 121 }
122 } 122 }
123 return _rms; 123 return _rms;
124 } 124 }
125 125
126 bool convert(unsigned int a_bins,double a_lo 126 bool convert(unsigned int a_bins,double a_lower_edge,double a_upper_edge){
127 if(m_histo) return true; 127 if(m_histo) return true;
128 m_histo = new histo::h1d(base_cloud::title 128 m_histo = new histo::h1d(base_cloud::title(),a_bins,a_lower_edge,a_upper_edge);
129 if(!m_histo) return false; 129 if(!m_histo) return false;
130 bool status = fill_histogram(*m_histo); 130 bool status = fill_histogram(*m_histo);
131 clear(); 131 clear();
132 return status; 132 return status;
133 } 133 }
134 134
135 bool convert(const std::vector<double>& a_ed 135 bool convert(const std::vector<double>& a_edges) {
136 if(m_histo) return true; 136 if(m_histo) return true;
137 m_histo = new histo::h1d(base_cloud::title 137 m_histo = new histo::h1d(base_cloud::title(),a_edges);
138 if(!m_histo) return false; 138 if(!m_histo) return false;
139 bool status = fill_histogram(*m_histo); 139 bool status = fill_histogram(*m_histo);
140 clear(); 140 clear();
141 return status; 141 return status;
142 } 142 }
143 143
144 const histo::h1d& histogram() const { 144 const histo::h1d& histogram() const {
145 if(!m_histo) const_cast<c1d&>(*this).conve 145 if(!m_histo) const_cast<c1d&>(*this).convert_to_histogram();
146 return *m_histo; 146 return *m_histo;
147 } 147 }
148 template <class HISTO> 148 template <class HISTO>
149 bool fill_histogram(HISTO& a_histo) const { 149 bool fill_histogram(HISTO& a_histo) const {
150 size_t number = m_xs.size(); 150 size_t number = m_xs.size();
151 for(size_t index=0;index<number;index++) { 151 for(size_t index=0;index<number;index++) {
152 if(!a_histo.fill(m_xs[index],m_ws[index] 152 if(!a_histo.fill(m_xs[index],m_ws[index])) return false;
153 } 153 }
154 return true; 154 return true;
155 } 155 }
156 bool set_conversion_parameters(unsigned int 156 bool set_conversion_parameters(unsigned int aCnvXnumber,
157 double aCnvXm 157 double aCnvXmin,double aCnvXmax){
158 m_cnv_x_num = aCnvXnumber; 158 m_cnv_x_num = aCnvXnumber;
159 m_cnv_x_min = aCnvXmin; 159 m_cnv_x_min = aCnvXmin;
160 m_cnv_x_max = aCnvXmax; 160 m_cnv_x_max = aCnvXmax;
161 return true; 161 return true;
162 } 162 }
163 163
164 public: 164 public:
165 c1d() 165 c1d()
166 :base_cloud(UNLIMITED()) 166 :base_cloud(UNLIMITED())
167 ,m_lower_x(0),m_upper_x(0) 167 ,m_lower_x(0),m_upper_x(0)
168 ,m_Sxw(0),m_Sx2w(0) 168 ,m_Sxw(0),m_Sx2w(0)
169 ,m_cnv_x_num(0),m_cnv_x_min(0),m_cnv_x_max(0 169 ,m_cnv_x_num(0),m_cnv_x_min(0),m_cnv_x_max(0),m_histo(0)
170 {} 170 {}
171 171
172 c1d(const std::string& a_title,int aLimit = 172 c1d(const std::string& a_title,int aLimit = base_cloud::UNLIMITED())
173 :base_cloud(aLimit) 173 :base_cloud(aLimit)
174 ,m_lower_x(0),m_upper_x(0) 174 ,m_lower_x(0),m_upper_x(0)
175 ,m_Sxw(0),m_Sx2w(0) 175 ,m_Sxw(0),m_Sx2w(0)
176 ,m_cnv_x_num(0),m_cnv_x_min(0),m_cnv_x_max(0 176 ,m_cnv_x_num(0),m_cnv_x_min(0),m_cnv_x_max(0),m_histo(0)
177 { 177 {
178 set_title(a_title); 178 set_title(a_title);
179 } 179 }
180 180
181 virtual ~c1d(){delete m_histo;} 181 virtual ~c1d(){delete m_histo;}
182 public: 182 public:
183 c1d(const c1d& a_from) 183 c1d(const c1d& a_from)
184 :base_cloud(a_from) 184 :base_cloud(a_from)
185 ,m_xs(a_from.m_xs) 185 ,m_xs(a_from.m_xs)
186 ,m_lower_x(a_from.m_lower_x) 186 ,m_lower_x(a_from.m_lower_x)
187 ,m_upper_x(a_from.m_upper_x) 187 ,m_upper_x(a_from.m_upper_x)
188 ,m_Sxw(a_from.m_Sxw) 188 ,m_Sxw(a_from.m_Sxw)
189 ,m_Sx2w(a_from.m_Sx2w) 189 ,m_Sx2w(a_from.m_Sx2w)
190 ,m_cnv_x_num(a_from.m_cnv_x_num) 190 ,m_cnv_x_num(a_from.m_cnv_x_num)
191 ,m_cnv_x_min(a_from.m_cnv_x_min) 191 ,m_cnv_x_min(a_from.m_cnv_x_min)
192 ,m_cnv_x_max(a_from.m_cnv_x_max) 192 ,m_cnv_x_max(a_from.m_cnv_x_max)
193 ,m_histo(0) 193 ,m_histo(0)
194 { 194 {
195 if(a_from.m_histo) { 195 if(a_from.m_histo) {
196 m_histo = new histo::h1d(*a_from.m_histo 196 m_histo = new histo::h1d(*a_from.m_histo);
197 } 197 }
198 } 198 }
199 199
200 c1d& operator=(const c1d& a_from){ 200 c1d& operator=(const c1d& a_from){
201 base_cloud::operator=(a_from); 201 base_cloud::operator=(a_from);
202 if(&a_from==this) return *this; 202 if(&a_from==this) return *this;
203 m_xs = a_from.m_xs; 203 m_xs = a_from.m_xs;
204 m_lower_x = a_from.m_lower_x; 204 m_lower_x = a_from.m_lower_x;
205 m_upper_x = a_from.m_upper_x; 205 m_upper_x = a_from.m_upper_x;
206 m_Sxw = a_from.m_Sxw; 206 m_Sxw = a_from.m_Sxw;
207 m_Sx2w = a_from.m_Sx2w; 207 m_Sx2w = a_from.m_Sx2w;
208 m_cnv_x_num = a_from.m_cnv_x_num; 208 m_cnv_x_num = a_from.m_cnv_x_num;
209 m_cnv_x_min = a_from.m_cnv_x_min; 209 m_cnv_x_min = a_from.m_cnv_x_min;
210 m_cnv_x_max = a_from.m_cnv_x_max; 210 m_cnv_x_max = a_from.m_cnv_x_max;
211 delete m_histo; 211 delete m_histo;
212 m_histo = 0; 212 m_histo = 0;
213 if(a_from.m_histo) { 213 if(a_from.m_histo) {
214 m_histo = new histo::h1d(*a_from.m_histo 214 m_histo = new histo::h1d(*a_from.m_histo);
215 } 215 }
216 return *this; 216 return *this;
217 } 217 }
218 218
219 public: //AIDA API 219 public: //AIDA API
220 double lowerEdge() const {return lower_edge( 220 double lowerEdge() const {return lower_edge();}
221 double upperEdge() const {return upper_edge( 221 double upperEdge() const {return upper_edge();}
222 template <class HISTO> 222 template <class HISTO>
223 bool fillHistogram(HISTO& a_histo) const {re 223 bool fillHistogram(HISTO& a_histo) const {return fill_histogram<HISTO>(a_histo);}
224 protected: 224 protected:
225 void clear(){ 225 void clear(){
226 m_lower_x = 0; 226 m_lower_x = 0;
227 m_upper_x = 0; 227 m_upper_x = 0;
228 m_Sw = 0; 228 m_Sw = 0;
229 m_Sxw = 0; 229 m_Sxw = 0;
230 m_Sx2w = 0; 230 m_Sx2w = 0;
231 m_xs.clear(); 231 m_xs.clear();
232 m_ws.clear(); 232 m_ws.clear();
233 } 233 }
234 234
235 protected: 235 protected:
236 std::vector<double> m_xs; 236 std::vector<double> m_xs;
237 double m_lower_x; 237 double m_lower_x;
238 double m_upper_x; 238 double m_upper_x;
239 double m_Sxw; 239 double m_Sxw;
240 double m_Sx2w; 240 double m_Sx2w;
241 // 241 //
242 unsigned int m_cnv_x_num; 242 unsigned int m_cnv_x_num;
243 double m_cnv_x_min; 243 double m_cnv_x_min;
244 double m_cnv_x_max; 244 double m_cnv_x_max;
245 histo::h1d* m_histo; 245 histo::h1d* m_histo;
246 }; 246 };
247 247
248 }} 248 }}
249 249
250 #endif 250 #endif