#PlotMuonsDepth.py
#
#PURPOSE:
#  To plot the fraction of charge lost vs the depth of the
#  muon in the detector when it ionized the atoms
#
from Return_CR_DbQuantities import *
from Djs_Iterstat import *
from Linefit_func import *
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
execfile('/afs/slac/u/ki/lances/python/ThesisPlotSettings.py')

PlotDir = '/nfs/slac/g/ki/ki04/lances/LSST/KPNO/Latex/Thesis/CosmicRayTraps/Figures/'
TvFlag  = 2
SaveImg = 1

#DBKeyStr   = '_4_10_09'
DBKeyStr   = ''
TextSize   = 14
XExtentMin = 15
YExtentMin = 15
CosEDepMin = 300
CosEDepMax = 10000
Temp       = 160
VSUB       = 20

if Temp == 150:
  QueryStr = ' AND (EXTENT  > ' +str(XExtentMin)+ ' OR ' + \
             '      YEXTENT > ' +str(YExtentMin)+ ') AND '+\
             ' MAXCOSEDEP > '   +str(CosEDepMin)+ ' AND ' +\
             ' MAXCOSEDEP < '   +str(CosEDepMax)+ ' AND ' +\
             ' LMBDA2/LMBDA1 < 0.1 AND CLASS = "Muon" AND '+\
             ' COSREAD > 20 AND COSREAD < 60 AND '+\
             ' COSID REGEXP "'+str(VSUB)+'.00V" '

elif Temp == 160:
  QueryStr = ' AND (EXTENT  > ' +str(XExtentMin)+ ' OR ' + \
             '      YEXTENT > ' +str(YExtentMin)+ ') AND '+\
             ' MAXCOSEDEP > '   +str(CosEDepMin)+ ' AND ' +\
             ' MAXCOSEDEP < '   +str(CosEDepMax)+ ' AND ' +\
             ' LMBDA2/LMBDA1 < 0.1 AND CLASS = "Muon" AND '+\
             ' COSREAD > 20 AND COSREAD < 80 AND '+\
             ' COSID REGEXP "'+str(VSUB)+'.00V" ' 
else:
  QueryStr = ' AND (EXTENT  > ' +str(XExtentMin)+ ' OR ' + \
             '      YEXTENT > ' +str(YExtentMin)+ ') AND '+\
             ' MAXCOSEDEP > '   +str(CosEDepMin)+ ' AND ' +\
             ' MAXCOSEDEP < '   +str(CosEDepMax)+ ' AND ' +\
             ' LMBDA2/LMBDA1 < 0.1 AND CLASS = "Muon" AND '+\
             ' COSREAD > 20 AND COSREAD < 80 AND '+\
             ' COSID REGEXP "'+str(VSUB)+'.00V" '


CosIDs, MaxCosEDeps, MaxQLosts, MaxQRises, TotCosEDeps, TotQLosts, TotQRises=\
Return_CR_DbQuantities('H2RG-001', 'ASIC', QueryStr, 2, Temp=Temp, \
                       DBKeyStr = DBKeyStr)

mplot.figure(0)
mplot.figure(1)

DepthsLostCum   = []
QLostsCum       = []
CosEDepsLostCum = []
DepthsRiseCum   = []
QRisesCum       = []
CosEDepsRiseCum = []

for CosID in CosIDs :
  QueryStr = ' AND COSID = "' + CosID +'"'+\
             ' AND (MINREADAC-COSREAD) > 0 ' +\
             ' AND (MINREADAC-COSREAD) > 0 ' +\
             ' AND COSENEDEP > 500 AND VSUB = '+str(VSUB)
  CosReads, CosEDeps, QLosts, QRises, Depths, MaxDarks = \
  Return_CR_DbQuantities('H2RG-001', 'ASIC', QueryStr, 3, Temp=Temp, \
                         DBKeyStr = DBKeyStr)
  if size(where(QLosts > 0)) > 0:
    if size(where(QRises > 0)) > 0:
      DepthsFit   = append(Depths[QLosts > 0], Depths[QRises > 0])
      QLostsFit   = append(QLosts[QLosts > 0], -QRises[QRises > 0])
      CosEDepsFit = append(CosEDeps[QLosts > 0], CosEDeps[QRises > 0])
    else:
      DepthsFit   = Depths[QLosts > 0] 
      QLostsFit   = QLosts[QLosts > 0]
      CosEDepsFit = CosEDeps[QLosts > 0]
    m, b, [Err_b, Err_m], CosRead = Linefit(DepthsFit, \
                                  QLostsFit/CosEDepsFit,\
                                  RejectCos=0)
  else:
    m = 0
  #Adjust for cases where the cosmic was not fit in the right direction
  if m < 0 :
    Depths = 100-Depths
  DepthsLostCum   = append(DepthsLostCum,   Depths[QLosts > 0])
  QLostsCum       = append(QLostsCum,       QLosts[QLosts > 0])
  CosEDepsLostCum = append(CosEDepsLostCum, CosEDeps[QLosts > 0])
  DepthsRiseCum   = append(DepthsRiseCum,   Depths[QRises > 0])
  QRisesCum       = append(QRisesCum,       QRises[QRises > 0])
  CosEDepsRiseCum = append(CosEDepsRiseCum, CosEDeps[QRises > 0])

