(**
                 FeynInit.m

(basic initialization routine)
**)

BeginPackage["FeynInit`","FeynArts`"]

FeynInit::usage = "
FeynInit[ Model->MOD, modeloptions ] reads the file ./Model/model.MOD \n
and creates the functions PossibleParticles, AntiParticle, ParticleList, \n
VertexList, AnalyticalCoupling and AnalyticalPropagator." 

Options[ FeynInit ] := Options[InsertFields`InsertFields]

Begin["`Private`"]

VertexList[ i___Integer,mod_List] := VL @@ 
		  ( Union @@ Join[ VertexList[i,#]& /@ mod ] )

AntiParticle[ mod_List ][ f_[a_,b___] ] := 
     Join[ First[  Select[ (AntiParticle[#][f[a]]& /@ mod) , 
	           FreeQ[#,AntiParticle]& ]  ] ,
	   f[b]  
	   ] /; Not[f===fi] 

AnalyticalCoupling[ mod_List ][ c__List ]:=
 First[  Select[ (AnalyticalCoupling[#][c]& /@ mod) , 
		 FreeQ[#,AnalyticalCoupling]& ]  ] 

FirstExtended[ {} ] = {};
FirstExtended[ x_ ] := First[x];

AnalyticalPropagator[ mod_List ][ type_ ][ c___ ]:=
 FirstExtended[  Select[ (AnalyticalPropagator[#][type][c]& /@ mod) , 
		 FreeQ[#,AnalyticalPropagator]& ]  ] 

FeynInit::initialized = "
the model `1` is already initialized ! (no initialization) "

FeynInit[ options___Rule ] :=
	 FeynInit[ options , ModelNameSymbol@@ localoptions ];

FeynInit[ options___Rule , modelname_ ] :=
   Module[ { localoptions,holdoptions,model,mess },
	 holdoptions = Options[FeynInit] ;
	 localoptions = SetOptions[FeynInit,options] ;
	 Options[ FeynInit ] = holdoptions ;
	 model = SymbolToModel[ modelname ] ; 
	 If[ MemberQ[ InitializedModels , modelname ],
	     Message[ FeynInit::initialized , modelname ];
	     Return[],
	     If[ ScreenMessages /. localoptions ,
		 $Output = {"stdout"} ,
		 $Output = {FeynArtsNullDevice}      ];
	     mess = Init[ model ,  modelname  , localoptions ];
	     $Output = {"stdout"};
             If[ mess===$Aborted,
                 mess,
	         InitializedModels = Append[ InitializedModels,
                                            modelname ]
              ]
	   ]
	] 

Init::filenotfound = "
Could not find file \"model.`1`\"  (no initialization). \n
Allowed Models are `2`. " 

Init[mod_, name_, opt_ ]:=
    Module[ {modelfile},
	  Unprotect[QMM];
          If[ QuarkMixing /. opt ,
	      QMM = QuarkMixingMatrix ,
	      QMM = IdentityMatrix[3] 
	    ];
	  Protect[QMM];
	  modelfile = Select[ $ModelFiles, 
			      StringMatchQ[#,StringJoin["*.",ToString[mod]]]&];
	  If[ modelfile === {},
	      Message[Init::filenotfound, mod , Models];
	      Return[ $Aborted ],
              Off[Syntax::newl];
	      Get @@ modelfile;
              On[Syntax::newl]
	     ];
	  Set[ PropList[ name ] , PropList[ mod ] ];
	  Set[ CoupList[ name ] , CoupList[ mod ] ];
	  CoupList[ name ] = Select[ CoupList[name],
				      Not[ SameQ[ #[[2]],PV[0] ] ]& ];
          If[ UnitaryGauge /. opt ,
              PropList[ name ] = DelForUnitGauge[ PropList[name] ];
              PropList[ name ] =  PropList[name] /. 
                  {( (1-1/GaugeXi[A]) PropagatorDenominator[ _,0 ])->0,
                   ( (1-1/GaugeXi[_]) PropagatorDenominator[ 
                           _,m_ /Sqrt[GaugeXi[_]] ])-> 1/ m^2  } ;
              CoupList[ name ] = DelForUnitGauge[ CoupList[name] ];
              Print[" ... Unitary Gauge: tree graphs only ! "]
             ];
          If[ Not[ LightFHCoupling/.opt ] ,
	      CoupList[ name ] = DeleteLFHC[ CoupList[ name ] ];
	      Print[" ... LightFermionHiggsCouplings deleted "]
	    ];
          If[ Not[ ElectronHCoupling/.opt ] ,
	      CoupList[ name ] = DeleteEHC[ CoupList[ name ] ];
	      Print[" ... ElectronHiggsCouplings deleted "]
	    ];
          If[ Not[ Generation1/.opt ] ,
	      PropList[ name ] = DeleteG1[ PropList[ name ] ];
	      CoupList[ name ] = DeleteG1[ CoupList[ name ] ];
	      Print[" ... Generation1 deleted "]
	    ];
          If[ Not[ Generation2/.opt ] ,
	      PropList[ name ] = DeleteG2[ PropList[ name ] ];
	      CoupList[ name ] = DeleteG2[ CoupList[ name ] ];
	      Print[" ... Generation2 deleted "]
	    ];
          If[ Not[ Generation3/.opt ] ,
	      PropList[ name ] = DeleteG3[ PropList[ name ] ];
	      CoupList[ name ] = DeleteG3[ CoupList[ name ] ];
	      Print[" ... Generation3 deleted "]
	    ];
	  CreateAntiParticles[ name ];
     CreateVertexList[ name ];
     CreateParticleList[ name ];
	  CreateMasses[ name ];
	  CreatePossibleParticles[ name ];
	  Unprotect[PV];
	  CreateAnalyticalCoupling[ name ];
	  CreateAnalyticalPropagator[ name ];
	  Protect[PV];
	  Null
         ]
			      (** deleting parts of the model **)

FreeList[ x_ , l_List ] := 
	Module[ {freelist,memberlist={}},
	       Do[ memberlist = Append[ memberlist , 
					Select[ x , 
						Not[ FreeQ[#,l[[i]]]]&
					      ]
				      ] ,
		   {i,Length[l]}
		  ];
		memberlist = Union@@ memberlist ;
               freelist = Select[ x , Not[ MemberQ[ memberlist,# ] ]& ];
	       Return[ freelist ]
	     ]

MemberList[ x_ , l_List ] := 
	Module[ {memberlist={}},
	       Do[ memberlist = Append[ memberlist , 
					Select[ x , 
						Not[ FreeQ[#,l[[i]]]]&
					      ]
				      ] ,
		   {i,Length[l]}
		  ];
		memberlist = Union@@ memberlist ;
		Return[memberlist]
	       ]

DeleteLFHC[ x_ ] := Module[{del,hold},
			  hold = FreeList[ x , {S[1],S[2],S[3],
						S[-1],S[-2],S[-3]} ];
			  del = MemberList[ x , {S[1],S[2],S[3],
						 S[-1].S[-2],S[-3]} ];
			  del = Union[FreeList[ del, {F[1],F[2],F[3],F[4],
						F[5],F[6],F[7],F[8],
						F[9],F[10],F[12],
			                        F[-1],F[-2],F[-3],F[-4],
						F[-5],F[-6],F[-7],F[-8],
						F[-9],F[-10],F[-12]} ],
				      MemberList[ del, {F[11],F[-11]}]];
			  Global`un= Union[ hold , del ] ;
			  Union[ hold , del ] /.
			  {ChiralityProjector[ 1] Global`MB :> 0,
			   ChiralityProjector[-1] Global`MB :> 0,
			   ChiralityProjector[ 1] Global`MS :> 0,
			   ChiralityProjector[-1] Global`MS :> 0,
			   ChiralityProjector[ 1] Global`MD :> 0,
			   ChiralityProjector[-1] Global`MD :> 0}
			]

