Page MenuHomeHEPForge
Files F10881573

No OneTemporary
Actions

View Options

diff --git a/src/EvtGenBase/EvtAmp.cpp b/src/EvtGenBase/EvtAmp.cpp
index 2d65b99..95e0e5f 100644
--- a/src/EvtGenBase/EvtAmp.cpp
+++ b/src/EvtGenBase/EvtAmp.cpp
@@ -1,551 +1,551 @@
//--------------------------------------------------------------------------
//
// 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) 1998 Caltech, UCSB
//
// Module: EvtAmp.cc
//
// Description: Class to manipulate the amplitudes in the decays.
//
// Modification history:
//
// RYD May 29, 1997 Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include <math.h>
#include <assert.h>
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtSpinDensity.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
using std::endl;
EvtAmp::EvtAmp(){
_ndaug=0;
_pstates=0;
_nontrivial=0;
}
EvtAmp::EvtAmp(const EvtAmp& amp){
int i;
_ndaug=amp._ndaug;
_pstates=amp._pstates;
for(i=0;i<_ndaug;i++){
dstates[i]=amp.dstates[i];
_dnontrivial[i]=amp._dnontrivial[i];
}
_nontrivial=amp._nontrivial;
int namp=1;
for(i=0;i<_nontrivial;i++){
_nstate[i]=amp._nstate[i];
namp*=_nstate[i];
}
for(i=0;i<namp;i++){
assert(i<125);
_amp[i]=amp._amp[i];
}
}
void EvtAmp::init(EvtId p,int ndaugs,EvtId *daug){
setNDaug(ndaugs);
int ichild;
int daug_states[100],parstates;
for(ichild=0;ichild<ndaugs;ichild++){
daug_states[ichild]=
EvtSpinType::getSpinStates(EvtPDL::getSpinType(daug[ichild]));
}
parstates=EvtSpinType::getSpinStates(EvtPDL::getSpinType(p));
setNState(parstates,daug_states);
}
void EvtAmp::setNDaug(int n){
_ndaug=n;
}
void EvtAmp::setNState(int parent_states,int *daug_states){
_nontrivial=0;
_pstates=parent_states;
if(_pstates>1) {
_nstate[_nontrivial]=_pstates;
_nontrivial++;
}
int i;
for(i=0;i<_ndaug;i++){
dstates[i]=daug_states[i];
_dnontrivial[i]=-1;
if(daug_states[i]>1) {
_nstate[_nontrivial]=daug_states[i];
_dnontrivial[i]=_nontrivial;
_nontrivial++;
}
}
if (_nontrivial>5) {
report(ERROR,"EvtGen") << "Too many nontrivial states in EvtAmp!"<<endl;
}
}
void EvtAmp::setAmp(int *ind, const EvtComplex& a){
int nstatepad = 1;
int position = ind[0];
for ( int i=1; i<_nontrivial; i++ ) {
nstatepad *= _nstate[i-1];
position += nstatepad*ind[i];
}
assert(position<125);
_amp[position] = a;
}
const EvtComplex& EvtAmp::getAmp(int *ind)const{
int nstatepad = 1;
int position = ind[0];
for ( int i=1; i<_nontrivial; i++ ) {
nstatepad *= _nstate[i-1];
position += nstatepad*ind[i];
}
return _amp[position];
}
EvtSpinDensity EvtAmp::getSpinDensity(){
EvtSpinDensity rho;
rho.setDim(_pstates);
EvtComplex temp;
int i,j,n;
if (_pstates==1) {
if (_nontrivial==0) {
rho.set(0,0,_amp[0]*conj(_amp[0]));
return rho;
}
n=1;
temp = EvtComplex(0.0);
for(i=0;i<_nontrivial;i++){
n*=_nstate[i];
}
for(i=0;i<n;i++){
temp+=_amp[i]*conj(_amp[i]);
}
rho.set(0,0,temp);;
return rho;
}
else{
for(i=0;i<_pstates;i++){
for(j=0;j<_pstates;j++){
temp = EvtComplex(0.0);
int kk;
int allloop = 1;
- for (kk=0;kk<(_nontrivial-1); kk++ ) {
+ for (kk=0;kk<_ndaug; kk++ ) {
allloop *= dstates[kk];
}
for (kk=0; kk<allloop; kk++) {
temp += _amp[_pstates*kk+i]*conj(_amp[_pstates*kk+j]);}
// if (_nontrivial>3){
//report(ERROR,"EvtGen") << "Can't handle so many states in EvtAmp!"<<endl;
//}
rho.set(i,j,temp);
}
}
return rho;
}
}
EvtSpinDensity EvtAmp::getBackwardSpinDensity(EvtSpinDensity *rho_list){
EvtSpinDensity rho;
rho.setDim(_pstates);
if (_pstates==1){
rho.set(0,0,EvtComplex(1.0,0.0));
return rho;
}
int k;
EvtAmp ampprime;
ampprime=(*this);
for(k=0;k<_ndaug;k++){
if (dstates[k]!=1){
ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
}
}
return ampprime.contract(0,(*this));
}
EvtSpinDensity EvtAmp::getForwardSpinDensity(EvtSpinDensity *rho_list,int i){
EvtSpinDensity rho;
rho.setDim(dstates[i]);
int k;
if (dstates[i]==1) {
rho.set(0,0,EvtComplex(1.0,0.0));
return rho;
}
EvtAmp ampprime;
ampprime=(*this);
if (_pstates!=1){
ampprime=ampprime.contract(0,rho_list[0]);
}
for(k=0;k<i;k++){
if (dstates[k]!=1){
ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
}
}
return ampprime.contract(_dnontrivial[i],(*this));
}
EvtAmp EvtAmp::contract(int k,const EvtSpinDensity& rho){
EvtAmp temp;
int i,j;
temp._ndaug=_ndaug;
temp._pstates=_pstates;
temp._nontrivial=_nontrivial;
for(i=0;i<_ndaug;i++){
temp.dstates[i]=dstates[i];
temp._dnontrivial[i]=_dnontrivial[i];
}
if (_nontrivial==0) {
report(ERROR,"EvtGen")<<"Should not be here EvtAmp!"<<endl;
}
for(i=0;i<_nontrivial;i++){
temp._nstate[i]=_nstate[i];
}
EvtComplex c;
int index[10];
for (i=0;i<10;i++) {
index[i] = 0;
}
int allloop = 1;
int indflag,ii;
for (i=0;i<_nontrivial;i++) {
allloop *= _nstate[i];
}
for ( i=0;i<allloop;i++) {
c = EvtComplex(0.0);
int tempint = index[k];
for (j=0;j<_nstate[k];j++) {
index[k] = j;
c+=rho.get(j,tempint)*getAmp(index);
}
index[k] = tempint;
temp.setAmp(index,c);
indflag = 0;
for ( ii=0;ii<_nontrivial;ii++) {
if ( indflag == 0 ) {
if ( index[ii] == (_nstate[ii]-1) ) {
index[ii] = 0;
}
else {
indflag = 1;
index[ii] += 1;
}
}
}
}
return temp;
}
EvtSpinDensity EvtAmp::contract(int k,const EvtAmp& amp2){
int i,j,l;
EvtComplex temp;
EvtSpinDensity rho;
rho.setDim(_nstate[k]);
int allloop = 1;
int indflag,ii;
for (i=0;i<_nontrivial;i++) {
allloop *= _nstate[i];
}
int index[10];
int index1[10];
// int l;
for(i=0;i<_nstate[k];i++){
for(j=0;j<_nstate[k];j++){
if (_nontrivial==0) {
report(ERROR,"EvtGen")<<"Should not be here1 EvtAmp!"<<endl;
rho.set(0,0,EvtComplex(1.0,0.0));
return rho;
}
for (ii=0;ii<10;ii++) {
index[ii] = 0;
index1[ii] = 0;
}
index[k] = i;
index1[k] = j;
temp = EvtComplex(0.0);
for ( l=0;l<int(allloop/_nstate[k]);l++) {
temp+=getAmp(index)*conj(amp2.getAmp(index1));
indflag = 0;
for ( ii=0;ii<_nontrivial;ii++) {
if ( ii!= k) {
if ( indflag == 0 ) {
if ( index[ii] == (_nstate[ii]-1) ) {
index[ii] = 0;
index1[ii] = 0;
}
else {
indflag = 1;
index[ii] += 1;
index1[ii] += 1;
}
}
}
}
}
rho.set(i,j,temp);
}
}
return rho;
}
EvtAmp EvtAmp::contract(int i, const EvtAmp& a1,const EvtAmp& a2){
//Do we need this method?
assert(a2._pstates>1&&a2._nontrivial==1);
assert(i<=a1._nontrivial);
EvtAmp tmp;
report(DEBUG,"EvtGen") << "EvtAmp::contract not written yet" << endl;
return tmp;
}
void EvtAmp::dump(){
int i,list[10];
report(DEBUG,"EvtGen") << "Number of daugthers:"<<_ndaug<<endl;
report(DEBUG,"EvtGen") << "Number of states of the parent:"<<_pstates<<endl;
report(DEBUG,"EvtGen") << "Number of states on daughters:";
for (i=0;i<_ndaug;i++){
report(DEBUG,"EvtGen") <<dstates[i]<<" ";
}
report(DEBUG,"EvtGen") << endl;
report(DEBUG,"EvtGen") << "Nontrivial index of daughters:";
for (i=0;i<_ndaug;i++){
report(DEBUG,"EvtGen") <<_dnontrivial[i]<<" ";
}
report(DEBUG,"EvtGen") <<endl;
report(DEBUG,"EvtGen") <<"number of nontrivial states:"<<_nontrivial<<endl;
report(DEBUG,"EvtGen") << "Nontrivial particles number of states:";
for (i=0;i<_nontrivial;i++){
report(DEBUG,"EvtGen") <<_nstate[i]<<" ";
}
report(DEBUG,"EvtGen") <<endl;
report(DEBUG,"EvtGen") <<"Amplitudes:"<<endl;
if (_nontrivial==0){
list[0] = 0;
report(DEBUG,"EvtGen") << getAmp(list) << endl;
}
int allloop[10];
allloop[0]=1;
for (i=0;i<_nontrivial;i++) {
if (i==0){
allloop[i] *= _nstate[i];
}
else{
allloop[i] = allloop[i-1]*_nstate[i];
}
}
int index = 0;
for (i=0;i<allloop[_nontrivial-1];i++) {
report(DEBUG,"EvtGen") << getAmp(list) << " ";
if ( i==allloop[index]-1 ) {
index ++;
report(DEBUG,"EvtGen") << endl;
}
}
report(DEBUG,"EvtGen") << "-----------------------------------"<<endl;
}
void EvtAmp::vertex(const EvtComplex& c){
int list[1];
list[0] = 0;
setAmp(list,c);
}
void EvtAmp::vertex(int i,const EvtComplex& c){
int list[1];
list[0] = i;
setAmp(list,c);
}
void EvtAmp::vertex(int i,int j,const EvtComplex& c){
int list[2];
list[0] = i;
list[1] = j;
setAmp(list,c);
}
void EvtAmp::vertex(int i,int j,int k,const EvtComplex& c){
int list[3];
list[0] = i;
list[1] = j;
list[2] = k;
setAmp(list,c);
}
void EvtAmp::vertex(int *i1,const EvtComplex& c){
setAmp(i1,c);
}
EvtAmp& EvtAmp::operator=(const EvtAmp& amp){
int i;
_ndaug=amp._ndaug;
_pstates=amp._pstates;
for(i=0;i<_ndaug;i++){
dstates[i]=amp.dstates[i];
_dnontrivial[i]=amp._dnontrivial[i];
}
_nontrivial=amp._nontrivial;
int namp=1;
for(i=0;i<_nontrivial;i++){
_nstate[i]=amp._nstate[i];
namp*=_nstate[i];
}
for(i=0;i<namp;i++){
assert(i<125);
_amp[i]=amp._amp[i];
}
return *this;
}
diff --git a/src/EvtGenBase/EvtSpinDensity.cpp b/src/EvtGenBase/EvtSpinDensity.cpp
index ec60482..c427613 100644
--- a/src/EvtGenBase/EvtSpinDensity.cpp
+++ b/src/EvtGenBase/EvtSpinDensity.cpp
@@ -1,206 +1,214 @@
//--------------------------------------------------------------------------
//
// 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) 1998 Caltech, UCSB
//
// Module: EvtSpinDensity.cc
//
// Description: Class to reperesent spindensity matrices.
//
// Modification history:
//
// RYD May 29,1997 Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <assert.h>
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtSpinDensity.hh"
#include "EvtGenBase/EvtReport.hh"
using std::endl;
using std::ostream;
EvtSpinDensity::EvtSpinDensity(const EvtSpinDensity& density){
dim=0;
rho=0;
int i,j;
setDim(density.dim);
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
rho[i][j]=density.rho[i][j];
}
}
}
EvtSpinDensity& EvtSpinDensity::operator=(const EvtSpinDensity& density){
int i,j;
setDim(density.dim);
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
rho[i][j]=density.rho[i][j];
}
}
return *this;
}
EvtSpinDensity::~EvtSpinDensity(){
if (dim!=0){
int i;
for(i=0;i<dim;i++) delete [] rho[i];
}
delete [] rho;
}
EvtSpinDensity::EvtSpinDensity(){
dim=0;
rho=0;
}
void EvtSpinDensity::setDim(int n){
if (dim==n) return;
if (dim!=0){
int i;
for(i=0;i<dim;i++) delete [] rho[i];
delete [] rho;
rho=0;
dim=0;
}
if (n==0) return;
dim=n;
rho=new EvtComplexPtr[n];
int i;
for(i=0;i<n;i++){
rho[i]=new EvtComplex[n];
}
}
int EvtSpinDensity::getDim() const {
return dim;
}
void EvtSpinDensity::set(int i,int j,const EvtComplex& rhoij){
assert(i<dim&&j<dim);
rho[i][j]=rhoij;
}
const EvtComplex& EvtSpinDensity::get(int i,int j)const{
assert(i<dim&&j<dim);
return rho[i][j];
}
void EvtSpinDensity::setDiag(int n){
setDim(n);
int i,j;
for(i=0;i<n;i++){
for(j=0;j<n;j++){
rho[i][j]=EvtComplex(0.0);
}
rho[i][i]=EvtComplex(1.0);
}
}
double EvtSpinDensity::normalizedProb(const EvtSpinDensity& d){
int i,j;
EvtComplex prob(0.0,0.0);
double norm=0.0;
if (dim!=d.dim) {
report(ERROR,"EvtGen")<<"Not matching dimensions in NormalizedProb"<<endl;
::abort();
}
for(i=0;i<dim;i++){
norm+=real(rho[i][i]);
for(j=0;j<dim;j++){
prob+=rho[i][j]*d.rho[i][j];
}
}
if (imag(prob)>0.00000001*real(prob)) {
report(ERROR,"EvtGen")<<"Imaginary probability:"<<prob<<" "<<norm<<endl;
}
if (real(prob)<0.0) {
report(ERROR,"EvtGen")<<"Negative probability:"<<prob<<" "<<norm<<endl;
}
return real(prob)/norm;
}
int EvtSpinDensity::check(){
if (dim<1) {
report(ERROR,"EvtGen")<<"dim="<<dim<<"in SpinDensity::Check"<<endl;
}
int i,j;
+ double trace(0.0);
+
+ for (i=0;i<dim;i++) {
+ trace += abs(rho[i][i]);
+ }
+
for(i=0;i<dim;i++){
if (real(rho[i][i])<0.0) return 0;
- if (imag(rho[i][i])*1000000.0>abs(rho[i][i])) {
+ if (imag(rho[i][i])*1000000.0>trace) {
+ report(INFO,"EvtGen") << *this << endl;
+ report(INFO,"EvtGen") << trace << endl;
report(INFO,"EvtGen") << "Failing 1"<<endl;
return 0;
}
}
for(i=0;i<dim;i++){
for(j=i+1;j<dim;j++){
if (fabs(real(rho[i][j]-rho[j][i]))>
0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
report(INFO,"EvtGen") << "Failing 2"<<endl;
return 0;
}
if (fabs(imag(rho[i][j]+rho[j][i]))>
0.00000001*(abs(rho[i][i])+abs(rho[j][j]))) {
report(INFO,"EvtGen") << "Failing 3"<<endl;
return 0;
}
}
}
return 1;
}
ostream& operator<<(ostream& s,const EvtSpinDensity& d){
int i,j;
s << endl;
s << "Dimension:"<<d.dim<<endl;
for (i=0;i<d.dim;i++){
for (j=0;j<d.dim;j++){
s << d.rho[i][j]<<" ";
}
s <<endl;
}
return s;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:33 AM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983071
Default Alt Text
(14 KB)

Event Timeline