LCDG4LcioHelper.cc

Go to the documentation of this file.

//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
//
// File: LCDG4LcioHelper.cc
// Module: lcdg4
//
// Purpose:  Prepare LCIO data objects for persistency, and
//      save them into an LCIO output file
//
// 20040805 - G.Lima - Created
//
//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/
#include "LCDG4LcioHelper.hh"
#include "LCDG4McPartManager.hh"
#include "G4ThreeVector.hh"
#include "G4Event.hh"
#include "G4Run.hh"

#include "IO/LCWriter.h"
#include "IMPL/LCRunHeaderImpl.h" 
#include "IMPL/LCEventImpl.h"
#include "IMPL/LCCollectionVec.h"
#include "UTIL/LCTOOLS.h"
using IMPL::LCEventImpl;
using IMPL::LCCollectionVec;
using UTIL::LCTOOLS;

#include "EVENT/LCIO.h"
#include "IMPL/SimCalorimeterHitImpl.h"
#include "IMPL/SimTrackerHitImpl.h"
using IMPL::SimCalorimeterHitImpl;
using IMPL::SimTrackerHitImpl;

#include "IOIMPL/LCFactory.h"
using IMPL::SimCalorimeterHitImpl;
using IMPL::SimTrackerHitImpl;
using IMPL::MCParticleImpl;
using EVENT::LCIO;
// using IO::LCWriter;
using IMPL::LCRunHeaderImpl;

#include <sstream>
#include <string>
#include <map>
#include <vector>
using std::stringstream;
using std::string;
using std::cout;
using std::endl;
using std::map;
using std::vector;

// initialize static pointer to first instance
// (only one instance should be created, but this is not inforced)
LCDG4LcioHelper* LCDG4LcioHelper::_me = new LCDG4LcioHelper();

// singleton access point
LCDG4LcioHelper* LCDG4LcioHelper::getInstance() {
  return _me;
}

// constructor
LCDG4LcioHelper::LCDG4LcioHelper()
  : _writer(0), _mcPartColl(0)
{
  if( _me==0 ) {
    cout << "*** LcioHelper constructor: _me="<< _me << endl;
    _writer = IOIMPL::LCFactory::getInstance()->createLCWriter();
    _mcPartColl = new map<int, MCParticleImpl*>();
  }
  else {
    // Inforce single instance
    this->~LCDG4LcioHelper();
  }
}

// destructor
LCDG4LcioHelper::~LCDG4LcioHelper() {
  if(_writer) {
    _writer->close();
    delete _writer;
    _writer = NULL;
  }
  if(_mcPartColl) {
    delete _mcPartColl;
    _mcPartColl = NULL;
  }
}

// File open
void LCDG4LcioHelper::openLcioFile(const string& filename) {
  G4cout << "attempting to open <" << filename << ">" << G4endl;
  if(_writer) _writer->open( filename, LCIO::WRITE_NEW );
}

//.. Create calorimeter-like hits for lcio output
SimCalorimeterHitImpl*
LCDG4LcioHelper::getCalHitLcio(const LCDG4CalHit& inhit) {
  SimCalorimeterHitImpl* outhit = new SimCalorimeterHitImpl();
  outhit->setCellID0( inhit.GetCellID() );
  outhit->setCellID1( 0 );

  const Hep3Vector inpos = inhit.GetPos()/mm;  
  float pos[3] = { inpos.x(), inpos.y(), inpos.z() };
  outhit->setPosition( pos );

  // loop over energy contributions
  LCDG4McPartManager* mcmgr = LCDG4McPartManager::GetPointer();
  unsigned int nmc = inhit.GetNMC();
  for(unsigned int imc=0; imc<nmc; ++imc) {
    float energy = inhit.GetEMC(imc)/GeV;
    float time = inhit.GetTimeMC(imc)/ns;
    // link to contributing MCParticle
    int trkid = inhit.GetTrkID(imc);
    G4int mctag = mcmgr->trkidToMctag(trkid);
    if(mctag != -1) {
      MCParticleImpl* ppart = _mcPartColl->find(mctag)->second;
      outhit->addMCParticleContribution( ppart, energy, time,
                                         ppart->getPDG() );
    }
    else {
      cout << "LcioHelper: mctag==-1?" << endl;
    }
  }
  return outhit;
}

