diff --git a/Analysis/DumpEvent/CMakeLists.txt b/Analysis/DumpEvent/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..3c5f697d2e4355b72a6ee88c16e271edba076353 --- /dev/null +++ b/Analysis/DumpEvent/CMakeLists.txt @@ -0,0 +1,19 @@ +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) diff --git a/Analysis/DumpEvent/src/DumpMCParticleAlg.cpp b/Analysis/DumpEvent/src/DumpMCParticleAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..56d68151489b3515c6089b9947b7a1ab7a78301a --- /dev/null +++ b/Analysis/DumpEvent/src/DumpMCParticleAlg.cpp @@ -0,0 +1,130 @@ +#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; +} diff --git a/Analysis/DumpEvent/src/DumpMCParticleAlg.h b/Analysis/DumpEvent/src/DumpMCParticleAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..bf9b16bce27ead67aa27b50db7fbb1a379dcfc73 --- /dev/null +++ b/Analysis/DumpEvent/src/DumpMCParticleAlg.h @@ -0,0 +1,51 @@ +#ifndef DumpMCParticleAlg_h +#define DumpMCParticleAlg_h 1 + +#include "k4FWCore/DataHandle.h" +#include "GaudiKernel/Algorithm.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/TrackCollection.h" + +#include "GaudiKernel/NTuple.h" + +class DumpMCParticleAlg : public Algorithm { + public: + // Constructor of this form must be provided + DumpMCParticleAlg( const std::string& name, ISvcLocator* pSvcLocator ); + + // Three mandatory member functions of any algorithm + StatusCode initialize() override; + StatusCode execute() override; + StatusCode finalize() override; + + private: + DataHandle<edm4hep::MCParticleCollection> _inMCColHdl{"MCParticle", Gaudi::DataHandle::Reader, this}; + + Gaudi::Property<double> m_field{this, "Field", 3.0}; + + NTuple::Tuple* m_tuple; + NTuple::Item<long> m_nParticles; + NTuple::Array<int> m_pdgID; + NTuple::Array<int> m_genStatus; + NTuple::Array<int> m_simStatus; + NTuple::Array<float> m_charge; + NTuple::Array<float> m_time; + NTuple::Array<double> m_mass; + NTuple::Array<double> m_vx; + NTuple::Array<double> m_vy; + NTuple::Array<double> m_vz; + NTuple::Array<float> m_px; + NTuple::Array<float> m_py; + NTuple::Array<float> m_pz; + NTuple::Array<float> m_d0; + NTuple::Array<float> m_phi0; + NTuple::Array<float> m_omega; + NTuple::Array<float> m_z0; + NTuple::Array<float> m_tanLambda; + + int _nEvt; + std::string m_thisName; +}; + +#endif diff --git a/Analysis/DumpEvent/src/DumpTrackAlg.cpp b/Analysis/DumpEvent/src/DumpTrackAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6a284e256481cc50a9c59a9060f6335695e86823 --- /dev/null +++ b/Analysis/DumpEvent/src/DumpTrackAlg.cpp @@ -0,0 +1,151 @@ +#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; +} diff --git a/Analysis/DumpEvent/src/DumpTrackAlg.h b/Analysis/DumpEvent/src/DumpTrackAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..91c526de7f1fe8d323a8dc6d4158050dbf53f348 --- /dev/null +++ b/Analysis/DumpEvent/src/DumpTrackAlg.h @@ -0,0 +1,56 @@ +#ifndef DumpTrackAlg_h +#define DumpTrackAlg_h 1 + +#include "k4FWCore/DataHandle.h" +#include "GaudiKernel/Algorithm.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/TrackCollection.h" + +#include "GaudiKernel/NTuple.h" + +class DumpTrackAlg : public Algorithm { + public: + // Constructor of this form must be provided + DumpTrackAlg( const std::string& name, ISvcLocator* pSvcLocator ); + + // Three mandatory member functions of any algorithm + StatusCode initialize() override; + StatusCode execute() override; + StatusCode finalize() override; + + private: + DataHandle<edm4hep::MCParticleCollection> _inMCColHdl{"MCParticle", Gaudi::DataHandle::Reader, this}; + DataHandle<edm4hep::TrackCollection> _inTrackColHdl{"SiTracks", Gaudi::DataHandle::Reader, this}; + + Gaudi::Property<double> m_field{this, "Field", 3.0}; + + NTuple::Tuple* m_tuple; + NTuple::Item<long> m_nTracks; + NTuple::Array<float> m_x; + NTuple::Array<float> m_y; + NTuple::Array<float> m_z; + NTuple::Array<float> m_px; + NTuple::Array<float> m_py; + NTuple::Array<float> m_pz; + NTuple::Array<float> m_d0; + NTuple::Array<float> m_phi0; + NTuple::Array<float> m_omega; + NTuple::Array<float> m_z0; + NTuple::Array<float> m_tanLambda; + NTuple::Array<float> m_sigma_d0; + NTuple::Array<float> m_sigma_phi0; + NTuple::Array<float> m_sigma_omega; + NTuple::Array<float> m_sigma_z0; + NTuple::Array<float> m_sigma_tanLambda; + NTuple::Array<int> m_nHitsVXD; + NTuple::Array<int> m_nHitsFTD; + NTuple::Array<int> m_nHitsSIT; + NTuple::Array<int> m_nHitsGAS; + NTuple::Array<int> m_nHitsSET; + + int _nEvt; + std::string m_thisName; +}; + +#endif