PhotosParticle.cxx

#include <vector>
#include <math.h>
#include "PhotosParticle.h"
#include "Log.h"
using std::vector;

namespace Photospp
{

bool PhotosParticle::hasDaughters()
{
        if(getDaughters().size()==0) return false;
        else                         return true;
}

PhotosParticle * PhotosParticle::findLastSelf()
{
        vector<PhotosParticle*> daughters = getDaughters();
        vector<PhotosParticle*>::iterator pcl_itr = daughters.begin();

        //get all daughters and look for stable with same pgd id
        for(;pcl_itr != daughters.end();pcl_itr++)
        {
                if((*pcl_itr)->getPdgID()==this->getPdgID())
                return (*pcl_itr)->findLastSelf();
        }

        return this;
}

vector<PhotosParticle*> PhotosParticle::findProductionMothers()
{
        vector<PhotosParticle*> mothers = getMothers();
        vector<PhotosParticle*>::iterator pcl_itr = mothers.begin();

        //get all mothers and check none have pdg id of this one
        for(;pcl_itr != mothers.end();pcl_itr++)
        {
                if((*pcl_itr)->getPdgID()==this->getPdgID())
                return (*pcl_itr)->findProductionMothers();
        }
        return mothers;
}

vector<PhotosParticle *> PhotosParticle::getDecayTree()
{
        vector<PhotosParticle *> particles;
        particles.push_back(this);
        vector<PhotosParticle *> daughters = getDaughters();
        for(int i=0;i<(int)daughters.size();i++)
        {
                // Check if we are the first mother of each daughters
                // If not - skip this daughter
                PhotosParticle *p = daughters.at(i);
                vector<PhotosParticle *> mothers = p->getMothers();
                if(mothers.size()>1 && mothers.at(0)->getBarcode()!=getBarcode()) continue;
                vector<PhotosParticle *> tree = p->getDecayTree();
                particles.insert(particles.end(),tree.begin(),tree.end());
        }
        return particles;
}

void PhotosParticle::boostDaughtersFromRestFrame(PhotosParticle * tau_momentum)
{
        if(!hasDaughters()) //if there are no daughters
        return;

        // get all daughters, granddaughters, etc. then rotate and boost them
        vector<PhotosParticle*> list = getAllDecayProducts();
        vector<PhotosParticle*>::iterator pcl_itr = list.begin();

        for(;pcl_itr != list.end();pcl_itr++)
        {
                (*pcl_itr)->boostFromRestFrame(tau_momentum);
        }

        //checkMomentumConservation();
}

void PhotosParticle::boostDaughtersToRestFrame(PhotosParticle * tau_momentum)
{
        if(!hasDaughters()) //if there are no daughters
        return;

        // get all daughters, granddaughters, etc. then rotate and boost them
        vector<PhotosParticle*> list = getAllDecayProducts();
        vector<PhotosParticle*>::iterator pcl_itr = list.begin();

        for(;pcl_itr != list.end();pcl_itr++)
        {
                (*pcl_itr)->boostToRestFrame(tau_momentum);
        }

        //checkMomentumConservation();
}


void PhotosParticle::boostToRestFrame(PhotosParticle * tau_momentum)
{
        double theta = tau_momentum->getRotationAngle(Y_AXIS);
        tau_momentum->rotate(Y_AXIS,theta);

        double phi = tau_momentum->getRotationAngle(X_AXIS);
        tau_momentum->rotate(Y_AXIS,-theta);

        //Now rotate coordinates to get boost in Z direction.
        rotate(Y_AXIS,theta);
        rotate(X_AXIS,phi);
        boostAlongZ(-1*tau_momentum->getP(),tau_momentum->getE());
        rotate(X_AXIS,-phi);
        rotate(Y_AXIS,-theta);
}

void PhotosParticle::boostFromRestFrame(PhotosParticle * tau_momentum)
{
        //get the rotation angles
        //and boost z

        double theta = tau_momentum->getRotationAngle(Y_AXIS);
        tau_momentum->rotate(Y_AXIS,theta);

        double phi = tau_momentum->getRotationAngle(X_AXIS);
        tau_momentum->rotate(Y_AXIS,-theta);

        //Now rotate coordinates to get boost in Z direction.
        rotate(Y_AXIS,theta);
        rotate(X_AXIS,phi);
        boostAlongZ(tau_momentum->getP(),tau_momentum->getE());
        rotate(X_AXIS,-phi);
        rotate(Y_AXIS,-theta);
}

double PhotosParticle::getRotationAngle(int axis, int second_axis)
{
        if(getP(second_axis)==0)
        {
                if(getP(axis)>0) return -M_PI/2.0;
                else             return  M_PI/2.0;
        }
        if(getP(second_axis)>0) return     -atan(getP(axis)/getP(second_axis));
        else                    return M_PI-atan(getP(axis)/getP(second_axis));

}

void PhotosParticle::boostAlongZ(double boost_pz, double boost_e)
{
        // Boost along the Z axis
        double m_tau=sqrt(boost_e*boost_e-boost_pz*boost_pz);

        double p=getPz();
        double e=getE();

        setPz((boost_e*p + boost_pz*e)/m_tau);
        setE((boost_pz*p + boost_e*e )/m_tau);
}

void PhotosParticle::rotate(int axis,double theta, int second_axis)
{
        double temp_px=getP(axis);
        double temp_pz=getP(second_axis);
        setP(axis,cos(theta)*temp_px + sin(theta)*temp_pz);
        setP(second_axis,-sin(theta)*temp_px + cos(theta)*temp_pz);
}

void PhotosParticle::rotateDaughters(int axis,double theta, int second_axis)
{
        if(!hasDaughters()) //if there are no daughters
        return;

        vector<PhotosParticle*> daughters = getDaughters();
        vector<PhotosParticle*>::iterator pcl_itr = daughters.begin();

        //get all daughters then rotate and boost them.
        for(;pcl_itr != daughters.end();pcl_itr++)
        {
                (*pcl_itr)->rotate(axis,theta,second_axis);
                (*pcl_itr)->rotateDaughters(axis,theta,second_axis);
        }
        //checkMomentumConservation();
}

double PhotosParticle::getVirtuality()
{
        double e_sq=getE()*getE();
        double p_sq=getP()*getP();

        if(e_sq>p_sq) return    sqrt(e_sq-p_sq);
        else          return -1*sqrt(p_sq-e_sq); //if it's negative
}

double PhotosParticle::getP()
{
        return sqrt(getPx()*getPx()+getPy()*getPy()+getPz()*getPz());
}

double PhotosParticle::getP(int axis)
{
        if(axis==X_AXIS) return getPx();
        if(axis==Y_AXIS) return getPy();
        if(axis==Z_AXIS) return getPz();
        return 0;
}

void PhotosParticle::setP(int axis, double p_component)
{
        if(axis==X_AXIS) setPx(p_component);
        if(axis==Y_AXIS) setPy(p_component);
        if(axis==Z_AXIS) setPz(p_component);
}

} // namespace Photospp