PhotonDefaultModel.cxx

Go to the documentation of this file.

#include "PhotonDefaultModel.h"
#include "PhotonTransportMaker.h"
#include "ScintPhoton.h"
#include "MessageService/MsgService.h"
#include "Digitization/DigiSignal.h"
#include "ChannelNoiseRates.h"
#include <cmath>

CVSID( "$Id: PhotonDefaultModel.cxx,v 1.30 2008/11/19 17:23:26 rhatcher Exp $" );

ClassImp(PhotonDefaultModel)

PhotonTransportMakerProxy<PhotonDefaultModel>
        gPhotonDefaultProxy("photonDefaultModel"); 


// Constructor
PhotonDefaultModel::PhotonDefaultModel() 
{
  Configure(PhotonConfiguration()); // Set up default values.
}

// Configure the model.
void 
PhotonDefaultModel::Configure( const Registry& r )
{
  // The other classses in the model need configuring:
  ScintPhoton::Configure(r);

  MSG("Photon",Msg::kDebug) << "Configuring DefaultModel." << std::endl;

  double tmpd;
  int tmpi;
  // FIXME!
  if (r.Get("OverallLightOutput",tmpd))   fOverallLightOutput  = tmpd;
  if (r.Get("BirksConstant",tmpd))        fBirksConstant  = tmpd;
  if (r.Get("ScintDecayTime",tmpd))       fScintDecayTime = tmpd;
  if (r.Get("ClearFibreIndex",tmpd))      fClearFibreIndex  = tmpd;
  if (r.Get("FibreRadius",tmpd))          fFibreRadius  = tmpd;
  if (r.Get("FibreFlurorDecayTime",tmpd)) fFibreFlurorDecayTime  = tmpd;

  if (r.Get("OuterCladdingIndex",tmpd))   fOuterCladdingIndex  = tmpd;
  if (r.Get("InnerCladdingIndex",tmpd))   fInnerCladdingIndex  = tmpd;
  if (r.Get("FibreIndex",tmpd))           fFibreIndex  = tmpd;
  if (r.Get("FibreVelocityFudge",tmpd))   fFibreVelocityFudge  = tmpd;
  if (r.Get("FibreCoreRadius",tmpd))      fFibreCoreRadius  = tmpd;
  if (r.Get("ReflectorFibreReflec",tmpd)) fReflectorFibreReflec  = tmpd;
  if (r.Get("QuantumEfficiency",tmpd))    fQuantumEfficiency = tmpd;

  if (r.Get("FibreAttenN1",tmpd))        fFibreAttenN1 = tmpd;
  if (r.Get("FibreAttenN2",tmpd))        fFibreAttenN2 = tmpd;
  if (r.Get("FibreAttenLength1",tmpd))   fFibreAttenLength1 = tmpd;
  if (r.Get("FibreAttenLength2",tmpd))   fFibreAttenLength2 = tmpd;
  if (r.Get("ClearFibreAttenN1",tmpd))   fClearFibreAttenN1 = tmpd;
  if (r.Get("ClearFibreAttenN2",tmpd))   fClearFibreAttenN2 = tmpd;
  if (r.Get("ClearFibreAttenLength1",tmpd)) fClearFibreAttenLength1 = tmpd;
  if (r.Get("ClearFibreAttenLength2",tmpd)) fClearFibreAttenLength2 = tmpd ;

  if (r.Get("DarkNoiseRate",tmpd))        fDarkNoiseRate = tmpd; 
  if (r.Get("GreenNoiseRate",tmpd))       fGreenNoiseRate = tmpd;
  if (r.Get("NoiseWindow",tmpd))          fNoiseWindow = tmpd;
  if (r.Get("PoissonNoiseMean",tmpd))     fPoissonNoiseMean = tmpd;
  if (r.Get("ExpTailFraction",tmpd))      fExpTailFraction = tmpd;
  if (r.Get("FDNoiseRate",tmpi))          fDetectorNoiseRate = tmpi;

  if (r.Get("UseSimpleNoiseModel",tmpi))  fUseSimpleNoiseModel = tmpi;

  // Set the prescale factors.
  SetBluePrescaleFactor(1.0); // No prescaling.
  SetWlsPrescaleFactor(0.2);  // We always accept a contained ray, no matter what.
  SetGreenPrescaleFactor(fQuantumEfficiency); // Phototube photocathode acceptance.


}