// Create tracker-like hits
SimTrackerHitImpl*
LCDG4LcioHelper::getTrkHitLcio(const LCDG4TrackerHit& inhit) {

  // pack hit info into tag
  unsigned int tag = inhit.GetLayer();
  tag |= (inhit.GetSystem() << 28);
  tag |= (inhit.GetNorthSouth() << 30);
  tag |= (inhit.GetBarEnd() << 31);

  SimTrackerHitImpl* outhit = new SimTrackerHitImpl();
  outhit->setCellID( tag );
  outhit->setdEdx( inhit.GetEdep()/GeV );
  outhit->setTime( inhit.GetTdep()/ns );

  const Hep3Vector inpos = inhit.GetPos()/mm;
  double pos[3] = { inpos.x(), inpos.y(), inpos.z() };
//   pos[0] = inhit.GetPos().x()/mm;
//   pos[1] = inhit.GetPos().y()/mm;
//   pos[2] = inhit.GetPos().z()/mm;
  outhit->setPosition( pos );

  // link to particle which produced hit
  int trkid = inhit.GetTrackID();

  G4int mctag = LCDG4McPartManager::GetPointer()->trkidToMctag(trkid);
  assert(mctag>=0);
//   if(mctag != -1) {
    MCParticleImpl* ppart = _mcPartColl->find(mctag)->second;
    outhit->setMCParticle( ppart );
//   }
  return outhit;
}

// Create MCParticleImpl objects for LCIO output
MCParticleImpl* LCDG4LcioHelper::getMCParticleLcio(LCDMcPart& inpart) {
  // create an MCParticleImpl
  MCParticleImpl* outpart = new MCParticleImpl();

  // set PDGID
  outpart->setPDG( inpart.GetParticleID() );

  // set simulation status
  int simstat = inpart.GetG4Status();
  switch(simstat) {
  case 0:
    cout<<"inpart: ALIVE, tag="<< inpart.GetParticleTag() << endl;
//     assert(false);
    break; // alive
  case 1: outpart->setDecayedInTracker(true); break; // decayed
  case 2: break;  // interacted
  case 3: outpart->setHasLeftDetector(true); break; // left
  case 4: outpart->setStopped(true); break; // stopped
  case 5: assert(false); break; // looped
  case 6: assert(false); break; // lost
  case 7: assert(false); break; // stuck
  case 8: break;  // documentation
  case 9: outpart->isDecayedInCalorimeter(); break;  // showered
  case 10: assert(false); break; // maxsteps
  default: cout<< "Problem: simstat="<< simstat << endl;
  }

  // set generator status
  outpart->setGeneratorStatus( inpart.GetGenStatus() );
  if( outpart->getGeneratorStatus()==4 ) {  // generated in simulation
    outpart->setGeneratorStatus( 0 );
    outpart->setCreatedInSimulation(true);
  }

  // set production vertex
  Hep3Vector pos = inpart.GetPosition()/mm;
  double daux[3] = { pos.x(), pos.y(), pos.z() };
  outpart->setVertex( daux );

  // set end point
  pos = inpart.GetTermPosition()/mm;
  daux[0] = pos.x();  daux[1] = pos.y();  daux[2] = pos.z();
  outpart->setEndpoint( daux );

  // set mass and momentum
  HepLorentzVector pmom = inpart.Get4Momentum()/GeV;
  float faux[3] = { pmom.px(), pmom.py(), pmom.pz() };
  outpart->setMomentum( faux );
  outpart->setMass( pmom.mag() );
  outpart->setCharge( inpart.GetCharge() );

  // Note: parentage info still not setup at this point
  return outpart;
}

// Build a MCParticle collection
void LCDG4LcioHelper::buildMCParticleColl(const vector<LCDMcPart*>& lcdMcpColl) {
  // make sure collection is cleared
  assert( _mcPartColl->size()==0 );

  // MCParticle* to tag is a handy map
  typedef std::map< MCParticleImpl*, int > PointerToIndexMap;

  vector<LCDMcPart*>::const_iterator mciter = lcdMcpColl.begin();
//   int mapindex = 0;  // tags start at one, no zero

//   // sio format requires that first element is parent of all other elements.
//   // Its ParticleID is set to 99999999 (HePEVT id).
//   // Sometimes the input stdhep file already has this HepEVT element, then
//   // keep it.  Otherwise, create a new one for the .sio output
//   //  if( (*mciter)->GetParticleID()!=99999999 ) { // no, it did not work...
//   MCParticleSIO* mcparthead = new MCParticleSIO(99999999,0,0);
//   mcparthead->setfsflag(300);
//   mcpartmap[mapindex] = mcparthead;
//   mapindex++;
//   //  }

  int mctag = 1;  // tags start at one, not zero
  for(mciter = lcdMcpColl.begin();
      mciter != lcdMcpColl.end();
      mciter++, mctag++) {
    // make sure LCDMcPart is valid
    LCDMcPart* pMcPart = *mciter;
    if(!pMcPart) continue;  // some McParts have been erased!

    // Create a new MCParticleImpl..
    MCParticleImpl* mcplcio = getMCParticleLcio(*pMcPart);
    //.. and save it into a map keyed by tag
    _mcPartColl->insert( std::pair<int,MCParticleImpl*>(mctag,mcplcio) );
    _mcp2tag[ mcplcio ] = mctag;
  }

  // Loop to set up parentage links
  buildParentageLinks(lcdMcpColl);
}

