/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Generator/BesEvtGen/BesEvtGen-00-03-58/src/EvtGen/EvtGenBase/EvtEvalHelAmp.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 2000      Caltech, UCSB
//
// Module: EvtHelAmp.cc
//
// Description: Decay model for implementation of generic 2 body
//              decay specified by the partial wave amplitudes
//
//
// Modification history:
//
//    fkw        February 2, 2001     changes to satisfy KCC
//    RYD       September 7, 2000       Module created
//    Ping RG   Apr.  15,2007     change to helicity angle(@Line=210) for sequential decays
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtEvalHelAmp.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtdFunction.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtHelSys.hh"

using std::endl;


EvtEvalHelAmp::~EvtEvalHelAmp() {

  //allocate memory
  delete [] _lambdaA2;
  delete [] _lambdaB2;
  delete [] _lambdaC2;

  int ia,ib,ic;
  for(ib=0;ib<_nB;ib++){
    delete [] _HBC[ib];
  }

  delete [] _HBC;


  for(ia=0;ia<_nA;ia++){
    delete [] _RA[ia];
  }
  delete [] _RA;

  for(ib=0;ib<_nB;ib++){
    delete [] _RB[ib];
  }
  delete [] _RB;

  for(ic=0;ic<_nC;ic++){
    delete [] _RC[ic];
  }
  delete [] _RC;

 
  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      delete [] _amp[ia][ib];
      delete [] _amp1[ia][ib];
      delete [] _amp3[ia][ib];
    }
    delete [] _amp[ia];
    delete [] _amp1[ia];
    delete [] _amp3[ia];
  }

  delete [] _amp;
  delete [] _amp1;
  delete [] _amp3;

}


EvtEvalHelAmp::EvtEvalHelAmp(EvtSpinType::spintype typeA,
                             EvtSpinType::spintype typeB,
                             EvtSpinType::spintype typeC,
                             EvtComplexPtrPtr HBC){

  //find out how many states each particle have
  _nA=EvtSpinType::getSpinStates(typeA);
  _nB=EvtSpinType::getSpinStates(typeB);
  _nC=EvtSpinType::getSpinStates(typeC);

  //find out what 2 times the spin is
  _JA2=EvtSpinType::getSpin2(typeA);
  _JB2=EvtSpinType::getSpin2(typeB);
  _JC2=EvtSpinType::getSpin2(typeC);


  //allocate memory
  _lambdaA2=new int[_nA];
  _lambdaB2=new int[_nB];
  _lambdaC2=new int[_nC];

  _HBC=new EvtComplexPtr[_nB];
  int ia,ib,ic;
  for(ib=0;ib<_nB;ib++){
    _HBC[ib]=new EvtComplex[_nC];
  }


  _RA=new EvtComplexPtr[_nA];
  for(ia=0;ia<_nA;ia++){
    _RA[ia]=new EvtComplex[_nA];
  }
  _RB=new EvtComplexPtr[_nB];
  for(ib=0;ib<_nB;ib++){
    _RB[ib]=new EvtComplex[_nB];
  }
  _RC=new EvtComplexPtr[_nC];
  for(ic=0;ic<_nC;ic++){
    _RC[ic]=new EvtComplex[_nC];
  }
  
  _amp=new EvtComplexPtrPtr[_nA];
  _amp1=new EvtComplexPtrPtr[_nA];
  _amp3=new EvtComplexPtrPtr[_nA];
  for(ia=0;ia<_nA;ia++){
    _amp[ia]=new EvtComplexPtr[_nB];
    _amp1[ia]=new EvtComplexPtr[_nB];
    _amp3[ia]=new EvtComplexPtr[_nB];
    for(ib=0;ib<_nB;ib++){
      _amp[ia][ib]=new EvtComplex[_nC];
      _amp1[ia][ib]=new EvtComplex[_nC];
      _amp3[ia][ib]=new EvtComplex[_nC];
    }
  }

  //find the allowed helicities (actually 2*times the helicity!)

  fillHelicity(_lambdaA2,_nA,_JA2);
  fillHelicity(_lambdaB2,_nB,_JB2);
  fillHelicity(_lambdaC2,_nC,_JC2);

  for(ib=0;ib<_nB;ib++){
    for(ic=0;ic<_nC;ic++){
      _HBC[ib][ic]=HBC[ib][ic];
    }
  }
}






double EvtEvalHelAmp::probMax(){

  double c=1.0/sqrt(4*EvtConst::pi/(_JA2+1));

  int ia,ib,ic;


  double theta;
  int itheta;

  double maxprob=0.0;

  for(itheta=-10;itheta<=10;itheta++){
    theta=acos(0.099999*itheta);
    for(ia=0;ia<_nA;ia++){
      double prob=0.0;
      for(ib=0;ib<_nB;ib++){
        for(ic=0;ic<_nC;ic++){
          _amp[ia][ib][ic]=0.0;
          if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
            _amp[ia][ib][ic]=c*_HBC[ib][ic]*
              EvtdFunction::d(_JA2,_lambdaA2[ia],
                              _lambdaB2[ib]-_lambdaC2[ic],theta);
            prob+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
          }
        }
      }
      
      prob*=sqrt(1.0*_nA);
      
      if (prob>maxprob) maxprob=prob;

    }
  }

  return maxprob;

}


