00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <iostream>
00023 #include <math.h>
00024 #include <assert.h>
00025 #include "EvtGenBase/EvtComplex.hh"
00026 #include "EvtGenBase/EvtSpinDensity.hh"
00027 #include "EvtGenBase/EvtAmp.hh"
00028 #include "EvtGenBase/EvtReport.hh"
00029 #include "EvtGenBase/EvtId.hh"
00030 #include "EvtGenBase/EvtPDL.hh"
00031 #include "EvtGenBase/EvtParticle.hh"
00032 using std::endl;
00033 using std::cout;
00034
00035
00036 EvtAmp::EvtAmp(){
00037 _ndaug=0;
00038 _pstates=0;
00039 _nontrivial=0;
00040 }
00041
00042
00043 EvtAmp::EvtAmp(const EvtAmp& amp){
00044
00045 int i;
00046
00047 _ndaug=amp._ndaug;
00048 _pstates=amp._pstates;
00049 for(i=0;i<_ndaug;i++){
00050 dstates[i]=amp.dstates[i];
00051 _dnontrivial[i]=amp._dnontrivial[i];
00052 }
00053 _nontrivial=amp._nontrivial;
00054
00055 int namp=_pstates;
00056
00057 for(i=0;i<_nontrivial;i++){
00058 _nstate[i]=amp._nstate[i];
00059 namp*=_nstate[i];
00060 }
00061
00062 for(i=0;i<namp;i++){
00063 _amp[i]=amp._amp[i];
00064 }
00065
00066 }
00067
00068
00069
00070 void EvtAmp::init(EvtId p,int ndaugs,EvtId *daug){
00071 setNDaug(ndaugs);
00072 int ichild;
00073 int daug_states[100],parstates;
00074 for(ichild=0;ichild<ndaugs;ichild++){
00075
00076 daug_states[ichild]=
00077 EvtSpinType::getSpinStates(EvtPDL::getSpinType(daug[ichild]));
00078
00079 }
00080
00081 parstates=EvtSpinType::getSpinStates(EvtPDL::getSpinType(p));
00082
00083 setNState(parstates,daug_states);
00084
00085 }
00086
00087
00088
00089
00090 void EvtAmp::setNDaug(int n){
00091 _ndaug=n;
00092 }
00093
00094 void EvtAmp::setNState(int parent_states,int *daug_states){
00095
00096 _nontrivial=0;
00097 _pstates=parent_states;
00098
00099 if(_pstates>1) {
00100 _nstate[_nontrivial]=_pstates;
00101 _nontrivial++;
00102 }
00103
00104 int i;
00105
00106 for(i=0;i<_ndaug;i++){
00107 dstates[i]=daug_states[i];
00108 _dnontrivial[i]=-1;
00109 if(daug_states[i]>1) {
00110 _nstate[_nontrivial]=daug_states[i];
00111 _dnontrivial[i]=_nontrivial;
00112 _nontrivial++;
00113 }
00114 }
00115
00116 if (_nontrivial>5) {
00117 report(ERROR,"EvtGen") << "Too many nontrivial states in EvtAmp!"<<endl;
00118 }
00119
00120 }
00121
00122 void EvtAmp::setAmp(int *ind, const EvtComplex& a){
00123
00124 int nstatepad = 1;
00125 int position = ind[0];
00126
00127 for ( int i=1; i<_nontrivial; i++ ) {
00128 nstatepad *= _nstate[i-1];
00129 position += nstatepad*ind[i];
00130 }
00131 _amp[position] = a;
00132
00133 }
00134
00135 const EvtComplex& EvtAmp::getAmp(int *ind)const{
00136
00137 int nstatepad = 1;
00138 int position = ind[0];
00139
00140 for ( int i=1; i<_nontrivial; i++ ) {
00141 nstatepad *= _nstate[i-1];
00142 position += nstatepad*ind[i];
00143 }
00144
00145 return _amp[position];
00146 }
00147
00148
00149 EvtSpinDensity EvtAmp::getSpinDensity(){
00150
00151 EvtSpinDensity rho;
00152 rho.SetDim(_pstates);
00153
00154 EvtComplex temp;
00155
00156 int i,j,n;
00157
00158 if (_pstates==1) {
00159
00160 if (_nontrivial==0) {
00161
00162 rho.Set(0,0,_amp[0]*conj(_amp[0]));
00163 return rho;
00164
00165 }
00166
00167 n=1;
00168
00169 temp = EvtComplex(0.0);
00170
00171 for(i=0;i<_nontrivial;i++){
00172 n*=_nstate[i];
00173 }
00174
00175 for(i=0;i<n;i++){
00176 temp+=_amp[i]*conj(_amp[i]);
00177 }
00178
00179 rho.Set(0,0,temp);;
00180
00181 return rho;
00182
00183 }
00184
00185 else{
00186
00187 for(i=0;i<_pstates;i++){
00188 for(j=0;j<_pstates;j++){
00189
00190 temp = EvtComplex(0.0);
00191
00192 int kk;
00193
00194 int allloop = 1;
00195 for (kk=0;kk<(_nontrivial-1); kk++ ) {
00196 allloop *= dstates[kk];
00197 }
00198
00199 for (kk=0; kk<allloop; kk++) {
00200 temp += _amp[_pstates*kk+i]*conj(_amp[_pstates*kk+j]);}
00201
00202
00203
00204
00205
00206 rho.Set(i,j,temp);
00207
00208 }
00209 }
00210 return rho;
00211 }
00212
00213 }
00214
00215
00216 EvtSpinDensity EvtAmp::getBackwardSpinDensity(EvtSpinDensity *rho_list){
00217
00218 EvtSpinDensity rho;
00219
00220 rho.SetDim(_pstates);
00221
00222 if (_pstates==1){
00223 rho.Set(0,0,EvtComplex(1.0,0.0));
00224 return rho;
00225 }
00226
00227 int k;
00228
00229 EvtAmp ampprime;
00230
00231 ampprime=(*this);
00232
00233 for(k=0;k<_ndaug;k++){
00234
00235 if (dstates[k]!=1){
00236 ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
00237 }
00238 }
00239
00240 return ampprime.contract(0,(*this));
00241
00242 }
00243
00244
00245 EvtSpinDensity EvtAmp::getForwardSpinDensity(EvtSpinDensity *rho_list,int i){
00246
00247 EvtSpinDensity rho;
00248
00249 rho.SetDim(dstates[i]);
00250
00251 int k;
00252
00253 if (dstates[i]==1) {
00254
00255 rho.Set(0,0,EvtComplex(1.0,0.0));
00256
00257 return rho;
00258
00259 }
00260
00261 EvtAmp ampprime;
00262
00263 ampprime=(*this);
00264
00265 if (_pstates!=1){
00266 ampprime=ampprime.contract(0,rho_list[0]);
00267 }
00268
00269 for(k=0;k<i;k++){
00270
00271 if (dstates[k]!=1){
00272 ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
00273 }
00274
00275 }
00276
00277 return ampprime.contract(_dnontrivial[i],(*this));
00278
00279 }
00280
00281 EvtAmp EvtAmp::contract(int k,const EvtSpinDensity& rho){
00282
00283 EvtAmp temp;
00284
00285 int i,j;
00286 temp._ndaug=_ndaug;
00287 temp._pstates=_pstates;
00288 temp._nontrivial=_nontrivial;
00289
00290 for(i=0;i<_ndaug;i++){
00291 temp.dstates[i]=dstates[i];
00292 temp._dnontrivial[i]=_dnontrivial[i];
00293 }
00294
00295 if (_nontrivial==0) {
00296 report(ERROR,"EvtGen")<<"Should not be here EvtAmp!"<<endl;
00297 }
00298
00299
00300 for(i=0;i<_nontrivial;i++){
00301 temp._nstate[i]=_nstate[i];
00302 }
00303
00304
00305 EvtComplex c;
00306
00307 int index[10];
00308 for (i=0;i<10;i++) {
00309 index[i] = 0;
00310 }
00311
00312 int allloop = 1;
00313 int indflag,ii;
00314 for (i=0;i<_nontrivial;i++) {
00315 allloop *= _nstate[i];
00316 }
00317
00318 for ( i=0;i<allloop;i++) {
00319
00320 c = EvtComplex(0.0);
00321 int tempint = index[k];
00322 for (j=0;j<_nstate[k];j++) {
00323 index[k] = j;
00324 c+=rho.Get(j,tempint)*getAmp(index);
00325 }
00326 index[k] = tempint;
00327
00328 temp.setAmp(index,c);
00329
00330 indflag = 0;
00331 for ( ii=0;ii<_nontrivial;ii++) {
00332 if ( indflag == 0 ) {
00333 if ( index[ii] == (_nstate[ii]-1) ) {
00334 index[ii] = 0;
00335 }
00336 else {
00337 indflag = 1;
00338 index[ii] += 1;
00339 }
00340 }
00341 }
00342
00343 }
00344 return temp;
00345
00346 }
00347
00348
00349 EvtSpinDensity EvtAmp::contract(int k,const EvtAmp& amp2){
00350
00351 int i,j,l;
00352
00353 EvtComplex temp;
00354 EvtSpinDensity rho;
00355
00356 rho.SetDim(_nstate[k]);
00357
00358 int allloop = 1;
00359 int indflag,ii;
00360 for (i=0;i<_nontrivial;i++) {
00361 allloop *= _nstate[i];
00362 }
00363
00364 int index[10];
00365 int index1[10];
00366
00367 for(i=0;i<_nstate[k];i++){
00368
00369 for(j=0;j<_nstate[k];j++){
00370 if (_nontrivial==0) {
00371 report(ERROR,"EvtGen")<<"Should not be here1 EvtAmp!"<<endl;
00372 rho.Set(0,0,EvtComplex(1.0,0.0));
00373 return rho;
00374 }
00375
00376 for (ii=0;ii<10;ii++) {
00377 index[ii] = 0;
00378 index1[ii] = 0;
00379 }
00380 index[k] = i;
00381 index1[k] = j;
00382
00383 temp = EvtComplex(0.0);
00384
00385 for ( l=0;l<int(allloop/_nstate[k]);l++) {
00386
00387 temp+=getAmp(index)*conj(amp2.getAmp(index1));
00388 indflag = 0;
00389 for ( ii=0;ii<_nontrivial;ii++) {
00390 if ( ii!= k) {
00391 if ( indflag == 0 ) {
00392 if ( index[ii] == (_nstate[ii]-1) ) {
00393 index[ii] = 0;
00394 index1[ii] = 0;
00395 }
00396 else {
00397 indflag = 1;
00398 index[ii] += 1;
00399 index1[ii] += 1;
00400 }
00401 }
00402 }
00403 }
00404 }
00405 rho.Set(i,j,temp);
00406
00407 }
00408 }
00409
00410 return rho;
00411 }
00412
00413
00414 EvtAmp EvtAmp::contract(int i, const EvtAmp& a1,const EvtAmp& a2){
00415
00416 assert(a2._pstates>1&&a2._nontrivial==1);
00417
00418 EvtAmp tmp;
00419 report(DEBUG,"EvtGen") << "EvtAmp::contract not written yet" << endl;
00420 return tmp;
00421
00422 }
00423
00424
00425 void EvtAmp::dump(){
00426
00427 int i,list[10];
00428
00429 report(DEBUG,"EvtGen") << "Number of daugthers:"<<_ndaug<<endl;
00430 report(DEBUG,"EvtGen") << "Number of states of the parent:"<<_pstates<<endl;
00431 report(DEBUG,"EvtGen") << "Number of states on daughters:";
00432 for (i=0;i<_ndaug;i++){
00433 report(DEBUG,"EvtGen") <<dstates[i]<<" ";
00434 }
00435 report(DEBUG,"EvtGen") << endl;
00436 report(DEBUG,"EvtGen") << "Nontrivial index of daughters:";
00437 for (i=0;i<_ndaug;i++){
00438 report(DEBUG,"EvtGen") <<_dnontrivial[i]<<" ";
00439 }
00440 report(DEBUG,"EvtGen") <<endl;
00441 report(DEBUG,"EvtGen") <<"number of nontrivial states:"<<_nontrivial<<endl;
00442 report(DEBUG,"EvtGen") << "Nontrivial particles number of states:";
00443 for (i=0;i<_nontrivial;i++){
00444 report(DEBUG,"EvtGen") <<_nstate[i]<<" ";
00445 }
00446 report(DEBUG,"EvtGen") <<endl;
00447 report(DEBUG,"EvtGen") <<"Amplitudes:"<<endl;
00448 if (_nontrivial==0){
00449 list[0] = 0;
00450 report(DEBUG,"EvtGen") << getAmp(list) << endl;
00451 }
00452
00453 int allloop[10];
00454 allloop[0]=1;
00455 for (i=0;i<_nontrivial;i++) {
00456 if (i==0){
00457 allloop[i] *= _nstate[i];
00458 }
00459 else{
00460 allloop[i] = allloop[i-1]*_nstate[i];
00461 }
00462 }
00463 int index = 0;
00464 for (i=0;i<allloop[_nontrivial-1];i++) {
00465 report(DEBUG,"EvtGen") << getAmp(list) << " ";
00466 if ( i==allloop[index]-1 ) {
00467 index ++;
00468 report(DEBUG,"EvtGen") << endl;
00469 }
00470 }
00471
00472 report(DEBUG,"EvtGen") << "-----------------------------------"<<endl;
00473
00474 }
00475
00476
00477 void EvtAmp::vertex(const EvtComplex& c){
00478 int list[1];
00479 list[0] = 0;
00480 setAmp(list,c);
00481 }
00482
00483 void EvtAmp::vertex(int i,const EvtComplex& c){
00484 int list[1];
00485 list[0] = i;
00486 setAmp(list,c);
00487 }
00488
00489 void EvtAmp::vertex(int i,int j,const EvtComplex& c){
00490 int list[2];
00491 list[0] = i;
00492 list[1] = j;
00493 setAmp(list,c);
00494 }
00495
00496 void EvtAmp::vertex(int i,int j,int k,const EvtComplex& c){
00497 int list[3];
00498 list[0] = i;
00499 list[1] = j;
00500 list[2] = k;
00501 setAmp(list,c);
00502 }
00503
00504 void EvtAmp::vertex(int *i1,const EvtComplex& c){
00505
00506 setAmp(i1,c);
00507 }
00508
00509
00510 EvtAmp& EvtAmp::operator=(const EvtAmp& amp){
00511
00512 int i;
00513
00514 _ndaug=amp._ndaug;
00515 _pstates=amp._pstates;
00516 for(i=0;i<_ndaug;i++){
00517 dstates[i]=amp.dstates[i];
00518 _dnontrivial[i]=amp._dnontrivial[i];
00519 }
00520 _nontrivial=amp._nontrivial;
00521
00522 int namp=_pstates;
00523
00524 for(i=0;i<_nontrivial;i++){
00525 _nstate[i]=amp._nstate[i];
00526 namp*=_nstate[i];
00527 }
00528
00529 for(i=0;i<namp;i++){
00530 _amp[i]=amp._amp[i];
00531 }
00532
00533 return *this;
00534 }
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545