Skip to content
Snippets Groups Projects
Commit db48fee6 authored by Bertram Kopf's avatar Bertram Kopf
Browse files

more generalized tensor amps

parent 4a7f7a67
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Examples/Tutorial/qft++Exams/TensorApp.cc
+ 8
− 1
View file @ db48fee6
......@@ -65,7 +65,7 @@ int main(int __argc,char *__argv[]){
// Tensor is a template class, here we'll work with double's and
// complex<double>'s...but you could use any data type that defines the
// proper operators. The constructor argument is the rank.
Tensor<double> x1(1),x2(2),x3(3);
Tensor<double> x1(1),x2(2),x3(3),x4(4);
Tensor<complex<double> > y1(1);
Vector4<double> v;
LeviCivitaTensor eps; // totally anti-symetric 4th rank tensor
......@@ -92,6 +92,9 @@ int main(int __argc,char *__argv[]){
x3 = x1 % x2; // tensor outer product
cout << "->x3: x1^{mu}x2^{nu,rho} (printing only defined for rank <= 2)"
<< endl;
x4 = x1 % x3;
cout << "->x4: x1^{mu}x3^{nu,rho,pi} (printing only defined for rank <= 2)"
<< endl;
// just make sure these don't throw errors
x3.Symmetric();
x3.AntiSymmetric();
......@@ -133,6 +136,10 @@ int main(int __argc,char *__argv[]){
cout << "(x1^{mu}x2^{nu,rho} + 2*x3^{mu,nu,rho})x3_{mu,nu,rho}:\n(((x1%x2) + 2*x3)|x3)\t->\t"
<< (((x1%x2) + 2*x3)|x3) << endl;
cout << "epsilon^{mu,nu,rho,pi} x2_{mu,nu}:\n(eps|x2)\t->\t" << (eps|x2) << endl;
cout << "x3(mu,nu,rho)x4(pi,mu,nu,rho):\nx4.Contract(x3, 3)->\t" << x4.Contract(x3, 3) << endl;
cout << "x3(mu,nu,rho)x4(pi,mu,nu,rho):\nx4.Permute(1,x4.Rank())\t x3|x4->\t" << (x4|x3) << endl;
cout << "x3(pi,mu,nu)x4(pi,mu,nu,rho):\nx3.Contract(x4, 3)->\t" << x3.Contract(x4, 3) << endl;
cout << "x3(pi,mu,nu)x4(pi,mu,nu,rho):\nx(x3|x4)->\t" << (x3|x4) << endl;
cout << "etc...as complicated as you want, it's still easy in the code."
<< endl;
......
PwaUtils/EvtDataBaseList.cc
+ 2
− 1
View file @ db48fee6
......@@ -76,7 +76,8 @@ void EvtDataBaseList::read4Vecs(EventList& evtList, std::vector<EvtData*>& theEv
int evtCount = 0;
while ((anEvent = evtList.nextEvent())){
if (evtCount>= maxEvts) break;
if (evtCount%10000 == 0) Info << "4vec calculation for event " << evtCount ; // << endmsg;
// if (evtCount%10000 == 0) Info << "4vec calculation for event " << evtCount ; // << endmsg;
if (evtCount%500 == 0) Info << "4vec calculation for event " << evtCount ; // << endmsg;
Vector4<double> V4_all_lab(0.,0.,0.,0.);
......
PwaUtils/IsobarTensorDecay.cc
+ 96
− 53
View file @ db48fee6
......@@ -54,6 +54,7 @@ IsobarTensorDecay::IsobarTensorDecay(boost::shared_ptr<const jpcRes> motherJPCPt
,_polDaughter1(PolVector(_daughter1JPCPtr->J))
,_polDaughter2(PolVector(_daughter2JPCPtr->J))
,_lctTensor(LeviCivitaTensor())
,_metricTensor(MetricTensor())
{
}
......@@ -137,26 +138,61 @@ void IsobarTensorDecay::fillWignerDs(std::map<std::string, Vector4<double> >& fs
DebugMsg << name() << endmsg;
std::vector< boost::shared_ptr<const JPCLS> >::iterator itJPCLS;
for(itJPCLS=theJPCLSAmps.begin(); itJPCLS!=theJPCLSAmps.end(); ++itJPCLS){
(*itJPCLS)->print(std::cout);
std::cout << std::endl;
// (*itJPCLS)->print(std::cout);
// std::cout << std::endl;
Spin L=(*itJPCLS)->L;
Spin S=(*itJPCLS)->S;
int s1s2S=spinDaughter1+spinDaughter2+(*itJPCLS)->S;
Tensor<complex<double> > leviPssTensor;
bool add_lctForChi=true;
if( s1s2S%2 ==0 ) add_lctForChi=false; //odd
if( s1s2S%2 ==0 ) { //odd
add_lctForChi=false;
leviPssTensor(0)=complex<double>(1.,0.);
}
else{
leviPssTensor=_lctTensor*mother_4Vec;
}
DebugMsg << "leviPssTensor: " << leviPssTensor << endmsg;
int SLms1=(*itJPCLS)->S+(*itJPCLS)->L+_motherJPCPtr->J;
Tensor<complex<double> > leviPlsTensor;
bool add_lctForTensor=true;
if( SLms1%2 ==0 ) add_lctForTensor=false; //odd
if( SLms1%2 ==0 ){ //odd
add_lctForTensor=false;
leviPlsTensor=complex<double>(1.,0.);
}
else{
leviPlsTensor=_lctTensor*mother_4Vec;
}
DebugMsg << "leviPlsTensor: " << leviPlsTensor << endmsg;
OrbitalTensor orbTensor(L);
orbTensor.SetP4(daughter1_4Vec, daughter2_4Vec);
DebugMsg << "orbTensor:\t" << orbTensor << endmsg;
for (Spin lamMother=-lamMotherMax; lamMother<=lamMotherMax; ++lamMother){
DebugMsg << "lamMother:\t" << lamMother << endmsg;
Tensor<complex<double> > epsilonMotherProject = _polMother(lamMother);
DebugMsg << "epsilonMotherProject:\t" << epsilonMotherProject << endmsg;
// //calculate ls part;
// Tensor<complex<double> > lsPartTensor;
// if(add_lctForTensor){
// if(epsilonMotherProject.Rank() >= orbTensor.Rank()){
// lsPartTensor= epsilonMotherProject*leviPlsTensor;
// lsPartTensor.Permute(1,lsPartTensor.Rank());
// lsPartTensor=lsPartTensor|orbTensor;
// }
// else{
// lsPartTensor= orbTensor*leviPlsTensor;
// lsPartTensor.Permute(1,lsPartTensor.Rank());
// lsPartTensor=lsPartTensor|epsilonMotherProject;
// }
// }
// for !add_lctForTensor calculation will be done in heli loops
//calculate chi
Tensor<complex<double> > chi12;
......@@ -164,66 +200,73 @@ void IsobarTensorDecay::fillWignerDs(std::map<std::string, Vector4<double> >& fs
PolVector part12PolVec(S);
part12PolVec.SetP4(mother_4Vec,mother_4Vec.M());
s12SpinProjector=part12PolVec.Projector();
DebugMsg << "s12SpinProjector(0,0,0,0):\t" << s12SpinProjector(0,0,0,0) << endmsg;
DebugMsg << "s12SpinProjector:\t" << s12SpinProjector << endmsg;
for (Spin lamDaughter1=-spinDaughter1; lamDaughter1 <= spinDaughter1; ++lamDaughter1){
DebugMsg << "lamDaughter1:\t" << lamDaughter1 << endmsg;
Tensor<complex<double> > epsilonDaughter1Project = _polDaughter1(lamDaughter1);
DebugMsg << "epsilonDaughter1Project:\t" << epsilonDaughter1Project << endmsg;
for (Spin lamDaughter2=-spinDaughter2; lamDaughter2<=spinDaughter2; ++lamDaughter2){
DebugMsg << "lamDaughter2:\t" << lamDaughter2 << endmsg;
Tensor<complex<double> > epsilonDaughter2Project=_polDaughter2(lamDaughter2);
DebugMsg << "epsilonDaughter2Project:\t" << epsilonDaughter2Project << endmsg;
if(spinDaughter2==0){
chi12= s12SpinProjector | conj(epsilonDaughter1Project);
// DebugMsg << "chi12 Tensor Amp (spinDaughter2==0):\t" << chi12 << endmsg;
}
Tensor<complex<double> > chiPart;
if(add_lctForChi){
chiPart=(conj(epsilonDaughter1Project)*leviPssTensor)*conj(epsilonDaughter2Project);
}
else{
Tensor<complex<double> > epsilonDaughter2Project=_polDaughter2(lamDaughter2);
if(spinDaughter1==0){
chi12 = s12SpinProjector | conj(epsilonDaughter2Project);
// DebugMsg << "chi12 Tensor Amp (spinDaughter1==0):\t" << chi12 << endmsg;
}
else{
if(add_lctForChi) chi12=s12SpinProjector | _lctTensor | (conj(epsilonDaughter1Project)%conj(epsilonDaughter2Project));
else chi12=s12SpinProjector | (conj(epsilonDaughter1Project)%conj(epsilonDaughter2Project));
}
}
if( S!=0 && (spinDaughter1!=0 || spinDaughter2!=0) ){
DebugMsg << "chi12 Tensor Amp:\t" << chi12 << endmsg;
Tensor<complex<double> > result;
if(add_lctForTensor){
// result = epsilonMotherProject | ( ( (_lctTensor * mother_4Vec ) | orbTensor) | chi12);
result = epsilonMotherProject | ( ( (_lctTensor * mother_4Vec ) * orbTensor) | chi12);
// result = (epsilonMotherProject | ( (_lctTensor * mother_4Vec ) | orbTensor) ) | chi12;
// result = epsilonMotherProject | ( (_lctTensor * mother_4Vec ) | (orbTensor * chi12) );
// result = epsilonMotherProject | ( (_lctTensor * mother_4Vec ) % (orbTensor | chi12) );
// result = epsilonMotherProject | ( (_lctTensor * mother_4Vec ) | (orbTensor | chi12));
// result = (epsilonMotherProject % mother_4Vec) | (orbTensor % chi12);
}
else{
if( orbTensor.Rank()+chi12.Rank() == epsilonMotherProject.Rank()) result = epsilonMotherProject | (orbTensor % chi12);
else if( orbTensor.Rank() == chi12.Rank() + epsilonMotherProject.Rank()) result = epsilonMotherProject | (orbTensor | chi12);
}
DebugMsg << "result Tensor Amp for S!=0 && (spinDaughter1!=0 || spinDaughter2!=0 " << _name << ":\t" << result << endmsg;
DebugMsg << "epsilonMotherProject " << epsilonMotherProject << endmsg;
DebugMsg << "orbTensor " << orbTensor << endmsg;
DebugMsg << "chi12 " << chi12 << endmsg;
evtData->ComplexDouble5SpinString[_name][L][S][lamMother][lamDaughter1][lamDaughter2]=result(0);
}
if(spinMother==0 && spinDaughter1==0 && spinDaughter2==0) evtData->ComplexDouble5SpinString[_name][L][S][lamMother][lamDaughter1][lamDaughter2]= complex<double> (1.,0.);
else if(S==0 && spinDaughter1==0 && spinDaughter2==0) {
Tensor<complex<double> > result = epsilonMotherProject | orbTensor;
DebugMsg << "result Tensor Amp for " << _name << ":\t" << result << endmsg;
evtData->ComplexDouble5SpinString[_name][L][S][lamMother][lamDaughter1][lamDaughter2]=result(0);
chiPart=conj(epsilonDaughter1Project)%conj(epsilonDaughter2Project);
}
DebugMsg << "chiPart:\t" << chiPart << endmsg;
chi12=s12SpinProjector|chiPart;
DebugMsg << "chi12:\t" << chi12 << endmsg;
//calculate ls part;
Tensor<complex<double> > lsPartTensor;
if(add_lctForTensor){
if(epsilonMotherProject.Rank() >= orbTensor.Rank()){
lsPartTensor= epsilonMotherProject*leviPlsTensor;
lsPartTensor.Permute(1,lsPartTensor.Rank());
lsPartTensor=lsPartTensor|orbTensor;
}
else{
lsPartTensor= orbTensor*leviPlsTensor;
lsPartTensor.Permute(1,lsPartTensor.Rank());
lsPartTensor=lsPartTensor|epsilonMotherProject;
}
}
if(!add_lctForTensor){
if (chi12.Rank()==fabs(epsilonMotherProject.Rank()-orbTensor.Rank())) lsPartTensor= epsilonMotherProject|orbTensor;
else if (chi12.Rank()==(epsilonMotherProject.Rank()+orbTensor.Rank())) lsPartTensor= epsilonMotherProject%orbTensor;
else if (epsilonMotherProject.Rank()==orbTensor.Rank()){
int noOfContractions=int((epsilonMotherProject.Rank()+orbTensor.Rank()-chi12.Rank())/2);
if (noOfContractions<=epsilonMotherProject.Rank() && noOfContractions>0) lsPartTensor= epsilonMotherProject.Contract(orbTensor, noOfContractions);
}
else{
Alert << "wrong ranks chi12.Rank()=" << chi12.Rank()
<< "\tepsilonMotherProject.Rank()=" << epsilonMotherProject.Rank()
<< "\torbTensor.Rank()=" << orbTensor.Rank()
<< endmsg;
exit(0);
}
}
DebugMsg << "lsPartTensor:\t" << lsPartTensor << endmsg;
Tensor<complex<double> > result=lsPartTensor|chi12;
DebugMsg << "result:\t" << result << endmsg;
if(result.Rank()>0){
Alert << "result.Rank()=" << result.Rank() << " > 0" << endmsg;
exit(0);
}
evtData->ComplexDouble5SpinString[_name][L][S][lamMother][lamDaughter1][lamDaughter2]=result(0);
}
}
}
......
qft++/tensor/Tensor.tcc
+ 1
− 1
View file @ db48fee6
......@@ -303,7 +303,7 @@ void Tensor<_Tp>::Print(std::ostream& __os) const {
}
else{
cout << "<Tensor::Print(ostream&)> Error! Can NOT print a Tensor with "
<< " Rank > 2." << endl;
<< " Rank "<< _rank << " > 2." << endl;
}
}
//_____________________________________________________________________________
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment