import ROOT,sys,os,math
sys.path.append('/home/mberger/fopiroot/fopiroot_dev/macro/tpc/FOPI/mberger')
pandapath=os.environ.get('PANDAPATH')
sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1')
import argparse
from functions import set_palette


parser=argparse.ArgumentParser(description='plot samples and some infos from datafile')
parser.add_argument('recofile',help='the reco file',type=str)

args=parser.parse_args()

ROOT.gROOT.ProcessLine(".x rootlogon.C") 
ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")')
ROOT.gROOT.LoadMacro("stlPYROOT.h+")
set_palette()

infile=ROOT.TFile(args.recofile,'read')
tree=infile.Get('cbmsim')
if tree==None:
        print 'no tree found'
        infile.Close()
        exit()
        
hsample=ROOT.TH1D('hsample','a sample',40,-20,20)
canvas=ROOT.TCanvas('canvas','sample',1000,1000)        
ev=-1
u='w'
for e in tree:
    ev+=1
    itr=-1
    for tr in e.CutCosmics:
        itr+=1
        for icl in range(tr.nhits()):
            hid=tr.GetHitID(icl)
            sphit=e.TpcSPHit.At(hid)
            clid=sphit.getClusterID()
            cluster=e.TpcCluster.At(clid)
            digimap=cluster.getDigiMap()
            idi=-1
            for dig in digimap:
                idi+=1
                digi=e.TpcDigi.At(dig.first)
                samplemap=digi.getSampleMap()
                isam=-1
                variance_w=0
                variance=0
                variance_wrtMean=0
                nsam=0
                sumsamamp=0
                mean=0
                for sam in samplemap:
                    sample=e.TpcSample.At(sam.first)
                    if type(sample)==ROOT.TObject:
                            continue
                    mean+=sample.t()
                    nsam+=1
                mean/=nsam
                for sam in samplemap:
                    sample=e.TpcSample.At(sam.first)
                    if type(sample)==ROOT.TObject:
                            continue
                    hsample.Fill(sample.t()-digi.t(),sample.amp())
                    variance_w+=sample.amp() * ( sample.t()-digi.t() )**2
                    variance+=( sample.t()-digi.t() )**2
                    variance_wrtMean+=( sample.t()-mean )**2
                    sumsamamp+=sample.amp()
                    
                variance_w/=sumsamamp
                variance/=nsam
                variance_wrtMean/=nsam
                std_dev_w=math.sqrt(variance_w)
                std_dev=math.sqrt(variance)
                std_dev_wrtMean=math.sqrt(variance_wrtMean)
                print '#samples:{0}, sum Amp={1}, digiT={2}'.format(nsam,sumsamamp,digi.t())
                print 'var_w={0}, std_dev_w={1}'.format(variance_w,std_dev_w)
                print 'var={0}, std_dev={1}'.format(variance,std_dev)
                print 'var_wrtMean={0}, std_dev_wrtMean={1}, mean={2}'.format(variance_wrtMean,std_dev_wrtMean,mean)
                print ''
                canvas.cd()
                hsample.Draw()
                canvas.Update()
                u=raw_input('quit=q')
                if u=='q':
                    exit()
                hsample.Reset()