G4HepMCInterface Class Reference

#include <G4HepMCInterface.h>

Inheritance diagram for G4HepMCInterface:

Public Member Functions
void	HepMC2G4 (HepMC::GenEvent *hepmcevt, G4Event *g4event)
void	Print (const HepMC::GenEvent *hepmcevt)
G4int	CheckType (const HepMC::GenEvent *hepmcevt)
void	Boost (HepMC::GenEvent *hepmcevt)
void	PrintHPlist ()
G4int	GetLogLevel ()
void	SetLogLevel (G4int level)
	G4HepMCInterface ()
virtual	~G4HepMCInterface ()
HepMC::GenEvent *	GetHepMCGenEvent () const
virtual void	GeneratePrimaryVertex (G4Event *anEvent)

Public Attributes
std::vector< G4HepMCParticle * >	HPlist

Protected Member Functions
virtual HepMC::GenEvent *	GenerateHepMCEvent ()

Protected Attributes
HepMC::GenEvent *	hepmcEvent
G4int	m_logLevel

Detailed Description

Definition at line 45 of file G4HepMCInterface.h.

Constructor & Destructor Documentation

G4HepMCInterface()

G4HepMCInterface::G4HepMCInterface

(

)

Definition at line 55 of file G4HepMCInterface.cpp.

  : hepmcEvent(0)
////////////////////////////////////
{
 
  IRealizationSvc *tmpReal;
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode status = svcLocator->service("RealizationSvc",tmpReal);
  if (!status.isSuccess())
  {
     std::cout << " Could not initialize Realization Service in G4HepMCInterface" << std::endl;
  } else {
    m_RealizationSvc=dynamic_cast<RealizationSvc*>(tmpReal);
  }
 
}

~G4HepMCInterface()

G4HepMCInterface::~G4HepMCInterface ( )

virtual

Definition at line 73 of file G4HepMCInterface.cpp.

{
  delete hepmcEvent;
}

Member Function Documentation

Boost()

void G4HepMCInterface::Boost

(

HepMC::GenEvent *

hepmcevt

)

Definition at line 154 of file G4HepMCInterface.cpp.

{
  //suppose relavtive vectoy is v(vx,vy,vz),position are (x,y,z,t)in CMS and(X,Y,Z,T)in Lab,the formulae for boost vertex position from Cms to Lab in natural units are following:
       // v2 = vx*vx + vy*vy + vz*vz;  
       // gamma = 1.0 / sqrt(1.0 - v2);
       // bp = vx*x + vy*y + vz*z;
       // X=x + gamma2*bp*vx + gamma*vx*t;
       // Y=y + gamma2*bp*vy + gamma*vy*t;
       // Z=z + gamma2*bp*vz + gamma*vz*t;
       // T=gamma*(t + bp);
  // the function of these formulae is the same as xvtx.boost(ThreeVector v)      
  //For the E_cms,need to loop and splice the momentum of all the final state particles
  G4LorentzVector ptot=0;
  double E_cms,p_Mag,e_Mass2,M_cms,theta;
  for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
      vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
    
    for(HepMC::GenVertex::particle_iterator
      vpitr = (*vitr)->particles_begin(HepMC::children);
    vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
      
      if((*vpitr)->status() !=1)continue;// Not final state particle
      
      G4LorentzVector p;
      p.setX( (*vpitr)-> momentum().x());
      p.setY( (*vpitr)-> momentum().y());
      p.setZ( (*vpitr)-> momentum().z());
      p.setT( (*vpitr)-> momentum().t());
      ptot = p + ptot; 
    }
  }
  E_cms=ptot.e()*GeV;
  
  //get G4Svc, allow user to access the beam momentum shift in JobOptions
    ISvcLocator* svcLocator = Gaudi::svcLocator();
  IG4Svc* tmpSvc;
  G4Svc* m_G4Svc;
  StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
  m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
  
 
  G4ThreeVector v=0;   //velocity of cms as seen from lab
  //for cms velocity 
  theta=11*mrad;     
  //theta=(m_G4Svc->GetBeamAngle())*mrad;
  //theta is half of the angle between the two beams,namely,is the angle between the beam and Z axis.
  M_cms=E_cms;  //for p_cms=0 in the CMS
  e_Mass2=electron_mass_c2*electron_mass_c2; //square of electron mass 
  p_Mag=sqrt((E_cms*E_cms-4*e_Mass2)/(2*(1-cos(pi-2*theta))));
  G4ThreeVector p_Positron(p_Mag*sin(theta),0,p_Mag*cos(theta));
  G4ThreeVector p_Electron(p_Mag*sin(pi-theta),0,p_Mag*cos(pi-theta));
  G4ThreeVector vee=p_Positron+p_Electron;
  G4ThreeVector beamshift(m_G4Svc->GetBeamShiftPx(),m_G4Svc->GetBeamShiftPy(),m_G4Svc->GetBeamShiftPz());
  if(m_G4Svc-> GetSetBeamShift()==true && fabs(M_cms-3686)<50.0) {vee = beamshift;} 
 
    double pee=vee.r(); 
    M_cms = sqrt(M_cms*M_cms + pee*pee);
    v=vee/M_cms;
 
  
  for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
      vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
    
    // Boost vertex position from cms to lab 
    G4LorentzVector xvtx;
    xvtx.setX( (*vitr)-> position().x());
    xvtx.setY( (*vitr)-> position().y());
    xvtx.setZ( (*vitr)-> position().z());
    xvtx.setT( (*vitr)-> position().t());
    xvtx.boost(v);
    (*vitr)->set_position(xvtx);
    
    // Boost the particles four momentum from cms to lab
    // the transform formulae are similary to the position boost       
    for (HepMC::GenVertex::particle_iterator 
       vpitr = (*vitr)->particles_begin(HepMC::children);
     vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
      G4LorentzVector p;
      p.setX( (*vpitr)-> momentum().x());
      p.setY( (*vpitr)-> momentum().y());
      p.setZ( (*vpitr)-> momentum().z());
      p.setT( (*vpitr)-> momentum().t());
      p=p.boost(v);
      (*vpitr)->set_momentum(p);
    }
  }
}

