Bug Summary

File:programs/Simulation/genr8/genr8.c
Location:line 335, column 13
Description:Value stored to 'pline' during its initialization is never read

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