PhotonTransport.cxx

Go to the documentation of this file.

// $Id: PhotonTransport.cxx,v 1.16 2009/05/19 18:19:17 kasahara Exp $
//
// FILL_IN: [Document your code!!]
//
// n.tagg1@physics.ox.ac.uk
#include "PhotonTransport.h"

#include <TRandom.h>
#include <TRandom3.h>

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro

#include "Digitization/DigiScintHit.h"
#include "Digitization/DigiPE.h"
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"

#include "PhotonTransportModule.h"
#include "PhotonTransportMaker.h"
#include "PhotonConfiguration.h"

#include "Plex/PlexHandle.h"

#include "Calibrator/Calibrator.h"

#include "Conventions/MinosMaterial.h"

#include "Util/LoadMinosPDG.h"

#include <TFile.h>
#include <TNtuple.h>
#include <TNtupleD.h>

#include <cmath>

JOBMODULE(PhotonTransport, "PhotonTransport",
          "Photon transport from energy deposition to photoelectroncs");
CVSID("$Id: PhotonTransport.cxx,v 1.16 2009/05/19 18:19:17 kasahara Exp $");


//......................................................................

PhotonTransport::PhotonTransport() :
  fName_PhotonComputer("photonDefaultModel"),   // These are NOT the true defaults.. See PhotonConfiguration
  fName_BluePhotonTracker("photonDefaultModel"),
  fName_WlsFibreModel("photonDefaultModel"),
  fName_GreenPhotonTracker("photonDefaultModel"),
  fName_NoiseMaker("photonDefaultModel"),
  fName_Afterpulser("photonDefaultModel"),
  fDebugFileName(""),
  fDebugMask(0xFFFF),
  fRandomSeed(0),
  fRandom(new TRandom3),
  fProtonRangeTable("PHOTONPROTONRANGE"),
  fElectronRangeTable("PHOTONELECTRONRANGE"),
  fPhotonComputer(0),
  fBluePhotonTracker(0),
  fWlsFibreModel(0),
  fGreenPhotonTracker(0),
  fNoiseMaker(0),
  fAfterpulser(0),
  fDebugFile(0),
  fDebugNtuple(0),
  fBadHitNtuple(0)
{

  LoadMinosPDG(); // initialize MINOS specific particles

  CreateModelObjects(); // Create the models.
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================  

}

//......................................................................

PhotonTransport::~PhotonTransport() 
{
//======================================================================
// FILL_IN: [Document your code!!]
//======================================================================
  if(fPhotonComputer) delete fPhotonComputer;
  if(fBluePhotonTracker) delete fBluePhotonTracker;
  if(fWlsFibreModel) delete fWlsFibreModel;
  if(fGreenPhotonTracker) delete fGreenPhotonTracker;
  if(fNoiseMaker) delete fNoiseMaker;
  if(fAfterpulser) delete fAfterpulser;
  if(fRandom) delete fRandom;
  if(fDebugFile){
    TDirectory* savedir = gDirectory;
    fDebugFile->cd();
    if ( fDebugNtuple ) fDebugNtuple->Write(NULL,TObject::kOverwrite);
    if ( fBadHitNtuple ) fBadHitNtuple->Write(NULL,TObject::kOverwrite);
    fDebugFile->Close();
    delete fDebugFile;
    savedir -> cd();
  }
}

//......................................................................

