21 max_bin_valid =
false;
30 SetLimits(xmin, xmax, ymin, ymax, Nbinsx, Nbinsy);
48 if(this->Nbinsx*this->Nbinsy != Nbinsx*Nbinsy){
55 if(!
hist)
hist =
new double[Nbinsx*Nbinsy];
61 this->Nbinsx = Nbinsx;
62 this->Nbinsy = Nbinsy;
64 bin_widthx = (xmax-xmin)/(
double)(Nbinsx-1);
65 bin_widthy = (ymax-
ymin)/(
double)(Nbinsy-1);
66 bin_size = bin_widthx<bin_widthy ? bin_widthx:bin_widthy;
76 for(
unsigned int i=0; i<Nbinsx*Nbinsy; i++)
hist[i]=0.0;
81 max_bin_content = 0.0;
89 double x = xmin + (0.5+(double)imax_binx)*bin_widthx;
90 double y =
ymin + (0.5+(double)imax_biny)*bin_widthy;
102 return max_bin_content;
110 return bin_widthx/
sqrt(12.0);
118 return bin_widthy/
sqrt(12.0);
142 double small_step = bin_size*1.0E-3;
144 for(
unsigned int i=0; i<points.size(); i++){
152 double beta = FindBeta(xmin,
ymin, xmax-xmin,
ymax-
ymin, pos, g);
161 FindIndexes(pos + small_step*g, ix, iy);
162 if(ix<0 || ix>=(
int)Nbinsx-1 || iy<0 || iy>=(
int)Nbinsy-1){
164 FindIndexes(pos + small_step*g, ix, iy);
167 if(ix<0 || ix>=(
int)Nbinsx-1 || iy<0 || iy>=(
int)Nbinsy-1){
172 unsigned int Niterations=0;
176 beta = FindBeta(xmin+(
double)ix*bin_widthx,
ymin+(
double)iy*bin_widthy, bin_widthx, bin_widthy, pos, g);
179 if(beta*beta > (bin_widthx*bin_widthx + bin_widthy*bin_widthy))
break;
182 if(ix<0 || ix>=(
int)Nbinsx || iy<0 || iy>=(
int)Nbinsy)
break;
183 int index = ix + iy*Nbinsx;
185 if(
hist[index]>max_bin_content){
193 FindIndexes(pos + small_step*g, ix, iy);
195 }
while(++Niterations<2*Nbinsx);
198 return GetMaxBinLocation();
214 int ixmin = (int)(((x-xmin) - 4.0*sigmax)/bin_widthx -0.5);
215 int ixmax = (int)(((x-xmin) + 4.0*sigmax)/bin_widthx +0.5);
216 int iymin = (int)(((y-
ymin) - 4.0*sigmay)/bin_widthy -0.5);
217 int iymax = (int)(((y-
ymin) + 4.0*sigmay)/bin_widthy +0.5);
220 if(ixmax>(
int)Nbinsx)ixmax=Nbinsx;
221 if(iymax>(
int)Nbinsy)iymax=Nbinsy;
226 for(
int i=ixmin; i<ixmax; i++){
227 double x_bin = xmin + (0.5+(double)i)*bin_widthx;
228 for(
int j=iymin; j<iymax; j++){
229 double y_bin =
ymin + (0.5+(double)j)*bin_widthy;
230 double k_x = (x - x_bin)/sigmax;
231 double k_y = (y - y_bin)/sigmay;
232 double gauss_x = exp(-k_x*k_x)/sigmax;
233 double gauss_y = exp(-k_y*k_y)/sigmay;
235 int index = i + j*Nbinsx;
261 if(houghs.size()<1)
return DVector2(0.0, 0.0);
263 unsigned int Nbinsx = houghs[0]->Nbinsx;
264 unsigned int Nbinsy = houghs[0]->Nbinsy;
266 unsigned int imax_binx=0, imax_biny=0;
267 double max_bin_content = 0.0;
268 for(
unsigned int i=0; i<Nbinsx; i++){
269 for(
unsigned int j=0; j<Nbinsy; j++){
270 unsigned int index = i + Nbinsx*j;
272 for(
unsigned int k=0; k<houghs.size(); k++){
273 tot += houghs[k]->hist[
index];
275 if(tot>max_bin_content){
276 max_bin_content = tot;
283 double x = houghs[0]->xmin + (0.5+(double)imax_binx)*houghs[0]->bin_widthx;
284 double y = houghs[0]->ymin + (0.5+(double)imax_biny)*houghs[0]->bin_widthy;
294 bool same_configuration =
true;
295 same_configuration &= hough->
Nbinsx==Nbinsx;
296 same_configuration &= hough->
Nbinsy==Nbinsy;
297 same_configuration &= hough->
xmin==xmin;
298 same_configuration &= hough->
ymin==
ymin;
299 same_configuration &= hough->
xmax==xmax;
300 same_configuration &= hough->
ymax==
ymax;
302 if(!same_configuration){
303 _DBG_<<
"WARNING: trying to add 2 DHoughFind objects with different dimensions!"<<endl;
305 <<
" Nbinsy="<<Nbinsy<<
":"<<hough->
Nbinsy
306 <<
" xmin="<<xmin<<
":"<<hough->
xmin
308 <<
" xmax="<<xmax<<
":"<<hough->
xmax
313 for(
unsigned int i=0; i<Nbinsx*Nbinsy; i++){
323 points.push_back(point);
340 for(
unsigned int i=0; i<points.size(); i++)this->points.push_back(points[i]);
362 string space(Nbinsx-10,
' ');
363 cout<<xmin<<space<<xmax<<endl;
366 for(
unsigned int i=0; i<Nbinsy; i++){
367 string row(Nbinsx,
' ');
368 for(
unsigned int j=0; j<Nbinsx; j++){
370 int index = j + i*Nbinsx;
371 unsigned int d = (
unsigned int)(floor(15.0*
hist[index]/max_bin_content));
374 row[j] = c[0] ==
'0' ?
'.':c[0];
376 cout<<row<<
" "<<
ymin+(double)i*bin_widthy<<endl;
388 TH2D *
h =
new TH2D(hname.c_str(),
"Copied from DHoughFind object",Nbinsx, xmin, xmax, Nbinsy,
ymin,
ymax);
389 for(
unsigned int i=0; i<Nbinsx; i++){
390 for(
unsigned int j=0; j<Nbinsy; j++){
391 h->SetBinContent(i, j,
hist[i+j*Nbinsx]);
void Add(const DHoughFind *hough)
DVector2 GetMaxBinLocation(void)
void AddPoints(const vector< DVector2 > &points)
sprintf(text,"Post KinFit Cut")
void AddPoint(const DVector2 &point)
static char index(char c)
TH2D * MakeIntoRootHist(string hname)
void SetLimits(double xmin, double xmax, double ymin, double ymax, unsigned int Nbinsx, unsigned int Nbinsy)
double GetMaxBinContent(void)
void Fill(double x, double sigmax, double y, double sigmay)