ThresholdFitAlgorithm.cpp

#include "ThresholdFitAlgorithm.h"
#include "FitAlgorithmMap.h"
#include "FitterMode.h"
#include "SctData/mVThresholdVariable.h"
#include "SctData/mVfromTrimTargetThresholdVariable.h"
#include "SctData/mVThresholdVariable.h"

#include "SctData/StandardDefects.h"
#include "SctData/Defect.h"
#include "SctData/ErfcFitObject.h"
#include <TH1.h>

using namespace SctData;

namespace SctFitter {
    
bool ThresholdFitAlgorithm::putInMap() throw() {
    bool first = FitAlgorithmMap::instance().setAlgorithm(mVThresholdVariable::instance().getVariableName(), auto_ptr<FitAlgorithm>(new ThresholdFitAlgorithm));
    bool second = FitAlgorithmMap::instance().setAlgorithm(mVfromTrimTargetThresholdVariable::instance().getVariableName(), auto_ptr<FitAlgorithm>(new ThresholdFitAlgorithm));
    return (first && second);
}
    
bool ThresholdFitAlgorithm::inMap = putInMap();

ThresholdFitAlgorithm::ThresholdFitAlgorithm() throw() {
    FitterMode mode = getMode();
    mode.fitNone(); 
    mode.fitChips(); 
    mode.fitChannels() ; 
    mode.doSummary(false); 
    setMode(mode);
}

auto_ptr<FitObject> ThresholdFitAlgorithm::getPrototype() const throw() {
    return auto_ptr<FitObject> (new ErfcFitObject());
}

void ThresholdFitAlgorithm::guessParameters(const TH1& hist, FitObject& fitOb) const throw (LogicError, MathsError) {
    fitOb.setParameter(0, hist.GetMaximum() ) ;
    fitOb.fixParameter(0, true);
    fitOb.setParameter(1, findLevel(hist,0.5,false)) ;                 // 50% level
    double p2 = findLevel(hist, 0.82, false) - fitOb.getParameter(1);    // level 1-sigma away from 50%
    if (p2<0.) p2=-p2;  
    fitOb.setParameter(2, p2 ) ;

    int lbin=findBin(hist, 0.99, false); if ((lbin-1)>=1) lbin--;
    int ubin=findBin(hist, 0.01, false); if ((ubin+1)<hist.GetNbinsX()) ubin++;
    fitOb.setVarMin(0 , hist.GetBinLowEdge( lbin ) );    // Lower edge of bin before first bin < 99%
    fitOb.setVarMax(0 , hist.GetBinLowEdge( ubin+1 ) );    // Upper edge of bin after first bin <1% 
}

void ThresholdFitAlgorithm::checkForDefects(const FitObject& fo, const ModuleElement& element, DefectList& defects) const {
  if (fo.getNDF() &&  fo.getChiSquared()/fo.getNDF() > StandardDefects::BADFIT.getParameter() ) {
      cout << "Chi-squared: " << fo.getChiSquared() << " NDF: " << fo.getNDF() << " so claiming is DEAD" << endl;
      defects.addDefect(Defect(StandardDefects::BADFIT,element));
  }
}

void ThresholdFitAlgorithm::checkForDefects(const TH1& hist, const ModuleElement& element, DefectList& defects) const throw (LogicError) {
    double maxVal = hist.GetMaximum();
    double minVal = hist.GetMinimum();
    
    //Prevent ROOT calculating this millions of times!
    TH1* h_nonconst = const_cast<TH1*>(&hist);
    h_nonconst->SetMaximum(maxVal);
    h_nonconst->SetMinimum(minVal);
    
    //Don't check for defects where the whole chip can be affected by a single channel,
    //which would already be marked as defective.  Examples are UNDER and OVER, where a single
    //UNDER or OVER channel would cause that chip to be UNDER or OVER
    bool minorDefectCheck = !(element.isChip() && getMode().fittingChannels());

    if ( maxVal < StandardDefects::DEAD.getParameter() ){
    defects.addDefect(Defect(StandardDefects::DEAD, element));
    } else if(minorDefectCheck && maxVal < StandardDefects::UNDER.getParameter() ){
    defects.addDefect(Defect(StandardDefects::UNDER,element));
    } else if(minorDefectCheck && maxVal > StandardDefects::OVER.getParameter() ){
    defects.addDefect(Defect(StandardDefects::OVER,element));
    }
    if( minVal > StandardDefects::STUCKON.getParameter() ){
    defects.addDefect(Defect(StandardDefects::STUCKON, element));
    } 
    
}

void ThresholdFitAlgorithm::createSummaryHistograms(FitScanResult& fits) const throw() {
    /*
    string name = getUniqueID() + "Mean";
    string title = getClassName() + "." + getUniqueID() + " Mean";
    TH1D* rooth = new TH1D(name.c_str(), title.c_str(), 30, 0, 20);
    if (!rooth) throw InvalidArgument("FitScanResult::createSummaryHistograms() couldn't make histogram");
    for (vector<FitObject*>::iterator is = m_channelFits.begin(); is!= m_channelFits.end(); ++is ) {
    rooth->Fill( (*is)->getParameter(1) );
    }
*/
}
    
}

---