NoiseOccupancyAlgorithm.cpp

#include "NoiseOccupancyAlgorithm.h"
#include "AnalysisAlgorithmMap.h"
#include "AnalysisService.h"
#include "SctData/Stat.h"
#include "SctData/NoiseOccupancyTestResult.h"
#include "SctData/OccupancyProjector.h"
#include "SctData/FitScanResult.h"
#include "SctData/FitObject.h"
#include "SctData/ResponseCurve.h"
#include "SctData/StandardDefects.h"
#include "SctData/ConfigurationVariable.h"
#include "Sct/UnsupportedOperationError.h"
#include <CommonWithDsp/ABCD/ABCDscans.h>
#include <TH1.h>
#include <TGraphAsymmErrors.h>
#include <string>

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

namespace SctAnalysis {
    
    NoiseOccupancyAlgorithm::~NoiseOccupancyAlgorithm() throw() {}

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

    bool NoiseOccupancyAlgorithm::inMap = AnalysisAlgorithmMap::instance().setAlgorithm("NoiseOccupancyTest", auto_ptr<AnalysisAlgorithm>(new NoiseOccupancyAlgorithm()));

    //----------------------------------------------------------------------
    
    void NoiseOccupancyAlgorithm::loadData() {
    loadAllRaws();
    loadAllFits();
    }
    
    bool NoiseOccupancyAlgorithm::canAnalyze() const {
    return hasAllFits() && hasAllRaws();
    }
    
    shared_ptr<SctData::TestResult> NoiseOccupancyAlgorithm::createTestResult() const {
    return shared_ptr<TestResult>(new NoiseOccupancyTestResult());
    }

    void NoiseOccupancyAlgorithm::analyze() {
    shared_ptr<NoiseOccupancyTestResult> result = dynamic_pointer_cast<NoiseOccupancyTestResult>(getTestResult());
    const SctData::ConfigurationVariable& scanVariable = getRaw(0)->getHeader().getVariable();
    result->setScanVariable(scanVariable);
    const ModuleConfiguration config=getRaw(0)->getConfiguration();
    
    m_scanVariable = SctData::ConfigurationVariableIOHelper::getTypeRep(scanVariable);
    
    //cout << " NoiseOccupancyAlgorithm::analyze()" << endl;
    //Loop over chips
    for (unsigned int ichip=0; ichip<nChipModule; ++ichip) {
        // get an raw result occupany projector 
        OccupancyProjector opr(*getRaw(0) );
        opr.vetoSeriousDefects(getFit(0)->getDefects());
        // get the chip configuration
        const ChipConfiguration& chipconfig=config.getChipConfiguration(ichip);
        // get the new defects
        DefectList& defects = result->getDefects();
        result->getChipResult(ichip) = analyzeChip(ichip, opr, chipconfig, defects);

        const float fittedOffset = getFit(0)->getChipFit(ichip).getParameter(1);
        switch(m_scanVariable) {
        case ST_VTHR:
        result->getChipResult(ichip).offset = fittedOffset;
        break;
        case ST_TTHR:
        result->getChipResult(ichip).offset = fittedOffset + config.getChipConfiguration(ichip).getTrimTarget()*2.5;
        break;
        case ST_QTHR: {
        auto_ptr<ResponseCurve> rc = ResponseCurveMap::getMap().get(config.getChipConfiguration(ichip));
        result->getChipResult(ichip).offset = rc->getFunction()->Eval(fittedOffset);
            }
        break;
        default: {
            ostringstream oss;
        oss << "Internal Error: NoiseOccupancy Algorithm doesn't know how to deal with scan of type : " << m_scanVariable;
        throw Sct::UnsupportedOperationError(oss.str(), __FILE__, __LINE__);
            }
        }
    }
    result->setPassed(true);
    }
    
