PhotonFullFibreModel.cxx

Go to the documentation of this file.

#include "PhotonFullFibreModel.h"
#include "PhotonTransportMaker.h"
#include "PhotonUtil.h"
#include "MessageService/MsgService.h"
#include "Conventions/Munits.h"
#include <cmath>

CVSID( "$Id: PhotonFullFibreModel.cxx,v 1.10 2006/12/01 19:59:25 rhatcher Exp $" );

ClassImp(PhotonFullFibreModel)

PhotonTransportMakerProxy<PhotonFullFibreModel>
        gPhotonFullFibreModelProxy("photonFullFibreModel"); 

// Constructor
PhotonFullFibreModel::PhotonFullFibreModel() :
  fFibreAbsorbTable("PHOTONFIBREABSORB"),
  fGreenSpectrumTable("PHOTONGREENSPECTRUM")
{
  Configure(PhotonConfiguration());
    
  if (fFibreAbsorbTable.GetNumRows()<1) {
    MSG("Photon",Msg::kError) << "Can't find data for table PHOTONFIBREABSORB. Aborting!\n";
    assert(0);
  }
  if (fGreenSpectrumTable.GetNumRows()<1) {
    MSG("Photon",Msg::kError) << "Can't find data for table PHOTONFIBREABSORB. Aborting!\n";
    assert(0);
  }
}

// Configure the model.
void 
PhotonFullFibreModel::Configure( const Registry& r )
{
  MSG("Photon",Msg::kDebug) << "Configuring FullFibreModel." << std::endl;
  double tmpd;
  if (r.Get("FibreRadius",tmpd))          fFibreRadius  = tmpd;
  if (r.Get("FibreCoreRadius",tmpd))      fFibreCoreRadius  = tmpd;
  if (r.Get("InnerCladdingRadius",tmpd))  fInnerCladdingRadius  = tmpd;

  if (r.Get("FibreFlurorDecayTime",tmpd)) fFibreFlurorDecayTime  = tmpd;

  if (r.Get("FibreIndex",tmpd))           fFibreIndex  = tmpd;
  if (r.Get("ScintIndex",tmpd))           fScintIndex  = tmpd;
  if (r.Get("OuterCladdingIndex",tmpd))   fOuterCladdingIndex  = tmpd;
  if (r.Get("InnerCladdingIndex",tmpd))   fInnerCladdingIndex  = tmpd;

  SetWlsPrescaleFactor(1.0); // No prescaling.
}

void 
PhotonFullFibreModel::Reset( const VldContext& cx )
{
  // Ensure tables are up-to-date.
  fFibreAbsorbTable.Reset(cx);
  fGreenSpectrumTable.Reset(cx);
}

Int_t
PhotonFullFibreModel:: FindBlueAbsorptionPosition( const UgliStripHandle& inStrip,
                                                   DigiPhoton*            inBluePhoton,
                                                   TVector3&              outPos)
{
  // Track the blue photon until it gets absorbed.
  // Return -ve number for a condition that shouldn't happen 
  // Return 0 if photon stopped
  // Return +ve number if photon escaped

  // Reset the blue photon coordinates so they're centered on our fibre.
  TVector3 x = inBluePhoton->GetPos();
  
  x -= TVector3(0,0,
                inStrip.GetHalfThickness() - fFibreRadius);
  fBlue_Start = fBlue_Stop = fBlue_OuterHit = fBlue_InnerHit = x;
  
  TVector3 v = inBluePhoton->GetDir();
  
  // Interact with the outer cladding.
  int interact  = PhotonUtil::InteractWithCylinder( x,v,
                                                    fFibreRadius,
                                                    fScintIndex,
                                                    fOuterCladdingIndex);
  if(interact!=0) { 
    // It reflected, not refracted. How did this happen?
    fBlue_OuterHit = x + v*0.0001;
    MSG("Photon",Msg::kWarning) << "Didn't refract into fibre. Huh? x=(" 
                                << x.x() << " " << x.y() << " " << x.z() << ") "
                                << "v=(" << v.x() << " " << v.y() << " "  << v.z() <<")\n";
                                
    return -2;
  } 
  // Interact with the inner cladding
  int inclad1 = PhotonUtil::InteractWithCylinder( x,v,
                                                   fInnerCladdingRadius,
                                                   fOuterCladdingIndex,                 
                                                   fInnerCladdingIndex );

  fBlue_OuterHit = x;
  if(inclad1!=0) { return 1; } // It got out.

  // Intersect with the core:
  int incore = PhotonUtil::InteractWithCylinder( x, v,
                                                 fFibreCoreRadius,
                                                 fInnerCladdingIndex,
                                                 fFibreIndex );
  fBlue_InnerHit = x;
  // Missed the core or reflected off it. Shucks.
  if(incore!=0) { return 2; };
  
  // Photon hit the core. See how far it could go.
  double d1, d2;
  incore = PhotonUtil::IntersectCylinder( x, v, fFibreCoreRadius, d1, d2 );
  
  // This shouldn't happen.. we've already said it hit the core.
  if(!incore) { return -3;}

  // Find wavelength of blue photon in nm.
  double wavelength = inBluePhoton->WaveLength()/Munits::nm;
  if(wavelength<=0) wavelength = 420; // Something reasonable.
  
  // Find absorption length of this photon in nm.
  double atten = fFibreAbsorbTable.Interpolate(wavelength);
  
  // find where the photon got absorbed.
  double d = d1 + fRandom->Exp(atten);

  // Compute the final postion.
  outPos = x + v*d; 
  fBlue_Stop = outPos;;
  if(d>d2) { return 3; } // the blue photon managed to escape through the other side.
  
  return 0;
}

