ADCoincFilter.cc

Go to the documentation of this file.

//*****************************
//
// AD Coincidence Tight for Filter
// chao@bnl.gov  2011-07-29
//
//*****************************

#include "ADCoincFilter.h"

#include "Event/HeaderObject.h"
#include "Event/RecHeader.h"
#include "Event/ReadoutHeader.h"
#include "DataSvc/IJobInfoSvc.h"
#include "DataUtilities/DybArchiveList.h"
#include "GaudiKernel/SystemOfUnits.h"

using namespace std;
using namespace Gaudi;
using namespace DayaBay;

ADCoincFilter::ADCoincFilter(const string& name, ISvcLocator* svcloc)
    : GaudiAlgorithm(name, svcloc),
    m_nTagged(0), m_nSeen(0), m_nKept(0), m_unPhysicalRun(false)
{
    declareProperty("RecLocation", m_recHeaderLoc="/Event/Rec/AdSimple",
            "the TES path of rec header location");
    declareProperty("TagADLocation", m_tagADLoc="/Event/Tag/Physics/CoincidenceTight",
            "the TES path of AD Coincidenc tag header location");
    declareProperty("TagMuonLocation", m_tagMuonLoc="/Event/Tag/Physics/CoincidenceTightRelated",
            "the TES path of Correlated Muon tag header location");
    declareProperty("TimeWindow", m_timeWindow=300*Units::microsecond,
        "AD coincidence time window");
    declareProperty("MuonVetoWindow", m_muonVetoWindow=10*Units::microsecond,
        "muon veto time window");
    declareProperty("MuonKeptWindow", m_muonKeptWindow=300*Units::microsecond,
        "muon kept time window");
    declareProperty("MaxCoincFraction", m_maxCoincFraction=0.5,
        "maximum allowed coincident trigger fraction");
}


ADCoincFilter::~ADCoincFilter()
{
}


StatusCode ADCoincFilter::initialize()
{
    info() << "initialize: time window: " << m_timeWindow/Units::microsecond << " microseconds" << endreq;    
    info() << "initialize: veto window: " << m_muonVetoWindow/Units::microsecond << " microseconds" << endreq;    

    m_jobInfoSvc = svc<IJobInfoSvc>("JobInfoSvc", true);
    if(!m_jobInfoSvc) {
        error() << "Failed to initialize JobInfoSvc" << endreq;
        return StatusCode::FAILURE;
    }
    
    // init archive event store
    StatusCode sc = service("EventDataArchiveSvc", m_archiveSvc);
    if(sc.isFailure()){
      error() << "Failed to access the Archive event store" << endreq;
      return sc;
    }
    
    // initialize
    m_inputHeaders[DetectorId::kAD1] = vector<const IHeader*>();
    m_inputHeaders[DetectorId::kAD2] = vector<const IHeader*>();
    m_inputHeaders[DetectorId::kAD3] = vector<const IHeader*>();
    m_inputHeaders[DetectorId::kAD4] = vector<const IHeader*>();
    
    return StatusCode::SUCCESS;
}