JobCResult PhotonTransport::Ana(const MomNavigator* mom )
{
  SimSnarlRecord* simrec = dynamic_cast<SimSnarlRecord*>
    (mom->GetFragment("SimSnarlRecord"));
  if ( !simrec ) {
    MSG("Photon",Msg::kWarning) << "No SimSnarlRecord in Mom" << endl;
    return JobCResult::kFailed;
  }
  
  const PhotonEventResult* per = dynamic_cast<const PhotonEventResult*>(simrec->FindComponent("PhotonEventResult"));
  if(per) per->Print();
  else MSG("Photon",Msg::kWarning) << "Can't find PhotonEventResult" << endl;
  // pass record to mom to own



  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

JobCResult PhotonTransport::Reco(MomNavigator* mom)
{
  return Ana(mom);
}

//......................................................................

JobCResult PhotonTransport::Get(MomNavigator* mom)
{
//======================================================================
// Does the work. 
// Creates photons at the scintHit vertex, then propagates them to the 
// phototube.
//======================================================================
  SimSnarlRecord* simsnarl = 0;
  TObject* tobj;
  TIter    fragiter = mom->FragmentIter();

  fPeTotal = 0;

  // Get the simsnarl.
  while( ( tobj = fragiter.Next() ) ) {
    simsnarl = dynamic_cast<SimSnarlRecord*>(tobj);
    if(simsnarl) break;
  }
  
  // Verify there IS a simsnarl.
  if(!simsnarl) {
    MSG("Photon",Msg::kError) << "No SimSnarl found. You must run RerootToTruthModule()!" << endl;
    return JobCResult::kFailed;
  }

  RecJobHistory& jobhist 
               = const_cast<RecJobHistory&>(simsnarl->GetJobHistory());
  jobhist.CreateJobRecord(RecJobHistory::kPhotonTransport);

  fSimHeader = simsnarl->GetSimSnarlHeader();
  if(fSimHeader ==0){
    MSG("Photon",Msg::kError) << "Cannot find SimSnarlHeader in SimSnarl." << endl;
    return JobCResult::kFailed;
  }

  // Set the random seed to match the event.
  // Also add some user-configurable randomness, in case they want to
  // generate the same event with a different outcome.
  fRandom->SetSeed( fSimHeader->GetRun() 
                    + fSimHeader->GetSnarl() 
                    + fRandomSeed );
  
  // We need a context to do table lookups
  VldContext simContext = fSimHeader->GetVldContext();

  // But.. we want our context to be "MC" not "Reroot". Ensure that's correct.
  fContext = VldContext(simContext.GetDetector(),SimFlag::kMC,simContext.GetTimeStamp());

  Calibrator::Instance().Reset(fContext);

  // Create the stat report.
  fEventResult.Reset();

  JobCResult res;

  // Get the DigiScintHits.
  fHitList = 
    dynamic_cast<const TObjArray*>(simsnarl->FindComponent(0,"DigiScintHits"));
  if(fHitList) {

    // Delete any pre-existing array of DigiPE.
    if( dynamic_cast<const TObjArray*>(simsnarl->FindTemporary(0,"DigiListPe")) != 0) {
      MSG("Photon",Msg::kError) << "Found pre-existing DigPE list. Aborting." << std::endl;
      return JobCResult::kFailed;
    } else {
      fPeList = new TObjArray(0);
      fPeList->SetName("DigiListPe");
      fPeList->SetOwner(true);
    }
    
    // Do the actual work.
    res = SimulateEvent();

  } else {

    MSG("Photon",Msg::kError) << "Can't find scint hit array.\n"; 
    return JobCResult::kFailed;  
    
  };

  simsnarl->AdoptTemporary(fPeList);

  // Add the event result.
  simsnarl->AdoptComponent(new PhotonEventResult(fEventResult));
  MSG("Photon",Msg::kSynopsis) << fEventResult;

  return res; // kNoDecision, kFailed, etc.
}

//......................................................................

const Registry& PhotonTransport::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================

  return PhotonConfiguration();
}

//......................................................................

void PhotonTransport::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================

  MSG("Photon",Msg::kDebug) << "PhotonTransport::Config" << endl;

  r.Get("randomSeed",fRandomSeed);

  const char* str;
  Bool_t modelChanged = false;
  if( r.Get("PhotonComputer",    str))
    if(fName_PhotonComputer != std::string(str))
      { modelChanged=true; fName_PhotonComputer = str; }

  if( r.Get("BluePhotonTracker", str))
    if(fName_BluePhotonTracker != std::string(str))
      { modelChanged=true; fName_BluePhotonTracker = str; }
  
  if( r.Get("WlsFibreModel",     str))
    if(fName_WlsFibreModel != std::string(str))
      { modelChanged=true; fName_WlsFibreModel = str; }
  
  if( r.Get("GreenPhotonTracker",str))
    if(fName_GreenPhotonTracker != std::string(str))
      { modelChanged=true; fName_GreenPhotonTracker = str; }

  if( r.Get("NoiseMaker",str))
    if(fName_NoiseMaker != std::string(str))
      { modelChanged=true; fName_NoiseMaker = str; }

  if( r.Get("Afterpulser",str))
    if(fName_Afterpulser != std::string(str))
      { modelChanged=true; fName_Afterpulser = str; }

  if(modelChanged) CreateModelObjects();

  fPhotonComputer->    Configure(r);
  fBluePhotonTracker-> Configure(r);
  fWlsFibreModel->     Configure(r);
  fGreenPhotonTracker->Configure(r);
  fNoiseMaker->        Configure(r);
  fAfterpulser->       Configure(r);

  r.Get("ElectronFudgeThreshEnergy",fElectronFudgeThreshEnergy);
  r.Get("ElectronFudgeThreshDist",  fElectronFudgeThreshDist);

  r.Get("DebugMask",fDebugMask);
  r.Get("DebugFileName",str);
  fDebugFileName = str;
  
}

//......................................................................

void PhotonTransport::BeginJob() 
{
//======================================================================
// Handle job module initialization
//======================================================================
  
  if(fDebugFileName.length()>0) {
    if (fDebugFile) delete fDebugFile;
      TDirectory* savedir = gDirectory;
      fDebugFile = new TFile(fDebugFileName.c_str(),"RECREATE");
      if ( fDebugMask & 1 ) fDebugNtuple = new TNtuple("pes","pes",
                "plane:strip:end:thit:tblue:tgreen:tpe:gcosx:gx");
      if ( fDebugMask & 2 ) fBadHitNtuple = new TNtupleD("badhit","badhit",
                "pE:pId:pln:stp:t0:xl0:yl0:zl0:xg0:yg0:zg0:dS:dE:failmode");
      savedir -> cd();
  }
}


void PhotonTransport::CreateModelObjects()
{
  if(fPhotonComputer) delete fPhotonComputer;
  if(fBluePhotonTracker) delete fBluePhotonTracker;
  if(fWlsFibreModel) delete fWlsFibreModel;
  if(fGreenPhotonTracker) delete fGreenPhotonTracker;
  if(fNoiseMaker)         delete fNoiseMaker;
  if(fAfterpulser)        delete fAfterpulser;
  
  fPhotonComputer     = PhotonTransportMaker::Create(fName_PhotonComputer);
  fBluePhotonTracker  = PhotonTransportMaker::Create(fName_BluePhotonTracker);
  fWlsFibreModel      = PhotonTransportMaker::Create(fName_WlsFibreModel);
  fGreenPhotonTracker = PhotonTransportMaker::Create(fName_GreenPhotonTracker);
  fNoiseMaker         = PhotonTransportMaker::Create(fName_NoiseMaker);
  fAfterpulser        = PhotonTransportMaker::Create(fName_Afterpulser);

  if(!fPhotonComputer)     { MSG("Photon",Msg::kFatal) << "FATAL: No such module: " << fName_PhotonComputer << endl; exit(1);};
  if(!fBluePhotonTracker)  { MSG("Photon",Msg::kFatal) << "FATAL: No such module: " << fName_BluePhotonTracker << endl; exit(1);};
  if(!fWlsFibreModel)      { MSG("Photon",Msg::kFatal) << "FATAL: No such module: " << fName_WlsFibreModel << endl; exit(1);};
  if(!fGreenPhotonTracker) { MSG("Photon",Msg::kFatal) << "FATAL: No such module: " << fName_GreenPhotonTracker << endl; exit(1);};
  if(!fNoiseMaker)         { MSG("Photon",Msg::kFatal) << "FATAL: No such module: " << fName_NoiseMaker << endl; exit(1);};
  if(!fAfterpulser)        { MSG("Photon",Msg::kFatal) << "FATAL: No such module: " << fName_Afterpulser << endl; exit(1);};


  fPhotonComputer->    SetRandom(fRandom);
  fBluePhotonTracker-> SetRandom(fRandom);
  fWlsFibreModel->     SetRandom(fRandom);
  fGreenPhotonTracker->SetRandom(fRandom);
  fNoiseMaker->        SetRandom(fRandom);
  fAfterpulser->       SetRandom(fRandom);
}



JobCResult PhotonTransport::ConfirmModelObjects() 
{  
  // Check for success.
  if(!fPhotonComputer) {
    MSG("Photon",Msg::kError) << "Photon Transport Model PhotonComputer = " 
                              << fName_PhotonComputer 
                              << " is not registed. Aborting!"
                              << endl;
    return JobCResult::kFatal;
  }

  if(!fBluePhotonTracker) {
    MSG("Photon",Msg::kError) << "Photon Transport Model BluePhotonTracker = " 
                              << fName_BluePhotonTracker
                              << " is not registed. Aborting!"
                              << endl;
    return JobCResult::kFatal;
  }

  if(!fWlsFibreModel) {
    MSG("Photon",Msg::kError) << "Photon Transport Model WlsFibreModel = " 
                              << fName_WlsFibreModel
                              << " is not registed. Aborting!"
                              << endl;
    return JobCResult::kFatal;
  }
  
  if(!fGreenPhotonTracker) {
    MSG("Photon",Msg::kError) << "Photon Transport Model GreenPhotonTracker = " 
                              << fName_GreenPhotonTracker
                              << " is not registed. Aborting!"
                              << endl;
    return JobCResult::kFatal;
  }

  if(!fNoiseMaker) {
    MSG("Photon",Msg::kError) << "Photon Transport Model NoiseMaker = " 
                              << fName_NoiseMaker
                              << " is not registed. Aborting!"
                              << endl;
    return JobCResult::kFatal;
  }

  if(!fAfterpulser) {
    MSG("Photon",Msg::kError) << "Photon Transport Model Afterpulser = " 
                              << fName_Afterpulser
                              << " is not registed. Aborting!"
                              << endl;
    return JobCResult::kFatal;
  }

  return JobCResult::kAOK;
}


