38 {
39
40
41
42
43
44
45
46
47
48
49
50
51 EvtScalarParticle* scalar_part;
52 EvtParticle* root_part;
53
54 scalar_part = new EvtScalarParticle;
55
56
58
59 EvtVector4R p_init;
60
62 scalar_part->
init( parent, p_init );
63 root_part = (EvtParticle*)scalar_part;
64
66
67 EvtParticle *daughter, *lep, *trino;
68
69 EvtAmp amp;
70
71 EvtId listdaug[3];
72 listdaug[0] = meson;
73 listdaug[1] = lepton;
74 listdaug[2] = nudaug;
75
76 amp.
init( parent, 3, listdaug );
77
79 daughter = root_part->
getDaug( 0 );
81 trino = root_part->
getDaug( 2 );
82
83
87
88
89
90 EvtSpinDensity rho;
92
94
95 double m = root_part->
mass();
96
97 EvtVector4R p4meson, p4lepton, p4nu, p4w;
98 double q2max;
99
100 double q2, elepton, plepton;
101 int i, j;
102 double erho, prho, costl;
103
104 double maxfoundprob = 0.0;
105 double prob = -10.0;
106 int massiter;
107
108 for ( massiter = 0; massiter < 3; massiter++ )
109 {
110
115 if ( massiter == 2 )
116 {
119 }
120
121 q2max = ( m -
mass[0] ) * ( m -
mass[0] );
122
123
124
125 for ( i = 0; i < 25; i++ )
126 {
127 q2 = ( ( i + 0.5 ) * q2max ) / 25.0;
128
129 erho = ( m * m +
mass[0] *
mass[0] - q2 ) / ( 2.0 * m );
130
131 prho = sqrt( erho * erho -
mass[0] *
mass[0] );
132
133 p4meson.
set( erho, 0.0, 0.0, -1.0 * prho );
134 p4w.
set( m - erho, 0.0, 0.0, prho );
135
136
137 elepton = ( q2 +
mass[1] *
mass[1] ) / ( 2.0 * sqrt( q2 ) );
138 plepton = sqrt( elepton * elepton -
mass[1] *
mass[1] );
139
140 double probctl[3];
141
142 for ( j = 0; j < 3; j++ )
143 {
144
145 costl = 0.99 * ( j - 1.0 );
146
147
148
149 p4lepton.
set( elepton, 0.0, plepton * sqrt( 1.0 - costl * costl ), plepton * costl );
150 p4nu.
set( plepton, 0.0, -1.0 * plepton * sqrt( 1.0 - costl * costl ),
151 -1.0 * plepton * costl );
152
153 EvtVector4R boost( ( m - erho ), 0.0, 0.0, 1.0 * prho );
154 p4lepton =
boostTo( p4lepton, boost );
156
157
158
159 daughter->
init( meson, p4meson );
160 lep->
init( lepton, p4lepton );
161 trino->
init( nudaug, p4nu );
162
163 CalcAmp( root_part, amp, FormFactors );
164
165
166
167
168
170
171 probctl[j] = prob;
172 }
173
174
175
176
177 double a = probctl[1];
178 double b = 0.5 * ( probctl[2] - probctl[0] );
179 double c = 0.5 * ( probctl[2] + probctl[0] ) - probctl[1];
180
181 prob = probctl[0];
182 if ( probctl[1] > prob ) prob = probctl[1];
183 if ( probctl[2] > prob ) prob = probctl[2];
184
185 if ( fabs( c ) > 1e-20 )
186 {
187 double ctlx = -0.5 * b / c;
188 if ( fabs( ctlx ) < 1.0 )
189 {
190 double probtmp = a + b * ctlx + c * ctlx * ctlx;
191 if ( probtmp > prob ) prob = probtmp;
192 }
193 }
194
195
196
197
198
199
200 if ( prob > maxfoundprob ) { maxfoundprob = prob; }
201 }
203 {
204
205 massiter = 4;
206 }
207 }
209
210 maxfoundprob *= 1.1;
211 return maxfoundprob;
212}
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)