1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
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 | |
49 | |
50 | |
51 | |
52 | |
53 | struct 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 | |
68 | |
69 | |
70 | |
71 | |
72 | |
73 | |
74 | int Debug = 0; |
75 | int Nprinted =0; |
76 | int PrintProduction=0; |
77 | int PrintRecoil=0; |
78 | double MassHighBW; |
79 | int UseName=0; |
80 | int FIRST_EVENT=1; |
81 | int PrintFlag=10; |
82 | int WriteAscii=0; |
83 | int runNo=9000; |
84 | int NFinalParts=0; |
85 | unsigned int RandomSeed=0; |
86 | int UseCurrentTimeForRandomSeed = TRUE1; |
87 | |
88 | |
89 | |
90 | |
91 | |
92 | |
93 | |
94 | |
95 | double rawthresh(struct particleMC_t *Isobar); |
96 | int decay(struct particleMC_t *Isobar); |
97 | int boost2lab(struct particleMC_t *Isobar); |
98 | int boostFamily(vector4_t *beta,struct particleMC_t *Isobar); |
99 | int boost(vector4_t *beta,vector4_t *vec); |
100 | int printParticle(struct particleMC_t *Isobar); |
101 | int printParticle(struct particleMC_t *Isobar); |
102 | vector4_t polarMake4v(double p, double theta, double phi, double mass); |
103 | double randm(double low, double high); |
104 | int printProduction(FILE *fp,struct particleMC_t *Isobar); |
105 | int printFinal(FILE *fp,struct particleMC_t *Isobar); |
106 | int printp2ascii(FILE *fp,struct particleMC_t *Isobar); |
107 | int setMass(struct particleMC_t *Isobar); |
108 | int initMass(struct particleMC_t *Isobar); |
109 | char *ParticleType(Particle_t p); |
110 | |
111 | |
112 | |
113 | |
114 | |
115 | |
116 | |
117 | |
118 | |
119 | int 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 | |
129 | fprintf(stderrstderr,"\t-P save flag= 11 & 01 events(default saves 11 & 10 events) \n"); |
130 | |
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 | |
144 | |
145 | |
146 | |
147 | |
148 | main(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 | |
169 | CM.child[0]= X; |
170 | CM.child[1]= Y; |
171 | |
172 | |
173 | if (argc == 1){ |
174 | PrintUsage(argv[0]); |
175 | exit (0); |
176 | } |
177 | else { |
178 | |
179 | |
180 | |
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 | |
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 | |
250 | |
251 | |
252 | |
253 | |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
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 | } |
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 | } |
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 | |
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 | } |
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 | |
325 | |
326 | |
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++){ |
332 | haveChildren=TRUE1; |
333 | isacomment=TRUE1; |
334 | while(isacomment==TRUE1){ |
335 | char *pline = fgets(line,sizeof(line),stdinstdin); |
336 | token=strtok(line," "); |
337 | if(!(*token == '%')) |
338 | isacomment=FALSE0; |
339 | } |
340 | if(!(*token == '*')) |
341 | part = atoi(token); |
342 | token=strtok(NULL((void*)0)," "); |
343 | if(!(*token == '*')){ |
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 == '*')) { |
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 == '*')) { |
358 | prnt = atoi(token); |
359 | particle[part].parent = &(particle[prnt]); |
360 | } |
361 | token=strtok(NULL((void*)0)," "); |
362 | if(!(*token == '*')) |
363 | |
364 | particle[part].particleID = atoi(token); |
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{ |
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 { |
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 | |
405 | |
406 | |
407 | |
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 | |
429 | |
430 | |
431 | |
432 | |
433 | |
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 | |
445 | MassHighBW = sqrt_s; |
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 | |
468 | |
469 | |
470 | |
471 | |
472 | |
473 | do{ |
474 | l0: imassc2=0; |
475 | if(!(X->width<0)) |
476 | do{ |
477 | initMass(X); |
478 | setMass(X); |
479 | |
480 | |
481 | |
482 | |
483 | |
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); |
490 | l1: imassc=setChildrenMass(X->child[i]); |
491 | if (Debug) fprintf(stderrstderr,"Return from setChildrenMass X... %d %d \n",i,imassc); |
492 | |
493 | |
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{ |
507 | fprintf(stderrstderr,"Cannot use a negative width!\n"); |
508 | exit(-1); |
509 | } |
510 | |
511 | |
512 | |
513 | if(!(Y->width<0)) |
514 | do{ |
515 | initMass(Y); |
516 | setMass(Y); |
517 | |
518 | |
519 | |
520 | |
521 | |
522 | |
523 | |
524 | for(i=0;i<Y->nchildren;i++) |
525 | { |
526 | if (Debug) fprintf(stderrstderr,"calling setChildrenMass Y... %d \n",i); |
527 | l2: imassc=setChildrenMass(Y->child[i]); |
528 | if (Debug) fprintf(stderrstderr,"Return from setChildrenMass Y... %d %d \n",i,imassc); |
529 | |
530 | |
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{ |
538 | fprintf(stderrstderr,"Cannot use a negative width!\n"); |
539 | exit(-1); |
540 | } |
541 | }while(sqrt_s < X->mass + Y->mass); |
542 | |
543 | |
544 | |
545 | |
546 | |
547 | |
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 | |
563 | |
564 | |
565 | |
566 | |
567 | |
568 | |
569 | |
570 | } else{ |
571 | |
572 | t_min=0.4; |
573 | t_max=10.0; |
| Value stored to 't_max' is never read |
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 | |
602 | |
603 | |
604 | |
605 | |
606 | |
607 | |
608 | |
609 | |
610 | |
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 | |
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; |
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 | |
640 | |
641 | |
642 | |
643 | |
644 | |
645 | |
646 | |
647 | |
648 | |
649 | |
650 | ngenerated++; |
651 | |
652 | if(lf > randm(0.0,lfmax) ){ |
653 | |
654 | naccepted++; |
655 | boost2lab(X); |
656 | boost2lab(Y); |
657 | |
658 | |
659 | |
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 | |
671 | |
672 | fprintf(fout,"%d %d %d\n",runNo,naccepted, NFinalParts); |
673 | |
674 | |
675 | |
676 | |
677 | |
678 | |
679 | |
680 | |
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 | } |
711 | } |
712 | |
713 | fprintf(stderrstderr, |
714 | "Max Lorentz Factor:%lf Events generated:%d Events acepted:%d\n\n", |
715 | lfmax,ngenerated,naccepted); |
716 | |
717 | |
718 | |
719 | fflush(fout); |
720 | fclose(fout); |
721 | |
722 | } |
723 | |
724 | |
725 | |
726 | |
727 | |
728 | |
729 | |
730 | |
731 | |
732 | |
733 | |
734 | |
735 | int 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 | |
758 | |
759 | |
760 | |
761 | |
762 | |
763 | |
764 | int setChildrenMass(struct particleMC_t *Isobar) |
765 | { |
766 | int i, imassc=0; |
767 | |
768 | |
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 | |
777 | l3: 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 | |
790 | if (Debug) fprintf(stderrstderr,"final call ... \n"); |
791 | |
792 | |
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 | |
809 | |
810 | |
811 | |
812 | |
813 | |
814 | |
815 | int 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] ? |
827 | rawthresh((Isobar->parent)->child[1]) : |
828 | rawthresh((Isobar->parent)->child[0]) ; |
829 | |
830 | |
831 | |
832 | hightail = (Isobar->parent->mass) - thresH2 ; |
833 | |
834 | |
835 | |
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 | |
855 | |
856 | |
857 | } |
858 | } |
859 | |
860 | |
861 | |
862 | |
863 | |
864 | |
865 | |
866 | |
867 | |
868 | int 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; |
877 | initMass(Isobar->child[i]); |
878 | } |
879 | } |
880 | |
881 | |
882 | |
883 | |
884 | |
885 | |
886 | |
887 | |
888 | |
889 | double 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) |
897 | rmassThresh += rawthresh(Isobar->child[i]); |
898 | else |
899 | rmassThresh += Isobar->child[i]->mass; |
900 | } |
901 | }else{ |
902 | rmassThresh=Isobar->mass; |
903 | if(Isobar->mass <0){ |
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 | |
915 | |
916 | |
917 | |
918 | |
919 | |
920 | |
921 | int 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 | |
948 | |
949 | |
950 | |
951 | |
952 | |
953 | |
954 | int 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 | |
970 | |
971 | |
972 | |
973 | |
974 | |
975 | |
976 | |
977 | |
978 | |
979 | int 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); |
989 | boostFamily(&beta,Isobar); |
990 | } |
991 | |
992 | |
993 | |
994 | |
995 | |
996 | |
997 | |
998 | |
999 | |
1000 | int 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 | |
1011 | |
1012 | |
1013 | |
1014 | |
1015 | int boost(vector4_t *beta,vector4_t *vec) |
1016 | { |
1017 | vector4_t temp; |
1018 | |
1019 | temp = lorentz(beta,vec); |
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 | |
1029 | |
1030 | |
1031 | |
1032 | int 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 | |
1046 | |
1047 | |
1048 | |
1049 | int 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 | |
1063 | |
1064 | |
1065 | int 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 | |
1084 | |
1085 | |
1086 | int 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 | |
1099 | |
1100 | |
1101 | int 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 | |
1115 | |
1116 | |
1117 | |
1118 | vector4_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 | |
1133 | |
1134 | |
1135 | double randm(double low, double high) |
1136 | { |
1137 | |
1138 | |
1139 | |
1140 | |
1141 | return ((high - low) * drand48() + low); |
1142 | } |
1143 | |
1144 | #if 0 |
1145 | char *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 | |
1279 | |
1280 | |
1281 | |
1282 | |
1283 | |
1284 | |
1285 | |