PROGRAM PoliMiPP
! Version 6.21
! Written January 7th, 2009
! By Brandon Grogan
! At The Oak Ridge National Laboratory
! Last Edited: January 8th, 2013

!==========================================================================================
! Variable declaration
!==========================================================================================

IMPLICIT NONE

! 1=History number, 2=particle number, 3=particle type, 4=interaction type, 5=target nucleus
! 6=cell, 7=Edep, 8=time(shakes), 9-11=x,y,z, 12=weight, 13=gen., 14=Num scat., 15=code, 16=Initial Energy
TYPE DatFile
  INTEGER(KIND=4) Col1, Col2, Col3, Col4, Col5, Col6
  REAL(KIND=4) Col7, Col8, Col9, Col10, Col11, Col12
  INTEGER(KIND=4) Col13, Col14, Col15
  ! REAL(KIND=4) Col16
END TYPE DatFile

TYPE DetectorPulse
  INTEGER(KIND=4) HistNo, DetNo, ParNo, Direct, SubHist, XTalk, HGamma, GenNo
  REAL(KIND=4) :: Time, PulseHeight
END TYPE DetectorPulse

TYPE PeaksFile
  INTEGER(Kind=4) DetNo, TotalPeak
  REAL(KIND=4) TotalMean
  INTEGER(KIND=4) DirectPeak
  REAL(KIND=4) DirectMean
  INTEGER(KIND=4) XTPeak
  REAL(KIND=4) XTMean
END TYPE PeaksFile

TYPE(DatFile), ALLOCATABLE :: ThisHist(:), DatFileInput(:)
TYPE(DatFile) :: HistSwap, TestInput

TYPE(DetectorPulse), ALLOCATABLE :: Pulses(:)
INTEGER(KIND=4) :: TotalPulses

! This array will temporarily store the data on the starting histories until the total
! number of histories in the .DAT file is known, at which time the data will be moved
! to the (smaller) DatHist array.
INTEGER(KIND=4), ALLOCATABLE :: TempHist(:,:)

! The DatHist array contains the history the history number, starting line number (in the
! DatFile array) of each history number, and the number of events for that history.
INTEGER(KIND=4), ALLOCATABLE :: DatHist(:,:), PulseHist(:,:), Detector(:)

! Records the number of detectors (NumDet) and detector cell numbers found in the .DAT file.
! CurrentDet holds the cell number of the detector in the current history that is being 
! manipulated by the program.
INTEGER(KIND=4) :: NumDet, CurrentDet, MaxDets
! NewDet indicates whether a detector cell number has been seen in the .DAT file previously.
LOGICAL :: StopAfterLoop, EndOfCycle, FirstCycle

! Integer variables used for DO loops and swapping
INTEGER(KIND=4) :: I, J, K, L, M, SubHistNo, NumSubHist

! TotalLines is the number of rows in the .DAT file.  TotalHist is the number of histories
! in the .DAT file.  MaxEvt is the maximum number of events recorded for a single history.
! MaxEvtHist is the history number with the maximum number of events.
INTEGER(KIND=4) :: TotalLines, TotalHist, MaxEvt, MaxEvtHist, NumHistP, MaxPulse, MaxPulseHist

! 'History' holds the history number of the line currently being read in the .DAT file.
! 'OldHistory' holds the history number of the previous line
INTEGER(KIND=4) :: History, OldHistory

! 'InputStatus' holds the IO status of a file read.  IO status < 0 indicates end of file.
INTEGER :: InputStatus, FileBaseSize
CHARACTER :: DatName*64, FileBase*64, DetectName*4, NPSText*11, OneLetter*1
CHARACTER :: CorrWindowText*11, DeadTimeText*11, DatFormatText*128, nThreshText*11
CHARACTER :: OutFile*64, ThisLine*80, ParamFile*80

! These variables store the time at the beginning and end of the program in order to 
! measure the execution time of the program.
REAL(Kind = 4) :: StartTime, FinishTime

! Dummy Variables
INTEGER(KIND=4) :: DummyInt
REAL(KIND=4) :: DummyReal
CHARACTER :: DummyChar*1

! DeadStart records the time that the detector deadtime window opens.  PulseStart records
! the time the pulse generation window opens.  PHThisEvent is the total light output
! attributed to a particular event in the .DAT file.  PHTotal is the total light output
! generated during the pulse generation window so far.
REAL(KIND = 4) :: DeadStart, PulseStart, PHThisEvent

REAL(KIND=4), ALLOCATABLE :: PHTotal(:)

! Input parameters.  These are currently hard-coded, but may be placed off-line later.
! Pulse generation time.  This is the amount of time it takes the detector to generate a 
! pulse once the initial energy is deposited in the detector cell.  If two events occur
! within pgentime of each other, their light output will combine.
REAL(KIND=4) :: pgentime
! The time after a pulse in a detector in which all further events are lost.
REAL(KIND=4) :: deadtime
! The light output (in MeVee) required to produce a pulse in the detector.
REAL(KIND=4) :: threshold, nthresh, sdthresh, samplethresh, deltathresh

! The minimum and maximum times over which pulses will be accepted.
REAL(KIND=4) :: CorrMax, CorrMin, TimeIncrement, WindMin, WindMax
INTEGER(KIND=4) :: NumTBins
! The cell number of the start detector which will be used for cross-correlation
! calculations
INTEGER(KIND=4) :: StartDet, StartDetRow, StopDetRow

! Used to hold the integer (rather than real) pulse time.
INTEGER(KIND=4) :: IntTime, TimeLag, CorrWindowOverflow, CorrWindowUnderflow

! Arrays used to hold the program output.
INTEGER(KIND=4), ALLOCATABLE :: Correlation(:,:,:), CrossCorr(:,:,:,:), EventScatter(:,:,:)

! Variables used for calculating multiplicities
INTEGER(KIND=4) :: nps
INTEGER(KIND=4), ALLOCATABLE :: Multiplicity(:,:)

TYPE(PeaksFile), ALLOCATABLE :: Peaks(:)
INTEGER(KIND=4) :: PeakWidth
INTEGER(KIND=4) :: PeakSumTotal, PeakSumDirect, PeakSumNoXT, PeakStart, PeakOld

! PH Spectrum parameters
REAL(KIND=4) :: PHSIncrement, PHSMax
INTEGER(KIND=4) :: PHSNumBins, PHSBin
INTEGER(KIND=4), ALLOCATABLE :: PHSpectrum(:,:,:)

! Variables Used for assigning the minimum particle number and subhist to a pulse if
! it is the aggregation of more than 1 event
INTEGER(KIND=4) :: MinParNo, MinSubHist, MinCollisions, MinGen, MinCode
REAL(KIND=4) :: MinWeight
REAL(KIND=4) :: NPSnorm

! These values are used to control the maximum number of lines and histories the program
! will read from the .DAT files.  If these values are too small, the entire .DAT file
! will not be processed.  If they are set too large, a huge .DAT file may overflow 
! system memory and cause a crash.  Values as large as MaxHist = 10000000 and 
! MaxLines = 24000000 have been tested.
INTEGER(KIND=4) :: MaxHist, MaxLines, LinesPerCycle, LinesThisCycle, HistThisCycle
INTEGER(KIND=4) :: PulsesThisCycle

! Added with version 4.0.  Timing parameters to take add timing uncertainty to the 
! source and detector signals.  Updated in version 5.0 to add a fixed amount of
! time to each history from a uniform distribution.  This will simulate a source
! with a finite pulse time width, e.g., MP-320 pulsed DT generator.
REAL(KIND=4) :: SrcFWHM, DetFWHM, SrcWidth, SrcStart, SrcDeltaT, DetDeltaT, PulseTime
! Added in version 5.6. These two time constants define the shape of the detector time
! response. Instead of sampling this uncertainty from a balanced (Gaussian) as was done
! in previous versions, this allows for the sampling of the detector time from a more
! realistic exponential sum with a rapid rise and a slower decay, giving it a long tail.
! The decay constant is identified as short to allow for the potential of a long decay
! constant (for liquid scintillators) to be added in the future.
REAL(KIND=4) :: DetRise, DetDecayShort

! These variables are used to simulate the energy resolution for gamma-ray spectroscopy.
! PEperMeV is the number of photoelectrons produced at the photocathode of the PMT.
! Assuming a 20% PMT efficiency, a typical values is 2300 for NaI.  For EJ-200 (plastic)
! scintillators the value is 1450 and for EJ-301 (liquid) scintillators the value is 1200.
! A comparable value for HPGe detectors is 2.55x10^6. 
REAL(KIND=4) :: PEperMeV, NumPE, NumPEFWHM, DeltaPE, HLineMin, HLineMax, E_Gamma
INTEGER(KIND=4) :: TotalHGammas
INTEGER(KIND=4) :: MinDetNum, MaxDetNum

! These are the coefficients of the neutron light curve.
REAL(KIND=4) :: Coeff0, Coeff1, Coeff2

! This coefficient is the empirical normalization used to match the simulated pixel
! profiles with measured values. SrcNorm will modify the nps value to nps = nps / SrcNorm.
! The default value is set to 1 (no change to nps).
REAL(KIND=4) :: SrcNorm

! The ".hits" file is an output with information about each pulse. HitsFile = 1 enables writing of the file.
INTEGER(KIND=4) :: HitsFile

110 FORMAT(I11,I5,I3,I5,I6,I4,F10.5,F10.3,F9.2,F8.2,F8.2,F7.3,I5,I6,I4)
130 FORMAT(I11, I4, I5, I4, I5, I3, I4, I4, F15.3, F12.6)

TotalHGammas = 0
MaxDets = 50

! These are default values for the user-definable parameters. 
! These values will only be used if they are not specified in the parameter file.
StartDet = 0
NPS = 0
SrcNorm = 1.0
CorrMin = 0.0
CorrMax = 256.0
NumTBins = 256
WindMin = 0.0
WindMax = 256.0
DeadTime = 50.0
nthresh = 1.0
sdthresh = 0.0
Coeff0 = 0.0
Coeff1 = 0.125
Coeff2 = 0.0364
pgentime = 1.0
PeakWidth = 5
PHSMax = 6.00
PHSNumBins = 4096
SrcWidth = 0.0
SrcFWHM = 1.0
DetFWHM = 1.0
DetRise = 0.0
DetDecayShort = 0.0
PEperMEV = 900
MinDetNum = 0
MaxDetNum = 99999
E_Gamma = 2.223250
MaxHist = 2000000000
MaxLines = 2000000000
HitsFile = 0
! End of user-definable defaults


LinesPerCycle = 5000000


!==========================================================================================
! Step 0: Read variables from command line, open .DAT file, initialize variables
!==========================================================================================

CALL CPU_TIME(StartTime)

WRITE(*,*) "PolimiPP, Version 6.21"
WRITE(*,*) "Modified 8 January 2013"

! Reads the parameter filename from the command line
CALL GetArg(1, ParamFile)


IF (LEN(TRIM(ParamFile)) == 0) THEN
  WRITE(*,*) "Enter the name of the PP Paramter File"
  READ (*,*) ParamFile
  CALL ReadParams()
  WRITE(*,*) "Opening file ", TRIM(DatName)
ELSEIF(TRIM(Paramfile) == "default" .OR. TRIM(Paramfile) == "DEFAULT") THEN
  WRITE(*,*) "Using default post-processor values."
  CALL GetArg(2, DatName)
  IF (LEN(TRIM(DatName)) == 0) THEN
    WRITE(*,*) "Enter the name of the .DAT File"
    READ (*,*) DatName
  END IF
ELSE
  CALL ReadParams()
  WRITE(*,*) "Opening file ", TRIM(DatName)
END IF

IF (StartDet == 0) THEN
  WRITE(*,*) "Starting Detector Cell set to ", StartDet
  WRITE(*,*) "Detector-detector cross-correlations will not be computed for this file."
  WRITE(*,*) " "
