175 {
176
177 static int iniflag=0;
178
180
182
183
184
185
186
187
188
189
190
192
193
194 int i,more;
196 int ndaugjs;
197 static int kf[20];
198 EvtId evtnumstable[20],evtnumparton[20];
199 int stableindex[20],partonindex[20];
200 int numstable;
201 int numparton;
202 static int km[20];
204
205 static double px[20],py[20],pz[20],e[20];
206
207
208 if (iniflag==0)
dectes_(&iniflag,&idtau,&ndaugjs,kf,km,px,py,pz,e);
209
211
213
214 do{
215
216 iniflag=iniflag+1;
217 dectes_(&iniflag,&idtau,&ndaugjs,kf,km,px,py,pz,e);
218
219 numstable=0;
220 numparton=0;
221
222 for(i=0;i<ndaugjs;i++){
223
224
226 report(
ERROR,
"EvtGen") <<
"Tauola returned particle:"<<kf[i]<<endl;
227 report(
ERROR,
"EvtGen") <<
"This can not be translated to evt number"<<endl;
228 report(
ERROR,
"EvtGen") <<
"and the decay will be rejected!"<<endl;
229 report(
ERROR,
"EvtGen") <<
"The decay was of particle:"<<idtau<<endl;
230
231 }
232
233
234 if (
abs(kf[i])<=6||kf[i]==21){
235 partonindex[numparton]=i;
237 numparton++;
238 }
239 else{
240 stableindex[numstable]=i;
242 numstable++;
243 }
244
245
246
247
248
249
250 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
251
252 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
253
254 }
255
256 p4[i].
set(e[i],px[i],py[i],pz[i]);
257
258
259 }
260
262
263
264 more=(channel!=-1);
265
266
268
269 }
while( more && (
count<10000) );
270
272 report(
INFO,
"EvtGen") <<
"Too many loops in EvtTauola!!!"<<endl;
274 for(i=0;i<numstable;i++){
276 }
277
278 }
279
280
281
282 if (numparton==0){
283
285 int ndaugFound=0;
286 for(i=0;i<numstable;i++){
287 p->
getDaug(i)->
init(evtnumstable[i],p4[stableindex[i]]);
288 ndaugFound++;
289 }
290 if ( ndaugFound == 0 ) {
291 report(
ERROR,
"EvtGen") <<
"Tauola has failed to do a decay ";
293 assert(0);
294 }
295
296 fixPolarizations(p);
297
298 return ;
299
300 }
301 else{
302
303
304
306
307 for(i=0;i<numparton;i++){
308 p4string+=p4[partonindex[i]];
309 }
310
311 int nprimary=1;
312 type[0]=STRNG;
313 for(i=0;i<numstable;i++){
314 if (km[stableindex[i]]==0){
315 type[nprimary++]=evtnumstable[i];
316 }
317 }
318
320
322
324
325 for(i=0;i<numparton;i++){
326 p4partons[i]=p4[partonindex[i]];
327 }
328
330
331
332
333 nprimary=1;
334
335 for(i=0;i<numstable;i++){
336
337 if (km[stableindex[i]]==0){
338 p->
getDaug(nprimary++)->
init(evtnumstable[i],p4[stableindex[i]]);
339 }
340 }
341
342
343 int nsecond=0;
344 for(i=0;i<numstable;i++){
345 if (km[stableindex[i]]!=0){
346 type[nsecond++]=evtnumstable[i];
347 }
348 }
349
350
352
353 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
354 -p4string.get(2),-p4string.get(3));
355
356 nsecond=0;
357 for(i=0;i<numstable;i++){
358 if (km[stableindex[i]]!=0){
359 p4[stableindex[i]]=
boostTo(p4[stableindex[i]],p4stringboost);
363 nsecond++;
364 }
365 }
366
367 if ( nsecond == 0 ) {
368 report(
ERROR,
"EvtGen") <<
"Jetset has failed to do a decay ";
370 assert(0);
371 }
372
373 fixPolarizations(p);
374
375 return ;
376
377 }
378
379}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
DOUBLE_PRECISION count[2]
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)