IsobarDecay.cc

// IsobarDecay class definition file. -*- C++ -*-
// Copyright 2012 Bertram Kopf

#include <getopt.h>
#include <fstream>
#include <algorithm>

#include "PwaUtils/IsobarDecay.hh"
#include "PwaUtils/IsobarDecayList.hh"
#include "pbarpUtils/pbarpEnv.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"

IsobarDecay::IsobarDecay(Particle* mother, Particle* daughter1, Particle* daughter2) :
  _mother(mother)
  ,_daughter1(daughter1)
  ,_daughter2(daughter2)
  ,_daughter1IsStable(true)
  ,_daughter2IsStable(true)
  ,_hasMotherPart(true)
  ,_withDyn(false)
  ,_motherJPCPtr(getJPCPtr(mother)) 
  ,_daughter1JPCPtr(getJPCPtr(daughter1))
  ,_daughter2JPCPtr(getJPCPtr(daughter2))
  ,_name(mother->name()+"To"+daughter1->name()+"_"+daughter2->name())
  ,_fitParamSuffix(_name)
  ,_massParamKey(mother->name())
{
  validJPCLS( _motherJPCPtr, _daughter1JPCPtr, _daughter2JPCPtr, _JPCLSDecAmps);
  _isoDecDaughter1=pbarpEnv::instance()->decayList()->decay(_daughter1);
  if(0 != _isoDecDaughter1){
    _daughter1IsStable=false;
    std::vector<Particle*> fsParticlesDaughter1=_isoDecDaughter1->finalStateParticles();
    _finalStateParticles.insert(_finalStateParticles.end(), fsParticlesDaughter1.begin(), fsParticlesDaughter1.end());
  }
  else _finalStateParticles.push_back(daughter1);
  _isoDecDaughter2=pbarpEnv::instance()->decayList()->decay(_daughter2);

  if(0 != _isoDecDaughter2){
    _daughter2IsStable=false;
    _finalStateParticlesDaughter2=_isoDecDaughter2->finalStateParticles();
    _finalStateParticles.insert(_finalStateParticles.end(), _finalStateParticlesDaughter2.begin(), _finalStateParticlesDaughter2.end());
  }
  else{
    _finalStateParticles.push_back(daughter2);
    _finalStateParticlesDaughter2.push_back(daughter2);
  }

  pawian::Collection::PtrLess thePtrLess;
  std::sort(_finalStateParticles.begin(), _finalStateParticles.end(), thePtrLess);
  _wignerDKey=FunctionUtils::particleListName(_finalStateParticlesDaughter2)+"_"+_motherJPCPtr->name()+FunctionUtils::particleListName(_finalStateParticles);

}

IsobarDecay::IsobarDecay(boost::shared_ptr<const jpcRes> motherJPCPtr, Particle* daughter1, Particle* daughter2, std::string motherName) :
  _mother(0)
  ,_daughter1(daughter1)
  ,_daughter2(daughter2)
  ,_daughter1IsStable(true)
  ,_daughter2IsStable(true)
  ,_hasMotherPart(false)
  ,_withDyn(false)
  ,_motherJPCPtr(motherJPCPtr)
  ,_daughter1JPCPtr(getJPCPtr(daughter1))
  ,_daughter2JPCPtr(getJPCPtr(daughter2))
  ,_name(motherName+"To"+daughter1->name()+"_"+daughter2->name())
  ,_fitParamSuffix(_name)
  ,_massParamKey("")
{
  validJPCLS( _motherJPCPtr, _daughter1JPCPtr, _daughter2JPCPtr, _JPCLSDecAmps);
  _isoDecDaughter1=pbarpEnv::instance()->decayList()->decay(_daughter1);


  if(0 != _isoDecDaughter1){
    _daughter1IsStable=false;
    std::vector<Particle*> fsParticlesDaughter1=_isoDecDaughter1->finalStateParticles();
    _finalStateParticles.insert(_finalStateParticles.end(), fsParticlesDaughter1.begin(), fsParticlesDaughter1.end());
  }
  else _finalStateParticles.push_back(daughter1);
  _isoDecDaughter2=pbarpEnv::instance()->decayList()->decay(_daughter2);

  if(0 != _isoDecDaughter2){
    _daughter2IsStable=false;
    _finalStateParticlesDaughter2=_isoDecDaughter2->finalStateParticles();
    _finalStateParticles.insert(_finalStateParticles.end(), _finalStateParticlesDaughter2.begin(), _finalStateParticlesDaughter2.end());
  }
  else{
    _finalStateParticles.push_back(daughter2);
    _finalStateParticlesDaughter2.push_back(daughter2);
  }
  pawian::Collection::PtrLess thePtrLess;
  std::sort(_finalStateParticles.begin(), _finalStateParticles.end(), thePtrLess);

  _wignerDKey=FunctionUtils::particleListName(_finalStateParticlesDaughter2)+"_"+motherName;
}

