19 even_contraction = even_contr_;
23 name_str=
new char[15];
24 sprintf(name_str,
"F_%lld_%lld",
lambda, nu);
26 for (
Int_t iLS=0; iLS<nLS; iLS++) {
27 if (LSampl[iLS]->Getdelta() == delta ||
28 LSampl[iLS]->Getdelta() == -delta ) {
34 cout <<
"Clebsch-Gordan is zero" << endl;
38 for (
Int_t iC=0; iC<Nterms; iC++) newPTR[iC] = LSt[iC];
39 newPTR[Nterms] =
new TLSContrib(LSampl[iLS], delta, SpinCouplFac);
50 if (flag !=
'i' && flag !=
'm'){
51 cerr <<
"TFhh::TFhh unknown flag "<<flag<<endl;
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;
62 name_str =
new char[15];
63 sprintf(name_str,
"%s[symm]", sFhh->
GetName());
75 name_str =
new char[15];
76 sprintf(name_str,
"%s[symm]", sFhh->
GetName());
85 cerr <<
"TFhh::TFhh(TFhh *, TFhh*): Something is wrong," << endl
86 <<
" source amplitudes for symmetrization do not match" << endl;
105 Int_t found_in_prevterm=0;
106 for (
Int_t j=0; j<prevterms; j++) {
107 if ( pLSt[j]->SameParameter(to_be_added) ) {
109 if (
i<Ns) pLSt[j]->
Add(to_be_added,
false);
110 else pLSt[j]->
Add(to_be_added,
true);
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);
125 if ( pLSt[
i]->GetNterms() != 0 ) non_zeros++;
134 if ( pLSt[
i]->GetNterms() != 0 ) {
148 if (! LSt[
i]->IsPureRelativistic()) {
150 for (
Int_t j=0; j<NNRterms; j++) {
151 if (NRLSt[j]->CheckJLS(LSt[
i])){
157 NRLSt[jfound]->Add(LSt[
i]);
162 for (
Int_t j=0; j<NNRterms-1; j++) {
170 cout << name_str <<
" (NR) = " << endl;
171 for (
Int_t j=0; j<NNRterms; j++) {
172 cout <<
"LS=" << NRLSt[j]->
GetL() << NRLSt[j]->GetS() <<
": ";
179 cout << name_str <<
" (NR) = " << endl;
180 for (
Int_t j=0; j<NNRterms; j++) {
181 cout <<
"LS=" << NRLSt[j]->GetL() << NRLSt[j]->GetS() <<
" => ";
190 cout << name_str <<
" =";
191 if (even_contraction) cout << endl;
192 else cout <<
" (iw)" << endl;
193 for (
Int_t iLSt=0; iLSt<Nterms; iLSt++) {
196 if (Nterms==0) cout <<
" 0" << endl;
204 cout <<
" Decay channel: " << JMother;
205 if (etaJ>0) cout <<
"+ -> ";
206 else cout <<
"- -> ";
208 if (eta1>0) cout <<
"+ ";
211 if (eta2>0) cout <<
"+";
217 Smin=SDecay1-SDecay2;
if (Smin<0) Smin=-Smin;
218 Smax=SDecay1+SDecay2;
221 cout <<
"possible S:";
230 Int_t intr_parity = etaJ*eta1*eta2;
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;
239 if (
debugJSS>=2) cout <<
"possible L:";
243 if (intr_parity<0) testL=1;
247 if (testL>=Lmin) NL++;
257 if (
debugJSS>=2) cout <<
" " << fL[NL];
265 Int_t even_contraction=1;
266 if ((SDecay1+SDecay2+testL-JMother) % 2) even_contraction=0;
269 if (even_contraction) cout <<
"contraction only with g~" << endl;
270 else cout <<
"contraction with eps*p and g~" << endl;
273 Int_t MaxPsiInternal=SDecay1;
274 if (SDecay2<SDecay1) MaxPsiInternal=SDecay2;
276 Int_t MaxPsiPhi = SDecay1+SDecay2;
277 if (JMother<MaxPsiPhi) MaxPsiPhi=JMother;
284 while (preloop==1 || preloop==0) {
287 LSAmplitudes =
new TLSAmpl*[NLSAmpl];
289 if (
debugJSS>=2) cout << endl <<
"*************" << endl;
292 for (
Int_t il=0; il<NL; il++) {
296 Int_t MaxPsiChi = SDecay1+SDecay2;
297 if (L<MaxPsiChi) MaxPsiChi=L;
299 Int_t MaxChiPhi = JMother;
300 if (L<JMother) MaxChiPhi=L;
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
311 << SDecay1+SDecay2 <<
" "
313 << JMother <<
"]" << endl;
314 Int_t totalRank=SDecay1+SDecay2+L+JMother;
316 if (JMother<S_L) MaxDelta=JMother;
318 Int_t IndexContractions=(totalRank-3)/2;
319 if (even_contraction) IndexContractions= totalRank / 2;
320 if (
debugJSS>=2) cout << IndexContractions
321 <<
" Lorentz contractions." << endl;
323 Int_t MaxContractionNumber=0;
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++) {
330 cout <<
"Checking " << PsiInternal <<
" " << cPsiChi <<
" "
331 << cChiPhi <<
" " << cPsiPhi;
332 if (PsiInternal+cPsiChi+cChiPhi+cPsiPhi != IndexContractions){
333 if (
debugJSS==3) cout <<
" C-" << endl;
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;
343 if (cPsiChi+cChiPhi != L) {
344 if (
debugJSS==3) cout <<
"L-" << endl;
348 if (cChiPhi+cPsiPhi != JMother) {
349 if (
debugJSS==3) cout <<
"J-" << endl;
355 if (L -cPsiChi-cChiPhi>1) {
356 if (
debugJSS==3) cout <<
"L-" << endl;
360 if (JMother-cPsiPhi-cChiPhi>1) {
361 if (
debugJSS==3) cout <<
"J-" << endl;
366 Int_t r_ome=SDecay1-PsiInternal;
367 Int_t r_eps=SDecay2-PsiInternal;
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;
377 if (r_ome>0) r_ome--;
382 if (r_eps>0) r_eps--;
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;
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) {
411 else if (
debugJSS==3) cout <<
"E+ OK" << endl;
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;
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 <<
"'";
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 <<
")";
449 if (preloop) cout <<
" delta:";
452 for (
Int_t delta=0; delta<=MaxDelta; delta++) {
454 if (
debugJSS>=2) cout <<
" " << delta;
458 if (
debugJSS>=2) cout <<
" Constructing LS-Amplitude "
462 for (SameContr=0; SameContr<NLSAmpl; SameContr++) {
463 if (LSAmplitudes[SameContr]->
464 CheckContraction(L, S_L, PsiInternal, cPsiChi,
466 cPhiOmega, cChiOmega,
467 cPhiEps, cChiEps))
break;
469 Int_t ContractionNumber=0;
470 if (SameContr<NLSAmpl)
471 ContractionNumber = LSAmplitudes[SameContr]->
474 ContractionNumber = MaxContractionNumber+1;
476 LSAmplitudes[NLSAmpl]=
477 new TLSAmpl(SDecay1, SDecay2, L, JMother,
479 PsiInternal, cPsiChi, cChiPhi, cPsiPhi,
480 cPhiOmega, cChiOmega, cPhiEps, cChiEps,
483 if (LSAmplitudes[NLSAmpl]->GetNterms() ) {
485 if (ContractionNumber > MaxContractionNumber)
486 MaxContractionNumber++;
489 delete LSAmplitudes[NLSAmpl];
493 if (
debugJSS>=2 && preloop) cout << endl;
504 cout << NLSAmpl <<
" LS-Amplitudes to be evaluated." << endl;
509 cout << NLSAmpl <<
" LS-Amplitudes found to be non-zero." << endl;
510 cout <<
"++++++++++++++++++++++++++++++++++++" << endl;
511 cout <<
"+++ Helicity-coupling amplitudes +++" << endl;
512 cout <<
"++++++++++++++++++++++++++++++++++++" << endl;
516 FhhAmpl =
new TFhh*[(SDecay1+1)*(SDecay2+1)];
519 for (
Int_t nu=-SDecay2; nu<=SDecay2; nu++) {
521 if (
lambda==0 && nu<0)
continue;
522 FhhAmpl[NFhhAmpl] =
new TFhh(JMother, SDecay1, SDecay2,
523 lambda, nu, NLSAmpl, LSAmplitudes,
525 if (FhhAmpl[NFhhAmpl]->GetNterms()) {
529 delete FhhAmpl[NFhhAmpl];
534 cout << NFhhAmpl <<
" non-zero helicity-coupling amplitudes" << endl;
537 FhhIdAmpl =
new TFhh*[(SDecay1+1)*(SDecay2+1)];
539 if ( SDecay1==SDecay2 && eta1==eta2 ) {
541 cout << endl <<
" for identical-particle decay:" << endl;
542 for (
Int_t ifhh=0; ifhh<NFhhAmpl; ifhh++) {
545 if (FhhAmpl[ifhh]->IsNuNu()) {
546 FhhIdAmpl[ifhh] =
new TFhh(FhhAmpl[ifhh],
'i');
548 else if (FhhAmpl[ifhh]->IsNuMinusNu() ) {
549 FhhIdAmpl[ifhh] =
new TFhh(FhhAmpl[ifhh],
'm');
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() ) {
557 FhhIdAmpl[ifhh] =
new TFhh(FhhAmpl[ifhh], FhhAmpl[jfhh]);
562 cerr <<
"?!?! No partner for amplitude "<< FhhAmpl[ifhh]->GetName()
571 cout << NFhhAmpl<<
" amplitudes: non-relativistic limit" << endl;
573 FhhAmpl[
i]->NonRelLimit();
576 cout <<
"Check non-relativistic G's" << endl;
578 FhhAmpl[
i]->PrintNRG();
586 sprintf(DecayName,
"%lld%lld%lld%c%c", JMother, SDecay1, SDecay2,
587 etaJ*eta1*eta2==-1?
'n':
'p',
' ');
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++){
597 << FhhAmpl[
i]->GetJ() <<
", "
598 << FhhAmpl[
i]->GetLambda() <<
", "
599 << FhhAmpl[
i]->GetNu() <<
", "
600 << FhhAmpl[
i]->GetEvenContr() <<
", // "
602 << FhhAmpl[
i]->GetLambda()
603 << FhhAmpl[
i]->GetNu()
604 <<
": J, lambda, nu, even_contr"
606 ofs <<
" " << FhhAmpl[
i]->GetNterms()
607 <<
", // number of contributions" << endl;
608 for (
int j=0; j<FhhAmpl[
i]->GetNterms(); j++) {
611 << lsa->
GetJ() <<
", "
612 << lsa->
GetL() <<
", "
613 << lsa->
GetS() <<
", "
617 << FhhAmpl[
i]->GetLambda()
618 << FhhAmpl[
i]->GetNu() <<
"-" << j
619 <<
": J, L, S, delta, #[g,f,h,...]" << endl;
625 <<
", // number of terms, squared norm factor sign/N/D" << endl;
633 <<
", // squared sign/N/D, exponents of g_s and g_sigma" << endl;