eic-opticks/u4/U4VPrimaryGenerator.h

0001 #pragma once
0002 /**
0003 U4VPrimaryGenerator.h
0004 =======================
0005
0006 Implemented header only to allow isolating use
0007 of MOCK_CURAND to the test executable only, and
0008 not the library.
0009
0010
0011 **/
0012
0013 struct NP ;
0014 struct sphoton ;
0015 class G4PrimaryVertex ;
0016 class G4Event ;
0017 #include "U4_API_EXPORT.hh"
0018
0019 struct U4_API U4VPrimaryGenerator
0020 {
0021     template<typename P>
0022     static void GetPhotonParam(
0023          G4ThreeVector& position_mm, G4double& time_ns,
0024          G4ThreeVector& direction,  G4double& wavelength_nm,
0025          G4ThreeVector& polarization, const P& p );
0026
0027     template<typename P>
0028     static G4PrimaryVertex* MakePrimaryVertexPhoton( const P& p);
0029
0030     static constexpr const char* _DEBUG_GENIDX = "U4VPrimaryGenerator__GeneratePrimaries_From_Photons_DEBUG_GENIDX" ;
0031     static void GeneratePrimaries_From_Photons(G4Event* event, const NP* ph );
0032
0033 };
0034
0035
0036 #include <cassert>
0037
0038 #include "G4Event.hh"
0039 #include "G4OpticalPhoton.hh"
0040 #include "G4PrimaryVertex.hh"
0041 #include "G4PrimaryParticle.hh"
0042 #include "G4PhysicalConstants.hh"
0043 #include "G4SystemOfUnits.hh"
0044
0045 #include "ssys.h"
0046 #include "scuda.h"
0047 #include "squad.h"
0048 #include "sphoton.h"
0049 #include "NP.hh"
0050
0051 /**
0052 U4VPrimaryGenerator::GetPhotonParam
0053 -------------------------------------
0054
0055 Photon parameters from Opticks photon types such as sphoton.h shuffled into Geant4 three vectors.
0056
0057 **/
0058
0059 template<typename P>
0060 inline void U4VPrimaryGenerator::GetPhotonParam(
0061      G4ThreeVector& position_mm, G4double& time_ns,
0062      G4ThreeVector& direction,  G4double& wavelength_nm,
0063      G4ThreeVector& polarization, const P& p )
0064 {
0065      position_mm.set(p.pos.x, p.pos.y, p.pos.z);
0066      time_ns = p.time ;
0067
0068      direction.set(p.mom.x, p.mom.y, p.mom.z );
0069      direction = direction.unit();
0070
0071      polarization.set(p.pol.x, p.pol.y, p.pol.z );
0072      wavelength_nm = p.wavelength ;
0073 }
0074
0075
0076 /**
0077 U4VPrimaryGenerator::MakePrimaryVertexPhoton
0078 ----------------------------------------------
0079
0080 Converts Opticks photon type P (eg sphoton.h) into G4PrimaryVertex
0081
0082 1. shuffle photon param from P into G4ThreeVector, G4double
0083 2. create G4PrimaryVertex and populate with the param obtained above
0084
0085 **/
0086
0087 template<typename P>
0088 inline G4PrimaryVertex* U4VPrimaryGenerator::MakePrimaryVertexPhoton( const P& p)
0089 {
0090     G4ThreeVector position_mm ;
0091     G4double time_ns  ;
0092     G4ThreeVector direction ;
0093     G4double wavelength_nm ;
0094     G4ThreeVector polarization;
0095
0096     GetPhotonParam( position_mm, time_ns, direction, wavelength_nm, polarization, p );
0097
0098     G4PrimaryVertex* vertex = new G4PrimaryVertex(position_mm, time_ns);
0099     G4double kineticEnergy = h_Planck*c_light/(wavelength_nm*nm) ;
0100     G4PrimaryParticle* particle = new G4PrimaryParticle(G4OpticalPhoton::Definition());
0101     particle->SetKineticEnergy( kineticEnergy );
0102     particle->SetMomentumDirection( direction );
0103     particle->SetPolarization(polarization);
0104
0105     vertex->SetPrimary(particle);
0106     return vertex ;
0107 }
0108
0109 /**
0110 U4VPrimaryGenerator::GeneratePrimaries
0111 ---------------------------------------
0112
0113 1. populates the G4Event argument with photon parameters from ph array tas G4PrimaryVertex
0114
0115 Notice that there are no G4Track in sight here, so there is no
0116 direct way to annotate the tracks with *spho* labels.
0117
0118 To do so would have to arrange to catch the tracks at labelling stage
0119 in U4Recorder::PreUserTrackingAction_Optical where could
0120 associate back to the originating array, that could be held in the
0121 SEvt as "input photons". This relies on Geant4 being consistent
0122 in the way PrimaryVertex become G4Track ? I expect it will work
0123 in purely optical case.
0124
0125 HMM: These *ph* are effectively input photons (even though generated from gensteps),
0126 should associate as such in the SEvt to retain access to these
0127
0128 **/
0129
0130
0131
0132 inline void U4VPrimaryGenerator::GeneratePrimaries_From_Photons(G4Event* event, const NP* ph )
0133 {
0134     int DEBUG_GENIDX = ssys::getenvint(_DEBUG_GENIDX, -1 ) ;
0135
0136     std::cout
0137          << "U4VPrimaryGenerator::GeneratePrimaries_From_Photons"
0138          << " ph " << ( ph ? ph->sstr() : "-" )
0139          << std::endl
0140          << " " << _DEBUG_GENIDX << " : " << DEBUG_GENIDX << " (when +ve, only generate tht photon idx)"
0141          << std::endl
0142          ;
0143
0144     if(ph == nullptr) std::cerr
0145          << "U4VPrimaryGenerator::GeneratePrimaries : FATAL : NO PHOTONS " << std::endl
0146          << "compile with MOCK_CURAND to use SGenerate.h curand on CPU" << std::endl
0147          ;
0148     if(ph == nullptr) return ;
0149
0150     //std::cout << "U4VPrimaryGenerator::GeneratePrimaries" << " ph " << ( ph ? ph->brief() : "-" ) << std::endl ;
0151
0152     if( ph->ebyte == 4 )
0153     {
0154         sphoton* pp = (sphoton*)ph->bytes() ;
0155         for(int i=0 ; i < ph->shape[0] ; i++)
0156         {
0157             if( DEBUG_GENIDX > -1 && i != DEBUG_GENIDX ) continue ;
0158             const sphoton& p = pp[i];
0159             G4PrimaryVertex* vertex = MakePrimaryVertexPhoton<sphoton>( p );
0160             event->AddPrimaryVertex(vertex);
0161         }
0162     }
0163     else if( ph->ebyte == 8 )
0164     {
0165         sphotond* pp = (sphotond*)ph->bytes() ;
0166         for(int i=0 ; i < ph->shape[0] ; i++)
0167         {
0168             if( DEBUG_GENIDX > -1 && i != DEBUG_GENIDX ) continue ;
0169             const sphotond& p = pp[i];
0170             G4PrimaryVertex* vertex = MakePrimaryVertexPhoton<sphotond>( p );
0171             event->AddPrimaryVertex(vertex);
0172         }
0173     }
0174 }
0175