Commit 1562a2dd authored by Chad Lantz's avatar Chad Lantz
Browse files

Corrected path length calculation and photon transit time

parent 15178ec6
......@@ -158,8 +158,9 @@ G4bool FiberSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
Rindex = RindexMPV->Value(aStep->GetTrack()->GetDynamicParticle()->GetTotalMomentum(), index);
//Use the refractive index to see if the photon will be transmitted at the top of the fiber
//instead of reflected back downward
// if( 3.1415/2.0 - atan( momentum.y() / ( 1-momentum.y() ) ) >= asin(1.0/Rindex) ){
if( momentum.y() > 0.0 ){
G4ThreeVector momDir = aStep->GetPreStepPoint()->GetMomentumDirection();
G4double angleFromY = atan( sqrt( 1 - pow(momDir.y(),2.0) )/momDir.y() );
if( angleFromY >= asin(1.0/Rindex) ){//if angle is greater than TIR critical angle
//Now check for total internal reflection
//There are potentially multiple processes for our post step point, so we get the vector
......@@ -173,25 +174,28 @@ G4bool FiberSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
//and roll a random number against it's chance of absorption. If it's not absorbed record it
G4MaterialPropertyVector* AbsMPV = MPT->GetProperty(kABSLENGTH);
G4double Absorption = AbsMPV->Value(aStep->GetTrack()->GetDynamicParticle()->GetTotalMomentum(), index);
G4ThreeVector momDir = aStep->GetPreStepPoint()->GetMomentumDirection();
//pathLength = (distance to top of fiber)/sin(polarAngle); polarAngle = atan(y/sqrt(x^2 + z^2));
G4double pathLength = (m_topOfVolume - pos.y())/sin( atan(momDir.y()/sqrt( pow(momDir.x(),2.0) + pow(momDir.z(),2.0) ) ) );
G4double pathLength = (m_topOfVolume - pos.y())/cos( angleFromY );
G4double absChance = 1 - exp(-pathLength/Absorption);
if(CLHEP::RandFlat::shoot(0.0,1.0) > absChance){//This check amounts to if(!absorbed)
//This check amounts to if(!absorbed)
if(CLHEP::RandFlat::shoot(0.0,1.0) > absChance){
//Calculate the time the photon would arrive at the top of the volume
G4double time = aStep->GetTrack()->GetGlobalTime() + Rindex*pathLength/CLHEP::c_light;
//Finally, record the ouput
if(REDUCED_TREE){
m_cherenkovVec->at(rodNum)++;
FillTimeVector( rodNum, aStep->GetTrack()->GetGlobalTime() );
FillTimeVector( rodNum, time );
m_nHits++;
return true;
}else{
FiberHit* newHit = new FiberHit();
newHit->setPos ( pos );
newHit->setOrigin ( aStep->GetTrack()->GetVertexPosition() );
newHit->setMomentum ( momentum );
newHit->setEnergy ( aStep->GetPreStepPoint()->GetTotalEnergy() );
newHit->setTime ( aStep->GetTrack()->GetGlobalTime() );
newHit->setTime ( time );
newHit->setRodNb ( rodNum );
fiberCollection->insert ( newHit );
......@@ -200,7 +204,7 @@ G4bool FiberSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
}//end if !absorbed
}//end if opProc && TotalInternalReflection
}//end for n_proc
}
}//end if > critical angle
//Kill all photons
aStep->GetTrack()->SetTrackStatus( fStopAndKill );
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment