#!/usr/bin/env python2

## Using a previously generated profiles of the collected charge vs
## distance from the MPPC to check the attenuation effect produced by
## SimDetectorResponse.

## The minimum probability for the fit to be considered good.
minProbability = 0.01

## The expected fraction of light going to the long wave length component
longComponent = 0.77
longComponentSigma = 0.1

## The expected "reflectivity".  This is modified from the nominal
## value since the bar simulation included non-linearity.
reflectivity = 0.9
reflectivitySigma = 0.1

import sys
import ROOT
import math

ROOT.gStyle.SetOptStat("")
ROOT.gStyle.SetOptFit(1111)

def OK(value,sigma, expected, delta):
    if 1.0 < abs(value-expected)/math.sqrt(delta*delta + sigma*sigma):
        return False
    return True
    
input = ROOT.TFile("SimDetectorResponse-800GPSMuon5GeVThroughP0D-p0dAtten.root","OLD")
if not input.IsOpen(): 
    print "MISSING INPUT"
    sys.exit(0)

p0dXAtten = input.Get("p0dXAtten")
p0dYAtten = input.Get("p0dYAtten")

# Functions to fit the bar attenuation.  The short (332mm) and long (4634mm)
# attenuation lengths are taken from the URegina measurements of the
# fiber attenuation.
p0dXFunc = ROOT.TF1("p0dXFunc","[0]*((1-[1])*exp(-x/332)+[1]*(exp(-x/4634)+[2]*exp((-2*2272+x)/4634)))")
p0dXFunc.SetParameter(0,30)
p0dXFunc.SetParameter(1,0.77)
p0dXFunc.SetParameter(2,0.80)

p0dYFunc = ROOT.TF1("p0dYFunc","[0]*((1-[1])*exp(-x/332)+[1]*(exp(-x/4634)+[2]*exp((-2*2133+x)/4634)))")
p0dYFunc.SetParameter(0,30)
p0dYFunc.SetParameter(1,0.77)
p0dYFunc.SetParameter(2,0.80)

p0dXAtten.Fit("p0dXFunc")
print p0dXFunc.GetParameter(0)
print p0dXFunc.GetParameter(1)
print p0dXFunc.GetParameter(2)
print p0dXFunc.GetProb()
p0dXAtten.SetMinimum(15)
p0dXAtten.Draw()
ROOT.gPad.Print("SimDetectorResponse-900GPSMuon5GeVThroughP0D-p0dXFunc.eps")

p0dYAtten.Fit("p0dYFunc")
print p0dYFunc.GetParameter(0)
print p0dYFunc.GetParameter(1)
print p0dYFunc.GetParameter(2)
print p0dYFunc.GetProb()
p0dYAtten.SetMinimum(15)
p0dYAtten.Draw()
ROOT.gPad.Print("SimDetectorResponse-900GPSMuon5GeVThroughP0D-p0dYFunc.eps")

fail = False
if not OK(p0dXFunc.GetParameter(1),p0dXFunc.GetParError(1), \
          longComponent, longComponentSigma):
    fail = True
    print "long component fraction wrong %f expected %f" \
        % (p0dXFunc.GetParameter(1), longComponent)
    
if not OK(p0dXFunc.GetParameter(2),p0dXFunc.GetParError(2), \
          reflectivity, reflectivitySigma):
    fail = True
    print "reflectivity wrong %f expected %f" \
        % (p0dXFunc.GetParameter(2), reflectivity)

if p0dXFunc.GetProb() < minProbability:
    fail = True
    print "bad X fit %f" % (p0dXFunc.GetProb())
    
if p0dYFunc.GetProb() < minProbability:
    fail = True
    print "bad Y fit %f" % (p0dYFunc.GetProb())

if fail: 
    print "FAIL"
    sys.exit(1)