In This Package:

CenterOfChargePosTool.cc

Go to the documentation of this file.

#include "CenterOfChargePosTool.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"
#include "DetHelpers/IPmtGeomInfoSvc.h"
#include "CLHEP/Vector/LorentzVector.h"

#include "TMath.h"

using namespace DayaBay;

CenterOfChargePosTool::CenterOfChargePosTool(const std::string& type,
                                             const std::string& name, 
                                             const IInterface* parent)
  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_calibDataSvc(0)
  , m_pmtGeomSvc(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");
    declareProperty("PmtGeomSvcName", m_pmtGeomSvcName="PmtGeomInfoSvc",
                    "Name of PMT geometry data service");
    declareProperty("ApplyScaling", m_applyScaling=true,
           "Apply empirical scaling factor to get a 'reconstructed' position");
}

CenterOfChargePosTool::~CenterOfChargePosTool(){}

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

  return StatusCode::SUCCESS;
}

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

StatusCode CenterOfChargePosTool::reconstruct(const CalibReadout& readout,
                                              RecTrigger& recTrigger)
{
  if( !readout.detector().isAD() ){
    debug() << "Not an AD readout; ignoring detector " 
            << readout.detector().detName() << endreq;
    recTrigger.setPositionStatus( 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.setPositionStatus( ReconStatus::kBadReadout );
    return StatusCode::FAILURE;
  }

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

  // Loop over hits and determine charge-weighted mean position
  double pos[3];
  double posSquare[3];
  for(int idx = 0; idx<3; idx++){ pos[idx]=0; posSquare[idx]=0; }
  double time = 0;
  double qSum = 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());
    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.setPositionStatus( ReconStatus::kBadReadout );
      return StatusCode::FAILURE;
    }
    if( pmtCalib->m_status != PmtCalibData::kGood ) { continue; }
    const CLHEP::Hep3Vector& pmtPos 
      = m_pmtGeomSvc->get(pmtId.fullPackedData())->localPosition();
    double peakAdc = channel.maxCharge();
    double firstTdc = channel.earliestTime();
    for(int idx = 0; idx<3; idx++){
      pos[idx] += peakAdc*pmtPos[idx];
      posSquare[idx] += peakAdc*pmtPos[idx]*pmtPos[idx];
    }
    time += firstTdc;
    qSum += peakAdc;
    nHits++;
  }

  if(nHits<1){
    recTrigger.setPositionStatus( ReconStatus::kNoHits );
    return StatusCode::SUCCESS;
  }
  double posStdDev = 0;
  for(int idx = 0; idx<3; idx++){
    posStdDev += posSquare[idx] - pos[idx]*pos[idx]/qSum;
    pos[idx] /= qSum;
  }
  posStdDev /= qSum;
  posStdDev = sqrt(posStdDev);
  time /= qSum;

  if(m_applyScaling){
    // Apply empirical rho, z scaling to vertex position
    // This is not very good, but is useful as a prefit
    double rho = sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
    double rhoNorm = rho/2340.;
    double zNorm = pos[2]/1750.;
    if(zNorm<0) zNorm *= -1;
    double rhoScaled = 1-(1-rhoNorm)/TMath::Power((1+rhoNorm), 3.13);
    double zScaled = 
      1-(1-zNorm)/TMath::Power((1+zNorm), 
                               10.0*(1-sqrt(rhoNorm)));
    //warning() << "Scaling by: " << rhoScaled << " , " << zScaled << endreq;
    pos[0] *= rhoScaled/rhoNorm;
    pos[1] *= rhoScaled/rhoNorm;
    pos[2] *= zScaled/zNorm;
  }

  CLHEP::HepLorentzVector position(pos[0], pos[1], pos[2], time);
  recTrigger.setPosition( position );
  // Set quality to charge-weighted standard deviation of hit positions
  recTrigger.setPositionQuality( posStdDev );
  recTrigger.setPositionStatus( ReconStatus::kGood );

  return StatusCode::SUCCESS;
}

---