Referenced by HepMC2G4().

CheckType()

G4int G4HepMCInterface::CheckType

(

const HepMC::GenEvent *

hepmcevt

)

Definition at line 129 of file G4HepMCInterface.cpp.

{
  G4int flag =0;
  for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
          vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
    for (HepMC::GenVertex::particle_iterator
             pitr= (*vitr)->particles_begin(HepMC::children);
                pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
      if((*pitr)->end_vertex())
        {flag = 1; break;}
    }
    if(flag) break;
  }
  if(m_logLevel <= 4)
  {
    if(flag == 0)
      G4cout<<G4endl<<"generator is GENBES!"<<G4endl;
    else
      G4cout<<G4endl<<"generator is not GENBES!"<<G4endl;
  }
  return flag;
}

Referenced by HepMC2G4().

GenerateHepMCEvent()

HepMC::GenEvent * G4HepMCInterface::GenerateHepMCEvent ( )

protectedvirtual

Reimplemented in BesHepMCInterface.

Definition at line 465 of file G4HepMCInterface.cpp.

{
  HepMC::GenEvent* aevent= new HepMC::GenEvent();
  return aevent;
}

Referenced by GeneratePrimaryVertex().

GeneratePrimaryVertex()

void G4HepMCInterface::GeneratePrimaryVertex ( G4Event *  anEvent )

virtual

Definition at line 473 of file G4HepMCInterface.cpp.

{
  // delete previous event object
  delete hepmcEvent;
 
  // generate next event
  hepmcEvent= GenerateHepMCEvent();
  if(! hepmcEvent) {
    G4Exception("G4HepMCInterface","GeneratePrimaryVertex",RunMustBeAborted,
              "HepMCInterface: no generated particles. run terminated..." );
    return;
  }
  HepMC2G4(hepmcEvent, anEvent);
}

GetHepMCGenEvent()

HepMC::GenEvent * G4HepMCInterface::GetHepMCGenEvent ( ) const

inline

Definition at line 94 of file G4HepMCInterface.h.

{ 
  return hepmcEvent; 
}

GetLogLevel()

G4int G4HepMCInterface::GetLogLevel ( )

inline

Definition at line 55 of file G4HepMCInterface.h.

{return m_logLevel;}

HepMC2G4()

void G4HepMCInterface::HepMC2G4	(	HepMC::GenEvent *	hepmcevt,
		G4Event *	g4event
	)

Definition at line 244 of file G4HepMCInterface.cpp.

