CHROMA
heavy_hadron_potentials_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Potential between 2 heavy hadrons : Detmold
3  * Correlators checked independentely by Savage
4  *
5 
6  * Includes Lambda_b's etc
7 
8  * Soon also will have the BBB "potential"
9 
10  */
11 
12 #include "util/ferm/antisymtensor.h"
13 #include "barQll_w.h"
14 #include "mesQl_w.h"
15 #include "heavy_hadron_potentials_w.h"
16 
17 namespace Chroma
18 {
19  //! Heavy hadron potentials for SU(2) isospin limit
20  /*!
21  * \ingroup hadron
22  *
23  * This routine is specific to Wilson fermions!
24  *
25 
26  * \param u gauge field (Read)
27  * \param quark1 quark propagator, src at 0 ( Read )
28  * \param quark2 quark propagator, src at R ( Read )
29  * \param src1 cartesian coordinates of one source "0"( Read )
30  * \param src2 cartesian coordinates of the other source "R" ( Read )
31  * \param phases object holds list of momenta and Fourier phases ( Read )
32  * \param xml xml file object ( Read )
33  * \param xml_group group name for xml data ( Read )
34  *
35  */
36 
37 
38  void QllQllPOT(const multi1d<LatticeColorMatrix>& u,
39  const LatticePropagator& quark1,
40  const LatticePropagator& quark2,
41  const multi1d<int>& src1,
42  const multi1d<int>& src2,
43  const SftMom& phases,
44  XMLWriter& xml,
45  const std::string& xml_group)
46  {
47  START_CODE();
48 
49  if ( Ns != 4 ) /* Code is specific to Ns=4 */
50  return;
51 
52  if (src1[Nd-1] != src2[Nd-1] ){
53  QDPIO::cerr<<"Sources must be at the same time slice\n";
54  return ;
55  }
56 
57  QDPIO::cout << "Sources at : " << src1[0] << " " << src1[1]
58  << " " << src1[2] << " " << src1[3] << std::endl
59  << " and : " << src2[0] << " " << src2[1]
60  << " " << src2[2] << " " << src2[3] << std::endl;
61 
62  int length = phases.numSubsets() ;
63  int Nt = length;
64  int num_mom = phases.numMom();
65 
66  //Spin matrices
67  SpinMatrix g_one = 1.0;
68  SpinMatrix Cg5 = g_one * Gamma(5);
69  SpinMatrix Cg3 = - g_one * Gamma(14);
70  SpinMatrix Cgminus = g_one * Gamma(11) + timesI(g_one * Gamma(8));
71  SpinMatrix Cgplus = g_one * Gamma(11) + timesMinusI(g_one * Gamma(8));
72  SpinMatrix G5 = g_one * Gamma(15);
73  SpinMatrix OnePlusSigma3 = g_one + timesMinusI(g_one*Gamma(3));
74  SpinMatrix OneMinusSigma3 = g_one + timesI(g_one*Gamma(3));
75  SpinMatrix Gup = g_one * Gamma(15) + timesMinusI((g_one*Gamma(3))*Gamma(15));
76  SpinMatrix Gdown = g_one * Gamma(15) + timesI((g_one*Gamma(3))*Gamma(15));
77  SpinMatrix PosEnergyProj = 0.5*((g_one + Gamma(8) * g_one));
78  SpinMatrix NegEnergyProj = 0.5*((g_one - Gamma(8) * g_one));
79  SpinMatrix OnePlusSigma3PEP = OnePlusSigma3 * PosEnergyProj;
80  SpinMatrix OneMinusSigma3PEP = OneMinusSigma3 * PosEnergyProj;
81  SpinMatrix G5PEP = G5 * PosEnergyProj;
82 
83 
84  // Make heavy quark propagators
85  LatticeColorMatrix Qprop1;
86  LatticeColorMatrix Qprop2;
87  HeavyQuarkProp(Qprop1,u,src1,length);
88  HeavyQuarkProp(Qprop2,u,src2,length);
89 
90  //peek propagators: reduce lattice problem to two site problem
91  multi1d<DPropagator> U1,U2,U3,U4;
92  multi1d<DPropagator> D1,D2,D3,D4;
93  multi1d<DPropagator> antiU1,antiU2,antiU3,antiU4;
94  multi1d<DPropagator> antiD1,antiD2,antiD3,antiD4;
95  multi1d<ColorMatrix> Q1, Q2, antiQ1, antiQ2;
96  Q1.resize(length);
97  Q2.resize(length);
98  antiQ1.resize(length);
99  antiQ2.resize(length);
100  U1.resize(length); U2.resize(length); U3.resize(length);U4.resize(length);
101  D1.resize(length); D2.resize(length); D3.resize(length);D4.resize(length);
102  antiU1.resize(length); antiU2.resize(length); antiU3.resize(length);antiU4.resize(length);
103  antiD1.resize(length); antiD2.resize(length); antiD3.resize(length);antiD4.resize(length);
104  U1=0;U2=0;U3=0;U4=0;
105  multi1d<int> currsrc1=src1;
106  multi1d<int> currsrc2=src2;
107  for (int t=0; t<length; t++)
108  {
109  currsrc1[Nd-1]=t;
110  currsrc2[Nd-1]=t;
111  Q1[t] = peekSite(Qprop1,currsrc1); // HQ prop from src1 == "0"
112  Q2[t] = peekSite(Qprop2,currsrc2); // HQ prop from src2 == "R"
113  antiQ1[t] = adj(Q1[t]); // antiHQ prop from src1 == "0"
114  antiQ2[t] = adj(Q2[t]); // antiHQ prop from src2 == "R"
115 
116  U1[t] = peekSite(quark1,currsrc1); // propagator from 0 to 0
117  U2[t] = peekSite(quark2,currsrc2); // propagator from R to R
118  U3[t] = peekSite(quark2,currsrc1); // propagator from R to 0
119  U4[t] = peekSite(quark1,currsrc2); // propagator from 0 to R
120 
121  antiU1[t]= Gamma(15) * adj(U1[t]) * Gamma(15);
122  antiU2[t]= Gamma(15) * adj(U2[t]) * Gamma(15);
123  antiU3[t]= Gamma(15) * adj(U3[t]) * Gamma(15);
124  antiU4[t]= Gamma(15) * adj(U4[t]) * Gamma(15);
125 
126  }
127 
128  D1 = U1; D2 = U2; D3 = U3; D4 = U4; // D quarks are the same as U quarks
129  antiD1 = antiU1; antiD2 = antiU2; antiD3 = antiU3; antiD4 = antiU4;
130 
131  // MESON BLOCKS
132 
133  QDPIO::cout<<"Making B^+ (bbar u) meson blocks\n";
134 
135  HeavyMesonBlock H_U4_G5_R(Nt,U4,G5,antiQ2,NegEnergyProj);
136  HeavyMesonBlock H_U2_G5_R(Nt,U2,G5,antiQ2,NegEnergyProj);
137  HeavyMesonBlock H_U3_G5_0(Nt,U3,G5,antiQ1,NegEnergyProj);
138  HeavyMesonBlock H_U1_G5_0(Nt,U1,G5,antiQ1,NegEnergyProj);
139 
140  HeavyMesonBlock H_U4_Gup_R(Nt,U4,Gup,antiQ2,NegEnergyProj);
141  HeavyMesonBlock H_U2_Gup_R(Nt,U2,Gup,antiQ2,NegEnergyProj);
142  HeavyMesonBlock H_U3_Gup_0(Nt,U3,Gup,antiQ1,NegEnergyProj);
143  HeavyMesonBlock H_U1_Gup_0(Nt,U1,Gup,antiQ1,NegEnergyProj);
144 
145  HeavyMesonBlock H_U4_Gdown_R(Nt,U4,Gdown,antiQ2,NegEnergyProj);
146  HeavyMesonBlock H_U2_Gdown_R(Nt,U2,Gdown,antiQ2,NegEnergyProj);
147  HeavyMesonBlock H_U3_Gdown_0(Nt,U3,Gdown,antiQ1,NegEnergyProj);
148  HeavyMesonBlock H_U1_Gdown_0(Nt,U1,Gdown,antiQ1,NegEnergyProj);
149 
150  QDPIO::cout<<"Making B^0 (bbar d) meson blocks\n";
151  HeavyMesonBlock H_D4_G5_R = H_U4_G5_R;
152  HeavyMesonBlock H_D2_G5_R = H_U2_G5_R;
153  HeavyMesonBlock H_D3_G5_0 = H_U3_G5_0;
154  HeavyMesonBlock H_D1_G5_0 = H_U1_G5_0;
155 
156  HeavyMesonBlock H_D4_Gup_R = H_U4_Gup_R;
157  HeavyMesonBlock H_D2_Gup_R = H_U2_Gup_R;
158  HeavyMesonBlock H_D3_Gup_0 = H_U3_Gup_0;
159  HeavyMesonBlock H_D1_Gup_0 = H_U1_Gup_0;
160 
161  HeavyMesonBlock H_D4_Gdown_R = H_U4_Gdown_R;
162  HeavyMesonBlock H_D2_Gdown_R = H_U2_Gdown_R;
163  HeavyMesonBlock H_D3_Gdown_0 = H_U3_Gdown_0;
164  HeavyMesonBlock H_D1_Gdown_0 = H_U1_Gdown_0;
165 
166  QDPIO::cout<<"Making Bbar^0 (b dbar) meson blocks\n";
167  HeavyMesonBlock Hbar_D4_G5_R(Nt,antiD4,G5,Q2,PosEnergyProj);
168  HeavyMesonBlock Hbar_D2_G5_R(Nt,antiD2,G5,Q2,PosEnergyProj);
169  HeavyMesonBlock Hbar_D3_G5_0(Nt,antiD3,G5,Q1,PosEnergyProj);
170  HeavyMesonBlock Hbar_D1_G5_0(Nt,antiD1,G5,Q1,PosEnergyProj);
171 
172  HeavyMesonBlock Hbar_D4_Gup_R(Nt,antiD4,Gup,Q2,PosEnergyProj);
173  HeavyMesonBlock Hbar_D2_Gup_R(Nt,antiD2,Gup,Q2,PosEnergyProj);
174  HeavyMesonBlock Hbar_D3_Gup_0(Nt,antiD3,Gup,Q1,PosEnergyProj);
175  HeavyMesonBlock Hbar_D1_Gup_0(Nt,antiD1,Gup,Q1,PosEnergyProj);
176 
177  HeavyMesonBlock Hbar_D4_Gdown_R(Nt,antiD4,Gdown,Q2,PosEnergyProj);
178  HeavyMesonBlock Hbar_D2_Gdown_R(Nt,antiD2,Gdown,Q2,PosEnergyProj);
179  HeavyMesonBlock Hbar_D3_Gdown_0(Nt,antiD3,Gdown,Q1,PosEnergyProj);
180  HeavyMesonBlock Hbar_D1_Gdown_0(Nt,antiD1,Gdown,Q1,PosEnergyProj);
181 
182 
183  // B meson
184  QDPIO::cout<<" Contracting B meson \n";
185  multi1d<DComplex> bmes, bmesup, bmesdown, bmes0, bmes0up, bmes0down, bmesR, bmesRup, bmesRdown;
186  bmes.resize(length); bmesup.resize(length); bmesdown.resize(length);
187  bmes0.resize(length); bmes0up.resize(length); bmes0down.resize(length);
188  bmesR.resize(length); bmesRup.resize(length); bmesRdown.resize(length);
189  bmes0 = bcontract(H_U1_G5_0,G5);
190  bmesR = bcontract(H_U2_G5_R,G5);
191  bmes0up = bcontract(H_U1_Gup_0,Gup);
192  bmes0down = bcontract(H_U1_Gdown_0,Gdown);
193  bmesRup = bcontract(H_U2_Gup_R,Gup);
194  bmesRdown = bcontract(H_U2_Gdown_R,Gdown);
195  for (int t=0; t<length;t++){
196  bmes[t] = 0.5 * (bmes0[t] + bmesR[t]);
197  bmesup[t] = 0.5 * (bmes0up[t] + bmesRup[t]);
198  bmesdown[t] = 0.5 * (bmes0down[t] + bmesRdown[t]);
199  };
200 
201  // B B I=I_3=1: Bu Bu
202  QDPIO::cout<<" Contracting Bu--Bu \n";
203  multi1d<DComplex> BuBuJ0m0,BuBuJ1m1,BuBuJ1m0,BuBuJ1mneg1;
204  // J=0
205  BuBuJ0m0.resize(length);
206  BuBuJ0m0 = m1contract(H_U1_Gup_0, H_U2_Gdown_R,
207  H_U3_Gup_0, H_U4_Gdown_R,
208  Gup, Gdown);
209  BuBuJ0m0 += m1contract(H_U1_Gdown_0, H_U2_Gup_R,
210  H_U3_Gdown_0, H_U4_Gup_R,
211  Gdown, Gup);
212  BuBuJ0m0 -= m1contract(H_U1_Gup_0, H_U2_Gdown_R,
213  H_U3_Gup_0, H_U4_Gdown_R,
214  Gdown, Gup);
215  BuBuJ0m0 -= m1contract(H_U1_Gdown_0, H_U2_Gup_R,
216  H_U3_Gdown_0, H_U4_Gup_R,
217  Gup, Gdown);
218  BuBuJ0m0 *= 0.5;
219 
220  //J=1
221  BuBuJ1m1.resize(length);
222  BuBuJ1m0.resize(length);
223  BuBuJ1mneg1.resize(length);
224 
225  //J=1, m=1
226  BuBuJ1m1 = m1contract(H_U1_Gup_0, H_U2_Gup_R,
227  H_U3_Gup_0, H_U4_Gup_R,
228  Gup, Gup);
229 
230  //J=1, m=0
231  BuBuJ1m0 = m1contract(H_U1_Gup_0, H_U2_Gdown_R,
232  H_U3_Gup_0, H_U4_Gdown_R,
233  Gup, Gdown);
234  BuBuJ1m0 += m1contract(H_U1_Gdown_0, H_U2_Gup_R,
235  H_U3_Gdown_0, H_U4_Gup_R,
236  Gdown, Gup);
237  BuBuJ1m0 += m1contract(H_U1_Gup_0, H_U2_Gdown_R,
238  H_U3_Gup_0, H_U4_Gdown_R,
239  Gdown, Gup);
240  BuBuJ1m0 += m1contract(H_U1_Gdown_0, H_U2_Gup_R,
241  H_U3_Gdown_0, H_U4_Gup_R,
242  Gup, Gdown);
243  BuBuJ1m0 *= 0.5;
244 
245  //J=1, m=-1
246  BuBuJ1mneg1 = m1contract(H_U1_Gdown_0, H_U2_Gdown_R,
247  H_U3_Gdown_0, H_U4_Gdown_R,
248  Gdown, Gdown);
249 
250  // B B I=0,1 with I_3=0: Bu Bd
251  // since there are two different particles, we might symmetrise over
252  // Bu(r) Bd(0) --> Bu(r) Bd(0) and Bd(r) Bu(0) --> Bd(r) Bu(0)
253  // but since the u and d quark propagators are the same this is
254  // pointless in this case
255 
256  QDPIO::cout<<" Contracting Bu--Bd \n";
257  multi1d<DComplex> BuBdJ0m0,BuBdJ1m1,BuBdJ1m0,BuBdJ1mneg1;
258  BuBdJ0m0.resize(length);
259  BuBdJ0m0 = m2contract(H_U1_Gup_0, H_D2_Gdown_R,
260  Gup, Gdown);
261  BuBdJ0m0 -= m2contract(H_U1_Gdown_0, H_D2_Gup_R,
262  Gup, Gdown);
263  BuBdJ0m0 += m2contract(H_U1_Gdown_0, H_D2_Gup_R,
264  Gdown, Gup);
265  BuBdJ0m0 -= m2contract(H_U1_Gup_0, H_D2_Gdown_R,
266  Gdown, Gup);
267  BuBdJ0m0 *= 0.5;
268 
269 
270  BuBdJ1m1.resize(length);
271  BuBdJ1m0.resize(length);
272  BuBdJ1mneg1.resize(length);
273 
274  SpinMatrix G5Gup = Gamma(15)* Gup;
275 
276  BuBdJ1m1 = m2contract(H_U1_Gup_0, H_D2_Gup_R,
277  Gup, Gup);
278 
279  BuBdJ1m0 = m2contract(H_U1_Gup_0, H_D2_Gdown_R,
280  Gup, Gdown);
281  BuBdJ1m0 += m2contract(H_U1_Gdown_0, H_D2_Gup_R,
282  Gup, Gdown);
283  BuBdJ1m0 += m2contract(H_U1_Gdown_0, H_D2_Gup_R,
284  Gdown, Gup);
285  BuBdJ1m0 += m2contract(H_U1_Gup_0, H_D2_Gdown_R,
286  Gdown, Gup);
287  BuBdJ1m0 *= 0.5;
288 
289 
290  BuBdJ1mneg1 = m2contract(H_U1_Gdown_0, H_D2_Gdown_R,
291  Gdown, Gdown);
292 
293  // B Bbar I=0,1 with I_3=0: Bu Bdbar
294  // since there are two different particles, we might symmetrise over
295  // Bu(r) Bd(0) --> Bu(r) Bd(0) and Bd(r) Bu(0) --> Bd(r) Bu(0)
296  // but since the u and d quark propagators are the same this is
297  // pointless in this case
298 
299  QDPIO::cout<<" Contracting Bu--Bdbar \n";
300  multi1d<DComplex> BuBdbarJ0m0,BuBdbarJ1m1,BuBdbarJ1m0,BuBdbarJ1mneg1;
301  BuBdbarJ0m0.resize(length);
302  BuBdbarJ0m0 = m2contract(H_U1_Gup_0, Hbar_D2_Gdown_R,
303  Gup, Gdown);
304  BuBdbarJ0m0 -= m2contract(H_U1_Gdown_0, Hbar_D2_Gup_R,
305  Gup, Gdown);
306  BuBdbarJ0m0 += m2contract(H_U1_Gdown_0, Hbar_D2_Gup_R,
307  Gdown, Gup);
308  BuBdbarJ0m0 -= m2contract(H_U1_Gup_0, Hbar_D2_Gdown_R,
309  Gdown, Gup);
310  BuBdbarJ0m0 *= 0.5;
311 
312 
313  BuBdbarJ1m1.resize(length);
314  BuBdbarJ1m0.resize(length);
315  BuBdbarJ1mneg1.resize(length);
316 
317  BuBdbarJ1m1 = m2contract(H_U1_Gup_0, Hbar_D2_Gup_R,
318  Gup, Gup);
319 
320 
321  BuBdbarJ1m0 = m2contract(H_U1_Gup_0, Hbar_D2_Gdown_R,
322  Gup, Gdown);
323  BuBdbarJ1m0 += m2contract(H_U1_Gdown_0, Hbar_D2_Gup_R,
324  Gup, Gdown);
325  BuBdbarJ1m0 += m2contract(H_U1_Gdown_0, Hbar_D2_Gup_R,
326  Gdown, Gup);
327  BuBdbarJ1m0 += m2contract(H_U1_Gup_0, Hbar_D2_Gdown_R,
328  Gdown, Gup);
329  BuBdbarJ1m0 *= 0.5;
330 
331 
332  BuBdbarJ1mneg1 = m2contract(H_U1_Gdown_0, Hbar_D2_Gdown_R,
333  Gdown, Gdown);
334 
335 
336  // BARYON BLOCKS
337 
338  // Add epsilon tensors onto ends of HQ propagators
339  multiNd<DComplex> HQBR, HQB0;
340  multiNd<DComplex> antiHQBR, antiHQB0;
341  multi1d<int> HQBarray;
342  HQBarray.resize(4);
343  HQBarray= Nc; HQBarray[0]=length;
344  HQBR.resize(HQBarray);
345  HQB0.resize(HQBarray);
346  HQB0 = HBQfunc(Q1);
347  HQBR = HBQfunc(Q2);
348  antiHQBR.resize(HQBarray);
349  antiHQB0.resize(HQBarray);
350  antiHQB0 = HBQfunc(antiQ1);
351  antiHQBR = HBQfunc(antiQ2);
352 
353  QDPIO::cout<<"Making lambdaB blocks\n";
354  // Lambda_b
355  // Note we will use spin to differentiate the Lambda_b from the \Sigma_B^0
356  // and simply take Bud to be the flavour structure
357  QllBlock B_U4_D2_R_Cg5(Nt,U4,D2,Cg5,HQBR);
358  QllBlock B_U3_D1_0_Cg5(Nt,U3,D1,Cg5,HQB0);
359  QllBlock B_U2_D4_R_Cg5(Nt,U2,D4,Cg5,HQBR);
360  QllBlock B_U1_D3_0_Cg5(Nt,U1,D3,Cg5,HQB0);
361  QllBlock B_U4_D4_R_Cg5(Nt,U4,D4,Cg5,HQBR);
362  QllBlock B_U3_D3_0_Cg5(Nt,U3,D3,Cg5,HQB0);
363  QllBlock B_U2_D2_R_Cg5(Nt,U2,D2,Cg5,HQBR);
364  QllBlock B_U1_D1_0_Cg5(Nt,U1,D1,Cg5,HQB0);
365 
366  // Lambda_b 2 pt
367  QDPIO::cout<<" Contracting Lambda_b \n";
368  multi1d<DComplex> lambdab,lambdaba,lambdabb;
369  lambdab.resize(length); lambdaba.resize(length); lambdabb.resize(length);
370  lambdaba = lambdabcontract(B_U1_D1_0_Cg5, Cg5);
371  lambdabb = lambdabcontract(B_U2_D2_R_Cg5, Cg5);
372  for (int t=0; t<length;t++){
373  lambdab[t]= 0.5 * ( lambdaba[t]+lambdabb[t]);
374  };
375 
376  // Lambda_b Lambda_b contractions
377  QDPIO::cout<<" Contracting Lambda_b--Lambda_b \n";
378  multi1d<DComplex> lambdablambdab;
379  lambdablambdab.resize(length);
380  lambdablambdab = c1contract(B_U1_D1_0_Cg5, B_U1_D3_0_Cg5,
381  B_U2_D2_R_Cg5, B_U2_D4_R_Cg5,
382  B_U3_D1_0_Cg5, B_U3_D3_0_Cg5,
383  B_U4_D2_R_Cg5, B_U4_D4_R_Cg5,
384  Cg5, Cg5);
385 
386  QDPIO::cout<<"Making sigmaB^+ blocks\n";
387  QllBlock B_U4_U2_R_Cg3(Nt,U4,U2,Cg3,HQBR);
388  QllBlock B_U3_U1_0_Cg3(Nt,U3,U1,Cg3,HQB0);
389  QllBlock B_U2_U4_R_Cg3(Nt,U2,U4,Cg3,HQBR);
390  QllBlock B_U1_U3_0_Cg3(Nt,U1,U3,Cg3,HQB0);
391  QllBlock B_U4_U4_R_Cg3(Nt,U4,U4,Cg3,HQBR);
392  QllBlock B_U3_U3_0_Cg3(Nt,U3,U3,Cg3,HQB0);
393  QllBlock B_U2_U2_R_Cg3(Nt,U2,U2,Cg3,HQBR);
394  QllBlock B_U1_U1_0_Cg3(Nt,U1,U1,Cg3,HQB0);
395 
396  QllBlock B_U4_U2_R_Cgplus(Nt,U4,U2,Cgplus,HQBR);
397  QllBlock B_U3_U1_0_Cgplus(Nt,U3,U1,Cgplus,HQB0);
398  QllBlock B_U2_U4_R_Cgplus(Nt,U2,U4,Cgplus,HQBR);
399  QllBlock B_U1_U3_0_Cgplus(Nt,U1,U3,Cgplus,HQB0);
400  QllBlock B_U4_U4_R_Cgplus(Nt,U4,U4,Cgplus,HQBR);
401  QllBlock B_U3_U3_0_Cgplus(Nt,U3,U3,Cgplus,HQB0);
402  QllBlock B_U2_U2_R_Cgplus(Nt,U2,U2,Cgplus,HQBR);
403  QllBlock B_U1_U1_0_Cgplus(Nt,U1,U1,Cgplus,HQB0);
404 
405  QllBlock B_U4_U2_R_Cgminus(Nt,U4,U2,Cgminus,HQBR);
406  QllBlock B_U3_U1_0_Cgminus(Nt,U3,U1,Cgminus,HQB0);
407  QllBlock B_U2_U4_R_Cgminus(Nt,U2,U4,Cgminus,HQBR);
408  QllBlock B_U1_U3_0_Cgminus(Nt,U1,U3,Cgminus,HQB0);
409  QllBlock B_U4_U4_R_Cgminus(Nt,U4,U4,Cgminus,HQBR);
410  QllBlock B_U3_U3_0_Cgminus(Nt,U3,U3,Cgminus,HQB0);
411  QllBlock B_U2_U2_R_Cgminus(Nt,U2,U2,Cgminus,HQBR);
412  QllBlock B_U1_U1_0_Cgminus(Nt,U1,U1,Cgminus,HQB0);
413 
414  QDPIO::cout<<"Making sigmaB^0 blocks\n";
415  QllBlock B_U4_D2_R_Cg3(Nt,U4,D2,Cg3,HQBR);
416  QllBlock B_U3_D1_0_Cg3(Nt,U3,D1,Cg3,HQB0);
417  QllBlock B_U2_D4_R_Cg3(Nt,U2,D4,Cg3,HQBR);
418  QllBlock B_U1_D3_0_Cg3(Nt,U1,D3,Cg3,HQB0);
419  QllBlock B_U4_D4_R_Cg3(Nt,U4,D4,Cg3,HQBR);
420  QllBlock B_U3_D3_0_Cg3(Nt,U3,D3,Cg3,HQB0);
421  QllBlock B_U2_D2_R_Cg3(Nt,U2,D2,Cg3,HQBR);
422  QllBlock B_U1_D1_0_Cg3(Nt,U1,D1,Cg3,HQB0);
423 
424  QllBlock B_U4_D2_R_Cgplus(Nt,U4,D2,Cgplus,HQBR);
425  QllBlock B_U3_D1_0_Cgplus(Nt,U3,D1,Cgplus,HQB0);
426  QllBlock B_U2_D4_R_Cgplus(Nt,U2,D4,Cgplus,HQBR);
427  QllBlock B_U1_D3_0_Cgplus(Nt,U1,D3,Cgplus,HQB0);
428  QllBlock B_U4_D4_R_Cgplus(Nt,U4,D4,Cgplus,HQBR);
429  QllBlock B_U3_D3_0_Cgplus(Nt,U3,D3,Cgplus,HQB0);
430  QllBlock B_U2_D2_R_Cgplus(Nt,U2,D2,Cgplus,HQBR);
431  QllBlock B_U1_D1_0_Cgplus(Nt,U1,D1,Cgplus,HQB0);
432 
433  QllBlock B_U4_D2_R_Cgminus(Nt,U4,D2,Cgminus,HQBR);
434  QllBlock B_U3_D1_0_Cgminus(Nt,U3,D1,Cgminus,HQB0);
435  QllBlock B_U2_D4_R_Cgminus(Nt,U2,D4,Cgminus,HQBR);
436  QllBlock B_U1_D3_0_Cgminus(Nt,U1,D3,Cgminus,HQB0);
437  QllBlock B_U4_D4_R_Cgminus(Nt,U4,D4,Cgminus,HQBR);
438  QllBlock B_U3_D3_0_Cgminus(Nt,U3,D3,Cgminus,HQB0);
439  QllBlock B_U2_D2_R_Cgminus(Nt,U2,D2,Cgminus,HQBR);
440  QllBlock B_U1_D1_0_Cgminus(Nt,U1,D1,Cgminus,HQB0);
441 
442 
443  QDPIO::cout<<"Making sigmaB^- blocks\n";
444  QllBlock B_D4_D2_R_Cg3(Nt,D4,D2,Cg3,HQBR);
445  QllBlock B_D3_D1_0_Cg3(Nt,D3,D1,Cg3,HQB0);
446  QllBlock B_D2_D4_R_Cg3(Nt,D2,D4,Cg3,HQBR);
447  QllBlock B_D1_D3_0_Cg3(Nt,D1,D3,Cg3,HQB0);
448  QllBlock B_D4_D4_R_Cg3(Nt,D4,D4,Cg3,HQBR);
449  QllBlock B_D3_D3_0_Cg3(Nt,D3,D3,Cg3,HQB0);
450  QllBlock B_D2_D2_R_Cg3(Nt,D2,D2,Cg3,HQBR);
451  QllBlock B_D1_D1_0_Cg3(Nt,D1,D1,Cg3,HQB0);
452 
453  QllBlock B_D4_D2_R_Cgplus(Nt,D4,D2,Cgplus,HQBR);
454  QllBlock B_D3_D1_0_Cgplus(Nt,D3,D1,Cgplus,HQB0);
455  QllBlock B_D2_D4_R_Cgplus(Nt,D2,D4,Cgplus,HQBR);
456  QllBlock B_D1_D3_0_Cgplus(Nt,D1,D3,Cgplus,HQB0);
457  QllBlock B_D4_D4_R_Cgplus(Nt,D4,D4,Cgplus,HQBR);
458  QllBlock B_D3_D3_0_Cgplus(Nt,D3,D3,Cgplus,HQB0);
459  QllBlock B_D2_D2_R_Cgplus(Nt,D2,D2,Cgplus,HQBR);
460  QllBlock B_D1_D1_0_Cgplus(Nt,D1,D1,Cgplus,HQB0);
461 
462  QllBlock B_D4_D2_R_Cgminus(Nt,D4,D2,Cgminus,HQBR);
463  QllBlock B_D3_D1_0_Cgminus(Nt,D3,D1,Cgminus,HQB0);
464  QllBlock B_D2_D4_R_Cgminus(Nt,D2,D4,Cgminus,HQBR);
465  QllBlock B_D1_D3_0_Cgminus(Nt,D1,D3,Cgminus,HQB0);
466  QllBlock B_D4_D4_R_Cgminus(Nt,D4,D4,Cgminus,HQBR);
467  QllBlock B_D3_D3_0_Cgminus(Nt,D3,D3,Cgminus,HQB0);
468  QllBlock B_D2_D2_R_Cgminus(Nt,D2,D2,Cgminus,HQBR);
469  QllBlock B_D1_D1_0_Cgminus(Nt,D1,D1,Cgminus,HQB0);
470 
471 
472  QDPIO::cout<<"Making anti-sigmaB^- blocks\n";
473  QllBlock Bbar_D4_D2_R_Cg3(Nt,antiD4,antiD2,Cg3,antiHQBR);
474  QllBlock Bbar_D3_D1_0_Cg3(Nt,antiD3,antiD1,Cg3,antiHQB0);
475  QllBlock Bbar_D2_D4_R_Cg3(Nt,antiD2,antiD4,Cg3,antiHQBR);
476  QllBlock Bbar_D1_D3_0_Cg3(Nt,antiD1,antiD3,Cg3,antiHQB0);
477  QllBlock Bbar_D4_D4_R_Cg3(Nt,antiD4,antiD4,Cg3,antiHQBR);
478  QllBlock Bbar_D3_D3_0_Cg3(Nt,antiD3,antiD3,Cg3,antiHQB0);
479  QllBlock Bbar_D2_D2_R_Cg3(Nt,antiD2,antiD2,Cg3,antiHQBR);
480  QllBlock Bbar_D1_D1_0_Cg3(Nt,antiD1,antiD1,Cg3,antiHQB0);
481 
482  QllBlock Bbar_D4_D2_R_Cgplus(Nt,antiD4,antiD2,Cgplus,antiHQBR);
483  QllBlock Bbar_D3_D1_0_Cgplus(Nt,antiD3,antiD1,Cgplus,antiHQB0);
484  QllBlock Bbar_D2_D4_R_Cgplus(Nt,antiD2,antiD4,Cgplus,antiHQBR);
485  QllBlock Bbar_D1_D3_0_Cgplus(Nt,antiD1,antiD3,Cgplus,antiHQB0);
486  QllBlock Bbar_D4_D4_R_Cgplus(Nt,antiD4,antiD4,Cgplus,antiHQBR);
487  QllBlock Bbar_D3_D3_0_Cgplus(Nt,antiD3,antiD3,Cgplus,antiHQB0);
488  QllBlock Bbar_D2_D2_R_Cgplus(Nt,antiD2,antiD2,Cgplus,antiHQBR);
489  QllBlock Bbar_D1_D1_0_Cgplus(Nt,antiD1,antiD1,Cgplus,antiHQB0);
490 
491  QllBlock Bbar_D4_D2_R_Cgminus(Nt,antiD4,antiD2,Cgminus,antiHQBR);
492  QllBlock Bbar_D3_D1_0_Cgminus(Nt,antiD3,antiD1,Cgminus,antiHQB0);
493  QllBlock Bbar_D2_D4_R_Cgminus(Nt,antiD2,antiD4,Cgminus,antiHQBR);
494  QllBlock Bbar_D1_D3_0_Cgminus(Nt,antiD1,antiD3,Cgminus,antiHQB0);
495  QllBlock Bbar_D4_D4_R_Cgminus(Nt,antiD4,antiD4,Cgminus,antiHQBR);
496  QllBlock Bbar_D3_D3_0_Cgminus(Nt,antiD3,antiD3,Cgminus,antiHQB0);
497  QllBlock Bbar_D2_D2_R_Cgminus(Nt,antiD2,antiD2,Cgminus,antiHQBR);
498  QllBlock Bbar_D1_D1_0_Cgminus(Nt,antiD1,antiD1,Cgminus,antiHQB0);
499 
500  // Sigma_b^+ and (=) \Sigma_b^- baryon 2pt
501  QDPIO::cout<<" Contracting Sigma_b^+ \n";
502  multi1d<DComplex> sigmabplusJ1m1, sigmabplusJ1m1a, sigmabplusJ1m1b;
503  multi1d<DComplex> sigmabplusJ1m0,sigmabplusJ1m0a,sigmabplusJ1m0b;
504  multi1d<DComplex> sigmabplusJ1mneg1,sigmabplusJ1mneg1a,sigmabplusJ1mneg1b;
505 
506  sigmabplusJ1m1.resize(length); sigmabplusJ1m1a.resize(length); sigmabplusJ1m1b.resize(length);
507  sigmabplusJ1m0.resize(length);sigmabplusJ1m0a.resize(length);sigmabplusJ1m0b.resize(length);
508  sigmabplusJ1mneg1.resize(length);sigmabplusJ1mneg1a.resize(length);sigmabplusJ1mneg1b.resize(length);
509 
510  sigmabplusJ1m1a = sigmabpluscontract(B_U1_U1_0_Cgplus, Cgplus);
511  sigmabplusJ1m1b = sigmabpluscontract(B_U2_U2_R_Cgplus, Cgplus);
512  sigmabplusJ1m0a = sigmabpluscontract(B_U1_U1_0_Cg3, Cg3);
513  sigmabplusJ1m0b = sigmabpluscontract(B_U2_U2_R_Cg3, Cg3);
514  sigmabplusJ1mneg1a = sigmabpluscontract(B_U1_U1_0_Cgminus, Cgminus);
515  sigmabplusJ1mneg1b = sigmabpluscontract(B_U2_U2_R_Cgminus, Cgminus);
516  for (int t=0; t<length;t++){
517  sigmabplusJ1m1[t] = 0.5 * (sigmabplusJ1m1a[t]+sigmabplusJ1m1b[t]);
518  sigmabplusJ1m0[t] = 0.5 * (sigmabplusJ1m0a[t]+sigmabplusJ1m0b[t]);
519  sigmabplusJ1mneg1[t] = 0.5 * (sigmabplusJ1mneg1a[t]+sigmabplusJ1mneg1b[t]);
520  };
521 
522 
523  // Sigma_b^0 2pt : average over two positions
524  QDPIO::cout<<" Contracting Sigma_b^0 \n";
525  multi1d<DComplex> sigmabzeroJ1m1, sigmabzeroJ1m0, sigmabzeroJ1mneg1;
526 
527  sigmabzeroJ1m1 = lambdabcontract(B_U1_D1_0_Cgplus, Cgplus);
528  sigmabzeroJ1m1 += lambdabcontract(B_U2_D2_R_Cgplus, Cgplus);
529  sigmabzeroJ1m1 *= 0.5;
530 
531  sigmabzeroJ1m0 = lambdabcontract(B_U1_D1_0_Cg3, Cg3);
532  sigmabzeroJ1m0 += lambdabcontract(B_U2_D2_R_Cg3, Cg3);
533  sigmabzeroJ1m0 *= 0.5;
534 
535  sigmabzeroJ1mneg1 = lambdabcontract(B_U1_D1_0_Cgminus, Cgminus);
536  sigmabzeroJ1mneg1 += lambdabcontract(B_U2_D2_R_Cgminus, Cgminus);
537  sigmabzeroJ1mneg1 *= 0.5;
538 
539 
540  // Sigma_b^+ \Sigma_b^+ contractions: I=I_3=2, J=0,1,2
541  QDPIO::cout<<" Contracting Sigma_b^+--Sigma_b^+ \n";
542  multi2d<DComplex> sigmabplussigmabplus;
543  sigmabplussigmabplus.resize(9,length); // Store all the spin states
544  sigmabplussigmabplus = c5J2corr(// Spin up blocks
545  B_U1_U1_0_Cgplus, B_U1_U3_0_Cgplus,
546  B_U2_U2_R_Cgplus, B_U2_U4_R_Cgplus,
547  B_U3_U1_0_Cgplus, B_U3_U3_0_Cgplus,
548  B_U4_U2_R_Cgplus, B_U4_U4_R_Cgplus,
549  // Spin zero blocks
550  B_U1_U1_0_Cg3, B_U1_U3_0_Cg3,
551  B_U2_U2_R_Cg3, B_U2_U4_R_Cg3,
552  B_U3_U1_0_Cg3, B_U3_U3_0_Cg3,
553  B_U4_U2_R_Cg3, B_U4_U4_R_Cg3,
554  // Spin down blocks
555  B_U1_U1_0_Cgminus, B_U1_U3_0_Cgminus,
556  B_U2_U2_R_Cgminus, B_U2_U4_R_Cgminus,
557  B_U3_U1_0_Cgminus, B_U3_U3_0_Cgminus,
558  B_U4_U2_R_Cgminus, B_U4_U4_R_Cgminus);
559 
560  // Sigma_b^+ \Sigma_b^0 contractions: I_3=1 (I=1,2), J=0,1,2
561  QDPIO::cout<<" Contracting Sigma_b^+--Sigma_b^0 \n";
562  multi2d<DComplex> sigmabplussigmabzero;
563  sigmabplussigmabzero.resize(9,length); // Store all the spin states
564  sigmabplussigmabzero = c4J2corr(// Spin up blocks
565  B_U1_D1_0_Cgplus, B_U2_U2_R_Cgplus,
566  B_U2_U4_R_Cgplus, B_U3_D1_0_Cgplus,
567  B_U4_U2_R_Cgplus,
568  // Spin zero blocks
569  B_U1_D1_0_Cg3, B_U2_U2_R_Cg3,
570  B_U2_U4_R_Cg3, B_U3_D1_0_Cg3,
571  B_U4_U2_R_Cg3,
572  // Spin down blocks
573  B_U1_D1_0_Cgminus, B_U2_U2_R_Cgminus,
574  B_U2_U4_R_Cgminus, B_U3_D1_0_Cgminus,
575  B_U4_U2_R_Cgminus);
576 
577 
578  // Sigma_b^+ \Sigma_b^- contractions: I_3=0 (I=0,1,2), J=0,1,2
579  QDPIO::cout<<" Contracting Sigma_b^+--Sigma_b^- \n";
580  multi2d<DComplex> sigmabplussigmabminus;
581  sigmabplussigmabminus.resize(9,length); // Store all the spin states
582  sigmabplussigmabminus = c6J2corr(B_D1_D1_0_Cgplus, B_U2_U2_R_Cgplus,
583  B_D1_D1_0_Cg3, B_U2_U2_R_Cg3,
584  B_D1_D1_0_Cgminus, B_U2_U2_R_Cgminus);
585 
586  // Sigma_b^+ \Lambda_b contractions: I=I_3=1, J=1
587  QDPIO::cout<<" Contracting Sigma_b^+--Lambda_b \n";
588  multi1d<DComplex> sigmabpluslambdabJ1m1, sigmabpluslambdabJ1m0, sigmabpluslambdabJ1mneg1;
589  sigmabpluslambdabJ1m1.resize(length);
590  sigmabpluslambdabJ1m0.resize(length);
591  sigmabpluslambdabJ1mneg1.resize(length);
592 
593  sigmabpluslambdabJ1m1 = c4contract(B_U1_D1_0_Cg5, B_U2_U2_R_Cgplus,
594  B_U2_U4_R_Cgplus, B_U3_D1_0_Cg5,
595  B_U4_U2_R_Cgplus,
596  Cgplus, Cg5);
597  sigmabpluslambdabJ1m0 = c4contract(B_U1_D1_0_Cg5, B_U2_U2_R_Cg3,
598  B_U2_U4_R_Cg3, B_U3_D1_0_Cg5,
599  B_U4_U2_R_Cg3,
600  Cg3, Cg5);
601  sigmabpluslambdabJ1mneg1 = c4contract(B_U1_D1_0_Cg5, B_U2_U2_R_Cgminus,
602  B_U2_U4_R_Cgminus, B_U3_D1_0_Cg5,
603  B_U4_U2_R_Cgminus,
604  Cgminus, Cg5);
605 
606  // Sigma_b^+ \Lambda_b transition contractions: I=I_3=1, J=1
607  // Since there are two different particles, we symmetrise over
608  // A(r) B(0) --> A(r) B(0) and B(r) A(0) --> B(r) A(0)
609  QDPIO::cout<<" Contracting Sigma_b^+--Lambda_b transition \n";
610  multi1d<DComplex> sigmabpluslambdabtransitionJ1m1, sigmabpluslambdabtransitionJ1m0, sigmabpluslambdabtransitionJ1mneg1;
611  sigmabpluslambdabtransitionJ1m1.resize(length);
612  sigmabpluslambdabtransitionJ1m0.resize(length);
613  sigmabpluslambdabtransitionJ1mneg1.resize(length);
614 
615  sigmabpluslambdabtransitionJ1m1 = c7contract(B_U2_D4_R_Cg5, B_U1_U3_0_Cgplus, B_U4_D4_R_Cg5,
616  B_U3_U1_0_Cgplus, B_U3_U3_0_Cgplus,
617  Cg5, Cgplus);
618  sigmabpluslambdabtransitionJ1m0 = c7contract(B_U2_D4_R_Cg5, B_U1_U3_0_Cg3, B_U4_D4_R_Cg5,
619  B_U3_U1_0_Cg3, B_U3_U3_0_Cg3,
620  Cg5, Cg3);
621  sigmabpluslambdabtransitionJ1mneg1 = c7contract(B_U2_D4_R_Cg5, B_U1_U3_0_Cgminus, B_U4_D4_R_Cg5,
622  B_U3_U1_0_Cgminus, B_U3_U3_0_Cgminus,
623  Cg5, Cgminus);
624 // sigmabpluslambdabtransitionJ1m1 = c7contract(B_U1_D3_0_Cg5, B_U2_U4_R_Cgplus, B_U3_D3_0_Cg5,
625 // B_U4_U2_R_Cgplus, B_U4_U4_R_Cgplus,
626 // Cg5, Cgplus);
627 // sigmabpluslambdabtransitionJ1m0 = c7contract(B_U1_D3_0_Cg5, B_U2_U4_R_Cg3, B_U3_D3_0_Cg5,
628 // B_U4_U2_R_Cg3, B_U4_U4_R_Cg3,
629 // Cg5, Cg3);
630 // sigmabpluslambdabtransitionJ1mneg1 = c7contract(B_U1_D3_0_Cg5, B_U2_U4_R_Cgminus, B_U3_D3_0_Cg5,
631 // B_U4_U2_R_Cgminus, B_U4_U4_R_Cgminus,
632 // Cg5, Cgminus);
633 
634  QDPIO::cout<<" Contracting Sigma_b^+--antiSigma_b^- \n";
635  // \overline{Sigma}_b^+ \Sigma_b^- contractions: I=I_3=1, J=0,1,2
636  // Since there are two different particles, we symmetrise over
637  // A(r) B(0) --> A(r) B(0) and B(r) A(0) --> B(r) A(0)
638  multi2d<DComplex> sigmabplussigmabarbminusa,sigmabplussigmabarbminusb,sigmabplussigmabarbminus;
639  sigmabplussigmabarbminusa.resize(9,length); // Store all the spin states
640  sigmabplussigmabarbminusb.resize(9,length); // Store all the spin states
641  sigmabplussigmabarbminus.resize(9,length); // Store all the spin states
642  sigmabplussigmabarbminusa = c6J2corr(Bbar_D1_D1_0_Cgplus, B_U2_U2_R_Cgplus,
643  Bbar_D1_D1_0_Cg3, B_U2_U2_R_Cg3,
644  Bbar_D1_D1_0_Cgminus, B_U2_U2_R_Cgminus);
645  sigmabplussigmabarbminusb = c6J2corr(B_U1_U1_0_Cgplus, Bbar_D2_D2_R_Cgplus,
646  B_U1_U1_0_Cg3, Bbar_D2_D2_R_Cg3,
647  B_U1_U1_0_Cgminus, Bbar_D2_D2_R_Cgminus);
648  for (int i=0;i<9;i++){ for (int t=0; t<length;t++){
649  sigmabplussigmabarbminus[i][t] = 0.5 * (sigmabplussigmabarbminusa[i][t] + sigmabplussigmabarbminusb[i][t]);
650  }}
651 
652  // MESON-BARYON CONTRACTIONS
653 
654  QDPIO::cout<<" Contracting Bu--Lambda_b \n";
655  // Bu \Lambda_b contractions: I=I_3=1/2, J=1/2, m=1/2
656  // Since there are two different particles, we symmetrise over
657  // A(r) B(0) --> A(r) B(0) and B(r) A(0) --> B(r) A(0)
658  multi1d<DComplex> bulambdabJhalfmhalf;
659  bulambdabJhalfmhalf.resize(length);
660 
661  bulambdabJhalfmhalf = d1contract(B_U1_D1_0_Cg5, B_U3_D1_0_Cg5,
662  H_U2_Gup_R, H_U4_Gup_R,
663  Gup, Cg5);
664  bulambdabJhalfmhalf += d1contract(B_U2_D2_R_Cg5, B_U4_D2_R_Cg5,
665  H_U1_Gup_0, H_U3_Gup_0,
666  Gup, Cg5);
667  bulambdabJhalfmhalf *= 0.5;
668 
669  multi1d<DComplex> bulambdabJhalfmneghalf;
670  bulambdabJhalfmneghalf.resize(length);
671  bulambdabJhalfmneghalf = d1contract(B_U1_D1_0_Cg5, B_U3_D1_0_Cg5,
672  H_U2_Gdown_R, H_U4_Gdown_R,
673  Gdown, Cg5);
674  bulambdabJhalfmneghalf += d1contract(B_U2_D2_R_Cg5, B_U4_D2_R_Cg5,
675  H_U1_Gdown_0, H_U3_Gdown_0,
676  Gdown, Cg5);
677  bulambdabJhalfmneghalf *= 0.5;
678 
679  QDPIO::cout<<" Contracting Bu--Sigma_b^+ \n";
680  // Bu \Sigma_b^+ contractions: I=I_3=3/2, J=1/2, 3/2
681  // Since there are two different particles, we symmetrise over
682  // A(r) B(0) --> A(r) B(0) and B(r) A(0) --> B(r) A(0)
683  multi2d<DComplex> busigmabplus,busigmabplusa,busigmabplusb;
684  busigmabplus.resize(6,length);
685  busigmabplusa.resize(6,length);
686  busigmabplusb.resize(6,length);
687  busigmabplusa = d2J32corr(// Spin up blocks
688  B_U1_U1_0_Cgplus, B_U3_U1_0_Cgplus, B_U1_U3_0_Cgplus,
689  H_U2_Gup_R, H_U4_Gup_R,
690  // Spin zero blocks
691  B_U1_U1_0_Cg3, B_U3_U1_0_Cg3, B_U1_U3_0_Cg3,
692  // Spin down blocks
693  B_U1_U1_0_Cgminus, B_U3_U1_0_Cgminus, B_U1_U3_0_Cgminus,
694  H_U2_Gdown_R, H_U4_Gdown_R);
695  busigmabplusb = d2J32corr(// Spin up blocks
696  B_U2_U2_R_Cgplus, B_U4_U2_R_Cgplus, B_U2_U4_R_Cgplus,
697  H_U1_Gup_0, H_U3_Gup_0,
698  // Spin zero blocks
699  B_U2_U2_R_Cg3, B_U4_U2_R_Cg3, B_U2_U4_R_Cg3,
700  // Spin down blocks
701  B_U2_U2_R_Cgminus, B_U4_U2_R_Cgminus, B_U2_U4_R_Cgminus,
702  H_U1_Gdown_0, H_U3_Gdown_0);
703  for (int i=0;i<6;i++){ for (int t=0; t<length;t++){
704  busigmabplus[i][t] = 0.5 *(busigmabplusa[i][t] + busigmabplusb[i][t]);
705  }}
706 
707  QDPIO::cout<<" Contracting Bd--Sigma_b^+ \n";
708  // Bd \Sigma_b^+ contractions: I_3=1/2 (I=1/2, 3/2), J=1/2, 3/2
709  // Since there are two different particles, we symmetrise over
710  // A(r) B(0) --> A(r) B(0) and B(r) A(0) --> B(r) A(0)
711  multi2d<DComplex> bdsigmabplus,bdsigmabplusa,bdsigmabplusb;
712  bdsigmabplus.resize(6,length);
713  bdsigmabplusa.resize(6,length);
714  bdsigmabplusb.resize(6,length);
715  bdsigmabplusa = d3J32corr(// Spin up blocks
716  B_U1_U1_0_Cgplus,
717  H_D2_Gup_R,
718  // Spin zero blocks
719  B_U1_U1_0_Cg3,
720  // Spin down blocks
721  B_U1_U1_0_Cgminus,
722  H_D2_Gdown_R);
723  bdsigmabplusb = d3J32corr(// Spin up blocks
724  B_U2_U2_R_Cgplus,
725  H_D1_Gup_0,
726  // Spin zero blocks
727  B_U2_U2_R_Cg3,
728  // Spin down blocks
729  B_U2_U2_R_Cgminus,
730  H_D1_Gdown_0);
731  for (int i=0;i<6;i++){ for (int t=0; t<length;t++){
732  bdsigmabplus[i][t] = 0.5 *(bdsigmabplusa[i][t] + bdsigmabplusb[i][t]);
733  }}
734 
735  QDPIO::cout<<" Contracting Bbard--Sigma_b^+ \n";
736  // Bbard \Sigma_b^+ (Q dbar Quu) contractions: I_3=1/2 (I=1/2, 3/2), J=1/2, 3/2
737  // Since there are two different particles, we symmetrise over
738  // A(r) B(0) --> A(r) B(0) and B(r) A(0) --> B(r) A(0)
739  multi2d<DComplex> bbardsigmabplus,bbardsigmabplusa,bbardsigmabplusb;
740  bbardsigmabplus.resize(6,length);
741  bbardsigmabplusa.resize(6,length);
742  bbardsigmabplusb.resize(6,length);
743  bbardsigmabplusa = d3J32corr(// Spin up blocks
744  B_U1_U1_0_Cgplus,
745  Hbar_D2_Gup_R,
746  // Spin zero blocks
747  B_U1_U1_0_Cg3,
748  // Spin down blocks
749  B_U1_U1_0_Cgminus,
750  Hbar_D2_Gdown_R);
751  bbardsigmabplusb = d3J32corr(// Spin up blocks
752  B_U2_U2_R_Cgplus,
753  Hbar_D1_Gup_0,
754  // Spin zero blocks
755  B_U2_U2_R_Cg3,
756  // Spin down blocks
757  B_U2_U2_R_Cgminus,
758  Hbar_D1_Gdown_0);
759  for (int i=0;i<6;i++){ for (int t=0; t<length;t++){
760  bbardsigmabplus[i][t] = 0.5 *(bbardsigmabplusa[i][t] + bbardsigmabplusb[i][t]);
761  }}
762 
763 
764 
765  push(xml, xml_group);
766  // Hadrons
767  write(xml, "Bmes",bmes);
768  write(xml, "Bmes0",bmes0);
769  write(xml, "BmesR", bmesR);
770  write(xml, "BmesJ0.5m0.5",bmesup);
771  write(xml, "Bmes0J0.5m0.5",bmes0up);
772  write(xml, "BmesRJ0.5m0.5", bmesRup);
773  write(xml, "BmesJ0.5mneg0.5",bmesdown);
774  write(xml, "Bmes0J0.5mneg0.5",bmes0down);
775  write(xml, "BmesRJ0.5mneg0.5", bmesRdown);
776  write(xml,"lambdab", lambdab);
777  write(xml,"lambdab0", lambdaba);
778  write(xml,"lambdabR", lambdabb);
779  write(xml,"sigmabplusJ1m1", sigmabplusJ1m1);
780  write(xml,"sigmabplusJ1m0", sigmabplusJ1m0);
781  write(xml,"sigmabplusJ1mneg1", sigmabplusJ1mneg1);
782  write(xml,"sigmabplusJ1m1_0", sigmabplusJ1m1a);
783  write(xml,"sigmabplusJ1m0_0", sigmabplusJ1m0a);
784  write(xml,"sigmabplusJ1mneg1_0", sigmabplusJ1mneg1a);
785  write(xml,"sigmabplusJ1m1_R", sigmabplusJ1m1b);
786  write(xml,"sigmabplusJ1m0_R", sigmabplusJ1m0b);
787  write(xml,"sigmabplusJ1mneg1_R", sigmabplusJ1mneg1b);
788  write(xml,"sigmabzeroJ1m1", sigmabzeroJ1m1);
789  write(xml,"sigmabzeroJ1m0", sigmabzeroJ1m0);
790  write(xml,"sigmabzeroJ1mneg1", sigmabzeroJ1mneg1);
791  // Meson-meson
792  write(xml, "Bu_Bu_J0m0", BuBuJ0m0);
793  write(xml, "Bu_Bu_J1m1", BuBuJ1m1);
794  write(xml, "Bu_Bu_J1m0",BuBuJ1m0);
795  write(xml, "Bu_Bu_J1mneg1", BuBuJ1mneg1);
796  write(xml, "Bu_Bd_J0m0", BuBdJ0m0);
797  write(xml, "Bu_Bd_J1m1", BuBdJ1m1);
798  write(xml, "Bu_Bd_J1m0", BuBdJ1m0);
799  write(xml, "Bu_Bd_J1mneg1", BuBdJ1mneg1);
800  write(xml, "Bu_Bdbar_J0m0", BuBdbarJ0m0);
801  write(xml, "Bu_Bdbar_J1m1", BuBdbarJ1m1);
802  write(xml, "Bu_Bdbar_J1m0", BuBdbarJ1m0);
803  write(xml, "Bu_Bdbar_J1mneg1", BuBdbarJ1mneg1);
804  // Baryon-baryon
805  write(xml, "lambdaB_lambdaB_J0m0", lambdablambdab);
806  write(xml, "sigmabplus_sigmabplus_J2m2", sigmabplussigmabplus[0]);
807  write(xml, "sigmabplus_sigmabplus_J2m1", sigmabplussigmabplus[1]);
808  write(xml, "sigmabplus_sigmabplus_J2m0", sigmabplussigmabplus[2]);
809  write(xml, "sigmabplus_sigmabplus_J2mneg1", sigmabplussigmabplus[3]);
810  write(xml, "sigmabplus_sigmabplus_J2mneg2", sigmabplussigmabplus[4]);
811  write(xml, "sigmabplus_sigmabplus_J1m1", sigmabplussigmabplus[5]);
812  write(xml, "sigmabplus_sigmabplus_J1m0", sigmabplussigmabplus[6]);
813  write(xml, "sigmabplus_sigmabplus_J1mneg1", sigmabplussigmabplus[7]);
814  write(xml, "sigmabplus_sigmabplus_J0m0", sigmabplussigmabplus[8]);
815 
816  write(xml, "sigmabplus_sigmabzero_J2m2", sigmabplussigmabzero[0]);
817  write(xml, "sigmabplus_sigmabzero_J2m1", sigmabplussigmabzero[1]);
818  write(xml, "sigmabplus_sigmabzero_J2m0", sigmabplussigmabzero[2]);
819  write(xml, "sigmabplus_sigmabzero_J2mneg1", sigmabplussigmabzero[3]);
820  write(xml, "sigmabplus_sigmabzero_J2mneg2", sigmabplussigmabzero[4]);
821  write(xml, "sigmabplus_sigmabzero_J1m1", sigmabplussigmabzero[5]);
822  write(xml, "sigmabplus_sigmabzero_J1m0", sigmabplussigmabzero[6]);
823  write(xml, "sigmabplus_sigmabzero_J1mneg1", sigmabplussigmabzero[7]);
824  write(xml, "sigmabplus_sigmabzero_J0m0", sigmabplussigmabzero[8]);
825 
826  write(xml, "sigmabplus_sigmabminus_J2m2", sigmabplussigmabminus[0]);
827  write(xml, "sigmabplus_sigmabminus_J2m1", sigmabplussigmabminus[1]);
828  write(xml, "sigmabplus_sigmabminus_J2m0", sigmabplussigmabminus[2]);
829  write(xml, "sigmabplus_sigmabminus_J2mneg1", sigmabplussigmabminus[3]);
830  write(xml, "sigmabplus_sigmabminus_J2mneg2", sigmabplussigmabminus[4]);
831  write(xml, "sigmabplus_sigmabminus_J1m1", sigmabplussigmabminus[5]);
832  write(xml, "sigmabplus_sigmabminus_J1m0", sigmabplussigmabminus[6]);
833  write(xml, "sigmabplus_sigmabminus_J1mneg1", sigmabplussigmabminus[7]);
834  write(xml, "sigmabplus_sigmabminus_J0m0", sigmabplussigmabminus[8]);
835 
836  write(xml, "sigmabplus_antisigmabminus_J2m2", sigmabplussigmabarbminus[0]);
837  write(xml, "sigmabplus_antisigmabminus_J2m1", sigmabplussigmabarbminus[1]);
838  write(xml, "sigmabplus_antisigmabminus_J2m0", sigmabplussigmabarbminus[2]);
839  write(xml, "sigmabplus_antisigmabminus_J2mneg1", sigmabplussigmabarbminus[3]);
840  write(xml, "sigmabplus_antisigmabminus_J2mneg2", sigmabplussigmabarbminus[4]);
841  write(xml, "sigmabplus_antisigmabminus_J1m1", sigmabplussigmabarbminus[5]);
842  write(xml, "sigmabplus_antisigmabminus_J1m0", sigmabplussigmabarbminus[6]);
843  write(xml, "sigmabplus_antisigmabminus_J1mneg1", sigmabplussigmabarbminus[7]);
844  write(xml, "sigmabplus_antisigmabminus_J0m0", sigmabplussigmabarbminus[8]);
845 
846  write(xml, "sigmabplus_lambda_J1m1", sigmabpluslambdabJ1m1);
847  write(xml, "sigmabplus_lambda_J1m0", sigmabpluslambdabJ1m0);
848  write(xml, "sigmabplus_lambda_J1mneg1", sigmabpluslambdabJ1mneg1);
849 
850  write(xml, "sigmabplus_lambda_transition_J1m1", sigmabpluslambdabtransitionJ1m1);
851  write(xml, "sigmabplus_lambda_transition_J1m0", sigmabpluslambdabtransitionJ1m0);
852  write(xml, "sigmabplus_lambda_transition_J1mneg1", sigmabpluslambdabtransitionJ1mneg1);
853 
854  // Meson-baryon
855  write(xml, "bu_lambdab_J0.5m0.5", bulambdabJhalfmhalf);
856  write(xml, "bu_lambdab_J0.5mneg0.5", bulambdabJhalfmneghalf);
857 
858  write(xml, "bu_sigmabplus_J1.5m1.5", busigmabplus[0]);
859  write(xml, "bu_sigmabplus_J1.5m0.5", busigmabplus[1]);
860  write(xml, "bu_sigmabplus_J1.5mneg0.5", busigmabplus[2]);
861  write(xml, "bu_sigmabplus_J1.5mneg1.5", busigmabplus[3]);
862  write(xml, "bu_sigmabplus_J0.5m0.5", busigmabplus[4]);
863  write(xml, "bu_sigmabplus_J0.5mneg0.5", busigmabplus[5]);
864 
865  write(xml, "bd_sigmabplus_J1.5m1.5", bdsigmabplus[0]);
866  write(xml, "bd_sigmabplus_J1.5m0.5", bdsigmabplus[1]);
867  write(xml, "bd_sigmabplus_J1.5mneg0.5", bdsigmabplus[2]);
868  write(xml, "bd_sigmabplus_J1.5mneg1.5", bdsigmabplus[3]);
869  write(xml, "bd_sigmabplus_J0.5m0.5", bdsigmabplus[4]);
870  write(xml, "bd_sigmabplus_J0.5mneg0.5", bdsigmabplus[5]);
871 
872  write(xml, "bbard_sigmabplus_J1.5m1.5", bbardsigmabplus[0]);
873  write(xml, "bbard_sigmabplus_J1.5m0.5", bbardsigmabplus[1]);
874  write(xml, "bbard_sigmabplus_J1.5mneg0.5", bbardsigmabplus[2]);
875  write(xml, "bbard_sigmabplus_J1.5mneg1.5", bbardsigmabplus[3]);
876  write(xml, "bbard_sigmabplus_J0.5m0.5", bbardsigmabplus[4]);
877  write(xml, "bbard_sigmabplus_J0.5mneg0.5", bbardsigmabplus[5]);
878 
879  pop(xml);
880 
881  END_CODE();
882  }
883 
884 
885 
886 
887 
888 
889  multiNd<DComplex> HBQfunc(const multi1d<ColorMatrix>& HQ)
890  {
891  int length = HQ.size();
892  multiNd<DComplex> result;
893  multi1d<int> Barray;
894  Barray.resize(4); Barray= Nc; Barray[0]=length;
895  result.resize(Barray);
896  result = 0;
897  DComplex thisres;
898  for (int t=0; t<length; t++){
899  Barray[0]=t;
900  for (int aa=0; aa<Nc; aa++){
901  Barray[1]=aa;
902  for (int bb=0; bb<Nc; bb++){
903  Barray[2]=bb;
904  for (int c=0; c<Nc; c++){
905  Barray[3]=c;
906  thisres = 0;
907  for (int cc=0;cc<Nc; cc++) // summed colour index
908  thisres += antiSymTensor3d(aa,bb,cc) * peekColor(HQ[t],cc,c);
909  result[Barray] = thisres;
910  }
911  }
912  }
913  }
914  return result;
915  }
916 
917  multiNd<DComplex> antiHBQfunc(const multi1d<ColorMatrix>& HQ)
918  {
919  int length = HQ.size();
920  multiNd<DComplex> result;
921  multi1d<int> Barray;
922  Barray.resize(4); Barray= Nc; Barray[0]=length;
923  result.resize(Barray);
924  result = 0;
925  DComplex thisres;
926  for (int t=0; t<length; t++){
927  Barray[0]=t;
928  for (int aa=0; aa<Nc; aa++){
929  Barray[1]=aa;
930  for (int bb=0; bb<Nc; bb++){
931  Barray[2]=bb;
932  for (int c=0; c<Nc; c++){
933  Barray[3]=c;
934  thisres = 0;
935  for (int cc=0;cc<Nc; cc++) // summed colour index
936  thisres += antiSymTensor3d(aa,bb,cc) * conj(peekColor(HQ[t],c,cc)); // conjugate transpose
937  result[Barray] = thisres;
938  }
939  }
940  }
941  }
942  return result;
943  }
944 
945 
946 
947  multi1d<DComplex> c1contract(const QllBlock& BzU1zD1z0zCg5, const QllBlock& BzU1zD3z0zCg5,
948  const QllBlock& BzU2zD2zRzCg5, const QllBlock& BzU2zD4zRzCg5,
949  const QllBlock& BzU3zD1z0zCg5, const QllBlock& BzU3zD3z0zCg5,
950  const QllBlock& BzU4zD2zRzCg5, const QllBlock& BzU4zD4zRzCg5,
951  const SpinMatrix& S1, const SpinMatrix& S2)
952  // Contractions for Lambda_b (R) \Lambda_b (0)
953  // spin matrix S1 is at R and S2 is at 0
954  //
955 
956  {
957  int length=BzU1zD1z0zCg5.length();
958  multi1d<DComplex> result;
959  result.resize(length);
960 
961  SpinMatrix S1tilde, S2tilde; // \gamma_0 S_i^\dagger\gamma_0
962 
963  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
964  S2tilde = Gamma(8) * adj(S2) * Gamma(8);
965 
966  DComplex tmpSpin, c2res, c1res, c3res;
967 
968  result = 0;
969  for (int t=0; t<length; t++){
970  for (int is=0; is<Nd; is++){ // summed spin index
971  for (int js=0; js<Nd; js++){ // summed spin index
972  for (int ks=0; ks<Nd; ks++){ // summed spin index
973  for (int ls=0; ls<Nd; ls++){ // summed spin index
974  tmpSpin = peekSpin(S1tilde,is,js) * peekSpin(S2tilde,ks,ls);
975  c1res=0;
976  for (int ic=0; ic<Nc; ic++){ // summed colour index
977  for (int kc=0; kc<Nc; kc++){ // summed colour index
978  for (int lc=0; lc<Nc; lc++){ // summed colour index
979  if (ic != kc && ic != lc && kc != lc)
980  {
981  c2res =0;
982  for (int jc=0; jc<Nc; jc++){ // summed colour index
983  for (int mc=0; mc<Nc; mc++){ // summed colour index
984  for (int nc=0; nc<Nc; nc++){ // summed colour index
985  c3res =0;
986  if (jc != mc && jc != nc && mc != nc)
987  {
988  c3res = -BzU1zD1z0zCg5(t,jc,mc,nc,ks,ls)*BzU2zD2zRzCg5(t,ic,kc,lc,is,js) +
989  BzU1zD3z0zCg5(t,jc,kc,nc,ks,js)*BzU2zD4zRzCg5(t,ic,mc,lc,is,ls) +
990  BzU3zD1z0zCg5(t,ic,mc,nc,is,ls)*BzU4zD2zRzCg5(t,jc,kc,lc,ks,js) -
991  BzU3zD3z0zCg5(t,ic,kc,nc,is,js)*BzU4zD4zRzCg5(t,jc,mc,lc,ks,ls);
992 
993 
994  c3res *= antiSymTensor3d(jc,mc,nc);
995  }
996  c2res += c3res;
997  } // colour nc
998  } //colour mc
999  } // colour jc
1000  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1001  } // if epsilon
1002  } // colour lc
1003  } // colour kc
1004  } // colour ic
1005  result[t] += tmpSpin * c1res;
1006  } // spin
1007  } // spin
1008  } // spin
1009  } // spin
1010  } // t
1011 
1012  return result;
1013  }
1014 
1015 
1016 
1017  multi1d<DComplex> c4contract(const QllBlock& BzU1zD1z0zCgjj, const QllBlock& BzU2zU2zRzCgii,
1018  const QllBlock& BzU2zU4zRzCgii, const QllBlock& BzU3zD1z0zCgjj,
1019  const QllBlock& BzU4zU2zRzCgii,
1020  const SpinMatrix& S1, const SpinMatrix& S2)
1021  // Contractions for Sigma_b^+ (R) \Sigma_b^0 (0)
1022  // and Sigma_b^+ (R) \Lambda_b (0)
1023  // spin matrix S1 is at R and S2 is at 0
1024  //
1025 
1026  {
1027  int length=BzU1zD1z0zCgjj.length();
1028  multi1d<DComplex> result;
1029  result.resize(length);
1030 
1031  SpinMatrix S1tilde, S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1032 
1033  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
1034  S2tilde = Gamma(8) * adj(S2) * Gamma(8);
1035 
1036  DComplex tmpSpin, c2res, c1res, c3res;
1037 
1038  result = 0;
1039  for (int t=0; t<length; t++){
1040  for (int is=0; is<Nd; is++){ // summed spin index
1041  for (int js=0; js<Nd; js++){ // summed spin index
1042  for (int ks=0; ks<Nd; ks++){ // summed spin index
1043  for (int ls=0; ls<Nd; ls++){ // summed spin index
1044  tmpSpin = peekSpin(S1tilde,is,js) * peekSpin(S2tilde,ks,ls);
1045  c1res=0;
1046  for (int ic=0; ic<Nc; ic++){ // summed colour index
1047  for (int kc=0; kc<Nc; kc++){ // summed colour index
1048  for (int lc=0; lc<Nc; lc++){ // summed colour index
1049  if (ic != kc && ic != lc && kc != lc)
1050  {
1051  c2res =0;
1052  for (int jc=0; jc<Nc; jc++){ // summed colour index
1053  for (int mc=0; mc<Nc; mc++){ // summed colour index
1054  for (int nc=0; nc<Nc; nc++){ // summed colour index
1055  c3res =0;
1056  if (jc != mc && jc != nc && mc != nc)
1057  {
1058  c3res = -BzU1zD1z0zCgjj(t,jc,mc,nc,ks,ls)*BzU2zU2zRzCgii(t,ic,kc,lc,is,js) +
1059  BzU1zD1z0zCgjj(t,jc,mc,nc,ks,ls)*BzU2zU2zRzCgii(t,kc,ic,lc,js,is) -
1060  BzU2zU4zRzCgii(t,kc,jc,lc,js,ks)*BzU3zD1z0zCgjj(t,ic,mc,nc,is,ls) +
1061  BzU2zU4zRzCgii(t,ic,jc,lc,is,ks)*BzU3zD1z0zCgjj(t,kc,mc,nc,js,ls) -
1062  BzU3zD1z0zCgjj(t,kc,mc,nc,js,ls)*BzU4zU2zRzCgii(t,jc,ic,lc,ks,is) +
1063  BzU3zD1z0zCgjj(t,ic,mc,nc,is,ls)*BzU4zU2zRzCgii(t,jc,kc,lc,ks,js);
1064 
1065  c3res *= antiSymTensor3d(jc,mc,nc);
1066  }
1067  c2res += c3res;
1068  } // colour nc
1069  } //colour mc
1070  } // colour jc
1071  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1072  } // if epsilon
1073  } // colour lc
1074  } // colour kc
1075  } // colour ic
1076  result[t] += tmpSpin * c1res;
1077  } // spin
1078  } // spin
1079  } // spin
1080  } // spin
1081  } // t
1082 
1083  return result;
1084  }
1085 
1086 
1087  multi1d<DComplex> c5contract(const QllBlock& BzU1zU1z0zCgjj, const QllBlock& BzU1zU3z0zCgjj,
1088  const QllBlock& BzU2zU2zRzCgii, const QllBlock& BzU2zU4zRzCgii,
1089  const QllBlock& BzU3zU1z0zCgjj, const QllBlock& BzU3zU3z0zCgjj,
1090  const QllBlock& BzU4zU2zRzCgii, const QllBlock& BzU4zU4zRzCgii,
1091  const SpinMatrix& S1, const SpinMatrix& S2)
1092  // Contractions for Sigma_b^+ \Sigma_b^+
1093  // spin matrix S1 is at R and S2 is at 0
1094  //
1095 
1096  {
1097  int length=BzU1zU1z0zCgjj.length();
1098  multi1d<DComplex> result;
1099  result.resize(length);
1100 
1101  SpinMatrix S1tilde, S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1102 
1103  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
1104  S2tilde = Gamma(8) * adj(S2) * Gamma(8);
1105 
1106  DComplex tmpSpin, c2res, c1res, c3res;
1107 
1108  result = 0;
1109  for (int t=0; t<length; t++){
1110  for (int is=0; is<Nd; is++){ // summed spin index
1111  for (int js=0; js<Nd; js++){ // summed spin index
1112  for (int ks=0; ks<Nd; ks++){ // summed spin index
1113  for (int ls=0; ls<Nd; ls++){ // summed spin index
1114  tmpSpin = peekSpin(S1tilde,is,js) * peekSpin(S2tilde,ks,ls);
1115  c1res=0;
1116  for (int ic=0; ic<Nc; ic++){ // summed colour index
1117  for (int kc=0; kc<Nc; kc++){ // summed colour index
1118  for (int lc=0; lc<Nc; lc++){ // summed colour index
1119  if (ic != kc && ic != lc && kc != lc)
1120  {
1121  c2res =0;
1122  for (int jc=0; jc<Nc; jc++){ // summed colour index
1123  for (int mc=0; mc<Nc; mc++){ // summed colour index
1124  for (int nc=0; nc<Nc; nc++){ // summed colour index
1125  c3res =0;
1126  if (jc != mc && jc != nc && mc != nc)
1127  {
1128  c3res = - BzU1zU1z0zCgjj(t,jc,mc,nc,ks,ls)*BzU2zU2zRzCgii(t,ic,kc,lc,is,js) +
1129  BzU1zU1z0zCgjj(t,mc,jc,nc,ls,ks)*BzU2zU2zRzCgii(t,ic,kc,lc,is,js) +
1130  BzU1zU1z0zCgjj(t,jc,mc,nc,ks,ls)*BzU2zU2zRzCgii(t,kc,ic,lc,js,is) -
1131  BzU1zU1z0zCgjj(t,mc,jc,nc,ls,ks)*BzU2zU2zRzCgii(t,kc,ic,lc,js,is) -
1132  BzU1zU3z0zCgjj(t,mc,kc,nc,ls,js)*BzU2zU4zRzCgii(t,ic,jc,lc,is,ks) +
1133  BzU1zU3z0zCgjj(t,jc,kc,nc,ks,js)*BzU2zU4zRzCgii(t,ic,mc,lc,is,ls) +
1134  BzU1zU3z0zCgjj(t,mc,ic,nc,ls,is)*BzU2zU4zRzCgii(t,kc,jc,lc,js,ks) -
1135  BzU1zU3z0zCgjj(t,jc,ic,nc,ks,is)*BzU2zU4zRzCgii(t,kc,mc,lc,js,ls) +
1136  BzU2zU4zRzCgii(t,kc,mc,lc,js,ls)*BzU3zU1z0zCgjj(t,ic,jc,nc,is,ks) -
1137  BzU2zU4zRzCgii(t,kc,jc,lc,js,ks)*BzU3zU1z0zCgjj(t,ic,mc,nc,is,ls) -
1138  BzU2zU4zRzCgii(t,ic,mc,lc,is,ls)*BzU3zU1z0zCgjj(t,kc,jc,nc,js,ks) +
1139  BzU2zU4zRzCgii(t,ic,jc,lc,is,ks)*BzU3zU1z0zCgjj(t,kc,mc,nc,js,ls) +
1140  BzU1zU3z0zCgjj(t,mc,kc,nc,ls,js)*BzU4zU2zRzCgii(t,jc,ic,lc,ks,is) -
1141  BzU3zU1z0zCgjj(t,kc,mc,nc,js,ls)*BzU4zU2zRzCgii(t,jc,ic,lc,ks,is) -
1142  BzU1zU3z0zCgjj(t,mc,ic,nc,ls,is)*BzU4zU2zRzCgii(t,jc,kc,lc,ks,js) +
1143  BzU3zU1z0zCgjj(t,ic,mc,nc,is,ls)*BzU4zU2zRzCgii(t,jc,kc,lc,ks,js) -
1144  BzU1zU3z0zCgjj(t,jc,kc,nc,ks,js)*BzU4zU2zRzCgii(t,mc,ic,lc,ls,is) +
1145  BzU3zU1z0zCgjj(t,kc,jc,nc,js,ks)*BzU4zU2zRzCgii(t,mc,ic,lc,ls,is) +
1146  BzU1zU3z0zCgjj(t,jc,ic,nc,ks,is)*BzU4zU2zRzCgii(t,mc,kc,lc,ls,js) -
1147  BzU3zU1z0zCgjj(t,ic,jc,nc,is,ks)*BzU4zU2zRzCgii(t,mc,kc,lc,ls,js) -
1148  BzU3zU3z0zCgjj(t,ic,kc,nc,is,js)*BzU4zU4zRzCgii(t,jc,mc,lc,ks,ls) +
1149  BzU3zU3z0zCgjj(t,kc,ic,nc,js,is)*BzU4zU4zRzCgii(t,jc,mc,lc,ks,ls) +
1150  BzU3zU3z0zCgjj(t,ic,kc,nc,is,js)*BzU4zU4zRzCgii(t,mc,jc,lc,ls,ks) -
1151  BzU3zU3z0zCgjj(t,kc,ic,nc,js,is)*BzU4zU4zRzCgii(t,mc,jc,lc,ls,ks);
1152 
1153 
1154  c3res *= antiSymTensor3d(jc,mc,nc);
1155  }
1156  c2res += c3res;
1157  } // colour nc
1158  } //colour mc
1159  } // colour jc
1160  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1161  } // if epsilon
1162  } // colour lc
1163  } // colour kc
1164  } // colour ic
1165  result[t] += tmpSpin * c1res;
1166  } // spin
1167  } // spin
1168  } // spin
1169  } // spin
1170  } // t
1171 
1172  return result;
1173  }
1174 
1175 
1176 
1177  multi1d<DComplex> c6contract(const QllBlock& BzD1zD1z0zCgjj, const QllBlock& BzU2zU2zRzCgii,
1178  const SpinMatrix& S1, const SpinMatrix& S2)
1179  // Contractions for Sigma_b^+ \Sigma_b^- and Sigma_b^+ \overline{\Sigma}_b^
1180  // spin matrix S1 is at R and S2 is at 0
1181  //
1182 
1183  {
1184  int length=BzD1zD1z0zCgjj.length();
1185  multi1d<DComplex> result;
1186  result.resize(length);
1187 
1188  SpinMatrix S1tilde, S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1189 
1190  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
1191  S2tilde = Gamma(8) * adj(S2) * Gamma(8);
1192 
1193  DComplex tmpSpin, c2res, c1res, c3res;
1194 
1195  result = 0;
1196  for (int t=0; t<length; t++){
1197  for (int is=0; is<Nd; is++){ // summed spin index
1198  for (int js=0; js<Nd; js++){ // summed spin index
1199  for (int ks=0; ks<Nd; ks++){ // summed spin index
1200  for (int ls=0; ls<Nd; ls++){ // summed spin index
1201  tmpSpin = peekSpin(S1tilde,is,js) * peekSpin(S2tilde,ks,ls);
1202  c1res=0;
1203  for (int ic=0; ic<Nc; ic++){ // summed colour index
1204  for (int kc=0; kc<Nc; kc++){ // summed colour index
1205  for (int lc=0; lc<Nc; lc++){ // summed colour index
1206  if (ic != kc && ic != lc && kc != lc)
1207  {
1208  c2res =0;
1209  for (int jc=0; jc<Nc; jc++){ // summed colour index
1210  for (int mc=0; mc<Nc; mc++){ // summed colour index
1211  for (int nc=0; nc<Nc; nc++){ // summed colour index
1212  c3res =0;
1213  if (jc != mc && jc != nc && mc != nc)
1214  {
1215  c3res = -(BzD1zD1z0zCgjj(t,jc,mc,nc,ks,ls)*BzU2zU2zRzCgii(t,ic,kc,lc,is,js)) +
1216  BzD1zD1z0zCgjj(t,mc,jc,nc,ls,ks)*BzU2zU2zRzCgii(t,ic,kc,lc,is,js) +
1217  BzD1zD1z0zCgjj(t,jc,mc,nc,ks,ls)*BzU2zU2zRzCgii(t,kc,ic,lc,js,is) -
1218  BzD1zD1z0zCgjj(t,mc,jc,nc,ls,ks)*BzU2zU2zRzCgii(t,kc,ic,lc,js,is);
1219 
1220  c3res *= antiSymTensor3d(jc,mc,nc);
1221  }
1222  c2res += c3res;
1223  } // colour nc
1224  } //colour mc
1225  } // colour jc
1226  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1227  } // if epsilon
1228  } // colour lc
1229  } // colour kc
1230  } // colour ic
1231  result[t] += tmpSpin * c1res;
1232  } // spin
1233  } // spin
1234  } // spin
1235  } // spin
1236  } // t
1237 
1238  return result;
1239  }
1240 
1241 
1242  multi1d<DComplex> c7contract(const QllBlock& BzU1zD3z0zCG5, const QllBlock& BzU2zU4zRzCGii, const QllBlock& BzU3zD3z0zCG5,
1243  const QllBlock& BzU4zU2zRzCGii, const QllBlock& BzU4zU4zRzCGii,
1244  const SpinMatrix& S1, const SpinMatrix& S2)
1245  // Contractions for Sigma_b^+ \Lambda_b transition
1246  // spin matrix S1 is CG5 and S2 is at CGi
1247  //
1248 
1249  {
1250  int length=BzU1zD3z0zCG5.length();
1251  multi1d<DComplex> result;
1252  result.resize(length);
1253 
1254  SpinMatrix S1tilde, S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1255 
1256  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
1257  S2tilde = Gamma(8) * adj(S2) * Gamma(8);
1258 
1259  DComplex tmpSpin, c2res, c1res, c3res;
1260 
1261  result = 0;
1262  for (int t=0; t<length; t++){
1263  for (int is=0; is<Nd; is++){ // summed spin index
1264  for (int js=0; js<Nd; js++){ // summed spin index
1265  for (int ks=0; ks<Nd; ks++){ // summed spin index
1266  for (int ls=0; ls<Nd; ls++){ // summed spin index
1267  tmpSpin = peekSpin(S1tilde,is,js) * peekSpin(S2tilde,ks,ls);
1268  c1res=0;
1269  for (int ic=0; ic<Nc; ic++){ // summed colour index
1270  for (int kc=0; kc<Nc; kc++){ // summed colour index
1271  for (int lc=0; lc<Nc; lc++){ // summed colour index
1272  if (ic != kc && ic != lc && kc != lc)
1273  {
1274  c2res =0;
1275  for (int jc=0; jc<Nc; jc++){ // summed colour index
1276  for (int mc=0; mc<Nc; mc++){ // summed colour index
1277  for (int nc=0; nc<Nc; nc++){ // summed colour index
1278  c3res =0;
1279  if (jc != mc && jc != nc && mc != nc)
1280  {
1281  c3res =
1282  BzU1zD3z0zCG5(t,mc,kc,nc,ls,js)*BzU2zU4zRzCGii(t,ic,jc,lc,is,ks) -
1283  BzU1zD3z0zCG5(t,jc,kc,nc,ks,js)*BzU2zU4zRzCGii(t,ic,mc,lc,is,ls) -
1284  BzU1zD3z0zCG5(t,mc,kc,nc,ls,js)*BzU4zU2zRzCGii(t,jc,ic,lc,ks,is) +
1285  BzU1zD3z0zCG5(t,jc,kc,nc,ks,js)*BzU4zU2zRzCGii(t,mc,ic,lc,ls,is) +
1286  BzU3zD3z0zCG5(t,ic,kc,nc,is,js)*BzU4zU4zRzCGii(t,jc,mc,lc,ks,ls) -
1287  BzU3zD3z0zCG5(t,ic,kc,nc,is,js)*BzU4zU4zRzCGii(t,mc,jc,lc,ls,ks);
1288 
1289  c3res *= antiSymTensor3d(jc,mc,nc);
1290  }
1291  c2res += c3res;
1292  } // colour nc
1293  } //colour mc
1294  } // colour jc
1295  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1296  } // if epsilon
1297  } // colour lc
1298  } // colour kc
1299  } // colour ic
1300  result[t] += tmpSpin * c1res;
1301  } // spin
1302  } // spin
1303  } // spin
1304  } // spin
1305  } // t
1306 
1307  return result;
1308  }
1309 
1310 
1311  multi1d<DComplex> d1contract(const QllBlock& BzU1zD1z0zCG5, const QllBlock& BzU3zD1z0zCG5,
1312  const HeavyMesonBlock& HzU2zRzG5, const HeavyMesonBlock& HzU4zRzG5,
1313  const SpinMatrix& mesonS1, const SpinMatrix& baryonS2)
1314  // Contractions for Bu (R) Lambda_b (0)
1315  // spin matrix S1 is at R and S2 is at 0
1316  //
1317 
1318  {
1319  int length=BzU1zD1z0zCG5.length();
1320  multi1d<DComplex> result;
1321  result.resize(length);
1322 
1323  SpinMatrix S1tilde,S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1324  S1tilde = mesonS1;
1325  // S1tilde = adj(mesonS1);
1326  S2tilde = Gamma(8) * adj(baryonS2) * Gamma(8);
1327 
1328  DComplex tmpSpin, c2res, c1res;
1329 
1330  result = 0;
1331  for (int t=0; t<length; t++){
1332  for (int is=0; is<Nd; is++){ // summed spin index
1333  for (int js=0; js<Nd; js++){ // summed spin index
1334  for (int ks=0; ks<Nd; ks++){ // summed spin index
1335  for (int ls=0; ls<Nd; ls++){ // summed spin index
1336  tmpSpin = peekSpin(S1tilde,ks,ls) * peekSpin(S2tilde,is,js);
1337  c1res=0;
1338  for (int ic=0; ic<Nc; ic++){ // summed colour index
1339  for (int kc=0; kc<Nc; kc++){ // summed colour index
1340  for (int lc=0; lc<Nc; lc++){ // summed colour index
1341  if (ic != kc && ic != lc && kc != lc)
1342  {
1343  c2res =0;
1344  for (int jc=0; jc<Nc; jc++){ // summed colour index
1345  c2res += -BzU1zD1z0zCG5(t,ic,kc,lc,is,js)*HzU2zRzG5(t,jc,jc,ls,ks) +
1346  BzU3zD1z0zCG5(t,jc,kc,lc,ks,js)*HzU4zRzG5(t,jc,ic,ls,is);
1347 
1348  } // colour jc
1349  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1350  } // if epsilon
1351  } // colour lc
1352  } // colour kc
1353  } // colour ic
1354  result[t] += tmpSpin * c1res;
1355  } // spin
1356  } // spin
1357  } // spin
1358  } // spin
1359  } // t
1360 
1361  return result;
1362  }
1363 
1364 
1365 
1366  multi1d<DComplex> d2contract(const QllBlock& BzU1zU1z0zCGi, const QllBlock& BzU3zU1z0zCGi, const QllBlock& BzU1zU3z0zCGi,
1367  const HeavyMesonBlock& HzU2zRzG5, const HeavyMesonBlock& HzU4zRzG5,
1368  const SpinMatrix& mesonS1, const SpinMatrix& baryonS2)
1369  // Contractions for Bu (R) Sigma_b ^+(0)
1370  // spin matrix S1 is at R and S2 is at 0
1371  //
1372 
1373  {
1374  int length=BzU1zU1z0zCGi.length();
1375  multi1d<DComplex> result;
1376  result.resize(length);
1377 
1378  SpinMatrix S1tilde,S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1379  S1tilde = mesonS1;
1380  S2tilde = Gamma(8) * adj(baryonS2) * Gamma(8);
1381 
1382  DComplex tmpSpin, c2res, c1res;
1383 
1384  result = 0;
1385  for (int t=0; t<length; t++){
1386  for (int is=0; is<Nd; is++){ // summed spin index
1387  for (int js=0; js<Nd; js++){ // summed spin index
1388  for (int ks=0; ks<Nd; ks++){ // summed spin index
1389  for (int ls=0; ls<Nd; ls++){ // summed spin index
1390  tmpSpin = peekSpin(S1tilde,ks,ls) * peekSpin(S2tilde,is,js);
1391  c1res=0;
1392  for (int ic=0; ic<Nc; ic++){ // summed colour index
1393  for (int kc=0; kc<Nc; kc++){ // summed colour index
1394  for (int lc=0; lc<Nc; lc++){ // summed colour index
1395  if (ic != kc && ic != lc && kc != lc)
1396  {
1397  c2res =0;
1398  for (int jc=0; jc<Nc; jc++){ // summed colour index
1399  c2res += -BzU1zU1z0zCGi(t,ic,kc,lc,is,js)*HzU2zRzG5(t,jc,jc,ls,ks) +
1400  BzU1zU1z0zCGi(t,kc,ic,lc,js,is)*HzU2zRzG5(t,jc,jc,ls,ks) -
1401  BzU1zU3z0zCGi(t,kc,jc,lc,js,ks)*HzU4zRzG5(t,jc,ic,ls,is) +
1402  BzU3zU1z0zCGi(t,jc,kc,lc,ks,js)*HzU4zRzG5(t,jc,ic,ls,is) +
1403  BzU1zU3z0zCGi(t,ic,jc,lc,is,ks)*HzU4zRzG5(t,jc,kc,ls,js) -
1404  BzU3zU1z0zCGi(t,jc,ic,lc,ks,is)*HzU4zRzG5(t,jc,kc,ls,js);
1405 
1406  } // colour jc
1407  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1408  } // if epsilon
1409  } // colour lc
1410  } // colour kc
1411  } // colour ic
1412  result[t] += tmpSpin * c1res;
1413  } // spin
1414  } // spin
1415  } // spin
1416  } // spin
1417  } // t
1418 
1419  return result;
1420  }
1421 
1422 
1423  multi1d<DComplex> d3contract(const QllBlock& BzU1zU1z0zCGi, const HeavyMesonBlock& HzD2zRzG5,
1424  const SpinMatrix& mesonS1, const SpinMatrix& baryonS2)
1425  // Contractions for Bd (R) Sigma_b ^+(0)
1426  // spin matrix S1 is at R and S2 is at 0
1427  //
1428 
1429  {
1430  int length=BzU1zU1z0zCGi.length();
1431  multi1d<DComplex> result;
1432  result.resize(length);
1433 
1434  SpinMatrix S1tilde,S2tilde; // \gamma_0 S_i^\dagger\gamma_0
1435  S1tilde = mesonS1;
1436  S2tilde = Gamma(8) * adj(baryonS2) * Gamma(8);
1437 
1438  DComplex tmpSpin, c2res, c1res;
1439 
1440  result = 0;
1441  for (int t=0; t<length; t++){
1442  for (int is=0; is<Nd; is++){ // summed spin index
1443  for (int js=0; js<Nd; js++){ // summed spin index
1444  for (int ks=0; ks<Nd; ks++){ // summed spin index
1445  for (int ls=0; ls<Nd; ls++){ // summed spin index
1446  tmpSpin = peekSpin(S1tilde,ks,ls) * peekSpin(S2tilde,is,js);
1447  c1res=0;
1448  for (int ic=0; ic<Nc; ic++){ // summed colour index
1449  for (int kc=0; kc<Nc; kc++){ // summed colour index
1450  for (int lc=0; lc<Nc; lc++){ // summed colour index
1451  if (ic != kc && ic != lc && kc != lc)
1452  {
1453  c2res =0;
1454  for (int jc=0; jc<Nc; jc++){ // summed colour index
1455  c2res += (BzU1zU1z0zCGi(t,kc,ic,lc,js,is)-BzU1zU1z0zCGi(t,ic,kc,lc,is,js))*HzD2zRzG5(t,jc,jc,ls,ks);
1456  } // colour jc
1457  c1res += c2res * antiSymTensor3d(ic,kc,lc);
1458  } // if epsilon
1459  } // colour lc
1460  } // colour kc
1461  } // colour ic
1462  result[t] += tmpSpin * c1res;
1463  } // spin
1464  } // spin
1465  } // spin
1466  } // spin
1467  } // t
1468 
1469  return result;
1470  }
1471 
1472 
1473 
1474 
1475  multi1d<DComplex> m1contract( const HeavyMesonBlock& HzU1z0zG5, const HeavyMesonBlock& HzU2zRzG5,
1476  const HeavyMesonBlock& HzU3z0zG5, const HeavyMesonBlock& HzU4zRzG5,
1477  const SpinMatrix& S1, const SpinMatrix& S2)
1478  // Contractions for B B potential
1479  // spin matrix S1 is at 0 and S2 is at R
1480  // blocks 1,3 are at R and the others at 0
1481 
1482  {
1483  int length=HzU3z0zG5.length();
1484  multi1d<DComplex> result;
1485  result.resize(length);
1486 
1487  DComplex tmpSpin, c1res;
1488 
1489  result = 0;
1490 
1491  for (int t=0; t<length; t++){
1492  for (int is=0; is<Nd; is++){
1493  for (int js=0; js<Nd; js++){
1494  for (int ks=0; ks<Nd; ks++){
1495  for (int ls=0; ls<Nd; ls++){
1496  tmpSpin = peekSpin(S1,is,js) * peekSpin(S2,ks,ls);
1497  c1res=0;
1498  for (int ic=0; ic<Nc; ic++){
1499  for (int jc=0; jc<Nc; jc++){
1500  c1res += HzU1z0zG5(t,ic,ic,js,is)*HzU2zRzG5(t,jc,jc,ls,ks);
1501  c1res -= HzU3z0zG5(t,ic,jc,js,ks)*HzU4zRzG5(t,jc,ic,ls,is);
1502  }
1503  }
1504  result[t] += tmpSpin * c1res;
1505  }
1506  }
1507  }
1508  }
1509  }
1510  return result;
1511  }
1512 
1513 
1514  multi1d<DComplex> m2contract( const HeavyMesonBlock& HzU1z0zG5, const HeavyMesonBlock& HzD2zRzG5,
1515  const SpinMatrix& S1, const SpinMatrix& S2)
1516  // Contractions for B Bs ( or Bu Bd in I=0) potential
1517  // spin matrix S1 is at 0 and S2 is at R
1518  // blocks 1 are at R and block 2 at 0
1519  {
1520  int length=HzD2zRzG5.length();
1521  multi1d<DComplex> result;
1522  result.resize(length);
1523 
1524  DComplex tmpSpin, c1res;
1525 
1526  result = 0;
1527  for (int t=0; t<length; t++){
1528  for (int is=0; is<Nd; is++){
1529  for (int js=0; js<Nd; js++){
1530  for (int ks=0; ks<Nd; ks++){
1531  for (int ls=0; ls<Nd; ls++){
1532  tmpSpin = peekSpin(S1,is,js) * peekSpin(S2,ks,ls);
1533  c1res=0;
1534  for (int ic=0; ic<Nc; ic++){
1535  for (int jc=0; jc<Nc; jc++){
1536  c1res += HzD2zRzG5(t,jc,jc,ls,ks)*HzU1z0zG5(t,ic,ic,js,is);
1537  }
1538  }
1539  result[t] += tmpSpin * c1res;
1540  }
1541  }
1542  }
1543  }
1544  }
1545  return result;
1546  }
1547 
1548 
1549 
1550 
1551 
1552  multi1d<DComplex> lambdabcontract( const QllBlock& BzU1zD1z0zCG5,
1553  const SpinMatrix& S1)
1554  // Contractions for lambda_B baryon
1555  {
1556  int length=BzU1zD1z0zCG5.length();
1557  multi1d<DComplex> result;
1558  result.resize(length);
1559 
1560  SpinMatrix S1tilde; // \gamma_0 S_i^\dagger\gamma_0
1561 
1562  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
1563 
1564  DComplex tmpSpin, c1res;
1565 
1566  result = 0;
1567  for (int t=0; t<length; t++){
1568  result[t] =0;
1569  for (int is=0; is<Nd; is++){
1570  for (int js=0; js<Nd; js++){
1571  tmpSpin = peekSpin(S1tilde,is,js);
1572  c1res=0;
1573  for (int ic=0; ic<Nc; ic++){
1574  for (int kc=0; kc<Nc; kc++){
1575  for (int lc=0; lc<Nc; lc++){
1576  if (ic != kc && ic != lc && kc != lc)
1577  {
1578  c1res += BzU1zD1z0zCG5(t, ic, kc, lc, is, js) * antiSymTensor3d(ic,kc,lc);
1579  } // if epsilon
1580  } // colour lc
1581  } // colour kc
1582  } // colour ic
1583  result[t] += tmpSpin * c1res;
1584  } // spin js
1585  } // spin is
1586  } //t
1587  return result;
1588  }
1589 
1590  multi1d<DComplex> sigmabpluscontract( const QllBlock& BzU1zU1z0zCGi,
1591  const SpinMatrix& S1)
1592  // Contractions for lambda_B baryon
1593  {
1594  int length=BzU1zU1z0zCGi.length();
1595  multi1d<DComplex> result;
1596  result.resize(length);
1597 
1598  SpinMatrix S1tilde; // \gamma_0 S_i^\dagger\gamma_0
1599 
1600  S1tilde = Gamma(8) * adj(S1) * Gamma(8);
1601 
1602  DComplex tmpSpin, c1res;
1603 
1604  result = 0;
1605  for (int t=0; t<length; t++){
1606  for (int is=0; is<Nd; is++){
1607  for (int js=0; js<Nd; js++){
1608  tmpSpin = peekSpin(S1tilde,is,js);
1609  c1res=0;
1610  for (int ic=0; ic<Nc; ic++){
1611  for (int kc=0; kc<Nc; kc++){
1612  for (int lc=0; lc<Nc; lc++){
1613  if (ic != kc && ic != lc && kc != lc)
1614  {
1615  c1res += (-BzU1zU1z0zCGi(t,ic,kc,lc,is,js) + BzU1zU1z0zCGi(t,kc,ic,lc,js,is))
1616  * antiSymTensor3d(ic,kc,lc);
1617  } // if epsilon
1618  } // colour lc
1619  } // colour kc
1620  } // colour ic
1621  result[t] += tmpSpin * c1res;
1622  } // spin js
1623  } // spin is
1624  } //t
1625  return result;
1626  }
1627 
1628  multi1d<DComplex> bcontract( const HeavyMesonBlock& H1,
1629  const SpinMatrix& S1)
1630  // Contractions for B meson
1631  {
1632  int length=H1.length();
1633  int size= H1.size();
1634  multi1d<DComplex> result;
1635  result.resize(length);
1636 
1637  DComplex tmpSpin, c1res;
1638 
1639  result = 0;
1640  for (int t=0; t<length; t++){
1641  for (int is=0; is<Nd; is++){
1642  for (int js=0; js<Nd; js++){
1643  tmpSpin = peekSpin(S1,is,js);
1644  c1res=0;
1645  for (int ic=0; ic<Nc; ic++){
1646  c1res += H1(t,ic,ic,js,is) ;
1647  }
1648  result[t] += tmpSpin * c1res;
1649  }
1650  }
1651  }
1652  return result;
1653  }
1654 
1655 
1656 
1657  multi2d<DComplex> c4J2corr(
1658  // Spin up blocks
1659  const QllBlock& BzU1zD1z0zCgplus, const QllBlock& BzU2zU2zRzCgplus,
1660  const QllBlock& BzU2zU4zRzCgplus, const QllBlock& BzU3zD1z0zCgplus,
1661  const QllBlock& BzU4zU2zRzCgplus,
1662  // Spin zero blocks
1663  const QllBlock& BzU1zD1z0zCg3, const QllBlock& BzU2zU2zRzCg3,
1664  const QllBlock& BzU2zU4zRzCg3, const QllBlock& BzU3zD1z0zCg3,
1665  const QllBlock& BzU4zU2zRzCg3,
1666  // Spin down blocks
1667  const QllBlock& BzU1zD1z0zCgminus, const QllBlock& BzU2zU2zRzCgminus,
1668  const QllBlock& BzU2zU4zRzCgminus, const QllBlock& BzU3zD1z0zCgminus,
1669  const QllBlock& BzU4zU2zRzCgminus)
1670  {
1671  // Calculate all 9 spin states of the J=1 \otimes J=1 baryon system
1672 
1673  SpinMatrix g_one = 1.0;
1674  SpinMatrix Cg5 = g_one * Gamma(5);
1675  SpinMatrix Cg3 = - g_one * Gamma(14);
1676  SpinMatrix Cgminus = g_one * Gamma(11) + timesI(g_one * Gamma(8));
1677  SpinMatrix Cgplus = g_one * Gamma(11) + timesMinusI(g_one * Gamma(8));
1678 
1679  int length = BzU1zD1z0zCgplus.length();
1680  multi2d<DComplex> c4J2result;
1681  c4J2result.resize(9,length);
1682 
1683  multi1d<DComplex> cJ2m2(length), cJ2m1(length), cJ2m0(length), cJ2mneg1(length), cJ2mneg2(length);
1684  multi1d<DComplex> cJ1m1(length), cJ1m0(length), cJ1mneg1(length), cJ0m0(length);
1685 
1686  // Below cIJ_KL correspons to spin
1687  // I at 0,0
1688  // J at r,0
1689  // K at x,0
1690  // L at x+r,0
1691 
1692 
1693 
1694  // m=2
1695  multi1d<DComplex> c11_11(length);
1696  c11_11 = c4contract(BzU1zD1z0zCgplus, BzU2zU2zRzCgplus,
1697  BzU2zU4zRzCgplus, BzU3zD1z0zCgplus,
1698  BzU4zU2zRzCgplus,
1699  Cgplus, Cgplus);
1700 
1701 
1702  // m=1
1703  multi1d<DComplex> c10_10(length), c10_01(length), c01_10(length), c01_01(length);
1704  c10_10 = c4contract(BzU1zD1z0zCgplus, BzU2zU2zRzCg3,
1705  BzU2zU4zRzCg3, BzU3zD1z0zCgplus,
1706  BzU4zU2zRzCg3,
1707  Cg3, Cgplus);
1708  c10_01 = c4contract(BzU1zD1z0zCgplus, BzU2zU2zRzCg3,
1709  BzU2zU4zRzCg3, BzU3zD1z0zCgplus,
1710  BzU4zU2zRzCg3,
1711  Cgplus, Cg3);
1712  c01_10 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCgplus,
1713  BzU2zU4zRzCgplus, BzU3zD1z0zCg3,
1714  BzU4zU2zRzCgplus,
1715  Cg3, Cgplus);
1716  c01_01 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCgplus,
1717  BzU2zU4zRzCgplus, BzU3zD1z0zCg3,
1718  BzU4zU2zRzCgplus,
1719  Cgplus, Cg3);
1720 
1721 
1722  // m=0
1723  multi1d<DComplex> c1m1_1m1(length), c1m1_00(length), c1m1_m11(length);
1724  multi1d<DComplex> c00_1m1(length), c00_00(length), c00_m11(length);
1725  multi1d<DComplex> cm11_1m1(length), cm11_00(length), cm11_m11(length);
1726  c1m1_1m1 = c4contract(BzU1zD1z0zCgplus, BzU2zU2zRzCgminus,
1727  BzU2zU4zRzCgminus, BzU3zD1z0zCgplus,
1728  BzU4zU2zRzCgminus,
1729  Cgminus, Cgplus);
1730  c1m1_00 = c4contract(BzU1zD1z0zCgplus, BzU2zU2zRzCgminus,
1731  BzU2zU4zRzCgminus, BzU3zD1z0zCgplus,
1732  BzU4zU2zRzCgminus,
1733  Cg3, Cg3);
1734  c1m1_m11 = c4contract(BzU1zD1z0zCgplus, BzU2zU2zRzCgminus,
1735  BzU2zU4zRzCgminus, BzU3zD1z0zCgplus,
1736  BzU4zU2zRzCgminus,
1737  Cgplus, Cgminus);
1738  c00_1m1 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCg3,
1739  BzU2zU4zRzCg3, BzU3zD1z0zCg3,
1740  BzU4zU2zRzCg3,
1741  Cgminus, Cgplus);
1742  c00_00 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCg3,
1743  BzU2zU4zRzCg3, BzU3zD1z0zCg3,
1744  BzU4zU2zRzCg3,
1745  Cg3, Cg3);
1746  c00_m11 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCg3,
1747  BzU2zU4zRzCg3, BzU3zD1z0zCg3,
1748  BzU4zU2zRzCg3,
1749  Cgplus, Cgminus);
1750  cm11_1m1 = c4contract(BzU1zD1z0zCgminus, BzU2zU2zRzCgplus,
1751  BzU2zU4zRzCgplus, BzU3zD1z0zCgminus,
1752  BzU4zU2zRzCgplus,
1753  Cgminus, Cgplus);
1754  cm11_00 = c4contract(BzU1zD1z0zCgminus, BzU2zU2zRzCgplus,
1755  BzU2zU4zRzCgplus, BzU3zD1z0zCgminus,
1756  BzU4zU2zRzCgplus,
1757  Cg3, Cg3);
1758  cm11_m11 = c4contract(BzU1zD1z0zCgminus, BzU2zU2zRzCgplus,
1759  BzU2zU4zRzCgplus, BzU3zD1z0zCgminus,
1760  BzU4zU2zRzCgplus,
1761  Cgplus, Cgminus);
1762 
1763 
1764 
1765  // m=-1
1766  multi1d<DComplex> cm10_m10(length), cm10_0m1(length), c0m1_m10(length), c0m1_0m1(length);
1767  cm10_m10 = c4contract(BzU1zD1z0zCgminus, BzU2zU2zRzCg3,
1768  BzU2zU4zRzCg3, BzU3zD1z0zCgminus,
1769  BzU4zU2zRzCg3,
1770  Cg3, Cgminus);
1771  cm10_0m1 = c4contract(BzU1zD1z0zCgminus, BzU2zU2zRzCg3,
1772  BzU2zU4zRzCg3, BzU3zD1z0zCgminus,
1773  BzU4zU2zRzCg3,
1774  Cgminus, Cg3);
1775  c0m1_m10 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCgminus,
1776  BzU2zU4zRzCgminus, BzU3zD1z0zCg3,
1777  BzU4zU2zRzCgminus,
1778  Cg3, Cgminus);
1779  c0m1_0m1 = c4contract(BzU1zD1z0zCg3, BzU2zU2zRzCgminus,
1780  BzU2zU4zRzCgminus, BzU3zD1z0zCg3,
1781  BzU4zU2zRzCgminus,
1782  Cgminus, Cg3);
1783 
1784 
1785  // m=-2
1786  multi1d<DComplex> cm1m1_m1m1(length);
1787  cm1m1_m1m1 = c4contract(BzU1zD1z0zCgminus, BzU2zU2zRzCgminus,
1788  BzU2zU4zRzCgminus, BzU3zD1z0zCgminus,
1789  BzU4zU2zRzCgminus,
1790  Cgminus, Cgminus);
1791 
1792 
1793 
1794 
1795 
1796 
1797  for (int t=0; t < length; t++) {
1798  cJ2m2[t] = c11_11[t];
1799  cJ2m1[t] = (c01_01[t] + c01_10[t] + c10_01[t] + c10_10[t]) /2.;
1800  cJ1m1[t] = (c01_01[t] - c01_10[t] - c10_01[t] + c10_10[t]) /2.;
1801  cJ2m0[t] = (cm11_m11[t] + 2 * cm11_00[t] + cm11_1m1[t]
1802  + 2 * c00_m11[t] + 4 * c00_00[t] + 2 * c00_1m1[t]
1803  + c1m1_m11[t] + 2 * c1m1_00[t] + c1m1_1m1[t])/6.;
1804  cJ1m0[t] = (cm11_m11[t] - cm11_1m1[t] - c1m1_m11[t] + c1m1_1m1[t])/2.;
1805  cJ0m0[t] = (cm11_m11[t] - cm11_00[t] + cm11_1m1[t]
1806  - c00_m11[t] + c00_00[t] - c00_1m1[t]
1807  + c1m1_m11[t] - c1m1_00[t] + c1m1_1m1[t])/3.;
1808  cJ2mneg1[t] = (c0m1_0m1[t] + c0m1_m10[t] + cm10_0m1[t] + cm10_m10[t])/2.;
1809  cJ1mneg1[t] = (c0m1_0m1[t] - c0m1_m10[t] - cm10_0m1[t] + cm10_m10[t])/2.;
1810  cJ2mneg2[t] = cm1m1_m1m1[t];
1811 
1812  c4J2result[0][t] = cJ2m2[t]; // J=2, m=2 state
1813  c4J2result[1][t] = cJ2m1[t]; // J=2, m=1 state
1814  c4J2result[2][t] = cJ2m0[t]; // J=2, m=0 state
1815  c4J2result[3][t] = cJ2mneg1[t]; // J=2, m=-1 state
1816  c4J2result[4][t] = cJ2mneg2[t]; // J=2, m=-2 state
1817  c4J2result[5][t] = cJ1m1[t];; // J=1, m=1 state
1818  c4J2result[6][t] = cJ1m0[t]; // J=1, m=0 state
1819  c4J2result[7][t] = cJ1mneg1[t]; // J=1, m=-1 state
1820  c4J2result[8][t] = cJ0m0[t]; // J=0, m=0 state
1821  }
1822 
1823  return c4J2result;
1824  }
1825 
1826 
1827 
1828 
1829 
1830 
1831  multi2d<DComplex> c5J2corr(
1832  // Spin up blocks
1833  const QllBlock& BzU1zU1z0zCgplus, const QllBlock& BzU1zU3z0zCgplus,
1834  const QllBlock& BzU2zU2zRzCgplus, const QllBlock& BzU2zU4zRzCgplus,
1835  const QllBlock& BzU3zU1z0zCgplus, const QllBlock& BzU3zU3z0zCgplus,
1836  const QllBlock& BzU4zU2zRzCgplus, const QllBlock& BzU4zU4zRzCgplus,
1837  // Spin zero blocks
1838  const QllBlock& BzU1zU1z0zCg3, const QllBlock& BzU1zU3z0zCg3,
1839  const QllBlock& BzU2zU2zRzCg3, const QllBlock& BzU2zU4zRzCg3,
1840  const QllBlock& BzU3zU1z0zCg3, const QllBlock& BzU3zU3z0zCg3,
1841  const QllBlock& BzU4zU2zRzCg3, const QllBlock& BzU4zU4zRzCg3,
1842  // Spin down blocks
1843  const QllBlock& BzU1zU1z0zCgminus, const QllBlock& BzU1zU3z0zCgminus,
1844  const QllBlock& BzU2zU2zRzCgminus, const QllBlock& BzU2zU4zRzCgminus,
1845  const QllBlock& BzU3zU1z0zCgminus, const QllBlock& BzU3zU3z0zCgminus,
1846  const QllBlock& BzU4zU2zRzCgminus, const QllBlock& BzU4zU4zRzCgminus)
1847  {
1848  // Calculate all 9 spin states of the J=1 \otimes J=1 baryon system
1849 
1850  SpinMatrix g_one = 1.0;
1851  SpinMatrix Cg5 = g_one * Gamma(5);
1852  SpinMatrix Cg3 = - g_one * Gamma(14);
1853  SpinMatrix Cgminus = g_one * Gamma(11) + timesI(g_one * Gamma(8));
1854  SpinMatrix Cgplus = g_one * Gamma(11) + timesMinusI(g_one * Gamma(8));
1855 
1856  int length = BzU1zU1z0zCgplus.length();
1857  multi2d<DComplex> c5J2result;
1858  c5J2result.resize(9,length);
1859 
1860  multi1d<DComplex> cJ2m2(length), cJ2m1(length), cJ2m0(length), cJ2mneg1(length), cJ2mneg2(length);
1861  multi1d<DComplex> cJ1m1(length), cJ1m0(length), cJ1mneg1(length), cJ0m0(length);
1862 
1863  // Below cIJ_KL correspons to spin
1864  // I at 0,0
1865  // J at r,0
1866  // K at x,0
1867  // L at x+r,0
1868 
1869  // m=2
1870  multi1d<DComplex> c11_11(length);
1871  c11_11 = c5contract( BzU1zU1z0zCgplus, BzU1zU3z0zCgplus,
1872  BzU2zU2zRzCgplus, BzU2zU4zRzCgplus,
1873  BzU3zU1z0zCgplus, BzU3zU3z0zCgplus,
1874  BzU4zU2zRzCgplus, BzU4zU4zRzCgplus,
1875  Cgplus, Cgplus);
1876 
1877 
1878  // m=1
1879  multi1d<DComplex> c10_10(length), c10_01(length), c01_10(length), c01_01(length);
1880  c10_10 = c5contract( BzU1zU1z0zCgplus, BzU1zU3z0zCgplus,
1881  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1882  BzU3zU1z0zCgplus, BzU3zU3z0zCgplus,
1883  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1884  Cg3, Cgplus);
1885  c10_01 = c5contract( BzU1zU1z0zCgplus, BzU1zU3z0zCgplus,
1886  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1887  BzU3zU1z0zCgplus, BzU3zU3z0zCgplus,
1888  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1889  Cgplus, Cg3);
1890  c01_10 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1891  BzU2zU2zRzCgplus, BzU2zU4zRzCgplus,
1892  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1893  BzU4zU2zRzCgplus, BzU4zU4zRzCgplus,
1894  Cg3, Cgplus);
1895  c01_01 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1896  BzU2zU2zRzCgplus, BzU2zU4zRzCgplus,
1897  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1898  BzU4zU2zRzCgplus, BzU4zU4zRzCgplus,
1899  Cgplus, Cg3);
1900 
1901 
1902  // m=0
1903  multi1d<DComplex> c1m1_1m1(length), c1m1_00(length), c1m1_m11(length);
1904  multi1d<DComplex> c00_1m1(length), c00_00(length), c00_m11(length);
1905  multi1d<DComplex> cm11_1m1(length), cm11_00(length), cm11_m11(length);
1906  c1m1_1m1 = c5contract( BzU1zU1z0zCgplus, BzU1zU3z0zCgplus,
1907  BzU2zU2zRzCgminus, BzU2zU4zRzCgminus,
1908  BzU3zU1z0zCgplus, BzU3zU3z0zCgplus,
1909  BzU4zU2zRzCgminus, BzU4zU4zRzCgminus,
1910  Cgminus, Cgplus);
1911  c1m1_00 = c5contract( BzU1zU1z0zCgplus, BzU1zU3z0zCgplus,
1912  BzU2zU2zRzCgminus, BzU2zU4zRzCgminus,
1913  BzU3zU1z0zCgplus, BzU3zU3z0zCgplus,
1914  BzU4zU2zRzCgminus, BzU4zU4zRzCgminus,
1915  Cg3, Cg3);
1916  c1m1_m11 = c5contract( BzU1zU1z0zCgplus, BzU1zU3z0zCgplus,
1917  BzU2zU2zRzCgminus, BzU2zU4zRzCgminus,
1918  BzU3zU1z0zCgplus, BzU3zU3z0zCgplus,
1919  BzU4zU2zRzCgminus, BzU4zU4zRzCgminus,
1920  Cgplus, Cgminus);
1921  c00_1m1 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1922  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1923  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1924  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1925  Cgminus, Cgplus);
1926  c00_00 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1927  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1928  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1929  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1930  Cg3, Cg3);
1931  c00_m11 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1932  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1933  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1934  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1935  Cgplus, Cgminus);
1936  cm11_1m1 = c5contract( BzU1zU1z0zCgminus, BzU1zU3z0zCgminus,
1937  BzU2zU2zRzCgplus, BzU2zU4zRzCgplus,
1938  BzU3zU1z0zCgminus, BzU3zU3z0zCgminus,
1939  BzU4zU2zRzCgplus, BzU4zU4zRzCgplus,
1940  Cgminus, Cgplus);
1941  cm11_00 = c5contract( BzU1zU1z0zCgminus, BzU1zU3z0zCgminus,
1942  BzU2zU2zRzCgplus, BzU2zU4zRzCgplus,
1943  BzU3zU1z0zCgminus, BzU3zU3z0zCgminus,
1944  BzU4zU2zRzCgplus, BzU4zU4zRzCgplus,
1945  Cg3, Cg3);
1946  cm11_m11 = c5contract( BzU1zU1z0zCgminus, BzU1zU3z0zCgminus,
1947  BzU2zU2zRzCgplus, BzU2zU4zRzCgplus,
1948  BzU3zU1z0zCgminus, BzU3zU3z0zCgminus,
1949  BzU4zU2zRzCgplus, BzU4zU4zRzCgplus,
1950  Cgplus, Cgminus);
1951 
1952 
1953 
1954  // m=-1
1955  multi1d<DComplex> cm10_m10(length), cm10_0m1(length), c0m1_m10(length), c0m1_0m1(length);
1956  cm10_m10 = c5contract( BzU1zU1z0zCgminus, BzU1zU3z0zCgminus,
1957  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1958  BzU3zU1z0zCgminus, BzU3zU3z0zCgminus,
1959  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1960  Cg3, Cgminus);
1961  cm10_0m1 = c5contract( BzU1zU1z0zCgminus, BzU1zU3z0zCgminus,
1962  BzU2zU2zRzCg3, BzU2zU4zRzCg3,
1963  BzU3zU1z0zCgminus, BzU3zU3z0zCgminus,
1964  BzU4zU2zRzCg3, BzU4zU4zRzCg3,
1965  Cgminus, Cg3);
1966  c0m1_m10 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1967  BzU2zU2zRzCgminus, BzU2zU4zRzCgminus,
1968  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1969  BzU4zU2zRzCgminus, BzU4zU4zRzCgminus,
1970  Cg3, Cgminus);
1971  c0m1_0m1 = c5contract( BzU1zU1z0zCg3, BzU1zU3z0zCg3,
1972  BzU2zU2zRzCgminus, BzU2zU4zRzCgminus,
1973  BzU3zU1z0zCg3, BzU3zU3z0zCg3,
1974  BzU4zU2zRzCgminus, BzU4zU4zRzCgminus,
1975  Cgminus, Cg3);
1976 
1977 
1978  // m=-2
1979  multi1d<DComplex> cm1m1_m1m1(length);
1980  cm1m1_m1m1 = c5contract( BzU1zU1z0zCgminus, BzU1zU3z0zCgminus,
1981  BzU2zU2zRzCgminus, BzU2zU4zRzCgminus,
1982  BzU3zU1z0zCgminus, BzU3zU3z0zCgminus,
1983  BzU4zU2zRzCgminus, BzU4zU4zRzCgminus,
1984  Cgminus, Cgminus);
1985 
1986 
1987 
1988 
1989 
1990 
1991  for (int t=0; t < length; t++) {
1992  cJ2m2[t] = c11_11[t];
1993  cJ2m1[t] = (c01_01[t] + c01_10[t] + c10_01[t] + c10_10[t]) /2.;
1994  cJ1m1[t] = (c01_01[t] - c01_10[t] - c10_01[t] + c10_10[t]) /2.;
1995  cJ2m0[t] = (cm11_m11[t] + 2 * cm11_00[t] + cm11_1m1[t]
1996  + 2 * c00_m11[t] + 4 * c00_00[t] + 2 * c00_1m1[t]
1997  + c1m1_m11[t] + 2 * c1m1_00[t] + c1m1_1m1[t])/6.;
1998  cJ1m0[t] = (cm11_m11[t] - cm11_1m1[t] - c1m1_m11[t] + c1m1_1m1[t])/2.;
1999  cJ0m0[t] = (cm11_m11[t] - cm11_00[t] + cm11_1m1[t]
2000  - c00_m11[t] + c00_00[t] - c00_1m1[t]
2001  + c1m1_m11[t] - c1m1_00[t] + c1m1_1m1[t])/3.;
2002  cJ2mneg1[t] = (c0m1_0m1[t] + c0m1_m10[t] + cm10_0m1[t] + cm10_m10[t])/2.;
2003  cJ1mneg1[t] = (c0m1_0m1[t] - c0m1_m10[t] - cm10_0m1[t] + cm10_m10[t])/2.;
2004  cJ2mneg2[t] = cm1m1_m1m1[t];
2005 
2006  c5J2result[0][t] = cJ2m2[t]; // J=2, m=2 state
2007  c5J2result[1][t] = cJ2m1[t]; // J=2, m=1 state
2008  c5J2result[2][t] = cJ2m0[t]; // J=2, m=0 state
2009  c5J2result[3][t] = cJ2mneg1[t]; // J=2, m=-1 state
2010  c5J2result[4][t] = cJ2mneg2[t]; // J=2, m=-2 state
2011  c5J2result[5][t] = cJ1m1[t];; // J=1, m=1 state
2012  c5J2result[6][t] = cJ1m0[t]; // J=1, m=0 state
2013  c5J2result[7][t] = cJ1mneg1[t]; // J=1, m=-1 state
2014  c5J2result[8][t] = cJ0m0[t]; // J=0, m=0 state
2015  }
2016 
2017  return c5J2result;
2018  }
2019 
2020 
2021 
2022 
2023  multi2d<DComplex> c6J2corr(
2024  // Spin up blocks
2025  const QllBlock& BzD1zD1z0zCgplus, const QllBlock& BzU2zU2zRzCgplus,
2026  // Spin zero blocks
2027  const QllBlock& BzD1zD1z0zCg3, const QllBlock& BzU2zU2zRzCg3,
2028  // Spin down blocks
2029  const QllBlock& BzD1zD1z0zCgminus, const QllBlock& BzU2zU2zRzCgminus)
2030  {
2031  // Calculate all 9 spin states of the J=1 \otimes J=1 baryon system
2032 
2033  SpinMatrix g_one = 1.0;
2034  SpinMatrix Cg5 = g_one * Gamma(5);
2035  SpinMatrix Cg3 = - g_one * Gamma(14);
2036  SpinMatrix Cgminus = g_one * Gamma(11) + timesI(g_one * Gamma(8));
2037  SpinMatrix Cgplus = g_one * Gamma(11) + timesMinusI(g_one * Gamma(8));
2038 
2039  int length = BzD1zD1z0zCgplus.length();
2040  multi2d<DComplex> c6J2result;
2041  c6J2result.resize(9,length);
2042 
2043  multi1d<DComplex> cJ2m2(length), cJ2m1(length), cJ2m0(length), cJ2mneg1(length), cJ2mneg2(length);
2044  multi1d<DComplex> cJ1m1(length), cJ1m0(length), cJ1mneg1(length), cJ0m0(length);
2045 
2046  // Below cIJ_KL correspons to spin
2047  // I at 0,0
2048  // J at r,0
2049  // K at x,0
2050  // L at x+r,0
2051 
2052  // m=2
2053  multi1d<DComplex> c11_11(length);
2054  c11_11 = c6contract( BzD1zD1z0zCgplus, BzU2zU2zRzCgplus,
2055  Cgplus, Cgplus);
2056 
2057  // m=1
2058  multi1d<DComplex> c10_10(length), c10_01(length), c01_10(length), c01_01(length);
2059  c10_10 = c6contract( BzD1zD1z0zCgplus, BzU2zU2zRzCg3,
2060  Cg3, Cgplus);
2061  c10_01 = c6contract( BzD1zD1z0zCgplus, BzU2zU2zRzCg3,
2062  Cgplus, Cg3);
2063  c01_10 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCgplus,
2064  Cg3, Cgplus);
2065  c01_01 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCgplus,
2066  Cgplus, Cg3);
2067 
2068  // m=0
2069  multi1d<DComplex> c1m1_1m1(length), c1m1_00(length), c1m1_m11(length);
2070  multi1d<DComplex> c00_1m1(length), c00_00(length), c00_m11(length);
2071  multi1d<DComplex> cm11_1m1(length), cm11_00(length), cm11_m11(length);
2072  c1m1_1m1 = c6contract( BzD1zD1z0zCgplus, BzU2zU2zRzCgminus,
2073  Cgminus, Cgplus);
2074  c1m1_00 = c6contract( BzD1zD1z0zCgplus, BzU2zU2zRzCgminus,
2075  Cg3, Cg3);
2076  c1m1_m11 = c6contract( BzD1zD1z0zCgplus, BzU2zU2zRzCgminus,
2077  Cgplus, Cgminus);
2078  c00_1m1 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCg3,
2079  Cgminus, Cgplus);
2080  c00_00 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCg3,
2081  Cg3, Cg3);
2082  c00_m11 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCg3,
2083  Cgplus, Cgminus);
2084  cm11_1m1 = c6contract( BzD1zD1z0zCgminus, BzU2zU2zRzCgplus,
2085  Cgminus, Cgplus);
2086  cm11_00 = c6contract( BzD1zD1z0zCgminus, BzU2zU2zRzCgplus,
2087  Cg3, Cg3);
2088  cm11_m11 = c6contract( BzD1zD1z0zCgminus, BzU2zU2zRzCgplus,
2089  Cgplus, Cgminus);
2090 
2091  // m=-1
2092  multi1d<DComplex> cm10_m10(length), cm10_0m1(length), c0m1_m10(length), c0m1_0m1(length);
2093  cm10_m10 = c6contract( BzD1zD1z0zCgminus, BzU2zU2zRzCg3,
2094  Cg3, Cgminus);
2095  cm10_0m1 = c6contract( BzD1zD1z0zCgminus, BzU2zU2zRzCg3,
2096  Cgminus, Cg3);
2097  c0m1_m10 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCgminus,
2098  Cg3, Cgminus);
2099  c0m1_0m1 = c6contract( BzD1zD1z0zCg3, BzU2zU2zRzCgminus,
2100  Cgminus, Cg3);
2101 
2102  // m=-2
2103  multi1d<DComplex> cm1m1_m1m1(length);
2104  cm1m1_m1m1 = c6contract( BzD1zD1z0zCgminus, BzU2zU2zRzCgminus,
2105  Cgminus, Cgminus);
2106 
2107 
2108 
2109 
2110 
2111 
2112  for (int t=0; t < length; t++) {
2113  cJ2m2[t] = c11_11[t];
2114  cJ2m1[t] = (c01_01[t] + c01_10[t] + c10_01[t] + c10_10[t]) /2.;
2115  cJ1m1[t] = (c01_01[t] - c01_10[t] - c10_01[t] + c10_10[t]) /2.;
2116  cJ2m0[t] = (cm11_m11[t] + 2 * cm11_00[t] + cm11_1m1[t]
2117  + 2 * c00_m11[t] + 4 * c00_00[t] + 2 * c00_1m1[t]
2118  + c1m1_m11[t] + 2 * c1m1_00[t] + c1m1_1m1[t])/6.;
2119  cJ1m0[t] = (cm11_m11[t] - cm11_1m1[t] - c1m1_m11[t] + c1m1_1m1[t])/2.;
2120  cJ0m0[t] = (cm11_m11[t] - cm11_00[t] + cm11_1m1[t]
2121  - c00_m11[t] + c00_00[t] - c00_1m1[t]
2122  + c1m1_m11[t] - c1m1_00[t] + c1m1_1m1[t])/3.;
2123  cJ2mneg1[t] = (c0m1_0m1[t] + c0m1_m10[t] + cm10_0m1[t] + cm10_m10[t])/2.;
2124  cJ1mneg1[t] = (c0m1_0m1[t] - c0m1_m10[t] - cm10_0m1[t] + cm10_m10[t])/2.;
2125  cJ2mneg2[t] = cm1m1_m1m1[t];
2126 
2127  c6J2result[0][t] = cJ2m2[t]; // J=2, m=2 state
2128  c6J2result[1][t] = cJ2m1[t]; // J=2, m=1 state
2129  c6J2result[2][t] = cJ2m0[t]; // J=2, m=0 state
2130  c6J2result[3][t] = cJ2mneg1[t]; // J=2, m=-1 state
2131  c6J2result[4][t] = cJ2mneg2[t]; // J=2, m=-2 state
2132  c6J2result[5][t] = cJ1m1[t];; // J=1, m=1 state
2133  c6J2result[6][t] = cJ1m0[t]; // J=1, m=0 state
2134  c6J2result[7][t] = cJ1mneg1[t]; // J=1, m=-1 state
2135  c6J2result[8][t] = cJ0m0[t]; // J=0, m=0 state
2136  }
2137 
2138  return c6J2result;
2139  }
2140 
2141 
2142  multi2d<DComplex> d2J32corr(
2143  // Spin up blocks
2144  const QllBlock& BzU1zU1z0zCgplus, const QllBlock& BzU3zU1z0zCgplus, const QllBlock& BzU1zU3z0zCgplus,
2145  const HeavyMesonBlock& HzU2zRzGup, const HeavyMesonBlock& HzU4zRzGup,
2146  // Spin zero blocks
2147  const QllBlock& BzU1zU1z0zCg3, const QllBlock& BzU3zU1z0zCg3, const QllBlock& BzU1zU3z0zCg3,
2148  // Spin down blocks
2149  const QllBlock& BzU1zU1z0zCgminus, const QllBlock& BzU3zU1z0zCgminus, const QllBlock& BzU1zU3z0zCgminus,
2150  const HeavyMesonBlock& HzU2zRzGdown, const HeavyMesonBlock& HzU4zRzGdown)
2151  {
2152  // Calculate all 6 spin states of the J=1 \otimes J=1/2 baryon - meson system
2153 
2154  SpinMatrix g_one = 1.0;
2155  SpinMatrix Cg5 = g_one * Gamma(5);
2156  SpinMatrix Cg3 = - g_one * Gamma(14);
2157  SpinMatrix Cgminus = g_one * Gamma(11) + timesI(g_one * Gamma(8));
2158  SpinMatrix Cgplus = g_one * Gamma(11) + timesMinusI(g_one * Gamma(8));
2159  SpinMatrix Gup = g_one * Gamma(15) + timesMinusI((g_one*Gamma(3))*Gamma(15));
2160  SpinMatrix Gdown = g_one * Gamma(15) + timesI((g_one*Gamma(3))*Gamma(15));
2161 
2162 
2163  int length = BzU1zU1z0zCgplus.length();
2164  multi2d<DComplex> d2J32result;
2165  d2J32result.resize(6,length);
2166 
2167  multi1d<DComplex> cJ32m32(length), cJ32m12(length), cJ32mneg12(length), cJ32mneg32(length);
2168  multi1d<DComplex> cJ12m12(length), cJ12mneg12(length);
2169 
2170  // Below cIJ_KL correspons to spin
2171  // I at 0,0
2172  // J at r,0
2173  // K at x,0
2174  // L at x+r,0
2175 
2176  // notation is spin of spin 1 first then of spin1/2 and I don't keep the /2
2177 
2178 
2179  // m=3/2
2180  multi1d<DComplex> c11_11(length);
2181  c11_11 = d2contract(BzU1zU1z0zCgplus, BzU3zU1z0zCgplus, BzU1zU3z0zCgplus,
2182  HzU2zRzGup, HzU4zRzGup,
2183  Gup, Cgplus);
2184 
2185  // m=1/2
2186  multi1d<DComplex> c01_01(length), c1m1_1m1(length), c1m1_01(length), c01_1m1(length);
2187  c01_01 = d2contract(BzU1zU1z0zCg3, BzU3zU1z0zCg3, BzU1zU3z0zCg3,
2188  HzU2zRzGup, HzU4zRzGup,
2189  Gup, Cg3);
2190  c01_1m1 = d2contract(BzU1zU1z0zCg3, BzU3zU1z0zCg3, BzU1zU3z0zCg3,
2191  HzU2zRzGup, HzU4zRzGup,
2192  Gdown, Cgplus);
2193  c1m1_01 = d2contract(BzU1zU1z0zCgplus, BzU3zU1z0zCgplus, BzU1zU3z0zCgplus,
2194  HzU2zRzGdown, HzU4zRzGdown,
2195  Gup, Cg3);
2196  c1m1_1m1 = d2contract(BzU1zU1z0zCgplus, BzU3zU1z0zCgplus, BzU1zU3z0zCgplus,
2197  HzU2zRzGdown, HzU4zRzGdown,
2198  Gdown, Cgplus);
2199 
2200  // m=-1/2
2201  multi1d<DComplex> c0m1_0m1(length), cm11_m11(length), cm11_0m1(length), c0m1_m11(length);
2202  c0m1_0m1 = d2contract(BzU1zU1z0zCg3, BzU3zU1z0zCg3, BzU1zU3z0zCg3,
2203  HzU2zRzGdown, HzU4zRzGdown,
2204  Gdown, Cg3);
2205  c0m1_m11 = d2contract(BzU1zU1z0zCg3, BzU3zU1z0zCg3, BzU1zU3z0zCg3,
2206  HzU2zRzGdown, HzU4zRzGdown,
2207  Gup, Cgminus);
2208  cm11_0m1 = d2contract(BzU1zU1z0zCgminus, BzU3zU1z0zCgminus, BzU1zU3z0zCgminus,
2209  HzU2zRzGup, HzU4zRzGup,
2210  Gdown, Cg3);
2211  cm11_m11 = d2contract(BzU1zU1z0zCgminus, BzU3zU1z0zCgminus, BzU1zU3z0zCgminus,
2212  HzU2zRzGup, HzU4zRzGup,
2213  Gup, Cgminus);
2214 
2215  // m=-3/2
2216  multi1d<DComplex> cm1m1_m1m1(length);
2217  cm1m1_m1m1 = d2contract(BzU1zU1z0zCgminus, BzU3zU1z0zCgminus, BzU1zU3z0zCgminus,
2218  HzU2zRzGdown, HzU4zRzGdown,
2219  Gdown, Cgminus);
2220 
2221 
2222  for (int t=0; t < length; t++) {
2223  // notation is spin of spin 1 first then of spin1/2 and I don't keep the /2
2224  cJ32m32[t] = c11_11[t];
2225  cJ32m12[t] = (c1m1_1m1[t] + sqrt(2.) * c1m1_01[t] + sqrt(2.) * c01_1m1[t] + 2.* c01_01[t]) /3.;
2226  cJ32mneg12[t] = (2.* c0m1_0m1[t] + sqrt(2.) * c0m1_m11[t] + sqrt(2.) * cm11_0m1[t] + cm11_m11[t]) /3.;
2227  cJ32mneg32[t] = cm1m1_m1m1[t];
2228  cJ12m12[t] = (2. * c1m1_1m1[t] - sqrt(2.) * c1m1_01[t] - sqrt(2.) * c01_1m1[t] + c01_01[t]) /3.;
2229  cJ12mneg12[t] = (c0m1_0m1[t] - sqrt(2.) * c0m1_m11[t] - sqrt(2.) * cm11_0m1[t] + 2. * cm11_m11[t]) /3.;
2230 
2231  d2J32result[0][t] = cJ32m32[t]; // J=3/2, m=3/2 state
2232  d2J32result[1][t] = cJ32m12[t]; // J=3/2, m=1/2 state
2233  d2J32result[2][t] = cJ32mneg12[t]; // J=3/2, m=-1/2 state
2234  d2J32result[3][t] = cJ32mneg32[t]; // J=3/2, m=-3/2 state
2235  d2J32result[4][t] = cJ12m12[t]; // J=1/2, m=1/2 state
2236  d2J32result[5][t] = cJ12mneg12[t];; // J=1/2, m=-1/2 state
2237  }
2238 
2239  return d2J32result;
2240  }
2241 
2242 
2243 
2244  multi2d<DComplex> d3J32corr(
2245  // Spin up blocks
2246  const QllBlock& BzU1zU1z0zCgplus,
2247  const HeavyMesonBlock& HzU2zRzGup,
2248  // Spin zero blocks
2249  const QllBlock& BzU1zU1z0zCg3,
2250  // Spin down blocks
2251  const QllBlock& BzU1zU1z0zCgminus,
2252  const HeavyMesonBlock& HzU2zRzGdown)
2253  {
2254  // Calculate all 6 spin states of the J=1 \otimes J=1/2 baryon - meson system
2255 
2256  SpinMatrix g_one = 1.0;
2257  SpinMatrix Cg5 = g_one * Gamma(5);
2258  SpinMatrix Cg3 = - g_one * Gamma(14);
2259  SpinMatrix Cgminus = g_one * Gamma(11) + timesI(g_one * Gamma(8));
2260  SpinMatrix Cgplus = g_one * Gamma(11) + timesMinusI(g_one * Gamma(8));
2261  SpinMatrix Gup = g_one * Gamma(15) + timesMinusI((g_one*Gamma(3))*Gamma(15));
2262  SpinMatrix Gdown = g_one * Gamma(15) + timesI((g_one*Gamma(3))*Gamma(15));
2263 
2264 
2265  int length = BzU1zU1z0zCgplus.length();
2266  multi2d<DComplex> d3J32result;
2267  d3J32result.resize(6,length);
2268 
2269  multi1d<DComplex> cJ32m32(length), cJ32m12(length), cJ32mneg12(length), cJ32mneg32(length);
2270  multi1d<DComplex> cJ12m12(length), cJ12mneg12(length);
2271 
2272  // Below cIJ_KL correspons to spin
2273  // I at 0,0
2274  // J at r,0
2275  // K at x,0
2276  // L at x+r,0
2277 
2278  // notation is spin of spin 1 first then of spin1/2 and I don't keep the /2
2279 
2280 
2281  // m=3/2
2282  multi1d<DComplex> c11_11(length);
2283  c11_11 = d3contract(BzU1zU1z0zCgplus,
2284  HzU2zRzGup,
2285  Gup, Cgplus);
2286 
2287  // m=1/2
2288  multi1d<DComplex> c01_01(length), c1m1_1m1(length), c1m1_01(length), c01_1m1(length);
2289  c01_01 = d3contract(BzU1zU1z0zCg3,
2290  HzU2zRzGup,
2291  Gup, Cg3);
2292  c01_1m1 = d3contract(BzU1zU1z0zCg3,
2293  HzU2zRzGup,
2294  Gdown, Cgplus);
2295  c1m1_01 = d3contract(BzU1zU1z0zCgplus,
2296  HzU2zRzGdown,
2297  Gup, Cg3);
2298  c1m1_1m1 = d3contract(BzU1zU1z0zCgplus,
2299  HzU2zRzGdown,
2300  Gdown, Cgplus);
2301 
2302  // m=-1/2
2303  multi1d<DComplex> c0m1_0m1(length), cm11_m11(length), cm11_0m1(length), c0m1_m11(length);
2304  c0m1_0m1 = d3contract(BzU1zU1z0zCg3,
2305  HzU2zRzGdown,
2306  Gdown, Cg3);
2307  c0m1_m11 = d3contract(BzU1zU1z0zCg3,
2308  HzU2zRzGdown,
2309  Gup, Cgminus);
2310  cm11_0m1 = d3contract(BzU1zU1z0zCgminus,
2311  HzU2zRzGup,
2312  Gdown, Cg3);
2313  cm11_m11 = d3contract(BzU1zU1z0zCgminus,
2314  HzU2zRzGup,
2315  Gup, Cgminus);
2316 
2317  // m=-3/2
2318  multi1d<DComplex> cm1m1_m1m1(length);
2319  cm1m1_m1m1 = d3contract(BzU1zU1z0zCgminus,
2320  HzU2zRzGdown,
2321  Gdown, Cgminus);
2322 
2323 
2324  for (int t=0; t < length; t++) {
2325  // notation is spin of spin 1 first then of spin1/2 and I don't keep the /2
2326  cJ32m32[t] = c11_11[t];
2327  cJ32m12[t] = (c1m1_1m1[t] + sqrt(2.) * c1m1_01[t] + sqrt(2.) * c01_1m1[t] + 2.* c01_01[t]) /3.;
2328  cJ32mneg12[t] = (2.* c0m1_0m1[t] + sqrt(2.) * c0m1_m11[t] + sqrt(2.) * cm11_0m1[t] + cm11_m11[t]) /3.;
2329  cJ32mneg32[t] = cm1m1_m1m1[t];
2330  cJ12m12[t] = (2. * c1m1_1m1[t] - sqrt(2.) * c1m1_01[t] - sqrt(2.) * c01_1m1[t] + c01_01[t]) /3.;
2331  cJ12mneg12[t] = (c0m1_0m1[t] - sqrt(2.) * c0m1_m11[t] - sqrt(2.) * cm11_0m1[t] + 2. * cm11_m11[t]) /3.;
2332 
2333  d3J32result[0][t] = cJ32m32[t]; // J=3/2, m=3/2 state
2334  d3J32result[1][t] = cJ32m12[t]; // J=3/2, m=1/2 state
2335  d3J32result[2][t] = cJ32mneg12[t]; // J=3/2, m=-1/2 state
2336  d3J32result[3][t] = cJ32mneg32[t]; // J=3/2, m=-3/2 state
2337  d3J32result[4][t] = cJ12m12[t]; // J=1/2, m=1/2 state
2338  d3J32result[5][t] = cJ12mneg12[t];; // J=1/2, m=-1/2 state
2339  }
2340 
2341  return d3J32result;
2342  }
2343 
2344 
2345 }
antisymtensor.h
Compute anti-symmetric tensors.
barQll_w.h
Heavy Baryon (Qll) 2-pt function : Orginos and Savage.
Chroma::SftMom
Fourier transform phase factor support.
Definition: sftmom.h:35
Chroma::SftMom::numSubsets
int numSubsets() const
Number of subsets - length in decay direction.
Definition: sftmom.h:63
Chroma::SftMom::numMom
int numMom() const
Number of momenta.
Definition: sftmom.h:60
Chroma::antiSymTensor3d
int antiSymTensor3d(int i, int j, int k)
Return 3d symmetric tensor.
Definition: antisymtensor.cc:54
Chroma::write
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
Definition: asqtad_fermact_params_s.cc:40
Chroma::HeavyQuarkProp
void HeavyQuarkProp(LatticeColorMatrix &Qprop, const multi1d< LatticeColorMatrix > &u, const multi1d< int > &src_coord, int length, int bc)
Heavy Quark Propagator.
Definition: barQll_w.cc:149
Chroma::QllQllPOT
void QllQllPOT(const multi1d< LatticeColorMatrix > &u, const LatticePropagator &quark1, const LatticePropagator &quark2, const multi1d< int > &src1, const multi1d< int > &src2, const SftMom &phases, XMLWriter &xml, const std::string &xml_group)
Heavy hadron potentials for SU(2) isospin limit.
Definition: heavy_hadron_potentials_w.cc:38
heavy_hadron_potentials_w.h
Potential between 2 heavy hadrons : Detmold.
mesQl_w.h
Heavy light meson (Qlbar) 2-pt function : Orginos and Savage.
t
int t
Definition: meslate.cc:37
Nd
Nd
Definition: meslate.cc:74
Chroma::BaryonSpinMats::Cg5
SpinMatrix Cg5()
C g_5 = C gamma_5 = Gamma(5)
Definition: barspinmat_w.cc:165
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::c1contract
multi1d< DComplex > c1contract(const QllBlock &BzU1zD1z0zCg5, const QllBlock &BzU1zD3z0zCg5, const QllBlock &BzU2zD2zRzCg5, const QllBlock &BzU2zD4zRzCg5, const QllBlock &BzU3zD1z0zCg5, const QllBlock &BzU3zD3z0zCg5, const QllBlock &BzU4zD2zRzCg5, const QllBlock &BzU4zD4zRzCg5, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:947
Chroma::c4contract
multi1d< DComplex > c4contract(const QllBlock &BzU1zD1z0zCgjj, const QllBlock &BzU2zU2zRzCgii, const QllBlock &BzU2zU4zRzCgii, const QllBlock &BzU3zD1z0zCgjj, const QllBlock &BzU4zU2zRzCgii, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:1017
Chroma::G5
int G5
Definition: pbg5p_w.cc:57
Chroma::c5contract
multi1d< DComplex > c5contract(const QllBlock &BzU1zU1z0zCgjj, const QllBlock &BzU1zU3z0zCgjj, const QllBlock &BzU2zU2zRzCgii, const QllBlock &BzU2zU4zRzCgii, const QllBlock &BzU3zU1z0zCgjj, const QllBlock &BzU3zU3z0zCgjj, const QllBlock &BzU4zU2zRzCgii, const QllBlock &BzU4zU4zRzCgii, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:1087
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::bcontract
multi1d< DComplex > bcontract(const HeavyMesonBlock &H1, const SpinMatrix &S1)
Definition: heavy_hadron_potentials_w.cc:1628
Chroma::push
push(xml_out,"Condensates")
Chroma::c
Double c
Definition: invbicg.cc:108
Chroma::c5J2corr
multi2d< DComplex > c5J2corr(const QllBlock &BzU1zU1z0zCgplus, const QllBlock &BzU1zU3z0zCgplus, const QllBlock &BzU2zU2zRzCgplus, const QllBlock &BzU2zU4zRzCgplus, const QllBlock &BzU3zU1z0zCgplus, const QllBlock &BzU3zU3z0zCgplus, const QllBlock &BzU4zU2zRzCgplus, const QllBlock &BzU4zU4zRzCgplus, const QllBlock &BzU1zU1z0zCg3, const QllBlock &BzU1zU3z0zCg3, const QllBlock &BzU2zU2zRzCg3, const QllBlock &BzU2zU4zRzCg3, const QllBlock &BzU3zU1z0zCg3, const QllBlock &BzU3zU3z0zCg3, const QllBlock &BzU4zU2zRzCg3, const QllBlock &BzU4zU4zRzCg3, const QllBlock &BzU1zU1z0zCgminus, const QllBlock &BzU1zU3z0zCgminus, const QllBlock &BzU2zU2zRzCgminus, const QllBlock &BzU2zU4zRzCgminus, const QllBlock &BzU3zU1z0zCgminus, const QllBlock &BzU3zU3z0zCgminus, const QllBlock &BzU4zU2zRzCgminus, const QllBlock &BzU4zU4zRzCgminus)
Definition: heavy_hadron_potentials_w.cc:1831
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::HBQfunc
multiNd< DComplex > HBQfunc(const multi1d< ColorMatrix > &HQ)
Definition: heavy_hadron_potentials_w.cc:889
Chroma::m2contract
multi1d< DComplex > m2contract(const HeavyMesonBlock &HzU1z0zG5, const HeavyMesonBlock &HzD2zRzG5, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:1514
Chroma::d2J32corr
multi2d< DComplex > d2J32corr(const QllBlock &BzU1zU1z0zCgplus, const QllBlock &BzU3zU1z0zCgplus, const QllBlock &BzU1zU3z0zCgplus, const HeavyMesonBlock &HzU2zRzGup, const HeavyMesonBlock &HzU4zRzGup, const QllBlock &BzU1zU1z0zCg3, const QllBlock &BzU3zU1z0zCg3, const QllBlock &BzU1zU3z0zCg3, const QllBlock &BzU1zU1z0zCgminus, const QllBlock &BzU3zU1z0zCgminus, const QllBlock &BzU1zU3z0zCgminus, const HeavyMesonBlock &HzU2zRzGdown, const HeavyMesonBlock &HzU4zRzGdown)
Definition: heavy_hadron_potentials_w.cc:2142
Chroma::c4J2corr
multi2d< DComplex > c4J2corr(const QllBlock &BzU1zD1z0zCgplus, const QllBlock &BzU2zU2zRzCgplus, const QllBlock &BzU2zU4zRzCgplus, const QllBlock &BzU3zD1z0zCgplus, const QllBlock &BzU4zU2zRzCgplus, const QllBlock &BzU1zD1z0zCg3, const QllBlock &BzU2zU2zRzCg3, const QllBlock &BzU2zU4zRzCg3, const QllBlock &BzU3zD1z0zCg3, const QllBlock &BzU4zU2zRzCg3, const QllBlock &BzU1zD1z0zCgminus, const QllBlock &BzU2zU2zRzCgminus, const QllBlock &BzU2zU4zRzCgminus, const QllBlock &BzU3zD1z0zCgminus, const QllBlock &BzU4zU2zRzCgminus)
Definition: heavy_hadron_potentials_w.cc:1657
Chroma::c6contract
multi1d< DComplex > c6contract(const QllBlock &BzD1zD1z0zCgjj, const QllBlock &BzU2zU2zRzCgii, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:1177
Chroma::c6J2corr
multi2d< DComplex > c6J2corr(const QllBlock &BzD1zD1z0zCgplus, const QllBlock &BzU2zU2zRzCgplus, const QllBlock &BzD1zD1z0zCg3, const QllBlock &BzU2zU2zRzCg3, const QllBlock &BzD1zD1z0zCgminus, const QllBlock &BzU2zU2zRzCgminus)
Definition: heavy_hadron_potentials_w.cc:2023
Chroma::sigmabpluscontract
multi1d< DComplex > sigmabpluscontract(const QllBlock &BzU1zU1z0zCGi, const SpinMatrix &S1)
Definition: heavy_hadron_potentials_w.cc:1590
Chroma::m1contract
multi1d< DComplex > m1contract(const HeavyMesonBlock &HzU1z0zG5, const HeavyMesonBlock &HzU2zRzG5, const HeavyMesonBlock &HzU3z0zG5, const HeavyMesonBlock &HzU4zRzG5, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:1475
Chroma::d3contract
multi1d< DComplex > d3contract(const QllBlock &BzU1zU1z0zCGi, const HeavyMesonBlock &HzD2zRzG5, const SpinMatrix &mesonS1, const SpinMatrix &baryonS2)
Definition: heavy_hadron_potentials_w.cc:1423
Chroma::pop
pop(xml_out)
Chroma::c7contract
multi1d< DComplex > c7contract(const QllBlock &BzU1zD3z0zCG5, const QllBlock &BzU2zU4zRzCGii, const QllBlock &BzU3zD3z0zCG5, const QllBlock &BzU4zU2zRzCGii, const QllBlock &BzU4zU4zRzCGii, const SpinMatrix &S1, const SpinMatrix &S2)
Definition: heavy_hadron_potentials_w.cc:1242
Chroma::antiHBQfunc
multiNd< DComplex > antiHBQfunc(const multi1d< ColorMatrix > &HQ)
Definition: heavy_hadron_potentials_w.cc:917
Chroma::d1contract
multi1d< DComplex > d1contract(const QllBlock &BzU1zD1z0zCG5, const QllBlock &BzU3zD1z0zCG5, const HeavyMesonBlock &HzU2zRzG5, const HeavyMesonBlock &HzU4zRzG5, const SpinMatrix &mesonS1, const SpinMatrix &baryonS2)
Definition: heavy_hadron_potentials_w.cc:1311
Chroma::lc
multi1d< LatticeComplex > lc(Ncb)
Chroma::START_CODE
START_CODE()
Chroma::d3J32corr
multi2d< DComplex > d3J32corr(const QllBlock &BzU1zU1z0zCgplus, const HeavyMesonBlock &HzU2zRzGup, const QllBlock &BzU1zU1z0zCg3, const QllBlock &BzU1zU1z0zCgminus, const HeavyMesonBlock &HzU2zRzGdown)
Definition: heavy_hadron_potentials_w.cc:2244
Chroma::d2contract
multi1d< DComplex > d2contract(const QllBlock &BzU1zU1z0zCGi, const QllBlock &BzU3zU1z0zCGi, const QllBlock &BzU1zU3z0zCGi, const HeavyMesonBlock &HzU2zRzG5, const HeavyMesonBlock &HzU4zRzG5, const SpinMatrix &mesonS1, const SpinMatrix &baryonS2)
Definition: heavy_hadron_potentials_w.cc:1366
Chroma::lambdabcontract
multi1d< DComplex > lambdabcontract(const QllBlock &BzU1zD1z0zCG5, const SpinMatrix &S1)
Definition: heavy_hadron_potentials_w.cc:1552
testing::internal::string
::std::string string
Definition: gtest.h:1979
Generated by doxygen 1.9.1