ThresholdFitAlgorithm2.cpp

#include "ThresholdFitAlgorithm2.h"
#include "FitAlgorithmMap.h"
#include "FitterMode.h"
#include "SctData/mVThresholdVariable.h"
#include "SctData/mVfromTrimTargetThresholdVariable.h"
#include "SctData/ErfcFitObject.h"
#include "SctData/StandardDefects.h"
#include <TH1.h>

using namespace SctData;

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

bool ThresholdFitAlgorithm2::putInMap() throw() {    
    inMap = FitAlgorithmMap::instance().setAlgorithm(mVThresholdVariable::instance().getVariableName(), auto_ptr<FitAlgorithm>(new ThresholdFitAlgorithm2));
    inMap |= FitAlgorithmMap::instance().setAlgorithm(mVfromTrimTargetThresholdVariable::instance().getVariableName(), auto_ptr<FitAlgorithm>(new ThresholdFitAlgorithm2));
    return inMap;
}
    
bool ThresholdFitAlgorithm2::inMap = false;

void ThresholdFitAlgorithm2::doFit(const TH1& hist, FitObject& fitObject,
                   const ModuleElement& element, DefectList& defects) const throw (LogicError) {
    
    double amp = hist.GetBinContent(1);
    unsigned int ampCount = 1;
    double mean = 0;
    double meanx2 = 0;
    double norm = 0;
    double last = hist.GetBinContent(1);
    bool flatTop = true;
    
    for (unsigned int i=2; i<=hist.GetNbinsX(); ++i) {
    double bin = hist.GetBinContent(i);
    double bincentre = hist.GetBinLowEdge(i);
    double error = hist.GetBinError(i);
    double diff = bin-last;
    norm += diff;
    if (flatTop && error < (amp/ampCount - bin)) {
        //cout << "diff>error: " << amp << " error: " << error << " last: "<< last << " bin: " << i << " bincentre" << bincentre << endl; 
        flatTop = false;
    }
    if (flatTop) {
        amp += bin;
        ++ampCount;
    }
    mean += diff * bincentre;
    meanx2 += diff*bincentre* bincentre;
    last = bin;
    }
 
    //cout << "amp: " << amp << " ampCount: " << ampCount << " mean: " << mean << " meanx2 " << meanx2 << endl;
    
    amp /= ampCount;
    //mean /= 1-hist.GetNbinsX();
    mean /= norm;
    //double sigma = sqrt(meanx2 /(1-hist.GetNbinsX()) - mean * mean);
    double sigma = sqrt(meanx2/norm - mean * mean);
    
    //cout << "amp: " << amp << " mean: " << mean << " sigma: " << sigma << endl;;
    
    
    fitObject.setParameter(0, amp);
    fitObject.setParameter(1, mean);
    fitObject.setParameter(2, sigma);
    fitObject.setParError(0, 0);
    fitObject.setParError(1, 0);
    fitObject.setParError(2, 0);    
}

    
auto_ptr<FitObject> ThresholdFitAlgorithm2::getPrototype() const throw() {
    return auto_ptr<FitObject> (new ErfcFitObject());
}
    
void ThresholdFitAlgorithm2::guessParameters(const TH1& hist, FitObject& fitOb) const throw (LogicError, MathsError) {
}


void ThresholdFitAlgorithm2::checkForDefects(const FitObject& fo, const ModuleElement& element, 
                         DefectList& defects) const {
  if (fo.getNDF() &&  fo.getChiSquared()/fo.getNDF() > 5000 ) {
       defects.addDefect(Defect(StandardDefects::DEAD,element));
  }
     
}

void ThresholdFitAlgorithm2::checkForDefects(const TH1& hist, const ModuleElement& element, 
                         DefectList& defects) const throw (LogicError) {
}

void ThresholdFitAlgorithm2::createSummaryHistograms(FitScanResult& fits) const throw() {
}
}

---