StatusCode ADCoincFilter::execute()
{   
    if (m_unPhysicalRun) {
        return StatusCode::SUCCESS; // looks like not a physical run, skip.
    }
    
    RecHeader* recHeader = get<RecHeader>(m_recHeaderLoc);
    if ( !recHeader ) {
        warning() << "Cannot find header at " << m_recHeaderLoc << endreq;
        return StatusCode::FAILURE;
    } 
    
    TimeStamp ts = recHeader->context().GetTimeStamp();
    DetectorId::DetectorId_t detectorId = recHeader->context().GetDetId();
    
    TimeStamp lastTS = ts;
    
    if(m_lastTimeStamp.find(detectorId) != m_lastTimeStamp.end()){
        lastTS = m_lastTimeStamp[detectorId];
    }       
    m_lastTimeStamp[detectorId] = ts;
    
    // info() << detectorId << " " << ts.GetNanoSec() << endreq;
    
    // only look for AD Coincidence
    if (   detectorId == DetectorId::kIWS 
        || detectorId == DetectorId::kOWS 
        || detectorId == DetectorId::kRPC) {
        return StatusCode::SUCCESS; 
    }
    
    TimeStamp dt = TimeStamp(ts);
    dt.Subtract(lastTS);
    if (dt.GetSeconds()*Units::second > m_timeWindow) {
        int multiplicity = m_inputHeaders[detectorId].size();
        if (multiplicity < 2) {
            // found singles, clear the inputHeaders
            m_inputHeaders[detectorId].clear();
        }
        else {
            // found a gap, save last conincidence's multiplets         
            HeaderObject *header = new HeaderObject();
            vector<const IHeader*>& inputHeaders = m_inputHeaders[detectorId];
            header->setInputHeaders(inputHeaders);
            header->setExecNumber(recHeader->execNumber());
            header->setContext(dynamic_cast<const HeaderObject*>(inputHeaders.at(0))->context());
            header->setEarliest(dynamic_cast<const HeaderObject*>(inputHeaders.at(0))->earliest());
            header->setLatest(dynamic_cast<const HeaderObject*>(inputHeaders.at(inputHeaders.size()-1))->latest());
            header->setJobId(m_jobInfoSvc->currentJobInfo()->jobId());

            put(header, m_tagADLoc);
            m_nTagged++;
            m_nKept += multiplicity;
            // info() << "found one: multiplicty: " << multiplicity << endreq;
            
            
            // create correlated muon tags
            HeaderObject *header2 = new HeaderObject();
            vector<const IHeader*> muonInputHeaders = vector<const IHeader*>();
            DybArchiveList* arcList = get<DybArchiveList*>(m_archiveSvc, "/Event/Readout/ReadoutHeader");
            if(!arcList){
                error() << "Failed to get archive of ReadoutHeader" << endreq;
                return StatusCode::FAILURE;
            }
            TimeStamp firstCoincTS = dynamic_cast<const HeaderObject*>(inputHeaders.at(0))->timeStamp();
            TimeStamp lastCoincTS  = dynamic_cast<const HeaderObject*>(inputHeaders.at(inputHeaders.size()-1))->timeStamp();
            for(unsigned int headerIdx=0; headerIdx<arcList->size(); headerIdx++){
                  ReadoutHeader* arcHeader = dynamic_cast<ReadoutHeader*>(arcList->dataObject(headerIdx));
                  if(!arcHeader) {
                      error() << "Invalid ReadoutHeader in AES!" << endreq;
                      return StatusCode::FAILURE;
                  }
                  DetectorId::DetectorId_t detId = arcHeader->context().GetDetId();
                  if ( detId == DetectorId::kIWS || detId == DetectorId::kOWS || detId == DetectorId::kRPC) {
                      TimeStamp arcTS = arcHeader->timeStamp();
                      TimeStamp dtLastCoinc = TimeStamp(arcTS);
                      dtLastCoinc.Subtract(lastCoincTS);
                      TimeStamp dtFirstCoinc = TimeStamp(arcTS);
                      dtFirstCoinc.Subtract(firstCoincTS);
                      if (dtLastCoinc.GetSeconds() > 0) { 
                          continue; 
                      }
                      else if (dtLastCoinc.GetSeconds()*Units::second > -m_muonKeptWindow ) {
                          // info() << "muon saved: " << detId << " " << arcTS.GetNanoSec() << endreq;
                          muonInputHeaders.push_back(arcHeader);
                      }
                      else { 
                          break; 
                      }
                  }
            }
            header2->setInputHeaders(muonInputHeaders); // muonHeaders may be empty
            header2->setExecNumber(recHeader->execNumber());
            header2->setContext(recHeader->context());
            header2->setEarliest(recHeader->earliest());
            header2->setLatest(recHeader->latest());
            header2->setJobId(m_jobInfoSvc->currentJobInfo()->jobId());
            put(header2, m_tagMuonLoc);
            m_nKeptMuon += muonInputHeaders.size();
            
            // clear the inputHeaders
            m_inputHeaders[detectorId].clear();
        }
    }
    
    // muon veto
    TimeStamp dt_IWS = TimeStamp(ts);
    if(m_lastTimeStamp.find(DetectorId::kIWS) != m_lastTimeStamp.end()){
        dt_IWS.Subtract(m_lastTimeStamp[DetectorId::kIWS]);
    }
    TimeStamp dt_OWS = TimeStamp(ts);
    if(m_lastTimeStamp.find(DetectorId::kOWS) != m_lastTimeStamp.end()){
        dt_OWS.Subtract(m_lastTimeStamp[DetectorId::kOWS]);
    }
    if (!(dt_IWS*Units::second < m_muonVetoWindow && dt_OWS*Units::second < m_muonVetoWindow)) {
        m_inputHeaders[detectorId].push_back(recHeader);    
    }
    
    // catch unphysical runs
    m_nSeen++;
    if (float(m_nKept)/m_nSeen > m_maxCoincFraction && m_nSeen > 1000) {
        m_unPhysicalRun = true;
        warning() << "saw " << m_nKept << " coincident triggers out of " << m_nSeen
            << " total. This is unphysical, disabling IBD Filter Tag" << endreq;
    }

    return StatusCode::SUCCESS;
}


StatusCode ADCoincFilter::finalize()
{
    info() << m_nTagged << " multiplets (" << m_nKept << " triggers) tagged at " << m_tagADLoc 
           << " out of " << m_nSeen << " triggers" << endreq;
    info() << m_nKeptMuon << " muon triggers kept"<< endreq;

    return StatusCode::SUCCESS;
}