NeutronSpTightTag.cc

Go to the documentation of this file.

/*****************************
A tight spallation neutron tag, for energy calibration, by FeiHong.
zhangfh@ihep.ac.cn 2011-07-18
*****************************/

#include <math.h>
#include "NeutronSpTightTag.h"
#include "Event/CalibReadoutHeader.h"
#include "Event/RecHeader.h"
#include "Event/JobInfo.h"
#include "Event/UserDataHeader.h"
#include "Conventions/DetectorId.h"
#include "Conventions/JobId.h"
#include "DataSvc/IJobInfoSvc.h"
#include "DataUtilities/DybArchiveList.h"
#include "GaudiKernel/SystemOfUnits.h"



using namespace std;
using namespace DayaBay;


NeutronSpTightTag::NeutronSpTightTag(const string& name, ISvcLocator* svcloc)
: GaudiAlgorithm(name, svcloc)
{
  declareProperty("location", m_recHeaderLoc="/Event/Rec/AdSimple", 
                  "the TES path of target event header location");
  declareProperty("NeutronSpTag_location", m_tagloc="/Event/Tag/Physics/NeutronSpTight", 
                  "tagPath: path for spallation neutron tag to put");
  declareProperty("first_muon_location", m_muonloc="/Event/Tag/Muon/MuonAll",
                  "path for first muon trigger tag");
  declareProperty("low_energy_cut", m_eCutLow=6*Gaudi::Units::MeV, 
                  "low energy cut, units MeV");
  declareProperty("high_energy_cut", m_eCutHigh=10*Gaudi::Units::MeV, 
                  "high energy cut, units MeV");
  declareProperty("zcut", m_zCut=1500*Gaudi::Units::millimeter, 
                  "vertex cut on z direction, units mm");
  declareProperty("rcut", m_rCut=1500*Gaudi::Units::millimeter, 
                  "vertex cut on r direction, units mm");
  declareProperty("time_window_start", m_timeCutLow=15*Gaudi::Units::microsecond, 
                  "time window start time, units micro second");
  declareProperty("time_window_end", m_timeCutHigh=100*Gaudi::Units::microsecond, 
                  "time window end time, units micro second");
}

NeutronSpTightTag::~NeutronSpTightTag()
{
}

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

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

  //Get Archive Svc
  StatusCode status = service("EventDataArchiveSvc", m_archiveSvc);
  if (status.isFailure()) {
    Error("Service [EventDataArchiveSvc] not found", status);
  }

  m_count=0;

  return StatusCode::SUCCESS;
}

StatusCode NeutronSpTightTag::execute()
{
  debug() << "execute() ______________________________ start" << endreq;
  
  // Get AdSimple
  RecHeader* recHeader = get<RecHeader>(m_recHeaderLoc);
  if(!recHeader) {
    error() << "Get no recHeader from TES!!" << endreq;
    return StatusCode::FAILURE;
  }

  DetectorId::DetectorId_t detectorId = recHeader->context().GetDetId();
  //Only in AD
  if(detectorId == DetectorId::kAD1 || detectorId == DetectorId::kAD2
      || detectorId == DetectorId::kAD3 || detectorId == DetectorId::kAD4)
    {
      double E = recHeader->recTrigger().energy()*Gaudi::Units::MeV;
      double Z = fabs(recHeader->recTrigger().position().z())*Gaudi::Units::millimeter;
      double R = sqrt((recHeader->recTrigger().position().x())*(recHeader->recTrigger().position().x())+(recHeader->recTrigger().position().y())*(recHeader->recTrigger().position().y()))*Gaudi::Units::millimeter;
      debug() << "Energy: " << E/Gaudi::Units::MeV << endreq;
      debug() << "Z: " << Z/Gaudi::Units::millimeter << endreq;
      debug() << "R: " << R/Gaudi::Units::millimeter << endreq;

      //Energy Cut & Vertex Cut
      if(E>m_eCutLow && E<m_eCutHigh && Z<m_zCut && R<m_rCut)
        {
          //Get muon list
          
          SmartDataPtr<DybArchiveList>  muonlist(m_archiveSvc, m_muonloc);
          if(muonlist)
            {
              debug() << "Number of muon in the ArchiveList: " << muonlist->size() << endreq;
            }
          else
            {
              debug() << "No MuonList in AES" << endreq;
              return StatusCode::SUCCESS;
            }
          
          DybArchiveList::const_iterator iter = muonlist->begin();
          for(; iter!=muonlist->end(); iter++)
            {
              bool getMuon = false;
              HeaderObject* muonTag = dynamic_cast<HeaderObject*>(*iter);
              const std::vector<const DayaBay::IHeader*>& muonTagList = muonTag->inputHeaders();
              vector<const IHeader*>::const_iterator mit = muonTagList.begin();
              for(;mit!=muonTagList.end();mit++)
                {
                  HeaderObject* muonTagHeader = dynamic_cast<HeaderObject*>(const_cast<IHeader*>(*mit));
                  DetectorId::DetectorId_t muonDetId = muonTagHeader->context().GetDetId();
                  TimeStamp muonTime(muonTagHeader->timeStamp());
                  TimeStamp currentTime(recHeader->timeStamp());
                  currentTime.Subtract(muonTime);
                  double deltaT = currentTime.GetSeconds()*Gaudi::Units::second;
                  if(deltaT>m_timeCutLow && deltaT<m_timeCutHigh && muonDetId == detectorId)
                    {
                      //Make a new Header as a tag
                      HeaderObject *SpalNTightTag = new HeaderObject();
                      // Use inputHeaders to record the RecHeader which is connect to this tag
                      vector<const IHeader*> inputHeaders;
                      inputHeaders.push_back(recHeader);
                      SpalNTightTag->setInputHeaders(inputHeaders);
                      // Copy from the recHeader
                      SpalNTightTag->setExecNumber(recHeader->execNumber());
                      SpalNTightTag->setContext(recHeader->context());
                      SpalNTightTag->setEarliest(recHeader->earliest());
                      SpalNTightTag->setLatest(recHeader->latest());
                      // Get job id
                      const DayaBay::JobId &m_currentJobId = m_jobInfoSvc->currentJobInfo()->jobId();
                      SpalNTightTag->setJobId(m_currentJobId);
                      
                      put(SpalNTightTag, m_tagloc);
                      m_count++;
                      debug() << "N Tagged: " << m_count<< endreq;
                      debug() << "DeltaT: " << deltaT/Gaudi::Units::second << endreq;
                      getMuon = true;
                      break;
                    }
                }
              if(getMuon==true) break;
            }
        }
    }
  debug() << "execute() ______________________________ end" << endreq;
  return StatusCode::SUCCESS;
}

StatusCode NeutronSpTightTag::finalize()
{
  debug() << "finalize()" << endreq;
  info() << "Summary: Tagger NeutronSpLooseTag : "<< m_count <<" events tagged at "<< m_tagloc << endreq;
  return StatusCode::SUCCESS;
}