ROOTPWA
TFhh.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <string>
4 #include "TFhh.h"
5 using namespace std;
6 
7 extern ClebschGordanBox box;
8 
10 
12  Int_t lambda_, Int_t nu_,
13  Int_t nLS, TLSAmpl* *LSampl,
14  Int_t even_contr_) {
15 
16  J=J_;
17  lambda= lambda_;
18  nu= nu_;
19  even_contraction = even_contr_;
20  Int_t delta=lambda-nu;
21  Nterms=0;
22 
23  name_str=new char[15];
24  sprintf(name_str, "F_%lld_%lld", lambda, nu);
25 
26  for (Int_t iLS=0; iLS<nLS; iLS++) {
27  if (LSampl[iLS]->Getdelta() == delta ||
28  LSampl[iLS]->Getdelta() == -delta ) {
29  // cout << "iLS=" << iLS << ", delta=" << delta << endl;
30  TFracNum *CG3S = box.GetCG(LSampl[iLS]->GetS(), S1, S2);
31  TFracNum SpinCouplFac = CG3S[CGIndex(S1, lambda, S2, -nu)];
32  if (SpinCouplFac==TFracNum_Zero) {
33  if (debugFhh==2)
34  cout << "Clebsch-Gordan is zero" << endl;
35  }
36  else {
37  TLSContrib* *newPTR = new TLSContrib*[Nterms+1];
38  for (Int_t iC=0; iC<Nterms; iC++) newPTR[iC] = LSt[iC];
39  newPTR[Nterms] = new TLSContrib(LSampl[iLS], delta, SpinCouplFac);
40  // delete[] LSt;
41  Nterms++;
42  LSt=newPTR;
43  }
44  }
45  }
46  if (debugFhh) Print();
47 }
48 
49 TFhh::TFhh(TFhh *sFhh, char flag) {
50  if (flag != 'i' && flag !='m'){
51  cerr << "TFhh::TFhh unknown flag "<<flag<<endl;
52  return;
53  }
54  if (debugFhh)
55  cout << "Initializing from single Amplitude" << endl;
56  if ( ( flag =='i' && (sFhh->GetJ()) % 2 ) ||
57  ( flag =='m' && (sFhh->GetJ()) % 2 == 0 ) ) {
58  cout << sFhh->GetName() << "[symm] = 0" << endl;
59  Nterms=0;
60  }
61  else {
62  name_str = new char[15];
63  sprintf(name_str, "%s[symm]", sFhh->GetName());
64  J = sFhh->GetJ();
65  lambda = sFhh->GetLambda();
66  nu = sFhh->GetNu();
67  even_contraction = sFhh->GetEvenContr();
68  Nterms = sFhh->GetNterms();
69  LSt = sFhh->GetLStPtr();
70  Print();
71  }
72 }
73 
74 TFhh::TFhh(TFhh *sFhh, TFhh *xFhh) {
75  name_str = new char[15];
76  sprintf(name_str, "%s[symm]", sFhh->GetName());
77  J = sFhh->GetJ();
78  lambda = sFhh->GetLambda();
79  nu = sFhh->GetNu();
80  even_contraction = sFhh->GetEvenContr();
81  Nterms=0;
82  LSt =0;
83  if (J != xFhh->GetJ() || even_contraction != xFhh->GetEvenContr() ||
84  lambda != xFhh->GetNu() || nu != xFhh->GetLambda() ){
85  cerr << "TFhh::TFhh(TFhh *, TFhh*): Something is wrong," << endl
86  << " source amplitudes for symmetrization do not match" << endl;
87  return;
88  }
89 
90  // Since some LS-terms cancel out, first a max-length array of
91  // LSContrib pointers is filled, and then squeezed to only
92  // non-zero contributions
93 
94  Int_t Ns = sFhh->GetNterms();
95  Int_t Nx = xFhh->GetNterms();
96  Nterms = Ns + Nx;
97  TLSContrib* *pLSt = new TLSContrib*[Nterms];
98 
99  Int_t prevterms=0;
100 
101  for (Int_t i=0; i<Nterms; i++) {
102  TLSContrib *to_be_added;
103  if (i<Ns) to_be_added = sFhh->GetLStPtr()[i];
104  else to_be_added = xFhh->GetLStPtr()[i-Ns];
105  Int_t found_in_prevterm=0;
106  for (Int_t j=0; j<prevterms; j++) {
107  if ( pLSt[j]->SameParameter(to_be_added) ) {
108  found_in_prevterm=1;
109  if (i<Ns) pLSt[j]->Add(to_be_added, false);
110  else pLSt[j]->Add(to_be_added, true);
111  }
112  }
113  if (!found_in_prevterm) {
114  if (i<Ns) pLSt[prevterms] = new TLSContrib(to_be_added, false);
115  else pLSt[prevterms] = new TLSContrib(to_be_added, true);
116  prevterms++;
117  }
118  }
119 
120  //
121  // Cancel zeros
122  //
123  Int_t non_zeros=0;
124  for (Int_t i=0; i<prevterms; i++) {
125  if ( pLSt[i]->GetNterms() != 0 ) non_zeros++;
126  }
127 
128  Nterms=non_zeros;
129 
130  if (Nterms) {
131  LSt = new TLSContrib*[Nterms];
132  Int_t j=0;
133  for (Int_t i=0; i<prevterms; i++) {
134  if ( pLSt[i]->GetNterms() != 0 ) {
135  LSt[j]=new TLSContrib(pLSt[i], false);
136  j++;
137  }
138  }
139  }
140  Print();
141 
142 }
143 
144 
146  NNRterms=0;
147  for (Int_t i=0; i<Nterms; i++){
148  if (! LSt[i]->IsPureRelativistic()) {
149  Int_t jfound=-1;
150  for (Int_t j=0; j<NNRterms; j++) {
151  if (NRLSt[j]->CheckJLS(LSt[i])){
152  jfound=j;
153  break;
154  }
155  }
156  if (jfound != -1) {
157  NRLSt[jfound]->Add(LSt[i]);
158  }
159  else {
160  NNRterms++;
161  TLSNonRel **newNRLS=new TLSNonRel*[NNRterms];
162  for (Int_t j=0; j<NNRterms-1; j++) {
163  newNRLS[j]=NRLSt[j];
164  }
165  newNRLS[NNRterms-1] = new TLSNonRel(LSt[i]);
166  NRLSt=newNRLS;
167  }
168  }
169  }
170  cout << name_str << " (NR) = " << endl;
171  for (Int_t j=0; j<NNRterms; j++) {
172  cout << "LS=" << NRLSt[j]->GetL() << NRLSt[j]->GetS() << ": ";
173  NRLSt[j]->Print();
174  }
175  return 0;
176 }
177 
179  cout << name_str << " (NR) = " << endl;
180  for (Int_t j=0; j<NNRterms; j++) {
181  cout << "LS=" << NRLSt[j]->GetL() << NRLSt[j]->GetS() << " => ";
182  NRLSt[j]->PrintG();
183  }
184  return 0;
185 }
186 
187 
188 Int_t
190  cout << name_str << " =";
191  if (even_contraction) cout << endl;
192  else cout << " (iw)" << endl;
193  for (Int_t iLSt=0; iLSt<Nterms; iLSt++) {
194  LSt[iLSt]->Print();
195  }
196  if (Nterms==0) cout << " 0" << endl;
197  return 0;
198 }
199 
201 
203  if (debugJSS>=2) {
204  cout << " Decay channel: " << JMother;
205  if (etaJ>0) cout << "+ -> ";
206  else cout << "- -> ";
207  cout << SDecay1;
208  if (eta1>0) cout << "+ ";
209  else cout << "- ";
210  cout << SDecay2;
211  if (eta2>0) cout << "+";
212  else cout << "-";
213  cout << endl;
214  }
215 
216  // range for coupled Spin
217  Smin=SDecay1-SDecay2; if (Smin<0) Smin=-Smin;
218  Smax=SDecay1+SDecay2;
219 
220  if (debugJSS>=2) {
221  cout << "possible S:";
222  Int_t iS=Smin;
223  while (iS<=Smax) {
224  cout << " " << iS;
225  iS++;
226  }
227  }
228  cout << endl;
229 
230  Int_t intr_parity = etaJ*eta1*eta2;
231 
232  Lmax=JMother+Smax;
233  Lmin=Lmax;
234  for (Int_t iS=Smin; iS<=Smax; iS++) {
235  Int_t Lm1=JMother-iS; if (Lm1<0) Lm1=-Lm1;
236  if (Lm1<Lmin) Lmin=Lm1;
237  }
238 
239  if (debugJSS>=2) cout << "possible L:";
240  Int_t NL=0;
241  Int_t *fL=0;
242  Int_t testL=0; // L even
243  if (intr_parity<0) testL=1; // L odd
244  Int_t tL=testL;
245 
246  while(testL<=Lmax) {
247  if (testL>=Lmin) NL++;
248  testL+=2;
249  }
250 
251  if (NL) {
252  fL=new Int_t[NL];
253  NL=0;
254  while(tL<=Lmax) {
255  if (tL>=Lmin) {
256  fL[NL]=tL;
257  if (debugJSS>=2) cout << " " << fL[NL];
258  NL++;
259  }
260  tL+=2;
261  }
262  }
263  if (debugJSS>=2) cout << endl;
264 
265  Int_t even_contraction=1;
266  if ((SDecay1+SDecay2+testL-JMother) % 2) even_contraction=0;
267 
268  if (debugJSS>=2) {
269  if (even_contraction) cout << "contraction only with g~" << endl;
270  else cout << "contraction with eps*p and g~" << endl;
271  }
272 
273  Int_t MaxPsiInternal=SDecay1;
274  if (SDecay2<SDecay1) MaxPsiInternal=SDecay2;
275 
276  Int_t MaxPsiPhi = SDecay1+SDecay2;
277  if (JMother<MaxPsiPhi) MaxPsiPhi=JMother;
278 
279  Int_t preloop=1;
280 
281  NLSAmpl=0;
282  LSAmplitudes = 0;
283 
284  while (preloop==1 || preloop==0) {
285 
286  if (preloop==0) {
287  LSAmplitudes = new TLSAmpl*[NLSAmpl];
288  NLSAmpl=0;
289  if (debugJSS>=2) cout << endl << "*************" << endl;
290  }
291 
292  for (Int_t il=0; il<NL; il++) {
293 
294  Int_t L = fL[il];
295 
296  Int_t MaxPsiChi = SDecay1+SDecay2;
297  if (L<MaxPsiChi) MaxPsiChi=L;
298 
299  Int_t MaxChiPhi = JMother;
300  if (L<JMother) MaxChiPhi=L;
301 
302  // possible S in this case:
303  Int_t SminL=JMother-L;
304  if (SminL<0) SminL=-SminL;
305  if (SminL<Smin) SminL=Smin;
306  Int_t SmaxL=JMother+L;
307  if (SmaxL>Smax) SmaxL=Smax;
308  for (Int_t S_L=SminL; S_L<=SmaxL; S_L++) {
309  if (debugJSS>=2) cout << "Amplitudes for L=" << L << " S="<< S_L
310  << " Rank scheme [ "
311  << SDecay1+SDecay2 << " "
312  << L << " "
313  << JMother << "]" << endl;
314  Int_t totalRank=SDecay1+SDecay2+L+JMother;
315  Int_t MaxDelta=S_L;
316  if (JMother<S_L) MaxDelta=JMother;
317 
318  Int_t IndexContractions=(totalRank-3)/2;
319  if (even_contraction) IndexContractions= totalRank / 2;
320  if (debugJSS>=2) cout << IndexContractions
321  << " Lorentz contractions." << endl;
322 
323  Int_t MaxContractionNumber=0;
324 
325  for (Int_t PsiInternal=0; PsiInternal<=MaxPsiInternal; PsiInternal++) {
326  for (Int_t cChiPhi=0; cChiPhi<=MaxChiPhi; cChiPhi++) {
327  for (Int_t cPsiChi=0; cPsiChi<=MaxPsiChi; cPsiChi++) {
328  for (Int_t cPsiPhi=0; cPsiPhi<=MaxPsiPhi; cPsiPhi++) {
329  if (debugJSS==3)
330  cout << "Checking " << PsiInternal << " " << cPsiChi << " "
331  << cChiPhi << " " << cPsiPhi; // << endl;
332  if (PsiInternal+cPsiChi+cChiPhi+cPsiPhi != IndexContractions){
333  if (debugJSS==3) cout << " C-" << endl;
334  continue;
335  }
336  else if (debugJSS==3) cout << " C+";
337  if (even_contraction) {
338  if (2*PsiInternal+cPsiChi+cPsiPhi != SDecay1+SDecay2) {
339  if (debugJSS==3) cout << "S-" << endl;
340  continue;
341  }
342  else if (debugJSS==3) cout << "S+";
343  if (cPsiChi+cChiPhi != L) {
344  if (debugJSS==3) cout << "L-" << endl;
345  continue;
346  }
347  else if (debugJSS==3) cout << "L+";
348  if (cChiPhi+cPsiPhi != JMother) {
349  if (debugJSS==3) cout << "J-" << endl;
350  continue;
351  }
352  else if (debugJSS==3) cout << "J+";
353  }
354  else {
355  if (L -cPsiChi-cChiPhi>1) {
356  if (debugJSS==3) cout << "L-" << endl;
357  continue;
358  }
359  else if (debugJSS==3) cout << "L+";
360  if (JMother-cPsiPhi-cChiPhi>1) {
361  if (debugJSS==3) cout << "J-" << endl;
362  continue;
363  }
364  else if (debugJSS==3) cout << "J+";
365  }
366  Int_t r_ome=SDecay1-PsiInternal;
367  Int_t r_eps=SDecay2-PsiInternal;
368 
369  if (!even_contraction) {
370  Int_t PsiRest=SDecay1+SDecay2-2*PsiInternal-cPsiChi-cPsiPhi;
371  if ( PsiRest<0 || PsiRest>2 ) {
372  if (debugJSS==3) cout << "R-" << endl;
373  continue;
374  }
375  else if (debugJSS==3) cout << "R+";
376  if (PsiRest==2){
377  if (r_ome>0) r_ome--;
378  else {
379  if (debugJSS==3) cout << "O-";
380  continue;
381  }
382  if (r_eps>0) r_eps--;
383  else {
384  if (debugJSS==3) cout << "E-";
385  continue;
386  }
387  }
388  }
389  //
390  // Original ordering:
391  //
392  //for (Int_t cPhiOmega=0; cPhiOmega<=r_ome; cPhiOmega++) {
393  // for (Int_t cChiOmega=0; cChiOmega<=r_ome-cPhiOmega;
394  // cChiOmega++) {
395  //
396  // For agreement with paper:
397  //
398  // error found 14.10.07 "r_eps" replaced with "r_ome"
399  if (debugJSS==3) cout << "{"<<r_ome<<"}";
400  for (Int_t cChiOmega=0; cChiOmega<=r_ome; cChiOmega++) {
401  for (Int_t cPhiOmega=0; cPhiOmega<=r_ome-cChiOmega;
402  cPhiOmega++) {
403  Int_t cPhiEps=cPsiPhi-cPhiOmega;
404  Int_t cChiEps=cPsiChi-cChiOmega;
405  if (debugJSS==3) cout << "[" << cPhiEps
406  << cChiEps<< r_eps<<"]";
407  if (cPhiEps<0 || cChiEps<0 ||
408  cPhiEps+cChiEps>r_eps) {
409  continue;
410  }
411  else if (debugJSS==3) cout << "E+ OK" << endl;
412  //
413  //
414  //
415  if (debugJSS>=2)
416  cout << "Checking PsiInt=" << PsiInternal
417  << " cPsiChi=" << cPsiChi
418  << " cChiPhi=" << cChiPhi
419  << " cPsiPhi=" << cPsiPhi << endl
420  << " cPhiOmega=" << cPhiOmega
421  << " cChiOmega=" << cChiOmega
422  << " cPhiEps=" << cPhiEps
423  << " cChiEps=" << cChiEps << endl;
424  Int_t cc=0;
425  if (debugJSS>=2) {
426  cout << "Contraction pattern ";
427  cc=cPsiPhi; while (cc--) cout << "#";
428  if (cPsiPhi) cout << "(";
429  cc=cPhiOmega; while (cc--) cout << "o";
430  cc=cPhiEps; while (cc--) cout << "e";
431  if (cPsiPhi) cout << ")";
432  cout<<" "<< SDecay1+SDecay2;
433  cc=PsiInternal; while (cc--) cout << "'";
434  cout << " ";
435  cc=cPsiChi; while (cc--) cout << "#";
436  if (cPsiChi) cout << "(";
437  cc=cChiOmega; while (cc--) cout << "o";
438  cc=cChiEps; while (cc--) cout << "e";
439  if (cPsiChi) cout << ")";
440  cout << " "<< L << " ";
441  cc=cChiPhi; while (cc--) cout << "#";
442  cout << " "<< JMother << " ";
443  cc=cPsiPhi; while (cc--) cout << "#";
444  if (cPsiPhi) cout << "(";
445  cc=cPhiOmega; while (cc--) cout << "o";
446  cc=cPhiEps; while (cc--) cout << "e";
447  if (cPsiPhi) cout << ")";
448 
449  if (preloop) cout << " delta:";
450  else cout << endl;
451  }
452  for (Int_t delta=0; delta<=MaxDelta; delta++) {
453  if (preloop) {
454  if (debugJSS>=2) cout << " " << delta;
455  NLSAmpl++;
456  }
457  else {
458  if (debugJSS>=2) cout << " Constructing LS-Amplitude "
459  << NLSAmpl << endl;
460 
461  Int_t SameContr=0;
462  for (SameContr=0; SameContr<NLSAmpl; SameContr++) {
463  if (LSAmplitudes[SameContr]->
464  CheckContraction(L, S_L, PsiInternal, cPsiChi,
465  cChiPhi, cPsiPhi,
466  cPhiOmega, cChiOmega,
467  cPhiEps, cChiEps)) break;
468  }
469  Int_t ContractionNumber=0;
470  if (SameContr<NLSAmpl)
471  ContractionNumber = LSAmplitudes[SameContr]->
472  GetContraction();
473  else
474  ContractionNumber = MaxContractionNumber+1;
475 
476  LSAmplitudes[NLSAmpl]=
477  new TLSAmpl(SDecay1, SDecay2, L, JMother,
478  delta, S_L,
479  PsiInternal, cPsiChi, cChiPhi, cPsiPhi,
480  cPhiOmega, cChiOmega, cPhiEps, cChiEps,
481  ContractionNumber);
482 
483  if (LSAmplitudes[NLSAmpl]->GetNterms() ) {
484  NLSAmpl++;
485  if (ContractionNumber > MaxContractionNumber)
486  MaxContractionNumber++;
487  }
488  else {
489  delete LSAmplitudes[NLSAmpl];
490  }
491  }
492  }
493  if (debugJSS>=2 && preloop) cout << endl;
494  }
495  }
496  }
497  }
498  }
499  }
500  }
501  }
502  if (preloop) {
503  if (debugJSS)
504  cout << NLSAmpl << " LS-Amplitudes to be evaluated." << endl;
505  }
506  preloop--;
507  }
508  if (debugJSS) {
509  cout << NLSAmpl << " LS-Amplitudes found to be non-zero." << endl;
510  cout << "++++++++++++++++++++++++++++++++++++" << endl;
511  cout << "+++ Helicity-coupling amplitudes +++" << endl;
512  cout << "++++++++++++++++++++++++++++++++++++" << endl;
513  }
514 
515  NFhhAmpl=0;
516  FhhAmpl = new TFhh*[(SDecay1+1)*(SDecay2+1)];
517 
518  for (Int_t lambda=0; lambda<=SDecay1; lambda++)
519  for (Int_t nu=-SDecay2; nu<=SDecay2; nu++) {
520  // for (Int_t nu=-SDecay1; nu<=SDecay2; nu++) { bug!!!! found 4.3.08
521  if (lambda==0 && nu<0) continue;
522  FhhAmpl[NFhhAmpl] = new TFhh(JMother, SDecay1, SDecay2,
523  lambda, nu, NLSAmpl, LSAmplitudes,
524  even_contraction);
525  if (FhhAmpl[NFhhAmpl]->GetNterms()) {
526  NFhhAmpl++;
527  }
528  else {
529  delete FhhAmpl[NFhhAmpl];
530  }
531  }
532 
533  if (debugJSS)
534  cout << NFhhAmpl << " non-zero helicity-coupling amplitudes" << endl;
535 
536  NFhhIdAmpl=0;
537  FhhIdAmpl = new TFhh*[(SDecay1+1)*(SDecay2+1)];
538 
539  if ( SDecay1==SDecay2 && eta1==eta2 ) {
540  if (debugJSS)
541  cout << endl << " for identical-particle decay:" << endl;
542  for (Int_t ifhh=0; ifhh<NFhhAmpl; ifhh++) {
543  FhhIdAmpl[ifhh]=0;
544  if (FhhAmpl[ifhh]) {
545  if (FhhAmpl[ifhh]->IsNuNu()) {
546  FhhIdAmpl[ifhh] = new TFhh(FhhAmpl[ifhh], 'i');
547  }
548  else if (FhhAmpl[ifhh]->IsNuMinusNu() ) {
549  FhhIdAmpl[ifhh] = new TFhh(FhhAmpl[ifhh], 'm');
550  }
551  else {
552  Int_t found_partner=0;
553  for (Int_t jfhh=0; jfhh<NFhhAmpl; jfhh++) {
554  if (FhhAmpl[ifhh]->GetLambda()==FhhAmpl[jfhh]->GetNu() &&
555  FhhAmpl[ifhh]->GetNu()==FhhAmpl[jfhh]->GetLambda() ) {
556  found_partner=1;
557  FhhIdAmpl[ifhh] = new TFhh(FhhAmpl[ifhh], FhhAmpl[jfhh]);
558  // ** continue here **
559  }
560  }
561  if (!found_partner)
562  cerr << "?!?! No partner for amplitude "<< FhhAmpl[ifhh]->GetName()
563  << endl;
564  }
565  }
566  }
567 
568 
569  }
570 
571  cout << NFhhAmpl<< " amplitudes: non-relativistic limit" << endl;
572  for (Int_t i=0; i<NFhhAmpl; i++){
573  FhhAmpl[i]->NonRelLimit();
574  }
575 
576  cout << "Check non-relativistic G's" << endl;
577  for (Int_t i=0; i<NFhhAmpl; i++){
578  FhhAmpl[i]->PrintNRG();
579  }
580 
581  return 0;
582 }
583 
585  char DecayName[10];
586  sprintf(DecayName,"%lld%lld%lld%c%c", JMother, SDecay1, SDecay2,
587  etaJ*eta1*eta2==-1?'n':'p',' ');
588  char ofname[20];
589  sprintf(ofname, "CalcAmpl-%s.h", DecayName);
590  ofstream ofs(ofname);
591  ofs << "// CalcAmpl output for " << DecayName << endl;
592  ofs << "const int FhhAmpl_" << DecayName << "[] = { " << endl;
593  ofs << " " << NFhhAmpl
594  << ", // number of Fhh amplitudes" << endl;
595  for (int i=0;i<NFhhAmpl; i++){
596  ofs << " "
597  << FhhAmpl[i]->GetJ() <<", "
598  << FhhAmpl[i]->GetLambda() <<", "
599  << FhhAmpl[i]->GetNu() <<", "
600  << FhhAmpl[i]->GetEvenContr() <<", // "
601  << "F"
602  << FhhAmpl[i]->GetLambda()
603  << FhhAmpl[i]->GetNu()
604  << ": J, lambda, nu, even_contr"
605  << endl;
606  ofs << " " << FhhAmpl[i]->GetNterms()
607  << ", // number of contributions" << endl;
608  for (int j=0; j<FhhAmpl[i]->GetNterms(); j++) {
609  TLSContrib* lsa = FhhAmpl[i]->GetLStPtr()[j];
610  ofs << " "
611  << lsa->GetJ() <<", "
612  << lsa->GetL() <<", "
613  << lsa->GetS() <<", "
614  << lsa->GetDelta()<<", "
615  << lsa->GetRunningNumber()
616  << ", // contr. F"
617  << FhhAmpl[i]->GetLambda()
618  << FhhAmpl[i]->GetNu() << "-" << j
619  << ": J, L, S, delta, #[g,f,h,...]" << endl;
620  ofs << " "
621  << lsa->GetNterms() << ", "
622  << lsa->GetNormFactor()->GetSign() << ", "
623  << lsa->GetNormFactor()->GetNumerator() << ", "
624  << lsa->GetNormFactor()->GetDenominator()
625  << ", // number of terms, squared norm factor sign/N/D" << endl;
626  for (int k=0; k<lsa->GetNterms(); k++) {
627  ofs << " "
628  << lsa->GetTermFracNum()[k].GetSign() << ", "
629  << lsa->GetTermFracNum()[k].GetNumerator() << ", "
630  << lsa->GetTermFracNum()[k].GetDenominator() << ", "
631  << lsa->GetTermg1pot()[k] << ", "
632  << lsa->GetTermg2pot()[k]
633  << ", // squared sign/N/D, exponents of g_s and g_sigma" << endl;
634  }
635  }
636  }
637  ofs << "};" << endl;
638 
639  return 0;
640 }