Page MenuHomeHEPForge

No OneTemporary

diff --git a/current_generator/include/helspin.frm b/current_generator/include/helspin.frm
index 26d02fd..921f34e 100644
--- a/current_generator/include/helspin.frm
+++ b/current_generator/include/helspin.frm
@@ -1,193 +1,263 @@
*/**
* \brief Procedures for dealing with helicity spinors
*
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
on shortstats;
-cf Conjugate,Current;
-cf angle,square;
+* Complex conjugate
+cfunction Conjugate;
+* Generalised current, see eq:j_gen in the developer manual
+cfunction Current;
+* Angle and square products
+* See eq:angle_product and eq:square_product in the developer manual
+cfunction angle,square;
+* Generalised angle and square products and
+* generalised negative-helicity currents
+* See eq:spinors_soinhel and eq:current_spinhel in the developer manual
#define CHAINS "AngleChain,SquareChain,SpinorChain"
-ctensor `CHAINS';
+ctensors `CHAINS';
* internal symbols are all uppercase
-v P,Q,K,L;
-cf DUM;
-ctensor CT;
-s X,Y;
-i MU1,...,MU50;
+vectors P,Q,K,L;
+cfunctions DUM;
+ctensors CT;
+symbols X,Y;
+indices MU1,...,MU50;
+* Compute all current contractions
+*
+* Example:
+*
+* vectors p1,...,p4;
+* indices mu;
+* local foo = Current(-1, p1,mu,p2)*Current(-1, p3,mu,p4);
+* #call ContractCurrents
+* print;
+* .end
+*
+* Example Output:
+*
+* foo = - 2*angle(p1,p3)*square(p2,p4);
+*
#procedure ContractCurrents
#call InsCurrent;
.sort
#call GetMomentaInChains($MOMENTA)
#if "`$MOMENTA'" != ""
#message Rewrite momenta `$MOMENTA' into spinors
#call ToSpinors(`$MOMENTA')
#endif
while(match(SpinorChain(P?, MU1?, Q?)*SpinorChain(K?, MU1?, L?)));
#call Fierz(MU1?)
endwhile;
#call SortArgs
id AngleChain(P?, Q?) = angle(P, Q);
id SquareChain(P?, Q?) = square(P, Q);
argument denom_;
id AngleChain(P?, Q?) = angle(P, Q);
id SquareChain(P?, Q?) = square(P, Q);
endargument;
#endprocedure
+* INTERNAL PROCEDURE
+* Find all momenta that appear inside
+* generalised currents, angle products, and square products
#procedure GetMomentaInChains(MOMENTA)
id CT?{`CHAINS'}(?a) = DUM(?a)*CT(?a);
#$TMP = 0;
#$P = 0;
argument DUM;
dropcoefficient;
$P = term_;
$TMP = $TMP + $P;
endargument;
id DUM(?a) = 1;
moduleoption local $P;
moduleoption sum $TMP;
.sort
#$TMP = dum_($TMP);
#inside $TMP
splitarg dum_;
argument dum_;
dropcoefficient;
endargument;
repeat id dum_(?a, MU1?index_, ?b) = dum_(?a, ?b);
id dum_(?a`MOMENTA') = 0;
#endinside
#endprocedure
+* INTERNAL PROCEDURE
+* Rewrite Current as SpinorChain
+* See eq:jplus_gen, eq:current_spinhel in the developer manual
#procedure InsCurrent
id Current(+1, ?a) = Current(-1, reverse_(?a));
* hack
* with ?a as pattern FORM has problems if the arguments are sums of momenta
#do N=3,50
id Current(-1, MU1?,...,MU`N'?) = SpinorChain(MU1,...,MU`N');
#enddo
#endprocedure
+* INTERNAL PROCEDURE
+* Sort arguments of angle and square products
+* See eq:angle_product, eq:square_product in the developer manual
#procedure SortArgs
#do F={AngleChain,SquareChain}
antisymmetrize `F':2;
#enddo
#endprocedure
+* INTERNAL PROCEDURE
+* Simplify expressions containing complex conjugates
+* See eq:square_product, eq:jplus_gen in the developer manual
#procedure DoConjugate
splitarg Conjugate;
repeat id Conjugate(X?, Y?, ?a) = Conjugate(X) + Conjugate(Y, ?a);
normalize Conjugate;
id Conjugate(AngleChain(P?, Q?)) = -SquareChain(P, Q);
id Conjugate(SquareChain(P?, Q?)) = -AngleChain(P, Q);
id Conjugate(SpinorChain(?a)) = DUM(reverse_(?a));
id DUM(?a) = SpinorChain(?a);
#endprocedure
+* INTERNAL PROCEDURE
+* Apply a Fierz transform to currents contracted over index `MU'
+* See eq:Fierz in the developer manual
#procedure Fierz(MU)
once SpinorChain(P?, `MU', Q?)*SpinorChain(K?, `MU', L?) = (
2*AngleChain(P, K)*SquareChain(L, Q)
);
#endprocedure
+* INTERNAL PROCEDURE
+* Rewrite a single momentum `P' inside a SpinorChain into spinors
+* See eq:p_in_current in the developer manual
+* See ToSpinors for an example
#procedure MomentumToSpinorsInSpinorChain(P)
id SpinorChain(P?, ?a, `P', ?b, Q?) = DUM(
mod_(nargs_(P, ?a),2), mod_(nargs_(?b, Q),2),
SpinorChain(P, ?a, `P', ?b, Q)
);
id DUM(1, 1, SpinorChain(P?, ?a, `P', ?b, Q?)) = (
AngleChain(P, ?a, `P')*SquareChain(`P', ?b, Q)
);
id DUM(0, 0, SpinorChain(P?, ?a, `P', ?b, Q?)) = (
SpinorChain(P, ?a, `P')*SpinorChain(`P', ?b, Q)
);
if(count(DUM,1)>0);
print "even number of arguments in spinor chain in %t";
exit;
endif;
#endprocedure
+* INTERNAL PROCEDURE
+* Rewrite a single momentum `P' inside an AngleChain into spinors
+* See eq:p_in_angle in the developer manual
+* See ToSpinors for an example
#procedure MomentumToSpinorsInAngleChain(P)
id AngleChain(P?, ?a, `P', ?b, Q?) = DUM(
mod_(nargs_(P, ?a),2), mod_(nargs_(?b, Q),2),
AngleChain(P, ?a, `P', ?b, Q)
);
id DUM(1, 0, AngleChain(P?, ?a, `P', ?b, Q?)) = (
AngleChain(P, ?a, `P')*Conjugate(SpinorChain(`P', ?b, Q))
);
id DUM(0, 1, AngleChain(P?, ?a, `P', ?b, Q?)) = (
SpinorChain(P, ?a, `P')*AngleChain(`P', ?b, Q)
);
if(count(DUM,1)>0);
print "odd number of arguments in angle chain in %t";
exit;
endif;
#endprocedure
+* INTERNAL PROCEDURE
+* Rewrite a single momentum `P' inside a SquareChain into spinors
+* See eq:p_in_square in the developer manual
+* See ToSpinors for an example
#procedure MomentumToSpinorsInSquareChain(P)
id SquareChain(P?, ?a, `P', ?b, Q?) = DUM(
mod_(nargs_(P, ?a),2), mod_(nargs_(?b, Q),2),
SquareChain(P, ?a, `P', ?b, Q)
);
id DUM(1, 0, SquareChain(P?, ?a, `P', ?b, Q?)) = (
SquareChain(P, ?a, `P')*SpinorChain(`P', ?b, Q)
);
id DUM(0, 1, SquareChain(P?, ?a, `P', ?b, Q?)) = (
Conjugate(SpinorChain(P, ?a, `P'))*SquareChain(`P', ?b, Q)
);
if(count(DUM,1)>0);
print "odd number of arguments in angle chain in %t";
exit;
endif;
#endprocedure
+* INTERNAL PROCEDURE
+* Rewrite a single momentum `P' into spinors
+* See ToSpinors for an example
#procedure MomentumToSpinors(P)
repeat;
#do CHAIN={`CHAINS'}
id `CHAIN'(?a, `P', `P', ?b) = 0;
#call MomentumToSpinorsIn`CHAIN'(`P')
#enddo
endrepeat;
#endprocedure
+* INTERNAL PROCEDURE
+* Successively rewrite all momenta `?PS' into spinors,
+* i.e. replace p = <p| [p| + |p> |p]
+*
+* Example:
+*
+* v p1,p2,p3;
+* l foo = SpinorChain(p1,p2,p3,p2,p1);
+* #call ToSpinors(p2,p3)
+* print;
+* .end
+*
+* Example Output:
+*
+* foo =
+* AngleChain(p1,p2)*AngleChain(p3,p2)*SquareChain(p2,p1)*SquareChain(p2,p3)
+*
#procedure ToSpinors(?PS)
#do P={`?PS'}
#call MomentumToSpinors(`P')
argument;
#call MomentumToSpinors(`P')
endargument;
factarg Conjugate;
chainout Conjugate;
id Conjugate(Conjugate(X?)) = X;
#enddo
#call DoConjugate
#endprocedure
diff --git a/current_generator/include/write.frm b/current_generator/include/write.frm
index 9ef23c2..5741f5f 100644
--- a/current_generator/include/write.frm
+++ b/current_generator/include/write.frm
@@ -1,123 +1,129 @@
+* Write start of C++ header to `OUTPUT'
#procedure WriteHeader(OUTPUT)
#write <`OUTPUT'> "#pragma once"
#write <`OUTPUT'> ""
#write <`OUTPUT'> "#include <complex>"
#write <`OUTPUT'> "#include \"HEJ/helicity.hh\""
#write <`OUTPUT'> "#include \"HEJ/LorentzVector.hh\""
#write <`OUTPUT'> ""
#write <`OUTPUT'> "namespace HEJ {"
#endprocedure
+* Write optimised expression to C++ header `OUTPUT'
#procedure WriteOptimised(OUTPUT,EXPRESSION,NUMHELICITIES,?MOMENTA)
#write <`OUTPUT'> " template<%"
#define FIRST "1"
#do i=1,`NUMHELICITIES'
#ifdef `FIRST'
#undefine FIRST
#write <`OUTPUT'> "Helicity%"
#else
#write <`OUTPUT'> ", Helicity%"
#endif
#enddo
#write <`OUTPUT'> ">"
#write <`OUTPUT'> " std::complex<double> `EXPRESSION'("
#call WriteMomenta(`?MOMENTA')
#write <`OUTPUT'> "\n );\n"
#call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',`NUMHELICITIES',`?MOMENTA')
#endprocedure
+*INTERNAL PROCEDURE
#procedure WriteOptimisedHelper(OUTPUT,EXPRESSION,NUMHELICITIES,?REST)
#if `NUMHELICITIES' > 0
#do H={+,-}
#call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',{`NUMHELICITIES'-1},`?REST',`H')
#enddo
#else
#define HELSTRING ""
#define TEMPLATEARGS ""
#define MOMENTA ""
#do ARG={`?REST'}
#if ("`ARG'" == "+") || ("`ARG'" == "-")
* arguments that define helicities
#redefine HELSTRING "`HELSTRING'`ARG'"
#if "`TEMPLATEARGS'" != ""
#redefine TEMPLATEARGS "`TEMPLATEARGS',"
#endif
#if "`ARG'" == "+"
#redefine TEMPLATEARGS "`TEMPLATEARGS'Helicity::plus"
#else
#redefine TEMPLATEARGS "`TEMPLATEARGS'Helicity::minus"
#endif
#else
* arguments that define momenta
#if "`MOMENTA'" != ""
#redefine MOMENTA "`MOMENTA',"
#endif
#redefine MOMENTA "`MOMENTA'`ARG'"
#endif
#enddo
#optimize [`EXPRESSION'`HELSTRING']
#write <`OUTPUT'> " template<>"
#write <`OUTPUT'> " inline std::complex<double> `EXPRESSION'<%"
* we use a loop here because otherwise FORM will insert line breaks
* if the string is too large
#define FIRST "1"
#do TEMPLATEARG={`TEMPLATEARGS'}
#ifdef `FIRST'
#undefine FIRST
#else
#write <`OUTPUT'> ", %"
#endif
#write <`OUTPUT'> "`TEMPLATEARG'%"
#enddo
#write <`OUTPUT'> ">("
#call WriteMomenta(`MOMENTA')
#write <`OUTPUT'> "\n ) {"
#if `optimmaxvar_' > 0
#write <`OUTPUT'> " std::complex<double> %"
#define FIRST "1"
#do i=1,`optimmaxvar_'
#ifdef `FIRST'
#undefine FIRST
#write <`OUTPUT'> "Z`i'_%"
#else
#write <`OUTPUT'> ", Z`i'_%"
#endif
#enddo
#write <`OUTPUT'> ";"
#endif
#write <`OUTPUT'> " %O"
#write <`OUTPUT'> " return %E;" [`EXPRESSION'`HELSTRING']
#write <`OUTPUT'> " }\n" [`EXPRESSION'`HELSTRING']
#clearoptimize
#endif
#endprocedure
+*INTERNAL PROCEDURE
+* Write momenta as C++ function arguments to `OUTPUT'
#procedure WriteMomenta(?MOMENTA)
#define FIRST "1"
#do P={`?MOMENTA'}
#ifdef `FIRST'
#undefine FIRST
#write <`OUTPUT'> " CLHEP::HepLorentzVector const & `P'%"
#else
#write <`OUTPUT'> ",\n CLHEP::HepLorentzVector const & `P'%"
#endif
#enddo
#endprocedure
+* Write end of C++ header to `OUTPUT'
#procedure WriteFooter(OUTPUT)
#write <`OUTPUT'> "}"
#endprocedure

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:58 PM (5 h, 40 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486755
Default Alt Text
(11 KB)

Event Timeline