// ESAF : Euso Simulation and Analysis Framework
// $Id: MedianFit.cc,v 1.8 2004/06/17 08:20:12 thea Exp $
// R. Pesce created Jan, 29 2004
#include <iostream>
#include "MedianFit.hh"
#include "TMath.h"
ClassImp(MedianFit)
// ctor
MedianFit::MedianFit( Int_t numpoints, vector<Double_t> x, vector<Double_t> y ) {
fNumPoints = numpoints;
fSlope = 0;
fOffset = 0;
fAbsoluteDeviation = 0;
fPointsX = x;
fPointsY = y;
Do();
}
// dtor
MedianFit::~MedianFit() {
fPointsX.clear();
fPointsY.clear();
}
// process of median fit
// adapted from "Numerical Recipes in C" pag.704
void MedianFit::Do() {
Double_t XSum = 0;
Double_t XSqSum = 0;
Double_t YSum = 0;
Double_t XYSum = 0;
for( Int_t i=0; i<fNumPoints; i++ ) {
XSum += fPointsX[i];
XSqSum += fPointsX[i] * fPointsX[i];
YSum += fPointsY[i];
XYSum += fPointsX[i] * fPointsY[i];
}
// method of least squares for a first guess for slope and offset
fOffset = (XSqSum * YSum - XSum * XYSum) / (fNumPoints * XSqSum - XSum * XSum);
fSlope = (fNumPoints * XYSum - XSum * YSum) / (fNumPoints * XSqSum - XSum * XSum);
Double_t chisquare = 0;
for( Int_t i=0; i<fNumPoints; i++ ) {
chisquare += (fPointsY[i] - (fOffset + fSlope * fPointsX[i] )) *
(fPointsY[i] - (fOffset + fSlope * fPointsX[i] )) ;
}
// standard deviation: idea of how big an interation step to take
Double_t sigmaslope = TMath::Sqrt(chisquare / fNumPoints * XSqSum - XSum * XSum);
Double_t b1 = fSlope;
Double_t r1 = RoFunc( b1 );
if ( sigmaslope > 0 ) {
Double_t b2 = fSlope + TMath::Sign(3. * sigmaslope, r1);
// gues bracket as 3-sigma away, in the downhill direction known from r1
Double_t r2 = RoFunc( b2 );
if ( b2 == b1 ) {
fAbsoluteDeviation /= fNumPoints;
return;
}
// bracketing
while( r1*r2 > 0 ) {
fSlope = b2 + 1.6 * ( b2 - b1 );
b1 = b2;
r1 = r2;
b2 = fSlope;
r2 = RoFunc( b2 );
}
sigmaslope *= 0.01;
// refine until error a neglegible number of sigmas
while( TMath::Abs(b2-b1) > sigmaslope ) {
fSlope = b1 + 0.5 * ( b2 - b1 );
// bisection
if ( ( fSlope == b1 ) || ( fSlope == b2 ) ) break;
Double_t r = RoFunc( fSlope );
if ( r*r1 >= 0 ) {
r1 = r;
b1 = fSlope;
} else {
r2 = r;
b2 = fSlope;
}
}
}
fAbsoluteDeviation /= fNumPoints;
}
// adapted from "Numerical Recipes in C pag.705"
Double_t MedianFit::RoFunc( Double_t b ) {
vector<Double_t> dummy;
for( Int_t i=0; i<fNumPoints; i++ ) {
dummy.push_back( fPointsY[i]- b * fPointsX[i] );
}
if ( ( fNumPoints % 2 ) == 0 ) {
fOffset = 0.5 * (Select( fNumPoints/2, dummy ) + Select( (fNumPoints/2) + 1, dummy ));
} else {
fOffset = Select( (fNumPoints+1)/2, dummy);
}
Float_t sum = 0.;
fAbsoluteDeviation = 0.;
for( Int_t i=0; i<fNumPoints; i++ ) {
Double_t d = fPointsY[i] - ( fOffset + b * fPointsX[i] );
fAbsoluteDeviation += TMath::Abs(d);
if ( fPointsY[i] != 0 )
d /= fPointsY[i];
if ( TMath::Abs(d) > 1.0e-7)
sum += ( d >= 0 ? fPointsX[i] : -fPointsX[i] );
}
return sum;
}
// find the k-th smallest element in a vector; used to find the median
// adapted from "Numerical Recipes in C" pag.342
Double_t MedianFit::Select( Int_t k, vector<Double_t> vec ) {
Int_t l = 1;
Int_t ir = vec.size();
for( ; ; ) {
if ( (ir - l) <= 1 ) { // active partition contains 1 or 2 elements
if ( (ir - l) == 1 && vec[ir] < vec[l] ) { // active partition contains 2 elements
Swap( vec[l], vec[ir] );
}
return vec[k];
} else {
// choose median of left, center and right elements as partitioning element a
// rearrange so that vec[l] <= vec[l+1] <= vec[ir]
Int_t middle = ( l + ir ) / 2;
Swap( vec[middle], vec[l+1] );
if ( vec[l] > vec[ir] ) {
Swap( vec[l], vec[ir] );
}
if ( vec[l+1] > vec[ir] ) {
Swap( vec[l+1], vec[ir] );
}
if ( vec[l] > vec[l+1] ) {
Swap( vec[l], vec[l+1] );
}
// initialize pointers for partitioning
Int_t i = l + 1;
Int_t j = ir;
Double_t a = vec[l+1]; // partitioning element
for( ; ; ) {
do { // scan up to find element > a
i++;
} while( vec[i] < a && i<= (Int_t)vec.size() );
do { // scan down to find element < a
j--;
} while( vec[j] > a && j>0 );
if ( j < i ) break; // pointers crossed; partitioning complete
Swap( vec[i], vec[j] );
}
vec[l+1] = vec[j]; // insert partitioning element
vec[j] = a;
if ( j >= k ) ir = j - 1; // keep active the partiction that contains the k-th element
if ( j <= k ) l = i;
}
}
}
// swap two variables
void MedianFit::Swap( Double_t w, Double_t z ) {
Double_t temp = w;
w = z;
z = temp;
}