upmns.py

import numpy as np
def neutrino_data():
    '''From arxiv:1405.7540 (table I)
    and asumming a Normal Hierarchy:
    Output:
    mnu1in: laightest neutrino mass
    Dms2: \Delta m^2_{12}
    Dma2: \Delta m^2_{13}
    ThetSol,ThetAtm,ThetRec: in radians
    '''
    Dms2=np.array([7.11e-5, 7.60e-5, 8.18e-5])*1e-18 # In GeV
    Dma2=np.array([2.30e-3, 2.48e-3, 2.65e-3])*1e-18 # In GeV
    #input real values:
    #
    ThetSol = np.array([0.278,  0.323,  0.375]) 
    ThetAtm = np.array([0.392,  0.567,  0.643])
    ThetRec = np.array([0.0177, 0.0234, 0.0294])
    mnu1in=1E-5*1E-9

    return mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec

def CasasIbarra(ranMnu=True,bestfit=True):
    #import numpy as np
    if ranMnu==True:
        bestfit=True
        
    mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec=neutrino_data()
    
    # Phases of the PMNS matrix
    phases1=np.random.uniform(0.,0.0*np.pi,3)       # cero value for all phases         
    
    delta=1.*(0 if ranMnu else phases1[0])
    eta1 =1.*(0 if ranMnu else phases1[1])
    eta2 =1.*(0 if ranMnu else phases1[2])

    mnu1=1.*(mnu1in      if bestfit else  10**((np.log10(2.5e-3)-np.log10(1e-9))*np.random.uniform(0,1)+np.log10(1e-9))*1e-9) 
    mnu2=1.*(np.sqrt(Dms2[1]+mnu1in**2)     if bestfit else  np.sqrt(np.random.uniform(Dms2[0],Dms2[2]) + mnu1**2)  ) 
    mnu3=1.*(np.sqrt(Dma2[1]+mnu1in**2)     if bestfit else  np.sqrt(np.random.uniform(Dma2[0],Dma2[2]) + mnu1**2)  )
    
    # Square root of left-handed neutrino mass matrix 
    DMnu = np.diag( [np.sqrt(mnu1),np.sqrt(mnu2),np.sqrt(mnu3)] )  

    #print "NEUTRINOS",mnu1,mnu2,mnu3

    # Mixing angles (up 3 sigma range)
    t12 = 1.*( np.arcsin(np.sqrt(ThetSol[1])) if bestfit else np.arcsin(np.sqrt(np.random.uniform(ThetSol[0],ThetSol[2]))))
    t23 = 1.*( np.arcsin(np.sqrt(ThetAtm[1])) if bestfit else np.arcsin(np.sqrt(np.random.uniform(ThetAtm[0],ThetAtm[2]))))
    t13 = 1.*( np.arcsin(np.sqrt(ThetRec[1])) if bestfit else np.arcsin(np.sqrt(np.random.uniform(ThetRec[0],ThetRec[2]))))         

    # Building PMNS matrix
    U12 = np.array([ [np.cos(t12), np.sin(t12),0], [-np.sin(t12), np.cos(t12),0], [0,0,1.0] ])
    U13 = np.array([ [np.cos(t13),0, np.sin(t13)* np.exp(-delta*1j)], [0,1.0,0], [-np.sin(t13)*np.exp(delta*1j),0, np.cos(t13)] ])
    U23 = np.array([ [1.0,0,0], [0, np.cos(t23), np.sin(t23)], [0, -np.sin(t23), np.cos(t23)] ])
    Uphases = np.array([ [np.exp(eta1*1j),0,0], [0, np.exp(eta2*1j),0], [0,0,1.0] ])
    U = np.dot(U23,np.dot(U13,np.dot(U12,Uphases)))
    return U

if __name__=='__main__':
      print( CasasIbarra(ranMnu=True)   )