// Set parentage relationships.  Assumes MCParticle collection is complete
void LCDG4LcioHelper::buildParentageLinks(const vector<LCDMcPart*>& lcdMcpColl) {
  // Loop over LCDMcParts and tags
  vector<LCDMcPart*>::const_iterator mciter = lcdMcpColl.begin();
  int mctag=1; // tags start at 1
  for(mciter = lcdMcpColl.begin();
      mciter != lcdMcpColl.end();
      mciter++, mctag++) {
    // make sure LCDMcPart is valid, as some of them have been erased
    LCDMcPart* plcdPart = *mciter;
    if(!plcdPart) continue;

    // Make sure there is a corresponding MCParticleImpl
    MCParticleImpl* pLcioPart = NULL;
    TagToMCParticleMap::iterator ipair = _mcPartColl->find(mctag);
    if( ipair!=_mcPartColl->end() ) pLcioPart = ipair->second;
    assert( pLcioPart );

    // Fetch parent info from LCDMcPart
    int par1,par2;
    plcdPart->GetMotherTags(par1,par2);

    // loop over parents
    int nparents = 0;
    int parmax=par1;
    if(par2>par1) parmax=par2;
    for(int ii=par1; ii<=parmax; ++ii) {
      MCParticleImpl* pLcioParent = NULL;
      TagToMCParticleMap::iterator ipair = _mcPartColl->find(ii);
      if( ipair!=_mcPartColl->end() ) pLcioParent = ipair->second;
      if(pLcioParent) {
        // associate particles
        pLcioPart->addParent( pLcioParent );
        ++nparents;
      }
//       if(ii==parmax && nparents==0) {
//      cout<< "***** No parent to be added, iParTag="<< ii << endl;
//       }
    }
  }
}

// Make sure that no daughter can come before its parent
// in the output LCIO particle table
void LCDG4LcioHelper::prepareMCParticleCollection(LCCollectionVec* newColl) {
  TagToMCParticleMap::iterator head = NULL;
  while( _mcPartColl->size()>0 ) {
    head = _mcPartColl->begin();
    G4cout<< "*** prepare: size="<< _mcPartColl->size()
          <<", head tag="<< _mcp2tag[head->second] << G4endl;
    insertAllFamily(newColl, head);
  }
}

// Recursive method used for reordering a la SIO
// It inserts the head into the new collection, calls itself for each
// daughter and deletes the head from the original map.
// It needs the head iterator in order to erase it from the original
// collection, otherwise just MCParticleImpl* would be necessary.
void LCDG4LcioHelper::insertAllFamily(LCCollectionVec* newColl,
                                      TagToMCParticleMap::iterator head) {
  // First of all, insert head
  G4cout<< " pushing tag="<< head->second << G4endl;
  G4cout<< " tag="<< _mcp2tag[ head->second ] << G4endl;
  newColl->push_back( head->second );
  // then loop over all daughters, inserting them recursively
  int ndau = head->second->getNumberOfDaughters();
  for(int i=0; i<ndau; ++i) {
    G4cout<< " ndau="<< ndau <<", i="<< i << G4endl;
    MCParticleImpl* dau
      = dynamic_cast<MCParticleImpl*>(head->second->getDaughter(i));
    int key = _mcp2tag[ dau ];
    insertAllFamily( newColl, _mcPartColl->find(key) );
  }
  // remove head from the original collection
  G4cout<< " erasing tag="<< _mcp2tag[head->second] << G4endl;
  _mcPartColl->erase(head);
}

void LCDG4LcioHelper::saveCollection(const string& collName,
                                     LCCollectionVec* coll) {
  std::pair<string,LCCollectionVec*> tmpPair(collName,coll);
  _savedColls.insert( tmpPair );
}

