---

FitAlgorithm.cpp

Go to the documentation of this file.

#include "FitAlgorithm.h"
#include "../src/FitStrategy.h"
#include "Fitter.h"
#include <boost/scoped_ptr.hpp>
#include "SctData/RawScanResult.h"
#include "SctData/FitScanResult.h"
#include "SctData/FitObject.h"
#include "SctData/ModuleDefect.h"
#include "SctData/OccupancyProjector.h"
#include <TH1.h>

using namespace SctData;
using boost::scoped_ptr;

namespace SctFitter {

FitAlgorithm::FitAlgorithm() throw() {}

const FitterMode& FitAlgorithm::getMode() const throw() {
    return mode;
}

FitterMode& FitAlgorithm::getMode() throw() {
    return mode;
}

void FitAlgorithm::setMode(const FitterMode& mode) throw() {
    this->mode = mode;
}

//Main algorithm
auto_ptr<FitScanResult> FitAlgorithm::doFit(const RawScanResult& raw) const throw(Sct::LogicError) {
    auto_ptr<FitScanResult> fitted (new FitScanResult(raw));
    OccupancyProjector occ(raw);
     
    auto_ptr<FitObject> o (getPrototype());
    if (mode.fittingChannels() ) {
        fitted->initializeChannelFits(*o);
        for (unsigned int ichannel=0; ichannel<fitted->getNChannelFits() ; ++ichannel) {
            scoped_ptr<TH1> projection (occ.getOccupancy("name", Channel(ichannel)));
            FitObject& fitObject = fitted->getChannelFit(ichannel);
        this->doFit(*projection, fitObject, Channel(ichannel),fitted->getDefects());
        }
    }
    occ.vetoSevereDefects(fitted->getDefects());

    if ( mode.fittingChips() ) {
        fitted->initializeChipFits(*o);
        for (unsigned int ichip=0; ichip<fitted->getNChipFits() ; ++ichip) {
            scoped_ptr<TH1> projection (occ.getOccupancy("name", Chip(ichip)));
            FitObject& fitObject = fitted->getChipFit(ichip);
        this->doFit(*projection, fitObject, Chip(ichip), fitted->getDefects());
        }
    }
    occ.vetoSevereDefects(fitted->getDefects());
    if ( mode.fittingLinks() ) {
        fitted->initializeLinkFits(*o);
        for (unsigned int ilink=0; ilink<fitted->getNLinkFits(); ++ilink) {
            scoped_ptr<TH1> projection (occ.getOccupancy("name", Link(ilink)));
            FitObject& fitObject = fitted->getLinkFit(ilink);
        this->doFit(*projection, fitObject, Link(ilink), fitted->getDefects());
        }
    }

    if ( mode.doingSummary() )
        createSummaryHistograms(*fitted);

    return fitted;
}


void FitAlgorithm::doFit(const TH1& hist, FitObject& fitObject, 
             const ModuleElement& element, ModuleDefectList& defects) const throw (LogicError) {
    Fitter& fitter = Fitter::instance();
    
    // check the projection for defects
    checkForDefects(hist, element, defects);
    if ( defects.severeDefectEncompassingElement(element) ) return;
    // try to guess the initial parameters of the fit
    try {
        guessParameters(hist, fitObject);
    } catch ( MathsError& e ) {
        fitter.incrementFitErrors();
    }

    // make a root TF1 of the apropriate type from that FitObject
    scoped_ptr<TF1> fit (fitObject.makeRootTF1());

    // get the strategy to do the fitting
    try {
    fitter.getFitStrategy().fitTH1(hist, *fit);
    } catch ( MathsError& e){
    fitter.incrementFitErrors();
    }

    fitter.incrementFitsDone();
    fitObject = *fit;                     // make the fitObject equal to the TF1.
}

//Utility methods

int FitAlgorithm::findBin(const TH1& h, const float fraction, const bool forward) throw(MathsError) {
    int bin=-1;                         // to return: the bin we are interested in
    float y = h.GetMaximum()*fraction;  // the level we want to match to
    if (forward) {                          // find last bin which is > y
        for(int i=1; i<=h.GetNbinsX(); i++) {
            if(h.GetBinContent(i)>y) {
                bin = i;
                break;
            }
        }
    } else {                            // find first bin which is > y
        for(int i=h.GetNbinsX(); i>0; i--) {
            if(h.GetBinContent(i)>y) {
                bin = i;
                break;
            }
        }
    }
    if (bin==-1) {
        throw MathsError("FitAlgorithm::findBin no valid bin", __FILE__, __LINE__);
    }
    return bin;
}

float FitAlgorithm::findLevel(const TH1& h, float fraction, bool forward) throw(MathsError) {
    int bin =  FitAlgorithm::findBin(h, fraction, forward);
    // We have found the bin, now linearly interpolate for value of x at fraction

    float x1 = h.GetBinCenter(bin-1);
    float x2 = h.GetBinCenter(bin+1);
    float y1 = h.GetBinContent(bin-1);
    float y2 = h.GetBinContent(bin+1);
    float y =  h.GetMaximum()*fraction;

    float dy = y2-y1;
    if(dy==0)
        dy=1; // protection against divide by 0

    return x1 + ((y - y1)*(x2-x1)/(dy));
}




}

---