142 {
143
144 static int iniflag = 0;
145
147
149
150
151
152
153
154
155
156
157
158 EvtVector4R p4[20];
159
160 int i, more;
162 int ndaugjs;
163 static int kf[20];
164 EvtId evtnumstable[20], evtnumparton[20];
165 int stableindex[20], partonindex[20];
166 int numstable;
167 int numparton;
168 static int km[20];
170
171 static double px[20], py[20], pz[20], e[20];
172
173
174
175 if ( iniflag == 0 )
dectes_( &iniflag, &idtau, &ndaugjs, kf, km, px, py, pz, e );
176
178
180
181 do {
182
183 iniflag = iniflag + 1;
184 dectes_( &iniflag, &idtau, &ndaugjs, kf, km, px, py, pz, e );
185
186 numstable = 0;
187 numparton = 0;
188
189 for ( i = 0; i < ndaugjs; i++ )
190 {
191
192
193
195 {
196 report(
ERROR,
"EvtGen" ) <<
"Tauola returned particle:" << kf[i] << endl;
197 report(
ERROR,
"EvtGen" ) <<
"This can not be translated to evt number" << endl;
198 report(
ERROR,
"EvtGen" ) <<
"and the decay will be rejected!" << endl;
199 report(
ERROR,
"EvtGen" ) <<
"The decay was of particle:" << idtau << endl;
200 }
201
202
203 if (
abs( kf[i] ) <= 6 || kf[i] == 21 )
204 {
205 partonindex[numparton] = i;
207 numparton++;
208 }
209 else
210 {
211 stableindex[numstable] = i;
213 numstable++;
214 }
215
216
217
218
219
220 if ( px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] >= e[i] * e[i] )
221 { e[i] = sqrt( px[i] * px[i] + py[i] * py[i] + pz[i] * pz[i] ) + 0.0000000001; }
222
223 p4[i].
set( e[i], px[i], py[i], pz[i] );
224 }
225
227
228 more = ( channel != -1 );
229
231
232 }
while ( more && (
count < 10000 ) );
233
235 {
236 report(
INFO,
"EvtGen" ) <<
"Too many loops in EvtTauola!!!" << endl;
238 for ( i = 0; i < numstable; i++ )
239 {
240 report(
INFO,
"EvtGen" ) <<
"Daug(" << i <<
")"
242 }
243 }
244
245 if ( numparton == 0 )
246 {
247
249 int ndaugFound = 0;
250 for ( i = 0; i < numstable; i++ )
251 {
252 p->
getDaug( i )->
init( evtnumstable[i], p4[stableindex[i]] );
253 ndaugFound++;
254 }
255 if ( ndaugFound == 0 )
256 {
257 report(
ERROR,
"EvtGen" ) <<
"Tauola has failed to do a decay ";
259 << endl;
260 assert( 0 );
261 }
262
263 fixPolarizations( p );
264
265 return;
266 }
267 else
268 {
269
270
271
272 EvtVector4R p4string( 0.0, 0.0, 0.0, 0.0 );
273
274 for ( i = 0; i < numparton; i++ ) { p4string += p4[partonindex[i]]; }
275
276 int nprimary = 1;
277 type[0] = STRNG;
278 for ( i = 0; i < numstable; i++ )
279 {
280 if ( km[stableindex[i]] == 0 ) { type[nprimary++] = evtnumstable[i]; }
281 }
282
284
286
287 EvtVector4R p4partons[10];
288
289 for ( i = 0; i < numparton; i++ ) { p4partons[i] = p4[partonindex[i]]; }
290
291 ( (EvtStringParticle*)p->
getDaug( 0 ) )->initPartons( numparton, p4partons, evtnumparton );
292
293 nprimary = 1;
294
295 for ( i = 0; i < numstable; i++ )
296 {
297
298 if ( km[stableindex[i]] == 0 )
299 { p->
getDaug( nprimary++ )->
init( evtnumstable[i], p4[stableindex[i]] ); }
300 }
301
302 int nsecond = 0;
303 for ( i = 0; i < numstable; i++ )
304 {
305 if ( km[stableindex[i]] != 0 ) { type[nsecond++] = evtnumstable[i]; }
306 }
307
309
310 EvtVector4R p4stringboost( p4string.get( 0 ), -p4string.get( 1 ), -p4string.get( 2 ),
311 -p4string.get( 3 ) );
312
313 nsecond = 0;
314 for ( i = 0; i < numstable; i++ )
315 {
316 if ( km[stableindex[i]] != 0 )
317 {
318 p4[stableindex[i]] =
boostTo( p4[stableindex[i]], p4stringboost );
322 nsecond++;
323 }
324 }
325
326 if ( nsecond == 0 )
327 {
328 report(
ERROR,
"EvtGen" ) <<
"Jetset has failed to do a decay ";
330 << endl;
331 assert( 0 );
332 }
333
334 fixPolarizations( p );
335
336 return;
337 }
338}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
DOUBLE_PRECISION count[3]
void dectes_(int *, int *, int *, int *, int *, double *, double *, double *, double *)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
static int getStdHep(EvtId id)
static EvtId evtIdFromStdHep(int stdhep)
static std::string name(EvtId i)
static EvtId getId(const std::string &name)
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void setDiagonalSpinDensity()
EvtParticle * getDaug(int i)
void deleteDaughters(bool keepChannel=false)
void set(int i, double d)