In This Package:

RawData2Tree.cc

Go to the documentation of this file.

#include "RawData2Tree.h"
#include "Event/CalibReadoutHeader.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/ReadoutHeader.h"
#include "Event/ReadoutPmtCrate.h"
#include "Event/RawEventHeader.h"
#include "Event/RawRom.h"
#include "Event/RawRomLtb.h"
#include "Event/RawPmtChannel.h"
#include "Event/GenHeader.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
#include "Event/RecHeader.h"
#include "DetHelpers/IPmtGeomInfoSvc.h"
#include "DataSvc/ICableSvc.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "Event/DaqLtb.h"
#include "Event/DaqLtbFrame.h"
#include <map>

using namespace std;
using namespace DayaBay;

RawData2Tree::RawData2Tree(const string& name, ISvcLocator* svcloc)
: GaudiAlgorithm(name, svcloc)
, m_cableSvc(0)
, m_pmtGeomSvc(0)
{
  declareProperty("PrintFreq", m_printFreq = -1, "print frequency for event information");
  declareProperty("CheckGen", m_checkGen = false, "check GenEvent");
  declareProperty("CheckReadout", m_checkReadout = false, "check ReadoutEvent");
  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");
  m_execNum = 0;
  m_ntuple0 = 0;
  m_ntuple1 = 0;
}

RawData2Tree::~RawData2Tree()
{
}

StatusCode RawData2Tree::initialize()
{
  debug() << "initialize()" << endreq;

  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName, true);

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

  // Book Ntuple
  StatusCode sc = bookNtuple();
  if(sc!=StatusCode::SUCCESS) {
    error() << "Can't book ntuple!" << endreq;
  }

  return sc;
}

StatusCode RawData2Tree::execute()
{
  debug() << "execute() ______________________________ start" << endreq;

  if(m_printFreq>0 && m_execNum%m_printFreq==0) cout << "---------- " << m_execNum << endl;

  //----------------------------------------------------------------------------------------
  //GenHeader

  if(m_checkGen) {

    GenHeader* gh = get<GenHeader>("/Event/Gen/GenHeader");

    if (gh == 0) {
      error() << " =======> Requested Object can not be accessable." << endreq;
      return StatusCode::FAILURE;
    }

    if (m_printFreq > 0 && m_execNum % m_printFreq == 0) {
      //m_log << MSG::INFO << "Gen Header: " << *gh << endreq;
      info() << "execNumber: " << gh->execNumber() << endreq;
      //gh->event()->print();
      info() << "signal_process_id= " << gh->event()->signal_process_id() << endreq;
      info() << "vertices_size= " << gh->event()->vertices_size() << endreq;
    }

    m_eventType = gh->event()->signal_process_id();
    m_nVertex = gh->event()->vertices_size();
    m_nParticle = gh->event()->particles_size();

    unsigned int iVertex = 0;
    for(HepMC::GenEvent::vertex_const_iterator it = gh->event()->vertices_begin();
        it != gh->event()->vertices_end();
        it++ ) {
      if(m_printFreq > 0 && m_execNum % m_printFreq == 0) {
        info() << "vertex= " << (*it)->position().x() 
              <<","<< (*it)->position().y()<<","<<(*it)->position().z()
              << endreq;
      }
      m_verX0[iVertex] = (*it)->position().x();
      m_verY0[iVertex] = (*it)->position().y();
      m_verZ0[iVertex] = (*it)->position().z();
      iVertex++;
    }

    unsigned int iParticle = 0;
    for(HepMC::GenEvent::particle_const_iterator it = gh->event()->particles_begin();
        it != gh->event()->particles_end();
        it++ ) {
      if(m_printFreq > 0 && m_execNum % m_printFreq == 0) {
        info() << "particle momentum= " 
              << (*it)->momentum().x() <<","
              << (*it)->momentum().y() <<","
              << (*it)->momentum().z() 
              << endreq;
      }
      m_verPx[iParticle] = (*it)->momentum().x();
      m_verPy[iParticle] = (*it)->momentum().y();
      m_verPz[iParticle] = (*it)->momentum().z();
      iParticle++;
    }

    m_ntuple0->write();

  }//GenHeader

  //----------------------------------------------------------------------------------------
  //ReadoutHeader

  if(m_checkReadout) {

    ReadoutHeader* roh = get<ReadoutHeader>("/Event/Readout/ReadoutHeader");
    if (roh == 0) {
      error() << " =======> Requested Object can not be accessable." << endreq;
      return StatusCode::FAILURE;
    }

    const DaqCrate* daqCrate = roh->daqCrate();
    if(!daqCrate){
      error() << "Failed to get DAQ crate from header" << endreq;
      return StatusCode::FAILURE;
    } 

    if(daqCrate) {
      const Detector& detector = daqCrate->detector();
      info() << "Detector = " << detector << endreq;

      if(detector.isAD() || detector.isWaterShield()) {
        processPmtCrate(roh,daqCrate);
      }
    } else if(roh->readout()) {
      processReadout(roh);
    } else {
      error() << "Failed to get DAQ crate or Readout from header" << endreq;
      return StatusCode::FAILURE;
    }

  }//ReadoutHeader

  m_execNum++;
  debug() << "execute() ______________________________ end" << endreq;
  return StatusCode::SUCCESS;
}

