In This Package:

FastQCtrTool.cc

Go to the documentation of this file.

#include "FastQCtrTool.h"

#include "Event/RecHeader.h"
#include "Event/RecTrigger.h"

#include "Event/CalibReadoutHeader.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/CalibReadoutPmtChannel.h"

#include "CLHEP/Vector/ThreeVector.h"
#include "DetHelpers/IPmtGeomInfoSvc.h"
#include "DataSvc/ICableSvc.h"
#include "Conventions/Detectors.h"
#include "Conventions/Reconstruction.h"

using namespace CLHEP;
using namespace DayaBay;

FastQCtrTool::FastQCtrTool(const std::string& type,
      const std::string& name,
      const IInterface* parent)
  :  GaudiTool(type, name, parent)
  ,  m_cableSvc(0)
  ,  m_pmtGeomSvc(0)
  ,  m_usedChannels()
{
  declareInterface< IReconTool >(this);
  declareProperty("CableSvcName", m_cableSvcName = "StaticCableSvc",
    "Name of service to map between detector, hardware, and electronic IDs");
  declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc",
    "Name of Pmt Geometry Information Service");
}

FastQCtrTool::~FastQCtrTool() {}

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

  // Get PmtGeomInfo Service
  m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);

  return StatusCode::SUCCESS;
}


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

StatusCode FastQCtrTool::reconstruct(const DayaBay::CalibReadout& calibReadout, 
                                     DayaBay::RecTrigger& recTrigger)
{
  debug() << "running reconstruct() " << endreq;

  if(!calibReadout.detector().isAD()){
    debug() << "Ignore non-AD detector: " << calibReadout.detector().detName()
            << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kNotProcessed);
    recTrigger.setPositionStatus( ReconStatus::kNotProcessed);
    return StatusCode::SUCCESS;
  }

  // Get CalibPmtCrate
  const DayaBay::CalibReadoutPmtCrate* calibCrate
    = dynamic_cast<const DayaBay::CalibReadoutPmtCrate*>(&calibReadout);
  
  // Check CalibPmtCrate
  if(calibCrate == NULL) {
    error() << "reconstruction():  no calib readouts." << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kBadReadout);
    recTrigger.setPositionStatus( ReconStatus::kBadReadout);
    return StatusCode::FAILURE;
  }
  
  // Get CalibPmtChannelReadouts
  DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts calib_chMap;
  calib_chMap = calibCrate->channelReadout();
  DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator cpcrIter 
    = calib_chMap.begin();

  debug() << "successfully get PmtChannelReadouts " << endreq;
  
  //  Context, ServiceMode for this data
  int task = 0;
  ServiceMode svcMode(calibReadout.header()->context(), task);

  debug() << "successfully get service mode " << endreq;

  // determine FeeChannels that used for reconstruction
  setUsedChannels(calibCrate, &recTrigger);

  // Calculate the charge center, loop over all calibrated channels
  double sumQ = 0.0, Xc, Yc, Zc;
  
  for(; cpcrIter!=calib_chMap.end(); ++cpcrIter) {

    DayaBay::AdPmtSensor adPmtId(cpcrIter->pmtSensorId().fullPackedData());
    if(m_usedChannels[adPmtId]==false) continue;

    // Get charge on this PMT
    double charge = cpcrIter->maxCharge();

    // Get Pmt Position from its id
    unsigned int pmtid = adPmtId.fullPackedData();
    if(adPmtId.site() == Site::kSAB) {
       pmtid = DayaBay::AdPmtSensor(adPmtId.ring(),adPmtId.column(),Site::kDayaBay,adPmtId.detectorId()).fullPackedData();
    }
    const CLHEP::Hep3Vector pos = m_pmtGeomSvc->get(pmtid)->localPosition();
    
    Xc   +=  charge * pos.x();
    Yc   +=  charge * pos.y();
    Zc   +=  charge * pos.z();
    sumQ +=  charge;
  }
  if(sumQ!=0) { Xc /= sumQ; Yc /= sumQ; Zc /= sumQ;}
  double position_unit= CLHEP::millimeter;
  debug() << "Xc, Yc, Zc (mm): " << Xc/position_unit
         << ", " << Yc/position_unit << ", " << Zc/position_unit << endreq;

  CLHEP::HepLorentzVector Qcenter(Xc, Yc, Zc, 0);
  recTrigger.setEnergy(sumQ/100.);
  recTrigger.setPosition(Qcenter);
  
  return StatusCode::SUCCESS;
}

StatusCode FastQCtrTool::setUsedChannels(
                               const DayaBay::CalibReadoutPmtCrate* calibCrate,
                               DayaBay::RecTrigger* recTrigger)
{
  m_usedChannels.clear();
  
  Context context = calibCrate->header()->context();
  int task = 0;
  ServiceMode svcMode(context, task);

  Site::Site_t site = calibCrate->detector().site();
  DetectorId::DetectorId_t  detid = calibCrate->detector().detectorId();

  // default is use all channels
  for(int local_id=0; local_id<192; local_id++) {
    int ring = local_id/24 + 1;
    int col  = local_id%24 + 1;

    AdPmtSensor adPmtId(ring, col, site, detid);
    m_usedChannels[adPmtId] = true;
    //FeeChannelId channelId = m_cableSvc->feeChannelId(adPmtId, svcMode);
    //m_usedChannels[channelId] = true;
  }
  // Other options, depending on CalibReadoutPmtCrate data

  //recTrigger->setUsedChannels(m_usedChannels);

  return StatusCode::SUCCESS;
}

---