ELSEIF (StartDet < 0) THEN
  WRITE(*,*) "Starting Detector value of < 0.  Cross-correlations will be computed for ALL "
  WRITE(*,*) "combinations of start and stop detectors."
  WRITE(*,*) " "
ELSE
  WRITE(*,*) "Cross-correlations will be computed using Detector Cell number ", StartDet, &
  & " as the starting detector."
  WRITE(*,*) " "
END IF


! Calculate the width of each time bin that will be used for recording correlations.
TimeIncrement = CorrMax / INT(NumTBins)
! Because the pulse times will be rounded into the nearest bin using integers, the 
! minimum and maximum time bins need to be be expanded by 1/2 of the bin width.
! Otherwise, they will be only half the width of the other bins.
CorrMin = CorrMin - 0.5 * TimeIncrement
CorrMax = CorrMax + 0.5 * TimeIncrement

! Calculate the width of the energy bins which will be used for recording the 
! pulse height spectrum information.
PHSIncrement = PHSMax / REAL(PHSNumBins)

! These 2 variables will be used to determine if a given pulse falls within the FWHM
! of the hydrogen capture gamma line.  The factor 2.35482/2 will results in an
! acceptance of +/- 1.18 standard deviations, or 74%, of the true capture gammas.
! This width can be adjusted to larger or smaller values as desired.
HLineMin = (E_Gamma * PEperMeV - 2.35482/2.*SQRT(E_Gamma * PEperMeV)) / PEperMeV
HLineMax = (E_Gamma * PEperMeV + 2.35482/2.*SQRT(E_Gamma * PEperMeV)) / PEperMeV
! WRITE(*,*) HLineMin, HLineMax

FileBase = " "

J = ICHAR(".")
DO I = 1,LEN(TRIM(DatName))
  IF (DatName(I:I) == ".") THEN
    EXIT
  ELSE 
    OneLetter = DatName(I:I)
    FileBase(I:I) = OneLetter
  END IF
END DO

! Opens .DAT file
OPEN(UNIT=1, FILE=TRIM(DatName), ACTION="READ", STATUS="OLD", POSITION="REWIND", IOSTAT = InputStatus)
IF (InputStatus > 0) STOP "Error opening specified .DAT file"



ALLOCATE(DatFileInput(LinesPerCycle+50000))
ALLOCATE(DatHist(LinesPerCycle,3)) ! This is extremely conservative - allows each history to have exactly one line.
ALLOCATE(Detector(MaxDets))




!==========================================================================================
! Step 1 : Process the raw .DAT file. The program will read in a certain number of lines
! at a time (a cycle) then process the data it has read so far. After the data has been 
! processed, the program will read more until the end of the file is reached or one of the
! other stop conditions (MaxHist, MaxLines) has been reached.
!==========================================================================================
StopAfterLoop = .FALSE.
EndOfCycle = .FALSE.
FirstCycle = .TRUE.
LinesThisCycle = 0
TotalLines = 0
HistThisCycle = 0
TotalHist = 0
OldHistory = 0
NumDet = 0
Detector = 0
MaxEvt = 0
! This loop runs through the .DAT file to record the number of events and histories.
! The loop will terminate when the end of the file is reached.
DO M=1,MaxLines
  LinesThisCycle = LinesThisCycle + 1
  TotalLines = TotalLines + 1

  ! Reads the history number of each line
  READ(1,*, IOSTAT=InputStatus) DatFileInput(LinesThisCycle)
  ! Stops the program if the .DAT file cannot be read.
  IF (InputStatus > 0) STOP "*** Error Reading DAT file ***  Unrecognized .DAT file format"


!  IF(M <= 10) WRITE(*,*) History, CurrentDet

  IF(InputStatus == 0) THEN
    History = DatFileInput(LinesThisCycle)%Col1
    CurrentDet = DatFileInput(LinesThisCycle)%Col6
    CALL AssignAndSortDetectorNumbers(CurrentDet)
  END IF

  ! The loop will stop reading data and process data given one of these 4 conditions is met.
  IF(InputStatus < 0 .OR. TotalLines >= MaxLines .OR. TotalHist >= MaxHist .OR. MOD(TotalLines, LinesPerCycle)==0) THEN
   IF (MOD(TotalLines, LinesPerCycle)==0) WRITE(*,*) TotalLines, " lines read"

   EndOfCycle = .TRUE. ! Stop reading data and process it at the end of this iteration of the 'M' loop.
   IF(InputStatus == 0) THEN 
     ! We first have to record the last line read and determine if it was a new history.
     IF (History .NE. OldHistory) THEN ! Start of a new history
       TotalHist = TotalHist + 1
       HistThisCycle = HistThisCycle + 1
       DatHist(HistThisCycle,1) = History
       DatHist(HistThisCycle,2) = LinesThisCycle
       DatHist(HistThisCycle,3) = 1
       OldHistory = History
     ELSE ! This is another event in a history that has already been started
       DatHist(HistThisCycle,3) = DatHist(HistThisCycle,3) + 1
       IF (DatHist(HistThisCycle,3) > MaxEvt) THEN
        MaxEvt = DatHist(HistThisCycle,3) ! Record the maximum number of events (lines) for any one history
        MaxEvtHist = History ! The history # with the most events (lines.)
       END IF
     END IF

     ! Now read ahead until reaching the end of this history.
     DO J = 1, 50000 
      LinesThisCycle = LinesThisCycle + 1
      TotalLines = TotalLines + 1
      IF (MOD(TotalLines, LinesPerCycle)==0) WRITE(*,*) TotalLines, " lines read"
       READ(1,*, IOSTAT=InputStatus) DatFileInput(LinesThisCycle)
       IF (InputStatus .NE. 0) THEN ! End of file was reached while looking for the end of the history.
      !   LinesThisCycle = LinesThisCycle + J
       !  TotalLines = TotalLines + J
         EXIT
       ELSE IF (DatFileInput(LinesThisCycle)%Col1 == History ) THEN ! New line is the same history
        CALL AssignAndSortDetectorNumbers(DatFileInput(LinesThisCycle)%Col6)
        DatHist(HistThisCycle,3) = DatHist(HistThisCycle,3) + 1
        IF (DatHist(HistThisCycle,3) > MaxEvt) THEN
         MaxEvt = DatHist(HistThisCycle,3) ! Record the maximum number of events (lines) for any one history
         MaxEvtHist = History ! The history # with the most events (lines.)
        END IF
        CYCLE
       ELSE IF (DatFileInput(LinesThisCycle)%Col1 .NE. History) THEN ! New line is a new history.
         BACKSPACE 1 ! Backs the .DAT file up by one line.
         TotalLines = TotalLines - 1 ! Don't include this line in calculations since it's a new history.
         LinesThisCycle = LinesThisCycle - 1 ! Don't include this line in calculations since it's a new history.
         EXIT
       ELSE
         STOP "The program has experienced an unexpected error and stopped." ! This is just here to catch bad logic.
       END IF
        
     END DO
   END IF

  ! End of the .DAT file.
   IF (InputStatus < 0) THEN
    WRITE(*,*) "Last History in .DAT file: ", History
    IF(nps == 0) THEN
      nps = History
      WRITE(*,*) "NPS set to ", NINT(REAL(History) / SrcNorm), " based on the last history number and a source "
      WRITE(*,*) "normalization of ", SrcNorm
    END IF
      StopAfterLoop = .TRUE.
      LinesThisCycle = LinesThisCycle - 1 ! -1 because no history was read, just an end of file.
      TotalLines = TotalLines - 1 ! -1 because no history was read, just an end of file.
   END IF 

   ! MaxLines or MaxHist have been read from the .DAT file
   IF(TotalLines >= MaxLines .OR. TotalHist >= MaxHist) THEN
     WRITE(*,*) "Warning, .DAT file is too large.  Only histories up to number "
     WRITE(*,*)  History, " were processed."
     nps = History
     WRITE(*,*) "NPS adjusted to ", NINT(REAL(History) / SrcNorm), " based on the last history number read "
     WRITE(*,*) "and a source normalization of ", SrcNorm
     StopAfterLoop = .TRUE.
   END IF

   ! LinesPerCycle lines have been read in. Unlike the other two, this one is a temporary break to 
   ! process data (and keep array sizes manageable.)


 ELSE ! None of the 4 stopping criteria have been met, so read the data normally.
  IF (MOD(TotalLines, LinesPerCycle)==0) WRITE(*,*) TotalLines, " lines read"
  IF (History .NE. OldHistory) THEN ! Start of a new history
   TotalHist = TotalHist + 1
   HistThisCycle = HistThisCycle + 1
   DatHist(HistThisCycle,1) = History
   DatHist(HistThisCycle,2) = LinesThisCycle
   DatHist(HistThisCycle,3) = 1
   OldHistory = History
  ELSE ! This is another event in a history that has already been started
   DatHist(HistThisCycle,3) = DatHist(HistThisCycle,3) + 1
   IF (DatHist(HistThisCycle,3) > MaxEvt) THEN
    MaxEvt = DatHist(HistThisCycle,3) ! Record the maximum number of events (lines) for any one history
    MaxEvtHist = History ! The history # with the most events (lines.)
   END IF
  END IF

 END IF ! This is the end of IF statement at the very beginning of the M loop




 IF(EndOfCycle) THEN ! Stop reading histories and process the data read during this cycle
  CALL ComputePulses(FirstCycle, StopAfterLoop)
  CALL ComputePHS(FirstCycle, StopAfterLoop)
  CALL ComputeMultiplicities(FirstCycle, StopAfterLoop)
  CALL ComputeCorrelations(FirstCycle, StopAfterLoop)
  CALL FindPeaks(FirstCycle, StopAfterLoop)
  CALL ComputeAntiC(FirstCycle, StopAfterLoop)


 ! Get ready for the next cycle now that the data has been processed.
  LinesThisCycle = 0
  HistThisCycle = 0
  FirstCycle = .False.
  EndOfCycle = .False.
 END IF

 IF (StopAfterLoop) THEN
   CALL ComputePeaks
   CALL WriteOutput
   EXIT
 ELSE
  CYCLE
 END IF

END DO ! This is the end of the main (M) loop of the program.

CLOSE(UNIT=1)

DEALLOCATE(DatHist)
DEALLOCATE(Pulses)

!DEALLOCATE (Correlation, CrossCorr, Multiplicity, Peaks, PHSpectrum)

CALL CPU_TIME(FinishTime)
WRITE (*, '(A, F8.3, A)') " The program took ", FinishTime - StartTime, " seconds to finish."

CONTAINS

!==========================================================================================
! This subroutine reads all of the parameters from the input file so that they don't
! have to be hard-coded or entered in at the command line.
!==========================================================================================
SUBROUTINE ReadParams()

OPEN(UNIT=101, FILE=TRIM(ParamFile), ACTION="READ", STATUS="OLD", IOSTAT=InputStatus)
IF (InputStatus > 0) STOP "Could not find specified input parameter file."

