ROOTPWA
ClebschGordanBox.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include "ClebschGordanBox.h"
4 using namespace std;
5 
6 Int_t debugCG=0;
7 
8 TFracNum* ClebschGordanInteger(Int_t J, Int_t J1, Int_t J2) {
9 
10  Int_t nom_ptr[1]={1};
11  TFracNum two(1, 0, nom_ptr, 0, 1);
12  TFracNum mtwo(1, 0, nom_ptr, 0, -1);
13 
14  Int_t max1=2*J1+1;
15  Int_t max2=2*J2+1;
16  TFracNum *CG=new TFracNum[max1*max2];
17 
18  CG[max2-1]=TFracNum(0,0,0,0,1);
19 
20  for (Int_t m1=-J1; m1<=J1; m1++) {
21  Int_t i1=m1+J1;
22  Int_t m2=0;
23  Int_t m=0;
24  Int_t i2=0;
25  Int_t i=0;
26  for (m2=J2-1; m2>=-J2; m2--) {
27  if (debugCG==2)
28  cout << " Setting m1="<<m1<<", m2="<<m2<<endl;
29  m=m1+m2;
30  i2=m2+J2;
31  i=i1*max2+i2;
32  if (m1+m2<-J||m1+m2>J){
33  CG[i]=TFracNum(0,0,0,0,0);
34  if (debugCG==2) {
35  cout << "#############################################"<<endl;
36  cout << "zero CG["<<i<<"]:" << endl;
37  CG[i].Print();
38  }
39  }
40  else if (m1==-J1) {
41  CG[i] = CG[i+1]*TFracNum((J-m)*(J+m+1), (J2-m2)*(J2+m2+1));
42  if (debugCG==2) {
43  cout << "#############################################"<<endl;
44  cout << " first row CG["<<i<<"]:" << endl;
45  CG[i].Print();
46  }
47  }
48  else if (m1+m2!=-J) {
49  Int_t Denom=(J1+m1)*(J1-m1+1);
50  TFracNum cgp=CG[i-max2]*TFracNum((J+m)*(J-m+1), Denom);
51  TFracNum cgm=CG[i-max2+1]*TFracNum((J2-m2)*(J2+m2+1), Denom);
52  if (debugCG==2) {
53  cout << "cgp:" <<endl; cgp.Print();
54  cout << "cgm:" <<endl; cgm.Print();
55  }
56 
57  TFracNum cgmixed=cgp*cgm;
58  if (cgmixed.Sqrt()) {
59  Bool_t flipsign=cgm>cgp;
60  cgp.Abs();
61  cgm.Abs();
62  CG[i]=cgp+cgm+mtwo*cgmixed;
63  if (flipsign) CG[i].FlipSign();
64  if (debugCG==2) {
65  cout << "#############################################"<<endl;
66  cout << " middle CG Works fine!!! CG["<<i<<"]:" << endl;
67  CG[i].Print();
68  }
69  }
70  else {
71  if (debugCG==2) {
72  cout << "J="<<J<<",J1="<<J1<<",J2="<<J2<<",m="
73  <<m<<",m1="<<m1<<",m2="<<m2<<endl;
74  cgp.Print(); cgm.Print();
75  }
76  return 0;
77  }
78  }
79  }
80  //
81  // special case maximum m2=J2
82  //
83  m2=J2;
84  m=m1+m2;
85  i2=m2+J2;
86  i=i1*max2+i2;
87  if (m1+m2>J) {
88  CG[i]=TFracNum(0,0,0,0,0);
89  if (debugCG==2) {
90  cout << "#############################################"<<endl;
91  cout << "first line Simple setting: CG["<<i<<"]:" << endl;
92  CG[i].Print();
93  }
94  }
95  else if (m1!=-J1) {
96  Int_t Denom=(J-m+1)*(J+m);
97  TFracNum cg1=CG[i-max2]*TFracNum((J1+m1)*(J1-m1+1), Denom);
98  TFracNum cg2=CG[i-1]*TFracNum((J2+m2)*(J2-m2+1), Denom);
99  TFracNum cgmixed=cg1*cg2;
100  if (debugCG==2) {
101  cout << "cg1:" <<endl; cg1.Print();
102  cout << "cg2:" <<endl; cg2.Print();
103  cout << "cgmixed:" <<endl; cgmixed.Print();
104  }
105  if (cgmixed.Sqrt()) {
106  cg1.Abs();
107  cg2.Abs();
108  CG[i]=cg1+cg2+two*cgmixed;
109  if (debugCG==2) {
110  cout << "#############################################"<<endl;
111  cout << " first line Works fine!!! CG["<<i<<"]:" << endl;
112  CG[i].Print();
113  }
114  }
115  else {
116  if (debugCG==2) {
117  cout << "J="<<J<<",J1="<<J1<<",J2="<<J2<<",m="
118  <<m<<",m1="<<m1<<",m2="<<m2<<endl;
119  cg1.Print(); cg2.Print();
120  }
121  return 0;
122  }
123  }
124  //
125  // special case m2=J-m1 (when CG(m1-1,m2)=0)
126  //
127  if (J<J1+J2 && m1!=-J1 && J+m1<=J2) {
128  m2=-J-m1;
129  m=-J;
130  i2=m2+J2;
131  i=i1*max2+i2;
132  Int_t Denom=(J2-m2)*(J2+m2+1);
133  TFracNum cg1=CG[i+1]*TFracNum((J-m)*(J+m+1), Denom);
134  TFracNum cg2=CG[i-max2+1]*TFracNum((J1+m1)*(J1-m1+1), Denom);
135  if (debugCG==2) {
136  cout << "cg1:" <<endl; cg1.Print();
137  cout << "cg2:" <<endl; cg2.Print();
138  }
139  TFracNum cgmixed=cg1*cg2;
140  if (cgmixed.Sqrt()) {
141  Bool_t flipsign=cg2>cg1;
142  cg1.Abs();
143  cg2.Abs();
144  CG[i]=cg1+cg2+mtwo*cgmixed;
145  if (flipsign) CG[i].FlipSign();
146  if (debugCG==2) {
147  cout << "#############################################"<<endl;
148  cout << "lower triangle Works fine!!! CG["<<i<<"]:" << endl;
149  CG[i].Print();
150  }
151  }
152  else {
153  if (debugCG==2) {
154  cout << "J="<<J<<",J1="<<J1<<",J2="<<J2<<",m="
155  <<m<<",m1="<<m1<<",m2="<<m2<<endl;
156  cg1.Print(); cg2.Print();
157  }
158  return 0;
159  }
160  }
161  }
162 
163  Int_t overall_sign=0;
164  TFracNum ZERO(0,0,0,0,0);
165 
166  for (Int_t i=0; i<max1*max2; i++) {
167  Int_t m1=i/max2-J1;
168  Int_t m2=i%max2-J2;
169  if (m1+m2==J && !(CG[i]==ZERO)) {
170  overall_sign=CG[i].GetSign();
171  }
172  }
173 
174  for (Int_t m=-J; m<=J; m++) {
175  TFracNum Sum(0,0,0,0,0);
176  for (Int_t i=0; i<max1*max2; i++) {
177  // i=i1*max2+i2 => i1=i/max2, i2=i%max2
178  // i1=m1+J1, i2=m2+J2 => m1=i/max2-J1, m2=i%max2-J2
179  Int_t m1=i/max2-J1;
180  Int_t m2=i%max2-J2;
181  if (m1+m2==m) {
182  TFracNum CGAbs=CG[i];
183  CGAbs.Abs();
184  Sum=Sum+CGAbs;
185  }
186  }
187  TFracNum SumInv=Sum;
188  SumInv.Invert();
189  for (Int_t i=0; i<max1*max2; i++) {
190  Int_t m1=i/max2-J1;
191  Int_t m2=i%max2-J2;
192  if (m1+m2==m) {
193  CG[i]=CG[i]*SumInv;
194  if (overall_sign==-1) CG[i].FlipSign();
195  }
196  }
197  }
198 
199  return CG;
200 }
201 
202 Int_t debugCGBox=0;
203 
204 ClebschGordanBox box;
205 
206 TFracNum*
207 ClebschGordanBox::GetCG(Int_t J, Int_t J1, Int_t J2) {
208  for (Int_t i=0; i<NCG; i++) {
209  if ( J==CGJ[i] && J1==CGJ1[i] && J2==CGJ2[i] )
210  return CG[i];
211  }
212  if (debugCGBox)
213  cout << "Registering Clebsch-Gordans for J="<<J
214  <<" J1="<<J1<<" J2="<<J2<< endl;
215  NCG++;
216  Int_t *newJ = new Int_t[NCG];
217  Int_t *newJ1 = new Int_t[NCG];
218  Int_t *newJ2 = new Int_t[NCG];
219  TFracNum* *newCG = new TFracNum*[NCG];
220  for (Int_t i=0; i<NCG-1; i++) {
221  newJ [i]=CGJ [i];
222  newJ1[i]=CGJ1[i];
223  newJ2[i]=CGJ2[i];
224  newCG[i]= CG[i];
225  }
226  newJ [NCG-1] = J;
227  newJ1[NCG-1] = J1;
228  newJ2[NCG-1] = J2;
229  newCG[NCG-1] = ClebschGordan(2*J, 2*J1, 2*J2);
230 
231  if (NCG-1) {
232  delete[] CGJ; // wrong delete corrected 27.10.08
233  delete[] CGJ1;
234  delete[] CGJ2;
235  delete[] CG;
236  }
237 
238  CGJ = newJ;
239  CGJ1 = newJ1;
240  CGJ2 = newJ2;
241  CG = newCG;
242 
243  return newCG[NCG-1];
244 }
245 
246 TFracNum*
247 ClebschGordanBox::ClebschGordan(Int_t twoJ, Int_t twoJ1, Int_t twoJ2) {
248 
249  if ( twoJ<0 || twoJ1<0 || twoJ2<0 ||
250  twoJ<twoJ1-twoJ2 || twoJ<twoJ2-twoJ1 || twoJ>twoJ1+twoJ2 ) {
251  cerr << "Clebsch-Gordan-Coefficient only for positive J's"
252  <<" with |J1-J2|<=J<=J1+J2" << endl;
253  }
254 
255  Int_t nom_ptr[1]={1};
256  TFracNum two(1, 0, nom_ptr, 0, 1);
257  TFracNum mtwo(1, 0, nom_ptr, 0, -1);
258 
259  Int_t max1=twoJ1+1;
260  Int_t max2=twoJ2+1;
261 
262  TFracNum *CG=new TFracNum[max1*max2];
263 
264  // first non-zero element
265  Int_t mintwom1=-twoJ1;
266  if (twoJ+twoJ2<twoJ1) mintwom1=-twoJ-twoJ2;
267  Int_t minIndex=(mintwom1+twoJ1)/2*max2 + max2-1;
268  CG[minIndex]=TFracNum(0,0,0,0,1);
269  if (debugCG==2) {
270  cout << "#############################################"<<endl;
271  cout << " first row CG["<<minIndex<<"]:" << endl;
272  CG[minIndex].Print();
273  }
274 
275  for (Int_t twom1=-twoJ1; twom1<=twoJ1; twom1+=2) {
276  Int_t twoi1=twom1+twoJ1;
277  Int_t twom2=0;
278  Int_t twom=0;
279  Int_t twoi2=0;
280  Int_t twoi=0;
281  Int_t i=0;
282  for (twom2=twoJ2-2; twom2>=-twoJ2; twom2-=2) {
283  if (debugCG==2)
284  cout << " Setting twom1="<<twom1<<", twom2="<<twom2<<endl;
285  twom=twom1+twom2;
286  twoi2=twom2+twoJ2;
287  twoi=twoi1*max2+twoi2;
288  i=twoi/2;
289  if ( twom1+twom2<-twoJ || twom1+twom2>twoJ ) {
290  CG[i]=TFracNum(0,0,0,0,0);
291  if (debugCG==2) {
292  cout << "#############################################"<<endl;
293  cout << "zero CG["<<i<<"]:" << endl;
294  CG[i].Print();
295  }
296  }
297  else if (twom1==mintwom1) {
298  CG[i] = CG[i+1]*TFracNum((twoJ-twom)*(twoJ+twom+2),
299  (twoJ2-twom2)*(twoJ2+twom2+2));
300  if (debugCG==2) {
301  cout << "#############################################"<<endl;
302  cout << " first row CG["<<i<<"]:" << endl;
303  CG[i].Print();
304  }
305  }
306  else if (twom1+twom2!=-twoJ) {
307  Int_t Denom=(twoJ1+twom1)*(twoJ1-twom1+2);
308  TFracNum cgp=CG[i-max2]*TFracNum((twoJ+twom)*(twoJ-twom+2), Denom);
309  TFracNum cgm=CG[i-max2+1]*TFracNum((twoJ2-twom2)*(twoJ2+twom2+2),
310  Denom);
311  if (debugCG==2) {
312  cout << "cgp:" <<endl; cgp.Print();
313  cout << "cgm:" <<endl; cgm.Print();
314  }
315 
316  TFracNum cgmixed=cgp*cgm;
317  if (cgmixed.Sqrt()) {
318  Bool_t flipsign=cgm>cgp;
319  cgp.Abs();
320  cgm.Abs();
321  CG[i]=cgp+cgm+mtwo*cgmixed;
322  if (flipsign) CG[i].FlipSign();
323  if (debugCG==2) {
324  cout << "#############################################"<<endl;
325  cout << " middle CG Works fine!!! CG["<<i<<"]:" << endl;
326  CG[i].Print();
327  }
328  }
329  else {
330  if (debugCG==2) {
331  cout << "2*J="<<twoJ<<",2*J1="<<twoJ1<<",2*J2="<<twoJ2<<",2*m="
332  <<twom<<",2*m1="<<twom1<<",2*m2="<<twom2<<endl;
333  cgp.Print(); cgm.Print();
334  }
335  return 0;
336  }
337  }
338  else {
339  CG[i] = CG[i-max2+1]*TFracNum(-(twoJ2-twom2)*(twoJ2+twom2+2),
340  (twoJ1+twom1)*(twoJ1-twom1+2));
341  if (debugCG==2) {
342  cout << "#############################################"<<endl;
343  cout << " first row CG["<<i<<"]:" << endl;
344  CG[i].Print();
345  }
346  }
347  }
348  //
349  // special case maximum m2=J2
350  //
351  twom2=twoJ2;
352  if (debugCG==2)
353  cout << " special setting twom1="<<twom1<<", twom2="<<twom2<<endl;
354  twom=twom1+twom2;
355  twoi2=twom2+twoJ2;
356  twoi=twoi1*max2+twoi2;
357  i=twoi/2;
358  if (twom1+twom2<-twoJ || twom1+twom2>twoJ) {
359  CG[i]=TFracNum(0,0,0,0,0);
360  if (debugCG==2) {
361  cout << "#############################################"<<endl;
362  cout << "first line Simple setting: CG["<<i<<"]:" << endl;
363  CG[i].Print();
364  }
365  }
366  else if (twom1>mintwom1 && twom1<=-mintwom1) {
367  Int_t Denom=(twoJ-twom+2)*(twoJ+twom);
368  TFracNum cg1=CG[i-max2]*TFracNum((twoJ1+twom1)*(twoJ1-twom1+2), Denom);
369  TFracNum cg2=CG[i-1]*TFracNum((twoJ2+twom2)*(twoJ2-twom2+2), Denom);
370  TFracNum cgmixed=cg1*cg2;
371  if (debugCG==2) {
372  cout << "cg1:" <<endl; cg1.Print();
373  cout << "cg2:" <<endl; cg2.Print();
374  cout << "cgmixed:" <<endl; cgmixed.Print();
375  }
376  if (cgmixed.Sqrt()) {
377  cg1.Abs();
378  cg2.Abs();
379  CG[i]=cg1+cg2+two*cgmixed;
380  if (debugCG==2) {
381  cout << "#############################################"<<endl;
382  cout << " first line Works fine!!! CG["<<i<<"]:" << endl;
383  CG[i].Print();
384  }
385  }
386  else {
387  if (debugCG==2) {
388  cout << "2*J="<<twoJ<<",2*J1="<<twoJ1<<",2*J2="<<twoJ2<<",2*m="
389  <<twom<<",2*m1="<<twom1<<",2*m2="<<twom2<<endl;
390  cg1.Print(); cg2.Print();
391  }
392  return 0;
393  }
394  }
395  //
396  // special case m2=J-m1 (when CG(m1-1,m2)=0)
397  //
398  if (twoJ<twoJ1+twoJ2 &&
399  twom1> mintwom1 &&
400  twom1<=-mintwom1 &&
401  twoJ+twom1<=twoJ2) {
402  twom2=-twoJ-twom1;
403  twom=-twoJ;
404  twoi2=twom2+twoJ2;
405  twoi=twoi1*max2+twoi2;
406  i=twoi/2;
407  Int_t Denom=(twoJ2-twom2)*(twoJ2+twom2+2);
408  TFracNum cg1=CG[i+1]*TFracNum((twoJ-twom)*(twoJ+twom+2), Denom);
409  TFracNum cg2=CG[i-max2+1]*TFracNum((twoJ1+twom1)*(twoJ1-twom1+2), Denom);
410  if (debugCG==2) {
411  cout << "cg1:" <<endl; cg1.Print();
412  cout << "cg2:" <<endl; cg2.Print();
413  }
414  TFracNum cgmixed=cg1*cg2;
415  if (cgmixed.Sqrt()) {
416  Bool_t flipsign=cg2>cg1;
417  cg1.Abs();
418  cg2.Abs();
419  CG[i]=cg1+cg2+mtwo*cgmixed;
420  if (flipsign) CG[i].FlipSign();
421  if (debugCG==2) {
422  cout << "#############################################"<<endl;
423  cout << "lower triangle Works fine!!! CG["<<i<<"]:" << endl;
424  CG[i].Print();
425  }
426  }
427  else {
428  if (debugCG==2) {
429  cout << "2*J="<<twoJ<<",2*J1="<<twoJ1<<",2*J2="<<twoJ2<<",2*m="
430  <<twom<<",2*m1="<<twom1<<",2*m2="<<twom2<<endl;
431  cg1.Print(); cg2.Print();
432  }
433  return 0;
434  }
435  }
436  }
437 
438  Int_t overall_sign=0;
439  TFracNum ZERO(0,0,0,0,0);
440 
441  for (Int_t i=0; i<max1*max2; i++) {
442  Int_t twom1=2*(i/max2)-twoJ1;
443  Int_t twom2=2*(i%max2)-twoJ2;
444  if (twom1+twom2==twoJ && !(CG[i]==ZERO)) {
445  overall_sign=CG[i].GetSign();
446  }
447  }
448 
449  for (Int_t twom=-twoJ; twom<=twoJ; twom+=2) {
450  TFracNum Sum(0,0,0,0,0);
451  for (Int_t i=0; i<max1*max2; i++) {
452  // i=i1*max2+i2 => i1=i/max2, i2=i%max2
453  // i1=m1+J1, i2=m2+J2 => m1=i/max2-J1, m2=i%max2-J2
454  Int_t twom1=2*(i/max2)-twoJ1;
455  Int_t twom2=2*(i%max2)-twoJ2;
456  if (twom1+twom2==twom) {
457  TFracNum CGAbs=CG[i];
458  CGAbs.Abs();
459  Sum=Sum+CGAbs;
460  if (debugCG==2)
461  cout << "2m="<<twom <<": adding " << CGAbs.FracString()<< " -> "
462  << Sum.FracString() << endl;
463  }
464  }
465  TFracNum SumInv=Sum;
466  SumInv.Invert();
467  for (Int_t i=0; i<max1*max2; i++) {
468  Int_t twom1=2*(i/max2)-twoJ1;
469  Int_t twom2=2*(i%max2)-twoJ2;
470  if (twom1+twom2==twom) {
471  CG[i]=CG[i]*SumInv;
472  if (overall_sign==-1) CG[i].FlipSign();
473  CG[i].SetINTs();
474  }
475  }
476  }
477  return CG;
478 }
479 
480 Int_t CGIndex(Int_t J1, Int_t m1, Int_t J2, Int_t m2) {
481  return (J1+m1)*(2*J2+1) + J2+m2;
482 }
Generated by   doxygen 1.8.1.2