StatusCode RawData2Tree::finalize()
{
  debug() << "finalize()" << endreq;
  return StatusCode::SUCCESS;
}
StatusCode RawData2Tree::bookNtuple()
{
  //  Book NTuple 0
  if(m_checkGen) {
    NTuplePtr nt0(ntupleSvc(), "FILE1/Gen");
    if ( nt0 ) m_ntuple0 = nt0;
    else {
      m_ntuple0 = ntupleSvc()->book("FILE1/Gen", CLID_ColumnWiseTuple, "Gen");
      if ( m_ntuple0 )    {
        m_ntuple0->addItem ("EventType",      m_eventType);
        m_ntuple0->addItem ("VertexSize",     m_nVertex, 0, 30);
        m_ntuple0->addItem ("ver_x0",  m_nVertex, m_verX0);
        m_ntuple0->addItem ("ver_y0",  m_nVertex, m_verY0);
        m_ntuple0->addItem ("ver_z0",  m_nVertex, m_verZ0);
        m_ntuple0->addItem ("ParticleSize",   m_nParticle, 0, 30);
        m_ntuple0->addItem ("ver_px",  m_nParticle, m_verPx);
        m_ntuple0->addItem ("ver_py",  m_nParticle, m_verPy);
        m_ntuple0->addItem ("ver_pz",  m_nParticle, m_verPz);
      }
      else { // Did not manage to book the NTuple....
        error() << "Can not book NTuple:" << long(m_ntuple0) << endreq;
        return StatusCode::FAILURE;
      }
    }
  }

  //  Book NTuple 1
  if(m_checkReadout) {
    NTuplePtr nt1(ntupleSvc(), "FILE1/Readout");
    if ( nt1 ) m_ntuple1 = nt1;
    else {
      m_ntuple1 = ntupleSvc()->book("FILE1/Readout", CLID_ColumnWiseTuple, "Readout");
      if ( m_ntuple1 )    {
        info() << "Start m_ntuple1 addItem successful!" << endreq;
        // Trigger number serial
        m_ntuple1->addItem("Run",m_run);
        m_ntuple1->addItem("Event",m_event);
        m_ntuple1->addItem("Det",m_det);
        m_ntuple1->addItem("TrigNum",m_trigNum);
        m_ntuple1->addItem("TrigSeconde",m_trigSecond);
        m_ntuple1->addItem("TrigNanoSec",m_trigNanoSec);
        m_ntuple1->addItem("TrigType",  m_trigType);
        m_ntuple1->addItem("NChannel",m_nChannel);
        m_ntuple1->addItem("NHit",m_nHit,0,1000);
        m_ntuple1->addItem("Board",m_nHit,m_board);
        m_ntuple1->addItem("Channel",m_nHit,m_channel);
        m_ntuple1->addItem("Ring",m_nHit,m_ring);
        m_ntuple1->addItem("Column",m_nHit,m_column);
        m_ntuple1->addItem("Wall",m_nHit,m_wall);
        m_ntuple1->addItem("Spot",m_nHit,m_spot);
        m_ntuple1->addItem("Inout",m_nHit,m_inout);
        m_ntuple1->addItem("Adc",m_nHit,m_adc);
        m_ntuple1->addItem("Tdc",m_nHit,m_tdc);
        m_ntuple1->addItem("Cycle",m_nHit,m_cycle);
        m_ntuple1->addItem("PreAdc",m_nHit,m_preAdc);
        m_ntuple1->addItem("Range",m_nHit,m_range);
        m_ntuple1->addItem("HitCount",m_nHit,m_hitCount);
        m_ntuple1->addItem("MultiTdc",m_nHit,m_multiTdc);
        m_ntuple1->addItem("NTrigger",m_nTrigger,0,20);
        m_ntuple1->addItem("TriggerSec",m_nTrigger,m_triggerSec);
        m_ntuple1->addItem("TriggerNano",m_nTrigger,m_triggerNano);
        m_ntuple1->addItem("TriggerHSum",m_nTrigger,m_triggerHSum);
        m_ntuple1->addItem("TriggerESumComp",m_nTrigger,m_triggerESumComp);
        m_ntuple1->addItem("TriggerESumADC",m_nTrigger,m_triggerESumAdc);
      }
      else { // Did not manage to book the NTuple....
        error() << "Can not book NTuple:" << long(m_ntuple1) << endmsg;
        return StatusCode::FAILURE;
      }
    }
  }

  return StatusCode::SUCCESS;
}

