NeutronSpLooseData.cc

Go to the documentation of this file.

/*
 * NeutronSpLooseData
 *
 * Store useful spallation neutron data
 *
 * Created by: zhangfh@ihep.ac.cn 2011/08/04
 */

#include "NeutronSpLooseData.h"

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

#include "GaudiKernel/SystemOfUnits.h"
#include "TMath.h"
#include "TString.h"

using namespace std;
using namespace DayaBay;


NeutronSpLooseData::NeutronSpLooseData(const std::string& name, ISvcLocator* pSvcLocator)
  : DybAlgorithm<DayaBay::UserDataHeader>(name,pSvcLocator)
{
  declareProperty("NeutronSpTag_location", m_neutronSpTagLoc="/Event/Tag/Physics/NeutronSpLoose", 
                  "Path for spallation neutron tag");
  declareProperty("first_muon_location", m_muonloc="/Event/Tag/Muon/MuonAll",
                  "path for first muon trigger tag");
}

NeutronSpLooseData::~NeutronSpLooseData()
{
}

StatusCode NeutronSpLooseData::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);
    return status;
  }

  m_count=0;

  return StatusCode::SUCCESS;
}

StatusCode NeutronSpLooseData::execute()
{
  debug() << "execute() ______________________________ start" << endreq;
  
  // Get tagged spallation neutron recHeader
  if(!exist<HeaderObject>(m_neutronSpTagLoc)) {
    return StatusCode::SUCCESS;
  }

  HeaderObject *neutronSpTag = get<HeaderObject>(m_neutronSpTagLoc);
  if(!neutronSpTag) {
    error() << "Get no neutronSpTag from TES." << endreq;
    return StatusCode::FAILURE;
  }
  
  const std::vector<const DayaBay::IHeader*>& recHeaderList = neutronSpTag->inputHeaders();

  // For NeutronSpTag, only one recHeader in each inputHeaders.
  const RecHeader* recHeader = dynamic_cast<const RecHeader*> (recHeaderList[0]);
  if(!recHeader) {
    error() << "Get no recHeader from NeutronSpTag!!" << endreq;
    return StatusCode::FAILURE;
  }

  // Create output user data header
  UserDataHeader* neutronSpData = MakeHeaderObject();

  // Copy from the recHeader
  neutronSpData->setExecNumber(recHeader->execNumber());
  neutronSpData->setContext(recHeader->context());
  neutronSpData->setEarliest(recHeader->earliest());
  neutronSpData->setLatest(recHeader->latest());


  // Get job id
  const DayaBay::JobId &m_currentJobId = m_jobInfoSvc->currentJobInfo()->jobId();
  neutronSpData->setJobId(m_currentJobId);

  float energy=recHeader->recTrigger().energy();
  float x=recHeader->recTrigger().position().x();
  float y=recHeader->recTrigger().position().y();
  float z=recHeader->recTrigger().position().z();
  TimeStamp trigTime(recHeader->timeStamp());
  int trigTimeSec=trigTime.GetSec();
  int trigTimeNanoSec=trigTime.GetNanoSec();
  int detectorId = recHeader->context().GetDetId();

  int muonTimeSec=0;
  int muonTimeNanoSec=0;
  float deltaT = 0;
  float muonPESum = 0;

  float nPESum=0;
  float nPEMax=0;
  float tRMS=0;
  //Get muon list         
  SmartDataPtr<DybArchiveList>  muonlist(m_archiveSvc, m_muonloc);
  if(muonlist)
    {
      debug() << "Number of muon in the ArchiveList: " << muonlist->size() << endreq;
    }
  else
    {
      error() << "No MuonList for NeutronSpTag event?!" << endreq;
      return StatusCode::FAILURE;
    }
  DybArchiveList::const_iterator iter = muonlist->begin();
  for(; iter!=muonlist->end(); iter++)
    {
      HeaderObject* muonTag = dynamic_cast<HeaderObject*>(*iter);
      //const std::vector<const DayaBay::IHeader*>& muonTagList = muonTag->inputHeaders();
      //HeaderObject* muonTagHeader = dynamic_cast<HeaderObject*>(const_cast<IHeader*>(muonTagList[0]));
      //const std::vector<const DayaBay::IHeader*>& muonIHeaders = muonTagHeader->inputHeaders();
      //info() << "muonIHeaders=" << muonIHeaders.size() << endreq;
      //UserDataHeader* muonClbStsHdr = dynamic_cast<UserDataHeader*>(dynamic_cast<HeaderObject*>(const_cast<IHeader*>(muonIHeaders[1])));

      // Make it compatible with both python and C++ muon tag
      // Find CalibReadoutHeader
      const std::vector<const IHeader*> &crhs = muonTag->findHeaders(CalibReadoutHeader::classID());
      HeaderObject* muonTagHeader = dynamic_cast<HeaderObject*>(const_cast<IHeader*>(crhs[0]));

      /*
      for(unsigned int i = 0; i < crhs.size(); i++) {
        IHeader* iHeader = const_cast<IHeader*>(crhs[i]);
        HeaderObject* header = dynamic_cast<HeaderObject*>(iHeader);
        debug() << header->name() << ", " << header->context() << endreq;
        if(header->context().GetDetId() == detectorId) {
          muonTagHeader = header;
          break;
        }
      }
      */

      // Find CalibStatsHeader
      const std::vector<const IHeader*> &calibstats = muonTag->findHeaders(UserDataHeader::classID());
      UserDataHeader* muonClbStsHdr = dynamic_cast<UserDataHeader*>(const_cast<IHeader*>(calibstats[0]));

      /*
      for(unsigned int i = 0; i < calibstats.size(); i++) {
        IHeader* iHeader = const_cast<IHeader*>(calibstats[i]);
        UserDataHeader* header = dynamic_cast<UserDataHeader*>(iHeader);
        debug() << header->name() << ", " << header->context() << endreq;
        if(header->name()!="/CalibStats") continue;
        if(header->context().GetDetId() == detectorId) {
          muonClbStsHdr = header;
          break;
        }
      }
      */

      if(muonTagHeader) {
        TimeStamp muonTime(muonTag->timeStamp());
        muonTimeSec=muonTime.GetSec();
        muonTimeNanoSec=muonTime.GetNanoSec();
        TimeStamp currentTime(recHeader->timeStamp());
        currentTime.Subtract(muonTime);
        deltaT = currentTime.GetSeconds();
        debug() << "DeltaT: " << deltaT << endreq;
      }
      if(muonClbStsHdr) {
        muonPESum = muonClbStsHdr->getFloat("nPESum");
        debug() << "nPESum: " << muonPESum << endreq;
      }

      break;
    }

  debug()<<"energy:"<<energy<<";x:"<<x<<";y:"<<y<<";z:"<<z<<";timeToMuon:"<<deltaT<<";detectorId:"<<detectorId<<endreq;

  //Get some calibstats data
  const std::vector<const IHeader*> &calibstats = neutronSpTag->findHeaders(UserDataHeader::classID());
  UserDataHeader* spNeutronClbStsHdr = 0;
  for(unsigned int i = 0; i < calibstats.size(); i++) {
    IHeader* iHeader = const_cast<IHeader*>(calibstats[i]);
    UserDataHeader* header = dynamic_cast<UserDataHeader*>(iHeader);
    debug() << header->name() << ", " << header->context() << endreq;
    if(header->name()!="/CalibStats") continue;
    if(header->context().GetDetId() == detectorId) {
      spNeutronClbStsHdr = header;
      break;
    }
  }
  if(spNeutronClbStsHdr) {
    nPESum = spNeutronClbStsHdr->getFloat("nPESum");
    nPEMax = spNeutronClbStsHdr->getFloat("nPEMax");
    tRMS   = spNeutronClbStsHdr->getFloat("tRMS");
    debug() << "nPESum: " << nPESum <<";nPEMax: "<<nPEMax<<";tRMS: "<<tRMS<< endreq;
  }

  // Fill user data
  neutronSpData->set("energy", energy);
  neutronSpData->set("x", x);
  neutronSpData->set("y", y);
  neutronSpData->set("z", z);
  neutronSpData->set("trigTimeSec", trigTimeSec);
  neutronSpData->set("trigTimeNanoSec", trigTimeNanoSec);
  neutronSpData->set("timeToMuon", deltaT);
  neutronSpData->set("detectorId", detectorId);
  neutronSpData->set("muonTimeSec", muonTimeSec);
  neutronSpData->set("muonTimeNanoSec", muonTimeNanoSec);
  neutronSpData->set("muonPESum", muonPESum);
  neutronSpData->set("nPESum", nPESum);
  neutronSpData->set("nPEMax", nPEMax);
  neutronSpData->set("tRMS", tRMS);
  m_count++;
  debug() << "execute() ______________________________ end" << endreq;
  return StatusCode::SUCCESS;
}

StatusCode NeutronSpLooseData::finalize()
{
  debug() << "finalize()" << endreq;
  info() << "Summary: UserData NeutronSpLooseData : "<< m_count <<" events data filled."<< endreq;
  return StatusCode::SUCCESS;
}