CHROMA
inline_heavyhadspec_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Inline construction of heavy hadron spectrum in SU(3)
3  *
4  */
5 
9 #include "meas/glue/mesplq.h"
10 #include "util/ft/sftmom.h"
11 #include "util/info/proginfo.h"
12 #include "io/param_io.h"
13 #include "io/qprop_io.h"
17 
18 namespace Chroma
19 {
20  namespace InlineHeavyHadSpecEnv
21  {
22  namespace
23  {
24  AbsInlineMeasurement* createMeasurement(XMLReader& xml_in,
25  const std::string& path)
26  {
27  return new InlineHeavyHadSpec(InlineHeavyHadSpecParams(xml_in, path));
28  }
29 
30  //! Local registration flag
31  bool registered = false;
32  }
33 
34  const std::string name = "HEAVY_HADRON_SPECTRUM";
35 
36  //! Register all the factories
37  bool registerAll()
38  {
39  bool success = true;
40  if (! registered)
41  {
42  success &= TheInlineMeasurementFactory::Instance().registerObject(name, createMeasurement);
43  registered = true;
44  }
45  return success;
46  }
47  }
48 
49 
50 
51  //! Reader for parameters
52  void read(XMLReader& xml, const std::string& path, InlineHeavyHadSpecParams::Param_t& param)
53  {
54  XMLReader paramtop(xml, path);
55 
56  int version;
57  read(paramtop, "version", version);
58 
59  switch (version)
60  {
61  case 1:
62  break;
63 
64  default:
65  QDPIO::cerr << "Input parameter version " << version << " unsupported." << std::endl;
66  QDP_abort(1);
67  }
68 
69  read(paramtop, "time_rev", param.time_rev);
70 
71  read(paramtop, "mom2_max", param.mom2_max);
72  read(paramtop, "avg_equiv_mom", param.avg_equiv_mom);
73  }
74 
75 
76  //! Writer for parameters
77  void write(XMLWriter& xml, const std::string& path, const InlineHeavyHadSpecParams::Param_t& param)
78  {
79  push(xml, path);
80 
81  int version = 1;
82  write(xml, "version", version);
83 
84  write(xml, "time_rev", param.time_rev);
85 
86  write(xml, "mom2_max", param.mom2_max);
87  write(xml, "avg_equiv_mom", param.avg_equiv_mom);
88 
89  pop(xml);
90  }
91 
92 
93  //! Propagator input
94  void read(XMLReader& xml, const std::string& path, InlineHeavyHadSpecParams::NamedObject_t::Props_t& input)
95  {
96  XMLReader inputtop(xml, path);
97 
98  read(inputtop, "first_id", input.first_id);
99  read(inputtop, "second_id", input.second_id);
100  }
101 
102  //! Propagator output
103  void write(XMLWriter& xml, const std::string& path, const InlineHeavyHadSpecParams::NamedObject_t::Props_t& input)
104  {
105  push(xml, path);
106 
107  write(xml, "first_id", input.first_id);
108  write(xml, "second_id", input.second_id);
109 
110  pop(xml);
111  }
112 
113 
114  //! Propagator input
115  void read(XMLReader& xml, const std::string& path, InlineHeavyHadSpecParams::NamedObject_t& input)
116  {
117  XMLReader inputtop(xml, path);
118 
119  read(inputtop, "gauge_id", input.gauge_id);
120  read(inputtop, "sink_pairs", input.sink_pairs);
121  }
122 
123  //! Propagator output
124  void write(XMLWriter& xml, const std::string& path, const InlineHeavyHadSpecParams::NamedObject_t& input)
125  {
126  push(xml, path);
127 
128  write(xml, "gauge_id", input.gauge_id);
129  write(xml, "sink_pairs", input.sink_pairs);
130 
131  pop(xml);
132  }
133 
134 
135  // Param stuff
137  {
138  frequency = 0;
139  }
140 
142  {
143  try
144  {
145  XMLReader paramtop(xml_in, path);
146 
147  if (paramtop.count("Frequency") == 1)
148  read(paramtop, "Frequency", frequency);
149  else
150  frequency = 1;
151 
152  // Parameters for source construction
153  read(paramtop, "Param", param);
154 
155  // Read in the output propagator/source configuration info
156  read(paramtop, "NamedObject", named_obj);
157 
158  // Possible alternate XML file pattern
159  if (paramtop.count("xml_file") != 0)
160  {
161  read(paramtop, "xml_file", xml_file);
162  }
163  }
164  catch(const std::string& e)
165  {
166  QDPIO::cerr << __func__ << ": Caught Exception reading XML: " << e << std::endl;
167  QDP_abort(1);
168  }
169  }
170 
171 
172  void
173  InlineHeavyHadSpecParams::write(XMLWriter& xml_out, const std::string& path)
174  {
175  push(xml_out, path);
176 
177  Chroma::write(xml_out, "Param", param);
178  Chroma::write(xml_out, "NamedObject", named_obj);
179  QDP::write(xml_out, "xml_file", xml_file);
180 
181  pop(xml_out);
182  }
183 
184 
185 
186  // Anonymous namespace
187  namespace
188  {
189  //! Useful structure holding sink props
190  struct SinkPropContainer_t
191  {
192  ForwardProp_t prop_header;
194  Real Mass;
195 
196  multi1d<int> bc;
197 
202  };
203 
204 
205  //! Useful structure holding sink props
206  struct AllSinkProps_t
207  {
208  SinkPropContainer_t sink_prop_1;
209  SinkPropContainer_t sink_prop_2;
210  };
211 
212 
213  //! Read a sink prop
214  void readSinkProp(SinkPropContainer_t& s, const std::string& id)
215  {
216  try
217  {
218  // Try a cast to see if it succeeds
219  const LatticePropagator& foo =
220  TheNamedObjMap::Instance().getData<LatticePropagator>(id);
221 
222  // Snarf the data into a copy
223  s.quark_propagator_id = id;
224 
225  // Snarf the prop info. This is will throw if the prop_id is not there
226  XMLReader prop_file_xml, prop_record_xml;
227  TheNamedObjMap::Instance().get(id).getFileXML(prop_file_xml);
228  TheNamedObjMap::Instance().get(id).getRecordXML(prop_record_xml);
229 
230  // Try to invert this record XML into a ChromaProp struct
231  // Also pull out the id of this source
232  {
233  std::string xpath;
234  read(prop_record_xml, "/SinkSmear", s.prop_header);
235 
236  read(prop_record_xml, "/SinkSmear/PropSource/Source/SourceType", s.source_type);
237  xpath = "/SinkSmear/PropSource/Source/Displacement/DisplacementType";
238  if (prop_record_xml.count(xpath) != 0)
239  read(prop_record_xml, xpath, s.source_disp_type);
240  else
241  s.source_disp_type = NoQuarkDisplacementEnv::getName();
242 
243  read(prop_record_xml, "/SinkSmear/PropSink/Sink/SinkType", s.sink_type);
244  xpath = "/SinkSmear/PropSink/Sink/Displacement/DisplacementType";
245  if (prop_record_xml.count(xpath) != 0)
246  read(prop_record_xml, xpath, s.sink_disp_type);
247  else
248  s.sink_disp_type = NoQuarkDisplacementEnv::getName();
249  }
250  }
251  catch( std::bad_cast )
252  {
253  QDPIO::cerr << InlineHeavyHadSpecEnv::name << ": caught dynamic cast error"
254  << std::endl;
255  QDP_abort(1);
256  }
257  catch (const std::string& e)
258  {
259  QDPIO::cerr << InlineHeavyHadSpecEnv::name << ": error message: " << e
260  << std::endl;
261  QDP_abort(1);
262  }
263 
264 
265  // Derived from input prop
266  // Hunt around to find the mass
267  // NOTE: this may be problematic in the future if actions are used with no
268  // clear def. of a Mass
269  QDPIO::cout << "Try action and mass" << std::endl;
270  s.Mass = getMass(s.prop_header.prop_header.fermact);
271 
272  // Only baryons care about boundaries
273  // Try to find them. If not present, assume dirichlet.
274  // This turns off any attempt to time reverse which is the
275  // only thing that the BC are affecting.
276  s.bc.resize(Nd);
277  s.bc = 0;
278 
279  try
280  {
281  s.bc = getFermActBoundary(s.prop_header.prop_header.fermact);
282  }
283  catch (const std::string& e)
284  {
285  QDPIO::cerr << InlineHeavyHadSpecEnv::name
286  << ": caught exception. No BC found in these headers. Will assume dirichlet: " << e
287  << std::endl;
288  }
289 
290  QDPIO::cout << "FermAct = " << s.prop_header.prop_header.fermact.id << std::endl;
291  QDPIO::cout << "Mass = " << s.Mass << std::endl;
292  }
293 
294 
295  //! Read all sinks
296  void readAllSinks(AllSinkProps_t& s,
297  InlineHeavyHadSpecParams::NamedObject_t::Props_t sink_pair)
298  {
299  QDPIO::cout << "Attempt to parse forward propagator = " << sink_pair.first_id << std::endl;
300  readSinkProp(s.sink_prop_1, sink_pair.first_id);
301  QDPIO::cout << "Forward propagator successfully parsed" << std::endl;
302 
303  QDPIO::cout << "Attempt to parse forward propagator = " << sink_pair.second_id << std::endl;
304  readSinkProp(s.sink_prop_2, sink_pair.second_id);
305  QDPIO::cout << "Forward propagator successfully parsed" << std::endl;
306  }
307 
308  } // namespace anonymous
309 
310 
311 
312  // Function call
313  void
314  InlineHeavyHadSpec::operator()(unsigned long update_no,
315  XMLWriter& xml_out)
316  {
317  // If xml file not empty, then use alternate
318  if (params.xml_file != "")
319  {
320  std::string xml_file = makeXMLFileName(params.xml_file, update_no);
321 
322  push(xml_out, "heavyhadspec");
323  write(xml_out, "update_no", update_no);
324  write(xml_out, "xml_file", xml_file);
325  pop(xml_out);
326 
327  XMLFileWriter xml(xml_file);
328  func(update_no, xml);
329  }
330  else
331  {
332  func(update_no, xml_out);
333  }
334  }
335 
336 
337  // Real work done here
338  void
339  InlineHeavyHadSpec::func(unsigned long update_no,
340  XMLWriter& xml_out)
341  {
342  START_CODE();
343 
344  StopWatch snoop;
345  snoop.reset();
346  snoop.start();
347 
348  // Test and grab a reference to the gauge field
349  XMLBufferWriter gauge_xml;
350  try
351  {
352  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
353  TheNamedObjMap::Instance().get(params.named_obj.gauge_id).getRecordXML(gauge_xml);
354  }
355  catch( std::bad_cast )
356  {
357  QDPIO::cerr << InlineHeavyHadSpecEnv::name << ": caught dynamic cast error"
358  << std::endl;
359  QDP_abort(1);
360  }
361  catch (const std::string& e)
362  {
363  QDPIO::cerr << InlineHeavyHadSpecEnv::name << ": std::map call failed: " << e
364  << std::endl;
365  QDP_abort(1);
366  }
367  const multi1d<LatticeColorMatrix>& u =
368  TheNamedObjMap::Instance().getData< multi1d<LatticeColorMatrix> >(params.named_obj.gauge_id);
369 
370  push(xml_out, "heavyhadspec");
371  write(xml_out, "update_no", update_no);
372 
373  QDPIO::cout << " HEAVYHADSPEC: SU(3) heavy hadron spectroscopy for Wilson-like fermions" << std::endl;
374  QDPIO::cout << std::endl << " Gauge group: SU(" << Nc << ")" << std::endl;
375  QDPIO::cout << " volume: " << Layout::lattSize()[0];
376  for (int i=1; i<Nd; ++i) {
377  QDPIO::cout << " x " << Layout::lattSize()[i];
378  }
379  QDPIO::cout << std::endl;
380 
381  proginfo(xml_out); // Print out basic program info
382 
383  // Write out the input
384  params.write(xml_out, "Input");
385 
386  // Write out the config info
387  write(xml_out, "Config_info", gauge_xml);
388 
389  push(xml_out, "Output_version");
390  write(xml_out, "out_version", 15);
391  pop(xml_out);
392 
393 
394  // First calculate some gauge invariant observables just for info.
395  MesPlq(xml_out, "Observables", u);
396 
397  // Keep an array of all the xml output buffers
398  push(xml_out, "Wilson_hadron_measurements");
399 
400  // Now loop over the various fermion pairs
401  for(int lpair=0; lpair < params.named_obj.sink_pairs.size(); ++lpair)
402  {
404 
405  push(xml_out, "elem");
406 
407  AllSinkProps_t all_sinks;
408  readAllSinks(all_sinks, named_obj);
409 
410  // Derived from input prop
411  multi1d<int> t_srce
412  = all_sinks.sink_prop_1.prop_header.source_header.getTSrce();
413  int j_decay = all_sinks.sink_prop_1.prop_header.source_header.j_decay;
414  int t0 = all_sinks.sink_prop_1.prop_header.source_header.t_source;
415 
416  // Sanity checks
417  {
418  if (all_sinks.sink_prop_2.prop_header.source_header.j_decay != j_decay)
419  {
420  QDPIO::cerr << "Error!! j_decay must be the same for all propagators " << std::endl;
421  QDP_abort(1);
422  }
423  if (all_sinks.sink_prop_2.prop_header.source_header.t_source !=
424  all_sinks.sink_prop_1.prop_header.source_header.t_source)
425  {
426  QDPIO::cerr << "Error!! t_source must be the same for all propagators " << std::endl;
427  QDP_abort(1);
428  }
429  if (all_sinks.sink_prop_1.source_type != all_sinks.sink_prop_2.source_type)
430  {
431  QDPIO::cerr << "Error!! source_type must be the same in a pair " << std::endl;
432  QDP_abort(1);
433  }
434  if (all_sinks.sink_prop_1.sink_type != all_sinks.sink_prop_2.sink_type)
435  {
436  QDPIO::cerr << "Error!! source_type must be the same in a pair " << std::endl;
437  QDP_abort(1);
438  }
439  }
440 
441  // Only baryons care about bc
442  int bc_spec = 0;
443  bc_spec = all_sinks.sink_prop_1.bc[j_decay] ;
444  if (all_sinks.sink_prop_2.bc[j_decay] != bc_spec)
445  {
446  QDPIO::cerr << "Error!! bc must be the same for all propagators " << std::endl;
447  QDP_abort(1);
448  }
449 
450 
451 
452  // Initialize the slow Fourier transform phases
454  j_decay);
455 
456  // Keep a copy of the phases with NO momenta
457  SftMom phases_nomom(0, true, j_decay);
458 
459  // Masses
460  write(xml_out, "Mass_1", all_sinks.sink_prop_1.Mass);
461  write(xml_out, "Mass_2", all_sinks.sink_prop_2.Mass);
462  write(xml_out, "t0", t0);
463 
464  // Save prop input
465  push(xml_out, "Forward_prop_headers");
466  write(xml_out, "First_forward_prop", all_sinks.sink_prop_1.prop_header);
467  write(xml_out, "Second_forward_prop", all_sinks.sink_prop_2.prop_header);
468  pop(xml_out);
469 
470  // Sanity check - write out the norm2 of the forward prop in the j_decay direction
471  // Use this for any possible verification
472  push(xml_out, "Forward_prop_correlator");
473  {
474  const LatticePropagator& sink_prop_1 =
475  TheNamedObjMap::Instance().getData<LatticePropagator>(all_sinks.sink_prop_1.quark_propagator_id);
476  const LatticePropagator& sink_prop_2 =
477  TheNamedObjMap::Instance().getData<LatticePropagator>(all_sinks.sink_prop_2.quark_propagator_id);
478 
479  write(xml_out, "forward_prop_corr_1", sumMulti(localNorm2(sink_prop_1), phases.getSet()));
480  write(xml_out, "forward_prop_corr_2", sumMulti(localNorm2(sink_prop_2), phases.getSet()));
481  }
482  pop(xml_out);
483 
484 
485  push(xml_out, "SourceSinkType");
486  {
487  QDPIO::cout << "Source_type_1 = " << all_sinks.sink_prop_1.source_type << std::endl;
488  QDPIO::cout << "Sink_type_1 = " << all_sinks.sink_prop_1.sink_type << std::endl;
489  QDPIO::cout << "Source_type_2 = " << all_sinks.sink_prop_2.source_type << std::endl;
490  QDPIO::cout << "Sink_type_2 = " << all_sinks.sink_prop_2.sink_type << std::endl;
491 
492  write(xml_out, "source_type_1", all_sinks.sink_prop_1.source_type);
493  write(xml_out, "source_disp_type_1", all_sinks.sink_prop_1.source_disp_type);
494  write(xml_out, "sink_type_1", all_sinks.sink_prop_1.sink_type);
495  write(xml_out, "sink_disp_type_1", all_sinks.sink_prop_1.sink_disp_type);
496 
497  write(xml_out, "source_type_2", all_sinks.sink_prop_2.source_type);
498  write(xml_out, "source_disp_type_2", all_sinks.sink_prop_2.source_disp_type);
499  write(xml_out, "sink_type_2", all_sinks.sink_prop_2.sink_type);
500  write(xml_out, "sink_disp_type_2", all_sinks.sink_prop_2.sink_disp_type);
501  }
502  pop(xml_out);
503 
504 
505  // References for use later
506  const LatticePropagator& sink_prop_1 =
507  TheNamedObjMap::Instance().getData<LatticePropagator>(all_sinks.sink_prop_1.quark_propagator_id);
508  const LatticePropagator& sink_prop_2 =
509  TheNamedObjMap::Instance().getData<LatticePropagator>(all_sinks.sink_prop_2.quark_propagator_id);
510 
511 
512  // Construct group name for output
513  std::string src_type;
514  if (all_sinks.sink_prop_1.source_type == "POINT_SOURCE")
515  src_type = "Point";
516  else if (all_sinks.sink_prop_1.source_type == "SF_POINT_SOURCE")
517  src_type = "Point";
518  else if (all_sinks.sink_prop_1.source_type == "SHELL_SOURCE")
519  src_type = "Shell";
520  else if (all_sinks.sink_prop_1.source_type == "SF_SHELL_SOURCE")
521  src_type = "Shell";
522  else if (all_sinks.sink_prop_1.source_type == "WALL_SOURCE")
523  src_type = "Wall";
524  else if (all_sinks.sink_prop_1.source_type == "SF_WALL_SOURCE")
525  src_type = "Wall";
526  else
527  {
528  QDPIO::cerr << "Unsupported source type = " << all_sinks.sink_prop_1.source_type << std::endl;
529  QDP_abort(1);
530  }
531 
532  std::string snk_type;
533  if (all_sinks.sink_prop_1.sink_type == "POINT_SINK")
534  snk_type = "Point";
535  else if (all_sinks.sink_prop_1.sink_type == "SHELL_SINK")
536  snk_type = "Shell";
537  else if (all_sinks.sink_prop_1.sink_type == "WALL_SINK")
538  snk_type = "Wall";
539  else
540  {
541  QDPIO::cerr << "Unsupported sink type = " << all_sinks.sink_prop_1.sink_type << std::endl;
542  QDP_abort(1);
543  }
544 
545  std::string source_sink_type = src_type + "_" + snk_type;
546  QDPIO::cout << "Source type = " << src_type << std::endl;
547  QDPIO::cout << "Sink type = " << snk_type << std::endl;
548 
549 
551  phases_nomom,xml_out,
552  source_sink_type + "_Heavy_Hadrons_SU3");
553 
554  pop(xml_out); // array element
555  }
556  pop(xml_out); // Wilson_spectroscopy
557  pop(xml_out); // heavyhadspec
558 
559  snoop.stop();
560  QDPIO::cout << InlineHeavyHadSpecEnv::name << ": total time = "
561  << snoop.getTimeInSeconds()
562  << " secs" << std::endl;
563 
564  QDPIO::cout << InlineHeavyHadSpecEnv::name << ": ran successfully" << std::endl;
565 
566  END_CODE();
567  }
568 
569 }
Inline measurement factory.
Inline measurement of hadron spectrum.
void operator()(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
InlineHeavyHadSpecParams params
void func(const unsigned long update_no, XMLWriter &xml_out)
Do the measurement.
Fourier transform phase factor support.
Definition: sftmom.h:35
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void static_light_su3(const multi1d< LatticeColorMatrix > &u, const LatticePropagator &quark1, const LatticePropagator &quark2, const multi1d< int > &src, const SftMom &phases, XMLWriter &xml, const std::string &xml_group)
Heavy hadron spectrum for SU(3) isospin limit.
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
std::string makeXMLFileName(std::string xml_file, unsigned long update_no)
Return a xml file name for inline measurements.
Potential between 2 heavy hadrons : Detmold Correlators checked independentely by Savage.
int bc_spec
multi1d< int > t_srce
std::string source_disp_type
multi1d< int > bc
std::string source_type
std::string quark_propagator_id
SinkPropContainer_t sink_prop_2
std::string sink_type
SinkPropContainer_t sink_prop_1
ForwardProp_t prop_header
std::string sink_disp_type
Inline hadron spectrum calculations.
Make xml file writer.
int j_decay
Definition: meslate.cc:22
Nd
Definition: meslate.cc:74
Named object function std::map.
static bool registered
Local registration flag.
bool registerAll()
Register all the factories.
std::string getName()
Return the name.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Real getMass(const GroupXML_t &fermact)
Given a fermion action in std::string form, return the Mass.
Definition: qprop_io.cc:16
static multi1d< LatticeColorMatrix > u
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > & Mass
Definition: pbg5p_w.cc:29
push(xml_out,"Condensates")
int i
Definition: pbg5p_w.cc:55
multi1d< int > getFermActBoundary(const GroupXML_t &fermact)
Given a fermion action in std::string form, return the boundary.
Definition: qprop_io.cc:61
void write(XMLWriter &xml, const std::string &path, const InlineHeavyHadSpecParams::NamedObject_t &input)
Propagator output.
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
START_CODE()
multi1d< LatticeFermion > s(Ncb)
::std::string string
Definition: gtest.h:1979
No quark displacement.
Various parameter structs and reader/writers.
Print out basic info about this program.
Routines associated with Chroma propagator IO.
Fourier transform phase factor support.
void write(XMLWriter &xml_out, const std::string &path)
struct Chroma::InlineHeavyHadSpecParams::Param_t param
struct Chroma::InlineHeavyHadSpecParams::NamedObject_t named_obj