In This Package:

FlasherIdTool.cc

Go to the documentation of this file.

#include "FlasherIdTool.h"

#include "StatisticsSvc/IStatisticsSvc.h"
#include "Event/ReadoutHeader.h"
#include "Event/DaqCrate.h"
#include "Event/DaqPmtCrate.h"
#include "Event/DaqPmtChannel.h"
#include "Conventions/Electronics.h"

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

#include "TH2.h"
#include "TMath.h"
#include "TList.h"
#include "TTree.h"

#include <sstream>


char timeFormat[] = "#splitline{%H:%M:%S}{%b %d %Y}";

// Another frustrating feature of ROOT: timezones are a pain.
// It doesn't matter if you set the china local time here,
// it will always be drawn in the local time of the
// machine which draws the histogram. :(

FlasherIdTool::FlasherIdTool(const std::string& type,
                             const std::string& name, 
                             const IInterface* parent) 
  : GaudiTool(type,name,parent),
    m_statsSvc(0),
    m_cableSvc(0),
    m_pmtInfoSvc(0),
    m_firstTriggerTime(0),
    m_detector(Site::kSAB, DetectorId::kAD1),
    m_nEvents(0),
    m_maxTdcCalibEvents(1000),
    m_topCalibTdcOffset(0),
    m_topCalibNTdcs(0),
    m_bottomCalibTdcOffset(0),
    m_bottomCalibNTdcs(0),
    m_nChannels(0),
    m_minDt(-MAXDT/2.)
{
  declareInterface<IReadoutProcessor>(this);
  declareProperty("FlasherThreshold", m_flasherThreshold=0.4, 
                  "Threshold for flasher discrimination");
  declareProperty("Detector", m_detectorName="SABAD1", 
                  "Detector to process with this tool");

  for(unsigned int board=0; board<NBOARDS; board++){
    m_boardTdcOffset[board]=0;
    m_boardNTdcs[board]=0;
  }
  unsigned int maxChanIdx = NBOARDS*NCONNECTORS;
  for(unsigned int chanIdx=0; chanIdx<maxChanIdx; chanIdx++){
    m_flashBoard[chanIdx]=0;
    m_flashConnector[chanIdx]=0;
    m_flashDeltaAdc[chanIdx]=0;
    m_flashTdc[chanIdx]=0;
    for(unsigned int chanIdx2=0; chanIdx2<maxChanIdx; chanIdx2++){
      m_distance[chanIdx][chanIdx2]=0;
    }
  }

  // Pre-calculate log-likelihood for dt
  // Use distribution that roughly approximates data
  double levyScale = 50.;
  double bgScale = 0.01;
  double probability[MAXDT];
  double maxProb=0;
  double sqrt2pi = TMath::Sqrt(2*TMath::Pi());
  for(unsigned int idx=0; idx < MAXDT; idx++){
    double levyX0 = 25.;
    double dt = (m_minDt + idx + levyX0)/levyScale;
    probability[idx] = bgScale;
    if(dt>0){
      probability[idx] += (1-bgScale)*(TMath::Exp(-0.5/dt)
                                       /(sqrt2pi*TMath::Power(dt,1.5)));
    }
    if(probability[idx]>maxProb) maxProb = probability[idx];
  }
  for(unsigned int idx=0; idx < MAXDT; idx++){
    m_logLikelihood[idx] = -2*TMath::Log( probability[idx]/maxProb );
    debug() << m_logLikelihood[idx] << endreq;
  }
}

FlasherIdTool::~FlasherIdTool()
{
}

