ScintPhoton.cxx

Go to the documentation of this file.

// A class to track photons in the scintillator of the
//MINOS far detector.  
//
#include <iostream>
#include <TMath.h>
#include <TRandom.h>

#include "ScintPhoton.h"

#include "Conventions/Munits.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "MessageService/MsgService.h"

CVSID( " $Id: ScintPhoton.cxx,v 1.8 2006/12/01 20:12:11 rhatcher Exp $" );

  // define the scintillator surfaces as follows
  //
  //                 I
  //        -------------------          ^ z
  //       |         O         |         |
  //    IV |          \        | II      |
  //       |           fiber   |         ---> y  
  //        -------------------     x out of page.
  //                III
  //
  // Surface V is positive end (towards us)
  // Surface VI is negative end (away from us)
  // Surface VII is the fibre cladding.
  //
  //coordinate system: (0,0,0) at center of strip
  //                    x-axis along long axis
  //                    z-axis up
  //                    y-axis to make right handed system

ClassImp(ScintPhoton)
 
// Static members.
Double_t ScintPhoton::fsOuterCladdingIndex = 1.42;  //index of fraction of outer cladding
Double_t ScintPhoton::fsInnerCladdingIndex = 1.49;  //index of fraction of outer cladding
Double_t ScintPhoton::fsFibreIndex = 1.59;      //index of refraction of fiber
Double_t ScintPhoton::fsFibreRadius = 0.6 * Munits::mm;  //radius of the WLS fiber (m)
Double_t ScintPhoton::fsScintIndex = 1.6;      //index of refraction of fiber
Double_t ScintPhoton::fsScintAttenLength = 0.33; // attenuation length in scintillator
Double_t ScintPhoton::fsTiO2ReflecDiffuse  = 0.93; // <tuned> diffuse reflection coefficient of TiO2
Double_t ScintPhoton::fsGrooveReflecCoeff = 0.5; // reflection coefficinet of the groove material
Double_t ScintPhoton::fsGrooveHalfWidth = 0.7 * Munits::mm; // 
Double_t ScintPhoton::fsReflectorScintReflec = 0.0; // Based upon guess from Leon
Int_t    ScintPhoton::fsLambertianReflection = 1; // 1= lambert, 0= isotropic.

const double kTooLong = 999. * Munits::meter; // "Fa, a long long way to run...."
                                              // (God help me, I can't believe I just put that in my comments.)

//----------------------------------------------------------------------
//default constructor 
ScintPhoton::ScintPhoton() : 
  DigiPhoton(),
  fRandom(gRandom),
  fStart(0.,0.,0.),
  fStartT(0.),
  fStartV(0.,0.,1.),
  fHalfHeight(0.005),
  fHalfLength(4.00),
  fHalfWidth(0.02083),
  fHitSurface(0),
  fInFibre(false),
  fSpecularReflected(false),
  fDiffuseReflected(false),
  fAttenuated(false),
  fAbsorbed(false),
  fGeomError(false),
  fBounces(0),
  fTrackLength(0.),
  fNormal(0,0,0),
  fScintAttenLength(fsScintAttenLength),
  fTiO2ReflecDiffuse(fsTiO2ReflecDiffuse)
{
}

//----------------------------------------------------------------------
//real constructor
ScintPhoton::ScintPhoton(TRandom* random,
                         const DigiScintHit* hit,
                         const UgliStripHandle& ush,
                         const TVector3& startPos, 
                         Double_t t,
                         const TVector3& startDir) : DigiPhoton()
{
  Reset( random, hit, ush, startPos, t, startDir );
}


//----------------------------------------------------------------------
//real constructor
ScintPhoton::ScintPhoton(TRandom* random,
                         const DigiScintHit* hit,
                         const UgliStripHandle& ush,
                         const TVector3& startPos, Double_t t) : DigiPhoton()
{
  TVector3 startDir(0,0,0);
  Reset( random, hit, ush, startPos, t, startDir );
}


//----------------------------------------------------------------------
//real constructor
ScintPhoton::ScintPhoton(TRandom* random,
                         const DigiScintHit* hit,
                         const UgliStripHandle& ush,
                         Double_t x, Double_t y, Double_t z, Double_t t,
                         Double_t u, Double_t v, Double_t w ) : DigiPhoton()
{
  TVector3 startPos(x,y,z);
  TVector3 startDir(u,v,w);
  Reset( random, hit, ush, startPos, t, startDir );
}