DeleteEHC[ x_ ] := Module[{del,hold},
			  hold = FreeList[ x , {S[1],S[2],S[3],
						S[-1],S[-2],S[-3]} ];
			  del = MemberList[ x , {S[1],S[2],S[3],
        					 S[-1].S[-2],S[-3]} ];
                          del = DeleteElecPosi[ del ];
                          Union[ hold, del ]
			]

LeaveQEDonly[ x_ ] := 
 FreeList[ x, {F[-12],F[-11],F[-10],F[-8],F[-7],F[-6],F[-4],F[-3],F[-2],
               F[2],F[3],F[4],F[6],F[7],F[8],F[10].F[11],F[12],
               S[-3],S[1],S[2],S[3],
               V[-3,_],V[2,_],V[3,_], 
               U[-4],U[-3],U[-2],U[-1],U[1],U[2],U[3],U[4]} ] 

DeleteElecPosi[ x_ ] := FreeList[ x , {F[-1],F[1]} ]

DelForUnitGauge[ x_ ] := 
FreeList[ x , {U[-4],U[-3],U[-2],U[-1],U[1],U[2],U[3],U[4],
	       S[-3],S[2],S[3]} ]

DeleteG1[ x_ ] :=
FreeList[ x , {F[-4],F[-3],F[-2],F[-1],F[1],F[2],F[3],F[4]} ]

