EvtVectorIsr.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtVectorIsr.cc
//
// Description: Model to simulate e+e- -> vector + gamma
//              This is implemented as a decay of the VPHO.
//
// Modification history:
//
//    Joe Izen        Dec  16, 2002         Fix cos_theta distribution - prevents boom at cos_theta=+/-1 
//    RYD/Adriano     June 16, 1998         Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenModels/EvtVectorIsr.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtConst.hh"
#include <string>
#include "EvtGenBase/EvtVector4C.hh"
 
EvtVectorIsr::~EvtVectorIsr() {}
 
void EvtVectorIsr::getName(std::string& model_name){
 
  model_name="VECTORISR";     
 
}
 
EvtDecayBase* EvtVectorIsr::clone(){
 
  return new EvtVectorIsr;
 
}
 
void EvtVectorIsr::init(){
 
  // check that there are 2 arguments
  checkNArg(2);
  checkNDaug(2);
 
  checkSpinParent(EvtSpinType::VECTOR);
  checkSpinDaughter(0,EvtSpinType::VECTOR);
  checkSpinDaughter(1,EvtSpinType::PHOTON);
 
  //copy the arguments into eaiser to remember names
 
  csfrmn=getArg(0);
  csbkmn=getArg(1);
 
 
}
 
void EvtVectorIsr::initProbMax(){
 
  noProbMax();
 
}
 
void EvtVectorIsr::decay( EvtParticle *p ){
 
  EvtParticle *phi;
  EvtParticle *gamma;
 
  //the elctron mass
  double electMass=EvtPDL::getMeanMass(EvtPDL::getId("e-"));
 
 
  //get pointers to the daughters set
  //get masses/initial phase space - will overwrite the
  //p4s below to get the kinematic distributions correct
  p->initializePhaseSpace(getNDaug(),getDaugs());
  phi=p->getDaug(0);
  gamma=p->getDaug(1);
 
  double wcm=p->mass();
  double beta=2.*electMass/wcm; //electMass/Ebeam = betagamma
  beta=sqrt(1. - beta*beta);   //sqrt (1 - (m/ebeam)**2)
 
  //gamma momentum in the parents restframe
  double pg=(wcm*wcm-phi->mass()*phi->mass())/(2*wcm);
 
//    //generate kinematics according to Bonneau-Martin article
//    //Nucl. Phys. B27 (1971) 381-397
 
//    double y=pow((1+csmn)/(1-csmn),EvtRandom::Flat(0.0,1.0));
//    double cs=kcs*(y-1)/(y+1);
 
  // For backward compatibility with .dec files before SP5, the backward cos limit for
  //the ISR photon is actually given as *minus* the actual limit. Sorry, this wouldn't be
  //my choice.  -Joe
 
  double ymax=log((1.+beta*csfrmn)/(1.-beta*csfrmn));
  double ymin=log((1.-beta*csbkmn)/(1.+beta*csbkmn));
 
  // photon theta distributed as  2*beta/(1-beta**2*cos(theta)**2)
  double y=(ymax-ymin)*EvtRandom::Flat(0.0,1.0) + ymin;
  double cs=exp(y);
  cs=(cs - 1.)/(cs + 1.)/beta;
  double sn=sqrt(1-cs*cs);
 
  double fi=EvtRandom::Flat(EvtConst::twoPi);
 
  //four-vector for the phi
  EvtVector4R p4phi(sqrt(phi->mass()*phi->mass()+pg*pg),pg*sn*cos(fi),
         pg*sn*sin(fi),pg*cs);
 
  EvtVector4R p4gamma(pg,-p4phi.get(1),-p4phi.get(2),-p4phi.get(3));
 
  //save momenta for particles
  phi->init( getDaug(0),p4phi);
  gamma->init( getDaug(1),p4gamma);
 
  //try setting the spin density matrix of the phi
 
 
  //first get the polarization vectors of the gamma and phi
 
  EvtVector4C phi0=phi->epsParent(0);
  EvtVector4C phi1=phi->epsParent(1);
  EvtVector4C phi2=phi->epsParent(2);
 
  EvtVector4C gamma0=gamma->epsParentPhoton(0);
  EvtVector4C gamma1=gamma->epsParentPhoton(1);
 
  EvtComplex r1p=phi0*gamma0;
  EvtComplex r2p=phi1*gamma0;
  EvtComplex r3p=phi2*gamma0;
 
 
  EvtComplex r1m=phi0*gamma1;
  EvtComplex r2m=phi1*gamma1;
  EvtComplex r3m=phi2*gamma1;
 
  EvtComplex rho33=r3p*conj(r3p)+r3m*conj(r3m);
  EvtComplex rho22=r2p*conj(r2p)+r2m*conj(r2m);
  EvtComplex rho11=r1p*conj(r1p)+r1m*conj(r1m);
 
  EvtComplex rho13=r3p*conj(r1p)+r3m*conj(r1m);
  EvtComplex rho12=r2p*conj(r1p)+r2m*conj(r1m);
  EvtComplex rho23=r3p*conj(r2p)+r3m*conj(r2m);
 
  EvtComplex rho31=conj(rho13);
  EvtComplex rho32=conj(rho23);
  EvtComplex rho21=conj(rho12);
 
 
  EvtSpinDensity rho;
  rho.SetDim(3);
 
  rho.Set(0,0,rho11);
  rho.Set(0,1,rho12);
  rho.Set(0,2,rho13);
  rho.Set(1,0,rho21);
  rho.Set(1,1,rho22);
  rho.Set(1,2,rho23);
  rho.Set(2,0,rho31);
  rho.Set(2,1,rho32);
  rho.Set(2,2,rho33);
 
 
  setDaughterSpinDensity(0);
  phi->setSpinDensityForward(rho);
 
 
  return ;
}