StatusCode FlasherIdTool::initialize()
{
  // Initialize the tool
  m_detector = DayaBay::Detector(m_detectorName);

  // Initialize the necessary services
  StatusCode sc = this->service("StatisticsSvc",m_statsSvc,true);
  if(sc.isFailure()){
    error() << "Failed to get StatisticsSvc" << endreq;
    return sc;
  }

  sc = this->service("StaticCableSvc",m_cableSvc,true);
  if(sc.isFailure()){
    error() << "Failed to get CableSvc" << endreq;
    return sc;
  }

  sc = this->service("PmtGeomInfoSvc",m_pmtInfoSvc,true);
  if(sc.isFailure()){
    error() << "Failed to get PmtGeomInfoSvc" << endreq;
    return sc;
  }

  return sc;
}

StatusCode FlasherIdTool::process(const DayaBay::ReadoutHeader* readoutHeader)
{
  // Check this readout for possible PMT flashers
  //
  const DayaBay::DaqCrate* daqCrate = readoutHeader->daqCrate();
  if(!daqCrate){
    error() << "Failed to get daq readout from header" << endreq;
    return StatusCode::FAILURE;
  }

  if(!daqCrate->detector().isAD() || daqCrate->detector()!=m_detector){
    // Not the detector we want, continue
    return StatusCode::SUCCESS;
  }

  // Convert to PMT crate readout
  const DayaBay::DaqPmtCrate* pmtCrate = daqCrate->asPmtCrate();

  if(!m_firstTriggerTime){
    m_firstTriggerTime = new TimeStamp(daqCrate->triggerTime());
    StatusCode sc = this->initializeDistanceTable(readoutHeader->context());
    if(!sc.isSuccess()) return sc;
  }

  // Find histograms
  int runNumber = daqCrate->runNumber();
  const DayaBay::Detector& detector = daqCrate->detector();

  TH2F* flasherRingVsColumn = 0;

  const std::vector<DayaBay::DaqPmtChannel*>& channels = pmtCrate->channelReadouts();
  
  std::vector<DayaBay::DaqPmtChannel*>::const_iterator channelIter, 
    channelEnd = channels.end();

  TimeStamp dtTriggerTime = daqCrate->triggerTime();
  dtTriggerTime.Subtract(m_lastTriggerTime);
  m_dt = dtTriggerTime.GetSeconds();

  // Stage 1: Collect TDC offset data for each board
  //  2" calib PMTs need an additional correction due to top/bottom cable length
  if( m_nEvents < m_maxTdcCalibEvents ){
    for(channelIter = channels.begin();
        channelIter!=channelEnd; 
        channelIter++){ 
      const DayaBay::DaqPmtChannel* channel = *channelIter;      
      if( channel->hitCount() < 1 ) continue;
      unsigned int board = channel->channelId().board();
      if( board!=17 ){
        m_boardTdcOffset[board] += channel->tdc(0);
        m_boardNTdcs[board]++;
      }else{
        // Deal with 2" calib PMT
        unsigned int connector = channel->channelId().connector();
        if(connector % 2 == 0){
          // Even PMTs on bottom
          m_bottomCalibTdcOffset += channel->tdc(0);
          m_bottomCalibNTdcs++;
        }else{
          // Odd PMTs on top
          m_topCalibTdcOffset += channel->tdc(0);
          m_topCalibNTdcs++;
        } 
      } 
    }
    m_nEvents++;
    return StatusCode::SUCCESS;
  }

  // Stage 2: Calculate current board TDC offset
  if(m_nEvents == m_maxTdcCalibEvents){
    // Calculate current TDC offsets
    for(unsigned int boardIdx = 0; boardIdx<NBOARDS; boardIdx++){
      if(m_boardNTdcs[boardIdx]>0){
        m_boardTdcOffset[boardIdx] /= m_boardNTdcs[boardIdx];
      }
    }
    if(m_topCalibNTdcs>0){
      m_topCalibTdcOffset /= m_topCalibNTdcs;
    }
    if(m_bottomCalibNTdcs>0){
      m_bottomCalibTdcOffset /= m_bottomCalibNTdcs;
    }
  }

  // Stage 3: Find first hit on each channel
  unsigned int nChan = 0;
  m_adcSum=0;
  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
    const DayaBay::DaqPmtChannel* channel = *channelIter;
    if( channel->hitCount() <1 ){
      continue;
    }   
    int board = channel->channelId().board();
    int connector = channel->channelId().connector();

    // Loop over hits, find first hit on channel
    unsigned int earliestTdc = channel->tdc( 0 );
    unsigned int adc = channel->adc( 0 );
    float preAdc = channel->preAdcAvg( 0 );
    bool isHighGainAdc = channel->isHighGainAdc( 0 );
    for(unsigned int hitIdx=1; hitIdx<channel->hitCount(); hitIdx++){
      if( channel->tdc( hitIdx ) > earliestTdc ){
        earliestTdc = channel->tdc( hitIdx );
        adc = channel->adc( hitIdx );
        preAdc = channel->preAdcAvg( hitIdx );
        isHighGainAdc = channel->isHighGainAdc( hitIdx );
      }
    }
    unsigned int chanIdx = this->channelIndex(board,connector);
    m_flashTdc[nChan] = earliestTdc - this->tdcOffset(board, connector);
    m_flashDeltaAdc[nChan] = adc-preAdc;
    if(!isHighGainAdc){
      m_flashDeltaAdc[nChan]*=19;
    }
    m_flashBoard[nChan] = board;
    m_flashConnector[nChan] = connector;
    m_flashRing[nChan] = m_ring[chanIdx];
    m_flashColumn[nChan] = m_column[chanIdx];
    m_adcSum += m_flashDeltaAdc[nChan];
    nChan++;
  }

  //Stage 4: calculate flasher discriminator for each PMT that was hit
  double tdcToNs = -(25./16.);
  double sqrt2 = sqrt(2.);
  for( unsigned int chanCount=0; chanCount<nChan; chanCount++ ){
    unsigned int flashBoard = m_flashBoard[chanCount];
    unsigned int flashConnector = m_flashConnector[chanCount];
    unsigned int flashChanIdx = this->channelIndex(flashBoard, flashConnector);
    // First check Q-discriminator
    double flashQDisc = 1 - m_flashDeltaAdc[chanCount]/m_adcSum;
    if(flashQDisc > m_flasherThreshold){
      // Skip this channel if Q-discriminator is already too high
      continue;
    }
    double flashDisc = 0;
    double flashTdc = m_flashTdc[chanCount];
    int nHit = 0;
    for(channelIter = channels.begin(); channelIter!=channelEnd;
        channelIter++) { 
      const DayaBay::DaqPmtChannel* channel = *channelIter;
      if( channel->hitCount() < 1 ) continue;
      int board = channel->channelId().board();
      int connector = channel->channelId().connector();
      unsigned int targetChanIdx = this->channelIndex(board, connector);
      if(targetChanIdx==flashChanIdx) continue;
      double tdcOff = this->tdcOffset(board, connector);
      for(unsigned int hitIdx=0; hitIdx<channel->hitCount(); hitIdx++){
        double tdcCorr = channel->tdc(hitIdx) - tdcOff;
        double dtHit = (tdcCorr - flashTdc) * tdcToNs;
        double timeResidual = (dtHit - m_distance[flashChanIdx][targetChanIdx]) 
          - this->transitTime(flashBoard);
        // Use integer ns resolution on likelihood to improve speed
        flashDisc += m_logLikelihood[int(timeResidual-m_minDt)];
        nHit++;
      }
    }
    if(nHit>0){
      flashDisc /= (nHit*m_logLikelihood[0]);
    }
    // Record flasher likelihood
    m_flashLikelihood[chanCount] = sqrt(flashDisc*flashDisc 
                                        + flashQDisc*flashQDisc)/sqrt2;
    if(m_flashLikelihood[chanCount] < m_flasherThreshold){
      TH2F* flasherQVsT = 
        dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                detector,
                                                flashBoard,
                                                flashConnector, 
                                                FLASHERQVST));
      flasherQVsT->Fill(flashDisc,flashQDisc);
      if(!flasherRingVsColumn){
        flasherRingVsColumn = 
          dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                  detector,
                                                  0, 0, 
                                                  FLASHERRINGVSCOLUMN));
      }
      flasherRingVsColumn->Fill(m_flashColumn[chanCount],
                                m_flashRing[chanCount]);
    }
  }

  m_nChannels = nChan;

  m_lastTriggerTime = daqCrate->triggerTime();

  m_nEvents++;
  return StatusCode::SUCCESS;
}

