ROOTPWA
TJwfTensor.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <string>
3 #include "TJwfTensor.h"
4 using namespace std;
5 
6 TTensorTerm::TTensorTerm(char name, Int_t RJ, Int_t* pzm_field,
7  TFracNum* prefac_){
8  Rome=0; ome_pzm=0;
9  Reps=0; eps_pzm=0;
10  Rchi=0; chi_pzm=0;
11  Rphi=0; phi_pzm=0;
12  gam_s_pot=0;
13  gam_sig_pot=0;
14  prefac=*prefac_;
15  if (name == 'o') {Rome=RJ; ome_pzm=pzm_field;}
16  if (name == 'e') {Reps=RJ; eps_pzm=pzm_field;}
17  if (name == 'c') {Rchi=RJ; chi_pzm=pzm_field;}
18  if (name == 'p') {Rphi=RJ; phi_pzm=pzm_field;}
19 }
20 
21 TTensorTerm::TTensorTerm(TTensorTerm *S, TTensorTerm *L,
22  Int_t contractions,
23  Int_t o_share, Int_t e_share,
24  char con_type) {
25 
26  if ( L->Rome || L->Reps ) return;
27  Rome = S->Rome;
28  Reps = S->Reps;
29  Rchi = S->Rchi + L->Rchi;
30  Rphi = S->Rphi + L->Rphi;
31  gam_s_pot = S->gam_s_pot +L->gam_s_pot;
32  gam_sig_pot = S->gam_sig_pot+L->gam_sig_pot;
33  prefac=S->prefac*L->prefac;
34 
35  Int_t ocount=o_share;
36  Int_t ecount=e_share;
37  for (Int_t con=0; con<contractions; con++) {
38  Int_t cp = 0;
39  if (con_type=='c') { cp=L->chi_pzm[Rchi-1]; Rchi--; }
40  else { cp=L->phi_pzm[Rphi-1]; Rphi--; }
41  Int_t oe = 0;
42  if (ocount) { // o is to be contracted
43  oe = S->ome_pzm[Rome-1];
44  Rome--;
45  }
46  else { // e is to be contracted
47  oe = S->eps_pzm[Reps-1];
48  ecount--;
49  Reps--;
50  }
51  if ( con_type == 'c' && ( (oe==1&&cp==-1) || (oe==-1&&cp==1) ) )
52  prefac.FlipSign();
53  else if
54  ( con_type == 'p' && ( (oe==1&&cp==1) || (oe==-1&&cp==-1) ) ) ;
55  else if (oe==0&&cp==0 ) {
56  if (ocount) gam_s_pot++;
57  else gam_sig_pot++;
58  }
59  else {
60  prefac=TFracNum(0,0,0,0,0);
61  return;
62  }
63  if (ocount) ocount--;
64  }
65 
66  ome_pzm = new Int_t[Rome];
67  for (Int_t i=0; i<Rome; i++) ome_pzm[i] = S->ome_pzm[i];
68  eps_pzm = new Int_t[Reps];
69  for (Int_t i=0; i<Reps; i++) eps_pzm[i] = S->eps_pzm[i];
70  chi_pzm = new Int_t[Rchi];
71  for (Int_t i=0; i<Rchi; i++) {
72  if (i < L->Rchi) chi_pzm[i] = L->chi_pzm[i];
73  else chi_pzm[i] = S->chi_pzm[i-(L->Rchi)];
74  }
75  phi_pzm = new Int_t[Rphi];
76  for (Int_t i=0; i<Rphi; i++){
77  if (i < L->Rphi) phi_pzm[i] = L->phi_pzm[i];
78  else phi_pzm[i] = S->phi_pzm[i-(L->Rphi)];
79  }
80 }
81 
82 Int_t
83 TTensorTerm::LJContraction(Int_t ncon, Int_t even) {
84 
85  for (Int_t con=0; con<ncon; con++) {
86  Int_t c=chi_pzm[Rchi-1];
87  Int_t p=phi_pzm[Rphi-1];
88  if ( c != p ) {
89  prefac=TFracNum(0,0,0,0,0);
90  return 0;
91  }
92  Rchi--;
93  Rphi--;
94  }
95 
96  Int_t trouble=0;
97 
98  if (even==0) {
99  if (Rome+Reps+Rchi+Rphi!=3) {
100  cerr << "TTensorTerm::LJContraction:"
101  <<" Contraction ended with wrong number of indices!!" <<endl;
102  return 0;
103  }
104  else { // eq. (5.21) - (5.23)
105  Int_t os=-100; if (Rome==1) os=ome_pzm[0];
106  Int_t es=-100; if (Reps==1) es=eps_pzm[0];
107  Int_t cs=-100; if (Rchi==1) cs=chi_pzm[0];
108  Int_t ps=-100; if (Rphi==1) ps=phi_pzm[0];
109  if ( Rome==0 ) {
110  Reps--; Rchi--; Rphi--;
111  if (es== 1 && cs==-1 && ps== 0 ) ;
112  else if (es==-1 && cs== 1 && ps== 0 ) prefac.FlipSign();
113  else if (es== 1 && cs== 0 && ps== 1 ) ;
114  else if (es==-1 && cs== 0 && ps==-1 ) prefac.FlipSign();
115  else if (es== 0 && cs== 1 && ps== 1 ){prefac.FlipSign();gam_sig_pot++;}
116  else if (es== 0 && cs==-1 && ps==-1 ){ gam_sig_pot++;}
117  else { prefac=TFracNum(0,0,0,0,0); return 1;}
118  }
119  else if ( Reps==0 ) {
120  Rome--; Rchi--; Rphi--;
121  if (os== 1 && cs==-1 && ps== 0 ) ;
122  else if (os==-1 && cs== 1 && ps== 0 ) prefac.FlipSign();
123  else if (os== 1 && cs== 0 && ps== 1 ) ;
124  else if (os==-1 && cs== 0 && ps==-1 ) prefac.FlipSign();
125  else if (os== 0 && cs== 1 && ps== 1 ){prefac.FlipSign();gam_s_pot++;}
126  else if (os== 0 && cs==-1 && ps==-1 ){ gam_s_pot++;}
127  else { prefac=TFracNum(0,0,0,0,0); return 1;}
128  }
129  else if ( Rchi==0 ) {
130  Rome--; Reps--; Rphi--;
131  if (os== 1 && es==-1 && ps== 0 ) ;
132  else if (os==-1 && es== 1 && ps== 0 ) prefac.FlipSign();
133  else if (os== 1 && es== 0 && ps== 1 ){ gam_sig_pot++;}
134  else if (os==-1 && es== 0 && ps==-1 ){prefac.FlipSign();gam_sig_pot++;}
135  else if (os== 0 && es== 1 && ps== 1 ){prefac.FlipSign();gam_s_pot++;}
136  else if (os== 0 && es==-1 && ps==-1 ){ gam_s_pot++;}
137  else { prefac=TFracNum(0,0,0,0,0); return 1;}
138  }
139  else if ( Rphi==0 ) {
140  Rome--; Reps--; Rchi--;
141  if (os== 1 && es==-1 && cs== 0 ) ;
142  else if (os==-1 && es== 1 && cs== 0 ) prefac.FlipSign();
143  else if (os== 1 && es== 0 && cs==-1 ){prefac.FlipSign();gam_sig_pot++;}
144  else if (os==-1 && es== 0 && cs== 1 ){ gam_sig_pot++;}
145  else if (os== 0 && es== 1 && cs==-1 ){ gam_s_pot++;}
146  else if (os== 0 && es==-1 && cs== 1 ){prefac.FlipSign();gam_s_pot++;}
147  else { prefac=TFracNum(0,0,0,0,0); return 1;}
148  }
149  else { trouble=5;}
150  }
151  }
152 
153  else {
154  if (Rome!=0 || Reps!=0 || Rchi!=0 || Rphi!=0 ) trouble=1;
155  }
156  if (trouble){
157  cerr << "TTensorTerm::LJContraction: "
158  << "Invalid espilon-contraction occurred "<< trouble << endl;
159  return 0;
160  }
161  return 1;
162 }
163 
164 Int_t
165 TTensorTerm::Multiply(char name, Int_t RJ, Int_t* pzm_field,
166  TFracNum* prefac_){
167 
168  prefac = prefac * *prefac_;
169 
170  Int_t Merr=0;
171  if (name == 'o') {if (Rome) Merr=1; else {Rome=RJ; ome_pzm=pzm_field;}}
172  if (name == 'e') {if (Reps) Merr=1; else {Reps=RJ; eps_pzm=pzm_field;}}
173  if (name == 'c') {if (Rchi) Merr=1; else {Rchi=RJ; chi_pzm=pzm_field;}}
174  if (name == 'p') {if (Rphi) Merr=1; else {Rphi=RJ; phi_pzm=pzm_field;}}
175  if (Merr) {
176  cerr << "TTensorTerm::Multiply: Each type can be multiplied only once!"
177  <<endl;
178  return 0;
179  }
180  return 1;
181 }
182 
183 Int_t
184 TTensorTerm::SpinInnerContraction(Int_t cPsiInt) {
185  Int_t res=0;
186  for (Int_t ic=0; ic<cPsiInt; ic++) {
187  res=0;
188  Int_t o = ome_pzm[Rome-1];
189  Int_t e = eps_pzm[Reps-1];
190  if ( (o==1&&e==-1) || (o==-1&&e==1) ) {
191  prefac.FlipSign();
192  res=1;
193  }
194  if (o==0&&e==0) {
195  gam_s_pot += 1;
196  gam_sig_pot += 1;
197  res=1;
198  }
199  Rome--;
200  Reps--;
201  }
202  return res;
203 }
204 
205 Int_t
206 TTensorTerm::SameStructure(TTensorTerm *other) {
207  if (Rome==other->Rome &&
208  Reps==other->Reps &&
209  Rchi==other->Rchi &&
210  Rphi==other->Rphi &&
211  gam_s_pot==other->gam_s_pot &&
212  gam_sig_pot==other->gam_sig_pot) return 1;
213  return 0;
214 }
215 
216 Int_t
217 TTensorTerm::AddTwoTerms(TTensorTerm *other) {
218  if (!SameStructure(other)) {
219  cerr << "NO NO NO these terms cannot be added!" << endl;
220  return 0;
221  }
222  else {
223  TFracNum *sum = prefac.SumSignedRoots(&(other->prefac));
224  if (sum) {
225  prefac= *sum;
226  return 1;
227  }
228  }
229  return 0;
230 }
231 
232 Int_t
233 TTensorTerm::Print(char flag) {
234  if (flag=='s') cout << prefac.FracStringSqrt()<<" ";
235  else cout << prefac.FracString()<<" ";
236  if (Rome) {
237  cout << "o(";
238  for (Int_t i=0; i<Rome; i++) {
239  if (ome_pzm[i]== 1) cout <<"+";
240  if (ome_pzm[i]== 0) cout <<"0";
241  if (ome_pzm[i]==-1) cout <<"-";
242  }
243  cout << ")";
244  }
245  if (Reps) {
246  cout << "e(";
247  for (Int_t i=0; i<Reps; i++) {
248  if (eps_pzm[i]== 1) cout <<"+";
249  if (eps_pzm[i]== 0) cout <<"0";
250  if (eps_pzm[i]==-1) cout <<"-";
251  }
252  cout << ")";
253  }
254  if (Rchi) {
255  cout << "c(";
256  for (Int_t i=0; i<Rchi; i++) {
257  if (chi_pzm[i]== 1) cout <<"+";
258  if (chi_pzm[i]== 0) cout <<"0";
259  if (chi_pzm[i]==-1) cout <<"-";
260  }
261  cout << ")";
262  }
263  if (Rphi) {
264  cout << "p(";
265  for (Int_t i=0; i<Rphi; i++) {
266  if (phi_pzm[i]== 1) cout <<"+";
267  if (phi_pzm[i]== 0) cout <<"0";
268  if (phi_pzm[i]==-1) cout <<"-";
269  }
270  cout << ")";
271  }
272  if (gam_s_pot) {
273  if (gam_s_pot==1) cout << " gs";
274  else cout << " gs^" << gam_s_pot;
275  }
276  if (gam_sig_pot) {
277  if (gam_sig_pot==1) cout << " gsig";
278  else cout << " gsig^" << gam_sig_pot;
279  }
280  return 0;
281 }
282 
283 //Int_t
284 //TTensorSum::Print(char flag='n'){ // CINT limitation for overloading
285 Int_t
286 TTensorSum::Print(char flag){
287  for (Int_t i=0; i<Nterms; i++) {
288  terms[i].Print(flag);
289  if (i<Nterms-1) cout << " ";
290  }
291  cout << endl;
292  return 0;
293 };
294 
295 Int_t
296 TTensorSum::AddTerm (TTensorTerm* addt) {
297  TTensorTerm* nt=new TTensorTerm[Nterms+1];
298  for (Int_t i=0; i<Nterms; i++)
299  nt[i]=terms[i];
300  nt[Nterms]= *addt;
301  delete[] terms;
302  terms=nt;
303  Nterms++;
304  return 0;
305 }
306 
307 Int_t
308 TTensorSum::SpinInnerContraction(Int_t cPsiInt) {
309  Int_t non_zero_terms=0;
310  Int_t val[Nterms];
311  for (Int_t i=0; i<Nterms; i++) {
312  val[i] = terms[i].SpinInnerContraction(cPsiInt);
313  if (val[i]!=0) non_zero_terms++;
314  }
315  if (non_zero_terms<Nterms) {
316  TTensorTerm* nt=new TTensorTerm[non_zero_terms];
317  Int_t j=0;
318  for (Int_t i=0; i<Nterms; i++)
319  if (val[i]) {
320  nt[j]=terms[i];
321  j++;
322  }
323  delete[] terms;
324  terms=nt;
325  Nterms=non_zero_terms;
326  }
327  return Nterms;
328 }
329 
330 TTensorSum*
331 TTensorSum::LSContraction(TTensorSum *L, Int_t contr,
332  Int_t co, Int_t ce, char con_type) {
333 
334  TTensorSum *tls = new TTensorSum();
335 
336  for (Int_t i=0; i<Nterms; i++) {
337  for (Int_t j=0; j<L->Nterms; j++) {
338  TTensorTerm *nt = new TTensorTerm(&(terms[i]), &(L->terms[j]),
339  contr, co, ce, con_type);
340  if ( nt->IsNonZero() ) tls->AddTerm(nt);
341  }
342  }
343 
344  return tls;
345 }
346 
347 TTensorSum*
348 TTensorSum::LJContraction(Int_t cChiPhi, Int_t even) {
349 
350  for (Int_t i=0; i<Nterms; i++) terms[i].LJContraction(cChiPhi, even);
351 
352  TTensorSum *tls = new TTensorSum();
353 
354  for (Int_t i=0; i<Nterms; i++) {
355  Int_t found_same=0;
356  for (Int_t j=0; j<tls->Nterms; j++) {
357  if (terms[i].SameStructure(&(tls->terms[j]))) {
358  if ( (tls->terms[j]).AddTwoTerms(&(terms[i])) ) {
359  found_same=1;
360  break;
361  }
362  }
363  }
364  if (!found_same) {
365  TTensorTerm *nt = new TTensorTerm(terms[i]);
366  tls->AddTerm(nt);
367  }
368  }
369 
370  TTensorSum *tls_nonzero = new TTensorSum();
371 
372  for (Int_t i=0; i<tls->Nterms; i++) {
373  if ( tls->terms[i].IsNonZero() ) {
374  TTensorTerm *nt = new TTensorTerm(tls->terms[i]);
375  tls_nonzero->AddTerm(nt);
376  }
377  }
378 
379  return tls_nonzero;
380 }
Generated by   doxygen 1.8.1.2