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