Co60Tag.cc

Go to the documentation of this file.

#include "Co60Tag.h"
#include "Event/ReadoutHeader.h"
#include "Event/CalibReadoutHeader.h"
#include "Event/RecHeader.h"
#include "Event/HeaderObject.h"
#include "Event/UserDataHeader.h"
#include "Event/JobInfo.h"
#include "Conventions/JobId.h"
#include "DataSvc/IJobInfoSvc.h"
#include "GaudiKernel/SystemOfUnits.h"


Co60Tag::Co60Tag(const std::string& name, ISvcLocator* svcloc)
: GaudiAlgorithm(name, svcloc)
{
  declareProperty("CalibStatsLocation", m_calibStatsLocation="/Event/Data/CalibStats",
                  "Location of CalibStats User data");
  declareProperty("RecHeaderLocation", m_recHeaderLocation="/Event/Rec/AdSimple",
                  "Location of RecHeader");
  declareProperty("MuonLooseLocation", m_muonLooseLocation="/Event/Tag/Physics/MuonLoose",
                  "Location of Muon Loose Tags");
  declareProperty("Co60Location", m_Co60Location="/Event/Tag/Calib/Co60",
                  "Location of Co60 Tag");

  declareProperty("VertexCut", m_vertexCut = 0.5 * Gaudi::Units::m,
                  "Radius of the spherical vertex cut");
  declareProperty("MuonVetoTime", m_muonVetoTime = 1.0 * Gaudi::Units::ms,
                  "Duration of muon veto after a \" loose muon\".");
}

Co60Tag::~Co60Tag()
{
}

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

  m_jobInfoSvc = svc<IJobInfoSvc>("JobInfoSvc",true);
  if(!m_jobInfoSvc) {
    error() << "Failed to initialize JobInfoSvc" << endreq;
    return StatusCode::FAILURE;
  }
  m_tagged = 0;
  m_freshStart = true;

  return StatusCode::SUCCESS;
}

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

  //----------------------------------------
  // Get all necessary headers
  //----------------------------------------

  const DayaBay::RecHeader * recHdr = 0;  
  DayaBay::UserDataHeader *calibStatsHdr = 0;  //
  
  // Get CalibStats Header
  if (!exist<DayaBay::UserDataHeader>(evtSvc(), m_calibStatsLocation)) {
    warning() << "Cannot find header at " << m_calibStatsLocation << endreq;
    return StatusCode::FAILURE;
  } else {
    calibStatsHdr = get<DayaBay::UserDataHeader>(m_calibStatsLocation);
  }

  // Get Rec Header
  if ( !exist<DayaBay::RecHeader>(evtSvc(), m_recHeaderLocation)) {
    warning() << "Cannot find header at " << m_recHeaderLocation << endreq;
    return StatusCode::FAILURE;
  } else {
    recHdr = get<DayaBay::RecHeader>(m_recHeaderLocation);
  } 
  const DayaBay::RecTrigger& recTrigger = recHdr->recTrigger();
       
  DetectorId::DetectorId_t detectorId = recHdr->context().GetDetId(); 
  TimeStamp triggerTime = recTrigger.triggerTime();

  //----------------------------------------
  // Actual checking
  //----------------------------------------

  // Fresh start! Assume it is muon event and skip.
  if(m_freshStart){
    debug() << "Fresh start! Assume it is muon event and skip." << endreq;
    m_lastMuonTime = triggerTime;
    m_freshStart = false;
    return StatusCode::SUCCESS;
  }

  // Check if it is a muon, if so, start veto and skip.
  bool isMuon = exist<DayaBay::HeaderObject>(evtSvc(), m_muonLooseLocation);
  if (isMuon){
    debug() << "Muon found! Start veto and skip." << endreq;
    m_lastMuonTime = triggerTime;
    return StatusCode::SUCCESS;
  }
  
  // Only consider AD triggers
  if (!(detectorId == DetectorId::kAD1 || detectorId == DetectorId::kAD2
   || detectorId == DetectorId::kAD3 || detectorId == DetectorId::kAD4)) {
    debug() << "Non-AD trigger. Skip." << endreq;
    return StatusCode::SUCCESS;
  }

  // Apply muon veto
  TimeStamp timeSinceLastMuon = TimeStamp(triggerTime);
  timeSinceLastMuon.Subtract(m_lastMuonTime);
  debug() << "Current trigger time: " << triggerTime.AsString() << endreq;
  debug() << "Last muon time      : " << m_lastMuonTime.AsString() << endreq;
  debug() << "Time since last muon: " << timeSinceLastMuon.GetSeconds() << endreq;
  if (timeSinceLastMuon.GetSeconds() * Gaudi::Units::s < m_muonVetoTime){
    //Muon veto on!
    debug() << "Muon veto on! Skip." << endreq;
    return StatusCode::SUCCESS;
  }

  // Skip if recon has failed
  float recES = recTrigger.energyStatus();
  float recPS = recTrigger.positionStatus();
  if (recES != 1  || recPS != 1){ 
    debug() << "Recon failed. Skip." << endreq;
    return StatusCode::SUCCESS;
  }

  // Apply spherical position cut  
  // FIXME: Source position is hard-coded to be AD center!
  float recX = recTrigger.position().x() / Gaudi::Units::mm;
  float recY = recTrigger.position().y() / Gaudi::Units::mm;
  float recZ = recTrigger.position().z() / Gaudi::Units::mm;
  float posX = 0;
  float posY = 0;
  float posZ = 0;
  float delX = recX - posX;
  float delY = recY - posY;
  float delZ = recZ - posZ;
  float delR2 = delX*delX + delY*delY + delZ*delZ;
  debug() << "Source position: ("<< posX << ", " << posY << ", " << posZ << ")" << endreq;
  debug() << "Recon position : ("<< recX << ", " << recY << ", " << recZ << ")" << endreq;
  if (delR2 > m_vertexCut * m_vertexCut){
    debug() << "Recon position lies outside of the vertex cut! Skip." << endreq;
    return StatusCode::SUCCESS;
  }

  //----------------------------------------
  // Got a candidate, let's tag it!
  //----------------------------------------

  // Tag as Co60 event
  DayaBay::HeaderObject *header = new DayaBay::HeaderObject();
  std::vector<const DayaBay::IHeader*> inputHeaders;
  inputHeaders.push_back(recHdr);
  inputHeaders.push_back(calibStatsHdr);
  header->setInputHeaders(inputHeaders);
  header->setExecNumber(recHdr->execNumber());
  header->setContext(recHdr->context());
  header->setEarliest(recHdr->earliest());
  header->setLatest(recHdr->latest());
  m_tagged++;

  const DayaBay::JobId &m_currentJobId = m_jobInfoSvc->currentJobInfo()->jobId();
  header->setJobId(m_currentJobId);

  put(header, m_Co60Location);

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

StatusCode Co60Tag::finalize()
{
  debug() << "finalize()" << endreq;
  switch(m_tagged){
    case 0: info () << "No trigger was tagged as Co60 event." << endreq; break;
    case 1: info () << "Only one trigger was tagged as Co60 event." << endreq; break;
    default: info () << m_tagged << " triggers were tagged as Co60 events." << endreq;
  }
  return StatusCode::SUCCESS;
}