1 | #include <stdio.h>
|
---|
2 | #include <stdlib.h>
|
---|
3 | #include <math.h>
|
---|
4 |
|
---|
5 | #define DEG2RAD 0.01745329252
|
---|
6 | #define TWOPI 6.283185307
|
---|
7 | #define MAGIC_HEIGHT 220000.
|
---|
8 |
|
---|
9 | /* EVTH from cerfile. See CORSIKA manual for explanations */
|
---|
10 | typedef struct
|
---|
11 | { char EVTH[4];
|
---|
12 | float EvtNumber;
|
---|
13 | float PrimaryID;
|
---|
14 | float Etotal;
|
---|
15 | float Thick0;
|
---|
16 | float FirstTarget;
|
---|
17 | float zFirstInt;
|
---|
18 | float p[3];
|
---|
19 | float Theta;
|
---|
20 | float Phi;
|
---|
21 | float NumRndSeq;
|
---|
22 | float RndData[10][3];
|
---|
23 | float RunNumber;
|
---|
24 | float DateRun;
|
---|
25 | float Corsika_version;
|
---|
26 | float NumObsLev;
|
---|
27 | float HeightLev[10];
|
---|
28 | float SlopeSpec;
|
---|
29 | float ELowLim;
|
---|
30 | float EUppLim;
|
---|
31 | float Ecutoffh;
|
---|
32 | float Ecutoffm;
|
---|
33 | float Ecutoffe;
|
---|
34 | float Ecutoffg;
|
---|
35 | float NFLAIN;
|
---|
36 | float NFLDIF;
|
---|
37 | float NFLPI0;
|
---|
38 | float NFLPIF;
|
---|
39 | float NFLCHE;
|
---|
40 | float NFRAGM;
|
---|
41 | float Bx;
|
---|
42 | float By;
|
---|
43 | float EGS4yn;
|
---|
44 | float NKGyn;
|
---|
45 | float GHEISHAyn;
|
---|
46 | float VENUSyn;
|
---|
47 | float CERENKOVyn;
|
---|
48 | float NEUTRINOyn;
|
---|
49 | float HORIZONTyn;
|
---|
50 | float COMPUTER;
|
---|
51 | float ThetaMin;
|
---|
52 | float ThetaMax;
|
---|
53 | float PhiMin;
|
---|
54 | float PhiMax;
|
---|
55 | float CBunchSize;
|
---|
56 | float CDetInX,CDetInY;
|
---|
57 | float CSpacInX,CSpacInY;
|
---|
58 | float CLenInX,CLenInY;
|
---|
59 | float COutput;
|
---|
60 | float AngleNorthX;
|
---|
61 | float MuonInfo;
|
---|
62 | float StepLength;
|
---|
63 | float CWaveLower;
|
---|
64 | float CWaveUpper;
|
---|
65 | float Multipl;
|
---|
66 | float CorePos[2][20];
|
---|
67 | float SIBYLL[2];
|
---|
68 | float QGSJET[2];
|
---|
69 | float DPMJET[2];
|
---|
70 | float VENUS_cross;
|
---|
71 | float mu_mult_scat;
|
---|
72 | float NKG_range;
|
---|
73 | float EFRCTHN[2];
|
---|
74 | float WMAX[2];
|
---|
75 | float rthin_rmax;
|
---|
76 | float viewcone_angles[2];
|
---|
77 | float dmmy2[119];
|
---|
78 | } CerEventHeader;
|
---|
79 |
|
---|
80 |
|
---|
81 | /* Corsika Photon block: */
|
---|
82 | typedef struct
|
---|
83 | { float w, /* photon wavelength (nm) */
|
---|
84 | x, y, /* (ground) impact point (cm) */
|
---|
85 | u, v, /* direction cosines */
|
---|
86 | t, /* arrival time (ns) */
|
---|
87 | h; /* production height (cm) */
|
---|
88 | } photon;
|
---|
89 |
|
---|
90 | /*****************************/
|
---|
91 |
|
---|
92 | main(short argc, char *argv[])
|
---|
93 | {
|
---|
94 | float theta, phi; // photons direction
|
---|
95 | double u, v, w; // direction cosines of photons
|
---|
96 | double source_x, source_y, source_z; // z is above sea level.
|
---|
97 | float xcenter = 0., ycenter = 0.;
|
---|
98 | double cphz, norm;
|
---|
99 | float minradius, maxradius; // area of photon impact points
|
---|
100 | long maxevents = 1;
|
---|
101 | long event, iphot;
|
---|
102 | FILE *out;
|
---|
103 |
|
---|
104 | float dummy[272];
|
---|
105 | char runh[5], rune[5], evte[5];
|
---|
106 | CerEventHeader event_head;
|
---|
107 | photon cph;
|
---|
108 |
|
---|
109 | float r, psi, maxpsi;
|
---|
110 |
|
---|
111 | /*****************************/
|
---|
112 |
|
---|
113 | if (argc < 4 || argc > 5)
|
---|
114 | {
|
---|
115 | printf ("\nUsage: cermaker source_x(cm) source_y(cm) source_z(cm) [events]\n\n");
|
---|
116 | exit(-1);
|
---|
117 | }
|
---|
118 |
|
---|
119 | /* Source position with respect to telescope! */
|
---|
120 |
|
---|
121 | sscanf (argv[1], "%lf", &source_x);
|
---|
122 | sscanf (argv[2], "%lf", &source_y);
|
---|
123 | sscanf (argv[3], "%lf", &source_z);
|
---|
124 |
|
---|
125 | source_z += MAGIC_HEIGHT;
|
---|
126 |
|
---|
127 | /*
|
---|
128 | sscanf (argv[4], "%f", &minradius);
|
---|
129 | sscanf (argv[5], "%f", &maxradius);
|
---|
130 | sscanf (argv[6], "%f", &maxpsi);
|
---|
131 | */
|
---|
132 |
|
---|
133 | minradius = 0.; maxradius = 900.; maxpsi = 360.;
|
---|
134 |
|
---|
135 | if (argc == 5)
|
---|
136 | sscanf (argv[4], "%ld", &maxevents);
|
---|
137 |
|
---|
138 |
|
---|
139 | memset ((char*)dummy,0,sizeof(dummy));
|
---|
140 | dummy[0] = 1;
|
---|
141 | dummy[2] = -1;
|
---|
142 | dummy[3] = 1;
|
---|
143 | dummy[4] = 220000;
|
---|
144 |
|
---|
145 |
|
---|
146 | memset ((char*)event_head.EVTH,0,sizeof(event_head));
|
---|
147 | event_head.EVTH[0] = 'E';
|
---|
148 | event_head.EVTH[1] = 'V';
|
---|
149 | event_head.EVTH[2] = 'T';
|
---|
150 | event_head.EVTH[3] = 'H';
|
---|
151 | event_head.PrimaryID = 1;
|
---|
152 | event_head.Etotal = 1.;
|
---|
153 | event_head.Theta = atan2(sqrt(source_x*source_x+source_y*source_y), source_z-MAGIC_HEIGHT);
|
---|
154 | event_head.Phi = atan2(-source_y,-source_x);
|
---|
155 | if (event_head.Phi < 0)
|
---|
156 | event_head.Phi += TWOPI;
|
---|
157 | event_head.NumObsLev = 1;
|
---|
158 | event_head.HeightLev[0] = 220000;
|
---|
159 | event_head.CERENKOVyn = 1;
|
---|
160 | event_head.Corsika_version = -1;
|
---|
161 | event_head.CBunchSize = 1;
|
---|
162 | event_head.CDetInX = 1;
|
---|
163 | event_head.CDetInY = 1;
|
---|
164 | event_head.COutput = 1;
|
---|
165 | event_head.CWaveLower = 290.;
|
---|
166 | event_head.CWaveUpper = 600.;
|
---|
167 | event_head.Multipl =1;
|
---|
168 | event_head.CorePos[0][0] = 0.;
|
---|
169 | event_head.CorePos[1][0] = 0.;
|
---|
170 |
|
---|
171 | memset ((char*)&(cph.w),0,sizeof(cph));
|
---|
172 | cph.w = 500;
|
---|
173 | cph.t = 1000;
|
---|
174 |
|
---|
175 | sprintf(evte, "EVTE");
|
---|
176 | sprintf(runh, "RUNH");
|
---|
177 | sprintf(rune, "RUNE");
|
---|
178 |
|
---|
179 | out = fopen("cer000001","w");
|
---|
180 | if (!out)
|
---|
181 | {
|
---|
182 | printf ("\nERROR: Cannot open output file cer000001\n\n");
|
---|
183 | exit(-1);
|
---|
184 | }
|
---|
185 |
|
---|
186 | fwrite (runh, strlen(runh), 1, out);
|
---|
187 | fwrite (dummy, sizeof(dummy), 1, out);
|
---|
188 |
|
---|
189 | /* Here one event always means 39 photons! */
|
---|
190 |
|
---|
191 | for (event = 0; event < maxevents; event++)
|
---|
192 | {
|
---|
193 | event_head.EvtNumber = event+1;
|
---|
194 | fwrite (event_head.EVTH, sizeof(event_head), 1, out);
|
---|
195 |
|
---|
196 | for (iphot = 0; iphot < 39; iphot++)
|
---|
197 | {
|
---|
198 | psi = (double)rand()/(double)RAND_MAX*maxpsi*DEG2RAD;
|
---|
199 | r = sqrt( (double)rand()/(double)RAND_MAX*
|
---|
200 | (maxradius*maxradius-minradius*minradius) +
|
---|
201 | minradius*minradius);
|
---|
202 |
|
---|
203 | cph.x = xcenter+r*cos(psi);
|
---|
204 | cph.y = ycenter+r*sin(psi);
|
---|
205 |
|
---|
206 | cphz = MAGIC_HEIGHT; /* Above sea level */
|
---|
207 |
|
---|
208 |
|
---|
209 | /* Photon direction: */
|
---|
210 | u = (double)cph.x - source_x;
|
---|
211 | v = (double)cph.y - source_y;
|
---|
212 | w = cphz - source_z;
|
---|
213 |
|
---|
214 | norm = sqrt(u*u+v*v+w*w);
|
---|
215 |
|
---|
216 | cph.u = (float)(u/norm);
|
---|
217 | cph.v = (float)(v/norm);
|
---|
218 | w = (float)(w/norm);
|
---|
219 |
|
---|
220 | cph.h = source_z;
|
---|
221 |
|
---|
222 | fwrite(&(cph.w),sizeof(cph),1,out);
|
---|
223 | }
|
---|
224 |
|
---|
225 | fwrite (evte, strlen(evte), 1, out);
|
---|
226 | fwrite (dummy, sizeof(dummy), 1, out);
|
---|
227 | }
|
---|
228 |
|
---|
229 | fwrite (rune, strlen(rune), 1, out);
|
---|
230 | fwrite (dummy, sizeof(dummy), 1, out);
|
---|
231 |
|
---|
232 | fclose(out);
|
---|
233 |
|
---|
234 | return;
|
---|
235 |
|
---|
236 | }
|
---|
237 |
|
---|