Real Mass   = prop_header.Mass;
Real Mass   = prop_header.FermActHandle->getMass();
int main(int argc, char **argv){  // Put the machine into a known state  QDP_initialize(&argc, &argv);  // Input parameter structure  Spectrum_input_t  input;  // Instantiate xml reader for DATA  XMLReader xml_in("DATA");  // Read data  read(xml_in, "/spectrum_w", input);  // Specify lattice size, shape, etc.  Layout::setLattSize(input.param.nrow);  Layout::create();  QDPIO::cout << " SPECTRUM_W: Spectroscopy for Wilson fermions" << endl;  /*   * Sanity checks   */  if (input.param.wvf_param.size() != input.prop.prop_files.size())  {    QDPIO::cerr << "wvf_param size inconsistent with prop_files size" << endl;    QDP_abort(1);  }  QDPIO::cout << endl << "     Gauge group: SU(" << Nc << ")" << endl;  QDPIO::cout << "     volume: " << input.param.nrow[0];  for (int i=1; i<Nd; ++i) {    QDPIO::cout << " x " << input.param.nrow[i];  }  QDPIO::cout << endl;  // Read in the configuration along with relevant information.  multi1d<LatticeColorMatrix> u(Nd);  XMLReader gauge_file_xml, gauge_xml;  switch (input.cfg.cfg_type)   {  case CFG_TYPE_SZIN :    readSzin(gauge_xml, u, input.cfg.cfg_file);    break;  case CFG_TYPE_SZINQIO:    readGauge(gauge_file_xml, gauge_xml, u, input.cfg.cfg_file, QDPIO_SERIAL);    break;  case CFG_TYPE_NERSC:    readArchiv(gauge_xml, u, input.cfg.cfg_file);    break;  default :    QDP_error_exit("Configuration type is unsupported.");  }  // Check if the gauge field configuration is unitarized  unitarityCheck(u);  // Instantiate XML writer for XMLDAT  XMLFileWriter xml_out("XMLDAT");  push(xml_out, "spectrum_w");  proginfo(xml_out);    // Print out basic program info  // Write out the input  write(xml_out, "Input", xml_in);  // Write out the config info  write(xml_out, "Config_info", gauge_xml);  push(xml_out, "Output_version");  write(xml_out, "out_version", 10);  pop(xml_out);  xml_out.flush();  // First calculate some gauge invariant observables just for info.  Double w_plaq, s_plaq, t_plaq, link;  multi1d<DComplex> pollp(Nd);  MesPlq(u, w_plaq, s_plaq, t_plaq, link);  for(int mu = 0; mu < Nd; ++mu)    polylp(u, pollp[mu], mu);  push(xml_out, "Observables");  write(xml_out, "w_plaq", w_plaq);  write(xml_out, "s_plaq", s_plaq);  write(xml_out, "t_plaq", t_plaq);  write(xml_out, "link", link);  write(xml_out, "pollp", pollp);  pop(xml_out);  xml_out.flush();  // Keep an array of all the xml output buffers  XMLArrayWriter xml_array(xml_out,input.prop.prop_files.size());  push(xml_array, "Wilson_hadron_measurements");  // Now loop over the various fermion masses  for (int loop=0; loop < input.prop.prop_files.size(); ++loop)  {    // Read the quark propagator and extract headers    LatticePropagator quark_propagator;    ChromaProp_t prop_header;    PropSource_t source_header;    {      XMLReader prop_file_xml, prop_record_xml;      readQprop(prop_file_xml, 		prop_record_xml, quark_propagator,		input.prop.prop_files[loop], QDPIO_SERIAL);      // Try to invert this record XML into a ChromaProp struct      // Also pull out the id of this source      try      {	read(prop_record_xml, "/Propagator/ForwardProp", prop_header);	read(prop_record_xml, "/Propagator/PropSource", source_header);      }      catch (const string& e)       {	QDPIO::cerr << "Error extracting forward_prop header: " << e << endl;	throw;      }    }    // Derived from input prop    int j_decay = source_header.j_decay;    Real Mass   = prop_header.Mass;    multi1d<int> boundary = prop_header.boundary;    multi1d<int> t_source = source_header.t_source;    // Flags    int t0      = t_source[j_decay];    int bc_spec = boundary[j_decay];    // Initialize the slow Fourier transform phases    SftMom phases(input.param.mom2_max, input.param.avg_equiv_mom, j_decay);    // Keep a copy of the phases with NO momenta    SftMom phases_nomom(0, true, j_decay);    // Next array element - name auto-written    push(xml_array);    write(xml_array, "loop", loop);    write(xml_array, "Mass_mes", Mass);    write(xml_array, "t_source", t_source);    // Save prop input    write(xml_array, "ForwardProp", prop_header);    write(xml_array, "PropSource", source_header);    // Sanity check - write out the norm2 of the forward prop in the j_decay direction    // Use this for any possible verification    {      multi1d<Double> forward_prop_corr = sumMulti(localNorm2(quark_propagator), 						   phases.getSet());      push(xml_array, "Forward_prop_correlator");      write(xml_array, "forward_prop_corr", forward_prop_corr);      pop(xml_array);    }    // Determine what kind of source to use    bool Pt_src = (source_header.source_type == SRC_TYPE_POINT_SOURCE) ? true : false;    bool Sl_src = (source_header.source_type == SRC_TYPE_SHELL_SOURCE) ? true : false;    bool Wl_src = (source_header.source_type == SRC_TYPE_WALL_SOURCE)  ? true : false;    // Do the mesons first    if (input.param.MesonP)     {      // Construct {Point|Shell}-Point mesons, if desired      if (input.param.Pt_snk)       {	if (Pt_src)	  mesons(quark_propagator, quark_propagator, phases, t0,		 xml_array, "Point_Point_Wilson_Mesons");        	if (Sl_src)	  mesons(quark_propagator, quark_propagator, phases, t0,		 xml_array, "Shell_Point_Wilson_Mesons");        	if (Wl_src)	  mesons(quark_propagator, quark_propagator, phases_nomom, t0,		 xml_array, "Wall_Point_Wilson_Mesons");      } // end if (Pt_snk)      // Convolute the quark propagator with the sink smearing function.      // Make a copy of the quark propagator and then overwrite it with      // the convolution.       if (input.param.Sl_snk)       {	LatticePropagator quark_prop_smr;	quark_prop_smr = quark_propagator;	sink_smear2(u, quark_prop_smr, 		    input.param.wvf_kind, 		    input.param.wvf_param[loop],		    input.param.wvfIntPar[loop], 		    j_decay);	if (Pt_src)	  mesons(quark_prop_smr, quark_prop_smr, phases, t0,		 xml_array, "Point_Shell_Wilson_Mesons");	if (Sl_src)	  mesons(quark_prop_smr, quark_prop_smr, phases, t0,		 xml_array, "Shell_Shell_Wilson_Mesons");	if (Wl_src)	  mesons(quark_prop_smr, quark_prop_smr, phases_nomom, t0,		 xml_array, "Wall_Shell_Wilson_Mesons");      } // end if (Sl_snk)      // Wall sink      if (input.param.Wl_snk)       {	LatticePropagator wall_quark_prop;	wall_qprop(wall_quark_prop, quark_propagator, phases_nomom);	if (Pt_src)	  mesons(wall_quark_prop, wall_quark_prop, phases_nomom, t0,		 xml_array, "Point_Wall_Wilson_Mesons");	if (Sl_src)	  mesons(wall_quark_prop, wall_quark_prop, phases_nomom, t0,		 xml_array, "Shell_Wall_Wilson_Mesons");	if (Wl_src)	  mesons(wall_quark_prop, wall_quark_prop, phases_nomom, t0,		 xml_array, "Wall_Wall_Wilson_Mesons");      } // end if (Wl_snk)    } // end if (MesonP)    // Do the currents next    if (input.param.CurrentP)     {      // Construct the rho vector-current and the pion axial current divergence      if (Pt_src)	curcor2(u, quark_propagator, quark_propagator, phases, 		t0, 3,		xml_array, "Point_Point_Meson_Currents");              if (Sl_src)	curcor2(u, quark_propagator, quark_propagator, phases, 		t0, 3,		xml_array, "Shell_Point_Meson_Currents");              if (Wl_src)	curcor2(u, quark_propagator, quark_propagator, phases_nomom, 		t0, 3,		xml_array, "Wall_Point_Meson_Currents");    } // end if (CurrentP)    // Do the baryons    if (input.param.BaryonP)     {      // Construct {Point|Shell}-Point mesons, if desired      if (input.param.Pt_snk)       {	if (Pt_src)	  baryon(quark_propagator, phases, 		 t0, bc_spec, input.param.time_rev, 		 xml_array, "Point_Point_Wilson_Baryons");        	if (Sl_src)	  baryon(quark_propagator, phases, 		 t0, bc_spec, input.param.time_rev, 		 xml_array, "Shell_Point_Wilson_Baryons");        	if (Wl_src)	  baryon(quark_propagator, phases_nomom, 		 t0, bc_spec, input.param.time_rev, 		 xml_array, "Wall_Point_Wilson_Baryons");      } // end if (Pt_snk)      // Convolute the quark propagator with the sink smearing function.      // Make a copy of the quark propagator and then overwrite it with      // the convolution.       if (input.param.Sl_snk)       {	LatticePropagator quark_prop_smr;	quark_prop_smr = quark_propagator;	sink_smear2(u, quark_prop_smr, 		    input.param.wvf_kind, 		    input.param.wvf_param[loop],		    input.param.wvfIntPar[loop], 		    j_decay);	if (Pt_src)	  baryon(quark_prop_smr, phases, 		 t0, bc_spec, input.param.time_rev, 		 xml_array, "Point_Shell_Wilson_Baryons");        	if (Sl_src)	  baryon(quark_prop_smr, phases, 		 t0, bc_spec, input.param.time_rev, 			 xml_array, "Shell_Shell_Wilson_Baryons");        	if (Wl_src)	  baryon(quark_prop_smr, phases_nomom, 		 t0, bc_spec, input.param.time_rev, 			 xml_array, "Wall_Shell_Wilson_Baryons");      } // end if (Sl_snk)      // Wall sink      if (input.param.Wl_snk)       {	LatticePropagator wall_quark_prop;	wall_qprop(wall_quark_prop, quark_propagator, phases_nomom);	if (Pt_src)	  baryon(wall_quark_prop, phases_nomom, 		 t0, bc_spec, input.param.time_rev, 		 xml_array, "Point_Wall_Wilson_Baryons");        	if (Sl_src)	  baryon(wall_quark_prop, phases_nomom, 		 t0, bc_spec, input.param.time_rev, 			 xml_array, "Shell_Wall_Wilson_Baryons");        	if (Wl_src)	  baryon(wall_quark_prop, phases_nomom, 		 t0, bc_spec, input.param.time_rev, 			 xml_array, "Wall_Wall_Wilson_Baryons");      } // end if (Wl_snk)    } // end if (BaryonP)    pop(xml_array);  // array element  } // end for(loop)  pop(xml_array);  // Wilson_spectroscopy  pop(xml_out);  // spectrum_w  xml_out.close();  xml_in.close();  // Time to bolt  QDP_finalize();  exit(0);}