LcioPrimaryGenerator.cc

Go to the documentation of this file.

// $Header: /cvs/lcd/slic/src/LcioPrimaryGenerator.cc,v 1.32 2007/04/27 01:54:35 jeremy Exp $

// slic
#include "LcioPrimaryGenerator.hh"
#include "LcioManager.hh"
#include "LcioMcpManager.hh"
#include "LcioMcpUtil.hh"
#include "StringUtil.hh"

// std
#include <sstream>

using IMPL::MCParticleImpl;
using EVENT::LCCollection;
using EVENT::LCIO;

namespace slic
{

  LcioPrimaryGenerator::LcioPrimaryGenerator(LcioManager* mgr)
    : Module( "LcioPrimaryGenerator" )
  {
    m_mgr = mgr;

    m_mcpManager = LcioMcpManager::instance();
  }

  LcioPrimaryGenerator::~LcioPrimaryGenerator()
  {}

  void LcioPrimaryGenerator::generatePrimaryVertexFromMcpCollection( LCCollection* mcpVec, G4Event* anEvent )
  {

#ifdef SLIC_LOG
    log() << LOG::debug << LOG::debug << "********** Generating Event from LCIO MCParticles **********" << LOG::done;
#endif

    assert( mcpVec );
    assert( mcpVec->getTypeName() == LCIO::MCPARTICLE );

    G4int nhep = mcpVec->getNumberOfElements();

    if ( nhep < 1 ) {
      return;
    }

#ifdef SLIC_LOG
    log() << LOG::debug << "Number of particles in Lcio input Mcp <" << nhep << ">." << LOG::done;
#endif

    for ( int i=0; i < nhep; i++ ) {

#ifdef SLIC_LOG
      log() << LOG::debug << "processing particle # <" << i << ">" << LOG::done;      
#endif

      MCParticleImpl* mcp = dynamic_cast<MCParticleImpl*>
        (mcpVec->getElementAt( i ) );

      assert( mcp );

#ifdef SLIC_LOG
      log() << LOG::debug << "pdg <" << mcp->getPDG() << ">" << LOG::done;
#endif

      G4int numPar = mcp->getParents().size();

      /*
       * Handle 4 cases, following Ron Cassell's LCSHEPEvtInterface from LCS package:
       * case 1: predecay
       * case 2: not predecay; "final state"
       * case 3: not predecay; "intermediate" or "documentation" TO BE TRACKED
       * case 4: not predecay; "intermediate" or "documentation" *NOT* TRACKED
       */
      // loop vars
      MCParticleImpl* parMcp = 0;
      G4PrimaryParticle* g4parent = 0;

      G4PrimaryParticle* thePrimary = 0;
      G4PrimaryVertex* theVertex = 0;

      G4bool isPreDecay = false;
      G4bool createPrimary = false;
      G4bool createVertex = false;

      // check if case 1 (PreDecay)
      if  ( numPar > 0 ) {
        parMcp = dynamic_cast<MCParticleImpl*>(mcp->getParents()[0]);
        g4parent = m_mcpManager->getMaps()->findPrimaryFromMcp( parMcp );

        if ( g4parent != 0 ) {

#ifdef SLIC_LOG
          log() << LOG::debug << "primary is preDecay" << LOG::done;
          std::istringstream ss;
          log() << LOG::debug << "g4parent <" << g4parent << ">" << LOG::done;
#endif

          isPreDecay = true;
        }
      }

      // case 1
      if ( isPreDecay ) {

#ifdef SLIC_LOG
        log() << LOG::debug << "case 1: predecay" << LOG::done;
#endif

        createPrimary = true;
        // no vertex
      }
      else {
        // case 2: final state
        if ( mcp->getGeneratorStatus() == 1 ) {

#ifdef SLIC_LOG
          log() << LOG::debug << "case 2: final state" << LOG::done;
#endif

          createPrimary = true;
          createVertex = true;
        }
        else {
          G4double dist = 0;

          // check case 3, if traveled > minDist
          if ( mcp->getDaughters().size() > 0 ) {

            IMPL::MCParticleImpl* firstDau = dynamic_cast<IMPL::MCParticleImpl*>
              ( mcp->getDaughters()[0] );

            dist = LcioMcpUtil::computeMcpDistance( mcp, firstDau );

#ifdef SLIC_LOG
            log() << LOG::debug << "mcpDistance=" << dist << LOG::done;
            log() << LOG::debug << "minTrackingDist=" << m_mcpManager->getMinimumTrackingDistance() << LOG::done; 
#endif

            // case 3
            if ( dist > m_mcpManager->getMinimumTrackingDistance() ) {

#ifdef SLIC_LOG
              log() << LOG::debug << "case 3: intermediate or doc TO BE TRACKED" << LOG::done;
#endif

              createPrimary = true;
              createVertex = true;
            }
            // case 4 = no-op
            // *This particle will not be tracked.*
#ifdef SLIC_LOG
            else {
              log() << LOG::debug << "case 4: intermediate or doc WILL NOT BE TRACKED" << LOG::done;
            }
#endif
          }

        }
      }

      // create a primary
      if ( createPrimary ) {

        thePrimary = createPrimaryParticleFromMcp( mcp );

        // set daughters for PreDecay
        if ( isPreDecay ) {

          assert( g4parent );

          // computation of proper_time from RC
          G4ThreeVector parMom = g4parent->GetMomentum();
          G4double E = sqrt(pow(g4parent->GetMass(), 2) + pow( parMom.x(), 2 ) + pow( parMom.y(), 2) + pow( parMom.z(), 2 ) );
          G4double proper_time = ( ( mcp->getTime() - parMcp->getTime() ) * g4parent->GetMass() ) / E;

#ifdef SLIC_LOG
          log() << LOG::debug << "parMom " << parMom << LOG::done;
          log() << LOG::debug << "e <" << E << ">" << LOG::done;
          log() << LOG::debug << "proper_time <" << proper_time << ">" << LOG::done;      
#endif

          g4parent->SetDaughter( thePrimary );
          g4parent->SetProperTime( proper_time );

#ifdef SLIC_LOG
          log() << LOG::debug << "mcp decay time <" << mcp->getTime() - parMcp->getTime() << ">" << LOG::done;    
#endif
        }
      }

      // create a vertex, add primary and set in event
      if ( createVertex ) {
        theVertex = createPrimaryVertexFromMcp( mcp );
        theVertex->SetPrimary( thePrimary );
        anEvent->AddPrimaryVertex( theVertex );
      }

      // insert mcp, primary pair into LcioManager's map (could be null)
      if ( thePrimary ) {

#ifdef SLIC_LOG
          log() << LOG::debug << "adding mcp <" << mcp << "> to primary <" << thePrimary << "> link" << LOG::done;       
#endif
        m_mcpManager->getMaps()->addMcpToPrimaryLink(mcp, thePrimary );
      }

#ifdef SLIC_LOG
      log() << LOG::debug << LOG::endl;
#endif
    }
  }

  G4PrimaryParticle* LcioPrimaryGenerator::createPrimaryParticleFromMcp(IMPL::MCParticleImpl* mcp)
  {
    assert( mcp );

    G4PrimaryParticle* primary = new G4PrimaryParticle(mcp->getPDG(),
                                                       mcp->getMomentum()[0] * GeV,
                                                       mcp->getMomentum()[1] * GeV,
                                                       mcp->getMomentum()[2] * GeV
                                                       );
    primary->SetMass( mcp->getMass() * GeV );

    return primary;
  }

  G4PrimaryVertex* LcioPrimaryGenerator::createPrimaryVertexFromMcp(IMPL::MCParticleImpl* mcp)
  {
    G4ThreeVector pos = G4ThreeVector( mcp->getVertex()[0],
                                       mcp->getVertex()[1],
                                       mcp->getVertex()[2]);

    return new G4PrimaryVertex(pos,
                               mcp->getTime() );

  }
}

---