import ROOT, glob, math, sys, os
from ROOT import std


# Histograms
histos = []

nBins = 1000

zBins = 80
zMin = 0
zMax = 75

ClSBins = 100
ClAmpBins = 100
ClAmpMax = 100000

pullMax = 15 #10

errMax = 0.2

# Digi --------------------------------------------------------

DigiResXvsZ = ROOT.TH2D("DigiResXvsZ", "", zBins, zMin, zMax, nBins, -3,3)
histos.append(DigiResXvsZ)
DigiResYvsZ = ROOT.TH2D("DigiResYvsZ", "", zBins, zMin, zMax, nBins, -3,3)
histos.append(DigiResYvsZ)
DigiResZvsZ = ROOT.TH2D("DigiResZvsZ", "", zBins, zMin, zMax, nBins, -3,3)
histos.append(DigiResZvsZ)

DigiResXvsclAmp = ROOT.TH2D("DigiResXvsclAmp", "", ClAmpBins, 0, ClAmpMax, nBins, -2,2)
histos.append(DigiResXvsclAmp)
DigiResYvsclAmp = ROOT.TH2D("DigiResYvsclAmp", "", ClAmpBins, 0, ClAmpMax, nBins, -2,2)
histos.append(DigiResYvsclAmp)
DigiResZvsclAmp = ROOT.TH2D("DigiResZvsclAmp", "", ClAmpBins, 0, ClAmpMax, nBins, -2,2)
histos.append(DigiResZvsclAmp)


# Cluster ------------------------------------------------------

ResXvsZ = ROOT.TH2D("ResXvsZ", "", zBins, zMin, zMax, nBins, -2,2)
histos.append(ResXvsZ)
PullXvsZ = ROOT.TH2D("PullXvsZ", "", zBins, zMin, zMax, nBins, -1.*pullMax, pullMax)
histos.append(PullXvsZ)
ErrXvsZ = ROOT.TH2D("ErrXvsZ", "", zBins, zMin, zMax, nBins, 0, errMax)
histos.append(ErrXvsZ)

ResYvsZ = ROOT.TH2D("ResYvsZ", "", zBins, zMin, zMax, nBins, -2,2)
histos.append(ResYvsZ)
PullYvsZ = ROOT.TH2D("PullYvsZ", "", zBins, zMin, zMax, nBins, -1.*pullMax, pullMax)
histos.append(PullYvsZ)

ResZvsZ = ROOT.TH2D("ResZvsZ", "", zBins, zMin, zMax, nBins, -2,2)
histos.append(ResZvsZ)
PullZvsZ = ROOT.TH2D("PullZvsZ", "", zBins, zMin, zMax, nBins, -1.*pullMax, pullMax)
histos.append(PullZvsZ)

ClSize3D = ROOT.TH1D("ClSize3D", "", ClSBins, 0, ClSBins)
histos.append(ClSize3D)





nCluster = 10000
clSizeMPV = 15
clSizeSigma = 5
diffT = 0.026767 # cm/sqrt(cm)


for i in range(nCluster) :
    clSize = int(ROOT.gRandom.Landau(clSizeMPV, clSizeSigma))
    if clSize < 1 :
        continue
    clPos = 0
    clAmp = 0
    z = ROOT.gRandom.Uniform(zMin, zMax)
    for n in range(clSize) :
        digiAmp = 1
        digiPos = ROOT.gRandom.Gaus(0, diffT*math.sqrt(z))
        clPos += digiPos
        clAmp += digiAmp
    clPos /= clAmp
    
    ResXvsZ.Fill(z, clPos)
    ClSize3D.Fill(clSize)
    




c1=ROOT.TCanvas()
c1.Divide(3,3)

c1.cd(1)
ResXvsZ.Draw("colz")

c1.cd(2)
ResXvsZ.FitSlicesY(0,0,-1,0,"QNRL")
ROOT.gDirectory.Get("ResXvsZ_2").Draw()

c1.cd(3)
ClSize3D.Draw()

input()