{
  //Print HepMC Event before boost
  if(m_logLevel <=4 )
    Print(hepmcevt);
 
  //get G4Svc
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  IG4Svc* tmpSvc;
  G4Svc* m_G4Svc;
  StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
  m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc); 
 
  //considering 22rmad beam clossing angle ,need to boost the Cms to Lab 
  // need to boost ?
  if(m_G4Svc->GetBoostLab()== true) 
    Boost(hepmcevt); 
  
  //Print HepMC Event after boost
  if(m_logLevel <=2  && m_G4Svc->GetBoostLab()== true)
    Print(hepmcevt);
 
  //*********************send particles to HPlist*************************//
  G4int flag = CheckType(hepmcevt);
  if(flag)
  {
    for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
        vitr != hepmcevt->vertices_end(); ++vitr ) 
    { // loop for vertex ...
      G4int vtxFlag=0;
      G4int pOut = (*vitr)->particles_out_size();
      HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
      G4int pdgcode= (*pitr)-> pdg_id();
      if(pOut == 1 && abs(pdgcode) == 11)
        vtxFlag=1;
      G4int barcodeVtx = (*vitr)->barcode();
 
      for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
           pitr != (*vitr)->particles_end(HepMC::children); ++pitr) 
      {// loop for particle from this vertex 
        G4int pdgcode = (*pitr)->pdg_id();
        if(vtxFlag==0)
        {
          if(pdgcode==-22) pdgcode=22;
          G4LorentzVector p;
          p.setX( (*pitr)-> momentum().x());
          p.setY( (*pitr)-> momentum().y());
          p.setZ( (*pitr)-> momentum().z());    
          p.setT( (*pitr)-> momentum().t());
          G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
          //G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode);
          //particle->Set4Momentum(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
          G4int ISTHEP = (*pitr)->status();
  
          G4int barcodeEndVtx = 0;
          if((*pitr)->end_vertex())
            barcodeEndVtx = (*pitr)->end_vertex()->barcode();
  
          G4HepMCParticle* hepmcParticle = new G4HepMCParticle(particle,ISTHEP,barcodeEndVtx);
          HPlist.push_back(hepmcParticle);
          
          for( size_t ii=0; ii<HPlist.size(); ii++ )
          {  
            if(barcodeVtx == HPlist[ii]->GetBarcodeEndVtx())
            {
              HPlist[ii]->GetTheParticle()->SetDaughter(particle); 
              hepmcParticle->Done();
              break; 
            }
          }
        }
      }
    }
  }
  else 
  {
    for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
                vitr != hepmcevt->vertices_end(); ++vitr )
    {
      for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
                     pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
      {
        G4int ISTHEP = (*pitr)->status();
        G4LorentzVector p;
        p.setX( (*pitr)-> momentum().x());
        p.setY( (*pitr)-> momentum().y());
        p.setZ( (*pitr)-> momentum().z());
        p.setT( (*pitr)-> momentum().t());
 
        G4int pdgcode = (*pitr)->pdg_id();
        G4int barcodeEndVtx = 0;
        if((*pitr)->end_vertex())
          barcodeEndVtx = (*pitr)->end_vertex()->barcode();
        G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
        //G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode);
        //particle->Set4Momentum(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
 
        G4HepMCParticle* hepmcParticle = new G4HepMCParticle(particle,ISTHEP,barcodeEndVtx);
        HPlist.push_back(hepmcParticle);
        if(ISTHEP>1)
          hepmcParticle->Done();
      }
    }
  }
 
  //print particles in HPlist
  if(m_logLevel<=4)
    PrintHPlist();
 
  //get time and position information from G4Svc
  G4double pmPosX,pmPosY,pmPosZ,pmTime;
  G4double tmpPosX,tmpPosY,tmpPosZ,tmpT;
  G4double beamPosX,beamPosY,beamPosZ,beamSizeX,beamSizeY,beamSizeZ;
 
  if(m_RealizationSvc->UseDBFlag() == false) {
    beamPosX = m_G4Svc->GetBeamPosX();
    beamPosY = m_G4Svc->GetBeamPosY();
    beamPosZ = m_G4Svc->GetBeamPosZ();
 
    beamSizeX = m_G4Svc->GetBeamSizeX();
    beamSizeY = m_G4Svc->GetBeamSizeY();
    beamSizeZ = m_G4Svc->GetBeamSizeZ()/sqrt(2);
  } 
  if(m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == true) {
    beamPosX = m_RealizationSvc->getBunchPosX();
    beamPosY = m_RealizationSvc->getBunchPosY();
    beamPosZ = m_RealizationSvc->getBunchPosZ();
 
    beamSizeX = m_RealizationSvc->getBunchSizeX();
    beamSizeY = m_RealizationSvc->getBunchSizeY();
    beamSizeZ = m_RealizationSvc->getBunchSizeZ();
  }
  if(m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == false) {
    beamPosX = m_G4Svc->GetBeamPosX();
    beamPosY = m_G4Svc->GetBeamPosY();
    beamPosZ = m_G4Svc->GetBeamPosZ();
 
    beamSizeX = m_G4Svc->GetBeamSizeX();
    beamSizeY = m_G4Svc->GetBeamSizeY();
    beamSizeZ = m_G4Svc->GetBeamSizeZ()/sqrt(2);
  }
 
  G4double gaussX = G4RandGauss::shoot();
  G4double gaussY = G4RandGauss::shoot();
  G4double gaussZ = G4RandGauss::shoot();
  G4double gaussT = G4RandGauss::shoot();
 
  G4double beamStartTime = m_G4Svc->GetBeamStartTime() * ns;
  G4double beamDeltaTime = m_G4Svc->GetBeamDeltaTime() * ns;
  G4double nBunch = m_G4Svc->GetNBunch();
  G4double bunchTimeSigma = m_G4Svc->GetBunchTimeSigma() * ns;
  
  G4double ran=G4UniformRand();
  G4double beamTime = bunchTimeSigma*G4RandGauss::shoot() + beamStartTime + beamDeltaTime*int(ran*nBunch); 
 
  tmpPosX = (beamPosX + beamSizeX*gaussX ) *mm;
  tmpPosY = (beamPosY + beamSizeY*gaussY ) *mm;
  tmpPosZ = (beamPosZ + beamSizeZ*gaussZ ) *mm;
  tmpT = (beamSizeZ * gaussT ) * mm/c_light +beamTime;
 
  G4LorentzVector tmpv(tmpPosX,tmpPosY,tmpPosZ,tmpT*c_light/mm);
 
  HepMC::GenEvent::vertex_const_iterator vitr0= hepmcevt->vertices_begin();
  G4LorentzVector xvtx0 ;
  xvtx0.setX( (*vitr0)-> position().x());
  xvtx0.setY( (*vitr0)-> position().y());
  xvtx0.setZ( (*vitr0)-> position().z());
  xvtx0.setT( (*vitr0)-> position().t());
  pmPosX = xvtx0.x()*mm + tmpPosX;
  pmPosY = xvtx0.y()*mm + tmpPosY;
  pmPosZ = xvtx0.z()*mm + tmpPosZ;
  pmTime = xvtx0.t()*mm/c_light + tmpT;
 
  if(m_logLevel<=4)
  {
    G4cout<<G4endl;
    G4cout<<xvtx0.x()<<" "<<xvtx0.y()<<" "<<xvtx0.z()<<" "
          <<beamSizeX*gaussX<<" "
          <<beamSizeY*gaussY<<" "
          <<beamSizeZ*gaussZ<<" "<<G4endl;
    G4cout<<xvtx0.t()* mm/c_light<<" "
          <<beamSizeZ * gaussT/sqrt(2)*mm/c_light<<" "
          <<beamTime<<G4endl;
  }
  for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
              vitr != hepmcevt->vertices_end(); ++vitr )
  {
    G4LorentzVector xvtx;
    xvtx.setX((*vitr)-> position().x());
    xvtx.setY((*vitr)-> position().y());
    xvtx.setZ((*vitr)-> position().z());
    xvtx.setT((*vitr)-> position().t());
    (*vitr)->set_position(xvtx+tmpv);
  }
  m_G4Svc->SetBeamTime(pmTime);
 
  G4PrimaryVertex* g4vtx= new G4PrimaryVertex(pmPosX,pmPosY,pmPosZ,pmTime);
 
  // put initial particles to the vertex
  for( size_t ii=0; ii<HPlist.size(); ii++ )
  {
    if( HPlist[ii]->GetISTHEP() > 0 ) // ISTHEP of daughters had been 
                                      // set to negative
    {
      G4PrimaryParticle* initialParticle = HPlist[ii]->GetTheParticle();
      g4vtx->SetPrimary( initialParticle );
    }
  }
  
  //clear G4HepMCInterface
  for(size_t iii=0;iii<HPlist.size();iii++)
  { delete HPlist[iii]; }
  HPlist.clear();
 
  g4event->AddPrimaryVertex(g4vtx);
} 