    ChipNOResult NoiseOccupancyAlgorithm::analyzeChip(const unsigned ichip, const OccupancyProjector& opr, const ChipConfiguration& chipconfig, DefectList& defects) {
    
    auto_ptr<TH1> occ = opr.getOccupancy("temp", ModuleElement::Chip(ichip)) ;
    auto_ptr<ResponseCurve> rc = ResponseCurveMap::getMap().get(chipconfig);
    shared_ptr<TF1> invFn = rc->getInverseFunction();
    
    //cout << "Config RC: " << (int)chipconfig.getRcFunctionIndex() << " params: " << chipconfig.getRcParam(0) << "  " << chipconfig.getRcParam(1) << "  " << chipconfig.getRcParam(2) << endl;
    //cout << "RC: " << rc->getIndex() << " params: " << rc->getFunction()->GetParameter(0) << "  " << rc->getFunction()->GetParameter(1) << "  " << rc->getFunction()->GetParameter(2) << endl;

    int ttbin=-1; //< trim target bin
    int fcbin=-1; //< 1 fC bin

    const float mVstep=2.5;  // number of mv per threshold step
    const float trimTargetmV=chipconfig.getTrimTarget()*mVstep;

    if (m_scanVariable==ST_VTHR){
      ttbin = occ->GetXaxis()->FindBin(chipconfig.getTrimTarget()*mVstep);
      fcbin = occ->GetXaxis()->FindBin(rc->getFunction()->Eval(1.));
    }else if(m_scanVariable==ST_TTHR){
      ttbin = occ->GetXaxis()->FindBin(0.);
      fcbin = occ->GetXaxis()->FindBin(rc->getFunction()->Eval(1.)-trimTargetmV);
    }else if(m_scanVariable==ST_QTHR){
      ttbin = occ->GetXaxis()->FindBin(invFn->Eval(trimTargetmV));
      fcbin = occ->GetXaxis()->FindBin(1.);
    }else{
      ostringstream oss;
      oss << "Internal Error: NoiseOccupancy Algorithm doesn't know how to deal with scan of type : " << m_scanVariable;
      throw Sct::UnsupportedOperationError(oss.str(), __FILE__, __LINE__);
    }

    //cout << "TrimTarget: " << chipconfig.getTrimTarget()*2.5 << " mV, bin: " << ttbin << " Occ: " << occ->GetBinContent(ttbin) << endl;
    //cout << "1fC: " << rc->getFunction()->Eval(1.) << " mV, bin: " << fcbin << " Occ: " << occ->GetBinContent(fcbin) << endl;
    
    // Construct TGraph of log occupancy against thr^2
    // project raw data from scan 0 of the test
    // Construction of TGraphAsymmetricErrors in this way because root is so 
    // inefficient at re-allocating memory for graphs
    double x[occ->GetNbinsX()];
    double xUpperErr[occ->GetNbinsX()];
    double xLowerErr[occ->GetNbinsX()];
    double yUpperErr[occ->GetNbinsX()];
    double yLowerErr[occ->GetNbinsX()];
    double y[occ->GetNbinsX()];
    unsigned int nBins = 0;
    
    double factor = 1;//11.3;  //warning - Should be 1.  Use 11.3 to simulate SCTDAQ
    
    for(int bin=1; bin<=occ->GetNbinsX(); ++bin){
        //We only want the "bottom" half of the threshold scan (occupancy < 50%)
        if (occ->GetBinContent(bin) > 0 && occ->GetBinContent(bin) < 0.5 ) {
        y[nBins] = log(occ->GetBinContent(bin));            
        yUpperErr[nBins] = log(occ->GetBinContent(bin) + factor*occ->GetBinError(bin)) - y[nBins];
        if (occ->GetBinContent(bin) - factor*occ->GetBinError(bin) > 0) {
            yLowerErr[nBins] = y[nBins] - log(occ->GetBinContent(bin) - factor*occ->GetBinError(bin));
        } else {
            yLowerErr[nBins] = yUpperErr[nBins];//occ->GetBinContent(bin);
            //yLowerErr[nBins] = 50;
            //cout << yLowerErr[nBins] << endl;         
        } 
        
        //Convert into fC:
        switch(m_scanVariable) {
        case ST_VTHR:
            x[nBins] = invFn->Eval(occ->GetBinCenter(bin));
            break;
        case ST_TTHR:
            x[nBins] = invFn->Eval(occ->GetBinCenter(bin)+trimTargetmV);
            break;
        case ST_QTHR:
            x[nBins] = occ->GetBinCenter(bin);
            break;
        default:
            ostringstream oss;
            oss << "Internal Error: NoiseOccupancy Algorithm doesn't know how to deal with scan of type : " << m_scanVariable;
            throw Sct::UnsupportedOperationError(oss.str(), __FILE__, __LINE__);
        }
        
        x[nBins] *= x[nBins];
        xUpperErr[nBins] = xLowerErr[nBins] = 0;            
                
        nBins++;
        }
    }

        auto_ptr<TF1> f = NoiseOccupancyTestResult::createFitFunction();

        //Check for a dead chip
        if (nBins == 0) {
            defects.addDefect(Defect(StandardDefects::DEAD, ModuleElement::Chip(ichip)));
            shared_ptr<TGraphAsymmErrors> g( new TGraphAsymmErrors());
            return ChipNOResult(g, shared_ptr<TF1>(f), 0, 0);
        }

    shared_ptr<TGraphAsymmErrors> g( new TGraphAsymmErrors(nBins, x, y, xLowerErr, xUpperErr, yLowerErr, yUpperErr));

    //Don't use FitStrategy because there is no asymm errors method
//  AnalysisService::instance().getFitStrategy().fitTGraphAsymmErrors(*g, *f);
    g->Fit(f.get(), "NQR");

    // Calculate mean and RMS NO at 1fC for channels in this chip
    Stats<double> occAtOnefC(nChannelChip);
    for (unsigned ichan=0; ichan < nChannelChip ; ++ichan ){
        unsigned channelInModule = ichip*nChannelChip + ichan;

        // get chip projection:
        char name[10];
        sprintf(name,"temp%d",channelInModule);
        auto_ptr<TH1> occChannel = opr.getOccupancy(name, ModuleElement::Channel(channelInModule)) ;
        
        // find 1 fC point - gets threshold in mV
        int bin = 0;
        switch(m_scanVariable) {
        case ST_VTHR: {
        double vthr = rc->getFunction()->Eval(1.);
        bin = occChannel->GetXaxis()->FindBin(vthr);        
            } break;
        case ST_TTHR: {
        double vthr = rc->getFunction()->Eval(1.);
        bin = occChannel->GetXaxis()->FindBin(vthr-trimTargetmV);       
            } break;
        case ST_QTHR:
        bin = occChannel->GetXaxis()->FindBin(1.0);
        break;
        default:
        ostringstream oss;
        oss << "Internal Error: NoiseOccupancy Algorithm doesn't know how to deal with scan of type : " << m_scanVariable;
        throw Sct::UnsupportedOperationError(oss.str(), __FILE__, __LINE__);
        }       
        
        //cout << "Thr bin="<<bin<<endl;
        occAtOnefC.modifyAt(ichan).value = occChannel->GetBinContent(bin);
        
        //cout << ichip << "\t " << ichan << "\t" << occAtOnefC.getAt(ichan).value << endl;
        
        // check for defects:
        if (occAtOnefC.getAt(ichan).value > StandardDefects::NO_HI.getParameter() ){
        defects.addDefect(Defect(StandardDefects::NO_HI, ModuleElement::Channel(channelInModule)));
                occAtOnefC.modifyAt(ichan).valid=false;
        }
    }
    double mean = occAtOnefC.mean();
    //cout << "NoiseOcc mean = " << mean << endl;
    double rms  = sqrt( occAtOnefC.var() );
    
    return ChipNOResult(g,shared_ptr<TF1>(f),mean,rms);
    }
} // end of namespace SctAnalysis

---