StatusCode FlasherIdTool::finalize()
{
  // Clean-up histogram shortcuts
  std::map<int,TH1**>::iterator histIter, histEnd = m_shortCuts.end();
  for(histIter = m_shortCuts.begin(); histIter!=histEnd; histIter++){
    delete [] (histIter->second);
    histIter->second = 0;
  }
  if( m_statsSvc ) m_statsSvc->release();
  if(m_firstTriggerTime){
    delete m_firstTriggerTime;
    m_firstTriggerTime = 0;
  }
  return StatusCode::SUCCESS;
}


StatusCode FlasherIdTool::initializeDistanceTable(const Context& context)
{
  // Pre-calculate PMT-to-PMT distances
  info() << "Init Dist A" << endreq;
  double lightspeed = 300.; // mm/ns
  ServiceMode svcMode(context,0);
  const std::vector<DayaBay::AdPmtSensor>& adPmts = 
    m_cableSvc->adPmtSensors(svcMode);
  std::vector<DayaBay::AdPmtSensor>::const_iterator pmtIdIter1, pmtIdIter2, 
    pmtIdEnd = adPmts.end();
  for(pmtIdIter1 = adPmts.begin(); pmtIdIter1!=pmtIdEnd; pmtIdIter1++){
    info() << "Init Dist b" << endreq;
    const DayaBay::AdPmtSensor& pmtId1 = *pmtIdIter1;
    info() << "  ring/column: " << pmtId1.ring() << "  " 
           << pmtId1.column() << endreq;
    if(pmtId1.ring()==0 && pmtId1.column()>6) continue;
    const DayaBay::FeeChannelId& channel1 = m_cableSvc->feeChannelId(pmtId1, 
                                                                     svcMode);
    unsigned int chanIdx1 = this->channelIndex(channel1.board(),
                                               channel1.connector());
    m_ring[chanIdx1] = pmtId1.ring();
    m_column[chanIdx1] = pmtId1.column();
    const CLHEP::Hep3Vector& pmtPos1 
      = m_pmtInfoSvc->get(pmtId1.fullPackedData())->localPosition();    
    for(pmtIdIter2 = adPmts.begin(); pmtIdIter2!=pmtIdEnd; pmtIdIter2++){
      //info() << "Init Dist c" << endreq;
      const DayaBay::AdPmtSensor& pmtId2 = *pmtIdIter2;
      //info() << "  ring/column: " << pmtId2.ring() << "  " 
      //     << pmtId2.column() << endreq;
      if(pmtId2.ring()==0 && pmtId2.column()>6) continue;
      const DayaBay::FeeChannelId& channel2 = m_cableSvc->feeChannelId(pmtId2,
                                                                       svcMode);
      //info() << "  board/connector: " << channel2.board() << "  " 
      //             << channel2.connector() << endreq;
      unsigned int chanIdx2 = this->channelIndex(channel2.board(),
                                                 channel2.connector());
      const CLHEP::Hep3Vector& pmtPos2 
        = m_pmtInfoSvc->get(pmtId2.fullPackedData())->localPosition();
      double distance = 0;
      double dr = pmtPos1.x()-pmtPos2.x();
      distance += dr*dr;
      dr = pmtPos1.y()-pmtPos2.y();
      distance += dr*dr;
      dr = pmtPos1.z()-pmtPos2.z();
      distance += dr*dr;
      distance = TMath::Sqrt(distance);
      //info() << "  distance " << distance << endreq;
      m_distance[chanIdx1][chanIdx2] = distance / lightspeed;
      //info() << "  done " << endreq;
    }
    //info() << "Init Dist d" << endreq;
  }
  return StatusCode::SUCCESS;
}

