Page MenuHomeHEPForge
Files F10881497

No OneTemporary
Actions

View Options

diff --git a/src/EvtGenModels/EvtPropSLPole.cpp b/src/EvtGenModels/EvtPropSLPole.cpp
index 7e1f02a..9204ff8 100644
--- a/src/EvtGenModels/EvtPropSLPole.cpp
+++ b/src/EvtGenModels/EvtPropSLPole.cpp
@@ -1,549 +1,549 @@
//--------------------------------------------------------------------------
//
// 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: EvtPropSLPole.cc
//
// Description: Routine to implement semileptonic decays according
// to light cone sum rules
//
// Modification history:
//
// DJL April 23, 1998 Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtPropSLPole.hh"
#include "EvtGenModels/EvtSLPoleFF.hh"
#include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh"
#include "EvtGenBase/EvtSemiLeptonicVectorAmp.hh"
#include "EvtGenBase/EvtSemiLeptonicTensorAmp.hh"
#include "EvtGenBase/EvtIntervalFlatPdf.hh"
#include "EvtGenBase/EvtScalarParticle.hh"
#include "EvtGenBase/EvtVectorParticle.hh"
#include "EvtGenBase/EvtTensorParticle.hh"
#include "EvtGenBase/EvtTwoBodyVertex.hh"
#include "EvtGenBase/EvtPropBreitWignerRel.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtAmpPdf.hh"
#include "EvtGenBase/EvtMassAmp.hh"
#include "EvtGenBase/EvtSpinType.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include <string>
std::string EvtPropSLPole::getName(){
return "PROPSLPOLE";
}
EvtDecayBase* EvtPropSLPole::clone(){
return new EvtPropSLPole;
}
void EvtPropSLPole::decay( EvtParticle *p ){
if(! _isProbMaxSet){
EvtId parnum,mesnum,lnum,nunum;
parnum = getParentId();
mesnum = getDaug(0);
lnum = getDaug(1);
nunum = getDaug(2);
double mymaxprob = calcMaxProb(parnum,mesnum,
lnum,nunum,SLPoleffmodel.get());
setProbMax(mymaxprob);
_isProbMaxSet = true;
}
double minKstMass = EvtPDL::getMinMass(p->getDaug(0)->getId());
double maxKstMass = EvtPDL::getMaxMass(p->getDaug(0)->getId());
EvtIntervalFlatPdf flat(minKstMass, maxKstMass);
EvtPdfGen<EvtPoint1D> gen(flat);
EvtPoint1D point = gen();
double massKst = point.value();
p->getDaug(0)->setMass(massKst);
p->initializePhaseSpace(getNDaug(),getDaugs());
// EvtVector4R p4meson = p->getDaug(0)->getP4();
calcamp->CalcAmp(p,_amp2,SLPoleffmodel.get());
EvtParticle *mesonPart = p->getDaug(0);
double meson_BWAmp = calBreitWigner(mesonPart, point);
int list[2];
list[0]=0; list[1]=0;
_amp2.vertex(0,0,_amp2.getAmp(list)*meson_BWAmp);
list[0]=0; list[1]=1;
_amp2.vertex(0,1,_amp2.getAmp(list)*meson_BWAmp);
list[0]=1; list[1]=0;
_amp2.vertex(1,0,_amp2.getAmp(list)*meson_BWAmp);
list[0]=1; list[1]=1;
_amp2.vertex(1,1,_amp2.getAmp(list)*meson_BWAmp);
list[0]=2; list[1]=0;
_amp2.vertex(2,0,_amp2.getAmp(list)*meson_BWAmp);
list[0]=2; list[1]=1;
_amp2.vertex(2,1,_amp2.getAmp(list)*meson_BWAmp);
return;
}
void EvtPropSLPole::initProbMax(){
_isProbMaxSet = false;
return;
}
void EvtPropSLPole::init(){
checkNDaug(3);
//We expect the parent to be a scalar
//and the daughters to be X lepton neutrino
checkSpinParent(EvtSpinType::SCALAR);
checkSpinDaughter(1,EvtSpinType::DIRAC);
checkSpinDaughter(2,EvtSpinType::NEUTRINO);
EvtSpinType::spintype mesontype=EvtPDL::getSpinType(getDaug(0));
SLPoleffmodel = std::make_unique<EvtSLPoleFF>(getNArg(),getArgs());
switch(mesontype) {
case EvtSpinType::SCALAR:
calcamp = std::make_unique<EvtSemiLeptonicScalarAmp>();
break;
case EvtSpinType::VECTOR:
calcamp = std::make_unique<EvtSemiLeptonicVectorAmp>();
break;
case EvtSpinType::TENSOR:
calcamp = std::make_unique<EvtSemiLeptonicTensorAmp>();
break;
default:
;
}
}
double EvtPropSLPole::calBreitWignerBasic(double maxMass){
if ( _width< 0.0001) return 1.0;
//its not flat - but generated according to a BW
double mMin=_massMin;
double mMax=_massMax;
if ( maxMass>-0.5 && maxMass< mMax) mMax=maxMass;
double massGood = EvtRandom::Flat(mMin, mMax);
double ampVal = sqrt(1.0/(pow(massGood-_mass, 2.0) + pow(_width, 2.0)/4.0));
return ampVal;
}
double EvtPropSLPole::calBreitWigner(EvtParticle *pmeson, EvtPoint1D point){
EvtId mesnum = pmeson->getId();
- double _mass = EvtPDL::getMeanMass(mesnum);
- double _width = EvtPDL::getWidth(mesnum);
- double _maxRange = EvtPDL::getMaxRange(mesnum);
+ _mass = EvtPDL::getMeanMass(mesnum);
+ _width = EvtPDL::getWidth(mesnum);
+ _maxRange = EvtPDL::getMaxRange(mesnum);
EvtSpinType::spintype mesontype=EvtPDL::getSpinType(mesnum);
_includeDecayFact=true;
_includeBirthFact=true;
_spin = mesontype;
_blatt = 3.0;
double maxdelta = 15.0*_width;
if ( _maxRange > 0.00001 ) {
_massMax=_mass+maxdelta;
_massMin=_mass-_maxRange;
}
else{
_massMax=_mass+maxdelta;
_massMin=_mass-15.0*_width;
}
_massMax=_mass+maxdelta;
if ( _massMin< 0. ) _massMin=0.;
EvtParticle* par=pmeson->getParent();
double maxMass=-1.;
if ( par != 0 ) {
if ( par->hasValidP4() ) maxMass=par->mass();
for ( size_t i=0;i<par->getNDaug();i++) {
EvtParticle *tDaug=par->getDaug(i);
if ( pmeson != tDaug )
maxMass-=EvtPDL::getMinMass(tDaug->getId());
}
}
EvtId *dauId=0;
double *dauMasses=0;
size_t nDaug = pmeson->getNDaug();
if ( nDaug > 0) {
dauId=new EvtId[nDaug];
dauMasses=new double[nDaug];
for (size_t j=0;j<nDaug;j++) {
dauId[j]=pmeson->getDaug(j)->getId();
dauMasses[j]=pmeson->getDaug(j)->mass();
}
}
EvtId *parId=0;
EvtId *othDaugId=0;
EvtParticle *tempPar=pmeson->getParent();
if (tempPar) {
parId=new EvtId(tempPar->getId());
if ( tempPar->getNDaug()==2 ) {
if ( tempPar->getDaug(0) == pmeson ) othDaugId=new EvtId(tempPar->getDaug(1)->getId());
else othDaugId=new EvtId(tempPar->getDaug(0)->getId());
}
}
if ( nDaug!=2) return calBreitWignerBasic(maxMass);
if ( _width< 0.00001) return 1.0;
//first figure out L - take the lowest allowed.
EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
int t1=EvtSpinType::getSpin2(spinD1);
int t2=EvtSpinType::getSpin2(spinD2);
int t3=EvtSpinType::getSpin2(_spin);
int Lmin=-10;
// allow for special cases.
if (Lmin<-1 ) {
//There are some things I don't know how to deal with
if ( t3>4) return calBreitWignerBasic(maxMass);
if ( t1>4) return calBreitWignerBasic(maxMass);
if ( t2>4) return calBreitWignerBasic(maxMass);
//figure the min and max allowwed "spins" for the daughters state
Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
if (Lmin<0) Lmin=0;
assert(Lmin==0||Lmin==2||Lmin==4);
}
//double massD1=EvtPDL::getMeanMass(dauId[0]);
//double massD2=EvtPDL::getMeanMass(dauId[1]);
double massD1=dauMasses[0];
double massD2=dauMasses[1];
// I'm not sure how to define the vertex factor here - so retreat to nonRel code.
if ( (massD1+massD2)> _mass ) return calBreitWignerBasic(maxMass);
//parent vertex factor not yet implemented
double massOthD=-10.;
double massParent=-10.;
int birthl=-10;
if ( othDaugId) {
EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
int tt1=EvtSpinType::getSpin2(spinOth);
int tt2=EvtSpinType::getSpin2(spinPar);
int tt3=EvtSpinType::getSpin2(_spin);
//figure the min and max allowwed "spins" for the daughters state
if ( (tt1<=4) && ( tt2<=4) ) {
birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
if (birthl<0) birthl=0;
massOthD=EvtPDL::getMeanMass(*othDaugId);
massParent=EvtPDL::getMeanMass(*parId);
}
}
double massM=_massMax;
if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
//special case... if the parent mass is _fixed_ we can do a little better
//and only for a two body decay as that seems to be where we have problems
// Define relativistic propagator amplitude
EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
vd.set_f(_blatt);
EvtPropBreitWignerRel bw(_mass,_width);
EvtMassAmp amp(bw,vd);
// if ( _fixMassForMax) amp.fixUpMassForMax();
// else std::cout << "problem problem\n";
if ( _includeDecayFact) {
amp.addDeathFact();
amp.addDeathFactFF();
}
if ( massParent>-1.) {
if ( _includeBirthFact ) {
EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
amp.setBirthVtx(vb);
amp.addBirthFact();
amp.addBirthFactFF();
}
}
EvtAmpPdf<EvtPoint1D> pdf(amp);
double ampVal = sqrt(pdf.evaluate(point));
if ( parId) delete parId;
if ( othDaugId) delete othDaugId;
if ( dauId) delete [] dauId;
if ( dauMasses) delete [] dauMasses;
return ampVal;
}
double EvtPropSLPole::calcMaxProb( EvtId parent, EvtId meson,
EvtId lepton, EvtId nudaug,
EvtSemiLeptonicFF *FormFactors ) {
//This routine takes the arguements parent, meson, and lepton
//number, and a form factor model, and returns a maximum
//probability for this semileptonic form factor model. A
//brute force method is used. The 2D cos theta lepton and
//q2 phase space is probed.
//Start by declaring a particle at rest.
//It only makes sense to have a scalar parent. For now.
//This should be generalized later.
EvtScalarParticle *scalar_part;
EvtParticle *root_part;
scalar_part=new EvtScalarParticle;
//cludge to avoid generating random numbers!
scalar_part->noLifeTime();
EvtVector4R p_init;
p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
scalar_part->init(parent,p_init);
root_part=(EvtParticle *)scalar_part;
// root_part->set_type(EvtSpinType::SCALAR);
root_part->setDiagonalSpinDensity();
EvtParticle *daughter, *lep, *trino;
EvtAmp amp;
EvtId listdaug[3];
listdaug[0] = meson;
listdaug[1] = lepton;
listdaug[2] = nudaug;
amp.init(parent,3,listdaug);
root_part->makeDaughters(3,listdaug);
daughter=root_part->getDaug(0);
lep=root_part->getDaug(1);
trino=root_part->getDaug(2);
EvtDecayBase *decayer;
decayer = EvtDecayTable::getInstance()->getDecayFunc(daughter);
if ( decayer ) {
daughter->makeDaughters(decayer->nRealDaughters(),decayer->getDaugs());
for(int ii=0; ii<decayer->nRealDaughters(); ii++){
daughter->getDaug(ii)->setMass(EvtPDL::getMeanMass(daughter->getDaug(ii)->getId()));
}
}
//cludge to avoid generating random numbers!
daughter->noLifeTime();
lep->noLifeTime();
trino->noLifeTime();
//Initial particle is unpolarized, well it is a scalar so it is
//trivial
EvtSpinDensity rho;
rho.setDiag(root_part->getSpinStates());
double mass[3];
double m = root_part->mass();
EvtVector4R p4meson, p4lepton, p4nu, p4w;
double q2max;
double q2, elepton, plepton;
int i,j;
double erho,prho,costl;
double maxfoundprob = 0.0;
double prob = -10.0;
int massiter;
for (massiter=0;massiter<3;massiter++){
mass[0] = EvtPDL::getMeanMass(meson);
mass[1] = EvtPDL::getMeanMass(lepton);
mass[2] = EvtPDL::getMeanMass(nudaug);
if ( massiter==1 ) {
mass[0] = EvtPDL::getMinMass(meson);
}
if ( massiter==2 ) {
mass[0] = EvtPDL::getMaxMass(meson);
if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
}
q2max = (m-mass[0])*(m-mass[0]);
//loop over q2
for (i=0;i<25;i++) {
q2 = ((i+0.5)*q2max)/25.0;
erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
prho = sqrt(erho*erho-mass[0]*mass[0]);
p4meson.set(erho,0.0,0.0,-1.0*prho);
p4w.set(m-erho,0.0,0.0,prho);
//This is in the W rest frame
elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
double probctl[3];
for (j=0;j<3;j++) {
costl = 0.99*(j - 1.0);
//These are in the W rest frame. Need to boost out into
//the B frame.
p4lepton.set(elepton,0.0,
plepton*sqrt(1.0-costl*costl),plepton*costl);
p4nu.set(plepton,0.0,
-1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
p4lepton=boostTo(p4lepton,boost);
p4nu=boostTo(p4nu,boost);
//Now initialize the daughters...
daughter->init(meson,p4meson);
lep->init(lepton,p4lepton);
trino->init(nudaug,p4nu);
calcamp->CalcAmp(root_part,amp,FormFactors);
EvtPoint1D *point = new EvtPoint1D(mass[0]);
double meson_BWAmp = calBreitWigner(daughter, *point);
int list[2];
list[0]=0; list[1]=0;
amp.vertex(0,0,amp.getAmp(list)*meson_BWAmp);
list[0]=0; list[1]=1;
amp.vertex(0,1,amp.getAmp(list)*meson_BWAmp);
list[0]=1; list[1]=0;
amp.vertex(1,0,amp.getAmp(list)*meson_BWAmp);
list[0]=1; list[1]=1;
amp.vertex(1,1,amp.getAmp(list)*meson_BWAmp);
list[0]=2; list[1]=0;
amp.vertex(2,0,amp.getAmp(list)*meson_BWAmp);
list[0]=2; list[1]=1;
amp.vertex(2,1,amp.getAmp(list)*meson_BWAmp);
//Now find the probability at this q2 and cos theta lepton point
//and compare to maxfoundprob.
//Do a little magic to get the probability!!
prob = rho.normalizedProb(amp.getSpinDensity());
probctl[j]=prob;
}
//probclt contains prob at ctl=-1,0,1.
//prob=a+b*ctl+c*ctl^2
double a=probctl[1];
double b=0.5*(probctl[2]-probctl[0]);
double c=0.5*(probctl[2]+probctl[0])-probctl[1];
prob=probctl[0];
if (probctl[1]>prob) prob=probctl[1];
if (probctl[2]>prob) prob=probctl[2];
if (fabs(c)>1e-20){
double ctlx=-0.5*b/c;
if (fabs(ctlx)<1.0){
double probtmp=a+b*ctlx+c*ctlx*ctlx;
if (probtmp>prob) prob=probtmp;
}
}
//EvtGenReport(EVTGEN_DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" "
// << probctl[0]<<" "
// << probctl[1]<<" "
// << probctl[2]<<endl;
if ( prob > maxfoundprob ) {
maxfoundprob = prob;
}
}
if ( EvtPDL::getWidth(meson) <= 0.0 ) {
//if the particle is narrow dont bother with changing the mass.
massiter = 4;
}
}
root_part->deleteTree();
maxfoundprob *=1.1;
return maxfoundprob;
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:27 AM (1 d, 17 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983032
Default Alt Text
(15 KB)

Event Timeline