void 
PhotonDefaultModel::Reset( const VldContext& )
{
  // Ensure tables are up-to-date.
  
  // VldContext should be set by now, so can set up noise rates

  // Set up the channel noise rates if in far detector
  if ( (fContext.GetDetector() == Detector::kFar) 
       && fUseSimpleNoiseModel==0 ) {
    DbiResultPtr<ChannelNoiseRates> resPtr(fContext);
    if (!resPtr.GetNumRows())
      MSG("Photon", Msg::kFatal) << "Got no DB rows with noise rates!" 
                                 << endl;
    
    for (UInt_t irow=0; irow < resPtr.GetNumRows(); ++irow) {
      const ChannelNoiseRates* cnr = resPtr.GetRow(irow);
      chanhitsmap[cnr->GetChAdd()+5329*cnr->GetCrate()]=cnr->GetNormRate();
    }
  }
}

// Computes the number of photons we'll need to track for the 
// hit. See the PhotonCount class for the various ways this
// can be calculated.

Bool_t 
PhotonDefaultModel::ComputePhotons( const DigiScintHit* inHit,
                                    Double_t inBluePrescale,
                                    Double_t inWlsPrescale,
                                    Double_t inGreenPrescale,
                                    PhotonCount& outCount )
{
  // Number of photons made:
  if(!inHit) return false;
  double dE = inHit->DE();
  double dx = inHit->DS();

  double dEdx = dE/dx;
  
  //get the number of photons
  //from the PDG 2002 Section 27 on Particle Detectors you get
  //about 1 photon/100eV of energy deposited.
  // .. so, 10photons/keV = 10^7 photons/GeV

  // From rough comparison with REROOT results, the total 
  // constant for this algorithm should be ~27500

  // Here's what I can account for:
  const double kLightYield = ( 1/100.  // photons/eV
                               * 1e9 // GeV/ev
                               * inBluePrescale
                               * inWlsPrescale
                               * inGreenPrescale
                               );

  // This fudge factor comes from a comparison of using GMINOS FLS digits
  // and this method.  It has no physical justification.
  const double kFudgeFactor = 2.0;
                        
  // Apply normalization and birks' constant.
  double meanBlue = ( kLightYield * kFudgeFactor * fOverallLightOutput 
                      * dE/(1. + fBirksConstant*dEdx)   );

  // Set to a poisson number about this mean.
  int nblue =  fRandom->Poisson(meanBlue);
    
  //MSG("Photon",Msg::kDebug) << "DefaultComputer: " << dE/Munits::MeV 
  //                        << " Mev -> " << meanBlue << " expected blue -> " 
  //                        << nblue << " blue photons.  " 
  //                        << inBluePrescale << "  " 
  //                        << inWlsPrescale << "  "
  //                        << inGreenPrescale << "  "
  //                        << endl;
  

  outCount.SetBlueTotal( nblue );
  return true;
}

