UsermDIY.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/BesCopyright
//      Copyright (A) 2006      Ping Rong-Gang
//
// Module: Modified DIY model, see EvtDeay.cc 
//
// Description: Routine to sample the decays using the amplitude specified by users
// This model allow use to specify the intermediate states
//
// Modification history:
//
// Ping R.-G.    4-25, 2010       Module created
//
//------------------------------------------------------------------------
#ifndef EVTDECAY_H
#define EVTDECAY_H
 
#include "EvtDecay.h"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtVector3R.hh"
#include "EvtGenBase/EvtHelSys.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGen.hh"
#include "EvtGenBase/EvtRandomEngine.hh"
#include "EvtGenBase/EvtDecayTable.hh"
 
#include <stdlib.h>
#include <string.h>
using std::cout;
using std::endl;
 
double EvtDecay::CalAmpsMDIY(EvtParticle* part ){
//---- read out the momentum for each particle, then using them to calculate the amplitude squared, for example:
/*********** using the following decay cards ****************
Decay J/psi
1 gamma eta_c  PHSP;
Enddecay
 
Decay eta_c
1 phi phi PHSP;
Enddecay
 
Decay phi
1 K+ K- PHSP;
Enddecay
 
End
*////////
//----------------------------------
    EvtVector4R pjsi,petac,pphi1,pphi2,pkp1,pkp2;
    pjsi = part->getP4();
    petac= part->getDaug(1)->getP4();  //etac momentum 
    pphi1= part->getDaug(1)->getDaug(0)->getP4();  //phi1 momentum
    pphi2= part->getDaug(1)->getDaug(1)->getP4();  //phi2 momentum
 
    pkp1 = part->getDaug(1)->getDaug(0)->getDaug(0)->getP4();  //K+ from phi1 
    pkp2 = part->getDaug(1)->getDaug(1)->getDaug(0)->getP4();  //K+ from phi2 
 
    EvtHelSys angles0(pjsi,petac);
    double theta0 = angles0.getHelAng(1);
    double phi0   = angles0.getHelAng(2);
 
    EvtHelSys angles2(pphi1,pkp1);
    double theta2 = angles2.getHelAng(1);
    double phi2   = angles2.getHelAng(2);
 
    EvtHelSys angles3(pphi2,pkp2);
    double theta3 = angles3.getHelAng(1);
    double phi3   = angles3.getHelAng(2);
 
    double amps=(3+cos(2*theta0)) * sin(theta2)*sin(theta2) * sin(theta3)*sin(theta3) * sin(phi2+phi3)*sin(phi2+phi3);
 
 
///////////======== don't touch follows =======================
    if(amps <=0){
    report(INFO,"EvtGen") << "Amplitude square of modified DIY should be positive, but found to be equal "<<amps<<endl;
    abort();
    } else {
    return amps;
    }
}
 
 
#endif