JobCResult PhotonTransport::SimulateEvent()
{
  JobCResult res = JobCResult::kPassed;
  res = ConfirmModelObjects();
  if(res.HaveError()) return res;

  fProtonRangeTable.Reset(fContext);
  fElectronRangeTable.Reset(fContext);

  // Set up the modules.
  fPhotonComputer->    DoReset(fContext);
  fBluePhotonTracker-> DoReset(fContext);
  fWlsFibreModel->     DoReset(fContext);
  fGreenPhotonTracker->DoReset(fContext);
  fNoiseMaker        ->DoReset(fContext);
  fAfterpulser       ->DoReset(fContext);

  TObject* tobj;
  TIter hitarrayIter(fHitList);
  while( (tobj = hitarrayIter.Next()) ) {
    DigiScintHit* scinthit = dynamic_cast<DigiScintHit*>(tobj);
    if(scinthit) {
      JobCResult hitRes = SimulateScintHit(scinthit); // Add up all the errors.
      res |= hitRes;
    }
  }  

  res |= AddNoise();

  res |= AddAfterpulsing();

  if((fEventResult.energyDiscardedGeom+
      fEventResult.energyDiscardedBad) > 1.0*Munits::MeV) 
    MSG("Photon",Msg::kError) 
      << "Run: " << fSimHeader->GetRun()
      << "  Event: " << fSimHeader->GetSnarl() 
      << " Threw out " 
      << fEventResult.energyDiscardedGeom/Munits::MeV << " MeV (bad geom) and "
      << fEventResult.energyDiscardedBad/Munits::MeV << " MeV (bad hits) of a total " 
      << fEventResult.totalHitEnergy/Munits::MeV << " MeV in the event" 
      << endl;
  
  fEventResult.totalStripsHit = fEventResult.fStripsHit.size();
  fEventResult.totalPixels    = fEventResult.fPixelsHit.size();

  return res;
}

enum EBoundaryStatus { kBoundOkay, 
                       kX1OffEnd, 
                       kX2OffEnd, 
                       kInsideBypass,
                       kY1OffWidth,
                       kY2OffWidth,
                       kZ1OffThick,
                       kZ2OffThick };

const char* kBoundaryString [] = { "OK",
                                   "X1 off end",
                                   "X2 off end",
                                   "Inside bypass",
                                   "Y1 off width",
                                   "Y2 off width",
                                   "Z1 off thick",
                                   "Z2 off thick" };

int PhotonTransport::VerifyScintHit(DigiScintHit* hit,UgliStripHandle ustrip)
{
  // Look at the scint hit and see if it's bogus.
  // Return 0 if it's OK
  // Return an integer corresponding to the coordinate violated if not.
  
  if(fabs(hit->X1()) > (ustrip.GetHalfLength()+0.1*Munits::cm)) return kX1OffEnd;
  if(fabs(hit->X2()) > (ustrip.GetHalfLength()+0.1*Munits::cm)) return kX2OffEnd;

  // Length checking: it may be on a two-part strip. Which end is it on?
  if(ustrip.WlsBypass()>0.) {
    float x01 = hit->X1() +  ustrip.GetHalfLength();
    float x02 = hit->X2() +  ustrip.GetHalfLength();
    float p_neg = ustrip.PartialLength(StripEnd::kNegative);
    float p_pos = ustrip.GetHalfLength()*2.0 - ustrip.PartialLength(StripEnd::kPositive);
    bool good = false;
    if((x01<p_neg) && (x02<p_neg)) good = true; // They're both on the neg side
    if((x01>p_pos) && (x02>p_pos)) good = true; // They're both on the pos side
     if(good==false) return kInsideBypass;
  }

  if(fabs(hit->Y1()) > ustrip.GetHalfWidth()) return kY1OffWidth;
  if(fabs(hit->Y2()) > ustrip.GetHalfWidth()) return kY2OffWidth;
  if(fabs(hit->Z1()) > ustrip.GetHalfThickness()) return kZ1OffThick;
  if(fabs(hit->Z2()) > ustrip.GetHalfThickness()) return kZ2OffThick;
  return 0;
}

