Geant4Output2EDM4hep.cpp

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : F.Gaede, DESY
//
//==========================================================================

#ifndef DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H
#define DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H

//  Framework include files
#include "DD4hep/Detector.h"
#include "DD4hep/VolumeManager.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4OutputAction.h"
#include "DDG4/Geant4SensDetAction.h"
#include "DDG4/EventParameters.h"

// Geant4 headers
#include "G4Threading.hh"
#include "G4AutoLock.hh"
#include <G4Version.hh>
#include <G4SystemOfUnits.hh>

// edm4hep include files
#include "edm4hep/EventHeaderCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
#include "edm4hep/CaloHitContributionCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "podio/EventStore.h"
#include "podio/ROOTWriter.h"

#include <typeinfo>
#include <iostream>
#include <ctime>


/// Namespace for the AIDA detector description toolkit
namespace dd4hep {

  class ComponentCast;

  /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
  namespace sim {

    template <class T=podio::EventStore> void EventParameters::extractParameters(T& event){
      auto& lcparameters = event.getEventMetaData();

      for(auto const& ival: this->intParameters()) {
        lcparameters.setValues(ival.first, ival.second);
      }
      for(auto const& ival: this->fltParameters()) {
        lcparameters.setValues(ival.first, ival.second);
      }
      for(auto const& ival: this->strParameters()) {
        lcparameters.setValues(ival.first, ival.second);
      }
    }

    class  Geant4ParticleMap;

    /// Base class to output Geant4 event data to EDM4hep
    /**
     *  \author  F.Gaede
     *  \version 1.0
     *  \ingroup DD4HEP_SIMULATION
     */
    class Geant4Output2EDM4hep : public Geant4OutputAction  {
    protected:
      podio::EventStore*  m_store;
      podio::ROOTWriter*  m_file;
      int              m_runNo;
      int              m_runNumberOffset;
      int              m_eventNumberOffset;
      std::map< std::string, std::string > m_runHeader;
      std::map< std::string, std::string > m_eventParametersInt;
      std::map< std::string, std::string > m_eventParametersFloat;
      std::map< std::string, std::string > m_eventParametersString;
      bool m_FirstEvent =  true  ;

      /// create the podio collections for the particles and hits
      void createCollections(OutputContext<G4Event>& ctxt) ;
      /// Data conversion interface for MC particles to EDM4hep format
      void saveParticles(Geant4ParticleMap* particles);
    public:
      /// Standard constructor
      Geant4Output2EDM4hep(Geant4Context* ctxt, const std::string& nam);
      /// Default destructor
      virtual ~Geant4Output2EDM4hep();
      /// Callback to store the Geant4 run information
      virtual void beginRun(const G4Run* run);
      /// Callback to store the Geant4 run information
      virtual void endRun(const G4Run* run);

      /// Callback to store the Geant4 run information
      virtual void saveRun(const G4Run* run);
      /// Callback to store the Geant4 event
      virtual void saveEvent( OutputContext<G4Event>& ctxt);
      /// Callback to store each Geant4 hit collection
      virtual void saveCollection( OutputContext<G4Event>& ctxt, G4VHitsCollection* collection);
      /// Commit data at end of filling procedure
      virtual void commit( OutputContext<G4Event>& ctxt);

      /// begin-of-event callback - creates EDM4hep event and adds it to the event context
      virtual void begin(const G4Event* event);
    protected:
      /// Fill event parameters in EDM4hep event
      template <typename T>
      void saveEventParameters(const std::map<std::string, std::string >& parameters);
    };
    
    /// Fill event parameters in EDM4hep event
    template <typename T>
    inline void Geant4Output2EDM4hep::saveEventParameters(const std::map<std::string, std::string >& parameters)  {
      for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter)  {
        T parameter;
        std::istringstream iss(iter->second);
        if ( (iss >> parameter).fail() )  {
          printout(FATAL,"saveEventParameters","+++ Event parameter %s: FAILED to convert to type :%s",iter->first.c_str(),typeid(T).name());
          continue;
        }
	auto& evtMD = m_store->getEventMetaData();
	evtMD.setValue(iter->first,parameter);
      }
    }

