Added absorption

MonteCarlo/src/FiberSD.cc

@@ -160,30 +160,44 @@ G4bool FiberSD::ProcessHits(G4Step* aStep,G4TouchableHistory*){
       //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 ){
+
         //Now check for total internal reflection
+        //There are potentially multiple processes for our post step point, so we get the vector
+        //and step through it looking for TIR
         G4ProcessVector* postStepDoItVector = opticalphoton->GetProcessManager()->GetPostStepProcessVector(typeDoIt);
         G4int n_proc = postStepDoItVector->entries();
-
         for(G4int i = 0; i < n_proc; ++i){
           G4OpBoundaryProcess* opProc = dynamic_cast<G4OpBoundaryProcess*>((*postStepDoItVector)[i]);
           if(opProc && opProc->GetStatus() == TotalInternalReflection){
-            if(REDUCED_TREE){
-              m_cherenkovVec->at(rodNum)++;
-              FillTimeVector( rodNum, aStep->GetTrack()->GetGlobalTime() );
-
-              m_nHits++;
-              return true;
-            }//end if REDUCED_TREE
-            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->setRodNb    ( rodNum );
-
-            fiberCollection->insert ( newHit );
-            m_nHits++;
+            //Here we're simulating absorption. We determine the full path length of the photon within the fiber
+            //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 absChance = exp(-pathLength/Absorption);
+            if(CLHEP::RandFlat::shoot(0.0,1.0) < absChance){//This check amounts to if(!absorbed)
+              //Finally, record the ouput
+              if(REDUCED_TREE){
+                m_cherenkovVec->at(rodNum)++;
+                FillTimeVector( rodNum, aStep->GetTrack()->GetGlobalTime() );
+
+                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->setRodNb    ( rodNum );
+
+                fiberCollection->insert ( newHit );
+                m_nHits++;
+              }//end if REDUCED_TREE
+            }//end if !absorbed
           }//end if opProc && TotalInternalReflection
         }//end for n_proc
       }