void RawData2Tree::processPmtCrate(const ReadoutHeader *roh, const DaqCrate *daqCrate)
{
  info() << "processPmtCrate"  << endl;
  ServiceMode svcMode( roh->context(),0 );

  // Convert to PMT crate readout
  const DaqPmtCrate* pmtCrate = daqCrate->asPmtCrate();
  if(!pmtCrate){
    error() << "Can't find pmt crate for AD or Water Shield" << endreq;
  } 

  m_run = pmtCrate->runNumber();
  m_event = pmtCrate->eventNumber();
  m_det = pmtCrate->detector().fullPackedData();
  m_trigNum = pmtCrate->localTriggerNumber();
  m_trigSecond = pmtCrate->triggerTime().GetSec();
  m_trigNanoSec = pmtCrate->triggerTime().GetNanoSec();
  m_trigType = pmtCrate->triggerType();

  const DaqPmtCrate::PmtChannelPtrList& channels
    = pmtCrate->channelReadouts();

  m_nChannel= channels.size();
  DaqPmtCrate::PmtChannelPtrList::const_iterator channelIter,
    channelEnd = channels.end();

  // Calculate number of hits in this readout
  int nHit = 0;
  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) {
    const DaqPmtChannel& channel = *(*channelIter);
    nHit += channel.hitCount();
  }
  m_nHit = nHit;
  info() << "Number of channel: " << channels.size() 
    << ", Number of hit: " << nHit << endreq;

  int index = -1;
  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) {
    const DaqPmtChannel& channel = *(*channelIter);

    int board = channel.channelId().board(); 
    int connector = channel.channelId().connector();
    int ring = -1;
    int clmn = -1;
    int wall = -1;
    int spot = -1;
    bool inout = false;
    if( pmtCrate->detector().isAD() ) {
      AdPmtSensor pmtId = m_cableSvc->adPmtSensor( channel.channelId(), svcMode );
      ring = pmtId.ring();
      clmn = pmtId.column();
    } else if( pmtCrate->detector().isWaterShield() ) {
      PoolPmtSensor pmtId = m_cableSvc->poolPmtSensor( channel.channelId(), svcMode );
      wall = pmtId.wallNumber();
      spot = pmtId.wallSpot();
      inout = pmtId.inwardFacing();
    }
    debug() << "(board,connector) = (" << board << "," << connector << "), " 
      << "(ring,column) = (" << ring << "," << clmn << ")" << endreq;

    unsigned int hitSize = channel.hitCount();
    for(unsigned int i=0; i<hitSize; i++) {
      index++;
      m_board[index] = board;
      m_channel[index] = connector;
      m_ring[index] = ring;
      m_column[index] = clmn;
      m_wall[index] = wall;
      m_spot[index] = spot;
      m_inout[index] = inout ? 1 : 0;
      m_adc[index] = channel.adc(i);
      m_tdc[index] = channel.tdc(i);
      m_cycle[index] = channel.peakCycle(i);
      m_preAdc[index] = channel.preAdcRaw(i);
      m_range[index] = channel.isHighGainAdc(i) ? 1 : 2;
      m_hitCount[index] = i;
      m_multiTdc[index] = hitSize;
      debug() << "tdc=" << m_tdc[index] << ", adc=" << m_adc[index] << endreq;
    }
  }

  const DaqLtb &ltb = pmtCrate->localTriggerBoard();
  //const LtbFramePtrList frames = ltb.frames();
  const vector<DayaBay::DaqLtbFrame*> &frames = ltb.frames();
  if(frames.size()>20) {
    warning() << "Two many trigger frames: " << frames.size() << endreq;
  }

  for(unsigned int i=0; i<frames.size(); i++) {
    m_triggerSec[i] = frames[i]->triggerTime().GetSec();
    m_triggerNano[i] = frames[i]->triggerTime().GetNanoSec();
    m_triggerHSum[i] = frames[i]->hitSum();
    m_triggerESumComp[i] = frames[i]->energySum();
    m_triggerESumAdc[i] = frames[i]->energySum();
  }

  m_ntuple1->write();
}

