---

NoiseOccupancyAlgorithm.cpp

Go to the documentation of this file.

#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/ModuleDefect.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));
    }

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

    //----------------------------------------------------------------------
    
    void NoiseOccupancyAlgorithm::init() throw(){
    if (getTestResult().get()==0){
        setTestResult(shared_ptr<TestResult>(new NoiseOccupancyTestResult() ));
        initializeTestResult();
    }
    }

    void NoiseOccupancyAlgorithm::canAddFitScanResult(const string& name) throw(Sct::LogicError, Sct::IoError){
    init();
    addFit(name);
    analyze();
    }

    void NoiseOccupancyAlgorithm::canAddRawScanResult(const string& name) throw(Sct::LogicError, Sct::IoError){
    init();
    addRaw(name);
    analyze();
    }

    void NoiseOccupancyAlgorithm::analyze() throw(Sct::LogicError){
    if (!getTestResult()->hasAllFits() || !getTestResult()->hasAllRaws() ) return;

    shared_ptr<NoiseOccupancyTestResult> result = dynamic_pointer_cast<NoiseOccupancyTestResult>(getTestResult());

    const ModuleConfiguration config=result->getRaw(0)->getConfiguration();
    
    //Loop over chips
    for (unsigned int ichip=0; ichip<nChipModule; ++ichip) {
        // get an raw result occupany projector 
        OccupancyProjector opr(*result->getRaw(0) );
        opr.vetoSevereDefects(result->getFit(0)->getDefects());
        // get the chip configuration
        const ChipConfiguration& chipconfig=config.getChipConfiguration(ichip);
        // get the new defects
        ModuleDefectList& defects = result->getDefects();
        result->getChipResult(ichip) = analyzeChip(ichip, opr, chipconfig, defects);
        result->getChipResult(ichip).offset = result->getFit(0)->getChipFit(ichip).getParameter(1);
    }
    result->setPassed(true);
    finish();
    }
    
    ChipNOResult NoiseOccupancyAlgorithm::analyzeChip(const unsigned ichip, const OccupancyProjector& opr, const ChipConfiguration& chipconfig, ModuleDefectList& defects) throw(LogicError){
    
    auto_ptr<TH1> occ = opr.getOccupancy("temp", Chip(ichip)) ;
    auto_ptr<ResponseCurve> rc = ResponseCurveMap::getMap().get(chipconfig);
    shared_ptr<TF1> invFn = rc->getInverseFunction();
    
    //cout << "Response Curve params: " << chipconfig.getRcParam(0) << "  " << chipconfig.getRcParam(1) << "  " << chipconfig.getRcParam(2) << 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()];
    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 thr in mV to thr in fC
        x[nBins] = invFn->Eval(occ->GetBinCenter(bin));
        x[nBins] *= x[nBins];
        xUpperErr[nBins] = xLowerErr[nBins] = 0;            
                
        nBins++;
        }
    }
    shared_ptr<TGraphAsymmErrors> g( new TGraphAsymmErrors(nBins, x, y, xLowerErr, xUpperErr, yLowerErr, yUpperErr));
    
    auto_ptr<TF1> f = NoiseOccupancyTestResult::createFitFunction();

//  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, Channel(channelInModule)) ;
        // find 1 fC point:
        double vthr = rc->getFunction()->Eval(1.);
        int bin = occChannel->GetXaxis()->FindBin(vthr);

//      cout << "Vthr=" << vthr << " 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 > ModuleDefect::NO_HI.getParameter() ){
        defects.addDefect(ModuleDefect::NO_HI,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

---