Ignore:
Timestamp:
07/24/02 15:36:10 (22 years ago)
Author:
bigongia
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r870 r1431  
    2424#include "header.h"
    2525
     26#define MAX(x,y) ((x)>(y)? (x) : (y))
     27#define MIN(x,y) ((x)>(y)? (y) : (x))
     28
    2629FILE *rflfile = NULL,           /*  reflector (output) file     */
    27      *cerfile = NULL;           /*  current cerfile             */
     30*cerfile = NULL;                /*  current cerfile             */
    2831
    2932cphoton *CPhotons = NULL;                 /*  cphoton array     */
    … …  
    4245static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile);
    4346
    44 FILE *chkf = NULL;     /**********************/
    45 long myloop;     /**********************/
    46 
    47 void main(void)
     47FILE *chkf = NULL;
     48long myloop;
     49
     50Event_end* evt_end;
     51
     52main(void)
    4853{   long event = 0L;            /*  event counter       */
    4954
    5055    /*  Read init & geometry parms, init files and vars */
    5156    init(NULL);
    52     chkf = fopen("check", "w");     /**********************/
     57
     58    /*    chkf = fopen("check", "w");   */
    5359
    5460    /*  Processing loop */
    5561    while(event < max_Events && GetEvent())
    56     {   if (ProcessEvent(cheadp, cerfile, rflfile)) event++;
    57         if (event % SHOW_ME == 0) Log(INFO_EVENT_LOG,
    58             cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal);  }
     62      {
     63        if (ProcessEvent(cheadp, cerfile, rflfile))
     64          event++;
     65        if (event % SHOW_ME == 0)
     66          Log(INFO_EVENT_LOG,
     67              cheadp->RunNumber, cheadp->EvtNumber, cheadp->Etotal); 
     68      }
    5969
    6070    /*  Writing final flags  */
    … …  
    6474    Log(RFLF_CLOSE_LOG);
    6575    fclose(rflfile);
    66     fclose(chkf);     /**********************/
     76    /*    fclose(chkf);   */
    6777
    6878    /*  Clean memory and exit  */
    … …  
    7080    Message(RFL__EXIT__MSG, QUOTE(PROGRAM), QUOTE(VERSION));
    7181
    72 }  /*  end of main  */
     82  } /*  end of main  */
    7383
    7484static int ProcessEvent(CerHeader *cheadp, FILE *cerfile, FILE *rflfile)
    … …  
    8292    extern void makeOmegaI(float theta, float phi);
    8393
     94    char pp[256];
    8495
    8596    /*  Various counters: phs = absphs + refphs[0..3] + cphs  */
    86     long phs,           /*  Number of incoming photons  */
    87          absphs,        /*  Photons absorbed            */
    88          refphs[4],     /*  Photons not reflected       */
    89          cphs;          /*  Number of cphotons          */
    90 
    91     FILE *tmpf=NULL;    /*  Temp fp to handle o/f       */
    92     size_t read;        /*  items read: != 1 on error   */
    93     int overflow=0,     /*  how many o/f on cphs        */
    94         ref_type,       /*  ret value from reflection   */
    95         ph;             /*  photon number               */
    96 
    97     float first = 1e8f, /*  Photon arrival times        */
    98           last  = 0;
     97    long phs,                   /*  Number of incoming photons  */
     98    absphs,                     /*  Photons absorbed            */
     99    refphs[4],                  /*  Photons not reflected       */
     100    cphs;                       /*  Number of cphotons          */
     101
     102    FILE *tmpf=NULL;            /*  Temp fp to handle o/f       */
     103    size_t read;                /*  items read: != 1 on error   */
     104    int overflow=0,             /*  how many o/f on cphs        */
     105    ref_type,                   /*  ret value from reflection   */
     106    ph;                         /*  photon number               */
     107
     108    float first = 1e8f,         /*  Photon arrival times        */
     109    last  = 0;
    99110
    100111    /*  photon quantities  */
    101     float wlen,         /*  photon wavelength   */
    102           theta;        /*  photon zenith angle */
     112    float wlen,                 /*  photon wavelength   */
     113    theta;                      /*  photon zenith angle */
    103114
    104115    /*  Reset counters and set fileptr  */
    … …  
    111122    /*  Calculate OmegaCT matrices  */
    112123    if (is_Fixed_Target)
    113     {   makeOmega(fixed_Theta, fixed_Phi);
     124      { makeOmega(fixed_Theta, fixed_Phi);
    114125        makeOmegaI(fixed_Theta, fixed_Phi);   }
    115126    else
    116     {   makeOmega(cheadp->Theta, cheadp->Phi);
     127      { makeOmega(cheadp->Theta, cheadp->Phi);
    117128        makeOmegaI(cheadp->Theta, cheadp->Phi);   }
    118129    memcpy( OmegaCT, Omega, 9*sizeof(float) );
    119130    memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
    120131
     132    /* ADDED AM May 2002: now we take into account the telescope
     133     * position chosen in the Corsika input card via the CERTEL
     134     * option, which must be supplied also in the reflector input
     135     * card with the option "telescope_position x y". Otherwise
     136     * x = y = 0 is assumed, which is the standard mode of
     137     * generation for MAGIC.
     138     * Here we change rheadp so that the CorePos written to the
     139     * .rfl file is referred to the telescope position. However,
     140     * I believe that the original "CorePos" vector points
     141     * from the core to the origin of coordinates of Corsika, and
     142     * therefore after the subtraction of (Telescope_x,Telescope_y)
     143     * the resulting vector CorePos  points from the core to the
     144     * telescope!
     145     */
     146
     147    rheadp->CorePos[0][0] += Telescope_x;
     148    rheadp->CorePos[1][0] += Telescope_y;
     149
    121150    /*  Loop on data blocks  */
    122     while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile))
    123       &&  strncmp((char *)Photons, "EVTE", 4))
    124 
    125     /*  Loop on phs inside block: exit when wlength is 0  */
    126     {   for (ph=0; ph<PH_IN_DATABLOCK; ph++)
    127         {   if (Photons[ph].w == 0) break;
     151
     152    while(1 == (read = fread(Photons, sizeof(Photons), 1, cerfile)))
     153      /*  Loop on phs inside block: exit when wlength is 0  */
     154      {
     155        /* If "event end" flag is found, read relevant quantities
     156         * from event end subblock (added June 2002, AM):
     157         */
     158
     159        if (strncmp((char *)Photons, "EVTE", 4) == 0 )
     160          {
     161            evt_end = (Event_end*) Photons;
     162            rheadp->longi_Nmax = evt_end->longi_Nmax;
     163            rheadp->longi_t0   = evt_end->longi_t0;
     164            rheadp->longi_tmax = evt_end->longi_tmax;
     165            rheadp->longi_a    = evt_end->longi_a;
     166            rheadp->longi_b    = evt_end->longi_b;
     167            rheadp->longi_c    = evt_end->longi_c;
     168            rheadp->longi_chi2 = evt_end->longi_chi2;
     169            break;
     170          }
     171
     172        for (ph=0; ph<PH_IN_DATABLOCK; ph++)
     173          {
     174            if (Photons[ph].w <= 0.) break;
    128175
    129176            CPhotons[cphs].w = Photons[ph].w;
    130177            Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
    131             Photons[ph].x -= cheadp->CorePos[0][0];
    132             Photons[ph].y -= cheadp->CorePos[1][0];
     178
     179            /* TEMPORARY FIX, AM Nov 2001: we found that sometimes the value
     180               stored in Photons[ph].w is not correct, and can result in a
     181               wavelength beyond 600 nm, which makes the program crash later.
     182               Now we force wlen to its expected range:
     183               */
     184
     185            wlen = MIN(MAX(290.,wlen),600.);
     186
     187
     188            /* ADDED AM May 2002: now we take into account the telescope
     189             * position chosen in the Corsika input card via the CERTEL
     190             * option, which must be supplied also in the reflector input
     191             * card with the option "telescope_position x y". Otherwise
     192             * x = y = 0 is assumed, which is the standard mode of
     193             * generation for MAGIC.
     194             */
     195
     196            Photons[ph].x -= cheadp->CorePos[0][0]+Telescope_x;
     197            Photons[ph].y -= cheadp->CorePos[1][0]+Telescope_y;
    133198
    134199            /*  Increment number of photons     */
    … …  
    137202            /*  Calculate some quantities  */
    138203            theta = (float) acos(sqrt(
    139                     1.f - Photons[ph].u*Photons[ph].u
    140                         - Photons[ph].v*Photons[ph].v));
     204                                      MAX(0., 1.f - Photons[ph].u*Photons[ph].u
     205                                      - Photons[ph].v*Photons[ph].v)));
    141206
    142207            /*  Check absorption  */
     208
     209
    143210            if (absorption(wlen, Photons[ph].h, theta))
    144                 absphs++;
     211              absphs++;
    145212
    146213            /*  Check reflection  */
    147214            else if (0 != (ref_type =
    148215                           ph2cph(&Photons[ph], &CPhotons[cphs])))
    149                 refphs[ref_type-1]++;
    150             else    /*  Photon passed   */
    151             {
    152             Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph,
    153                    Photons[ph].x,Photons[ph].y,
    154                    Photons[ph].u,Photons[ph].v,theta);
    155             Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x,
    156                    CPhotons[cphs].y);
     216              refphs[ref_type-1]++;
     217            else                /*      Photon passed   */
     218              {
     219                Debug("Ph %d\t%f\t%f\t%f\t%f\t%f\n",ph,
     220                      Photons[ph].x,Photons[ph].y,
     221                      Photons[ph].u,Photons[ph].v,theta);
     222                Debug("CPh %d\t%d\t%f\t%f\n\n",cphs,ph,CPhotons[cphs].x,
     223                      CPhotons[cphs].y);
    157224
    158225
    … …  
    162229
    163230                /*  Update cphs */
    164                 if (++cphs == NR_OF_CPHOTONS)   /*  Overflow  */
    165                 {
     231                if (++cphs == NR_OF_CPHOTONS) /*  Overflow  */
     232                  {
    166233                    /*  if it is the first o/f open a tempfile  */
    167234                    if (overflow++ == 0)
    168                     {   Log("Spooling... ");
     235                      { Log("Spooling... ");
    169236                        tmpf = tmpfile();
    170237                        if (tmpf == NULL) FatalError(TEMP_ERROR_FTL);  }
    … …  
    176243                    cphs = 0;
    177244
    178                 }  /*  if overflow  */
    179             }   /*  else (=Photon passed)  */
    180         }  /*  end of loop inside datablock    */
    181     }   /*  end of loop on datablocks           */
     245                  } /*  if overflow  */
     246              } /*  else (=Photon passed)  */
     247          } /*  end of loop inside datablock    */
     248      } /*  end of loop on datablocks           */
    182249
    183250    /*  Check if there was an error while reading cerfile  */
    184251    if (read != 1) fseek(rflfile, writep, SEEK_SET);
    185252
    186     else    /*  no error: write the new event   */
    187 
    188     {   /*  Write "start of event" flag  */
     253    else                        /*  no error: write the new event       */
     254
     255      {                         /*  Write "start of event" flag  */
    189256        fwrite(FLAG_START_OF_EVENT, SIZE_OF_FLAGS, 1, rflfile);
    190257
    … …  
    202269        rheadp->CPhotons     = (long) overflow * NR_OF_CPHOTONS + cphs;
    203270        fwrite(rheadp, sizeof(RflHeader), 1, rflfile);
    204         for (myloop=0; myloop<sizeof(RflHeader)/4; myloop++)
    205             fprintf(chkf, "%e ", *((float *)rheadp+myloop));
    206         fputc('\n', chkf);
     271
     272/*        for (myloop=0; myloop<sizeof(RflHeader)/4; myloop++)
     273 *        fprintf(chkf, "%e ", *((float *)rheadp+myloop));
     274 *        fputc('\n', chkf);
     275 */
    207276
    208277        /*  If there was an overflow, append data from tempfile   */
    209278        if (overflow)
    210         {
     279          {
    211280            /*  Unload data from CPhotons  */
    212281            fwrite(CPhotons, sizeof(cphoton), cphs, tmpf);
    213282
    214283            /*  Transfer data  */
    215             fseek(tmpf, 0L, SEEK_SET);  /*  Start from the beginning  */
     284            fseek(tmpf, 0L, SEEK_SET); /*  Start from the beginning  */
    216285            while (overflow--)
    217             {   fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
     286              { fread (CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, tmpf);
    218287                fwrite(CPhotons, sizeof(cphoton), NR_OF_CPHOTONS, rflfile); }
    219288
    … …  
    236305
    237306    return read == 1;
    238 }  /*  end of ProcessEvent  */
     307  } /*  end of ProcessEvent  */
    239308
    240309static int GetEvent(void)
    241310{   int found = FALSE,          /*  event found         */
    242         isWrong = FALSE;        /*  cerfile is wrong    */
    243     static int newFile = TRUE; /*  if TRUE, check if cerfile is valid  */
    244 
    245     do
    246     {   /* In case the just-opened file is a valid cerfile,
    247             starting with a RUNH, loop until a valid event is found:
    248             id est with: first_Event <= EvtNumber <= last_Event
    249             and             low_Ecut <=  Etotal   <= high_Ecut
    250             If the search was successful, "found" is set to TRUE.
    251             If there are reading errors, "isWrong" is set to TRUE. */
     311      isWrong = FALSE;          /*  cerfile is wrong    */
     312      static int newFile = TRUE; /*  if TRUE, check if cerfile is valid */
     313
     314      do
     315        { /* In case the just-opened file is a valid cerfile,
     316             starting with a RUNH, loop until a valid event is found:
     317             id est with: first_Event <= EvtNumber <= last_Event
     318             and             low_Ecut <=  Etotal   <= high_Ecut
     319             If the search was successful, "found" is set to TRUE.
     320             If there are reading errors, "isWrong" is set to TRUE. */
    252321
    253322        if (newFile
    … …  
    308377
    309378FILE *GetNextFile(char *cername)
    310 {   FILE *f = NULL;             /*  return value (cerfile ptr)  */
     379{   FILE *inputfile = NULL;             /*  return value (cerfile ptr)  */
    311380    char *value_ptr;            /*  ptr at parm value   */
    312381    extern char line[];         /*  white chars (init)  */
    … …  
    324393
    325394        /*  If you found a line with some meaning, try to open the file  */
    326         if ((f=fopen(value_ptr, "r")) == NULL)
     395        if ((inputfile=fopen(value_ptr, "r")) == NULL)
    327396            Message(CERF_NFND__MSG, value_ptr);
    328397
    … …  
    346415                last_Event  = 1000000;  }}}
    347416
    348     while (f == NULL);  /*  Loop until a readable file is found  */
     417    while (inputfile == NULL);  /*  Loop until a readable file is found  */
    349418
    350419    /*  If a new cerfile is met, write "START_OF_RUN"  */
    351     if (f) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
    352 
    353     return f;
     420    if (inputfile) fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
     421
     422    return inputfile;
    354423}   /*  end of GetNextFile  */
    355424
    … …  
    378447    {   free(CPhotons);         Log(VECT_FREE__LOG, "CPhotons");  }
    379448}   /*  end of clean  */
     449
Note: See TracChangeset for help on using the changeset viewer.