#Try to bin by mean
NumBins = 10
BinDis  = 100/NumBins
ZDis    = arange(NumBins)*BinDis
ZAve    = zeros(NumBins, dtype=float32)
for i, Z in zip(arange(NumBins-1), ZDis):
  print i
  QLostZs  = QLostsCum[(where((DepthsLostCum > ZDis[i]) & \
                             (DepthsLostCum < ZDis[i+1])))[0]]
  QDepsLZs = CosEDepsLostCum[(where((DepthsLostCum > ZDis[i]) & \
                                   (DepthsLostCum < ZDis[i+1])))[0]]
  QRiseZs  = QRisesCum[(where((DepthsRiseCum > ZDis[i]) & \
                             (DepthsRiseCum < ZDis[i+1])))[0]]
  QDepsRZs = CosEDepsRiseCum[(where((DepthsRiseCum > ZDis[i]) & \
                                   (DepthsRiseCum < ZDis[i+1])))[0]]

  QRatios = append(QLostZs/QDepsLZs, QRiseZs/QDepsRZs)
  print mean(QRatios) 
  ZAve[i] = mean(QRatios)

  if TvFlag == 1:
    mplot.figure(0)
    mplot.clf()
    mplot.plot(Depths, QLosts/CosEDeps,'ro')
    mplot.figure(1)
    mplot.clf()
    mplot.plot(Depths, QRises, 'bo')
    mplot.figure(2)
    mplot.clf()
    mplot.plot(CosEDeps, QLosts, 'go')
    pdb.set_trace() 
    ylim(0,.5)

if TvFlag == 2:
  mplot.figure(1)
  mplot.clf()
  mplot.plot(DepthsRiseCum, QRisesCum, 'bo')
  ylim(0,1000)
  mplot.figure(2)
  mplot.clf()

  mplot.plot(CosEDepsLostCum, QLostsCum, 'go')
  mplot.figure(0)
  mplot.clf()
  mplot.plot(1.05*DepthsLostCum, QLostsCum/CosEDepsLostCum,'ro')
  mplot.plot(1.05*DepthsRiseCum, -QRisesCum/CosEDepsRiseCum, 'ro')
  xlim(0,100)
  ylim(-.1,.4)
  mplot.xlabel(r'\textbf{Distance from Front Side ($\mu m$)}')
  mplot.ylabel(r'$Q_{lost}/Q_{dep}$')
  mplot.title(r'Muons with Horizontal Track Length $\ge$ ' + str(18*XExtentMin) + \
              ' $\mu $m')
  #Figure Text
  mplot.figtext(.55,.30,r'\textbf{\# Muons \ : ' + str(size(CosIDs))+'} ',\
                  fontsize=TextSize)
  TempStr = "%.1f" % (round(float(Temp)/10)*10)
  mplot.figtext(.58,.23,r'\textbf{Temp \ \ : } ',\
                  fontsize=TextSize)
  mplot.figtext(.71,.24,r'\textbf{'+str(TempStr) +'\ K}', \
                  fontsize=TextSize, color='k')
  VSUBStr = "%5.1f" % (round(float(VSUB)))
  mplot.figtext(.58,.16, r'\textbf{V$_{SUB}$\ \ \ \ \ : '+ VSUBStr+'\ V}',
                  fontsize=TextSize)


if SaveImg == 1:
  mplot.savefig(PlotDir+'QLostVsDepth_'+str(Temp)+'K_'+str(VSUB)+'V.eps')
  mplot.savefig(PlotDir+'QLostVsDepth_'+str(Temp)+'K_'+str(VSUB)+'V.png')