//----------------------------------------------------------------------
//Resetter.
void ScintPhoton::Reset(TRandom* random,
                        const DigiScintHit* hit,
                        const UgliStripHandle& ush,
                        const TVector3& startPos, Double_t t,
                        const TVector3& startDir)
{
  
  if(random) fRandom = random;
  else       fRandom = gRandom;

  SetParentHit(hit);

  fUgliStrip = ush;
  if(fUgliStrip.IsValid()) {
    fHalfHeight = fUgliStrip.GetHalfThickness();
    fHalfLength = fUgliStrip.GetHalfLength();
    fHalfWidth  = fUgliStrip.GetHalfWidth();
  } else {
    // Hardwired defaults, just in case.
    fHalfHeight=(0.005);
    fHalfLength=(4.00);
    fHalfWidth=(0.02083);
  }

  fHitSurface = 0;
  fInFibre = false;
  fSpecularReflected = false;
  fDiffuseReflected = false;
  fAttenuated = false;
  fAbsorbed = false; 
  fGeomError = false; 
  fBounces = 0;        // 
  fTrackLength = 0;

  fStart = startPos;
  fStartT = t;

  if(startDir.Mag2()==0) {
    // Roll a starting direction.
    double cosTheta = 1. - 2.*fRandom->Rndm();
    double sinTheta = TMath::Sin(TMath::ACos(cosTheta));
    double phi = 2.*TMath::Pi()*fRandom->Rndm();
    
    fStartV.SetXYZ( sinTheta*TMath::Cos(phi), 
                    sinTheta*TMath::Sin(phi),
                    cosTheta );
  } 

  // Set up the DigiPhoton.
  fPos =  fStart ;
  fDir =  fStartV;
  fT   =  fStartT;

  // Configuration.
  fScintAttenLength =fsScintAttenLength;
  fTiO2ReflecDiffuse =fsTiO2ReflecDiffuse;
}



//---------------------------------------------------------------------
ScintPhoton::~ScintPhoton()
{
}

void ScintPhoton::Configure( const Registry& r )
{
  //
  // Override any defaults with external configuration parameters.
  //
  Double_t tmpd;
  if (r.Get("ScintIndex",tmpd))   fsScintIndex  = tmpd;
  if (r.Get("FibreIndex",tmpd))   fsFibreIndex  = tmpd;
  if (r.Get("OuterCladdingIndex",tmpd))       fsOuterCladdingIndex  = tmpd;
  if (r.Get("InnerCladdingIndex",tmpd))       fsInnerCladdingIndex  = tmpd;
  
  if (r.Get("FibreRadius",tmpd))          fsFibreRadius  = tmpd;
  if (r.Get("ScintAttenLength",tmpd))     fsScintAttenLength  = tmpd;
  if (r.Get("TiO2ReflecDiffuse",tmpd))   fsTiO2ReflecDiffuse   = tmpd;
  if (r.Get("GrooveReflecCoeff",tmpd))    fsGrooveReflecCoeff  = tmpd;
  if (r.Get("GrooveHalfWidth",tmpd))      fsGrooveHalfWidth  = tmpd;
  if (r.Get("ReflectorScintReflec",tmpd))      fsReflectorScintReflec  = tmpd;   

  Int_t tmpi;
  if (r.Get("LambertianReflection",tmpi))      fsLambertianReflection  = tmpi;   
}



Bool_t ScintPhoton::PropagateOneLeg()
{
  // Propagate a ScintPhoton through the scintillator
  // to the next surface and evaluate if it survives.

  // Return true if everything is OK.

  fAttenuated = false;
  fSpecularReflected = false;
  fDiffuseReflected = false;
  fAttenuated = false;     // Caught in scintillator.
  fAbsorbed = false;       // Caught on wall.
  fGeomError = false;      // Too tricky for our stupid little code


  bool bounceOk = FindNextSurface(); // Propagate to next boundary.
  if(!bounceOk) {
    MSG("Photon",Msg::kWarning) << "Bounce failed." << endl;
    fGeomError = true;
    return false;
  }
    
  // See if the photon survived this leg.
  Attenuate();
  Move();
  
  if(fAttenuated) return true;

  // Interact with the surface.
  InteractAtSurface();
  
  // See if the photon made it to the scintillator
  if(fInFibre) return true;
  if(fAbsorbed) return true;
 
  // It got reflected, so compute the next bounce direction.
  if(fSpecularReflected || fDiffuseReflected) {
    FindReflectedDirection();
    return true;
  }
  
  return false;
}


