Bug Summary

File:programs/Simulation/genr8/genr8.c
Location:line 768, column 12
Description:The right operand of '<' is a garbage value

Annotated Source Code

1 /********************************************************
2 *
3 * Usage: genr8 <options> < input.gen
4 * Use "genr8 -h" for help with options.
5 ********************************************************
6 * * Generate t-channel
7 * genr8.c * monte carlo events.
8 * *
9 ********************************************************
10 *
11 * created by: Paul M Eugenio
12 * Carnegie Mellon University
13 * 25-Mar-98
14 *
15 * minor modifications to avoid infinite loop in n_omega_pi0_pi+ generator
16 * garth huber, 04.04.21
17 ******************************************************** */
18
19#include <stdio.h>
20#include <unistd.h>
21#include <math.h>
22#include <string.h>
23#include <stdlib.h>
24#include <sys/types.h>
25#include <sys/signal.h>
26#include <sys/stat.h>
27#include <time.h>
28
29#include <genkin.h>
30#include <particleType.h>
31
32#define TRUE1 1
33#define FALSE0 0
34#define CONV(180.0/3.14159265358979323846) (180.0/M_PI3.14159265358979323846)
35#define BUFSIZE100000 100000
36#define T_CUT10.0 10.0
37#define RECOIL0.938 0.938
38
39#define RESTFRAME-1 -1
40#define PARENTFRAME+1 +1
41
42#define PRODUCTION_PARTICLE1 1
43
44
45
46
47/***********************/
48/* STRUCTURES */
49/***********************/
50
51
52
53struct particleMC_t{
54 int flag;
55 int nchildren;
56 int charge;
57 double mass;
58 double bookmass;
59 double width;
60 Particle_t particleID;
61 vector4_t p;
62 struct particleMC_t *parent, *child[2];
63} ;
64
65
66/******************************************************************
67 * GLOBAL VARIABLES
68 * NOTE: Please start the global name
69 * with one capital letter!!!!!!
70 * Use all lower case for local
71 * variable names. Thank You 8^)
72 **********************************************************************/
73
74int Debug = 0;
75int Nprinted =0;
76int PrintProduction=0;
77int PrintRecoil=0;
78double MassHighBW;
79int UseName=0;
80int FIRST_EVENT=1;
81int PrintFlag=10;
82int WriteAscii=0;
83int runNo=9000;
84int NFinalParts=0;
85unsigned int RandomSeed=0;
86int UseCurrentTimeForRandomSeed = TRUE1;
87/***********************/
88/* Declarations */
89/***********************/
90/*
91 * These functions are coded after main().
92 */
93
94
95double rawthresh(struct particleMC_t *Isobar);
96void decay(struct particleMC_t *Isobar);
97void boost2lab(struct particleMC_t *Isobar);
98void boostFamily(vector4_t *beta,struct particleMC_t *Isobar);
99void boost(vector4_t *beta,vector4_t *vec);
100void printParticle(struct particleMC_t *Isobar);
101vector4_t polarMake4v(double p, double theta, double phi, double mass);
102double randm(double low, double high);
103void printProduction(FILE *fp,struct particleMC_t *Isobar);
104void printFinal(FILE *fp,struct particleMC_t *Isobar);
105void printp2ascii(FILE *fp,struct particleMC_t *Isobar);
106void setMass(struct particleMC_t *Isobar);
107void initMass(struct particleMC_t *Isobar);
108char *ParticleType(Particle_t p);
109void checkFamily(struct particleMC_t *Isobar);
110int setChildrenMass(struct particleMC_t *Isobar);
111void printFamily(struct particleMC_t *Isobar);
112void lorentzFactor(double *lf,struct particleMC_t *Isobar);
113
114/*
115 ***********************
116 * *
117 * PrintUsage() *
118 * *
119 ***********************
120 */
121
122void PrintUsage(char *processName)
123{
124
125 fprintf(stderrstderr,"%s usage: [-A<name>] < infile \n",processName);
126 fprintf(stderrstderr,"\t-d debug flag\n");
127 fprintf(stderrstderr,"\t-n Use a particle name and not its ID number (ascii only) \n");
128 fprintf(stderrstderr,"\t-M<max> Process first max events\n");
129 fprintf(stderrstderr,"\t-l<lfevents> Determine the lorentz factor with this many number of events (default is 10000)\n");
130 fprintf(stderrstderr,"\t-r<runNo> default runNo is 9000. \n");
131 /*fprintf(stderr,"\t-o<name> The output file \n");*/
132 fprintf(stderrstderr,"\t-P save flag= 11 & 01 events(default saves 11 & 10 events) \n");
133 /* fprintf(stderr,"\t-R Save recoiling baryon information. \n"); */
134
135 fprintf(stderrstderr,"\t-A<filename> Save in ascii format. \n");
136 fprintf(stderrstderr,"\t-s<seed> Set random number seed to <seed>. \n");
137 fprintf(stderrstderr,"\t (default is to set using current time + pid) \n");
138 fprintf(stderrstderr,"\t-h Print this help message\n\n");
139
140}
141
142
143/*
144 ***********************
145 * *
146 * Main() *
147 * *
148 ***********************
149 */
150
151int main(int argc,char **argv)
152{
153 char *argptr,*token,line[2056];
154 int i,npart=0,ngenerated=0,naccepted=0, imassc, imassc2;
155 int nv4,max=10,part=0,chld1=-1,chld2=-1,prnt=-1,lfevents=10000;
156 FILE *fout=stdoutstdout;
157 struct particleMC_t particle[20],beam,target,CM;
158 //struct particleMC_t recoil;
159 struct particleMC_t *X,*Y;
160 vector4_t beta,v4[2];
161 // vector4_t initBeam4;
162 double t,expt_max,expt,expt_min,sqrt_s,t_min=0;
163 double CMenergy, t_max,slope=5.0;
164 double X_momentum, X_energy,xmass,ymass;
165 // double X_threshold ;
166 double costheta,theta,phi,lf,lfmax=0;
167 int isacomment=TRUE1,haveChildren=TRUE1;
168
169 Y= &(particle[0]);
170 Y->parent = &CM;
171 X= &(particle[1]);
172 X->parent = &CM;
173 // recoil.parent = &CM;
174 CM.child[0]= X;
175 CM.child[1]= Y;
176 /* CM.child[1]= &recoil; */
177
178 if (argc == 1){
1
Assuming 'argc' is not equal to 1
→
2
←
Taking false branch
→
179 PrintUsage(argv[0]);
180 exit (0);
181 }
182 else {
183
184 /*
185 * Read command line options.
186 */
187
188 for (i=1; i<argc; i++) {
3
←
Loop condition is true. Entering loop body
→
4
←
Assuming 'i' is >= 'argc'
→
5
←
Loop condition is false. Execution continues on line 246
→
189 argptr = argv[i];
190 if ((*argptr == '-') && (strlen(argptr) > 1)) {
191 argptr++;
192 switch (*argptr) {
193 case 'd':
194 Debug =1;
195 break;
196 case 'h':
197 PrintUsage(argv[0]);
198 exit(0);
199 break;
200 case 'n':
201 UseName =1;
202 break;
203 case 'A':
204 WriteAscii=1;
205 fout = fopen(++argptr,"w");
206 fprintf(stderrstderr,"Opening file %s for output. \n",argptr);
207 break;
208 case 'R':
209 fprintf(stderrstderr,"Printing recoil information.\n");
210 PrintRecoil=1;
211 break;
212 case 'P':
213 fprintf(stderrstderr,"Printing eta and pizeros and not gammas.\n");
214 PrintProduction=1;
215 break;
216 case 'l':
217 lfevents = atoi(++argptr);
218 fprintf(stderrstderr,"Using %d events to determine the lorentz factor\n",lfevents);
219 break;
220 case 'M':
221 max = atoi(++argptr);
222 fprintf(stderrstderr,"Maximum number of events: %d\n",max);
223 break;
224 case 'r':
225 runNo = atoi(++argptr);
226 fprintf(stderrstderr,"Using runNo: %d\n",runNo);
227 break;
228 case 's':
229 RandomSeed = atoi(++argptr);
230 UseCurrentTimeForRandomSeed = FALSE0;
231 break;
232 default:
233 fprintf(stderrstderr,"Unrecognized argument -%s\n\n",argptr);
234 PrintUsage(argv[0]);
235 exit(-1);
236 break;
237
238 }
239 }
240 }
241 }
242
243 /*
244 * Seed the random number generator.
245 */
246 if(UseCurrentTimeForRandomSeed){
6
←
Assuming 'UseCurrentTimeForRandomSeed' is 0
→
7
←
Taking false branch
→
247 RandomSeed=time(NULL((void*)0));
248 RandomSeed += getpid();
249 }
250 printf("Setting random number seed to: %d\n",RandomSeed);
251 srand48(RandomSeed);
252
253 /*
254 * Now read the input.gen file
255 * from the stdin.
256 *
257 *
258 * Any line starting with a "%"
259 * is a comment line and is ignored.
260 *
261 *
262 */
263
264
265/* Fill particle information */
266
267 isacomment=TRUE1;
268 while(isacomment==TRUE1){
8
←
Loop condition is true. Entering loop body
→
11
←
Loop condition is true. Entering loop body
→
14
←
Loop condition is true. Entering loop body
→
18
←
Loop condition is false. Execution continues on line 277
→
269 char *pline;
270 pline = fgets(line,sizeof(line),stdinstdin);
271 if (pline!=NULL((void*)0)){
9
←
Assuming 'pline' is equal to null
→
10
←
Taking false branch
→
12
←
Assuming 'pline' is equal to null
→
13
←
Taking false branch
→
15
←
Assuming 'pline' is not equal to null
→
16
←
Taking true branch
→
272 token=strtok(line," ");
273 if(!(*token == '%'))
17
←
Taking true branch
→
274 isacomment=FALSE0;
275 }
276 } /* get beam information */
277 beam.p.space.x = atof(token);
278 token=strtok(NULL((void*)0)," ");
279 beam.p.space.y = atof(token);
280 token=strtok(NULL((void*)0)," ");
281 beam.p.space.z = atof(token);
282 token=strtok(NULL((void*)0)," ");
283 beam.mass = atof(token);
284 fprintf(stderrstderr,"Reading: \tbeamp.x \tbeamp.y \tbeamp.z \tbeamMass\n");
285 fprintf(stderrstderr,"Found: \t\t%lf \t%lf \t%lf \t%lf \n",
286 beam.p.space.x, beam.p.space.y, beam.p.space.z, beam.mass);
287
288
289 isacomment=TRUE1;
290 while(isacomment==TRUE1){
19
←
Loop condition is true. Entering loop body
→
22
←
Loop condition is true. Entering loop body
→
25
←
Loop condition is true. Entering loop body
→
29
←
Loop condition is false. Execution continues on line 299
→
291 char *pline;
292 pline = fgets(line,sizeof(line),stdinstdin);
293 if (pline!=NULL((void*)0)){
20
←
Assuming 'pline' is equal to null
→
21
←
Taking false branch
→
23
←
Assuming 'pline' is equal to null
→
24
←
Taking false branch
→
26
←
Assuming 'pline' is not equal to null
→
27
←
Taking true branch
→
294 token=strtok(line," ");
295 if(!(*token == '%'))
28
←
Taking true branch
→
296 isacomment=FALSE0;
297 }
298 } /* get target information */
299 target.p.space.x = atof(token);
300 token=strtok(NULL((void*)0)," ");
301 target.p.space.y = atof(token);
302 token=strtok(NULL((void*)0)," ");
303 target.p.space.z = atof(token);
304 token=strtok(NULL((void*)0)," ");
305 target.mass = atof(token);
306 fprintf(stderrstderr,"Reading: \ttargetp.x \ttargetp.y \ttargetp.z \ttargetMass\n");
307 fprintf(stderrstderr,"Found: \t\t%lf \t%lf \t%lf \t%lf \n",
308 target.p.space.x, target.p.space.y, target.p.space.z, target.mass);
309
310
311 isacomment=TRUE1;
312 while(isacomment==TRUE1){
30
←
Loop condition is true. Entering loop body
→
33
←
Loop condition is true. Entering loop body
→
36
←
Loop condition is true. Entering loop body
→
40
←
Loop condition is false. Execution continues on line 322
→
313 char *pline;
314 pline = fgets(line,sizeof(line),stdinstdin);
315 if (pline!=NULL((void*)0)){
31
←
Assuming 'pline' is equal to null
→
32
←
Taking false branch
→
34
←
Assuming 'pline' is equal to null
→
35
←
Taking false branch
→
37
←
Assuming 'pline' is not equal to null
→
38
←
Taking true branch
→
316 token=strtok(line," ");
317 if(!(*token == '%'))
39
←
Taking true branch
→
318 isacomment=FALSE0;
319 }
320 }
321 /* get the t-channel slope */
322 slope=atof(token);
323 fprintf(stderrstderr,"Reading: t-channelSlope\n");
324 fprintf(stderrstderr,"Found: \t%lf \n",slope);
325
326 isacomment=TRUE1;
327 while(isacomment==TRUE1){
41
←
Loop condition is true. Entering loop body
→
44
←
Loop condition is true. Entering loop body
→
47
←
Loop condition is true. Entering loop body
→
51
←
Loop condition is false. Execution continues on line 336
→
328 char *pline;
329 pline = fgets(line,sizeof(line),stdinstdin);
330 if (pline!=NULL((void*)0)){
42
←
Assuming 'pline' is equal to null
→
43
←
Taking false branch
→
45
←
Assuming 'pline' is equal to null
→
46
←
Taking false branch
→
48
←
Assuming 'pline' is not equal to null
→
49
←
Taking true branch
→
331 token=strtok(line," ");
332 if(!(*token == '%'))
50
←
Taking true branch
→
333 isacomment=FALSE0;
334 }
335 } /* get the number of particles to read in below */
336 npart = atoi(token);
337 fprintf(stderrstderr,"Reading: number of particles need to describe the decay\n");
338 fprintf(stderrstderr,"Found: \t%d \n",npart);
339
340 /*
341 * read all particles needed
342 * to decsribe an isobar decay
343 * of the resonance (X)
344 */
345
346 fprintf(stderrstderr,"Reading: \tpart# \tchld1# \tchld2# \tprnt# \tId \tnchld \tmass \t\twidth \t\tchrg \tflag \n");
347
348 for(i=0;i<npart;i++){ /* read particle information*/
52
←
Assuming 'i' is >= 'npart'
→
53
←
Loop condition is false. Execution continues on line 433
→
349 haveChildren=TRUE1;
350 isacomment=TRUE1;
351 while(isacomment==TRUE1){
352 char *pline;
353 pline = fgets(line,sizeof(line),stdinstdin);
354 if (pline!=NULL((void*)0)){
355 token=strtok(line," ");
356 if(!(*token == '%'))
357 isacomment=FALSE0;
358 }
359 }
360 if(!(*token == '*')) /* "*" means it's unknown at this time */
361 part = atoi(token);
362 token=strtok(NULL((void*)0)," ");
363 if(!(*token == '*')){ /* set pointer to child */
364 chld1 = atoi(token);
365 particle[part].child[0] = &(particle[chld1]);
366 }
367 else
368 haveChildren = FALSE0;
369 token=strtok(NULL((void*)0)," ");
370 if(!(*token == '*')) {/* set pointer to child */
371 chld2 = atoi(token);
372 particle[part].child[1] = &(particle[chld2]);
373 }
374 else
375 haveChildren = FALSE0;
376 token=strtok(NULL((void*)0)," ");
377 if(!(*token == '*')) {/* set pointer to parent */
378 prnt = atoi(token);
379 particle[part].parent = &(particle[prnt]);
380 }
381 token=strtok(NULL((void*)0)," ");
382 if(!(*token == '*'))
383 /* particle[part].particleID = part; */
384 particle[part].particleID = atoi(token);/* see particleType.h */
385 token=strtok(NULL((void*)0)," ");
386 if(!(*token == '*')){
387 particle[part].nchildren = atoi(token);
388 if(haveChildren==FALSE0 && particle[part].nchildren>0 ){
389 fprintf(stderrstderr,
390 "If a particle has children then it must point to them!\n");
391 exit(-1);
392 }
393 }
394 token=strtok(NULL((void*)0)," ");
395 if(!(*token == '*'))
396 particle[part].bookmass = atof(token);
397 else{ /* get a list of the particle that need a mass generated */
398 fprintf(stderrstderr,"Every Particle needs a mass\n");
399 exit(-1);
400 }
401 token=strtok(NULL((void*)0)," ");
402 if(!(*token == '*')){
403 particle[part].width = atof(token);
404
405 }
406 else {/* for a fixed mass use a zero width */
407 fprintf(stderrstderr,"Every Particle needs a width\n");
408 exit(-1);
409 }
410 token=strtok(NULL((void*)0)," ");
411 if(!(*token == '*'))
412 particle[part].charge = atoi(token);
413 token=strtok(NULL((void*)0)," ");
414 if(!(*token == '*')){
415 particle[part].flag = atoi(token);
416 if(PrintProduction==1){
417 if(particle[part].flag==11 || particle[part].flag==01)
418 NFinalParts++;
419 }else{
420 if(particle[part].flag==11 || particle[part].flag==10)
421 NFinalParts++;
422 }
423 }
424 /* flag 00 = isobar or resonace
425 * flag 01 = production particle that decays i.e. eta, pizero ..
426 * flag 11 = production particle that does not decay i.e. piplus,...
427 * flag 10 = final state particle not in production i.e. gamma
428 */
429 fprintf(stderrstderr,
430 "Found: \t\t%d \t%d \t%d \t%d \t%d \t%d \t%lf \t%lf \t%d \t%d\n",
431 part,chld1,chld2,prnt,particle[part].particleID,particle[part].nchildren,particle[part].bookmass,particle[part].width,particle[part].charge,particle[part].flag);
432 }
433 isacomment=TRUE1;
434 while(isacomment==TRUE1){
54
←
Loop condition is true. Entering loop body
→
58
←
Loop condition is false. Execution continues on line 443
→
435 char *pline;
436 pline = fgets(line,sizeof(line),stdinstdin);
437 if (pline!=NULL((void*)0)){
55
←
Assuming 'pline' is not equal to null
→
56
←
Taking true branch
→
438 token=strtok(line," ");
439 if(!(*token == '%'))
57
←
Taking true branch
→
440 isacomment=FALSE0;
441 }
442 }
443 if(!(*token == '!')){
59
←
Taking false branch
→
444 fprintf(stderrstderr,"Failed to find EOI---- Check Input File\n");
445 exit(-1);
446 }
447
448 checkFamily(X);
60
←
Calling 'checkFamily'
→
449 fprintf(stderrstderr,"Found EOI---- Input File appears Fine.\n");
450
451 /* We are now done reading the input information */
452
453 /*
454 * The beam and target are in the lab frame.
455 * Put them in the overall center of momentum (CM) frame
456 * and calculate |t| & recoil angles.
457 */
458
459 if(X->nchildren == 0) X->mass = X->bookmass;
460 if(Y->nchildren == 0) Y->mass = Y->bookmass;
461
462 target.p.t = energy(target.mass,&(target.p.space));
463 beam.p.t = energy(beam.mass,&(beam.p.space));
464 //initBeam4.t= beam.p.t; initBeam4.space.x= beam.p.space.x;
465 //initBeam4.space.y= beam.p.space.y; initBeam4.space.z= beam.p.space.z;
466 sqrt_s = sqrt( SQ(beam.mass) +SQ(target.mass) + 2.0*beam.p.t * target.p.t);
467 /*MassHighBW = sqrt_s - recoil.mass; */
468 MassHighBW = sqrt_s; /* see do loop below */
469
470 v4[0]= beam.p;
471 v4[1]= target.p;
472 nv4=2;
473
474 CM.mass = sqrt_s;
475 CM.p = Sum4vec(v4,nv4);
476 beta = get_beta(&(CM.p),RESTFRAME-1);
477 boost(&beta,&(beam.p));
478 boost(&beta,&(target.p));
479
480 initMass(X);
481 initMass(Y);
482 CMenergy = beam.p.t + target.p.t;
483
484 while(naccepted <max){
485
486
487 if (Debug) fprintf(stderrstderr,"In main do loop ... %d \n",naccepted);
488
489 /*
490 * Generate the resonance
491 * in the CM frame, and
492 * fill the four vectors
493 * for both X and the recoil.
494 */
495
496 do{
497l0: imassc2=0;
498 if(!(X->width<0))
499 do{/*use BreitWigner--phasespace distribution */
500 initMass(X);
501 setMass(X);
502
503 /*
504 * set the children mass to the book mass or
505 * distribute the by a Breit-Wigner. If the isobar
506 * mass is unknown it's mass remains unknow at this time.
507 */
508
509 if(Debug) fprintf(stderrstderr,"looping over nchildren = %d \n",X->nchildren);
510 for(i=0;i<X->nchildren;i++)
511 {
512 if(Debug) fprintf(stderrstderr,"calling setChildrenMass X... %d \n",i);
513l1: imassc=setChildrenMass(X->child[i]);
514 if (Debug) fprintf(stderrstderr,"Return from setChildrenMass X... %d %d \n",i,imassc);
515
516 /* if the daughters of child[i] are more massive than child[i], generate masses */
517 if (imassc!=0) {
518 if (Debug) fprintf(stderrstderr,"Need new masses for X %d %d %f \n",i,imassc,(X->child[i])->mass);
519 imassc2=imassc2+1;
520 if (imassc2<1000)
521 goto l1;
522 else
523 goto l0;
524 }
525 }
526 }while((X->mass > MassHighBW) || ( X->nchildren==0 ? FALSE0 :
527 (X->mass < ( (X->child[0])->mass + (X->child[1])->mass))) );
528
529 else{/* there's an error.. */
530 fprintf(stderrstderr,"Cannot use a negative width!\n");
531 exit(-1);
532 }
533 /*
534 * Now do loop it for Y
535 */
536 if(!(Y->width<0))
537 do{/*use BreitWigner--phasespace distribution */
538 initMass(Y);
539 setMass(Y);
540
541 /*
542 * set the children mass to the book mass or
543 * distribute the by a Breit-Wigner. If the isobar
544 * mass is unknown it's mass remains unknown at this time.
545 */
546
547 for(i=0;i<Y->nchildren;i++)
548 {
549 if (Debug) fprintf(stderrstderr,"calling setChildrenMass Y... %d \n",i);
550l2: imassc=setChildrenMass(Y->child[i]);
551 if (Debug) fprintf(stderrstderr,"Return from setChildrenMass Y... %d %d \n",i,imassc);
552
553 /* if the daughters of child[i] are more massive than child[i], generate masses */
554 if (imassc!=0) {
555 if (Debug) fprintf(stderrstderr,"Need new masses for Y %d %d %f \n",i,imassc,(Y->child[i])->mass);
556 goto l2; }
557 }
558 }while((Y->mass > MassHighBW) || ( Y->nchildren==0 ? FALSE0 :
559 (Y->mass < ( (Y->child[0])->mass + (Y->child[1])->mass)) ) );
560 else{/* there's an error.. */
561 fprintf(stderrstderr,"Cannot use a negative width!\n");
562 exit(-1);
563 }
564 }while(sqrt_s < X->mass + Y->mass);
565 /*
566 xmass=rawthresh(X);
567 ymass=rawthresh(Y);
568 *
569 * fprintf(stderr," xmass= %lf ymass= %lf X->mass= %lf Y->mass= %lf\n",
570 * xmass,ymass, X->mass , Y->mass);
571 */
572 xmass = X->mass;
573 ymass = Y->mass;
574
575 //X_threshold = 0;
576 X_momentum = CMmomentum( CMenergy, X->mass, Y->mass);
577 X_energy = sqrt( (X->mass)*(X->mass) + X_momentum*X_momentum);
578
579 if(Y->nchildren ==0){
580
581 t_min = -( SQ( (SQ(beam.mass) -SQ(xmass) -SQ(target.mass) +SQ(ymass))/(2.0*sqrt_s))
582 -SQ(v3mag(&(beam.p.space)) - X_momentum ));
583 t_max = -( SQ( (SQ(beam.mass) -SQ(xmass) -SQ(target.mass) +SQ(ymass))/(2.0*sqrt_s))
584 -SQ(v3mag(&(beam.p.space)) + X_momentum ));
585 /*
586 *fprintf(stderr,
587 "beam.mass= %lf xmass= %lf target.mass=%lf ymass= %lf sqrt_s= %lf beam.p= %lf X->p= %lf X_momentum= %lf\n",
588 beam.mass,xmass,target.mass,ymass,sqrt_s,
589 v3mag(&(beam.p.space)),v3mag(&(X->p.space)), X_momentum);
590 */
591
592 /*fprintf(stderr,"t_min: %lf t_max: %lf\n", t_min,t_max);
593 */
594 } else{ /* it's some baryon pseudo t process */
595
596 t_min=0.4;
597 t_max=10.0;
598 t_min = -( SQ( (SQ(beam.mass) -SQ(xmass) -SQ(target.mass) +SQ(ymass))/(2.0*sqrt_s))
599 -SQ(v3mag(&(beam.p.space)) - X_momentum ));
600 t_max = -( SQ( (SQ(beam.mass) -SQ(xmass) -SQ(target.mass) +SQ(ymass))/(2.0*sqrt_s))
601 -SQ(v3mag(&(beam.p.space)) + X_momentum ));
602
603 }
604 expt_max = exp(-slope * t_max);
605 expt_min = exp(-slope * t_min);
606
607 do{
608
609 expt = randm(expt_max,expt_min);
610
611 t= -log(expt)/slope;
612 costheta = ( beam.p.t * X_energy -
613 0.5*(t + (beam.mass)*(beam.mass) + (X->mass)*(X->mass))
614 )/( v3mag(&(beam.p.space))*X_momentum ) ;
615
616 }while(fabs(costheta)>1.0 );
617
618 theta = acos(costheta);
619 phi = randm(-1*M_PI3.14159265358979323846,M_PI3.14159265358979323846);
620
621 X->p = polarMake4v(X_momentum,theta,phi,X->mass);
622 Y->p=polarMake4v(X_momentum,(M_PI3.14159265358979323846-theta),(M_PI3.14159265358979323846+phi),Y->mass);
623
624 /*
625 * Now decay X -> children -> grandchildren -> and so forth
626 *
627 * Note: all particles are generated in their parent's rest frame.
628 */
629
630 /*
631 if(Debug)
632 fprintf(stderr,"before decay\n");
633 if(Debug)
634 printFamily(X);
635 */
636
637 decay(X);
638 decay(Y);
639 if(Debug) {
640 fprintf(stderrstderr,"X after decay\n");
641 printFamily(X);
642 fprintf(stderrstderr,"Y after decay\n");
643 printFamily(Y);
644 }
645
646
647 /*
648 * Compute Lorentz Factor (used for phasespace weighting)
649 */
650 lf=v3mag(&(X->p.space));
651 lorentzFactor(&lf,X);
652 lorentzFactor(&lf,Y);
653 if (Debug) fprintf(stderrstderr,"lorentz factor information: %f ... %f ... \n",lf,lfmax);
654 if(lfevents-->0){
655 lfmax = lf >lfmax ? lf : lfmax; /* find the largest value */
656 if( (lfevents % 10) == 0 ) {
657 if ( lfevents <= 100 || (lfevents % 1000) == 0 )
658 fprintf(stderrstderr,"Calculating Lorentz Factor: %d \r",lfevents);
659 }
660 }
661 else{
662
663 if (Debug) fprintf(stderrstderr," inside loop: lorentz factor information: %f ... %f ... \n",lf,lfmax);
664 /*
665 * Now generate the events weighted by phasespace
666 * (the maximum Lorentz factor).
667 *
668 * Since each particle is in its parent's rest frame,
669 * it must be boosted through each parent's -> parent's-> ...
670 * rest frame to the lab frame.
671 */
672 /*
673 fprintf(stderr,"expt_min: %lf \t expt_max: %lf\n",expt_min,expt_max);
674 fprintf(stderr,"t_min: %lf \t t: %lf \t t_max: %lf\n",t_min,t,t_max);
675 */
676 ngenerated++;
677 /* fprintf(stderr,"ngen: %d \tlf: %lf \tlfmax: %lf\n",ngenerated,lf,lfmax);*/
678 if(lf > randm(0.0,lfmax) ){ /* phasespace distribution */
679
680 naccepted++;
681 boost2lab(X);
682 boost2lab(Y);
683 /*initBeam4; use lab frame beam */
684 /*
685 * We have a complete event. Now save it!
686 */
687
688 if(Debug) {
689 fprintf(stderrstderr,"X after boost2lab\n");
690 printFamily(X);
691 fprintf(stderrstderr,"Y after boost2lab\n");
692 printFamily(Y);
693 }
694
695 Nprinted =0;
696 /* event header information
697 fprintf(fout,"RunNo %d EventNo %d\n",runNo,naccepted);*/
698 fprintf(fout,"%d %d %d\n",runNo,naccepted, NFinalParts);
699
700 /*
701 * Print out the production
702 * or the final state particles.
703 *
704 if(WriteEsr)
705 WriteItape(&CM,&initBeam4);
706 * Remove old BNL-E852 dependence
707 */
708
709 if(PrintProduction){
710 if(X->nchildren==0)
711 printp2ascii(fout,X);
712 else
713 printProduction(fout,X);
714 if(Y->nchildren==0)
715 printp2ascii(fout,Y);
716 else
717 printProduction(fout,Y);
718 }
719 else{
720 if(X->nchildren==0)
721 printp2ascii(fout,X);
722 else
723 printFinal(fout,X);
724 if(Y->nchildren==0 && Y->flag/10 == 1 )
725 printp2ascii(fout,Y);
726 else
727 printFinal(fout,Y);
728 }
729
730
731 }
732 if(!(ngenerated % 100)) {
733 if( ngenerated <= 1000 || !(ngenerated % 10000))
734 fprintf(stderrstderr,"Events generated: %d Events accepted: %d \r",
735 ngenerated,naccepted);
736 }
737 if(Debug) fprintf(stderrstderr,"End of event\n");
738 } /* end of else{ */
739 }/* end of while */
740
741 fprintf(stderrstderr,
742 "Max Lorentz Factor:%lf Events generated:%d Events accepted:%d\n\n",
743 lfmax,ngenerated,naccepted);
744 /*
745 * Close the output file.
746 */
747 fflush(fout);
748 fclose(fout);
749
750 return 0;
751}/* end of main */
752
753
754/********************
755 *
756 * checkFamily()
757 * Testing the input file.
758 *
759 * If I have children
760 * then I should be
761 * my children's
762 * parent.
763 *******************/
764void checkFamily(struct particleMC_t *Isobar)
765{
766 int i;
767
768 for(i=0;i<Isobar->nchildren; i++){
61
←
The right operand of '<' is a garbage value
769 if(Isobar->nchildren != 2){
770 fprintf(stderrstderr,"Error in input file: Sorry, only 0 or 2 children are allowed.\n");
771 exit(-1);
772 }
773 if(Isobar != Isobar->child[i]->parent){
774 fprintf(stderrstderr,"Error in input file: Parent to children mismatch\n");
775 exit(-1);
776 }
777 checkFamily(Isobar->child[i]);
778 }
779
780}
781
782
783
784/********************
785 *
786 * setChildrenMass()
787 *
788 * Sets all children masses
789 * to bookmass or to a
790 * Breit-Wigner mass
791 *********************/
792
793int setChildrenMass(struct particleMC_t *Isobar)
794{
795 int i, imassc=0;
796
797/* fprintf(stderr,"In setChildrenMass ... %f \n",Isobar->mass);
798*/
799
800 initMass(Isobar);
801 setMass(Isobar);
802 for(i=0;i < Isobar->nchildren;i++){
803 if (Debug) fprintf(stderrstderr,"In loop ... %d %d %f \n",i,Isobar->nchildren,(Isobar->child[i])->mass);
804
805 /* Generate masses of all of the daughters */
806l3: imassc=setChildrenMass(Isobar->child[i]);
807
808 if (imassc!=0) {
809 if (Debug) fprintf(stderrstderr,"Need new masses for Isobar %d %d %f \n",i,imassc,(Isobar->child[i])->mass);
810 goto l3; }
811 }
812
813 if(Isobar->nchildren !=0)
814 {
815 if((Isobar->mass) < ((Isobar->child[0])->mass)+((Isobar->child[1])->mass))
816 {
817
818 /* If the daughters are more massive than the parent, set the return code and exit */
819 if (Debug) fprintf(stderrstderr,"final call ... \n");
820
821 /* setChildrenMass(Isobar); */
822 imassc=imassc+1;
823 if (Debug) fprintf(stderrstderr," %f %f %f \n",Isobar->mass,(Isobar->child[0])->mass,(Isobar->child[1])->mass);
824 } else {imassc=0;}
825 }
826
827 if (Debug) fprintf(stderrstderr,"Leaving setChildrenMass ... %f %d \n",Isobar->mass,imassc);
828
829 return imassc;
830}
831
832
833
834
835/********************
836 *
837 * setMass()
838 *
839 * Sets the particle mass
840 * using a
841 * Breit-Wigner distribution.
842 *********************/
843
844void setMass(struct particleMC_t *Isobar)
845{
846 double n,height,thresH2,lowtail,hightail,hcut,lcut;
847
848 if(Isobar->width > 0){
849
850 lowtail = rawthresh(Isobar);
851
852
853
854 thresH2 =
855 Isobar == (Isobar->parent)->child[0] ? /* is the 1st child me? */
856 rawthresh((Isobar->parent)->child[1]) : /* if true */
857 rawthresh((Isobar->parent)->child[0]) ;/* if false */
858 /*
859 hightail = ((Isobar->parent)->mass);
860 */
861 hightail = (Isobar->parent->mass) - thresH2 ;
862
863
864 /* cut off the tails */
865 hcut= Isobar->bookmass + 4.0*Isobar->width ;
866 lcut= Isobar->bookmass - 4.0*Isobar->width ;
867
868 if(hightail> hcut)
869 hightail=hcut;
870 if(lowtail < lcut)
871 lowtail= lcut;
872
873
874 do{
875 n=randm(0.0,0.9999);
876 Isobar->mass = randm( lowtail , hightail);
877
878 height= SQ((Isobar->bookmass)*(Isobar->width))/
879 ( SQ(SQ(Isobar->bookmass) - SQ(Isobar->mass)) +
880 SQ((Isobar->bookmass) * (Isobar->width) ));
881 }while(n > height );
882 /*
883 fprintf(stderr,"bookmass is %lf: low: %lf high: %lf bwmass: %lf\n",
884 Isobar->bookmass, lowtail,hightail,Isobar->mass);
885 */
886 }
887}
888/********************
889 *
890 * initMass()
891 *
892 * Sets the particle mass
893 * to bookmass or UNKNOWN
894 *
895 *********************/
896
897void initMass(struct particleMC_t *Isobar)
898{
899 int i;
900
901 for(i=0;i<Isobar->nchildren;i++){
902 if(Isobar->child[i]->width == 0.0)
903 Isobar->child[i]->mass = Isobar->child[i]->bookmass;
904 else
905 Isobar->child[i]->mass = -1.0; /* UNKNOWN */
906 initMass(Isobar->child[i]);
907 }
908}
909/********************
910 *
911 * rawthresh()
912 *
913 * Calculates the mass
914 * threshold for the
915 * isobar.
916 *******************/
917
918double rawthresh(struct particleMC_t *Isobar)
919{
920 int i;
921 double rmassThresh=0.0;
922
923 if(Isobar->nchildren){
924 for(i=0; i < Isobar->nchildren;i++){
925 if (Isobar->child[i]->mass <0)/* it is not known now */
926 rmassThresh += rawthresh(Isobar->child[i]);
927 else /* it is known now */
928 rmassThresh += Isobar->child[i]->mass;
929 }
930 }else{
931 rmassThresh=Isobar->mass;
932 if(Isobar->mass <0){ /* error */
933 fprintf(stderrstderr,"Error!: Isobar->mass <0 for Isobar with no children\nExit\n");
934 exit(-1);
935 }
936 }
937 return rmassThresh;
938}
939
940
941/**********************************
942 *
943 * decay(Isobar)
944 *
945 * Decay the isobar into its children
946 * and then repeat to decay each
947 * child isobar.
948 *************************************/
949
950void decay(struct particleMC_t *Isobar)
951{
952 int i;
953 double breakup_p,theta,phi;
954
955
956
957 if(Isobar->nchildren>0)
958 {
959 breakup_p = CMmomentum(Isobar->mass,
960 Isobar->child[0]->mass,
961 Isobar->child[1]->mass);
962 theta = acos(randm(-0.9999, 0.9999));
963 phi = randm(-1*M_PI3.14159265358979323846,M_PI3.14159265358979323846);
964 Isobar->child[0]->p =
965 polarMake4v(breakup_p,theta,phi,Isobar->child[0]->mass);
966 Isobar->child[1]->p =
967 polarMake4v(breakup_p,(M_PI3.14159265358979323846 - theta),(M_PI3.14159265358979323846 + phi),Isobar->child[1]->mass);
968 for(i=0;i<Isobar->nchildren;i++)
969 decay(Isobar->child[i]);
970 }
971}
972
973
974/********************************
975 *
976 * lorentzFactor()
977 *
978 * Returns the multiplication of
979 * all of the break-up momenta.
980 *
981 *********************************/
982
983void lorentzFactor(double *lf,struct particleMC_t *Isobar)
984{
985 int i;
986
987 if(!(Isobar->nchildren == 0 )){
988 *lf *= v3mag(&(Isobar->child[0]->p.space));
989 for(i=0;i<Isobar->nchildren;i++)
990 lorentzFactor(lf,Isobar->child[i]);
991 }
992
993}
994
995
996/********************************
997 *
998 * boost2lab()
999 *
1000 * o Starting w/ final state particles
1001 * boost to parent's frame.
1002 *
1003 * o Then boost parent & children to
1004 * parent's parent's frame and repeat
1005 * to the lab frame.
1006 *
1007 *********************************/
1008void boost2lab(struct particleMC_t *Isobar)
1009{
1010 int i;
1011 vector4_t beta;
1012
1013 for(i=0;i<Isobar->nchildren;i++)
1014
1015 boost2lab(Isobar->child[i]);
1016
1017 beta = get_beta(&(Isobar->parent->p),PARENTFRAME+1);/* see kinematics.c */
1018 boostFamily(&beta,Isobar);
1019}
1020
1021/********************************
1022 *
1023 * boostFamily()
1024 *
1025 * Boost particle and all children,
1026 * children's children, ...
1027 *
1028 *********************************/
1029void boostFamily(vector4_t *beta,struct particleMC_t *Isobar)
1030{
1031 int j;
1032 boost(beta,&(Isobar->p));
1033 for(j=0;j<Isobar->nchildren;j++)
1034 boostFamily(beta, Isobar->child[j]);
1035}
1036
1037/********************************
1038 *
1039 * boost()
1040 *
1041 * Boost a four vector.
1042 *
1043 *********************************/
1044void boost(vector4_t *beta,vector4_t *vec)
1045{
1046 vector4_t temp;
1047
1048 temp = lorentz(beta,vec);/* see kinematics.c */
1049 vec->t = temp.t;
1050 vec->space.x = temp.space.x;
1051 vec->space.y = temp.space.y;
1052 vec->space.z = temp.space.z;
1053}
1054
1055/********************************
1056 *
1057 * printProduction()
1058 *
1059 * Print out production particles.
1060 *******************************/
1061void printProduction(FILE *fp,struct particleMC_t *Isobar)
1062{
1063 int i;
1064
1065 for(i=0;i<Isobar->nchildren;i++){
1066 if((Isobar->child[i]->flag%10 ) == 1)
1067 printp2ascii(fp,Isobar->child[i]);
1068 printProduction(fp,Isobar->child[i]);
1069 }
1070}
1071
1072/********************************
1073 *
1074 * printFinal()
1075 *
1076 * Print out final state particles
1077 *******************************/
1078void printFinal(FILE *fp,struct particleMC_t *Isobar)
1079{
1080 int i;
1081
1082 for(i=0;i<Isobar->nchildren;i++){
1083 if((Isobar->child[i]->flag/10 ) == 1)
1084 printp2ascii(fp,Isobar->child[i]);
1085 printFinal(fp,Isobar->child[i]);
1086 }
1087}
1088
1089/********************************
1090 *
1091 * printp2ascii()
1092 *
1093 *******************************/
1094void printp2ascii(FILE *fp,struct particleMC_t *Isobar)
1095{
1096 Nprinted++;
1097 if(UseName)
1098 fprintf(fp,"%d %s %lf\n",Nprinted,ParticleType(Isobar->particleID),Isobar->mass);
1099 else
1100 fprintf(fp,"%d %d %lf\n",Nprinted,Isobar->particleID,Isobar->mass);
1101
1102 fprintf(fp," %d %lf %lf %lf %lf\n",Isobar->charge,
1103 Isobar->p.space.x,
1104 Isobar->p.space.y,
1105 Isobar->p.space.z,
1106 Isobar->p.t);
1107
1108}
1109
1110/********************************
1111 *
1112 * printFamily()
1113 *
1114 *******************************/
1115void printFamily(struct particleMC_t *Isobar)
1116{
1117 int j;
1118 printParticle(Isobar);
1119 for(j=0;j<Isobar->nchildren;j++)
1120 printFamily(Isobar->child[j]);
1121}
1122
1123
1124
1125/********************************
1126 *
1127 * printParticle()
1128 *
1129 *******************************/
1130void printParticle(struct particleMC_t *Isobar)
1131{
1132 fprintf(stderrstderr,"Particle ID %s with %d children\n",
1133 ParticleType(Isobar->particleID),Isobar->nchildren);
1134 fprintf(stderrstderr,"four momentum (E,p): %lf %lf %lf %lf\n\n",
1135 Isobar->p.t,
1136 Isobar->p.space.x,
1137 Isobar->p.space.y,
1138 Isobar->p.space.z);
1139}
1140
1141/******************************************************
1142 *
1143 * polarMake4v()
1144 *
1145 * make a four vector given (p,theta,phi) and it's mass
1146 ********************************************************/
1147vector4_t polarMake4v(double p, double theta, double phi, double mass)
1148{
1149 vector4_t temp;
1150
1151 temp.t = sqrt( SQ(mass) + SQ(p));
1152 temp.space.z = p*cos(theta);
1153 temp.space.x = p*sin(theta)*cos(phi);
1154 temp.space.y = p*sin(theta)*sin(phi);
1155
1156 return temp;
1157}
1158
1159/********************************
1160 *
1161 * randm(double low, double high)
1162 *
1163 *******************************/
1164double randm(double low, double high)
1165{
1166 /* Seed the random number generator using:
1167 * int now = time(NULL);
1168 * srand48(now);
1169 */
1170 return ((high - low) * drand48() + low);
1171}
1172
1173#if 0
1174char *ParticleType(Particle_t p)
1175{
1176 static char ret[20];
1177 switch (p) {
1178 case Unknown:
1179 strcpy(ret,"unknown");
1180 break;
1181 case Gamma:
1182 strcpy(ret,"gamma");
1183 break;
1184 case Positron:
1185 strcpy(ret,"positron");
1186 break;
1187 case Electron:
1188 strcpy(ret,"electron");
1189 break;
1190 case Neutrino:
1191 strcpy(ret,"neutrino");
1192 break;
1193 case MuonPlus:
1194 strcpy(ret,"mu+");
1195 break;
1196 case MuonMinus:
1197 strcpy(ret,"mu-");
1198 break;
1199 case Pi0:
1200 strcpy(ret,"pi0");
1201 break;
1202 case PiPlus:
1203 strcpy(ret,"pi+");
1204 break;
1205 case PiMinus:
1206 strcpy(ret,"pi-");
1207 break;
1208 case KLong:
1209 strcpy(ret,"KL");
1210 break;
1211 case KPlus:
1212 strcpy(ret,"K+");
1213 break;
1214 case KMinus:
1215 strcpy(ret,"K-");
1216 break;
1217 case Neutron:
1218 strcpy(ret,"neutron");
1219 break;
1220 case Proton:
1221 strcpy(ret,"proton");
1222 break;
1223 case AntiProton:
1224 strcpy(ret,"pbar");
1225 break;
1226 case KShort:
1227 strcpy(ret,"Ks");
1228 break;
1229 case Eta:
1230 strcpy(ret,"eta");
1231 break;
1232 case Lambda:
1233 strcpy(ret,"lambda");
1234 break;
1235 case SigmaPlus:
1236 strcpy(ret,"sigma+");
1237 break;
1238 case Sigma0:
1239 strcpy(ret,"sigma0");
1240 break;
1241 case SigmaMinus:
1242 strcpy(ret,"sigma-");
1243 break;
1244 case Xi0:
1245 strcpy(ret,"Xi0");
1246 break;
1247 case XiMinus:
1248 strcpy(ret,"Xi-");
1249 break;
1250 case OmegaMinus:
1251 strcpy(ret,"omega-");
1252 break;
1253 case AntiNeutron:
1254 strcpy(ret,"nbar");
1255 break;
1256
1257 case AntiLambda:
1258 strcpy(ret,"lambdabar");
1259 break;
1260 case AntiSigmaMinus:
1261 strcpy(ret,"sigmabar-");
1262 break;
1263 case AntiSigma0:
1264 strcpy(ret,"sigmabar0");
1265 break;
1266 case AntiSigmaPlus:
1267 strcpy(ret,"sigmabar+");
1268 break;
1269 case AntiXi0:
1270 strcpy(ret,"Xibar0");
1271 break;
1272 case AntiXiPlus:
1273 strcpy(ret,"Xibar+");
1274 break;
1275 case AntiOmegaPlus:
1276 strcpy(ret,"omegabar+");
1277 break;
1278 case Rho0:
1279 strcpy(ret,"rho0");
1280 break;
1281 case RhoPlus:
1282 strcpy(ret,"rho+");
1283 break;
1284 case RhoMinus:
1285 strcpy(ret,"rho;");
1286 break;
1287 case omega:
1288 strcpy(ret,"omega");
1289 break;
1290 case EtaPrime:
1291 strcpy(ret,"etaprime");
1292 break;
1293 case phiMeson:
1294 strcpy(ret,"phi");
1295 break;
1296 default:
1297 sprintf(ret,"type(%d)",(int)p);
1298 break;
1299 }
1300 return(ret);
1301}
1302
1303#endif
1304/*
1305 ***********************
1306 * *
1307 * END OF FILE *
1308 * *
1309 ***********************
1310 */
1311
1312
1313
1314