TH1* FlasherIdTool::getOrMakeHist(int run, const DayaBay::Detector& detector,
                                  int board, int connector, int histogram)
{
  std::map<int,TH1**>::iterator histIter = m_shortCuts.find(run);
  TH1** histograms = 0;
  if( histIter == m_shortCuts.end() ){
    // Initialize histogram shortcuts
    histograms = new TH1*[MAXHISTS];
    for(int i=0; i<MAXHISTS; i++) histograms[i] = 0;
    m_shortCuts[run] = histograms;
  }else{
    // Found run
    histograms = histIter->second;
  }
  
  TH1* hist = 
    histograms[board * NCONNECTORS * NHISTOGRAMS
               + connector * NHISTOGRAMS
               + histogram];
  if(!hist){
    // Make this histogram
    std::string histName;
    if(detector.detectorId()){
      if(board){
        // Make channel histograms
        switch(histogram){
        case FLASHERQVST:
          // Plot flasher charge vs. time discriminator
          histName = "flasherQVsT";
          hist = new TH2F(histName.c_str(),
                          "Flasher Q vs. T Discriminator",
                          200,0,1,200,0,1);
          hist->GetXaxis()->SetTitle("T-discriminator");
          hist->GetYaxis()->SetTitle("Q-discriminator");
          break;
        default:
          error() << "Unknown Histogram: " << histogram << endreq;
          return 0;
        }
      }else{
        // Make detector histogram
        switch(histogram){
        case FLASHERRINGVSCOLUMN:
          // Plot flasher ring vs. column
          histName = "flasherRingVsColumn";
          hist = new TH2F(histName.c_str(),
                          "Flasher Ring vs. Column",
                          49,0.25,24.75,19,-0.75,8.75);
          hist->GetXaxis()->SetTitle("PMT Column");
          hist->GetYaxis()->SetTitle("PMT Ring");
          hist->SetStats(0);              
          break;
        default:
          error() << "Unknown Histogram: " << histogram << endreq;
          return 0;
        }
      }
    }
    info() << "Making histogram: " << histName << endreq;
    m_statsSvc->put( this->getPath(run, detector.detName().c_str(), board, 
                                   connector, histName.c_str()), 
                     hist );
    histograms[board * NCONNECTORS * NHISTOGRAMS
               + connector * NHISTOGRAMS
               + histogram] = hist;
  }

  return hist;
}

unsigned int FlasherIdTool::channelIndex(unsigned int board, 
                                         unsigned int connector)
{
  // Return a counting index for channels
  return board * NCONNECTORS + connector;
}

double FlasherIdTool::tdcOffset(int board, int connector)
{
  double tdcOffset = m_boardTdcOffset[board];
  if(board==17){
    if(connector %2 == 0){
      tdcOffset = m_bottomCalibTdcOffset;
    }else{
      tdcOffset = m_topCalibTdcOffset;
    }
  }
  return tdcOffset;
}

double FlasherIdTool::transitTime(unsigned int board) 
{
  // Return transit time
  if(board==17) return 29.;
  return 55.;
}

std::string FlasherIdTool::getPath(int run, const char* detectorName,
                                   int board, int connector, 
                                   const char* histName)
{
  // Construct histogram path in statistics service
  std::stringstream path;
  path << "/file1/diagnostics/run_" << std::setfill('0') << std::setw(7) 
       << run;
  if(detectorName){
    path << "/detector_" << detectorName;
    if(board){
      path << "/channel_board" << std::setfill('0') << std::setw(2) 
           << board << "_connector" << std::setfill('0') << std::setw(2)
           << connector;
    }
  }
  path << "/" << histName;
  return path.str();
}

---