---

GParticle.cxx

Go to the documentation of this file.

// $Id: GParticle.cxx,v 1.4 2000/05/05 04:38:33 burnett Exp $
//
// This file is part of Gismo 2
//
//
// This is the implementation of the GParticle class.
// and Random::theGenerator


#include "gismo/GParticle.h"

#include "DecayEntry.h"
#include "gismo/PDT.h"
#include "gismo/Units.h"
#include "facilities/error.h"
#include "Random.h"
#include <algorithm>
#include <strstream>
inline static double sqr(double x){return x*x;}

GParticle::GParticle(const PData * data, float mass)
 : pData(data)
 , m_parent(0)
 , m_first_child(0)
 , m_next_sibling(0)
 , _time(0)
{
    setMass(mass);
}

GParticle::GParticle(long id,float mass)
 : m_parent(0)
 , m_first_child(0)
 , m_next_sibling(0)
 , _time(0)
{
    pData = _thePDT->lookUp(id);
    if (pData == 0)
        FATAL("cannot create particle!");
    setMass(mass);
}

GParticle::GParticle(const char* name, float mass)
 : m_parent(0)
 , m_first_child(0)
 , m_next_sibling(0)
 , _time(0)
{
    pData = _thePDT->lookUp(name);
    if (pData == 0)
        FATAL("cannot create particle!");
    setMass(mass);
}
// copy constructor -- do not copy children, make original "parent"
GParticle::GParticle(const GParticle& p)
 : HepLorentzVector((HepLorentzVector&)p)
 , pData(p.pData)
 , m_parent(&p)
 , m_first_child(0)
 , m_next_sibling(0)
 , _mass(p._mass)
 , _r(p._r)
 , _time(p._time)
{}

// set the mass: this is private, only done from constructor with momentum zero
void GParticle::setMass(float mass)
{  // negative is flag to generate with breitWigner
        _mass = mass<0 ? Random::breitWigner(pData->mass(), pData->width(),
                                                        pData->widthCut())
                                         : mass ;
        setE(_mass);
}

// set the mometum, adjusting the energy
GParticle&
GParticle::setMomentum(const Vector & p)
{
  setVect(p);
  setE(sqrt(vect().mag2() + _mass*_mass));
  return *this;
}

// set the energy, scaling the momentum to preserve the mass
GParticle&
GParticle::setEnergy(float e)
{
  e = std::max(e,_mass);
  float scale = sqrt((e*e - _mass*_mass)/vect().mag2());
  setE( e );
  setVect(vect()*scale);
  return *this;
}


float GParticle::GenerateTimeOfDecay()
{
    return  (pData->lifeTime() > 0 && pData->lifeTime()< 1e8)
          ? Random::exponential(pData->lifeTime())
          : pData->lifeTime();
}

GParticle::~GParticle()
{
    removeChildren();
}

void GParticle::setParent(GParticle* p){m_parent = p;}

void GParticle::decayOnce()
// Summary -------------------------------------------------------------
//
//    Decay the particle, creating daughters and attaching them
//
// Parameters
//
//   none
//
// Remarks
//
// End -----------------------------------------------------------------
{
    const DecayEntry* de = pData->decayInfo();
    if (de == 0) return;

    int anti = pData->idCode() < 0;

    // choose the decay

    float r = Random::flat(pData->sumBR());
    while ((r -= de->branchingFraction)>0 && de->next != 0) de=de->next;

    // create the secondary particles at rest for the decayer
    if (de->numChildren > 0)    {

        float Q = -1.0; // the failure code
        for(unsigned ntries = 0; ntries < MAX_TRIES; ntries++)  {
            Q = mass();

            for( int i=0; i < de->numChildren; i++ ){

                PData *childData = de->children[i];
                if( anti ) childData = childData->antiParticle;
                GParticle* child = addChild(childData); // this is virtual to allow addition of
                                                        // supclass
                Q -= child->mass();
            }

            // check that it is possible: maybe need to recreate?
            if (Q < 0.){
                std::ostrstream serr;
                serr<< "Decay kinematically forbidden:" <<de->name()<<'\n';
                serr<< "trying again\n";
                WARNING(serr.str());
                removeChildren();;
            }
            else
            {
                break;
            }
        }
        if (Q<0.)FATAL( "Failed, aborting...");
    }

    // boost to rest frame: first save current velocity;
    Vector beta(px()/e(),py()/e(),pz()/e());
    Vector oldMomentum( momentum() );  //save original (possible round-off when restore)
    setMomentum(Vector(0,0,0));

    // apply matrix elements, etc.
    if (! de->decay( this))     {
        // cannot decay: must reverse
        WARNING("*** Error in decay\n");
        removeChildren();
        return;
    }

    // boost secondaries from rest frame to lab frame

    if (beta.mag()!=0.) boost(beta);
    setMomentum( oldMomentum );  //restore

}