// Creates a single scintillator photon and tracks it to a fibre hit.
// Returns FALSE if the photon doesn't make it into the fibre.
// inStrip is the strip geometry.
// inBlueStart is the blue photon start position in strip coords
// inBlueTime is the time the photon starts, in Munits
// outBlueWavelength is the returned wavelength of the blue photon
//   (which may be zero if this model doesn't bother to track it)
// outHitPos is the X-position along the fibre of the hit
// outHitTime is the time the blue photon hit the fibre.
Bool_t 
PhotonDefaultModel::ScintPhotonToFibreHit( const DigiScintHit*    hit, 
                                           const UgliStripHandle& inStrip,
                                           DigiPhoton*            &outPhoton )
{
  outPhoton = NULL;
  // This default method uses the ScintPhoton class to do the dirty work.

  
  if(!hit) return false;

  // Make the photon.
  ScintPhoton*  photon = new ScintPhoton();  

  // Find the location/time for the blue photon.

  Double_t delta = fRandom->Rndm();
  TVector3 x1(hit->X1(),hit->Y1(),hit->Z1());
  TVector3 x2(hit->X2(),hit->Y2(),hit->Z2());

  TVector3 p = x1 + delta*(x2-x1);
  Double_t t = hit->T1() + delta * (hit->T2() - hit->T1())
    + fRandom->Exp(fScintDecayTime);
  

  // The point of this loop is to ensure that we actually manage to 
  // simulate a photon successfully.. it doesn't just disappear, but
  // is successfully tracked to SOMEWHERE.
  bool ok;
  do {
    TVector3 dir(0,0,0); // let the scint photon do it.
    photon->Reset(fRandom,
                  hit,
                  inStrip,
                  p, t,
                  dir);
    
    photon->Propagate();
    ok = ! photon->GeomError(); // Geometry errors are bad.
  } while(!ok);

  // Successful photon.. may or may not have hit green.
  outPhoton = (DigiPhoton*) photon;

  // Simple check for intra-strip light leakage. Not completely correct; this SHOULD
  // be part of the photon tracking code.. but simple and straightforward.
  if(inStrip.WlsBypass() > 0.) {
    float x = outPhoton->X();
    // 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;
    }
  }

  return (photon->InFibre()); 
}


// Makes a green photon from a blue photon.
// Returns true if a green photon is made and trapped, false if lost.
// inStrip is the strip geometry,
// inBlueWavelength is the wavelength of the blue photon (may be 0)
// inHitPos is the position along the fibre.
// inHitTime is the time of the hit
// outGreenTime is the returned time of the photon start(after decay)
// outGreenCosAngle is the X-direction cosine of the emitted photon 
//     (which can be used to calculate distance traveled in the green fibre)
// outGreenWavelength is the wavelength of the green photon
//     (which may be 0 if this model doesn't bother to calculate it)
Bool_t 
PhotonDefaultModel::FibreHitToGreenPhoton( const UgliStripHandle& inStrip,
                                           StripEnd::StripEnd_t   inDirection,
                                           DigiPhoton*            inBluePhoton,
                                           DigiPhoton*            &outGreenPhoton
                                           )
{
  // Simple model. No wavelength dependencies.
  // Return a contained green photon.
  // This code currently does not do counting correctly; it just keeps
  // simulating until it gets a green photon that is contained in the fibre.
  // This implicitly assumes there is no correllation between the blue photon
  // and the green photon in any way: one blue photon is as good as another.
  // This is both computationally efficient and pretty true.
  // The blue photon is assumed to uniformly in the green fibre.


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

  // Loop until we have a valid solution.
  TVector3 dir;
  TVector3 x0;
  while(true) {

    // Find the new direction of the green photon.
    // Go in the forward or backward direction, depending on the 
    // given parameter.
    
    // Without loss of generality, rotate the system coordinates about the X-axis
    // such that the photon start position it on the Z axis between 0 and R.
    
    // Chose a position.
    // FIXME: really, the position should be chosen by 
   // computing the depth the blue photon reaches in the fibre.
    // But frankly, who cares?

    // Chose radial position as volume-weighted.
    double R = fFibreCoreRadius;
    double r = sqrt(fRandom->Rndm())*R;
    x0.SetXYZ(inBluePhoton->X(),0,r);

    // Choose a direction.
    double cosTheta  = 2.*fRandom->Rndm()-1.0;
    double sinTheta = sqrt(1.0-cosTheta*cosTheta); //TMath::Sin(TMath::ACos(cosTheta));
    double phi = 2.*TMath::Pi()*fRandom->Rndm();
    dir.SetXYZ( cosTheta, 
                sinTheta*TMath::Cos(phi),
                sinTheta*TMath::Sin(phi) );
    
    // Find the position of first reflection:
    double xv = r*dir.z();         // xv = x0.y()*dir.y() + x0.z()*dir.z();
    double v2 = sinTheta*sinTheta; //  double v2 = dir.y()*dir.y() + dir.z()*dir.z();
    if(v2==0) return false;
    
    double lambda = (-xv +  sqrt(xv*xv + v2*(R*R-r*r)) )/v2;
    
    TVector3 xhit = x0 + lambda*dir;
    
    // The surface normal at this place.
    TVector3 n(0, -xhit.y()/R, -xhit.z()/R);

    // Find angle of incidence.
    double cosAngle =  n.Dot(dir);
    double sinAngle = sqrt(1.-cosAngle*cosAngle); //sin(acos(cosAngle));
    if(sinAngle > (fOuterCladdingIndex/fFibreIndex)) {
      // Solution is good!
      break;
    } else {
      // Return if we want to do counting properly, or if
      // we have some sort of correlation between green and blue.
      // Otherwise, why waste a perfectly good blue photon? Just
      // keep trying till it bloody well DOES go down the fibre.
      // return false;
    }
  }

  // Choose photon direction:
  if(inStrip.IsMirrored(StripEnd::kPositive)||inStrip.IsMirrored(StripEnd::kNegative)) {
    // If the strip is mirrored, we _don't_ want to make an a priori decision
    // since either direction can make a hit on the correct side.

  } else {

    // Let's ensure it's going in the demanded 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;    
    }
  }
  
  // There's another (unlikely) condition: the photon bounces around exactly sideways
  // and so never gets to the end of the tube.
  if(dir.X() ==0) return false;

  // Otherwise.. we have a winner.

  // Set up the time correctly now.
  Double_t decayTime = fRandom->Exp(fFibreFlurorDecayTime);

  outGreenPhoton = new DigiPhoton( inBluePhoton->ParentHit(),
                                   0, // Wavelength is generic.
                                   inBluePhoton->T() + decayTime,
                                   x0,
                                   dir
                                   );
                                   

  return true;
}


