40 {
41
42
43
44
45
46
47
48
49
50
51
52
55
57
58
60
62
64 scalar_part->
init(parent,p_init);
66
68
70
72
74 listdaug[0] = meson;
75 listdaug[1] = lepton;
76 listdaug[2] = nudaug;
77
78 amp.
init(parent,3,listdaug);
79
84
85
89
90
91
92
95
97
98 double m = root_part->
mass();
99
101 double q2max;
102
103 double q2, elepton, plepton;
104 int i,j;
105 double erho,prho,costl;
106
107 double maxfoundprob = 0.0;
108 double prob = -10.0;
109 int massiter;
110
111 for (massiter=0;massiter<3;massiter++){
112
116 if ( massiter==1 ) {
118 }
119 if ( massiter==2 ) {
122 }
123
125
126
127
128 for (i=0;i<25;i++) {
129 q2 = ((i+0.5)*q2max)/25.0;
130
131 erho = ( m*m +
mass[0]*
mass[0] - q2 )/(2.0*m);
132
133 prho = sqrt(erho*erho-
mass[0]*
mass[0]);
134
135 p4meson.
set(erho,0.0,0.0,-1.0*prho);
136 p4w.
set(m-erho,0.0,0.0,prho);
137
138
139 elepton = (q2+
mass[1]*
mass[1])/(2.0*sqrt(q2));
140 plepton = sqrt(elepton*elepton-
mass[1]*
mass[1]);
141
142 double probctl[3];
143
144 for (j=0;j<3;j++) {
145
146 costl = 0.99*(j - 1.0);
147
148
149
150 p4lepton.
set(elepton,0.0,
151 plepton*sqrt(1.0-costl*costl),plepton*costl);
152 p4nu.
set(plepton,0.0,
153 -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
154
156 p4lepton=
boostTo(p4lepton,boost);
158
159
160
161 daughter->
init(meson,p4meson);
162 lep->
init(lepton,p4lepton);
163 trino->
init(nudaug,p4nu);
164
165 CalcAmp(root_part,amp,FormFactors);
166
167
168
169
170
172
173 probctl[j]=prob;
174 }
175
176
177
178
179 double a=probctl[1];
180 double b=0.5*(probctl[2]-probctl[0]);
181 double c=0.5*(probctl[2]+probctl[0])-probctl[1];
182
183 prob=probctl[0];
184 if (probctl[1]>prob) prob=probctl[1];
185 if (probctl[2]>prob) prob=probctl[2];
186
187 if (fabs(c)>1e-20){
188 double ctlx=-0.5*b/c;
189 if (fabs(ctlx)<1.0){
190 double probtmp=a+b*ctlx+c*ctlx*ctlx;
191 if (probtmp>prob) prob=probtmp;
192 }
193
194 }
195
196
197
198
199
200
201 if ( prob > maxfoundprob ) {
202 maxfoundprob = prob;
203 }
204
205 }
207
208 massiter = 4;
209 }
210
211 }
213
214 maxfoundprob *=1.1;
215 return maxfoundprob;
216
217}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
void init(EvtId p, int ndaug, EvtId *daug)
EvtSpinDensity getSpinDensity()
static double getWidth(EvtId i)
static double getMeanMass(EvtId i)
static double getMinMass(EvtId i)
static double getMaxMass(EvtId i)
static double getMass(EvtId i)
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
int getSpinStates() const
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
void init(EvtId part_n, double e, double px, double py, double pz)
virtual void CalcAmp(EvtParticle *parent, EvtAmp &, EvtSemiLeptonicFF *FormFactors)
double NormalizedProb(const EvtSpinDensity &d)
void set(int i, double d)