void EvtEvalHelAmp::evalAmp( EvtParticle *p, EvtAmp& amp){

  //find theta and phi of the first daughter
  
  EvtVector4R pB=p->getDaug(0)->getP4();
  EvtVector4R pC=p->getDaug(1)->getP4();
  EvtVector4R pP=pB+pC;

  EvtHelSys angles(pP,pB);
  double theta=angles.getHelAng(1);
  double phi  =angles.getHelAng(2);
//  double theta=acos(pB.get(3)/pB.d3mag());  //pingrg ,2007,4,15
//  double phi=atan2(pB.get(2),pB.get(1));

  double c=sqrt((_JA2+1)/(4*EvtConst::pi));

  int ia,ib,ic;

  double prob1=0.0;

  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      for(ic=0;ic<_nC;ic++){
        _amp[ia][ib][ic]=0.0;
        if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
          double dfun=EvtdFunction::d(_JA2,_lambdaA2[ia],
                                      _lambdaB2[ib]-_lambdaC2[ic],theta);

          _amp[ia][ib][ic]=c*_HBC[ib][ic]*
            exp(EvtComplex(0.0,phi*0.5*(_lambdaA2[ia]-_lambdaB2[ib]+
                                        _lambdaC2[ic])))*dfun;
        }
        prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
      }
    }
  }

  setUpRotationMatrices(p,theta,phi);

  applyRotationMatrices();

  double prob2=0.0;

  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      for(ic=0;ic<_nC;ic++){
        prob2+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
        if (_nA==1){
          if (_nB==1){
            if (_nC==1){
              amp.vertex(_amp[ia][ib][ic]);
            }
            else{
              amp.vertex(ic,_amp[ia][ib][ic]);
            }
          }
          else{
            if (_nC==1){
              amp.vertex(ib,_amp[ia][ib][ic]);
            }
            else{
              amp.vertex(ib,ic,_amp[ia][ib][ic]);
              //report(INFO,"EvtGen") << "ib,ic:"<<ib<<" "<<ic<<" "<<_amp[ia][ib][ic]<<endl;
            }
          }
        }else{
          if (_nB==1){
            if (_nC==1){
              amp.vertex(ia,_amp[ia][ib][ic]);
            }
            else{
              amp.vertex(ia,ic,_amp[ia][ib][ic]);
            }
          }
          else{
            if (_nC==1){
              amp.vertex(ia,ib,_amp[ia][ib][ic]);
            }
            else{
              amp.vertex(ia,ib,ic,_amp[ia][ib][ic]);
            }
          }
        }
      }
    }
  }

  if (fabs(prob1-prob2)>0.000001*prob1){
    report(INFO,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl;
    ::abort();
  }

  return ;

}


void EvtEvalHelAmp::fillHelicity(int* lambda2,int n,int J2){
  
  int i;
  
  //photon is special case!
  if (n==2&&J2==2) {
    lambda2[0]=2;
    lambda2[1]=-2;
    return;
  }
  
  assert(n==J2+1);

  for(i=0;i<n;i++){
    lambda2[i]=n-i*2-1;
  }

  return;

}


void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){

  switch(_JA2){

  case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:

    {

      EvtSpinDensity R=p->rotateToHelicityBasis();

      
      int i,j,n;
      
      n=R.GetDim();
      
      assert(n==_nA);
        
      
      for(i=0;i<n;i++){
        for(j=0;j<n;j++){
          _RA[i][j]=R.Get(i,j);
        }
      }

    }

    break;

  default:
    report(ERROR,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl;
    ::abort();
  }
  
  
  switch(_JB2){


  case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:

    {
      
      int i,j,n;

      EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi);
      
      n=R.GetDim();
      
      assert(n==_nB);
        
      //report(INFO,"EvtGen") << "RB:"<< R <<endl; 
     
      for(i=0;i<n;i++){
        for(j=0;j<n;j++){
          _RB[i][j]=conj(R.Get(i,j));
        }
      }

    }

    break;

  default:
    report(ERROR,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl;
    ::abort();
  }
  
  switch(_JC2){

  case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:

    {

      int i,j,n;

      EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi);
            
      n=R.GetDim();

      //report(INFO,"EvtGen") << "RC:"<< R <<endl;

      assert(n==_nC);

      for(i=0;i<n;i++){
        for(j=0;j<n;j++){
          _RC[i][j]=conj(R.Get(i,j));
        }
      }

    }

    break;

  default:
    report(ERROR,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl;
    ::abort();
  }
  
  

}


void EvtEvalHelAmp::applyRotationMatrices(){

  int ia,ib,ic,i;
  
  EvtComplex temp;



  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      for(ic=0;ic<_nC;ic++){
        temp=0;
        for(i=0;i<_nC;i++){
          temp+=_RC[i][ic]*_amp[ia][ib][i];
        }
        _amp1[ia][ib][ic]=temp;
      }
    }
  }



  for(ia=0;ia<_nA;ia++){
    for(ic=0;ic<_nC;ic++){
      for(ib=0;ib<_nB;ib++){
        temp=0;
        for(i=0;i<_nB;i++){
          temp+=_RB[i][ib]*_amp1[ia][i][ic];
        }
        _amp3[ia][ib][ic]=temp;
      }
    }
  }
  


  for(ib=0;ib<_nB;ib++){
    for(ic=0;ic<_nC;ic++){
      for(ia=0;ia<_nA;ia++){
        temp=0;
        for(i=0;i<_nA;i++){
          temp+=_RA[i][ia]*_amp3[i][ib][ic];
        }
        _amp[ia][ib][ic]=temp;
      }
    }
  }


}

---