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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#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 < 4 || argc > 5)
  {
    printf ("\nUsage: cermaker source_x(cm) source_y(cm) source_z(cm) [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);
*/

minradius = 0.; maxradius = 900.; maxpsi = 360.;

if (argc == 5)
  sscanf (argv[4], "%ld", &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;

}