// Creates a green photon from a fibre hit, and tracks to the PMT.
// Returns true if green photon made it, false if attenuated away.
// Takes the blue wavelength (may not be used), position, and time,
// and fills a DigiPE object.
Bool_t 
PhotonDefaultModel::GreenPhotonToPe( const UgliStripHandle& inStrip,
                                     DigiPhoton*         inGreenPhoton,
                                     DigiPE*             &outPe,
                                     StripEnd::StripEnd_t &outEnd)
{
  outPe = NULL;
  if(!inGreenPhoton) return false;

  // First, figure out which way we're going.
  StripEnd::StripEnd_t dir = 
    (inGreenPhoton->GetCosX()>0) ? StripEnd::kPositive : StripEnd::kNegative;

  StripEnd::StripEnd_t readout = dir; // Assume we're reading out the end it's going.

  // Is the end we've gotten to mirrored? 
  if(inStrip.IsMirrored(dir)) {
    // Ok, then we're going to need to flip it.
    // Does this photon survive the flip?
    if(fRandom->Rndm() > fReflectorFibreReflec) {
      return false;
    }

    // We survived the reflection.
    readout = (dir==StripEnd::kPositive) ? StripEnd::kNegative : StripEnd::kPositive;
  }

  // The photon is reflecting it's way down the strip.
  // Figure out how far it has to go in each medium.

  double distToCenter = (inGreenPhoton->GetCosX()>0) ? (-inGreenPhoton->X())
    : (inGreenPhoton->X());

  Double_t greenDist =
    distToCenter // Distance to center of strip.
    +inStrip.GetHalfLength()         // Distance to end of strip;
    +inStrip.WlsPigtail(readout); // Pigtail on readout end.

  //MSG("Photon",Msg::kVerbose) << "Distance to center of strip: " << distToCenter 
  //                          << "  end= " << readout 
  //                          << "  halflen=" << inStrip.GetHalfLength()
  //                          << "  x="<<inGreenPhoton->X() 
  //                          << "  dx=" << inGreenPhoton->GetCosX() << endl;
 
  // Deal with WLS bypass, if any
  double bypassExtra = 0;
  if(inStrip.WlsBypass() >0) {    
    bypassExtra = ( inStrip.PartialLength(StripEnd::kNegative) + 
                    inStrip.PartialLength(StripEnd::kPositive) +
                    inStrip.WlsBypass() - 
                    inStrip.GetHalfLength()*2.0 );

    // place in the middle of the break;
    float x_break =0.5* (inStrip.PartialLength(StripEnd::kNegative) 
                         - inStrip.PartialLength(StripEnd::kPositive) );

    if(inGreenPhoton->X() < x_break) {
      // If the photon is at -ve end and going +ve, it goes through the bypass
      if(dir==StripEnd::kPositive)  greenDist+= bypassExtra;
    } else {
      // If the photon is at +ve end and going -ve, it goes through the bypass
      if(dir==StripEnd::kNegative)  greenDist+= bypassExtra;
    }
  }

  if(readout!=dir) {
    // i.e. we're going towards the mirror end but being read out on the non-mirror end
    greenDist+= inStrip.GetHalfLength()*2.0; // Gotta double back...
    greenDist+= bypassExtra;                 // ... through the ENTIRE green fibre
    greenDist+= inStrip.WlsPigtail(dir)*2.0; // And we have to go through the pigtail, if any
                                             // (N.B. This is for caldet, which has pigtails before the reflector)
  }


  // Find the distance down the clear strip.
  Double_t clearDist = inStrip.ClearFiber(readout);

  // But... we don't travel in a straight line!
  // We bounce around the fibre a lot. This makes the path length longer.
  greenDist /= fabs(inGreenPhoton->GetCosX());
  clearDist /= fabs(inGreenPhoton->GetCosX());
  
  //MSG("Photon",Msg::kVerbose) << "After pathlength correction: greenDist = " << greenDist
  //                          << "  clearDist = " << clearDist << endl;

  // Compute the probablility we just lose it.
  // FIXME: could use full treatment, and at the very least
  // shouldn't be hard-coded.
  // This treatment comes straight out of GMINOS: init/scintfiber_postffr.F
  double norm_g = 1.0/(fFibreAttenN1+fFibreAttenN2);

  // 1 is no attenuation, 0 is complete loss
  // So, a small number is a small attenuation prob.
  double greenAtten = norm_g * ( fFibreAttenN1 * exp(-greenDist/fFibreAttenLength1) +
                                 fFibreAttenN2 * exp(-greenDist/fFibreAttenLength2) );

  //MSG("Photon",Msg::kVerbose) << "Green atten prob = " << greenAtten << endl;
  if(fRandom->Rndm() > greenAtten) return false;

  // Clear fibres:
  // Again, these numbers are from the GMINOS defaults.
  double norm_c =  1.0/(fClearFibreAttenN1+fClearFibreAttenN2);

  // 1 is no attenuation, 0 is complete loss
  double clearAtten = norm_c * ( fClearFibreAttenN1 * exp(-clearDist/fClearFibreAttenLength1) +
                                 fClearFibreAttenN2 * exp(-clearDist/fClearFibreAttenLength2) );
  
  //MSG("Photon",Msg::kVerbose) << "Clear atten prob = " << clearAtten << endl;
  if(fRandom->Rndm() > clearAtten) return false;

  // ASSUME: no loss at connector interfaces

  // PMT: ASSUME:
  // FIXME: Assume 100% quantum efficiency, or efficiency dealt with
  // in photon computation.

  // Build the resultant digipe.
  // Find the time of arrival.
  double time = inGreenPhoton->T()
    + greenDist * fFibreIndex * fFibreVelocityFudge / Munits::c_light
    + clearDist * fClearFibreIndex * fFibreVelocityFudge / Munits::c_light;

  PlexStripEndId seid =  inStrip.GetSEId();
  seid.SetEnd(readout);

  PlexHandle plex(fContext);
  PlexPixelSpotId psid = plex.GetPixelSpotId(seid);

  outPe = new DigiPE( time, psid, inGreenPhoton->ParentHit() );
  outEnd = readout;

  return true;  
}

