Index: trunk/omega/src/fusion.mli =================================================================== --- trunk/omega/src/fusion.mli (revision 8498) +++ trunk/omega/src/fusion.mli (revision 8499) @@ -1,381 +1,381 @@ (* fusion.mli -- Copyright (C) 1999-2021 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 T = sig val options : Options.t (* JRR's implementation of Majoranas needs a special case. *) val vintage : bool (* Wavefunctions are an abstract data type, containing a momentum~[p] and additional quantum numbers, collected in~[flavor]. *) type wf val conjugate : wf -> wf (* Obviously, [flavor] is not restricted to the physical notion of flavor, but can carry spin, color, etc. *) type flavor val flavor : wf -> flavor type flavor_sans_color val flavor_sans_color : wf -> flavor_sans_color (* Momenta are represented by an abstract datatype (defined in~[Momentum]) that is optimized for performance. They can be accessed either abstractly or as lists of indices of the external momenta. These indices are assigned sequentially by [amplitude] below. *) type p val momentum : wf -> p val momentum_list : wf -> int list (* At tree level, the wave functions are uniquely specified by [flavor] and momentum. If loops are included, we need to distinguish among orders. Also, if we build a result from an incomplete sum of diagrams, we need to add a distinguishing mark. At the moment, we assume that a [string] that can be attached to the symbol suffices. *) val wf_tag : wf -> string option (* Coupling constants *) type constant (* and right hand sides of assignments. The latter are formed from a sign from Fermi statistics, a coupling (constand and Lorentz structure) and wave functions. *) type coupling type rhs type 'a children val sign : rhs -> int val coupling : rhs -> constant Coupling.t val coupling_tag : rhs -> string option type exclusions val no_exclusions : exclusions (* In renormalized perturbation theory, couplings come in different orders of the loop expansion. Be prepared: [val order : rhs -> int] *) (* \begin{dubious} This is here only for the benefit of [Target] and shall become [val children : rhs -> wf children] later \ldots \end{dubious} *) val children : rhs -> wf list (* Fusions come in two types: fusions of wave functions to off-shell wave functions: \begin{equation*} \phi(p+q) = \phi(p)\phi(q) \end{equation*} *) type fusion val lhs : fusion -> wf val rhs : fusion -> rhs list (* and products at the keystones: \begin{equation*} \phi(-p-q)\cdot\phi(p)\phi(q) \end{equation*} *) type braket val bra : braket -> wf val ket : braket -> rhs list (* [amplitude goldstones incoming outgoing] calculates the amplitude for scattering of [incoming] to [outgoing]. If [goldstones] is true, also non-propagating off-shell Goldstone amplitudes are included to allow the checking of Slavnov-Taylor identities. *) type amplitude type amplitude_sans_color type selectors val amplitudes : bool -> exclusions -> selectors -> flavor_sans_color list -> flavor_sans_color list -> amplitude list val amplitude_sans_color : bool -> exclusions -> selectors -> flavor_sans_color list -> flavor_sans_color list -> amplitude_sans_color val dependencies : amplitude -> wf -> (wf, coupling) Tree2.t (* We should be precise regarding the semantics of the following functions, since modules implementating [Target] must not make any mistakes interpreting the return values. Instead of calculating the amplitude \begin{subequations} \begin{equation} \label{eq:physical-amplitude} \Braket{f_3,p_3,f_4,p_4,\ldots|T|f_1,p_1,f_2,p_2} \end{equation} directly, O'Mega calculates the---equivalent, but more symmetrical---crossed amplitude \begin{equation} \Braket{\bar f_1,-p_1,\bar f_2,-p_2,f_3,p_3,f_4,p_4,\ldots|T|0} \end{equation} Internally, all flavors are represented by their charge conjugates \begin{equation} \label{eq:internal-amplitude} A(f_1,-p_1,f_2,-p_2,\bar f_3,p_3,\bar f_4,p_4,\ldots) \end{equation} \end{subequations} The correspondence of vertex and term in the lagrangian \begin{equation} \parbox{26\unitlength}{% \fmfframe(5,3)(5,3){% \begin{fmfgraph*}(15,20) \fmfleft{v} \fmfright{p,A,e} \fmflabel{$\mathrm{e}^-$}{e} \fmflabel{$\mathrm{e}^+$}{p} \fmflabel{$\mathrm{A}$}{A} \fmf{fermion}{p,v,e} \fmf{photon}{A,v} \fmfdot{v} \end{fmfgraph*}}}: \bar\psi\fmslash{A}\psi \end{equation} suggests to denote the \emph{outgoing} particle by the flavor of the \emph{anti}particle and the \emph{outgoing} \emph{anti}particle by the flavor of the particle, since this choice allows to represent the vertex by a triple \begin{equation} \bar\psi\fmslash{A}\psi: (\mathrm{e}^+,A,\mathrm{e}^-) \end{equation} which is more intuitive than the alternative $(\mathrm{e}^-,A,\mathrm{e}^+)$. Also, when thinking in terms of building wavefunctions from the outside in, the outgoing \emph{antiparticle} is represented by a \emph{particle} propagator and vice versa\footnote{Even if this choice will appear slightly counter-intuitive on the [Target] side, one must keep in mind that much more people are expected to prepare [Model]s.}. [incoming] and [outgoing] are the physical flavors as in~(\ref{eq:physical-amplitude}) *) val incoming : amplitude -> flavor list val outgoing : amplitude -> flavor list (* [externals] are flavors and momenta as in~(\ref{eq:internal-amplitude}) *) val externals : amplitude -> wf list val variables : amplitude -> wf list val fusions : amplitude -> fusion list val brakets : amplitude -> braket list val on_shell : amplitude -> (wf -> bool) val is_gauss : amplitude -> (wf -> bool) val constraints : amplitude -> string option val symmetry : amplitude -> int val allowed : amplitude -> bool (* \thocwmodulesubsection{Performance Hacks} *) val initialize_cache : string -> unit val set_cache_name : string -> unit (* \thocwmodulesubsection{Diagnostics} *) val check_charges : unit -> flavor_sans_color list list val count_fusions : amplitude -> int val count_propagators : amplitude -> int val count_diagrams : amplitude -> int val forest : wf -> amplitude -> ((wf * coupling option, wf) Tree.t) list val poles : amplitude -> wf list list val s_channel : amplitude -> wf list val tower_to_dot : out_channel -> amplitude -> unit val amplitude_to_dot : out_channel -> amplitude -> unit (* \thocwmodulesubsection{WHIZARD} *) val phase_space_channels : out_channel -> amplitude_sans_color -> unit val phase_space_channels_flipped : out_channel -> amplitude_sans_color -> unit end (* There is more than one way to make fusions. *) module type Maker = functor (P : Momentum.T) -> functor (M : Model.T) -> T with type p = P.t and type flavor = Colorize.It(M).flavor and type flavor_sans_color = M.flavor and type constant = M.constant and type selectors = Cascade.Make(M)(P).selectors (*i If we want or need to expose [Make], here's how to do it: module type Stat = sig type flavor type stat exception Impossible val stat : flavor -> int -> stat val stat_fuse : stat -> stat -> flavor -> stat val stat_sign : stat -> int end module type Stat_Maker = functor (M : Model.T) -> Stat with type flavor = M.flavor module Make : functor (PT : Tuple.Poly) (Stat : Stat_Maker) (T : Topology.T with type 'a children = 'a PT.t) -> Maker i*) (* Straightforward Dirac fermions vs. slightly more complicated Majorana fermions: *) exception Majorana module Binary : Maker module Binary_Majorana : Maker module Mixed23 : Maker module Mixed23_Majorana : Maker module Nary : functor (B : Tuple.Bound) -> Maker module Nary_Majorana : functor (B : Tuple.Bound) -> Maker (* We can also proceed \'a la~\cite{HELAC:2000}. Empirically, this will use slightly~($O(10\%)$) fewer fusions than the symmetric factorization. Our implementation uses significantly~($O(50\%)$) fewer fusions than reported by~\cite{HELAC:2000}. Our pruning of the DAG might be responsible for this. *) module Helac : functor (B : Tuple.Bound) -> Maker -(* [module Helac_Majorana : functor (B : Tuple.Bound) -> Maker] *) +module Helac_Majorana : functor (B : Tuple.Bound) -> Maker (* \thocwmodulesection{Multiple Amplitudes} *) module type Multi = sig exception Mismatch val options : Options.t type flavor type process = flavor list * flavor list type amplitude type fusion type wf type exclusions val no_exclusions : exclusions type selectors type amplitudes (* Construct all possible color flow amplitudes for a given process. *) val amplitudes : bool -> int option -> exclusions -> selectors -> process list -> amplitudes val empty : amplitudes (* Precompute the vertex table cache. *) val initialize_cache : string -> unit val set_cache_name : string -> unit (* The list of all combinations of incoming and outgoing particles with a nonvanishing scattering amplitude. *) val flavors : amplitudes -> process list (* The list of all combinations of incoming and outgoing particles that don't lead to any color flow with non vanishing scattering amplitude. *) val vanishing_flavors : amplitudes -> process list (* The list of all color flows with a nonvanishing scattering amplitude. *) val color_flows : amplitudes -> Color.Flow.t list (* The list of all valid helicity combinations. *) val helicities : amplitudes -> (int list * int list) list (* The list of all amplitudes. *) val processes : amplitudes -> amplitude list (* [(process_table a).(f).(c)] returns the amplitude for the [f]th allowed flavor combination and the [c]th allowed color flow as an [amplitude option]. *) val process_table : amplitudes -> amplitude option array array (* The list of all non redundant fusions together with the amplitudes they came from. *) val fusions : amplitudes -> (fusion * amplitude) list (* If there's more than external flavor state, the wavefunctions are \emph{not} uniquely specified by [flavor] and [Momentum.t]. This function can be used to determine how many variables must be allocated. *) val multiplicity : amplitudes -> wf -> int (* This function can be used to disambiguate wavefunctions with the same combination of [flavor] and [Momentum.t]. *) val dictionary : amplitudes -> amplitude -> wf -> int (* [(color_factors a).(c1).(c2)] power of~$N_C$ for the given product of color flows. *) val color_factors : amplitudes -> Color.Flow.factor array array (* A description of optional diagram selectors. *) val constraints : amplitudes -> string option end module type Multi_Maker = functor (Fusion_Maker : Maker) -> functor (P : Momentum.T) -> functor (M : Model.T) -> Multi with type flavor = M.flavor and type amplitude = Fusion_Maker(P)(M).amplitude and type fusion = Fusion_Maker(P)(M).fusion and type wf = Fusion_Maker(P)(M).wf and type selectors = Fusion_Maker(P)(M).selectors module Multi : Multi_Maker (* \thocwmodulesection{Tags} *) (* It appears that there are useful applications for tagging couplings and wave functions, e.\,g.~skeleton expansion and diagram selections. We can abstract this in a [Tags] signature: *) module type Tags = sig type wf type coupling type 'a children val null_wf : wf val null_coupling : coupling val fuse : coupling -> wf children -> wf val wf_to_string : wf -> string option val coupling_to_string : coupling -> string option end module type Tagger = functor (PT : Tuple.Poly) -> Tags with type 'a children = 'a PT.t module type Tagged_Maker = functor (Tagger : Tagger) -> functor (P : Momentum.T) -> functor (M : Model.T) -> T with type p = P.t and type flavor = Colorize.It(M).flavor and type flavor_sans_color = M.flavor and type constant = M.constant module Tagged_Binary : Tagged_Maker (*i * Local Variables: * mode:caml * indent-tabs-mode:nil * page-delimiter:"^(\\* .*\n" * End: i*) Index: trunk/omega/src/omega.ml =================================================================== --- trunk/omega/src/omega.ml (revision 8498) +++ trunk/omega/src/omega.ml (revision 8499) @@ -1,694 +1,699 @@ (* omega.ml -- Copyright (C) 1999-2021 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. *) let (<<) f g x = f (g x) let (>>) f g x = g (f x) module P = Momentum.Default module P_Whizard = Momentum.DefaultW module type T = sig val main : unit -> unit type flavor val diagrams : flavor -> flavor -> flavor list -> ((flavor * Momentum.Default.t) * (flavor * Momentum.Default.t, flavor * Momentum.Default.t) Tree.t) list end module Make (Fusion_Maker : Fusion.Maker) (Target_Maker : Target.Maker) (M : Model.T) = struct module CM = Colorize.It(M) type flavor = M.flavor module Proc = Process.Make(M) (* \begin{dubious} We must have initialized the vertices \emph{before} applying [Fusion_Maker], at least if we want to continue using the vertex cache! \end{dubious} *) (* \begin{dubious} NB: this causes the constant initializers in [Fusion_Maker] more than once. Such side effects must be avoided if the initializers involve expensive computations. \emph{Relying on the fact that the functor will be called only once is not a good idea!} \end{dubious} *) module F = Fusion_Maker(P)(M) module CF = Fusion.Multi(Fusion_Maker)(P)(M) module T = Target_Maker(Fusion_Maker)(P)(M) module W = Whizard.Make(Fusion_Maker)(P)(P_Whizard)(M) module C = Cascade.Make(M)(P) module VSet = Set.Make (struct type t = F.constant Coupling.t let compare = compare end) (* For the phase space, we need asymmetric DAGs. - HACK: since we will not use this to compute amplitudes, there's + Since we will not use this to compute amplitudes, there's no need to supply the proper statistics module and we may - assume Dirac fermions. - - HACK: for the phase space, we should be able to work on the - uncolored model. *) + always use Majorana fermions to be as general as possible. + In principle, we could expose in [Fusion.T] the [Fusion.Stat_Maker] + used by [Fusion_Maker] to construct it, but that is just not + worth the effort. + + \begin{dubious} + For the phase space, we should be able to work on the + uncolored model. + \end{dubious} *) module PHS = - Fusion.Helac(struct let max_arity () = pred (M.max_degree ()) end)(P)(M) + Fusion.Helac_Majorana(struct let max_arity () = pred (M.max_degree ()) end)(P)(M) (* Form a ['a list] from a ['a option array], containing the elements that are not [None] in order. *) let opt_array_to_list a = let rec opt_array_to_list' acc i a = if i < 0 then acc else begin match a.(i) with | None -> opt_array_to_list' acc (pred i) a | Some x -> opt_array_to_list' (x :: acc) (pred i) a end in opt_array_to_list' [] (Array.length a - 1) a (* Return a list of [CF.amplitude list]s, corresponig to the diagrams for a specific color flow for each flavor combination. *) let amplitudes_by_flavor amplitudes = List.map opt_array_to_list (Array.to_list (CF.process_table amplitudes)) (* \begin{dubious} If we plan to distiguish different couplings later on, we can no long map all instances of [coupling option] in the tree to [None]. In this case, we will need to normalize different fusion orders [Coupling.fuse2], [Coupling.fuse3] or [Coupling.fusen], because they would otherwise lead to inequivalent diagrams. Unfortunately, this stuff packaged deep in [Fusion.Tagged_Coupling]. \end{dubious} *) (*i let strip_fuse' = function | Coupling.V3 (v, f, c) -> Coupling.V3 (v, Coupling.F12, c) | Coupling.V4 (v, f, c) -> Coupling.V4 (v, Coupling.F123, c) | Coupling.Vn (v, f, c) -> Coupling.Vn (v, [], c) let strip_fuse = function | Some c -> Some (strip_fuse' c) | None -> None i*) (* \begin{dubious} The [Tree.canonicalize] below should be necessary to remove topologically equivalent duplicates. \end{dubious} *) (* Take a [CF.amplitude list] assumed to correspond to the same external states after stripping the color and return a pair of the list of external particles and the corresponding Feynman diagrams without color. *) let wf1 amplitude = match F.externals amplitude with | wf :: _ -> wf | [] -> failwith "Omega.forest_sans_color: no external particles" let uniq l = ThoList.uniq (List.sort compare l) let forest_sans_color = function | amplitude :: _ as amplitudes -> let externals = F.externals amplitude in let prune_color wf = (F.flavor_sans_color wf, F.momentum_list wf) in let prune_color_and_couplings (wf, c) = (prune_color wf, None) in (List.map prune_color externals, uniq (List.map (fun t -> Tree.canonicalize (Tree.map prune_color_and_couplings prune_color t)) (ThoList.flatmap (fun a -> F.forest (wf1 a) a) amplitudes))) | [] -> ([], []) let dag_sans_color = function | amplitude :: _ as amplitudes -> let prune a = a in List.map prune amplitudes | [] -> [] let p2s p = if p >= 0 && p <= 9 then string_of_int p else if p <= 36 then String.make 1 (Char.chr (Char.code 'A' + p - 10)) else "_" let format_p wf = String.concat "" (List.map p2s (F.momentum_list wf)) let variable wf = M.flavor_to_string (F.flavor_sans_color wf) ^ "[" ^ format_p wf ^ "]" let variable' wf = CM.flavor_to_TeX (F.flavor wf) ^ "(" ^ format_p wf ^ ")" let feynmf_style propagator color = { Tree.style = begin match propagator with | Coupling.Prop_Feynman | Coupling.Prop_Gauge _ -> begin match color with | Color.AdjSUN _ -> Some ("gluon", "") | _ -> Some ("boson", "") end | Coupling.Prop_Col_Feynman -> Some ("gluon", "") | Coupling.Prop_Unitarity | Coupling.Prop_Rxi _ -> Some ("dbl_wiggly", "") | Coupling.Prop_Spinor | Coupling.Prop_ConjSpinor -> Some ("fermion", "") | _ -> None end; Tree.rev = begin match propagator with | Coupling.Prop_Spinor -> true | Coupling.Prop_ConjSpinor -> false | _ -> false end; Tree.label = None; Tree.tension = None } let header incoming outgoing = "$ " ^ String.concat " " (List.map (CM.flavor_to_TeX << F.flavor) incoming) ^ " \\to " ^ String.concat " " (List.map (CM.flavor_to_TeX << CM.conjugate << F.flavor) outgoing) ^ " $" let header_sans_color incoming outgoing = "$ " ^ String.concat " " (List.map (M.flavor_to_TeX << fst) incoming) ^ " \\to " ^ String.concat " " (List.map (M.flavor_to_TeX << M.conjugate << fst) outgoing) ^ " $" let diagram incoming tree = let fmf wf = let f = F.flavor wf in feynmf_style (CM.propagator f) (CM.color f) in Tree.map (fun (n, _) -> let n' = fmf n in if List.mem n incoming then { n' with Tree.rev = not n'.Tree.rev } else n') (fun l -> if List.mem l incoming then l else F.conjugate l) tree let diagram_sans_color incoming (tree) = let fmf (f, p) = feynmf_style (M.propagator f) (M.color f) in Tree.map (fun (n, c) -> let n' = fmf n in if List.mem n incoming then { n' with Tree.rev = not n'.Tree.rev } else n') (fun (f, p) -> if List.mem (f, p) incoming then (f, p) else (M.conjugate f, p)) tree let feynmf_set amplitude = match F.externals amplitude with | wf1 :: wf2 :: wfs -> let incoming = [wf1; wf2] in { Tree.header = header incoming wfs; Tree.incoming = incoming; Tree.diagrams = List.map (diagram incoming) (F.forest wf1 amplitude) } | _ -> failwith "less than two external particles" let feynmf_set_sans_color (externals, trees) = match externals with | wf1 :: wf2 :: wfs -> let incoming = [wf1; wf2] in { Tree.header = header_sans_color incoming wfs; Tree.incoming = incoming; Tree.diagrams = List.map (diagram_sans_color incoming) trees } | _ -> failwith "less than two external particles" let feynmf_set_sans_color_empty (externals, trees) = match externals with | wf1 :: wf2 :: wfs -> let incoming = [wf1; wf2] in { Tree.header = header_sans_color incoming wfs; Tree.incoming = incoming; Tree.diagrams = [] } | _ -> failwith "less than two external particles" let uncolored_colored amplitudes = { Tree.outer = feynmf_set_sans_color (forest_sans_color amplitudes); Tree.inner = List.map feynmf_set amplitudes } let uncolored_only amplitudes = { Tree.outer = feynmf_set_sans_color (forest_sans_color amplitudes); Tree.inner = [] } let colored_only amplitudes = { Tree.outer = feynmf_set_sans_color_empty (forest_sans_color amplitudes); Tree.inner = List.map feynmf_set amplitudes } let momentum_to_TeX (_, p) = String.concat "" (List.map p2s p) let wf_to_TeX (f, _ as wf) = M.flavor_to_TeX f ^ "(" ^ momentum_to_TeX wf ^ ")" let amplitudes_to_feynmf latex name amplitudes = Tree.feynmf_sets_wrapped latex name wf_to_TeX momentum_to_TeX variable' format_p (List.map uncolored_colored (amplitudes_by_flavor amplitudes)) let amplitudes_to_feynmf_sans_color latex name amplitudes = Tree.feynmf_sets_wrapped latex name wf_to_TeX momentum_to_TeX variable' format_p (List.map uncolored_only (amplitudes_by_flavor amplitudes)) let amplitudes_to_feynmf_color_only latex name amplitudes = Tree.feynmf_sets_wrapped latex name wf_to_TeX momentum_to_TeX variable' format_p (List.map colored_only (amplitudes_by_flavor amplitudes)) let debug (str, descr, opt, var) = [ "-warning:" ^ str, Arg.Unit (fun () -> var := (opt, false):: !var), " check " ^ descr ^ " and print warning on error"; "-error:" ^ str, Arg.Unit (fun () -> var := (opt, true):: !var), " check " ^ descr ^ " and terminate on error" ] let rec include_goldstones = function | [] -> false | (T.Gauge, _) :: _ -> true | _ :: rest -> include_goldstones rest let read_lines_rev file = let ic = open_in file in let rev_lines = ref [] in let rec slurp () = rev_lines := input_line ic :: !rev_lines; slurp () in try slurp () with | End_of_file -> close_in ic; !rev_lines let read_lines file = List.rev (read_lines_rev file) type cache_mode = | Cache_Default | Cache_Initialize of string let cache_option = ref Cache_Default let unphysical_polarization = ref None (* \thocwmodulesection{Main Program} *) let main () = (* Delay evaluation of [M.external_flavors ()]! *) let usage () = "usage: " ^ Sys.argv.(0) ^ " [options] [" ^ String.concat "|" (List.map M.flavor_to_string (ThoList.flatmap snd (M.external_flavors ()))) ^ "]" and rev_scatterings = ref [] and rev_decays = ref [] and cascades = ref [] and checks = ref [] and output_file = ref None and print_forest = ref false and template = ref false and diagrams_all = ref None and diagrams_sans_color = ref None and diagrams_color_only = ref None and diagrams_LaTeX = ref false and quiet = ref false and write = ref true and params = ref false and poles = ref false and dag_out = ref None and dag0_out = ref None and phase_space_out = ref None in Options.parse (Options.cmdline "-target:" T.options @ Options.cmdline "-model:" M.options @ Options.cmdline "-fusion:" CF.options @ ThoList.flatmap debug ["a", "arguments", T.All, checks; "n", "# of input arguments", T.Arguments, checks; "m", "input momenta", T.Momenta, checks; "g", "internal Ward identities", T.Gauge, checks] @ [("-o", Arg.String (fun s -> output_file := Some s), "file write to given file instead of /dev/stdout"); ("-scatter", Arg.String (fun s -> rev_scatterings := s :: !rev_scatterings), "expr in1 in2 -> out1 out2 ..."); ("-scatter_file", Arg.String (fun s -> rev_scatterings := read_lines_rev s @ !rev_scatterings), "name each line: in1 in2 -> out1 out2 ..."); ("-decay", Arg.String (fun s -> rev_decays := s :: !rev_decays), "expr in -> out1 out2 ..."); ("-decay_file", Arg.String (fun s -> rev_decays := read_lines_rev s @ !rev_decays), "name each line: in -> out1 out2 ..."); ("-cascade", Arg.String (fun s -> cascades := s :: !cascades), "expr select diagrams"); ("-initialize", Arg.String (fun s -> cache_option := Cache_Initialize s), "dir precompute lookup tables and store them in directory"); ("-unphysical", Arg.Int (fun i -> unphysical_polarization := Some i), "n use unphysical polarization for n-th particle / test WIs"); ("-template", Arg.Set template, " write a template for handwritten amplitudes"); ("-forest", Arg.Set print_forest, " Diagrammatic expansion"); ("-diagrams", Arg.String (fun s -> diagrams_sans_color := Some s), "file produce FeynMP output for Feynman diagrams"); ("-diagrams:c", Arg.String (fun s -> diagrams_color_only := Some s), "file produce FeynMP output for color flow diagrams"); ("-diagrams:C", Arg.String (fun s -> diagrams_all := Some s), "file produce FeynMP output for Feynman and color flow diagrams"); ("-diagrams_LaTeX", Arg.Set diagrams_LaTeX, " enclose FeynMP output in LaTeX wrapper"); ("-quiet", Arg.Set quiet, " don't print a summary"); ("-summary", Arg.Clear write, " print only a summary"); ("-params", Arg.Set params, " print the model parameters"); ("-poles", Arg.Set poles, " print the Monte Carlo poles"); ("-dag", Arg.String (fun s -> dag_out := Some s), " print minimal DAG"); ("-full_dag", Arg.String (fun s -> dag0_out := Some s), " print complete DAG"); ("-phase_space", Arg.String (fun s -> phase_space_out := Some s), " print minimal DAG for phase space")]) (*i ("-T", Arg.Int Topology.Binary.debug_triplet, ""); ("-P", Arg.Int Topology.Binary.debug_partition, "")]) i*) (fun _ -> prerr_endline (usage ()); exit 1) usage; let cmdline = String.concat " " (List.map ThoString.quote (Array.to_list Sys.argv)) in let output_channel, close_output_channel = match !output_file with | None -> (stdout, fun () -> ()) | Some name -> let oc = open_out name in (oc, fun () -> close_out oc) in let processes = try ThoList.uniq (List.sort compare (match List.rev !rev_scatterings, List.rev !rev_decays with | [], [] -> [] | scatterings, [] -> Proc.expand_scatterings (List.map Proc.parse_scattering scatterings) | [], decays -> Proc.expand_decays (List.map Proc.parse_decay decays) | scatterings, decays -> invalid_arg "mixed scattering and decay!")) with | Invalid_argument s -> begin Printf.eprintf "O'Mega: invalid process specification: %s!\n" s; flush stderr; [] end in (* \begin{dubious} This is still crude. Eventually, we want to catch \emph{all} exceptions and write an empty (but compilable) amplitude unless one of the special options is selected. \end{dubious} *) begin match processes, !cache_option, !params with | [], Cache_Initialize dir, false -> F.initialize_cache dir; exit 0 | _, _, true -> if !write then T.parameters_to_channel output_channel; exit 0 | [], _, false -> if !write then T.amplitudes_to_channel cmdline output_channel !checks CF.empty; exit 0 | _, _, false -> let selectors = let fin, fout = List.hd processes in C.to_selectors (C.of_string_list (List.length fin + List.length fout) !cascades) in let amplitudes = try begin match F.check_charges () with | [] -> () | violators -> let violator_strings = String.concat ", " (List.map (fun flist -> "(" ^ String.concat "," (List.map M.flavor_to_string flist) ^ ")") violators) in failwith ("charge violating vertices: " ^ violator_strings) end; CF.amplitudes (include_goldstones !checks) !unphysical_polarization CF.no_exclusions selectors processes with | Fusion.Majorana -> begin Printf.eprintf "O'Mega: found Majorana fermions, switching representation!\n"; flush stderr; close_output_channel (); Arg.current := 0; raise Fusion.Majorana end | exc -> begin Printf.eprintf "O'Mega: exception %s in amplitude construction!\n" (Printexc.to_string exc); flush stderr; CF.empty; end in if !write then T.amplitudes_to_channel cmdline output_channel !checks amplitudes; if not !quiet then begin List.iter (fun amplitude -> Printf.eprintf "SUMMARY: %d fusions, %d propagators" (F.count_fusions amplitude) (F.count_propagators amplitude); flush stderr; Printf.eprintf ", %d diagrams" (F.count_diagrams amplitude); Printf.eprintf "\n") (CF.processes amplitudes); let couplings = List.fold_left (fun acc p -> let fusions = ThoList.flatmap F.rhs (F.fusions p) and brakets = ThoList.flatmap F.ket (F.brakets p) in let couplings = VSet.of_list (List.map F.coupling (fusions @ brakets)) in VSet.union acc couplings) VSet.empty (CF.processes amplitudes) in Printf.eprintf "SUMMARY: %d vertices\n" (VSet.cardinal couplings); let ufo_couplings = VSet.fold (fun v acc -> match v with | Coupling.Vn (Coupling.UFO (_, v, _, _, _), _, _) -> Sets.String.add v acc | _ -> acc) couplings Sets.String.empty in if not (Sets.String.is_empty ufo_couplings) then Printf.eprintf "SUMMARY: %d UFO vertices: %s\n" (Sets.String.cardinal ufo_couplings) (String.concat ", " (Sets.String.elements ufo_couplings)) end; if !poles then begin List.iter (fun amplitude -> W.write output_channel "omega" (W.merge (W.trees amplitude))) (CF.processes amplitudes) end; begin match !dag0_out with | Some name -> let ch = open_out name in List.iter (F.tower_to_dot ch) (CF.processes amplitudes); close_out ch | None -> () end; begin match !dag_out with | Some name -> let ch = open_out name in List.iter (F.amplitude_to_dot ch) (CF.processes amplitudes); close_out ch | None -> () end; begin match !phase_space_out with | Some name -> let ch = open_out name in begin try List.iter (fun (fin, fout) -> Printf.fprintf ch "%s -> %s ::\n" (String.concat " " (List.map M.flavor_to_string fin)) (String.concat " " (List.map M.flavor_to_string fout)); match fin with | [] -> failwith "Omega(): phase space: no incoming particles" | [f] -> PHS.phase_space_channels ch (PHS.amplitude_sans_color false PHS.no_exclusions selectors fin fout) | [f1; f2] -> PHS.phase_space_channels ch (PHS.amplitude_sans_color false PHS.no_exclusions selectors fin fout); PHS.phase_space_channels_flipped ch (PHS.amplitude_sans_color false PHS.no_exclusions selectors [f2; f1] fout) | _ -> failwith "Omega(): phase space: 3 or more incoming particles") processes; close_out ch with | exc -> begin close_out ch; Printf.eprintf "O'Mega: exception %s in phase space construction!\n" (Printexc.to_string exc); flush stderr end end | None -> () end; if !print_forest then List.iter (fun amplitude -> List.iter (fun t -> Printf.eprintf "%s\n" (Tree.to_string (Tree.map (fun (wf, _) -> variable wf) (fun _ -> "") t))) (F.forest (List.hd (F.externals amplitude)) amplitude)) (CF.processes amplitudes); begin match !diagrams_all with | Some name -> amplitudes_to_feynmf !diagrams_LaTeX name amplitudes | None -> () end; begin match !diagrams_sans_color with | Some name -> amplitudes_to_feynmf_sans_color !diagrams_LaTeX name amplitudes | None -> () end; begin match !diagrams_color_only with | Some name -> amplitudes_to_feynmf_color_only !diagrams_LaTeX name amplitudes | None -> () end; close_output_channel (); exit 0 end (* \begin{dubious} This was only intended for debugging O'Giga \ldots \end{dubious} *) let decode wf = (F.flavor wf, (F.momentum wf : Momentum.Default.t)) let diagrams in1 in2 out = match F.amplitudes false F.no_exclusions C.no_cascades [in1; in2] out with | a :: _ -> let wf1 = List.hd (F.externals a) and wf2 = List.hd (List.tl (F.externals a)) in let wf2 = decode wf2 in List.map (fun t -> (wf2, Tree.map (fun (wf, _) -> decode wf) decode t)) (F.forest wf1 a) | [] -> [] let diagrams in1 in2 out = failwith "Omega().diagrams: disabled" end Index: trunk/share/tests/functional_tests/ref-output/bjet_cluster.ref =================================================================== --- trunk/share/tests/functional_tests/ref-output/bjet_cluster.ref (revision 8498) +++ trunk/share/tests/functional_tests/ref-output/bjet_cluster.ref (revision 8499) @@ -1,111 +1,111 @@ ?openmp_logging = false ?vis_history = false ?integration_timer = false ?pacify = true ?omega_write_phs_output = true | Switching to model 'SM', scheme 'default' SM.alphas => 1.18000E-01 ?alphas_is_fixed = false ?alphas_from_mz = true ?alphas_from_lambda_qcd = false alphas_nf = 5 SM.ms => 0.00000E+00 SM.mc => 0.00000E+00 [user variable] lightjet = PDG(2, -2, 1, -1, 3, -3, 4, -4, 21) [user variable] jet = PDG(2, -2, 1, -1, 3, -3, 4, -4, 21, 5, -5) $phs_method = "fast_wood" $restrictions = "!H" | Process library 'bjet_cluster_lib': recorded process 'bjet_cluster_p1' seed = 1234 sqrts = 1.00000E+03 jet_algorithm = 2 jet_r = 5.00000E-01 | Integrate: current process library needs compilation | Process library 'bjet_cluster_lib': compiling ... | Process library 'bjet_cluster_lib': writing makefile | Process library 'bjet_cluster_lib': removing old files | Process library 'bjet_cluster_lib': writing driver | Process library 'bjet_cluster_lib': creating source code | Process library 'bjet_cluster_lib': compiling sources | Process library 'bjet_cluster_lib': linking | Process library 'bjet_cluster_lib': loading | Process library 'bjet_cluster_lib': ... success. | Integrate: compilation done | QCD alpha: using a running strong coupling | RNG: Initializing TAO random-number generator | RNG: Setting seed for random-number generator to 1234 | Initializing integration for process bjet_cluster_p1: | Beam structure: [any particles] | Beam data (collision): | e- (mass = 5.1099700E-04 GeV) | e+ (mass = 5.1099700E-04 GeV) | sqrts = 1.000000000000E+03 GeV | Phase space: generating configuration ... | Phase space: ... success. | Phase space: writing configuration file 'bjet_cluster_p1.i1.phs' | ------------------------------------------------------------------------ | Process [scattering]: 'bjet_cluster_p1' | Library name = 'bjet_cluster_lib' | Process index = 1 | Process components: | 1: 'bjet_cluster_p1_i1': e-, e+ => b, bbar, u:ubar:d:dbar:s:sbar:c:cbar:gl, u:ubar:d:dbar:s:sbar:c:cbar:gl [omega] | ------------------------------------------------------------------------ | Phase space: 68 channels, 8 dimensions | Phase space: found 68 channels, collected in 12 groves. | Phase space: Using 104 equivalences between channels. | Phase space: wood | Applying user-defined cuts. | Using user-defined general scale. | Starting integration for process 'bjet_cluster_p1' | Integrate: iterations = 1:760:"gw" | Integrator: 12 chains, 68 channels, 8 dimensions | Integrator: Using VAMP channel equivalences | Integrator: 760 initial calls, 20 bins, stratified = T | Integrator: VAMP |=============================================================================| | It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] | |=============================================================================| - 1 748 3.021E+01 2.23E+01 73.79 20.18 9.5 + 1 748 2.862E+01 2.22E+01 77.61 21.23 9.7 |-----------------------------------------------------------------------------| - 1 748 3.021E+01 2.23E+01 73.79 20.18 9.5 + 1 748 2.862E+01 2.22E+01 77.61 21.23 9.7 |=============================================================================| | QCD alpha: using a running strong coupling | RNG: Initializing TAO random-number generator | RNG: Setting seed for random-number generator to 1235 | Initializing integration for process bjet_cluster_p1: | Beam structure: [any particles] | Beam data (collision): | e- (mass = 5.1099700E-04 GeV) | e+ (mass = 5.1099700E-04 GeV) | sqrts = 1.000000000000E+03 GeV | Phase space: generating configuration ... | Phase space: ... success. | Phase space: writing configuration file 'bjet_cluster_p1.i1.phs' | ------------------------------------------------------------------------ | Process [scattering]: 'bjet_cluster_p1' | Library name = 'bjet_cluster_lib' | Process index = 1 | Process components: | 1: 'bjet_cluster_p1_i1': e-, e+ => b, bbar, u:ubar:d:dbar:s:sbar:c:cbar:gl, u:ubar:d:dbar:s:sbar:c:cbar:gl [omega] | ------------------------------------------------------------------------ | Phase space: 68 channels, 8 dimensions | Phase space: found 68 channels, collected in 12 groves. | Phase space: Using 104 equivalences between channels. | Phase space: wood | Applying user-defined cuts. | Using user-defined general scale. | Starting integration for process 'bjet_cluster_p1' | Integrate: iterations = 1:760:"gw" | Integrator: 12 chains, 68 channels, 8 dimensions | Integrator: Using VAMP channel equivalences | Integrator: 760 initial calls, 20 bins, stratified = T | Integrator: VAMP |=============================================================================| | It Calls Integral[fb] Error[fb] Err[%] Acc Eff[%] Chi2 N[It] | |=============================================================================| - 1 748 2.365E+01 4.13E+00 17.46 4.78 13.6 + 1 748 2.517E+01 4.46E+00 17.72 4.85 13.0 |-----------------------------------------------------------------------------| - 1 748 2.365E+01 4.13E+00 17.46 4.78 13.6 + 1 748 2.517E+01 4.46E+00 17.72 4.85 13.0 |=============================================================================| | WHIZARD run finished. |=============================================================================|