OccupancyProjector.cpp

#include "Sct/SctParameters.h"
#include "../ScanResultWriter/dataTypes.h"
#include "OccupancyProjector.h"
#include "DefectList.h"
#include "TH1.h"
#include "TH2.h"

namespace SctData{
    using namespace Sct;

    OccupancyProjector::OccupancyProjector() throw() : raw(0) {
    unsetVetos();
    }

    OccupancyProjector::OccupancyProjector(const RawScanResult& r) throw() : raw(&r) {
    unsetVetos();
    }

    OccupancyProjector::~OccupancyProjector() throw() {;}
    
    void OccupancyProjector::vetoAllDefects(const DefectList& l) throw(){
    defects=&l; m_vetoAllDefects=true;
    }
    void OccupancyProjector::vetoSeriousDefects(const DefectList& l) throw(){
    defects=&l;
    }
    //void OccupancyProjector::setRawScanResult(const RawScanResult& r) throw(){
    //  raw=&r; unsetVetos();
    //}
    void OccupancyProjector::vetoMaskedChannels(const ModuleConfiguration& c) throw(){
    config=&c;
    }
    void OccupancyProjector::unsetVetos() throw(){
    config=0; defects=0; m_vetoAllDefects = false;
    }

    auto_ptr<TH1> OccupancyProjector::getOccupancy(const char* name, const ModuleElement& element) const throw (LogicError) {
    if (!raw) throw IllegalStateError("OccupancyProjector::getOccupancy no raw result", __FILE__, __LINE__);

    if ( raw->getDataType() == SR_DT_RAWHIST) throw IllegalStateError("Occupany project data in SR_DT_RAWHIST format????", __FILE__, __LINE__);

    double* bins = raw->getPoints().getEdgesAscending();
    auto_ptr<TH1> occ ( new TH1D(name, "Raw data projection", raw->getPoints().getNPoints(), bins) );
    delete[] bins;
    
    // project data in ascending order in scan point.
    bool ascending=raw->getPoints().ascending();

    //Explicit binomial errors.
    for (unsigned int j=0; j<raw->getPoints().getNPoints(); ++j){
      unsigned ipt=ascending ? j : raw->getPoints().getNPoints()-1-j;
      Stat_t n=raw->getPoints().getNEvents(ipt); 
      double z=0;
      unsigned int nChans = 0;
        
        //Do projection, taking into account that we might want data on both links
        for (unsigned int i=element.getFirst(); i<=element.getLast(); ++i) {
        if (config && config->channelIsMasked(i)) {
            //cout << "Channel: " << i << " masked, so skipping" << endl;
            continue;
        }
        if (defects) {
            DefectSeverity s = defects->defectSeverityAffectingElement(ModuleElement::Channel(i));
            if (m_vetoAllDefects && s != NONE ) continue;
            else if (s >= SERIOUS) {
            //if (nChans>1 && j==0) 
            //cout << "Skipping channel " << i << endl;
            continue;
            }
        }
        // Make sure to get correct scan data TH2:
        z += raw->getScanData(i/nChannelLink).GetBinContent(i%nChannelLink+1, j+1);
        nChans++;
        }               
        
        if (n==0 || nChans==0) {
        // no triggers.
        //cout << "No triggers?! " << n << " Chans: " << nChans << endl;
        occ->SetBinContent(j+1, z);
        occ->SetBinError(j+1, 0);
        continue;
        }
        n*=nChans;                         // number of triggers corrected for #channels
        
        double err2=(z+1)*(n-z+1)/(n*n*n);
        double occupancy=z/n;
        
        //if (nChans>1 && j==0) cout << "Number of channels: " << nChans << " Occ: " << occupancy << " tot: " << z << endl;
        
        occ->SetBinContent(j+1, occupancy);
        occ->SetBinError(j+1,sqrt(err2));
    }       
    return occ;
    }
    
}// end of namepace SctData

---