TEvtGen/EvtGenBase/EvtAmp.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtAmp.cc
  12 //
  13 // Description: Class to manipulate the amplitudes in the decays.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     May 29, 1997         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <math.h>
  24 #include <assert.h>
  25 #include "EvtGenBase/EvtComplex.hh"
  26 #include "EvtGenBase/EvtSpinDensity.hh"
  27 #include "EvtGenBase/EvtAmp.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include "EvtGenBase/EvtId.hh"
  30 #include "EvtGenBase/EvtPDL.hh"
  31 #include "EvtGenBase/EvtParticle.hh"
  32 using std::endl;
  33
  34
  35
  36 EvtAmp::EvtAmp(){
  37   _ndaug=0;
  38   _pstates=0;
  39   _nontrivial=0;
  40 }
  41
  42
  43 EvtAmp::EvtAmp(const EvtAmp& amp){
  44
  45   int i;
  46
  47   _ndaug=amp._ndaug;
  48   _pstates=amp._pstates;
  49   for(i=0;i<_ndaug;i++){
  50     dstates[i]=amp.dstates[i];
  51     _dnontrivial[i]=amp._dnontrivial[i];
  52   }
  53   _nontrivial=amp._nontrivial;
  54
  55   int namp=1;
  56
  57   for(i=0;i<_nontrivial;i++){
  58     _nstate[i]=amp._nstate[i];
  59     namp*=_nstate[i];
  60   }
  61
  62   for(i=0;i<namp;i++){
  63     assert(i<125);
  64     _amp[i]=amp._amp[i];
  65   }
  66
  67 }
  68
  69
  70
  71 void EvtAmp::init(EvtId p,int ndaugs,EvtId *daug){
  72   setNDaug(ndaugs);
  73   int ichild;
  74   int daug_states[100],parstates;
  75   for(ichild=0;ichild<ndaugs;ichild++){
  76
  77     daug_states[ichild]=
  78       EvtSpinType::getSpinStates(EvtPDL::getSpinType(daug[ichild]));
  79
  80   }
  81
  82   parstates=EvtSpinType::getSpinStates(EvtPDL::getSpinType(p));
  83
  84   setNState(parstates,daug_states);
  85
  86 }
  87
  88
  89
  90
  91 void EvtAmp::setNDaug(int n){
  92   _ndaug=n;
  93 }
  94
  95 void EvtAmp::setNState(int parent_states,int *daug_states){
  96
  97   _nontrivial=0;
  98   _pstates=parent_states;
  99
 100   if(_pstates>1) {
 101      _nstate[_nontrivial]=_pstates;
 102      _nontrivial++;
 103   }
 104
 105   int i;
 106
 107   for(i=0;i<_ndaug;i++){
 108     dstates[i]=daug_states[i];
 109     _dnontrivial[i]=-1;
 110     if(daug_states[i]>1) {
 111       _nstate[_nontrivial]=daug_states[i];
 112       _dnontrivial[i]=_nontrivial;
 113       _nontrivial++;
 114     }
 115   }
 116
 117   if (_nontrivial>5) {
 118     report(ERROR,"EvtGen") << "Too many nontrivial states in EvtAmp!"<<endl;
 119   }
 120
 121 }
 122
 123 void EvtAmp::setAmp(int *ind, const EvtComplex& a){
 124
 125   int nstatepad = 1;
 126   int position = ind[0];
 127
 128   for ( int i=1; i<_nontrivial; i++ ) {
 129     nstatepad *= _nstate[i-1];
 130     position += nstatepad*ind[i];
 131   }
 132   assert(position<125);
 133   _amp[position] = a;
 134
 135 }
 136
 137 const EvtComplex& EvtAmp::getAmp(int *ind)const{
 138
 139   int nstatepad = 1;
 140   int position = ind[0];
 141
 142   for ( int i=1; i<_nontrivial; i++ ) {
 143     nstatepad *= _nstate[i-1];
 144     position += nstatepad*ind[i];
 145   }
 146
 147   return _amp[position];
 148 }
 149
 150
 151 EvtSpinDensity EvtAmp::getSpinDensity(){
 152
 153   EvtSpinDensity rho;
 154   rho.setDim(_pstates);
 155
 156   EvtComplex temp;
 157
 158   int i,j,n;
 159
 160   if (_pstates==1) {
 161
 162     if (_nontrivial==0) {
 163
 164        rho.set(0,0,_amp[0]*conj(_amp[0]));
 165        return rho;
 166
 167     }
 168
 169     n=1;
 170
 171     temp = EvtComplex(0.0);
 172
 173     for(i=0;i<_nontrivial;i++){
 174       n*=_nstate[i];
 175     }
 176
 177     for(i=0;i<n;i++){
 178       temp+=_amp[i]*conj(_amp[i]);
 179     }
 180
 181     rho.set(0,0,temp);;
 182
 183     return rho;
 184
 185   }
 186
 187   else{
 188
 189     for(i=0;i<_pstates;i++){
 190       for(j=0;j<_pstates;j++){
 191
 192         temp = EvtComplex(0.0);
 193
 194         int kk;
 195
 196         int allloop = 1;
 197         for (kk=0;kk<(_nontrivial-1); kk++ ) {
 198           allloop *= dstates[kk];
 199         }
 200
 201         for (kk=0; kk<allloop; kk++) {
 202           temp += _amp[_pstates*kk+i]*conj(_amp[_pstates*kk+j]);}
 203
 204         //        if (_nontrivial>3){
 205         //report(ERROR,"EvtGen") << "Can't handle so many states in EvtAmp!"<<endl;
 206         //}
 207
 208         rho.set(i,j,temp);
 209
 210       }
 211     }
 212     return rho;
 213   }
 214
 215 }
 216
 217
 218 EvtSpinDensity EvtAmp::getBackwardSpinDensity(EvtSpinDensity *rho_list){
 219
 220   EvtSpinDensity rho;
 221
 222   rho.setDim(_pstates);
 223
 224   if (_pstates==1){
 225     rho.set(0,0,EvtComplex(1.0,0.0));
 226     return rho;
 227   }
 228
 229   int k;
 230
 231   EvtAmp ampprime;
 232
 233   ampprime=(*this);
 234
 235   for(k=0;k<_ndaug;k++){
 236
 237     if (dstates[k]!=1){
 238       ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
 239     }
 240   }
 241
 242   return ampprime.contract(0,(*this));
 243
 244 }
 245
 246
 247 EvtSpinDensity EvtAmp::getForwardSpinDensity(EvtSpinDensity *rho_list,int i){
 248
 249   EvtSpinDensity rho;
 250
 251   rho.setDim(dstates[i]);
 252
 253   int k;
 254
 255   if (dstates[i]==1) {
 256
 257     rho.set(0,0,EvtComplex(1.0,0.0));
 258
 259     return rho;
 260
 261   }
 262
 263   EvtAmp ampprime;
 264
 265   ampprime=(*this);
 266
 267   if (_pstates!=1){
 268     ampprime=ampprime.contract(0,rho_list[0]);
 269   }
 270
 271   for(k=0;k<i;k++){
 272
 273     if (dstates[k]!=1){
 274       ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
 275     }
 276
 277   }
 278
 279   return ampprime.contract(_dnontrivial[i],(*this));
 280
 281 }
 282
 283 EvtAmp EvtAmp::contract(int k,const EvtSpinDensity& rho){
 284
 285   EvtAmp temp;
 286
 287   int i,j;
 288   temp._ndaug=_ndaug;
 289   temp._pstates=_pstates;
 290   temp._nontrivial=_nontrivial;
 291
 292   for(i=0;i<_ndaug;i++){
 293     temp.dstates[i]=dstates[i];
 294     temp._dnontrivial[i]=_dnontrivial[i];
 295   }
 296
 297   if (_nontrivial==0) {
 298     report(ERROR,"EvtGen")<<"Should not be here EvtAmp!"<<endl;
 299   }
 300
 301
 302   for(i=0;i<_nontrivial;i++){
 303     temp._nstate[i]=_nstate[i];
 304   }
 305
 306
 307   EvtComplex c;
 308
 309   int index[10];
 310   for (i=0;i<10;i++) {
 311      index[i] = 0;
 312   }
 313
 314   int allloop = 1;
 315   int indflag,ii;
 316   for (i=0;i<_nontrivial;i++) {
 317      allloop *= _nstate[i];
 318   }
 319
 320   for ( i=0;i<allloop;i++) {
 321
 322      c = EvtComplex(0.0);
 323      int tempint = index[k];
 324      for (j=0;j<_nstate[k];j++) {
 325        index[k] = j;
 326        c+=rho.get(j,tempint)*getAmp(index);
 327      }
 328      index[k] = tempint;
 329
 330      temp.setAmp(index,c);
 331
 332      indflag = 0;
 333      for ( ii=0;ii<_nontrivial;ii++) {
 334        if ( indflag == 0 ) {
 335          if ( index[ii] == (_nstate[ii]-1) ) {
 336            index[ii] = 0;
 337          }
 338          else {
 339            indflag = 1;
 340            index[ii] += 1;
 341          }
 342        }
 343      }
 344
 345   }
 346   return temp;
 347
 348 }
 349
 350
 351 EvtSpinDensity EvtAmp::contract(int k,const EvtAmp& amp2){
 352
 353   int i,j,l;
 354
 355   EvtComplex temp;
 356   EvtSpinDensity rho;
 357
 358   rho.setDim(_nstate[k]);
 359
 360   int allloop = 1;
 361   int indflag,ii;
 362   for (i=0;i<_nontrivial;i++) {
 363      allloop *= _nstate[i];
 364   }
 365
 366   int index[10];
 367   int index1[10];
 368   //  int l;
 369   for(i=0;i<_nstate[k];i++){
 370
 371     for(j=0;j<_nstate[k];j++){
 372       if (_nontrivial==0) {
 373         report(ERROR,"EvtGen")<<"Should not be here1 EvtAmp!"<<endl;
 374         rho.set(0,0,EvtComplex(1.0,0.0));
 375         return rho;
 376       }
 377
 378       for (ii=0;ii<10;ii++) {
 379         index[ii] = 0;
 380         index1[ii] = 0;
 381       }
 382       index[k] = i;
 383       index1[k] = j;
 384
 385       temp = EvtComplex(0.0);
 386
 387       for ( l=0;l<int(allloop/_nstate[k]);l++) {
 388
 389         temp+=getAmp(index)*conj(amp2.getAmp(index1));
 390         indflag = 0;
 391         for ( ii=0;ii<_nontrivial;ii++) {
 392           if ( ii!= k) {
 393             if ( indflag == 0 ) {
 394               if ( index[ii] == (_nstate[ii]-1) ) {
 395                 index[ii] = 0;
 396                 index1[ii] = 0;
 397               }
 398               else {
 399                 indflag = 1;
 400                 index[ii] += 1;
 401                 index1[ii] += 1;
 402               }
 403             }
 404           }
 405         }
 406       }
 407       rho.set(i,j,temp);
 408
 409     }
 410   }
 411
 412   return rho;
 413 }
 414
 415
 416 EvtAmp EvtAmp::contract(int i, const EvtAmp& a1,const EvtAmp& a2){
 417
 418   //Do we need this method?
 419
 420   assert(a2._pstates>1&&a2._nontrivial==1);
 421   assert(i<=a1._nontrivial);
 422
 423   EvtAmp tmp;
 424   report(DEBUG,"EvtGen") << "EvtAmp::contract not written yet" << endl;
 425   return tmp;
 426
 427 }
 428
 429
 430 void EvtAmp::dump(){
 431
 432   int i,list[10];
 433
 434   report(DEBUG,"EvtGen") << "Number of daugthers:"<<_ndaug<<endl;
 435   report(DEBUG,"EvtGen") << "Number of states of the parent:"<<_pstates<<endl;
 436   report(DEBUG,"EvtGen") << "Number of states on daughters:";
 437   for (i=0;i<_ndaug;i++){
 438     report(DEBUG,"EvtGen") <<dstates[i]<<" ";
 439   }
 440   report(DEBUG,"EvtGen") << endl;
 441   report(DEBUG,"EvtGen") << "Nontrivial index of  daughters:";
 442   for (i=0;i<_ndaug;i++){
 443     report(DEBUG,"EvtGen") <<_dnontrivial[i]<<" ";
 444   }
 445   report(DEBUG,"EvtGen") <<endl;
 446   report(DEBUG,"EvtGen") <<"number of nontrivial states:"<<_nontrivial<<endl;
 447   report(DEBUG,"EvtGen") << "Nontrivial particles number of states:";
 448   for (i=0;i<_nontrivial;i++){
 449     report(DEBUG,"EvtGen") <<_nstate[i]<<" ";
 450   }
 451   report(DEBUG,"EvtGen") <<endl;
 452   report(DEBUG,"EvtGen") <<"Amplitudes:"<<endl;
 453   if (_nontrivial==0){
 454     list[0] = 0;
 455     report(DEBUG,"EvtGen") << getAmp(list) << endl;
 456   }
 457
 458   int allloop[10];
 459   allloop[0]=1;
 460   for (i=0;i<_nontrivial;i++) {
 461     if (i==0){
 462       allloop[i] *= _nstate[i];
 463     }
 464     else{
 465       allloop[i] = allloop[i-1]*_nstate[i];
 466     }
 467   }
 468   int index = 0;
 469   for (i=0;i<allloop[_nontrivial-1];i++) {
 470     report(DEBUG,"EvtGen") << getAmp(list) << " ";
 471     if ( i==allloop[index]-1 ) {
 472       index ++;
 473       report(DEBUG,"EvtGen") << endl;
 474     }
 475   }
 476
 477   report(DEBUG,"EvtGen") << "-----------------------------------"<<endl;
 478
 479 }
 480
 481
 482 void EvtAmp::vertex(const EvtComplex& c){
 483    int list[1];
 484    list[0] = 0;
 485    setAmp(list,c);
 486 }
 487
 488 void EvtAmp::vertex(int i,const EvtComplex& c){
 489    int list[1];
 490    list[0] = i;
 491    setAmp(list,c);
 492 }
 493
 494 void EvtAmp::vertex(int i,int j,const EvtComplex& c){
 495    int list[2];
 496    list[0] = i;
 497    list[1] = j;
 498    setAmp(list,c);
 499 }
 500
 501 void EvtAmp::vertex(int i,int j,int k,const EvtComplex& c){
 502    int list[3];
 503    list[0] = i;
 504    list[1] = j;
 505    list[2] = k;
 506    setAmp(list,c);
 507 }
 508
 509 void EvtAmp::vertex(int *i1,const EvtComplex& c){
 510
 511    setAmp(i1,c);
 512 }
 513
 514
 515 EvtAmp& EvtAmp::operator=(const EvtAmp& amp){
 516
 517   int i;
 518
 519   _ndaug=amp._ndaug;
 520   _pstates=amp._pstates;
 521   for(i=0;i<_ndaug;i++){
 522     dstates[i]=amp.dstates[i];
 523     _dnontrivial[i]=amp._dnontrivial[i];
 524   }
 525   _nontrivial=amp._nontrivial;
 526
 527   int namp=1;
 528
 529   for(i=0;i<_nontrivial;i++){
 530     _nstate[i]=amp._nstate[i];
 531     namp*=_nstate[i];
 532   }
 533
 534   for(i=0;i<namp;i++){
 535     assert(i<125);
 536     _amp[i]=amp._amp[i];
 537   }
 538
 539   return *this;
 540 }
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551