eic-opticks/u4/U4.cc

0001 #include <iomanip>
0002 #include <iostream>
0003 #include <cassert>
0004 #include <csignal>
0005 #include <cstdlib>
0006
0007 #include "G4VParticleChange.hh"
0008 #include "G4SystemOfUnits.hh"
0009 #include "G4PhysicalConstants.hh"
0010 #include "G4Track.hh"
0011 #include "G4OpticalPhoton.hh"
0012 #include "G4Event.hh"
0013
0014 #include "SEvt.hh"
0015 #include "scuda.h"
0016 #include "squad.h"
0017 #include "sphoton.h"
0018 #include "OpticksGenstep.h"
0019
0020 #include "SPath.hh"
0021 #include "SStr.hh"
0022 #include "NP.hh"
0023 #include "SLOG.hh"
0024 #include "sscint.h"
0025 #include "scerenkov.h"
0026 #include "sgs.h"
0027
0028 #ifdef WITH_CUSTOM4
0029 #include "C4GS.h"
0030 #include "C4Pho.h"
0031 #include "C4TrackInfo.h"
0032 #else
0033 #include "spho.h"
0034 #include "STrackInfo.h"
0035 #endif
0036
0037 #include "U4.hh"
0038
0039 const plog::Severity U4::LEVEL = SLOG::EnvLevel("U4", "DEBUG");
0040
0041
0042 /**
0043 U4 private state
0044 --------------------
0045
0046 Here are holding the state of the genstep collection in translation-unit-local static variables.
0047
0048 An alternative more o.o. approach would be to use a U4Private/U4Impl struct
0049 that a U4 instance holds a pointer to and passes along messages to.
0050 That is like the PImpl pattern : pointer to implementation.
0051
0052 * https://www.geeksforgeeks.org/pimpl-idiom-in-c-with-examples/
0053 * https://www.cppstories.com/2018/01/pimpl/
0054
0055 HMM: perhapa this state belongs better within SEvt together with the full gensteps ?
0056
0057 **/
0058
0059
0060 #ifdef WITH_CUSTOM4
0061 static C4GS gs = {} ;            // updated by eg U4::CollectGenstep_DsG4Scintillation_r4695 prior to each photon generation loop
0062 static C4Pho ancestor = {} ;     // updated by U4::GenPhotonAncestor prior to the photon generation loop(s)
0063 static C4Pho pho = {} ;          // updated by U4::GenPhotonBegin at start of photon generation loop
0064 static C4Pho secondary = {} ;    // updated by U4::GenPhotonEnd   at end of photon generation loop
0065 #else
0066 static sgs gs = {} ;            // updated by eg U4::CollectGenstep_DsG4Scintillation_r4695 prior to each photon generation loop
0067 static spho ancestor = {} ;     // updated by U4::GenPhotonAncestor prior to the photon generation loop(s)
0068 static spho pho = {} ;          // updated by U4::GenPhotonBegin at start of photon generation loop
0069 static spho secondary = {} ;    // updated by U4::GenPhotonEnd   at end of photon generation loop
0070 #endif
0071
0072 static bool dump = false ;
0073
0074 /**
0075 MakeGenstep... Hidden Functions only usable from this translation unit
0076 ------------------------------------------------------------------------
0077
0078 Using hidden static non-member functions allows keeping Opticks types like quad6 out of U4 header
0079
0080 **/
0081
0082 static quad6 MakeGenstep_DsG4Scintillation_r4695(
0083      const G4Track* aTrack,
0084      const G4Step* aStep,
0085      G4int    numPhotons,
0086      G4int    scnt,
0087      G4double ScintillationTime
0088     )
0089 {
0090     G4StepPoint* pPreStepPoint  = aStep->GetPreStepPoint();
0091     G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
0092
0093     G4ThreeVector x0 = pPreStepPoint->GetPosition();
0094     G4double      t0 = pPreStepPoint->GetGlobalTime();
0095     G4ThreeVector deltaPosition = aStep->GetDeltaPosition() ;
0096     G4double meanVelocity = (pPreStepPoint->GetVelocity()+pPostStepPoint->GetVelocity())/2. ;
0097
0098     const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
0099     const G4Material* aMaterial = aTrack->GetMaterial();
0100
0101     quad6 _gs ;
0102     _gs.zero() ;
0103
0104     sscint* gs = (sscint*)(&_gs) ;   // warning: dereferencing type-punned pointer will break strict-aliasing rules
0105
0106     gs->gentype = OpticksGenstep_DsG4Scintillation_r4695 ;
0107     gs->trackid = aTrack->GetTrackID() ;
0108     gs->matline = aMaterial->GetIndex() + SEvt::G4_INDEX_OFFSET ;  // offset signals that a mapping must be done in SEvt::setGenstep
0109     gs->numphoton = numPhotons ;
0110
0111     // note that gs->matline is not currently used for scintillation,
0112     // but done here as check of SEvt::addGenstep mtindex to mtline mapping
0113
0114     gs->pos.x = x0.x() ;
0115     gs->pos.y = x0.y() ;
0116     gs->pos.z = x0.z() ;
0117     gs->time = t0 ;
0118
0119     gs->DeltaPosition.x = deltaPosition.x() ;
0120     gs->DeltaPosition.y = deltaPosition.y() ;
0121     gs->DeltaPosition.z = deltaPosition.z() ;
0122     gs->step_length = aStep->GetStepLength() ;
0123
0124     gs->code = aParticle->GetDefinition()->GetPDGEncoding() ;
0125     gs->charge = aParticle->GetDefinition()->GetPDGCharge() ;
0126     gs->weight = aTrack->GetWeight() ;
0127     gs->meanVelocity = meanVelocity ;
0128
0129     gs->scnt = scnt ;
0130     gs->f41 = 0.f ;
0131     gs->f42 = 0.f ;
0132     gs->f43 = 0.f ;
0133
0134     gs->ScintillationTime = ScintillationTime ;
0135     gs->f51 = 0.f ;
0136     gs->f52 = 0.f ;
0137     gs->f53 = 0.f ;
0138
0139     return _gs ;
0140 }
0141
0142
0143 const char* U4::CollectGenstep_DsG4Scintillation_r4695_DISABLE = "U4__CollectGenstep_DsG4Scintillation_r4695_DISABLE" ;
0144 const char* U4::CollectGenstep_DsG4Scintillation_r4695_ZEROPHO = "U4__CollectGenstep_DsG4Scintillation_r4695_ZEROPHO" ;
0145
0146
0147 /**
0148
0149 U4::CollectGenstep_DsG4Scintillation_r4695
0150 -------------------------------------------
0151
0152 1. makes quad6 gs_ collecting parameters needed to generate photons on GPU
0153 2. adds gs_ to SEvt using SEvt::AddGenstep which collects into a vectors of quad6 for each SEvt
0154    instance. For standard GPU running there is just one GPU instance, for GPU CPU comparison
0155    running that there will be two SEvt instances.
0156
0157 **/
0158
0159 void U4::CollectGenstep_DsG4Scintillation_r4695(
0160          const G4Track* aTrack,
0161          const G4Step* aStep,
0162          G4int    numPhotons,
0163          G4int    scnt,
0164          G4double ScintillationTime
0165     )
0166 {
0167     if(getenv(CollectGenstep_DsG4Scintillation_r4695_DISABLE))
0168     {
0169         LOG(error) << CollectGenstep_DsG4Scintillation_r4695_DISABLE ;
0170         return ;
0171     }
0172     if(getenv(CollectGenstep_DsG4Scintillation_r4695_ZEROPHO))
0173     {
0174         LOG(error) << CollectGenstep_DsG4Scintillation_r4695_ZEROPHO ;
0175         numPhotons = 0 ;
0176     }
0177
0178
0179
0180     quad6 gs_ = MakeGenstep_DsG4Scintillation_r4695( aTrack, aStep, numPhotons, scnt, ScintillationTime);
0181
0182 #ifdef WITH_CUSTOM4
0183     sgs _gs = SEvt::AddGenstep(gs_);    // returns sgs struct which is a simple 4 int label
0184     gs = C4GS::Make(_gs.index, _gs.photons, _gs.offset, _gs.gentype );
0185 #else
0186     gs = SEvt::AddGenstep(gs_);    // returns sgs struct which is a simple 4 int label
0187 #endif
0188     // gs is private static genstep label
0189
0190     //if(dump) std::cout << "U4::CollectGenstep_DsG4Scintillation_r4695 " << gs.desc() << std::endl ;
0191     LOG(LEVEL) << gs.desc();
0192 }
0193
0194
0195
0196 static quad6 MakeGenstep_G4Cerenkov_modified(
0197     const G4Track* aTrack,
0198     const G4Step* aStep,
0199     G4int    numPhotons,
0200     G4double    betaInverse,
0201     G4double    pmin,
0202     G4double    pmax,
0203     G4double    maxCos,
0204
0205     G4double    maxSin2,
0206     G4double    meanNumberOfPhotons1,
0207     G4double    meanNumberOfPhotons2
0208     )
0209 {
0210     G4StepPoint* pPreStepPoint  = aStep->GetPreStepPoint();
0211     G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
0212
0213     G4ThreeVector x0 = pPreStepPoint->GetPosition();
0214     G4double      t0 = pPreStepPoint->GetGlobalTime();
0215     G4ThreeVector deltaPosition = aStep->GetDeltaPosition() ;
0216
0217     const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
0218
0219     G4double Wmin_nm = h_Planck*c_light/pmax/nm ;
0220     G4double Wmax_nm = h_Planck*c_light/pmin/nm ;
0221
0222     const G4Material* aMaterial = aTrack->GetMaterial();
0223
0224     quad6 _gs ;
0225     _gs.zero() ;
0226
0227     scerenkov* gs = (scerenkov*)(&_gs) ;
0228
0229     gs->gentype = OpticksGenstep_G4Cerenkov_modified ;
0230     gs->trackid = aTrack->GetTrackID() ;
0231     gs->matline = aMaterial->GetIndex() + SEvt::G4_INDEX_OFFSET ;  // offset signals that a mapping must be done in SEvt::setGenstep
0232     gs->numphoton = numPhotons ;
0233
0234     gs->pos.x = x0.x() ;
0235     gs->pos.y = x0.y() ;
0236     gs->pos.z = x0.z() ;
0237     gs->time = t0 ;
0238
0239     gs->DeltaPosition.x = deltaPosition.x() ;
0240     gs->DeltaPosition.y = deltaPosition.y() ;
0241     gs->DeltaPosition.z = deltaPosition.z() ;
0242     gs->step_length = aStep->GetStepLength() ;
0243
0244     gs->code = aParticle->GetDefinition()->GetPDGEncoding() ;
0245     gs->charge = aParticle->GetDefinition()->GetPDGCharge() ;
0246     gs->weight = aTrack->GetWeight() ;
0247     gs->preVelocity = pPreStepPoint->GetVelocity() ;
0248
0249     gs->BetaInverse = betaInverse ;
0250     gs->Wmin = Wmin_nm ;
0251     gs->Wmax = Wmax_nm   ;
0252     gs->maxCos = maxCos ;
0253
0254     gs->maxSin2 = maxSin2 ;
0255     gs->MeanNumberOfPhotons1 = meanNumberOfPhotons1 ;
0256     gs->MeanNumberOfPhotons2 = meanNumberOfPhotons2 ;
0257     gs->postVelocity = pPostStepPoint->GetVelocity() ;
0258
0259     return _gs ;
0260 }
0261
0262 const char* U4::CollectGenstep_G4Cerenkov_modified_DISABLE = "U4__CollectGenstep_G4Cerenkov_modified_DISABLE" ;
0263 const char* U4::CollectGenstep_G4Cerenkov_modified_ZEROPHO = "U4__CollectGenstep_G4Cerenkov_modified_ZEROPHO" ;
0264
0265 /**
0266 U4::CollectGenstep_G4Cerenkov_modified
0267 ---------------------------------------
0268
0269 1. makes quad6 gs_ collecting parameters needed to generate photons on GPU
0270 2. adds gs_ to SEvt using SEvt::AddGenstep which collects into a vectors of quad6 for each SEvt
0271    instance. For standard GPU running there is just one GPU instance, for GPU CPU comparison
0272    running that there will be two SEvt instances.
0273
0274
0275 **/
0276
0277
0278 void U4::CollectGenstep_G4Cerenkov_modified(
0279     const G4Track* aTrack,
0280     const G4Step* aStep,
0281     G4int    numPhotons,
0282     G4double    betaInverse,
0283     G4double    pmin,
0284     G4double    pmax,
0285     G4double    maxCos,
0286     G4double    maxSin2,
0287     G4double    meanNumberOfPhotons1,
0288     G4double    meanNumberOfPhotons2)
0289 {
0290
0291     if(getenv(CollectGenstep_G4Cerenkov_modified_DISABLE))
0292     {
0293         LOG(error) << CollectGenstep_G4Cerenkov_modified_DISABLE ;
0294         return ;
0295     }
0296     if(getenv(CollectGenstep_G4Cerenkov_modified_ZEROPHO))
0297     {
0298         LOG(error) << CollectGenstep_G4Cerenkov_modified_ZEROPHO ;
0299         numPhotons = 0 ;
0300     }
0301
0302
0303
0304
0305     quad6 gs_ = MakeGenstep_G4Cerenkov_modified( aTrack, aStep, numPhotons, betaInverse, pmin, pmax, maxCos, maxSin2, meanNumberOfPhotons1, meanNumberOfPhotons2 );
0306
0307 #ifdef WITH_CUSTOM4
0308     sgs _gs = SEvt::AddGenstep(gs_);    // returns sgs struct which is a simple 4 int label
0309     gs = C4GS::Make(_gs.index, _gs.photons, _gs.offset , _gs.gentype );
0310 #else
0311     gs = SEvt::AddGenstep(gs_);    // returns sgs struct which is a simple 4 int label
0312 #endif
0313     // gs is primate static genstep label
0314     // TODO: avoid the duplication betweek C and S with common SetGenstep private method
0315
0316     if(dump) std::cout << "U4::CollectGenstep_G4Cerenkov_modified " << gs.desc() << std::endl ;
0317     LOG(LEVEL) << gs.desc();
0318 }
0319
0320
0321
0322
0323
0324
0325
0326 /**
0327 U4::GenPhotonAncestor
0328 ----------------------
0329
0330 NB this is called prior to generation loops to get the ancestor spho.h label
0331
0332 This label is needed for BOTH Scintillation and Cerenkov in order for photon G4Track
0333 labelling done by U4::GenPhotonEnd to work.
0334
0335 When the track has no user info the ancestor is set to spho::Placeholder label {-1,-1,-1,-1}.
0336
0337 If the call to U4::GenPhotonAncestor is omitted from the Scintillation OR Cerenkov
0338 PostStepDoIt then the ancestor will default to {0,0,0,0} : that will cause
0339 unexpected labels.
0340
0341 **/
0342
0343 void U4::GenPhotonAncestor( const G4Track* aTrack )
0344 {
0345 #ifdef WITH_CUSTOM4_OLD
0346     ancestor = C4TrackInfo<C4Pho>::Get(aTrack) ;
0347 #elif WITH_CUSTOM4
0348     ancestor = C4TrackInfo::Get(aTrack) ;
0349 #else
0350     ancestor = STrackInfo::Get(aTrack) ;
0351 #endif
0352     if(dump) std::cout << "U4::GenPhotonAncestor " << ancestor.desc() << std::endl ;
0353     LOG(LEVEL) << ancestor.desc() ;
0354 }
0355
0356 /**
0357 U4::GenPhotonBegin
0358 -------------------
0359
0360 This is called from head of Scintillation and Cerenkov generation loops.
0361 It updates the private "pho" spho.h label which carries
0362 genstep index, photon index within genstep, photon identity, reemission index
0363
0364 **/
0365
0366 void U4::GenPhotonBegin( int genloop_idx )
0367 {
0368     assert(genloop_idx > -1);
0369     pho = gs.MakePho(genloop_idx, ancestor);
0370
0371     int align_id = ancestor.isPlaceholder() ? gs.offset + genloop_idx : ancestor.id ;
0372
0373     bool align_id_expect = pho.id == align_id ;
0374     assert( align_id_expect );
0375     if(!align_id_expect) std::raise(SIGINT);
0376
0377 #ifdef DEBUG
0378     if(dump) std::cout
0379         << "U4::GenPhotonBegin"
0380         << " genloop_idx " << std::setw(6) << genloop_idx
0381         << " gs.offset " << std::setw(6) << gs.offset
0382         << " pho.id " << std::setw(6) << pho.id
0383         << std::endl
0384         ;
0385 #endif
0386 }
0387
0388 /**
0389 U4::GenPhotonEnd
0390 ------------------
0391
0392 This is called from tail of Scintillation and Cerenkov generation loops
0393 following instanciation of the secondary track.
0394
0395 Sets spho label into the secondary track using U4PhotonInfo::Set
0396
0397 **/
0398
0399 void U4::GenPhotonEnd( int genloop_idx, G4Track* aSecondaryTrack )
0400 {
0401     assert(genloop_idx > -1);
0402     secondary = gs.MakePho(genloop_idx, ancestor) ;
0403     assert( secondary.isIdentical(pho) );  // make sure paired U4::GenPhotonBegin and U4::GenPhotonEnd
0404
0405 #ifdef DEBUG
0406     if(dump) std::cout << "U4::GenPhotonEnd " << secondary.desc() << std::endl ;
0407 #endif
0408
0409 #ifdef WITH_CUSTOM4_OLD
0410     C4TrackInfo<C4Pho>::Set(aSecondaryTrack, secondary );
0411 #elif WITH_CUSTOM4
0412     C4TrackInfo::Set(aSecondaryTrack, secondary );
0413 #else
0414     STrackInfo::Set(aSecondaryTrack, secondary );
0415 #endif
0416
0417 }
0418
0419
0420
0421 /**
0422 U4::CollectOpticalSecondaries
0423 ------------------------------
0424
0425 Optional and as yet incomplete conversion of G4VParticleChange into NP array of photons
0426
0427 **/
0428
0429 NP* U4::CollectOpticalSecondaries(const G4VParticleChange* pc )
0430 {
0431     G4int num = pc->GetNumberOfSecondaries();
0432
0433     std::cout << "U4::CollectOpticalSecondaries num " << num << std::endl ;
0434
0435     NP* p = NP::Make<float>(num, 4, 4);
0436     sphoton* pp = (sphoton*)p->bytes() ;
0437
0438     for(int i=0 ; i < num ; i++)
0439     {
0440         G4Track* track =  pc->GetSecondary(i) ;
0441         assert( track->GetParticleDefinition() == G4OpticalPhoton::Definition() );
0442         const G4DynamicParticle* ph = track->GetDynamicParticle() ;
0443         const G4ThreeVector& pmom = ph->GetMomentumDirection() ;
0444         const G4ThreeVector& ppol = ph->GetPolarization() ;
0445         sphoton& sp = pp[i] ;
0446
0447         // position ?
0448
0449         sp.mom.x = pmom.x();
0450         sp.mom.y = pmom.y();
0451         sp.mom.z = pmom.z();
0452
0453         sp.pol.x = ppol.x();
0454         sp.pol.y = ppol.y();
0455         sp.pol.z = ppol.z();
0456     }
0457     return p ;
0458 }
0459
0460
0461 void U4::GenPhotonSecondaries( const G4Track* , const G4VParticleChange* )
0462 {
0463     // do nothing
0464 }
0465
0466