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

added tensor amps; so far only working for spin0 decay particles

parent df48eafb
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 812 additions and 5 deletions
Examples/pbarp/pbarpReactionApp.cc
+ 2
− 0
View file @ fe4f5432
...@@ -52,6 +52,7 @@ ...@@ -52,6 +52,7 @@
#include "pbarpUtils/pbarpBaseLh.hh" #include "pbarpUtils/pbarpBaseLh.hh"
#include "pbarpUtils/pbarpHeliLh.hh" #include "pbarpUtils/pbarpHeliLh.hh"
#include "pbarpUtils/pbarpCanoLh.hh" #include "pbarpUtils/pbarpCanoLh.hh"
#include "pbarpUtils/pbarpTensorLh.hh"
#include "Event/EventReaderDefault.hh" #include "Event/EventReaderDefault.hh"
...@@ -106,6 +107,7 @@ int main(int __argc,char *__argv[]){ ...@@ -106,6 +107,7 @@ int main(int __argc,char *__argv[]){
boost::shared_ptr<AbsLh> theLhPtr; boost::shared_ptr<AbsLh> theLhPtr;
if(prodFormalism=="Cano") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpCanoLh()); if(prodFormalism=="Cano") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpCanoLh());
else if(prodFormalism=="Heli") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpHeliLh()); else if(prodFormalism=="Heli") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpHeliLh());
else if(prodFormalism=="Tensor") theLhPtr=boost::shared_ptr<AbsLh>(new pbarpTensorLh());
else { else {
Alert << "prodFormalism\t" << prodFormalism << "\tdoesn't exist!!!" << endmsg; Alert << "prodFormalism\t" << prodFormalism << "\tdoesn't exist!!!" << endmsg;
exit(1); exit(1);
......
PwaUtils/AbsEnv.cc
+ 2
− 0
View file @ fe4f5432
...@@ -33,6 +33,7 @@ ...@@ -33,6 +33,7 @@
#include "PwaUtils/AbsDecayList.hh" #include "PwaUtils/AbsDecayList.hh"
#include "PwaUtils/IsobarLSDecay.hh" #include "PwaUtils/IsobarLSDecay.hh"
#include "PwaUtils/IsobarHeliDecay.hh" #include "PwaUtils/IsobarHeliDecay.hh"
#include "PwaUtils/IsobarTensorDecay.hh"
#include "PwaUtils/OmegaTo3PiLSDecay.hh" #include "PwaUtils/OmegaTo3PiLSDecay.hh"
#include "PwaUtils/OmegaTo3PiTensorDecay.hh" #include "PwaUtils/OmegaTo3PiTensorDecay.hh"
#include "PwaUtils/ParserBase.hh" #include "PwaUtils/ParserBase.hh"
...@@ -146,6 +147,7 @@ void AbsEnv::setup(ParserBase* theParser){ ...@@ -146,6 +147,7 @@ void AbsEnv::setup(ParserBase* theParser){
if(daughterParticles.size()==2){ if(daughterParticles.size()==2){
if (usedSystem=="Heli") tmpDec= boost::shared_ptr<AbsDecay>(new IsobarHeliDecay(motherParticle, daughterParticles[0], daughterParticles[1], this)); if (usedSystem=="Heli") tmpDec= boost::shared_ptr<AbsDecay>(new IsobarHeliDecay(motherParticle, daughterParticles[0], daughterParticles[1], this));
else if (usedSystem=="Cano") tmpDec= boost::shared_ptr<AbsDecay>(new IsobarLSDecay(motherParticle, daughterParticles[0], daughterParticles[1], this)); else if (usedSystem=="Cano") tmpDec= boost::shared_ptr<AbsDecay>(new IsobarLSDecay(motherParticle, daughterParticles[0], daughterParticles[1], this));
else if (usedSystem=="Tensor") tmpDec= boost::shared_ptr<AbsDecay>(new IsobarTensorDecay(motherParticle, daughterParticles[0], daughterParticles[1], this));
else { else {
Alert << "used decay system\t" << usedSystem << "\tnot supported!!!\n" << endmsg; Alert << "used decay system\t" << usedSystem << "\tnot supported!!!\n" << endmsg;
exit(1); exit(1);
......
PwaUtils/EvtDataBaseList.hh
+ 2
− 0
View file @ fe4f5432
...@@ -41,6 +41,7 @@ typedef std::map<std::string, Vector4<double> > mapString4Vec; ...@@ -41,6 +41,7 @@ typedef std::map<std::string, Vector4<double> > mapString4Vec;
typedef std::map<std::string, map<Spin,map<Spin,map<Spin,complex<double> > > > > mapStringSpinComplex; typedef std::map<std::string, map<Spin,map<Spin,map<Spin,complex<double> > > > > mapStringSpinComplex;
typedef std::map<std::string, double> mapStringDouble; typedef std::map<std::string, double> mapStringDouble;
typedef std::map<std::string, map<Spin,map<Spin, complex<double> > > > mapStringComplDouble; typedef std::map<std::string, map<Spin,map<Spin, complex<double> > > > mapStringComplDouble;
typedef std::map<std::string, map<Spin,map<Spin,map<Spin, map<Spin, map<Spin,complex<double> > > > > > > mapString5SpinsComplex;
struct EvtData { struct EvtData {
mapInt4Vec FourVecsProd; mapInt4Vec FourVecsProd;
...@@ -51,6 +52,7 @@ struct EvtData { ...@@ -51,6 +52,7 @@ struct EvtData {
mapStringSpinComplex WignerDsString; mapStringSpinComplex WignerDsString;
mapStringDouble DoubleString; mapStringDouble DoubleString;
mapStringComplDouble ComplexDoubleString; mapStringComplDouble ComplexDoubleString;
mapString5SpinsComplex ComplexDouble5SpinString;
double evtWeight; double evtWeight;
int evtNo; int evtNo;
}; };
......
PwaUtils/IsobarTensorDecay.cc 0 → 100644
+ 209
− 0
View file @ fe4f5432
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
// IsobarTensorDecay class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#include <getopt.h>
#include <fstream>
#include <algorithm>
#include "PwaUtils/IsobarTensorDecay.hh"
#include "PwaUtils/AbsEnv.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "ErrLogger/ErrLogger.hh"
#include "Particle/Particle.hh"
#include "Utils/PawianCollectionUtils.hh"
#include "Utils/FunctionUtils.hh"
#include "PwaUtils/KinUtils.hh"
#include "PwaUtils/EvtDataBaseList.hh"
IsobarTensorDecay::IsobarTensorDecay(Particle* mother, Particle* daughter1, Particle* daughter2, AbsEnv* theEnv) :
IsobarLSDecay(mother, daughter1, daughter2, theEnv)
,_polMother(PolVector(_motherJPCPtr->J))
,_polDaughter1(PolVector(_daughter1JPCPtr->J))
,_polDaughter2(PolVector(_daughter2JPCPtr->J))
,_lctTensor(LeviCivitaTensor())
{
}
IsobarTensorDecay::IsobarTensorDecay(boost::shared_ptr<const jpcRes> motherJPCPtr, Particle* daughter1, Particle* daughter2, AbsEnv* theEnv, std::string motherName) :
IsobarLSDecay(motherJPCPtr, daughter1, daughter2, theEnv, motherName)
,_polMother(PolVector(_motherJPCPtr->J))
,_polDaughter1(PolVector(_daughter1JPCPtr->J))
,_polDaughter2(PolVector(_daughter2JPCPtr->J))
,_lctTensor(LeviCivitaTensor())
{
}
IsobarTensorDecay::~IsobarTensorDecay(){
}
void IsobarTensorDecay::print(std::ostream& os) const{
os << "\nJPCLS tensor amplitudes for decay\t" << _name << ":\n";
os << "suffix for fit parameter name:\t" << _fitParamSuffix << "\n";
std::vector< boost::shared_ptr<const JPCLS> >::const_iterator it;
for (it = _JPCLSDecAmps.begin(); it!= _JPCLSDecAmps.end(); ++it){
(*it)->print(os);
os << "\n";
}
AbsDecay::print(os);
os << "\n";
}
void IsobarTensorDecay::fillWignerDs(std::map<std::string, Vector4<double> >& fsMap, EvtData* evtData){
int evtNo=evtData->evtNo;
std::map<int, bool>::const_iterator it = _alreadyFilledMap.find(evtNo);
if(it!=_alreadyFilledMap.end() && it->second) return; //already filled
if (!_daughter1IsStable) _absDecDaughter1->fillWignerDs(fsMap, evtData);
if (!_daughter2IsStable) _absDecDaughter2->fillWignerDs(fsMap, evtData);
Vector4<double> all4Vec(0.,0.,0.,0.);
//fill all4Vec
std::map<std::string , Vector4<double> >::iterator itMap;
for(itMap=fsMap.begin(); itMap!=fsMap.end(); ++itMap){
all4Vec+=itMap->second;
}
//fill daughter1 and daughter2 4Vec
Vector4<double> daughter1_4Vec(0.,0.,0.,0.);
std::vector<Particle*> finalStatePart1;
if(0!=_absDecDaughter1) finalStatePart1=_absDecDaughter1->finalStateParticles();
else finalStatePart1.push_back(_daughter1);
std::vector<Particle*>::iterator itParticle;
for (itParticle= finalStatePart1.begin(); itParticle != finalStatePart1.end(); ++itParticle){
itMap=fsMap.find( (*itParticle)->name() );
daughter1_4Vec+=itMap->second;
}
// DebugMsg << "daughter1_4Vec:\t" << _daughter1->name() << "\t" << daughter1_4Vec << endmsg;
Vector4<double> daughter2_4Vec(0.,0.,0.,0.);
std::vector<Particle*> finalStatePart2;
if(0!=_absDecDaughter2) finalStatePart2=_absDecDaughter2->finalStateParticles();
else finalStatePart2.push_back(_daughter2);
for (itParticle= finalStatePart2.begin(); itParticle != finalStatePart2.end(); ++itParticle){
itMap=fsMap.find( (*itParticle)->name() );
daughter2_4Vec+=itMap->second;
}
// DebugMsg << "daughter2_4Vec:\t" << _daughter2->name() << "\t" << daughter2_4Vec << endmsg;
Vector4<double> mother_4Vec=daughter1_4Vec+daughter2_4Vec;
Spin spinMother=_motherJPCPtr->J;
Spin spinDaughter1=_daughter1JPCPtr->J;
Spin spinDaughter2=_daughter2JPCPtr->J;
if(spinMother>0){
_polMother.SetP4(mother_4Vec,mother_4Vec.M());
}
// _polMother.SetP4(mother_4Vec,mother_4Vec.M());
// DebugMsg << "_polMother:\t" << _polMother << endmsg;
if(spinDaughter1>0) _polDaughter1.SetP4(daughter1_4Vec, daughter1_4Vec.M());
if(spinDaughter2>0) _polDaughter2.SetP4(daughter2_4Vec, daughter2_4Vec.M());
Spin lam12Max=spinDaughter1+spinDaughter2;
if(lam12Max>spinMother) lam12Max=spinMother;
Spin lamMotherMax=spinMother;
if (!_hasMotherPart && spinMother>1) lamMotherMax=1; //attention: this is only valid for pbar p or e+ e- reactions; must be moved to individual specific classes
std::vector< boost::shared_ptr<const JPCLS> > theJPCLSAmps=JPCLSAmps();
// 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;
Spin L=(*itJPCLS)->L;
Spin S=(*itJPCLS)->S;
int s1s2S=spinDaughter1+spinDaughter2+(*itJPCLS)->S;
// bool add_lctForChi=true;
// if( s1s2S%2 ==0 ) add_lctForChi=false; //odd
// int SLms1=(*itJPCLS)->S+(*itJPCLS)->L+_motherJPCPtr->J;
// bool add_lctForTensor=true;
// if( SLms1%2 ==0 ) add_lctForTensor=false; //odd
OrbitalTensor orbTensor(L);
if (L>0) orbTensor.SetP4(daughter1_4Vec, daughter2_4Vec);
for (Spin lamMother=-lamMotherMax; lamMother<=lamMotherMax; ++lamMother){
DebugMsg << "lamMother:\t" << lamMother << endmsg;
Tensor<complex<double> > epsilonMotherProject(0);
if(spinMother>0) epsilonMotherProject=_polMother(lamMother);
// //calculate chi
// // Tensor<complex<double> > chi12;
// // Tensor<complex<double> > s12SpinProjector(0);
// // s12SpinProjector(0)=complex<double>(1.,0);
// // if (S>0) PolVector::Projector(S, 2*S, mother_4Vec, mother_4Vec.M(), s12SpinProjector);
for (Spin lamDaughter1=-spinDaughter1; lamDaughter1 <= spinDaughter1; ++lamDaughter1){
DebugMsg << "lamDaughter1:\t" << lamDaughter1 << endmsg;
// // Tensor<complex<double> > epsilonDaughter1Project(0);
// // epsilonDaughter1Project(0)=complex<double>(1.,0);
// // if(spinDaughter1>0) epsilonDaughter1Project=_polDaughter1(lamDaughter1);
for (Spin lamDaughter2=-spinDaughter2; lamDaughter2<=spinDaughter2; ++lamDaughter2){
DebugMsg << "lamDaughter2:\t" << lamDaughter2 << endmsg;
// // Tensor<complex<double> > epsilonDaughter2Project(0);
// // epsilonDaughter2Project(0)=complex<double>(1.,0);
// // if(spinDaughter2>0) epsilonDaughter2Project=_polDaughter2(lamDaughter2);
// // // chi12=conj(epsilonDaughter2Project)%conj(epsilonDaughter2Project);
// // if(add_lctForChi) chi12=epsilonMotherProject| _lctTensor | (conj(epsilonDaughter2Project)%conj(epsilonDaughter2Project));
// //chi12=epsilonMotherProject| (conj(epsilonDaughter2Project)%conj(epsilonDaughter2Project));
// // DebugMsg << "chi12 Tensor Amp:\t" << chi12 << endmsg;
// // Tensor<complex<double> > result=orbTensor|chi12;
// // if(add_lctForTensor) result=_lctTensor_lctTensor|result;
// // Tensor<complex<double> > result;
if(spinMother==0 && spinDaughter1 && spinDaughter1) evtData->ComplexDouble5SpinString[_name][L][S][lamMother][lamDaughter1][lamDaughter2]= complex<double> (1.,0.);
else if(S==0 && spinDaughter1 && spinDaughter1) {
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);
}
// // else{
// // Tensor<complex<double> > result;
// // if(add_lctForTensor) result=epsilonMotherProject | _lctTensor | orbTensor| chi12;
// // else result=epsilonMotherProject | orbTensor | chi12;
// // }
// // DebugMsg << "result Tensor Amp:\t" << result << endmsg;
}
}
}
}
}
PwaUtils/IsobarTensorDecay.hh 0 → 100644
+ 62
− 0
View file @ fe4f5432
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
// IsobarTensorDecay class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#pragma once
#include <iostream>
#include <vector>
#include <complex>
#include <map>
#include <vector>
#include <string>
#include <sstream>
#include <boost/shared_ptr.hpp>
#include "PwaUtils/IsobarLSDecay.hh"
#include "PwaUtils/DataUtils.hh"
#include "Utils/PawianCollectionUtils.hh"
class Particle;
class EvtData;
class AbsEnv;
class IsobarTensorDecay : public IsobarLSDecay{
public:
IsobarTensorDecay(Particle* mother, Particle* daughter1, Particle* daughter2, AbsEnv* theEnv);
IsobarTensorDecay(boost::shared_ptr<const jpcRes> motherJPCPtr, Particle* daughter1, Particle* daughter2, AbsEnv* theEnv, std::string motherName="pbarp");
virtual ~IsobarTensorDecay();
virtual void fillWignerDs(std::map<std::string , Vector4<double> >& fsMap, EvtData* evtData);
virtual void print(std::ostream& os) const;
virtual std::string type() {return "IsobarTensorDecay";}
protected:
PolVector _polMother;
PolVector _polDaughter1;
PolVector _polDaughter2;
LeviCivitaTensor _lctTensor;
};
PwaUtils/TensorDecAmps.cc 0 → 100644
+ 253
− 0
View file @ fe4f5432
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
// TensorDecAmps class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#include <getopt.h>
#include <fstream>
#include <string>
#include <mutex>
#include "PwaUtils/TensorDecAmps.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "ErrLogger/ErrLogger.hh"
#include "PwaUtils/DataUtils.hh"
#include "PwaUtils/IsobarTensorDecay.hh"
#include "Particle/Particle.hh"
TensorDecAmps::TensorDecAmps(boost::shared_ptr<IsobarTensorDecay> theDec) :
AbsXdecAmp(theDec)
,_JPCLSs(theDec->JPCLSAmps())
,_factorMag(1.)
{
Particle* daughter1=_decay->daughter1Part();
Particle* daughter2=_decay->daughter2Part();
if(_JPCLSs.size()>0) _factorMag=1./sqrt(_JPCLSs.size());
}
TensorDecAmps::TensorDecAmps(boost::shared_ptr<AbsDecay> theDec) :
AbsXdecAmp(theDec)
{
Particle* daughter1=_decay->daughter1Part();
Particle* daughter2=_decay->daughter2Part();
if(_JPCLSs.size()>0) _factorMag=1./sqrt(_JPCLSs.size());
}
TensorDecAmps::~TensorDecAmps()
{
}
complex<double> TensorDecAmps::XdecPartAmp(Spin lamX, Spin lamDec, short fixDaughterNr, EvtData* theData, Spin lamFs, AbsXdecAmp* grandmaAmp){
Spin lam1Min=-_Jdaughter1;
Spin lam1Max= _Jdaughter1;
Spin lam2Min=-_Jdaughter2;
Spin lam2Max=_Jdaughter2;
if(fixDaughterNr == 1){
lam1Min = lam1Max = lamDec;
}
else if(fixDaughterNr == 2){
lam2Min = lam2Max = lamDec;
}
else{
Alert << "Invalid fixDaughterNr in XdecPartAmp." << endmsg;
}
if(_enabledlamFsDaughter1){
lam1Min=lamFs;
lam1Max=lamFs;
}
else if(_enabledlamFsDaughter2){
lam2Min=lamFs;
lam2Max=lamFs;
}
complex<double> result=lsLoop(lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, false);
return result;
}
complex<double> TensorDecAmps::XdecAmp(Spin lamX, EvtData* theData, Spin lamFs, AbsXdecAmp* grandmaAmp){
complex<double> result(0.,0.);
if( fabs(lamX) > _JPCPtr->J) return result;
int evtNo=theData->evtNo;
if ( _cacheAmps && !_recalculate){
result=_cachedAmpMap[evtNo][lamX][lamFs];
result*=_absDyn->eval(theData, grandmaAmp);
return result;
}
// Spin lam1Min=-_Jdaughter1;
Spin lam1Min=-_Jdaughter1;
Spin lam1Max= _Jdaughter1;
Spin lam2Min=-_Jdaughter2;
Spin lam2Max=_Jdaughter2;
if(_enabledlamFsDaughter1){
lam1Min=lamFs;
lam1Max=lamFs;
}
else if(_enabledlamFsDaughter2){
lam2Min=lamFs;
lam2Max=lamFs;
}
result=lsLoop(lamX, theData, lam1Min, lam1Max, lam2Min, lam2Max, true, lamFs );
if ( _cacheAmps){
theMutex.lock();
_cachedAmpMap[evtNo][lamX][lamFs]=result;
theMutex.unlock();
}
result*=_absDyn->eval(theData, grandmaAmp);
return result;
}
complex<double> TensorDecAmps::lsLoop(Spin lamX, EvtData* theData, Spin lam1Min, Spin lam1Max, Spin lam2Min, Spin lam2Max, bool withDecs, Spin lamFs ){
complex<double> result(0.,0.);
std::vector< boost::shared_ptr<const JPCLS> >::iterator it;
for (it=_JPCLSs.begin(); it!=_JPCLSs.end(); ++it){
double theMag=_currentParamMags[*it];
double thePhi=_currentParamPhis[*it];
complex<double> expi(cos(thePhi), sin(thePhi));
for(Spin lambda1=lam1Min; lambda1<=lam1Max; ++lambda1){
for(Spin lambda2=lam2Min; lambda2<=lam2Max; ++lambda2){
Spin lambda = lambda1-lambda2;
if( fabs(lambda)>(*it)->J || fabs(lambda)>(*it)->S) continue;
complex<double> amp = theMag*expi*theData->ComplexDouble5SpinString[_name][(*it)->L][(*it)->S][lamX][lambda1][lambda2];
// Info << "theData->ComplexDouble5SpinString[" << _name << "]["
// << (*it)->L << "][" << (*it)->S << "]["
// << lamX <<"][" << lambda1 << "][" << lambda2 << "]:\t"
// << theData->ComplexDouble5SpinString[_name][(*it)->L][(*it)->S][lamX][lambda1][lambda2] << endmsg;
// Info << "amp:\t" << amp << endmsg;
if(withDecs) amp *=daughterAmp(lambda1, lambda2, theData, lamFs);
result+=amp;
}
}
}
return result;
}
void TensorDecAmps::getDefaultParams(fitParams& fitVal, fitParams& fitErr){
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > currentMagValMap;
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > currentPhiValMap;
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > currentMagErrMap;
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > currentPhiErrMap;
std::vector< boost::shared_ptr<const JPCLS> >::const_iterator itLS;
for(itLS=_JPCLSs.begin(); itLS!=_JPCLSs.end(); ++itLS){
currentMagValMap[*itLS]=_factorMag;
currentPhiValMap[*itLS]=0.;
currentMagErrMap[*itLS]=_factorMag;
currentPhiErrMap[*itLS]=0.3;
}
fitVal.Mags[_key]=currentMagValMap;
fitVal.Phis[_key]=currentPhiValMap;
fitErr.Mags[_key]=currentMagErrMap;
fitErr.Phis[_key]=currentPhiErrMap;
_absDyn->getDefaultParams(fitVal, fitErr);
if(!_daughter1IsStable) _decAmpDaughter1->getDefaultParams(fitVal, fitErr);
if(!_daughter2IsStable) _decAmpDaughter2->getDefaultParams(fitVal, fitErr);
}
void TensorDecAmps::print(std::ostream& os) const{
return; //dummy
}
bool TensorDecAmps::checkRecalculation(fitParams& theParamVal){
_recalculate=false;
if(_absDyn->checkRecalculation(theParamVal)) _recalculate=true;
if(!_daughter1IsStable) {
if(_decAmpDaughter1->checkRecalculation(theParamVal)) _recalculate=true;
}
if(!_daughter2IsStable){
if(_decAmpDaughter2->checkRecalculation(theParamVal)) _recalculate=true;
}
if(!_recalculate){
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >& magMap=theParamVal.Mags[_key];
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >& phiMap=theParamVal.Phis[_key];
std::vector< boost::shared_ptr<const JPCLS> >::iterator it;
for (it=_JPCLSs.begin(); it!=_JPCLSs.end(); ++it){
double theMag=magMap[*it];
double thePhi=phiMap[*it];
if(!CheckDoubleEquality(theMag, _currentParamMags[*it])){
_recalculate=true;
return _recalculate;
}
if(!CheckDoubleEquality(thePhi, _currentParamPhis[*it])){
_recalculate=true;
return _recalculate;
}
}
}
return _recalculate;
}
void TensorDecAmps::updateFitParams(fitParams& theParamVal){
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >& magMap=theParamVal.Mags[_key];
std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >& phiMap=theParamVal.Phis[_key];
std::vector< boost::shared_ptr<const JPCLS> >::iterator it;
for (it=_JPCLSs.begin(); it!=_JPCLSs.end(); ++it){
double theMag=magMap[*it];
double thePhi=phiMap[*it];
_currentParamMags[*it]=theMag;
_currentParamPhis[*it]=thePhi;
}
_absDyn->updateFitParams(theParamVal);
if(!_daughter1IsStable) _decAmpDaughter1->updateFitParams(theParamVal);
if(!_daughter2IsStable) _decAmpDaughter2->updateFitParams(theParamVal);
}
PwaUtils/TensorDecAmps.hh 0 → 100644
+ 82
− 0
View file @ fe4f5432
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
// TensorDecAmps class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#pragma once
#include <iostream>
#include <vector>
#include <complex>
#include <map>
#include <string>
#include <cassert>
#include <boost/shared_ptr.hpp>
#include "PwaUtils/AbsXdecAmp.hh"
class IsobarTensorDecay;
class AbsDecay;
class TensorDecAmps : public AbsXdecAmp{
public:
// create/copy/destroy:
///Constructor
TensorDecAmps(boost::shared_ptr<IsobarTensorDecay> theDec);
TensorDecAmps(boost::shared_ptr<AbsDecay> theDec);
/** Destructor */
virtual ~TensorDecAmps();
// Getters:
virtual complex<double> XdecAmp(Spin lamX, EvtData* theData, Spin lamFs, AbsXdecAmp* grandmaAmp);
virtual complex<double> XdecPartAmp(Spin lamX, Spin lamDec, short fixDaughterNr,
EvtData* theData, Spin lamFs, AbsXdecAmp* grandmaAmp);
virtual void print(std::ostream& os) const;
std::vector< boost::shared_ptr<const JPCLS> >& jpclsVec() {return _JPCLSs;}
virtual void getDefaultParams(fitParams& fitVal, fitParams& fitErr);
virtual bool checkRecalculation(fitParams& theParamVal);
virtual void updateFitParams(fitParams& theParamVal);
protected:
std::vector< boost::shared_ptr<const JPCLS> > _JPCLSs;
double _factorMag;
virtual complex<double> lsLoop(Spin lamX, EvtData* theData, Spin lam1Min, Spin lam1Max, Spin lam2Min, Spin lam2Max, bool withDecs, Spin lamFs=0 );
private:
};
PwaUtils/XdecAmpRegistry.cc
+ 6
− 0
View file @ fe4f5432
...@@ -32,11 +32,13 @@ ...@@ -32,11 +32,13 @@
#include "PwaUtils/AbsDecay.hh" #include "PwaUtils/AbsDecay.hh"
#include "PwaUtils/IsobarLSDecay.hh" #include "PwaUtils/IsobarLSDecay.hh"
#include "PwaUtils/IsobarHeliDecay.hh" #include "PwaUtils/IsobarHeliDecay.hh"
#include "PwaUtils/IsobarTensorDecay.hh"
#include "PwaUtils/OmegaTo3PiLSDecay.hh" #include "PwaUtils/OmegaTo3PiLSDecay.hh"
#include "PwaUtils/OmegaTo3PiTensorDecay.hh" #include "PwaUtils/OmegaTo3PiTensorDecay.hh"
#include "PwaUtils/AbsXdecAmp.hh" #include "PwaUtils/AbsXdecAmp.hh"
#include "PwaUtils/LSDecAmps.hh" #include "PwaUtils/LSDecAmps.hh"
#include "PwaUtils/HeliDecAmps.hh" #include "PwaUtils/HeliDecAmps.hh"
#include "PwaUtils/TensorDecAmps.hh"
#include "PwaUtils/LSOmegaTo3PiDecAmps.hh" #include "PwaUtils/LSOmegaTo3PiDecAmps.hh"
#include "PwaUtils/TensorOmegaTo3PiDecAmps.hh" #include "PwaUtils/TensorOmegaTo3PiDecAmps.hh"
#include "ErrLogger/ErrLogger.hh" #include "ErrLogger/ErrLogger.hh"
...@@ -74,6 +76,10 @@ boost::shared_ptr<AbsXdecAmp> XdecAmpRegistry::getXdecAmp(boost::shared_ptr<AbsD ...@@ -74,6 +76,10 @@ boost::shared_ptr<AbsXdecAmp> XdecAmpRegistry::getXdecAmp(boost::shared_ptr<AbsD
boost::shared_ptr<IsobarHeliDecay> decLamLam = boost::dynamic_pointer_cast<IsobarHeliDecay>(theAbsXDec); boost::shared_ptr<IsobarHeliDecay> decLamLam = boost::dynamic_pointer_cast<IsobarHeliDecay>(theAbsXDec);
result=boost::shared_ptr<AbsXdecAmp>(new HeliDecAmps(decLamLam)); result=boost::shared_ptr<AbsXdecAmp>(new HeliDecAmps(decLamLam));
} }
else if(theAbsXDec->type()=="IsobarTensorDecay"){
boost::shared_ptr<IsobarTensorDecay> decTensor = boost::dynamic_pointer_cast<IsobarTensorDecay>(theAbsXDec);
result=boost::shared_ptr<AbsXdecAmp>(new TensorDecAmps(decTensor));
}
else if(theAbsXDec->type()=="OmegaTo3PiLSDecay"){ else if(theAbsXDec->type()=="OmegaTo3PiLSDecay"){
boost::shared_ptr<OmegaTo3PiLSDecay> decOmega = boost::dynamic_pointer_cast<OmegaTo3PiLSDecay>(theAbsXDec); boost::shared_ptr<OmegaTo3PiLSDecay> decOmega = boost::dynamic_pointer_cast<OmegaTo3PiLSDecay>(theAbsXDec);
result=boost::shared_ptr<AbsXdecAmp>(new LSOmegaTo3PiDecAmps(decOmega)); result=boost::shared_ptr<AbsXdecAmp>(new LSOmegaTo3PiDecAmps(decOmega));
......
pbarpUtils/pbarpEnv.cc
+ 18
− 4
View file @ fe4f5432
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include "PwaUtils/AbsDecay.hh" #include "PwaUtils/AbsDecay.hh"
#include "PwaUtils/AbsDecayList.hh" #include "PwaUtils/AbsDecayList.hh"
#include "PwaUtils/IsobarLSDecay.hh" #include "PwaUtils/IsobarLSDecay.hh"
#include "PwaUtils/IsobarTensorDecay.hh"
//#include "PwaUtils/IsobarDecayList.hh" //#include "PwaUtils/IsobarDecayList.hh"
#include "pbarpUtils/pbarpReaction.hh" #include "pbarpUtils/pbarpReaction.hh"
//#include "pbarpUtils/pbarpEventList.hh" //#include "pbarpUtils/pbarpEventList.hh"
...@@ -92,10 +93,23 @@ void pbarpEnv::setup(pbarpParser* thePbarpParser){ ...@@ -92,10 +93,23 @@ void pbarpEnv::setup(pbarpParser* thePbarpParser){
_pbarpReaction=boost::shared_ptr<pbarpReaction>(new pbarpReaction(_producedParticlePairs, _lmax)); _pbarpReaction=boost::shared_ptr<pbarpReaction>(new pbarpReaction(_producedParticlePairs, _lmax));
//fill prodDecayList //fill prodDecayList
std::vector< boost::shared_ptr<IsobarLSDecay> > prodDecs= _pbarpReaction->productionDecays(); if(thePbarpParser->productionFormalism()=="Cano"){
std::vector< boost::shared_ptr<IsobarLSDecay> >::iterator itDec; std::vector< boost::shared_ptr<IsobarLSDecay> > prodDecs= _pbarpReaction->productionDecays();
for (itDec=prodDecs.begin(); itDec!=prodDecs.end(); ++itDec){ std::vector< boost::shared_ptr<IsobarLSDecay> >::iterator itDec;
_prodDecList->addDecay(*itDec); for (itDec=prodDecs.begin(); itDec!=prodDecs.end(); ++itDec){
_prodDecList->addDecay(*itDec);
}
}
else if(thePbarpParser->productionFormalism()=="Tensor"){
std::vector< boost::shared_ptr<IsobarTensorDecay> > prodDecs= _pbarpReaction->productionTensorDecays();
std::vector< boost::shared_ptr<IsobarTensorDecay> >::iterator itDec;
for (itDec=prodDecs.begin(); itDec!=prodDecs.end(); ++itDec){
_prodDecList->addDecay(*itDec);
}
}
else{
Alert <<"production formalism\t" << thePbarpParser->productionFormalism() << "\t is not supported!!!" << endmsg;
exit(0);
} }
//set suffixes //set suffixes
......
pbarpUtils/pbarpReaction.cc
+ 5
− 0
View file @ fe4f5432
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include "pbarpUtils/pbarpEnv.hh" #include "pbarpUtils/pbarpEnv.hh"
#include "PwaUtils/IsobarLSDecay.hh" #include "PwaUtils/IsobarLSDecay.hh"
#include "PwaUtils/IsobarHeliDecay.hh" #include "PwaUtils/IsobarHeliDecay.hh"
#include "PwaUtils/IsobarTensorDecay.hh"
#include "qft++/relativistic-quantum-mechanics/Utils.hh" #include "qft++/relativistic-quantum-mechanics/Utils.hh"
#include "ErrLogger/ErrLogger.hh" #include "ErrLogger/ErrLogger.hh"
#include "Particle/Particle.hh" #include "Particle/Particle.hh"
...@@ -70,6 +71,10 @@ pbarpReaction::pbarpReaction(std::vector<std::pair<Particle*, Particle*> >& prod ...@@ -70,6 +71,10 @@ pbarpReaction::pbarpReaction(std::vector<std::pair<Particle*, Particle*> >& prod
CheckLmaxLimit((std::string&)(*itPartPairs).first->name(), *itJPC) && CheckLmaxLimit((std::string&)(*itPartPairs).first->name(), *itJPC) &&
CheckLmaxLimit((std::string&)(*itPartPairs).second->name(), *itJPC)){ CheckLmaxLimit((std::string&)(*itPartPairs).second->name(), *itJPC)){
_prodDecs.push_back(currentDec); _prodDecs.push_back(currentDec);
boost::shared_ptr<IsobarTensorDecay> currentTensorDec(new IsobarTensorDecay( (*itJPC),itPartPairs->first, itPartPairs->second, pbarpEnv::instance(), decName));
_prodTensorDecs.push_back(currentTensorDec);
boost::shared_ptr<IsobarHeliDecay> currentHeliDec(new IsobarHeliDecay( (*itJPC),itPartPairs->first, itPartPairs->second, pbarpEnv::instance(), decName)); boost::shared_ptr<IsobarHeliDecay> currentHeliDec(new IsobarHeliDecay( (*itJPC),itPartPairs->first, itPartPairs->second, pbarpEnv::instance(), decName));
_prodHeliDecs.push_back(currentHeliDec); _prodHeliDecs.push_back(currentHeliDec);
acceptJPC=true; acceptJPC=true;
......
pbarpUtils/pbarpReaction.hh
+ 4
− 1
View file @ fe4f5432
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
class Particle; class Particle;
class IsobarLSDecay; class IsobarLSDecay;
class IsobarHeliDecay; class IsobarHeliDecay;
class IsobarTensorDecay;
class pbarpReaction { class pbarpReaction {
...@@ -51,7 +52,8 @@ public: ...@@ -51,7 +52,8 @@ public:
virtual void print(std::ostream& os) const; virtual void print(std::ostream& os) const;
std::vector< boost::shared_ptr<IsobarLSDecay> >& productionDecays() {return _prodDecs;} std::vector< boost::shared_ptr<IsobarLSDecay> >& productionDecays() {return _prodDecs;}
std::vector< boost::shared_ptr<IsobarHeliDecay> >& productionHeliDecays() {return _prodHeliDecs;} std::vector< boost::shared_ptr<IsobarHeliDecay> >& productionHeliDecays() {return _prodHeliDecs;}
std::vector< boost::shared_ptr<IsobarTensorDecay> >& productionTensorDecays() {return _prodTensorDecs;}
std::vector< boost::shared_ptr<const jpcRes> >& jpcStates() {return _pbarpJPCs;} std::vector< boost::shared_ptr<const jpcRes> >& jpcStates() {return _pbarpJPCs;}
std::vector< boost::shared_ptr<const JPCLS> >& jpclsStates() {return _pbarpJPCLSs;} std::vector< boost::shared_ptr<const JPCLS> >& jpclsStates() {return _pbarpJPCLSs;}
std::vector< boost::shared_ptr<const JPCLS> >& jpclsSingletStates() {return _pbarpJPCLSsinglet;} std::vector< boost::shared_ptr<const JPCLS> >& jpclsSingletStates() {return _pbarpJPCLSsinglet;}
...@@ -72,6 +74,7 @@ private: ...@@ -72,6 +74,7 @@ private:
std::vector< boost::shared_ptr<IsobarLSDecay> > _prodDecs; std::vector< boost::shared_ptr<IsobarLSDecay> > _prodDecs;
std::vector< boost::shared_ptr<IsobarHeliDecay> > _prodHeliDecs; std::vector< boost::shared_ptr<IsobarHeliDecay> > _prodHeliDecs;
std::vector< boost::shared_ptr<IsobarTensorDecay> > _prodTensorDecs;
std::map< boost::shared_ptr<const jpcRes>, Spin, pawian::Collection::SharedPtrLess > _minLMap; std::map< boost::shared_ptr<const jpcRes>, Spin, pawian::Collection::SharedPtrLess > _minLMap;
......
pbarpUtils/pbarpTensorLh.cc 0 → 100644
+ 95
− 0
View file @ fe4f5432
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
// pbarpTensorLh class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#include <getopt.h>
#include <fstream>
#include <string>
#include "pbarpUtils/pbarpTensorLh.hh"
#include "pbarpUtils/pbarpEnv.hh"
#include "pbarpUtils/pbarpReaction.hh"
#include "PwaUtils/LSDecAmps.hh"
#include "PwaUtils/EvtDataBaseList.hh"
#include "PwaUtils/AbsXdecAmp.hh"
#include "PwaUtils/AbsDecay.hh"
#include "PwaUtils/IsobarTensorDecay.hh"
#include "PwaUtils/FitParamsBase.hh"
#include "PwaUtils/XdecAmpRegistry.hh"
#include "Particle/Particle.hh"
#include "ErrLogger/ErrLogger.hh"
#include <boost/bind.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
// pbarpTensorLh::pbarpTensorLh(boost::shared_ptr<const EvtDataBaseList> theEvtList) :
// pbarpBaseLh(theEvtList)
// {
// initialize();
// }
pbarpTensorLh::pbarpTensorLh() :
pbarpBaseLh()
{
initialize();
}
pbarpTensorLh::~pbarpTensorLh()
{;
}
void pbarpTensorLh::print(std::ostream& os) const{
}
void pbarpTensorLh::initialize(){
std::vector< boost::shared_ptr<IsobarTensorDecay> > theDecs = _pbarpReactionPtr->productionTensorDecays();
std::vector< boost::shared_ptr<IsobarTensorDecay> >::iterator it;
for (it=theDecs.begin(); it!=theDecs.end(); ++it){
boost::shared_ptr<AbsDecay> currentDec( (*it).get() );
boost::shared_ptr<AbsXdecAmp> currentAmp=XdecAmpRegistry::instance()->getXdecAmp(currentDec);
_decAmps.push_back(currentAmp);
}
std::vector< boost::shared_ptr<const JPCLS> > jpclsSingletStates=_pbarpReactionPtr->jpclsSingletStates();
fillMap(jpclsSingletStates, _decAmps, _decAmpsSinglet);
std::vector< boost::shared_ptr<const JPCLS> > jpclsTriplet0States=_pbarpReactionPtr->jpclsTriplet0States();
fillMap(jpclsTriplet0States, _decAmps, _decAmpsTriplet0);
std::vector< boost::shared_ptr<const JPCLS> > jpclsTripletp1States=_pbarpReactionPtr->jpclsTripletp1States();
fillMap(jpclsTripletp1States, _decAmps, _decAmpsTripletp1);
std::vector< boost::shared_ptr<const JPCLS> > jpclsTripletm1States=_pbarpReactionPtr->jpclsTripletm1States();
fillMap(jpclsTripletm1States, _decAmps, _decAmpsTripletm1);
}
pbarpUtils/pbarpTensorLh.hh 0 → 100644
+ 72
− 0
View file @ fe4f5432
//************************************************************************//
// //
// Copyright 2013 Bertram Kopf (bertram@ep1.rub.de) //
// Julian Pychy (julian@ep1.rub.de) //
// - Ruhr-Universität Bochum //
// //
// This file is part of Pawian. //
// //
// Pawian is free software: you can redistribute it and/or modify //
// it under the terms of the GNU General Public License as published by //
// the Free Software Foundation, either version 3 of the License, or //
// (at your option) any later version. //
// //
// Pawian is distributed in the hope that it will be useful, //
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
// GNU General Public License for more details. //
// //
// You should have received a copy of the GNU General Public License //
// along with Pawian. If not, see <http://www.gnu.org/licenses/>. //
// //
//************************************************************************//
// pbarpTensorLh class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#pragma once
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <complex>
#include <boost/shared_ptr.hpp>
#include <boost/function.hpp>
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
#include "pbarpUtils/pbarpBaseLh.hh"
#include "PwaUtils/DataUtils.hh"
#include "Minuit2/MnUserParameters.h"
class AbsXdecAmp;
class pbarpReaction;
class LSDecAmps;
class pbarpTensorLh : public pbarpBaseLh {
public:
// pbarpTensorLh(boost::shared_ptr<const EvtDataBaseList>);
pbarpTensorLh();
virtual ~pbarpTensorLh();
virtual AbsLh* clone_() const {
AbsLh* theClone=new pbarpTensorLh();
theClone->setDataVec(_evtDataVec);
theClone->setMcVec(_evtMCVec);
return theClone;
}
virtual void print(std::ostream& os) const;
protected:
private:
void initialize();
};
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