#! /usr/bin/env python

import re
import math

def grabValueFromAlpgenDict(textLine):
    #extract the number from line of the form
    #blah = 223 or blah = 223d0
    return (textLine.split('=')[1]).split('d')[0]

def add(x,y): return x+y

def alpgenValue(lines,alpgenLabel,defaultValue):
    bob=[x for x in lines if len(x.split())==4 if x.split()[3]==alpgenLabel]
    if len(bob)>0:
        return float(bob[0].split()[1])
    else:
        return defaultValue

def generateFome(bookkeepingfilename): #return a FOME dictionnary from the bookeeping file
    result={}
    result['process']='unspecified'
    result['decay']='unspecified'
    f=open(bookkeepingfilename,'r')
    lines=f.readlines()
    
    bob=[x for x in lines if re.compile("BeamEnergy").match(x)]
    result['root-s']=2*float(grabValueFromAlpgenDict(bob[0]))

    bob=[x for x in lines if re.compile("unwtd events").search(x)]
    evt=[int(x.split()[0]) for x in bob]
    result['nevts']=reduce(add,evt,0)

    bob=[x for x in lines if re.compile("NWarmup").match(x)]
    result['itmx1']=int(grabValueFromAlpgenDict(bob[0]))
    result['itmx2']=0

    bob=[x for x in lines if re.compile("NeventsWarmup").match(x)]
    result['ncal1']=int(grabValueFromAlpgenDict(bob[0]))

    bob=[x for x in lines if re.compile("NEventAfterWarmup").match(x)]
    result['ncal2']=int(grabValueFromAlpgenDict(bob[0]))

    result['ranseed1']=str([int(grabValueFromAlpgenDict(x)) for x in lines if re.compile("ranseed1").match(x)])
    result['ranseed2']=str([int(grabValueFromAlpgenDict(x)) for x in lines if re.compile("ranseed2").match(x)])
    result['ranseed3']=str([int(grabValueFromAlpgenDict(x)) for x in lines if re.compile("ranseed3").match(x)])
    result['ranseed4']=str([int(grabValueFromAlpgenDict(x)) for x in lines if re.compile("ranseed4").match(x)])


    bob=[x for x in lines if re.compile("Crosssection").search(x)]
    xsect=[float(x.split()[0]) for x in bob if float(x.split()[1])!=0 ]
    dxsect=[float(x.split()[1]) for x in bob if float(x.split()[1])!=0]
    xsectweight=[1/(x*x) for x in dxsect]
    totweight=reduce(add,xsectweight,0.0)
    tmpxsect=[xsect[i]*xsectweight[i] for i in range(len(xsect))]
    finalxsect=reduce(add,tmpxsect,0.0)/totweight
    finaldxsect=1/math.sqrt(totweight)
    #print bob
    #print xsect
    #print dxsect
    #print xsectweight,'sum =',totweight
    #print tmpxsect
    #print finalxsect,' +- ',finaldxsect,' (pb)'
    result['comments']=''
    bob=[int(x.split()[1]) for x in lines if x.split()[0]=="pythia_generated="]
    if len(bob)==0:
        ngenerated=1
    else:
        ngenerated=bob[len(bob)-1]
        result['comments']=result['comments']+' pythia_generated='+str(ngenerated)
    bob=[int(x.split()[1]) for x in lines if x.split()[0]=="pythia_tried="]
    if len(bob)==0:
        ntried=1
    else:
        ntried=bob[len(bob)-1]
        result['comments']=result['comments']+' pythia_tried='+str(ntried)
    if ntried==0:
        ntried=1
        ngenerated=1
    #I note  S=unweighted alpgen efficiency
    #       dS=unweighted alpgen efficiency error
    #       Ng=number of matched (generated) parton shower events
    #       Nt=number of tried events by parton shower
    # matching efficiency is e=Ng/Nt
    # The variance on the number of matched events is given by sqrt(Nt*e*(1-e))
    # This results in x-sections variance = variance/int_lum
    # where int_lum=Nt/S
    # so the variance on the x-section due to matching efficiency variance is
    # S*sqrt(e(1-e)/Nt)
    # The total error is the quadratic sum of the x-section error in
    # alpgen unweigted cross-section and matching efficiency
    # so error = sqrt(dS**2+S**2*e(1-e)/Nt)
    matcheff=float(ngenerated)/float(ntried)
    result['xsect-lo']=finalxsect*matcheff
    result['d-xsect-lo']=math.sqrt(finaldxsect**2+(finalxsect**2)*matcheff*(1-matcheff)/(ntried))
    result['xsect-nlo']=-1.0
    result['d-xsect-nlo']=-1.0

    result['po']='LO'
    result['gupi-wgt']=3
    result['mx-wgt-f']=1.0
    
    bob=[x for x in lines if re.compile("matchingMode").match(x)]
    if len(bob) > 0:
        result['match-algo']='MLM'
        result['in-ex-cl']=bob[0].split('=')[1].split()[0]
        marc=[x for x in lines if re.compile("EtClus").match(x)]
        if len(marc) > 0:
            result['matching_pt']=float(grabValueFromAlpgenDict(marc[0]))
        else:
            result['matching_pt']=alpgenValue(lines,'ptjmin',-1.0)
        marc=[x for x in lines if re.compile("RClus").match(x)]
        if len(marc) > 0:
            result['matching_dr']=float(grabValueFromAlpgenDict(marc[0]))
        else:
            result['matching_dr']=alpgenValue(lines,'drjmin',-1.0)
        cluopt=int(alpgenValue(lines,'cluopt',-1))
        if cluopt == 1:
            result['matching_opt']='pt'
        elif cluopt == 2:
            result['matching_opt']='mt'
        else:
            result['matching_opt']=str(cluopt)
        #algo = JCM2 if delta R =0.2, JCM7 if delta R=0.7
        tendr=10*result['matching_dr']
        algodr=int(tendr)
        if abs(tendr-algodr) <0.01 :
            result['jet-algo']='JCM'+str(algodr)
        else:
            result['jet-algo']='JCM dr='+str(result['matching_dr'])
    else:
        result['match-algo']='NONE'
        result['in-ex-cl']='NONE'
        result['matching_pt']=-1.0
        result['matching_dr']=-1.0
        result['matching_opt']='NONE'
        result['jet-algo']='NONE'

    result['part-kt-max']=''
    result['part-kt-min']=''


    nparton=[int(x.split()[2]) for x in lines if x.split()[0]=='LightPartons' if x.split()[1]=='=']
    nparton=nparton+[int(x.split()[2]) for x in lines if x.split()[0]=='c' if x.split()[1]=='=']
    nparton=nparton+[int(x.split()[2]) for x in lines if x.split()[0]=='b' if x.split()[1]=='=']
    nparton=nparton+[int(x.split()[2]) for x in lines if x.split()[0]=='t' if x.split()[1]=='=']
    result['npartons']=reduce(add,nparton,0)


    bob=[x for x in lines if re.compile("pdfset").match(x)]
    result['pdfcode']=bob[0].split()[2]
    result['pdfsource']='NATIVE'

    result['rscal-opt']=int(alpgenValue(lines,'iqopt',-1))
    result['fscal-opt']=result['rscal-opt']
    result['rscal-val']=alpgenValue(lines,'qfac',-1)
    result['fscal-val']=result['rscal-val']
    result['fact-scheme']='MSbar'
    result['ewopt']=3

    
    result['mt']=alpgenValue(lines,'mt',174.3) #174.3 is alpgen default
    result['mb']=alpgenValue(lines,'mb',4.7)
    result['mc']=alpgenValue(lines,'mc',0.0)
    result['mhl']=alpgenValue(lines,'mh',-1.0)
    result['mha']=-1.0
    result['mhp']=-1.0
    result['mhh']=-1.0


    bob=[x for x in lines if re.compile("AlpgenExec").match(x)]
    if bob[0].split('=')[1].split()[0]=='vbjetgen':
        #result['z-w']=True
        result['z-w']=1
    else:
        #result['z-w']=False
        result['z-w']=0


    ihvy=int(alpgenValue(lines,'ihvy',-1))
    if ihvy==4:
        result['pt_min_hf']=alpgenValue(lines,'ptcmin',-1.0)
        result['y_max_hf']=alpgenValue(lines,'etacmax',-1000.0)
        result['y_min_hf']=-alpgenValue(lines,'etacmax',-1000.0)
        result['dr-hf-hf']=alpgenValue(lines,'drcmin',-1.0)
    elif ihvy==5:
        result['pt_min_hf']=alpgenValue(lines,'ptbmin',-1.0)
        result['y_max_hf']=alpgenValue(lines,'etabmax',-1000.0)
        result['y_min_hf']=-alpgenValue(lines,'etabmax',-1000.0)
        result['dr-hf-hf']=alpgenValue(lines,'drbmin',-1.0)
    else:
        result['pt_min_hf']=-1.0
        result['y_max_hf']=-1000.0
        result['y_min_hf']=1000.0
        result['dr-hf-hf']=-1.0
    result['pt_min_lf']=alpgenValue(lines,'ptjmin',-1.0)
    result['y_max_lf']=alpgenValue(lines,'etajmax',-1000.0)
    result['y_min_lf']=-alpgenValue(lines,'etajmax',-1000.0)


    result['mass-max-34']=alpgenValue(lines,'mllmax',1960.0)
    result['mass-min-34']=alpgenValue(lines,'mllmin',0.0)
    result['mass-max-56']=1960.0
    result['mass-min-56']=0.0

    result['dr-hf-lf']=alpgenValue(lines,'drjmin',-1.0)
    result['dr-lf-lf']=alpgenValue(lines,'drjmin',-1.0)
    result['dr-lep-hf']=alpgenValue(lines,'drlmin',-1.0)
    result['dr-lep-lep']=alpgenValue(lines,'drlmin',-1.0)

    result['met']=alpgenValue(lines,'metmin',0.0)
    result['y-max-lep1']=alpgenValue(lines,'etalmax',-1000.0)
    result['y-max-lep2']=result['y-max-lep1']
    result['y-min-lep1']=-alpgenValue(lines,'etalmax',-1000.0)
    result['y-min-lep2']=result['y-min-lep1']
    result['r-lep-cone']=-1.0

    result['pt-min-gam1']=alpgenValue(lines,'ptphmin',-1.0)
    result['pt-min-gam2']=result['pt-min-gam1']
    result['y-max-gam1']=alpgenValue(lines,'etaphmax',-1000.0)
    result['y-max-gam2']=result['y-max-gam1']
    result['y-min-gam1']=-alpgenValue(lines,'etaphmax',-1000.0)
    result['y-min-gam2']=result['y-min-gam1']
    result['r-gam-cone']=-1.0
    result['dr-gam-lep']=alpgenValue(lines,'drphlmin',-1.0)
    result['dr-gam-lf']=alpgenValue(lines,'drphjmin',-1.0)
    result['ptfr-gam-cone']=0.0

    
    cosvma=alpgenValue(lines,'cosvma',-1000.0)
    if (cosvma < -500):
        #result['mod-couplings']=False
        result['mod-couplings']=0
        result['mod-coupl-val']=''
    else:
        #result['mod-couplings']=True
        result['mod-couplings']=1
        result['mod-coupl-val']='top-W coupling: cosvma*(V-A)+sinvma*(V+A) with cosvma='+str(cosvma)

    irap=int(alpgenValue(lines,'irapgap',0))
    if irap == 1:
        result['comments']=result['comments']+' etagap='+str(alpgenValue(lines,'etagap',-1))+' ptcen='+str(alpgenValue(lines,'ptcen',-1))

    return result


if __name__ == "__main__":
    print repr(generateFome('alpgen.bookkeeping'))