void RawData2Tree::processReadout(const ReadoutHeader *roh)
{
  info() << "processReadout"  << endl;
  ServiceMode svcMode( roh->context(),0 );

  Readout* readout = const_cast<Readout*>(roh->readout());

  ReadoutPmtCrate* pmtCrate = dynamic_cast<DayaBay::ReadoutPmtCrate* > ( readout );

  Detector det = pmtCrate->detector();

  m_det = det.fullPackedData();
  m_trigNum = readout->triggerNumber();
  m_trigSecond = readout->triggerTime().GetSec();
  m_trigNanoSec= readout->triggerTime().GetNanoSec();
  m_trigType = readout->triggerType();

  const ReadoutPmtCrate::PmtChannelReadouts &channels = pmtCrate->channelReadout();
  m_nChannel = channels.size();
  ReadoutPmtCrate::PmtChannelReadouts::const_iterator ci, ci_end = channels.end();
  
  // Calculate number of hits in this readout
  int nHit = 0;
  for(ci = channels.begin(); ci!=ci_end; ci++) {
    nHit += ci->second.size();
  }
  m_nHit = nHit;
  info() << "Number of channel: " << channels.size() 
    << ", Number of hit: " << m_nHit << endreq;

  int index = -1;
  for(ci = channels.begin(); ci!=ci_end; ci++) {

    const ReadoutPmtChannel &channel = ci->second;

    int board = channel.channelId().board(); 
    int connector = channel.channelId().connector();
    int ring = -1;
    int clmn = -1;
    int wall = -1;
    int spot = -1;
    bool inout = false;
    if( readout->detector().isAD() ) {
      AdPmtSensor pmtId = m_cableSvc->adPmtSensor( channel.channelId(), svcMode );
      ring = pmtId.ring();
      clmn = pmtId.column();
    } else if( readout->detector().isWaterShield() ) {
      PoolPmtSensor pmtId = m_cableSvc->poolPmtSensor( channel.channelId(), svcMode );
      wall = pmtId.wallNumber();
      spot = pmtId.wallSpot();
      inout = pmtId.inwardFacing();
    }
    debug() << "(board,connector) = (" << board << "," << connector << "), " 
      << "(ring,column) = (" << ring << "," << clmn << ")" << endreq;

    unsigned int hitSize = channel.size();
    for(unsigned int i = 0; i<hitSize; i++) {
      index++;
      m_board[index] = board;
      m_channel[index] = connector;
      m_ring[index] = ring;
      m_column[index] = clmn;
      m_wall[index] = wall;
      m_spot[index] = spot;
      m_inout[index] = inout ? 1 : 0;
      m_adc[index] = channel.adc(i);
      m_tdc[index] = channel.tdc(i);
      m_cycle[index] = channel.adcCycle(i);
      m_preAdc[index] = channel.pedestal(i);
      m_range[index] = channel.adcRange(i);
      m_hitCount[index] = channel.tdcHitCount(i);
      m_multiTdc[index] = hitSize;
      debug() << "index=" << index << ", tdc=" << m_tdc[index] << ", adc=" << m_adc[index] << endreq;
    }
  } // end of channels

  m_ntuple1->write();
}

---