RxThresholdAlgorithm.cpp

#include "RxThresholdAlgorithm.h"
#include "AnalysisAlgorithmMap.h"
#include "SctData/RxThresholdTestResult.h"
#include "SctData/RawScanResult.h"
#include "SctData/FitObject.h"
#include "SctData/TopHatFitObject.h"
#include "SctData/DefectList.h"
#include "SctData/StandardDefects.h"
#include <string>
#include "TH2.h"

using namespace std;
using namespace boost;
using namespace SctData;
using namespace Sct;

namespace SctAnalysis {
    
bool RxThresholdAlgorithm::inMap = AnalysisAlgorithmMap::instance().setAlgorithm("RxThresholdTest", auto_ptr<AnalysisAlgorithm>(new RxThresholdAlgorithm()));

shared_ptr<AnalysisAlgorithm> RxThresholdAlgorithm::clone(const TestData& testData, const string& moduleName) const throw(){
    return shared_ptr<AnalysisAlgorithm>(new RxThresholdAlgorithm(testData, moduleName, *this));
}

shared_ptr<TestResult> RxThresholdAlgorithm::createTestResult() const {
    return shared_ptr<RxThresholdTestResult> (new RxThresholdTestResult()) ;
}

void RxThresholdAlgorithm::loadData() {
    loadAllRaws();
}

bool RxThresholdAlgorithm::canAnalyze() const {
    return hasAllRaws();
}

void RxThresholdAlgorithm::analyze() {
    RxThresholdTestResult& result = *dynamic_pointer_cast<RxThresholdTestResult>(getTestResult());
    const RawScanResult& raw = *getRaw(0);
    result.setScanVariable(raw.getHeader().getVariable());  
    
    const SctData::ScanPoints& points=raw.getPoints();
    bool ascending = points.ascending();
    bool pass = true;
    result.setNOptima(nLinkModule);
    for (short unsigned int ilink=0; ilink<nLinkModule; ++ilink ) {
      const TH2& data = raw.getScanData(ilink);
      const unsigned nx = data.GetNbinsX();
      const unsigned ny = data.GetNbinsY();
      double* projection = new double[ny];
      int start=-1;
      int stop =-1;
      for (unsigned iy=0; iy<ny; ++iy){
    const unsigned ipt = ascending ? iy : ny-iy-1;
    projection[ipt]=0;
    for (unsigned ix=0; ix<nx; ++ix){
      projection[iy]+=data.GetBinContent(ix+1,ipt+1);
    }
        double compare = points.getNEvents(ipt)*nx*0.75;
    if (ipt>0 && projection[ipt]<compare && projection[ipt-1]>compare) {
      std::cout << ipt << "=" << points[ipt] << " " 
            << projection[ipt-1] << " -> " << projection[ipt]
                << " : " << compare << std::endl;
      if (start!=-1) pass=false;
      start=ipt;
    }
        compare = points.getNEvents(ipt)*nx*0.25;
    if (ipt>0 && projection[ipt]<compare && projection[ipt-1]>compare) {
      std::cout << ipt << "=" << points[ipt] << " " 
            << projection[ipt-1] << " -> " << projection[ipt]
                << " : " << compare << std::endl;     
      if (stop!=-1) pass=false;
      stop=ipt;
    }
      }
      
      cout << "start = " << start << ", stop = " << stop << endl;

      double optimum=-1;
      const double fraction= 0.5;
      const unsigned zero  = ascending ? 0 : ny-1;
      const unsigned end   = ascending ? ny-1 : 0;
      
      if (start==-1&&stop==-1){
    optimum=points[zero]+fraction*(points[end]-points[zero]);
      } else if (start==-1){
    optimum=points[zero]+fraction*(points[stop]-points[zero]);
      } else if (stop==-1){
    optimum=points[start]+fraction*(points[end]-points[start]);
      }else{
    optimum=points[start]+fraction*(points[stop]-points[start]);
      }
      cout << "optimum=" << optimum << endl;
      result.setOptimum(ilink, optimum);
    }
    result.setPassed(pass);
}


}

---