Index: trunk/omega/src/algebra.ml =================================================================== --- trunk/omega/src/algebra.ml (revision 8848) +++ trunk/omega/src/algebra.ml (revision 8849) @@ -1,803 +1,807 @@ (* algebra.ml -- Copyright (C) 1999-2022 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) (* Avoid refering to [Pervasives.compare], because [Pervasives] will become [Stdlib.Pervasives] in O'Caml 4.07 and [Stdlib] in O'Caml 4.08. *) let pcompare = compare module type Test = sig val suite : OUnit.test end (* The terms will be small and there's no need to be fancy and/or efficient. It's more important to have a unique representation. *) module PM = Pmap.List (* \thocwmodulesection{Coefficients} *) (* For our algebra, we need coefficient rings. *) module type CRing = sig type t val null : t val unit : t val mul : t -> t -> t val add : t -> t -> t val sub : t -> t -> t val neg : t -> t val to_string : t -> string end (* And rational numbers provide a particularly important example: *) module type Rational = sig include CRing val is_null : t -> bool val is_unit : t -> bool val is_positive : t -> bool val is_negative : t -> bool val is_integer : t -> bool val make : int -> int -> t val abs : t -> t val inv : t -> t val div : t -> t -> t val pow : t -> int -> t val sum : t list -> t val to_ratio : t -> int * int val to_float : t -> float val to_integer : t -> int module Test : Test end (* \thocwmodulesection{Naive Rational Arithmetic} *) (* \begin{dubious} This \emph{is} dangerous and will overflow even for simple applications. The production code will have to be linked to a library for large integer arithmetic. \end{dubious} *) (* Anyway, here's Euclid's algorithm: *) let rec gcd i1 i2 = if i2 = 0 then abs i1 else gcd i2 (i1 mod i2) let lcm i1 i2 = (i1 / gcd i1 i2) * i2 module Small_Rational : Rational = struct type t = int * int let is_null (n, _) = (n = 0) let is_unit (n, d) = (n <> 0) && (n = d) let is_positive (n, d) = n * d > 0 let is_negative (n, d) = n * d < 0 let is_integer (n, d) = (gcd n d = d) let null = (0, 1) let unit = (1, 1) let make n d = let c = gcd n d in (n / c, d / c) let abs (n, d) = (abs n, abs d) let inv (n, d) = (d, n) let mul (n1, d1) (n2, d2) = make (n1 * n2) (d1 * d2) let div q1 q2 = mul q1 (inv q2) let add (n1, d1) (n2, d2) = make (n1 * d2 + n2 * d1) (d1 * d2) let sub (n1, d1) (n2, d2) = make (n1 * d2 - n2 * d1) (d1 * d2) let neg (n, d) = (- n, d) let rec pow q p = if p = 0 then unit else if p < 0 then pow (inv q) (-p) else mul q (pow q (pred p)) let sum qs = List.fold_right add qs null let to_ratio (n, d) = if d < 0 then (-n, -d) else (n, d) let to_float (n, d) = float n /. float d let to_string (n, d) = if d = 1 then Printf.sprintf "%d" n else let n, d = to_ratio (n, d) in Printf.sprintf "(%d/%d)" n d let to_integer (n, d) = if is_integer (n, d) then n else invalid_arg "Algebra.Small_Rational.to_integer" module Test = struct open OUnit let equal z1 z2 = is_null (sub z1 z2) let assert_equal_rational z1 z2 = assert_equal ~printer:to_string ~cmp:equal z1 z2 let suite_mul = "mul" >::: [ "1*1=1" >:: (fun () -> assert_equal_rational (mul unit unit) unit) ] let suite = "Algebra.Small_Rational" >::: [suite_mul] end end module Q = Small_Rational (* \thocwmodulesection{Rational Complex Numbers} *) module type QComplex = sig type q type t val make : q -> q -> t val null : t val unit : t val real : t -> q val imag : t -> q val conj : t -> t val neg : t -> t val add : t -> t -> t val sub : t -> t -> t val mul : t -> t -> t val inv : t -> t val div : t -> t -> t val pow : t -> int -> t val sum : t list -> t val is_null : t -> bool val is_unit : t -> bool val is_positive : t -> bool val is_negative : t -> bool val is_integer : t -> bool val is_real : t -> bool val to_string : t -> string module Test : Test end module QComplex (Q : Rational) : QComplex with type q = Q.t = struct type q = Q.t type t = { re : q; im : q } let make re im = { re; im } let null = { re = Q.null; im = Q.null } let unit = { re = Q.unit; im = Q.null } let real z = z.re let imag z = z.im let conj z = { re = z.re; im = Q.neg z.im } let neg z = { re = Q.neg z.re; im = Q.neg z.im } let add z1 z2 = { re = Q.add z1.re z2.re; im = Q.add z1.im z2.im } let sub z1 z2 = { re = Q.sub z1.re z2.re; im = Q.sub z1.im z2.im } let sum qs = List.fold_right add qs null (* Save one multiplication with respect to the standard formula \begin{equation} (x+iy)(u+iv) = \lbrack xu-yv\rbrack + i\lbrack(x+u)(y+v)-xu-yv\rbrack\, \end{equation} at the expense of one addition and two subtractions. *) let mul z1 z2 = let re12 = Q.mul z1.re z2.re and im12 = Q.mul z1.im z2.im in { re = Q.sub re12 im12; im = Q.sub (Q.sub (Q.mul (Q.add z1.re z1.im) (Q.add z2.re z2.im)) re12) im12 } let inv z = let modulus = Q.add (Q.mul z.re z.re) (Q.mul z.im z.im) in { re = Q.div z.re modulus; im = Q.div (Q.neg z.im) modulus } let div n d = mul (inv d) n let rec pow q p = if p = 0 then unit else if p < 0 then pow (inv q) (-p) else mul q (pow q (pred p)) let is_real q = Q.is_null q.im let test_real test q = is_real q && test q.re let is_null = test_real Q.is_null let is_unit = test_real Q.is_unit let is_positive = test_real Q.is_positive let is_negative = test_real Q.is_negative let is_integer = test_real Q.is_integer let q_to_string q = (if Q.is_negative q then "-" else " ") ^ Q.to_string (Q.abs q) let to_string z = if Q.is_null z.im then q_to_string z.re else if Q.is_null z.re then if Q.is_unit z.im then " I" else if Q.is_unit (Q.neg z.im) then "-I" else q_to_string z.im ^ "*I" else Printf.sprintf "(%s%s*I)" (Q.to_string z.re) (q_to_string z.im) module Test = struct open OUnit let equal z1 z2 = is_null (sub z1 z2) let assert_equal_complex z1 z2 = assert_equal ~printer:to_string ~cmp:equal z1 z2 let suite_mul = "mul" >::: [ "1*1=1" >:: (fun () -> assert_equal_complex (mul unit unit) unit) ] let suite = "Algebra.QComplex" >::: [suite_mul] end end module QC = QComplex(Q) (* \thocwmodulesection{Laurent Polynomials} *) module type Laurent = sig type c type t val null : t val is_null : t -> bool val unit : t val atom : c -> int -> t val const : c -> t val scale : c -> t -> t + val neg : t -> t val add : t -> t -> t val diff : t -> t -> t val sum : t list -> t val mul : t -> t -> t val product : t list -> t val pow : int -> t -> t val eval : c -> t -> c val compare : t -> t -> int val to_string : string -> t -> string val pp : Format.formatter -> t -> unit module Test : Test end module Laurent : Laurent with type c = QC.t = struct module IMap = Map.Make (struct type t = int let compare i1 i2 = pcompare i2 i1 end) type c = QC.t let qc_minus_one = QC.neg QC.unit type t = c IMap.t let null = IMap.empty let is_null l = IMap.for_all (fun _ -> QC.is_null) l let atom qc n = if qc = QC.null then null else IMap.singleton n qc let const z = atom z 0 let unit = const QC.unit let add1 n qc l = try let qc' = QC.add qc (IMap.find n l) in if qc' = QC.null then IMap.remove n l else IMap.add n qc' l with | Not_found -> IMap.add n qc l let add l1 l2 = IMap.fold add1 l1 l2 let sum = function | [] -> null | [l] -> l | l :: l_list -> List.fold_left add l l_list let scale qc l = IMap.map (QC.mul qc) l + let neg l = + IMap.map QC.neg l + let diff l1 l2 = add l1 (scale qc_minus_one l2) (* cf.~[Product.fold2_rev] *) let fold2 f l1 l2 acc = IMap.fold (fun n1 qc1 acc1 -> IMap.fold (fun n2 qc2 acc2 -> f n1 qc1 n2 qc2 acc2) l2 acc1) l1 acc let mul l1 l2 = fold2 (fun n1 qc1 n2 qc2 acc -> add1 (n1 + n2) (QC.mul qc1 qc2) acc) l1 l2 null let product = function | [] -> unit | [l] -> l | l :: l_list -> List.fold_left mul l l_list let poly_pow multiply one inverse n x = let rec pow' i x' acc = if i < 1 then acc else pow' (pred i) x' (multiply x' acc) in if n < 0 then let x' = inverse x in pow' (pred (-n)) x' x' else if n = 0 then one else pow' (pred n) x x let qc_pow n z = poly_pow QC.mul QC.unit QC.inv n z let pow n l = poly_pow mul unit (fun _ -> invalid_arg "Algebra.Laurent.pow") n l let q_to_string q = (if Q.is_positive q then "+" else "-") ^ Q.to_string (Q.abs q) let qc_to_string z = let r = QC.real z and i = QC.imag z in if Q.is_null i then q_to_string r else if Q.is_null r then if Q.is_unit i then "+I" else if Q.is_unit (Q.neg i) then "-I" else q_to_string i ^ "*I" else Printf.sprintf "(%s%s*I)" (Q.to_string r) (q_to_string i) let to_string1 name (n, qc) = if n = 0 then qc_to_string qc else if n = 1 then if QC.is_unit qc then name else if qc = qc_minus_one then "-" ^ name else Printf.sprintf "%s*%s" (qc_to_string qc) name else if n = -1 then Printf.sprintf "%s/%s" (qc_to_string qc) name else if n > 1 then if QC.is_unit qc then Printf.sprintf "%s^%d" name n else if qc = qc_minus_one then Printf.sprintf "-%s^%d" name n else Printf.sprintf "%s*%s^%d" (qc_to_string qc) name n else Printf.sprintf "%s/%s^%d" (qc_to_string qc) name (-n) let to_string name l = match IMap.bindings l with | [] -> "0" | l -> String.concat "" (List.map (to_string1 name) l) let pp fmt l = Format.fprintf fmt "%s" (to_string "N" l) let eval v l = IMap.fold (fun n qc acc -> QC.add (QC.mul qc (qc_pow n v)) acc) l QC.null let compare l1 l2 = pcompare (List.sort pcompare (IMap.bindings l1)) (List.sort pcompare (IMap.bindings l2)) let compare l1 l2 = IMap.compare pcompare l1 l2 module Test = struct open OUnit let equal l1 l2 = compare l1 l2 = 0 let assert_equal_laurent l1 l2 = assert_equal ~printer:(to_string "N") ~cmp:equal l1 l2 let suite_mul = "mul" >::: [ "(1+N)(1-N)=1-N^2" >:: (fun () -> assert_equal_laurent (sum [unit; atom (QC.neg QC.unit) 2]) (product [sum [unit; atom QC.unit 1]; sum [unit; atom (QC.neg QC.unit) 1]])); "(1+N)(1-1/N)=N-1/N" >:: (fun () -> assert_equal_laurent (sum [atom QC.unit 1; atom (QC.neg QC.unit) (-1)]) (product [sum [unit; atom QC.unit 1]; sum [unit; atom (QC.neg QC.unit) (-1)]])); ] let suite = "Algebra.Laurent" >::: [suite_mul] end end (* \thocwmodulesection{Expressions: Terms, Rings and Linear Combinations} *) (* The tensor algebra will be spanned by an abelian monoid: *) module type Term = sig type 'a t val unit : unit -> 'a t val is_unit : 'a t -> bool val atom : 'a -> 'a t val power : int -> 'a t -> 'a t val mul : 'a t -> 'a t -> 'a t val map : ('a -> 'b) -> 'a t -> 'b t val to_string : ('a -> string) -> 'a t -> string val derive : ('a -> 'b option) -> 'a t -> ('b * int * 'a t) list val product : 'a t list -> 'a t val atoms : 'a t -> 'a list end module type Ring = sig module C : Rational type 'a t val null : unit -> 'a t val unit : unit -> 'a t val is_null : 'a t -> bool val is_unit : 'a t -> bool val atom : 'a -> 'a t val scale : C.t -> 'a t -> 'a t val add : 'a t -> 'a t -> 'a t val sub : 'a t -> 'a t -> 'a t val mul : 'a t -> 'a t -> 'a t val neg : 'a t -> 'a t val derive_inner : ('a -> 'a t) -> 'a t -> 'a t (* this? *) val derive_inner' : ('a -> 'a t option) -> 'a t -> 'a t (* or that? *) val derive_outer : ('a -> 'b option) -> 'a t -> ('b * 'a t) list val sum : 'a t list -> 'a t val product : 'a t list -> 'a t val atoms : 'a t -> 'a list val to_string : ('a -> string) -> 'a t -> string end module type Linear = sig module C : Ring type ('a, 'c) t val null : unit -> ('a, 'c) t val atom : 'a -> ('a, 'c) t val singleton : 'c C.t -> 'a -> ('a, 'c) t val scale : 'c C.t -> ('a, 'c) t -> ('a, 'c) t val add : ('a, 'c) t -> ('a, 'c) t -> ('a, 'c) t val sub : ('a, 'c) t -> ('a, 'c) t -> ('a, 'c) t val partial : ('c -> ('a, 'c) t) -> 'c C.t -> ('a, 'c) t val linear : (('a, 'c) t * 'c C.t) list -> ('a, 'c) t val map : ('a -> 'c C.t -> ('b, 'd) t) -> ('a, 'c) t -> ('b, 'd) t val sum : ('a, 'c) t list -> ('a, 'c) t val atoms : ('a, 'c) t -> 'a list * 'c list val to_string : ('a -> string) -> ('c -> string) -> ('a, 'c) t -> string end module Term : Term = struct module M = PM type 'a t = ('a, int) M.t let unit () = M.empty let is_unit = M.is_empty let atom f = M.singleton f 1 let power p x = M.map (( * ) p) x let insert1 binop f p term = let p' = binop (try M.find compare f term with Not_found -> 0) p in if p' = 0 then M.remove compare f term else M.add compare f p' term let mul1 f p term = insert1 (+) f p term let mul x y = M.fold mul1 x y let map f term = M.fold (fun t -> mul1 (f t)) term M.empty let to_string fmt term = String.concat "*" (M.fold (fun f p acc -> (if p = 0 then "1" else if p = 1 then fmt f else "[" ^ fmt f ^ "]^" ^ string_of_int p) :: acc) term []) let derive derive1 x = M.fold (fun f p dx -> if p <> 0 then match derive1 f with | Some df -> (df, p, mul1 f (pred p) (M.remove compare f x)) :: dx | None -> dx else dx) x [] let product factors = List.fold_left mul (unit ()) factors let atoms t = List.map fst (PM.elements t) end module Make_Ring (C : Rational) (T : Term) : Ring = struct module C = C let one = C.unit module M = PM type 'a t = ('a T.t, C.t) M.t let null () = M.empty let is_null = M.is_empty let power t p = M.singleton t p let unit () = power (T.unit ()) one let is_unit t = unit () = t (* \begin{dubious} The following should be correct too, but produces to many false positives instead! What's going on? \end{dubious} *) let broken__is_unit t = match M.elements t with | [(t, p)] -> T.is_unit t || C.is_null p | _ -> false let atom t = power (T.atom t) one let scale c x = M.map (C.mul c) x let insert1 binop t c sum = let c' = binop (try M.find compare t sum with Not_found -> C.null) c in if C.is_null c' then M.remove compare t sum else M.add compare t c' sum let add x y = M.fold (insert1 C.add) x y let sub x y = M.fold (insert1 C.sub) y x (* One might be tempted to use [Product.outer_self M.fold] instead, but this would require us to combine~[tx] and~[cx] to~[(tx, cx)]. *) let fold2 f x y = M.fold (fun tx cx -> M.fold (f tx cx) y) x let mul x y = fold2 (fun tx cx ty cy -> insert1 C.add (T.mul tx ty) (C.mul cx cy)) x y (null ()) let neg x = sub (null ()) x let neg x = scale (C.neg C.unit) x (* Multiply the [derivatives] by [c] and add the result to [dx]. *) let add_derivatives derivatives c dx = List.fold_left (fun acc (df, dt_c, dt_t) -> add (mul df (power dt_t (C.mul c (C.make dt_c 1)))) acc) dx derivatives let derive_inner derive1 x = M.fold (fun t -> add_derivatives (T.derive (fun f -> Some (derive1 f)) t)) x (null ()) let derive_inner' derive1 x = M.fold (fun t -> add_derivatives (T.derive derive1 t)) x (null ()) let collect_derivatives derivatives c dx = List.fold_left (fun acc (df, dt_c, dt_t) -> (df, power dt_t (C.mul c (C.make dt_c 1))) :: acc) dx derivatives let derive_outer derive1 x = M.fold (fun t -> collect_derivatives (T.derive derive1 t)) x [] let sum terms = List.fold_left add (null ()) terms let product factors = List.fold_left mul (unit ()) factors let atoms t = ThoList.uniq (List.sort compare (ThoList.flatmap (fun (t, _) -> T.atoms t) (PM.elements t))) let to_string fmt sum = "(" ^ String.concat " + " (M.fold (fun t c acc -> if C.is_null c then acc else if C.is_unit c then T.to_string fmt t :: acc else if C.is_unit (C.neg c) then ("(-" ^ T.to_string fmt t ^ ")") :: acc else (C.to_string c ^ "*[" ^ T.to_string fmt t ^ "]") :: acc) sum []) ^ ")" end module Make_Linear (C : Ring) : Linear with module C = C = struct module C = C module M = PM type ('a, 'c) t = ('a, 'c C.t) M.t let null () = M.empty let is_null = M.is_empty let atom a = M.singleton a (C.unit ()) let singleton c a = M.singleton a c let scale c x = M.map (C.mul c) x let insert1 binop t c sum = let c' = binop (try M.find compare t sum with Not_found -> C.null ()) c in if C.is_null c' then M.remove compare t sum else M.add compare t c' sum let add x y = M.fold (insert1 C.add) x y let sub x y = M.fold (insert1 C.sub) y x let map f t = M.fold (fun a c -> add (f a c)) t M.empty let sum terms = List.fold_left add (null ()) terms let linear terms = List.fold_left (fun acc (a, c) -> add (scale c a) acc) (null ()) terms let partial derive t = let d t' = let dt' = derive t' in if is_null dt' then None else Some dt' in linear (C.derive_outer d t) let atoms t = let a, c = List.split (PM.elements t) in (a, ThoList.uniq (List.sort compare (ThoList.flatmap C.atoms c))) let to_string fmt cfmt sum = "(" ^ String.concat " + " (M.fold (fun t c acc -> if C.is_null c then acc else if C.is_unit c then fmt t :: acc else if C.is_unit (C.neg c) then ("(-" ^ fmt t ^ ")") :: acc else (C.to_string cfmt c ^ "*" ^ fmt t) :: acc) sum []) ^ ")" end Index: trunk/omega/src/color.mli =================================================================== --- trunk/omega/src/color.mli (revision 8848) +++ trunk/omega/src/color.mli (revision 8849) @@ -1,362 +1,378 @@ (* color.mli -- Copyright (C) 1999-2022 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) module type Test = sig val suite : OUnit.test val suite_long : OUnit.test end (* \thocwmodulesection{Quantum Numbers} *) (* Color is not necessarily the~$\textrm{SU}(3)$ of QCD. Conceptually, it can be any \emph{unbroken} symmetry (\emph{broken} symmetries correspond to [Model.flavor]). In order to keep the group theory simple, we confine ourselves to the fundamental and adjoint representation of a single~$\textrm{SU}(N_C)$ for the moment. Therefore, particles are either color singlets or live in the defining representation of $\textrm{SU}(N_C)$: [SUN]$(|N_C|)$, its conjugate [SUN]$(-|N_C|)$ or in the adjoint representation of $\textrm{SU}(N_C)$: [AdjSUN]$(N_C)$. *) type t = Singlet | SUN of int | AdjSUN of int val conjugate : t -> t val compare : t -> t -> int (* \thocwmodulesection{Color Flows} *) (* This computes the color flow as used by WHIZARD: *) module type Flow = sig type color type t = color list * color list val rank : t -> int val of_list : int list -> color val ghost : unit -> color val to_lists : t -> int list list val in_to_lists : t -> int list list val out_to_lists : t -> int list list val ghost_flags : t -> bool list val in_ghost_flags : t -> bool list val out_ghost_flags : t -> bool list (* A factor is a list of powers \begin{equation} \sum_{i} \left( \frac{\ocwlowerid{num}_i}{\ocwlowerid{den}_i} \right)^{\ocwlowerid{power}_i} \end{equation} *) type power = { num : int; den : int; power : int } type factor = power list val factor : t -> t -> factor val zero : factor module Test : Test end module Flow : Flow (* \thocwmodulesection{Vertex Color Flows} *) (* \begin{dubious} The following is (still work-in-progress) infrastructure for translating UFO style color factors into color flows. \end{dubious} *) (* \begin{dubious} It might be beneficial, to use the color flow representation here. This will simplify the colorizer at the price of some complexity in [UFO] or here. \end{dubious} *) (* The datatypes [Arrow.free] and [Arrow.factor] will be used as building blocks for [Birdtracks.t] below. *) module type Arrow = sig (* For fundamental and adjoint representations, the endpoints of arrows are uniquely specified by a vertex (which will be represented by a number). For representations with more than one outgoing or incoming arrow, we need an additional index. This is abrcated in the [endpoint] type. *) type endpoint (* Endpoints can be the the tip or tail of an arrow or a ghost. Currently, we use the types for illustration only, but we might eventually try to make them abstract for additional safety.. *) type tip = endpoint type tail = endpoint type ghost = endpoint (* The position of the endpoint is encoded as an integer, which can be mapped, if necessary. *) val position : endpoint -> int val relocate : (int -> int) -> endpoint -> endpoint (* An [Arrow.t] is either a genuine arrow or a ghost \ldots *) type ('tail, 'tip, 'ghost) t = | Arrow of 'tail * 'tip | Ghost of 'ghost + | Epsilon of 'tip list + | Epsilon_bar of 'tail list (* {}\ldots and we distuish [free] arrows that must not contain summation indices from [factor]s that may. Indices are opaque. [('tail, 'tip, 'ghost) t] is polymorphic so that we can use richer ['tail], ['tip] and ['ghost] in [factor]. *) type free = (tail, tip, ghost) t type factor (* For debugging, logging, etc. *) val free_to_string : free -> string val factor_to_string : factor -> string (* Change the [endpoint]s in a [free] arrow. *) val map : (endpoint -> endpoint) -> free -> free (* Turn the [endpoint]s satisfying the predicate into a left or right hand side summation index. Left and right refer to the two factors in a product and we must only match arrows with [endpoint]s in both factors, not double lines on either side. Typically, the predicate will be set up to select only the summation indices that appear on both sides.*) val to_left_factor : (endpoint -> bool) -> free -> factor val to_right_factor : (endpoint -> bool) -> free -> factor (* The incomplete inverse [of_factor] raises an exception if there are remaining summation indices. [is_free] can be used to check first. *) val of_factor : factor -> free val is_free : factor -> bool (* Return all the endpoints of the arrow that have a [position] encoded as a negative integer. These are treated as summation indices in our applications. *) val negatives : free -> endpoint list (* We will need to test whether an arrow represents a ghost. *) val is_ghost : free -> bool - (* Merging two arrows can give a variety of results: *) + (* An arrow looping back to itself. *) + val is_tadpole : factor -> bool + + (* Merging two arrows can give a variety of results. Note that + we return the determinant resulting from merging an~$\epsilon$ + and an~$\bar\epsilon$ rather than the list of [Arrow]s + with permuted tips to avoid having to pass the relative signs. + These will be handled by [Birdtracks] below. *) type merge = | Match of factor (* a tip fits the other's tail: make one arrow out of two *) + | Determinant of factor list list * factor list list (* even and odd parts of $\bar\epsilon_{kj_1j_2\ldots}\epsilon_{ki_1i_2\ldots}$ *) | Ghost_Match (* two matching ghosts *) | Loop_Match (* both tips fit both tails: drop the arrows *) | Mismatch (* ghost meets arrow: error *) | No_Match (* nothing to be done *) val merge : factor -> factor -> merge (* Break up an arrow [tee a (i => j) -> [i => a; a => j]], i.\,e.~insert - a gluon. *) + a gluon. Returns an empty list for a ghost and raises an exception + for~$\epsilon$ and~$\bar\epsilon$. *) val tee : int -> free -> free list -(* [dir i j arrow] returns the direction of the arrow relative to [j => i] *) +(* [dir i j arrow] returns the direction of the arrow relative to [j => i]. + Returns 0 for a ghost and raises an exception for~$\epsilon$ + and~$\bar\epsilon$. *) val dir : int -> int -> free -> int (* It's intuitive to use infix operators to construct the lines. *) val single : endpoint -> endpoint -> free val double : endpoint -> endpoint -> free list val ghost : endpoint -> free module Infix : sig (* [single i j] or [i => j] creates a single line from [i] to [j] and [i ==> j] is a shorthard for [[i => j]]. *) val (=>) : int -> int -> free val (==>) : int -> int -> free list (* [double i j] or [i <=> j] creates a double line from [i] to [j] and back. *) val (<=>) : int -> int -> free list (* Single lines with subindices at the tip and/or tail *) val (>=>) : int * int -> int -> free val (=>>) : int -> int * int -> free val (>=>>) : int * int -> int * int -> free (* [?? i] creates a ghost at [i]. *) val (??) : int -> free (* NB: I wanted to use [~~] instead of [??], but ocamlweb can't handle operators starting with [~] in the index properly. *) end + val epsilon : int list -> free + val epsilon_bar : int list -> free + (* [chain [1;2;3]] is a shorthand for [[1 => 2; 2 => 3]] and [cycle [1;2;3]] for [[1 => 2; 2 => 3; 3 => 1]]. Other lists and edge cases are handled in the natural way. *) val chain : int list -> free list val cycle : int list -> free list module Test : Test (* Pretty printer for the toplevel. *) val pp_free : Format.formatter -> free -> unit val pp_factor : Format.formatter -> factor -> unit end module Arrow : Arrow (* Possible color flows for a single propagator, as currently supported by WHIZARD. *) module type Propagator = sig type cf_in = int type cf_out = int type t = W | I of cf_in | O of cf_out | IO of cf_in * cf_out | G val to_string : t -> string end module Propagator : Propagator (* Implement birdtracks operations as generally as possible. Below, the signature will be extended with group specific generators for $\mathrm{SU}(N_C)$ and $\mathrm{U}(N_C)$ and even $N_C=3$. *) module type Birdtracks = sig type t (* Strip out redundancies. *) val canonicalize : t -> t (* Debugging, logging, etc. *) val to_string : t -> string (* Test for trivial color flows that are just a number. *) val trivial : t -> bool (* Test for vanishing coefficients. *) val is_null : t -> bool (* Purely numeric factors, implemented as Laurent polynomials (cf.~[Algebra.Laurent] in~$N_C$ with complex rational coefficients. *) val const : Algebra.Laurent.t -> t val null : t (* $0$ *) val one : t (* $1$ *) val two : t (* $2$ *) val half : t (* $1/2$ *) val third : t (* $1/3$ *) val minus : t (* $-1$ *) val int : int -> t (* $n$ *) val fraction : int -> t (* $1/n$ *) val nc : t (* $N_C$ *) val over_nc : t (* $1/N_C$ *) val imag : t (* $\ii$ *) (* Shorthand: $\{(c_i,p_i)\}_i\to \sum_i c_i (N_C)^{p_i}$*) val ints : (int * int) list -> t val scale : Algebra.QC.t -> t -> t val sum : t list -> t val diff : t -> t -> t val times : t -> t -> t val multiply : t list -> t (* For convenience, here are infix versions of the above operations. *) module Infix : sig val ( +++ ) : t -> t -> t val ( --- ) : t -> t -> t val ( *** ) : t -> t -> t end (* We can compute the $f_{abc}$ and $d_{abc}$ invariant tensors from the generators of an arbitrary representation: \begin{subequations} \begin{align} f_{a_1a_2a_3} &= - \ii \tr\left(T_{a_1}\left\lbrack T_{a_2},T_{a_3}\right\rbrack_-\right) = - \ii \tr\left(T_{a_1}T_{a_2}T_{a_3}\right) + \ii \tr\left(T_{a_1}T_{a_3}T_{a_2}\right) \\ d_{a_1a_2a_3} &= \tr\left(T_{a_1}\left\lbrack T_{a_2},T_{a_3}\right\rbrack_+\right) = \tr\left(T_{a_1}T_{a_2}T_{a_3}\right) + \tr\left(T_{a_1}T_{a_3}T_{a_2}\right)\, \end{align} \end{subequations} assuming the normalization $ \tr(T_aT_b) = \delta_{ab}$. NB: this uses the summation indices $-1$, $-2$ and $-3$. Therefore it \emph{must not} appear unevaluated more than once in a product! *) val f_of_rep : (int -> int -> int -> t) -> int -> int -> int -> t val d_of_rep : (int -> int -> int -> t) -> int -> int -> int -> t (* Rename the indices of endpoints in a birdtrack. *) val relocate : (int -> int) -> t -> t (* [fuse nc vertex children] use the color flows in the [vertex] to combine the color flows in the incoming [children] and return the color flows for outgoing particle together with their weights. *) val fuse : int -> t -> Propagator.t list -> (Algebra.QC.t * Propagator.t) list module Test : Test (* Pretty printer for the toplevel. *) val pp : Format.formatter -> t -> unit end module Birdtracks : Birdtracks module type SU3 = sig include Birdtracks val delta3 : int -> int -> t val delta8 : int -> int -> t val delta8_loop : int -> int -> t val gluon : int -> int -> t val delta6 : int -> int -> t val delta10 : int -> int -> t val t : int -> int -> int -> t val f : int -> int -> int -> t val d : int -> int -> int -> t - val epsilon : int -> int -> int -> t - val epsilonbar : int -> int -> int -> t + val epsilon : int list -> t + val epsilon_bar : int list -> t val t8 : int -> int -> int -> t val t6 : int -> int -> int -> t val t10 : int -> int -> int -> t val k6 : int -> int -> int -> t val k6bar : int -> int -> int -> t val delta_of_tableau : int Young.tableau -> int -> int -> t val t_of_tableau : int Young.tableau -> int -> int -> int -> t end module SU3 : SU3 module Vertex : SU3 (* \begin{dubious} This must not be used, because it has not yet been updated to the correctly symmetrized version! \end{dubious} *) module U3 : SU3 Index: trunk/omega/src/algebra.mli =================================================================== --- trunk/omega/src/algebra.mli (revision 8848) +++ trunk/omega/src/algebra.mli (revision 8849) @@ -1,295 +1,296 @@ (* algebra.mli -- Copyright (C) 1999-2022 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) module type Test = sig val suite : OUnit.test end (* \thocwmodulesection{Coefficients} *) (* For our algebra, we need coefficient rings. *) module type CRing = sig type t val null : t val unit : t val mul : t -> t -> t val add : t -> t -> t val sub : t -> t -> t val neg : t -> t val to_string : t -> string end (* And rational numbers provide a particularly important example: *) module type Rational = sig include CRing val is_null : t -> bool val is_unit : t -> bool val is_positive : t -> bool val is_negative : t -> bool val is_integer : t -> bool val make : int -> int -> t val abs : t -> t val inv : t -> t val div : t -> t -> t val pow : t -> int -> t val sum : t list -> t val to_ratio : t -> int * int val to_float : t -> float val to_integer : t -> int module Test : Test end (* \thocwmodulesection{Naive Rational Arithmetic} *) (* \begin{dubious} This \emph{is} dangerous and will overflow even for simple applications. The production code will have to be linked to a library for large integer arithmetic. \end{dubious} *) module Small_Rational : Rational module Q : Rational (* \thocwmodulesection{Rational Complex Numbers} *) module type QComplex = sig type q type t val make : q -> q -> t val null : t val unit : t val real : t -> q val imag : t -> q val conj : t -> t val neg : t -> t val add : t -> t -> t val sub : t -> t -> t val mul : t -> t -> t val inv : t -> t val div : t -> t -> t val pow : t -> int -> t val sum : t list -> t val is_null : t -> bool val is_unit : t -> bool val is_positive : t -> bool val is_negative : t -> bool val is_integer : t -> bool val is_real : t -> bool val to_string : t -> string module Test : Test end module QComplex : functor (Q' : Rational) -> QComplex with type q = Q'.t module QC : QComplex with type q = Q.t (* \thocwmodulesection{Laurent Polynomials} *) (* Polynomials, including negative powers, in one variable. In our applications, the variable~$x$ will often be~$N_C$, the number of colors \begin{equation} \sum_n c_n N_C^n \end{equation} *) module type Laurent = sig (* The type of coefficients. In the implementation below, it is [QComplex.t]: complex numbers with rational real and imaginary parts. *) type c type t (* Elementary constructors *) val null : t val is_null : t -> bool val unit : t (* [atom c n] constructs a term $c x^n$, where $x$ denotes the variable. *) val atom : c -> int -> t (* Shortcut: [const c = atom c 0] *) val const : c -> t (* Elementary arithmetic *) val scale : c -> t -> t + val neg : t -> t val add : t -> t -> t val diff : t -> t -> t val sum : t list -> t val mul : t -> t -> t val product : t list -> t val pow : int -> t -> t (* [eval c p] evaluates the polynomial [p] by substituting the constant [c] for the variable. *) val eval : c -> t -> c (* A total ordering. Does not correspond to any mathematical order. *) val compare : t -> t -> int (* Logging, debugging and toplevel integration. *) val to_string : string -> t -> string val pp : Format.formatter -> t -> unit module Test : Test end (* \begin{dubious} Could (should?) be functorialized over [QComplex]. We had to wait until we upgraded our O'Caml requirements to 4.02, but that has been done. \end{dubious} *) module Laurent : Laurent with type c = QC.t (* \thocwmodulesection{Expressions: Terms, Rings and Linear Combinations} *) (* The tensor algebra will be spanned by an abelian monoid: *) module type Term = sig type 'a t val unit : unit -> 'a t val is_unit : 'a t -> bool val atom : 'a -> 'a t val power : int -> 'a t -> 'a t val mul : 'a t -> 'a t -> 'a t val map : ('a -> 'b) -> 'a t -> 'b t val to_string : ('a -> string) -> 'a t -> string (* The derivative of a term is \emph{not} a term, but a sum of terms instead: \begin{equation} D (f_1^{p_1}f_2^{p_2}\cdots f_n^{p_n}) = \sum_i (Df_i) p_i f_1^{p_1}f_2^{p_2}\cdots f_i^{p_i-1} \cdots f_n^{p_n} \end{equation} The function returns the sum as a list of triples $(Df_i,p_i, f_1^{p_1}f_2^{p_2}\cdots f_i^{p_i-1} \cdots f_n^{p_n})$. Summing the terms is left to the calling module and the $Df_i$ are \emph{not} guaranteed to be different. NB: The function implementating the inner derivative, is supposed to return~[Some]~$Df_i$ and [None], iff~$Df_i$ vanishes. *) val derive : ('a -> 'b option) -> 'a t -> ('b * int * 'a t) list (* convenience function *) val product : 'a t list -> 'a t val atoms : 'a t -> 'a list end module type Ring = sig module C : Rational type 'a t val null : unit -> 'a t val unit : unit -> 'a t val is_null : 'a t -> bool val is_unit : 'a t -> bool val atom : 'a -> 'a t val scale : C.t -> 'a t -> 'a t val add : 'a t -> 'a t -> 'a t val sub : 'a t -> 'a t -> 'a t val mul : 'a t -> 'a t -> 'a t val neg : 'a t -> 'a t (* Again \begin{equation} D (f_1^{p_1}f_2^{p_2}\cdots f_n^{p_n}) = \sum_i (Df_i) p_i f_1^{p_1}f_2^{p_2}\cdots f_i^{p_i-1} \cdots f_n^{p_n} \end{equation} but, iff~$Df_i$ can be identified with a~$f'$, we know how to perform the sum. *) val derive_inner : ('a -> 'a t) -> 'a t -> 'a t (* this? *) val derive_inner' : ('a -> 'a t option) -> 'a t -> 'a t (* or that? *) (* Below, we will need partial derivatives that lead out of the ring: [derive_outer derive_atom term] returns a list of partial derivatives ['b] with non-zero coefficients ['a t]: *) val derive_outer : ('a -> 'b option) -> 'a t -> ('b * 'a t) list (* convenience functions *) val sum : 'a t list -> 'a t val product : 'a t list -> 'a t (* The list of all generators appearing in an expression: *) val atoms : 'a t -> 'a list val to_string : ('a -> string) -> 'a t -> string end module type Linear = sig module C : Ring type ('a, 'c) t val null : unit -> ('a, 'c) t val atom : 'a -> ('a, 'c) t val singleton : 'c C.t -> 'a -> ('a, 'c) t val scale : 'c C.t -> ('a, 'c) t -> ('a, 'c) t val add : ('a, 'c) t -> ('a, 'c) t -> ('a, 'c) t val sub : ('a, 'c) t -> ('a, 'c) t -> ('a, 'c) t (* A partial derivative w.\,r.\,t.~a vector maps from a coefficient ring to the dual vector space. *) val partial : ('c -> ('a, 'c) t) -> 'c C.t -> ('a, 'c) t (* A linear combination of vectors \begin{equation} \text{[linear]} \lbrack (v_1, c_1); (v_2, c_2); \ldots; (v_n, c_n)\rbrack = \sum_{i=1}^{n} c_i\cdot v_i \end{equation} *) val linear : (('a, 'c) t * 'c C.t) list -> ('a, 'c) t (* Some convenience functions *) val map : ('a -> 'c C.t -> ('b, 'd) t) -> ('a, 'c) t -> ('b, 'd) t val sum : ('a, 'c) t list -> ('a, 'c) t (* The list of all generators and the list of all generators of coefficients appearing in an expression: *) val atoms : ('a, 'c) t -> 'a list * 'c list val to_string : ('a -> string) -> ('c -> string) -> ('a, 'c) t -> string end module Term : Term module Make_Ring (C : Rational) (T : Term) : Ring module Make_Linear (C : Ring) : Linear with module C = C Index: trunk/omega/src/omega.tex =================================================================== --- trunk/omega/src/omega.tex (revision 8848) +++ trunk/omega/src/omega.tex (revision 8849) @@ -1,1200 +1,1201 @@ % omega.tex -- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \NeedsTeXFormat{LaTeX2e} \RequirePackage{ifpdf} \ifpdf \documentclass[a4paper,notitlepage,chapters]{flex} \usepackage{type1cm} \usepackage[pdftex,colorlinks]{hyperref} \usepackage[pdftex]{graphicx,feynmp,emp} \DeclareGraphicsRule{*}{mps}{*}{} \else \documentclass[a4paper,notitlepage,chapters]{flex} \usepackage[T1]{fontenc} % \usepackage[hypertex]{hyperref} \usepackage{graphicx,feynmp,emp} \fi \usepackage{verbatim,array,amsmath,amssymb} \usepackage{url,thophys,thohacks} +\usepackage{pgf} \usepackage{ytableau} \setlength{\unitlength}{1mm} \empaddtoTeX{\usepackage{amsmath,amssymb}} \empaddtoTeX{\usepackage{thophys,thohacks}} \empaddtoprelude{input graph;} \empaddtoprelude{input boxes;} \IfFileExists{geometry.sty}% {\usepackage{geometry}% \geometry{a4paper,margin=2cm}}% {\relax} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% This should be part of flex.cls and/or thopp.sty \makeatletter \@ifundefined{frontmatter}% {\def\frontmatter{\pagenumbering{roman}}% \def\mainmatter{\cleardoublepage\pagenumbering{arabic}}} {} \makeatother %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% \makeatletter %%% %%% Italic figure captions to separate them visually from the text %%% %%% (this should be supported by flex.cls): %%% \makeatletter %%% \@secpenalty=-1000 %%% \def\fps@figure{t} %%% \def\fps@table{b} %%% \long\def\@makecaption#1#2{% %%% \vskip\abovecaptionskip %%% \sbox\@tempboxa{#1: \textit{#2}}% %%% \ifdim\wd\@tempboxa>\hsize %%% #1: \textit{#2}\par %%% \else %%% \global\@minipagefalse %%% \hb@xt@\hsize{\hfil\box\@tempboxa\hfil}% %%% \fi %%% \vskip\belowcaptionskip} %%% \makeatother \widowpenalty=4000 \clubpenalty=4000 \displaywidowpenalty=4000 %%% \pagestyle{headings} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \allowdisplaybreaks \renewcommand{\topfraction}{0.8} \renewcommand{\bottomfraction}{0.8} \renewcommand{\textfraction}{0.2} \setlength{\abovecaptionskip}{.5\baselineskip} \setlength{\belowcaptionskip}{\baselineskip} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% allow VERY overfull hboxes \setlength{\hfuzz}{5cm} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \usepackage{noweb} %%% \usepackage{nocondmac} \setlength{\nwmarginglue}{1em} \noweboptions{smallcode,noidentxref}%%%{webnumbering} %%% Saving paper: \def\nwendcode{\endtrivlist\endgroup} \nwcodepenalty=0 \let\nwdocspar\relax %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newcommand{\ttfilename}[1]{\texttt{\detokenize{#1}}} \usepackage[noweb,bypages]{ocamlweb} \empaddtoTeX{\usepackage[noweb,bypages]{ocamlweb}} \renewcommand{\ocwinterface}[1]{\section{Interface of \ocwupperid{#1}}} \renewcommand{\ocwmodule}[1]{\section{Implementation of \ocwupperid{#1}}} \renewcommand{\ocwinterfacepart}{\relax} \renewcommand{\ocwcodepart}{\relax} \renewcommand{\ocwbeginindex}{\begin{theindex}} \newcommand{\thocwmodulesection}[1]{\subsection{#1}} \newcommand{\thocwmodulesubsection}[1]{\subsubsection{#1}} \newcommand{\thocwmoduleparagraph}[1]{\paragraph{#1}} \renewcommand{\ocwindent}[1]{\noindent\ignorespaces} \renewcommand{\ocwbegincode}{\renewcommand{\ocwindent}[1]{\noindent\kern##1}} \renewcommand{\ocwendcode}{\renewcommand{\ocwindent}[1]{\noindent\ignorespaces}} \renewcommand{\ocweol}{\setlength\parskip{0pt}\par} \makeatletter \renewcommand{\@oddfoot}{\reset@font\hfil\thepage\hfil} \let\@evenfoot\@oddfoot \def\@evenhead{\leftmark{} \hrulefill}% \def\@oddhead{\hrulefill{} \rightmark}% \let\@mkboth\markboth \renewcommand{\chaptermark}[1]{\markboth{\hfil}{\hfil}}% \renewcommand{\sectionmark}[1]{\markboth{#1}{#1}} \renewcommand{\chapter}{% \clearpage\global\@topnum\z@\@afterindentfalse \secdef\@chapter\@schapter} \makeatother \newcommand{\signature}[1]{% \InputIfFileExists{#1.interface}{}% {\begin{dubious}\textit{Interface \ttfilename{#1.mli} unavailable!}\end{dubious}}} \newcommand{\application}[1]{% \InputIfFileExists{#1.implementation}{}% {\begin{dubious}\textit{Application \ttfilename{#1.ml} unavailable!}\end{dubious}}} \newcommand{\module}[1]{% \label{mod:#1}% \InputIfFileExists{#1.interface}{}% {\begin{dubious}\textit{Interface \ttfilename{#1.mli} unavailable!}\end{dubious}}% \InputIfFileExists{#1.implementation}{}% {\begin{dubious}\textit{Implementation \ttfilename{#1.ml} unavailable!}\end{dubious}}} \newcommand{\lexer}[1]{\application{#1_lexer}} \renewcommand{\ocwlexmodule}[1]{\relax} \newcommand{\parser}[1]{\application{#1_parser}} \renewcommand{\ocwyaccmodule}[1]{\relax} \newcommand{\thocwincludegraphics}[2]{\includegraphics[#1]{#2}} \ifpdf \newcommand{\thocwdefref}[1]{\textbf{\pageref{#1}}}% \newcommand{\thocwuseref}[1]{\textrm{\pageref{#1}}}% \renewcommand{\ocwrefindexentry}[5]% {\item #1,\quad\let\ref\thocwdefref{#4}, used: \let\ref\thocwuseref{#5}} \fi \newcommand{\thocwmakebox}[4]{\makebox(#1,#2)[#3]{#4}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newenvironment{modules}[1]% {\begin{list}{}% {\setlength{\leftmargin}{3em}% \setlength{\rightmargin}{2em}% \setlength{\itemindent}{-1em}% \setlength{\listparindent}{0pt}% %%%\setlength{\itemsep}{0pt}% \settowidth{\labelwidth}{\textbf{\ocwupperid{#1}:}}% \renewcommand{\makelabel}[1]{\ocwupperid{##1:}}}}% {\end{list}} \newenvironment{JR}% {\begin{dubious}\textit{JR sez' (regarding the Majorana Feynman rules):}} {\textit{(JR's probably right, but I need to check myself \ldots)} \end{dubious}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \DeclareMathOperator{\tr}{tr} \newcommand{\dd}{\mathrm{d}} \newcommand{\ii}{\mathrm{i}} \newcommand{\ee}{\mathrm{e}} \renewcommand{\Re}{\text{Re}} \renewcommand{\Im}{\text{Im}} \newcommand{\ketbra}[2]{\ket{#1}\!\bra{#2}} \newcommand{\Ketbra}[2]{\Ket{#1}\!\Bra{#2}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \makeindex \begin{document} \begin{fmffile}{\jobname pics} \fmfset{arrow_ang}{10} \fmfset{curly_len}{2mm} \fmfset{wiggly_len}{3mm} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \fmfcmd{% numeric joindiameter; joindiameter := 7thick;} \fmfcmd{% vardef sideways_at (expr d, p, frac) = save len; len = length p; (point frac*len of p) shifted ((d,0) rotated (90 + angle direction frac*len of p)) enddef; secondarydef p sideways d = for frac = 0 step 0.01 until 0.99: sideways_at (d, p, frac) .. endfor sideways_at (d, p, 1) enddef; secondarydef p choptail d = subpath (ypart (fullcircle scaled d shifted (point 0 of p) intersectiontimes p), infinity) of p enddef; secondarydef p choptip d = reverse ((reverse p) choptail d) enddef; secondarydef p pointtail d = fullcircle scaled d shifted (point 0 of p) intersectionpoint p enddef; secondarydef p pointtip d = (reverse p) pointtail d enddef; secondarydef pa join pb = pa choptip joindiameter .. pb choptail joindiameter enddef; vardef cyclejoin (expr p) = subpath (0.5*length p, infinity) of p join subpath (0, 0.5*length p) of p .. cycle enddef;} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \fmfcmd{% style_def double_line_arrow expr p = save pi, po; path pi, po; pi = reverse (p sideways thick); po = p sideways -thick; cdraw pi; cdraw po; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_line_arrow_beg expr p = save pi, po, pc; path pi, po, pc; pc = p choptail 7thick; pi = reverse (pc sideways thick); po = pc sideways -thick; cdraw pi .. p pointtail 5thick .. po; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_line_arrow_end expr p = save pi, po, pc; path pi, po, pc; pc = p choptip 7thick; pi = reverse (pc sideways thick); po = pc sideways -thick; cdraw po .. p pointtip 5thick .. pi; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_line_arrow_both expr p = save pi, po, pc; path pi, po, pc; pc = p choptip 7thick choptail 7thick; pi = reverse (pc sideways thick); po = pc sideways -thick; cdraw po .. p pointtip 5thick .. pi .. p pointtail 5thick .. cycle; cfill (arrow pi); cfill (arrow po); enddef;} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \fmfcmd{vardef middir (expr p, ang) = dir (angle direction length(p)/2 of p + ang) enddef;} \fmfcmd{style_def arrow_left expr p = shrink (.7); cfill (arrow p shifted (4thick * middir (p, 90))); endshrink enddef;} \fmfcmd{style_def arrow_right expr p = shrink (.7); cfill (arrow p shifted (4thick * middir (p, -90))); endshrink enddef;} \fmfcmd{style_def warrow_left expr p = shrink (.7); cfill (arrow p shifted (8thick * middir (p, 90))); endshrink enddef;} \fmfcmd{style_def warrow_right expr p = shrink (.7); cfill (arrow p shifted (8thick * middir (p, -90))); endshrink enddef;} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newcommand{\threeexternal}[3]{% \fmfsurround{d1,e1,d2,e2,d3,e3}% \fmfv{label=$#1$,label.ang=0}{e1}% \fmfv{label=$#2$,label.ang=180}{e2}% \fmfv{label=$#3$,label.ang=0}{e3}} \newcommand{\Threeexternal}[3]{% \fmfsurround{d1,e1,d3,e3,d2,e2}% \fmfv{label=$#1$,label.ang=0}{e1}% \fmfv{label=$#2$,label.ang=0}{e2}% \fmfv{label=$#3$,label.ang=180}{e3}} \newcommand{\Fourexternal}[4]{% \fmfsurround{d2,e2,d1,e1,d4,e4,d3,e3}% \fmfv{label=$#1$,label.ang=180}{e1}% \fmfv{label=$#2$,label.ang=0}{e2}% \fmfv{label=$#3$,label.ang=0}{e3}% \fmfv{label=$#4$,label.ang=180}{e4}} \newcommand{\Fiveexternal}[5]{% \fmfsurround{d2,e2,d1,e1,d5,e5,d4,e4,d3,e3}% \fmfv{label=$#1$,label.ang=180}{e1}% \fmfv{label=$#2$,label.ang=0}{e2}% \fmfv{label=$#3$,label.ang=0}{e3}% \fmfv{label=$#4$,label.ang=0}{e4}% \fmfv{label=$#5$,label.ang=180}{e5}} \newcommand{\twoincoming}{% \fmfdot{v}% \fmffreeze% \fmf{warrow_right}{e1,v}% \fmf{warrow_right}{e2,v}% \fmf{warrow_right}{v,e3}} \newcommand{\threeincoming}{% \fmfdot{v}% \fmffreeze% \fmf{warrow_right}{e1,v}% \fmf{warrow_right}{e2,v}% \fmf{warrow_right}{e3,v}} \newcommand{\threeoutgoing}{% \fmfdot{v}% \fmffreeze% \fmf{warrow_right}{v,e1}% \fmf{warrow_right}{v,e2}% \fmf{warrow_right}{v,e3}} \newcommand{\fouroutgoing}{% \threeoutgoing% \fmf{warrow_right}{v,e4}} \newcommand{\fiveoutgoing}{% \fouroutgoing% \fmf{warrow_right}{v,e5}} \newcommand{\setupthreegluons}{% \fmftop{g3} \fmfbottom{g1,g2} \fmf{phantom}{v,g1} \fmf{phantom}{v,g2} \fmf{phantom}{v,g3} \fmffreeze \fmfipair{v,g[],a[],b[]} \fmfiset{g1}{vloc (__g1)} \fmfiset{g2}{vloc (__g2)} \fmfiset{g3}{vloc (__g3)} \fmfiset{v}{vloc (__v)} \fmfiset{a1}{g1 shifted (-3thin,0)} \fmfiset{b1}{g1 shifted (+1thin,-2thin)} \fmfiset{a2}{g2 shifted (0,-3thin)} \fmfiset{b2}{g2 shifted (0,+3thin)} \fmfiset{a3}{g3 shifted (+1thin,+2thin)} \fmfiset{b3}{g3 shifted (-3thin,0)}} \begin{empfile} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \frontmatter \title{ O'Mega:\\ Optimal~Monte-Carlo\\ Event~Generation~Amplitudes} \author{% Thorsten Ohl\thanks{% \texttt{ohl@physik.uni-wuerzburg.de}, \texttt{http://physik.uni-wuerzburg.de/ohl}}\\ \hfil\\ Institut f\"ur Theoretische~Physik und Astrophysik\\ Julius-Maximilians-Universit\"at~W\"urzburg\\ Emil-Hilb-Weg 22, 97074~W\"urzburg, Germany\\ \hfil\\ J\"urgen Reuter\thanks{\texttt{juergen.reuter@desy.de}}\\ \hfil\\ DESY Theory Group, Notkestr. 85, 22603 Hamburg, Germany\\ \hfil\\ Wolfgang Kilian${}^{c,}$\thanks{\texttt{kilian@physik.uni-siegen.de}}\\ \hfil\\ Theoretische Physik 1\\ Universit\"at Siegen\\ Walter-Flex-Str.~3, 57068 Siegen, Germany\\ \hfil\\ with contributions from Christian Speckner${}^{d,}$\thanks{\texttt{cnspeckn@googlemail.com}}\\ as well as Christian Schwinn et al.} \date{\textbf{unpublished draft, printed \timestamp}} \maketitle \begin{abstract} \ldots \end{abstract} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage \begin{quote} Copyright \textcopyright~1999-2017 by \begin{itemize} \item Wolfgang~Kilian ~\texttt{} \item Thorsten~Ohl~\texttt{} \item J\"urgen~Reuter~\texttt{} \end{itemize} \end{quote} \begin{quote} WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. \end{quote} \begin{quote} WHIZARD is distributed in the hope that it will be useful, but \emph{without any warranty}; without even the implied warranty of \emph{merchantability} or \emph{fitness for a particular purpose}. See the GNU General Public License for more details. \end{quote} \begin{quote} You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. \end{quote} \setcounter{tocdepth}{2} \tableofcontents \mainmatter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Introduction} \label{sec:intro} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Complexity} \label{sec:complexity} \begin{table} \begin{center} \begin{tabular}{r|r|r} $n$ & $P(n)$& $F(n)$ \\\hline 4 & 3 & 3 \\ 5 & 10 & 15 \\ 6 & 25 & 105 \\ 7 & 56 & 945 \\ 8 & 119 & 10395 \\ 9 & 246 & 135135 \\ 10 & 501 & 2027025 \\ 11 & 1012 & 34459425 \\ 12 & 2035 & 654729075 \\ 13 & 4082 & 13749310575 \\ 14 & 8177 & 316234143225 \\ 15 & 16368 & 7905853580625 \\ 16 & 32751 & 213458046676875 \end{tabular} \end{center} \caption{\label{tab:P(n),F(n)} The number of $\phi^3$ Feynman diagrams~$F(n)$ and independent poles~$P(n)$.} \end{table} There are \begin{equation} P(n) = \frac{2^n-2}{2} - n = 2^{n-1} - n - 1 \end{equation} independent internal momenta in a $n$-particle scattering amplitude~\cite{ALPHA:1997}. This grows much slower than the number \begin{equation} F(n) = (2n-5)!! = (2n-5)\cdot(2n-7)\cdot\ldots\cdot3\cdot1 \end{equation} of tree Feynman diagrams in vanilla $\phi^3$ (see table~\ref{tab:P(n),F(n)}). There are no known corresponding expressions for theories with more than one particle type. However, empirical evidence from numerical studies~\cite{ALPHA:1997,HELAC:2000} as well as explicit counting results from O'Mega suggest \begin{equation} P^*(n) \propto 10^{n/2} \end{equation} while he factorial growth of the number of Feynman diagrams remains unchecked, of course. The number of independent momenta in an amplitude is a better measure for the complexity of the amplitude than the number of Feynman diagrams, since there can be substantial cancellations among the latter. Therefore it should be possible to express the scattering amplitude more compactly than by a sum over Feynman diagrams. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Ancestors} \label{sec:ancestors} Some of the ideas that O'Mega is based on can be traced back to HELAS~\cite{HELAS}. HELAS builts Feynman amplitudes by recursively forming off-shell `wave functions' from joining external lines with other external lines or off-shell `wave functions'. The program Madgraph~\cite{MADGRAPH:1994} automatically generates Feynman diagrams and writes a Fortran program corresponding to their sum. The amplitudes are calculated by calls to HELAS~\cite{HELAS}. Madgraph uses one straightforward optimization: no statement is written more than once. Since each statement corresponds to a collection of trees, this optimization is very effective for up to four particles in the final state. However, since the amplitudes are given as a sum of Feynman diagrams, this optimization can, by design, \emph{not} remove the factorial growth and is substantially weaker than the algorithms of~\cite{ALPHA:1997,HELAC:2000} and the algorithm of O'Mega for more particles in the final state. Then ALPHA~\cite{ALPHA:1997} (see also the slightly modified variant~\cite{HELAC:2000}) provided a numerical algorithm for calculating scattering amplitudes and it could be shown empirically, that the calculational costs are rising with a power instead of factorially. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Architecture} \label{sec:architecture} \begin{figure} \begin{center} \includegraphics[width=\textwidth]{modules} %includegraphics[height=.8\textheight]{modules} \end{center} \caption{\label{fig:modules}% Module dependencies in O'Mega.} %% The diamond shaped nodes are abstract signatures defininng functor %% domains and co-domains. The rectangular boxes are modules and %% functors and oval boxes are examples for applications. \end{figure} \subsection{General purpose libraries} Functions that are not specific to O'Mega and could be part of the O'Caml standard library \begin{modules}{} \item[ThoList] (mostly) simple convenience functions for lists that are missing from the standard library module \ocwupperid{List} (section~\ref{sec:tholist}, p.~\pageref{sec:tholist}) \item[Product] effcient tensor products for lists and sets (section~\ref{sec:product}, p.~\pageref{sec:product}) \item[Combinatorics] combinatorical formulae, sets of subsets, etc. (section~\ref{sec:combinatorics}, p.~\pageref{sec:combinatorics}) \end{modules} \subsection{O'Mega} The non-trivial algorithms that constitute O'Mega: \begin{modules}{} \item[DAG] Directed Acyclical Graphs (section~\ref{sec:DAG}, p.~\pageref{sec:DAG}) \item[Topology] unusual enumerations of unflavored tree diagrams (section~\ref{sec:topology}, p.~\pageref{sec:topology}) \item[Momentum] finite sums of external momenta (section~\ref{sec:momentum}, p.~\pageref{sec:momentum}) \item[Fusion] off shell wave functions (section~\ref{sec:fusion}, p.~\pageref{sec:fusion}) \item[Omega] functor constructing an application from a model and a target (section~\ref{sec:omega}, p.~\pageref{sec:omega}) \end{modules} \subsection{Abstract interfaces} The domains and co-domains of functors (section~\ref{sec:coupling}, p.~\pageref{sec:coupling}) \begin{modules}{} \item[Coupling] all possible couplings (not comprensive yet) \item[Model] physical models \item[Target] target programming languages \end{modules} \subsection{Models} (section~\ref{sec:models}, p.~\pageref{sec:models}) \begin{modules}{} \item[Modellib_SM.QED] Quantum Electrodynamics \item[Modellib_SM.QCD] Quantum Chromodynamics (not complete yet) \item[Modellib_SM.SM] Minimal Standard Model (not complete yet) \end{modules} etc. \subsection{Targets} Any programming language that supports arithmetic and a textual representation of programs can be targeted by O'Caml. The implementations translate the abstract expressions derived by \ocwupperid{Fusion} to expressions in the target (section~\ref{sec:targets}, p.~\pageref{sec:targets}). \begin{modules}{} \item[Targets.Fortran] Fortran95 language implementation, calling subroutines \end{modules} Other targets could come in the future: \texttt{C}, \texttt{C++}, O'Caml itself, symbolic manipulation languages, etc. \subsection{Applications} (section~\ref{sec:omega}, p.~\pageref{sec:omega}) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{The Big To Do Lists} \label{sec:TODO} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Required} All features required for leading order physics applications are in place. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Useful} \begin{enumerate} \item select allowed helicity combinations for massless fermions \item Weyl-Van der Waerden spinors \item speed up helicity sums by using discrete symmetries \item general triple and quartic vector couplings \item diagnostics: count corresponding Feynman diagrams more efficiently for more than ten external lines \item recognize potential cascade decays ($\tau$, $b$, etc.) \begin{itemize} \item warn the user to add additional \item kill fusions (at runtime), that contribute to a cascade \end{itemize} \item complete standard model in $R_\xi$-gauge \item groves (the simple method of cloned generations works) \end{enumerate} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Future Features} \begin{enumerate} \item investigate if unpolarized squared matrix elements can be calculated faster as traces of densitiy matrices. Unfortunately, the answer apears to be \emph{no} for fermions and \emph{up to a constant factor} for massive vectors. Since the number of fusions in the amplitude grows like~$10^{n/2}$, the number of fusions in the squared matrix element grows like~$10^n$. On the other hand, there are $2^{\#\text{fermions}+\#\text{massless vectors}} \cdot3^{\#\text{massive vectors}}$ terms in the helicity sum, which grows \emph{slower} than~$10^{n/2}$. The constant factor is probably also not favorable. However, there will certainly be asymptotic gains for sums over gauge (and other) multiplets, like color sums. \item compile Feynman rules from Lagrangians \item evaluate amplitues in O'Caml by compiling it to three address code for a virtual machine \begin{flushleft} \ocwkw{type}~$\ocwlowerid{mem}~=~\ocwlowerid{scalar}~$\ocwbt{array}~$% \times{}~\ocwlowerid{spinor}~$\ocwbt{array}~$% \times{}~\ocwlowerid{spinor}~$\ocwbt{array}~$% \times{}~\ocwlowerid{vector}~$\ocwbt{array}\\ \ocwkw{type}~$\ocwlowerid{instr}~=$\\ \qquad|~$\ocwupperid{VSS}~$\ocwkw{of}~\ocwbt{int}~$% \times{}~$\ocwbt{int}~$\times{}~$\ocwbt{int}\\ \qquad|~$\ocwupperid{SVS}~$\ocwkw{of}~\ocwbt{int}~$% \times{}~$\ocwbt{int}~$\times{}~$\ocwbt{int}\\ \qquad|~$\ocwupperid{AVA}~$\ocwkw{of}~\ocwbt{int}~$% \times{}~$\ocwbt{int}~$\times{}~$\ocwbt{int}\\ \qquad\ldots \end{flushleft} this could be as fast as~\cite{ALPHA:1997} or~\cite{HELAC:2000}. \item a virtual machine will be useful for for other target as well, because native code appears to become to large for most compilers for more than ten external particles. Bytecode might even be faster due to improved cache locality. \item use the virtual machine in O'Giga \end{enumerate} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Science Fiction} \begin{enumerate} \item numerical and symbolical loop calculations with \textsc{O'Tera: O'Mega Tool for Evaluating Renormalized Amplitudes} \end{enumerate} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Tuples and Polytuples} \label{sec:tuple} \module{tuple} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Topologies} \label{sec:topology} \module{topology} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Directed Acyclical Graphs} \label{sec:DAG} \module{DAG} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Momenta} \label{sec:momentum} \module{momentum} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Cascades} \label{sec:cascades} \module{cascade_syntax} \section{Lexer} \lexer{cascade} \section{Parser} \parser{cascade} \module{cascade} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Color} \label{sec:color} \module{color} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Fusions} \label{sec:fusion} \module{fusion} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Lorentz Representations, Couplings, Models and Targets} \label{sec:coupling} \signature{coupling} \signature{model} \module{dirac} \module{vertex} \signature{target} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Conserved Quantum Numbers} \label{sec:charges} \module{charges} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Colorization} \label{sec:colorize} \module{colorize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Processes} \label{sec:process} \module{process} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Model Files} \label{sec:model-files} \module{vertex_syntax} \section{Lexer} \lexer{vertex} \section{Parser} \parser{vertex} \module{vertex} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{UFO Models} \label{sec:ufo} \module{UFOx_syntax} \section{Expression Lexer} \lexer{UFOx} \section{Expression Parser} \parser{UFOx} \module{UFOx} \module{UFO_syntax} \section{Lexer} \lexer{UFO} \section{Parser} \parser{UFO} \module{UFO_Lorentz} \module{UFO} \section{Targets} \module{UFO_targets} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Hardcoded Targets} \label{sec:targets} \module{format_Fortran} \module{targets} \module{targets_Kmatrix} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Phase Space} \label{sec:phasespace} \module{phasespace} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Whizard} \label{sec:whizard} Talk to~\cite{Kilian:WHIZARD}. \module{whizard} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Applications} \label{sec:omega} \section{Sample} {\small\verbatiminput{sample.prc}} \module{omega} %application{omega_Phi3} %application{omega_Phi3h} %application{omega_Phi4} %application{omega_Phi4h} \application{omega_QED} %application{omega_QCD} %application{omega_SM3} %application{omega_SM3_ac} \application{omega_SM} \application{omega_SYM} %application{omega_SM_ac} %application{f90Maj_SM} %application{f90Maj_SM4} %application{omega_MSSM} %application{omega_MSSM_g} %application{omega_SM_Rxi} %application{omega_SM_clones} %application{omega_THDM} %application{omega_SMh} %application{omega_SM4h} %application{helas_QED} %application{helas_QCD} %application{helas_SM} %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% \chapter{O'Giga: O'Mega Graphical Interface for Generation and Analysis} %%% \label{sec:ogiga} %%% {\itshape NB: The code in this chapter \emph{must} be compiled with %%% \verb+-labels+, since \verb+lablgtk+ doesn't appear to work in classic mode.} %%% \begin{dubious} %%% Keep in mind that \texttt{ocamlweb} doesn't work properly with %%% O'Caml~3 yet. The colons in label declarations are typeset with %%% erroneous white space. %%% \end{dubious} %%% %%% \application{ogiga} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter*{Acknowledgements} We thank Mauro Moretti for fruitful discussions of the ALPHA algorithm~\cite{ALPHA:1997}, that inspired our solution of the double counting problem. We thank Wolfgang Kilian for providing the WHIZARD environment that turns our numbers into real events with unit weight. Thanks to the ECFA/DESY workshops and their participants for providing a showcase. Thanks to Edward Boos for discussions in Kaluza-Klein gravitons. This research is supported by Bundesministerium f\"ur Bildung und Forschung, Germany, (05\,HT9RDA) and Deutsche Forschungsgemeinschaft (MA\,676/6-1). Thanks to the Caml and Objective Caml teams from INRIA for the development and the lean and mean implementation of a programming language that does not insult the programmer's intelligence. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{thebibliography}{10} \bibitem{ALPHA:1997} F. Caravaglios, M. Moretti, Z.{} Phys.{} \textbf{C74} (1997) 291. \bibitem{HELAC:2000} A. Kanaki, C. Papadopoulos, DEMO-HEP-2000/01, hep-ph/0002082, February 2000. \bibitem{Ler97} Xavier Leroy, \textit{The Objective Caml system, documentation and user's guide}, Technical Report, INRIA, 1997. \bibitem{Okasaki:1998:book} Chris Okasaki, \textit{Purely Functional Data Structures}, Cambridge University Press, 1998. \bibitem{HELAS} H. Murayama, I. Watanabe, K. Hagiwara, KEK Report 91-11, January 1992. \bibitem{MADGRAPH:1994} T. Stelzer, W.F. Long, Comput.{} Phys.{} Commun.{} \textbf{81} (1994) 357. \bibitem{Denner:Majorana} A. Denner, H. Eck, O. Hahn and J. K\"ublbeck, Phys.{} Lett.{} \textbf{B291} (1992) 278; Nucl.{} Phys.{} \textbf{B387} (1992) 467. \bibitem{Barger/etal:1992:color} V.~Barger, A.~L.~Stange, R.~J.~N.~Phillips, Phys.~Rev.~\textbf{D45}, (1992) 1751. \bibitem{Ohl:LOTR} T. Ohl, \textit{Lord of the Rings}, (Computer algebra library for O'Caml, unpublished). \bibitem{Ohl:bocages} T. Ohl, \textit{Bocages}, (Feynman diagram library for O'Caml, unpublished). \bibitem{Kilian:WHIZARD} W. Kilian, \textit{\texttt{WHIZARD}}, University of Karlsruhe, 2000. \bibitem{Boos/Ohl:groves} E.\,E. Boos, T. Ohl, Phys.\ Rev.\ Lett.\ \textbf{83} (1999) 480. \bibitem{Han/Lykken/Zhang:1999:Kaluza-Klein} T.~Han, J.~D.~Lykken and R.~Zhang, %``On Kaluza-Klein states from large extra dimensions,'' Phys.{} Rev.{} \textbf{D59} (1999) 105006 [hep-ph/9811350]. %%CITATION = HEP-PH 9811350;%% \bibitem{PTVF92} William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, \textit{Numerical Recipes: The Art of Scientific Computing}, Second Edition, Cambridge University Press, 1992. \bibitem{Cvi76} P.~Cvitanovi\'c, % author={Predrag Cvitanovi\'c}, % title={Group Theory for {Feynman} Diagrams in Non-{Abelian} % Gauge Theories}, Phys.{} Rev.{} \textbf{D14} (1976) 1536. %%%\bibitem{Kleiss/etal:Color-Monte-Carlo} %%% \begin{dubious} %%% ``\texttt{Kleiss/etal:Color-Monte-Carlo}'' %%% \end{dubious} %\cite{Kilian:2012pz} \bibitem{Kilian:2012pz} W.~Kilian, T.~Ohl, J.~Reuter and C.~Speckner, %``QCD in the Color-Flow Representation,'' JHEP \textbf{1210} (2012) 022 [arXiv:1206.3700 [hep-ph]]. %%CITATION = doi:10.1007/JHEP10(2012)022;%% %37 citations counted in INSPIRE as of 23 Apr 2019 %\cite{Degrande:2011ua} \bibitem{Degrande:2011ua} C.~Degrande, C.~Duhr, B.~Fuks, D.~Grellscheid, O.~Mattelaer and T.~Reiter, %``UFO - The Universal FeynRules Output,'' Comput.{} Phys.{} Commun.{} \textbf{183} (2012), 1201-1214 doi:10.1016/j.cpc.2012.01.022 [arXiv:1108.2040 [hep-ph]]. \end{thebibliography} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \appendix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Autotools} \label{sec:autotools} \module{config} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Textual Options} \label{sec:options} \module{options} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Progress Reports} \label{sec:progress} \module{progress} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{More on Filenames} \label{sec:thoFilename} \module{thoFilename} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Cache Files} \label{sec:cache} \module{cache} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{More On Lists} \label{sec:tholist} \module{thoList} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{More On Arrays} \label{sec:thoarray} \module{thoArray} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{More On Strings} \label{sec:thostring} \module{thoString} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Polymorphic Maps} \label{sec:pmap} From~\cite{Ohl:LOTR}. \module{pmap} \module{partial} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Tries} \label{sec:trie} From~\cite{Okasaki:1998:book}, extended for~\cite{Ohl:LOTR}. \module{trie} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Tensor Products} \label{sec:product} From~\cite{Ohl:LOTR}. \module{product} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{(Fiber) Bundles} \label{sec:bundle} \module{bundle} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Power Sets} \label{sec:powSet} \module{powSet} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Combinatorics} \label{sec:combinatorics} \module{combinatorics} \module{permutation} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Partitions} \label{sec:partition} \module{partition} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Young Diagrams and Tableaux} \label{sec:young} \module{young} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Trees} \label{sec:tree} From~\cite{Ohl:bocages}: Trees with one root admit a straightforward recursive definition \begin{equation} \label{eq:trees} T(N,L) = L \cup N\times T(N,L)\times T(N,L) \end{equation} that is very well adapted to mathematical reasoning. Such recursive definitions are useful because they allow us to prove properties of elements by induction \begin{multline} \label{eq:tree-induction} \forall l\in L: p(l) \land (\forall n\in N: \forall t_1,t_2\in T(N,L): p(t_1) \land p(t_2) \Rightarrow p(n\times t_1\times t_2)) \\ \Longrightarrow \forall t\in T(N,L): p(t) \end{multline} i.\,e.~establishing a property for all leaves and showing that a node automatically satisfies the property if it is true for all children proves the property for \emph{all} trees. This induction is of course modelled after standard mathematical induction \begin{equation} p(1) \land (\forall n\in \mathbf{N}: p(n) \Rightarrow p(n+1)) \Longrightarrow \forall n\in \mathbf{N}: p(n) \end{equation} The recursive definition~(\ref{eq:trees}) is mirrored by the two tree construction functions\footnote{To make the introduction more accessible to non-experts, I avoid the `curried' notation for functions with multiple arguments and use tuples instead. The actual implementation takes advantage of curried functions, however. Experts can read $\alpha\to\beta\to\gamma$ for $\alpha\times\beta\to\gamma$.} \begin{subequations} \begin{align} \ocwlowerid{leaf}:\;& \nu\times\lambda \to(\nu,\lambda) T \\ \ocwlowerid{node}:\;& \nu\times(\nu,\lambda)T \times(\nu,\lambda)T \to(\nu,\lambda)T \end{align} \end{subequations} Renaming leaves and nodes leaves the structure of the tree invariant. Therefore, morphisms~$L\to L'$ and~$N\to N'$ of the sets of leaves and nodes induce natural homomorphisms~$T(N,L)\to T(N',L')$ of trees \begin{equation} \ocwlowerid{map}:\; (\nu\to\nu')\times(\lambda\to\lambda') \times(\nu,\lambda)T \to(\nu',\lambda') T \end{equation} The homomorphisms constructed by \ocwlowerid{map} are trivial, but ubiquitous. More interesting are the morphisms \begin{equation} \begin{aligned} \ocwlowerid{fold}:\;& (\nu\times\lambda\to\alpha) \times(\nu\times\alpha\times\alpha\to\alpha) \times(\nu,\lambda)T \to\alpha \\ & (f_1,f_2,l\in L) \mapsto f_1(l) \\ & (f_1,f_2,(n,t_1,t_2)) \mapsto f_2(n,\ocwlowerid{fold}(f_1,f_2,t_1), \ocwlowerid{fold}(f_1,f_2,t_2)) \end{aligned} \end{equation} and \begin{equation} \begin{aligned} \ocwlowerid{fan}:\;& (\nu\times\lambda\to\{\alpha\}) \times(\nu\times\alpha\times\alpha\to\{\alpha\}) \times(\nu,\lambda)T \to\{\alpha\} \\ & (f_1,f_2,l\in L) \mapsto f_1(l) \\ & (f_1,f_2,(n,t_1,t_2)) \mapsto f_2(n, \ocwlowerid{fold}(f_1,f_2,t_1) \otimes\ocwlowerid{fold}(f_1,f_2,t_2)) \end{aligned} \end{equation} where the tensor product notation means that~$f_2$ is applied to all combinations of list members in the argument: \begin{equation} \phi(\{x\}\otimes \{y\}) = \left\{ \phi(x,y) | x\in\{x\} \land y\in\{y\} \right\} \end{equation} But note that due to the recursive nature of trees, \ocwlowerid{fan} is \emph{not} a morphism from $T(N,L)$ to $T(N\otimes N,L)$.\par If we identify singleton sets with their members, \ocwlowerid{fold} could be viewed as a special case of \ocwlowerid{fan}, but that is probably more confusing than helpful. Also, using the special case~$\alpha=(\nu',\lambda')T$, the homomorphism \ocwlowerid{map} can be expressed in terms of \ocwlowerid{fold} and the constructors \begin{equation} \begin{aligned} \ocwlowerid{map}:\;& (\nu\to\nu')\times(\lambda\to\lambda') \times(\nu,\lambda)T \to(\nu',\lambda')T \\ &(f,g,t) \mapsto \ocwlowerid{fold} (\ocwlowerid{leaf}\circ (f\times g), \ocwlowerid{node}\circ (f\times\ocwlowerid{id} \times\ocwlowerid{id}), t) \end{aligned} \end{equation} \ocwlowerid{fold} is much more versatile than \ocwlowerid{map}, because it can be used with constructors for other tree representations to translate among different representations. The target type can also be a mathematical expression. This is used extensively below for evaluating Feynman diagrams.\par Using \ocwlowerid{fan} with~$\alpha=(\nu',\lambda')T$ can be used to construct a multitude of homomorphic trees. In fact, below it will be used extensively to construct all Feynman diagrams~$\{(\nu,\{p_1,\ldots,p_n\})T\}$ of a given topology~$t\in (\emptyset,\{1,\ldots,n\})T$. \begin{dubious} The physicist in me guesses that there is another morphism of trees that is related to \ocwlowerid{fan} like a Lie-algebra is related to the it's Lie-group. I have not been able to pin it down, but I guess that it is a generalization of \ocwlowerid{grow} below. \end{dubious} \module{tree} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Dependency Trees} \label{sec:tree2} \module{tree2} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Consistency Checks} \label{sec:count} \application{count} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Complex Numbers} \label{sec:complex} \module{complex} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Algebra} \label{sec:algebra} \module{algebra} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Simple Linear Algebra} \label{sec:linalg} \module{linalg} %application{test_linalg} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Partial Maps} \label{sec:partial} \module{partial} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Talk To The WHiZard \ldots} \label{sec:whizard_tool} Talk to~\cite{Kilian:WHIZARD}. \begin{dubious} Temporarily disabled, until, we implement some conditional weaving\ldots \end{dubious} %application{whizard_tool} %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% \chapter{Widget Library and Class Hierarchy for O'Giga} %%% \label{sec:thogtk} %%% {\itshape NB: The code in this chapter \emph{must} be compiled with %%% \verb+-labels+, since \verb+lablgtk+ doesn't appear to work in classic mode.} %%% \begin{dubious} %%% Keep in mind that \texttt{ocamlweb} doesn't work properly with %%% O'Caml~3 yet. The colons in label declarations are typeset with %%% erroneous white space. %%% \end{dubious} %%% %%% \section{Architecture} %%% In \texttt{lablgtk}, O'Caml objects are typically constructed in %%% parallel to constructors for \texttt{GTK+} widgets. The objects %%% provide inheritance and all that, while the constructors implement the %%% semantics. %%% %%% \subsection{Inheritance vs.~Aggregation} %%% We have two mechanisms for creating new widgets: inheritance and %%% aggregation. Inheritance makes it easy to extend a given widget with %%% new methods or to combine orthogonal widgets (\emph{multiple %%% inheritance}). Aggregation is more suitable for combining %%% non-orthogonal widgets (e.\,g.~multiple instances of the same widget). %%% %%% The problem with inheritance in \texttt{lablgtk} is, that it is a %%% \emph{bad} idea to implement the semantics in the objects. In a %%% multi-level inheritance hierarchy, O'Caml can evaluate class functions %%% more than once. Since functions accessing \texttt{GTK+} change the %%% state of \texttt{GTK+}, we could accidentally violate invariants. %%% Therefore inheritance forces us to use the two-tiered approach of %%% \texttt{lablgtk} ourselves. It is not really complicated, but tedious %%% and it appears to be a good idea to use aggregation whenever in doubt. %%% %%% Nevertheless, there are examples (like %%% \ocwupperid{ThoGButton.mutable\_button} below, where just one new %%% method is added), that cry out for inheritance for the benefit of the %%% application developer. %%% %%% \module{thoGWindow} %%% \module{thoGButton} %%% \module{thoGMenu} %%% \module{thoGDraw} %%% %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% \chapter{O'Mega Virtual Machine} %%% \label{sec:ovm} %%% \module{OVM} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{\texttt{Fortran} Libraries} \label{sec:fortran} \input{omegalib} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{raggedright} \ifpdf \chapter{Index} \let\origtwocolumn\twocolumn \def\twocolumn[#1]{\origtwocolumn}% This index has been generated automatically and might not be 100\%ly accurate. In particular, hyperlinks have been observed to be off by one page. \fi \input{index.tex} \end{raggedright} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \end{empfile} \end{fmffile} \end{document} \endinput Local Variables: mode:latex indent-tabs-mode:nil page-delimiter:"^%%%%%.*\n" End: Index: trunk/omega/src/color.ml =================================================================== --- trunk/omega/src/color.ml (revision 8848) +++ trunk/omega/src/color.ml (revision 8849) @@ -1,3212 +1,3586 @@ (* color.ml -- Copyright (C) 1999-2022 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) (* Avoid refering to [Pervasives.compare], because [Pervasives] will become [Stdlib.Pervasives] in O'Caml 4.07 and [Stdlib] in O'Caml 4.08. *) let pcompare = compare module type Test = sig val suite : OUnit.test val suite_long : OUnit.test end (* \thocwmodulesection{Quantum Numbers} *) type t = | Singlet | SUN of int | AdjSUN of int let conjugate = function | Singlet -> Singlet | SUN n -> SUN (-n) | AdjSUN n -> AdjSUN n let compare c1 c2 = match c1, c2 with | Singlet, Singlet -> 0 | Singlet, _ -> -1 | _, Singlet -> 1 | SUN n, SUN n' -> compare n n' | SUN _, AdjSUN _ -> -1 | AdjSUN _, SUN _ -> 1 | AdjSUN n, AdjSUN n' -> compare n n' module type Line = sig type t val conj : t -> t val equal : t -> t -> bool val to_string : t -> string end module type Cycles = sig type line type t = (line * line) list (* Contract the graph by connecting lines and return the number of cycles together with the contracted graph. \begin{dubious} The semantics of the contracted graph is not yet 100\%ly fixed. \end{dubious} *) val contract : t -> int * t (* The same as [contract], but returns only the number of cycles and raises [Open_line] when not all lines are closed. *) val count : t -> int exception Open_line (* Mainly for debugging \ldots *) val to_string : t -> string end module Cycles (L : Line) : Cycles with type line = L.t = struct type line = L.t type t = (line * line) list exception Open_line (* NB: The following algorithm for counting the cycles is quadratic since it performs nested scans of the lists. If this was a serious problem one could replace the lists of pairs by a [Map] and replace one power by a logarithm. *) let rec find_fst c_final c1 disc seen = function | [] -> ((L.conj c_final, c1) :: disc, List.rev seen) | (c1', c2') as c12' :: rest -> if L.equal c1 c1' then find_snd c_final (L.conj c2') disc [] (List.rev_append seen rest) else find_fst c_final c1 disc (c12' :: seen) rest and find_snd c_final c2 disc seen = function | [] -> ((L.conj c_final, L.conj c2) :: disc, List.rev seen) | (c1', c2') as c12' :: rest-> if L.equal c2' c2 then begin if L.equal c1' c_final then (disc, List.rev_append seen rest) else find_fst c_final (L.conj c1') disc [] (List.rev_append seen rest) end else find_snd c_final c2 disc (c12' :: seen) rest let consume = function | [] -> ([], []) | (c1, c2) :: rest -> find_snd (L.conj c1) (L.conj c2) [] [] rest let contract lines = let rec contract' acc disc = function | [] -> (acc, List.rev disc) | rest -> begin match consume rest with | [], rest' -> contract' (succ acc) disc rest' | disc', rest' -> contract' acc (List.rev_append disc' disc) rest' end in contract' 0 [] lines let count lines = match contract lines with | n, [] -> n | n, _ -> raise Open_line let to_string lines = String.concat "" (List.map (fun (c1, c2) -> "[" ^ L.to_string c1 ^ "," ^ L.to_string c2 ^ "]") lines) end (* \thocwmodulesection{Color Flows} *) module type Flow = sig type color type t = color list * color list val rank : t -> int val of_list : int list -> color val ghost : unit -> color val to_lists : t -> int list list val in_to_lists : t -> int list list val out_to_lists : t -> int list list val ghost_flags : t -> bool list val in_ghost_flags : t -> bool list val out_ghost_flags : t -> bool list type power = { num : int; den : int; power : int } type factor = power list val factor : t -> t -> factor val zero : factor module Test : Test end module Flow : Flow = struct (* All [int]s are non-zero! *) type color = | N of int | N_bar of int | SUN of int * int | Singlet | Ghost (* Incoming and outgoing, since we need to cross the incoming states. *) type t = color list * color list let rank cflow = 2 (* \thocwmodulesubsection{Constructors} *) let ghost () = Ghost let of_list = function | [0; 0] -> Singlet | [c; 0] -> N c | [0; c] -> N_bar c | [c1; c2] -> SUN (c1, c2) | _ -> invalid_arg "Color.Flow.of_list: num_lines != 2" let to_list = function | N c -> [c; 0] | N_bar c -> [0; c] | SUN (c1, c2) -> [c1; c2] | Singlet -> [0; 0] | Ghost -> [0; 0] let to_lists (cfin, cfout) = (List.map to_list cfin) @ (List.map to_list cfout) let in_to_lists (cfin, _) = List.map to_list cfin let out_to_lists (_, cfout) = List.map to_list cfout let ghost_flag = function | N _ | N_bar _ | SUN (_, _) | Singlet -> false | Ghost -> true let ghost_flags (cfin, cfout) = (List.map ghost_flag cfin) @ (List.map ghost_flag cfout) let in_ghost_flags (cfin, _) = List.map ghost_flag cfin let out_ghost_flags (_, cfout) = List.map ghost_flag cfout (* \thocwmodulesubsection{Evaluation} *) type power = { num : int; den : int; power : int } type factor = power list let zero = [] let count_ghosts1 colors = List.fold_left (fun acc -> function Ghost -> succ acc | _ -> acc) 0 colors let count_ghosts (fin, fout) = count_ghosts1 fin + count_ghosts1 fout type 'a square = | Square of 'a | Mismatch let conjugate = function | N c -> N_bar (-c) | N_bar c -> N (-c) | SUN (c1, c2) -> SUN (-c2, -c1) | Singlet -> Singlet | Ghost -> Ghost let cross_in (cin, cout) = cin @ (List.map conjugate cout) let cross_out (cin, cout) = (List.map conjugate cin) @ cout module C = Cycles (struct type t = int let conj = (~-) let equal = (=) let to_string = string_of_int end) (* Match lines in the color flows [f1] and [f2] after crossing the incoming states. This will be used to compute squared diagrams in [square] and [square2] below. *) let match_lines match1 match2 f1 f2 = let rec match_lines' acc f1' f2' = match f1', f2' with (* If we encounter an empty list, we're done --- unless the lengths don't match (which should never happen!): *) | [], [] -> Square (List.rev acc) | _ :: _, [] | [], _ :: _ -> Mismatch (* Handle matching \ldots *) | Ghost :: rest1, Ghost :: rest2 | Singlet :: rest1, Singlet :: rest2 -> match_lines' acc rest1 rest2 (* \ldots{} and mismatched ghosts and singlet gluons: *) | Ghost :: _, Singlet :: _ | Singlet :: _, Ghost :: _ -> Mismatch (* Ghosts and singlet gluons can't match anything else *) | (Ghost | Singlet) :: _, (N _ | N_bar _ | SUN (_, _)) :: _ | (N _ | N_bar _ | SUN (_, _)) :: _, (Ghost | Singlet) :: _ -> Mismatch (* Handle matching \ldots *) | N_bar c1 :: rest1, N_bar c2 :: rest2 | N c1 :: rest1, N c2 :: rest2 -> match_lines' (match1 c1 c2 acc) rest1 rest2 (* \ldots{} and mismatched $N$ or $\bar N$ states: *) | N _ :: _, N_bar _ :: _ | N_bar _ :: _, N _ :: _ -> Mismatch (* The $N$ and $\bar N$ don't match non-singlet gluons: *) | (N _ | N_bar _) :: _, SUN (_, _) :: _ | SUN (_, _) :: _, (N _ | N_bar _) :: _ -> Mismatch (* Now we're down to non-singlet gluons: *) | SUN (c1, c1') :: rest1, SUN (c2, c2') :: rest2 -> match_lines' (match2 c1 c1' c2 c2' acc) rest1 rest2 in match_lines' [] (cross_out f1) (cross_out f2) (* NB: in WHIZARD versions before 3.0, the code for [match_lines] contained a bug in the pattern matching of [Singlet], [N], [N_bar] and [SUN] states, because they all were represented as [SUN (c1, c2)], only distinguished by the numeric conditions [c1 = 0] and/or [c2 = 0]. This prevented the use of exhaustiveness checking and introduced a subtle dependence on the pattern order. *) let square f1 f2 = match_lines (fun c1 c2 pairs -> (c1, c2) :: pairs) (fun c1 c1' c2 c2' pairs -> (c1', c2') :: (c1, c2) :: pairs) f1 f2 (*i let square f1 f2 = let ll2s ll = String.concat "; " (List.map (ThoList.to_string string_of_int) ll) and lp2s lp = String.concat "; " (List.map (fun (c1, c2) -> string_of_int c1 ^ ", " ^ string_of_int c2) lp) in Printf.eprintf "square ([%s], [%s]) ([%s], [%s]) = " (ll2s (in_to_lists f1)) (ll2s (out_to_lists f1)) (ll2s (in_to_lists f2)) (ll2s (out_to_lists f2)); let res = square f1 f2 in begin match res with | Mismatch -> Printf.eprintf "Mismatch!\n" | Square f12 -> Printf.eprintf "Square [%s]\n" (lp2s f12) end; res i*) (* In addition to counting closed color loops, we also need to count closed gluon loops. Fortunately, we can use the same algorithm on a different data type, provided it doesn't require all lines to be closed. *) module C2 = Cycles (struct type t = int * int let conj (c1, c2) = (- c2, - c1) let equal (c1, c2) (c1', c2') = c1 = c1' && c2 = c2' let to_string (c1, c2) = "(" ^ string_of_int c1 ^ "," ^ string_of_int c2 ^ ")" end) let square2 f1 f2 = match_lines (fun c1 c2 pairs -> pairs) (fun c1 c1' c2 c2' pairs -> ((c1, c1'), (c2, c2')) :: pairs) f1 f2 (* $\ocwlowerid{int\_power}: n\, p \to n^p$ for integers is missing from [Pervasives]! *) let int_power n p = let rec int_power' acc i = if i < 0 then invalid_arg "int_power" else if i = 0 then acc else int_power' (n * acc) (pred i) in int_power' 1 p (* Instead of implementing a full fledged algebraic evaluator, let's simply expand the binomial by hand: \begin{equation} \left(\frac{N_C^2-2}{N_C^2}\right)^n = \sum_{i=0}^n \binom{n}{i} (-2)^i N_C^{-2i} \end{equation} *) (* NB: Any result of [square] other than [Mismatch] guarantees [count_ghosts f1 = count_ghosts f2]. *) let factor f1 f2 = match square f1 f2, square2 f1 f2 with | Mismatch, _ | _, Mismatch -> [] | Square f12, Square f12' -> let num_cycles = C.count f12 and num_cycles2, disc = C2.contract f12' and num_ghosts = count_ghosts f1 in (*i Printf.eprintf "f12 = %s -> #loops = %d\n" (C.to_string f12) num_cycles; Printf.eprintf "f12' = %s -> #loops = %d, disc = %s\n" (C2.to_string f12') num_cycles2 (C2.to_string disc); flush stderr; i*) List.map (fun i -> let parity = if num_ghosts mod 2 = 0 then 1 else -1 and power = num_cycles - num_ghosts in let coeff = int_power (-2) i * Combinatorics.binomial num_cycles2 i and power2 = - 2 * i in { num = parity * coeff; den = 1; power = power + power2 }) (ThoList.range 0 num_cycles2) module Test : Test = struct open OUnit (* Here and elsewhere, we have to resist the temptation to define these tests as functions with an additional argument [()] in the hope to avoid having to package them into an explicit thunk [fun () -> eq v1 v2] in order to delay evaluation. It turns out that the runtime would then sometimes evaluate the argument [v1] or [v2] even \emph{before} the test is run. For pure functions, there is no difference, but the compiler appears to treat explicit thunks specially. \begin{dubious} I haven't yet managed to construct a small demonstrator to find out in which circumstances the premature evaluation happens. \end{dubious} *) let suite_square = "square" >::: [ "square ([], []) ([], [])" >:: (fun () -> assert_equal (Square []) (square ([], []) ([], []))); "square ([3], [3; 0]) ([3], [3; 0])" >:: (fun () -> assert_equal (Square [(-1, -1); (1, 1)]) (square ([N 1], [N 1; Singlet]) ([N 1], [N 1; Singlet]))); "square ([0], [3; -3]) ([0], [3; -3])" >:: (fun () -> assert_equal (Square [(1, 1); (-1, -1)]) (square ([Singlet], [N 1; N_bar (-1)]) ([Singlet], [N 1; N_bar (-1)]))); "square ([3], [3; 0]) ([0], [3; -3])" >:: (fun () -> assert_equal Mismatch (square ([N 1], [N 1; Singlet]) ([Singlet], [N 1; N_bar (-1)]))); "square ([3; 8], [3]) ([3; 8], [3])" >:: (fun () -> assert_equal (Square [-1, -1; 1, 1; -2, -2; 2, 2]) (square ([N 1; SUN (2, -1)], [N 2]) ([N 1; SUN (2, -1)], [N 2]))) ] let suite = "Color.Flow" >::: [suite_square] let suite_long = "Color.Flow long" >::: [] end end (* later: *) module General_Flow = struct type color = | Lines of int list | Ghost of int type t = color list * color list let rank_default = 2 (* Standard model *) let rank cflow = try begin match List.hd cflow with | Lines lines -> List.length lines | Ghost n_lines -> n_lines end with | _ -> rank_default end (* \thocwmodulesection{Vertex Color Flows} *) (* \newcommand{\setupFourAmp}{% \fmfleft{i1,i2} \fmfright{o1,o2} \fmf{phantom}{i1,v1,i2} \fmf{phantom}{o2,v2,o1} \fmf{phantom}{v1,v2} \fmffreeze} \fmfcmd{% numeric joindiameter; joindiameter := 7thick;} \fmfcmd{% vardef sideways_at (expr d, p, frac) = save len; len = length p; (point frac*len of p) shifted ((d,0) rotated (90 + angle direction frac*len of p)) enddef; secondarydef p sideways d = for frac = 0 step 0.01 until 0.99: sideways_at (d, p, frac) .. endfor sideways_at (d, p, 1) enddef; secondarydef p choptail d = subpath (ypart (fullcircle scaled d shifted (point 0 of p) intersectiontimes p), infinity) of p enddef; secondarydef p choptip d = reverse ((reverse p) choptail d) enddef; secondarydef p pointtail d = fullcircle scaled d shifted (point 0 of p) intersectionpoint p enddef; secondarydef p pointtip d = (reverse p) pointtail d enddef; secondarydef pa join pb = pa choptip joindiameter .. pb choptail joindiameter enddef; vardef cyclejoin (expr p) = subpath (0.5*length p, infinity) of p join subpath (0, 0.5*length p) of p .. cycle enddef;} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \fmfcmd{% style_def double_line_arrow expr p = save pi, po; path pi, po; pi = reverse (p sideways thick); po = p sideways -thick; cdraw pi; cdraw po; cfill (arrow (subpath (0, 0.9 length pi) of pi)); cfill (arrow (subpath (0, 0.9 length po) of po)); enddef;} \fmfcmd{% style_def double_line_arrow_beg expr p = save pi, po, pc; path pi, po, pc; pc = p choptail 7thick; pi = reverse (pc sideways thick); po = pc sideways -thick; cdraw pi .. p pointtail 5thick .. po; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_line_arrow_end expr p = save pi, po, pc; path pi, po, pc; pc = p choptip 7thick; pi = reverse (pc sideways thick); po = pc sideways -thick; cdraw po .. p pointtip 5thick .. pi; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_line_arrow_both expr p = save pi, po, pc; path pi, po, pc; pc = p choptip 7thick choptail 7thick; pi = reverse (pc sideways thick); po = pc sideways -thick; cdraw po .. p pointtip 5thick .. pi .. p pointtail 5thick .. cycle; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_arrow_parallel expr p = save pi, po; path pi, po; pi = p sideways thick; po = p sideways -thick; save li, lo; li = length pi; lo = length po; cdraw pi; cdraw po; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_arrow_crossed_beg expr p = save lp; lp = length p; save pi, po; path pi, po; pi = p sideways thick; po = p sideways -thick; save li, lo; li = length pi; lo = length po; cdraw subpath (0, 0.1 li) of pi .. subpath (0.3 lo, lo) of po; cdraw subpath (0, 0.1 lo) of po .. subpath (0.3 li, li) of pi; cfill (arrow pi); cfill (arrow po); enddef;} \fmfcmd{% style_def double_arrow_crossed_end expr p = save lp; lp = length p; save pi, po; path pi, po; pi = p sideways thick; po = p sideways -thick; save li, lo; li = length pi; lo = length po; cdraw subpath (0, 0.7 li) of pi .. subpath (0.9 lo, lo) of po; cdraw subpath (0, 0.7 lo) of po .. subpath (0.9 li, li) of pi; cfill (arrow pi); cfill (arrow po); enddef;} *) module Q = Algebra.Q module QC = Algebra.QC module type Arrow = sig type endpoint type tip = endpoint type tail = endpoint type ghost = endpoint val position : endpoint -> int val relocate : (int -> int) -> endpoint -> endpoint type ('tail, 'tip, 'ghost) t = | Arrow of 'tail * 'tip | Ghost of 'ghost + | Epsilon of 'tip list + | Epsilon_bar of 'tail list type free = (tail, tip, ghost) t type factor val free_to_string : free -> string val factor_to_string : factor -> string val map : (endpoint -> endpoint) -> free -> free val to_left_factor : (endpoint -> bool) -> free -> factor val to_right_factor : (endpoint -> bool) -> free -> factor val of_factor : factor -> free val is_free : factor -> bool val negatives : free -> endpoint list val is_ghost : free -> bool + val is_tadpole : factor -> bool type merge = | Match of factor + | Determinant of factor list list * factor list list | Ghost_Match | Loop_Match | Mismatch | No_Match val merge : factor -> factor -> merge val tee : int -> free -> free list val dir : int -> int -> free -> int val single : endpoint -> endpoint -> free val double : endpoint -> endpoint -> free list val ghost : endpoint -> free module Infix : sig val (=>) : int -> int -> free val (==>) : int -> int -> free list val (<=>) : int -> int -> free list val (>=>) : int * int -> int -> free val (=>>) : int -> int * int -> free val (>=>>) : int * int -> int * int -> free val (??) : int -> free end + val epsilon : int list -> free + val epsilon_bar : int list -> free val chain : int list -> free list val cycle : int list -> free list module Test : Test val pp_free : Format.formatter -> free -> unit val pp_factor : Format.formatter -> factor -> unit end module Arrow : Arrow = struct type endpoint = | I of int | M of int * int let position = function | I i -> i | M (i, _) -> i let relocate f = function | I i -> I (f i) | M (i, n) -> M (f i, n) type tip = endpoint type tail = endpoint type ghost = endpoint (* Note that in the case of double lines for the adjoint representation the \emph{same} [endpoint] appears twice: once as a [tip] and once as a [tail]. If we want to multiply two factors by merging arrows with matching [tip] and [tail], we must make sure that the [tip] is from one factor and the [tail] from the other factor. *) (* The [Free] variant contains positive indices as well as negative indices that don't appear on both sides and will be summed in a later product. [SumL] and [SumR] indices appear on both sides. *) type 'a index = | Free of 'a | SumL of 'a | SumR of 'a + let is_free_index = function + | Free _ -> true + | SumL _ | SumR _ -> false + type ('tail, 'tip, 'ghost) t = | Arrow of 'tail * 'tip | Ghost of 'ghost + | Epsilon of 'tip list + | Epsilon_bar of 'tail list type free = (tail, tip, ghost) t type factor = (tail index, tip index, ghost index) t let endpoint_to_string = function | I i -> string_of_int i | M (i, n) -> Printf.sprintf "%d.%d" i n let index_to_string = function | Free i -> endpoint_to_string i | SumL i -> endpoint_to_string i ^ "L" | SumR i -> endpoint_to_string i ^ "R" let to_string i2s = function | Arrow (tail, tip) -> Printf.sprintf "%s>%s" (i2s tail) (i2s tip) | Ghost ghost -> Printf.sprintf "{%s}" (i2s ghost) + | Epsilon tips -> Printf.sprintf ">>>%s" (ThoList.to_string i2s tips) + | Epsilon_bar tails -> Printf.sprintf "<<<%s" (ThoList.to_string i2s tails) let free_to_string = to_string endpoint_to_string let factor_to_string = to_string index_to_string let index_matches i1 i2 = match i1, i2 with | SumL i1, SumR i2 | SumR i1, SumL i2 -> i1 = i2 | _ -> false let map f = function | Arrow (tail, tip) -> Arrow (f tail, f tip) | Ghost ghost -> Ghost (f ghost) + | Epsilon tips -> Epsilon (List.map f tips) + | Epsilon_bar tails -> Epsilon_bar (List.map f tails) let free_index = function | Free i -> i | SumL i -> invalid_arg "Color.Arrow.free_index: leftover LHS summation" | SumR i -> invalid_arg "Color.Arrow.free_index: leftover RHS summation" let to_left_index is_sum i = if is_sum i then SumL i else Free i let to_right_index is_sum i = if is_sum i then SumR i else Free i let to_left_factor is_sum = map (to_left_index is_sum) let to_right_factor is_sum = map (to_right_index is_sum) let of_factor = map free_index let negatives = function | Arrow (tail, tip) -> if position tail < 0 then if position tip < 0 then [tail; tip] else [tail] else if position tip < 0 then [tip] else [] | Ghost ghost -> if position ghost < 0 then [ghost] else [] + | Epsilon tips -> List.filter (fun tip -> position tip < 0) tips + | Epsilon_bar tails -> List.filter (fun tail -> position tail < 0) tails let is_free = function | Arrow (Free _, Free _) | Ghost (Free _) -> true - | _ -> false + | Arrow (_, _) | Ghost _ -> false + | Epsilon tips -> List.for_all is_free_index tips + | Epsilon_bar tails -> List.for_all is_free_index tails let is_ghost = function | Ghost _ -> true | Arrow _ -> false - + | Epsilon _ | Epsilon_bar _ -> false + let single tail tip = Arrow (tail, tip) let double a b = if a = b then [single a b] else [single a b; single b a] let ghost g = Ghost g + module Infix = + struct + let ( => ) i j = single (I i) (I j) + let ( ==> ) i j = [i => j] + let ( <=> ) i j = double (I i) (I j) + let ( >=> ) (i, n) j = single (M (i, n)) (I j) + let ( =>> ) i (j, m) = single (I i) (M (j, m)) + let ( >=>> ) (i, n) (j, m) = single (M (i, n)) (M (j, m)) + let ( ?? ) i = ghost (I i) + end + + open Infix + +(* Split [a_list] at the first element equal to [a] according + to [eq]. Return the reversed first part and the rest as a + pair and wrap it in [Some]. Return [None] if there is no match. *) + let take_first_match_opt ?(eq=(=)) a a_list = + let rec take_first_match_opt' rev_head = function + | [] -> None + | elt :: tail -> + if eq elt a then + Some (rev_head, tail) + else + take_first_match_opt' (elt :: rev_head) tail in + take_first_match_opt' [] a_list + +(* Split [a_list] and [b_list] at the first element equal according + to [eq]. Return the reversed first part and the rest of each + as a pair of pairs wrap it in [Some]. + Return [None] if there is no match. *) + let take_first_matching_pair_opt ?(eq=(=)) a_list b_list = + let rec take_first_matching_pair_opt' rev_a_head = function + | [] -> None + | a :: a_tail -> + begin match take_first_match_opt ~eq a b_list with + | Some (rev_b_head, b_tail) -> + Some ((rev_a_head, a_tail), (rev_b_head, b_tail)) + | None -> + take_first_matching_pair_opt' (a :: rev_a_head) a_tail + end in + take_first_matching_pair_opt' [] a_list + +(* Replace the first occurence of an element equal to [a] according + to [eq] in [a_list] by [a'] and wrap the new list in [Some]. + Return [None] if there is no match. *) + let replace_first_opt ?(eq=(=)) a a' a_list = + match take_first_match_opt ~eq a a_list with + | Some (rev_head, tail) -> Some (List.rev_append rev_head (a' :: tail)) + | None -> None + let tee a = function | Arrow (tail, tip) -> [Arrow (tail, I a); Arrow (I a, tip)] - | Ghost _ -> [] + | Ghost _ as g -> [g] + | Epsilon _ -> invalid_arg "Arrow.tee not defined for Epsilon" + | Epsilon_bar _ -> invalid_arg "Arrow.tee not defined for Epsilon_bar" let dir i j = function | Arrow (tail, tip) -> let tail = position tail and tip = position tip in if tip = i && tail = j then - 1 + 1 else if tip = j && tail = i then -1 else invalid_arg "Arrow.dir" - | Ghost _ -> 0 + | Ghost _ | Epsilon _ | Epsilon_bar _ -> 0 type merge = | Match of factor + | Determinant of factor list list * factor list list | Ghost_Match | Loop_Match | Mismatch | No_Match - let merge arrow1 arrow2 = +(* When computing + \begin{equation} + \epsilon_{ki_1i_2\cdots i_n} \bar\epsilon^{kj_1j_2\cdots j_n} + = \sum_{\sigma} (-1)^{\epsilon(\sigma)} + \delta_{i_1}^{\sigma(j_1)} + \delta_{i_2}^{\sigma(j_2)} + \cdots + \delta_{i_n}^{\sigma(j_n)}\,, + \end{equation} + we must keep track of the position of summation indices. + We can use the fact that cyclic permutations are even for + $\epsilon$-tensors with an odd number of indices, corresponding + to $n$ even and odd otherwise. *) + + let fuse_epsilons tails tips = + match take_first_matching_pair_opt ~eq:index_matches tails tips with + | None -> No_Match + | Some ((rev_tails_head, tails_tail), (rev_tips_head, tips_tail)) -> + let tails = tails_tail @ List.rev rev_tails_head + and tips = tips_tail @ List.rev rev_tips_head in + let num_tails = List.length tails + and num_tips = List.length tips in + if num_tails <> num_tips then + invalid_arg + (Printf.sprintf + "Color.Arrow.fuse_epsilons: length mismatch %d <> %d" + (succ num_tails) (succ num_tips)) + else + let is_odd n = n mod 2 <> 0 in + let flip = + is_odd num_tips && + is_odd (List.length rev_tails_head - List.length rev_tips_head) in + let even_tips = Combinatorics.permute_even tips + and odd_tips = Combinatorics.permute_odd tips in + let even = List.rev_map (List.rev_map2 single tails) even_tips + and odd = List.rev_map (List.rev_map2 single tails) odd_tips in + if flip then + Determinant (odd, even) + else + Determinant (even, odd) + + let merge' arrow1 arrow2 = match arrow1, arrow2 with | Ghost g1, Ghost g2 -> if index_matches g1 g2 then Ghost_Match else No_Match | Arrow (tail, tip), Ghost g | Ghost g, Arrow (tail, tip) -> if index_matches g tail || index_matches g tip then Mismatch else No_Match | Arrow (tail, tip), Arrow (tail', tip') -> if index_matches tip tail' then if index_matches tip' tail then Loop_Match else Match (Arrow (tail, tip')) else if index_matches tip' tail then Match (Arrow (tail', tip)) else No_Match + | Arrow (tail, tip), Epsilon tips | Epsilon tips, Arrow (tail, tip) -> + begin match replace_first_opt ~eq:index_matches tail tip tips with + | None -> No_Match + | Some tips -> Match (Epsilon tips) + end + | Arrow (tail, tip), Epsilon_bar tails | Epsilon_bar tails, Arrow (tail, tip) -> + begin match replace_first_opt ~eq:index_matches tip tail tails with + | None -> No_Match + | Some tails -> Match (Epsilon_bar tails) + end + | Epsilon tips, Ghost g | Ghost g, Epsilon tips -> + if List.exists (index_matches g) tips then + Mismatch + else + No_Match + | Epsilon_bar tails, Ghost g | Ghost g, Epsilon_bar tails -> + if List.exists (index_matches g) tails then + Mismatch + else + No_Match + | Epsilon _, Epsilon _ | Epsilon_bar _, Epsilon_bar _ -> + No_Match + | Epsilon tips, Epsilon_bar tails | Epsilon_bar tails, Epsilon tips -> + fuse_epsilons tails tips - module Infix = - struct - let ( => ) i j = single (I i) (I j) - let ( ==> ) i j = [i => j] - let ( <=> ) i j = double (I i) (I j) - let ( >=> ) (i, n) j = single (M (i, n)) (I j) - let ( =>> ) i (j, m) = single (I i) (M (j, m)) - let ( >=>> ) (i, n) (j, m) = single (M (i, n)) (M (j, m)) - let ( ?? ) i = ghost (I i) - end +(* As an optimization, don't attempt to merge neither of the arrows + contains a summation index and return immediately. *) - open Infix + let merge arrow1 arrow2 = + if is_free arrow1 || is_free arrow2 then + No_Match + else + merge' arrow1 arrow2 + + let merge_to_string = function + | Match factor -> + Printf.sprintf "Match (%s)" (factor_to_string factor) + | Determinant (even, odd) -> + Printf.sprintf + "Det (%s, %s)" + (ThoList.to_string (ThoList.to_string factor_to_string) even) + (ThoList.to_string (ThoList.to_string factor_to_string) odd) + | Ghost_Match -> "Ghost" + | Loop_Match -> "Loop" + | Mismatch -> "Mismatch" + | No_Match -> "No_Match" + + let logging_merge arrow1 arrow2 = + let result = merge arrow1 arrow2 in + Printf.eprintf + "merge %s with %s ==> %s\n" + (factor_to_string arrow1) + (factor_to_string arrow2) + (merge_to_string result); + result + + let is_tadpole = function + | Arrow (tail, tip) -> + index_matches tail tip + | _ -> false + +(*i + let merge = logging_merge +i*) + let epsilon = function + | [] -> invalid_arg "Color.Arrow.epsilon []" + | [_] -> invalid_arg "Color.Arrow.epsilon lone index" + | tips -> + Epsilon (List.map (fun tip -> I tip) tips) + + let epsilon_bar = function + | [] -> invalid_arg "Color.Arrow.epsilon []" + | [_] -> invalid_arg "Color.Arrow.epsilon lone index" + | tails -> + Epsilon_bar (List.map (fun tail -> I tail) tails) (* Composite Arrows. *) let rec chain = function | [] -> [] | [a] -> [a => a] | [a; b] -> [a => b] | a :: (b :: _ as rest) -> (a => b) :: chain rest let rec cycle' a = function | [] -> [a => a] | [b] -> [b => a] | b :: (c :: _ as rest) -> (b => c) :: cycle' a rest let cycle = function | [] -> [] | a :: _ as a_list -> cycle' a a_list module Test : Test = struct open OUnit let suite_chain = "chain" >::: - - [ "chain []" >:: - (fun () -> - assert_equal [] (chain [])); - - "chain [1]" >:: - (fun () -> - assert_equal [1 => 1] (chain [1])); - - "chain [1;2]" >:: - (fun () -> - assert_equal [1 => 2] (chain [1; 2])); - - "chain [1;2;3]" >:: - (fun () -> - assert_equal [1 => 2; 2 => 3] (chain [1; 2; 3])); - - "chain [1;2;3;4]" >:: - (fun () -> - assert_equal [1 => 2; 2 => 3; 3 => 4] (chain [1; 2; 3; 4])) ] + [ "[]" >:: (fun () -> assert_equal [] (chain [])); + "[1]" >:: (fun () -> assert_equal [1 => 1] (chain [1])); + "[1;2]" >:: (fun () -> assert_equal [1 => 2] (chain [1; 2])); + "[1;2;3]" >:: (fun () -> assert_equal [1 => 2; 2 => 3] (chain [1; 2; 3])); + "[1;2;3;4]" >:: (fun () -> assert_equal [1 => 2; 2 => 3; 3 => 4] (chain [1; 2; 3; 4])) ] let suite_cycle = "cycle" >::: + [ "[]" >:: (fun () -> assert_equal [] (cycle [])); + "[1]" >:: (fun () -> assert_equal [1 => 1] (cycle [1])); + "[1;2]" >:: (fun () -> assert_equal [1 => 2; 2 => 1] (cycle [1; 2])); + "[1;2;3]" >:: (fun () -> assert_equal [1 => 2; 2 => 3; 3 => 1] (cycle [1; 2; 3])); + + "[1;2;3;4]" >:: (fun () -> assert_equal [1 => 2; 2 => 3; 3 => 4; 4 => 1] (cycle [1; 2; 3; 4])) ] + + let suite_take = + "take" >::: + [ "1 []" >:: (fun () -> assert_equal None (take_first_match_opt 1 [])); + "1 [1]" >:: (fun () -> assert_equal (Some ([], [])) (take_first_match_opt 1 [1])); + "1 [2;3;4]" >:: (fun () -> assert_equal None (take_first_match_opt 1 [2;3;4])); + "1 [1;2;3]" >:: (fun () -> assert_equal (Some ([], [2;3])) (take_first_match_opt 1 [1;2;3])); + "2 [1;2;3]" >:: (fun () -> assert_equal (Some ([1], [3])) (take_first_match_opt 2 [1;2;3])); + "3 [1;2;3]" >:: (fun () -> assert_equal (Some ([2;1], [])) (take_first_match_opt 3 [1;2;3])) ] + + let suite_take2 = + "take2" >::: + [ "[] []" >:: + (fun () -> assert_equal None (take_first_matching_pair_opt [] [])); + + "[] [1;2;3]" >:: + (fun () -> assert_equal None (take_first_matching_pair_opt [] [1;2;3])); + + "[1] [2;3;4]" >:: + (fun () -> assert_equal None (take_first_matching_pair_opt [1] [2;3;4])); + + "[2;3;4] [1]" >:: + (fun () -> assert_equal None (take_first_matching_pair_opt [2;3;4] [1])); + + "[1;2;3] [4;5;6;7]" >:: + (fun () -> assert_equal None (take_first_matching_pair_opt [1;2;3] [4;5;6;7])); + + "[1] [1;2;3]" >:: + (fun () -> + assert_equal + (Some (([],[]), ([],[2;3]))) + (take_first_matching_pair_opt [1] [1;2;3])); + + "[1;2;3] [1;20;30]" >:: + (fun () -> + assert_equal + (Some (([],[2;3]), ([],[20;30]))) + (take_first_matching_pair_opt [1;2;3] [1;20;30])); + + "[1;2;3;4;5;6] [10;20;4;30;40]" >:: + (fun () -> + assert_equal + (Some (([3;2;1],[5;6]), ([20;10],[30;40]))) + (take_first_matching_pair_opt [1;2;3;4;5;6] [10;20;4;30;40])) ] + + let suite_replace = + "replace" >::: + [ "1 10 []" >:: (fun () -> assert_equal None (replace_first_opt 1 2 [])); + "1 10 [1]" >:: (fun () -> assert_equal (Some [10]) (replace_first_opt 1 10 [1])); + "1 [2;3;4]" >:: (fun () -> assert_equal None (replace_first_opt 1 10 [2;3;4])); + "1 [1;2;3]" >:: (fun () -> assert_equal (Some [10;2;3]) (replace_first_opt 1 10 [1;2;3])); + "2 [1;2;3]" >:: (fun () -> assert_equal (Some [1;10;3]) (replace_first_opt 2 10 [1;2;3])); + "3 [1;2;3]" >:: (fun () -> assert_equal (Some [1;2;10]) (replace_first_opt 3 10 [1;2;3])) ] + + let determinant_to_string = function + | Determinant (even, odd) -> + Printf.sprintf + "Determinant (even = %s, odd = %s)" + (ThoList.to_string (ThoList.to_string factor_to_string) even) + (ThoList.to_string (ThoList.to_string factor_to_string) odd) + | _ -> "not a Determinant" + + let make_determinant even odd = + let make_free_single (tail, tip) = + single (Free (I tail)) (Free (I tip)) in + Determinant + (List.map (List.map make_free_single) even, + List.map (List.map make_free_single) odd) + + let canonicalize_determinant = function + | Determinant (even, odd) -> + Determinant + (List.sort pcompare (List.map (List.sort pcompare) even), + List.sort pcompare (List.map (List.sort pcompare) odd)) + | other -> other + + let merge_epsilon_pair eps eps_bar = + merge + (to_left_factor (fun i -> position i < 0) (epsilon eps)) + (to_right_factor (fun i -> position i < 0) (epsilon_bar eps_bar)) + + let make_even tails tips = + List.rev_map + (List.rev_map2 (fun tail tip -> (tail, tip)) tails) + (Combinatorics.permute_even tips) + + let make_odd tails tips = + List.rev_map + (List.rev_map2 (fun tail tip -> (tail, tip)) tails) + (Combinatorics.permute_odd tips) + + let assert_eps_aux even odd eps eps_bar = + assert_equal ~printer:determinant_to_string + (canonicalize_determinant (make_determinant even odd)) + (canonicalize_determinant (merge_epsilon_pair eps eps_bar)) + + let assert_eps unit eps eps_bar = + let tips, tails = List.split unit in + let even = make_even tails tips + and odd = make_odd tails tips in + assert_eps_aux even odd eps eps_bar + + (* A single arrow needs special treatment to get the + sign to the proper place.*) + let assert_eps1 odd (tip, tail) eps eps_bar = + if odd then + assert_eps_aux [] [[(tail,tip)]] eps eps_bar + else + assert_eps_aux [[(tail,tip)]] [] eps eps_bar - [ "cycle []" >:: - (fun () -> - assert_equal [] (cycle [])); - - "cycle [1]" >:: - (fun () -> - assert_equal [1 => 1] (cycle [1])); - - "cycle [1;2]" >:: - (fun () -> - assert_equal [1 => 2; 2 => 1] (cycle [1; 2])); + let suite_fuse_epsilons = + "fuse_epsilons" >::: - "cycle [1;2;3]" >:: - (fun () -> - assert_equal [1 => 2; 2 => 3; 3 => 1] (cycle [1; 2; 3])); - - "cycle [1;2;3;4]" >:: - (fun () -> - assert_equal - [1 => 2; 2 => 3; 3 => 4; 4 => 1] - (cycle [1; 2; 3; 4])) ] + [ "1a*2a" >:: (fun () -> assert_eps1 false (1,2) [1;-9] [2;-9]); + "a1*a2" >:: (fun () -> assert_eps1 false (1,2) [-9;1] [-9;2]); + "1a*a2" >:: (fun () -> assert_eps1 true (1,2) [1;-9] [-9;2]); + "a1*2a" >:: (fun () -> assert_eps1 true (1,2) [-9;1] [2;-9]); + + "13a*24a" >:: (fun () -> assert_eps [(1,2);(3,4)] [1;3;-9] [2;4;-9]); + "1a3*24a" >:: (fun () -> assert_eps [(1,4);(3,2)] [1;-9;3] [2;4;-9]); + "a13*2a4" >:: (fun () -> assert_eps [(1,4);(3,2)] [-9;1;3] [2;-9;4]); + "1a3*2a4" >:: (fun () -> assert_eps [(1,2);(3,4)] [1;-9;3] [2;-9;4]); + + "135a*246a" >:: (fun () -> assert_eps [(1,2);(3,4);(5,6)] [1;3;5;-9] [2;4;6;-9]); + "315a*246a" >:: (fun () -> assert_eps [(3,2);(1,4);(5,6)] [3;1;5;-9] [2;4;6;-9]); + "5a13*246a" >:: (fun () -> assert_eps [(1,2);(3,4);(5,6)] [5;-9;1;3] [2;4;6;-9]); + "a135*2a46" >:: (fun () -> assert_eps [(3,2);(1,4);(5,6)] [-9;1;3;5] [2;-9;4;6]) ] let suite = "Color.Arrow" >::: [suite_chain; - suite_cycle] + suite_cycle; + suite_take; + suite_take2; + suite_replace; + suite_fuse_epsilons] let suite_long = "Color.Arrow long" >::: [] end let pp_free fmt f = Format.fprintf fmt "%s" (free_to_string f) let pp_factor fmt f = Format.fprintf fmt "%s" (factor_to_string f) end module type Propagator = sig type cf_in = int type cf_out = int type t = W | I of cf_in | O of cf_out | IO of cf_in * cf_out | G val to_string : t -> string end module Propagator : Propagator = struct type cf_in = int type cf_out = int type t = W | I of cf_in | O of cf_out | IO of cf_in * cf_out | G let to_string = function | W -> "W" | I cf -> Printf.sprintf "I(%d)" cf | O cf' -> Printf.sprintf "O(%d)" cf' | IO (cf, cf') -> Printf.sprintf "IO(%d,%d)" cf cf' | G -> "G" end module type LP = sig val rationals : (Algebra.Q.t * int) list -> Algebra.Laurent.t val ints : (int * int) list -> Algebra.Laurent.t val rational : Algebra.Q.t -> Algebra.Laurent.t val int : int -> Algebra.Laurent.t val fraction : int -> Algebra.Laurent.t val imag : int -> Algebra.Laurent.t val nc : int -> Algebra.Laurent.t val over_nc : int -> Algebra.Laurent.t end module LP : LP = struct module L = Algebra.Laurent (* Rationals from integers. *) let q_int n = Q.make n 1 let q_fraction n = Q.make 1 n (* Complex rationals: *) let qc_rational q = QC.make q Q.null let qc_int n = qc_rational (q_int n) let qc_fraction n = qc_rational (q_fraction n) let qc_imag n = QC.make Q.null (q_int n) (* Laurent polynomials: *) let of_pairs f pairs = L.sum (List.map (fun (coeff, power) -> L.atom (f coeff) power) pairs) let rationals = of_pairs qc_rational let ints = of_pairs qc_int let rational q = rationals [(q, 0)] let int n = ints [(n, 0)] let fraction n = L.const (qc_fraction n) let imag n = L.const (qc_imag n) let nc n = ints [(n, 1)] let over_nc n = ints [(n, -1)] end module type Birdtracks = sig type t val canonicalize : t -> t val to_string : t -> string val trivial : t -> bool val is_null : t -> bool val const : Algebra.Laurent.t -> t val null : t val one : t val two : t val half : t val third : t val minus : t val int : int -> t val fraction : int -> t val nc : t val over_nc : t val imag : t val ints : (int * int) list -> t val scale : QC.t -> t -> t val sum : t list -> t val diff : t -> t -> t val times : t -> t -> t val multiply : t list -> t module Infix : sig val ( +++ ) : t -> t -> t val ( --- ) : t -> t -> t val ( *** ) : t -> t -> t end val f_of_rep : (int -> int -> int -> t) -> int -> int -> int -> t val d_of_rep : (int -> int -> int -> t) -> int -> int -> int -> t val relocate : (int -> int) -> t -> t val fuse : int -> t -> Propagator.t list -> (QC.t * Propagator.t) list module Test : Test val pp : Format.formatter -> t -> unit end module Birdtracks = struct module A = Arrow open A.Infix module P = Propagator module L = Algebra.Laurent type connection = L.t * A.free list type t = connection list let trivial = function | [] -> true | [(coeff, [])] -> coeff = L.unit | _ -> false (* Rationals from integers. *) let q_int n = Q.make n 1 let q_fraction n = Q.make 1 n (* Complex rationals: *) let qc_rational q = QC.make q Q.null let qc_int n = qc_rational (q_int n) let qc_fraction n = qc_rational (q_fraction n) let qc_imag n = QC.make Q.null (q_int n) (* Laurent polynomials: *) let laurent_of_pairs f pairs = L.sum (List.map (fun (coeff, power) -> L.atom (f coeff) power) pairs) let l_rationals = laurent_of_pairs qc_rational let l_ints = laurent_of_pairs qc_int let l_rational q = l_rationals [(q, 0)] let l_int n = l_ints [(n, 0)] let l_fraction n = L.const (qc_fraction n) let l_imag n = L.const (qc_imag n) let l_nc n = l_ints [(n, 1)] let l_over_nc n = l_ints [(n, -1)] (* Expressions *) let const c = [c, []] let ints pairs = const (LP.ints pairs) let null = const L.null let half = const (LP.fraction 2) let third = const (LP.fraction 3) let fraction n = const (LP.fraction n) let one = const (LP.int 1) let two = const (LP.int 2) let minus = const (LP.int (-1)) let int n = const (LP.int n) let nc = const (LP.nc 1) let over_nc = const (LP.ints [(1, -1)]) let imag = const (LP.imag 1) module AMap = Pmap.Tree let find_arrows_opt arrows map = try Some (AMap.find pcompare arrows map) with Not_found -> None let canonicalize1 (coeff, io_list) = (coeff, List.sort pcompare io_list) let canonicalize terms = let map = List.fold_left (fun acc term -> let coeff, arrows = canonicalize1 term in if L.is_null coeff then acc else match find_arrows_opt arrows acc with | None -> AMap.add pcompare arrows coeff acc | Some coeff' -> let coeff'' = L.add coeff coeff' in if L.is_null coeff'' then AMap.remove pcompare arrows acc else AMap.add pcompare arrows coeff'' acc) AMap.empty terms in if AMap.is_empty map then null else AMap.fold (fun arrows coeff acc -> (coeff, arrows) :: acc) map [] let arrows_to_string_aux f arrows = ThoList.to_string f arrows let to_string1_aux f (coeff, arrows) = Printf.sprintf "(%s) * %s" (L.to_string "N" coeff) (arrows_to_string_aux f arrows) - let to_string1_opt_aux f = function - | None -> "None" - | Some v -> to_string1_aux f v - let to_string_raw_aux f v = ThoList.to_string (to_string1_aux f) v let to_string_aux f v = to_string_raw_aux f (canonicalize v) let factor_arrows_to_string = arrows_to_string_aux A.factor_to_string let factor_to_string1 = to_string1_aux A.factor_to_string - let factor_to_string1_opt = to_string1_opt_aux A.factor_to_string let factor_to_string_raw = to_string_raw_aux A.factor_to_string let factor_to_string = to_string_aux A.factor_to_string let arrows_to_string = arrows_to_string_aux A.free_to_string let to_string1 = to_string1_aux A.free_to_string - let to_string1_opt = to_string1_opt_aux A.free_to_string let to_string_raw = to_string_raw_aux A.free_to_string let to_string = to_string_aux A.free_to_string let pp fmt v = Format.fprintf fmt "%s" (to_string v) let is_null v = List.for_all (fun (c, _) -> L.is_null c) (canonicalize v) let is_white = function | P.W -> true | _ -> false let relocate1 f (c, v) = (c, List.map (A.map (A.relocate f)) v) let relocate f = List.map (relocate1 f) + (* Only for documentiation: a [term] is a list of arrows with + a coefficient. *) + type term = L.t * A.factor list + + (* Avoid the recursion, if there is no summation index in [arrow]. + If [arrow] loops back to itself, replace it by a factor of~$N_C$. *) + let rec add_arrow : A.factor -> term -> term list = + fun arrow (coeff, arrows) -> + if A.is_free arrow then + [(coeff, arrow :: arrows)] + else if A.is_tadpole arrow then + [(L.mul (LP.nc 1) coeff, arrows)] + else + add_arrow' coeff [] arrow arrows + (* Add one [arrow] to a list of arrows, updating [coeff] - if necessary. Accumulate already processed arrows in [acc]. - Returns [None] if there is a mismatch (a gluon meeting - a ghost), [Some (coeff', arrows')] otherwise. *) - let rec add_arrow' arrow (coeff, acc) = function + if necessary. Accumulate already processed arrows in [seen]. + Returns an empty list if there is a mismatch (a gluon meeting + a ghost) and a list of pairs consisting of a coefficient and a + list of arrows otherwise. There can be more than one pair, + because matching $\epsilon$ and $\bar\epsilon$ results + in a sum over permutations. *) + + and add_arrow' : L.t -> A.factor list -> A.factor -> A.factor list -> term list = + fun coeff seen arrow -> function | [] -> (* visited all [arrows]: no opportunities for further matches *) - Some (coeff, arrow :: acc) + [(coeff, arrow :: seen)] | arrow' :: arrows' -> begin match A.merge arrow arrow' with | A.Mismatch -> - None + [] | A.Ghost_Match -> (* replace matching ghosts by $-1/N_C$ *) - Some (L.mul (LP.over_nc (-1)) coeff, List.rev_append acc arrows') + [(L.mul (LP.over_nc (-1)) coeff, List.rev_append seen arrows')] | A.Loop_Match -> (* replace a loop by $N_C$ *) - Some (L.mul (LP.nc 1) coeff, List.rev_append acc arrows') + [(L.mul (LP.nc 1) coeff, List.rev_append seen arrows')] | A.Match arrow'' -> (* two arrows have been merged into one *) if A.is_free arrow'' then (* no opportunities for further matches *) - Some (coeff, arrow'' :: List.rev_append acc arrows') + [(coeff, arrow'' :: List.rev_append seen arrows')] else (* the new [arrow''] ist not yet saturated, try again: *) - add_arrow' arrow'' (coeff, acc) arrows' + add_arrow' coeff seen arrow'' arrows' + | A.Determinant (even, odd) -> + add_determinant seen even odd (coeff, arrows') | A.No_Match -> (* recurse to the remaining arrows *) - add_arrow' arrow (coeff, arrow' :: acc) arrows' + add_arrow' coeff (arrow' :: seen) arrow arrows' end - (* Avoid the recursion, if there is no summation index in [arrow]. *) - let add_arrow arrow (coeff, arrows) = - if A.is_free arrow then - Some (coeff, arrow :: arrows) - else - add_arrow' arrow (coeff, []) arrows + and add_determinant : A.factor list -> A.factor list list -> A.factor list list -> term -> term list = + fun seen even odd (coeff, arrows as term) -> + distribute seen even term (distribute seen odd (L.neg coeff, arrows) []) + + and distribute : A.factor list -> A.factor list list-> term -> term list -> term list = + fun seen permutations term terms -> + List.fold_left + (fun acc permutation -> splice_arrows seen permutation term :: acc) + terms permutations + + and splice_arrows : A.factor list -> A.factor list -> term -> term = + fun seen arrows term -> + let coeff', arrows' = add_arrow_list arrows term in + (coeff', List.rev_append seen arrows') + + (* \begin{dubious} + Here we would like to use the type system to prove + that the two failing cases can't happen. + In real life they can't happen, because [arrow] is + never [A.Epsilon] \ldots + \end{dubious} *) + and add_arrow_list : A.factor list -> term -> term = + fun arrows term -> + match arrows with + | [] -> term + | arrow :: rest -> + begin match add_arrow arrow term with + | [term] -> add_arrow_list rest term + | [] -> failwith "add_arrow_list: unexpected empty list" + | _ -> failwith "add_arrow_list: unexpected multi element list" + end + +(*i and add_arrow_list arrows (coeff, acc as term) = + let result = add_arrow_list_raw arrows term in + Printf.eprintf + "add_arrow_list (%s) * %s %s ==> %s\n" + (factor_arrows_to_string arrows) + (L.to_string "N" coeff) + (factor_arrows_to_string acc) + (factor_to_string1 result); + result +i*) let logging_add_arrow arrow (coeff, arrows) = let result = add_arrow arrow (coeff, arrows) in Printf.eprintf "add_arrow %s to %s ==> %s\n" (A.factor_to_string arrow) (factor_to_string1 (coeff, arrows)) - (factor_to_string1_opt result); + (ThoList.to_string factor_to_string1 result); result - (* We can reject the contributions with unsaturated summation indices +(*i + let add_arrow = logging_add_arrow +i*) + + (* [add_arrows term arrows] add the [arrows] to [term] by calling + [add_arrow] for each one. Return an empty list if there are + leftover summation indices in the end. *) + + (* NB: we can reject the contributions with unsaturated summation indices from Ghost contributions to~$T_a$ only \emph{after} adding all arrows that might saturate an open index. *) - let add_arrows factor1 arrows2 = - let rec add_arrows' (_, arrows as acc) = function - | [] -> - if List.for_all A.is_free arrows then - Some acc - else - None - | arrow :: arrows -> - begin match add_arrow arrow acc with - | None -> None - | Some acc' -> add_arrows' acc' arrows - end in - add_arrows' factor1 arrows2 + let rec add_arrows : term -> A.factor list -> term list = + fun (_, acc_arrows as acc) -> function + | [] -> + if List.for_all A.is_free acc_arrows then + [acc] + else + [] + | arrow :: arrows -> + ThoList.flatmap (fun term -> add_arrows term arrows) (add_arrow arrow acc) let logging_add_arrows factor1 arrows2 = let result = add_arrows factor1 arrows2 in Printf.eprintf "add_arrows %s to %s ==> %s\n" (factor_to_string1 factor1) (factor_arrows_to_string arrows2) - (factor_to_string1_opt result); + (ThoList.to_string factor_to_string1 result); result +(*i + let add_arrows = logging_add_arrows +i*) (* Note that a negative index might be summed only later in a sequence of binary products and must therefore be treated as free in this product. Therefore, we have to classify the indices as summation indices \emph{not only} based on their sign, but in addition based on whether they appear in both factors. Only then can we reject surviving ghosts. *) module ESet = Set.Make (struct type t = A.endpoint let compare = pcompare end) let negatives arrows = List.fold_left (fun acc arrow -> List.fold_left (fun acc' i -> ESet.add i acc') acc (A.negatives arrow)) ESet.empty arrows let times1 (coeff1, arrows1) (coeff2, arrows2) = let summations = ESet.inter (negatives arrows1) (negatives arrows2) in let is_sum i = ESet.mem i summations in let arrows1' = List.map (A.to_left_factor is_sum) arrows1 and arrows2' = List.map (A.to_right_factor is_sum) arrows2 in - match add_arrows (coeff1, arrows1') arrows2' with - | None -> None - | Some (coeff1, arrows) -> - Some (L.mul coeff1 coeff2, List.map A.of_factor arrows) + List.map + (fun (coeff1, arrows) -> + (L.mul coeff1 coeff2, List.map A.of_factor arrows)) + (add_arrows (coeff1, arrows1') arrows2') let logging_times1 factor1 factor2 = let result = times1 factor1 factor2 in Printf.eprintf "%s times1 %s ==> %s\n" (to_string1 factor1) (to_string1 factor2) - (to_string1_opt result); + (ThoList.to_string to_string1 result); result let sum terms = canonicalize (List.concat terms) let times term term' = - canonicalize (Product.list2_opt times1 term term') + canonicalize + (Product.fold2 + (fun x y -> List.rev_append (times1 x y)) + term term' []) (* \begin{dubious} Is that more efficient than the following implementation? \end{dubious} *) let rec multiply1' acc = function - | [] -> Some acc + | [] -> [acc] | factor :: factors -> - begin match times1 acc factor with - | None -> None - | Some acc' -> multiply1' acc' factors - end + List.fold_right multiply1' (times1 acc factor) factors let multiply1 = function - | [] -> Some (L.unit, []) - | [factor] -> Some factor + | [] -> [(L.unit, [])] + | [factor] -> [factor] | factor :: factors -> multiply1' factor factors - let multiply termss = - canonicalize (Product.list_opt multiply1 termss) + let multiply terms = + canonicalize + (Product.fold (fun x -> List.rev_append (multiply1 x)) terms []) (* \begin{dubious} Isn't that the more straightforward implementation? \end{dubious} *) let multiply = function | [] -> [] | term :: terms -> canonicalize (List.fold_left times term terms) let scale1 q (coeff, arrows) = (L.scale q coeff, arrows) let scale q = List.map (scale1 q) let diff term1 term2 = canonicalize (List.rev_append term1 (scale (qc_int (-1)) term2)) module Infix = struct let ( +++ ) term term' = sum [term; term'] let ( --- ) = diff let ( *** ) = times end open Infix (* Compute $ \tr(r(T_a) r(T_b) r(T_c)) $. NB: this uses the summation indices $-1$, $-2$ and $-3$. Therefore it \emph{must not} appear unevaluated more than once in a product! *) let trace3 r a b c = r a (-1) (-2) *** r b (-2) (-3) *** r c (-3) (-1) let f_of_rep r a b c = minus *** imag *** (trace3 r a b c --- trace3 r a c b) (* $ d_{abc} = \tr(r(T_a) [r(T_b), r(T_c)]_+) $ *) let d_of_rep r a b c = trace3 r a b c +++ trace3 r a c b (* \thocwmodulesubsection{Feynman Rules} *) module IMap = Map.Make (struct type t = int let compare = pcompare end) let line_map lines = let _, map = List.fold_left (fun (i, acc) line -> (succ i, match line with | P.W -> acc | _ -> IMap.add i line acc)) (1, IMap.empty) lines in map (*i Redundant since ocaml 4.05 let find_opt i map = try Some (IMap.find i map) with Not_found -> None i*) let lines_to_string lines = match IMap.bindings lines with | [] -> "W" | lines -> String.concat " " (List.map (fun (i, c) -> Printf.sprintf "%s@%d" (P.to_string c) i) lines) let clear = IMap.remove let add_in i cf lines = match IMap.find_opt i lines with | Some (P.O cf') -> IMap.add i (P.IO (cf, cf')) lines | _ -> IMap.add i (P.I cf) lines let add_out i cf' lines = match IMap.find_opt i lines with | Some (P.I cf) -> IMap.add i (P.IO (cf, cf')) lines | _ -> IMap.add i (P.O cf') lines let add_ghost i lines = IMap.add i P.G lines let connect1 n arrow lines = match arrow with | A.Ghost g -> let g = A.position g in if g = n then Some (add_ghost n lines) else begin match IMap.find_opt g lines with | Some P.G -> Some (clear g lines) | _ -> None end | A.Arrow (i, o) -> let i = A.position i and o = A.position o in if o = n then - match IMap.find_opt i lines with + begin match IMap.find_opt i lines with | Some (P.I cfi) -> Some (add_in o cfi (clear i lines)) | Some (P.IO (cfi, cfi')) -> Some (add_in o cfi (add_out i cfi' lines)) | _ -> None + end else if i = n then - match IMap.find_opt o lines with + begin match IMap.find_opt o lines with | Some (P.O cfo') -> Some (add_out i cfo' (clear o lines)) | Some (P.IO (cfo, cfo')) -> Some (add_out i cfo' (add_in o cfo lines)) | _ -> None + end else - match IMap.find_opt i lines, IMap.find_opt o lines with + begin match IMap.find_opt i lines, IMap.find_opt o lines with | Some (P.I cfi), Some (P.O cfo') when cfi = cfo' -> Some (clear o (clear i lines)) | Some (P.I cfi), Some (P.IO (cfo, cfo')) when cfi = cfo'-> Some (add_in o cfo (clear i lines)) | Some (P.IO (cfi, cfi')), Some (P.O cfo') when cfi = cfo' -> Some (add_out i cfi' (clear o lines)) | Some (P.IO (cfi, cfi')), Some (P.IO (cfo, cfo')) when cfi = cfo' -> Some (add_in o cfo (add_out i cfi' lines)) | _ -> None + end + | A.Epsilon _ -> + failwith "Birdtracks.connect not yet defined for Epsilon" + | A.Epsilon_bar _ -> + failwith "Birdtracks.connect not yet defined for Epsilon_bar" let connect connections lines = let n = succ (List.length lines) and lines = line_map lines in let rec connect' acc = function | arrow :: arrows -> begin match connect1 n arrow acc with | None -> None | Some acc -> connect' acc arrows end | [] -> Some acc in match connect' lines connections with | None -> None | Some acc -> begin match IMap.bindings acc with | [] -> Some P.W | [(i, cf)] when i = n -> Some cf | _ -> None end let fuse1 nc lines (c, vertex) = match connect vertex lines with | None -> [] | Some cf -> [(L.eval (qc_int nc) c, cf)] let fuse nc vertex lines = match vertex with | [] -> if List.for_all is_white lines then [(QC.unit, P.W)] else [] | vertex -> ThoList.flatmap (fuse1 nc lines) vertex module Test : Test = struct open OUnit let vertices_equal v1 v2 = - match v1, v2 with - | None, None -> true - | Some v1, Some v2 -> (canonicalize1 v1) = (canonicalize1 v2) - | _ -> false + (canonicalize v1) = (canonicalize v2) let eq v1 v2 = - assert_equal ~printer:to_string1_opt ~cmp:vertices_equal v1 v2 + assert_equal ~printer:(ThoList.to_string to_string1) ~cmp:vertices_equal v1 v2 let suite_times1 = "times1" >::: [ "merge two" >:: (fun () -> eq - (Some (L.unit, 1 ==> 2)) + [(L.unit, 1 ==> 2)] (times1 (L.unit, 1 ==> -1) (L.unit, -1 ==> 2))); "merge two exchanged" >:: (fun () -> eq - (Some (L.unit, 1 ==> 2)) + [(L.unit, 1 ==> 2)] (times1 (L.unit, -1 ==> 2) (L.unit, 1 ==> -1))); "ghost1" >:: (fun () -> eq - (Some (l_over_nc (-1), 1 ==> 2)) + [(l_over_nc (-1), 1 ==> 2)] (times1 (L.unit, [-1 => 2; ?? (-3)]) (L.unit, [ 1 => -1; ?? (-3)]))); "ghost2" >:: (fun () -> eq - None + [] (times1 (L.unit, [ 1 => -1; ?? (-3)]) (L.unit, [-1 => 2; -3 => -4; -4 => -3]))); "ghost2 exchanged" >:: (fun () -> eq - None + [] (times1 (L.unit, [-1 => 2; -3 => -4; -4 => -3]) (L.unit, [ 1 => -1; ?? (-3)]))) ] let suite_canonicalize = "canonicalize" >::: [ ] let line_option_to_string = function | None -> "no match" | Some line -> P.to_string line let test_connect_msg vertex formatter (expected, result) = Format.fprintf formatter "[%s]: expected %s, got %s" (arrows_to_string vertex) (line_option_to_string expected) (line_option_to_string result) let test_connect expected lines vertex = assert_equal ~printer:line_option_to_string expected (connect vertex lines) let test_connect_permutations expected lines vertex = List.iter (fun v -> assert_equal ~pp_diff:(test_connect_msg v) expected (connect v lines)) (Combinatorics.permute vertex) let suite_connect = "connect" >::: [ "delta" >:: (fun () -> test_connect_permutations (Some (P.I 1)) [ P.I 1; P.W ] ( 1 ==> 3 )); "f: 1->3->2->1" >:: (fun () -> test_connect_permutations (Some (P.IO (1, 3))) [P.IO (1, 2); P.IO (2, 3)] (A.cycle [1; 3; 2])); "f: 1->2->3->1" >:: (fun () -> test_connect_permutations (Some (P.IO (1, 2))) [P.IO (3, 2); P.IO (1, 3)] (A.cycle [1; 2; 3])) ] let suite = "Color.Birdtracks" >::: [suite_times1; suite_canonicalize; suite_connect] let suite_long = "Color.Birdtracks long" >::: [] end let vertices_equal v1 v2 = is_null (v1 --- v2) let assert_zero_vertex v = OUnit.assert_equal ~printer:to_string ~cmp:vertices_equal null v (* As an extra protection agains vacuous tests, we make sure that the LHS does not vanish. *) let eq v1 v2 = OUnit.assert_bool "LHS = 0" (not (is_null v1)); OUnit.assert_equal ~printer:to_string ~cmp:vertices_equal v1 v2 end (* \thocwmodulesection{$\mathrm{SU}(N_C)$} We're computing with a general $N_C$, but [epsilon] and [epsilonbar] make only sense for $N_C=3$. Also some of the terminology alludes to $N_C=3$: triplet, sextet, octet. *) (* Using the normalization~$\tr(T_{a}T_{b}) = \delta_{ab}$, we can check the selfconsistency of the completeness relation \begin{equation} T_{a}^{i_1j_1} T_{a}^{i_2j_2} = \left( \delta^{i_1j_2} \delta^{i_2j_1} - \frac{1}{N_C} \delta^{i_1j_1} \delta^{j_1j_2}\right) \end{equation} as \begin{multline} T_{a}^{i_1j_1} T_{a}^{i_2j_2} = \tr\left(T_{a_1}T_{a_2}\right) T_{a_1}^{i_1j_1} T_{a_2}^{i_2j_2} = T_{a_1}^{l_1l_2} T_{a_2}^{l_2l_1} T_{a_1}^{i_1j_1} T_{a_2}^{i_2j_2} \\ = \left( \delta^{l_1j_1} \delta^{i_1l_2} - \frac{1}{N_C} \delta^{l_1l_2} \delta^{i_1j_1}\right) \left( \delta^{l_2j_2} \delta^{i_2l_1} - \frac{1}{N_C} \delta^{l_2l_1} \delta^{i_2j_2}\right) = \left( \delta^{i_1j_2} \delta^{i_2j_1} - \frac{1}{N_C} \delta^{i_1i_2} \delta^{j_2j_1}\right) \end{multline} With \begin{equation} \label{eq:f=tr(TTT)'} \ii f_{a_1a_2a_3} = \tr\left(T_{a_1}\left\lbrack T_{a_2},T_{a_3}\right\rbrack\right) = \tr\left(T_{a_1}T_{a_2}T_{a_3}\right) - \tr\left(T_{a_1}T_{a_3}T_{a_2}\right) \end{equation} and \begin{multline} \tr\left(T_{a_1}T_{a_2}T_{a_3}\right) T_{a_1}^{i_1j_1} T_{a_2}^{i_2j_2} T_{a_3}^{i_3j_3} = T_{a_1}^{l_1l_2} T_{a_2}^{l_2l_3} T_{a_3}^{l_3l_1} T_{a_1}^{i_1j_1} T_{a_2}^{i_2j_2} T_{a_3}^{i_3j_3} = \\ \left( \delta^{l_1j_1} \delta^{i_1l_2} - \frac{1}{N_C} \delta^{l_1l_2} \delta^{i_1j_1}\right) \left( \delta^{l_2j_2} \delta^{i_2l_3} - \frac{1}{N_C} \delta^{l_2l_3} \delta^{i_2j_2}\right) \left( \delta^{l_3j_3} \delta^{i_3l_1} - \frac{1}{N_C} \delta^{l_3l_1} \delta^{i_3j_3}\right) \end{multline} we find the decomposition \begin{equation} \label{eq:fTTT'} \ii f_{a_1a_2a_3} T_{a_1}^{i_1j_1}T_{a_2}^{i_2j_2}T_{a_3}^{i_3j_3} = \delta^{i_1j_2}\delta^{i_2j_3}\delta^{i_3j_1} - \delta^{i_1j_3}\delta^{i_3j_2}\delta^{i_2j_1}\,. \end{equation} *) (* Indeed, \begin{verbatim} symbol nc; Dimension nc; vector i1, i2, i3, j1, j2, j3; index l1, l2, l3; local [TT] = ( j1(l1) * i1(l2) - d_(l1,l2) * i1.j1 / nc ) * ( j2(l2) * i2(l1) - d_(l2,l1) * i2.j2 / nc ); #procedure TTT(sign) local [TTT`sign'] = ( j1(l1) * i1(l2) - d_(l1,l2) * i1.j1 / nc ) * ( j2(l2) * i2(l3) - d_(l2,l3) * i2.j2 / nc ) * ( j3(l3) * i3(l1) - d_(l3,l1) * i3.j3 / nc ) `sign' ( j1(l1) * i1(l2) - d_(l1,l2) * i1.j1 / nc ) * ( j3(l2) * i3(l3) - d_(l2,l3) * i3.j3 / nc ) * ( j2(l3) * i2(l1) - d_(l3,l1) * i2.j2 / nc ); #endprocedure #call TTT(-) #call TTT(+) bracket nc; print; .sort .end \end{verbatim} gives \begin{verbatim} [TT] = + nc^-1 * ( - i1.j1*i2.j2 ) + i1.j2*i2.j1; [TTT-] = + i1.j2*i2.j3*i3.j1 - i1.j3*i2.j1*i3.j2; [TTT+] = + nc^-2 * ( 4*i1.j1*i2.j2*i3.j3 ) + nc^-1 * ( - 2*i1.j1*i2.j3*i3.j2 - 2*i1.j2*i2.j1*i3.j3 - 2*i1.j3*i2.j2*i3.j1 ) + i1.j2*i2.j3*i3.j1 + i1.j3*i2.j1*i3.j2; \end{verbatim} *) module type SU3 = sig include Birdtracks val delta3 : int -> int -> t val delta8 : int -> int -> t val delta8_loop : int -> int -> t val gluon : int -> int -> t val delta6 : int -> int -> t val delta10 : int -> int -> t val t : int -> int -> int -> t val f : int -> int -> int -> t val d : int -> int -> int -> t - val epsilon : int -> int -> int -> t - val epsilonbar : int -> int -> int -> t + val epsilon : int list -> t + val epsilon_bar : int list -> t val t8 : int -> int -> int -> t val t6 : int -> int -> int -> t val t10 : int -> int -> int -> t val k6 : int -> int -> int -> t val k6bar : int -> int -> int -> t val delta_of_tableau : int Young.tableau -> int -> int -> t val t_of_tableau : int Young.tableau -> int -> int -> int -> t end module SU3 : SU3 = struct module A = Arrow open Arrow.Infix module B = Birdtracks type t = B.t let canonicalize = B.canonicalize let to_string = B.to_string let pp = B.pp let trivial = B.trivial let is_null = B.is_null let null = B.null let const = B.const let one = B.one let two = B.two let int = B.int let half = B.half let third = B.third let fraction = B.fraction let nc = B.nc let over_nc = B.over_nc let minus = B.minus let imag = B.imag let ints = B.ints let sum = B.sum let diff = B.diff let scale = B.scale let times = B.times let multiply = B.multiply let relocate = B.relocate let fuse = B.fuse let f_of_rep = B.f_of_rep let d_of_rep = B.d_of_rep module Infix = B.Infix (* \thocwmodulesubsection{Fundamental and Adjoint Representation} *) let delta3 i j = [(LP.int 1, j ==> i)] let delta8 a b = [(LP.int 1, a <=> b)] (* If the~$\delta_{ab}$ originates from a~$\tr(T_aT_b)$, like an effective~$gg\to H$ coupling, it makes a difference in the color flow basis and we must write the full expression~(6.2) from~\cite{Kilian:2012pz} including the ghosts instead. Note that the sign for the terms with one ghost has not been spelled out in that reference. *) let delta8_loop a b = [(LP.int 1, a <=> b); (LP.int (-1), [a => a; ?? b]); (LP.int (-1), [?? a; b => b]); (LP.nc 1, [?? a; ?? b])] (* The following can be used for computing polarization sums (eventually, this could make the [Flow] module redundant). Note that we have $-N_C$ instead of $-1/N_C$ in the ghost contribution here, because two factors of $-1/N_C$ will be produced by [add_arrow] below, when contracting two ghost indices. Indeed, with this definition we can maintain [multiply [delta8 1 (-1); gluon (-1) (-2); delta8 (-2) 2] = delta8 1 2]. *) let ghost a b = [ (LP.nc (-1), [?? a; ?? b])] let gluon a b = delta8 a b @ ghost a b (* Note that the arrow is directed from the second to the first index, opposite to our color flow paper~\cite{Kilian:2012pz}. Fortunately, this is just a matter of conventions. \begin{subequations} \begin{align} \parbox{28\unitlength}{% \fmfframe(4,4)(4,4){% \begin{fmfgraph*}(20,20) \fmfleft{f1,f2} \fmfright{g} \fmfv{label=$i$}{f2} \fmfv{label=$j$}{f1} \fmfv{label=$a$}{g} \fmf{fermion}{f1,v} \fmf{fermion}{v,f2} \fmf{gluon}{v,g} \end{fmfgraph*}}} &\Longrightarrow \parbox{28\unitlength}{% \fmfframe(4,4)(4,4){% \begin{fmfgraph*}(20,20) \fmfleft{f1,f2} \fmfright{g} \fmfv{label=$i$}{f2} \fmfv{label=$j$}{f1} \fmfv{label=$a$}{g} \fmf{phantom}{f1,v} \fmf{phantom}{v,f2} \fmf{phantom}{v,g} \fmffreeze \fmfi{phantom_arrow}{vpath (__v, __g) sideways -thick} \fmfi{phantom_arrow}{(reverse vpath (__v, __g)) sideways -thick} \fmfi{phantom_arrow}{vpath (__f1, __v)} \fmfi{phantom_arrow}{vpath (__v, __f2)} \fmfi{plain}{% (vpath (__f1, __v) join (vpath (__v, __g)) sideways -thick)} \fmfi{plain}{% ((reverse vpath (__g, __v) sideways -thick) join vpath (__v, __f2))} \end{fmfgraph*}}} \parbox{28\unitlength}{% \fmfframe(4,4)(4,4){% \begin{fmfgraph*}(20,20) \fmfleft{f1,f2} \fmfright{g} \fmfv{label=$i$}{f1} \fmfv{label=$j$}{f2} \fmfv{label=$a$}{g} \fmf{fermion}{f1,v} \fmf{fermion}{v,f2} \fmf{dots}{v,g} \end{fmfgraph*}}}\\ T_a^{ij} \qquad\quad &\Longrightarrow \qquad\quad \delta^{ia}\delta^{aj} \qquad\qquad\qquad - \delta^{ij} \end{align} \end{subequations} *) let t a i j = [ (LP.int 1, [j => a; a => i]); (LP.int (-1), [j => i; ?? a]) ] (* Note that while we expect $\tr(T_a)=T_a^{ii}=0$, the evaluation of the expression [t 1 (-1) (-1)] will stop at [ [ -1 => 1; 1 => -1 ] --- [ -1 => -1; ?? 1 ] ], because the summation index appears in a single term. However, a naive further evaluation would get stuck at [ [ 1 => 1 ] --- nc *** [ ?? 1 ] ]. Fortunately, traces of single generators are never needed in our applications. We just have to resist the temptation to use them in unit tests. *) (* \begin{equation} \parbox{29\unitlength}{% \fmfframe(2,2)(2,2){% \begin{fmfgraph*}(25,25) \fmfleft{g1,g2} \fmfright{g3} \fmfv{label=$a$}{g1} \fmfv{label=$b$}{g2} \fmfv{label=$c$}{g3} \fmf{gluon}{g1,v} \fmf{gluon}{g2,v} \fmf{gluon}{g3,v} \end{fmfgraph*}}} \qquad\Longrightarrow \parbox{29\unitlength}{% \fmfframe(2,2)(2,2){% \begin{fmfgraph*}(25,25) \fmfleft{g1,g2} \fmfright{g3} \fmfv{label=$a$}{g1} \fmfv{label=$b$}{g2} \fmfv{label=$c$}{g3} \fmf{phantom}{g1,v} \fmf{phantom}{g2,v} \fmf{phantom}{g3,v} \fmffreeze \fmfi{plain}{(vpath(__g1,__v) join (reverse vpath(__g2,__v))) sideways thick} \fmfi{plain}{(vpath(__g2,__v) join (reverse vpath(__g3,__v))) sideways thick} \fmfi{plain}{(vpath(__g3,__v) join (reverse vpath(__g1,__v))) sideways thick} \fmfi{phantom_arrow}{vpath (__g1, __v) sideways thick} \fmfi{phantom_arrow}{vpath (__g2, __v) sideways thick} \fmfi{phantom_arrow}{vpath (__g3, __v) sideways thick} \fmfi{phantom_arrow}{(reverse vpath (__g1, __v)) sideways thick} \fmfi{phantom_arrow}{(reverse vpath (__g2, __v)) sideways thick} \fmfi{phantom_arrow}{(reverse vpath (__g3, __v)) sideways thick} \end{fmfgraph*}}} \qquad \parbox{29\unitlength}{% \fmfframe(2,2)(2,2){% \begin{fmfgraph*}(25,25) \fmfleft{g1,g2} \fmfright{g3} \fmfv{label=$a$}{g1} \fmfv{label=$b$}{g2} \fmfv{label=$c$}{g3} \fmf{phantom}{g1,v} \fmf{phantom}{g2,v} \fmf{phantom}{g3,v} \fmffreeze \fmfi{plain}{(vpath(__g1,__v) join (reverse vpath(__g3,__v))) sideways thick} \fmfi{plain}{(vpath(__g2,__v) join (reverse vpath(__g1,__v))) sideways thick} \fmfi{plain}{(vpath(__g3,__v) join (reverse vpath(__g2,__v))) sideways thick} \fmfi{phantom_arrow}{vpath (__g1, __v) sideways thick} \fmfi{phantom_arrow}{vpath (__g2, __v) sideways thick} \fmfi{phantom_arrow}{vpath (__g3, __v) sideways thick} \fmfi{phantom_arrow}{(reverse vpath (__g1, __v)) sideways thick} \fmfi{phantom_arrow}{(reverse vpath (__g2, __v)) sideways thick} \fmfi{phantom_arrow}{(reverse vpath (__g3, __v)) sideways thick} \end{fmfgraph*}}} \end{equation} *) let f a b c = [ (LP.imag ( 1), A.cycle [a; b; c]); (LP.imag (-1), A.cycle [a; c; b]) ] (* The generator in the adjoint representation $T_a^{bc}=-\ii f_{abc}$: *) let t8 a b c = Birdtracks.Infix.( minus *** imag *** f a b c ) (* This $d_{abc}$ is now compatible with~(6.11) in our color flow paper~\cite{Kilian:2012pz}. The signs had been wrong in earlier versions of the code to match the missing sign in the ghost contribution to the generator~$T_a^{ij}$ above. *) let d a b c = [ (LP.int 1, A.cycle [a; b; c]); (LP.int 1, A.cycle [a; c; b]); (LP.int (-2), (a <=> b) @ [?? c]); (LP.int (-2), (b <=> c) @ [?? a]); (LP.int (-2), (c <=> a) @ [?? b]); (LP.int 2, [a => a; ?? b; ?? c]); (LP.int 2, [?? a; b => b; ?? c]); (LP.int 2, [?? a; ?? b; c => c]); (LP.nc (-2), [?? a; ?? b; ?? c]) ] (* \thocwmodulesubsection{Decomposed Tensor Product Representations} *) let pass_through m n incoming outgoing = List.rev_map2 (fun i o -> (m, i) >=>> (n, o)) incoming outgoing let delta_of_permutations n permutations k l = let incoming = ThoList.range 0 (pred n) and normalization = List.length permutations in List.rev_map (fun (eps, outgoing) -> (LP.fraction (eps * normalization), pass_through l k incoming outgoing)) permutations let totally_symmetric n = List.map (fun p -> (1, p)) (Combinatorics.permute (ThoList.range 0 (pred n))) let totally_antisymmetric n = (Combinatorics.permute_signed (ThoList.range 0 (pred n))) let delta_S n k l = delta_of_permutations n (totally_symmetric n) k l let delta_A n k l = delta_of_permutations n (totally_antisymmetric n) k l let delta6 = delta_S 2 let delta10 = delta_S 3 let delta15 = delta_S 4 let delta3bar = delta_A 2 (* Mixed symmetries, as in section 9.4 of the birdtracks book. *) module IM = Partial.Make (struct type t = int let compare = pcompare end) module P = Permutation.Default (* Map the elements of [original] to [permuted] in [all], with [all] a list of $n$ integers from $0$ to $n-1$ in order, and use the resulting list to define a permutation. E.\,g.~[permute_partial [1;3] [3;1] [0;1;2;3;4]] will define a permutation that transposes the second and fourth element in a 5 element list. *) let permute_partial original permuted all = P.of_list (List.map (IM.auto (IM.of_lists original permuted)) all) let apply1 (sign, indices) (eps, p) = (eps * sign, P.list p indices) let apply signed_permutations signed_indices = List.rev_map (apply1 signed_indices) signed_permutations let apply_list signed_permutations signed_indices = ThoList.flatmap (apply signed_permutations) signed_indices let symmetrizer_of_permutations n original signed_permutations = let incoming = ThoList.range 0 (pred n) in List.rev_map (fun (eps, permuted) -> (eps, permute_partial original permuted incoming)) signed_permutations let symmetrizer n indices = symmetrizer_of_permutations n indices (List.rev_map (fun p -> (1, p)) (Combinatorics.permute indices)) let anti_symmetrizer n indices = symmetrizer_of_permutations n indices (Combinatorics.permute_signed indices) let symmetrize n elements indices = apply_list (symmetrizer n elements) indices let anti_symmetrize n elements indices = apply_list (anti_symmetrizer n elements) indices let id n = [(1, ThoList.range 0 (pred n))] (* \begin{dubious} We can avoid the recursion here, if we use [Combinatorics.permute_tensor_signed] in [symmetrizer] above. \end{dubious} *) let rec apply_tableau f n tableau indices = match tableau with | [] | [_] :: _ -> indices | cells :: rest -> apply_tableau f n rest (f n cells indices) (* \begin{dubious} Here we should at a sanity test for [tableau]: all integers should be consecutive starting from 0 with no duplicates. In additions the rows must not grow in length. \end{dubious} *) let delta_of_tableau tableau i j = let n = Young.num_cells_tableau tableau and num, den = Young.normalization (Young.diagram_of_tableau tableau) and rows = tableau - and cols = Young.transpose_tableau tableau in + and cols = Young.conjugate_tableau tableau in let permutations = apply_tableau symmetrize n rows (apply_tableau anti_symmetrize n cols (id n)) in Birdtracks.Infix.( int num *** fraction den *** delta_of_permutations n permutations i j ) let incomplete tensor = failwith ("Color.Vertex: " ^ tensor ^ " not supported yet!") let experimental tensor = Printf.eprintf "Color.Vertex: %s support still experimental and untested!\n" tensor -(* \begin{dubious} - Can we avoid nonlocality of the $\epsilon_{ijk}$ reduction, - as described in the revision of our color flow paper, - by simply using $\bar N\otimes_A \bar N$ instead of~$N$ on one - of the lines? - - This should work trivially, if we could always pick one flavor - appearing in the $\epsilon_{ijk}$ for this conversion, but this - is not guaranteed. - - As a hack, we could choose the color triplet bosons for - the $\bar N\otimes_A \bar N$ treatment, - as long as we can expect only $\epsilon_{ijk} \psi_i\psi_j\phi_k$ - couplings. This would take care of the RPV MSSM. - \end{dubious} *) + let distinct integers = + let rec distinct' seen = function + | [] -> true + | i :: rest -> + if Sets.Int.mem i seen then + false + else + distinct' (Sets.Int.add i seen) rest in + distinct' Sets.Int.empty integers + + (* All lines start here: they point towards the vertex. *) + let epsilon tips = + if distinct tips then + [(LP.int 1, [Arrow.epsilon tips])] + else + null (* All lines end here: they point away from the vertex. *) - let epsilon i j k = incomplete "epsilon-tensor" + let epsilon_bar tails = + if distinct tails then + [(LP.int 1, [Arrow.epsilon_bar tails])] + else + null - (* All lines start here: they point towards the vertex. *) - let epsilonbar i j k = incomplete "epsilon-tensor" (* In order to get the correct $N_C$ dependence of quadratic Casimir operators, the arrows in the vertex must have the same permutation symmetry as the propagator. This is demonstrated by the unit tests involving Casimir operators on page \pageref{pg:casimir-tests} below. These tests also provide a check of our normalization. The implementation takes a propagator and uses [Arrow.tee] to - replace one arrow by the pair of arrows correspondig to the + replace one arrow by the pair of arrows corresponding to the insertion of a gluon. This is repeated for each arrow. The normalization remains unchanged from the propagator. A minus sign is added for antiparallel arrows, since the conjugate representation is~$-T^*_a$. To this, we add the diagrams with a gluon connected to one arrow. Since these are identical, only one diagram multiplied by the difference of the number of parallel and antiparallel arrows is added. *) let insert_gluon a k l (norm, arrows) = let rec insert_gluon' acc left = function | [] -> acc | arrow :: right -> insert_gluon' ((Algebra.Laurent.mul (LP.int (A.dir k l arrow)) norm, List.rev_append left ((A.tee a arrow) @ right)) :: acc) (arrow :: left) right in insert_gluon' [] [] arrows let t_of_delta delta a k l = match delta k l with | [] -> [] | (_, arrows) :: _ as delta_kl -> let n = List.fold_left (fun acc arrow -> acc + A.dir k l arrow) 0 arrows in let ghosts = List.rev_map (fun (norm, arrows) -> (Algebra.Laurent.mul (LP.int (-n)) norm, ?? a :: arrows)) delta_kl in List.fold_left (fun acc arrows -> insert_gluon a k l arrows @ acc) ghosts delta_kl let t_of_delta delta a k l = canonicalize (t_of_delta delta a k l) let t_S n a k l = t_of_delta (delta_S n) a k l let t_A n a k l = t_of_delta (delta_A n) a k l let t6 = t_S 2 let t10 = t_S 3 let t15 = t_S 4 let t3bar = t_A 2 (* Equivalent definition: *) let t8' a b c = t_of_delta delta8 a b c let t_of_tableau tableau a k l = t_of_delta (delta_of_tableau tableau) a k l (* \begin{dubious} Check the following for a real live UFO file! \end{dubious} *) (* In the UFO paper, the Clebsh-Gordan is defined as~$K^{(6),ij}_{\hphantom{(6),ij}m}$. Therefore, keeping our convention for the generators~$T_{a\hphantom{(6),j}i}^{(6),j}$, the must arrows \emph{end} at~$m$. *) let k6 m i j = experimental "k6"; [ (LP.int 1, [i =>> (m, 0); j =>> (m, 1)]); (LP.int 1, [i =>> (m, 1); j =>> (m, 0)]) ] (* The arrow are reversed for~$\bar K^{(6),m}_{\hphantom{(6),m}ij}$ and \emph{start} at~$m$. *) let k6bar m i j = experimental "k6bar"; [ (LP.int 1, [(m, 0) >=> i; (m, 1) >=> j]); (LP.int 1, [(m, 1) >=> i; (m, 0) >=> j]) ] (* \thocwmodulesubsection{Unit Tests} *) module Test : Test = struct open OUnit module L = Algebra.Laurent module B = Birdtracks open Birdtracks open Birdtracks.Infix let exorcise vertex = List.filter (fun (_, arrows) -> not (List.exists A.is_ghost arrows)) vertex let eqx v1 v2 = eq (exorcise v1) (exorcise v2) (* \thocwmodulesubsection{Trivia} *) let suite_sum = "sum" >::: [ "atoms" >:: (fun () -> eq (two *** delta3 1 2) (delta3 1 2 +++ delta3 1 2)) ] let suite_diff = "diff" >::: [ "atoms" >:: (fun () -> eq (delta3 3 4) (delta3 1 2 +++ delta3 3 4 --- delta3 1 2)) ] let suite_times = "times" >::: [ "reorder components t1*t2" >:: (* trivial $T_a^{ik}T_a^{kj}=T_a^{kj}T_a^{ik}$ *) (fun () -> let t1 = t (-1) 1 (-2) and t2 = t (-1) (-2) 2 in eq (t1 *** t2) (t2 *** t1)); "reorder components tr(t1*t2)" >:: (* trivial $T_a^{ij}T_a^{ji}=T_a^{ji}T_a^{ij}$ *) (fun () -> let t1 = t 1 (-1) (-2) and t2 = t 2 (-2) (-1) in eq (t1 *** t2) (t2 *** t1)); "reorderings" >:: (fun () -> let v1 = [(L.unit, [ 1 => -2; -2 => -1; -1 => 1])] and v2 = [(L.unit, [-1 => 2; 2 => -2; -2 => -1])] and v' = [(L.unit, [ 1 => 1; 2 => 2])] in - eq v' (v1 *** v2)) ] + eq v' (v1 *** v2)); + + "eps*epsbar" >:: + (fun () -> + eq + (delta3 1 2 *** delta3 3 4 --- delta3 1 4 *** delta3 3 2) + (epsilon [-1; 1; 3] *** epsilon_bar [-1; 2; 4])); + + "eps*epsbar cyclic 1" >:: + (fun () -> + eq + (delta3 1 2 *** delta3 3 4 --- delta3 1 4 *** delta3 3 2) + (epsilon [3; -1; 1] *** epsilon_bar [-1; 2; 4])); + + "eps*epsbar cyclic 2" >:: + (fun () -> + eq + (delta3 1 2 *** delta3 3 4 --- delta3 1 4 *** delta3 3 2) + (epsilon [-1; 1; 3] *** epsilon_bar [4; -1; 2])); + + "eps*epsbar 2" >:: + (fun () -> + eq + (const (LP.ints [ (1, 1); (-1,0) ]) *** delta3 1 2) + (epsilon [-1; -2; 1] *** epsilon_bar [-1; -2; 2])); + + "eps*epsbar 3" >:: + (fun () -> + eq + (const (LP.ints [ (1, 2); (-1,1) ])) + (epsilon [-1; -2; -3] *** epsilon_bar [-1; -2; -3])) ] (* \thocwmodulesubsection{Propagators} *) (* Verify the normalization of the propagators by making sure that $D^{ij}D^{jk}=D^{ik}$ *) let projection_id rep_d = eq (rep_d 1 2) (rep_d 1 (-1) *** rep_d (-1) 2) let orthogonality d d' = assert_zero_vertex (d 1 (-1) *** d' (-1) 2) (* Pass every arrow straight through, without (anti-)symmetrization. *) let delta_unsymmetrized n k l = delta_of_permutations n [(1, ThoList.range 0 (pred n))] k l let completeness n tableaux = eq (delta_unsymmetrized n 1 2) (sum (List.map (fun t -> delta_of_tableau t 1 2) tableaux)) (* The following names are of historical origin. From the time, when we didn't have full support for Young tableaux and implemented figure 9.1 from the birdtrack book. \ytableausetup{centertableaux,smalltableaux} \begin{equation} \ytableaushort{01,2} \end{equation} *) let delta_SAS i j = delta_of_tableau [[0;1];[2]] i j (* \begin{equation} \ytableaushort{02,1} \end{equation} *) let delta_ASA i j = delta_of_tableau [[0;2];[1]] i j let suite_propagators = "propagators" >::: [ "D*D=D" >:: (fun () -> projection_id delta3); "D8*D8=D8" >:: (fun () -> projection_id delta8); "G*G=G" >:: (fun () -> projection_id gluon); "D6*D6=D6" >:: (fun () -> projection_id delta6); "D10*D10=D10" >:: (fun () -> projection_id delta10); "D15*D15=D15" >:: (fun () -> projection_id delta15); "D3bar*D3bar=D3bar" >:: (fun () -> projection_id delta3bar); "D6*D3bar=0" >:: (fun () -> orthogonality delta6 delta3bar); "D_A3*D_A3=D_A3" >:: (fun () -> projection_id (delta_A 3)); "D10*D_A3=0" >:: (fun () -> orthogonality delta10 (delta_A 3)); "D_SAS*D_SAS=D_SAS" >:: (fun () -> projection_id delta_SAS); "D_ASA*D_ASA=D_ASA" >:: (fun () -> projection_id delta_ASA); "D_SAS*D_S3=0" >:: (fun () -> orthogonality delta_SAS (delta_S 3)); "D_SAS*D_A3=0" >:: (fun () -> orthogonality delta_SAS (delta_A 3)); "D_SAS*D_ASA=0" >:: (fun () -> orthogonality delta_SAS delta_ASA); "D_ASA*D_SAS=0" >:: (fun () -> orthogonality delta_ASA delta_SAS); "D_ASA*D_S3=0" >:: (fun () -> orthogonality delta_ASA (delta_S 3)); "D_ASA*D_A3=0" >:: (fun () -> orthogonality delta_ASA (delta_A 3)); "DU*DU=DU" >:: (fun () -> projection_id (delta_unsymmetrized 3)); "S3=[0123]" >:: (fun () -> eq (delta_S 4 1 2) (delta_of_tableau [[0;1;2;3]] 1 2)); "A3=[0,1,2,3]" >:: (fun () -> eq (delta_A 4 1 2) (delta_of_tableau [[0];[1];[2];[3]] 1 2)); "[0123]*[012,3]=0" >:: (fun () -> orthogonality (delta_of_tableau [[0;1;2;3]]) (delta_of_tableau [[0;1;2];[3]])); "[0123]*[01,23]=0" >:: (fun () -> orthogonality (delta_of_tableau [[0;1;2;3]]) (delta_of_tableau [[0;1];[2;3]])); "[012,3]*[012,3]=[012,3]" >:: (fun () -> projection_id (delta_of_tableau [[0;1;2];[3]])); (* \ytableausetup{centertableaux,smalltableaux} \begin{equation} \ytableaushort{01} + \ytableaushort{0,1} \end{equation} *) "completeness 2" >:: (fun () -> completeness 2 [ [[0;1]]; [[0];[1]] ]) ; "completeness 2'" >:: (fun () -> eq (delta_unsymmetrized 2 1 2) (delta_S 2 1 2 +++ delta_A 2 1 2)); (* The normalization factors are written for illustration. They are added by [delta_of_tableau] automatically. \ytableausetup{centertableaux,smalltableaux} \begin{equation} \ytableaushort{012} + \frac{4}{3}\cdot\ytableaushort{01,2} + \frac{4}{3}\cdot\ytableaushort{02,1} + \ytableaushort{0,1,2} \end{equation} *) "completeness 3" >:: (fun () -> completeness 3 [ [[0;1;2]]; [[0;1];[2]]; [[0;2];[1]]; [[0];[1];[2]] ]); "completeness 3'" >:: (fun () -> eq (delta_unsymmetrized 3 1 2) (delta_S 3 1 2 +++ delta_SAS 1 2 +++ delta_ASA 1 2 +++ delta_A 3 1 2)); (* \ytableausetup{centertableaux,smalltableaux} \begin{equation} \ytableaushort{0123} + \frac{3}{2}\cdot\ytableaushort{012,3} + \frac{3}{2}\cdot\ytableaushort{013,2} + \frac{3}{2}\cdot\ytableaushort{023,1} + \frac{4}{3}\cdot\ytableaushort{01,23} + \frac{4}{3}\cdot\ytableaushort{02,13} + \frac{3}{2}\cdot\ytableaushort{01,2,3} + \frac{3}{2}\cdot\ytableaushort{02,1,3} + \frac{3}{2}\cdot\ytableaushort{03,1,2} + \ytableaushort{0,1,2,3} \end{equation} *) "completeness 4" >:: (fun () -> completeness 4 [ [[0;1;2;3]]; [[0;1;2];[3]]; [[0;1;3];[2]]; [[0;2;3];[1]]; [[0;1];[2;3]]; [[0;2];[1;3]]; [[0;1];[2];[3]]; [[0;2];[1];[3]]; [[0;3];[1];[2]]; [[0];[1];[2];[3]] ]) ] (* \thocwmodulesubsection{Normalization} *) let suite_normalization = "normalization" >::: [ "tr(t*t)" >:: (* $\tr(T_aT_b)=\delta_{ab} + \text{ghosts}$ *) (fun () -> eq (delta8_loop 1 2) (t 1 (-1) (-2) *** t 2 (-2) (-1))); "tr(t*t) sans ghosts" >:: (* $\tr(T_aT_b)=\delta_{ab}$ *) (fun () -> eqx (delta8 1 2) (t 1 (-1) (-2) *** t 2 (-2) (-1))); (* The additional ghostly terms were unexpected, but arises like~(6.2) in our color flow paper~\cite{Kilian:2012pz}. *) "t*t*t" >:: (* $T_aT_bT_a=-T_b/N_C + \ldots$ *) (fun () -> eq (minus *** over_nc *** t 1 2 3 +++ [(LP.int 1, [1 => 1; 3 => 2]); (LP.nc (-1), [3 => 2; ?? 1])]) (t (-1) 2 (-2) *** t 1 (-2) (-3) *** t (-1) (-3) 3)); (* As expected, these ghostly terms cancel in the summed squares \begin{equation} \tr(T_aT_bT_aT_cT_bT_c) = \tr(T_bT_b)/N_C^2 = \delta_{bb}/N_C^2 = (N_C^2-1) / N_C^2 = 1 - 1 / N_C^2 \end{equation} *) "sum((t*t*t)^2)" >:: (fun () -> eq (ints [(1, 0); (-1, -2)]) (t (-1) (-11) (-12) *** t (-2) (-12) (-13) *** t (-1) (-13) (-14) *** t (-3) (-14) (-15) *** t (-2) (-15) (-16) *** t (-3) (-16) (-11))); "d*d" >:: (fun () -> eqx [ (LP.ints [(2, 1); (-8,-1)], 1 <=> 2); (LP.ints [(2, 0); ( 4,-2)], [1=>1; 2=>2]) ] (d 1 (-1) (-2) *** d 2 (-2) (-1))) ] (* As proposed in our color flow paper~\cite{Kilian:2012pz}, we can get the correct (anti-)symmetrized generators by sandwiching the following unsymmetrized generators between the corresponding (anti-)symmetrized projectors. Therefore, the unsymmetrized generators work as long as they're used in Feynman diagrams, where they are connected by propagators that contain (anti-)symmetrized projectors. They even work in the Lie algebra relations and give the correct normalization there. They fail however for more general color algebra expressions that can appear in UFO files. In particular, the Casimir operators come out really wrong. *) let t_unsymmetrized n k l = t_of_delta (delta_unsymmetrized n) k l (* The following trivial vertices are \emph{not} used anymore, since they don't get the normalization of the Ward identities right. For the quadratic casimir operators, they always produce a result proportional to~$C_F=C_2(S_1)$. This can be understood because they correspond to a fundamental representation with spectators. (Anti-)symmetrizing by sandwiching with projectors almost works, but they must be multiplied by hand by the number of arrows to get the normalization right. They're here just for documenting what doesn't work. *) let t_trivial n a k l = let sterile = List.map (fun i -> (l, i) >=>> (k, i)) (ThoList.range 1 (pred n)) in [ (LP.int ( 1), ((l, 0) >=> a) :: (a =>> (k, 0)) :: sterile); (LP.int (-1), (?? a) :: ((l, 0) >=>> (k, 0)) :: sterile) ] let t6_trivial = t_trivial 2 let t10_trivial = t_trivial 3 let t15_trivial = t_trivial 4 let t_SAS = t_of_delta delta_SAS let t_ASA = t_of_delta delta_ASA let symmetrization ?rep_ts rep_tu rep_d = let rep_ts = match rep_ts with | None -> rep_tu | Some rep_t -> rep_t in eq (rep_ts 1 2 3) (gluon 1 (-1) *** rep_d 2 (-2) *** rep_tu (-1) (-2) (-3) *** rep_d (-3) 3) let suite_symmetrization = "symmetrization" >::: [ "t6" >:: (fun () -> symmetrization t6 delta6); "t10" >:: (fun () -> symmetrization t10 delta10); "t15" >:: (fun () -> symmetrization t15 delta15); "t3bar" >:: (fun () -> symmetrization t3bar delta3bar); "t_SAS" >:: (fun () -> symmetrization t_SAS delta_SAS); "t_ASA" >:: (fun () -> symmetrization t_ASA delta_ASA); "t6'" >:: (fun () -> symmetrization ~rep_ts:t6 (t_unsymmetrized 2) delta6); "t10'" >:: (fun () -> symmetrization ~rep_ts:t10 (t_unsymmetrized 3) delta10); "t15'" >:: (fun () -> symmetrization ~rep_ts:t15 (t_unsymmetrized 4) delta15); "t6''" >:: (fun () -> eq (t6 1 2 3) (int 2 *** delta6 2 (-1) *** t6_trivial 1 (-1) (-2) *** delta6 (-2) 3)); "t10''" >:: (fun () -> eq (t10 1 2 3) (int 3 *** delta10 2 (-1) *** t10_trivial 1 (-1) (-2) *** delta10 (-2) 3)); "t15''" >:: (fun () -> eq (t15 1 2 3) (int 4 *** delta15 2 (-1) *** t15_trivial 1 (-1) (-2) *** delta15 (-2) 3)) ] (* \thocwmodulesubsection{Traces} *) (* Compute (anti-)commutators of generators in the representation~$r$, i.\,e.~$[r(t_a)r(t_b)]_{ij}\mp[r(t_b)r(t_a)]_{ij}$, using [isum<0] as summation index in the matrix products. *) let commutator rep_t i_sum a b i j = multiply [rep_t a i i_sum; rep_t b i_sum j] --- multiply [rep_t b i i_sum; rep_t a i_sum j] let anti_commutator rep_t i_sum a b i j = multiply [rep_t a i i_sum; rep_t b i_sum j] +++ multiply [rep_t b i i_sum; rep_t a i_sum j] (* Trace of the product of three generators in the representation~$r$, i.\,e.~$\tr_r(r(t_a)r(t_b)r(t_c))$, using $-1,-2,-3$ as summation indices in the matrix products. *) let trace3 rep_t a b c = rep_t a (-1) (-2) *** rep_t b (-2) (-3) *** rep_t c (-3) (-1) let loop3 a b c = [ (LP.int 1, A.cycle (List.rev [a; b; c])); (LP.int (-1), (a <=> b) @ [?? c]); (LP.int (-1), (b <=> c) @ [?? a]); (LP.int (-1), (c <=> a) @ [?? b]); (LP.int 1, [a => a; ?? b; ?? c]); (LP.int 1, [?? a; b => b; ?? c]); (LP.int 1, [?? a; ?? b; c => c]); (LP.nc (-1), [?? a; ?? b; ?? c]) ] let suite_trace = "trace" >::: [ "tr(ttt)" >:: (fun () -> eq (trace3 t 1 2 3) (loop3 1 2 3)); "tr(ttt) cyclic 1" >:: (* $\tr(T_aT_bT_c)=\tr(T_bT_cT_a)$ *) (fun () -> eq (trace3 t 1 2 3) (trace3 t 2 3 1)); "tr(ttt) cyclic 2" >:: (* $\tr(T_aT_bT_c)=\tr(T_cT_aT_b)$ *) (fun () -> eq (trace3 t 1 2 3) (trace3 t 3 1 2)); (* \begin{dubious} Do we expect this? \end{dubious} *) "tr(tttt)" >:: (* $\tr(T_aT_bT_cT_d)=\ldots$ *) (fun () -> eqx [(LP.int 1, A.cycle [4; 3; 2; 1])] (t 1 (-1) (-2) *** t 2 (-2) (-3) *** t 3 (-3) (-4) *** t 4 (-4) (-1))) ] let suite_ghosts = "ghosts" >::: [ "H->gg" >:: (fun () -> eq (delta8_loop 1 2) (t 1 (-1) (-2) *** t 2 (-2) (-1))); "H->ggg f" >:: (fun () -> eq (imag *** f 1 2 3) (trace3 t 1 2 3 --- trace3 t 1 3 2)); "H->ggg d" >:: (fun () -> eq (d 1 2 3) (trace3 t 1 2 3 +++ trace3 t 1 3 2)); "H->ggg f'" >:: (fun () -> eq (imag *** f 1 2 3) (t 1 (-3) (-2) *** commutator t (-1) 2 3 (-2) (-3))); "H->ggg d'" >:: (fun () -> eq (d 1 2 3) (t 1 (-3) (-2) *** anti_commutator t (-1) 2 3 (-2) (-3))); "H->ggg cyclic'" >:: (fun () -> let trace a b c = t a (-3) (-2) *** commutator t (-1) b c (-2) (-3) in eq (trace 1 2 3) (trace 2 3 1)) ] let ff a1 a2 a3 a4 = [ (LP.int (-1), A.cycle [a1; a2; a3; a4]); (LP.int ( 1), A.cycle [a2; a1; a3; a4]); (LP.int ( 1), A.cycle [a1; a2; a4; a3]); (LP.int (-1), A.cycle [a2; a1; a4; a3]) ] let tf j i a b = [ (LP.imag ( 1), A.chain [i; a; b; j]); (LP.imag (-1), A.chain [i; b; a; j]) ] let suite_ff = "f*f" >::: [ "1" >:: (fun () -> eq (ff 1 2 3 4) (f (-1) 1 2 *** f (-1) 3 4)); "2" >:: (fun () -> eq (ff 1 2 3 4) (f (-1) 1 2 *** f 3 4 (-1))); "3" >:: (fun () -> eq (ff 1 2 3 4) (f (-1) 1 2 *** f 4 (-1) 3)) ] let suite_tf = "t*f" >::: [ "1" >:: (fun () -> eq (tf 1 2 3 4) (t (-1) 1 2 *** f (-1) 3 4)) ] (* \thocwmodulesubsection{Completeness Relation} *) (* Check the completeness relation corresponding to $q\bar q$-scattering: \begin{equation} \parbox{38\unitlength}{% \fmfframe(4,2)(4,4){% \begin{fmfgraph*}(30,20) \setupFourAmp \fmflabel{$i$}{i2} \fmflabel{$j$}{i1} \fmflabel{$k$}{o1} \fmflabel{$l$}{o2} \fmf{fermion}{i1,v1,i2} \fmf{fermion}{o2,v2,o1} \fmf{gluon}{v1,v2} \end{fmfgraph*}}} = \parbox{38\unitlength}{% \fmfframe(4,2)(4,4){% \begin{fmfgraph*}(30,20) \setupFourAmp \fmflabel{$i$}{i2} \fmflabel{$j$}{i1} \fmflabel{$k$}{o1} \fmflabel{$l$}{o2} \fmfi{phantom_arrow}{vpath (__i1, __v1)} \fmfi{phantom_arrow}{vpath (__v1, __v2) sideways -thick} \fmfi{phantom_arrow}{vpath (__v2, __o1)} \fmfi{phantom_arrow}{vpath (__o2, __v2)} \fmfi{phantom_arrow}{reverse vpath (__v1, __v2) sideways -thick} \fmfi{phantom_arrow}{vpath (__v1, __i2)} \fmfi{plain}{vpath (__i1, __v1) join (vpath (__v1, __v2) sideways -thick) join vpath (__v2, __o1)} \fmfi{plain}{vpath (__o2, __v2) join (reverse vpath (__v1, __v2) sideways -thick) join vpath (__v1, __i2)} \end{fmfgraph*}}} + \parbox{38\unitlength}{% \fmfframe(4,2)(4,4){% \begin{fmfgraph*}(30,20) \setupFourAmp \fmflabel{$i$}{i2} \fmflabel{$j$}{i1} \fmflabel{$k$}{o1} \fmflabel{$l$}{o2} \fmfi{phantom_arrow}{vpath (__i1, __v1)} \fmfi{phantom_arrow}{vpath (__v2, __o1)} \fmfi{phantom_arrow}{vpath (__o2, __v2)} \fmfi{phantom_arrow}{vpath (__v1, __i2)} \fmfi{plain}{vpath (__i1, __v1) join vpath (__v1, __i2)} \fmfi{plain}{vpath (__o2, __v2) join vpath (__v2, __o1)} \fmfi{dots,label=$-1/N_C$}{vpath (__v1, __v2)} \end{fmfgraph*}}} \end{equation} *) (* $T_{a}^{ij} T_{a}^{kl}$ *) let tt i j k l = t (-1) i j *** t (-1) k l (* $ \delta^{il}\delta^{kj} - \delta^{ij}\delta^{kl}/N_C$ *) let tt_expected i j k l = [ (LP.int 1, [l => i; j => k]); (LP.over_nc (-1), [j => i; l => k]) ] let suite_tt = "t*t" >::: [ "1" >:: (* $T_{a}^{ij} T_{a}^{kl} = \delta^{il}\delta^{kj} - \delta^{ij}\delta^{kl}/N_C$ *) (fun () -> eq (tt_expected 1 2 3 4) (tt 1 2 3 4)) ] (* \thocwmodulesubsection{Lie Algebra} *) (* Check the commutation relations $[T_a,T_b]=\ii f_{abc} T_c$ in various representations. *) let lie_algebra_id rep_t = let lhs = imag *** f 1 2 (-1) *** t (-1) 3 4 and rhs = commutator t (-1) 1 2 3 4 in eq lhs rhs (* Check the normalization of the structure consistants $\mathcal{N} f_{abc} = - \ii \tr(T_a[T_b,T_c])$ *) let f_of_rep_id norm rep_t = let lhs = norm *** f 1 2 3 and rhs = f_of_rep rep_t 1 2 3 in eq lhs rhs (* \begin{dubious} Are the normalization factors for the traces of the higher dimensional representations correct? \end{dubious} *) (* \begin{dubious} The traces don't work for the symmetrized generators that we need elsewhere! \end{dubious} *) let suite_lie = "Lie algebra relations" >::: [ "[t,t]=ift" >:: (fun () -> lie_algebra_id t); "[t8,t8]=ift8" >:: (fun () -> lie_algebra_id t8); "[t6,t6]=ift6" >:: (fun () -> lie_algebra_id t6); "[t10,t10]=ift10" >:: (fun () -> lie_algebra_id t10); "[t15,t15]=ift15" >:: (fun () -> lie_algebra_id t15); "[t3bar,t3bar]=ift3bar" >:: (fun () -> lie_algebra_id t3bar); "[tSAS,tSAS]=iftSAS" >:: (fun () -> lie_algebra_id t_SAS); "[tASA,tASA]=iftASA" >:: (fun () -> lie_algebra_id t_ASA); "[t6,t6]=ift6'" >:: (fun () -> lie_algebra_id (t_unsymmetrized 2)); "[t10,t10]=ift10'" >:: (fun () -> lie_algebra_id (t_unsymmetrized 3)); "[t15,t15]=ift15'" >:: (fun () -> lie_algebra_id (t_unsymmetrized 4)); "[t6,t6]=ift6''" >:: (fun () -> lie_algebra_id t6_trivial); "[t10,t10]=ift10''" >:: (fun () -> lie_algebra_id t10_trivial); "[t15,t15]=ift15''" >:: (fun () -> lie_algebra_id t15_trivial); "if = tr(t[t,t])" >:: (fun () -> f_of_rep_id one t); "2n*if = tr(t8[t8,t8])" >:: (fun () -> f_of_rep_id (two *** nc) t8); "n*if = tr(t6[t6,t6])" >:: (fun () -> f_of_rep_id nc t6_trivial); "n^2*if = tr(t10[t10,t10])" >:: (fun () -> f_of_rep_id (nc *** nc) t10_trivial); "n^3*if = tr(t15[t15,t15])" >:: (fun () -> f_of_rep_id (nc *** nc *** nc) t15_trivial) ] (* \thocwmodulesubsection{Ward Identities} *) (* Testing the color part of basic Ward identities is essentially the same as testing the Lie algebra equations above, but with generators sandwiched between propagators, as in Feynman diagrams, where the relative signs come from the kinematic part of the diagrams after applying the equations of motion.. *) (* First the diagram with the three gluon vertex $\ii f_{abc} D_{cd}^{\text{gluon}} D^{ik} T_d^{kl} D^{lj}$ *) let ward_ft rep_t rep_d a b i j = imag *** f a b (-11) *** gluon (-11) (-12) *** rep_d i (-1) *** rep_t (-12) (-1) (-2) *** rep_d (-2) j (* then one diagram with two gauge couplings $D^{ik} T_c^{kl} D^{lm} T_c^{mn} D^{nj}$ *) let ward_tt1 rep_t rep_d a b i j = rep_d i (-1) *** rep_t a (-1) (-2) *** rep_d (-2) (-3) *** rep_t b (-3) (-4) *** rep_d (-4) j (* finally the difference of exchanged orders: $D^{ik} T_a^{kl} D^{lm} T_b^{mn} D^{nj} -D^{ik} T_b^{kl} D^{lm} T_a^{mn} D^{nj}$ *) let ward_tt rep_t rep_d a b i j = ward_tt1 rep_t rep_d a b i j --- ward_tt1 rep_t rep_d b a i j (* \begin{dubious} The optional [~fudge] factor was used for debugging normalizations. \end{dubious} *) let ward_id ?(fudge=one) rep_t rep_d = let lhs = ward_ft rep_t rep_d 1 2 3 4 and rhs = ward_tt rep_t rep_d 1 2 3 4 in eq lhs (fudge *** rhs) let suite_ward = "Ward identities" >::: [ "fund." >:: (fun () -> ward_id t delta3); "adj." >:: (fun () -> ward_id t8 delta8); "S2" >:: (fun () -> ward_id t6 delta6); "S3" >:: (fun () -> ward_id t10 delta10); "A2" >:: (fun () -> ward_id t3bar delta3bar); "A3" >:: (fun () -> ward_id (t_A 3) (delta_A 3)); "SAS" >:: (fun () -> ward_id t_SAS delta_SAS); "ASA" >:: (fun () -> ward_id t_ASA delta_ASA); "S2'" >:: (fun () -> ward_id ~fudge:two t6_trivial delta6); "S3'" >:: (fun () -> ward_id ~fudge:(int 3) t10_trivial delta10) ] let suite_ward_long = "Ward identities" >::: [ "S4" >:: (fun () -> ward_id t15 delta15); "S4'" >:: (fun () -> ward_id ~fudge:(int 4) t15_trivial delta15) ] (* \thocwmodulesubsection{Jacobi Identities} *) (* $T_aT_bT_c$ *) let prod3 rep_t a b c i j = rep_t a i (-1) *** rep_t b (-1) (-2) *** rep_t c (-2) j (* $[T_a,[T_b,T_c]]$ *) let jacobi1 rep_t a b c i j = (prod3 rep_t a b c i j --- prod3 rep_t a c b i j) --- (prod3 rep_t b c a i j --- prod3 rep_t c b a i j) (* sum of cyclic permutations of $[T_a,[T_b,T_c]]$ *) let jacobi rep_t = sum [jacobi1 rep_t 1 2 3 4 5; jacobi1 rep_t 2 3 1 4 5; jacobi1 rep_t 3 1 2 4 5] let jacobi_id rep_t = assert_zero_vertex (jacobi rep_t) let suite_jacobi = "Jacobi identities" >::: [ "fund." >:: (fun () -> jacobi_id t); "adj." >:: (fun () -> jacobi_id f); "S2" >:: (fun () -> jacobi_id t6); "S3" >:: (fun () -> jacobi_id t10); "A2" >:: (fun () -> jacobi_id (t_A 2)); "A3" >:: (fun () -> jacobi_id (t_A 3)); "SAS" >:: (fun () -> jacobi_id t_SAS); "ASA" >:: (fun () -> jacobi_id t_ASA); "S2'" >:: (fun () -> jacobi_id t6_trivial); "S3'" >:: (fun () -> jacobi_id t10_trivial) ] let suite_jacobi_long = "Jacobi identities" >::: [ "S4" >:: (fun () -> jacobi_id t15); "S4'" >:: (fun () -> jacobi_id t15_trivial) ] (* \thocwmodulesubsection{Casimir Operators} \label{pg:casimir-tests} *) (* We can read of the eigenvalues of the Casimir operators for the adjoint, totally symmetric and totally antisymmetric representations of~$\mathrm{SU}(N)$ from table~II of \texttt{hep-ph/0611341} \begin{subequations} \begin{align} C_2(\text{adj}) &= 2N \\ C_2(S_n) &= \frac{n(N-1)(N+n)}{N} \\ C_2(A_n) &= \frac{n(N-n)(N+1)}{N} \end{align} \end{subequations} adjusted for our normalization. Also from \texttt{arxiv:1912.13302} \begin{equation} C_3(S_1) =(N^2-1)(N^2-4)/N^2=\frac{N_C^4-5N_C^2+4}{N_C^2} \end{equation} *) (* Building blocks $n/N_C$ and $N_C+n$ *) let n_over_nc n = const (LP.ints [ (n, -1) ]) let nc_plus n = const (LP.ints [ (1, 1); (n,0) ]) (* $C_2(S_n) = n/N_C(N_C-1)(N_C+n)$ *) let c2_S n = n_over_nc n *** nc_plus (-1) *** nc_plus n (* $C_2(A_n) = n/N_C(N_C-n)(N_C+1)$ *) let c2_A n = n_over_nc n *** nc_plus (-n) *** nc_plus 1 let casimir_tt i j = c2_S 1 *** delta3 i j let casimir_t6t6 i j = c2_S 2 *** delta6 i j let casimir_t10t10 i j = c2_S 3 *** delta10 i j let casimir_t15t15 i j = c2_S 4 *** delta15 i j let casimir_t3bart3bar i j = c2_A 2 *** delta3bar i j let casimir_tA3tA3 i j = c2_A 3 *** delta_A 3 i j (* $C_2(\text{adj})=2N_C$ *) let ca = LP.ints [(2, 1)] let casimir_ff a b = [(ca, 1 <=> 2); (LP.int (-2), [1=>1; 2=>2])] (* $C_3(S_1)=N_C^2-5+4/N_C^2$ *) let c3f = LP.ints [(1, 2); (-5, 0); (4, -2)] let casimir_ttt i j = const c3f *** delta3 i j let suite_casimir = "Casimir operators" >::: [ "t*t" >:: (fun () -> eq (casimir_tt 1 2) (t (-1) 1 (-2) *** t (-1) (-2) 2)); "t*t*t" >:: (fun () -> eq (casimir_ttt 1 2) (d (-1) (-2) (-3) *** t (-1) 1 (-4) *** t (-2) (-4) (-5) *** t (-3) (-5) 2)); "f*f" >:: (fun () -> eq (casimir_ff 1 2) (minus *** f (-1) 1 (-2) *** f (-1) (-2) 2)); "t6*t6" >:: (fun () -> eq (casimir_t6t6 1 2) (t6 (-1) 1 (-2) *** t6 (-1) (-2) 2)); "t3bar*t3bar" >:: (fun () -> eq (casimir_t3bart3bar 1 2) (t3bar (-1) 1 (-2) *** t3bar (-1) (-2) 2)); "tA3*tA3" >:: (fun () -> eq (casimir_tA3tA3 1 2) (t_A 3 (-1) 1 (-2) *** t_A 3 (-1) (-2) 2)); "t_SAS*t_SAS" >:: (fun () -> eq (const (LP.ints [(3,1); (-9,-1)]) *** delta_SAS 1 2) (t_SAS (-1) 1 (-2) *** t_SAS (-1) (-2) 2)); "t_ASA*t_ASA" >:: (fun () -> eq (const (LP.ints [(3,1); (-9,-1)]) *** delta_ASA 1 2) (t_ASA (-1) 1 (-2) *** t_ASA (-1) (-2) 2)); "t10*t10" >:: (fun () -> eq (casimir_t10t10 1 2) (t10 (-1) 1 (-2) *** t10 (-1) (-2) 2)) ] let suite_casimir_long = "Casimir operators" >::: [ "t15*t15" >:: (fun () -> eq (casimir_t15t15 1 2) (t15 (-1) 1 (-2) *** t15 (-1) (-2) 2)) ] (* \thocwmodulesubsection{Color Sums} *) let suite_colorsums = "(squared) color sums" >::: [ "gluon normalization" >:: (fun () -> eq (delta8 1 2) (delta8 1 (-1) *** gluon (-1) (-2) *** delta8 (-2) 2)); "f*f" >:: (fun () -> let sum_ff = multiply [ f (-11) (-12) (-13); f (-21) (-22) (-23); gluon (-11) (-21); gluon (-12) (-22); gluon (-13) (-23) ] and expected = ints [(2, 3); (-2, 1)] in eq expected sum_ff); "d*d" >:: (fun () -> let sum_dd = multiply [ d (-11) (-12) (-13); d (-21) (-22) (-23); gluon (-11) (-21); gluon (-12) (-22); gluon (-13) (-23) ] and expected = ints [(2, 3); (-10, 1); (8, -1)] in eq expected sum_dd); "f*d" >:: (fun () -> let sum_fd = multiply [ f (-11) (-12) (-13); d (-21) (-22) (-23); gluon (-11) (-21); gluon (-12) (-22); gluon (-13) (-23) ] in assert_zero_vertex sum_fd); "Hgg" >:: (fun () -> let sum_hgg = multiply [ delta8_loop (-11) (-12); delta8_loop (-21) (-22); gluon (-11) (-21); gluon (-12) (-22) ] and expected = ints [(1, 2); (-1, 0)] in eq expected sum_hgg) ] let suite = "Color.SU3" >::: [suite_sum; suite_diff; suite_times; suite_normalization; suite_symmetrization; suite_ghosts; suite_propagators; suite_trace; suite_ff; suite_tf; suite_tt; suite_lie; suite_ward; suite_jacobi; suite_casimir; suite_colorsums] let suite_long = "Color.SU3 long" >::: [suite_ward_long; suite_jacobi_long; suite_casimir_long] end end (* \thocwmodulesection{$\mathrm{U}(N_C)$} *) (* \begin{dubious} This must not be used, because it has not yet been updated to the correctly symmetrized version! \end{dubious} *) module U3 : SU3 = struct module A = Arrow open Arrow.Infix module B = Birdtracks type t = B.t let canonicalize = B.canonicalize let to_string = B.to_string let pp = B.pp let trivial = B.trivial let is_null = B.is_null let null = B.null let const = B.const let one = B.one let two = B.two let int = B.int let half = B.half let third = B.third let fraction = B.fraction let nc = B.nc let over_nc = B.over_nc let minus = B.minus let imag = B.imag let ints = B.ints let sum = B.sum let diff = B.diff let scale = B.scale let times = B.times let multiply = B.multiply let relocate = B.relocate let fuse = B.fuse let f_of_rep = B.f_of_rep let d_of_rep = B.d_of_rep module Infix = B.Infix let delta3 i j = [(LP.int 1, j ==> i)] let delta8 a b = [(LP.int 1, a <=> b)] let delta8_loop = delta8 let gluon a b = delta8 a b let delta6 n m = [ (LP.fraction 2, [(m, 0) >=>> (n, 0); (m, 1) >=>> (n, 1)]); (LP.fraction 2, [(m, 0) >=>> (n, 1); (m, 1) >=>> (n, 0)]) ] let triples = [(0, 1, 2); (1, 2, 0); (2, 0, 1); (2, 1, 0); (0, 2, 1); (1, 0, 2)] let delta10 n m = List.map (fun (i, j, k) -> (LP.fraction 6, [(m, 0) >=>> (n, i); (m, 1) >=>> (n, j); (m, 2) >=>> (n, k)])) triples let t a i j = [ (LP.int 1, [j => a; a => i]) ] let f a b c = [ (LP.imag ( 1), A.cycle [a; b; c]); (LP.imag (-1), A.cycle [a; c; b]) ] let t8 a b c = Birdtracks.Infix.( minus *** imag *** f a b c ) let d a b c = [ (LP.int 1, A.cycle [a; b; c]); (LP.int 1, A.cycle [a; c; b]) ] let incomplete tensor = failwith ("Color.Vertex: " ^ tensor ^ " not supported yet!") let experimental tensor = Printf.eprintf "Color.Vertex: %s support still experimental and untested!\n" tensor - let epsilon i j k = incomplete "epsilon-tensor" - let epsilonbar i j k = incomplete "epsilon-tensor" + let epsilon tips = incomplete "epsilon-tensor" + let epsilon_bar tails = incomplete "epsilon-tensor" let t6 a m n = [ (LP.int ( 1), [(n, 0) >=> a; a =>> (m, 0); (n, 1) >=>> (m, 1)]); (LP.int ( 1), [(n, 1) >=> a; a =>> (m, 0); (n, 0) >=>> (m, 1)]) ] let t10 a m n = [ (LP.int ( 1), [(n, 0) >=> a; a =>> (m, 0); (n, 1) >=>> (m, 1); (n, 2) >=>> (m, 2)]); (LP.int (-1), [(n, 0) >=>> (m, 0); (n, 1) >=>> (m, 1); (n, 2) >=>> (m, 2)]) ] let k6 m i j = experimental "k6-tensor"; [ (LP.int 1, [(m, 0) >=> i; (m, 1) >=> j]); (LP.int 1, [(m, 1) >=> i; (m, 0) >=> j]) ] let k6bar m i j = experimental "k6-tensor"; [ (LP.int 1, [i =>> (m, 0); j =>> (m, 1)]); (LP.int 1, [i =>> (m, 1); j =>> (m, 0)]) ] let delta_of_tableau t i j = incomplete "delta_of_tableau" let t_of_tableau tableau a k l = incomplete "t_of_tableau" (* \thocwmodulesubsection{Unit Tests} *) module Test : Test = struct open OUnit open Birdtracks open Infix let suite_lie = "Lie algebra relations" >::: [ "if = tr(t[t,t])" >:: (fun () -> eq (f 1 2 3) (f_of_rep t 1 2 3)) ] (* $N_C=N_C^2/N_C$ *) let cf = LP.ints [(1, 1)] let casimir_tt i j = [(cf, i ==> j)] let suite_casimir = "Casimir operators" >::: [ "t*t" >:: (fun () -> eq (casimir_tt 2 1) (t (-1) (-2) 2 *** t (-1) 1 (-2))) ] let suite = "Color.U3" >::: [suite_lie; suite_casimir] let suite_long = "Color.U3 long" >::: [] end end module Vertex = SU3 Index: trunk/omega/src/UFO.ml =================================================================== --- trunk/omega/src/UFO.ml (revision 8848) +++ trunk/omega/src/UFO.ml (revision 8849) @@ -1,2944 +1,2944 @@ (* UFO.ml -- Copyright (C) 1999-2022 by Wolfgang Kilian Thorsten Ohl Juergen Reuter with contributions from Christian Speckner WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) (* Unfortunately, \texttt{ocamlweb} will not typeset all multi character operators nicely. E.\,g.~\verb+f @< g+ comes out as [f @< g]. *) let (<*>) f g x = f (g x) let (<**>) f g x y = f (g x y) module SMap = Map.Make (struct type t = string let compare = compare end) module SSet = Sets.String module CMap = Map.Make (struct type t = string let compare = ThoString.compare_caseless end) module CSet = Sets.String_Caseless let error_in_string text start_pos end_pos = let i = start_pos.Lexing.pos_cnum and j = end_pos.Lexing.pos_cnum in String.sub text i (j - i) let error_in_file name start_pos end_pos = Printf.sprintf "%s:%d.%d-%d.%d" name start_pos.Lexing.pos_lnum (start_pos.Lexing.pos_cnum - start_pos.Lexing.pos_bol) end_pos.Lexing.pos_lnum (end_pos.Lexing.pos_cnum - end_pos.Lexing.pos_bol) let parse_string text = try UFO_parser.file UFO_lexer.token (UFO_lexer.init_position "" (Lexing.from_string text)) with | UFO_tools.Lexical_Error (msg, start_pos, end_pos) -> invalid_arg (Printf.sprintf "lexical error (%s) at: `%s'" msg (error_in_string text start_pos end_pos)) | UFO_syntax.Syntax_Error (msg, start_pos, end_pos) -> invalid_arg (Printf.sprintf "syntax error (%s) at: `%s'" msg (error_in_string text start_pos end_pos)) | Parsing.Parse_error -> invalid_arg ("parse error: " ^ text) exception File_missing of string let parse_file name = let ic = try open_in name with | Sys_error msg as exc -> if msg = name ^ ": No such file or directory" then raise (File_missing name) else raise exc in let result = begin try UFO_parser.file UFO_lexer.token (UFO_lexer.init_position name (Lexing.from_channel ic)) with | UFO_tools.Lexical_Error (msg, start_pos, end_pos) -> begin close_in ic; invalid_arg (Printf.sprintf "%s: lexical error (%s)" (error_in_file name start_pos end_pos) msg) end | UFO_syntax.Syntax_Error (msg, start_pos, end_pos) -> begin close_in ic; invalid_arg (Printf.sprintf "%s: syntax error (%s)" (error_in_file name start_pos end_pos) msg) end | Parsing.Parse_error -> begin close_in ic; invalid_arg ("parse error: " ^ name) end end in close_in ic; result (* These are the contents of the Python files after lexical analysis as context-free variable declarations, before any semantic interpretation. *) module type Files = sig type t = private { particles : UFO_syntax.t; couplings : UFO_syntax.t; coupling_orders : UFO_syntax.t; vertices : UFO_syntax.t; lorentz : UFO_syntax.t; parameters : UFO_syntax.t; propagators : UFO_syntax.t; decays : UFO_syntax.t } val parse_directory : string -> t end module Files : Files = struct type t = { particles : UFO_syntax.t; couplings : UFO_syntax.t; coupling_orders : UFO_syntax.t; vertices : UFO_syntax.t; lorentz : UFO_syntax.t; parameters : UFO_syntax.t; propagators : UFO_syntax.t; decays : UFO_syntax.t } let parse_directory dir = let filename stem = Filename.concat dir (stem ^ ".py") in let parse stem = parse_file (filename stem) in let parse_optional stem = try parse stem with File_missing _ -> [] in { particles = parse "particles"; couplings = parse "couplings"; coupling_orders = parse_optional "coupling_orders"; vertices = parse "vertices"; lorentz = parse "lorentz"; parameters = parse "parameters"; propagators = parse_optional "propagators"; decays = parse_optional "decays" } end let dump_file pfx f = List.iter (fun s -> print_endline (pfx ^ ": " ^ s)) (UFO_syntax.to_strings f) type charge = | Q_Integer of int | Q_Fraction of int * int let charge_to_string = function | Q_Integer i -> Printf.sprintf "%d" i | Q_Fraction (n, d) -> Printf.sprintf "%d/%d" n d module S = UFO_syntax let find_attrib name attribs = try (List.find (fun a -> name = a.S.a_name) attribs).S.a_value with | Not_found -> failwith ("UFO.find_attrib: \"" ^ name ^ "\" not found") let find_attrib name attribs = (List.find (fun a -> name = a.S.a_name) attribs).S.a_value let name_to_string ?strip name = let stripped = begin match strip, List.rev name with | Some pfx, head :: tail -> if pfx = head then tail else failwith ("UFO.name_to_string: expected prefix '" ^ pfx ^ "', got '" ^ head ^ "'") | _, name -> name end in String.concat "." stripped let name_attrib ?strip name attribs = match find_attrib name attribs with | S.Name n -> name_to_string ?strip n | _ -> invalid_arg ("UFO.name_attrib: " ^ name) let integer_attrib name attribs = match find_attrib name attribs with | S.Integer i -> i | _ -> invalid_arg ("UFO.integer_attrib: " ^ name) let charge_attrib name attribs = match find_attrib name attribs with | S.Integer i -> Q_Integer i | S.Fraction (n, d) -> Q_Fraction (n, d) | _ -> invalid_arg ("UFO.charge_attrib: " ^ name) let string_attrib name attribs = match find_attrib name attribs with | S.String s -> s | _ -> invalid_arg ("UFO.string_attrib: " ^ name) let string_expr_attrib name attribs = match find_attrib name attribs with | S.Name n -> [S.Macro n] | S.String s -> [S.Literal s] | S.String_Expr e -> e | _ -> invalid_arg ("UFO.string_expr_attrib: " ^ name) let boolean_attrib name attribs = try match ThoString.lowercase (name_attrib name attribs) with | "true" -> true | "false" -> false | _ -> invalid_arg ("UFO.boolean_attrib: " ^ name) with | Not_found -> false type value = | Integer of int | Fraction of int * int | Float of float | Expr of UFOx.Expr.t | Name of string list let map_expr f default = function | Integer _ | Fraction (_, _) | Float _ | Name _ -> default | Expr e -> f e let variables = map_expr UFOx.Expr.variables CSet.empty let functions = map_expr UFOx.Expr.functions CSet.empty let add_to_set_in_map key element map = let set = try CMap.find key map with Not_found -> CSet.empty in CMap.add key (CSet.add element set) map (* Add all variables in [value] to the [map] from variables to the names in which they appear, indicating that [name] depends on these variables. *) let dependency name value map = CSet.fold (fun variable acc -> add_to_set_in_map variable name acc) (variables value) map let dependencies name_value_list = List.fold_left (fun acc (name, value) -> dependency name value acc) CMap.empty name_value_list let dependency_to_string (variable, appearences) = Printf.sprintf "%s -> {%s}" variable (String.concat ", " (CSet.elements appearences)) let dependencies_to_strings map = List.map dependency_to_string (CMap.bindings map) let expr_to_string = UFOx.Value.to_string <*> UFOx.Value.of_expr let value_to_string = function | Integer i -> Printf.sprintf "%d" i | Fraction (n, d) -> Printf.sprintf "%d/%d" n d | Float x -> string_of_float x | Expr e -> "'" ^ expr_to_string e ^ "'" | Name n -> name_to_string n let value_to_expr substitutions = function | Integer i -> Printf.sprintf "%d" i | Fraction (n, d) -> Printf.sprintf "%d/%d" n d | Float x -> string_of_float x | Expr e -> expr_to_string (substitutions e) | Name n -> name_to_string n let value_to_coupling substitutions atom = function | Integer i -> Coupling.Integer i | Fraction (n, d) -> Coupling.Quot (Coupling.Integer n, Coupling.Integer d) | Float x -> Coupling.Float x | Expr e -> UFOx.Value.to_coupling atom (UFOx.Value.of_expr (substitutions e)) | Name n -> failwith "UFO.value_to_coupling: Name not supported yet!" let value_to_numeric = function | Integer i -> Printf.sprintf "%d" i | Fraction (n, d) -> Printf.sprintf "%g" (float n /. float d) | Float x -> Printf.sprintf "%g" x | Expr e -> invalid_arg ("UFO.value_to_numeric: expr = " ^ (expr_to_string e)) | Name n -> invalid_arg ("UFO.value_to_numeric: name = " ^ name_to_string n) let value_to_float = function | Integer i -> float i | Fraction (n, d) -> float n /. float d | Float x -> x | Expr e -> invalid_arg ("UFO.value_to_float: string = " ^ (expr_to_string e)) | Name n -> invalid_arg ("UFO.value_to_float: name = " ^ name_to_string n) let value_attrib name attribs = match find_attrib name attribs with | S.Integer i -> Integer i | S.Fraction (n, d) -> Fraction (n, d) | S.Float x -> Float x | S.String s -> Expr (UFOx.Expr.of_string s) | S.Name n -> Name n | _ -> invalid_arg ("UFO.value_attrib: " ^ name) let string_list_attrib name attribs = match find_attrib name attribs with | S.String_List l -> l | _ -> invalid_arg ("UFO.string_list_attrib: " ^ name) let name_list_attrib ~strip name attribs = match find_attrib name attribs with | S.Name_List l -> List.map (name_to_string ~strip) l | _ -> invalid_arg ("UFO.name_list_attrib: " ^ name) let integer_list_attrib name attribs = match find_attrib name attribs with | S.Integer_List l -> l | _ -> invalid_arg ("UFO.integer_list_attrib: " ^ name) let order_dictionary_attrib name attribs = match find_attrib name attribs with | S.Order_Dictionary d -> d | _ -> invalid_arg ("UFO.order_dictionary_attrib: " ^ name) let coupling_dictionary_attrib ~strip name attribs = match find_attrib name attribs with | S.Coupling_Dictionary d -> List.map (fun (i, j, c) -> (i, j, name_to_string ~strip c)) d | _ -> invalid_arg ("UFO.coupling_dictionary_attrib: " ^ name) let decay_dictionary_attrib name attribs = match find_attrib name attribs with | S.Decay_Dictionary d -> List.map (fun (p, w) -> (List.map List.hd p, w)) d | _ -> invalid_arg ("UFO.decay_dictionary_attrib: " ^ name) (*i The following doesn't typecheck in applications, even with type annotations ... let attrib_handlers : type attribs value. string -> string -> attribs -> ((string -> attribs -> value) -> string -> value) * ((string -> attribs -> value) -> string -> value -> value) = fun kind symbol attribs -> let required query name = try query name attribs with | Not_found -> invalid_arg (Printf.sprintf "fatal UFO error: mandatory attribute `%s' missing for %s `%s'!" name kind symbol) and optional query name default = try query name attribs with | Not_found -> default in (required, optional) i*) let required_handler kind symbol attribs query name = try query name attribs with | Not_found -> invalid_arg (Printf.sprintf "fatal UFO error: mandatory attribute `%s' missing for %s `%s'!" name kind symbol) let optional_handler attribs query name default = try query name attribs with | Not_found -> default (* The UFO paper~\cite{Degrande:2011ua} is not clear on the question whether the \texttt{name} attribute of an instance must match its Python name. While the examples appear to imply this, there are examples of UFO files in the wild that violate this constraint. *) let warn_symbol_name file symbol name = if name <> symbol then Printf.eprintf "UFO: warning: symbol '%s' <> name '%s' in %s.py: \ while legal in UFO, it is unusual and can cause problems!\n" symbol name file let valid_fortran_id kind name = if not (ThoString.valid_fortran_id name) then invalid_arg (Printf.sprintf "fatal UFO error: the %s `%s' is not a valid fortran id!" kind name) let map_to_alist map = SMap.fold (fun key value acc -> (key, value) :: acc) map [] let keys map = SMap.fold (fun key _ acc -> key :: acc) map [] let keys_caseless map = CMap.fold (fun key _ acc -> key :: acc) map [] let values map = SMap.fold (fun _ value acc -> value :: acc) map [] module SKey = struct type t = string let hash = Hashtbl.hash let equal = (=) end module SHash = Hashtbl.Make (SKey) module type Particle = sig type t = private { pdg_code : int; name : string; antiname : string; spin : UFOx.Lorentz.r; color : UFOx.Color.r; mass : string; width : string; propagator : string option; texname : string; antitexname : string; charge : charge; ghost_number : int; lepton_number : int; y : charge; goldstone : bool; propagating : bool; (* NOT HANDLED YET! *) line : string option; (* NOT HANDLED YET! *) is_anti : bool } val of_file : S.t -> t SMap.t val to_string : string -> t -> string val conjugate : t -> t val force_spinor : t -> t val force_conjspinor : t -> t val force_majorana : t -> t val is_majorana : t -> bool val is_ghost : t -> bool val is_goldstone : t -> bool val is_physical : t -> bool val filter : (t -> bool) -> t SMap.t -> t SMap.t end module Particle : Particle = struct type t = { pdg_code : int; name : string; antiname : string; spin : UFOx.Lorentz.r; color : UFOx.Color.r; mass : string; width : string; propagator : string option; texname : string; antitexname : string; charge : charge; ghost_number : int; lepton_number : int; y : charge; goldstone : bool; propagating : bool; (* NOT HANDLED YET! *) line : string option; (* NOT HANDLED YET! *) is_anti : bool } let to_string symbol p = Printf.sprintf "particle: %s => [pdg = %d, name = '%s'/'%s', \ spin = %s, color = %s, \ mass = %s, width = %s,%s \ Q = %s, G = %d, L = %d, Y = %s, \ TeX = '%s'/'%s'%s]" symbol p.pdg_code p.name p.antiname (UFOx.Lorentz.rep_to_string p.spin) (UFOx.Color.rep_to_string p.color) p.mass p.width (match p.propagator with | None -> "" | Some p -> " propagator = " ^ p ^ ",") (charge_to_string p.charge) p.ghost_number p.lepton_number (charge_to_string p.y) p.texname p.antitexname (if p.goldstone then ", GB" else "") let conjugate_charge = function | Q_Integer i -> Q_Integer (-i) | Q_Fraction (n, d) -> Q_Fraction (-n, d) let is_neutral p = (p.name = p.antiname) (* We \emph{must not} mess with [pdg_code] and [color] if the particle is neutral! *) let conjugate p = if is_neutral p then p else { pdg_code = - p.pdg_code; name = p.antiname; antiname = p.name; spin = UFOx.Lorentz.rep_conjugate p.spin; color = UFOx.Color.rep_conjugate p.color; mass = p.mass; width = p.width; propagator = p.propagator; texname = p.antitexname; antitexname = p.texname; charge = conjugate_charge p.charge; ghost_number = - p.ghost_number; lepton_number = - p.lepton_number; y = conjugate_charge p.y; goldstone = p.goldstone; propagating = p.propagating; line = p.line; is_anti = not p.is_anti } let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Particle" ], attribs -> let required query name = required_handler "particle" symbol attribs query name and optional query name default = optional_handler attribs query name default in let name = required string_attrib "name" and antiname = required string_attrib "antiname" in let neutral = (name = antiname) in let pdg_code = required integer_attrib "pdg_code" in SMap.add symbol { (* The required attributes per UFO docs. *) pdg_code; name; antiname; spin = UFOx.Lorentz.rep_of_int neutral (required integer_attrib "spin"); color = UFOx.Color.rep_of_int neutral (required integer_attrib "color"); mass = required (name_attrib ~strip:"Param") "mass"; width = required (name_attrib ~strip:"Param") "width"; texname = required string_attrib "texname"; antitexname = required string_attrib "antitexname"; charge = required charge_attrib "charge"; (* The optional attributes per UFO docs. *) ghost_number = optional integer_attrib "GhostNumber" 0; lepton_number = optional integer_attrib "LeptonNumber" 0; y = optional charge_attrib "Y" (Q_Integer 0); goldstone = optional boolean_attrib "goldstone" false; propagating = optional boolean_attrib "propagating" true; line = (try Some (name_attrib "line" attribs) with _ -> None); (* Undocumented extensions. *) propagator = (try Some (name_attrib ~strip:"Prop" "propagator" attribs) with _ -> None); (* O'Mega extensions. *) (* Instead of ``first come is particle'' rely on a negative PDG code to identify antiparticles. *) is_anti = pdg_code < 0 } map | [ "anti"; p ], [] -> begin try SMap.add symbol (conjugate (SMap.find p map)) map with | Not_found -> invalid_arg ("Particle.of_file: " ^ p ^ ".anti() not yet defined!") end | _ -> invalid_arg ("Particle.of_file: " ^ name_to_string d.S.kind) let of_file particles = List.fold_left of_file1 SMap.empty particles let is_spinor p = match UFOx.Lorentz.omega p.spin with | Coupling.Spinor | Coupling.ConjSpinor | Coupling.Majorana -> true | _ -> false (* \begin{dubious} TODO: this is a bit of a hack: try to expose the type [UFOx.Lorentz_Atom'.r] instead. \end{dubious} *) let force_spinor p = if is_spinor p then { p with spin = UFOx.Lorentz.rep_of_int false 2 } else p let force_conjspinor p = if is_spinor p then { p with spin = UFOx.Lorentz.rep_of_int false (-2) } else p let force_majorana p = if is_spinor p then { p with spin = UFOx.Lorentz.rep_of_int true 2 } else p let is_majorana p = match UFOx.Lorentz.omega p.spin with | Coupling.Majorana | Coupling.Vectorspinor | Coupling.Maj_Ghost -> true | _ -> false let is_ghost p = p.ghost_number <> 0 let is_goldstone p = p.goldstone let is_physical p = not (is_ghost p || is_goldstone p) let filter predicate map = SMap.filter (fun symbol p -> predicate p) map end module type UFO_Coupling = sig type t = private { name : string; value : UFOx.Expr.t; order : (string * int) list } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module UFO_Coupling : UFO_Coupling = struct type t = { name : string; value : UFOx.Expr.t; order : (string * int) list } let order_to_string orders = String.concat ", " (List.map (fun (s, i) -> Printf.sprintf "'%s':%d" s i) orders) let to_string symbol c = Printf.sprintf "coupling: %s => [name = '%s', value = '%s', order = [%s]]" symbol c.name (expr_to_string c.value) (order_to_string c.order) let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Coupling" ], attribs -> let required query name = required_handler "coupling" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "couplings" symbol name; valid_fortran_id "coupling" name; SMap.add symbol { name; value = UFOx.Expr.of_string (required string_attrib "value"); order = required order_dictionary_attrib "order" } map | _ -> invalid_arg ("UFO_Coupling.of_file: " ^ name_to_string d.S.kind) let of_file couplings = List.fold_left of_file1 SMap.empty couplings end module type Coupling_Order = sig type t = private { name : string; expansion_order : int; hierarchy : int } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Coupling_Order : Coupling_Order = struct type t = { name : string; expansion_order : int; hierarchy : int } let to_string symbol c = Printf.sprintf "coupling_order: %s => [name = '%s', \ expansion_order = '%d', \ hierarchy = %d]" symbol c.name c.expansion_order c.hierarchy let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "CouplingOrder" ], attribs -> let required query name = required_handler "coupling order" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "coupling_orders" symbol name; SMap.add symbol { name; expansion_order = required integer_attrib "expansion_order"; hierarchy = required integer_attrib "hierarchy" } map | _ -> invalid_arg ("Coupling_order.of_file: " ^ name_to_string d.S.kind) let of_file coupling_orders = List.fold_left of_file1 SMap.empty coupling_orders end module type Lorentz_UFO = sig (* If the \texttt{name} attribute of a \texttt{Lorentz} object does \emph{not} match the the name of the object, we need the latter for weeding out unused Lorentz structures (see [Vertex.contains] below). Therefore, we keep it around. *) type t = private { name : string; symbol : string; spins : int list; structure : UFOx.Lorentz.t } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Lorentz_UFO : Lorentz_UFO = struct type t = { name : string; symbol : string; spins : int list; structure : UFOx.Lorentz.t } let to_string symbol l = Printf.sprintf "lorentz: %s => [name = '%s', spins = [%s], \ structure = %s]" symbol l.name (String.concat ", " (List.map string_of_int l.spins)) (UFOx.Lorentz.to_string l.structure) let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Lorentz" ], attribs -> let required query name = required_handler "lorentz" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "lorentz" symbol name; valid_fortran_id "lorentz" symbol; SMap.add symbol { name; symbol; spins = required integer_list_attrib "spins"; structure = UFOx.Lorentz.of_string (required string_attrib "structure") } map | _ -> invalid_arg ("Lorentz.of_file: " ^ name_to_string d.S.kind) let of_file lorentz = List.fold_left of_file1 SMap.empty lorentz end module type Vertex = sig type lcc = private (* Lorentz-color-coupling *) { lorentz : string; color : UFOx.Color.t; coupling : string } type t = private { name : string; particles : string array; lcc : lcc list } val of_file : Particle.t SMap.t -> S.t -> t SMap.t val to_string : string -> t -> string val to_string_expanded : Lorentz_UFO.t SMap.t -> UFO_Coupling.t SMap.t -> t -> string val contains : Particle.t SMap.t -> (Particle.t -> bool) -> t -> bool val filter : (t -> bool) -> t SMap.t -> t SMap.t end module Vertex : Vertex = struct type lcc = { lorentz : string; color : UFOx.Color.t; coupling : string } type t = { name : string; particles : string array; lcc : lcc list } let to_string symbol c = Printf.sprintf "vertex: %s => [name = '%s', particles = [%s], \ lorentz-color-couplings = [%s]" symbol c.name (String.concat ", " (Array.to_list c.particles)) (String.concat ", " (List.map (fun lcc -> Printf.sprintf "%s * %s * %s" lcc.coupling lcc.lorentz (UFOx.Color.to_string lcc.color)) c.lcc)) let to_string_expanded lorentz couplings c = let expand_lorentz s = try UFOx.Lorentz.to_string (SMap.find s lorentz).Lorentz_UFO.structure with | Not_found -> "?" in Printf.sprintf "expanded: [%s] -> { lorentz-color-couplings = [%s] }" (String.concat ", " (Array.to_list c.particles)) (String.concat ", " (List.map (fun lcc -> Printf.sprintf "%s * %s * %s" lcc.coupling (expand_lorentz lcc.lorentz) (UFOx.Color.to_string lcc.color)) c.lcc)) let contains particles predicate v = let p = v.particles in let rec contains' i = if i < 0 then false else if predicate (SMap.find p.(i) particles) then true else contains' (pred i) in contains' (Array.length p - 1) let force_adj_identity1 adj_indices = function | UFOx.Color_Atom.Identity (a, b) as atom -> begin match List.mem a adj_indices, List.mem b adj_indices with | true, true -> UFOx.Color_Atom.Identity8 (a, b) | false, false -> atom | true, false | false, true -> invalid_arg "force_adj_identity: mixed representations!" end | atom -> atom let force_adj_identity adj_indices tensor = UFOx.Color.map_atoms (force_adj_identity1 adj_indices) tensor let find_adj_indices map particles = let adj_indices = ref [] in Array.iteri (fun i p -> (* We must pattern match against the O'Mega representation, because [UFOx.Color.r] is abstract. *) match UFOx.Color.omega (SMap.find p map).Particle.color with | Color.AdjSUN _ -> adj_indices := succ i :: !adj_indices | _ -> ()) particles; !adj_indices let classify_color_indices map particles = let fund_indices = ref [] and conj_indices = ref [] and adj_indices = ref [] in Array.iteri (fun i p -> (* We must pattern match against the O'Mega representation, because [UFOx.Color.r] is abstract. *) match UFOx.Color.omega (SMap.find p map).Particle.color with | Color.SUN n -> if n > 0 then fund_indices := succ i :: !fund_indices else if n < 0 then conj_indices := succ i :: !conj_indices else failwith "classify_color_indices: SU(0)" | Color.AdjSUN n -> if n <> 0 then adj_indices := succ i :: !adj_indices else failwith "classify_color_indices: SU(0)" | _ -> ()) particles; (!fund_indices, !conj_indices, !adj_indices) (* FIXME: would have expected the opposite order \ldots *) let force_identity1 (fund_indices, conj_indices, adj_indices) = function | UFOx.Color_Atom.Identity (a, b) as atom -> if List.mem a fund_indices then begin if List.mem b conj_indices then UFOx.Color_Atom.Identity (b, a) else invalid_arg "force_adj_identity: mixed representations!" end else if List.mem a conj_indices then begin if List.mem b fund_indices then UFOx.Color_Atom.Identity (a, b) else invalid_arg "force_adj_identity: mixed representations!" end else if List.mem a adj_indices then begin if List.mem b adj_indices then UFOx.Color_Atom.Identity8 (a, b) else invalid_arg "force_adj_identity: mixed representations!" end else atom | atom -> atom let force_identity indices tensor = UFOx.Color.map_atoms (force_identity1 indices) tensor (* Here we don't have the Lorentz structures available yet. Thus we set [fermion_lines = []] for now and correct this later. *) let of_file1 particle_map map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Vertex" ], attribs -> let required query name = required_handler "vertex" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "vertices" symbol name; let particles = Array.of_list (required (name_list_attrib ~strip:"P") "particles") in let color = let indices = classify_color_indices particle_map particles in Array.of_list (List.map (force_identity indices <*> UFOx.Color.of_string) (required string_list_attrib "color")) and lorentz = Array.of_list (required (name_list_attrib ~strip:"L") "lorentz") and couplings_alist = required (coupling_dictionary_attrib ~strip:"C") "couplings" in let lcc = List.map (fun (i, j, c) -> { lorentz = lorentz.(j); color = color.(i); coupling = c }) couplings_alist in SMap.add symbol { name; particles; lcc } map | _ -> invalid_arg ("Vertex.of_file: " ^ name_to_string d.S.kind) let of_file particles vertices = List.fold_left (of_file1 particles) SMap.empty vertices let filter predicate map = SMap.filter (fun symbol p -> predicate p) map end module type Parameter = sig type nature = private Internal | External type ptype = private Real | Complex type t = private { name : string; nature : nature; ptype : ptype; value : value; texname : string; lhablock : string option; lhacode : int list option; sequence : int } val of_file : S.t -> t SMap.t val to_string : string -> t -> string val missing : string -> t end module Parameter : Parameter = struct type nature = Internal | External let nature_to_string = function | Internal -> "internal" | External -> "external" let nature_of_string = function | "internal" -> Internal | "external" -> External | s -> invalid_arg ("Parameter.nature_of_string: " ^ s) type ptype = Real | Complex let ptype_to_string = function | Real -> "real" | Complex -> "complex" let ptype_of_string = function | "real" -> Real | "complex" -> Complex | s -> invalid_arg ("Parameter.ptype_of_string: " ^ s) type t = { name : string; nature : nature; ptype : ptype; value : value; texname : string; lhablock : string option; lhacode : int list option; sequence : int } let to_string symbol p = Printf.sprintf "parameter: %s => [#%d, name = '%s', nature = %s, type = %s, \ value = %s, texname = '%s', \ lhablock = %s, lhacode = [%s]]" symbol p.sequence p.name (nature_to_string p.nature) (ptype_to_string p.ptype) (value_to_string p.value) p.texname (match p.lhablock with None -> "???" | Some s -> s) (match p.lhacode with | None -> "" | Some c -> String.concat ", " (List.map string_of_int c)) let of_file1 (map, n) d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Parameter" ], attribs -> let required query name = required_handler "particle" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "parameters" symbol name; valid_fortran_id "parameter" name; (SMap.add symbol { name; nature = nature_of_string (required string_attrib "nature"); ptype = ptype_of_string (required string_attrib "type"); value = required value_attrib "value"; texname = required string_attrib "texname"; lhablock = (try Some (string_attrib "lhablock" attribs) with Not_found -> None); lhacode = (try Some (integer_list_attrib "lhacode" attribs) with Not_found -> None); sequence = n } map, succ n) | _ -> invalid_arg ("Parameter.of_file: " ^ name_to_string d.S.kind) let of_file parameters = let map, _ = List.fold_left of_file1 (SMap.empty, 0) parameters in map let missing name = { name; nature = External; ptype = Real; value = Integer 0; texname = Printf.sprintf "\\texttt{%s}" name; lhablock = None; lhacode = None; sequence = 0 } end (* Macros are encoded as a special [S.declaration] with [S.kind = "$"]. This is slightly hackish, but general enough and the overhead of a special union type is probably not worth the effort. *) module type Macro = sig type t val empty : t (* The domains and codomains are still a bit too much ad hoc, but it does the job. *) val define : t -> string -> S.value -> t val expand_string : t -> string -> S.value val expand_expr : t -> S.string_atom list -> string (* Only for documentation: *) val expand_atom : t -> S.string_atom -> string end module Macro : Macro = struct type t = S.value SMap.t let empty = SMap.empty let define macros name expansion = SMap.add name expansion macros let expand_string macros name = SMap.find name macros let rec expand_atom macros = function | S.Literal s -> s | S.Macro [name] -> begin try begin match SMap.find name macros with | S.String s -> s | S.String_Expr expr -> expand_expr macros expr | _ -> invalid_arg ("expand_atom: not a string: " ^ name) end with | Not_found -> invalid_arg ("expand_atom: not found: " ^ name) end | S.Macro [] -> invalid_arg "expand_atom: empty" | S.Macro name -> invalid_arg ("expand_atom: compound name: " ^ String.concat "." name) and expand_expr macros expr = String.concat "" (List.map (expand_atom macros) expr) end module type Propagator_UFO = sig type t = (* private *) { name : string; numerator : UFOx.Lorentz.t; denominator : UFOx.Lorentz.t } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Propagator_UFO : Propagator_UFO = struct type t = { name : string; numerator : UFOx.Lorentz.t; denominator : UFOx.Lorentz.t } let to_string symbol p = Printf.sprintf "propagator: %s => [name = '%s', numerator = '%s', \ denominator = '%s']" symbol p.name (UFOx.Lorentz.to_string p.numerator) (UFOx.Lorentz.to_string p.denominator) (* The \texttt{denominator} attribute is optional and there is a default (cf.~\texttt{arXiv:1308.1668}) *) let default_denominator = "P('mu', id) * P('mu', id) \ - Mass(id) * Mass(id) \ + complex(0,1) * Mass(id) * Width(id)" let of_string_with_error_correction symbol num_or_den s = try UFOx.Lorentz.of_string s with | Invalid_argument msg -> begin let fixed = s ^ ")" in try let tensor = UFOx.Lorentz.of_string fixed in Printf.eprintf "UFO.Propagator.of_string: added missing closing parenthesis \ in %s of %s: \"%s\"\n" num_or_den symbol s; tensor with | Invalid_argument _ -> invalid_arg (Printf.sprintf "UFO.Propagator.of_string: %s of %s: %s in \"%s\"\n" num_or_den symbol msg fixed) end let of_file1 (macros, map) d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Propagator" ], attribs -> let required query name = required_handler "particle" symbol attribs query name and optional query name default = optional_handler attribs query name default in let name = required string_attrib "name" in warn_symbol_name "propagators" symbol name; let num_string_expr = required string_expr_attrib "numerator" and den_string = begin match optional find_attrib "denominator" (S.String default_denominator) with | S.String s -> s | S.Name [n] -> begin match Macro.expand_string macros n with | S.String s -> s | _ -> invalid_arg "Propagator.denominator" end | _ -> invalid_arg "Propagator.denominator: " end in let num_string = Macro.expand_expr macros num_string_expr in let numerator = of_string_with_error_correction symbol "numerator" num_string and denominator = of_string_with_error_correction symbol "denominator" den_string in (macros, SMap.add symbol { name; numerator; denominator } map) | [ "$" ], [ macro ] -> begin match macro.S.a_value with | S.String _ as s -> (Macro.define macros symbol s, map); | S.String_Expr expr -> let expanded = S.String (Macro.expand_expr macros expr) in (Macro.define macros symbol expanded, map) | _ -> invalid_arg ("Propagator:of_file: not a string " ^ symbol) end | [ "$" ], [] -> invalid_arg ("Propagator:of_file: empty declaration " ^ symbol) | [ "$" ], _ -> invalid_arg ("Propagator:of_file: multiple declaration " ^ symbol) | _ -> invalid_arg ("Propagator:of_file: " ^ name_to_string d.S.kind) let of_file propagators = let _, propagators' = List.fold_left of_file1 (Macro.empty, SMap.empty) propagators in propagators' end module type Decay = sig type t = private { name : string; particle : string; widths : (string list * string) list } val of_file : S.t -> t SMap.t val to_string : string -> t -> string end module Decay : Decay = struct type t = { name : string; particle : string; widths : (string list * string) list } let width_to_string ws = String.concat ", " (List.map (fun (ps, w) -> "(" ^ String.concat ", " ps ^ ") -> '" ^ w ^ "'") ws) let to_string symbol d = Printf.sprintf "decay: %s => [name = '%s', particle = '%s', widths = [%s]]" symbol d.name d.particle (width_to_string d.widths) let of_file1 map d = let symbol = d.S.name in match d.S.kind, d.S.attribs with | [ "Decay" ], attribs -> let required query name = required_handler "particle" symbol attribs query name in let name = required string_attrib "name" in warn_symbol_name "decays" symbol name; SMap.add symbol { name; particle = required (name_attrib ~strip:"P") "particle"; widths = required decay_dictionary_attrib "partial_widths" } map | _ -> invalid_arg ("Decay.of_file: " ^ name_to_string d.S.kind) let of_file decays = List.fold_left of_file1 SMap.empty decays end (* We can read the spinor representations off the vertices to check for consistency. *) (* \begin{dubious} Note that we have to conjugate the representations! \end{dubious} *) let collect_spinor_reps_of_vertex particles lorentz v sets = List.fold_left (fun sets' lcc -> let l = (SMap.find lcc.Vertex.lorentz lorentz).Lorentz_UFO.structure in List.fold_left (fun (spinors, conj_spinors as sets'') (i, rep) -> let p = v.Vertex.particles.(pred i) in match UFOx.Lorentz.omega rep with | Coupling.ConjSpinor -> (SSet.add p spinors, conj_spinors) | Coupling.Spinor -> (spinors, SSet.add p conj_spinors) | _ -> sets'') sets' (UFOx.Lorentz.classify_indices l)) sets v.Vertex.lcc let collect_spinor_reps_of_vertices particles lorentz vertices = SMap.fold (fun _ v -> collect_spinor_reps_of_vertex particles lorentz v) vertices (SSet.empty, SSet.empty) let lorentz_reps_of_vertex particles v = ThoList.alist_of_list ~predicate:(not <*> UFOx.Lorentz.rep_trivial) ~offset:1 (List.map (fun p -> (* Why do we need to conjugate??? *) UFOx.Lorentz.rep_conjugate (SMap.find p particles).Particle.spin) (Array.to_list v.Vertex.particles)) let rep_compatible rep_vertex rep_particle = let open UFOx.Lorentz in let open Coupling in match omega rep_vertex, omega rep_particle with | (Spinor | ConjSpinor), Majorana -> true | r1, r2 -> r1 = r2 let reps_compatible reps_vertex reps_particles = List.for_all2 (fun (iv, rv) (ip, rp) -> iv = ip && rep_compatible rv rp) reps_vertex reps_particles let check_lorentz_reps_of_vertex particles lorentz v = let reps_particles = List.sort compare (lorentz_reps_of_vertex particles v) in List.iter (fun lcc -> let l = (SMap.find lcc.Vertex.lorentz lorentz).Lorentz_UFO.structure in let reps_vertex = List.sort compare (UFOx.Lorentz.classify_indices l) in if not (reps_compatible reps_vertex reps_particles) then begin Printf.eprintf "%s <> %s [%s]\n" (UFOx.Index.classes_to_string UFOx.Lorentz.rep_to_string reps_particles) (UFOx.Index.classes_to_string UFOx.Lorentz.rep_to_string reps_vertex) v.Vertex.name (* [(Vertex.to_string v.Vertex.name v)] *); (* [invalid_arg "check_lorentz_reps_of_vertex"] *) () end) v.Vertex.lcc let color_reps_of_vertex particles v = ThoList.alist_of_list ~predicate:(not <*> UFOx.Color.rep_trivial) ~offset:1 (List.map (fun p -> (SMap.find p particles).Particle.color) (Array.to_list v.Vertex.particles)) let check_color_reps_of_vertex particles v = let reps_particles = List.sort compare (color_reps_of_vertex particles v) in List.iter (fun lcc -> let reps_vertex = List.sort compare (UFOx.Color.classify_indices lcc.Vertex.color) in if reps_vertex <> reps_particles then begin Printf.printf "%s <> %s\n" (UFOx.Index.classes_to_string UFOx.Color.rep_to_string reps_particles) (UFOx.Index.classes_to_string UFOx.Color.rep_to_string reps_vertex); invalid_arg "check_color_reps_of_vertex" end) v.Vertex.lcc module P = Permutation.Default module type Lorentz = sig type spins = private | Unused | Unique of Coupling.lorentz array | Ambiguous of Coupling.lorentz array SMap.t type t = private { name : string; n : int; spins : spins; structure : UFO_Lorentz.t; fermion_lines : Coupling.fermion_lines; variables : string list } val required_charge_conjugates : t -> t list val permute : P.t -> t -> t val of_lorentz_UFO : Particle.t SMap.t -> Vertex.t SMap.t -> Lorentz_UFO.t SMap.t -> t SMap.t val lorentz_to_string : Coupling.lorentz -> string val to_string : string -> t -> string end module Lorentz : Lorentz = struct let rec lorentz_to_string = function | Coupling.Scalar -> "Scalar" | Coupling.Spinor -> "Spinor" | Coupling.ConjSpinor -> "ConjSpinor" | Coupling.Majorana -> "Majorana" | Coupling.Maj_Ghost -> "Maj_Ghost" | Coupling.Vector -> "Vector" | Coupling.Massive_Vector -> "Massive_Vector" | Coupling.Vectorspinor -> "Vectorspinor" | Coupling.Tensor_1 -> "Tensor_1" | Coupling.Tensor_2 -> "Tensor_2" | Coupling.BRS l -> "BRS(" ^ lorentz_to_string l ^ ")" (* Unlike UFO, O'Mega distinguishes bewteen spinors and conjugate spinors. However, we can inspect the particles in the vertices in which a Lorentz structure is used to determine the correct quantum numbers. Most model files in the real world contain unused Lorentz structures. This is not a problem, we can just ignore them. *) type spins = | Unused | Unique of Coupling.lorentz array | Ambiguous of Coupling.lorentz array SMap.t (* \begin{dubious} Use [UFO_targets.Fortran.fusion_name] below in order to avoid communication problems. Or even move away from strings alltogether. \end{dubious} *) type t = { name : string; n : int; spins : spins; structure : UFO_Lorentz.t; fermion_lines : Coupling.fermion_lines; variables : string list } (* Add one charge conjugated fermion lines. *) let charge_conjugate1 l (ket, bra as fermion_line) = { name = l.name ^ Printf.sprintf "_c%x%x" ket bra; n = l.n; spins = l.spins; structure = UFO_Lorentz.charge_conjugate fermion_line l.structure; fermion_lines = l.fermion_lines; variables = l.variables } (* Add several charge conjugated fermion lines. *) let charge_conjugate l fermion_lines = List.fold_left charge_conjugate1 l fermion_lines (*i let all_charge_conjugates l = List.map (charge_conjugate l) (ThoList.power l.fermion_lines) i*) (* Add all combinations of charge conjugated fermion lines that don't leave the fusion. *) let required_charge_conjugates l = let saturated_fermion_lines = List.filter (fun (ket, bra) -> ket != 1 && bra != 1) l.fermion_lines in List.map (charge_conjugate l) (ThoList.power saturated_fermion_lines) let permute_spins p = function | Unused -> Unused | Unique s -> Unique (P.array p s) | Ambiguous map -> Ambiguous (SMap.map (P.array p) map) (* Note that we apply the \emph{inverse} permutation to the indices in order to match the permutation of the particles/spins. *) let permute_structure n p (l, f) = let permuted = P.array (P.inverse p) (Array.init n succ) in let permute_index i = if i > 0 then UFOx.Index.map_position (fun pos -> permuted.(pred pos)) i else i in (UFO_Lorentz.map_indices permute_index l, UFO_Lorentz.map_fermion_lines permute_index f) let permute p l = let structure, fermion_lines = permute_structure l.n p (l.structure, l.fermion_lines) in { name = l.name ^ "_p" ^ P.to_string (P.inverse p); n = l.n; spins = permute_spins p l.spins; structure; fermion_lines; variables = l.variables } let omega_lorentz_reps n alist = let reps = Array.make n Coupling.Scalar in List.iter (fun (i, rep) -> reps.(pred i) <- UFOx.Lorentz.omega rep) alist; reps let contained lorentz vertex = List.exists (fun lcc1 -> lcc1.Vertex.lorentz = lorentz.Lorentz_UFO.symbol) vertex.Vertex.lcc (* Find all vertices in with the Lorentz structure [lorentz] is used and build a map from those vertices to the O'Mega Lorentz representations inferred from UFO's Lorentz structure and the [particles] involved. Then scan the bindings and check that we have inferred the same Lorentz representation from all vertices. *) let lorentz_reps_of_structure particles vertices lorentz = let uses = SMap.fold (fun name v acc -> if contained lorentz v then SMap.add name (omega_lorentz_reps (Array.length v.Vertex.particles) (lorentz_reps_of_vertex particles v)) acc else acc) vertices SMap.empty in let variants = ThoList.uniq (List.sort compare (List.map snd (SMap.bindings uses))) in match variants with | [] -> Unused | [s] -> Unique s | _ -> Printf.eprintf "UFO.Lorentz.lorentz_reps_of_structure: AMBIGUOUS!\n"; List.iter (fun variant -> Printf.eprintf "UFO.Lorentz.lorentz_reps_of_structure: %s\n" (ThoList.to_string lorentz_to_string (Array.to_list variant))) variants; Ambiguous uses let of_lorentz_tensor spins lorentz = match spins with | Unique s -> begin try Some (UFO_Lorentz.parse (Array.to_list s) lorentz) with | Failure msg -> begin prerr_endline msg; Some (UFO_Lorentz.dummy) end end | Unused -> Printf.eprintf "UFO.Lorentz: stripping unused structure %s\n" (UFOx.Lorentz.to_string lorentz); None | Ambiguous _ -> invalid_arg "UFO.Lorentz.of_lorentz_tensor: Ambiguous" (* NB: if the \texttt{name} attribute of a \texttt{Lorentz} object does \emph{not} match the the name of the object, the former has a better chance to correspond to a valid Fortran name. Therefore we use it. *) let of_lorentz_UFO particles vertices lorentz_UFO = SMap.fold (fun name l acc -> let spins = lorentz_reps_of_structure particles vertices l in match of_lorentz_tensor spins l.Lorentz_UFO.structure with | None -> acc | Some structure -> SMap.add name { name = l.Lorentz_UFO.symbol; n = List.length l.Lorentz_UFO.spins; spins; structure; fermion_lines = UFO_Lorentz.fermion_lines structure; variables = UFOx.Lorentz.variables l.Lorentz_UFO.structure } acc) lorentz_UFO SMap.empty let to_string symbol l = Printf.sprintf "lorentz: %s => [name = '%s', spins = %s, \ structure = %s, fermion_lines = %s]" symbol l.name (match l.spins with | Unique s -> "[" ^ String.concat ", " (List.map lorentz_to_string (Array.to_list s)) ^ "]" | Ambiguous _ -> "AMBIGUOUS!" | Unused -> "UNUSED!") (UFO_Lorentz.to_string l.structure) (UFO_Lorentz.fermion_lines_to_string l.fermion_lines) end (* According to arxiv:1308:1668, there should not be a factor of~$i$ in the numerators of propagators, but the (unused) \texttt{propagators.py} in most models violate this rule! *) let divide_propagators_by_i = ref false module type Propagator = sig type t = (* private *) { name : string; spins : Coupling.lorentz * Coupling.lorentz; numerator : UFO_Lorentz.t; denominator : UFO_Lorentz.t; variables : string list } val of_propagator_UFO : ?majorana:bool -> Propagator_UFO.t -> t val of_propagators_UFO : ?majorana:bool -> Propagator_UFO.t SMap.t -> t SMap.t val transpose : t -> t val to_string : string -> t -> string end module Propagator : Propagator = struct type t = (* private *) { name : string; spins : Coupling.lorentz * Coupling.lorentz; numerator : UFO_Lorentz.t; denominator : UFO_Lorentz.t; variables : string list } let lorentz_rep_at rep_classes i = try UFOx.Lorentz.omega (List.assoc i rep_classes) with | Not_found -> Coupling.Scalar let imaginary = Algebra.QC.make Algebra.Q.null Algebra.Q.unit let scalars = [Coupling.Scalar; Coupling.Scalar] (* If~$51$ and~$52$ show up as indices, we must map $(1,51)\to(1001,2001)$ and $(2,52)\to(1002,2002)$, as per the UFO conventions for Lorentz structures. *) (* \begin{dubious} This does not work yet, because [UFOx.Lorentz.map_indices] affects also the position argument of [P], [Mass] and [Width]. \end{dubious} *) let contains_51_52 tensor = List.exists (fun (i, _) -> i = 51 || i = 52) (UFOx.Lorentz.classify_indices tensor) let remap_51_52 = function | 1 -> 1001 | 51 -> 2001 | 2 -> 1002 | 52 -> 2002 | i -> i let canonicalize_51_52 tensor = if contains_51_52 tensor then UFOx.Lorentz.rename_indices remap_51_52 tensor else tensor let force_majorana = function | Coupling.Spinor | Coupling.ConjSpinor -> Coupling.Majorana | s -> s let string_list_union l1 l2 = Sets.String.elements (Sets.String.union (Sets.String.of_list l1) (Sets.String.of_list l2)) (* In the current conventions, the factor of~$i$ is not included: *) let of_propagator_UFO ?(majorana=false) p = let numerator = canonicalize_51_52 p.Propagator_UFO.numerator in let lorentz_reps = UFOx.Lorentz.classify_indices numerator in let spin1 = lorentz_rep_at lorentz_reps 1 and spin2 = lorentz_rep_at lorentz_reps 2 in let numerator_sans_i = if !divide_propagators_by_i then UFOx.Lorentz.map_coeff (fun q -> Algebra.QC.div q imaginary) numerator else numerator in { name = p.Propagator_UFO.name; spins = if majorana then (force_majorana spin1, force_majorana spin2) else (spin1, spin2); numerator = UFO_Lorentz.parse ~allow_denominator:true [spin1; spin2] numerator_sans_i; denominator = UFO_Lorentz.parse scalars p.Propagator_UFO.denominator; variables = string_list_union (UFOx.Lorentz.variables p.Propagator_UFO.denominator) (UFOx.Lorentz.variables numerator_sans_i) } let of_propagators_UFO ?majorana propagators_UFO = SMap.fold (fun name p acc -> SMap.add name (of_propagator_UFO ?majorana p) acc) propagators_UFO SMap.empty let permute12 = function | 1 -> 2 | 2 -> 1 | n -> n let transpose_positions t = UFOx.Index.map_position permute12 t let transpose p = { name = p.name; spins = (snd p.spins, fst p.spins); numerator = UFO_Lorentz.map_indices transpose_positions p.numerator; denominator = p.denominator; variables = p.variables } let to_string symbol p = Printf.sprintf "propagator: %s => [name = '%s', spin = '(%s, %s)', numerator/I = '%s', \ denominator = '%s']" symbol p.name (Lorentz.lorentz_to_string (fst p.spins)) (Lorentz.lorentz_to_string (snd p.spins)) (UFO_Lorentz.to_string p.numerator) (UFO_Lorentz.to_string p.denominator) end type t = { particles : Particle.t SMap.t; particle_array : Particle.t array; (* for diagnostics *) couplings : UFO_Coupling.t SMap.t; coupling_orders : Coupling_Order.t SMap.t; vertices : Vertex.t SMap.t; lorentz_UFO : Lorentz_UFO.t SMap.t; lorentz : Lorentz.t SMap.t; parameters : Parameter.t SMap.t; propagators_UFO : Propagator_UFO.t SMap.t; propagators : Propagator.t SMap.t; decays : Decay.t SMap.t; nc : int } let use_majorana_spinors = ref false let fallback_to_majorana_if_necessary particles vertices lorentz_UFO = let majoranas = SMap.fold (fun p particle acc -> if Particle.is_majorana particle then SSet.add p acc else acc) particles SSet.empty in let spinors, conj_spinors = collect_spinor_reps_of_vertices particles lorentz_UFO vertices in let ambiguous = SSet.diff (SSet.inter spinors conj_spinors) majoranas in let no_majoranas = SSet.is_empty majoranas and no_ambiguities = SSet.is_empty ambiguous in if no_majoranas && no_ambiguities && not !use_majorana_spinors then (SMap.mapi (fun p particle -> if SSet.mem p spinors then Particle.force_spinor particle else if SSet.mem p conj_spinors then Particle.force_conjspinor particle else particle) particles, false) else begin if !use_majorana_spinors then Printf.eprintf "O'Mega: Majorana fermions requested.\n"; if not no_majoranas then Printf.eprintf "O'Mega: found Majorana fermions!\n"; if not no_ambiguities then Printf.eprintf "O'Mega: found ambiguous spinor representations for %s!\n" (String.concat ", " (SSet.elements ambiguous)); Printf.eprintf "O'Mega: falling back to the Majorana representation for all fermions.\n"; (SMap.map Particle.force_majorana particles, true) end let nc_of_particles particles = let nc_set = List.fold_left (fun nc_set (_, p) -> match UFOx.Color.omega p.Particle.color with | Color.Singlet -> nc_set | Color.SUN nc -> Sets.Int.add (abs nc) nc_set | Color.AdjSUN nc -> Sets.Int.add (abs nc) nc_set) Sets.Int.empty (SMap.bindings particles) in match Sets.Int.elements nc_set with | [] -> 0 | [n] -> n | nc_list -> invalid_arg ("UFO.Model: more than one value of N_C: " ^ String.concat ", " (List.map string_of_int nc_list)) let of_file u = let particles = Particle.of_file u.Files.particles in let vertices = Vertex.of_file particles u.Files.vertices and lorentz_UFO = Lorentz_UFO.of_file u.Files.lorentz and propagators_UFO = Propagator_UFO.of_file u.Files.propagators in let particles, majorana = fallback_to_majorana_if_necessary particles vertices lorentz_UFO in let particle_array = Array.of_list (values particles) and lorentz = Lorentz.of_lorentz_UFO particles vertices lorentz_UFO and propagators = Propagator.of_propagators_UFO ~majorana propagators_UFO in let model = { particles; particle_array; couplings = UFO_Coupling.of_file u.Files.couplings; coupling_orders = Coupling_Order.of_file u.Files.coupling_orders; vertices; lorentz_UFO; lorentz; parameters = Parameter.of_file u.Files.parameters; propagators_UFO; propagators; decays = Decay.of_file u.Files.decays; nc = nc_of_particles particles } in SMap.iter (fun _ v -> check_color_reps_of_vertex model.particles v; check_lorentz_reps_of_vertex model.particles model.lorentz_UFO v) model.vertices; model let parse_directory dir = of_file (Files.parse_directory dir) let dump model = Printf.printf "NC = %d\n" model.nc; SMap.iter (print_endline <**> Particle.to_string) model.particles; SMap.iter (print_endline <**> UFO_Coupling.to_string) model.couplings; SMap.iter (print_endline <**> Coupling_Order.to_string) model.coupling_orders; (* [SMap.iter (print_endline <**> Vertex.to_string) model.vertices;] *) SMap.iter (fun symbol v -> (print_endline <**> Vertex.to_string) symbol v; print_endline (Vertex.to_string_expanded model.lorentz_UFO model.couplings v)) model.vertices; SMap.iter (print_endline <**> Lorentz_UFO.to_string) model.lorentz_UFO; SMap.iter (print_endline <**> Lorentz.to_string) model.lorentz; SMap.iter (print_endline <**> Parameter.to_string) model.parameters; SMap.iter (print_endline <**> Propagator_UFO.to_string) model.propagators_UFO; SMap.iter (print_endline <**> Propagator.to_string) model.propagators; SMap.iter (print_endline <**> Decay.to_string) model.decays; SMap.iter (fun symbol d -> List.iter (fun (_, w) -> ignore (UFOx.Expr.of_string w)) d.Decay.widths) model.decays exception Unhandled of string let unhandled s = raise (Unhandled s) module Model = struct (* NB: we could use [type flavor = Particle.t], but that would be very inefficient, because we will use [flavor] as a key for maps below. *) type flavor = int type constant = string type gauge = unit module M = Modeltools.Mutable (struct type f = flavor type g = gauge type c = constant end) let flavors = M.flavors let external_flavors = M.external_flavors let external_flavors = M.external_flavors let lorentz = M.lorentz let color = M.color let nc = M.nc let propagator = M.propagator let width = M.width let goldstone = M.goldstone let conjugate = M.conjugate let fermion = M.fermion let vertices = M.vertices let fuse2 = M.fuse2 let fuse3 = M.fuse3 let fuse = M.fuse let max_degree = M.max_degree let parameters = M.parameters let flavor_of_string = M.flavor_of_string let flavor_to_string = M.flavor_to_string let flavor_to_TeX = M.flavor_to_TeX let flavor_symbol = M.flavor_symbol let gauge_symbol = M.gauge_symbol let pdg = M.pdg let mass_symbol = M.mass_symbol let width_symbol = M.width_symbol let constant_symbol = M.constant_symbol module Ch = M.Ch let charges = M.charges let rec fermion_of_lorentz = function | Coupling.Spinor -> 1 | Coupling.ConjSpinor -> -1 | Coupling.Majorana -> 2 | Coupling.Maj_Ghost -> 2 | Coupling.Vectorspinor -> 1 | Coupling.Vector | Coupling.Massive_Vector -> 0 | Coupling.Scalar | Coupling.Tensor_1 | Coupling.Tensor_2 -> 0 | Coupling.BRS f -> fermion_of_lorentz f module Q = Algebra.Q module QC = Algebra.QC let dummy_tensor3 = Coupling.Scalar_Scalar_Scalar 1 let dummy_tensor4 = Coupling.Scalar4 1 let triplet p = (p.(0), p.(1), p.(2)) let quartet p = (p.(0), p.(1), p.(2), p.(3)) let half_times q1 q2 = Q.mul (Q.make 1 2) (Q.mul q1 q2) let name g = g.UFO_Coupling.name let fractional_coupling g r = let g = name g in match Q.to_ratio r with | 0, _ -> "0.0_default" | 1, 1 -> g | -1, 1 -> Printf.sprintf "(-%s)" g | n, 1 -> Printf.sprintf "(%d*%s)" n g | 1, d -> Printf.sprintf "(%s/%d)" g d | -1, d -> Printf.sprintf "(-%s/%d)" g d | n, d -> Printf.sprintf "(%d*%s/%d)" n g d let lorentz_of_symbol model symbol = try SMap.find symbol model.lorentz with | Not_found -> invalid_arg ("lorentz_of_symbol: " ^ symbol) let lorentz_UFO_of_symbol model symbol = try SMap.find symbol model.lorentz_UFO with | Not_found -> invalid_arg ("lorentz_UFO_of_symbol: " ^ symbol) let coupling_of_symbol model symbol = try SMap.find symbol model.couplings with | Not_found -> invalid_arg ("coupling_of_symbol: " ^ symbol) let spin_triplet model name = match (lorentz_of_symbol model name).Lorentz.spins with | Lorentz.Unique [|s0; s1; s2|] -> (s0, s1, s2) | Lorentz.Unique _ -> invalid_arg "spin_triplet: wrong number of spins" | Lorentz.Unused -> invalid_arg "spin_triplet: Unused" | Lorentz.Ambiguous _ -> invalid_arg "spin_triplet: Ambiguous" let spin_quartet model name = match (lorentz_of_symbol model name).Lorentz.spins with | Lorentz.Unique [|s0; s1; s2; s3|] -> (s0, s1, s2, s3) | Lorentz.Unique _ -> invalid_arg "spin_quartet: wrong number of spins" | Lorentz.Unused -> invalid_arg "spin_quartet: Unused" | Lorentz.Ambiguous _ -> invalid_arg "spin_quartet: Ambiguous" let spin_multiplet model name = match (lorentz_of_symbol model name).Lorentz.spins with | Lorentz.Unique sarray -> sarray | Lorentz.Unused -> invalid_arg "spin_multiplet: Unused" | Lorentz.Ambiguous _ -> invalid_arg "spin_multiplet: Ambiguous" (* If we have reason to belive that a $\delta_{ab}$-vertex is an effective $\tr(T_aT_b)$-vertex generated at loop level, like~$gg\to H\ldots$ in the SM, we should interpret it as such and use the expression~(6.2) from~\cite{Kilian:2012pz}. *) (* AFAIK, there is no way to distinguish these cases directly in a UFO file. Instead we rely in a heuristic, in which each massless color octet vector particle or ghost is a gluon and colorless scalars are potential Higgses. *) let is_massless p = match ThoString.uppercase p.Particle.mass with | "ZERO" -> true | _ -> false let is_gluon model f = let p = model.particle_array.(f) in match UFOx.Color.omega p.Particle.color, UFOx.Lorentz.omega p.Particle.spin with | Color.AdjSUN _, Coupling.Vector -> is_massless p | Color.AdjSUN _, Coupling.Scalar -> if p.Particle.ghost_number <> 0 then is_massless p else false | _ -> false let is_color_singlet model f = let p = model.particle_array.(f) in match UFOx.Color.omega p.Particle.color with | Color.Singlet -> true | _ -> false let is_higgs_gluon_vertex model p adjoints = if Array.length p > List.length adjoints then List.for_all (fun (i, p) -> if List.mem i adjoints then is_gluon model p else is_color_singlet model p) (ThoList.enumerate 1 (Array.to_list p)) else false let delta8_heuristics model p a b = if is_higgs_gluon_vertex model p [a; b] then Color.Vertex.delta8_loop a b else Color.Vertex.delta8 a b let verbatim_higgs_glue = ref false let translate_color_atom model p = function | UFOx.Color_Atom.Identity (i, j) -> Color.Vertex.delta3 j i | UFOx.Color_Atom.Identity8 (a, b) -> if !verbatim_higgs_glue then Color.Vertex.delta8 a b else delta8_heuristics model p a b | UFOx.Color_Atom.T (a, i, j) -> Color.Vertex.t a i j | UFOx.Color_Atom.F (a, b, c) -> Color.Vertex.f a b c | UFOx.Color_Atom.D (a, b, c) -> Color.Vertex.d a b c - | UFOx.Color_Atom.Epsilon (i, j, k) -> Color.Vertex.epsilon i j k - | UFOx.Color_Atom.EpsilonBar (i, j, k) -> Color.Vertex.epsilonbar i j k + | UFOx.Color_Atom.Epsilon (i, j, k) -> Color.Vertex.epsilon [i; j; k] + | UFOx.Color_Atom.EpsilonBar (i, j, k) -> Color.Vertex.epsilon_bar [i; j; k] | UFOx.Color_Atom.T6 (a, i, j) -> Color.Vertex.t6 a i j | UFOx.Color_Atom.K6 (i, j, k) -> Color.Vertex.k6 i j k | UFOx.Color_Atom.K6Bar (i, j, k) -> Color.Vertex.k6bar i j k let translate_color_term model p = function | [], q -> Color.Vertex.scale q Color.Vertex.one | [atom], q -> Color.Vertex.scale q (translate_color_atom model p atom) | atoms, q -> let atoms = List.map (translate_color_atom model p) atoms in Color.Vertex.scale q (Color.Vertex.multiply atoms) let translate_color model p terms = match terms with | [] -> invalid_arg "translate_color: empty" | [ term ] -> translate_color_term model p term | terms -> Color.Vertex.sum (List.map (translate_color_term model p) terms) let translate_coupling_1 model p lcc = let l = lcc.Vertex.lorentz in let s = Array.to_list (spin_multiplet model l) and fl = (SMap.find l model.lorentz).Lorentz.fermion_lines and c = name (coupling_of_symbol model lcc.Vertex.coupling) in match lcc.Vertex.color with | UFOx.Color.Linear color -> let col = translate_color model p color in (Array.to_list p, Coupling.UFO (QC.unit, l, s, fl, col), c) | UFOx.Color.Ratios _ as color -> invalid_arg ("UFO.Model.translate_coupling: invalid color structure" ^ UFOx.Color.to_string color) let translate_coupling model p lcc = List.map (translate_coupling_1 model p) lcc let long_flavors = ref false module type Lookup = sig type f = private { flavors : flavor list; flavor_of_string : string -> flavor; flavor_of_symbol : string -> flavor; particle : flavor -> Particle.t; flavor_symbol : flavor -> string; conjugate : flavor -> flavor } type flavor_format = | Long | Decimal | Hexadecimal val flavor_format : flavor_format ref val of_model : t -> f end module Lookup : Lookup = struct type f = { flavors : flavor list; flavor_of_string : string -> flavor; flavor_of_symbol : string -> flavor; particle : flavor -> Particle.t; flavor_symbol : flavor -> string; conjugate : flavor -> flavor } type flavor_format = | Long | Decimal | Hexadecimal let flavor_format = ref Hexadecimal (*i let match_pdf_code p1 p2 = p1.Particle.pdg_code = p2.Particle.pdg_code i*) let conjugate_of_particle_array particles = Array.init (Array.length particles) (fun i -> let f' = Particle.conjugate particles.(i) in match ThoArray.match_all f' particles with | [i'] -> i' | [] -> invalid_arg ("no charge conjugate: " ^ f'.Particle.name) | _ -> invalid_arg ("multiple charge conjugates: " ^ f'.Particle.name)) let invert_flavor_array a = let table = SHash.create 37 in Array.iteri (fun i s -> SHash.add table s i) a; (fun name -> try SHash.find table name with | Not_found -> invalid_arg ("not found: " ^ name)) let digits base n = let rec digits' acc n = if n < 1 then acc else digits' (succ acc) (n / base) in if n < 0 then digits' 1 (-n) else if n = 0 then 1 else digits' 0 n let of_model model = let particle_array = Array.of_list (values model.particles) in let conjugate_array = conjugate_of_particle_array particle_array and name_array = Array.map (fun f -> f.Particle.name) particle_array and symbol_array = Array.of_list (keys model.particles) in let flavor_symbol f = begin match !flavor_format with | Long -> symbol_array.(f) | Decimal -> let w = digits 10 (Array.length particle_array - 1) in Printf.sprintf "%0*d" w f | Hexadecimal -> let w = digits 16 (Array.length particle_array - 1) in Printf.sprintf "%0*X" w f end in { flavors = ThoList.range 0 (Array.length particle_array - 1); flavor_of_string = invert_flavor_array name_array; flavor_of_symbol = invert_flavor_array symbol_array; particle = Array.get particle_array; flavor_symbol = flavor_symbol; conjugate = Array.get conjugate_array } end (* \begin{dubious} We appear to need to conjugate all flavors. Why??? \end{dubious} *) let translate_vertices model tables = let vn = List.fold_left (fun acc v -> let p = Array.map tables.Lookup.flavor_of_symbol v.Vertex.particles and lcc = v.Vertex.lcc in let p = Array.map conjugate p in (* FIXME: why? *) translate_coupling model p lcc @ acc) [] (values model.vertices) in ([], [], vn) let propagator_of_lorentz = function | Coupling.Scalar -> Coupling.Prop_Scalar | Coupling.Spinor -> Coupling.Prop_Spinor | Coupling.ConjSpinor -> Coupling.Prop_ConjSpinor | Coupling.Majorana -> Coupling.Prop_Majorana | Coupling.Maj_Ghost -> invalid_arg "UFO.Model.propagator_of_lorentz: SUSY ghosts do not propagate" | Coupling.Vector -> Coupling.Prop_Feynman | Coupling.Massive_Vector -> Coupling.Prop_Unitarity | Coupling.Tensor_2 -> Coupling.Prop_Tensor_2 | Coupling.Vectorspinor -> invalid_arg "UFO.Model.propagator_of_lorentz: Vectorspinor" | Coupling.Tensor_1 -> invalid_arg "UFO.Model.propagator_of_lorentz: Tensor_1" | Coupling.BRS _ -> invalid_arg "UFO.Model.propagator_of_lorentz: no BRST" let filter_unphysical model = let physical_particles = Particle.filter Particle.is_physical model.particles in let physical_particle_array = Array.of_list (values physical_particles) in let physical_vertices = Vertex.filter (not <*> (Vertex.contains model.particles (not <*> Particle.is_physical))) model.vertices in { model with particles = physical_particles; particle_array = physical_particle_array; vertices = physical_vertices } let whizard_constants = SSet.of_list [ "ZERO" ] let filter_constants parameters = List.filter (fun p -> not (SSet.mem (ThoString.uppercase p.Parameter.name) whizard_constants)) parameters let add_name set parameter = CSet.add parameter.Parameter.name set let hardcoded_parameters = CSet.of_list ["cmath.pi"] let missing_parameters input derived couplings = let input_parameters = List.fold_left add_name hardcoded_parameters input in let all_parameters = List.fold_left add_name input_parameters derived in let derived_dependencies = dependencies (List.map (fun p -> (p.Parameter.name, p.Parameter.value)) derived) in let coupling_dependencies = dependencies (List.map (fun p -> (p.UFO_Coupling.name, Expr p.UFO_Coupling.value)) (values couplings)) in let missing_input = CMap.filter (fun parameter derived_parameters -> not (CSet.mem parameter all_parameters)) derived_dependencies and missing = CMap.filter (fun parameter couplings -> not (CSet.mem parameter all_parameters)) coupling_dependencies in CMap.iter (fun parameter derived_parameters -> Printf.eprintf "UFO warning: undefined input parameter %s appears in derived \ parameters {%s}: will be added to the list of input parameters!\n" parameter (String.concat "; " (CSet.elements derived_parameters))) missing_input; CMap.iter (fun parameter couplings -> Printf.eprintf "UFO warning: undefined parameter %s appears in couplings {%s}: \ will be added to the list of input parameters!\n" parameter (String.concat "; " (CSet.elements couplings))) missing; keys_caseless missing_input @ keys_caseless missing let classify_parameters model = let compare_parameters p1 p2 = compare p1.Parameter.sequence p2.Parameter.sequence in let input, derived = List.fold_left (fun (input, derived) p -> match p.Parameter.nature with | Parameter.Internal -> (input, p :: derived) | Parameter.External -> begin match p.Parameter.ptype with | Parameter.Real -> () | Parameter.Complex -> Printf.eprintf "UFO warning: invalid complex declaration of input \ parameter `%s' ignored!\n" p.Parameter.name end; (p :: input, derived)) ([], []) (filter_constants (values model.parameters)) in let additional = missing_parameters input derived model.couplings in (List.sort compare_parameters input @ List.map Parameter.missing additional, List.sort compare_parameters derived) (*i List.iter (fun line -> Printf.eprintf "par: %s\n" line) (dependencies_to_strings derived_dependencies); List.iter (fun line -> Printf.eprintf "coupling: %s\n" line) (dependencies_to_strings coupling_dependencies); i*) let translate_name map name = try SMap.find name map with Not_found -> name let translate_input map p = (translate_name map p.Parameter.name, value_to_float p.Parameter.value) let alpha_s_half e = UFOx.Expr.substitute "aS" (UFOx.Expr.half "aS") e let alpha_s_half_etc map e = UFOx.Expr.rename (map_to_alist map) (alpha_s_half e) let translate_derived map p = let make_atom s = s in let c = make_atom (translate_name map p.Parameter.name) and v = value_to_coupling (alpha_s_half_etc map) make_atom p.Parameter.value in match p.Parameter.ptype with | Parameter.Real -> (Coupling.Real c, v) | Parameter.Complex -> (Coupling.Complex c, v) let translate_coupling_constant map c = let make_atom s = s in (Coupling.Complex c.UFO_Coupling.name, Coupling.Quot (value_to_coupling (alpha_s_half_etc map) make_atom (Expr c.UFO_Coupling.value), Coupling.I)) module Lowercase_Parameters = struct type elt = string type base = string let compare_elt = compare let compare_base = compare let pi = ThoString.lowercase end module Lowercase_Bundle = Bundle.Make (Lowercase_Parameters) let coupling_names model = SMap.fold (fun _ c acc -> c.UFO_Coupling.name :: acc) model.couplings [] let parameter_names model = SMap.fold (fun _ c acc -> c.Parameter.name :: acc) model.parameters [] let ambiguous_parameters model = let all_names = List.rev_append (coupling_names model) (parameter_names model) in let lc_bundle = Lowercase_Bundle.of_list all_names in let lc_set = List.fold_left (fun acc s -> SSet.add s acc) SSet.empty (Lowercase_Bundle.base lc_bundle) and ambiguities = List.filter (fun (_, names) -> List.length names > 1) (Lowercase_Bundle.fibers lc_bundle) in (lc_set, ambiguities) let disambiguate1 lc_set name = let rec disambiguate1' i = let name' = Printf.sprintf "%s_%d" name i in let lc_name' = ThoString.lowercase name' in if SSet.mem lc_name' lc_set then disambiguate1' (succ i) else (SSet.add lc_name' lc_set, name') in disambiguate1' 1 let disambiguate lc_set names = let _, replacements = List.fold_left (fun (lc_set', acc) name -> let lc_set'', name' = disambiguate1 lc_set' name in (lc_set'', SMap.add name name' acc)) (lc_set, SMap.empty) names in replacements let omegalib_names = ["u"; "ubar"; "v"; "vbar"; "eps"] let translate_parameters model = let lc_set, ambiguities = ambiguous_parameters model in let replacements = disambiguate lc_set (ThoList.flatmap snd ambiguities) in SMap.iter (Printf.eprintf "warning: case sensitive parameter names: renaming '%s' -> '%s'\n") replacements; let replacements = List.fold_left (fun acc name -> SMap.add name ("UFO_" ^ name) acc) replacements omegalib_names in let input_parameters, derived_parameters = classify_parameters model and couplings = values model.couplings in { Coupling.input = List.map (translate_input replacements) input_parameters; Coupling.derived = List.map (translate_derived replacements) derived_parameters @ List.map (translate_coupling_constant replacements) couplings; Coupling.derived_arrays = [] } (* UFO requires us to look up the mass parameter to distinguish between massless and massive vectors. TODO: this is a candidate for another lookup table. *) let lorentz_of_particle p = match UFOx.Lorentz.omega p.Particle.spin with | Coupling.Vector -> begin match ThoString.uppercase p.Particle.mass with | "ZERO" -> Coupling.Vector | _ -> Coupling.Massive_Vector end | s -> s type state = { directory : string; model : t } let initialized = ref None let is_initialized_from dir = match !initialized with | None -> false | Some state -> dir = state.directory let dump_raw = ref false let init dir = let model = filter_unphysical (parse_directory dir) in if !dump_raw then dump model; let tables = Lookup.of_model model in let vertices () = translate_vertices model tables in let particle f = tables.Lookup.particle f in let lorentz f = lorentz_of_particle (particle f) in let propagator f = let p = particle f in match p.Particle.propagator with | None -> propagator_of_lorentz (lorentz_of_particle p) | Some s -> Coupling.Prop_UFO s in let gauge_symbol () = "?GAUGE?" in let constant_symbol s = s in let parameters = translate_parameters model in M.setup ~color:(fun f -> UFOx.Color.omega (particle f).Particle.color) ~nc:(fun () -> model.nc) ~pdg:(fun f -> (particle f).Particle.pdg_code) ~lorentz ~propagator ~width:(fun f -> Coupling.Constant) ~goldstone:(fun f -> None) ~conjugate:tables.Lookup.conjugate ~fermion:(fun f -> fermion_of_lorentz (lorentz f)) ~vertices ~flavors:[("All Flavors", tables.Lookup.flavors)] ~parameters:(fun () -> parameters) ~flavor_of_string:tables.Lookup.flavor_of_string ~flavor_to_string:(fun f -> (particle f).Particle.name) ~flavor_to_TeX:(fun f -> (particle f).Particle.texname) ~flavor_symbol:tables.Lookup.flavor_symbol ~gauge_symbol ~mass_symbol:(fun f -> (particle f).Particle.mass) ~width_symbol:(fun f -> (particle f).Particle.width) ~constant_symbol; initialized := Some { directory = dir; model = model } let ufo_directory = ref Config.default_UFO_dir let load () = if is_initialized_from !ufo_directory then () else init !ufo_directory let include_all_fusions = ref false (* In case of Majorana spinors, also generate all combinations of charge conjugated fermion lines. The naming convention is to append \texttt{\_c}$nm$ if the $\gamma$-matrices of the fermion line $n\to m$ has been charge conjugated (this could become impractical for too many fermions at a vertex, but shouldn't matter in real life). *) (* Here we alway generate \emph{all} charge conjugations, because we treat \emph{all} fermions as Majorana fermion, if there is at least one Majorana fermion in the model! *) let is_majorana = function | Coupling.Majorana | Coupling.Vectorspinor | Coupling.Maj_Ghost -> true | _ -> false let name_spins_structure spins l = (l.Lorentz.name, spins, l.Lorentz.structure) let fusions_of_model ?only model = let include_fusion = match !include_all_fusions, only with | true, _ | false, None -> (fun name -> true) | false, Some names -> (fun name -> SSet.mem name names) in SMap.fold (fun name l acc -> if include_fusion name then List.fold_left (fun acc p -> let l' = Lorentz.permute p l in match l'.Lorentz.spins with | Lorentz.Unused -> acc | Lorentz.Unique spins -> if Array.exists is_majorana spins then List.map (name_spins_structure spins) (Lorentz.required_charge_conjugates l') @ acc else name_spins_structure spins l' :: acc | Lorentz.Ambiguous _ -> failwith "fusions: Lorentz.Ambiguous") [] (Permutation.Default.cyclic l.Lorentz.n) @ acc else acc) model.lorentz [] let fusions ?only () = match !initialized with | None -> [] | Some { model = model } -> fusions_of_model ?only model let propagators_of_model ?only model = let include_propagator = match !include_all_fusions, only with | true, _ | false, None -> (fun name -> true) | false, Some names -> (fun name -> SSet.mem name names) in SMap.fold (fun name p acc -> if include_propagator name then (name, p) :: acc else acc) model.propagators [] let propagators ?only () = match !initialized with | None -> [] | Some { model = model } -> propagators_of_model ?only model let include_hadrons = ref true let ufo_majorana_warnings = [ "***************************************************"; "* *"; "* CAVEAT: *"; "* *"; "* These amplitudes have been computed for a *"; "* UFO model containing Majorana fermions. *"; "* This version of O'Mega contains some known *"; "* bugs for this case. It was released early at *"; "* the request of the Linear Collider community. *"; "* *"; "* These amplitudes MUST NOT be used for *"; "* publications without prior consulation *"; "* with the WHIZARD authors !!! *"; "* *"; "***************************************************" ] let caveats () = if !use_majorana_spinors then ufo_majorana_warnings else [] module Whizard : sig val write : unit -> unit end = struct let write_header dir = Printf.printf "# WHIZARD Model file derived from UFO directory\n"; Printf.printf "# '%s'\n\n" dir; List.iter (fun s -> Printf.printf "# %s\n" s) (M.caveats ()); Printf.printf "model \"%s\"\n\n" (Filename.basename dir) let write_input_parameters parameters = let open Parameter in Printf.printf "# Independent (input) Parameters\n"; List.iter (fun p -> Printf.printf "parameter %s = %s" p.name (value_to_numeric p.value); begin match p.lhablock, p.lhacode with | None, None -> () | Some name, Some (index :: indices) -> Printf.printf " slha_entry %s %d" name index; List.iter (fun i -> Printf.printf " %d" i) indices | Some name, None -> Printf.eprintf "UFO: parameter %s: slhablock %s without slhacode\n" p.name name | Some name, Some [] -> Printf.eprintf "UFO: parameter %s: slhablock %s with empty slhacode\n" p.name name | None, Some _ -> Printf.eprintf "UFO: parameter %s: slhacode without slhablock\n" p.name end; Printf.printf "\n") parameters; Printf.printf "\n" let write_derived_parameters parameters = let open Parameter in Printf.printf "# Dependent (derived) Parameters\n"; List.iter (fun p -> Printf.printf "derived %s = %s\n" p.name (value_to_expr alpha_s_half p.value)) parameters let write_particles particles = let open Particle in Printf.printf "# Particles\n"; Printf.printf "# NB: hypercharge assignments appear to be unreliable\n"; Printf.printf "# therefore we can't infer the isospin\n"; Printf.printf "# NB: parton-, gauge- & handedness are unavailable\n"; List.iter (fun p -> if not p.is_anti then begin Printf.printf "particle \"%s\" %d ### parton? gauge? left?\n" p.name p.pdg_code; Printf.printf " spin %s charge %s color %s ### isospin?\n" (UFOx.Lorentz.rep_to_string_whizard p.spin) (charge_to_string p.charge) (UFOx.Color.rep_to_string_whizard p.color); Printf.printf " name \"%s\"\n" p.name; if p.antiname <> p.name then Printf.printf " anti \"%s\"\n" p.antiname; Printf.printf " tex_name \"%s\"\n" p.texname; if p.antiname <> p.name then Printf.printf " tex_anti \"%s\"\n" p.antitexname; Printf.printf " mass %s width %s\n\n" p.mass p.width end) (values particles); Printf.printf "\n" let write_hadrons () = Printf.printf "# Hadrons (protons and beam remnants)\n"; Printf.printf "# NB: these are NOT part of the UFO model\n"; Printf.printf "# but added for WHIZARD's convenience!\n"; Printf.printf "particle PROTON 2212\n"; Printf.printf " spin 1/2 charge 1\n"; Printf.printf " name p \"p+\"\n"; Printf.printf " anti pbar \"p-\"\n"; Printf.printf "particle HADRON_REMNANT 90\n"; Printf.printf " name hr\n"; Printf.printf " tex_name \"had_r\"\n"; Printf.printf "particle HADRON_REMNANT_SINGLET 91\n"; Printf.printf " name hr1\n"; Printf.printf " tex_name \"had_r^{(1)}\"\n"; Printf.printf "particle HADRON_REMNANT_TRIPLET 92\n"; Printf.printf " color 3\n"; Printf.printf " name hr3\n"; Printf.printf " tex_name \"had_r^{(3)}\"\n"; Printf.printf " anti hr3bar\n"; Printf.printf " tex_anti \"had_r^{(\\bar 3)}\"\n"; Printf.printf "particle HADRON_REMNANT_OCTET 93\n"; Printf.printf " color 8\n"; Printf.printf " name hr8\n"; Printf.printf " tex_name \"had_r^{(8)}\"\n"; Printf.printf "\n" let vertex_to_string model v = String.concat " " (List.map (fun s -> "\"" ^ (SMap.find s model.particles).Particle.name ^ "\"") (Array.to_list v.Vertex.particles)) let write_vertices3 model vertices = Printf.printf "# Vertices (for phasespace generation only)\n"; Printf.printf "# NB: particles should be sorted increasing in mass.\n"; Printf.printf "# This is NOT implemented yet!\n"; List.iter (fun v -> if Array.length v.Vertex.particles = 3 then Printf.printf "vertex %s\n" (vertex_to_string model v)) (values vertices); Printf.printf "\n" let write_vertices_higher model vertices = Printf.printf "# Higher Order Vertices (ignored by phasespace generation)\n"; List.iter (fun v -> if Array.length v.Vertex.particles <> 3 then Printf.printf "# vertex %s\n" (vertex_to_string model v)) (values vertices); Printf.printf "\n" let write_vertices model vertices = write_vertices3 model vertices; write_vertices_higher model vertices let write () = match !initialized with | None -> failwith "UFO.Whizard.write: UFO model not initialized" | Some { directory = dir; model = model } -> let input_parameters, derived_parameters = classify_parameters model in write_header dir; write_input_parameters input_parameters; write_derived_parameters derived_parameters; write_particles model.particles; if !include_hadrons then write_hadrons (); write_vertices model model.vertices; exit 0 end let options = Options.create [ ("UFO_dir", Arg.String (fun name -> ufo_directory := name), "UFO model directory (default: " ^ !ufo_directory ^ ")"); ("Majorana", Arg.Set use_majorana_spinors, "use Majorana spinors (must come _before_ exec!)"); ("divide_propagators_by_i", Arg.Set divide_propagators_by_i, "divide propagators by I (pre 2013 FeynRules convention)"); ("verbatim_Hg", Arg.Set verbatim_higgs_glue, "don't correct the color flows for effective Higgs Gluon couplings"); ("write_WHIZARD", Arg.Unit Whizard.write, "write the WHIZARD model file (required once per model)"); ("long_flavors", Arg.Unit (fun () -> Lookup.flavor_format := Lookup.Long), "write use the UFO flavor names instead of integers"); ("dump", Arg.Set dump_raw, "dump UFO model for debugging the parser (must come _before_ exec!)"); ("all_fusions", Arg.Set include_all_fusions, "include all fusions in the fortran module"); ("no_hadrons", Arg.Clear include_hadrons, "don't add any particle not in the UFO file"); ("add_hadrons", Arg.Set include_hadrons, "add protons and beam remants for WHIZARD"); ("exec", Arg.Unit load, "load the UFO model files (required _before_ using particles names)"); ("help", Arg.Unit (fun () -> prerr_endline "..."), "print information on the model")] end module type Fortran_Target = sig val fuse : Algebra.QC.t -> string -> Coupling.lorentzn -> Coupling.fermion_lines -> string -> string list -> string list -> Coupling.fusen -> unit val lorentz_module : ?only:SSet.t -> ?name:string -> ?fortran_module:string -> ?parameter_module:string -> Format_Fortran.formatter -> unit -> unit end module Targets = struct module Fortran : Fortran_Target = struct open Format_Fortran let fuse = UFO_targets.Fortran.fuse let lorentz_functions ff fusions () = List.iter (fun (name, s, l) -> UFO_targets.Fortran.lorentz ff name s l) fusions let propagator_functions ff parameter_module propagators () = List.iter (fun (name, p) -> UFO_targets.Fortran.propagator ff name parameter_module p.Propagator.variables p.Propagator.spins p.Propagator.numerator p.Propagator.denominator) propagators let lorentz_module ?only ?(name="omega_amplitude_ufo") ?(fortran_module="omega95") ?(parameter_module="parameter_module") ff () = let printf fmt = fprintf ff fmt and nl = pp_newline ff in printf "module %s" name; nl (); printf " use kinds"; nl (); printf " use %s" fortran_module; nl (); printf " implicit none"; nl (); printf " private"; nl (); let fusions = Model.fusions ?only () and propagators = Model.propagators () in List.iter (fun (name, _, _) -> printf " public :: %s" name; nl ()) fusions; List.iter (fun (name, _) -> printf " public :: pr_U_%s" name; nl ()) propagators; UFO_targets.Fortran.eps4_g4_g44_decl ff (); UFO_targets.Fortran.eps4_g4_g44_init ff (); printf "contains"; nl (); UFO_targets.Fortran.inner_product_functions ff (); lorentz_functions ff fusions (); propagator_functions ff parameter_module propagators (); printf "end module %s" name; nl (); pp_flush ff () end end module type Test = sig val suite : OUnit.test end module Test : Test = struct open OUnit let lexer s = UFO_lexer.token (UFO_lexer.init_position "" (Lexing.from_string s)) let suite_lexer_escapes = "escapes" >::: [ "single-quote" >:: (fun () -> assert_equal (UFO_parser.STRING "a'b'c") (lexer "'a\\'b\\'c'")); "unterminated" >:: (fun () -> assert_raises End_of_file (fun () -> lexer "'a\\'b\\'c")) ] let suite_lexer = "lexer" >::: [suite_lexer_escapes] let suite = "UFO" >::: [suite_lexer] end Index: trunk/omega/src/young.mli =================================================================== --- trunk/omega/src/young.mli (revision 8848) +++ trunk/omega/src/young.mli (revision 8849) @@ -1,145 +1,147 @@ (* young.mli -- Copyright (C) 2022- by Wolfgang Kilian Thorsten Ohl Juergen Reuter WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) (* Caveat: the following are not optimized for large Young diagrams and tableaux. They are straightforward implementations of the definitions, since we are unlikely to meet large diagrams. To make matters worse, native integer arithmetic will overflow already for diagrams with more than 20 cells. Since the [Num] library has been removed from the O'Caml distribution with version 4.06, we can not use it as a shortcut. Requiring Whizard/O'Mega users to install [Num] or its successor [Zarith] is probably not worth the effort. *) (* \ytableausetup{centertableaux,smalltableaux} *) (* \thocwmodulesection{Young Diagrams} *) (* Young diagrams can be represented by a non-increasing list of positive integers, corresponding to the number of boxes in each row: \begin{equation} \ydiagram{5,4,4,2} \Longleftrightarrow \lbrack 5;4;4;2 \rbrack \end{equation} *) type diagram = int list (* Check that the diagram is valid, i.\,e.~the number of boxes is non-increasing from top to bottom. *) val valid_diagram : diagram -> bool (* Count the number of cells. *) val num_cells_diagram : diagram -> int -(* Transpose a diagram: +(* Conjugate a diagram: \begin{equation} \ydiagram{5,4,4,2} \mapsto \ydiagram{4,4,3,3,1} \end{equation} *) -val transpose_diagram : diagram -> diagram +val conjugate_diagram : diagram -> diagram (* The product of all the ``hook lengths'' in the diagram, e.\,g. \begin{equation} \ydiagram{5,4,4,2} \mapsto \ytableaushort{87541,6532,5421,21} \mapsto 8 \cdot 7 \cdot 6 \cdot 5^3 \cdot 4^2 \cdot 3 \cdot 2^3 = 16128000 \end{equation} where the intermediate step is only for illustration and does not represent a Young tableau! *) val hook_lengths_product : diagram -> int -(* Dimension of the representation of~$S_n$ described by the diagram +(* Number of standard tableaux corresponding to the diagram. + Also, the dimension of the representation of~$S_n$ described + by this diagram \begin{equation} d = \frac{n!}{\prod_{i=1}^n h_i} \end{equation} with~$n$ the number of cells and~$h_i$ the hook length of the $i$th cell. *) -val dim_rep_Sn : diagram -> int +val num_standard_tableaux : diagram -> int (* Normalization of the projector on the representation of $\mathrm{GL(N)}$ described by the diagram \begin{equation} \alpha = \frac{\prod_{R} |R|!\prod_{C} |C|!}{\prod_{i=1}^n h_i} \end{equation} with~$|R|$ and~$|C|$ the lengths of the row~$R$ and column~$C$, respectively. Returned as a pair of numerator and denominator, because it is not guaranteed to be integer. *) val normalization : diagram -> int * int (* \thocwmodulesection{Young Tableaux} *) (* There is an obvious representation as a list of lists: \begin{equation} \ytableaushort{023,14} \Longleftrightarrow \lbrack \lbrack 0; 2; 3 \rbrack; \lbrack 1; 4 \rbrack \rbrack \end{equation} *) type 'a tableau = 'a list list (* Ignoring the contents of the cells of a Young tableau produces a unique corresponding Young diagram. \begin{equation} \ytableaushort{023,14} \mapsto \ydiagram{3,2} \end{equation} *) val diagram_of_tableau : 'a tableau -> diagram (* The number of columns must be non-increasing. Obviously, [valid_tableau] is the composition of [diagram_of_tableau] and [valid_diagram].*) val valid_tableau : 'a tableau -> bool (* A tableau is called \textit{semistandard}, iff the entries don't increase along rows and strictly increase along columns. - Therefore, the transpose of a semistandard tableau is \emph{not} + Therefore, the conjugate of a semistandard tableau is \emph{not} necessarily semistandard. *) val semistandard_tableau : 'a tableau -> bool (* A tableau is called \textit{standard}, iff it is semistandard and the entries are an uninterrupted sequence of natural numbers. If the optional [offset] is specified, it must match the smallest of these numbers. Some authors expect [offset=1], but we want to be able to start from 0 as well. - The transpose of a standard tableau is again a standard tableau. *) + The conjugate of a standard tableau is again a standard tableau. *) val standard_tableau : ?offset:int -> int tableau -> bool (* The contents of the cells and their number. *) val cells_tableau : 'a tableau -> 'a list val num_cells_tableau : 'a tableau -> int -(* Transpose a Young tableau +(* Conjugate a Young tableau \begin{equation} \ytableaushort{023,14} \mapsto \ytableaushort{01,24,3} \end{equation} *) -val transpose_tableau : 'a tableau -> 'a tableau +val conjugate_tableau : 'a tableau -> 'a tableau (* \thocwmodulesection{Unit Tests} *) module type Test = sig val suite : OUnit.test val suite_long : OUnit.test end module Test : Test Index: trunk/omega/src/young.ml =================================================================== --- trunk/omega/src/young.ml (revision 8848) +++ trunk/omega/src/young.ml (revision 8849) @@ -1,276 +1,276 @@ (* young.ml -- Copyright (C) 2022- by Wolfgang Kilian Thorsten Ohl Juergen Reuter WHIZARD is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. WHIZARD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *) (* Avoid refering to [Pervasives.compare], because [Pervasives] will become [Stdlib.Pervasives] in O'Caml 4.07 and [Stdlib] in O'Caml 4.08. *) let pcompare = compare type diagram = int list type 'a tableau = 'a list list (* Not exposed. Just for documentation. *) type 'a table = 'a option array array (* The following three are candidates for [ThoList]. *) let rec sum = function | [] -> 0 | n :: rest -> n + sum rest let rec product = function | [] -> 1 | n :: rest -> n * product rest (* Test a predicate for each pair of consecutive elements of a list. Trivially true for empty and one-element lists. *) let rec for_all_pairs predicate = function | [] | [_] -> true | a1 :: (a2 :: _ as a_list) -> if not (predicate a1 a2) then false else for_all_pairs predicate a_list let decreasing l = for_all_pairs (fun a1 a2 -> pcompare a1 a2 > 0) l let increasing l = for_all_pairs (fun a1 a2 -> pcompare a1 a2 < 0) l let non_increasing l = for_all_pairs (fun a1 a2 -> pcompare a1 a2 >= 0) l let non_decreasing l = for_all_pairs (fun a1 a2 -> pcompare a1 a2 <= 0) l let valid_diagram = non_increasing let diagram_rows d = List.length d let diagram_columns = function | [] -> 0 | nc :: _ -> nc let take_column d = let rec take_column' len acc = function | [] -> (len, List.rev acc) | cols :: rest -> if cols <= 1 then take_column' (succ len) acc rest else take_column' (succ len) (pred cols :: acc) rest in take_column' 0 [] d -let transpose_diagram_new d = - let rec transpose_diagram' rows = +let conjugate_diagram_new d = + let rec conjugate_diagram' rows = match take_column rows with | n, [] -> [n] - | n, rest -> n :: transpose_diagram' rest in - transpose_diagram' d + | n, rest -> n :: conjugate_diagram' rest in + conjugate_diagram' d let tableau_rows t = List.length t let tableau_columns = function | [] -> 0 | row :: _ -> List.length row let num_cells_diagram d = sum d let cells_tableau t = List.flatten t let num_cells_tableau t = List.fold_left (fun acc row -> acc + List.length row) 0 t let diagram_of_tableau t = List.map List.length t let tableau_of_diagram cell d = List.map (ThoList.clone cell) d (* Note that the first index counts the rows and the second the columns! *) let array_of_tableau t = let nr = tableau_rows t and nc = tableau_columns t in let a = Array.make_matrix nr nc None in List.iteri (fun ir -> List.iteri (fun ic cell -> a.(ir).(ic) <- Some cell)) t; a let transpose_array a = let nr = Array.length a in if nr <= 0 then invalid_arg "Young.transpose_array" else let nc = Array.length a.(0) in let a' = Array.make_matrix nc nr None in for ic = 0 to pred nc do for ir = 0 to pred nr do a'.(ic).(ir) <- a.(ir).(ic) done done; a' let list_of_array_row a = let n = Array.length a in let rec list_of_array_row' ic = if ic >= n then [] else match a.(ic) with | None -> [] | Some cell -> cell :: list_of_array_row' (succ ic) in list_of_array_row' 0 let tableau_of_array a = Array.fold_right (fun row acc -> list_of_array_row row :: acc) a [] -let transpose_tableau t = +let conjugate_tableau t = array_of_tableau t |> transpose_array |> tableau_of_array -let transpose_diagram d = - tableau_of_diagram () d |> transpose_tableau |> diagram_of_tableau +let conjugate_diagram d = + tableau_of_diagram () d |> conjugate_tableau |> diagram_of_tableau let valid_tableau t = valid_diagram (diagram_of_tableau t) let semistandard_tableau t = let rows = t - and columns = transpose_tableau t in + and columns = conjugate_tableau t in valid_tableau t && List.for_all non_decreasing rows && List.for_all increasing columns let standard_tableau ?offset t = match List.sort pcompare (cells_tableau t) with | [] -> true | cell :: _ as cell_list -> (match offset with None -> true | Some o -> cell = o) && for_all_pairs (fun c1 c2 -> c2 = c1 + 1) cell_list && semistandard_tableau t let hook_lengths_table d = let nr = diagram_rows d and nc = diagram_columns d in if min nr nc <= 0 then invalid_arg "Young.hook_lengths_table" else let a = array_of_tableau (tableau_of_diagram 0 d) in let cols = Array.of_list d and rows = transpose_array a |> tableau_of_array |> diagram_of_tableau |> Array.of_list in for ir = 0 to pred nr do for ic = 0 to pred cols.(ir) do a.(ir).(ic) <- Some (rows.(ic) - ir + cols.(ir) - ic - 1) done done; a (* \begin{dubious} The following products and factorials can easily overflow, even if the final ratio is a smallish number. We can avoid this by representing them as lists of factors (or maps from factors to powers). The ratio can be computed by first cancelling all common factors and multiplying the remaining factors at the very end. \end{dubious} *) let hook_lengths_product d = let nr = diagram_rows d and nc = diagram_columns d in if min nr nc <= 0 then 0 else let cols = Array.of_list d - and rows = Array.of_list (transpose_diagram d) in + and rows = Array.of_list (conjugate_diagram d) in let n = ref 1 in for ir = 0 to pred nr do for ic = 0 to pred cols.(ir) do n := !n * (rows.(ic) - ir + cols.(ir) - ic - 1) done done; !n -let dim_rep_Sn d = +let num_standard_tableaux d = let num = Combinatorics.factorial (num_cells_diagram d) and den = hook_lengths_product d in if num mod den <> 0 then - failwith "Young.dim_rep_Sn" + failwith "Young.num_standard_tableaux" else num / den -(* Note that [hook_lengths_product] calls [transpose_diagram] +(* Note that [hook_lengths_product] calls [conjugate_diagram] and this calls it again. This is wasteful, but probably no big deal for our applications. *) let normalization d = let num = - product (List.map Combinatorics.factorial (d @ transpose_diagram d)) + product (List.map Combinatorics.factorial (d @ conjugate_diagram d)) and den = hook_lengths_product d in (num, den) module type Test = sig val suite : OUnit.test val suite_long : OUnit.test end module Test = struct open OUnit let random_int ratio = truncate (Random.float ratio +. 0.5) let random_diagram ?(ratio=1.0) rows = let rec random_diagram' acc row cols = if row >= rows then acc else let cols' = cols + random_int ratio in random_diagram' (cols' :: acc) (succ row) cols' in random_diagram' [] 0 (1 + random_int ratio) let suite_hook_lengths_product = "hook_lengths_product" >::: [ "[4;3;2]" >:: (fun () -> assert_equal 2160 (hook_lengths_product [4; 3; 2])) ] - let suite_dim_rep_Sn = - "dim_rep_Sn" >::: + let suite_num_standard_tableaux = + "num_standard_tableaux" >::: [ "[4;3;2]" >:: - (fun () -> assert_equal 168 (dim_rep_Sn [4; 3; 2])) ] + (fun () -> assert_equal 168 (num_standard_tableaux [4; 3; 2])) ] let suite_normalization = "normalization" >::: [ "[2;1]" >:: (fun () -> assert_equal (4, 3) (normalization [2; 1])) ] let suite = "Young" >::: [suite_hook_lengths_product; - suite_dim_rep_Sn; + suite_num_standard_tableaux; suite_normalization] let suite_long = "Young long" >::: [] end