1 //--------------------------------------------------------------------------
4 // This software is part of the EvtGen package developed jointly
5 // for the BaBar and CLEO collaborations. If you use all or part
6 // of it, please give an appropriate acknowledgement.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtSemiLeptonicAmp.cc
13 // Description: Base class for semileptonic decays
15 // Modification history:
17 // DJL April 17,1998 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
22 #include "EvtGenBase/EvtPatches.hh"
23 #include "EvtGenBase/EvtParticle.hh"
24 #include "EvtGenBase/EvtGenKine.hh"
25 #include "EvtGenBase/EvtPDL.hh"
26 #include "EvtGenBase/EvtReport.hh"
27 #include "EvtGenBase/EvtVector4C.hh"
28 #include "EvtGenBase/EvtTensor4C.hh"
29 #include "EvtGenBase/EvtDiracSpinor.hh"
30 #include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh"
31 #include "EvtGenBase/EvtId.hh"
32 #include "EvtGenBase/EvtAmp.hh"
33 #include "EvtGenBase/EvtScalarParticle.hh"
34 #include "EvtGenBase/EvtVectorParticle.hh"
35 #include "EvtGenBase/EvtTensorParticle.hh"
37 double EvtSemiLeptonicAmp::CalcMaxProb( EvtId parent, EvtId meson,
38 EvtId lepton, EvtId nudaug,
39 EvtSemiLeptonicFF *FormFactors ) {
41 //This routine takes the arguements parent, meson, and lepton
42 //number, and a form factor model, and returns a maximum
43 //probability for this semileptonic form factor model. A
44 //brute force method is used. The 2D cos theta lepton and
45 //q2 phase space is probed.
47 //Start by declaring a particle at rest.
49 //It only makes sense to have a scalar parent. For now.
50 //This should be generalized later.
52 EvtScalarParticle *scalar_part;
53 EvtParticle *root_part;
55 scalar_part=new EvtScalarParticle;
57 //cludge to avoid generating random numbers!
58 scalar_part->noLifeTime();
62 p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
63 scalar_part->init(parent,p_init);
64 root_part=(EvtParticle *)scalar_part;
65 root_part->setDiagonalSpinDensity();
67 EvtParticle *daughter, *lep, *trino;
76 amp.init(parent,3,listdaug);
78 root_part->makeDaughters(3,listdaug);
79 daughter=root_part->getDaug(0);
80 lep=root_part->getDaug(1);
81 trino=root_part->getDaug(2);
83 //cludge to avoid generating random numbers!
84 daughter->noLifeTime();
89 //Initial particle is unpolarized, well it is a scalar so it is
92 rho.setDiag(root_part->getSpinStates());
96 double m = root_part->mass();
98 EvtVector4R p4meson, p4lepton, p4nu, p4w;
101 double q2, elepton, plepton;
103 double erho,prho,costl;
105 double maxfoundprob = 0.0;
109 for (massiter=0;massiter<3;massiter++){
111 mass[0] = EvtPDL::getMeanMass(meson);
112 mass[1] = EvtPDL::getMeanMass(lepton);
113 mass[2] = EvtPDL::getMeanMass(nudaug);
115 mass[0] = EvtPDL::getMinMass(meson);
118 mass[0] = EvtPDL::getMaxMass(meson);
119 if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
122 q2max = (m-mass[0])*(m-mass[0]);
127 q2 = ((i+0.5)*q2max)/25.0;
129 erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
131 prho = sqrt(erho*erho-mass[0]*mass[0]);
133 p4meson.set(erho,0.0,0.0,-1.0*prho);
134 p4w.set(m-erho,0.0,0.0,prho);
136 //This is in the W rest frame
137 elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
138 plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
144 costl = 0.99*(j - 1.0);
146 //These are in the W rest frame. Need to boost out into
148 p4lepton.set(elepton,0.0,
149 plepton*sqrt(1.0-costl*costl),plepton*costl);
150 p4nu.set(plepton,0.0,
151 -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
153 EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
154 p4lepton=boostTo(p4lepton,boost);
155 p4nu=boostTo(p4nu,boost);
157 //Now initialize the daughters...
159 daughter->init(meson,p4meson);
160 lep->init(lepton,p4lepton);
161 trino->init(nudaug,p4nu);
163 CalcAmp(root_part,amp,FormFactors);
165 //Now find the probability at this q2 and cos theta lepton point
166 //and compare to maxfoundprob.
168 //Do a little magic to get the probability!!
169 prob = rho.normalizedProb(amp.getSpinDensity());
174 //probclt contains prob at ctl=-1,0,1.
175 //prob=a+b*ctl+c*ctl^2
178 double b=0.5*(probctl[2]-probctl[0]);
179 double c=0.5*(probctl[2]+probctl[0])-probctl[1];
182 if (probctl[1]>prob) prob=probctl[1];
183 if (probctl[2]>prob) prob=probctl[2];
186 double ctlx=-0.5*b/c;
188 double probtmp=a+b*ctlx+c*ctlx*ctlx;
189 if (probtmp>prob) prob=probtmp;
194 if ( prob > maxfoundprob ) {
199 if ( EvtPDL::getWidth(meson) <= 0.0 ) {
200 //if the particle is narrow dont bother with changing the mass.
205 root_part->deleteTree();