#!/usr/bin/env python

#
#  Python verion of NP03 detector simulation program.
#
#  The files contain a series of events generated by the GEANT4 simulation program.
#  The simulation is an infinite block of either Silicon or Germanium.  Photons
#  are generated at the origin pointing along the positive x axis.  The simulation
#  generates the interactions in the following particle cascade.  We are going to
#  make two simplifying assumptions to avoid the computing becoming too complicted.
#  (1) The electron and positrons deposit their energy all in one location as soon
#      as they are generated.  This is almost true, they usually travel a few
#      tens of microns to deposit their energy.
#  (2) We neglect tertiary bremstrahlung - i.e. the process where electrons emit
#      photons (which move a large distance) before making electrons again.  This
#      is not bad either - in reality <1% of events have a bremstrahlung.  In the
#      files, the energy of the bremstrahlung photon is added in to the electron
#      which made it.
#  This means that the events in the files almost always contain a single chain of
#  electrons which have been knocked off by the primary photon.  The exception to
#  this is pair production.
#
#  There are four versions of the input files simulated with different parameters
#  the variations in their file names are given here in brackets.
#  Two of the files are simulated with Silicon (Si), and two with Germanium (Ge).
#  Two start with monoenergetic 'Cesium' photons E=0.6617MeV (Cs) and the other
#  two have photons with random energies uniformly from 0.030MeV to 3.030MeV (All)
#
#  The next clump of comments document the input file format in too much detail
#  to begin with.  You don't need to know all the detail to start, probably.  The
#  program loops over the events and reads them in.  Better start by looking at the
#  comments in the loop which gives the meaning of the variables which have been read.
#
#  Here is an example of an event.  It contains a series of 'E' lines followed by
#  one 'M' line which indicates the end of the event.  The 'E' lines represent a
#  particle within one interaction.  The 3 numbers in the 1st bracket are (x,y,z)
#  of the end of the first step the particle made (we think it is where it first
#  interacts?). The next 4 numbers are also in brackets and are not used (except
#  that the first 3 are a repeat of the x,y,z of the parent particle).  The final six
#  bits of information are partice-type particle-energy ignore-next partID parentID
#  interact-code.  The one to ignore is the energy loss inside the GEANT step.  This
#  is not useful without more info from GEANT because it can step the particle
#  more than once.
#  Energies are in MeV.  Distances are in mm (must check the mm).
#  The interact-code is compt=compton scatter, conv=pair production (aka conversion)
#  phot=photoelectric.  For the photon, it gives the first interaction which happened.
#  The interact-code of the electron is usually eIoni=Electron ionisation.
#
#E (103.595025,0.000000,0.000000) (0.000000,0.000000,0.000000,0.000000) gamma 0.661700 0.000000 1 0 compt
#E (110.054005,129.214988,-29.651278) (103.595025,0.000000,0.000000,103.595025) e- 0.052832 0.052832 11 1 eIoni
#E (104.175232,126.425472,-27.431369) (103.595025,0.000000,0.000000,103.595025) e- 0.002926 0.002926 10 1 eIoni
#E (104.103995,123.438379,-24.163798) (103.595025,0.000000,0.000000,103.595025) e- 0.014981 0.014981 9 1 eIoni
#E (104.992130,124.195952,-25.239895) (103.595025,0.000000,0.000000,103.595025) e- 0.028011 0.028011 8 1 eIoni
#E (100.363148,119.641707,-21.950251) (103.595025,0.000000,0.000000,103.595025) e- 0.058091 0.058091 7 1 eIoni
#E (104.220751,127.052374,-33.705091) (103.595025,0.000000,0.000000,103.595025) e- 0.140492 0.140492 6 1 eIoni
#E (114.121669,96.487046,-29.259889) (103.595025,0.000000,0.000000,103.595025) e- 0.018915 0.018915 5 1 eIoni
#E (143.836794,59.362828,-7.262470) (103.595025,0.000000,0.000000,103.595025) e- 0.054741 0.054741 4 1 eIoni
#E (143.944857,31.562242,-3.126021) (103.595025,0.000000,0.000000,103.595025) e- 0.145172 0.145172 3 1 eIoni
#E (103.624575,-0.036994,0.003667) (103.595025,0.000000,0.000000,103.595025) e- 0.143694 0.143694 2 1 eIoni
#M 0 0.661700 gamma compt 0.659855 0.997212 -0.001845

# The numbers on the 'M' line are: event-number first-particle-energy first-particle-type
# first-interaction-of-first-particle sum-of-all-secondary-energy.  The last two give the
# discrepancy between the first-particle-energy and the sum-of-all-secondary-energy,
# the first as a fraction, the second as a difference.  The difference seems to be
# the more consistent, perhaps there is some atomic binding energy or something.
# If there was a pair production, the rest mass energy of the e+e- is not included in
# the total.

import sys
import math

def main():
  # Check correct arguments
  if len(sys.argv) != 3:
      print "usage: %s <in-text-file> <out-text-file>" % sys.argv[0]
      sys.exit(1)

  infile = sys.argv[1]
  outfile = sys.argv[2]

  # Open input text file (this is printout from the MC with 'vent:' lines added
  if infile == "-": inf = sys.stdin
  else:             inf = open(infile,'rU')
  if not inf:
      print "Input file %d cannot be opened" % outfile
      sys.exit(1)

  # Open output file
  if outfile == "-": ouf = sys.stdout
  else:              ouf = open(outfile,'wU')
  if not ouf:
      print "Output file %d cannot be opened" % outfile
      sys.exit(1)

  # ***Place to modify*** Reserve space for histogramming
  NBINS = 300
  henerall  = [0]*NBINS    # Makes an array with 100 elements and puts zero in each
  henerpeak = [0]*NBINS
  htener    = [0]*NBINS

  # Main loop
  tener = 0.     # Variable for total energy detected
  while (1):    # Infinite loop, exit with break statement below
      line = inf.readline()    # read one line from file
      if line == "": break     # readline only returns empty at end of file
      line = line.rstrip('\n') # get rid of the carriage return character

      # An E line.   Unpack and add energy to total if wanted
      if line[0] == 'E':
         (text, Sxyz, Spxyzr, Spart, Sener, Sdener, Sid, Spid, Sinter) = line.split()
         (Sx,Sy,Sz) = Sxyz.lstrip('(').rstrip(')').split(',')   # Splits the first xyz bracket
         (Spx,Spy,Spz,Spr) = Spxyzr.lstrip('(').rstrip(')').split(',')   # Splits the second xyzr bracket

         # The above variables are all still character strings (start with capital S)
         # Convert them to floats and integers (use lowercase s)
         sx = float(Sx)   # x position of step (mm)
         sy = float(Sy)   # y position of step (mm)
         sz = float(Sz)   # z position of step (mm)
         spx = float(Spx) # x position of photon in it's step (mm)
         spy = float(Spy) # y position of photon in it's step (mm)
         spz = float(Spz) # z position of photon in it's step (mm)
         spart = Spart    # Particle name 'gamma', 'e-' or 'e+'
         sener = float(Sener) # Particle energy (MeV)
         sid = int(Sid)   # Particle ID (1=primary photon)
         spid = int(Spid) # Particle ID of parent
         sinter = Sinter  # Inteaction-code (see above)

         # ***Place to modify***  Decide whether to add in the energy
         if (Spart == "e-" or Spart == "e+"):  # Only e+ and e-, not photons
            if math.sqrt(sy*sy+sz*sz) < 30.:   # radius 3cm
               if sx > 0. and sx < 100.:       # 10cm long
                  tener = tener + sener

      # End of 'E' line processing  (python just ends indent to end if or loop!)

      # An 'M' line.  Process the event, and reset total for next
      if line[0] == 'M':
         (text,Mevno,Mener,Mpart,Minter,Menert,Mratio,Mdifference) = line.split()
         mevno = int(Mevno)      # Event number
         mener = float(Mener)    # Energy of photon
         mpart = Mpart           # Particle type of primary (should be 'photon')
         minter = Minter         # First interaction of primary usually 'conv', 'compt' or 'phot'
         menert = Menert         # Total energy of secondaries

         # ***Place to modify***  Histogram filling
         # Make an integer which is the bin number from 0 to 99 for E=0 to 3
         # Put events above 3MeV in bin 99
         ib = int(mener/3.0*float(NBINS))
         if ib < 0:  ib = 0
         if ib >= NBINS: ib = NBINS-1

         henerall[ib] = henerall[ib]+1
         d = tener - mener     # Difference in energy (nearly 0 if on peak)
         if d > -0.02 and d < 0.02: henerpeak[ib] = henerpeak[ib]+1

         # Similar bin calculation to above for tener
         ibt = int(tener/3.0*float(NBINS))
         if ibt < 0: ibt = 0
         if ibt >= NBINS: ibt = NBINS-1
         htener[ibt] = htener[ibt] + 1

         # Reset variables which are accumulted in the E lines for next event
         tener = 0.

         # Print occasionally so we know it is running
         if mevno%10000 == 0: print "Just done event ",mevno

      # End of 'M' line processing (python just ends indent to end if or loops)

  # End of loop.

  # Now write a table from histograms to output file
  for ib in range(NBINS):
     e = ib*3.0/float(NBINS)     # Invert the ib calculation above to get energy at bottom edge of bin
     print >>ouf,e,henerall[ib],henerpeak[ib],htener[ib]

  inf.close()
  ouf.close()

if __name__=='__main__':main()