JobCResult PhotonTransport::SimulateScintHit(DigiScintHit* hit)
{
  // Does all the real work.

  if(!hit) {
    MSG("Photon",Msg::kWarning) << "SimulateScintHit got null ScintHit." << endl;
    return JobCResult::kWarning; // Null hit.. shouldn't happen.
  }
  // Set up ugli.
  UgliGeomHandle ugli = UgliGeomHandle(fContext);

  // Get the goods on this strip.
  UgliStripHandle ustrip = ugli.GetStripHandle(hit->StripEndId());

  // Check to see if this strip is valid.
  // This error crops up from PTSim where it generates veto shield hits, but Ugli can't cope
  if(! ustrip.IsValid()) {
    MAXMSG("Photon",Msg::kWarning,10) << "Throwing out hit with no Ugli info: " << hit->StripEndId().AsString() << std::endl;
    hit -> SetFailBit(DigiScintHit::kInvalidStrip);
    return JobCResult::kWarning;
  }

  // Check for some validity things on the hit. Fix them up if neccessary.
  if(hit->DE()<0.) {
    MAXMSG("Photon",Msg::kWarning,5) 
      << "DigiScintHit has negative energy (" 
      << hit->DE()/Munits::GeV 
      << " GeV). Throwing it out." << endl;    
    hit -> SetFailBit(DigiScintHit::kNegEnergy);
    return JobCResult::kWarning;
  }

  if(hit->DE()==0.) {
    MAXMSG("Photon",Msg::kWarning,5) 
      << "DigiScintHit has zero energy. Throwing it out." << endl;
    hit -> SetFailBit(DigiScintHit::kNullEnergy);
    return JobCResult::kWarning;
  }

  fEventResult.totalHits +=1;
  fEventResult.totalHitEnergy += hit->DE();
  fEventResult.fStripsHit.insert(hit->StripEndId().Build18BitPlnStripKey());

  DigiScintHit* computerHit = hit; // A pointer used only for the Computer.

  // Verify that this hit is on a strip that can eventually read out somewhere.
  // NB. Up 0-3 and Vp 64-69 in the ND have scintillator but no readout.
  PlexHandle plex(fContext);
  PlexStripEndId seid1(hit->StripEndId()); seid1.SetEnd(StripEnd::kPositive);
  PlexStripEndId seid2(hit->StripEndId()); seid2.SetEnd(StripEnd::kNegative);
  int numEndsReadOut = 0;
  if(plex.GetPixelSpotId(seid1).IsValid()) numEndsReadOut++;
  if(plex.GetPixelSpotId(seid2).IsValid()) numEndsReadOut++;
  if(numEndsReadOut==0) return JobCResult::kPassed;  // Everything is fine, but don't bother with this hit

  // Verify the geometry in this scint hit makes sense.
  int badBoundary = VerifyScintHit(hit,ustrip);

  if(badBoundary>0) {
    MAXMSG("Photon",Msg::kWarning,20) 
      << "DigiScintHit geometry inconsistency: is not contained in strip " << endl 
      << "   Boundary violation: " << badBoundary 
      << "  " << kBoundaryString[badBoundary] << endl
      << "   Strip: " <<hit->StripEndId().AsString() 
      << " dE = " << hit->DE()/Munits::MeV << " MeV  "
      << " dx = " << hit->DS()/Munits::cm << " cm  " << endl
      << "   ScintHit coords: " 
      << Form("(%.2f,%.2f,%.2f) cm ",hit->X1()/Munits::cm,hit->Y1()/Munits::cm,hit->Z1()/Munits::cm)
      << "to" 
      << Form("(%.2f,%.2f,%.2f) cm",hit->X2()/Munits::cm,hit->Y2()/Munits::cm,hit->Z2()/Munits::cm)
      << endl
      << "   Ugli Geometry: " 
      << Form("(+/-%.2f, +/-%.2f,+/-%.2f) cm",ustrip.GetHalfLength()/Munits::cm,ustrip.GetHalfWidth()/Munits::cm,ustrip.GetHalfThickness()/Munits::cm)
      << ".. skipping it." << endl;
    if ( fBadHitNtuple ) {
      Int_t failmode = 1; // => failed geometry
      FillBadHitNtuple(hit,failmode);
    }
    hit -> SetFailBit(DigiScintHit::kBadGeom);
    fEventResult.hitsDiscardedGeom +=1;
    fEventResult.energyDiscardedGeom += hit->DE();
    return JobCResult::kWarning;
  }     

  if(hit->DS()<=0.) {
    // Zero range! 
    // Try to do a fix.

    if(hit->ParticleId()==2212) {
      // It's a proton. Attempt to do make our own range determination.
      double range = fProtonRangeTable.Interpolate(hit->ParticleKineticEnergy()/Munits::MeV);
      MinosMaterial::StdMaterial_t fStdMatSc=MinosMaterial::ePolystyreneMinos;
      double fRoSc=MinosMaterial::Density(fStdMatSc);//Was 1.06;

      double dx = range / fRoSc * Munits::cm; 
      
      MSG("Photon",Msg::kDebug) << "dx==0 for a proton. Table lookup gives range of " 
                                << dx/Munits::mm << " mm for KE of " 
                                << hit->ParticleKineticEnergy()/Munits::MeV << " MeV" << endl;
      
      if(dx<=0.0) {
        MSG("Photon",Msg::kDebug) << "dx<=0 for a proton, after table lookup!";
        if ( fBadHitNtuple ) {
          Int_t failmode = 3; // => proton dS is zero even after table lookup 
          FillBadHitNtuple(hit,failmode);
        }
        hit -> SetFailBit(DigiScintHit::kNullPath);
        fEventResult.hitsDiscardedBad+=1;
        fEventResult.energyDiscardedBad+=hit->DE();
        return JobCResult::kWarning;
      }

      // Make a new, temporary ScintHit that the Computer can use
      // that has an adjusted track length.
      static DigiScintHit* tempHit = 0;
      if(tempHit) delete tempHit;
      tempHit = new DigiScintHit(hit->TrackId(),hit->ParticleId(),
                                 hit->ParticleEnergy(), hit->StripEndId(),
                                 hit->T1(), hit->X1(), hit->Y1(), hit->Z1(),
                                 hit->T2(), hit->X2(), hit->Y2(), hit->Z2(),
                                 dx, hit->DE());
      computerHit = tempHit;

    } else {

      MAXMSG("Photon",Msg::kWarning,5) 
        << "DigiScintHit has zero path length. Throwing it out." << endl;

      if ( fBadHitNtuple ) {
        Int_t failmode = 2; // => zero pathlength and not proton
        FillBadHitNtuple(hit,failmode);
      }

      hit -> SetFailBit(DigiScintHit::kNullPath);
      fEventResult.hitsDiscardedBad+=1;
      fEventResult.energyDiscardedBad+=hit->DE();
      return JobCResult::kAOK;
    }
  }

  // Hack for low-energy electrons.
  if(abs(hit->ParticleId())==11) {
    // Is it below threshold?
    if((hit->ParticleKineticEnergy() < fElectronFudgeThreshEnergy) ||
       (hit->DS()                    < fElectronFudgeThreshDist) ) {
      // We're in the fudge regime. Do it.

      // Find the start/stop energy.
      double E1 = hit->ParticleKineticEnergy();
      double E2 = E1 - hit->DE();
      if(E2<0) E2 = 0;
      double dE = E1-E2;

      // Compute the range from the tables.      
      double Range1 = fElectronRangeTable.Interpolate(E1);
      double Range2 = fElectronRangeTable.Interpolate(E2);
      double dX = Range1-Range2;

      // Make a new, temporary ScintHit that the Computer can use
      // that has an adjusted track length.
      static DigiScintHit* tempHit = 0;
      if(tempHit) delete tempHit;
      tempHit = new DigiScintHit(hit->TrackId(),hit->ParticleId(),
                                 hit->ParticleEnergy(), hit->StripEndId(),
                                 hit->T1(), hit->X1(), hit->Y1(), hit->Z1(),
                                 hit->T2(), hit->X2(), hit->Y2(), hit->Z2(),
                                 dX, dE);
      computerHit = tempHit;
    }
  }

  // First, ask the photon computer for how many hits we're going to demand.
  PhotonCount demand;
  
  // Set a 'null' demand value.
  demand.SetBlueTotal(0);
  

  // Ask the computer.
  // Note that we give it what each of the modules says is it's 
  // own prescale factor. For instance, the green tracker will usually
  // be prescaled by 0.13, or 13% for the photocathode quantum efficiency.
  // This indicates to the computer that it only needs to generate 
  // 13% of the photons it would otherwise need.  The GreenTracker, in 
  // turn, scales up the quantum efficiency by a factor of 1/0.13.
  // This is used as a trick to speed up the simulation by several factors.
  Bool_t hitOK = fPhotonComputer->
    ComputePhotons(computerHit,
                   fBluePhotonTracker->GetBluePrescaleFactor(),
                   fWlsFibreModel->GetWlsPrescaleFactor(),
                   fGreenPhotonTracker->GetGreenPrescaleFactor(),
                   demand);

  if(!hitOK) {
    if ( fBadHitNtuple ) {
      Int_t failmode = 4;  // failed by PhotonComputer
      FillBadHitNtuple(hit,failmode);
    }
    hit -> SetFailBit(DigiScintHit::kPhotonComputer);
    fEventResult.hitsDiscardedBad+=1;
    fEventResult.energyDiscardedBad+=hit->DE();
    return JobCResult::kAOK;
  }
  // The photon computer will have asked us for one of:
  // 1) A total number of blue photons
  // 2) A number of blue photons to be made for the +ve end readout
  //   and -ve end readout
  // 3) A total number of green photons to be made
  // 4) A number of green photons to be made (+,-)
  // 5) A total number of photoelectrons (usually used for mirrored strips)
  // 6) A number of final photoelectrons to be mad (+,-)

  // Our job now is to keep simulating photons until we reach 
  // the desired number, whichever it is.
  // The unused demand numbers are automatically set to -1, so we can simply 
  // start and zero and continue until something matches.

  Int_t nBlue_Total=0;
  Int_t nBlue_Pos=0;
  Int_t nBlue_Neg=0;
  Int_t nGreen_Total=0;
  Int_t nGreen_Pos=0;
  Int_t nGreen_Neg=0;
  Int_t nPe_Total = 0;
  Int_t nPe_Pos=0;
  Int_t nPe_Neg=0;
  
  DigiPhoton* bluePhoton =0;
  DigiPhoton* greenPhoton =0;
  DigiPE*     pe =0;
  
  bool done = false;
  while(!done) {

    // --------------------------------------
    // Clean up from last pass.
    if(bluePhoton) delete bluePhoton;
    if(greenPhoton) delete greenPhoton;
    if(pe) delete pe;
    bluePhoton = 0;
    greenPhoton = 0;
    pe = 0;

    
    // --------------------------------------
    // Are we done?

    // Deal with finish condition (1) above (the simplest):
    // if we have made enough blue photons, this is the last one.
    if(nBlue_Total == demand.GetBlue_Total()) done = true; // We reached demand.
        
    // Finish condition (2):
    if( (nBlue_Pos == demand.GetBlue_Pos()) && 
        (nBlue_Neg == demand.GetBlue_Neg()) ) done = true; // We reached demand.

    // Deal with finish condition (3) above
    // if we have made enough green photons.
    if(demand.GetGreen_Total() == nGreen_Total) done = true; // We reached demand.
    
    // Finish condition (4):
    if( (nGreen_Pos == demand.GetGreen_Pos()) && 
        (nGreen_Neg == demand.GetGreen_Neg()) ) done = true; // We reached demand.
    
    // Finish condition (5): we got enough PE.
    if(nPe_Total == demand.GetPe_Total()) done = true;
    
    // Condition (6)
    if( (nPe_Pos == demand.GetPe_Pos()) && 
        (nPe_Neg == demand.GetPe_Neg()) ) done = true; // We reached demand.
    
    if(done) break;


    // ----------------------------------------------------------
    // Choose which direction this photon will go.

    // By default, we don't know which end this thing should go to.
    // i.e. cases (1) and (3) above.
    StripEnd::StripEnd_t direction = StripEnd::kUnknown; 

    // Deal with cases (2),(4),and(5) to choose direction:
    if(demand.IsDirectional()) {
      // If we don't have enough positives, go that way.
      // If we do, go the other way.
      if( (nBlue_Pos  < demand.GetBlue_Pos()) ||
          (nGreen_Pos < demand.GetGreen_Pos()) ||
          (nPe_Pos    < demand.GetPe_Pos()) ) {
        direction = StripEnd::kPositive;
      } else {
        direction = StripEnd::kNegative;
      }
    }

    //---------------------------------------
    // Step I::Simulate scintillator photons
    Bool_t blue_hit = fBluePhotonTracker->ScintPhotonToFibreHit(hit,
                                                                ustrip,
                                                                bluePhoton);

    if(!bluePhoton) continue; // Didn't actually make a blue photon.. tracking error or somesuch.
    
    nBlue_Total++; // We made a blue photon.
    if(direction==StripEnd::kPositive) nBlue_Pos++; // Blue +ve going photon
    if(direction==StripEnd::kNegative) nBlue_Neg++; // Blue -ve going photon

    // Don't bother with the rest of the steps if it didn't even hit the fibre.
    if(!blue_hit) continue;

    // ------------------------------------------------
    // Step II:Convert from blue photon to a fibre hit.
    Bool_t fibre_hit = 
      fWlsFibreModel->FibreHitToGreenPhoton(ustrip,
                                            direction,
                                            bluePhoton,
                                            greenPhoton);

    if(!fibre_hit) continue; // Didn't successfully make a green photon..
    
    nGreen_Total++;
    if(direction==StripEnd::kPositive) nGreen_Pos++; // Green +ve going photon
    if(direction==StripEnd::kNegative) nGreen_Neg++; // Green -ve going photon
    
    
    // ------------------------------------------------
    // Step III: Track the green photon to the phototube.

    Bool_t greenSurvived = 
      fGreenPhotonTracker->GreenPhotonToPe( ustrip,
                                            greenPhoton,
                                            pe,
                                            direction);
    
    
    if((!greenSurvived)||(!pe)) continue; // green photon didn't make it to the photocathode.

    nPe_Total++;
    if(direction==StripEnd::kPositive) nPe_Pos++; // Blue +ve going photon
    if(direction==StripEnd::kNegative) nPe_Neg++; // Blue -ve going photon

    
    
    // Finally! We have made a photoelectron. Tack it onto the list...
    // if it's a valid pixel.
    if(pe->GetPixelSpotId().IsValid()) {
      fEventResult.totalPE += 1;
      fEventResult.fPixelsHit.insert(pe->GetPixelSpotId().GetUniquePixelEncodedValue());
      fPeList->Add(pe);

      if(fDebugNtuple)
        fDebugNtuple->Fill(hit->Plane(),hit->Strip(),(float)direction,
                           hit->T1(), bluePhoton->T(), 
                           greenPhoton->T(), pe->GetTime(),
                           greenPhoton->GetCosX(),
                           greenPhoton->X()
                           );

      pe = 0; // So it doesn't get deleted.
    }
  };

  MSG("Photon",Msg::kDebug) << "Hit: " 
                              << hit->DE()/Munits::MeV << " MeV -> " 
                              << nBlue_Total << " blue -> "
                              << nGreen_Total << " green -> "
                              << nPe_Total << " pe"
                              << endl;

  MSG("Photon",Msg::kVerbose) << "Hit result:" << endl 
                            << "Blue photons:    Total ( " << nBlue_Total << " / " << demand.GetBlue_Total() << " ) "
                            << " +ve ( "<< nBlue_Pos << " / " << demand.GetBlue_Pos() << " ) "
                            << " -ve ( "<< nBlue_Neg << " / " << demand.GetBlue_Neg() << " ) " << endl
                            << "Green photons:   Total ( " << nGreen_Total << " / " << demand.GetGreen_Total() << " ) "
                            << " +ve ( "<< nGreen_Pos << " / " << demand.GetGreen_Pos() << " ) "
                            << " -ve ( "<< nGreen_Neg << " / " << demand.GetGreen_Neg() << " ) " << endl
                            << "Photoelectrons:  Total ( " << nPe_Total << " / " << demand.GetPe_Total() << " ) "
                            << " +ve ( "<< nPe_Pos << " / " << demand.GetPe_Pos() << " ) "
                            << " -ve ( "<< nPe_Neg << " / " << demand.GetPe_Neg() << " ) " << endl;
    
  fEventResult.bluePhotons+=nBlue_Total;
  fEventResult.greenPhotons+=nGreen_Total;
  
  fEventResult.bluePhotons_nonprescaled += (Int_t)(nBlue_Total/(fBluePhotonTracker->GetBluePrescaleFactor()
                                                          *fWlsFibreModel->GetWlsPrescaleFactor()
                                                          *fGreenPhotonTracker->GetGreenPrescaleFactor() ) );
  fEventResult.greenPhotons_nonprescaled += (Int_t)(nGreen_Total/fGreenPhotonTracker->GetGreenPrescaleFactor());


  // Cleanup.
  if(bluePhoton) delete bluePhoton;
  if(greenPhoton) delete greenPhoton;
  if(pe) delete pe;
 
  return JobCResult::kPassed;
}