    /// Fill event parameters in EDM4hep event - std::string specialization
    template <>
    inline void Geant4Output2EDM4hep::saveEventParameters<std::string>(const std::map<std::string, std::string >& parameters)  {
      for(std::map<std::string, std::string >::const_iterator iter = parameters.begin(), endIter = parameters.end() ; iter != endIter ; ++iter)  {
	auto& evtMD = m_store->getEventMetaData();
	evtMD.setValue(iter->first,iter->second);
      }
    }

  }    // End namespace sim
}      // End namespace dd4hep
#endif // DD4HEP_DDG4_GEANT4OUTPUT2EDM4hep_H

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : F.Gaede, DESY
//
//==========================================================================

// Framework include files
#include "DD4hep/InstanceCount.h"
#include "DD4hep/Detector.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4DataConversion.h"
#include "DDG4/Geant4Context.h"
#include "DDG4/Geant4Particle.h"
#include "DDG4/Geant4Data.h"
#include "DDG4/Geant4Action.h"

//#include "DDG4/Geant4Output2EDM4hep.h"
#include "G4ParticleDefinition.hh"
#include "G4VProcess.hh"
#include "G4Event.hh"
#include "G4Run.hh"


using namespace dd4hep::sim;
using namespace dd4hep;
using namespace std;
namespace {
  G4Mutex action_mutex=G4MUTEX_INITIALIZER;
}

#include "DDG4/Factories.h"
DECLARE_GEANT4ACTION(Geant4Output2EDM4hep)

/// Standard constructor
Geant4Output2EDM4hep::Geant4Output2EDM4hep(Geant4Context* ctxt, const string& nam)
: Geant4OutputAction(ctxt,nam), m_store(0), m_file(0), m_runNo(0), m_runNumberOffset(0), m_eventNumberOffset(0)
{
  declareProperty("RunHeader", m_runHeader);
  declareProperty("EventParametersInt",    m_eventParametersInt);
  declareProperty("EventParametersFloat",  m_eventParametersFloat);
  declareProperty("EventParametersString", m_eventParametersString);
  declareProperty("RunNumberOffset", m_runNumberOffset);
  declareProperty("EventNumberOffset", m_eventNumberOffset);
  printout( INFO, "Geant4Output2EDM4hep" ," instantiated ..." ) ;
  InstanceCount::increment(this);
}

/// Default destructor
Geant4Output2EDM4hep::~Geant4Output2EDM4hep()  {
  G4AutoLock protection_lock(&action_mutex);
  if ( m_file )  {
    m_file->finish();
    detail::deletePtr(m_file);
  }
  if (nullptr != m_store) {
    delete m_store;
  }
  InstanceCount::decrement(this);
}

// Callback to store the Geant4 run information
void Geant4Output2EDM4hep::beginRun(const G4Run* run)  {
  G4AutoLock protection_lock(&action_mutex);
  if ( 0 == m_file && !m_output.empty() )   {
    m_store = new podio::EventStore ;
    m_file = new podio::ROOTWriter(m_output, m_store);

    printout( INFO, "Geant4Output2EDM4hep" ," opened %s for output", m_output.c_str() ) ;
  }
  
  saveRun(run);
}

/// Callback to store the Geant4 run information
void Geant4Output2EDM4hep::endRun(const G4Run* /*run*/)  {
  // saveRun(run);
}

/// Commit data at end of filling procedure
void Geant4Output2EDM4hep::commit( OutputContext<G4Event>& /* ctxt */)   {
  if ( m_file )   {
    G4AutoLock protection_lock(&action_mutex);
    m_file->writeEvent();
    m_store->clearCollections();
    return;
  }
  except("+++ Failed to write output file. [Stream is not open]");
}