void   
PhotonDefaultModel::MakeNoise( std::vector<DigiPE*> &peList,
                               Double_t t_start,
                               Double_t t_end )
{


  // So far, this only works for the far detector, so if in ND (or the
  // use specifically requested it), use old method:
  if ( (fContext.GetDetector() != Detector::kFar) 
       || fUseSimpleNoiseModel ) {
    MakeNoiseOld(peList, t_start, t_end);
    return;
  }
  
  // Work out our window.
  double t1 = t_start - fNoiseWindow*0.5;
  double t2 = t_end   + fNoiseWindow*0.5;
  double dt = t2-t1;

  MSG("Photon",Msg::kDebug) << "Computing noise for window at" << t1 
                           << " s for " << (t2-t1)/Munits::ns << " ns."
                           << std::endl;

  PlexHandle plex(fContext);
  
  const std::vector<PlexStripEndId>& stripEnds = plex.GetAllStripEnds();
  
  // How many fibre hits would we expect?
  Double_t fibreHits = dt * fDetectorNoiseRate;
  // Now find the number that DID happen...
  Int_t nGreenHits = fRandom->Poisson(fibreHits);
  MSG("Photon",Msg::kDebug) 
    << "Creating " << nGreenHits << " green WLS noise hits.\n";
  
  // Create the hits.
  for(int i=0; i<nGreenHits; i++) {
    // Pick a strip end at random.
    Int_t whichSeid = fRandom->Integer((int)stripEnds.size());
    // Decide whether or not to put this hit on this channel based on
    // its normalized noise rate
    while ( fRandom->Uniform(0,1) > 
            GetStripEndRate(stripEnds[whichSeid] , plex) )
      whichSeid = fRandom->Integer((int)stripEnds.size());

    PlexPixelSpotId psid = plex.GetPixelSpotId(stripEnds[whichSeid]);
    if(psid.IsValid()) {
      Int_t npe = 0;
      // Discard the 0 bin
      while (npe==0) npe=fRandom->Poisson(fPoissonNoiseMean);
      Double_t hitTime = fRandom->Uniform(t1,t2);
      for (int j=0; j<npe; ++j) {
        DigiPE* pe = new DigiPE( hitTime,
                                 psid,
                                 DigiSignal::kDarkNoise ); // No Scint Hit associated with the PE.
        peList.push_back(pe);
      }
    }
  }
  // Do the same for the exponential tail
  Int_t nExpHits = fRandom->Poisson(fibreHits*fExpTailFraction);
  for (Int_t i=0; i<nExpHits; ++i) {
    // Pick a strip end at random.
    Int_t whichSeid = fRandom->Integer((int)stripEnds.size());
    // Decide whether or not to put this hit on this channel based on
    // its normalized noise rate
    while ( fRandom->Uniform(0,1) > GetStripEndRate(stripEnds[whichSeid] , plex))
      whichSeid = fRandom->Integer((int)stripEnds.size());
    
    PlexPixelSpotId psid = plex.GetPixelSpotId(stripEnds[whichSeid]);

    Int_t npe = 4+TMath::Nint(fRandom->Exp(0.7));
    Double_t hitTime = fRandom->Uniform(t1,t2);
    for (int j=0; j<npe; ++j) {
      DigiPE* pe = new DigiPE( hitTime,
                               psid,
                               DigiSignal::kFibreLight ); // No Scint Hit associated with the PE.
      peList.push_back(pe);
    }
  }
}

