switch (baryons)
XMLArrayWriter xml_sink_mom(xml_bar,num_mom); push(xml_sink_mom, "momenta");  for(int sink_mom_num = 0; sink_mom_num < num_mom; ++sink_mom_num)
void baryon(LatticePropagator& quark_propagator,             const SftMom& phases,	    multi3d<DComplex>& barprop){  START_CODE("baryon");  // Length of lattice in decay direction  int length = phases.numSubsets() ;  if ( Ns != 4 || Nc != 3 )		/* Code is specific to Ns=4 and Nc=3. */    return;  // Setup the return stuff  const int num_baryons = 21;  int num_mom = phases.numMom();  barprop.resize(num_baryons,num_mom,length);  SpinMatrix Cgm;  SpinMatrix Cg4m;  SpinMatrix CgmNR;  SpinMatrix g_one = 1.0;  SpinMatrix g_tmp1;  // C = Gamma(10)  g_tmp1 = 0.5 * (Gamma(2) * g_one  +  timesI(Gamma(1) * g_one));  Cgm = Gamma(10) * g_tmp1;  Cg4m = Gamma(10) * (Gamma(8) * g_tmp1);  CgmNR = Cgm - Cg4m;  // S_proj = (1 + \Sigma_3)*(1 + gamma_4) / 2   //        = (1 + Gamma(8) - i G(3) - i G(11)) / 2  SpinMatrix S_proj =     0.5*((g_one + Gamma(8) * g_one) - timesI(Gamma(3) * g_one  +  Gamma(11) * g_one));  // S_proj_unpol = (1/2)(1 + gamma_4)  SpinMatrix S_proj_unpol = 0.5 * (g_one + (g_one * Gamma(8)));  // C g_5 NR = (1/2)*C gamma_5 * ( 1 + g_4 )  SpinMatrix Cg5NR = Gamma(5) * S_proj_unpol;  // C NR = (1/2)*C * ( 1 + g_4 )  SpinMatrix CNR = Gamma(10) * S_proj_unpol;  LatticeComplex b_prop;  // Loop over baryons  for(int baryons = 0; baryons < num_baryons; ++baryons)  {    LatticePropagator di_quark;    switch (baryons)    {    case 0:      // Proton_1; use also for Lambda_1!      // |P_1, s_z=1/2> = (d C gamma_5 u) "u_up", see comments at top      // C gamma_5 = Gamma(5)      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      di_quark = quarkContract13(quark_propagator * Gamma(5),				 Gamma(5) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	     + trace(S_proj * traceColor(quark_propagator * di_quark));      break;		      case 1:      // Lambda_1 = 3*Proton_1 (for compatibility with heavy-light routine)      // |L_1, s_z=1/2> = 2*(u C gamma_5 d) "s_up" + (s C gamma_5 d) "u_up"      //                  + (u C gamma_5 s) "d_up" , see comments at top         // C gamma_5 = Gamma(5)      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      b_prop *= 3.0;      break;    case 2:      // Delta^+_1      // |D_1, s_z=3/2> = 2*(d C gamma_- u) "u_up" + (u C gamma_- u) "d_up"      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      di_quark = quarkContract13(quark_propagator * Cgm, 				 Cgm * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 3:      // Proton_2; use also for Lambda_2!      // |P_2, s_z=1/2> = (d C gamma_4 gamma_5 u) "u_up"       // C gamma_5 gamma_4 = - Gamma(13)      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      di_quark = quarkContract13(quark_propagator * Gamma(13),				 Gamma(13) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	     + trace(S_proj * traceColor(quark_propagator * di_quark));      break;    case 4:      // Lambda_2 = 3*Proton_2 (for compatibility with heavy-light routine)      // |L_2, s_z=1/2> = 2*(u C gamma_4 gamma_5 d) "s_up"      //                  + (s C gamma_4 gamma_5 d) "u_up"      //                  + (u C gamma_4 gamma_5 s) "d_up"      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      b_prop *= 3.0;      break;    case 5:      // Sigma^{*+}_2      // |D_2, s_z=3/2> = 2*(d C gamma_4 gamma_- u) "u_up"       //                  + (u C gamma_4 gamma_- u) "d_up"       // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      di_quark = quarkContract13(quark_propagator * Cg4m,				 Cg4m * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));      // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 6:      // Proton^+_3; use also for Lambda_3!      // |P_3, s_z=1/2> = (d C (1/2)(1 + gamma_4) gamma_5 u) "u_up"       // C gamma_5 - C gamma_5 gamma_4 = Gamma(5) + Gamma(13)      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      di_quark = quarkContract13(quark_propagator * Cg5NR,  				 Cg5NR * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	     + trace(S_proj * traceColor(quark_propagator * di_quark));      break;    case 7:      // Lambda_3 = 3*Proton_3 (for compatibility with heavy-light routine)      // |L_3, s_z=1/2> = 2*(u C (1/2)(1 + gamma_4) gamma_5 d) "s_up"      //                  + (s C (1/2)(1 + gamma_4) gamma_5 d) "u_up"      //                  + (u C (1/2)(1 + gamma_4) gamma_5 s) "d_up"      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      b_prop *= 3.0;      break;    case 8:      // Sigma^{*+}_3      // |D_3, s_z=3/2> = 2*(d C (1/2)(1 + gamma_4) gamma_- d) u) "u_up"      //                  + (u C (1/2)(1 + gamma_4) gamma_- d) u) "d_up"      // Polarized:      // S_proj = T = (1 + \Sigma_3)*(1 + gamma_4) / 2       //            = (1 + Gamma(8) - i G(3) - i G(11)) / 2      di_quark = quarkContract13(quark_propagator * CgmNR,				 CgmNR * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));      // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 9:      // Proton_4 -- but unpolarised ; use also for Lambda_4!      // |P_4, s_z=1/2> = (d C gamma_5 u) "u_up", see comments at top      // C gamma_5 = Gamma(5)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(5),				 Gamma(5) * quark_propagator);      b_prop = trace(S_proj_unpol * traceColor(quark_propagator * traceSpin(di_quark)))	     + trace(S_proj_unpol * traceColor(quark_propagator * di_quark));      break;    case 10:      // Proton_5; use also for Lambda_5!      // |P_5, s_z=1/2> = (d C gamma_4 gamma_5 u) "u_up", see comments at top      // C gamma_5 gamma_4 = - Gamma(13)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(13),				 Gamma(13) * quark_propagator);      b_prop = trace(S_proj_unpol * traceColor(quark_propagator * traceSpin(di_quark)))	     + trace(S_proj_unpol * traceColor(quark_propagator * di_quark));      break;        case 11:      // Proton^+_6; use also for Lambda_6!      // |P_6, s_z=1/2> = (d C (1/2)(1 + gamma_4) gamma_5 u) "u_up", see comments at top      // C gamma_5 = Gamma(5)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Cg5NR,				 Cg5NR * quark_propagator);      b_prop = trace(S_proj_unpol * traceColor(quark_propagator * traceSpin(di_quark)))	     + trace(S_proj_unpol * traceColor(quark_propagator * di_quark));      break;    case 12:      // Delta_x^+_4 -- unpolarised with explicit gamma_k interpolation      // |D_4, s_z=3/2> = 2*(d C gamma_1 u) "u_up" + (u C gamma_1 u) "d_up"      // C gamma_1 = Gamma(10) * Gamma(1) = Gamma(11)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(11),				 Gamma(11) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 13:      // Delta_y^+_4 -- unpolarised with explicit gamma_k interpolation      // |D_4, s_z=3/2> = 2*(d C gamma_2 u) "u_up" + (u C gamma_2 u) "d_up"      // C gamma_2 = Gamma(10) * Gamma(2) = Gamma(8)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(8),				 Gamma(8) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 14:      // Delta_z^+_4 -- unpolarised with explicit gamma_k interpolation      // |D_4, s_z=3/2> = 2*(d C gamma_3 u) "u_up" + (u C gamma_3 u) "d_up"      // C gamma_3 = Gamma(10) * Gamma(3) = Gamma(11)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(11),				 Gamma(11) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 15:      // Delta_x^+_5 -- unpolarised with explicit gamma_k interpolation      // |D_5, s_z=3/2> = 2*(d C gamma_4 gamma_1 u) "u_up" + (u C gamma_4 gamma_1 u) "d_up"      // C gamma_4 gamma_1 = Gamma(10) * Gamma(8) * Gamma(1) = Gamma(3)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(3),				 Gamma(3) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 16:      // Delta_y^+_5 -- unpolarised with explicit gamma_k interpolation      // |D_4, s_z=3/2> = 2*(d C gamma_4 gamma_2 u) "u_up" + (u C gamma_4 gamma_2 u) "d_up"      // C gamma_4 gamma_2 = Gamma(10) * Gamma(8) * Gamma(2) = Gamma(0)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(0),				 Gamma(0) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 17:      // Delta_z^+_5 -- unpolarised with explicit gamma_k interpolation      // |D_4, s_z=3/2> = 2*(d C gamma_4 gamma_3 u) "u_up" + (u C gamma_4 gamma_3 u) "d_up"      // C gamma_4 gamma_3 = Gamma(10) * Gamma(8) * Gamma(3) = Gamma(1)      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13(quark_propagator * Gamma(1),				 Gamma(1) * quark_propagator);      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 18:      // Delta_x^+_6 -- unpolarised NR with explicit gamma_k interpolation      // |D_6, s_z=3/2> = 2*(d C gamma_1 (1/2)(1 + gamma_4) d) u) "u_up"      //                  + (u C gamma_1 (1/2)(1 + gamma_4) d) u) "d_up"      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13((quark_propagator*Gamma(1))*CNR,				 Gamma(1)*(CNR*quark_propagator));      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 19:      // Delta_y^+_6 -- unpolarised NR with explicit gamma_k interpolation      // |D_6, s_z=3/2> = 2*(d C gamma_2 (1/2)(1 + gamma_4) d) u) "u_up"      //                  + (u C gamma_2 (1/2)(1 + gamma_4) d) u) "d_up"      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13((quark_propagator*Gamma(2))*CNR,				 Gamma(2)*(CNR*quark_propagator));      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    case 20:      // Delta_z^+_6 -- unpolarised NR with explicit gamma_k interpolation      // |D_6, s_z=3/2> = 2*(d C gamma_3 (1/2)(1 + gamma_4) d) u) "u_up"      //                  + (u C gamma_3 (1/2)(1 + gamma_4) d) u) "d_up"      // Unpolarized:      // S_proj_unpol = T = (1/2)(1 + gamma_4)      di_quark = quarkContract13((quark_propagator*Gamma(4))*CNR,				 Gamma(4)*(CNR*quark_propagator));      b_prop = trace(S_proj * traceColor(quark_propagator * traceSpin(di_quark)))	   + 2*trace(S_proj * traceColor(quark_propagator * di_quark));            // Multiply by 3 for compatibility with heavy-light routine      b_prop *= 3.0;      break;    default:      QDP_error_exit("Unknown baryon: baryons=%d",baryons);    }    // Project onto zero and if desired non-zero momentum    multi2d<DComplex> hsum;    hsum = phases.sft(b_prop);    for(int sink_mom_num=0; sink_mom_num < num_mom; ++sink_mom_num)       for(int t = 0; t < length; ++t)      {	// NOTE: there is NO  1/2  multiplying hsum	barprop[baryons][sink_mom_num][t] = hsum[sink_mom_num][t];      }  } // end loop over baryons  END_CODE("baryon");}