TEvtGen/EvtGenBase/EvtEvalHelAmp.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) 2000      Caltech, UCSB
  10 //
  11 // Module: EvtHelAmp.cc
  12 //
  13 // Description: Decay model for implementation of generic 2 body
  14 //              decay specified by the partial wave amplitudes
  15 //
  16 //
  17 // Modification history:
  18 //
  19 //    fkw        February 2, 2001     changes to satisfy KCC
  20 //    RYD       September 7, 2000       Module created
  21 //
  22 //------------------------------------------------------------------------
  23 //
  24 #include "EvtGenBase/EvtPatches.hh"
  25 #include <stdlib.h>
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtVector4C.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30 #include "EvtGenBase/EvtEvalHelAmp.hh"
  31 #include "EvtGenBase/EvtId.hh"
  32 #include "EvtGenBase/EvtConst.hh"
  33 #include "EvtGenBase/EvtdFunction.hh"
  34 #include "EvtGenBase/EvtAmp.hh"
  35 using std::endl;
  36
  37
  38 EvtEvalHelAmp::~EvtEvalHelAmp() {
  39
  40   //deallocate memory
  41   delete [] _lambdaA2;
  42   delete [] _lambdaB2;
  43   delete [] _lambdaC2;
  44
  45   int ia,ib,ic;
  46   for(ib=0;ib<_nB;ib++){
  47     delete [] _HBC[ib];
  48   }
  49
  50   delete [] _HBC;
  51
  52
  53   for(ia=0;ia<_nA;ia++){
  54     delete [] _RA[ia];
  55   }
  56   delete [] _RA;
  57
  58   for(ib=0;ib<_nB;ib++){
  59     delete [] _RB[ib];
  60   }
  61   delete [] _RB;
  62
  63   for(ic=0;ic<_nC;ic++){
  64     delete [] _RC[ic];
  65   }
  66   delete [] _RC;
  67
  68
  69   for(ia=0;ia<_nA;ia++){
  70     for(ib=0;ib<_nB;ib++){
  71       delete [] _amp[ia][ib];
  72       delete [] _amp1[ia][ib];
  73       delete [] _amp3[ia][ib];
  74     }
  75     delete [] _amp[ia];
  76     delete [] _amp1[ia];
  77     delete [] _amp3[ia];
  78   }
  79
  80   delete [] _amp;
  81   delete [] _amp1;
  82   delete [] _amp3;
  83
  84 }
  85
  86
  87 EvtEvalHelAmp::EvtEvalHelAmp(EvtId idA,
  88                              EvtId idB,
  89                              EvtId idC,
  90                              EvtComplexPtrPtr HBC){
  91
  92
  93   EvtSpinType::spintype typeA=EvtPDL::getSpinType(idA);
  94   EvtSpinType::spintype typeB=EvtPDL::getSpinType(idB);
  95   EvtSpinType::spintype typeC=EvtPDL::getSpinType(idC);
  96
  97   //find out how many states each particle have
  98   _nA=EvtSpinType::getSpinStates(typeA);
  99   _nB=EvtSpinType::getSpinStates(typeB);
 100   _nC=EvtSpinType::getSpinStates(typeC);
 101
 102   //find out what 2 times the spin is
 103   _JA2=EvtSpinType::getSpin2(typeA);
 104   _JB2=EvtSpinType::getSpin2(typeB);
 105   _JC2=EvtSpinType::getSpin2(typeC);
 106
 107
 108   //allocate memory
 109   _lambdaA2=new int[_nA];
 110   _lambdaB2=new int[_nB];
 111   _lambdaC2=new int[_nC];
 112
 113   _HBC=new EvtComplexPtr[_nB];
 114   int ia,ib,ic;
 115   for(ib=0;ib<_nB;ib++){
 116     _HBC[ib]=new EvtComplex[_nC];
 117   }
 118
 119
 120   _RA=new EvtComplexPtr[_nA];
 121   for(ia=0;ia<_nA;ia++){
 122     _RA[ia]=new EvtComplex[_nA];
 123   }
 124   _RB=new EvtComplexPtr[_nB];
 125   for(ib=0;ib<_nB;ib++){
 126     _RB[ib]=new EvtComplex[_nB];
 127   }
 128   _RC=new EvtComplexPtr[_nC];
 129   for(ic=0;ic<_nC;ic++){
 130     _RC[ic]=new EvtComplex[_nC];
 131   }
 132
 133   _amp=new EvtComplexPtrPtr[_nA];
 134   _amp1=new EvtComplexPtrPtr[_nA];
 135   _amp3=new EvtComplexPtrPtr[_nA];
 136   for(ia=0;ia<_nA;ia++){
 137     _amp[ia]=new EvtComplexPtr[_nB];
 138     _amp1[ia]=new EvtComplexPtr[_nB];
 139     _amp3[ia]=new EvtComplexPtr[_nB];
 140     for(ib=0;ib<_nB;ib++){
 141       _amp[ia][ib]=new EvtComplex[_nC];
 142       _amp1[ia][ib]=new EvtComplex[_nC];
 143       _amp3[ia][ib]=new EvtComplex[_nC];
 144     }
 145   }
 146
 147   //find the allowed helicities (actually 2*times the helicity!)
 148
 149   fillHelicity(_lambdaA2,_nA,_JA2,idA);
 150   fillHelicity(_lambdaB2,_nB,_JB2,idB);
 151   fillHelicity(_lambdaC2,_nC,_JC2,idC);
 152
 153   for(ib=0;ib<_nB;ib++){
 154     for(ic=0;ic<_nC;ic++){
 155       _HBC[ib][ic]=HBC[ib][ic];
 156     }
 157   }
 158 }
 159
 160
 161
 162
 163
 164
 165 double EvtEvalHelAmp::probMax(){
 166
 167   double c=1.0/sqrt(4*EvtConst::pi/(_JA2+1));
 168
 169   int ia,ib,ic;
 170
 171
 172   double theta;
 173   int itheta;
 174
 175   double maxprob=0.0;
 176
 177   for(itheta=-10;itheta<=10;itheta++){
 178     theta=acos(0.099999*itheta);
 179     for(ia=0;ia<_nA;ia++){
 180       double prob=0.0;
 181       for(ib=0;ib<_nB;ib++){
 182         for(ic=0;ic<_nC;ic++){
 183           _amp[ia][ib][ic]=0.0;
 184           if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
 185             _amp[ia][ib][ic]=c*_HBC[ib][ic]*
 186               EvtdFunction::d(_JA2,_lambdaA2[ia],
 187                               _lambdaB2[ib]-_lambdaC2[ic],theta);
 188             prob+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
 189           }
 190         }
 191       }
 192
 193       prob*=sqrt(1.0*_nA);
 194
 195       if (prob>maxprob) maxprob=prob;
 196
 197     }
 198   }
 199
 200   return maxprob;
 201
 202 }
 203
 204
 205 void EvtEvalHelAmp::evalAmp( EvtParticle *p, EvtAmp& amp){
 206
 207   //find theta and phi of the first daughter
 208
 209   EvtVector4R pB=p->getDaug(0)->getP4();
 210
 211   double theta=acos(pB.get(3)/pB.d3mag());
 212   double phi=atan2(pB.get(2),pB.get(1));
 213
 214   double c=sqrt((_JA2+1)/(4*EvtConst::pi));
 215
 216   int ia,ib,ic;
 217
 218   double prob1=0.0;
 219
 220   for(ia=0;ia<_nA;ia++){
 221     for(ib=0;ib<_nB;ib++){
 222       for(ic=0;ic<_nC;ic++){
 223         _amp[ia][ib][ic]=0.0;
 224         if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
 225           double dfun=EvtdFunction::d(_JA2,_lambdaA2[ia],
 226                                       _lambdaB2[ib]-_lambdaC2[ic],theta);
 227
 228           _amp[ia][ib][ic]=c*_HBC[ib][ic]*
 229             exp(EvtComplex(0.0,phi*0.5*(_lambdaA2[ia]-_lambdaB2[ib]+
 230                                         _lambdaC2[ic])))*dfun;
 231         }
 232         prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
 233       }
 234     }
 235   }
 236
 237   setUpRotationMatrices(p,theta,phi);
 238
 239   applyRotationMatrices();
 240
 241   double prob2=0.0;
 242
 243   for(ia=0;ia<_nA;ia++){
 244     for(ib=0;ib<_nB;ib++){
 245       for(ic=0;ic<_nC;ic++){
 246         prob2+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
 247         if (_nA==1){
 248           if (_nB==1){
 249             if (_nC==1){
 250               amp.vertex(_amp[ia][ib][ic]);
 251             }
 252             else{
 253               amp.vertex(ic,_amp[ia][ib][ic]);
 254             }
 255           }
 256           else{
 257             if (_nC==1){
 258               amp.vertex(ib,_amp[ia][ib][ic]);
 259             }
 260             else{
 261               amp.vertex(ib,ic,_amp[ia][ib][ic]);
 262             }
 263           }
 264         }else{
 265           if (_nB==1){
 266             if (_nC==1){
 267               amp.vertex(ia,_amp[ia][ib][ic]);
 268             }
 269             else{
 270               amp.vertex(ia,ic,_amp[ia][ib][ic]);
 271             }
 272           }
 273           else{
 274             if (_nC==1){
 275               amp.vertex(ia,ib,_amp[ia][ib][ic]);
 276             }
 277             else{
 278               amp.vertex(ia,ib,ic,_amp[ia][ib][ic]);
 279             }
 280           }
 281         }
 282       }
 283     }
 284   }
 285
 286   if (fabs(prob1-prob2)>0.000001*prob1){
 287     report(INFO,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl;
 288     ::abort();
 289   }
 290
 291   return ;
 292
 293 }
 294
 295
 296 void EvtEvalHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){
 297
 298   int i;
 299
 300   //photon is special case!
 301   if (n==2&&J2==2) {
 302     lambda2[0]=2;
 303     lambda2[1]=-2;
 304     return;
 305   }
 306
 307   //and so is the neutrino!
 308   if (n==1&&J2==1) {
 309     if (EvtPDL::getStdHep(id)>0){
 310         //particle i.e. lefthanded
 311         lambda2[0]=-1;
 312     }else{
 313         //anti particle i.e. righthanded
 314         lambda2[0]=1;
 315     }
 316     return;
 317   }
 318
 319
 320   assert(n==J2+1);
 321
 322   for(i=0;i<n;i++){
 323     lambda2[i]=n-i*2-1;
 324   }
 325
 326   return;
 327
 328 }
 329
 330
 331 void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){
 332
 333   switch(_JA2){
 334
 335   case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
 336
 337     {
 338
 339       EvtSpinDensity R=p->rotateToHelicityBasis();
 340
 341
 342       int i,j,n;
 343
 344       n=R.getDim();
 345
 346       assert(n==_nA);
 347
 348
 349       for(i=0;i<n;i++){
 350         for(j=0;j<n;j++){
 351           _RA[i][j]=R.get(i,j);
 352         }
 353       }
 354
 355     }
 356
 357     break;
 358
 359   default:
 360     report(ERROR,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl;
 361     ::abort();
 362   }
 363
 364
 365   switch(_JB2){
 366
 367
 368   case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
 369
 370     {
 371
 372       int i,j,n;
 373
 374       EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi);
 375
 376       n=R.getDim();
 377
 378       assert(n==_nB);
 379
 380       for(i=0;i<n;i++){
 381         for(j=0;j<n;j++){
 382           _RB[i][j]=conj(R.get(i,j));
 383         }
 384       }
 385
 386     }
 387
 388     break;
 389
 390   default:
 391     report(ERROR,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl;
 392     ::abort();
 393   }
 394
 395   switch(_JC2){
 396
 397   case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
 398
 399     {
 400
 401       int i,j,n;
 402
 403       EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi);
 404
 405       n=R.getDim();
 406
 407       assert(n==_nC);
 408
 409       for(i=0;i<n;i++){
 410         for(j=0;j<n;j++){
 411           _RC[i][j]=conj(R.get(i,j));
 412         }
 413       }
 414
 415     }
 416
 417     break;
 418
 419   default:
 420     report(ERROR,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl;
 421     ::abort();
 422   }
 423
 424
 425
 426 }
 427
 428
 429 void EvtEvalHelAmp::applyRotationMatrices(){
 430
 431   int ia,ib,ic,i;
 432
 433   EvtComplex temp;
 434
 435
 436
 437   for(ia=0;ia<_nA;ia++){
 438     for(ib=0;ib<_nB;ib++){
 439       for(ic=0;ic<_nC;ic++){
 440         temp=0;
 441         for(i=0;i<_nC;i++){
 442           temp+=_RC[i][ic]*_amp[ia][ib][i];
 443         }
 444         _amp1[ia][ib][ic]=temp;
 445       }
 446     }
 447   }
 448
 449
 450
 451   for(ia=0;ia<_nA;ia++){
 452     for(ic=0;ic<_nC;ic++){
 453       for(ib=0;ib<_nB;ib++){
 454         temp=0;
 455         for(i=0;i<_nB;i++){
 456           temp+=_RB[i][ib]*_amp1[ia][i][ic];
 457         }
 458         _amp3[ia][ib][ic]=temp;
 459       }
 460     }
 461   }
 462
 463
 464
 465   for(ib=0;ib<_nB;ib++){
 466     for(ic=0;ic<_nC;ic++){
 467       for(ia=0;ia<_nA;ia++){
 468         temp=0;
 469         for(i=0;i<_nA;i++){
 470           temp+=_RA[i][ia]*_amp3[i][ib][ic];
 471         }
 472         _amp[ia][ib][ic]=temp;
 473       }
 474     }
 475   }
 476
 477
 478 }
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490