Float_t 
PhotonDefaultModel::GetStripEndRate(PlexStripEndId seid, 
                                             PlexHandle ph)
{
  // We could perhaps speed this up by storing the rate by *strip end*
  // in the config file, instead of the rate by channel (saves asking
  // the plex)
  RawChannelId rcid = ph.GetRawChannelId(seid);
  UShort_t chadd = rcid.GetChAdd();
  Int_t crate = rcid.GetCrate();
  MSG("Photon",Msg::kDebug) << "chanhitsmap[" << chadd+5329*crate 
                            << "] = " << chanhitsmap[chadd+5329*crate] 
                            << endl;
  return chanhitsmap[chadd+5329*crate];
}

void   
PhotonDefaultModel::MakeNoiseOld( std::vector<DigiPE*> &peList,
                               Double_t t_start,
                               Double_t t_end )
{

  // Work out our window.
  double t1 = t_start - fNoiseWindow*0.5;
  double t2 = t_end   + fNoiseWindow*0.5;
  double dt = t2-t1;

  MSG("Photon",Msg::kDebug) << "Computing noise using simple noise"
                            << " model for window at" << t1 
                            << " s for " << (t2-t1)/Munits::ns << " ns."
                            << std::endl;

  PlexHandle plex(fContext);
  // Get a list of PMTs from the Plex.
  const std::vector<PlexPixelSpotId>& tubes = plex.GetAllTubes();
  
  // How many PMT hits would we expect?
  Double_t pmtHits = (double)tubes.size() * dt * fDarkNoiseRate;
  // Now find the number that DID happen...
  Int_t nPmtHits = fRandom->Poisson(pmtHits);
  MSG("Photon",Msg::kDebug) << "Creating " << nPmtHits 
                            << " dark noise hits.\n";
  
  // Create the hits.
  for(int i=0; i<nPmtHits; i++) {
    // Pick a tube at random.
    Int_t whichPmt = fRandom->Integer((int)tubes.size());
    // Pick a pixelspot.
    PlexPixelSpotId psid = tubes[whichPmt];
    if(psid.GetElecType()==ElecType::kQIE) {
      psid.SetPixel(fRandom->Integer(64));
      psid.SetSpot(1);      
    }
    if(psid.GetElecType()==ElecType::kVA) {
      psid.SetPixel(fRandom->Integer(16));
      psid.SetSpot(fRandom->Integer(8)+1);      
    }
    if(psid.IsValid()) {
      DigiPE* pe = new DigiPE( fRandom->Uniform(t1,t2),
                               psid,
                               DigiSignal::kDarkNoise ); // No Scint Hit associated with the PE.
      peList.push_back(pe);
    }
  }

  // Get list of strip ends from the Plex.  
  const std::vector<PlexStripEndId>& stripEnds = plex.GetAllStripEnds();

  // How many fibre hits would we expect?
  Double_t fibreHits = (double)stripEnds.size() * dt * fGreenNoiseRate;
  // Now find the number that DID happen...
  Int_t nGreenHits = fRandom->Poisson(fibreHits);
  MSG("Photon",Msg::kDebug) << "Creating " << nGreenHits 
                            << " green WLS noise hits.\n";
  
  // Create the hits.
  for(int i=0; i<nGreenHits; i++) {
    // Pick a strip end at random.
    Int_t whichSeid = fRandom->Integer((int)stripEnds.size());
    
    PlexPixelSpotId psid = plex.GetPixelSpotId(stripEnds[whichSeid]);

    if(psid.IsValid()) {  
      DigiPE* pe =  new DigiPE( fRandom->Uniform(t1,t2),
                                psid,
                                DigiSignal::kFibreLight ); // No Scint Hit associated with the PE.
      peList.push_back(pe);
    }
  }
}


