NPtGainAlgorithm.cpp

#include "NPtGainAlgorithm.h"
#include "AnalysisAlgorithmMap.h"
#include "AnalysisService.h"
#include "Sct/SctParameters.h"
#include "SctData/NPtGainTestResult.h"
#include "SctData/FitScanResult.h"
#include "SctData/FitObject.h"
#include "SctData/ResponseCurve.h"
#include "SctData/StandardDefects.h"
#include "SctData/ModuleConfiguration.h"
#include "SctData/ChipConfiguration.h"
#include "Sct/SctNames.h"
#include "Sct/ConfigurationException.h"
#include "CalibrationController/IS/TestData.h"
#include "sctConf/configuration.h"
#include <boost/shared_ptr.hpp>

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

using std::auto_ptr;
using std::ostringstream;

namespace SctAnalysis {

unsigned NPtGainAlgorithm::NPtGainAlgorithm::s_configurationWarnings=0;

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

NPtGainAlgorithm::NPtGainAlgorithm(shared_ptr<const TestData> testData, const string& moduleName, const AnalysisAlgorithm& alg) throw() : AnalysisAlgorithm(testData, moduleName, alg) , m_mur_type(SctConfiguration::UNKNOWN), m_endcap_type(SctConfiguration::ENDCAP_UNKNOWN){
  bool error=false;
  string error_message="unknown problem";
  const unsigned max_warnings=5;

 try{
    boost::shared_ptr<SctConfiguration::Configuration> config=AnalysisService::instance().getConfiguration();
    if (!config.get()) throw ConfigurationException("Configuration not available", __FILE__, __LINE__);
    unsigned mur=0, mod=0, quadrant=0, number=0;
    int disk=0;
    config->translateFromSN(moduleName, mur, mod);
    m_mur_type=config->getMURType(mur);
    if (m_mur_type==SctConfiguration::ENDCAP){
      config->translateToEndcap(mur, mod, disk, quadrant, number);
      m_endcap_type= SctConfiguration::getEndcapType(disk, number);
    }
  }catch(std::exception& e){
    error=true;
    error_message=e.what();
  }catch(CORBA::Exception& e){
    error=true;
    error_message="Corba exception : ";
    error_message+=e._name();
    error_message+=" ... probably could not find configuration service";
  }catch(...){
    error=true;
  }
  if (error){
    // only send configuration warning the first few times!
    if (s_configurationWarnings<max_warnings){
      ostringstream msg;
      if (s_configurationWarnings==max_warnings-1){
    msg << "NO MORE WARNINGS OF THIS TYPE WILL BE PRINTED:\n";
      }
      msg << "NPtGainAlgorithm (AnalysisService) could not find geometry-dependant expected noise levels for test " 
      << testData->testName
      << " run:" << testData->runNumber << " scan:" << testData->startScanNumber
      << " module:" << moduleName 
      << " Problem in configuration look-up: " 
      << error_message;
      Sct::SctNames::Mrs() << MRS_TEXT(msg.str()) << "CONFIG_LOOKUP_FAIL" << MRS_WARNING << ENDM;
      s_configurationWarnings++;
    }
  }
}
    
shared_ptr<ResponseCurve> NPtGainAlgorithm::getResponseCurve(unsigned int nPts) throw(LogicError){
    if (nPts <= 3) {
    return shared_ptr<ResponseCurve> (new LinearResponseCurve());
    } else {
    return shared_ptr<ResponseCurve> (new ExponentialResponseCurve());
    }   
}
    
void NPtGainAlgorithm::setResponseCurve(auto_ptr<ResponseCurve> rc) throw() {
    static shared_ptr<ResponseCurve> s_responseCurve=shared_ptr<ResponseCurve>(rc.release());
}

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

shared_ptr<SctData::TestResult> NPtGainAlgorithm::createTestResult() const {
    shared_ptr<NPtGainTestResult> result (new NPtGainTestResult());
    result->setChipDataSize(nChipModule);
    result->setChannelDataSize(nChannelModule);
    
    return result;
}

void NPtGainAlgorithm::loadData() {
    loadAllFits();
}

bool NPtGainAlgorithm::canAnalyze() const {
    return hasAllFits();
}

void NPtGainAlgorithm::analyze() {
    shared_ptr<NPtGainTestResult> result = dynamic_pointer_cast<NPtGainTestResult> ( getTestResult() );
    result->setScanVariable(getFit(0)->getHeader().getVariable());  
    try{
    result->setSpecialScanPointValue(2.0);
    } catch (Sct::Throwable& e){
    e.sendToMrs(MRS_ERROR);
    }
    
    //Reset counters
    nOnePointTrimmed = 0;
    nTwoPointTrimmed = 0;
    
    //Do chips
    for (unsigned int i=0; i<nChipModule; ++i) {
        setupGraph(i, ModuleElement::Chip(i), &FitScanResult::getChipFit, *result, result->getChipData(i), true);
        doFit(i, &FitScanResult::getChipFit, *result, result->getChipData(i));
    }

    //Do channels
    for (unsigned int i=0; i<nChannelModule; ++i) {
    ModuleElement element = ModuleElement::Channel(i);
        setupGraph(i, element, &FitScanResult::getChannelFit, *result, result->getChannelData(i), true);
    //Note - make sure we do this after we have setup the graph - that way we don't cause any problems with missing data
    //And we can check later if necessary
    if (getFit(0)->getConfiguration().channelIsMasked(i)) continue; 
    
        doFit(i, &FitScanResult::getChannelFit, *result, result->getChannelData(i));
        doDefect(element, result->getDefects(), 
         result->getChannelData(i), result->getChipData(i/nChannelChip));
    }
    
    doSlopes(*result);

    //mergeDefects();
    
    //Check for defects - fail if more than 15 defective channels or more than 7 consecutive defective channels
    const DefectList::DefectCollection& defects = result->getDefects().getAllDefects();
    bool passed = true;
    unsigned int totalChannels = 0;
    for (DefectList::DefectCollection::const_iterator i=defects.begin(); i!=defects.end(); ++i) {
        // DODGY defects don't count towards failure of module
        if (i->getPrototype().getSeverity() > DODGY) {
            totalChannels += i->getModuleElement().getNChannels();
            if (i->getModuleElement().getNChannels()>7) {
                passed = false;
                break;
            }
        }
    }
    if (totalChannels > 15) passed = false;
    
    //Add comment about number of channels trimmed
    if (nOnePointTrimmed != 0) {
    ostringstream oss;
    oss << nOnePointTrimmed << " channels had 1 point trimmed from fit";
    result->addComment(oss.str());
    }
    if (nTwoPointTrimmed != 0) {
    ostringstream oss;
    oss << nTwoPointTrimmed << " channels had 2 (or more) points trimmed from fit";
    result->addComment(oss.str());
    }
    
    result->setPassed(passed);
}
 
//Setup the graphs - also prune points
void NPtGainAlgorithm::setupGraph(unsigned int id, const ModuleElement& element, getFitFunction fitFunc, NPtGainTestResult& test, NPtGainTestResultData& testData, bool trimPoints) throw(LogicError) {
    
    unsigned int lastPoint = test.getNScans()-1;
    //Filter out points above 5fC if they are broad and we have more than 5 points
    if (trimPoints && test.getNScans() > 5) {
    double sigmaCut = 0;
    for (unsigned int i=0; i<test.getNScans(); ++i) {
        sigmaCut+= (getFit(i).get()->*fitFunc)(id).getParameter("Sigma");
    }
    sigmaCut *= 1.5/test.getNScans();   //We want 1.5*mean sigma
    //cout << "Sigma Cut set to: " << sigmaCut << endl;
    for (unsigned int i=test.getNScans() - 1; i>4; --i) {
        //cout << "Point: " << i << " sigma: " << (getFit(i).get()->*fitFunc)(id).getParameter("Sigma") << " TestPoint: " << test.getTestPointAt(i) << endl;
        const DefectList& defects = getFit(i)->getDefects();
        if (test.getTestPointAt(i) > 5 && ((getFit(i).get()->*fitFunc)(id).getParameter("Sigma") > sigmaCut || 
                           defects.defectSeverityAffectingElement(element) >= SERIOUS)) {
        lastPoint = i-1;
        }
    }
    
    //Note we should only record trimed points for chips... 
    if (lastPoint != test.getNScans() - 1) {
        if (element.isChip()) {
        ostringstream oss;
        oss << "Points trimed from response curve for chip: " << id << ". Last point " << lastPoint << " with charge " << test.getTestPointAt(lastPoint);
        test.addComment(oss.str());
        } else {    //Element is a channel
        if (lastPoint+1-test.getNScans() == 1) ++nOnePointTrimmed;
        else ++nTwoPointTrimmed;
        }
    }
    }
    
    mergeDefects(test, element, lastPoint);

    shared_ptr<TGraph> g (new TGraph(lastPoint+1));
    if (!g)
        throw InvariantViolatedError("NPtGainTestResult::setupGraph couldn't make TGraph", __FILE__, __LINE__);
    
    
    for (unsigned int j = 0; j <= lastPoint; ++j) {
        // Cal charge register is set in RodDaq/Dsp/Sct/Code/ABCDConfig.c
        // find out which chip from the module element
        unsigned ichip=element.getFirst()/nChannelChip;
    const ChipConfiguration& chipconfig = getFit(j)->getConfiguration().getChipConfiguration(ichip);
    // get the cal charge register value
        const unsigned char cal = chipconfig.getCalCharge();
    // convert this to fC, taking into account the calibration factor:
        const float qcal= static_cast<float>(cal)*0.0625*chipconfig.getCalFactor();
        g->SetPoint(j, qcal, (getFit(j).get()->*fitFunc)(id).getParameter("Mean"));
    //if (element.getFirst()%nChannelChip==0){
    //  std::cout << "Chip " << element.getFirst()/nChannelChip << " Point (" << j << ") value requested = " << test.getTestPointAt(j)
    //      << " register = " << (int)cal 
    //      << " step = " << 0.0625*chipconfig.getCalFactor()
    //      << " actual value = " << qcal << std::endl; 
    //}
    }
    testData.graph = g;

    shared_ptr<ResponseCurve> r=getResponseCurve(test.getNScans());
    testData.rc = r;    
}

void NPtGainAlgorithm::mergeDefects(NPtGainTestResult& test, const ModuleElement& element, unsigned int lastPoint) {
    for (unsigned int i=0; i<lastPoint; ++i) {
    auto_ptr<DefectList> defects = getFit(i)->getDefects().getDefectsEncompassingElement(element);
    const DefectList::DefectCollection& allDefects = defects->getAllDefects();
    for (DefectList::DefectCollection::const_iterator it=allDefects.begin(); it!=allDefects.end(); ++it) {
        test.getDefects().addDefect(Defect(it->getPrototype(), element));
    }
    }
}

void NPtGainAlgorithm::doFit(unsigned int id, getFitFunction fitFunc, NPtGainTestResult& test,
                              NPtGainTestResultData& testData) throw(LogicError) {
    //testData.graph->Fit(&testData.rc->getFunction(), "NQ", "", 0, 10);
    // replaced below with generalised fitting, 17/7/03 Alan Barr

    testData.rc->getFunction()->SetRange(0,10);
    AnalysisService::instance().getFitStrategy().fitTGraph(*testData.graph, *testData.rc->getFunction() );
    
    testData.gain = testData.rc->getGain(test.getSpecialScanPointValue());

    shared_ptr<const FitScanResult> specialFit = getFit(test.getSpecialScanIndex());
    //6250 is a magic number that converts noise to ENC
    if (testData.gain != 0) testData.noise = 6250 * (specialFit.get()->*fitFunc)(id).getParameter("Sigma") / testData.gain;
    else  testData.noise = 0;
    testData.offset= testData.rc->getFunction()->Eval(0.);
}

void NPtGainAlgorithm::doDefect(const ModuleElement& element, DefectList& defects, 
                const NPtGainTestResultData& data,
                const NPtGainTestResultData& comparisonData) throw(LogicError) {
    
    if (data.gain < StandardDefects::VLO_GAIN.getParameter() * comparisonData.gain){
        defects.addDefect(Defect(StandardDefects::VLO_GAIN, element));
    }
    if (data.gain < StandardDefects::LO_GAIN.getParameter() * comparisonData.gain){
        defects.addDefect(Defect(StandardDefects::LO_GAIN, element));
    } else if (data.gain > StandardDefects::HI_GAIN.getParameter() * comparisonData.gain){
        defects.addDefect(Defect(StandardDefects::HI_GAIN, element));
    }

    if (data.offset < StandardDefects::LO_OFFSET.getParameter()) {
        defects.addDefect(Defect(StandardDefects::LO_OFFSET, element));
    } else if (data.offset > StandardDefects::HI_OFFSET.getParameter()) {
        defects.addDefect(Defect(StandardDefects::HI_OFFSET, element));
    }

    bool short_module = ( m_mur_type==SctConfiguration::ENDCAP && 
              ( m_endcap_type==SctConfiguration::ENDCAP_SHORT_MIDDLE
                || m_endcap_type==SctConfiguration::ENDCAP_INNER) );
    
    if (data.noise < StandardDefects::UNBONDED.getParameter()) {
        defects.addDefect(Defect(StandardDefects::UNBONDED, element));
    } else if (data.noise < StandardDefects::PARTBONDED.getParameter() && !short_module) {
        defects.addDefect(Defect(StandardDefects::PARTBONDED, element));
    } else if (data.noise > StandardDefects::NOISY.getParameter() * comparisonData.noise) {
        defects.addDefect(Defect(StandardDefects::NOISY, element));
    }
  
    if (data.noise > StandardDefects::V_NOISY.getParameter() * comparisonData.noise) {
        defects.addDefect(Defect(StandardDefects::V_NOISY, element));
    }
}