void LCDG4LcioHelper::writeEvent(const G4Event* g4evt) {

  // we need to use the implementation classes here 
  LCEventImpl*  lcevt = new LCEventImpl();
  lcevt->setEventNumber( g4evt->GetEventID() );
  lcevt->setRunNumber( _runHdr->getRunNumber() );
  lcevt->setDetectorName( _runHdr->getDetectorName() );

  //***** MCParticle collection

  // Create collection for writing into LCIO file
  LCCollectionVec* mcVec = new LCCollectionVec( LCIO::MCPARTICLE );
  TagToMCParticleMap::const_iterator ipart;

  // Setup LCCollectionVec for MCParticle table
//   prepareMCParticleCollection(mcVec);
  for( ipart = _mcPartColl->begin();
       ipart != _mcPartColl->end(); ++ipart ) {
    mcVec->push_back( ipart->second );
  }

  // Save particles into LCIO event
  lcevt->addCollection(mcVec,"MCParticle");

  // other collections to be saved
  map<string,LCCollectionVec*>::const_iterator icoll;
  for( icoll=_savedColls.begin();
       icoll!=_savedColls.end();
       ++icoll ) {
    lcevt->addCollection( icoll->second, icoll->first);
  }

  // write the event to the file
  _writer->writeEvent( lcevt );

//   // dump the event to the screen
//   G4cout<< "Dumping LCIO event:"<< G4endl;
//   LCTOOLS::dumpEvent( lcevt );

  // ------------ IMPORTANT ------------- !
  // we created the event so we need to delete it ...
  // the LCEvent destructor takes care of deleting all the hits and particles
  delete lcevt;
  reset();
}

void LCDG4LcioHelper::beginRun(const G4Run* aRun, const string& detname) {

  // Create the LCIO run header
  _runHdr = new LCRunHeaderImpl();
  _runHdr->setRunNumber( aRun->GetRunID() );

  // Need a way to pass the detector name, description and
  // subdetectors from construction class. ????
  std::string detName(detname);
  _runHdr->setDetectorName( detName );

  stringstream description;
  description << " Run: " << aRun->GetRunID()
              << ", just for testing lcio with G4";
  _runHdr->setDescription( description.str() );

  _runHdr->addActiveSubdetector( string("EMcal") );
  _runHdr->addActiveSubdetector( string("HADcal") );
  _runHdr->addActiveSubdetector( string("LUMcal") );
  _runHdr->addActiveSubdetector( string("MUONcal") );
  _runHdr->addActiveSubdetector( string("CENt") );
  _runHdr->addActiveSubdetector( string("VXDt") );

//
// Write the run header
//
  _writer->writeRunHeader( _runHdr );
}

void  LCDG4LcioHelper::endRun() {
  if(_runHdr) {
    delete _runHdr;
    _runHdr=0;
  }
}

void  LCDG4LcioHelper::reset() {
  _mcPartColl->clear();
  _savedColls.clear();
  _mcp2tag.clear();
}

void LCDG4LcioHelper::dumpParticleLists(const vector<LCDMcPart*>& lcdMcpColl) const {
  //*** Collection dumps for debugging
  G4cout<< "MCParticle list sizes: "<< lcdMcpColl.size()
        <<' '<<_mcPartColl->size() << G4endl;

  for(unsigned int itag=0; itag<lcdMcpColl.size(); ++itag) {
    //== dump LCDMcPart from tag
    LCDMcPart* plcdmc = *(lcdMcpColl.begin() + itag);
    if(plcdmc) {
      // parents
      G4cout<< "---LCDMcPart: itag="<< itag << G4endl;
      plcdmc->printParents();
      plcdmc->printDaughters();
    }

    //== dump MCParticle from tag
    TagToMCParticleMap::const_iterator imcp = _mcPartColl->find(itag);
    if( imcp==_mcPartColl->end() ) continue;
    MCParticleImpl* pmcp = imcp->second;
    // parents
    G4int npars=pmcp->getNumberOfParents();
    G4cout<< "MCParticle: "<< npars <<" parents=";
    for(int ii=0; ii<npars; ++ii) {
      MCParticleImpl* pp= dynamic_cast<MCParticleImpl*>(pmcp->getParent(ii));
      if(pp) G4cout<<" "<< _mcp2tag.find(pp)->second;
    }
    G4cout << endl;
    // daughters
    G4int ndaus=pmcp->getNumberOfDaughters();
    G4cout<< "MCParticle: "<< ndaus <<" daughters=";
    for(int ii=0; ii<ndaus; ++ii) {
      MCParticleImpl* pp= dynamic_cast<MCParticleImpl*>(pmcp->getDaughter(ii));
      if(pp) G4cout<<" "<< _mcp2tag.find(pp)->second;
    }
    G4cout << endl;
  }
}

---