FitAlgorithm.cpp

#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/OccupancyProjector.h"
#include <TH1.h>
#include <TF1.h>

using namespace SctData;
using boost::scoped_ptr;

namespace SctFitter {

FitAlgorithm::FitAlgorithm() throw() {}

const FitterMode& FitAlgorithm::getMode() const throw() {
  boost::recursive_mutex::scoped_lock lock (getMutex());
  return mode;
}

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

void FitAlgorithm::setMode(const FitterMode& mode) throw() {
  boost::recursive_mutex::scoped_lock lock (getMutex());
  this->mode = mode;
}

//Main algorithm
auto_ptr<FitScanResult> FitAlgorithm::doFit(const RawScanResult& raw) const throw(Sct::LogicError) {
  boost::recursive_mutex::scoped_lock lock (getMutex());
  auto_ptr<FitScanResult> fitted (new FitScanResult(raw));
  OccupancyProjector occ(raw);
     
    bool debug = (Fitter::instance().getFitStrategy().getOptions().find('Q')==string::npos && Fitter::instance().getFitStrategy().getOptions().find('q')==string::npos );
    
    //const ScanPoints& points = raw.getPoints();
    //for (unsigned int i=0; i<points.getNPoints(); ++i)
    //cout << "Point " << i << " NEvents: " << points.getNEvents(i) << endl;
    
    auto_ptr<FitObject> o (getPrototype());
    if (mode.fittingChannels() ) {
        fitted->initializeChannelFits(*o);
        for (unsigned int ichannel=0; ichannel<fitted->getNChannelFits() ; ++ichannel) {
        if (debug) 
        cout << "Fitting channel " << ichannel << endl << endl;
            scoped_ptr<TH1> projection (occ.getOccupancy("name", ModuleElement::Channel(ichannel)));
            FitObject& fitObject = fitted->getChannelFit(ichannel);
        this->doFit(*projection, fitObject, ModuleElement::Channel(ichannel),fitted->getDefects());
        }
    }
    occ.vetoSeriousDefects(fitted->getDefects());
    occ.vetoMaskedChannels(raw.getConfiguration());

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

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

    return fitted;
}


void FitAlgorithm::doFit(const TH1& hist, FitObject& fitObject, 
             const ModuleElement& element, DefectList& defects) const throw (LogicError) {
    Fitter& fitter = Fitter::instance();
    
    // check the projection for defects
    checkForDefects(hist, element, defects);
    if ( defects.defectSeverityEncompassingElement(element) >= SERIOUS) 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.
    checkForDefects(fitObject, element, defects);
}

//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));
}




}

---