Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/Makefile (revision 1601) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/Makefile (revision 1601) @@ -0,0 +1,4 @@ +cermaker: cermaker.o + cc -o cermaker -lm cermaker.o +.c.o: + cc -g -c cermaker.c Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/cermaker.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/cermaker.c (revision 1601) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/tester/cermaker.c (revision 1601) @@ -0,0 +1,232 @@ +#include +#include +#include + +#define DEG2RAD 0.01745329252 +#define TWOPI 6.283185307 +#define MAGIC_HEIGHT 220000. + +/* EVTH from cerfile. See CORSIKA manual for explanations */ +typedef struct +{ char EVTH[4]; + float EvtNumber; + float PrimaryID; + float Etotal; + float Thick0; + float FirstTarget; + float zFirstInt; + float p[3]; + float Theta; + float Phi; + float NumRndSeq; + float RndData[10][3]; + float RunNumber; + float DateRun; + float Corsika_version; + float NumObsLev; + float HeightLev[10]; + float SlopeSpec; + float ELowLim; + float EUppLim; + float Ecutoffh; + float Ecutoffm; + float Ecutoffe; + float Ecutoffg; + float NFLAIN; + float NFLDIF; + float NFLPI0; + float NFLPIF; + float NFLCHE; + float NFRAGM; + float Bx; + float By; + float EGS4yn; + float NKGyn; + float GHEISHAyn; + float VENUSyn; + float CERENKOVyn; + float NEUTRINOyn; + float HORIZONTyn; + float COMPUTER; + float ThetaMin; + float ThetaMax; + float PhiMin; + float PhiMax; + float CBunchSize; + float CDetInX,CDetInY; + float CSpacInX,CSpacInY; + float CLenInX,CLenInY; + float COutput; + float AngleNorthX; + float MuonInfo; + float StepLength; + float CWaveLower; + float CWaveUpper; + float Multipl; + float CorePos[2][20]; + float SIBYLL[2]; + float QGSJET[2]; + float DPMJET[2]; + float VENUS_cross; + float mu_mult_scat; + float NKG_range; + float EFRCTHN[2]; + float WMAX[2]; + float rthin_rmax; + float viewcone_angles[2]; + float dmmy2[119]; +} CerEventHeader; + + +/* Corsika Photon block: */ +typedef struct +{ float w, /* photon wavelength (nm) */ + x, y, /* (ground) impact point (cm) */ + u, v, /* direction cosines */ + t, /* arrival time (ns) */ + h; /* production height (cm) */ +} photon; + +/*****************************/ + +main(short argc, char *argv[]) +{ +float theta, phi; // photons direction +double u, v, w; // direction cosines of photons +double source_x, source_y, source_z; // z is above sea level. +float xcenter = 0., ycenter = 0.; +double cphz, norm; +float minradius, maxradius; // area of photon impact points +long maxevents = 1; +long event, iphot; +FILE *out; + +float dummy[272]; +char runh[5], rune[5], evte[5]; +CerEventHeader event_head; +photon cph; + +float r, psi, maxpsi; + +/*****************************/ + +if (argc < 7 || argc > 8) + { + printf ("\nUsage: cermaker source_x(cm) source_y(cm) source_z(cm) minradius(cm) maxradius(cm) maxpsi(deg) [events]\n\n"); + exit(-1); + } + +/* Source position with respect to telescope! */ + +sscanf (argv[1], "%lf", &source_x); +sscanf (argv[2], "%lf", &source_y); +sscanf (argv[3], "%lf", &source_z); + +source_z += MAGIC_HEIGHT; + +sscanf (argv[4], "%f", &minradius); +sscanf (argv[5], "%f", &maxradius); +sscanf (argv[6], "%f", &maxpsi); + +if (argc == 8) + sscanf (argv[7], "%d", &maxevents); + +memset ((char*)dummy,0,sizeof(dummy)); +dummy[0] = 1; +dummy[2] = -1; +dummy[3] = 1; +dummy[4] = 220000; + + +memset ((char*)event_head.EVTH,0,sizeof(event_head)); +event_head.EVTH[0] = 'E'; +event_head.EVTH[1] = 'V'; +event_head.EVTH[2] = 'T'; +event_head.EVTH[3] = 'H'; +event_head.PrimaryID = 1; +event_head.Etotal = 1.; +event_head.Theta = atan2(sqrt(source_x*source_x+source_y*source_y), source_z-MAGIC_HEIGHT); +event_head.Phi = atan2(-source_y,-source_x); +if (event_head.Phi < 0) + event_head.Phi += TWOPI; +event_head.NumObsLev = 1; +event_head.HeightLev[0] = 220000; +event_head.CERENKOVyn = 1; +event_head.Corsika_version = -1; +event_head.CBunchSize = 1; +event_head.CDetInX = 1; +event_head.CDetInY = 1; +event_head.COutput = 1; +event_head.CWaveLower = 290.; +event_head.CWaveUpper = 600.; +event_head.Multipl =1; +event_head.CorePos[0][0] = 0.; +event_head.CorePos[1][0] = 0.; + +memset ((char*)&(cph.w),0,sizeof(cph)); +cph.w = 500; +cph.t = 1000; + +sprintf(evte, "EVTE"); +sprintf(runh, "RUNH"); +sprintf(rune, "RUNE"); + +out = fopen("cer000001","w"); +if (!out) + { + printf ("\nERROR: Cannot open output file cer000001\n\n"); + exit(-1); + } + +fwrite (runh, strlen(runh), 1, out); +fwrite (dummy, sizeof(dummy), 1, out); + +/* Here one event always means 39 photons! */ + +for (event = 0; event < maxevents; event++) + { + event_head.EvtNumber = event+1; + fwrite (event_head.EVTH, sizeof(event_head), 1, out); + + for (iphot = 0; iphot < 39; iphot++) + { + psi = (double)rand()/(double)RAND_MAX*maxpsi*DEG2RAD; + r = sqrt( (double)rand()/(double)RAND_MAX* + (maxradius*maxradius-minradius*minradius) + + minradius*minradius); + + cph.x = xcenter+r*cos(psi); + cph.y = ycenter+r*sin(psi); + + cphz = MAGIC_HEIGHT; /* Above sea level */ + + + /* Photon direction: */ + u = (double)cph.x - source_x; + v = (double)cph.y - source_y; + w = cphz - source_z; + + norm = sqrt(u*u+v*v+w*w); + + cph.u = (float)(u/norm); + cph.v = (float)(v/norm); + w = (float)(w/norm); + + cph.h = source_z; + + fwrite(&(cph.w),sizeof(cph),1,out); + } + + fwrite (evte, strlen(evte), 1, out); + fwrite (dummy, sizeof(dummy), 1, out); + } + +fwrite (rune, strlen(rune), 1, out); +fwrite (dummy, sizeof(dummy), 1, out); + +fclose(out); + +return; + +} +