eic-opticks/u4/U4StepPoint.cc

0001 /**
0002 U4StepPoint.cc
0003 ===============
0004
0005 Boundary class changes need to match in all the below::
0006
0007     U4OpBoundaryProcess.h
0008     U4Physics.cc
0009     U4Recorder.cc
0010     U4StepPoint.cc
0011
0012 **/
0013
0014
0015 #include <sstream>
0016
0017 #include "G4StepPoint.hh"
0018 #include "G4VProcess.hh"
0019 #include "G4ThreeVector.hh"
0020 #include "G4SystemOfUnits.hh"
0021 #include "G4PhysicalConstants.hh"
0022
0023 #include "SLOG.hh"
0024 #include "SSys.hh"
0025 #include "OpticksPhoton.h"
0026 #include "OpticksPhoton.hh"
0027 #include "scuda.h"
0028 #include "squad.h"
0029 #include "sphoton.h"
0030 #include "SFastSimOpticalModel.hh"
0031
0032 #include "U4StepStatus.h"
0033 #include "U4OpBoundaryProcess.h"
0034 #include "U4OpBoundaryProcessStatus.h"
0035 #include "U4Physics.hh"
0036 #include "U4StepPoint.hh"
0037
0038
0039 const plog::Severity U4StepPoint::LEVEL = SLOG::EnvLevel("U4StepPoint", "DEBUG");
0040 const char* U4StepPoint::OpFastSim_ = SSys::getenvvar("U4StepPoint_OpFastSim", "fast_sim_man" );
0041
0042
0043 /**
0044 U4StepPoint::Update
0045 ---------------------
0046
0047 * cf CWriter::writeStepPoint_
0048
0049 **/
0050
0051 void U4StepPoint::Update(sphoton& photon, const G4StepPoint* point)  // static
0052 {
0053     const G4ThreeVector& pos = point->GetPosition();
0054     const G4ThreeVector& mom = point->GetMomentumDirection();
0055     const G4ThreeVector& pol = point->GetPolarization();
0056
0057     G4double time = point->GetGlobalTime();
0058     G4double energy = point->GetKineticEnergy();
0059     G4double wavelength = h_Planck*c_light/energy ;
0060
0061     photon.pos.x = pos.x();
0062     photon.pos.y = pos.y();
0063     photon.pos.z = pos.z();
0064     photon.time  = time/ns ;
0065
0066     photon.mom.x = mom.x();
0067     photon.mom.y = mom.y();
0068     photon.mom.z = mom.z();
0069     //photon.iindex = 0u ;
0070
0071     photon.pol.x = pol.x();
0072     photon.pol.y = pol.y();
0073     photon.pol.z = pol.z();
0074     photon.wavelength = wavelength/nm ;
0075 }
0076
0077 std::string U4StepPoint::DescPositionTime(const G4StepPoint* point )
0078 {
0079     const G4ThreeVector& pos = point->GetPosition();
0080     G4double time = point->GetGlobalTime();
0081     std::stringstream ss ;
0082
0083     ss << "U4StepPoint::DescPositionTime ("
0084        << " " << std::setw(10) << std::fixed << std::setw(10) << std::setprecision(3) << pos.x()
0085        << " " << std::setw(10) << std::fixed << std::setw(10) << std::setprecision(3) << pos.y()
0086        << " " << std::setw(10) << std::fixed << std::setw(10) << std::setprecision(3) << pos.z()
0087        << " " << std::setw(10) << std::fixed << std::setw(10) << std::setprecision(3) << time/ns
0088        << ")"
0089        ;
0090
0091     std::string s = ss.str();
0092     return s ;
0093 }
0094
0095 const char* U4StepPoint::ProcessName(const G4StepPoint* point) // static
0096 {
0097     const G4VProcess* process = point->GetProcessDefinedStep() ;
0098     if(process == nullptr) return nullptr ;
0099     const G4String& processName = process->GetProcessName() ;
0100     return processName.c_str();
0101 }
0102
0103 unsigned U4StepPoint::ProcessDefinedStepType(const G4StepPoint* point) // static
0104 {
0105     const G4VProcess* process = point->GetProcessDefinedStep() ;
0106     if(process == nullptr) return U4StepPoint_NoProc ;
0107     const G4String& processName = process->GetProcessName() ;
0108     return ProcessDefinedStepType( processName.c_str() );
0109 }
0110
0111 unsigned U4StepPoint::ProcessDefinedStepType(const char* name) // static
0112 {
0113     unsigned type = U4StepPoint_Undefined ;
0114     if(strcmp(name, NoProc_) == 0 )         type = U4StepPoint_NoProc ;
0115     if(strcmp(name, Transportation_) == 0 ) type = U4StepPoint_Transportation ;
0116     if(strcmp(name, OpRayleigh_) == 0)      type = U4StepPoint_OpRayleigh ;
0117     if(strcmp(name, OpAbsorption_) == 0)    type = U4StepPoint_OpAbsorption ;
0118     if(strcmp(name, OpFastSim_) == 0)       type = U4StepPoint_OpFastSim ;
0119     if(strcmp(name, OTHER_) == 0)           type = U4StepPoint_OTHER ;
0120     return type ;
0121 }
0122
0123 const char* U4StepPoint::ProcessDefinedStepTypeName( unsigned type ) // static
0124 {
0125     const char* s = nullptr ;
0126     switch(type)
0127     {
0128         case U4StepPoint_Undefined:      s = Undefined_      ; break ;
0129         case U4StepPoint_NoProc:         s = NoProc_         ; break ;
0130         case U4StepPoint_Transportation: s = Transportation_ ; break ;
0131         case U4StepPoint_OpRayleigh:     s = OpRayleigh_     ; break ;
0132         case U4StepPoint_OpAbsorption:   s = OpAbsorption_   ; break ;
0133         case U4StepPoint_OpFastSim:      s = OpFastSim_      ; break ;
0134         default:                         s = OTHER_          ; break ;
0135     }
0136     return s ;
0137 }
0138
0139
0140 template <typename T>
0141 bool U4StepPoint::IsTransportationBoundary(const G4StepPoint* point)
0142 {
0143     G4StepStatus status = point->GetStepStatus()  ;
0144     unsigned proc = ProcessDefinedStepType(point);
0145     return status == fGeomBoundary && proc == U4StepPoint_Transportation ;
0146 }
0147
0148
0149
0150 /**
0151 U4StepPoint::Flag
0152 ------------------
0153
0154 Adapted from cfg4/OpStatus.cc:OpStatus::OpPointFlag
0155
0156 Q: Why does this never return BULK_REEMIT ?
0157 A: As that always starts as BULK_ABSORB and then gets changed into BULK_REEMIT
0158    by subsequent history rewriting when another track comes along with
0159    the corresponding ancestry.
0160
0161 **/
0162
0163 template <typename T>
0164 unsigned U4StepPoint::Flag(const G4StepPoint* point, bool warn, bool& tir )
0165 {
0166     G4StepStatus status = point->GetStepStatus()  ;
0167     unsigned proc = ProcessDefinedStepType(point);
0168     unsigned flag = 0 ;
0169     tir = false ;
0170
0171     if( status == fPostStepDoItProc && proc == U4StepPoint_OpAbsorption )
0172     {
0173         flag = BULK_ABSORB ;
0174     }
0175     else if( status == fPostStepDoItProc && proc == U4StepPoint_OpRayleigh )
0176     {
0177         flag = BULK_SCATTER ;
0178     }
0179     else if( status == fGeomBoundary && proc == U4StepPoint_Transportation )
0180     {
0181         unsigned bstat = U4OpBoundaryProcess::GetStatus<T>();
0182         T* proc = U4OpBoundaryProcess::Get<T>();
0183
0184         LOG_IF( error, bstat == Undefined )
0185             << " U4OpBoundaryProcess::GetStatus<T>() : Undefined "
0186             << std::endl
0187             << " U4OpBoundaryProcess::Get<T>() " << ( proc ? "YES" : "NO " )
0188             << std::endl
0189             << " U4Physics::Switches() "
0190             << std::endl
0191             << U4Physics::Switches()
0192             ;
0193
0194         tir = bstat == TotalInternalReflection ;
0195
0196         flag = BoundaryFlag(bstat) ;   // BT BR NA SA SD SR DR
0197
0198         LOG_IF(error, flag == 0 )   // NAN_ABORT are common from StepTooSmall
0199             << " UNEXPECTED BoundaryFlag ZERO  "
0200             << std::endl
0201             << " flag " << flag
0202             << " OpticksPhoton::Flag(flag) " << OpticksPhoton::Flag(flag)
0203             << std::endl
0204             << " bstat " << bstat
0205             << " U4OpBoundaryProcessStatus::Name(bstat) " << U4OpBoundaryProcessStatus::Name(bstat)
0206             ;
0207     }
0208     else if( status == fGeomBoundary && proc == U4StepPoint_OpFastSim )
0209     {
0210         flag = DEFER_FSTRACKINFO ; // signal that must get FastSim DoIt status from trackinfo label
0211     }
0212     else if( status == fWorldBoundary && proc == U4StepPoint_Transportation )
0213     {
0214         flag = MISS ;
0215     }
0216     else
0217     {
0218         if(warn)
0219         {
0220             const char* procName = ProcessName(point);
0221             LOG(error)
0222                 << " failed to define flag for StepPoint "
0223                 << " G4StepStatus " << U4StepStatus::Name(status)
0224                 << " proc " << ProcessDefinedStepTypeName(proc )
0225                 << " procName " << ( procName ? procName : "-" )
0226                 ;
0227         }
0228
0229     }
0230     return flag ;
0231 }
0232
0233 /**
0234 U4StepPoint::BoundaryFlag
0235 --------------------------
0236
0237 Canonical stack::
0238
0239     U4StepPoint::BoundaryFlag
0240     U4StepPoint::Flag
0241     U4Recorder::UserSteppingAction_Optical
0242     U4Recorder::UserSteppingAction
0243
0244 **/
0245
0246 unsigned U4StepPoint::BoundaryFlag(unsigned status) // BT BR NA SA SD SR DR
0247 {
0248     unsigned flag = 0 ;
0249     switch(status)
0250     {
0251         case FresnelRefraction:
0252         case SameMaterial:
0253         case Transmission:
0254                                flag=BOUNDARY_TRANSMIT;
0255                                break;
0256         case TotalInternalReflection:
0257         case       FresnelReflection:
0258                                flag=BOUNDARY_REFLECT;
0259                                break;
0260         case StepTooSmall:
0261                                flag=NAN_ABORT;
0262                                break;
0263         case Absorption:
0264                                flag=SURFACE_ABSORB ;
0265                                break;
0266         case Detection:
0267                                flag=SURFACE_DETECT ;
0268                                break;
0269         case SpikeReflection:
0270                                flag=SURFACE_SREFLECT ;
0271                                break;
0272         case LobeReflection:
0273         case LambertianReflection:
0274                                flag=SURFACE_DREFLECT ;
0275                                break;
0276         case NoRINDEX:
0277                                flag=SURFACE_ABSORB ;
0278                                //flag=NAN_ABORT;
0279                                break;
0280         case Undefined:
0281         case BackScattering:
0282         case NotAtBoundary:
0283         case PolishedLumirrorAirReflection:
0284         case PolishedLumirrorGlueReflection:
0285         case PolishedAirReflection:
0286         case PolishedTeflonAirReflection:
0287         case PolishedTiOAirReflection:
0288         case PolishedTyvekAirReflection:
0289         case PolishedVM2000AirReflection:
0290         case PolishedVM2000GlueReflection:
0291         case EtchedLumirrorAirReflection:
0292         case EtchedLumirrorGlueReflection:
0293         case EtchedAirReflection:
0294         case EtchedTeflonAirReflection:
0295         case EtchedTiOAirReflection:
0296         case EtchedTyvekAirReflection:
0297         case EtchedVM2000AirReflection:
0298         case EtchedVM2000GlueReflection:
0299         case GroundLumirrorAirReflection:
0300         case GroundLumirrorGlueReflection:
0301         case GroundAirReflection:
0302         case GroundTeflonAirReflection:
0303         case GroundTiOAirReflection:
0304         case GroundTyvekAirReflection:
0305         case GroundVM2000AirReflection:
0306         case GroundVM2000GlueReflection:
0307         case Dichroic:
0308                                flag=0;   // leads to bad flag asserts
0309                                break;
0310     }
0311     return flag ;
0312 }
0313
0314
0315 template <typename T>
0316 std::string U4StepPoint::Desc(const G4StepPoint* point)
0317 {
0318     G4StepStatus status = point->GetStepStatus()  ;
0319     const char* statusName = U4StepStatus::Name(status);
0320
0321     unsigned proc = ProcessDefinedStepType(point);
0322     const char* procName = ProcessDefinedStepTypeName(proc);
0323     const char* procNameRaw = ProcessName(point);
0324
0325     unsigned bstat = U4OpBoundaryProcess::GetStatus<T>();
0326     const char* bstatName = U4OpBoundaryProcessStatus::Name(bstat);
0327
0328     bool warn = false ;
0329     bool is_tir = false ;
0330     unsigned flag = Flag<T>(point, warn, is_tir);
0331     const char* flagName = OpticksPhoton::Flag(flag);
0332
0333     std::stringstream ss ;
0334     ss << "U4StepPoint::Desc"
0335        << std::endl
0336        << " proc " << proc
0337        << " procName " << procName
0338        << " procNameRaw " << ( procNameRaw ? procNameRaw : "-" )
0339        << std::endl
0340        << " status " << status
0341        << " statusName " << statusName
0342        << std::endl
0343        << " bstat " << bstat
0344        << " bstatName " << bstatName
0345        << " is_tir " << is_tir
0346        << std::endl
0347        << " flag " << flag
0348        << " flagName " << flagName
0349        ;
0350     std::string s = ss.str();
0351     return s ;
0352 }
0353
0354 #if defined(WITH_CUSTOM4)
0355 template unsigned U4StepPoint::Flag<C4OpBoundaryProcess>(const G4StepPoint*, bool, bool& );
0356 template std::string U4StepPoint::Desc<C4OpBoundaryProcess>(const G4StepPoint* point);
0357 #elif defined(WITH_PMTSIM)
0358 template unsigned U4StepPoint::Flag<CustomG4OpBoundaryProcess>(const G4StepPoint*, bool, bool& );
0359 template std::string U4StepPoint::Desc<CustomG4OpBoundaryProcess>(const G4StepPoint* point);
0360 #elif defined(WITH_INSTRUMENTED_DEBUG)
0361 template unsigned U4StepPoint::Flag<InstrumentedG4OpBoundaryProcess>(const G4StepPoint*, bool, bool& );
0362 template std::string U4StepPoint::Desc<InstrumentedG4OpBoundaryProcess>(const G4StepPoint* point);
0363 #else
0364 template unsigned U4StepPoint::Flag<G4OpBoundaryProcess>(const G4StepPoint*, bool, bool& );
0365 template std::string U4StepPoint::Desc<G4OpBoundaryProcess>(const G4StepPoint* point);
0366 #endif
0367
0368