Line data Source code
1 : //--------------------------------------------------------------------------
2 : //
3 : // Environment:
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.
7 : //
8 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 1998 Caltech, UCSB
10 : //
11 : // Module: EvtSemiLeptonicAmp.cc
12 : //
13 : // Description: Base class for semileptonic decays
14 : //
15 : // Modification history:
16 : //
17 : // DJL April 17,1998 Module created
18 : //
19 : //------------------------------------------------------------------------
20 : //
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"
36 :
37 : double EvtSemiLeptonicAmp::CalcMaxProb( EvtId parent, EvtId meson,
38 : EvtId lepton, EvtId nudaug,
39 : EvtSemiLeptonicFF *FormFactors ) {
40 :
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.
46 :
47 : //Start by declaring a particle at rest.
48 :
49 : //It only makes sense to have a scalar parent. For now.
50 : //This should be generalized later.
51 :
52 : EvtScalarParticle *scalar_part;
53 : EvtParticle *root_part;
54 :
55 0 : scalar_part=new EvtScalarParticle;
56 :
57 : //cludge to avoid generating random numbers!
58 0 : scalar_part->noLifeTime();
59 :
60 0 : EvtVector4R p_init;
61 :
62 0 : p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
63 0 : scalar_part->init(parent,p_init);
64 : root_part=(EvtParticle *)scalar_part;
65 0 : root_part->setDiagonalSpinDensity();
66 :
67 : EvtParticle *daughter, *lep, *trino;
68 :
69 0 : EvtAmp amp;
70 :
71 0 : EvtId listdaug[3];
72 0 : listdaug[0] = meson;
73 0 : listdaug[1] = lepton;
74 0 : listdaug[2] = nudaug;
75 :
76 0 : amp.init(parent,3,listdaug);
77 :
78 0 : root_part->makeDaughters(3,listdaug);
79 0 : daughter=root_part->getDaug(0);
80 0 : lep=root_part->getDaug(1);
81 0 : trino=root_part->getDaug(2);
82 :
83 : //cludge to avoid generating random numbers!
84 0 : daughter->noLifeTime();
85 0 : lep->noLifeTime();
86 0 : trino->noLifeTime();
87 :
88 :
89 : //Initial particle is unpolarized, well it is a scalar so it is
90 : //trivial
91 0 : EvtSpinDensity rho;
92 0 : rho.setDiag(root_part->getSpinStates());
93 :
94 : double mass[3];
95 :
96 0 : double m = root_part->mass();
97 :
98 0 : EvtVector4R p4meson, p4lepton, p4nu, p4w;
99 : double q2max;
100 :
101 : double q2, elepton, plepton;
102 : int i,j;
103 : double erho,prho,costl;
104 :
105 : double maxfoundprob = 0.0;
106 : double prob = -10.0;
107 : int massiter;
108 :
109 0 : for (massiter=0;massiter<3;massiter++){
110 :
111 0 : mass[0] = EvtPDL::getMeanMass(meson);
112 0 : mass[1] = EvtPDL::getMeanMass(lepton);
113 0 : mass[2] = EvtPDL::getMeanMass(nudaug);
114 0 : if ( massiter==1 ) {
115 0 : mass[0] = EvtPDL::getMinMass(meson);
116 0 : }
117 0 : if ( massiter==2 ) {
118 0 : mass[0] = EvtPDL::getMaxMass(meson);
119 0 : if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
120 : }
121 :
122 0 : q2max = (m-mass[0])*(m-mass[0]);
123 :
124 : //loop over q2
125 :
126 0 : for (i=0;i<25;i++) {
127 0 : q2 = ((i+0.5)*q2max)/25.0;
128 :
129 0 : erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
130 :
131 0 : prho = sqrt(erho*erho-mass[0]*mass[0]);
132 :
133 0 : p4meson.set(erho,0.0,0.0,-1.0*prho);
134 0 : p4w.set(m-erho,0.0,0.0,prho);
135 :
136 : //This is in the W rest frame
137 0 : elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
138 0 : plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
139 :
140 0 : double probctl[3];
141 :
142 0 : for (j=0;j<3;j++) {
143 :
144 0 : costl = 0.99*(j - 1.0);
145 :
146 : //These are in the W rest frame. Need to boost out into
147 : //the B frame.
148 0 : p4lepton.set(elepton,0.0,
149 0 : plepton*sqrt(1.0-costl*costl),plepton*costl);
150 0 : p4nu.set(plepton,0.0,
151 0 : -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
152 :
153 0 : EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
154 0 : p4lepton=boostTo(p4lepton,boost);
155 0 : p4nu=boostTo(p4nu,boost);
156 :
157 : //Now initialize the daughters...
158 :
159 0 : daughter->init(meson,p4meson);
160 0 : lep->init(lepton,p4lepton);
161 0 : trino->init(nudaug,p4nu);
162 :
163 0 : CalcAmp(root_part,amp,FormFactors);
164 :
165 : //Now find the probability at this q2 and cos theta lepton point
166 : //and compare to maxfoundprob.
167 :
168 : //Do a little magic to get the probability!!
169 0 : prob = rho.normalizedProb(amp.getSpinDensity());
170 :
171 0 : probctl[j]=prob;
172 0 : }
173 :
174 : //probclt contains prob at ctl=-1,0,1.
175 : //prob=a+b*ctl+c*ctl^2
176 :
177 0 : double a=probctl[1];
178 0 : double b=0.5*(probctl[2]-probctl[0]);
179 0 : double c=0.5*(probctl[2]+probctl[0])-probctl[1];
180 :
181 : prob=probctl[0];
182 0 : if (probctl[1]>prob) prob=probctl[1];
183 0 : if (probctl[2]>prob) prob=probctl[2];
184 :
185 0 : if (fabs(c)>1e-20){
186 0 : double ctlx=-0.5*b/c;
187 0 : if (fabs(ctlx)<1.0){
188 0 : double probtmp=a+b*ctlx+c*ctlx*ctlx;
189 0 : if (probtmp>prob) prob=probtmp;
190 0 : }
191 :
192 0 : }
193 :
194 0 : if ( prob > maxfoundprob ) {
195 : maxfoundprob = prob;
196 0 : }
197 :
198 0 : }
199 0 : if ( EvtPDL::getWidth(meson) <= 0.0 ) {
200 : //if the particle is narrow dont bother with changing the mass.
201 : massiter = 4;
202 0 : }
203 :
204 : }
205 0 : root_part->deleteTree();
206 :
207 0 : maxfoundprob *=1.1;
208 : return maxfoundprob;
209 :
210 0 : }
211 :
|