DO I = 1, 10000
  READ(101, '(A)', IOSTAT=InputStatus) ThisLine
  IF (InputStatus < 0) EXIT !End of parameter file

  IF (ThisLine(1:7) == "DATNAME") THEN
  DO J = 8, 80
    IF (ThisLine(J:J) /= " ") THEN
      DatName = TRIM(ThisLine(J:80))
      EXIT
    END IF
  END DO
  ELSE IF (ThisLine(1:8) == "STARTDET") THEN 
    READ(ThisLine(9:80), *) StartDet
  ELSE IF (ThisLine(1:3) == "NPS") THEN 
    READ(ThisLine(4:80), *) NPS
  ELSE IF (ThisLine(1:7) == "SRCNORM") THEN 
    READ(ThisLine(8:80), *) SrcNorm
  ELSE IF (ThisLine(1:7) == "CORRMIN") THEN 
    READ(ThisLine(8:80), *) CorrMin
  ELSE IF (ThisLine(1:7) == "CORRMAX") THEN 
    READ(ThisLine(8:80), *) CorrMax
  ELSE IF (ThisLine(1:8) == "NUMTBINS") THEN 
    READ(ThisLine(9:80), *) NumTBins
  ELSE IF (ThisLine(1:8) == "DEADTIME") THEN 
    READ(ThisLine(9:80), *) DeadTime
  ELSE IF (ThisLine(1:7) == "NTHRESH") THEN 
    READ(ThisLine(8:80), *) nthresh
  ELSE IF (ThisLine(1:8) == "SDTHRESH") THEN 
    READ(ThisLine(9:80), *) sdthresh
  ELSE IF (ThisLine(1:6) == "COEFF0") THEN 
    READ(ThisLine(7:80), *) Coeff0
  ELSE IF (ThisLine(1:6) == "COEFF1") THEN 
    READ(ThisLine(7:80), *) Coeff1
  ELSE IF (ThisLine(1:6) == "COEFF2") THEN 
    READ(ThisLine(7:80), *) Coeff2
  ELSE IF (ThisLine(1:8) == "PGENTIME") THEN 
    READ(ThisLine(9:80), *) pgentime
  ELSE IF (ThisLine(1:9) == "PEAKWIDTH") THEN 
    READ(ThisLine(10:80), *) PeakWidth
  ELSE IF (ThisLine(1:6) == "PHSMAX") THEN 
    READ(ThisLine(7:80), *) PHSMax
  ELSE IF (ThisLine(1:7) == "PHSBINS") THEN 
    READ(ThisLine(8:80), *) PHSNumBins
  ELSE IF (ThisLine(1:6) == "SRCWTH") THEN 
    READ(ThisLine(7:80), *) SrcWidth
  ELSE IF (ThisLine(1:7) == "SRCFWHM") THEN 
    READ(ThisLine(8:80), *) SrcFWHM
  ELSE IF (ThisLine(1:7) == "DETFWHM") THEN 
    READ(ThisLine(8:80), *) DetFWHM
  ELSE IF (ThisLine(1:7) == "DETRISE") THEN 
    READ(ThisLine(8:80), *) DetRise
  ELSE IF (ThisLine(1:9) == "DETDECAYS") THEN 
    READ(ThisLine(10:80), *) DetDecayShort
  ELSE IF (ThisLine(1:8) == "PEPERMEV") THEN 
    READ(ThisLine(9:80), *) PEperMEV
  ELSE IF (ThisLine(1:9) == "MINDETNUM") THEN 
    READ(ThisLine(10:80), *) MinDetNum
  ELSE IF (ThisLine(1:9) == "MAXDETNUM") THEN 
    READ(ThisLine(10:80), *) MaxDetNum
  ELSE IF (ThisLine(1:6) == "EGAMMA") THEN 
    READ(ThisLine(7:80), *) E_Gamma
  ELSE IF (ThisLine(1:7) == "MAXHIST") THEN 
    READ(ThisLine(8:80), *) MaxHist
  ELSE IF (ThisLine(1:8) == "MAXLINES") THEN 
    READ(ThisLine(9:80), *) MaxLines
  ELSE IF (ThisLine(1:8) == "HITSFILE") THEN 
    READ(ThisLine(9:80), *) HitsFile
  ELSE IF (ThisLine(1:7) == "WINDMIN") THEN 
    READ(ThisLine(8:80), *) WindMin
  ELSE IF (ThisLine(1:7) == "WINDMAX") THEN 
    READ(ThisLine(8:80), *) WindMax
  END IF
END DO

CLOSE(UNIT=101)

END SUBROUTINE ReadParams

!==========================================================================================
! Subroutine PulseHeight converts the energy deposited in a collision (MeV) to light
! output in the detector (in MeVee)
!==========================================================================================
SUBROUTINE PulseHeight(Zprojectile, ZRx, Ztarget, ZEin, Zpheight)

INTEGER(KIND=4), INTENT(IN) :: Zprojectile, Ztarget, ZRx
REAL(KIND=4), INTENT(IN) :: ZEin

REAL(KIND=4), INTENT(OUT) :: Zpheight

! Neutron elastic scattering on hydrogen
IF (Zprojectile == 1 .AND. Ztarget == 1001 .AND. ZRx == -99) THEN
  Zpheight = Coeff2 * ZEin ** 2 + Coeff1 * ZEin + Coeff0
  RETURN
! Helium-3 Detector
ELSEIF (Zprojectile == 1 .AND. Ztarget == 2003 .AND. ZRx == 0) THEN
  Zpheight = ZEin
  RETURN
! Boron-10 detector
ELSEIF (Zprojectile == 1 .AND. Ztarget == 5010 .AND. ZRx == 0) THEN
  Zpheight = ZEin
  RETURN
! Inelastic scattering on carbon.
ELSEIF (Ztarget == 6000 .AND. Zprojectile == 1) THEN
  Zpheight = 0.02 * ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 92235 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 92238 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94238 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94239 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94240 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94241 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94242 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 95241 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 2 .AND. Ztarget == 1) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 2 .AND. Ztarget == 6) THEN
  Zpheight = ZEin
  RETURN
! Sodium (in an NaI scintillator)
ELSEIF (Zprojectile == 2 .AND. Ztarget == 11) THEN
  Zpheight = ZEin
  RETURN
! Iodine (in an NaI detector)
ELSEIF (Zprojectile == 2 .AND. Ztarget == 53) THEN
  Zpheight = ZEin
  RETURN
! Germanium Detector
ELSEIF (Zprojectile == 2 .AND. Ztarget == 32) THEN
  Zpheight = ZEin
  RETURN
ELSE
  Zpheight = 0.
  RETURN
END IF

END SUBROUTINE PulseHeight

!==========================================================================================
! Subroutine New_PulseHeight converts the energy deposited in a collision (MeV) to light
! output in the detector (in MeVee).  This new (and hopefully improved) version was added
! in Version 5.
!==========================================================================================
SUBROUTINE New_PulseHeight(Zprojectile, ZRx, Ztarget, ZEin, Zpheight)

INTEGER(KIND=4), INTENT(IN) :: Zprojectile, Ztarget, ZRx
REAL(KIND=4), INTENT(IN) :: ZEin

REAL(KIND=4), INTENT(OUT) :: Zpheight

! Neutron elastic scattering on hydrogen
IF (Zprojectile == 1 .AND. Ztarget == 1001 .AND. ZRx == -99) THEN
  Zpheight = Coeff2 * ZEin ** 2 + Coeff1 * ZEin + Coeff0
  RETURN
! (n,gamma) on hydrogen
ELSEIF (Zprojectile == 1 .AND. Ztarget == 1001 .AND. ZRx == 0) THEN
  Zpheight = Coeff2 * ZEin ** 2 + Coeff1 * ZEin + Coeff0
  RETURN
! Absorption in a He-3 detector.
ELSEIF (Zprojectile == 1 .AND. Ztarget == 2003 .AND. ZRx == 0) THEN
  Zpheight = ZEin
  RETURN
! Absorption in a B-10 detector.
ELSEIF (Zprojectile == 1 .AND. Ztarget == 5010 .AND. ZRx == 0) THEN
  Zpheight = ZEin
  RETURN
! Absorption on carbon.  Could be either (n,gamma) or (n,alpha) reaction.
ELSEIF (Ztarget == 6000 .AND. Zprojectile == 1 .AND. ZRx == 0) THEN
  ! Assume the reaction is (n,gamma) because the (n,alpha) xs is ~0 at this energy.
  IF (ZEin < 1.3) THEN
   Zpheight = 0.00006174 * ZEin ** 2 + 0.0047094 * ZEin + 0.0014204
  ! (n,alpha) >> (n,gamma) above this energy deposition level, so assume everything is
  ! (n,alpha) above this level.
  ELSE
   ! This is the light generated by the Be-9 recoil.  I'm just assuming it generates the 
   ! same amount of light as a carbon nucleus.  This is probably an underestimation, but
   ! the alpha should produce the majority of the light, so the difference should be
   ! negligible.
   Zpheight = 0.00006174 * (ZEin*4./13.) ** 2 + 0.0047094 * (ZEin*4./13.) + 0.0014204
   ! Now add the alpha contribution to the light.
   ZPheight = Zpheight + 0.014432*(ZEin*9./13.)**2 - 0.013526*(ZEin*4./13.) + 0.0096963
  END IF
  RETURN
! Elastic Scattering on carbon.
ELSEIF (Ztarget == 6000 .AND. Zprojectile == 1 .AND. ZRx == -99) THEN
  Zpheight = 0.00006174 * ZEin ** 2 + 0.0047094 * ZEin + 0.0014204
  RETURN
! Inelastic scattering on carbon.
ELSEIF (Ztarget == 6000 .AND. Zprojectile == 1 .AND. ABS(ZRx) == 1) THEN
  Zpheight = 0.00006174 * ZEin ** 2 + 0.0047094 * ZEin + 0.0014204
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 92235 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 92238 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94238 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94239 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94240 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94241 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 94242 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 1 .AND. Ztarget == 95241 .AND. ZRx == 19) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 2 .AND. Ztarget == 1) THEN
  Zpheight = ZEin
  RETURN
ELSEIF (Zprojectile == 2 .AND. Ztarget == 6) THEN
  Zpheight = ZEin
  RETURN
! Sodium (in an NaI scintillator)
ELSEIF (Zprojectile == 2 .AND. Ztarget == 11) THEN
  Zpheight = ZEin
  RETURN
! Iodine (in an NaI detector)
ELSEIF (Zprojectile == 2 .AND. Ztarget == 53) THEN
  Zpheight = ZEin
  RETURN
! Germanium Detector
ELSEIF (Zprojectile == 2 .AND. Ztarget == 32) THEN
  Zpheight = ZEin
  RETURN
ELSE
  Zpheight = 0.
  RETURN
END IF

END SUBROUTINE New_PulseHeight

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Added with Version 4.0.  This subroutine samples from a Gaussian distribution
! with a mean of zero and a FWHM provided by the user.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SUBROUTINE SAMPLE_GAUSSIAN(ZFWHM,ZRV)
REAL(KIND=4), INTENT(IN) :: ZFWHM
REAL(KIND=4), INTENT(OUT) :: ZRV
REAL(KIND=4) :: Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8, Z9, Z10, Z11, Z12

CALL RANDOM_NUMBER(Z1)
CALL RANDOM_NUMBER(Z2)
CALL RANDOM_NUMBER(Z3)
CALL RANDOM_NUMBER(Z4)
CALL RANDOM_NUMBER(Z5)
CALL RANDOM_NUMBER(Z6)
CALL RANDOM_NUMBER(Z7)
CALL RANDOM_NUMBER(Z8)
CALL RANDOM_NUMBER(Z9)
CALL RANDOM_NUMBER(Z10)
CALL RANDOM_NUMBER(Z11)
CALL RANDOM_NUMBER(Z12)

ZRV = (Z1+Z2+Z3+Z4+Z5+Z6+Z7+Z8+Z9+Z10+Z11+Z12) - 6.0
ZRV = ZRV * ZFWHM / 2.35482
RETURN
END SUBROUTINE

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Added with Version 5.0.  This subroutine samples from a uniform distribution
! with a width provided via the input file.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SUBROUTINE PulseWidth(ZSourceWidth,ZSourceStart)
REAL(KIND=4), INTENT(IN) :: ZSourceWidth
REAL(KIND=4), INTENT(OUT) :: ZSourceStart
REAL(KIND=4) :: ZZ1

CALL RANDOM_NUMBER(ZZ1)

ZSourceStart = ZZ1 * ZSourceWidth
RETURN