Referenced by G4SvcRunManager::GenerateEvent(), and GeneratePrimaryVertex().

Print()

void G4HepMCInterface::Print

(

const HepMC::GenEvent *

hepmcevt

)

Definition at line 78 of file G4HepMCInterface.cpp.

{
  G4int i=0;
  for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
      vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
    G4cout<<G4endl<<"vertex:"<<i<<"  barcode:"<<(*vitr)->barcode()<<" "; 
    i++;
    //if((*vitr)->mother())
    //  G4cout<<"mother particle: "<<(*vitr)->mother()-> pdg_id()<<"  ";
    //if((*vitr)->secondMother())
    //  G4cout<<" second mother: "<<(*vitr)->secondMother()-> pdg_id()<<G4endl;
    G4LorentzVector xvtx;
    xvtx.setX( (*vitr)-> position().x());
    xvtx.setY( (*vitr)-> position().y());
    xvtx.setZ( (*vitr)-> position().z());
    xvtx.setT( (*vitr)-> position().t());    
    G4cout<<"x:"<<xvtx.x()<<" y:"<<xvtx.y()<<" z:"<<xvtx.z()
          <<" t:"<<xvtx.t()*mm/c_light<<G4endl;
    
    G4int j=0;
    for (HepMC::GenVertex::particle_iterator 
       pitr= (*vitr)->particles_begin(HepMC::children);
     pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
      G4int pdgcode= (*pitr)-> pdg_id();
      G4LorentzVector p;
      p.setX( (*pitr)-> momentum().x());
      p.setY( (*pitr)-> momentum().y());
      p.setZ( (*pitr)-> momentum().z());
      p.setT( (*pitr)-> momentum().t());
      G4cout<<"particle:"<<j<<" pdgcode:"<<pdgcode<<" "; 
      G4cout<<"p:"<<p.x()<<" "<<p.y()<<" "<<p.z()<<" ";
      if((*pitr)->end_vertex()) 
        G4cout<<"endVtx:"<<(*pitr)->end_vertex()->barcode()<<" ";
      G4cout<<"status:"<<(*pitr)->status()<<G4endl;
      j++;
    }
  }
}

