#!/usr/bin/env python2

## Example skeleton for validation check of the previously filled histograms 
## of G4 trajectories from IValidationAnalyzer.

import sys
import ROOT
import math

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

def OK(value,sigma, expected, delta):
    if 3.0 < abs(value-expected)/math.sqrt(delta*delta + sigma*sigma):
        return False
    return True

## The minimum entries for the histogram to be considered good.             
minEntries = 30

## The minimum mean start time of muons
minStartTime = 1
 
## The minimum mean stopping Z position
minStopZ = 6000
## The maximum mean stopping Z position
maxStopZ = 8000

## The minimum mean stopping X position
minStopX = 1500
## The maximum mean stopping X position
maxStopX = 2500

## The maximum mean stopping Y position
maxStopY = 150

## Momentum ranges
minStartPx = -50
maxStartPx = 50
minStartPz = 0
maxStartPz = 150
minStartP = 0
maxStartP = 200

file = "IcedustValidate-830_1.root"
#if sys.argv[1:]: #option to specify different input file as argument                            
#    file = sys.argv[1]

print file

input = ROOT.TFile(file,"OLD")

if not input.IsOpen(): 
    print "MISSING INPUT"
    sys.exit(0)

input.cd()
muonStartTime = input.Get("G4Trajectories_IValidationAnalyzer/hStartTime/hStartTime_13")
muonStartPx = input.Get("G4Trajectories_IValidationAnalyzer/hMomentumX/hMomentumX_13")
muonStartPz = input.Get("G4Trajectories_IValidationAnalyzer/hMomentumZ/hMomentumZ_13")
muonStartP = input.Get("G4Trajectories_IValidationAnalyzer/hMomentum/hMomentum_13")
muonStopZ = input.Get("G4Trajectories_IValidationAnalyzer/hPositionZ/hPositionZ_13")
muonStopY = input.Get("G4Trajectories_IValidationAnalyzer/hPositionY/hPositionY_13")
muonStopX = input.Get("G4Trajectories_IValidationAnalyzer/hPositionX/hPositionX_13")

## Save files of the plots   
muonStartTime.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStartTime.eps")
muonStopZ.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStopZ.eps")
muonStopY.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStopY.eps")
muonStopX.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStopX.eps")
muonStartPx.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStartPx.eps")
muonStartPz.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStartPz.eps")
muonStartP.Draw()
ROOT.gPad.Print("IcedustValidate-910MuonBeamValidationAnalysis-muonStartP.eps")

fail = False

#Check number of entries    
if muonStartTime.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStartTime.GetEntries())    

#Check time distribution    
if muonStartTime.GetMean() < minStartTime:
    fail = True
    print "Too low mean start time %f" % (muonStartTime.GetMean())  

#Check number of entries                                                                                                    
if muonStartPx.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStartPx.GetEntries())

#Check Px distribution                              
if muonStartPx.GetMean() < minStartPx or muonStartPx.GetMean() > maxStartPx:
    fail = True
    print "Mean start Px of muons out of expected range %f" % (muonStartPx.GetMean())

#Check number of entries              
if muonStartPz.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStartPz.GetEntries())

#Check Pz distribution                                                  
if muonStartPz.GetMean() < minStartPz or muonStartPz.GetMean() > maxStartPz:
    fail = True
    print "Mean start Pz of muons out of expected range %f" % (muonStartPz.GetMean())

#Check number of entries                    
if muonStartP.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStartP.GetEntries())

#Check P distribution           
if muonStartP.GetMean() < minStartP or muonStartP.GetMean() > maxStartP:
    fail = True
    print "Mean start P of muons out of expected range %f" % (muonStartP.GetMean())

#Check number of entries    
if muonStopZ.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStopZ.GetEntries())    

#Check Z distribution    
if muonStopZ.GetMean() < minStopZ or muonStopZ.GetMean() > maxStopZ :
    fail = True
    print "Mean Z position of muons out of expected range  %f" % (muonStopZ.GetMean()) 

#Check number of entries    
if muonStopX.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStopX.GetEntries())    

#Check X distribution    
if muonStopX.GetMean() < minStopX or muonStopX.GetMean() > maxStopX :
    fail = True
    print "Mean X position of muons out of expected range  %f" % (muonStopX.GetMean()) 

#Check number of entries    
if muonStopY.GetEntries() < minEntries:
    fail = True
    print "Too low number of entries %f" % (muonStopY.GetEntries())    

#Check Y distribution    
if abs(muonStopY.GetMean()) > maxStopY:
    fail = True
    print "Mean Y position of muons out of expected range %f" % (muonStopY.GetMean()) 

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