END SUBROUTINE

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! Added with Version 5.6.  This subroutine samples the detector timing 
! uncertainty from an exponential sum consisting of a rise component and a
! decay component. This distribution tends to have a long tail (skewed right.)
! The time is adjusted so that the peak of the distribution corresponds to a 
! DeltaT of 0.
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
SUBROUTINE FIND_DETECTOR_TIME(TauRise,TauDecay1,StartTime)
REAL(KIND=4), INTENT(OUT) :: StartTime
INTEGER(KIND=4) :: I1
REAL(KIND=4) :: Y1, Y2, YMaxTime, YPDFVal
REAL(KIND=4), INTENT(IN) :: TAURise, TauDecay1
! This value must be larger than the largest possible of the pdf being
! sampled or else the start times will be sampled incorrectly.
REAL(KIND=4) :: PDFMax, PDFPeakTime

PDFPeakTime = (LOG(TauDecay1) - LOG(TauRise)) / (1./TauRise - 1./TauDecay1)
PDFMax = EXP(-PDFPeakTime/TauDecay1) - EXP(-PDFPeakTime/TauRise)
! The PDF is sampled 7 decay times past the peak, which should capture all PDF values
! larger than 0.1% of the peak value.
YMaxTime = PDFPeakTime + 7.0 * TauDecay1

DO I1 = 1, 1000000
 CALL RANDOM_NUMBER(Y1)
 CALL RANDOM_NUMBER(Y2)
 Y1 = Y1 * YMaxTime
 YPDFVal = (EXP(-Y1/TauDecay1) - EXP(-Y1/TauRise))

 IF ((YPDFVal / PDFMax) >= Y2) THEN
  ! DeltaT is adjusted so that the peak time represents a DeltaT of 0.
  StartTime = Y1 - PDFPeakTime
  EXIT
 ELSE
  CYCLE
 END IF
END DO

RETURN

END SUBROUTINE

SUBROUTINE AssignAndSortDetectorNumbers(TempCellNum)
INTEGER(KIND=4), INTENT(IN) :: TempCellNum
INTEGER(KIND=4) :: SwapInt, MinRowLoc(1), II
LOGICAL :: NewDet

 ! This is the very first valid detector cell in the .DAT file.
 IF(NumDet == 0 .AND. TempCellNum >= MinDetNum .AND. TempCellNum <= MaxDetNum) THEN 
  NumDet = 1
  Detector(1) = TempCellNum
  RETURN
 END IF

 NewDet = .TRUE.
 ! Determine if this is a new (valid) detector cell, and if so add it to the 
 IF (TempCellNum >= MinDetNum .AND. TempCellNum <= MaxDetNum) THEN
  DO II = 1, NumDet
    IF (TempCellNum == Detector(II)) THEN
      NewDet = .FALSE.
      EXIT   
    END IF
  END DO
 ELSE
  NewDet = .FALSE.
 END IF

 IF(NewDet) THEN ! If this is a new detector, resort the detector cell numbers from low to high
  NumDet = NumDet + 1 ! New detector found, so increment the number of detectors by 1.
  Detector(NumDet) = TempCellNum
!  WRITE(*,*) NumDet, TempCellNum, Detector(NumDet), NewDet
  IF (NumDet > MaxDets) THEN
    WRITE(*,*) "Error! Number of detector cells in the .DAT file exceeds the maximum value of"
    WRITE(*,*) MaxDets, " detector cells. Recompile with a larger MaxDets value if required."
    STOP
  END IF
  ! Sorts the detector cell numbers from low to high.
  DO II = 1, NumDet-1
  MinRowLoc = II + MINLOC(Detector(II+1:NumDet))
   IF (Detector(II) .GT. Detector(MinRowLoc(1))) THEN
     SwapInt = Detector(II)
     Detector(II) = Detector(MinRowLoc(1))
     Detector(MinRowLoc(1)) = SwapInt
   END IF
  END DO

  ! Find the start detector number (if applicable) now that the detectors have been sorted.
  StartDetRow = 0
  DO II = 1, NumDet
    IF (Detector(II) == StartDet) StartDetRow = II
    EXIT
  END DO
 END IF


RETURN
END SUBROUTINE AssignAndSortDetectorNumbers

