00001 #include "ThresholdFitAlgorithm2.h"
00002 #include "FitAlgorithmMap.h"
00003 #include "FitterMode.h"
00004 #include "SctData/mVThresholdVariable.h"
00005 #include "SctData/mVfromTrimTargetThresholdVariable.h"
00006 #include "SctData/ErfcFitObject.h"
00007 #include "SctData/StandardDefects.h"
00008 #include <TH1.h>
00009
00010 using namespace SctData;
00011
00012 namespace SctFitter {
00013
00014
00015 ThresholdFitAlgorithm2::ThresholdFitAlgorithm2() throw() {
00016 FitterMode mode = getMode();
00017 mode.fitNone();
00018 mode.fitChips();
00019 mode.fitChannels() ;
00020 mode.doSummary(false);
00021 setMode(mode);
00022 }
00023
00024 bool ThresholdFitAlgorithm2::putInMap() throw() {
00025 return FitAlgorithmMap::instance().setAlgorithm("ThresholdAnalytic", auto_ptr<FitAlgorithm>(new ThresholdFitAlgorithm2));
00026 }
00027
00028 bool ThresholdFitAlgorithm2::inMap = putInMap();
00029
00030 void ThresholdFitAlgorithm2::doFit(const TH1& hist, FitObject& fitObject,
00031 const ModuleElement& element, DefectList& defects) const throw (LogicError) {
00032
00033 double amp = hist.GetBinContent(1);
00034 unsigned int ampCount = 1;
00035 double mean = 0;
00036 double meanx2 = 0;
00037 double norm = 0;
00038 double last = hist.GetBinContent(1);
00039 bool flatTop = true;
00040
00041 for (unsigned int i=2; i<=hist.GetNbinsX(); ++i) {
00042 double bin = hist.GetBinContent(i);
00043 double bincentre = hist.GetBinLowEdge(i);
00044 double error = hist.GetBinError(i);
00045 double diff = bin-last;
00046 norm += diff;
00047 if (flatTop && error < (amp/ampCount - bin)) {
00048
00049 flatTop = false;
00050 }
00051 if (flatTop) {
00052 amp += bin;
00053 ++ampCount;
00054 }
00055 mean += diff * bincentre;
00056 meanx2 += diff*bincentre* bincentre;
00057 last = bin;
00058 }
00059
00060
00061
00062 amp /= ampCount;
00063
00064 mean /= norm;
00065
00066 double sigma = sqrt(meanx2/norm - mean * mean);
00067
00068
00069
00070
00071 fitObject.setParameter(0, amp);
00072 fitObject.setParameter(1, mean);
00073 fitObject.setParameter(2, sigma);
00074 fitObject.setParError(0, 0);
00075 fitObject.setParError(1, 0);
00076 fitObject.setParError(2, 0);
00077 }
00078
00079
00080 auto_ptr<FitObject> ThresholdFitAlgorithm2::getPrototype() const throw() {
00081 return auto_ptr<FitObject> (new ErfcFitObject());
00082 }
00083
00084 void ThresholdFitAlgorithm2::guessParameters(const TH1& hist, FitObject& fitOb) const throw (LogicError, MathsError) {
00085 }
00086
00087
00088 void ThresholdFitAlgorithm2::checkForDefects(const FitObject& fo, const ModuleElement& element,
00089 DefectList& defects) const {
00090 if (fo.getNDF() && fo.getChiSquared()/fo.getNDF() > StandardDefects::BADFIT.getParameter() ) {
00091 cout << "Chi-squared: " << fo.getChiSquared() << " NDF: " << fo.getNDF() << endl;
00092 defects.addDefect(Defect(StandardDefects::BADFIT,element));
00093 }
00094
00095 }
00096
00097 void ThresholdFitAlgorithm2::checkForDefects(const TH1& hist, const ModuleElement& element,
00098 DefectList& defects) const throw (LogicError) {
00099 }
00100
00101 void ThresholdFitAlgorithm2::createSummaryHistograms(FitScanResult& fits) const throw() {
00102 }
00103 }