/// Callback to store the Geant4 run information
void Geant4Output2EDM4hep::saveRun(const G4Run* /*run*/)  {
  //G4AutoLock protection_lock(&action_mutex);

  printout( WARNING, "Geant4Output2EDM4hep" ,"saveRun(): RunHeader not implemented in EDM4hep, nothing written ..." ) ;

  // --- write an edm4hep::RunHeader ---------
  //FIXME: need a suitable Runheader object in EDM4hep

//  edm4hep::LCRunHeaderImpl* rh =  new edm4hep::LCRunHeaderImpl;
//  for (std::map< std::string, std::string >::iterator it = m_runHeader.begin(); it != m_runHeader.end(); ++it) {
//    rh->parameters().setValue( it->first, it->second );
//  }
//  m_runNo = m_runNumberOffset > 0 ? m_runNumberOffset + run->GetRunID() : run->GetRunID();
//  rh->parameters().setValue("GEANT4Version", G4Version);
//  rh->parameters().setValue("DD4HEPVersion", versionString());
//  rh->setRunNumber(m_runNo);
//  rh->setDetectorName(context()->detectorDescription().header().name());
//  m_file->writeRunHeader(rh);
}

void Geant4Output2EDM4hep::begin(const G4Event* /* event */)  {
}

/// Data conversion interface for MC particles to EDM4hep format
void Geant4Output2EDM4hep::saveParticles(Geant4ParticleMap* particles)    {
  typedef detail::ReferenceBitMask<const int> PropertyMask;
  typedef Geant4ParticleMap::ParticleMap ParticleMap;
  const ParticleMap& pm = particles->particleMap;
  size_t nparts = pm.size();
  edm4hep::MCParticleCollection* mcpc =
    const_cast<edm4hep::MCParticleCollection*>(
      &m_store->get<edm4hep::MCParticleCollection>("MCParticles"));

  if ( nparts > 0 )  {
    size_t cnt = 0;
    map<int,int> p_ids;
    vector<const Geant4Particle*> p_part(pm.size(),0);
    vector<edm4hep::MCParticle> p_edm4hep(pm.size());
    // First create the particles
    for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end();++i, ++cnt)   {
      int id = (*i).first;
      const Geant4ParticleHandle p = (*i).second;
      PropertyMask mask(p->status);
      //      std::cout << " ********** mcp status : 0x" << std::hex << p->status << ", mask.isSet(G4PARTICLE_GEN_STABLE) x" << std::dec << mask.isSet(G4PARTICLE_GEN_STABLE)  <<std::endl ;
      const G4ParticleDefinition* def = p.definition();
      edm4hep::MCParticle mcp = edm4hep::MCParticle();
      mcp.setPDG(p->pdgID);

      float ps_fa[3] = {float(p->psx/CLHEP::GeV),float(p->psy/CLHEP::GeV),float(p->psz/CLHEP::GeV)};
      mcp.setMomentum( ps_fa );

      float pe_fa[3] = {float(p->pex/CLHEP::GeV),float(p->pey/CLHEP::GeV),float(p->pez/CLHEP::GeV)};
      mcp.setMomentumAtEndpoint( pe_fa );

      double vs_fa[3] = { p->vsx/CLHEP::mm, p->vsy/CLHEP::mm, p->vsz/CLHEP::mm } ;
      mcp.setVertex( vs_fa );

      double ve_fa[3] = { p->vex/CLHEP::mm, p->vey/CLHEP::mm, p->vez/CLHEP::mm } ;
      mcp.setEndpoint( ve_fa );

      mcp.setTime(p->time/CLHEP::ns);
      mcp.setMass(p->mass/CLHEP::GeV);
      mcp.setCharge(def ? def->GetPDGCharge() : 0); // Charge(e+) = 1 !

      // Set generator status
      mcp.setGeneratorStatus(0);
      if( p->genStatus ) {
        mcp.setGeneratorStatus( p->genStatus ) ;
      } else {
        if ( mask.isSet(G4PARTICLE_GEN_STABLE) )             mcp.setGeneratorStatus(1);
        else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) )       mcp.setGeneratorStatus(2);
        else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) mcp.setGeneratorStatus(3);
        else if ( mask.isSet(G4PARTICLE_GEN_BEAM) )          mcp.setGeneratorStatus(4);
        else if ( mask.isSet(G4PARTICLE_GEN_OTHER) )         mcp.setGeneratorStatus(9);
      }

      // Set simulation status
      mcp.setCreatedInSimulation(         mask.isSet(G4PARTICLE_SIM_CREATED) );
      mcp.setBackscatter(                 mask.isSet(G4PARTICLE_SIM_BACKSCATTER) );
      mcp.setVertexIsNotEndpointOfParent( mask.isSet(G4PARTICLE_SIM_PARENT_RADIATED) );
      mcp.setDecayedInTracker(            mask.isSet(G4PARTICLE_SIM_DECAY_TRACKER) );
      mcp.setDecayedInCalorimeter(        mask.isSet(G4PARTICLE_SIM_DECAY_CALO) );
      mcp.setHasLeftDetector(             mask.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) );
      mcp.setStopped(                     mask.isSet(G4PARTICLE_SIM_STOPPED) );
      mcp.setOverlay(                     false );

      //fg: if simstatus !=0 we have to set the generator status to 0:
      if( mcp.isCreatedInSimulation() )
        mcp.setGeneratorStatus( 0 )  ;

      mcp.setSpin(p->spin);
      mcp.setColorFlow(p->colorFlow);

      mcpc->push_back(mcp);
      p_ids[id] = cnt;
      p_part[cnt] = p;
      p_edm4hep[cnt] = mcp;
    }

    // Now establish parent-daughter relationships
    for(size_t i=0, n=p_ids.size(); i<n; ++i)   {
      map<int,int>::iterator k;
      const Geant4Particle* p = p_part[i];
      edm4hep::MCParticle q = p_edm4hep[i];
      const Geant4Particle::Particles& dau = p->daughters;
      for(Geant4Particle::Particles::const_iterator j=dau.begin(); j!=dau.end(); ++j)  {
        int idau = *j;
        if ( (k=p_ids.find(idau)) == p_ids.end() )  {  // Error!!!
          printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
          continue;
        }
        int iqdau = (*k).second;
        edm4hep::MCParticle qdau = p_edm4hep[iqdau];
        qdau.addToParents(q);
      }
      const Geant4Particle::Particles& par = p->parents;
      for(Geant4Particle::Particles::const_iterator j=par.begin(); j!=par.end(); ++j)  {
        int ipar = *j; // A parent ID iof -1 means NO parent, because a base of 0 is perfectly leagal!
        if ( ipar>=0 && (k=p_ids.find(ipar)) == p_ids.end() )  {  // Error!!!
          printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
          continue;
        }
        int iqpar = (*k).second;
        edm4hep::MCParticle qpar = p_edm4hep[iqpar];
        q.addToParents(qpar);
      }
    }
  }
}