SUBROUTINE ComputePulses(ResetFlag, TerminateFlag)
LOGICAL,INTENT(IN) :: ResetFlag, TerminateFlag
!==========================================================================================
! Step 2 : Extract pulses for each history.  This also includes the addition of the source
! pulse time width, source time uncertainty, and the detector time and energy uncertainties.
!==========================================================================================
IF(ResetFlag) THEN
  ALLOCATE(Pulses(LinesPerCycle))
  TotalPulses = 0
  ! Open the .hits file if specified.
  IF(HitsFile == 1) THEN
    OPEN(UNIT=8,FILE=TRIM(FileBase)//".hits",ACTION="WRITE",STATUS="REPLACE")
  END IF
END IF

ALLOCATE(ThisHist(MaxEvt))
ALLOCATE(PHTotal(MaxEvt))
PHTotal = 0.

PulsesThisCycle = 0

! Loop through the history file.  
DO I = 1, HistThisCycle
  
  ! Determines the time of the alpha event by adding the pulse width and timing uncertainty.
  CALL PulseWidth(SrcWidth, SrcStart)
  CALL SAMPLE_GAUSSIAN(SrcFWHM, SrcDeltaT)

  ! Reads all of the events for history 'I' into the ThisHist array.
  DO J = 1, DatHist(I,3)
  ThisHist(J) = DatFileInput(DatHist(I,2)+J-1)

    ! Converts the time from shakes (10^-8 s.) to ns.
    ThisHist(J)%Col8 = 10. * ThisHist(J)%Col8

  END DO ! J Loop

  !==========================================================================================
  ! Step 2a : Sort Histories ascending by detector cell and then by time within the detector
  ! cell.  If there is a tie for both detector cell and time, the event with the lowest
  ! collision number and then code is brought up. 
  !==========================================================================================

  IF(DatHist(I,3) > 1) THEN
  DO J = 1, DatHist(I,3)-1
    DO K = J, DatHist(I,3)
     IF (ThisHist(J)%Col6 > ThisHist(K)%Col6) THEN
       HistSwap = ThisHist(J)
       ThisHist(J) = ThisHist(K)
       ThisHist(K) = HistSwap
     ELSEIF(ThisHist(J)%Col6 == ThisHist(K)%Col6 .AND. &
     & ThisHist(J)%Col8 > ThisHist(K)%Col8) THEN
       HistSwap = ThisHist(J)
       ThisHist(J) = ThisHist(K)
       ThisHist(K) = HistSwap
     ELSEIF(ThisHist(J)%Col6 == ThisHist(K)%Col6 .AND. &
     & ThisHist(J)%Col8 == ThisHist(K)%Col8 .AND. ThisHist(J)%Col14 > ThisHist(K)%Col14) THEN
       HistSwap = ThisHist(J)
       ThisHist(J) = ThisHist(K)
       ThisHist(K) = HistSwap
     ELSEIF(ThisHist(J)%Col6 == ThisHist(K)%Col6 .AND. &
     & ThisHist(J)%Col8 == ThisHist(K)%Col8 .AND. ThisHist(J)%Col14 == ThisHist(K)%Col14 &
     & .AND. ThisHist(J)%Col15 > ThisHist(K)%Col15) THEN
       HistSwap = ThisHist(J)
       ThisHist(J) = ThisHist(K)
       ThisHist(K) = HistSwap
     END IF
    END DO ! K Loop
  END DO ! J Loop
  END IF

  DO J = 1, DatHist(I,3)
    ! Assigns a light output of 0 to any pulses outside of the minimum and maximum detector numbers.
    ! This will prevent any pulses (and thus multiplicities, cross-correlations, etc.) from being
    ! recorded for those detectors.
    IF (ThisHist(J)%Col6 < MinDetNum .OR. ThisHist(J)%Col6 > MaxDetNum) THEN
      PHTotal(J) = 0.
      CYCLE
    ELSE
      CALL New_PulseHeight(ThisHist(J)%Col3, ThisHist(J)%Col4, ThisHist(J)%Col5, ThisHist(J)%Col7, PHThisEvent)
      PHTotal(J) = PHThisEvent
	END IF

  END DO  ! J Loop
  
  CurrentDet = 0

  DO J = 1, DatHist(I,3)

    ! Cycles if no light was deposited by this event.
    IF (PHTotal(J) <= 0.) CYCLE

	! New detector cell. Reset dead time and pulse generation time.
    IF (ThisHist(J)%Col6 /= CurrentDet) THEN
      ! Set the start of the deadtime and pulse generation time to a negative number large
      ! enough that the first event will not be in either window.
      DeadStart = -2.*deadtime
      PulseStart = -2.*pgentime
      CurrentDet = ThisHist(J)%Col6
    END IF

	! If the detector cell is dead because of a previous pulse, ignore this one.
    IF (ThisHist(J)%Col8 < DeadStart + deadtime) CYCLE

    ! These values are used to identify what type of pulse this is in the event that multiple
    ! events get summed together to generate it.
    MinSubHist = ThisHist(J)%Col2
    MinParNo = ThisHist(J)%Col3
    MinGen = ThisHist(J)%Col13
    MinCollisions = ThisHist(J)%Col14
    MinCode = ThisHist(J)%Col15
    MinWeight = ThisHist(J)%Col12

    ! All events that occur within the pulse generation time are summed together to
    ! determine the total light output in the (potential) pulse.
    DO K = J+1, DatHist(I,3)
     IF (ThisHist(K)%Col6 /= CurrentDet) EXIT
     IF (ThisHist(K)%Col8 > ThisHist(J)%Col8 + pgentime) EXIT
     PHTotal(J) = PHTotal(J) + PHTotal(K)

    END DO ! K Loop

    ! Added in Version 5.  This section accounts for the statistical variation in the 
    ! pulse height due to random fluctuations in the number of photons produced in the
    ! detector and the number of photoelectrons produced in the photomultiplier tube.
    NumPE = PEperMeV * PHTotal(J)
    NumPEFWHM = 2.35482 * SQRT(NumPE)
    CALL SAMPLE_GAUSSIAN(NumPEFWHM,DeltaPE)
    NumPE = NumPE + DeltaPE
    PHTotal(J) = NumPE / PEperMeV

    ! Added in version 5.6. Threshold is no longer treated as a simple binary. Now, it is
    ! sampled from a distribution to account for the fact that events in different parts of
    ! the detector have varying probabilities of having enough photons reach the PMT to
    ! overcome the threshold

    ! This determines the light output threshold in MeVee based on the neutron energy threshold
    ! input by the user.
    CALL SAMPLE_GAUSSIAN(2.35482*sdthresh,deltathresh)
    samplethresh = nthresh + deltathresh
    threshold = Coeff2*samplethresh**2 + Coeff1*samplethresh + Coeff0
    IF (PHTotal(J) >= threshold) THEN

      ! Add pulse time width, source time uncertainty, and detector time uncertainty 
      ! to the pulse time.
      IF (DetRise > 0. .AND. DetDecayShort > DetRise) THEN
        CALL FIND_DETECTOR_TIME(DetRise, DetDecayShort, DetDeltaT)
      ELSE
        CALL SAMPLE_GAUSSIAN(DetFWHM, DetDeltaT)
      END IF
      PulseTime = ThisHist(J)%Col8 + SrcStart + SrcDeltaT + DetDeltaT

      ! Check to see if the new pulse time falls within the correlation window.  Cycles the J
      ! loop if it doesn't.
      IF (PulseTime < CorrMin .OR. PulseTime > CorrMax) CYCLE

      PulsesThisCycle = PulsesThisCycle + 1
      TotalPulses = TotalPulses + 1
      Pulses(PulsesThisCycle)%HistNo = ThisHist(J)%Col1
      Pulses(PulsesThisCycle)%DetNo = ThisHist(J)%Col6
      Pulses(PulsesThisCycle)%ParNo = MinParNo
      Pulses(PulsesThisCycle)%Time = PulseTime
      Pulses(PulsesThisCycle)%PulseHeight = PHTotal(J)
      Pulses(PulsesThisCycle)%SubHist = MinSubHist
      Pulses(PulsesThisCycle)%GenNo = MinGen


      ! Check to see if the event that generated the pulse was a directly transmitted DT neutron.
      ! SubHistory = 1, Particle Type = 1 (neutron), Weight > 0.95, Generation = 0, 
      ! # Collisions = 0, Code = 0
      IF(MinSubHist == 1 .AND. MinParNo == 1 .AND. MinWeight >= 0.95 &
      & .AND. MinGen == 0 .AND. MinCollisions == 0 .AND. MinCode == 0) THEN
        Pulses(PulsesThisCycle)%Direct = 1
      ELSE
        Pulses(PulsesThisCycle)%Direct = 0
      END IF

      ! Check to see if the pulse height falls within the FWHM of the hydrogen absorption gamma
      ! ray (2.22 MeV).  Added in version 5.
      IF (PHTotal(J) >= HLineMin .AND. PHTotal(J) <= HLineMax) THEN
       Pulses(PulsesThisCycle)%HGamma = 1
      ELSE
       Pulses(PulsesThisCycle)%HGamma = 0
      END IF

      DeadStart = ThisHist(J)%Col8
      PulseStart = -2.*pgentime

      IF(HitsFile == 1) THEN
        WRITE(8,'(I12,I6,ES13.5)') Pulses(PulsesThisCycle)%HistNo, Pulses(PulsesThisCycle)%DetNo, Pulses(PulsesThisCycle)%Time
      END IF

      CYCLE
    ELSE
      PulseStart = ThisHist(J)%Col8
    END IF

  END DO ! J Loop

  PHTotal = 0.

END DO ! I Loop

DO I = 1, PulsesThisCycle
  Pulses(I)%XTalk = 0
!  WRITE(8,130, ADVANCE="NO") Pulses(I)
!  WRITE(8,'(I12)') INT(Pulses(I)%Time/TimeIncrement)
  ! This loop changes the detector number value from the actual cell number to the line
  ! referencing that cell number in the Detector() array.  This will make matching the 
  ! detector cell faster in the correlation section.
  DO J = 1, NumDet
    IF (Pulses(I)%DetNo == Detector(J)) THEN
      Pulses(I)%DetNo = J
      EXIT
    END IF
  END DO
END DO

DEALLOCATE(ThisHist)
DEALLOCATE(PHTotal)

IF(TerminateFlag) THEN
  CLOSE(UNIT=8)
END IF

RETURN
END SUBROUTINE ComputePulses




SUBROUTINE ComputePHS(ResetFlag, TerminateFlag)
LOGICAL,INTENT(IN) :: ResetFlag, TerminateFlag
INTEGER(KIND=4), ALLOCATABLE :: TempPHSpectrum(:,:,:)

!==========================================================================================
! Step 3 : Compute the Pulse Height Spectra.
!==========================================================================================

! Allocate the main array for the Pulse Height Spectrum. If this is not the first cycle,
! there may be more detectors now, so the PHSpectrum array is rebinned into the TempPHSpectrum
! array, and then copied back.
IF(ResetFlag) THEN
  ALLOCATE(PHSpectrum(-1:PHSNumBins,0:MaxDets,3))
  MaxPulse = 0
  PHSpectrum = 0
  DO I = 1, 3
    PHSpectrum(-1,0,I) = nps
    DO J = 1, NumDet
      PHSpectrum(-1,J,I) = Detector(J)
    END DO
  END DO
ELSE
  ALLOCATE(TempPHSpectrum(-1:PHSNumBins,0:MaxDets,3))
  DO I = 1, 3
    TempPHSpectrum(-1,0,I) = nps
    DO J = 1, NumDet
      DO K = 1, NumDet
        IF(PHSpectrum(-1,K,I) == Detector(J)) THEN
          TempPHSpectrum(:,J,:) = PHSpectrum(:,K,:)
        END IF
      END DO
    END DO
  END DO
  PHSpectrum = TempPHSpectrum
  DEALLOCATE(TempPHSpectrum)
  DEALLOCATE(PulseHist)
END IF
ALLOCATE(TempHist(PulsesThisCycle,2))

NumHistP = 0
OldHistory = 0




DO I = 1, PulsesThisCycle

PHSBin = INT(Pulses(I)%PulseHeight / PHSIncrement)+1
IF (PHSBin > PHSNumBins) PHSBin = PHSNumBins
PHSpectrum(PHSBin,Pulses(I)%DetNo,1) = PHSpectrum(PHSBin,Pulses(I)%DetNo,1) + 1
IF (Pulses(I)%ParNo == 1) THEN
  PHSpectrum(PHSBin,Pulses(I)%DetNo,2) = PHSpectrum(PHSBin,Pulses(I)%DetNo,2) + 1 
ELSE IF (Pulses(I)%ParNo == 2) THEN
  PHSpectrum(PHSBin,Pulses(I)%DetNo,3) = PHSpectrum(PHSBin,Pulses(I)%DetNo,3) + 1 
END IF

! qqqqq
! IF (Pulses(I)%Direct == 1 .AND. Pulses(I)%Time > 50) WRITE(*,*) Pulses(I)%HistNo


  IF (Pulses(I)%HistNo == OldHistory) THEN
    CYCLE
  ELSE
    OldHistory = Pulses(I)%HistNo
    NumHistP = NumHistP + 1
    TempHist(NumHistP,1) = Pulses(I)%HistNo
    TempHist(NumHistP,2) = I
  END IF
END DO

ALLOCATE(PulseHist(NumHistP,3))

DO I = 1, NumHistP
  PulseHist(I,1) = TempHist(I,1) ! History number of a history with a pulse
  PulseHist(I,2) = TempHist(I,2) ! Starting line number of that history in Pulses array.
  IF (I < NumHistP) THEN
    PulseHist(I,3) = TempHist(I+1,2) - TempHist(I,2)
  ELSE
    PulseHist(I,3) = PulsesThisCycle+1 - TempHist(I,2)
  END IF
  IF (PulseHist(I,3) > MaxPulse) THEN
    MaxPulse = PulseHist(I,3)
    MaxPulseHist = PulseHist(I,1)
  END IF
END DO

DEALLOCATE(TempHist)
RETURN
END SUBROUTINE ComputePHS



SUBROUTINE ComputeMultiplicities(ResetFlag, TerminateFlag)

!==========================================================================================
! Step 4 : Extract multiplicities from pulse data.
!==========================================================================================
LOGICAL,INTENT(IN) :: ResetFlag, TerminateFlag
INTEGER(KIND=4) :: TotalMult, TotalNeutronMult, NumGen0, NumGen1, TotalGen0Mult
INTEGER(KIND=4) :: NumPulses, TotalGen1Mult, TotalHGammaMult, NumNeutronPulses, NumHGammas
INTEGER(KIND=4), ALLOCATABLE :: TempMultiplicity(:,:)

! First time through - allocate multiplicity array.
IF(ResetFlag) THEN
  ALLOCATE(Multiplicity(0:MaxPulse,5))
! First time through the multiplicity loop - initialize counting variables
  Multiplicity = 0
! 2nd+ Time through. Because MaxPulse may be larger now, I transfer the old values into a temporary array,
! then reallocate multiplicity with the new MaxPulse value and copy the old values back.
ELSE
  ALLOCATE(TempMultiplicity(0:MaxPulse,5))
  TempMultiplicity = 0
  DO I = 0, UBOUND(Multiplicity,1)
    TempMultiplicity(I,:) = Multiplicity(I,:)
  END DO
  DEALLOCATE(Multiplicity)
  ALLOCATE(Multiplicity(0:MaxPulse,5))
  Multiplicity = TempMultiplicity
  DEALLOCATE(TempMultiplicity)
END IF

WRITE(*,*) "Calculating Multiplicities"

! Multiplicity structure: Rows = number of source triggered multiplicities; Columns - 
! Column 1 = total multiplicities; Column 2 = neutron multiplicities; Column 3 =
! hydrogen gamma ray multiplicities; Column 4 = Generation 0; Column 5 = Generation 1+

! Go through the pulse file and record the number of pulses created for each history in
! the Multiplicity array.
DO I = 1, NumHistP

  NumPulses = 0
  NumNeutronPulses = 0
  NumHGammas = 0
  NumGen0 = 0
  NumGen1 = 0
  DO J = 1, PulseHist(I,3)
    IF (Pulses(PulseHist(I,2)+J-1)%Time < WindMin .OR. Pulses(PulseHist(I,2)+J-1)%Time > WindMax) CYCLE
    NumPulses = NumPulses + 1
    IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 1) NumNeutronPulses = NumNeutronPulses + 1
    IF (Pulses(PulseHist(I,2)+J-1)%HGamma == 1) THEN
     NumHGammas = NumHGammas + 1
     TotalHGammas = TotalHGammas + 1
    END IF
    IF (Pulses(PulseHist(I,2)+J-1)%GenNo == 0 .AND. Pulses(PulseHist(I,2)+J-1)%ParNo == 1) NumGen0 = NumGen0 + 1
    IF (Pulses(PulseHist(I,2)+J-1)%GenNo > 0 .AND. Pulses(PulseHist(I,2)+J-1)%ParNo == 1) NumGen1 = NumGen1 + 1
  END DO
  Multiplicity(NumPulses,1) = Multiplicity(NumPulses,1) + 1
  Multiplicity(NumNeutronPulses, 2) = Multiplicity(NumNeutronPulses, 2) + 1
  Multiplicity(NumHGammas, 3) = Multiplicity(NumHGammas, 3) + 1
  Multiplicity(NumGen0, 4) = Multiplicity(NumGen0, 4) + 1
  Multiplicity(NumGen1, 5) = Multiplicity(NumGen1, 5) + 1

END DO



IF(TerminateFlag) THEN
  ! Assigns to 'TotalMult', et al. the total number of histories which have >= 1 pulse.
  TotalMult = 0
  TotalNeutronMult = 0
  TotalHGammaMult = 0
  TotalGen0Mult = 0
  TotalGen1Mult = 0
  DO I = 1, MaxPulse
    TotalMult = TotalMult + Multiplicity(I, 1)
    TotalNeutronMult = TotalNeutronMult + Multiplicity(I, 2) 
    TotalHGammaMult = TotalHGammaMult + Multiplicity(I,3)
    TotalGen0Mult = TotalGen0Mult + Multiplicity(I,4)
    TotalGen1Mult = TotalGen1Mult + Multiplicity(I,5)
  END DO
  ! Total histories minus all histories with >= 1 pulse = no. of histories with no pulses.
  Multiplicity(0, 1) = nps - TotalMult
  Multiplicity(0, 2) = nps - TotalNeutronMult
  Multiplicity(0, 3) = nps - TotalHGammaMult
  Multiplicity(0, 4) = nps - TotalGen0Mult
  Multiplicity(0, 5) = nps - TotalGen1Mult
END IF

RETURN
END SUBROUTINE ComputeMultiplicities



SUBROUTINE ComputeCorrelations(ResetFlag, TerminateFlag)
LOGICAL,INTENT(IN) :: ResetFlag, TerminateFlag
INTEGER(KIND=4), ALLOCATABLE :: TempCorrelation(:,:,:), TempCrossCorr(:,:,:,:)

!==========================================================================================
! Step 5 : Calculate source-detector correlations
!==========================================================================================
130 FORMAT(I11, I4, I5, I4, I5, I3, I4, I4, F15.3, F12.6)

IF(ResetFlag) THEN
  CorrWindowOverflow = 0
  CorrWindowUnderflow = 0
  ! Correlation Structure - Rows = Time, Columns = Det #, Panes = Total, n, gamma, direct, 
  ! no x-talk, H_Gamma, Gen. 0, Gen. 1+
  ALLOCATE(Correlation(-1:NumTBins,0:MaxDets,8))
  Correlation = 0
