OpeTool.cpp

#include "OpeTool.h"
#include "SctData/RawScanResult.h"
#include "Sct/LogicErrors.h"
#include <boost/shared_ptr.hpp>

#include "TH2.h"
#include "TH1.h"
#include "TMath.h"
#include "TFile.h"

using namespace SctData;
using namespace Sct;
using boost::shared_ptr;

namespace SctAnalysis {

  OpeTool::OpeTool(){}

  OpeTool::~OpeTool(){}

  boost::shared_ptr<SctData::OpeResult> OpeTool::analyzeModule(const RawScanResult& raw){

    unsigned nOPE = raw.nOccupancyPerEvent();

    shared_ptr<OpeResult> result;
    if (nOPE==0) return result;
    
    result = shared_ptr<OpeResult>(new OpeResult(raw.getUniqueID()));

    //For each chip... 
    for (unsigned int ichip=0; ichip<nOPE; ++ichip) {
      const TH2D& chipHist = raw.getOccupancyPerEvent(ichip);
      const int biny = chipHist.GetNbinsY();

      std::ostringstream newhistname;
      newhistname << raw.getClassName() << "." << raw.getUniqueID() 
          << "_Chip" << ichip << "_Ope_Badness"; 

      std::ostringstream newtitle;
      newtitle << "Chip" << ichip << "_Ope_Badness"; 

      const char* c_name=newhistname.str().c_str();
      const char* c_title=newtitle.str().c_str();
      //  we need to make a non-const number of bins :(
      int NonConstBinY=biny;
      shared_ptr<TH1> badnessHist(new TH1D(c_name, c_title, NonConstBinY, chipHist.GetYaxis()->GetXbins()->GetArray()));

      // We want only bad OPEs over a threshold of -20. Below -20, OPE is dominated by s-curve wiggles...
      // The modules will run at around 1fC which corresponds to a threshold of 0 so ignore anything below -20. 
      // Lets find which bin -20 corresponds to. 
      double maxabove20 = 0;
      int badbin = 0;
      // loop over y bins making projections
      for(int iy=1; iy<biny+1; ++iy){
    std::ostringstream projectionName;
    projectionName << raw.getClassName() << "." << raw.getUniqueID() 
               << "_Chip" << ichip << "_bin" << iy; 
    shared_ptr<TH1D> proj(chipHist.ProjectionX(projectionName.str().c_str(),iy,iy,""));
    double mybadness = fractionVariance(*proj);
    badnessHist->SetBinContent(iy, mybadness);

    if ((badnessHist->GetBinCenter(iy))>-20.0) {
      if (mybadness > maxabove20) {
        maxabove20 = mybadness;
        badbin = iy;
      }
    }
      }

      SctData::ChipOpeResult& chipResult=result->getChipResult(ichip);

      // FILL IN THAT HISTOGRAM
      chipResult.hist=badnessHist;

      chipResult.badmax= maxabove20; 
      chipResult.threshold= badnessHist->GetBinCenter(badbin);

    }

    // Return result. The defect will be decided upon in the NoiseOccupancyAlgorithm. 
    return result;
  }

  double OpeTool::fractionVariance(const TH1& data) {

    double range = data.GetXaxis()->GetXmax() - data.GetXaxis()->GetXmin();
    double prob = (data.GetMean())/range;
    // theoretically calculate the RMS from the probability dist:
    double calcRMS = sqrt(range*prob*(1-prob));
    // measure the RMS from data
    double measRMS = data.GetRMS();
    // return the fraction of the measured and calculated RMS
    // std::cout << "OPE: measuredRMS :" << measRMS << " calcRMS :" << calcRMS << std::endl;    
    if (calcRMS ==0) return 1;
    else return measRMS/calcRMS;
  }
  
}

---