//************************************************************************//
// //
// 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/>. //
// //
//************************************************************************//
// gammapBaseLh class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf
#include <getopt.h>
#include <fstream>
#include <string>
#include <math.h>
#include "gammapUtils/gammapBaseLh.hh"
#include "gammapUtils/gammapReaction.hh"
#include "gammapUtils/GammapChannelEnv.hh"
#include "PwaUtils/GlobalEnv.hh"
#include "PwaUtils/EvtDataBaseList.hh"
#include "PwaUtils/AbsXdecAmp.hh"
#include "PwaUtils/AbsDecay.hh"
#include "PwaUtils/IsobarLSDecay.hh"
#include "FitParams/FitParColBase.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>
gammapBaseLh::gammapBaseLh(ChannelID channelID) :
AbsLh(channelID)
{
initialize();
}
gammapBaseLh::~gammapBaseLh()
{;
}
complex<double> gammapBaseLh::calcSpinDensity(Spin M1, Spin M2, std::string& nameDec, EvtData* theData){
complex<double> result(0.,0.);
return result;
}
complex<double> gammapBaseLh::calcProdPartAmp(Spin lamX, Spin lamDec, std::string nameDec, EvtData* theData,
std::map <std::shared_ptr<const JPCLS>,
std::vector< std::shared_ptr<AbsXdecAmp> >,
pawian::Collection::SharedPtrLess > pbarpAmps){
complex<double> resultAmp(0.,0.);
return resultAmp;
}
// double gammapBaseLh::calcEvtIntensity(EvtData* theData, fitParCol& theParamVal){
// double result=0.;
// std::vector<std::shared_ptr<AbsXdecAmp> >::iterator itDec;
// std::map< Spin, std::map< Spin, std::map< Spin, complex<double> > > > mapPinGamiPout;
// for (Spin mzPin=-1./2; mzPin<=1./2; mzPin++){
// for (Spin mzGam=-1.; mzGam<=1.; mzGam+=2){
// for (Spin lamPout=-1./2; lamPout<=1./2; lamPout++){
// mapPinGamiPout[mzPin][mzGam][lamPout]=complex<double> (0.,0.);
// }
// }
// }
// std::map< std::shared_ptr<const jpcRes>, std::vector< std::shared_ptr<const JPCLS> >, pawian::Collection::SharedPtrLess>::iterator itjpc;
// std::vector< std::shared_ptr<const JPCLS> >::iterator itJPClj;
// for(itDec=_decAmps.begin(); itDec !=_decAmps.end(); ++itDec){
// std::shared_ptr<const jpcRes> currentJPC=(*itDec)->jpcPtr();
// double iso12Val=_currentParamJPCIsos12.at(currentJPC);
// double iso32Val=_currentParamJPCIsos32.at(currentJPC);
// double isoFactor=iso12Val;
// if( int(2.*(*itDec)->absDec()->motherIGJPC()->I)==3) isoFactor=iso32Val;
// double theMag=_currentParamMags.at(currentJPC);
// double thePhi=0.;
// if( 2.*currentJPC->J > 1.) thePhi=_currentParamPhis.at(currentJPC);
// else if(currentJPC->P == 1){ // 1/2+ M1
// thePhi=M_PI/2.;
// }
// complex<double> current_amp_prod(0.,0.);
// std::vector< std::shared_ptr<const JPCLS> > currentJPCljVec = _jpcToJPCljMap.at(currentJPC);
// for(itJPClj=currentJPCljVec.begin(); itJPClj!=currentJPCljVec.end(); ++itJPClj){
// int lmp=(*itJPClj)->L;
// Spin jmp=(*itJPClj)->S;
// Spin J=(*itJPClj)->J;
// double elMagFactor=1.;
// if(currentJPC->P == pow(-1.,(int)(J + 1./2))){ // 1/2-,3/2+,5/2-,...
// if(jmp==J+Spin(1./2)){ //el. multipole
// elMagFactor=cos(thePhi);
// }
// else elMagFactor=sin(thePhi); //mag multipole
// }
// else{ // 1/2+,3/2-,5/2+,...
// if(jmp==J+Spin(1./2)){ //mag. multipole
// elMagFactor=sin(thePhi);
// }
// else elMagFactor=cos(thePhi); //el. multipole
// }
// for (Spin mzPin=-1./2; mzPin<=1./2; mzPin++){
// for (Spin mzGam=-1.; mzGam<=1.; mzGam+=2.){
// Spin mz = mzGam + mzPin; // spin-j projection
// if( J < fabs(mz) || jmp < fabs(mzGam) ) continue;
// complex<double> current_amp_prod=theMag*elMagFactor*isoFactor*sqrt(2*lmp+1)*Clebsch(1, mzGam,lmp,0,jmp,mzGam)*Clebsch(jmp, mzGam,1/2., mzPin,J,mz);
// for (Spin lamPout=-1./2; lamPout<=1./2; lamPout++){
// mapPinGamiPout.at(mzPin).at(mzGam).at(lamPout)+=current_amp_prod*(*itDec)->XdecAmp(mz, theData, lamPout);
// }
// }
// }
// }
// }
// std::map< Spin, std::map< Spin, std::map< Spin, complex<double> > > >::iterator itAmpMap;
// for(itAmpMap=mapPinGamiPout.begin(); itAmpMap!=mapPinGamiPout.end(); ++itAmpMap){
// std::map< Spin, std::map< Spin, complex<double> > >::iterator itAmpMap2;
// for(itAmpMap2=itAmpMap->second.begin(); itAmpMap2!=itAmpMap->second.end(); ++itAmpMap2){
// std::map< Spin, complex<double> >::iterator itAmpMap3;
// for(itAmpMap3=itAmpMap2->second.begin(); itAmpMap3!=itAmpMap2->second.end(); ++itAmpMap3){
// result+=norm(itAmpMap3->second);
// }
// }
// }
// if(_usePhasespace) result+=theParamVal.otherParams[_phasespaceKey];
// result *= theParamVal.otherParams.at(_channelScaleParam);
// return result;
// }
// void gammapBaseLh::getDefaultParams(fitParCol& fitVal, fitParCol& fitErr){
// AbsLh::getDefaultParams(fitVal, fitErr);
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess>::iterator itIso;
// for(itIso=_currentParamJPCIsos12.begin(); itIso!=_currentParamJPCIsos12.end(); ++itIso){
// fitVal.otherParams["Iso12"+(*itIso).first->name()+"Range01"]=(*itIso).second;
// fitErr.otherParams["Iso12"+(*itIso).first->name()+"Range01"]=0.5;
// }
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess > currentMagValMap;
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess > currentPhiValMap;
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess > currentMagErrMap;
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess > currentPhiErrMap;
// double magFactor=1./sqrt(_jpcToIGJPCMap.size());
// for ( std::map< std::shared_ptr<const jpcRes>, std::vector< std::shared_ptr<const IGJPC> >, pawian::Collection::SharedPtrLess>::iterator it = _jpcToIGJPCMap.begin(); it!=_jpcToIGJPCMap.end(); ++it){
// currentMagValMap[it->first] = magFactor;
// currentMagErrMap[it->first] = magFactor;
// if(it->first->J > 1./2){
// currentPhiValMap[it->first] = 0.;
// currentPhiErrMap[it->first] = 0.3;
// }
// }
// fitVal.MagsJPC["gammap"]=currentMagValMap;
// fitVal.PhisJPC["gammap"]=currentPhiValMap;
// fitErr.MagsJPC["gammap"]=currentMagErrMap;
// fitErr.PhisJPC["gammap"]=currentPhiErrMap;
// }
void gammapBaseLh::fillDefaultParams(std::shared_ptr<AbsPawianParameters> fitPar){
AbsLh::fillDefaultParams(fitPar);
std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess>::iterator itIso;
for(itIso=_currentParamJPCIsos12.begin(); itIso!=_currentParamJPCIsos12.end(); ++itIso){
fitPar->Add("Iso12"+(*itIso).first->name(), (*itIso).second, 0.5);
fitPar->SetLimits("Iso12"+(*itIso).first->name(), 0., 1.);
}
// double magFactor=1./sqrt(_jpcToIGJPCMap.size());
//attention: missing things have to be added !!!!
}
void gammapBaseLh::fillIsos(){
//first look for iso0 decay amplitudes and fill it in the map
std::vector< std::shared_ptr<AbsXdecAmp> >::iterator it;
for(it= _decAmps.begin(); it!= _decAmps.end(); ++it){
if( int(2.*(*it)->absDec()->motherIGJPC()->I) ==1){ //I=1/2
_iso12DecAmps.push_back(*it);
}
if(int(2.*(*it)->absDec()->motherIGJPC()->I) ==3){
_iso32DecAmps.push_back(*it);
}
}
for(it= _iso12DecAmps.begin(); it!= _iso12DecAmps.end(); ++it){
_currentParamJPCIsos12[(*it)->jpcPtr()]=1.;
_currentParamJPCIsos32[(*it)->jpcPtr()]=0.;
}
for(it= _iso32DecAmps.begin(); it!= _iso32DecAmps.end(); ++it){
_currentParamJPCIsos32[(*it)->jpcPtr()]=1.;
_currentParamJPCIsos12[(*it)->jpcPtr()]=0.;
}
//now look if iso12 and iso32 exisist
for(it= _iso12DecAmps.begin(); it!= _iso12DecAmps.end(); ++it){
// find corresponding iso32 partner
std::vector< std::shared_ptr<AbsXdecAmp> >::iterator itIso1;
for(itIso1= _iso32DecAmps.begin(); itIso1!= _iso32DecAmps.end(); ++itIso1){
if( (*(*it)->jpcPtr()) == (*(*itIso1)->jpcPtr()) ){
_currentParamJPCIsos12[(*it)->jpcPtr()]=sqrt(.5);
_currentParamJPCIsos32[(*it)->jpcPtr()]=sqrt(.5);
}
}
}
//now look if iso12 and iso32 exisist for the same amplitude
for(it= _iso12DecAmps.begin(); it!= _iso12DecAmps.end(); ++it){
std::vector< std::shared_ptr<AbsXdecAmp> >::iterator it2;
for(it2= _iso32DecAmps.begin(); it2!= _iso32DecAmps.end(); ++it2){
if((*it)->jpcDecsName()==(*it2)->jpcDecsName()){
//found
_currentParamJPCIsos12[(*it)->jpcPtr()]=sqrt(.7);
_currentParamJPCIsos32[(*it)->jpcPtr()]=sqrt(.3);
}
}
}
}
// void gammapBaseLh::updateFitParams(fitParCol& theParamVal){
// AbsLh::updateFitParams(theParamVal);
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess>::iterator itIso;
// for(itIso=_currentParamJPCIsos12.begin(); itIso!=_currentParamJPCIsos12.end(); ++ itIso){
// double theVal=theParamVal.otherParams.at("Iso12"+(*itIso).first->name()+"Range01");
// (*itIso).second=theVal;
// _currentParamJPCIsos32[(*itIso).first]=sqrt(1.-theVal*theVal);
// }
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess >& magMap=theParamVal.MagsJPC.at("gammap");
// std::map< std::shared_ptr<const jpcRes>, double, pawian::Collection::SharedPtrLess >& phiMap=theParamVal.PhisJPC.at("gammap");
// std::map< std::shared_ptr<const jpcRes>, std::vector< std::shared_ptr<const IGJPC> >, pawian::Collection::SharedPtrLess>::iterator it;
// for (it= _jpcToIGJPCMap.begin(); it!=_jpcToIGJPCMap.end(); ++it){
// double theMag=magMap.at(it->first);
// _currentParamMags[it->first]=theMag;
// if(it->first->J > 1./2){
// double thePhi=phiMap.at(it->first);
// _currentParamPhis[it->first]=thePhi;
// }
// }
// }
void gammapBaseLh::print(std::ostream& os) const{
}
void gammapBaseLh::initialize(){
// check that only one non spin0 final particle exists
std::vector<Particle*> fsParticles=std::static_pointer_cast<GammapChannelEnv>(GlobalEnv::instance()->GammapChannel(_channelID))->finalStateParticles();
std::vector<Particle*>::iterator itParticle;
bool spinParticleFound=false;
for (itParticle=fsParticles.begin(); itParticle != fsParticles.end(); ++itParticle){
int current2J = (*itParticle)->twoJ();
if(current2J>0.){
if(spinParticleFound){
Alert << "only one non spin0 final particle allowed!!!" << endmsg;
exit(0);
}
spinParticleFound=true;
}
}
_gammapReactionPtr = std::static_pointer_cast<GammapChannelEnv>(GlobalEnv::instance()->GammapChannel(_channelID))->reaction();
_jpcToIGJPCMap = _gammapReactionPtr->jpcToIGJPCMap();
_jpcToJPCljMap = _gammapReactionPtr->jpcToJPCljMap();
std::vector< std::shared_ptr<IsobarLSDecay> > theDecs = _gammapReactionPtr->productionDecays();
std::vector< std::shared_ptr<IsobarLSDecay> >::iterator it;
for (it=theDecs.begin(); it!=theDecs.end(); ++it){
std::shared_ptr<AbsXdecAmp> currentAmp=XdecAmpRegistry::instance()->getXdecAmp(_channelID, (*it)->absDecPtr());
_decAmps.push_back(currentAmp);
}
fillIsos();
}