algorithm to dump event infos to TTree format

Analysis/DumpEvent/CMakeLists.txt

+gaudi_add_module(DumpEvent
+                 SOURCES src/DumpMCParticleAlg.cpp
+		         #src/DumpHitAlg.cpp
+			 src/DumpTrackAlg.cpp
+			 #src/DumpCalorimeterAlg.cpp
+		 LINK DataHelperLib
+                      Gaudi::GaudiKernel
+                      EDM4HEP::edm4hep
+                      ${ROOT_LIBRARIES}
+                      ${CLHEP_LIBRARIES}
+		      ${DD4hep_COMPONENT_LIBRARIES}
+                      DetInterface
+)
+
+install(TARGETS DumpEvent
+  EXPORT CEPCSWTargets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  COMPONENT dev)

Analysis/DumpEvent/src/DumpMCParticleAlg.cpp

+#include "DumpMCParticleAlg.h"
+
+#include "GaudiKernel/DataObject.h"
+#include "GaudiKernel/IHistogramSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/SmartDataPtr.h"
+
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/HelixClass.h"
+
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <math.h>
+
+DECLARE_COMPONENT( DumpMCParticleAlg )
+
+//------------------------------------------------------------------------------
+DumpMCParticleAlg::DumpMCParticleAlg( const std::string& name, ISvcLocator* pSvcLocator )
+    : Algorithm( name, pSvcLocator ) {
+  declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
+
+  m_thisName = name;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpMCParticleAlg::initialize(){
+  info() << "Booking Ntuple" << endmsg;
+
+  NTuplePtr nt1(ntupleSvc(), "MyTuples/MC");
+  if ( !nt1 ) {
+    m_tuple = ntupleSvc()->book("MyTuples/MC",CLID_ColumnWiseTuple,"MC truth");
+    if ( 0 != m_tuple ) {
+      m_tuple->addItem ("nmc",        m_nParticles, 0, 1000 ).ignore();
+      m_tuple->addIndexedItem ("pdg",        m_nParticles, m_pdgID ).ignore();
+      m_tuple->addIndexedItem ("genStatus",  m_nParticles, m_genStatus ).ignore();
+      m_tuple->addIndexedItem ("simStatus",  m_nParticles, m_simStatus ).ignore();
+      m_tuple->addIndexedItem ("charge",     m_nParticles, m_charge ).ignore();
+      m_tuple->addIndexedItem ("time",       m_nParticles, m_time ).ignore();
+      m_tuple->addIndexedItem ("mass",       m_nParticles, m_mass ).ignore();
+      m_tuple->addIndexedItem ("vx",         m_nParticles, m_vx ).ignore();
+      m_tuple->addIndexedItem ("vy",         m_nParticles, m_vy ).ignore();
+      m_tuple->addIndexedItem ("vz",         m_nParticles, m_vz ).ignore();
+      m_tuple->addIndexedItem ("px",         m_nParticles, m_px ).ignore();
+      m_tuple->addIndexedItem ("py",         m_nParticles, m_py ).ignore();
+      m_tuple->addIndexedItem ("pz",         m_nParticles, m_pz ).ignore();
+      m_tuple->addIndexedItem ("d0",         m_nParticles, m_d0 ).ignore();
+      m_tuple->addIndexedItem ("phi0",       m_nParticles, m_phi0 ).ignore();
+      m_tuple->addIndexedItem ("omega",      m_nParticles, m_omega ).ignore();
+      m_tuple->addIndexedItem ("z0",         m_nParticles, m_z0 ).ignore();
+      m_tuple->addIndexedItem ("tanLambda",  m_nParticles, m_tanLambda ).ignore();
+    }
+    else { // did not manage to book the N tuple....
+      fatal() << "Cannot bool MyTuples/MC " << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+  else{
+    m_tuple = nt1;
+  }
+
+  auto geomSvc = service<IGeomSvc>("GeomSvc");
+  if(geomSvc){
+    const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
+    m_field = field.z()/dd4hep::tesla;
+    info() << "Magnetic field will obtain from GeomSvc = " << m_field << " tesla" << endmsg; 
+  }
+  else{
+    info() << "Failed to find GeomSvc ..." << endmsg;
+    info() << "Magnetic field will use what input through python option for this algorithm namse as Field, now " << m_field << " tesla" << endmsg;
+  }
+
+  _nEvt = 0;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpMCParticleAlg::execute(){
+  const edm4hep::MCParticleCollection* mcCols = nullptr;
+  try {
+    mcCols = _inMCColHdl.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "Collection " << _inMCColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+  }
+
+  if(mcCols){
+    m_nParticles = 0;
+    for(auto particle : *mcCols){
+      m_pdgID[m_nParticles] = particle.getPDG();
+      m_genStatus[m_nParticles] = particle.getGeneratorStatus();
+      m_simStatus[m_nParticles] = particle.getSimulatorStatus();
+      m_charge[m_nParticles] = particle.getCharge();
+      m_time[m_nParticles] = particle.getTime();
+      m_mass[m_nParticles] = particle.getMass();
+      const edm4hep::Vector3d& vertex = particle.getVertex();
+      m_vx[m_nParticles] = vertex.x;
+      m_vy[m_nParticles] = vertex.y;
+      m_vz[m_nParticles] = vertex.z;
+      const edm4hep::Vector3f& momentum = particle.getMomentum();
+      m_px[m_nParticles] = momentum.x;
+      m_py[m_nParticles] = momentum.y;
+      m_pz[m_nParticles] = momentum.z;
+
+      HelixClass helix;
+      float posV[3] = {vertex.x,vertex.y,vertex.z};
+      float momV[3] = {momentum.x,momentum.y,momentum.z};
+      helix.Initialize_VP(posV,momV,particle.getCharge(),m_field);
+      float phiMC = helix.getPhi0();
+      if(phiMC>CLHEP::pi) phiMC = phiMC - CLHEP::twopi;
+      m_phi0[m_nParticles] = phiMC;
+      m_d0[m_nParticles] = helix.getD0();
+      m_omega[m_nParticles] = helix.getOmega();
+      m_z0[m_nParticles] = helix.getZ0();
+      m_tanLambda[m_nParticles] = helix.getTanLambda();
+      m_nParticles++;
+    }
+    debug() << "MCParticle: " << m_nParticles <<endmsg;
+  }
+  m_tuple->write();
+  _nEvt++;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpMCParticleAlg::finalize(){
+  debug() << "Finalizing..." << endmsg;
+
+  return StatusCode::SUCCESS;
+}

Analysis/DumpEvent/src/DumpTrackAlg.cpp

+#include "DumpTrackAlg.h"
+
+#include "GaudiKernel/DataObject.h"
+#include "GaudiKernel/IHistogramSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/SmartDataPtr.h"
+
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/HelixClass.h"
+
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <math.h>
+
+DECLARE_COMPONENT( DumpTrackAlg )
+
+//------------------------------------------------------------------------------
+DumpTrackAlg::DumpTrackAlg( const std::string& name, ISvcLocator* pSvcLocator )
+    : Algorithm( name, pSvcLocator ) {
+  declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
+  declareProperty("TrackCollection", _inTrackColHdl, "Handle of the Input Track collection from CEPCSW");
+
+  m_thisName = name;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpTrackAlg::initialize(){
+  info() << "Booking Ntuple" << endmsg;
+
+  NTuplePtr nt1(ntupleSvc(), "MyTuples/Track"+m_thisName);
+  if ( !nt1 ) {
+    m_tuple = ntupleSvc()->book("MyTuples/Track"+m_thisName,CLID_ColumnWiseTuple,"Tracking result");
+    if ( 0 != m_tuple ) {
+      m_tuple->addItem ("ntrk",      m_nTracks, 0, 500 ).ignore();
+      m_tuple->addIndexedItem ("x",         m_nTracks, m_x ).ignore();
+      m_tuple->addIndexedItem ("y",         m_nTracks, m_y ).ignore();
+      m_tuple->addIndexedItem ("z",         m_nTracks, m_z ).ignore();
+      m_tuple->addIndexedItem ("px",        m_nTracks, m_px ).ignore();
+      m_tuple->addIndexedItem ("py",        m_nTracks, m_py ).ignore();
+      m_tuple->addIndexedItem ("pz",        m_nTracks, m_pz ).ignore();
+      m_tuple->addIndexedItem ("d0",        m_nTracks, m_d0 ).ignore();
+      m_tuple->addIndexedItem ("phi0",      m_nTracks, m_phi0 ).ignore();
+      m_tuple->addIndexedItem ("omega",     m_nTracks, m_omega ).ignore();
+      m_tuple->addIndexedItem ("z0",        m_nTracks, m_z0 ).ignore();
+      m_tuple->addIndexedItem ("tanLambda", m_nTracks, m_tanLambda ).ignore();
+      m_tuple->addIndexedItem ("sigma_d0",        m_nTracks, m_sigma_d0 ).ignore();
+      m_tuple->addIndexedItem ("sigma_phi0",      m_nTracks, m_sigma_phi0 ).ignore();
+      m_tuple->addIndexedItem ("sigma_omega",     m_nTracks, m_sigma_omega ).ignore();
+      m_tuple->addIndexedItem ("sigma_z0",        m_nTracks, m_sigma_z0 ).ignore();
+      m_tuple->addIndexedItem ("sigma_tanLambda", m_nTracks, m_sigma_tanLambda ).ignore();
+      m_tuple->addIndexedItem ("nvxd", m_nTracks, m_nHitsVXD ).ignore();
+      m_tuple->addIndexedItem ("nftd", m_nTracks, m_nHitsFTD ).ignore();
+      m_tuple->addIndexedItem ("nsit", m_nTracks, m_nHitsSIT ).ignore();
+      m_tuple->addIndexedItem ("ngas", m_nTracks, m_nHitsGAS ).ignore();
+      m_tuple->addIndexedItem ("nset", m_nTracks, m_nHitsSET ).ignore();
+    }
+    else { // did not manage to book the N tuple....
+      fatal() << "Cannot book MyTuples/Track"+m_thisName <<endmsg; 
+      return StatusCode::FAILURE;
+    }
+  }
+  else{
+    m_tuple = nt1;
+  }
+
+  auto geomSvc = service<IGeomSvc>("GeomSvc");
+  if(geomSvc){
+    const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
+    m_field = field.z()/dd4hep::tesla;
+    info() << "Magnetic field will obtain from GeomSvc = " << m_field << " tesla" << endmsg;
+  }
+  else{
+    info() << "Failed to find GeomSvc ..." << endmsg;
+    info() << "Magnetic field will use what input through python option for this algorithm namse as Field, now " << m_field << " tesla" << endmsg;
+  }
+
+  _nEvt = 0;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpTrackAlg::execute(){
+  const edm4hep::TrackCollection* trackCols = nullptr;
+  try {
+    trackCols = _inTrackColHdl.get();
+  }
+  catch ( GaudiException &e ) {
+    debug() << "Collection " << _inTrackColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+  }
+
+  if(trackCols){
+    m_nTracks = 0;
+    for(auto track : *trackCols){
+      // since possible more than one location=1 TrackState (not deleted in reconstruction), always use last one 
+      for(std::vector<edm4hep::TrackState>::const_iterator it=track.trackStates_end()-1; it!=track.trackStates_begin()-1; it--){
+	edm4hep::TrackState trackState = *it;
+	if(trackState.location!=1)continue;
+	m_d0[m_nTracks]              = trackState.D0;
+	m_phi0[m_nTracks]            = trackState.phi;
+	m_omega[m_nTracks]           = trackState.omega;
+	m_z0[m_nTracks]              = trackState.Z0;
+	m_tanLambda[m_nTracks]       = trackState.tanLambda;
+	m_sigma_d0[m_nTracks]        = std::sqrt(trackState.covMatrix[0]);
+	m_sigma_phi0[m_nTracks]      = std::sqrt(trackState.covMatrix[2]);
+	m_sigma_omega[m_nTracks]     = std::sqrt(trackState.covMatrix[5]);
+	m_sigma_z0[m_nTracks]        = std::sqrt(trackState.covMatrix[9]);
+        m_sigma_tanLambda[m_nTracks] = std::sqrt(trackState.covMatrix[14]);
+	HelixClass helix_fit;
+        helix_fit.Initialize_Canonical(trackState.phi,trackState.D0,trackState.Z0,trackState.omega,trackState.tanLambda,m_field);
+	m_x[m_nTracks] = helix_fit.getReferencePoint()[0];
+	m_y[m_nTracks] = helix_fit.getReferencePoint()[1];
+	m_z[m_nTracks] = helix_fit.getReferencePoint()[2];
+        m_px[m_nTracks] = helix_fit.getMomentum()[0];
+        m_py[m_nTracks] = helix_fit.getMomentum()[1];
+        m_pz[m_nTracks] = helix_fit.getMomentum()[2];
+	//info() << "size = " << track.subDetectorHitNumbers_size() << endmsg;
+	//for(int ii=0;ii<track.subDetectorHitNumbers_size();ii++){
+	//  std::cout << track.getSubDetectorHitNumbers(ii) << " ";
+	//}
+	//std::cout << std::endl;
+	if(track.subDetectorHitNumbers_size()>=5){
+	  m_nHitsVXD[m_nTracks] = track.getSubDetectorHitNumbers(0);
+	  m_nHitsFTD[m_nTracks] = track.getSubDetectorHitNumbers(1);
+	  m_nHitsSIT[m_nTracks] = track.getSubDetectorHitNumbers(2);
+	  m_nHitsGAS[m_nTracks] = track.getSubDetectorHitNumbers(3);
+	  m_nHitsSET[m_nTracks] = track.getSubDetectorHitNumbers(4);
+	}
+	else{
+	  m_nHitsVXD[m_nTracks] = 0;
+          m_nHitsSIT[m_nTracks] = 0;
+          m_nHitsSET[m_nTracks] = 0;
+          m_nHitsFTD[m_nTracks] = 0;
+          m_nHitsGAS[m_nTracks] = 0;
+	}
+	m_nTracks++;
+	break;
+      }
+    }
+  }
+  m_tuple->write();
+  _nEvt++;
+  return StatusCode::SUCCESS;
+}
+
+//------------------------------------------------------------------------------
+StatusCode DumpTrackAlg::finalize(){
+  debug() << "Finalizing..." << endmsg;
+
+  return StatusCode::SUCCESS;
+}