! Input detector cell numbers and time steps into correlation and cross-correlation arrays.
  DO I = 1, 8
    Correlation(-1,0,I) = NINT(REAL(nps) / SrcNorm)
    DO J = 1, NumDet
      Correlation(-1,J,I) = Detector(J)
    END DO
    DO J = 0, NumTBins
      Correlation(J,0,I) = J
    END DO
  END DO

  ! CrossCorr Structure - Rows = Time Lag, Columns = Stop Det #, Panes = Total, nn, gg, np, 
  ! pn, H_Gamma-H_Gamma, Gen0-Gen0, Induced-Induced;    4th Dimension = Start Det #.
  IF (StartDet < 0) THEN
    ALLOCATE(CrossCorr(-NumTBins-1:NumTBins+1,0:MaxDets,8,MaxDets))
    ALLOCATE(EventScatter(0:NumTBins,0:NumTBins,1))
    EventScatter = 0
  ELSEIF(StartDet > 0) THEN
    ALLOCATE(CrossCorr(-NumTBins-1:NumTBins+1,0:MaxDets,8,1))
    ALLOCATE(EventScatter(0:NumTBins,0:NumTBins,1))
    EventScatter = 0
  END IF

  IF(ALLOCATED(CrossCorr)) CrossCorr=0

  IF(StartDet .NE. 0) THEN
   DO I = 1, 8
    CrossCorr(-NumTBins-1,0,I,1) = NINT(REAL(nps) / SrcNorm)
    DO J = 1, NumDet
      IF (StartDet < 0) THEN
       DO K = 1, NumDet
        CrossCorr(-NumTBins-1,J,I,K) = Detector(J)
        CrossCorr(NumTBins+1,J,I,K) = Detector(K)
       END DO
      ELSE
       CrossCorr(-NumTBins-1,J,I,1) = Detector(J)
      END IF
    END DO
    DO J = -NumTBins, NumTBins
     IF (StartDet < 0) THEN
      DO K = 1, NumDet
       CrossCorr(J,0,I,K) = J
      END DO
     ELSE
      CrossCorr(J,0,I,1) = J
     END IF
    END DO
   END DO
  END IF


ELSE
  ALLOCATE(TempCorrelation(-1:NumTBins,0:MaxDets,8))
  TempCorrelation = 0
  DO I = 1, 8
    TempCorrelation(-1,0,I) = NINT(REAL(nps) / SrcNorm)
    DO J = 1, NumDet
      TempCorrelation(-1,J,I) = Detector(J)
    END DO
    DO J = 0, NumTBins
      TempCorrelation(J,0,I) = J
    END DO
  END DO
  DO I = 1, NumDet
    DO J = 1, NumDet
      IF(TempCorrelation(-1,I,1) == Correlation(-1,J,1)) TempCorrelation(:,I,:) = Correlation(:,J,:)
    END DO
  END DO
  Correlation = TempCorrelation
  DEALLOCATE(TempCorrelation)

  IF (StartDet < 0) THEN
    ALLOCATE(TempCrossCorr(-NumTBins-1:NumTBins+1,0:MaxDets,8,MaxDets))
  ELSEIF(StartDet > 0) THEN
    ALLOCATE(TempCrossCorr(-NumTBins-1:NumTBins+1,0:MaxDets,8,1))
  END IF

  IF(ALLOCATED(TempCrossCorr)) TempCrossCorr=0

  IF(StartDet .NE. 0) THEN
   DO I = 1, 8
    TempCrossCorr(-NumTBins-1,0,I,1) = NINT(REAL(nps) / SrcNorm)
    DO J = 1, NumDet
      IF (StartDet < 0) THEN
       DO K = 1, NumDet
        TempCrossCorr(-NumTBins-1,J,I,K) = Detector(J)
        TempCrossCorr(NumTBins+1,J,I,K) = Detector(K)
       END DO
      ELSE
       TempCrossCorr(-NumTBins-1,J,I,1) = Detector(J)
      END IF
    END DO
    DO J = -NumTBins, NumTBins
     IF (StartDet < 0) THEN
      DO K = 1, NumDet
       TempCrossCorr(J,0,I,K) = J
      END DO
     ELSE
      TempCrossCorr(J,0,I,1) = J
     END IF
    END DO
   END DO
  END IF

  IF(StartDet < 0) THEN
    DO I = 1, NumDet
     DO J = 1, NumDet
      DO K = 1, NumDet
       DO L = 1, NumDet
         IF(CrossCorr(-NumTBins-1,I,1,1) == TempCrossCorr(-NumTBins-1,K,1,1) .AND. &
         & CrossCorr(NumTBins+1,1,1,J) == TempCrossCorr(NumTBins+1,1,1,L)) THEN
           TempCrossCorr(:,K,:,L) = CrossCorr(:,I,:,J)
         END IF
       END DO
      END DO
     END DO
    END DO
  ELSEIF(StartDet > 0) THEN
    DO I = 1, NumDet
     DO K = 1, NumDet
       IF(CrossCorr(-NumTBins-1,I,1,1) == TempCrossCorr(-NumTBins-1,K,1,1)) THEN
         TempCrossCorr(:,K,:,1) = CrossCorr(:,I,:,1)
       END IF
     END DO
    END DO
  END IF

  IF(ALLOCATED(CrossCorr) .AND. ALLOCATED(TempCrossCorr)) THEN
    CrossCorr = TempCrossCorr
  END IF


  IF(ALLOCATED(TempCrossCorr)) DEALLOCATE(TempCrossCorr)

END IF
! qqqqq
! OPEN (UNIT=99, FILE="error.out", STATUS="REPLACE", ACTION="WRITE")

WRITE(*,*) "Calculating Source-Detector Correlations"

DO I = 1, NumHistP
  DO J = 1, PulseHist(I,3)
    ! IntTime is now the time bin of the pulse rather than (necessarily) the time.
    IntTime = NINT(Pulses(PulseHist(I,2)+J-1)%Time / TimeIncrement)
    IF (IntTime > NumTBins) THEN
      CorrWindowOverflow = CorrWindowOverflow + 1
      CYCLE
!    ESLEIF (IntTime < INT(CorrMin/TimeIncrement)) THEN
!      CorrWindowUnderflow = CorrWindowUnderflow + 1
!      CYCLE
    END IF
    CurrentDet = Pulses(PulseHist(I,2)+J-1)%DetNo
    Correlation(IntTime,CurrentDet,1) = Correlation(IntTime,CurrentDet,1) + 1
    IF (Pulses(PulseHist(I,2)+J-1)%Direct == 1) THEN
      Correlation(IntTime,CurrentDet,4) = Correlation(IntTime,CurrentDet,4) + 1
    END IF
    IF (Pulses(PulseHist(I,2)+J-1)%HGamma == 1) THEN
      Correlation(IntTime,CurrentDet,6) = Correlation(IntTime,CurrentDet,6) + 1
    END IF
    IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 1) THEN
      Correlation(IntTime,CurrentDet,2) = Correlation(IntTime,CurrentDet,2) + 1
      ! Check to see if the neutron is generation 0 or generation 1+ (induced)
      IF (Pulses(PulseHist(I,2)+J-1)%GenNo == 0) THEN
        Correlation(IntTime,CurrentDet,7) = Correlation(IntTime,CurrentDet,7) + 1
      ELSE
        Correlation(IntTime,CurrentDet,8) = Correlation(IntTime,CurrentDet,8) + 1
      END IF
    ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 2) THEN
      Correlation(IntTime,CurrentDet,3) = Correlation(IntTime,CurrentDet,3) + 1
    ELSE
      WRITE(*,*) "Error!  Particle type other than a photon or a neutron encountered &
      &in .DAT file, history number ", PulseHist(I,1), ".  Check .DAT file."
      WRITE (*, 130) Pulses(PulseHist(I,2)+J-1)
    END IF

!==========================================================================================
! Step 5a : Calculate detector-detector correlations if start detector is specified
!==========================================================================================

   IF (StartDet > 0) THEN

    ! Check cross-correlations and fill in the cross-correlation arrays.
    IF (CurrentDet == StartDetRow) THEN
      DO K = 1, PulseHist(I,3)
        IF (K == J) CYCLE
        StopDetRow = Pulses(PulseHist(I,2)+K-1)%DetNo

        ! Version 5 - because of the change in the time bin structure so that the bins
        ! no longer (necessarily) match the times attached, TimeLag is now the number of
        ! bins between the two pulses rather than the time between pulses.
        TimeLag = NINT((Pulses(PulseHist(I,2)+K-1)%Time - &
        & Pulses(PulseHist(I,2)+J-1)%Time) / TimeIncrement)

        CrossCorr(TimeLag,StopDetRow,1,1) = CrossCorr(TimeLag,StopDetRow,1,1) + 1
        IF(TimeLag >= 0) EventScatter(IntTime,TimeLag,1) = EventScatter(IntTime,TimeLag,1) + 1
        ! qqqqq
        ! IF (ABS(TimeLag) <= 5) WRITE(99,'(I11)') PulseHist(I,1)
        IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 1 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 1) THEN
          CrossCorr(TimeLag,StopDetRow,2,1) = CrossCorr(TimeLag,StopDetRow,2,1) + 1
          IF (Pulses(PulseHist(I,2)+J-1)%GenNo == 0 .AND. &
          & Pulses(PulseHist(I,2)+K-1)%GenNo == 0) THEN
            CrossCorr(TimeLag,StopDetRow,7,1) = CrossCorr(TimeLag,StopDetRow,7,1) + 1
          ELSE IF (Pulses(PulseHist(I,2)+J-1)%GenNo > 0 .AND. &
          & Pulses(PulseHist(I,2)+K-1)%GenNo > 0) THEN
            CrossCorr(TimeLag,StopDetRow,8,1) = CrossCorr(TimeLag,StopDetRow,8,1) + 1
          END IF
        ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 2 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 2) THEN
          CrossCorr(TimeLag,StopDetRow,3,1) = CrossCorr(TimeLag,StopDetRow,3,1) + 1
        ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 1 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 2) THEN
          CrossCorr(TimeLag,StopDetRow,4,1) = CrossCorr(TimeLag,StopDetRow,4,1) + 1
        ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 2 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 1) THEN
          CrossCorr(TimeLag,StopDetRow,5,1) = CrossCorr(TimeLag,StopDetRow,5,1) + 1
        END IF
        IF (Pulses(PulseHist(I,2)+J-1)%HGamma == 1 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%HGamma == 1) THEN
          CrossCorr(TimeLag,StopDetRow,6,1) = CrossCorr(TimeLag,StopDetRow,6,1) + 1
        END IF
       END DO ! K Loop (Cross-Correlations)
    END IF
   ELSEIF (StartDet < 0) THEN
      DO K = 1, PulseHist(I,3)
        IF (K == J) CYCLE
        StopDetRow = Pulses(PulseHist(I,2)+K-1)%DetNo

        ! Version 5 - because of the change in the time bin structure so that the bins
        ! no longer (necessarily) match the times attached, TimeLag is now the number of
        ! bins between the two pulses rather than the time between pulses.
        TimeLag = NINT((Pulses(PulseHist(I,2)+K-1)%Time - &
        & Pulses(PulseHist(I,2)+J-1)%Time) / TimeIncrement)

        CrossCorr(TimeLag,StopDetRow,1,CurrentDet) = CrossCorr(TimeLag,StopDetRow,1,CurrentDet) + 1
        IF(IntTime > NumTBins .OR. TimeLag > NumTBins) WRITE(*,*) MAX(IntTime,TimeLag)
        IF(TimeLag >= 0) EventScatter(IntTime,TimeLag,1) = EventScatter(IntTime,TimeLag,1) + 1
        ! qqqqq
        ! IF (ABS(TimeLag) <= 5) WRITE(99,'(I11)') PulseHist(I,1)
        IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 1 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 1) THEN
          CrossCorr(TimeLag,StopDetRow,2,CurrentDet) = &
          & CrossCorr(TimeLag,StopDetRow,2,CurrentDet) + 1
          IF (Pulses(PulseHist(I,2)+J-1)%GenNo == 0 .AND. &
          & Pulses(PulseHist(I,2)+K-1)%GenNo == 0) THEN
            CrossCorr(TimeLag,StopDetRow,7,CurrentDet) = &
            & CrossCorr(TimeLag,StopDetRow,7,CurrentDet) + 1
          ELSE IF (Pulses(PulseHist(I,2)+J-1)%GenNo > 0 .AND. &
          & Pulses(PulseHist(I,2)+K-1)%GenNo > 0) THEN
            CrossCorr(TimeLag,StopDetRow,8,CurrentDet) = &
            & CrossCorr(TimeLag,StopDetRow,8,CurrentDet) + 1
          END IF
        ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 2 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 2) THEN
          CrossCorr(TimeLag,StopDetRow,3,CurrentDet) = &
          & CrossCorr(TimeLag,StopDetRow,3,CurrentDet) + 1
        ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 1 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 2) THEN
          CrossCorr(TimeLag,StopDetRow,4,CurrentDet) = &
          & CrossCorr(TimeLag,StopDetRow,4,CurrentDet) + 1
        ELSE IF (Pulses(PulseHist(I,2)+J-1)%ParNo == 2 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%ParNo == 1) THEN
          CrossCorr(TimeLag,StopDetRow,5,CurrentDet) = &
          & CrossCorr(TimeLag,StopDetRow,5,CurrentDet) + 1
        END IF
        IF (Pulses(PulseHist(I,2)+J-1)%HGamma == 1 .AND. &
        & Pulses(PulseHist(I,2)+K-1)%HGamma == 1) THEN
          CrossCorr(TimeLag,StopDetRow,6,CurrentDet) = &
          & CrossCorr(TimeLag,StopDetRow,6,CurrentDet) + 1
        END IF
       END DO ! K Loop (Cross-Correlations)
   END IF
  END DO
