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

bool FitAlgorithm::debug() const{
  return Fitter::instance().debug();
}

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);
     
    auto_ptr<FitObject> o (getPrototype());
    if (mode.fittingChannels() ) {
        fitted->initializeChannelFits(*o);
        for (unsigned int ichannel=0; ichannel<fitted->getNChannelFits() ; ++ichannel) {
        if (debug()) {
          std::cout << "Fitting channel " << ichannel << std::endl << std::endl;
        }
        std::ostringstream histname;
        histname << raw.getHeader().getUniqueID() << ".chan." << ichannel;
            scoped_ptr<TH1> projection (occ.getOccupancy(histname.str().c_str(),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()) {
          std::cout << "Fitting chip " << ichip << std::endl << std::endl;
          std::cout << "Existing defects: \n" << *fitted->getDefects().getDefectsAffectingElement(ModuleElement::Chip(ichip)) << std::endl;
        }
        std::ostringstream histname;
        histname << raw.getHeader().getUniqueID() << ".chip." << ichip;
            scoped_ptr<TH1> projection (occ.getOccupancy(histname.str().c_str(),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()) std::cout << "Fitting link " << ilink << std::endl << std::endl;
        std::ostringstream histname;
        histname << raw.getHeader().getUniqueID() << ".lk." << ilink;
            scoped_ptr<TH1> projection (occ.getOccupancy(histname.str().c_str(),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();
    
    if (debug()) std::cout << "Doing fit on [" << hist.GetName() << "]" << std::endl;
    
    // check the projection for defects
    checkForDefects(hist, element, defects);
    if ( defects.defectSeverityEncompassingElement(element) >= SERIOUS) {
      if (debug()) {
    std::cout << "Skipping because of serious defect\n" << defects << std::endl;
      }
      return;
    }
    // try to guess the initial parameters of the fit
    try {
        guessParameters(hist, fitObject);
    if (debug()) std::cout << "Fit object after guess parameters :\n" << fitObject << std::endl;
    } catch ( MathsError& e ) {
        fitter.incrementFitErrors();
    if ( debug()) std::cout << e.what();
    }

    // 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.internal_getFitStrategy().fitTH1(hist, *fit);
    } catch ( MathsError& e){
    fitter.incrementFitErrors();
    if ( debug()) std::cout << e.what();
    }

    if (debug()) std::cout << "Fit object after fit :\n" << fitObject << std::endl;;
    
    fitter.incrementFitsDone();
    fitObject = *fit;                     // make the fitObject equal to the TF1.
    checkForDefects(fitObject, element, defects);

    if (debug()) std::cout << defects;
}

//Utility methods

int FitAlgorithm::findBin(const TH1& h, const float fraction, const bool forward, const bool relative) throw(MathsError) {
    int bin=-1;                         // to return: the bin we are interested in
    float y = fraction;  // the level we want to match to
    if (relative) y*= h.GetMaximum();
    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));
}




}

---