Total internal reflection process not working as expected

Hi,

I have optical physics turned on in my simulation and am propagating Cherenkov light up cladded/buffered fibers, and I seem to get a significant amount of light escaping the fibers that should be captured.

Is there any reason to expect an optical photon moving from core to cladding with an angle of incidence greater than the critical angle not to be total internally reflected? Because I see individual photons sometimes entering a “captured TIR state” where it proceeds to bounce back-and-forth between core/cladding indefinitely and other times where it escapes and continues on to the cladding/buffer interface. The latter case happens more often.

Some details on the simulation: To get the angle of incidence relative to the surface normal of the fiber layers I’m using the method @vyeroshe gave here: https://geant4-forum.web.cern.ch/t/angle-of-incidence-for-optical-photon/461/5. That seems to work well, as all the captured TIR photons that reach the top of my fibers have an angle of incidence greater than the critical angle. But as I described above, there are other times where the original track won’t enter this captured TIR mode, even though its angle of incidence appears greater than the critical angle. Is there another optical physics process happening here that I’m missing or some other reason not to expect every photon satisfying the TIR-condition to be captured? For the core and cladding, the only material properties I define are the indices of refraction and absorption lengths. And there’s no optical surface or skin along this boundary.

I’m attaching a printout with an example of what I’m seeing. The first track (6911361) has an incidence angle greater than the critical angle (82.148 in this case), and in the first step shown, it travels from core to cladding with an angle greater than the critical angle. The second track (6965460) does a similar thing, except you can see it becomes captured in the next step.

Thanks in advance

Another observation: the TIR appears to begin in the cladding to core step. So I can identify a captured photon by asking whether the incidence angle repeats indefinitely, starting from an initial cladding to core step (outer to inner).

This is surprising to me. I expected to identify captured light based off a core prestep and cladding poststep. Is there something I’m missing (either physical or inherent to Geant) that makes this not true?

I’m trying to understand this because there’s a significant amount of light escaping my fibers that seems to satisfy the TIR condition.

Well, according to the linked discussion, @vyeroshe 's method didn’t work for them, either.

I don’t have much in the way of concrete suggestions. We’d have to step through the code to see what is going on. Maybe the starting point would be compare your calculation of the angle of incidence to that in Geant4 here:
https://geant4.kek.jp/lxr/source/processes/optical/src/G4OpBoundaryProcess.cc#L1095
and see if they agree. (You could start with adding a printout at that location.)

Two other thoughts:

  1. What is the boundary process status?
  2. I don’t know what this code does.
    Geant4/processes/optical/src/G4OpBoundaryProcess.cc
    What happens if you enable it?

@dsawkey thanks for the feedback.

Yeah, so I did a few things to clarify what was going on. I changed the verboseLevel to 2 in G4OpBoundaryProcess and seeing theStatus from BoundaryProcessVerbose() was especially helpful. I also changed G4OpBoundaryProcess::GetIncidentAngle() from private to public (why isn’t this public by default!?).

What I found was that there are a lot of subtleties the user has to be aware of. First of all, @vyeroshe 's method does work well and almost always matches Geant’s internal calculation from GetIncidentAngle(). In my cladded/buffered fiber application, the primary exception was at fiber ends or other places where the direction of the normal changed. It also obviously delivers a different answer anywhere where microfacets are defined (as the calculation assumes polished surfaces). The other detail to take care of is where precisely you are calling each function from within the stepping process and which step point you’re probing. The presteppoint and poststeppoint predictably give complementary angles during total internal reflection, but I’ve also been using a FastSimulationImplementation that probes the track before the step is made. At this point, the poststep angle is undefined and theStatus often returns a stepTooShort, since the step size is still 0 at that point in time. If you compare theStatus and incident angle for that same step to something in a sensitive detector class (which are invoked after the step), you’ll see differences. You may also see differences in the first optical step. Depending when you’re calling GetIncidentAngle() it may still show the angle from the previous track’s final step. So in my case, when calling it from the FastSimulation class, the direct calculation gave better results for that first step.

Hopefully this helps others in the future. And to developers, it really would be helpful to make G4OpBoundaryProcess::GetIncidentAngle() public in future versions. It’s a valuable cross check and is useful for various detector R&D applications.

thanks for the awesome information.

It’s not clear to me how, or what, to implement. GetIncidentAngle() returns the angle of the particle track relative to a facet. In one boundary interaction, there can be multiple reflections off different facets, if the surface roughness is high enough. I suppose GetIncidentAngle() is only useful to the user if the surface roughness is zero.

Hello,

I’m working on the same project as the OP, though I’d like to simplify what we’re asking for. This might be better in a new thread, but I’ll check here to start.

Rather than getting the low level information and calculating total internal reflection criteria within the SD, I’d just like to perform a check to see if the TIR process was called. Unfortunately, this is presenting a new problem in that the only process I can detect within the SD is optical absorption using the following code:

if( aStep->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition() ){
  if( aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessType() == fOptical ){
    G4cout << aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessSubType() << G4endl;
  }
}