END DO

! qqqqq
! CLOSE(UNIT=99)

RETURN
END SUBROUTINE ComputeCorrelations



SUBROUTINE FindPeaks(ResetFlag, TerminateFlag)
LOGICAL,INTENT(IN) :: ResetFlag, TerminateFlag

!==========================================================================================
! Step 6 : Determine the Fast Neutron Time Window.  This is determined by finding the
! largest source-detector correlations in the entire measurement.  All detectors have the
! same time window.
!==========================================================================================
  PeakOld = 0
  PeakStart = 0
  PeakSumTotal = 0


! Version 5 - !!!WARNING!!! - peak value is now calculated in TIME BINS rather than ns.
WRITE(*,*) "Calculating Peak Values"

DO I = 1, NumDet
  DO J = 0, NumTBins
    PeakSumTotal = Correlation(J,I,1)
    IF (PeakSumTotal >= PeakOld) THEN
      PeakOld = PeakSumTotal
      PeakStart = J
    END IF
  END DO
END DO
  PeakStart = MAX(PeakStart - PeakWidth / 2, 0)

IF(TerminateFlag) THEN
  WRITE(*,*) "Fast Neutron Time Window: ", REAL(PeakStart)*TimeIncrement, " to ", &
  & REAL(PeakStart+Peakwidth-1)*TimeIncrement
END IF

RETURN
END SUBROUTINE FindPeaks

SUBROUTINE ComputeAntiC(ResetFlag, TerminateFlag)
LOGICAL,INTENT(IN) :: ResetFlag, TerminateFlag
REAL(KIND=4) :: ZZ

!==========================================================================================
! Step 7 : Now that the fast neutron time window is known, mark the cross talk between
! detectors.  Pulses are sorted by detector and then by time, so if more than one 
! correlation occurs in the fast time window for a given history, the one with the larger
! time is marked as cross-talk regardless of detector cell.
!==========================================================================================

WRITE(*,*) "Calculating Cross-Talk"

DO I = 1, NumHistP
  IF (PulseHist(I,3) == 1) THEN
! Revised in V6. This whole idea just seemed silly. If there's only 1 count, it's not cross-talk.
     ! This marks pulses which are the only pulse in the history .AND. outside of the
     ! peak window as cross-talk so that
     ! they won't appear in the NoXTalk values when they are extracted.
!     IF (Pulses(PulseHist(I,2))%Time/TimeIncrement < REAL(PeakStart) .OR. &
!     & Pulses(PulseHist(I,2))%Time/TimeIncrement >= REAL(PeakStart + PeakWidth)) THEN
!        Pulses(PulseHist(I,2))%XTalk = 1
!     END IF
     CYCLE
  END IF
  DO J = 1, PulseHist(I,3)
    ! Cycle if Pulse J not in time window or if it has already been marked as XTalk
    ! Note that the window is 1 smaller in integer math than in Real math because of 
    ! rounding.
    IF (Pulses(PulseHist(I,2)+J-1)%Time/TimeIncrement < REAL(PeakStart) .OR. &
    & Pulses(PulseHist(I,2)+J-1)%Time/TimeIncrement >= REAL(PeakStart + PeakWidth)) THEN
      Pulses(PulseHist(I,2)+J-1)%XTalk = 1
      CYCLE
    END IF
    IF (Pulses(PulseHist(I,2)+J-1)%XTalk == 1) CYCLE
    DO K = J+1, PulseHist(I,3)
      IF (Pulses(PulseHist(I,2)+K-1)%Time/TimeIncrement < REAL(PeakStart) .OR. &
      & Pulses(PulseHist(I,2)+K-1)%Time/TimeIncrement >= REAL(PeakStart+PeakWidth)) CYCLE
      ! If we have reached this point, both Pulse J and K are in the time window.  The
      ! one with the higher time gets marked a XTalk.  PoliMi only gives time to 0.1 ns,
      ! so it's possible (but very unlikely) the times could be equal.  If so, the 
      ! detector that gets counted as XTalk is chosen randomly.
      IF (Pulses(PulseHist(I,2)+J-1)%Time/TimeIncrement > Pulses(PulseHist(I,2)+K-1)%Time) THEN
        Pulses(PulseHist(I,2)+J-1)%XTalk = 1
      ELSEIF(Pulses(PulseHist(I,2)+J-1)%Time/TimeIncrement == Pulses(PulseHist(I,2)+K-1)%Time) THEN
        CALL RANDOM_NUMBER(ZZ)
        IF(ZZ < 0.5) THEN
          Pulses(PulseHist(I,2)+J-1)%XTalk = 1
        ELSE
          Pulses(PulseHist(I,2)+K-1)%XTalk = 1
        END IF
      ELSE
        Pulses(PulseHist(I,2)+K-1)%XTalk = 1
      END IF
    END DO
  END DO
END DO

! Unfortunately, in order to find the fast neutron window to remove the cross-talk, I had
! to compute correlations and now I have to go through again to compute the no cross-talk
! correlations.  These no cross-talk correlations should equal the directs plus scatter
! in the object being imaged inside of the fast neutron window.  If other features are in
! the geometry as well, such as the fancy detectors or the detector arm, scattered
! neutrons and induced gammas from those objects can contribute as well.  No cross-talk
! is computed outside of the peaks window, so those values should be exactly the same
! as the measured.

DO I = 1, NumHistP
  DO J = 1, PulseHist(I,3)
    ! IntTime is the time bin rather than (necessarily) the actual time.
    IntTime = INT(Pulses(PulseHist(I,2)+J-1)%Time/TimeIncrement)
    IF (ABS(IntTime) > NumTBins) CYCLE
    CurrentDet = Pulses(PulseHist(I,2)+J-1)%DetNo
    IF (Pulses(PulseHist(I,2)+J-1)%XTalk == 0) THEN
      Correlation(IntTime,CurrentDet,5) = Correlation(IntTime,CurrentDet,5) + 1
    END IF
  END DO
END DO

RETURN
END SUBROUTINE ComputeAntiC



SUBROUTINE ComputePeaks
!==========================================================================================
! Step 8 : Determine the number of counts in the fast neutron peaks.  These values are
! used for neutron radiography.
!==========================================================================================
ALLOCATE(Peaks(1:MaxDets))

PeakSumTotal = 0
PeakSumDirect = 0
PeakSumNoXT = 0

DO I = 1, NumDet
  Peaks(I)%DetNo = Detector(I)
  DO J = PeakStart, PeakStart + PeakWidth - 1
    PeakSumTotal = PeakSumTotal + Correlation(J,I,1)
    PeakSumDirect = PeakSumDirect + Correlation(J,I,4)
    PeakSumNoXT = PeakSumNoXT + Correlation(J,I,5)
  END DO

  Peaks(I)%TotalPeak = PeakSumTotal
  Peaks(I)%DirectPeak = PeakSumDirect
  Peaks(I)%XTPeak = PeakSumNoXT

  Peaks(I)%TotalMean = 0.
  Peaks(I)%DirectMean = 0.
  Peaks(I)%XTMean = 0.

  DO J = PeakStart, PeakStart + PeakWidth - 1

    IF (PeakSumTotal > 0) THEN
      Peaks(I)%TotalMean = Peaks(I)%TotalMean + REAL(J * Correlation(J,I,1))
    END IF

    IF (PeakSumDirect > 0) THEN
      Peaks(I)%DirectMean = Peaks(I)%DirectMean + REAL(J * Correlation(J,I,4))
    END IF

    IF (PeakSumNoXT > 0) THEN
      Peaks(I)%XTMean = Peaks(I)%XTMean + REAL(J * Correlation(J,I,5))
    END IF
  END DO

IF (PeakSumTotal > 0) THEN
  Peaks(I)%TotalMean = Peaks(I)%TotalMean / REAL(Peaks(I)%TotalPeak)
END IF
IF (PeakSumDirect > 0) THEN
  Peaks(I)%DirectMean = Peaks(I)%DirectMean / REAL(Peaks(I)%DirectPeak)
END IF
IF (PeakSumNoXT > 0) THEN
  Peaks(I)%XTMean = Peaks(I)%XTMean / REAL(Peaks(I)%XTPeak)
END IF

PeakSumTotal = 0
PeakSumDirect = 0
PeakSumNoXT = 0

END DO

RETURN
END SUBROUTINE ComputePeaks




SUBROUTINE WriteOutput
!==========================================================================================
! Step 9 : Print Output to screen and file(s).
!==========================================================================================

WRITE(*,*) " "
WRITE(*,*) "The .DAT file contains ", TotalLines," lines and ", TotalHist, " histories."
WRITE(*,*) " "
WRITE(*,*) "The largest history is number ", MaxEvtHist, ".  It has ", MaxEvt, " events."
WRITE(*,*) " "
WRITE(*,*) "The .DAT file has records for ", NumDet, " detectors."
WRITE(*,*) " "
WRITE(*,*) "A total of ", TotalPulses, " pulses were recorded."
WRITE(*,*) " "
WRITE(*,*) "The largest history is number ", MaxPulseHist, ".  It has ", MaxPulse, " pulses."
WRITE(*,*) " "
WRITE (*,*) TotalHGammas, " pulses were detected which correspond to the gamma energy of interest."
WRITE (*,*) "These pulses have an energy between ", HLineMin," and ", HLineMax, " MeV."
IF (CorrWindowOverflow > 0) THEN
  WRITE(*, '(I10, A)') CorrWindowOverflow, " pulses were discarded because they occured after &
  & the end of the correlation window.  Consider increasing the size of the window if these &
  & losses were not intended."
  WRITE(*,*) " "
END IF

IF (CorrWindowUnderflow > 0) THEN
  PRINT '(I10, A)', CorrWindowUnderflow, " pulses were discarded because they occurred before &
  & the start of the correlation window."
END IF