/// Callback to store the Geant4 event
void Geant4Output2EDM4hep::saveEvent(OutputContext<G4Event>& ctxt)  {

  if( m_FirstEvent ){
    createCollections( ctxt ) ;
    m_FirstEvent = false ;
  }

  EventParameters* parameters = context()->event().extension<EventParameters>(false);
  int runNumber(0), eventNumber(0);
  const int eventNumberOffset(m_eventNumberOffset > 0 ? m_eventNumberOffset : 0);
  const int runNumberOffset(m_runNumberOffset > 0 ? m_runNumberOffset : 0);
  // Get event number, run number and parameters from extension ...
  if ( parameters ) {
    runNumber = parameters->runNumber() + runNumberOffset;
    eventNumber = parameters->eventNumber() + eventNumberOffset;
    parameters->extractParameters(*m_store);
  } else { // ... or from DD4hep framework
    runNumber = m_runNo + runNumberOffset;
    eventNumber = ctxt.context->GetEventID() + eventNumberOffset;
  }
  printout(INFO,"Geant4Output2EDM4hep","+++ Saving EDM4hep event %d run %d.", eventNumber, runNumber);

  // this does not compile as create() is we only get a const ref - need to review PODIO EventStore API
  // auto& evtHCol = m_store->get<edm4hep::EventHeaderCollection>("EventHeader") ;
  // auto evtHdr = evtHCol.create() ;
  auto* evtHCol = const_cast<edm4hep::EventHeaderCollection*>(&m_store->get<edm4hep::EventHeaderCollection>("EventHeader") );
  auto evtHdr = evtHCol->create() ;

  evtHdr.setRunNumber(runNumber);
  evtHdr.setEventNumber(eventNumber);
//not implemented in EDM4hep ?  evtHdr.setDetectorName(context()->detectorDescription().header().name());
  evtHdr.setTimeStamp( std::time(nullptr) ) ;

  saveEventParameters<int>(m_eventParametersInt);
  saveEventParameters<float>(m_eventParametersFloat);
  saveEventParameters<std::string>(m_eventParametersString);

  Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>(false);
  if ( part_map )   {
    print("+++ Saving %d EDM4hep particles....",int(part_map->particleMap.size()));
    if ( part_map->particleMap.size() > 0 )  {
      saveParticles(part_map);
    }
  }
}