DeleteG2[ x_ ] := 
FreeList[ x , {F[-8],F[-7],F[-6],F[-6],F[5],F[6],F[7],F[8]} ]

DeleteG3[ x_ ] := 
FreeList[ x , {F[-12],F[-11],F[-10],F[-9],F[9],F[10],F[11],F[12]} ]

                              (** creating AntiParticle-function **)

MakeAP[ mod_ ][ Global`Prop[_][ fa_,fb_, ___ ] ] := Null 

MakeAP[ mod_ ][ Global`Prop[in][ fa_,fb_, ___ ] ] := 
        { Set[ AntiParticle[mod][Take[fa,1]] , Take[fb,1] ],
          Set[ AntiParticle[mod][Take[fb,1]] , Take[fa,1] ]   }

CreateAntiParticles[ mod_ ] :=
             MakeAP[mod][ #1[[1]] ]& /@ PropList[mod]

                                      (** creating VertexList **)
        
MakeVertexList[ Global`Coup[{fi1_,mom1_},{fi2_,mom2_},
			    {fi3_,mom3_}]              ] :=
            Vertex[ Take[fi1,1],Take[fi2,1],Take[fi3,1] ]

MakeVertexList[ Global`Coup[{fi1_,mom1_},{fi2_,mom2_},
		            {fi3_,mom3_},{fi4_,mom4_}] ] :=
            Vertex[ Take[fi1,1],Take[fi2,1],Take[fi3,1],Take[fi4,1] ]

CreateVertexList[ mod_ ] := 
      Module[{},
	    SetDelayed[ VertexList[i_Integer][mod],
			Set[ VertexList[i,mod], 
			     Select[ VertexList[mod],Length[#]===i& ]
			   ]
		       ];
            Set[ VertexList[mod] , 
                 FeynArts`VL@@ (
		 MakeVertexList[#[[1]]]& /@ CoupList[mod]) ]
	   ]

				     (** creating ParticleList **)

ExtractParticle[ Global`Prop[_][ fa_,fb_, ___ ] ] := {}

ExtractParticle[ Global`Prop[in][ fa_,fb_, ___ ] ] := 
			{Take[fa,1] , Take[fb,1]}

CreateParticleList[ mod_ ] := 
   Set[ ParticleList[mod] , 
        FeynArts`PL@@( Union[Flatten[ExtractParticle[#[[1]]]& /@ 
		       PropList[mod]]] )                  ]

            (** auxiliary functions for CreatePossibleParticles **)

c31[ part_,mod_ ] := 
     Set[ Ok[mod][ 
          Vertex[ part , zzz[1] , zzz[1] ]  ],
          True                                             ]

c41[ part_,mod_ ] := 
     Set[ Ok[mod][
          Vertex[part,zzz[1],zzz[1],zzz[1]]  ],
          True                                             ]

c32[ v_[part1_,part2_],mod_ ] := 
     Set[ Ok[mod][
          Vertex[part1,part2,zzz[1]         ] ],
          True                                             ]

c42[ v_[part1_,part2_],mod_ ] := 
     Set[ Ok[mod][
          Vertex[part1,part2,zzz[1],zzz[1]] ],
          True                                             ]

c33[ v_[part1_,part2_,part3_],mod_ ] := 
     Set[ Ok[mod][
          Vertex[part1,part2,part3] ],
          True                                             ]

c43[ v_[part1_,part2_,part3_],mod_ ] := 
     Set[ Ok[mod][
          Vertex[part1,part2,part3,zzz[1]] ],
          True                                             ]

c44[ v_[part1_,part2_,part3_,part4_],mod_ ] := 
     Set[ Ok[mod][
          Vertex[part1,part2,part3,part4] ],
          True                                             ]

pselect[mod_][vv_] := Select[ ParticleList[mod],
			      ( Ok[mod][Append[vv,#]]==True) & ]

o31[ part_,mod_ ] := 
     Set[ OkP[mod][
          Vertex[ part , zzz[1] , zzz[1] ]  ],
          pselect[mod][Vertex[part,zzz[1]]]                     ]

o41[ part_,mod_ ] := 
     Set[ OkP[mod][
          Vertex[part,zzz[1],zzz[1],zzz[1]]  ],
          pselect[mod][Vertex[part,zzz[1],zzz[1]]]               ]

o32[ v_[part1_,part2_],mod_ ] := 
     Set[ OkP[mod][
          Vertex[part1,part2,zzz[1]]      ],
          pselect[mod][Vertex[part1,part2]]             ]

o42[ v_[part1_,part2_],mod_ ] := 
     Set[ OkP[mod][
          Vertex[part1,part2,zzz[1],zzz[1]] ],
          pselect[mod][Vertex[part1,part2,zzz[1]]]                     ]

o43[ v_[part1_,part2_,part3_],mod_ ] := 
     Set[ OkP[mod][
          Vertex[part1,part2,part3,zzz[1]]  ],
          pselect[mod][Vertex[part1,part2,part3]]                     ]


h32[vx_]:=FeynArts`VL[ Take[vx,2],Take[vx,-2],Drop[vx,{2}] ]
h42[vx_]:=FeynArts`VL[ Take[vx,2],Take[vx,-2],Take[vx,{2,3}],
              {vx[[1]],vx[[4]]},{vx[[1]],vx[[3]]},{vx[[2]],vx[[4]]} ]
h43[vx_]:=FeynArts`VL[ Drop[vx,{1}],Drop[vx,{2}],
                  Drop[vx,{3}],Drop[vx,{4}] ]

                               (** creating PossibleParticles **)

CreatePossibleParticles[mod_] :=
 Module[{modvert3,modvert4,part31,part41,part32,part42,part43},
       SetDelayed[ FeynArts`Ok[FeynArts`Vertex[__],mod] , False ];
       modvert3=VertexList[3][mod];
       modvert4=VertexList[4][mod];
       Map[ c33[#,mod]&,modvert3 ];
       Map[ c44[#,mod]&,modvert4 ];
       part31=Union[  Flatten[ Apply[Vertex,modvert3] ]  ];
       part41=Union[  Flatten[ Apply[Vertex,modvert4] ]  ];
       Map[ c31[#,mod]&,part31 ];
       Map[ c41[#,mod]&,part41 ];
       part32=Union[  Flatten[ Map[h32,modvert3] ]  ];
       part42=Union[  Flatten[ Map[h42,modvert4] ]  ];
       part43=Union[  Flatten[ Map[h43,modvert4] ]  ];
       Map[ c32[#,mod]&,part32 ];
       Map[ c42[#,mod]&,part42 ];
       Map[ c43[#,mod]&,part43 ];
       SetDelayed[ OkP[mod][_],FeynArts`PL[] ];
       Set[  OkP[mod][Vertex[zzz[1],zzz[1],zzz[1]]],       
	     Apply[ FeynArts`PL,Apply[Union,VertexList[3][mod]] ]  ]; 
       Set[  OkP[mod][Vertex[zzz[1],zzz[1],zzz[1],zzz[1]]], 
	     Apply[ FeynArts`PL,Apply[Union,VertexList[4][mod]] ]  ]; 
       Print[ " |.    |"  ];
       Map[ o31[#,mod]&,part31 ];
       Print[ " |..   |"  ];
       Map[ o41[#,mod]&,part41 ];
       Print[ " |...  |"  ];
       Map[ o32[#,mod]&,part32 ];
       Print[ " |.... |"  ];
       Map[ o42[#,mod]&,part42 ];
       Print[ " |.....|"  ];
       Map[ o43[#,mod]&,part43 ];
       Print[" ... (done)"];
     ]

SetDelayed[ PossibleParticles[ mod_ ][ vx_ ] ,
            OkP[ mod ][ Apply[ Vertex,vx ]/.{
	                v_[fi[i_]]->zzz[1],fi[i_]->zzz[1]}    ] ];

                         (** create AnalyticalCoupling - functions **)

CreateAnalyticalCoupling[ mod_ ] := 
    Module[ { cl } ,
	   cl = CoupList[ mod ]/.(
                    Global`Coup[ sth__ ] :> Global`Coup@@ Sort[{sth}] );
	   cl = cl /. Global`Coup->AnalyticalCoupling[mod] ;
	   cl /. Equal->SetDelayed;
	   Null
        ]
                       (** create AnalyticalPropagator - functions **)
                                      
CreateAnalyticalPropagator[ mod_ ] :=
    Module[ { pl } ,
	   pl=PropList[mod] /. Global`Prop->AnalyticalPropagator[mod] ;
	   pl /. Equal->SetDelayed;
	   Null
         ]

CreateMasses[ mod_ ] :=
    Module[ { pl } ,
      pl=Select[ PropList[mod],!FreeQ[#,in]& ] /.
                PV[a___ * PropagatorDenominator[mom_,mas_]]:> mas  ;
	   pl=Join[ pl /.  Global`Prop[in][f1_[a1_,a2___],f2_[a3_,a4___],___]->
                      Mass[f1[a1]]  ,
	            pl /.  Global`Prop[in][f1_[a1_,a2___],f2_[a3_,a4___],___]->
                      Mass[f2[a3]]  ];
	   pl /. Equal->SetDelayed;
	   Null
         ]

End[]

EndPackage[]