IsobarDecay::~IsobarDecay(){
}

void IsobarDecay::fillWignerDs(std::map<std::string, Vector4<double> >& fsMap, EvtData* evtData){
  if (!_daughter1IsStable) _isoDecDaughter1->fillWignerDs(fsMap, evtData);
  if (!_daughter2IsStable) _isoDecDaughter2->fillWignerDs(fsMap, evtData);
  
  Vector4<double> all4Vec(0.,0.,0.,0.);
  Vector4<double> mother4Vec(0.,0.,0.,0.);
  Vector4<double> daughter2_4Vec(0.,0.,0.,0.);

  //fill all4Vec
  std::map<std::string, Vector4<double>, pawian::Collection::PtrLess>::iterator itMap;
  for(itMap=fsMap.begin(); itMap!=fsMap.end(); ++itMap){
    all4Vec+=itMap->second;
  }

  //fill mother4Vec
  std::vector<Particle*>::iterator itP;
  for(itP = _finalStateParticles.begin(); itP != _finalStateParticles.end(); ++itP){
    itMap=fsMap.find((*itP)->name());
    mother4Vec+=itMap->second;
  }

  //fill mother4Vec
  for(itP=_finalStateParticlesDaughter2.begin(); itP!=_finalStateParticlesDaughter2.end(); ++itP){
    itMap=fsMap.find((*itP)->name());
    daughter2_4Vec+=itMap->second;
  }

  Vector4<double> daughter2HelMother(0.,0.,0.,0.);
  if(_hasMotherPart){
    daughter2HelMother=helicityVec(all4Vec, mother4Vec, daughter2_4Vec);
  }
  else{
    daughter2HelMother=daughter2_4Vec;
    daughter2HelMother.Boost(mother4Vec);
  }

  Spin spinMother=_motherJPCPtr->J;
  Spin spinDaughter1=_daughter1JPCPtr->J;
  Spin spinDaughter2=_daughter2JPCPtr->J;
  Spin lam12Max=spinDaughter1+spinDaughter2;
  if(lam12Max>spinMother) lam12Max=spinMother;

  Spin lamMotherMax=spinMother;
  if (!_hasMotherPart && spinMother>1) lamMotherMax=1;

 
  for (Spin lamMother=-lamMotherMax; lamMother<=lamMotherMax; ++lamMother){
    for (Spin lam12=-lam12Max; lam12<=lam12Max; ++lam12){
      double thePhi=0.;
      if(_hasMotherPart) thePhi=daughter2HelMother.Phi();
      evtData->WignerDsString[_wignerDKey][spinMother][lamMother][lam12]=Wigner_D(thePhi,daughter2HelMother.Theta(),0,spinMother,lamMother,lam12);
    }
  }
}

void IsobarDecay::print(std::ostream& os) const{
  os << "\nJPCLS 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";
  }
  
  if(!_daughter1IsStable){
    os << "with further decay:";
    _isoDecDaughter1->print(os);
  }
  if(!_daughter2IsStable){
    os << "with further decay:";
    _isoDecDaughter2->print(os);
  }  

  os << "\n";
}