Implemented ContributingDiagrams[]

meta/EDM.m

@@ -209,7 +209,141 @@ If you add new kinds of vertices (e.g for new diagram types):
        Return[code];
        ];
 
+<<<<<<< 379db5394ffb4373e1ea0bcc6a72e5a0ce7ee357
 (* Find all diagrams of the type type_, testing all corresponding combinations of fields *)
+=======
+(************************ Begin helper routines *******************************)
+
+(* Return the name of the SARAH particle family containing the muon *)
+GetMuonFamily[] := If[TreeMasses`GetDimension[SARAH`Electron] =!= 1,
+                        SARAH`Electron,
+                        Cases[SARAH`ParticleDefinitions[FlexibleSUSY`FSEigenstates],
+                              {p_, {Description -> "Muon", ___}} -> p, 1][[1]]
+                        ];
+(* If the muon has a generation index, return it, otherwise return Null.
+ e.g. CMSSMNoFV has no muon generation index *)
+GetMuonIndex[] := If[TreeMasses`GetDimension[SARAH`Electron] =!= 1, 2, Null];
+
+GetPhoton[] := SARAH`Photon;
+
+IsLorentzIndex[index_] := StringMatchQ[ToString @ index, "lt" ~~ __];
+
+StripLorentzIndices[p_Symbol] := p;
+StripLorentzIndices[SARAH`bar[p_]] := SARAH`bar[StripLorentzIndices[p]];
+StripLorentzIndices[Susyno`LieGroups`conj[p_]] := Susyno`LieGroups`conj[StripLorentzIndices[p]];
+StripLorentzIndices[p_] := Module[{remainingIndices},
+                                  remainingIndices = Select[p[[1]], (!IsLorentzIndex[#] &)];
+                                  If[Length[remainingIndices] === 0, Head[p],
+                                     Head[p][remainingIndices]]
+                                  ];
+SetAttributes[StripLorentzIndices, {Listable}];
+
+(* Takes a SARAH particle and returns its antiparticle *)
+AntiParticle[SARAH`bar[p_]] := p;
+AntiParticle[Susyno`LieGroups`conj[p_]] := p;
+AntiParticle[p_] := Module[{pNoIndices = Vertices`StripFieldIndices[p]},
+                           If[IsScalar[pNoIndices] || IsVector[pNoIndices],
+                              Susyno`LieGroups`conj[p],
+                              SARAH`bar[p]]];
+SetAttributes[AntiParticle, {Listable}];
+
+(* Return a string corresponding to the c++ class name of the particle.
+ Note that "bar" and "conj" get turned into anti<...>::type! *)
+ParticleToCXXName[p_] := SymbolName[p];
+ParticleToCXXName[SARAH`bar[p_]] := "anti<" <> SymbolName[p] <> ">::type";
+ParticleToCXXName[Susyno`LieGroups`conj[p_]] := "anti<" <> SymbolName[p] <> ">::type";
+
+(* Return a string corresponding to the name of the particle.
+ Note that "bar" and "conj" are left as they are! *)
+ParticleToSARAHString[p_] := SymbolName[p];
+ParticleToSARAHString[SARAH`bar[p_]] := "bar" <> SymbolName[p];
+ParticleToSARAHString[Susyno`LieGroups`conj[p_]] := "conj" <> SymbolName[p];
+
+(* Change Symbol[indexNameNUMBEROLD] to Symbol[indexNameNUMBERNEW] *)
+ChangeIndexNumber[index_Symbol, number_Integer] := Symbol @ StringReplace[ToString[index],
+                       Shortest[name__] ~~ NumberString :> name ~~ ToString[number]];
+
+(* Returns the rules to change p[{iN,jN,kN}] to p[{iM,jM,kM}] *)
+IndexReplacementRulesForNewParticleIndexNumber[particle_, number_Integer] :=
+    ((# -> ChangeIndexNumber[#, number] &) /@ Vertices`FieldIndexList[particle]);
+
+(* Returns the rules that are needed to change order of a prticle list.
+ i.e turn {p1[{i1,j1,k1}], p2[{l2,m2,n2}]} to {p2[{l1,m1,n1}], p1[{i2,j2,k2}]}
+ For some reason SARAH expects the indices to always be in that order... *)
+IndexReplacementRulesForParticleReordering[particles_List, ordering_] :=
+Module[{indices = Table[i, {i, Length[particles]}], index, fieldIndexList},
+       Flatten[(IndexReplacementRulesForNewParticleIndexNumber[particles[[ordering[[#]]]], #] &) /@ indices]];
+
+(* Perform the actual transformation initiated by IndexReplacementRulesForParticleReordering[] *)
+OrderParticles[particles_List, ordering_] := Module[{indexRules},
+       indexRules = IndexReplacementRulesForParticleReordering[particles, ordering];
+       Return[particles[[ordering]] /. indexRules];
+                                                    ];
+
+(* Essentially the same as above, but using the obtained rules,
+ transform a whole vertex expression, as returned by SARAH`Vertex[] *)
+OrderVertex[vertex_, ordering_] :=
+    Module[{indexRules, particles, expr, newVertex},
+           indexRules = IndexReplacementRulesForParticleReordering[vertex[[1]], ordering];
+
+           particles = vertex[[1]][[ordering]];
+           expr = vertex[[2;;]];
+
+           newVertex = (Join[{particles}, expr] /. indexRules);
+           Return[newVertex];
+          ];
+
+(* MemoizingVertex[] works just like SARAH`Vertex[], but it caches the results *)
+(* MemoizingVertex[] only works when __no__ indices are specified!!! *)
+(* Use of memoization gives ~30% speedup for the MSSM! *)
+memoizedVertices = {};
+MemoizingVertex[particles_List, options : OptionsPattern[SARAH`Vertex]] :=
+    Module[{memo, ordering, orderedParticles},
+           (* First we sort the particles *)
+           ordering = Ordering[particles];
+           orderedParticles = particles[[ordering]];
+
+           memo = Select[memoizedVertices, MatchesMemoizedVertex[orderedParticles], 1];
+           If[memo =!= {}, memo = memo[[1]],
+              (* Create a new entry *)
+              memo = SARAH`Vertex[orderedParticles, options];
+              AppendTo[memoizedVertices, memo];];
+
+           (* Now return the particles to their original order *)
+           memo = OrderVertex[memo, Ordering[ordering]];
+           Return[memo]];
+
+MatchesMemoizedVertex[particles_List][vertex_] := MatchQ[particles, Vertices`StripFieldIndices /@ vertex[[1]]];
+
+(* Test whether a SARAH Vertex[] result is nonzero (exact) *)
+IsNonZeroVertex[v_] := MemberQ[v[[2 ;;]][[All, 1]], Except[0]];
+
+(* Returns the name of the coupling function that FlexibleSUSY generates for
+ a specific vertex in a canonical order! *)
+NameOfCouplingFunction[particles_List] :=
+((* FIXME: Not upwards compatible if naming conventions change *)
+ "Cp" <> StringJoin @ (ParticleToSARAHString /@ Sort[particles]));
+
+(********************** End helper routines **************************)
+
+(* The different vertex types that are supported.
+ They have the same names as their c++ counterparts. *)
+vertexTypes = {
+    SingleComponentedVertex,
+    LeftAndRightComponentedVertex
+};
+
+(* The different diagram types that should be taken into consideration *)
+(* They need to be called DIAGRAMTYPENAME[_Integer]! See CreateDiagramClasses[] below. *)
+(* There is no bounds check done on the integers, so they have to fit
+ into a standard c++ unsigned (!) int *)
+contributingDiagramTypes = {
+    OneLoopDiagram[0],
+    OneLoopDiagram[1]
+};
+
+(* Find all diagrams of the type type_, testing all corresponding combinations of particles *)
+>>>>>>> Implemented ContributingDiagrams[]
 (* IMPORTANT: Return value should have the format
  {{edmField1, {Diagram[DIAGRAMTYPENAME[_Integer], Fields___], Diagram[...], ...}},
   {edmField2, {...}},