(* *********************************************************** *)
(*                                                             *)
(**                    InsertFields.m                         **)
(*                                                             *)
(* by Hagen Eck and Sepp Kueblbeck (last change: 06/91)        *)
(*                                                             *)
(* *********************************************************** *)

BeginPackage["InsertFields`","FeynArts`"]

Needs["Utilities`Compare`"];
Needs["Utilities`FeynInit`"];

InsertFields::usage = "
InsertFields[topology, extli, options] constructs all Feynman- \n
diagrams of a Topology or a TopologyList with external fields  \n
given in `extli'. extli must be of the  form: \n
{ code(particle_1), ..., code(particle_i) } -> \n
      { code(particle_i+1), ..., code(particle_n) }" 

Begin["`Private`"] 

Options[ InsertFields ] = { Model             -> {SM},
			    Generation1	      -> True,
			    Generation2       -> True,
			    Generation3       -> True,
			    ElectronHCoupling -> True,
			    LightFHCoupling   -> True,
			    QuarkMixing       -> False,
			    UnitaryGauge      -> False,
			    RemoveEmptyTops   -> True,
			    ScreenMessages    -> True,
			    ProcessName       -> Automatic
			  }
                                         (** initialization **)

AppendFieldVar[ pr_[a_,b_,c___],i_Integer ] := pr[ a,b,fi[i],c ]

PropagatorSort[a_] :=  If[ Length[a]===2,
		           Sort[a],
		           Append[ Sort[Take[a,2]],a[[3]] ]
                         ]

ToType[ Incoming ] = inc;
ToType[ Outgoing ] = out;

InitTopology[ mod_ ] [ tt:Topology[_][___] , extli_ ] := 
 Module[{to},
       to  =  Sort[ (PropagatorSort /@ tt) /. ex->AAA
		  ] ;
       Apply[  Head[to],
	       Append[  Array[ AppendFieldVar[to[[#]],#]&,Length[to] ],
                        Insertions[Graph[]]
                     ]        
            ] //.{
  t_[ pp___,
      Propagator[AAA][a___,e[i_Integer],b___,fi[j_],c___],
      qq___,
      Insertions[Graph[v___]]    ]                        :>
  t[   pp,
       Propagator[ToType[Head[extli[[i]]]]][ a,e[i],b,fi[j],c ],
       qq,
       Insertions[ Graph[v,fi[j]->extli[[i]] /.
		 {Incoming->Identity,Outgoing->AntiParticle[mod]}] ] ]}
 ]

(* Insert compatible particles in 1 propagator for 1 set of rules  *)

Ins11[ vert2_ , ru_ , mod_ , i_Integer ] := 
 Module[{int,vx,newvertex,possiblevx1,possiblevx2},
       vx  = vert2 /. Apply[List,ru] ;
       possiblevx1 = Union @@ (PossibleParticles[#][vx[[1]]]& /@ mod);
       possiblevx2 = Union @@ (PossibleParticles[#][vx[[2]]]& /@ mod);
       int = Intersection[ AntiParticle[mod] /@ possiblevx1, possiblevx2 ];
       If[ vx[[1]] === vx[[2]],
           int = Select[ int, 
                         ( newvertex = vx[[1]]/.fi[i]-># ;
                           (!FreeQ[newvertex,fi]) ||
                           (MemberQ[ VertexList[mod], newvertex] )
                         )&
                       ]
         ]; 
       Apply[ Insertions,
              Array[ Append[ru,fi[i]->int[[#]]]&,Length[int] ]  
            ]     
      ]; 

FindFields[ver_,pr_[ver_,_   ,c_]]:=Vertex[ APa[c]   ];
FindFields[ver_,pr_[_   ,ver_,c_]]:=Vertex[        c ]; 
FindFields[ver_,pr_[ver_,ver_,c_]]:=Vertex[ APa[c],c ];
FindFields[ver_,pr_[_   ,_   ,c_]]:=Vertex[          ];

        (** insert the first i propagators (recursive definition) **)

Ins1[ tt : Topology[ _ ][ ___ ] , mod_ , i_Integer , extlength_Integer ] := 
 Module[ { t,rules,props,vert2 } ,
  If[ i === extlength ,
   tt ,
   t     = Ins1[ tt , mod , i-1 , extlength ] ;
   rules = Last[ t ];
   props = Take[ t , Length[t]-1 ] ;
   vert2 = {Flatten[ Apply[ Vertex,
                            Map[FindFields[props[[i,1]],#]&,props] ]],
            Flatten[ Apply[ Vertex,
	                    Map[FindFields[props[[i,2]],#]&,props] ]] 
           } /. APa->AntiParticle[mod] ;
   Append[ props,
           Flatten[Union[  Map[ Ins11[vert2,#,mod,i]&,rules ]  ]]
         ]
  ]
 ];

InsertFields::badparticle = "
Particle `1` does not live in models `2`. "

InsertFields::extnumber = "
You cannot fit `1` external particle(s) onto a `2` leg topology. "

                              (** supports equation form of extli **)

InsertFields[ tt_, Rule[ a_,b_ ] , Opt___Rule ] :=
 Module[ {insresult},
   insresult = InsertFields[  tt, 
                              Join[ List@@(Incoming/@a), 
                                    List@@(Outgoing/@b) ], 
                              Opt                               ];
   Append[ Head[insresult], Process->Rule[a,b] ] @@ insresult
 ] /; Head[a]===Head[b]

                              (** insert fields into 1 topology **)

InsertFields[ t:Topology[ _ ][ ___ ] , extli_List, modOpt___Rule  ] := 
 Module[ { modelname , localoptions , holdoptions ,
	  momli , ttt , tt , to , too , i=1 , topextlength } ,
  If[ Length[{modOpt}] > 0,
      holdoptions = Options[ InsertFields ];
      localoptions = SetOptions[ InsertFields, modOpt  ];
      Options[ InsertFields ] = holdoptions, 
      localoptions = Options[ InsertFields ]
    ];
  topextlength = Count[ t,Propagator[ex][___] ];
  If[ Not[ Length[extli]===topextlength ],
      Message[InsertFields::extnumber,Length[extli],topextlength];
      Return[$Aborted] ] ;
  If[ (modelname = CheckModel[localoptions])===$Aborted,Return[$Aborted] ];
  While[ i<=Length[extli],
	 If[ !MemberQ[ Union @@ (ParticleList[#]& /@ modelname),extli[[i,1]] ],
             Message[InsertFields::badparticle,
		     extli[[i,1]],
		     Model/.localoptions];
	     Return[$Aborted] ];
         i++
       ];	     
  ttt = InitTopology[ modelname ][ t , extli ] ;
  If[  Length[ Select[Drop[ttt,-1],FreeQ[#,e]&] ]==0,
       ttt=If[ MemberQ[VertexList[modelname],
		       Vertex @@ ( Array[Last[ttt][[1,#,2]]&,
	                                 Length[ Last[ttt][[1]] ]
                                        ] )
                      ]===False,
               ttt/.Insertions[__]->Insertions[],
	       ttt
             ]
    ];
  tt    = Append[  Map[ Take[#,3]&,Take[ttt,Length[t]] ],Last[ttt]  ] ;
  momli = Map   [ Drop[#,3]&,Take[ttt,Length[t]] ] ;
  to    = Ins1  [ tt,modelname,Length[t],Length[extli] ] ;
  Print   [  Length[ Last[to] ]," insertions before Compare "  ] ;
  too   = InsertionsCompare[ to,modelname ] ;
  Print   [  Length[ too ]     ," insertions after  Compare "  ] ;
  Print   [] ;
  too = Apply[  Topology,
          Append[ Array[ Head[tt[[#]]]@@ (
			    Join[List@@tt[[#]],List@@momli[[#]]] )&,
			 Length[tt]-1 ],
(*                  Map[Sort,too]               *)
                  Map[Identity,too]
       ] 	]  ;
  localoptions = localoptions/. Rule[ProcessName,Automatic]:>Rule[
            ProcessName,AutoProcessName[extli] ];
  (TopologyList@@ Select[localoptions,(#[[1]]=!=ScreenMessages &&
				       #[[1]]=!=RemoveEmptyTops   )& ]
	) @ too
 ]
                    (** insert fields into a list of topologies **)

InsertFields[ topl_TopologyList , 
              extli_List , modOpt___Rule ] :=
 Module[{result,holdoptions,localoptions,modelname,mod,
         i=1,topextlength,deltalength},
  If[ Length[{modOpt}] > 0,
      holdoptions = Options[ InsertFields ];
      localoptions = SetOptions[ InsertFields, modOpt  ];
      Options[ InsertFields ] = holdoptions, 
      localoptions = Options[ InsertFields ]
    ];
  topextlength = Count[ topl[[1]],Propagator[ex][___] ];
  If[ Not[ Length[extli]===topextlength ],
      Message[InsertFields::extnumber,Length[extli],topextlength];
      Return[$Aborted] ] ;
  If[ (modelname = CheckModel[localoptions])===$Aborted,Return[$Aborted] ];
  While[ i<=Length[extli],
	 If[ !MemberQ[ Union @@ (ParticleList[#]& /@ modelname),extli[[i,1]] ],
             Message[InsertFields::badparticle,
		     extli[[i,1]],
		     Model/.localoptions];
	     Return[$Aborted] ];
         i++
       ];	    
  mod = Select[ localoptions, FreeQ[#,ScreenMessages]& ];
  If[ (ProcessName/.mod)===Automatic,
      mod = mod/.{Automatic->AutoProcessName[extli]}
    ];
  result = Apply[  TopologyList@@mod,
  	    	   Array[ (Print[#,". topology (out of ",
		      	         Length[topl],") being inserted "];
                          InsertFields[topl[[#]],extli,modOpt][[1]])&,
	       	          Length[topl] ]  ];
              (*  this one discards the `empty' insertions :   *)
  Print[" "];
  If[ (RemoveEmptyTops/.mod),
       le=Length[result];
       deltalength = le - Length[  result = Select[ result , 
                                                    FreeQ[#,Insertions[]]& ] ];
       If[ deltalength>0,
           Print["Out of ",le,If[ le===1," top. "," tops. "],
                 deltalength,  	 
                 " have no inss. and are discarded" ]
         ]
    ] ;
  Print[ "total number of insertions: " , 
	 Plus  @@  ( Length[ Last[#] ]& /@ result ) ];
				 (* switch on messages again *)
  SetOptions[FeynInit`FeynInit,ScreenMessages->True];
  $Output = {"stdout"};
  Return[ result ]
 ]

AutoProcessName[ x_List ] :=
 Module[{name},
   Off[General::spell1];
   name =
    ToExpression[
    StringJoin @@ ( Flatten[ x //.{
           H_[ n_Integer?Positive ]:> (ToString[H] <> ToString[n]),
      H_[ n_Integer ]         :> ("a" <> ToString[H] <> ToString[Abs[n]]) ,
      H_[ I]    :> (ToString[H] <> "I"),
      H_[-I]         :> ("a" <> ToString[H] <> "I"),
      Incoming->List,
           Outgoing->List
                                   } ] ) ];
   On[General::spell1];
   name
 ]
                                   (** compare all insertions **)

LongInsertion[ mod_ ][ b___ ] := 
  Array[ List[b][[#,1]] -> Sort[{ List[b][[#,2]],
	                          AntiParticle[mod][ List[b][[#,2]] ]
			       }] [[1]] &,
	 Length[List[b]] 
       ]

InsertionsCompare[ Topology[___][___,Insertions[]],mod_ ] := 
		    Insertions[]

InsertionsCompare[ t:Topology[1][___],mod_ ] :=  
           Last[t]  /.  { Topology[1] -> Topology , Graph -> Graph[1] }

InsertionsCompare[ to : Topology[ _Integer ][ __ ] , mod_ ] := 
 Module[ { top,ru,inli,lili,fiinli,all,cases,pos,ar } ,
  top    = Drop [ to,{-1} ] ;
  ru     = Last [ to ] /. Graph->LongInsertion[mod] ;
  inli   = Select [ top,!FreeQ[#,in]& ] ;
  lili   = Select [ top/.{inc->in,out->in},FreeQ[#,in]& ] ;
  fiinli = Array [ inli[[#,3]]&,Length[inli] ] ;
  filili = Array [ lili[[#,3]]&,Length[lili] ] ;
  all    = Array [  Join[ fiinli/.ru[[#]],
	       	          Sort[filili/.ru[[#]]] ] & , 
                    Length[ru] ] ;
  cases  = Union[  all  ] ;
  pos    = Array[  Position[ all,cases[[#]] ]&,Length[cases]  ];
  ar     = 
    Array[ 
     Function[z,
              Append[top,
                     Insertions @@
			Array[ Last[to][[ Flatten[ pos[[z]] ][[#]] ]]&,
	                       Length[ Flatten[pos[[z]]] ] 
	                     ]  
             ]      ],
          Length[pos]
	 ];
  Join  @@  ( DisjointInsertionsCompare[ #,mod ]& /@ ar ) 
  ]

        (** comparing disjoint subsets using Utilities`Compare **)

DisjointInsertionsCompare[ to:T_[ __ ] , mod_ ] := 
  Insertions @@ (ExtractInsertion /@ 
       Compare[ Apply[ TopologyList,
                       Array[ Drop[to,{-1}] /.
                                  ( Last[to][[#]] /. Graph -> List ) & ,
	                      Length[ Last[to] ]
                     ]      ] ,
                mod 
	      ])

ExtractInsertion[ t:Topology[ i_Integer ][ ___ ]] :=
        Graph[i] @@ Array[ fi[#] -> t[[#,3]] &, Length[t] ]
 
               (** making insertions for drawing a topology **)

LineSpec[ TopDummy[1] ] = { straight, none, " " }

InsertFields[ tt_TopologyList,nrin_Integer,nrout_Integer ] :=
TopologyList[ Model->{Topology} , 
              ProcessName->ToExpression["T"] ,
              Process->Topology ]@@ (
               InsertFields[ #,nrin,nrout ]& /@ tt)

InsertFields[ tt:Topology[sym_][ props__ ],nri_,nro_ ] :=
   Module[ {rettop},
           rettop = Sort[ (PropagatorSort /@ tt) /. ex->xxt
		        ] /. xxt->ex ;
           rettop = rettop /. IncOutSpecification[ nri,nro ];
           rettop = Array[ AppendFieldVar[rettop[[#]],#]&,Length[tt] ];
           rettop = Append[ rettop,
                   Insertions[
                    Graph[sym]@@ 
                     Array[ fi[#]->TopDummy[1] & , Length[tt] ] ]
                   ];
          Return[ Topology@@ rettop ]
        ]

IncOutSpecification[ nin_, nou_ ] :=
   Join[ Array[ RuleDelayed[
                       Propagator[ex][e[#],vert__] ,
                       Propagator[inc][e[#],vert]   ]& , nin ],
         Array[ RuleDelayed[
                       Propagator[ex][e[#],vert__] ,
                       Propagator[out][e[#],vert]   ]& , nou , nin+1 ]
      ]


End[]

EndPackage[]

