PhotonFastBlueModel.cxx

Go to the documentation of this file.

#include "PhotonFastBlueModel.h"
#include "PhotonTransportMaker.h"
#include "MessageService/MsgService.h"
#include "Plex/PlexHandle.h"
#include <cmath>

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

ClassImp(PhotonFastBlueModel)

PhotonTransportMakerProxy<PhotonFastBlueModel>
        gPhotonFastBlueModelProxy("photonFastBlueModel"); 

// Constructor
PhotonFastBlueModel::PhotonFastBlueModel() :
  fFastBlueSpectrumTable("PHOTONFASTBLUESPECTRUM")
{
  Configure(PhotonConfiguration()); // Set up default values.
  
  if (fFastBlueSpectrumTable.GetNumRows()<1) {
    MSG("Photon",Msg::kError) << "Can't find data for table PHOTONFASTBLUESPECTRUM Aborting!\n";
    assert(0);
  }
  SetBluePrescaleFactor(0.20676); // Lots of prescaling! This model is very efficient.
}

// Configure the model.
void 
PhotonFastBlueModel::Configure( const Registry& r )
{
  MSG("Photon",Msg::kDebug) << "Configuring FastBlueModelModel." << std::endl;
  double tmpd;
  if (r.Get("ScintDecayTime",tmpd))          fScintDecayTime = tmpd;
  if (r.Get("FibreRadius",tmpd))             fFibreRadius  = tmpd;
  if (r.Get("ReflectorScintReflec",tmpd))    fReflectorScintReflec  = tmpd;
  SetBluePrescaleFactor(0.20676); // Lots of prescaling! This model is very efficient.
}

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


Double_t PhotonFastBlueModel::GenTime(void)
{ 
  // choose a number from a double-exponential.
  /* Old.
 const double kLambda1 = 0.909973;
 const double kLambda2 = 0.012212;
 const double kArea1 = 377.615909;
 const double kArea2 = 124.790295;
  */

  const double kLambda1 = 0.846807;
  const double kLambda2 = 0.067829;
  const double kArea1 = 410.350066;
  const double kArea2 = 28.261422;

  float tt;
  if( fRandom->Rndm()*(kArea1+kArea2) < kArea1 ) 
    tt = fRandom->Exp(kLambda1);
  else 
    tt = fRandom->Exp(kLambda2);
  
  return tt * Munits::ns;
}

Double_t PhotonFastBlueModel::GenDx(void)
{
  /* OLD.
  const double kLorentzWidth1 = 0.020726;
  const double kArea1 = 5.726562;
  const double kLorentzWidth2 = 0.000385;
  const double kArea2 = 2.729155;
  const double kLambda3 = 0.051335;
  const double kArea3 = 0.836376;
  */
  
  const double kLorentzWidth1 = 0.017853;
  const double kArea1 = 8.178052;
  const double kLorentzWidth2 = 0.001309;
  const double kArea2 = 0.864215;
  const double kLambda3 = 1.000000;
  const double kArea3 = 0.000000;

  const double kArea = kArea1+kArea2+kArea3;

  double choose = fRandom->Rndm() * kArea;

  if( choose < kArea1 ) {
    return tan(fRandom->Rndm() * TMath::Pi()) * kLorentzWidth1;
  
  } else if ( choose < (kArea2+kArea1) ) {
    return tan(fRandom->Rndm() * TMath::Pi()) * kLorentzWidth2;
  
  } else {
    double absdx = fRandom->Exp(kLambda3);
    if(fRandom->Rndm()>0.5) return -absdx;
    return absdx;
  }
  return 0; // for compiler happiness
}

Bool_t 
PhotonFastBlueModel::ScintPhotonToFibreHit( const DigiScintHit*    hit, 
                                            const UgliStripHandle& inStrip,
                                            DigiPhoton*            &outPhoton )
{
  outPhoton = NULL;

  // Find the location/time for the start of the blue photon.
  double along = fRandom->Rndm();
  double tstart = hit->T1() + along*(hit->T2() - hit->T1()) + fRandom->Exp(fScintDecayTime);
  double xstart = hit->X1() + along*(hit->X2() - hit->X1());
  
  double dt, dx;
  do {
    dt = GenTime();
    dx = GenDx();
  } while (dt*Munits::c_light/1.6 < fabs(dx)); // Throw out unphysical solutions.
  double t = tstart + dt;
  double x = xstart + dx;

  double wavelength = fFastBlueSpectrumTable.Pick(fRandom->Rndm()) * Munits::nm;

  // So, we have our chosen x, t, lambda.
  // Make an arbitrary direction and position relative to the fibre.

  // Deal with the ends of the strip.
  if ( fabs(x) > inStrip.GetHalfLength() ) {

    // Which end did it disappare through?
    StripEnd::StripEnd_t end = (x>0)?StripEnd::kPositive:StripEnd::kNegative;
    if(inStrip.IsMirrored(end)) {
      if(fRandom->Rndm()>fReflectorScintReflec) {
        return false; // Photon didn't get reflected    
      } else {
        // Photon bounced off of reflector.
        if(x>0) 
          x = inStrip.GetHalfLength() - x;
        else
          x = -inStrip.GetHalfLength() -x;
      }

    } else {
      // end wasn't mirrored.
      return false; // Photon escapes.
    }   
  }

  if(inStrip.WlsBypass() > 0.) {
    // Deal with the case where the hit was near the edge of a two-part strip.
    // See if the hit is contained or not.
    if( ( (x + inStrip.GetHalfLength()) > inStrip.PartialLength(StripEnd::kNegative) ) 
        && ( (inStrip.GetHalfLength() - x) > inStrip.PartialLength(StripEnd::kPositive) )
        ) {
      // The photon got leaked out into the coil hole bypass
      return false;
    }

  }

  outPhoton = new DigiPhoton(hit,
                             wavelength, t,
                             x, 0,  inStrip.GetHalfThickness()-fFibreRadius*2.0,
                             0, 0,  1.0);
  
  
  return true;
}

---