HPS-MC
 
Loading...
Searching...
No Matches
beam_coords.cc
Go to the documentation of this file.
1#include <stdio.h>
2#include <stdlib.h>
3#include <math.h>
4#include <string.h>
5#include <stdhep_util.hh>
6
7#include <gsl/gsl_rng.h>
8#include <gsl/gsl_randist.h>
9
10#include <unistd.h>
11
12void rotate_entry(stdhep_entry *entry, double theta_x, double theta_y)
13{
14 double px = entry->phep[0];
15 double py = entry->phep[1];
16 double pz = entry->phep[2];
17 double vx = entry->vhep[0];
18 double vy = entry->vhep[1];
19 double vz = entry->vhep[2];
20
21 // rotate about x-axis (rotation in y-z)
22 entry->phep[1] = py*cos(theta_x) - pz*sin(theta_x);
23 entry->phep[2] = py*sin(theta_x) + pz*cos(theta_x);
24
25 py = entry->phep[1];
26 pz = entry->phep[2];
27
28 // rotate about y-axis (rotation in z-x)
29 entry->phep[0] = px*cos(theta_y) + pz*sin(theta_y);
30 entry->phep[2] = pz*cos(theta_y) - px*sin(theta_y);
31
32 // rotate about x-axis (rotation in y-z)
33 entry->vhep[1] = vy*cos(theta_x) - vz*sin(theta_x);
34 entry->vhep[2] = vy*sin(theta_x) + vz*cos(theta_x);
35
36 vy = entry->vhep[1];
37 vz = entry->vhep[2];
38
39 // rotate about y-axis (rotation in z-x)
40 entry->vhep[0] = vx*cos(theta_y) + vz*sin(theta_y);
41 entry->vhep[2] = vz*cos(theta_y) - vx*sin(theta_y);
42}
43
48int main(int argc, char** argv)
49{
50 int nevhep;
51 vector<stdhep_entry> new_event;
52
53 int rseed = 0;
54
55 double theta_x = 0;
56 double theta_y = 0.0305;
57 double theta_z = 0;
58 double sigma_x = 0;
59 double sigma_y = 0;
60 double target_x = 0;
61 double target_y = 0;
62 double target_z = 0;
63
64 int c;
65
66 while ((c = getopt(argc,argv,"hs:x:y:u:v:w:X:Y:Z:")) != -1)
67 switch (c)
68 {
69 case 'h':
70 printf("-h: print this help\n");
71 printf("-s: RNG seed\n");
72 printf("-x: beam sigma_x in mm; default = 0\n");
73 printf("-y: beam sigma_y in mm; default = 0\n");
74 printf("-u: beam x rotation in radians; default = 0\n");
75 printf("-v: beam y rotation in radians; default = 0.0305 rad\n");
76 printf("-w: beam z rotation in radians; default = 0\n");
77 printf("-X: target X in mm; default = 0\n");
78 printf("-Y: target Y in mm; default = 0\n");
79 printf("-Z: target Z in mm; default = 0\n");
80 return(0);
81 break;
82 case 's':
83 rseed = atoi(optarg);
84 break;
85 case 'u':
86 theta_x = atof(optarg);
87 break;
88 case 'v':
89 theta_y = atof(optarg);
90 break;
91 case 'w':
92 theta_z = atof(optarg);
93 break;
94 case 'x':
95 sigma_x = atof(optarg);
96 break;
97 case 'y':
98 sigma_y = atof(optarg);
99 break;
100 case 'X':
101 target_x = atof(optarg);
102 break;
103 case 'Y':
104 target_y = atof(optarg);
105 break;
106 case 'Z':
107 target_z = atof(optarg);
108 break;
109 case '?':
110 printf("Invalid option or missing option argument; -h to list options\n");
111 return(1);
112 default:
113 abort();
114 }
115
116 if (argc-optind < 2)
117 {
118 printf("<input stdhep filename> <output stdhep filename>\n");
119 return 1;
120 }
121
123 const gsl_rng_type * T;
124 gsl_rng * r;
125 gsl_rng_env_setup();
126
127 T = gsl_rng_mt19937;
128 r = gsl_rng_alloc (T);
129 gsl_rng_set(r, rseed);
130
131 int n_events;
132 int istream = 0;
133 int ostream = 1;
134
135 n_events = open_read(argv[optind], istream);
136
137 open_write(argv[optind+1], ostream, n_events);
138
139 printf("Rotating by %f radians in X, %f radians in Y, and %f radians in Z;\nbeam size %f mm in X and %f mm in Y;\ntarget shifted by (X,Y,Z)=(%f,%f,%f) mm\n", theta_x, theta_y, theta_z, sigma_x, sigma_y, target_x, target_y, target_z);
140
141 while (true) {
142 if (!read_next(istream)) {
143 close_read(istream);
144 close_write(ostream);
145 return(0);
146 }
147 nevhep = read_stdhep(&new_event);
148
149 double shift_x = 0.0, shift_y = 0.0;
150 if (sigma_x > 0) shift_x = sigma_x * gsl_ran_gaussian(r, sigma_x);
151 if (sigma_y > 0) shift_y = sigma_y * gsl_ran_gaussian(r, sigma_y);
152
153 double temp_x, temp_y;
154 temp_x = shift_x * cos(theta_z) - shift_y * sin(theta_z);
155 temp_y = shift_y * cos(theta_z) + shift_x * sin(theta_z);
156
157 shift_x = temp_x;
158 shift_y = temp_y;
159
160 for (int i = 0; i < new_event.size(); i++) {
161 rotate_entry(&(new_event[i]), theta_x, theta_y);
162 new_event[i].vhep[2] += target_z;
163 new_event[i].vhep[0] += (shift_x + target_x);
164 new_event[i].vhep[1] += (shift_y + target_y);
165 }
166
167 write_stdhep(&new_event, nevhep);
168 write_file(ostream);
169 }
170}
171
int main(int argc, char **argv)
void rotate_entry(stdhep_entry *entry, double theta_x, double theta_y)
void open_write(char *filename, int ostream, int n_events)
int open_read(char *filename, int istream, int n_events)
void write_stdhep(vector< stdhep_entry > *new_event, int nevhep)
int read_stdhep(vector< stdhep_entry > *new_event)
bool read_next(int istream)
void close_read(int istream)
void write_file(int ostream)
void close_write(int ostream)
Definition stdhep_util.hh:7
double vhep[4]
double phep[5]