EvtPi0Dalitz.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: EvtPi0Dalitz.cc
//
// Description: pi0 -> e+ e- gamma 
//
// Modification history:
//
//    DJL/RYD    June 30, 1998          Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include <fstream>
#include <stdio.h>
#include <string>
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtPi0Dalitz.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtTensor4C.hh"
using std::fstream;
 
EvtPi0Dalitz::~EvtPi0Dalitz() {}
 
void EvtPi0Dalitz::getName(std::string& model_name){
 
    model_name="PI0_DALITZ";
 
}
 
EvtDecayBase* EvtPi0Dalitz::clone(){
 
  return new EvtPi0Dalitz;
 
}
 
 
void EvtPi0Dalitz::initProbMax(){
 
  setProbMax(3.5);
 
}
 
 
void EvtPi0Dalitz::init(){
 
  // check that there are 0 arguments
  checkNArg(0);
  checkNDaug(3);
 
    
  checkSpinParent(EvtSpinType::SCALAR);
 
  checkSpinDaughter(0,EvtSpinType::DIRAC);
  checkSpinDaughter(1,EvtSpinType::DIRAC);
  checkSpinDaughter(2,EvtSpinType::PHOTON);
 
}
 
 
void EvtPi0Dalitz::decay( EvtParticle *p){
 
  EvtParticle *ep, *em, *gamma;
  setWeight(p->initializePhaseSpace(getNDaug(),getDaugs(),0.00000002,0,1));
 
  ep=p->getDaug(0);
  em=p->getDaug(1);
  gamma=p->getDaug(2);
 
  // the next four lines generates events with a weight such that
  // the efficiency for selecting them is good. The parameter below of
  // 0.1 is the size of the peak at low q^2 (in arbitrary units).
  // The value of 0.1 is appropriate for muons. 
  // when you use this remember to remove the cut on q^2!
   
 
  //ep em invariant mass^2
  double m2=(ep->getP4()+em->getP4()).mass2();
  EvtVector4R q=ep->getP4()+em->getP4();
  //Just use the prob summed over spins...
 
  EvtTensor4C w,v;
 
  v=2.0*(gamma->getP4()*q)*directProd(q,gamma->getP4()) 
    - (gamma->getP4()*q)*(gamma->getP4()*q)*EvtTensor4C::g()
    -m2*directProd(gamma->getP4(),gamma->getP4());
 
  w=4.0*( directProd(ep->getP4(),em->getP4()) + directProd(em->getP4(),ep->getP4())
       -EvtTensor4C::g()*(ep->getP4()*em->getP4()-ep->getP4().mass2()));
 
  double prob=(real(cont(v,w)))/(m2*m2);
  prob *=(1.0/( (0.768*0.768-m2)*(0.768*0.768-m2)
           +0.768*0.768*0.151*0.151));
 
  //  report(INFO,"EvtGen") << "prob is "<<prob<<endl;
  setProb(prob);
 
  return;
}