JobCResult PhotonTransport::AddNoise()
{

  // Find the times of the hits to get an idea what the window is.
  const double kBig = 1e99;
  Double_t tstart = kBig;
  Double_t tend   = -kBig;
  TObject* tobj;
  TIter hitarrayIter(fHitList);
  while( (tobj = hitarrayIter.Next()) ) {
    DigiScintHit* scinthit = dynamic_cast<DigiScintHit*>(tobj);
    if(scinthit) {
      if(scinthit->T1()<tstart) tstart = scinthit->T1();
      if(scinthit->T2()>tend)   tend   = scinthit->T2();
    }
  }  

  // Couldn't find a reasonable time to put the noise in:
  if(tstart== kBig) return JobCResult::kFailed;
  if(tend  ==-kBig) return JobCResult::kFailed;

  std::vector<DigiPE*> peNoiseList;

  MSG("Photon",Msg::kDebug) << "Calling NoiseMaker  t=(" << tstart << ", " << tend << ")\n";
  fNoiseMaker->MakeNoise(peNoiseList,tstart,tend);

  for(UInt_t i=0;i<peNoiseList.size();i++) {
    fPeList->Add(peNoiseList[i]);
  }
  fEventResult.noisePE += peNoiseList.size();
  fEventResult.totalPE += peNoiseList.size();


  return JobCResult::kPassed;
}