Bool_t ScintPhoton::Propagate()
{
  // Propagate a ScintPhoton through the scintillator
  // until it is absorbed, attenuated, or hits the fibre.
  // Return false if the photon fails for some reason: indication
  // that the caller should not count this trial.

  fBounces = 0;
  while(true) {

    PropagateOneLeg();
    if(fGeomError) return false;
    if(fAttenuated) return true;
    if(fAbsorbed) return true;
    if(fInFibre) return true;

    // We got bounced.

    fBounces++;
  }
  
  return false;
}


void ScintPhoton::Print(Option_t* ) const
{
  cout << "ScintPhoton: (Units in cm)" << endl;
  cout << "  Started:  " << fStart.X() << "\t" << fStart.Y() << "\t" << fStart.Z() << endl;
  cout << "  StartDir: " << fStartV.X() << "\t" << fStartV.Y() << "\t" << fStartV.Z() << endl;

  cout << "  Pos:     " << X() << "\t" << Y() << "\t" << Z() << endl;
  cout << "  Dir After:  " << fDir.X() << "\t" << fDir.Y() << "\t" << fDir.Z() << endl;
  cout << "  Hit surface: " << fHitSurface << "\t Track length: " << fTrackLength << endl;
  if(fInFibre) cout << " In fibre. ";
  if(fSpecularReflected) cout << " SpecReflec. ";
  if(fDiffuseReflected) cout << " DiffuseReflec. ";
  if(fAttenuated) cout << " Attenuated. ";
  if(fAbsorbed) cout << " Absorbed. ";
  if(fGeomError) cout << " GeomError. ";
  cout << endl;
}


//---------------------------------------------------------------------
Bool_t ScintPhoton::FindNextSurface()
{
  //the photon cannot start off in the fiber or the groove as they dont
  //scintillate - for now let the groove be the same as the fiber - it 
  //basically is anyway.
  /*
    // let's ignore this.
    // it's a small correction... 
    // instead, we'll simply only allow incoming fibre hits, not outgoing
    // fibre hits.
  if((TMath::Power(fStartY,2.)
      +TMath::Power((fStartZ-(fHalfHeight-fsFibreRadius)),2.))<=TMath::Power(fsFibreRadius,2.)
     ){
    fInFibre = false;
    return false;
  }
  else if(TMath::Abs(fStartY)<=fsFibreRadius && fStartZ>=fHalfHeight-2.*fsFibreRadius){
    fInFibre = false;
    return false;
  }
  */

  //find which surface the photon would hit
  fTrackLength = kTooLong; // The distance to the nearest surface.
  fHitSurface = 0;
  double lengths[8];
  lengths[1] = DistanceToSurfaceI();
  lengths[2] = DistanceToSurfaceII();
  lengths[3] = DistanceToSurfaceIII();
  lengths[4] = DistanceToSurfaceIV();
  lengths[5] = DistanceToSurfaceV();
  lengths[6] = DistanceToSurfaceVI();
  lengths[7] = DistanceToFibre();
  
  for(int isurface = 1; isurface<=7; isurface++) {
    if( ( lengths[isurface] < fTrackLength ) && (lengths[isurface]>=0) ){
      fHitSurface = isurface;
      fTrackLength = lengths[isurface];
    }
  }
  
  if(fHitSurface == 0) {
    MSG("Photon",Msg::kWarning) << "Photon hit no surface!" << endl;
    fGeomError = true;
    return false;
  }
 
  return true;
}


Bool_t ScintPhoton::Attenuate()
{
  // Compute if photon was absorbed on this leg.

  // Figure out how far we go before being absorbed.
  // Doing it this way has the nice property that we
  // actually know where it stopped.
  double distance = fRandom->Exp(fScintAttenLength);

  if(distance < fTrackLength) {
    // Photon didn't make it.
    fAttenuated = true;
    // compute how far the photon DID make it.
    fTrackLength = distance;
  }

  return fAttenuated;
}


Bool_t ScintPhoton::Move()
{
  fPos += fDir*fTrackLength;
  fT += TimeOfFlight(fTrackLength);
  return true;
}