Referenced by HepMC2G4().

PrintHPlist()

void G4HepMCInterface::PrintHPlist

(

)

Definition at line 117 of file G4HepMCInterface.cpp.

{
  G4cout<<G4endl;
  G4cout<<"particles sent to Geant4: "<<G4endl;
  for( size_t ii=0; ii<HPlist.size(); ii++ )
  {  
    G4int pdgcode = HPlist[ii]->GetTheParticle()->GetPDGcode(); 
    G4cout<<pdgcode<<" ";
  }
  G4cout<<G4endl;
}

Referenced by HepMC2G4().

SetLogLevel()

void G4HepMCInterface::SetLogLevel ( G4int  level )

inline

Definition at line 56 of file G4HepMCInterface.h.

{m_logLevel = level;}

Referenced by G4SvcRunManager::GenerateEvent().

Member Data Documentation

hepmcEvent

HepMC::GenEvent* G4HepMCInterface::hepmcEvent

protected

Definition at line 66 of file G4HepMCInterface.h.

Referenced by GeneratePrimaryVertex(), GetHepMCGenEvent(), and ~G4HepMCInterface().

HPlist

std::vector<G4HepMCParticle*> G4HepMCInterface::HPlist

Definition at line 84 of file G4HepMCInterface.h.

Referenced by HepMC2G4(), and PrintHPlist().

m_logLevel

G4int G4HepMCInterface::m_logLevel

protected

Definition at line 71 of file G4HepMCInterface.h.

Referenced by CheckType(), GetLogLevel(), HepMC2G4(), and SetLogLevel().

---

The documentation for this class was generated from the following files:

• G4HepMCInterface.h
• G4HepMCInterface.cpp