void
GParticle::removeChildren()
{
    GParticle* next_child = m_first_child;
    while( next_child) {
        GParticle* tc = next_child->m_next_sibling;
        delete next_child;
        next_child = tc;
    }
    m_first_child = 0;
}

GParticle*
GParticle::addChild(GParticle * child)
// used by a client to set a single "child"
{
    if( ! m_first_child) m_first_child = child;
    else {
        GParticle* this_child = m_first_child;
        while(this_child-> m_next_sibling !=0 ) this_child= this_child->m_next_sibling;
        this_child->m_next_sibling = child;
    }
    child->setParent(this);
    child->setPosition(_r);
    child->setT(_time);
    return child;
}
// this is virtual, so that a sub class will add its own kind
GParticle *
GParticle::addChild(const PData* childData, float mass)
{
    return addChild(new GParticle(childData, mass));
}
int
GParticle::numChildren()const
{
    int n = 0;
    GParticle* the_child = m_first_child;
    while(the_child != 0){
        n++; the_child = the_child->m_next_sibling;
    }
    return n;

}

const GParticle*
GParticle::parent()const{return m_parent;}

GParticle*
GParticle::child(int n)const
{
    GParticle* the_child = m_first_child;
    while(n-- && the_child !=0 ) the_child = the_child->m_next_sibling;
    return the_child;
}


void
GParticle::boost(const Vector& b)
{
   HepLorentzVector::boost(b);

   for(int i=0; i<numChildren(); i++)
      child(i)->boost(b);
}


static int level;
void GParticle::printAll(std::ostream& cout) const
{
        if( !parent() ) level=0;
        printOn(cout);    // this is virtual!
        level++;
        for(int i=0; i<numChildren(); i++)
                child(i)->printAll(cout);
        level--;
}
void GParticle::printOn(std::ostream& cout) const
{
        cout << '\n';
        for(int i=0;i<level;i++)
                cout << "  ";
        cout << name() << ' ';
        cout << *(HepLorentzVector*)this;
        cout << " at: " << _r;
}

void
GParticle::scatterBy(double theta, double phi)
{
    double sinthe = sin(theta),
        sinps2=sqr(px())+sqr(py()),
        magnitude = sinps2 + sqr(pz());
    if( magnitude==0 ) return;
    sinps2 /=magnitude;
    magnitude = sqrt(magnitude);
    double cosTheta = cos(theta);

    //   if sinps2 is small, no rotation is needed

    if (sinps2 < 1.0e-20)
    {
        //small polar angle case
        if(pz() < 0.) cosTheta = -cosTheta;
        setVect(Vector(sinthe*cos(phi), sinthe*sin(phi), cosTheta ));
    }
    else
    {
        //large polar angle case
         float a=px()/magnitude,
                 b=py()/magnitude,
                 c=pz()/magnitude,
                 sinpsi=sqrt(sinps2),
                 us=sinthe*cos(phi),
                 vs=sinthe*sin(phi),
                 sindel=b/sinpsi,
                 cosdel=a/sinpsi;
         setVect(Vector(c*cosdel*us - sindel*vs + a*cosTheta ,
                        c*sindel*us + cosdel*vs + b*cosTheta ,
                        -sinpsi*us              + c*cosTheta ));
     }
     setVect(vect()*magnitude);
}
// ----------------------------------------------------------
PDT*  GParticle::_thePDT = PDT::createDefault();  // pdt used to create particles
const unsigned GParticle::MAX_TRIES=100;
//don't really want these
HepBoolean GParticle::hasPostionInfo() const{return true;}
float GParticle::timeOfDecay()const{return -1;}

---