Double_t ScintPhoton::TimeOfFlight( Double_t distance )
{
  // Time for a photon to travel a distance in scintillator.
  return distance * fsScintIndex / Munits::c_light;
}


Bool_t ScintPhoton::InteractAtSurface()
{
  // 
  // Deals with whatever process happens at the surface.
  //
  // Sets the interaction flags and the surface normal.
  // Returns false if tracking error.

  fNormal.SetXYZ(0,0,0);
  
  // Fibre interaction:

  if(fHitSurface == 7) {
  // If the photon hit the fibre cladding, find out if it reflected or got absorbed.
    
    // The normal to the fibre surface is the norm of (y,z)-(fy,fz)
    // Surface normals at interaction point.
    fNormal.SetXYZ(0, 
                   Y(),  
                   Z()- (fHalfHeight-fsFibreRadius));
    fNormal.SetMag(1.0);

    // Dot product of direction vector and surface normal:
    double n_dot_v = fNormal.Dot(fDir);

    //the condition for reflection back into the scintillator is
    //if incident angle is > sin(n_2/n_1) where n_1 is the index of
    //refraction of the scintillator and n_2 is that for the fiber outer cladding
    //because the index for the inner cladding is greater than that for the outer
    //cladding, there can be no total reflection between outer and inner.
    
    // -ve sign because n points outwards.
  
    double incidentAngle = TMath::ACos(- n_dot_v);
    if(incidentAngle<TMath::ASin(fsOuterCladdingIndex/fsScintIndex)){
        
      //the photon made it into the fiber and that should be enough
      //for it to wavelength shift to a green photon - the mean free 
      //path of a photon in the fiber is of order microns and the 
      //WLS core of the fiber is 1.128 mm in diameter, or about 3 orders
      //of magnitude greater than the mean free path
      fInFibre = true;      

    } else {
      fSpecularReflected = true;      
    }

    return true;
  } // end surface == 7;


  // Was it a wall?
  if( ( fHitSurface >=1 ) && (fHitSurface <=4) ) {
    if( fHitSurface == 1 ) {
      // hit the top
      fNormal.SetXYZ(0,0,-1);
      
      // Check to see if it hit the groove.
      if( (TMath::Abs(Y()) < fsGrooveHalfWidth) ) { // In the groove. 
        if(fRandom->Rndm()>fsGrooveReflecCoeff)
          fSpecularReflected = true;
        else
          fAbsorbed = true;      
        return true;
      }
    } else if ( fHitSurface == 2 ) {
      fNormal.SetXYZ(0,-1,0);
    } else if ( fHitSurface == 3 ) {
      fNormal.SetXYZ(0,0,1);
    } else if ( fHitSurface == 4 ) {
      fNormal.SetXYZ(0,1,0);
    }

    // It hit TiO2 wall.
    double rand = fRandom->Rndm();
    if(rand < fTiO2ReflecDiffuse) {
      fDiffuseReflected = true;
    } else {
      fAbsorbed = true;
    }

    return true;
  }

  
  // Check to see if it hit an end.
  if( (fHitSurface == 5) || (fHitSurface == 6)) {

    // It hit the end. Was the end mirrored?
    //if( fHitSurface == 5 ):
    StripEnd::StripEnd_t end = StripEnd::kPositive;
    fNormal.SetXYZ(-1,0,0);   

    if( fHitSurface == 6 ) {
      fNormal.SetXYZ(1,0,0);
      end = StripEnd::kNegative;      
    }

    if( fUgliStrip.IsMirrored(end) ) {
      if(fRandom->Rndm()<=fsReflectorScintReflec) 
        fSpecularReflected = true;
      else 
        fAbsorbed = true;
    } else {
      fAbsorbed = true; // Flies off into the darkness of the pigtail.
    }
    
    return true;
  }

  MSG("Photon",Msg::kError) << "Error in InteractAtSurface: surface unrecognized." << endl;
  return false;
}

