SourceXtractorPlusPlus  0.16
Please provide a description of the project.
BackgroundCell.cpp
Go to the documentation of this file.
1 
17 /*
18  * Created on Jan 05, 2015
19  * @author: mkuemmel@usm.lmu.de
20  *
21  * Date: $Date$
22  * Revision: $Revision$
23  * Author: $Author$
24  */
25 #include <cmath>
26 #include <limits>
27 #include "ElementsKernel/Exception.h" // for Elements Exception
31 
32 namespace SourceXtractor {
33 
34 BackgroundCell::BackgroundCell(const PIXTYPE* cellData, const size_t ndata, const PIXTYPE* cellWeight, const PIXTYPE weightThresh)
35 {
36  itsCellData=cellData;
37  itsCellWeight=cellWeight;
38  itsNdata=ndata;
39  itsWeightThresh=weightThresh;
40 
41  // compute some basic statistics on the
42  // cell data and weights: mean, sigma (with a clipping), npix
43  if (cellWeight){
45  itsHasWeight=true;
48  }
49  else{
51  itsHasWeight=false;
53  }
54 }
55 
57 {
58  delete itsHisto;
59  itsHisto=NULL;
60  delete itsWeightHisto;
61  itsWeightHisto=NULL;
62 }
63 
65 {
66  PIXTYPE pixVal=0.0;
67 
68  meanVal = itsHisto->itsMean;
69  sigmaVal = itsHisto->itsSigma;
70 
71  // check if something can be done at all
72  if (itsHisto->itsMean<=-BIG && itsHisto->itsSigma<=-BIG){
73  return;
74  }
75 
76  // go over the data
77  for (size_t index=0; index<itsNdata; index++) {
78 
79  // take the current value
80  pixVal = itsCellData[index];
81 
82  // add to the histogram
83  itsHisto->addDatum(pixVal);
84  }
85 
86  // determine the mean value and sigma
87  itsHisto->getBackGuess(meanVal, sigmaVal);
88 
89  return;
90 }
91 
92 void BackgroundCell::getBackgroundValues(PIXTYPE& meanVal, PIXTYPE& sigmaVal, PIXTYPE& whtMeanVal, PIXTYPE& whtSigmaVal)
93 {
94  PIXTYPE pixVal=0.0;
95  PIXTYPE whtVal=0.0;
96 
97  // transfer all mean ans sigma values
98  meanVal = itsHisto->itsMean;
99  sigmaVal = itsHisto->itsSigma;
100  if (itsHasWeight){
101  whtMeanVal = itsWeightHisto->itsMean;
102  whtSigmaVal = itsWeightHisto->itsSigma;
103  }
104 
105  // check if something can be done at all
106  if (itsHisto->itsMean<=-BIG && itsHisto->itsSigma<=-BIG){
107  return;
108  }
109 
110  if (itsHasWeight){
111 
112  // go over the data
113  for (size_t index=0; index<itsNdata; index++) {
114 
115  // take the current values
116  pixVal = itsCellData[index];
117  whtVal = itsCellWeight[index];
118  if (whtVal<itsWeightThresh){
119  // add to the histograms
120  itsHisto->addDatum(pixVal);
121  itsWeightHisto->addDatum(whtVal);
122  }
123  }
124  }
125  else {
126  // go over the data
127  for (size_t index=0; index<itsNdata; index++) {
128 
129  // take the current value
130  pixVal = itsCellData[index];
131 
132  // add to the histogram
133  itsHisto->addDatum(pixVal);
134  }
135  }
136 
137  // determine the mean value and sigma
138  itsHisto->getBackGuess(meanVal, sigmaVal);
139 
140  if (itsHasWeight){
141  itsWeightHisto->getBackGuess(whtMeanVal, whtSigmaVal);
142  }
143 
144  return;
145 }
146 
148 {
149  BackgroundHistogram* theHisto=NULL;
150 
151  // check if something can be done at all
152  if (itsMean<=-BIG && itsSigma<=-BIG)
153  {
154  bckVal = itsMean;
155  sigmaVal = itsSigma;
156  return;
157  }
158  else
159  {
160  bckVal = itsMean;
161  sigmaVal = itsSigma;
162  }
163 
164  // create the histogram and fill it
166  //theHisto->fillInData(itsCellData, itsNdata);
167  theHisto->getBackGuess(bckVal, sigmaVal);
168 
169  // release the memory
170  delete theHisto;
171 
172  return;
173 }
174 
175 void BackgroundCell::getStats(const PIXTYPE* cellData, const size_t& ndata, double& mean, double& sigma, size_t& statNData)
176 {
177  //size_t npix=0;
178  PIXTYPE lcut=0;
179  PIXTYPE hcut=0;
180  double pixVal=0;
181 
182  // go over the data and do the sum ups
183  mean=0.0;
184  sigma=0.0;
185  for (size_t index=0; index<ndata; index++)
186  {
187  pixVal =cellData[index];
188  if (pixVal>-BIG)
189  {
190  mean += pixVal;
191  sigma += pixVal*pixVal;
192  statNData +=1;
193  }
194  }
195 
196  // check whether there is enough data to continue
197  if ((float)statNData < (float)(ndata*BACK_MINGOODFRAC))
198  {
199  mean = -BIG;
200  sigma = -BIG;
201  return;
202  }
203 
204  // compute mean and sigma of all data
205  if (statNData == 0) {
206  throw Elements::Exception() << "Can not compute meaningful stats with statNData=" << statNData << "!";
207  }
208  mean /= (double)statNData;
209  sigma = sigma/statNData - mean*mean;
210  sigma = sigma > std::numeric_limits<double>::epsilon() ? sqrt(sigma):0.0;
211 
212  // define the range for the second round
213  lcut = (PIXTYPE)(mean - BACK_FKAPPA*sigma);
214  hcut = (PIXTYPE)(mean + BACK_FKAPPA*sigma);
215 
216  // go over the data, apply
217  // the cuts and do the sum ups
218  mean = 0.0;
219  sigma = 0.0;
220  statNData = 0;
221  for (size_t index=0; index<ndata; index++)
222  {
223  pixVal =cellData[index];
224  if (pixVal>=lcut && pixVal <= hcut)
225  {
226  mean += pixVal;
227  sigma += pixVal*pixVal;
228  statNData +=1;
229  }
230  }
231 
232  // compute mean and sigma
233  // in the restricted range
234  if (statNData == 0) {
235  throw Elements::Exception() << "Can not compute meaningful stats with statNData=" << statNData << "!";
236  }
237 
238  mean /= (double)statNData;
239  sigma = sigma/statNData - mean*mean;
240  sigma = sigma>0.0 ? sqrt(sigma):0.0;
241 
242  return;
243 }
244 
245 void BackgroundCell::getStatsWeight(const PIXTYPE* cellData, const size_t& ndata, const PIXTYPE* cellWeight, const PIXTYPE weightThresh, double& mean, double& sigma, size_t& statNData, double& weightMean, double& weightSigma, size_t& statNWeight)
246 {
247  double lcut=0;
248  double hcut=0;
249  double weightLcut=0;
250  double weightHcut=0;
251  double pixVal=0;
252  double weightVal=0;
253  int weightSigma_zero=0;
254 
255  // go over the data and make the sums
256  // for mean and sigma
257  weightMean = 0.0;
258  weightSigma = 0.0;
259  mean=0.0;
260  sigma=0.0;
261  for (size_t index=0; index<ndata; index++)
262  {
263  pixVal = static_cast<double>(cellData[index]);
264  weightVal = static_cast<double>(cellWeight[index]);
265  if (weightVal < static_cast<double>(weightThresh) && pixVal > -BIG)
266  {
267  mean += pixVal;
268  sigma += pixVal*pixVal;
269  weightMean += weightVal;
270  weightSigma += weightVal*weightVal;
271  statNData++;
272 
273  }
274  }
275 
276  // check whether there is enough data to continue
277  if (static_cast<float>(statNData) < static_cast<float>(ndata*BACK_MINGOODFRAC))
278  {
279  mean = -BIG;
280  sigma = -BIG;
281  weightMean = -BIG;
282  weightSigma = -BIG;
283  return;
284  }
285 
286  // compute mean and sigma of all data
287  mean /= static_cast<double>(statNData);
288  sigma = sigma/static_cast<double>(statNData) - mean*mean;
289  sigma = sigma > std::numeric_limits<double>::epsilon() ? sqrt(sigma):0.0;
290  weightMean /= static_cast<double>(statNData);
291  weightSigma = weightSigma/static_cast<double>(statNData) - weightMean*weightMean;
292  weightSigma = weightSigma > std::numeric_limits<double>::epsilon() ? sqrt(weightSigma) : 0.0;
293  if (weightSigma==0.0 && statNData>0)
294  weightSigma_zero=1;
295 
296  // define the ranges for the second round
297  lcut = mean - BACK_FKAPPA*sigma;
298  hcut = mean + BACK_FKAPPA*sigma;
299  weightLcut = weightMean - BACK_FKAPPA*weightSigma;
300  weightHcut = weightMean + BACK_FKAPPA*weightSigma;
301 
302  // go over the data and make the sums
303  // for mean and sigma
304  weightMean = 0.0;
305  weightSigma = 0.0;
306  mean=0.0;
307  sigma=0.0;
308  statNData=0;
309  statNWeight=0;
310  for (size_t index=0; index<ndata; index++)
311  {
312  pixVal = static_cast<double>(cellData[index]);
313  weightVal = static_cast<double>(cellWeight[index]);
314  if (weightVal < weightThresh && pixVal <=hcut && pixVal>=lcut)
315  {
316  mean += pixVal;
317  sigma += pixVal*pixVal;
318  statNData++;
319 
320  if ((weightVal<=weightHcut && weightVal>=weightLcut) || weightSigma_zero)
321  {
322  weightMean += weightVal;
323  weightSigma += weightVal*weightVal;
324  statNWeight++;
325  }
326  }
327  }
328 
329  // compute mean and sigma of the data
330  // within the cuts
331  if (statNData == 0) {
332  throw Elements::Exception() << "Can not compute meaningful data stats with statNData=" << statNData;
333  }
334  mean /= static_cast<double>(statNData);
335  sigma = sigma/statNData - mean*mean;
336  sigma = sigma>0.0 ? sqrt(sigma):0.0;
337 
338  // compute the mean and sigma of the weights
339  // within the cuts
340  if (statNWeight == 0) {
341  throw Elements::Exception() << "Can not compute meaningful weight stats with statNWeight=" << statNWeight;
342  }
343  weightMean /= static_cast<double>(statNWeight);
344  weightSigma = weightSigma/static_cast<double>(statNWeight) - weightMean*weightMean;
345  weightSigma = weightSigma>0.0 ? sqrt(weightSigma):0.0;
346 }
347 } // end of namespace SourceXtractor
#define BIG
#define BACK_MINGOODFRAC
#define BACK_FKAPPA
void getStats(const PIXTYPE *cellData, const size_t &ndata, double &mean, double &sigma, size_t &statNData)
BackgroundHistogram * itsHisto
void getStatsWeight(const PIXTYPE *cellData, const size_t &ndata, const PIXTYPE *cellWeight, const PIXTYPE weightThresh, double &mean, double &sigma, size_t &statNData, double &weightMean, double &weightSigma, size_t &statNWeight)
void getBackgroundValuesOld(PIXTYPE &meanVal, PIXTYPE &sigmaVal)
BackgroundHistogram * itsWeightHisto
void getBackgroundValues(PIXTYPE &meanVal, PIXTYPE &sigmaVal)
BackgroundCell(const PIXTYPE *cellData, const size_t ndata, const PIXTYPE *cellWeight=NULL, const PIXTYPE weightThresh=BIG)
void getBackGuess(PIXTYPE &bckVal, PIXTYPE &sigmaVal)
void addDatum(const PIXTYPE &pixVal)
T epsilon(T... args)
T sqrt(T... args)