/// Callback to store each Geant4 hit collection
void Geant4Output2EDM4hep::saveCollection(OutputContext<G4Event>& /*ctxt*/, G4VHitsCollection* collection)  {

  size_t nhits = collection->GetSize();
  std::string colName = collection->GetName();

  printout(DEBUG,"Geant4Output2EDM4hep","+++ Saving EDM4hep collection %s with %d entries.",
     colName.c_str(),int(nhits));

  Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(collection);
  if( coll == nullptr ){
    printout(ERROR, "Geant4Output2EDM4hep" , " no Geant4HitCollection:  %s ", colName.c_str() );
    return ;
  }

  Geant4ParticleMap* pm = context()->event().extension<Geant4ParticleMap>(false);

  edm4hep::MCParticleCollection* mcpc =
    const_cast<edm4hep::MCParticleCollection*>(
      &m_store->get<edm4hep::MCParticleCollection>("MCParticles"));

  //-------------------------------------------------------------------
  if( typeid( Geant4Tracker::Hit ) == coll->type().type  ){

    edm4hep::SimTrackerHitCollection* sthc =
      const_cast<edm4hep::SimTrackerHitCollection*>(&m_store->get<edm4hep::SimTrackerHitCollection>(colName));

    for(unsigned i=0 ; i < nhits ; ++i){
      auto sth = sthc->create() ;

      const Geant4Tracker::Hit* hit = coll->hit(i);
      const Geant4Tracker::Hit::Contribution& t = hit->truth;
      int trackID = pm->particleID(t.trackID);

      auto mcp = mcpc->at(trackID);

      sth.setCellID( hit->cellID ) ;
      sth.setEDep(hit->energyDeposit/CLHEP::GeV);
      sth.setPathLength(hit->length/CLHEP::mm);
      sth.setTime(hit->truth.time/CLHEP::ns);
      sth.setMCParticle(mcp);
      sth.setPosition({hit->position.x()/CLHEP::mm,
           hit->position.y()/CLHEP::mm,
           hit->position.z()/CLHEP::mm});
      sth.setMomentum(edm4hep::Vector3f(hit->momentum.x()/CLHEP::GeV,
          hit->momentum.y()/CLHEP::GeV,
          hit->momentum.z()/CLHEP::GeV ));

      auto particleIt = pm->particles().find(trackID);
      if( ( particleIt != pm->particles().end()) ){
        // if the original track ID of the particle is not the same as the
        // original track ID of the hit it was produced by an MCParticle that
        // is no longer stored
        sth.setProducedBySecondary( (particleIt->second->originalG4ID != t.trackID) );
      }
    }
  //-------------------------------------------------------------------
  }
  else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type ){

    Geant4Sensitive*       sd      = coll->sensitive();
    int hit_creation_mode = sd->hitCreationMode();

    edm4hep::SimCalorimeterHitCollection* sCaloHitColl =
      const_cast<edm4hep::SimCalorimeterHitCollection*>(
  &m_store->get<edm4hep::SimCalorimeterHitCollection>(colName));

    colName += "Contributions"  ;

    edm4hep::CaloHitContributionCollection* sCaloHitContColl =
      const_cast<edm4hep::CaloHitContributionCollection*>(
  &m_store->get<edm4hep::CaloHitContributionCollection>(colName));


    for(unsigned i=0 ; i < nhits ; ++i){
      auto sch = sCaloHitColl->create() ;

      const Geant4Calorimeter::Hit* hit = coll->hit(i);

      sch.setCellID( hit->cellID );
      sch.setPosition({
           float(hit->position.x()/CLHEP::mm),
           float(hit->position.y()/CLHEP::mm),
           float(hit->position.z()/CLHEP::mm)});
      sch.setEnergy( hit->energyDeposit );

      // now add the individual step contributions
      for(Geant4HitData::Contributions::const_iterator ci=hit->truth.begin();
        ci!=hit->truth.end(); ++ci){

        auto sCaloHitCont = sCaloHitContColl->create();
        sch.addToContributions( sCaloHitCont );

        const Geant4HitData::Contribution& c = *ci;
        int trackID = pm->particleID(c.trackID);
        auto mcp = mcpc->at(trackID);

        sCaloHitCont.setEnergy( c.deposit/CLHEP::GeV );
        sCaloHitCont.setTime( c.time/CLHEP::ns );
        sCaloHitCont.setParticle( mcp );

        if ( hit_creation_mode == Geant4Sensitive::DETAILED_MODE )     {
          sCaloHitCont.setPDG( c.pdgID );
          sCaloHitCont.setStepPosition( edm4hep::Vector3f(
            c.x/CLHEP::mm,
            c.y/CLHEP::mm, 
            c.z/CLHEP::mm) );
        }
      }
    }
  //-------------------------------------------------------------------
  } else {

    printout(ERROR, "Geant4Output2EDM4hep" , " unknown type in Geant4HitCollection  %s ",
       coll->type().type.name() );
  }
}


