Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 2 // See the file tools.license for terms. 3 4 #ifndef tools_histo_h3 5 #define tools_histo_h3 6 7 #include "b3" 8 9 namespace tools { 10 namespace histo { 11 12 template <class TC,class TO,class TN,class TW, 13 class h3 : public b3<TC,TO,TN,TW,TH> { 14 typedef b3<TC,TO,TN,TW,TH> parent; 15 public: 16 typedef histo_data<TC,TO,TN,TW> hd_t; 17 typedef typename parent::bn_t bn_t; 18 protected: 19 virtual TH get_bin_height(TO a_offset) const 20 return parent::m_bin_Sw[a_offset]; 21 } 22 public: 23 24 virtual TH bin_error(int aI,int aJ,int aK) c 25 TO offset; 26 if(!parent::_find_offset(aI,aJ,aK,offset)) 27 return ::sqrt(parent::m_bin_Sw2[offset]); 28 } 29 30 public: 31 bool multiply(TW a_factor){return parent::ba 32 bool scale(TW a_factor) {return multiply(a_f 33 34 void copy_from_data(const hd_t& a_from) {par 35 hd_t get_histo_data() const {return *this;} 36 37 bool reset() { 38 parent::base_reset(); 39 return true; 40 } 41 42 bool fill(TC aX,TC aY,TC aZ,TW aWeight = 1) 43 if(parent::m_dimension!=3) return false; 44 45 bn_t ibin,jbin,kbin; 46 if(!parent::m_axes[0].coord_to_absolute_in 47 if(!parent::m_axes[1].coord_to_absolute_in 48 if(!parent::m_axes[2].coord_to_absolute_in 49 TO offset = ibin + jbin * parent::m_axes[1 50 51 parent::m_bin_entries[offset]++; 52 parent::m_bin_Sw[offset] += aWeight; 53 parent::m_bin_Sw2[offset] += aWeight * aWe 54 55 TC xw = aX * aWeight; 56 TC x2w = aX * xw; 57 parent::m_bin_Sxw[offset][0] += xw; 58 parent::m_bin_Sx2w[offset][0] += x2w; 59 60 TC yw = aY * aWeight; 61 TC y2w = aY * yw; 62 parent::m_bin_Sxw[offset][1] += yw; 63 parent::m_bin_Sx2w[offset][1] += y2w; 64 65 TC zw = aZ * aWeight; 66 TC z2w = aZ * zw; 67 parent::m_bin_Sxw[offset][2] += zw; 68 parent::m_bin_Sx2w[offset][2] += z2w; 69 70 bool inRange = true; 71 if(ibin==0) inRange = false; 72 else if(ibin==(parent::m_axes[0].m_number_ 73 74 if(jbin==0) inRange = false; 75 else if(jbin==(parent::m_axes[1].m_number_ 76 77 if(kbin==0) inRange = false; 78 else if(kbin==(parent::m_axes[2].m_number_ 79 80 parent::m_all_entries++; 81 if(inRange) { 82 parent::m_in_range_plane_Sxyw[0] += aX * 83 parent::m_in_range_plane_Sxyw[1] += aY * 84 parent::m_in_range_plane_Sxyw[2] += aZ * 85 86 // fast getters : 87 parent::m_in_range_entries++; 88 parent::m_in_range_Sw += aWeight; 89 parent::m_in_range_Sw2 += aWeight*aWeigh 90 91 parent::m_in_range_Sxw[0] += xw; 92 parent::m_in_range_Sx2w[0] += x2w; 93 94 parent::m_in_range_Sxw[1] += yw; 95 parent::m_in_range_Sx2w[1] += y2w; 96 97 parent::m_in_range_Sxw[2] += zw; 98 parent::m_in_range_Sx2w[2] += z2w; 99 } 100 101 return true; 102 } 103 104 bool set_bin_content(bn_t a_ibin,bn_t a_jbin 105 TN a_entries,TW a_Sw,TW 106 TC a_Sxw,TC a_Sx2w,TC a 107 if(parent::m_dimension!=3) return false; 108 if(a_ibin>(parent::m_axes[0].m_number_of_b 109 if(a_jbin>(parent::m_axes[1].m_number_of_b 110 if(a_kbin>(parent::m_axes[2].m_number_of_b 111 112 bool inRange = true; 113 if(a_ibin==0) inRange = false; 114 else if(a_ibin==(parent::m_axes[0].m_numbe 115 if(a_jbin==0) inRange = false; 116 else if(a_jbin==(parent::m_axes[1].m_numbe 117 if(a_kbin==0) inRange = false; 118 else if(a_kbin==(parent::m_axes[2].m_numbe 119 120 TO offset = a_ibin + a_jbin * parent::m_ax 121 122 parent::m_all_entries -= parent::m_bin_ent 123 if(inRange) { 124 parent::m_in_range_entries -= parent::m_ 125 parent::m_in_range_Sw -= parent::m_bin_S 126 parent::m_in_range_Sw2 -= parent::m_bin_ 127 parent::m_in_range_Sxw[0] -= parent::m_b 128 parent::m_in_range_Sx2w[0] -= parent::m_ 129 parent::m_in_range_Sxw[1] -= parent::m_b 130 parent::m_in_range_Sx2w[1] -= parent::m_ 131 parent::m_in_range_Sxw[2] -= parent::m_b 132 parent::m_in_range_Sx2w[2] -= parent::m_ 133 } 134 135 parent::m_bin_entries[offset] = a_entries; 136 parent::m_bin_Sw[offset] = a_Sw; 137 parent::m_bin_Sw2[offset] = a_Sw2; 138 139 parent::m_bin_Sxw[offset][0] = a_Sxw; 140 parent::m_bin_Sx2w[offset][0] = a_Sx2w; 141 parent::m_bin_Sxw[offset][1] = a_Syw; 142 parent::m_bin_Sx2w[offset][1] = a_Sy2w; 143 parent::m_bin_Sxw[offset][2] = a_Szw; 144 parent::m_bin_Sx2w[offset][2] = a_Sz2w; 145 146 parent::m_all_entries += a_entries; 147 if(inRange) { 148 //parent::m_in_range_plane_Sxyw[0,1,2] ? 149 150 parent::m_in_range_entries += a_entries; 151 parent::m_in_range_Sw += a_Sw; 152 parent::m_in_range_Sw2 += a_Sw2; 153 154 parent::m_in_range_Sxw[0] += a_Sxw; 155 parent::m_in_range_Sx2w[0] += a_Sx2w; 156 parent::m_in_range_Sxw[1] += a_Syw; 157 parent::m_in_range_Sx2w[1] += a_Sy2w; 158 parent::m_in_range_Sxw[2] += a_Szw; 159 parent::m_in_range_Sx2w[2] += a_Sz2w; 160 } 161 162 return true; 163 } 164 165 bool get_bin_content(bn_t a_ibin,bn_t a_jbin 166 TN& a_entries,TW& a_Sw, 167 TC& a_Sxw,TC& a_Sx2w, 168 TC& a_Syw,TC& a_Sy2w, 169 TC& a_Szw,TC& a_Sz2w) { 170 if(parent::m_dimension!=3) { 171 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 172 a_Sxw = 0;a_Sx2w = 0; 173 a_Syw = 0;a_Sy2w = 0; 174 a_Szw = 0;a_Sz2w = 0; 175 return false; 176 } 177 if(a_ibin>(parent::m_axes[0].m_number_of_b 178 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 179 a_Sxw = 0;a_Sx2w = 0; 180 a_Syw = 0;a_Sy2w = 0; 181 a_Szw = 0;a_Sz2w = 0; 182 return false; 183 } 184 if(a_jbin>(parent::m_axes[1].m_number_of_b 185 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 186 a_Sxw = 0;a_Sx2w = 0; 187 a_Syw = 0;a_Sy2w = 0; 188 a_Szw = 0;a_Sz2w = 0; 189 return false; 190 } 191 if(a_kbin>(parent::m_axes[2].m_number_of_b 192 a_entries = 0;a_Sw = 0;a_Sw2 = 0; 193 a_Sxw = 0;a_Sx2w = 0; 194 a_Syw = 0;a_Sy2w = 0; 195 a_Szw = 0;a_Sz2w = 0; 196 return false; 197 } 198 199 TO offset = a_ibin + a_jbin * parent::m_ax 200 201 a_entries = parent::m_bin_entries[offset]; 202 a_Sw = parent::m_bin_Sw[offset]; 203 a_Sw2 = parent::m_bin_Sw2[offset]; 204 205 a_Sxw = parent::m_bin_Sxw[offset][0]; 206 a_Sx2w = parent::m_bin_Sx2w[offset][0]; 207 a_Syw = parent::m_bin_Sxw[offset][1]; 208 a_Sy2w = parent::m_bin_Sx2w[offset][1]; 209 a_Szw = parent::m_bin_Sxw[offset][2]; 210 a_Sz2w = parent::m_bin_Sx2w[offset][2]; 211 212 return true; 213 } 214 215 bool add(const h3& a_histo){ 216 parent::base_add(a_histo); 217 return true; 218 } 219 bool subtract(const h3& a_histo){ 220 parent::base_subtract(a_histo); 221 return true; 222 } 223 224 bool multiply(const h3& a_histo) { 225 return parent::base_multiply(a_histo); 226 } 227 228 bool divide(const h3& a_histo) { 229 return parent::base_divide(a_histo); 230 } 231 232 bool equals_TH(const h3& a_from,const TW& a_ 233 if(!parent::equals_TH(a_from,a_prec,a_fabs 234 return true; 235 } 236 237 void not_a_profile() const {} 238 239 public: //CERN-ROOT API 240 bool Fill(TC aX,TC aY,TC aZ,TW aWeight = 1) 241 242 public: 243 /* 244 // Slices : 245 h2d* slice_xy(int aKbeg,int aKend) const; 246 h2d* slice_yz(int aIbeg,int aIend) const; 247 h2d* slice_xz(int aJbeg,int aJend) const; 248 */ 249 public: 250 h3(const std::string& a_title, 251 bn_t aXnumber,TC aXmin,TC aXmax, 252 bn_t aYnumber,TC aYmin,TC aYmax, 253 bn_t aZnumber,TC aZmin,TC aZmax) 254 :parent(a_title,aXnumber,aXmin,aXmax, 255 aYnumber,aYmin,aYmax, 256 aZnumber,aZmin,aZmax) 257 {} 258 259 h3(const std::string& a_title, 260 const std::vector<TC>& a_edges_x, 261 const std::vector<TC>& a_edges_y, 262 const std::vector<TC>& a_edges_z) 263 :parent(a_title,a_edges_x,a_edges_y,a_edges_ 264 {} 265 266 virtual ~h3(){} 267 public: 268 h3(const h3& a_from): parent(a_from){} 269 h3& operator=(const h3& a_from){ 270 parent::operator=(a_from); 271 return *this; 272 } 273 }; 274 275 }} 276 277 #endif 278 279 280 281