MarkSpaceRatioAlgorithm.cpp

#include "MarkSpaceRatioAlgorithm.h"
#include "AnalysisAlgorithmMap.h"
#include "SctData/MarkSpaceRatioTestResult.h"
#include "SctData/RawScanResult.h"
#include <cfloat>
#include "TH2.h"
#include "Sct/SctNames.h"
#include <set>

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

namespace SctAnalysis {
    
    bool MarkSpaceRatioAlgorithm::inMap = AnalysisAlgorithmMap::instance().setAlgorithm("MarkSpaceRatioTest", auto_ptr<AnalysisAlgorithm>(new MarkSpaceRatioAlgorithm()));
    
    shared_ptr<AnalysisAlgorithm> MarkSpaceRatioAlgorithm::clone(shared_ptr<const TestData> testData, const string& moduleName) const throw() {
    return shared_ptr<AnalysisAlgorithm>(new MarkSpaceRatioAlgorithm(testData, moduleName, *this));
    }

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

    
    bool MarkSpaceRatioAlgorithm::canAnalyze() const{
    return hasAllRaws();
    }
    
    shared_ptr<TestResult> MarkSpaceRatioAlgorithm::createTestResult() const {
    return shared_ptr<MarkSpaceRatioTestResult> (new MarkSpaceRatioTestResult());
    }
  
    void MarkSpaceRatioAlgorithm::analyze() {
      MarkSpaceRatioTestResult & result =  *dynamic_pointer_cast<MarkSpaceRatioTestResult> ( getTestResult() );
      const unsigned nscans=result.getNScans();

      //std::cout << "MarkSpaceRatio::analyze with " << nscans << " scans" << std::endl;

      // check validity
      if (nscans<4 || nscans%2!=0){
    std::ostringstream msg;
    msg << "MarkSpaceRatio should have even number of scans (>=4) here we have " << nscans;
                
    throw Sct::LogicError(msg.str(), __FILE__, __LINE__);
      }
      for (unsigned iscan=0; iscan<nscans/2; ++iscan){
    MarkSpaceRatioTestResult::MsrPoint p;
    boost::shared_ptr<const RawScanResult> raw1 = getRaw(iscan);
    boost::shared_ptr<const RawScanResult> raw2 = getRaw(iscan+nscans/2);

    // link0
    p.occ0=getAverageOccupancy(*raw1,0);
    p.occ0flip=getAverageOccupancy(*raw2,0);

    // link1
    p.occ1=getAverageOccupancy(*raw1,1);
    p.occ1flip=getAverageOccupancy(*raw2,1);

    p.msr=static_cast<short int>(floor(result.getTestPointAt(iscan)+0.5));
    result.addPoint(p);
    //std::cout << "Adding point " << p.msr << " " 
    //    << p.occ0 << "->" << p.occ0flip << "\t" 
    //    << p.occ1 << "->" << p.occ1flip << "\t"
    //    << endl;
      }
      fixOptimum(result);
    }
    
    void MarkSpaceRatioAlgorithm::fixOptimum(MarkSpaceRatioTestResult& r){
      double c[2]={0.,0.};
      double m[2]={0.,0.};
      double opt[2]={0.,0.};
      bool good[2]={false, false};

      for (unsigned ilink=0; ilink<2; ++ilink){
    float xbar=0.;
    float ybar=0.;
    float x2bar=0.;
    float y2bar=0.;
    float xybar=0.;
    
    const unsigned& n=r.getNPoints();
    for (unsigned ipoint=0; ipoint<n; ++ipoint){
      const MarkSpaceRatioTestResult::MsrPoint& p = r.getPoint(ipoint);
      
      const double diff=p.getOccupancy(ilink) - p.getOccupancyFlipped(ilink);
      
      ybar  += diff;
      y2bar += diff*diff;
      xbar  += p.msr;
      x2bar  += p.msr*p.msr;
      xybar  += p.msr*diff;
    }

    const float det=n*x2bar-xbar*xbar;
    if (det>0.){
      c[ilink] = (x2bar*ybar - xbar*xybar) / det;
      m[ilink] = (n*xybar - xbar*ybar) / det;; 
    }
    if ( m[ilink] !=0.) { 
      opt[ilink] = -c[ilink]/m[ilink];
      good[ilink]=true;
    }
      }
      
      // use mean weighted by gradient for optimum:
      double numerator=0.;
      double denominator=0.;
      unsigned ngood=0;
      for (unsigned ilink=0; ilink<2; ++ilink){
    if (good[ilink]){
      ngood++;
      numerator+=opt[ilink];
    }
      }

      if(ngood==0){
    r.setPassed(false);
    return;
      }

      double optimum=numerator/ngood;
      r.setOptimum(optimum);
      if (optimum<-1. || optimum>33.) {
    r.setPassed(false);
      } else if (optimum<3. || optimum > 29){
    r.setProblem(true);
      }else{
        r.setPassed(true);
      }
    }

    double MarkSpaceRatioAlgorithm::getAverageOccupancy(const RawScanResult& raw, unsigned short ilink){
      double contents=0.;
      double trigs=0.;
      const TH2& hist=  raw.getScanData(ilink);
      const ScanPoints& points=raw.getPoints();

      for (unsigned iy=1;iy<=hist.GetNbinsY(); ++iy){
    for (unsigned ix=4;ix<=hist.GetNbinsX(); ++ix){
      contents += hist.GetBinContent(ix,iy);
    }
    unsigned ipt = points.ascending() ? iy-1 : points.getNPoints()-iy;
    trigs += points.getNEvents(ipt)*hist.GetNbinsX();
      }
      if (trigs<0.0001) return 0.; 
      return contents/trigs;
    }

} // end of namespace SctAnalysis

---