void Geant4Output2EDM4hep::createCollections(OutputContext<G4Event>& ctxt){

  m_store->create<edm4hep::MCParticleCollection>("MCParticles");
  m_file->registerForWrite("MCParticles");
  printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection MCParticles" );

  m_store->create<edm4hep::EventHeaderCollection>("EventHeader") ;
  m_file->registerForWrite("EventHeader");
  printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection EventHeader" );


  const G4Event* evt = ctxt.context ;
  G4HCofThisEvent* hce = evt->GetHCofThisEvent();
  int nCol = hce->GetNumberOfCollections();

  for (int i = 0; i < nCol; ++i) {
    G4VHitsCollection* hc = hce->GetHC(i);
    std::string colName =  hc->GetName() ;
    Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(hc);
    if( coll == nullptr ){
      printout(WARNING, "Geant4Output2EDM4hep" , " no Geant4HitCollection:  %s ", colName.c_str() );
      continue ;
    }

    if( typeid( Geant4Tracker::Hit ) == coll->type().type  ){

      m_store->create<edm4hep::SimTrackerHitCollection>(colName);
      m_file->registerForWrite(colName);
      printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );
    }
    else if( typeid( Geant4Calorimeter::Hit ) == coll->type().type ){

      m_store->create<edm4hep::SimCalorimeterHitCollection>(colName);
      m_file->registerForWrite(colName);
      printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );

      colName += "Contributions"  ;
      m_store->create<edm4hep::CaloHitContributionCollection>(colName);
      m_file->registerForWrite(colName);
      printout(DEBUG,"Geant4Output2EDM4hep","+++ created collection %s",colName.c_str() );

    } else {

      printout(WARNING, "Geant4Output2EDM4hep" ,
         " unknown type in Geant4HitCollection  %s ", coll->type().type.name() );
    }
  }


}