JobCResult PhotonTransport::AddAfterpulsing()
{

  // Find the times of the hits to get an idea what the window is.
  std::vector<DigiPE*> apList;
  std::vector<DigiPE*> peList;

  MSG("Photon",Msg::kDebug) << "Calling Afterpulser." << std::endl;

  peList.resize(fPeList->GetEntries(),0);
  for(UInt_t i=0; i<peList.size(); i++) {
    peList[i] = dynamic_cast<DigiPE*>((*fPeList)[i]);
  }

  fAfterpulser->MakeAfterpulses(peList,apList);

  for(UInt_t i=0;i<apList.size();i++) {
    fPeList->Add(apList[i]);
  }
  fEventResult.afterpulsePE += apList.size();
  fEventResult.totalPE      += apList.size();


  return JobCResult::kPassed;
}


void PhotonTransport::Print(Option_t* ) const
{
  MSG("Photon",kInfo) << "PhotonTransport Configuration:" << endl;
  if(fPhotonComputer) fPhotonComputer->Print("c");
  if(fBluePhotonTracker) fBluePhotonTracker->Print("b");
  if(fWlsFibreModel) fWlsFibreModel->Print("w");
  if(fGreenPhotonTracker) fWlsFibreModel->Print("g");
  if(fNoiseMaker) fNoiseMaker->Print("n");
  if(fAfterpulser) fAfterpulser->Print("n");
}


void PhotonTransport::FillBadHitNtuple(const DigiScintHit* hit, 
                                       Int_t failmode) {
  
  // Convert local to global coordinates
  UgliGeomHandle ugli = UgliGeomHandle(fContext);
  const UgliStripHandle& uglistp = ugli.GetStripHandle(hit->StripEndId());
  TVector3 vlocal0(hit->X1(),hit->Y1(),hit->Z1());
  const TVector3& vglobal0 = uglistp.LocalToGlobal(vlocal0);

  fBadHitNtuple -> Fill(hit->ParticleEnergy(),
                        hit->ParticleId(),hit->Plane(),hit->Strip(),
                        hit->T1(),hit->X1(),hit->Y1(),hit->Z1(),
                        (Float_t)(vglobal0.X()),(Float_t)vglobal0.Y(),
                        (Float_t)vglobal0.Z(),
                        hit->DS(),hit->DE(),failmode);

  return;
}

---