In This Package:

QsumEnergyTool.cc

Go to the documentation of this file.

#include "QsumEnergyTool.h"

#include "Event/CalibReadout.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/CalibReadoutPmtChannel.h"
#include "Event/RecTrigger.h"
#include "Conventions/Reconstruction.h"
#include "Conventions/Detectors.h"

#include "DataSvc/ICableSvc.h"
#include "DataSvc/ICalibDataSvc.h"

using namespace DayaBay;

QsumEnergyTool::QsumEnergyTool(const std::string& type,
                               const std::string& name, 
                               const IInterface* parent)
  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_calibDataSvc(0)
{
    declareInterface< IReconTool >(this) ;   
    declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc",
                    "Name of service to map between detector, hardware, and electronic IDs");
    declareProperty("CalibDataSvcName", m_calibDataSvcName="StaticCalibDataSvc",
                    "Name of calibration data service");
}

QsumEnergyTool::~QsumEnergyTool(){}

StatusCode QsumEnergyTool::initialize()
{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
  m_calibDataSvc = svc<ICalibDataSvc>(m_calibDataSvcName,true);

  return StatusCode::SUCCESS;
}

StatusCode QsumEnergyTool::finalize()
{
  return StatusCode::SUCCESS;
}

StatusCode QsumEnergyTool::reconstruct(const CalibReadout& readout,
                                       RecTrigger& recTrigger)
{
  if( !readout.detector().isAD() ){
    debug() << "Not an AD readout; ignoring detector " 
            << readout.detector().detName() << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kNotProcessed );
    return StatusCode::SUCCESS;
  }

  const CalibReadoutPmtCrate* pmtReadout 
    = dynamic_cast<const CalibReadoutPmtCrate*>(&readout);
  if(!pmtReadout){
    error() << "Incorrect type of readout crate for detector "
            << readout.detector().detName() << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kBadReadout );
    return StatusCode::FAILURE;
  }

  // Context for this data
  int task = 0;
  ServiceMode svcMode(readout.header()->context(), task);

  // Loop over hits and add up charge
  double qSum = 0;
  double qSquareSum = 0;
  int nHits = 0;
  CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
    chanEnd = pmtReadout->channelReadout().end();
  for(chanIter=pmtReadout->channelReadout().begin(); chanIter != chanEnd; 
      chanIter++){
    const CalibReadoutPmtChannel& channel = *chanIter;
    AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
    //AdPmtSensor pmtId = m_cableSvc->adPmtSensor(channel.channelId(), svcMode);
    if( pmtId.ring()<1 ){ continue; }  // Ignore calibration pmts
    const PmtCalibData* pmtCalib = m_calibDataSvc->pmtCalibData(pmtId, svcMode);
    if( !pmtCalib ){ 
      error() << "No calibration data for pmt ID: " << pmtId << endreq;
      recTrigger.setEnergyStatus( ReconStatus::kBadReadout );
      return StatusCode::FAILURE;
    }
    if( pmtCalib->m_status != PmtCalibData::kGood ) { continue; }
    double peakAdc = channel.maxCharge();
    qSum += peakAdc;
    qSquareSum += peakAdc*peakAdc;
    nHits++;
  }

  if(nHits<1){
    recTrigger.setEnergyStatus( ReconStatus::kNoHits );
    return StatusCode::SUCCESS;
  }

  double pePerMeV = 116.0; // FIXME: add to calibration data service
  double energy = qSum / pePerMeV;
  recTrigger.setEnergy( energy );
  // Set quality to standard deviation of hits
  recTrigger.setEnergyQuality(sqrt((qSquareSum - qSum*qSum/nHits)/(nHits-1)));
  recTrigger.setEnergyStatus( ReconStatus::kGood );

  return StatusCode::SUCCESS;
}

---