Changeset 1535 for trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c

Timestamp:

10/09/02 19:15:28 (22 years ago)

Author:

bigongia

Message:

Version 0.6. Changed output format: added run header, changed event header,
             added ascii parameter files attached at the end of every output
             file to keep all info.

File:

• trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c (modified) (14 diffs)

trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c

@@ -13 +13 @@
  *
  * ----------------------------------------------------------------------
+ *
+ *  Modified: October 2002
+ *  Author:   Abelardo Moralejo
+ *  Version 0.6: introduced a run header and changed event header to keep
+ *  all information from CORSIKA. In addition, now the magic.def, axisdev.dat
+ *  and reflectivity.dat are attached at the end of the reflector output file 
+ *  to keep also all the parameters with which the program is run.
+ *
+ * ----------------------------------------------------------------------
  */
 
@@ -36 +45 @@
 static long writep = 0L;        /*  write ptr (on rflfile)      */
 
-extern CerHeader *cheadp;       /*  var inited in header.c      */
-extern RflHeader *rheadp;       /*  var inited in header.c      */
+extern CerEventHeader *cheadp;  /*  var inited in header.c      */
+extern RflEventHeader *rheadp;  /*  var inited in header.c      */
 
 /*  Prototypes  */
@@ -43 +52 @@
 FILE *GetNextFile(char *cername);
 static int GetEvent(void);
-static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile);
+static int ProcessEvent(CerEventHeader *cheadp, FILE *cerfile, FILE *rflfile);
+void attach(FILE *f1, FILE *f2);
+
 
 FILE *chkf = NULL;
@@ -51 +62 @@
 
 main(void)
-{   long event = 0L;            /*  event counter       */
+{
+    extern char axisdev_filename[256], reflectivity_filename[256];
+    extern char geo_filename[256];
+    FILE *dummy;
+    long event = 0L;            /*  event counter       */
 
     /*  Read init & geometry parms, init files and vars */
@@ -71 +86 @@
     fwrite(FLAG_END_OF_FILE, SIZE_OF_FLAGS, 1, rflfile);
 
-    /*  Close file  */
+    /* Attach parameter files used for running the program: */
+    dummy = fopen (geo_filename,"r");
+    attach (rflfile, dummy);
+    fclose(dummy);
+
+    dummy = fopen (axisdev_filename,"r");
+    attach (rflfile, dummy);
+    fclose(dummy);
+
+    dummy = fopen (reflectivity_filename,"r");
+    attach (rflfile, dummy);
+    fclose(dummy);
+
+    /*  Close reflector output file  */
     Log(RFLF_CLOSE_LOG);
     fclose(rflfile);
@@ -82 +110 @@
   } /*  end of main  */
 
-static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile)
+static int ProcessEvent(CerEventHeader *cheadp, FILE *cerfile, FILE *rflfile)
 {   extern int absorption(float wlen, float height, float theta);
     extern int ph2cph(photon *ph, cphoton *cph);
@@ -92 +120 @@
     extern void makeOmegaI(float theta, float phi);
 
-    char pp[256];
 
     /*  Various counters: phs = absphs + refphs[0..3] + cphs  */
@@ -172 +199 @@
         for (ph=0; ph<PH_IN_DATABLOCK; ph++)
           {
-            if (Photons[ph].w <= 0.) break;
+            /* Added July 2002, AM: check integrity of photon info:
+               Sometimes we found NaN values inside cerXXXXX file.
+             */
+            if (isnan(Photons[ph].w) || isnan(Photons[ph].x) ||
+                isnan(Photons[ph].y) || isnan(Photons[ph].u) ||
+                isnan(Photons[ph].v) || isnan(Photons[ph].t) ||
+                isnan(Photons[ph].h))
+              {
+                printf("Warning: skipped one photon because its data contained Not-a-Number!\n");
+                continue;
+              }
+
+            if (! (Photons[ph].w > 0.) ) break;
 
             CPhotons[cphs].w = Photons[ph].w;
             Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
 
-            /* TEMPORARY FIX, AM Nov 2001: we found that sometimes the value
+            /* AM Nov 2001: we found that sometimes the value
                stored in Photons[ph].w is not correct, and can result in a 
                wavelength beyond 600 nm, which makes the program crash later.
-               Now we force wlen to its expected range:
-               */
-
-            wlen = MIN(MAX(290.,wlen),600.);
-
+
+               AM July 2002: we now simply skip the photon if the wavelength
+               is not in the expected range, which we now take from the corsika
+               event header (just in case we would change it in the future).
+            */
+
+               if (wlen < cheadp->CWaveLower || wlen > cheadp->CWaveUpper)
+                 {
+                   printf("Warning: skipped one photon with strange wavelength: %f nm\n", wlen);
+                   continue;
+                 }
 
             /* ADDED AM May 2002: now we take into account the telescope
@@ -215 +260 @@
                            ph2cph(&Photons[ph], &CPhotons[cphs])))
               refphs[ref_type-1]++; 
+
             else                /*      Photon passed   */
               {
@@ -260 +306 @@
         rheadp->TimeLast  = last;
 
-        /*  Update RflHeader with info on ph/cph nrs and write it  */
+        /*  Update RflEventHeader with info on ph/cph nrs and write it  */
         rheadp->CORSIKAPhs   = phs;
         rheadp->AtmAbsPhs    = absphs;
@@ -268 +314 @@
         rheadp->OutOfChamPhs = refphs[3];
         rheadp->CPhotons     = (long) overflow * NR_OF_CPHOTONS + cphs;
-        fwrite(rheadp, sizeof(RflHeader), 1, rflfile);
-
-/*        for (myloop=0; myloop<sizeof(RflHeader)/4; myloop++)
+        fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile);
+
+/*        for (myloop=0; myloop<sizeof(RflEventHeader)/4; myloop++)
  *        fprintf(chkf, "%e ", *((float *)rheadp+myloop)); 
  *        fputc('\n', chkf);
@@ -284 +330 @@
             fseek(tmpf, 0L, SEEK_SET); /*  Start from the beginning  */
             while (overflow--)
-              { fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
-                fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); }
+              { 
+                fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
+                fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); 
+              }
 
             /*  Reload data in CPhotons    */
@@ -291 +339 @@
 
             /*  Close (and delete) temp file  */
-            fclose(tmpf);  }
+            fclose(tmpf);  
+          }
 
         /*  Write (remaining) cphotons  */