void ScintPhoton::FindReflectedDirection()
{
  // New direction after a bounce.
  
  // FIXME: specular reflection only.
  double n_dot_v = fNormal.Dot(fDir);
  
  if(fSpecularReflected) {
    fDir = fDir - 2.0*n_dot_v*fNormal;
    return;
  } 

  if(fDiffuseReflected) {
    // Find a second vector, b, at right angles to the normal.
    // Cheat here: we know all normals are either along X or perp to X
    
    TVector3 a = fNormal.Orthogonal();
    TVector3 b = fNormal.Cross(a);

    double cosTheta;
    if(fsLambertianReflection==0) {
      cosTheta = fRandom->Rndm(); // 0 to 1
    } else {
      cosTheta = sqrt(fRandom->Rndm());
    }
    double phi = 2.*TMath::Pi()*fRandom->Rndm(); // 0 to 2pi
    double sinTheta = TMath::Sqrt(1.0-cosTheta*cosTheta);
    double sinPhi = TMath::Sin(phi);
    double cosPhi = TMath::Cos(phi);
    
    fDir = fNormal*cosTheta 
      + a*sinTheta*cosPhi 
      + b*sinTheta*sinPhi;
    if(fDir.Mag2()==0) {
      cout << "mag0: " << fNormal.X() << "\t" << fNormal.Y() << "\t" << fNormal.Z() << endl;
      cout << "mag0: " << a.x() << "\t" << a.y() << "\t" << a.z() << endl;
      cout << "mag0: " << b.x() << "\t" << b.y() << "\t" << b.z() << endl;
    }
  }
}



Double_t ScintPhoton::DistanceToSurfaceI()
{
  // Find the distance from pos to surface I.
  // Return kTooLong if no interesection.
  if(fDir.Z()<=0) return kTooLong;
  double lambda = ( fHalfHeight - fPos.Z() )/fDir.Z();
  return lambda;
}

Double_t ScintPhoton::DistanceToSurfaceII()
{
  // Find the distance from pos to surface II.
  // Return kTooLong if no interesection.
  if(fDir.Y()<=0) return kTooLong;
  double lambda = ( fHalfWidth - fPos.Y() )/fDir.Y();
  return lambda;
}

Double_t ScintPhoton::DistanceToSurfaceIII()
{
  // Find the distance from pos to surface III.
  // Return kTooLong if no interesection.
  if(fDir.Z()>=0) return kTooLong;
  double lambda = ( -fHalfHeight - fPos.Z() )/fDir.Z();
  return lambda;
}


Double_t ScintPhoton::DistanceToSurfaceIV()
{
  // Find the distance from pos to surface IV.
  // Return kTooLong if no interesection.
  if(fDir.Y()>=0) return kTooLong;
  double lambda = ( -fHalfWidth - fPos.Y() )/fDir.Y();
  return lambda;
}


Double_t ScintPhoton::DistanceToSurfaceV()
{
  // Find the distance from pos to surface V.
  // Return kTooLong if no interesection.
  if(fDir.X() <= 0.) return kTooLong;
  double lambda = ( fHalfLength - fPos.X() )/fDir.X();
  return lambda;
}

Double_t ScintPhoton::DistanceToSurfaceVI()
{
  // Find the distance from fPos to surface VI.
  // Return kTooLong if no interesection.
  if(fDir.X() >= 0.) return kTooLong;
  double lambda = ( -fHalfLength - fPos.X() )/fDir.X();
  return lambda;
}


Double_t ScintPhoton::DistanceToFibre()
{
  // Find if the photon will intersect with the fibre,
  // and if so, where. Return kTooLong if it doesn't intersect.

  // the fibre position is (y=0, z=fHalfHeight-fsFibreRadius)
  double dy = fPos.Y();
  double dz = fPos.Z() - (fHalfHeight-fsFibreRadius);

  // Solve for : | (y,z) + lambda (vy,vz) - (0,h-r) |= r
  // where h is the height, and r is the radius.
  
  double v2  = (fDir.Y()*fDir.Y() + fDir.Z()*fDir.Z());
  if(v2==0) return kTooLong; // Going exactly in x-direction.

  double vdotdx = (fDir.Y() * dy) + (fDir.Z() * dz);
  double insqrt =  (vdotdx*vdotdx) - 
    (v2)*((dy*dy) + (dz*dz) - (fsFibreRadius*fsFibreRadius));

  // The track does not intersect the fibre ever.
  if(insqrt<0) return kTooLong;

  // the first solution is in the forward direction: take it.
  double soln1 = (-vdotdx - TMath::Sqrt(insqrt))/v2;             
  if(soln1>=0) return soln1;
  
  // only the first solution can be valid; the second solution
  // must apply only if we're already in the fibre.

  return kTooLong;
}

---