void PhotonDefaultModel::MakeAfterpulses( const std::vector<DigiPE*>& /*inPeList*/,
                                             std::vector<DigiPE*>& /*outPeList*/ )
{

  return;
}




void PhotonDefaultModel::Print(Option_t*option) const
{
  if(option[0]=='c') {
    MSG("Photon",Msg::kInfo) << "-<PhotonComputer>     PhotonDefaultModel:-----" << endl;
    MSG("Photon",Msg::kInfo) << "  OverallLightOutput:   " << fOverallLightOutput << endl;
    MSG("Photon",Msg::kInfo) << "  Birk's Constant:      " << fBirksConstant << endl;   
  }
  if(option[0]=='b') {
    MSG("Photon",Msg::kInfo) << "-<BluePhotonTracker>  PhotonDefaultModel:-----" << endl;
    //FIXME
  }
  if(option[0]=='w') {
    MSG("Photon",Msg::kInfo) << "-<WlsFibreModel>      PhotonDefaultModel:-----" << endl;  
    // FIXME
  }
  if(option[0]=='g') {
    MSG("Photon",Msg::kInfo) << "-<GreenPhotonTracker> PhotonDefaultModel:-----" << endl;  
    MSG("Photon",Msg::kInfo) << "  ReflectorFibreReflec: " << fReflectorFibreReflec << endl; 
    MSG("Photon",Msg::kInfo) << "  Quantum Efficiency:   " << fQuantumEfficiency    << endl; 
    
  }
  if(option[0]=='n') {
    MSG("Photon",Msg::kInfo) << "-<NoiseMaker>          PhotonDefaultModel:-----" << endl;  
    //FIXME
  }
}

---