+
         fwrite(CPhotons, sizeof(cphoton), cphs, rflfile);
 
@@ -308 +358 @@
 
 static int GetEvent(void)
-{   int found = FALSE,          /*  event found         */
-      isWrong = FALSE;          /*  cerfile is wrong    */
-      static int newFile = TRUE; /*  if TRUE, check if cerfile is valid */
-
-      do
-        { /* In case the just-opened file is a valid cerfile,
-             starting with a RUNH, loop until a valid event is found:
-             id est with: first_Event <= EvtNumber <= last_Event
-             and             low_Ecut <=  Etotal   <= high_Ecut
-             If the search was successful, "found" is set to TRUE.
-             If there are reading errors, "isWrong" is set to TRUE. */
+{   
+  int found = FALSE,          /*  event found           */
+  isWrong = FALSE;            /*  cerfile is wrong      */
+  static int newFile = TRUE;  /*  if TRUE, check if cerfile is valid    */
+
+  CerRunHeader cRunHead;
+  
+  do
+    {
+      /* In case the just-opened file is a valid cerfile,
+         starting with a RUNH, loop until a valid event is found:
+         id est with: first_Event <= EvtNumber <= last_Event
+         and             low_Ecut <=  Etotal   <= high_Ecut
+         If the search was successful, "found" is set to TRUE.
+         If there are reading errors, "isWrong" is set to TRUE. */
 
         if (newFile
-         && (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile)
+         && (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile)
              || strncmp(cheadp->EVTH, "RUNH", 4)))
         {   isWrong = TRUE;   }
 
         else do
-        {   /*  Push fileptr in case it is a "EVTH"  */
-            readp = ftell(cerfile);
-
-            /*  Error: exit loop  */
-            if (1 != fread(cheadp, sizeof(CerHeader), 1, cerfile))
+        {
+          if (newFile)
+            {
+              /* Write Reflector "run header" (one per cer file!): */
+              memcpy(&cRunHead, cheadp, sizeof(CerEventHeader));
+              fwrite(&cRunHead, sizeof(CerRunHeader), 1, rflfile);
+            }
+
+          /*  Push fileptr in case it is a "EVTH"  */
+          readp = ftell(cerfile);
+
+          /*  Error: exit loop  */
+          if (1 != fread(cheadp, sizeof(CerEventHeader), 1, cerfile))
             {   isWrong = TRUE;    break;   }
 
-            /*  Ok: set found at TRUE and exit loop  */
-            if (strncmp(cheadp->EVTH, "EVTH", 4) == 0
-             && first_Event <= (long)cheadp->EvtNumber
-             &&                (long)cheadp->EvtNumber <= last_Event
-             &&    low_Ecut <= cheadp->Etotal
-             &&                cheadp->Etotal <= high_Ecut)
+          /*  Ok: set found at TRUE and exit loop  */
+          if (strncmp(cheadp->EVTH, "EVTH", 4) == 0
+              && first_Event <= (long)cheadp->EvtNumber
+              &&                       (long)cheadp->EvtNumber <= last_Event
+              &&    low_Ecut <= cheadp->Etotal
+              &&                       cheadp->Etotal <= high_Ecut)
             {   found = TRUE;   }
 
-            /*  File is finished: exit loop  */
-            else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0)
-                break;
+          /*  File is finished: exit loop  */
+          else if (strncmp(cheadp->EVTH, "RUNE", 4) == 0)
+            break;
 
         } while(found == FALSE);