casasibarra.py

def casasibarra(di,ranMnu=False):
    """
    di.keys()-> ['MH0','MA0','Mtr01','Mtr02','Mtr03',]
    """
    MH0=di['MH0'];MA0=di['MA0'];Mtr01=di['Mtr01'];Mtr02=di['Mtr02'];Mtr03=di['Mtr03']

    mnu1in,Dms2,Dma2,ThetSol,ThetAtm,ThetRec=neutrino_data()
    #Nupa r=       min,    best_fit, max

    #Thesol=np.array([0.278, 0.323, 0.375])
    #Theatm=np.array([0.278, 0.323, 0.375])

    #print Thetas()
    #break
    
    #Yout=casasibarra(dp)
    #"""requiered input parameters
    #   dp with keys: 'MH0',.....
    # writing the neutrino mass seesaw type
    Mtr01t=np.abs(Mtr0tilde(Mtr01,MH0,MA0))
    Mtr02t=np.abs(Mtr0tilde(Mtr02,MH0,MA0))
    Mtr03t=np.abs(Mtr0tilde(Mtr03,MH0,MA0))
    
    
    DMR=[ [np.sqrt( Mtr01t),0,0],[0,np.sqrt(Mtr02t),0],[0,0,np.sqrt(Mtr03t)] ]   
    
    #phases of the PMNS matrix
    #phases1=np.random.uniform(0,0*pi,3) # All the phases are set to be zero
    phases1=np.random.uniform(0,2*pi,3)                
    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 ranMnu 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 ranMnu else  np.sqrt(random.uniform(Dms2[0],Dms2[2]) + mnu1**2)                      ) 
    mnu3=1.*(np.sqrt(Dma2[1]+mnu1in**2)     if ranMnu else  np.sqrt(random.uniform(Dma2[0],Dma2[2]) + mnu1**2)                      ) 
    
    #light neutrino masses only for an estimation 
    #mnu1=0
    #mnu2=sqrt(8.2e-5*1e-18+mnu1**2)
    #mnu3=sqrt(2.74e-3*1e-18+mnu1**2)
    
    #Square root of left-handed nuetrino mass matrix 
    DMnu=[ [np.sqrt(mnu1),0,0],[0,np.sqrt(mnu2),0],[0,0,np.sqrt(mnu3)] ]
    
    #mixing angles using 3 sigma data from arxiv:1405.7540 (table I)                            #and asumming a Normal Hierarchy'''



    t12 = 1.*( np.arcsin(np.sqrt(ThetSol[1])) if ranMnu else np.arcsin(np.sqrt(random.uniform(ThetSol[0],ThetSol[2]))))
    t23 = 1.*( np.arcsin(np.sqrt(ThetAtm[1])) if ranMnu else np.arcsin(np.sqrt(random.uniform(ThetAtm[0],ThetAtm[2]))))
    t13 = 1.*( np.arcsin(np.sqrt(ThetRec[1])) if ranMnu else np.arcsin(np.sqrt(random.uniform(ThetRec[0],ThetRec[2]))))


    #
    c23 = np.cos(t23)
    s23 = np.sin(t23)
    c13 = np.cos(t13)
    s13 = np.sin(t13)
    c12 = np.cos(t12)
    s12 = np.sin(t12)
    tan12 = np.tan(t12)
    
    
    #Building PMNS matrix
    U12 = array([ [np.cos(t12),np.sin(t12),0], [-np.sin(t12),np.cos(t12),0], [0,0,1.0] ])
    U13 = array([ [np.cos(t13),0,np.sin(t13)*np.exp(-delta*1j)], [0,1.0,0], [-np.sin(t13)*exp(delta*1j),0,np.cos(t13)] ])
    U23 = array([ [1.0,0,0], [0,np.cos(t23),np.sin(t23)], [0,-np.sin(t23),np.cos(t23)] ])
    Uphases = array([ [np.exp(eta1*1j),0,0], [0,np.exp(eta2*1j),0], [0,0,1.0] ])
    U=dot(U23,dot(U13,dot(U12,Uphases)))

    #Building R matrix of the Casas-Ibarra parametrization
    
    phases2=np.random.uniform(0,2*pi,3) #real part of mixing angles of R
    

    b12 = 1.*(0 if ranMnu else phases2[0])
    b13 = 1.*(0 if ranMnu else phases2[1])
    b23 = 1.*(0 if ranMnu else phases2[2])
    
  
    # R 
    R12 = array([ [np.cos(b12),np.sin(b12),0], [-np.sin(b12),np.cos(b12),0], [0,0,1.0] ])
    R13 = array([ [np.cos(b13),0,np.sin(b13)], [0,1.0,0], [-np.sin(b13),0,np.cos(b13)] ])
    R23 = array([ [1.0,0,0], [0,np.cos(b23),np.sin(b23)], [0,-np.sin(b23),np.cos(b23)] ])
    
    R=dot(R23,dot(R13,R12))
    #print np.dot(R, R.T)   
    #Yukawa matrix of the Casas-Ibarra parametrization
    yuk2=dot(DMR,dot(R,dot(DMnu,transpose(conjugate(U)))))

    return yuk2
        



def leptophilic(dp,theta,phi,precision=0.9):
    """
       Normalized yukawa couplings are expressed as:
       ye_N  = np.sin(phi)*np.cos(theta) +- error
       ymu_N = np.sin(phi)*np.sin(theta) +- error
       ytau-N= np.cos(phi)
       Angles defined in radians
       examples:
       1. e-philic phi=pi/2 and theta=0
       2. mu-philic phi=pi/2 and theta=pi/2 
       normalized electron yukawa equal to precision
       DEBUG: implement minimization 
    """
    for mloop in range(0,10000000):
        yuk  = casasibarra(dp,ranMnu=False) 
        #Yukawa sector associated to the first generation of sigma
        yff=np.abs(yuk[0,:])
        #Renormalized Yukawa
        Yuk1R = yff[0]/np .sqrt(yff[0]**2 + yff[1]**2 +yff[2]**2)
        Yuk2R = yff[1]/np.sqrt(yff[0]**2 + yff[1]**2 +yff[2]**2)
        Yuk3R = yff[2]/np.sqrt(yff[0]**2 + yff[1]**2 +yff[2]**2)
        if ((theta < np.pi/(2*100)) or (theta > 0.99*np.pi/2) or (phi < np.pi/(2*100)) or (phi > 0.99*np.pi/2) ):
            error =5.0*precision
        else:
            error=precision
        if ((np.sin(phi)*np.cos(theta) - error < Yuk1R < np.sin(phi)*np.cos(theta) + error) and
             (np.sin(phi)*np.sin(theta) - error < Yuk2R < np.sin(phi)*np.sin(theta) + error)):
               a=yuk, Yuk1R, Yuk2R, Yuk3R,error
   
               return a # Modifications has been done here!