DigiPhoton*
PhotonFullFibreModel::MakeGreenPhoton(   DigiPhoton*            inBluePhoton,
                                         const TVector3&        inPos,
                                         StripEnd::StripEnd_t   inDirection
                                         )
{
  // Choose a direction.
  double cosTheta  = 2.*fRandom->Rndm()-1.0;
  double sinTheta = TMath::Sin(TMath::ACos(cosTheta));
  double phi = 2.*TMath::Pi()*fRandom->Rndm();
  TVector3 dir( cosTheta, 
                sinTheta*TMath::Cos(phi),
                sinTheta*TMath::Sin(phi) );
  
  // Let's ensure it's going in the correct direction:
  switch(inDirection) {
  case StripEnd::kPositive:
    dir.SetX(fabs(dir.X()));
    break;
    
  case StripEnd::kNegative:
    dir.SetX(-fabs(dir.X()));
    break;
    
  default:
    // Do nothing.
    break;    
  }

  // Choose wavelength.
  double wavelength = fGreenSpectrumTable.Pick(fRandom->Rndm())*Munits::nm;
  
  double time = inBluePhoton->T() + fRandom->Exp(fFibreFlurorDecayTime);

  DigiPhoton* green =  new DigiPhoton( inBluePhoton->ParentHit(),
                                       wavelength,
                                       time,
                                       inPos,
                                       dir
                                       );
  return green;
}

Bool_t 
PhotonFullFibreModel::IsGreenTrapped( DigiPhoton* inGreen )
{
  TVector3 x   = inGreen->GetPos();
  TVector3 dir = inGreen->GetDir();

  // Special cases:
  if(dir.x()==0) return false; // Going sideways.
  if(fabs(dir.x())>=1.) return true; // Going straight down the pipe.

  fGreen_Start = x;

  int interact;
  // Interact with the inner/outer cladding interface.
  interact = PhotonUtil::InteractWithCylinder( x,dir,
                                               fFibreCoreRadius,
                                               fFibreIndex,
                                               fInnerCladdingIndex);
  if(interact==-1) {
    MSG("Photon",Msg::kError) << "Didn't hit inner wall IsGreenTrapped(). Huh?!?" << endl;
  };
  fGreen_InnerHit = x;
  if(interact==1) {
    // It reflected off the inner surface. Bonus.
    fGreen_OuterHit = x;
    return true;
  }


  // Interact with the outer cladding/scint interface. 
  interact = PhotonUtil::InteractWithCylinder( x,dir,
                                               fInnerCladdingRadius,
                                               fInnerCladdingIndex,
                                               fOuterCladdingIndex);  
  if(interact==-1) {
    MSG("Photon",Msg::kError) << "Didn't hit outer wall IsGreenTrapped(). Huh?!?" << endl;
  };
  fGreen_OuterHit = x;

  if(interact==0) return false; // Refracted out.
  return true; // Reflected in.
}

Bool_t 
PhotonFullFibreModel::FibreHitToGreenPhoton( const UgliStripHandle& inStrip,
                                             StripEnd::StripEnd_t   inDirection,
                                             DigiPhoton*            inBluePhoton,
                                             DigiPhoton*            &outGreenPhoton
                                             )
{
  // Return a contained green photon.
  // This contains the full photon simulation.
  // Blue photons are tracked through the cladding and into the core
  // where they are absorbed by the correct attenuation.
  fLastResult = -1;

  outGreenPhoton = NULL;
  if(!inBluePhoton) { return false; };

  TVector3 flashPos;
  fLastResult = FindBlueAbsorptionPosition( inStrip,
                                            inBluePhoton,
                                            flashPos );

  if(fLastResult !=0 ) return false;

  // Loop until we have a valid solution.

  // Find the new direction of the green photon.
  // Go in the forward or backward direction, depending on the 
  // given parameter.
  
  if(outGreenPhoton) delete outGreenPhoton;
  outGreenPhoton = MakeGreenPhoton( inBluePhoton, 
                                    flashPos,
                                    inDirection);

  if(!outGreenPhoton) return false; // Huh?

  if(IsGreenTrapped(outGreenPhoton)) return true;
  
  // else:
  return false; // No trapped photon!
}

---