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

BeginPackage["CreateFeynAmp`","FeynArts`"]

Needs["Utilities`FeynInit`"];
Needs["Utilities`AddMomenta`"];
Needs["Utilities`FeynGrassman`"];
Needs["Utilities`VertexHandling`"];
Needs["Utilities`FeynMultiplication`"];
Needs["Utilities`OutputRules`"];

CreateFeynAmp::usage ="
CreateFeynAmp[ insertedTopologyList, options ] creates a list of \n
Feynman amplitudes (`FeynAmpList') of all the insertions created \n
by InsertFields (Head: TopologyList[modeloptions,name]."

FeynmanGauge::usage ="
FeynmanGauge is an option of CreateFeynAmp that only applies to 
tHooft gauge (UnitaryGauge->False). Set to True (default) 
it sets all GaugeXi[i]=1 with i=A,Z,W. "

Truncated::usage ="
Truncated is an option of CreateFeynAmp that determines wether to 
insert the wave functions of external fields (default=True)."
 
$DDimensions::usage ="
$DDimensions is an option for CreateFeynAmp to determine the prefactor 
of a Feynman amplitude."

Begin["`Private`"]

Options[ CreateFeynAmp ] = {
             FeynmanGauge   -> True,
			    $DDimensions   -> False,
			    Truncated	   -> False,
			    ScreenMessages -> True
			   }
Unprotect[FeynAmp];
Format[ FeynAmp[ GraphName[name___],mom___,ampl_ ] ] := 
	"FeynAmp"[name][ "Integral"[mom] ampl/.ReadRules ]  
Protect[FeynAmp];

                           (** omit wave funcions               **) 

TruncationRules = {
                   PolarizationVector[__] -> 1 ,
                   QuarkSpinor[__] -> 1 ,
                   LeptonSpinor[__] -> 1 
                  }
                                 (** times 3 for the colors **)

QuarkFree[x_] := And@@ ( FreeQ[x,#]& /@ { mu,md,mc,ms,mt,mb,
					  Global`MU,Global`MD,
					  Global`MC,Global`MS,
                 Global`MT,Global`MB } )

ColorFactorRule = {
SpinorChain[ sth__ ] :> 3 SpinorChain[ sth ] /; (
	    FreeQ[ {sth}, LeptonSpinor[__] ] 
	 && FreeQ[ {sth}, QuarkSpinor[__] ]
	 && FreeQ[ {sth}, GellMannMatrix ]
    && Not[ QuarkFree[ {sth} ] ] )
}
                                    (** exported function **)

CreateFeynAmp[ Topology[mod___Rule][props__,verslist_Insertions] ,
	       optional___Rule ] := 
   Module[ {localoptions,holdoptions,workoptions,holdworkoptions,
	   reamp,ru,check,momtop,newmodels} ,
         holdoptions = Options[ FeynInit ];
         localoptions = SetOptions[ FeynInit, mod ,
                   Rule[ScreenMessages,ScreenMessages/.{optional}] ];
         options[ FeynInit ] = holdoptions ; 
	 If[ (modelname = CheckModel[localoptions])===$Aborted,
	     Return[$Aborted] ];
	 Set[ WorkingModel, modelname ];
	 check = Length[ {props}[[1]] ];
	 If[ check < 4 ,
	     momtop = AddMomenta[ Topology[props,verslist ] ] ;
	     momtop = Drop[momtop,-1] ,
	     momtop = Topology[ props ]
	   ];
    holdworkoptions=Options[CreateFeynAmp]; 
    workoptions = SetOptions[CreateFeynAmp,optional];
    Options[CreateFeynAmp]=holdworkoptions;
	 ru = verslist /.  IndicesRule  ;
	 reamp = CreateAmplitude[modelname][
		  CrInit[momtop,#,workoptions] , workoptions
		 ]& /@ ru;
	 Return[ FeynAmpList@@reamp ]
	]

ExtendedProcess[ proc_,mtop_] :=
   Module[{mrules,flatproc},
    mrules=Select[ List @@ ( mtop/.{ Propagator[inc][e[i_Integer],_,_,mom_]:>
                                     Rule[e[i],mom] ,
                                     Propagator[out][e[i_Integer],_,_,mom_]:>
                                     Rule[e[i],-mom] } 
                           ) , FreeQ[#,fi]& ];
    header=Head[ proc[[1]] ];
    header @@ Table[ {proc[[1,i]],e[i]/.mrules,Mass[proc[[1,i]]]},
                     {i,1,Length[proc[[1]]]} ]    ->
    header @@ Table[ {proc[[2,i]],e[i+Length[proc[[1]]]]/.mrules,
                                               Mass[proc[[2,i]]]},
                     {i,1,Length[proc[[2]]]} ]    
   ]

CreateFeynAmp[ TopologyList[mod__Rule ][ tops__Topology ] ,
	       options___Rule ] :=
   Module[ {moddd,back,temp,len,name,modelname,n=1,
            topoptions,localoptions,holdoptions,result} ,
	      len = Length[{tops}] ;
         back = TopologyList[] ;
                      (** Add momenta and masses to `Process' **)
         moddd={mod};
         moddd=moddd/.(Rule[Process,proc_]:>
                   Rule[Process,
                        ExtendedProcess[proc,
                                        AddMomenta[{tops}[[1]]] ]  ] );
         name = Cases[ moddd , Rule[ProcessName,x_] -> x ][[1]];
         holdoptions=Options[CreateFeynAmp];
         localoptions=SetOptions[CreateFeynAmp,options];
         Options[CreateFeynAmp]=holdoptions;
         topoptions=Select[ Select[moddd,Not[#[[1]] === Process]& ],
                            FreeQ[#,Process]& ];
         If[ (modelname = CheckModel[  Append[topoptions,
                                           Rule[ScreenMessages,
	       	                           ScreenMessages/.localoptions]] ])===
			        $Aborted,
              Return[Aborted] ];
	      Do[ 
            Print[thestring[ i,len,Length[Last[{tops}[[i]]]] ]];
            If[ {options}==={},
	         temp = CreateFeynAmp[ (Topology[##]& @@ topoptions) 
                                            @@ {tops}[[i]]         ],
                 temp = CreateFeynAmp[ (Topology[##]& @@ topoptions)
                                            @@ {tops}[[i]],options ] ];
            Off[General::spell1];
            temp = Array[ Prepend[ temp[[#]] ,
	       GraphName[ name,
	                  ToExpression[StringJoin["T",ToString[i]]],
			  ToExpression[StringJoin["I",ToString[#]]],
			  ToExpression[StringJoin["N",ToString[n++]]]
                        ] ]& , Length[temp]
                   ];
             On[ General::spell1];
	     back = Append[ back , temp ] ,
	     {i,len}
	   ] ;
         $Output={"stdout"};
         Off[General::spell1];
         result =  ( (FeynAmpList @@ moddd) @@ (
                    Join@@ (back/.FeynAmpList->List) ) )//.LastRules ; 
         On[General::spell1];
         Return[ result ]
	]

                      (** supply with li's and gi's 
			      (Lorentz and gluon indices) **)


AddIndex[ Graph[s_][ elem__ ] , lnum_Integer , gnum_Integer ] := 
   Module[{vugwithout,firstwithout,with,newl=lnum,newg=gnum},
	  vugwithout = Flatten[ Union[ Position[ {elem},Rule[fi[_],V[_]] ],
	                               Position[ {elem},Rule[fi[_],U[5]] ],
	                               Position[ {elem},Rule[fi[_],U[-5]] ],
	                               Position[ {elem},Rule[fi[_],G[_]] ] ] ];
	  If[ vugwithout==={},
	      Return[ Graph[s][elem] ],
              firstwithout = {elem}[[  Min[ vugwithout ]  ]];
              with = firstwithout /. {V[a_]:>V[ a,li[lnum+1],li[newl=lnum+2] ],
                                      U[ 5]:>U[ 5,gi[gnum+1],gi[newg=gnum+2] ],
                                      U[-5]:>U[-5,gi[gnum+1],gi[newg=gnum+2] ],
				      G[a_]:>G[ a,li[lnum+1],li[newl=lnum+2],
				    	          gi[gnum+1],gi[newg=gnum+2] ]};
              AddIndex[ Graph[s][elem]/.(firstwithout->with),newl,newg ]
         ] ]

IndicesRule = Graph[s_][elem__] :> AddIndex[ Graph[s][elem] ,0 ,0 ]

                                        (** optical message **)

thestring[ somuch_ , outof_ ,actlen_ ] := 
   Module[ {backstring} ,
	 backstring="";
	 Do[ backstring=StringJoin[backstring,"."],{somuch}];
	 Do[ backstring=StringJoin[backstring," "],{j,somuch+1,outof}];
	 StringJoin[ backstring , "|(",ToString[actlen],")" ]
	]

    (** extract sym-fac , insert version , search grassman fields **)

CrInit[ acTop:Topology[ ps__ ] , actvers:Graph[ sym_ ][ ru__ ] ,
	opt_List ] :=
   Module[ {retop}  ,
	 retop =  Topology[sym][ps];
	 If[ $DDimensions/.opt/.Options[CreateFeynAmp],
	     retop = retop/. symruleddim ,
	     retop = retop/. symrule4dim 
	   ];
	 retop = retop /. {ru } ;
	 TopToGrassman[ retop ]
	]

                     (** constructing the global prefactor ... **)
                     (** for truncated the -I is eliminated 
                         when calling the truncation rules     **)
symruleddim = 
 {
  Topology[sym_][ props__ ]:>
  Topology[ -I (1/(2 Pi)^$D Mu^(4-$D))^LoopNr[{props}] / sym ][ props ]
 }

symrule4dim = 
 {
  Topology[sym_][ props__ ]:>
  Topology[ -I (1/(2 Pi)^4)^LoopNr[{props}] / sym ][ props ]
 }
                     (** ... using Eulers relation **)

 LoopNr[ ps_List ] := 
    Module[ {all,three,four,extern} ,
          all = Union[ Flatten[ ps /. Propagator[_]->List ] ];
	  extern = Count[ all , e[_] ];
	  three = Count[ all , v3[_] ];
	  four = Count[ all , v4[_] ];
	  Return[ 1/2 (three-extern) + four + 1 ]
	 ]

           (** Vertex -> AnalyticalCoupling : automatic insertion **)

VertInsert[mod_][ Topology[s_][elem__] ] :=
 Topology[s][elem] /.{ Vertex[___,vv_List,___]:>
		       AnalyticalCoupling[mod]@@ Sort[vv] }

                           (** put propagator to nice form **)

ppOne[ h_[ nr_ , one_ , two_ , three_ , four_ ] ] := h[nr,one,three];
ppTwo[ h_[ nr_ , one_ , two_ , three_ , four_ ] ] := h[nr,two,four];
ppOne[ h_[ nr_ , one_ , two_ ] ] := h[nr,one];
ppTwo[ h_[ nr_ , one_ , two_ ] ] := h[nr,two];
ppOne[ x_ ] := x;
ppTwo[ x_ ] := x;

PropOrder[mod_][ Propagator[type_][_,_,part_,mom_] ] :=
    Module[{tt,from,to},
          from = ppOne[part] ;
          to   = ppTwo[part] ;
	  If[ Or[ SameQ[inc,type],SameQ[out,type] ] , 
                                      (** if the momentum of the outgoing
                                          particle is incoming (now it is
                                          an outgoing ki), then the 
                                          rule out->1 has to be changed into
                                          out->-1 !!                      **)
              tt = type/.{inc->1,out->1};
              If[ Head[ AnalyticalPropagator[mod][ex][to,mom] ] === PV   ,
                  AnalyticalPropagator[mod][ ex ][  to , 
                                                    Expand[tt*mom] ] ,
                  AnalyticalPropagator[mod][ ex ][  AntiParticle[mod][from] ,
						                                  Expand[-tt*mom] ]
		          ],
	      tt = in ;
(* The direction of the in-propagator is only relevant for the fermions.
   As these are formed against the flow of fermion number, part is always
   the antifermion. So to apply our Feynman rules to to and AP[from] we just
   reverse the propagator and change the momentum sign. 
   (For us the second rule would be enough, but someone might enter 
    Feynman rules from fermion to antifermion (or particle to antiparticle)) *)

	      If[ Head[ AnalyticalPropagator[mod][tt][AntiParticle[mod][from],
						    to,mom] ] === PV   ,
                  AnalyticalPropagator[mod][ tt ][ 
				   AntiParticle[mod][from] , to ,
				   mom                  ],
                  AnalyticalPropagator[mod][ tt ][ 
				   to , AntiParticle[mod][from] ,
				   Expand[-mom] ]                  
            ] 
        ] 
 	 ]

      (** Propagator -> AnalyticalPropagator : automatic insertion **)

PropInsert[mod_][ tt:Topology[s_][elem__] ] := tt /.( 
		   Propagator[ty_][ff__] :>
		   PropOrder[mod][Propagator[ty][ff]] )

               (** resolve powers of PropagatorDenominator, 
                   collect PropagatorDenominator factors, expand them,
                   order them in FeynAmpDenominators              **)

allnonpd[list_]:=Select[ list[[1]],(FreeQ[#,PropagatorDenominator]||
                               FreeQ[#,q])& ];
allpd[list_]:=( ((Expand[#,PropagatorDenominator]& [ Select[ list[[1]],
                                  !(FreeQ[#,PropagatorDenominator]||
                                    FreeQ[#,q]                      )& ]
                       ]) /. rule1) /.rule2 ) //.rule3 ;
rule1= Power[PropagatorDenominator[a___],n_Integer?Positive]  :>
              FeynAmpDenominator @@ Array[prden[a]&,n];
rule2 = PropagatorDenominator[a__] :> FeynAmpDenominator[ prden[a] ];
rule3 = FeynAmpDenominator[ a___]  FeynAmpDenominator[ b___ ] :>
              FeynAmpDenominator[ a , b ] ;
CollectPD[list_] := 
FeynAmp[ allnonpd[list] * ( allpd[list] /. 
                          prden -> PropagatorDenominator ) ] 

               (** expand the FeynAmpDenominator, 
		   PropagatorDenominator and FourVector arguments **)

lv[ mom_ , args___ ] := FourVector[ Expand[mom] , args ]

df[ mom_, mass_ ] := If[ MemberQ[ mom , -q[1] ] ||
			  Head[ First[ mom ] ]===Times ||
			  MemberQ[ {-q[1],-q[2],-q[3]},mom ] ,
			  PropagatorDenominator[ Expand[-mom] , mass ],
			  PropagatorDenominator[ Expand[mom]  , mass ] 
	     	       ]/;!(mom===0);
df[ 0 , mass_ ] := PropagatorDenominator[ 0, mass ];

dn[ allfactors___ ] := Sort[ FeynAmpDenominator[ allfactors ] ]

                        (** construct integration information **)

IntegrationInfo[ faexpr_ ] := Module[{searchq,tolist},
	 tolist = Join[Transpose[{FeynArtsSymbols,
			      Table[List,{Length[FeynArtsSymbols]}]}],
			    Transpose[{Results,
			      Table[List,{Length[Results]}]}]
			   ];
	 tolist =  Apply[Rule,#]& /@ tolist;
	 searchq = Flatten[(List@@(List/@faexpr))/.{
                            Plus->List,Times->List}/.tolist ];
	 searchq = Union[ Cases[ searchq , q[_] ] ];
	 If[ searchq === {} ,
	     Return[ MomentumIntegration[] ],
	     Return[ MomentumIntegration @@ searchq ]
	   ]
	 ]

				 (** substituting dummy by appropr. gi[n] **)

DummyTreat[ FeynAmp[amp__] ] := DummyTreat[ FeynAmp[amp],
				  Last[   Sort[ Cases[{amp},
						      gi[_],
						      Infinity] ]  ][[1]] ]

DummyTreat[ FeynAmp[amp__],high_Integer ] := 
     Module[ {newamp,thisdummypair,thisgipair},
     If[ FreeQ[ {amp},dummy ],
	 Return[ FeynAmp[amp] ],
         newamp=FeynAmp[amp] /. ( SU3F[a___, dummy, b___] *
				  SU3F[c___, dummy, d___] :>
		                  Prod[SU3F[a, dummy, b], SU3F[c, dummy, d]] );
         thisdummypair = Cases[newamp,Prod[___],Infinity][[1]];
         thisgipair = thisdummypair /. dummy->gi[high+1];
         DummyTreat[ (newamp/.thisdummypair:>thisgipair)/.Prod->Times , high+1 ]
       ] ]

DummyRule =  FeynAmp[a___] :> If[ FreeQ[{a},dummy],
				  FeynAmp[a],
				  DummyTreat[ FeynAmp[a] ]  ];
OuterIndexRule = {
   Propagator[inc][fr_,to_,h_[ a_,b_,c_,d_,e_ ],mom_ ] :>
    Propagator[inc][fr,to,h[a,c,c,e,e],mom],
   Propagator[inc][fr_,to_,h_[ a_,b_,c_],mom_ ] :>
    Propagator[inc][fr,to,h[a,c,c],mom],
   Propagator[out][fr_,to_,h_[ a_,b_,c_,d_,e_ ],mom_ ] :>
    Propagator[out][fr,to,h[a,c,c,e,e],mom],
   Propagator[out][fr_,to_,h_[ a_,b_,c_ ],mom_ ] :>
    Propagator[out][fr,to,h[a,c,c],mom]
}

                                 (** putting it together **)

CreateAmplitude[ mod_ ][ acTop:Topology[ sym_ ][ ps__ ] , opt_ ] :=
       Module[
             { theamp,xirule,dummy1,dummy2 },
             If[ (FeynmanGauge/.opt) , 
		       xirule=GaugeXi[_]->1, xirule=dummy1->dummy2 ];
	          theamp = acTop/.OuterIndexRule;
             theamp = VertexHandling[ theamp ];
	          theamp = VertInsert[mod][ theamp ];
             theamp = PropInsert[mod][ theamp ] /.xirule ;
             theamp = FeynMultiplication[ theamp ];
             theamp = CollectPD[ theamp ];
	          theamp = theamp/.{FourVector->lv,
			             PropagatorDenominator->df,
			             FeynAmpDenominator->dn};
             theamp = theamp /. ColorFactorRule ;
	     If[ Truncated/.opt/.Options[CreateFeynAmp] ,
                             (** eliminate factor -I from global prefactor  **)
          theamp = theamp/.( Last[theamp]:>(I Last[theamp]) ) ;
   		 theamp = theamp/.TruncationRules 
	       ];
	     theamp = Prepend[ theamp, IntegrationInfo@@theamp ] ; 
             theamp = theamp /. DummyRule    ;
             Off[ General::spell1 ];
             theamp =  theamp//.FeynCalcRules ;   
             On[ General::spell1 ];
	     Return[theamp]
            ]

End[]

EndPackage[]
