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