  void NPtGainAlgorithm::doSlopes(NPtGainTestResult& r){
    shared_ptr<TGraph> noise  = r.getNoiseGraph();
    shared_ptr<TGraph> gain   = r.getGainGraph();
    shared_ptr<TGraph> offset = r.getOffsetGraph();

    //std::cout << "Noise" << std::endl;
    doSlopes(noise,  r.noiseSlope);
    
    //std::cout << "Gain" << std::endl;
    doSlopes(gain,   r.gainSlope);

    //std::cout << "Offset" << std::endl;
    doSlopes(offset, r.offsetSlope);
    for (unsigned ichip=0; ichip<nChipModule; ++ichip){
      if (abs(r.noiseSlope[ichip]) > 0.7 || r.noiseSlope[ichip] < -1.2){
    r.getDefects().addDefect(Defect(StandardDefects::NOISE_SLOPE, ModuleElement::Chip(ichip)));
      }
      if (abs(r.gainSlope[ichip])>StandardDefects::GAIN_SLOPE.getParameter()){
    r.getDefects().addDefect(Defect(StandardDefects::GAIN_SLOPE, ModuleElement::Chip(ichip)));
      }
      if (abs(r.offsetSlope[ichip])>StandardDefects::OFFSET_SLOPE.getParameter()){
    r.getDefects().addDefect(Defect(StandardDefects::OFFSET_SLOPE, ModuleElement::Chip(ichip)));
      }

    }

  }

  void NPtGainAlgorithm::doSlopes(shared_ptr<TGraph> graph, Sct::RangedVector<float>& result){
    result.resize(nChipModule);
    // linear least squares:
    for (unsigned ichip=0; ichip<nChipModule; ++ichip){
      double xsum=0.;
      double x2sum=0.;
      double ysum=0.;
      double xysum=0.;
      unsigned n=0;
      for (unsigned ichan = 0; ichan<nChannelChip; ++ichan) {
    double x=0;
    double y=0;
    int bin=ichip*nChannelChip + ichan;
    graph->GetPoint(bin, x, y);
    if (y>0.0001){
      xsum  += x;
      ysum  += y;
      x2sum += x*x;
      xysum += x*y;
      ++n;
    }
      }
      const double det = ( n*x2sum ) - ( xsum*xsum );
      if (det>0) {
    result[ichip]= (n*xysum - xsum*ysum) / det;
      }else{
    result[ichip]=0.;
      }
    }
    }

}

---