source: trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/cermaker.c@ 1602

Last change on this file since 1602 was 1602, checked in by bigongia, 22 years ago
*** empty log message ***
File size: 5.3 KB
Line 
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 */
10typedef 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: */
82typedef 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
92main(short argc, char *argv[])
93{
94float theta, phi; // photons direction
95double u, v, w; // direction cosines of photons
96double source_x, source_y, source_z; // z is above sea level.
97float xcenter = 0., ycenter = 0.;
98double cphz, norm;
99float minradius, maxradius; // area of photon impact points
100long maxevents = 1;
101long event, iphot;
102FILE *out;
103
104float dummy[272];
105char runh[5], rune[5], evte[5];
106CerEventHeader event_head;
107photon cph;
108
109float r, psi, maxpsi;
110
111/*****************************/
112
113if (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
121sscanf (argv[1], "%lf", &source_x);
122sscanf (argv[2], "%lf", &source_y);
123sscanf (argv[3], "%lf", &source_z);
124
125source_z += MAGIC_HEIGHT;
126
127/*
128sscanf (argv[4], "%f", &minradius);
129sscanf (argv[5], "%f", &maxradius);
130sscanf (argv[6], "%f", &maxpsi);
131*/
132
133minradius = 0.; maxradius = 900.; maxpsi = 360.;
134
135if (argc == 5)
136 sscanf (argv[4], "%ld", &maxevents);
137
138
139memset ((char*)dummy,0,sizeof(dummy));
140dummy[0] = 1;
141dummy[2] = -1;
142dummy[3] = 1;
143dummy[4] = 220000;
144
145
146memset ((char*)event_head.EVTH,0,sizeof(event_head));
147event_head.EVTH[0] = 'E';
148event_head.EVTH[1] = 'V';
149event_head.EVTH[2] = 'T';
150event_head.EVTH[3] = 'H';
151event_head.PrimaryID = 1;
152event_head.Etotal = 1.;
153event_head.Theta = atan2(sqrt(source_x*source_x+source_y*source_y), source_z-MAGIC_HEIGHT);
154event_head.Phi = atan2(-source_y,-source_x);
155if (event_head.Phi < 0)
156 event_head.Phi += TWOPI;
157event_head.NumObsLev = 1;
158event_head.HeightLev[0] = 220000;
159event_head.CERENKOVyn = 1;
160event_head.Corsika_version = -1;
161event_head.CBunchSize = 1;
162event_head.CDetInX = 1;
163event_head.CDetInY = 1;
164event_head.COutput = 1;
165event_head.CWaveLower = 290.;
166event_head.CWaveUpper = 600.;
167event_head.Multipl =1;
168event_head.CorePos[0][0] = 0.;
169event_head.CorePos[1][0] = 0.;
170
171memset ((char*)&(cph.w),0,sizeof(cph));
172cph.w = 500;
173cph.t = 1000;
174
175sprintf(evte, "EVTE");
176sprintf(runh, "RUNH");
177sprintf(rune, "RUNE");
178
179out = fopen("cer000001","w");
180if (!out)
181 {
182 printf ("\nERROR: Cannot open output file cer000001\n\n");
183 exit(-1);
184 }
185
186fwrite (runh, strlen(runh), 1, out);
187fwrite (dummy, sizeof(dummy), 1, out);
188
189/* Here one event always means 39 photons! */
190
191for (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
229fwrite (rune, strlen(rune), 1, out);
230fwrite (dummy, sizeof(dummy), 1, out);
231
232fclose(out);
233
234return;
235
236}
237
Note: See TracBrowser for help on using the repository browser.