Though if I set /process/optical/verbose 2 I can see that the process sub type is frequently fOpBoundary and total internal reflection does occur. Ideally I’d like my code to look like:

if( aStep->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition() ){
  if( aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessType() == fOptical ){
    if( aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessSubType() == fOpBoundary){
      G4OpBoundaryProcess *opProc = (G4OpBoundaryProcess*)aStep->GetPostStepPoint()->GetProcess()
      if( opProc->GetStatus() == G4OpBoundaryProcessStatus::TotalInternalReflection){
        //Do stuff
      }
    }
  } 
}

Thank you in advance

Edit:

After looking a little closer at the output with /process/optical/verbose 2, I noticed that total internal reflection gets tagged as absorption. This seems like a bug

G4VDiscreteProcess::PostStepGetPhysicalInteractionLength() - [ OpAbsorption]
 Particle type - opticalphoton
   mass:        0[GeV]
   charge:      0[e]
   Direction x: -0.65093, y: 0.0306664, z: 0.758518
   Total Momentum = 3.32434e-09[GeV]
   Momentum: -2.16392e-09[GeV], y: 1.01945e-10[GeV], z: 2.52157e-09[GeV]
   Total Energy   = 3.32434e-09[GeV]
   Kinetic Energy = 3.32434e-09[GeV]
 MagneticMoment  [MeV/T]: 0
   ProperTime     = 0[ns]
 in Material  Quartz
InteractionLength= 1100.53[cm]
 Photon at Boundary!
 thePrePV:  ZDC1_Core215
 thePostPV: ZDC1_Air_Physical
 Old Momentum Direction: (-0.65093,0.0306664,0.758518)
 Old Polarization:       (0.75763,0.0891851,0.646563)
 New Momentum Direction: (-0.992035,0.0306664,-0.122175)
 New Polarization:       (-0.12446,-0.0891851,0.988208)
 *** TotalInternalReflection ***

Indeed, the sensitive detector only detects absorbed photons. If you put your “do stuff” code into UserSteppingAction it will do what you want. Add a check to make sure the volume is correct.

For the edit: The output for OpAbsorption ends at “interaction length = 1100 cm”. The rest is from OpBoundary.

Thank you for the clarification.

I moved my code into UserSteppingAction as you suggested, but I’m still only seeing the absorption sub process.

Do you mean all the boundary interactions are absorption? If so, what is the material property REFLECTION set to?

Or is your code giving you the OpAbsorption process happening in the bulk material?

I’m not looking specifically for boundary interactions in my code, only that the track is within my sensitive volume (glass fiber). However, for all steps within that volume GetProcessSubType() returns fOpAbsorption. I suspect that if a boundary process were to occur GetProcessSubType() should return fOpBoundary.

My materials don’t have reflection set. I’m currently looking at simplified fibers, so the boundary is always glass->air. I didn’t think setting a value for reflection on those materials would be necessary since reflectivity at that boundary is governed by refractive index and incident angle.

If you’re not seeing the boundary process at all, something is wrong with your code! Compare what you do to OpNovice2 example.

The problem is that I’m seeing the boundary process when I set /process/optical/verbose 2. So I know it’s happening, but I can’t identify when it’s happening from within my code. I’m also using the same seed every time, so I’m 100% sure total internal reflection occurs even though my “do stuff” code doesn’t identify it. So the only issues with my code that I can think of are the checks before “do stuff”. Am I performing those checks correctly, or am I making some assumption that’s not accurate?

Update:

@dsawkey Thank you for your help. I’m leaving a brief summary in case it may help someone in the future.

I compared with OpNovice2 as you suggested and found the inaccurate assumption I was making. I assumed there would only be a single optical process associated with the post step point and I was using G4StepPoint::GetProcessDefinedStep() to find it. I thought this was the case because I was started at the doxygen documentation for G4Step and looked for what information was available from that class. It seems that a single step can have multiple optical processes associated with it and they are accessed via G4ProcessManager::GetPostStepProcessVector().

The code below is what I’ll be using. Taken almost verbatim from OpNovice2 and trimmed down for my purpose. It is usable in the SD, not just UserSteppingAction.

      static G4ParticleDefinition* opticalphoton =
      G4OpticalPhoton::OpticalPhotonDefinition();
      if( aStep->GetTrack()->GetDefinition() == opticalphoton){
        if( aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary){
          G4ProcessVector* postStepDoItVector =
            opticalphoton->GetProcessManager()->GetPostStepProcessVector(typeDoIt);
          G4int n_proc = postStepDoItVector->entries();
          for(G4int i = 0; i < n_proc; ++i){
            G4VProcess* currentProcess = (*postStepDoItVector)[i];
            G4OpBoundaryProcess* opProc = dynamic_cast<G4OpBoundaryProcess*>(currentProcess);
            if(opProc){
              if(opProc->GetStatus() == TotalInternalReflection){
                //Do stuff
              }//end if TotalInternalReflection
            }//end if opProc
          }//end for n_proc
        }//end if fGeomBoundary
      }//end if optical photon