! Writes the correlation output to files.
OPEN(UNIT=11, FILE=TRIM(FileBase)//".total.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=12, FILE=TRIM(FileBase)//".neutron.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=13, FILE=TRIM(FileBase)//".gamma.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=14, FILE=TRIM(FileBase)//".direct.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=15, FILE=TRIM(FileBase)//".noXT.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=16, FILE=TRIM(FileBase)//".HGamma.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=17, FILE=TRIM(FileBase)//".gen0.corr", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=18, FILE=TRIM(FileBase)//".induced.corr", STATUS="REPLACE", ACTION="WRITE")

WRITE(11,'(I13,1000I10)') (Correlation(-1,J,1),J=0,NumDet)
WRITE(12,'(I13,1000I10)') (Correlation(-1,J,2),J=0,NumDet)
WRITE(13,'(I13,1000I10)') (Correlation(-1,J,3),J=0,NumDet)
WRITE(14,'(I13,1000I10)') (Correlation(-1,J,4),J=0,NumDet)
WRITE(15,'(I13,1000I10)') (Correlation(-1,J,5),J=0,NumDet)
WRITE(16,'(I13,1000I10)') (Correlation(-1,J,6),J=0,NumDet)
WRITE(17,'(I13,1000I10)') (Correlation(-1,J,7),J=0,NumDet)
WRITE(18,'(I13,1000I10)') (Correlation(-1,J,8),J=0,NumDet)

DO I = 0, NumTBins
  WRITE(11,'(ES13.5,1000I10)') REAL(Correlation(I,0,1))*TimeIncrement, &
  & (Correlation(I,J,1),J=1,NumDet)
  WRITE(12,'(ES13.5,1000I10)') REAL(Correlation(I,0,2))*TimeIncrement, &
  & (Correlation(I,J,2),J=1,NumDet)
  WRITE(13,'(ES13.5,1000I10)') REAL(Correlation(I,0,3))*TimeIncrement, &
  & (Correlation(I,J,3),J=1,NumDet)
  WRITE(14,'(ES13.5,1000I10)') REAL(Correlation(I,0,4))*TimeIncrement, &
  & (Correlation(I,J,4),J=1,NumDet)
  WRITE(15,'(ES13.5,1000I10)') REAL(Correlation(I,0,5))*TimeIncrement, &
  & (Correlation(I,J,5),J=1,NumDet)
  WRITE(16,'(ES13.5,1000I10)') REAL(Correlation(I,0,6))*TimeIncrement, &
  & (Correlation(I,J,6),J=1,NumDet)
  WRITE(17,'(ES13.5,1000I10)') REAL(Correlation(I,0,7))*TimeIncrement, &
  & (Correlation(I,J,7),J=1,NumDet)
  WRITE(18,'(ES13.5,1000I10)') REAL(Correlation(I,0,8))*TimeIncrement, &
  & (Correlation(I,J,8),J=1,NumDet)
END DO
CLOSE(UNIT=11)
CLOSE(UNIT=12)
CLOSE(UNIT=13)
CLOSE(UNIT=14)
CLOSE(UNIT=15)
CLOSE(UNIT=16)
CLOSE(UNIT=17)
CLOSE(UNIT=18)

IF (StartDet /= 0) THEN
 ! Writes the cross-correlation output to files.
 OPEN(UNIT=21, FILE=TRIM(FileBase)//".total.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=22, FILE=TRIM(FileBase)//".nn.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=23, FILE=TRIM(FileBase)//".pp.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=24, FILE=TRIM(FileBase)//".np.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=25, FILE=TRIM(FileBase)//".pn.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=26, FILE=TRIM(FileBase)//".HH.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=27, FILE=TRIM(FileBase)//".gen0.cc", STATUS="REPLACE", ACTION="WRITE")
 OPEN(UNIT=28, FILE=TRIM(FileBase)//".induced.cc", STATUS="REPLACE", ACTION="WRITE")
END IF

IF (StartDet > 0) THEN
 WRITE(21,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,1,1)
 WRITE(22,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,2,1)
 WRITE(23,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,3,1)
 WRITE(24,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,4,1)
 WRITE(25,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,5,1)
 WRITE(26,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,6,1)
 WRITE(27,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,7,1)
 WRITE(28,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,8,1)

 WRITE(21,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(22,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(23,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(24,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(25,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(26,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(27,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)
 WRITE(28,'(1000(I6,A,I3))') (StartDet,"-",CrossCorr(-NumTBins-1,J,1,1),J=1,NumDet)

 DO I = -NumTBins, NumTBins
  WRITE(21,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,1,1))*TimeIncrement, &
  & (CrossCorr(I,J,1,1),J=1,NumDet)
  WRITE(22,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,2,1))*TimeIncrement, &
  & (CrossCorr(I,J,2,1),J=1,NumDet)
  WRITE(23,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,3,1))*TimeIncrement, &
  & (CrossCorr(I,J,3,1),J=1,NumDet)
  WRITE(24,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,4,1))*TimeIncrement, &
  & (CrossCorr(I,J,4,1),J=1,NumDet)
  WRITE(25,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,5,1))*TimeIncrement, &
  & (CrossCorr(I,J,5,1),J=1,NumDet)
  WRITE(26,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,6,1))*TimeIncrement, &
  & (CrossCorr(I,J,6,1),J=1,NumDet)
  WRITE(27,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,7,1))*TimeIncrement, &
  & (CrossCorr(I,J,7,1),J=1,NumDet)
  WRITE(28,'(ES13.5,1000I10)') REAL(CrossCorr(I,0,8,1))*TimeIncrement, &
  & (CrossCorr(I,J,8,1),J=1,NumDet)
 END DO

ELSEIF (StartDet < 0) THEN
 WRITE(21,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,1,1)
 WRITE(22,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,2,1)
 WRITE(23,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,3,1)
 WRITE(24,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,4,1)
 WRITE(25,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,5,1)
 WRITE(26,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,6,1)
 WRITE(27,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,7,1)
 WRITE(28,'(I13)',ADVANCE="NO") CrossCorr(-NumTBins-1,0,8,1)

 DO K = 1, NumDet
  WRITE(21,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(22,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(23,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(24,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(25,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(26,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(27,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
  WRITE(28,'(1000(I6,A,I3))',ADVANCE="NO") (Detector(K),"-",CrossCorr(-NumTBins-1,J,1,K),J=1,NumDet)
 END DO

 WRITE(21,*)
 WRITE(22,*)
 WRITE(23,*)
 WRITE(24,*)
 WRITE(25,*)
 WRITE(26,*)
 WRITE(27,*)
 WRITE(28,*)

 DO I = -NumTBins, NumTBins
  WRITE(21,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,1,1))*TimeIncrement
  WRITE(22,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,2,1))*TimeIncrement
  WRITE(23,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,3,1))*TimeIncrement
  WRITE(24,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,4,1))*TimeIncrement
  WRITE(25,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,5,1))*TimeIncrement
  WRITE(26,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,6,1))*TimeIncrement
  WRITE(27,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,7,1))*TimeIncrement
  WRITE(28,'(ES13.5)',ADVANCE="NO") REAL(CrossCorr(I,0,8,1))*TimeIncrement

  DO K = 1, NumDet
   WRITE(21,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,1,K),J=1,NumDet)
   WRITE(22,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,2,K),J=1,NumDet)
   WRITE(23,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,3,K),J=1,NumDet)
   WRITE(24,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,4,K),J=1,NumDet)
   WRITE(25,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,5,K),J=1,NumDet)
   WRITE(26,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,6,K),J=1,NumDet)
   WRITE(27,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,7,K),J=1,NumDet)
   WRITE(28,'(1000I10)',ADVANCE="NO") (CrossCorr(I,J,8,K),J=1,NumDet)
  END DO

  WRITE(21,*)
  WRITE(22,*)
  WRITE(23,*)
  WRITE(24,*)
  WRITE(25,*)
  WRITE(26,*)
  WRITE(27,*)
  WRITE(28,*)

 END DO
END IF

IF (StartDet /= 0) THEN
 CLOSE(UNIT=21)
 CLOSE(UNIT=22)
 CLOSE(UNIT=23)
 CLOSE(UNIT=24)
 CLOSE(UNIT=25)
 CLOSE(UNIT=26)
 CLOSE(UNIT=27)
 CLOSE(UNIT=28)
END IF

! Writes the Multiplicity array to file.
OPEN(UNIT=31, FILE=TRIM(FileBase)//".multip", STATUS="REPLACE", ACTION="WRITE")
WRITE(31,'(A)') "  N      Total   Neutrons   H_Gammas  Gen. Zero    Induced"
DO I = 0, MaxPulse
  WRITE (31, '(I3, 5I11)') I, (Multiplicity(I,J), J=1,5)
END DO
CLOSE(UNIT=31)

OPEN(UNIT=41, FILE=TRIM(FileBase)//".peaks", STATUS="REPLACE", ACTION="WRITE")
  WRITE (41,'(I12, A)') NINT(REAL(nps) / SrcNorm), &
  & "  Total       Mean(T)     Direct      Mean(D)     No XTalk    Mean(N)"
DO I = 1, NumDet
  WRITE (41, '(2I12, ES12.4, I12, ES12.4, I12, ES12.4)') Peaks(I)
END DO
CLOSE(UNIT=41)

! Writes the pulse height spectra to file.
OPEN(UNIT=51, FILE=TRIM(FileBase)//".total.ph", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=52, FILE=TRIM(FileBase)//".neutron.ph", STATUS="REPLACE", ACTION="WRITE")
OPEN(UNIT=53, FILE=TRIM(FileBase)//".gamma.ph", STATUS="REPLACE", ACTION="WRITE")

DO I = 0, NumDet - 1
  WRITE(51,'(I11)', ADVANCE="NO") PHSpectrum(-1,I,1)
  WRITE(52,'(I11)', ADVANCE="NO") PHSpectrum(-1,I,2)
  WRITE(53,'(I11)', ADVANCE="NO") PHSpectrum(-1,I,3)
END DO
  WRITE(51,'(I11)') PHSpectrum(-1,NumDet,1)
  WRITE(52,'(I11)') PHSpectrum(-1,NumDet,2)
  WRITE(53,'(I11)') PHSpectrum(-1,NumDet,3)

DO I = 1, PHSNumBins
  WRITE(51,'(F11.5)', ADVANCE="NO") REAL(I) * PHSIncrement
  WRITE(52,'(F11.5)', ADVANCE="NO") REAL(I) * PHSIncrement
  WRITE(53,'(F11.5)', ADVANCE="NO") REAL(I) * PHSIncrement
  DO J = 1, NumDet - 1
  WRITE(51,'(I11)', ADVANCE="NO") PHSpectrum(I,J,1)
  WRITE(52,'(I11)', ADVANCE="NO") PHSpectrum(I,J,2)
  WRITE(53,'(I11)', ADVANCE="NO") PHSpectrum(I,J,3)
  END DO
  WRITE(51,'(I11)') PHSpectrum(I,NumDet,1)
  WRITE(52,'(I11)') PHSpectrum(I,NumDet,2)
  WRITE(53,'(I11)') PHSpectrum(I,NumDet,3)

END DO

CLOSE(UNIT=51)
CLOSE(UNIT=52)
CLOSE(UNIT=53)

IF(StartDet /= 0) THEN
  OPEN(UNIT=61, FILE=TRIM(FileBase)//".scatter.txt", STATUS="REPLACE", ACTION="WRITE")
  WRITE(61,'(A13,10000ES13.5)') "Time Bin     ", (REAL(I,4) / TimeIncrement, I = 0, NumTBins)
  DO I = 0, NumTBins
    WRITE(61,'(10000ES13.5)') REAL(I,4) / TimeIncrement, (EventScatter(J,I,1)/(REAL(nps,4) / SrcNorm),J=0,NumTBins)
  END DO
  CLOSE(UNIT=61)
END IF

RETURN
END SUBROUTINE WriteOutput

END PROGRAM