Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

MedianFit - source file

// 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;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31