Changeset 1614


Ignore:
Timestamp:
11/14/02 21:39:04 (22 years ago)
Author:
bigongia
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Simulation/Detector/ReflectorII
Files:
14 added
3 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile

    r1535 r1614  
    1919#
    2020# $RCSfile: Makefile,v $
    21 # $Revision: 1.4 $
     21# $Revision: 1.5 $
    2222# $Author: bigongia $
    23 # $Date: 2002-10-09 18:15:27 $
     23# $Date: 2002-11-14 21:39:01 $
    2424#
    2525##################################################################

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

    r1431 r1614  
    99#define RandomNumber  ranf()
    1010
    11 /*  Local declarations */
    12 static int atmModel=0;          /*  current atm. model  */
     11/* AM Nov 2002: atmModel is now an external variable: */
     12int atmModel=0;   /* current atm. model */
    1313
    1414/*  Function declarations  */

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/config.mk.linux

    r1535 r1614  
    1818#---------------------------------------------------------------
    1919# $RCSfile: config.mk.linux,v $
    20 # $Revision: 1.2 $
     20# $Revision: 1.3 $
    2121# $Author: bigongia $
    22 # $Date: 2002-10-09 18:15:27 $
     22# $Date: 2002-11-14 21:39:01 $
    2323##################################################################
    2424# @maintitle

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

    r1535 r1614  
    2424float   *Reflectivity[2];       /*  reflectivity table          */
    2525float   *AxisDeviation[2];      /*  axis deviation table        */
    26 float   *ct_Focal;              /*  focal length table          */
    2726
    2827/*  Prototypes  */
    … …  
    3029static void ReadReflectivity(char *datname);
    3130static void ReadAxisDev(char *datname);
    32 static void ReadFocals(void);
    3331
    3432static void ReadMirrorTable(FILE *geofile)
    … …  
    5048                         &ct_data[i].xn,            &ct_data[i].yn,     &ct_data[i].zn))
    5149          break;
     50
    5251      }
    5352    if (i < ct_NMirrors)
    … …  
    121120}   /*  end of ReadAxisDev  */
    122121
    123 static void ReadFocals(void)
    124 {   int loop;
    125 
    126     if ((ct_Focal = (float *) malloc(sizeof(float) * ct_NMirrors)) == NULL)
    127         FatalError(FOCL_FEW___FTL, ct_NMirrors);
    128 
    129     for (loop=0; loop<ct_NMirrors; loop++)
    130         ct_Focal[loop] = ct_data[loop].f;
    131 }   /*  end of ReadFocals  */
    132 
    133122void GeometrySwitch(FILE *geofile)
    134123{   char *value_ptr = NULL;     /*  ptr at parm value   */
    … …  
    218207
    219208
    220     ReadFocals();
    221209}   /*  end of GeometrySwitch  */
    222210

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/header.h

    r1535 r1614  
    22#define __RFL_HEADER__
    33
    4 /*  RUNH from cerfile. See CORSIKA manual for explanations */
     4/*  RUNH for reflector file. Identical to Corsika's Run Header (see CORSIKA
     5 *  manual for explanations), plus some other info like the wobble mode
     6 *  which has been used. We have put this in the unused variables.
     7 */
    58typedef struct
    69{   char RUNH[4];
    … …  
    2124    /* Physical constants and interaction flags (see CORSIKA manual): */
    2225    float C[50];
    23     float dummy1[20]; /* not used */
     26    float wobble_mode;      /* Indicates wobble mode with which
     27                               reflector has been run */
     28    float atmospheric_model; /* Indicates atmospheric model used in
     29                                absorption simulation. 0 = no atmosphere,
     30                                1 = atm_90percent, 2 = atm_isothermal,
     31                                3 = atm_corsika.
     32                                */
     33    float dummy1[18]; /* not used */
    2434    float CKA[40];
    2535    float CETA[5];
    … …  
    3040    float CATM[5];
    3141    float NFL[4];
    32 }   CerRunHeader;
    33 
     42}   RflRunHeader;
    3443
    3544/*  EVTH from cerfile. See CORSIKA manual for explanations */
    … …  
    239248    float       OutOfChamPhs;   /*  Photons outside the camera          */
    240249    float       CPhotons;       /*  Photons reaching the camera         */
    241 
    242     float       dummy[10];      /* not used */
     250   
     251    /* Now follow the fraction of photons reaching the camera produced by
     252     * electrons, muons and other particles respectively:
     253     */
     254    float       elec_cph_fraction;
     255    float       muon_cph_fraction;
     256    float       other_cph_fraction;
     257
     258    float       dummy[7];      /* not used */
     259
    243260}   RflEventHeader;
    244261

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

    r1431 r1614  
    2323float Telescope_x = 0.;
    2424float Telescope_y = 0.;         /* Telescope coordinates (cm) */
     25
     26int wobble_position = 0;
    2527
    2628int   is_Random_Pointing=FALSE; /*  random pointing?    */

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h

    r1535 r1614  
    4343extern float Telescope_x, Telescope_y;  /* Telescope coordinates (cm) */
    4444
     45extern int   wobble_position;
     46
    4547extern int   is_Random_Pointing;        /*  random pointing?    */
    4648extern float Random_Pointing_MaxDist;   /*  in metres           */
    … …  
    7274extern float  *Reflectivity[];          /*  ptr to refl. table          */
    7375extern float  *AxisDeviation[];         /*  ptr to axisdev. table       */
    74 extern float  *ct_Focal;                /*  ptr to focal data           */
    7576
    7677/*  Photon structures  */

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/input.card

    r1432 r1614  
    1 reflector 0.5
     1reflector 0.6
    22#
    33# Sample parameters file
    … …  
    1010telescope_position 0. 0.
    1111#
     12wobble_mode 0
    1213max_events   1000000000
    1314#energy_cuts   800. 1000.
    1415#
    15 ct_file        ../Data/magic.def
     16ct_file        mypath/Data/magic.def
    1617#
    17 #reflectivity_file /mypath/Data/reflectivity.dat
     18reflectivity_file /mypath/Data/reflectivity.dat
    1819#
    19 #axisdev_file /mypath/Data/axisdev.dat
     20axisdev_file /mypath/Data/axisdev.dat
    2021#
    2122output_file   /myout/my.rfl

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

    r1535 r1614  
    128128                Message(MAX__EVTS__MSG, max_Events=atol(value_ptr));
    129129                break;
     130
     131                /* AM, Nov 2002, added wobble option: */
     132              case wobble_mode:
     133                wobble_position = (int) atoi(value_ptr);
     134                if ( (wobble_position > 1 || wobble_position < -1) &&
     135                    wobble_position != 3)
     136                  {
     137                    Message(WOBBLE_POS_ERR);
     138                    exit(-1);
     139                  }
     140                else
     141                  {
     142                    Message(WOBBLE_POS_MSG, wobble_position);
     143                    break;
     144                  }
    130145              case energy_cuts:
    131146                low_Ecut = (float) atof(value_ptr);

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.h

    r1431 r1614  
    1313T(telescope_position),/* position towards which CT is pointing */    \
    1414T(max_events),      /* maximum number of event to read */          \
     15T(wobble_mode),      /* wobble mode option */          \
    1516T(energy_cuts),     /* lowest/highest energy allowed */            \
    1617T(seeds),           /* seeds for random number generation */       \
    … …  
    4748#define MAX__EVTS__MSG          /*  max. nr. of evts.   */ \
    4849    "Processing at most %ld events.\n"
     50
     51#define WOBBLE_POS_ERR "wobble mode can only be  -1, 0, 1 or 3\n"
     52#define WOBBLE_POS_MSG "Using wobble position %d\n"
     53
    4954#define ENRG_LIMIT_ERR          /*  no parms            */ \
    5055    "Error while parsing \"energy_cuts\" fields.\n" \
    … …  
    6469    "No specified input file could be opened.\n" \
    6570    " *** Please specify some *valid* filename after the cer_files directive.\n"
    66 
    67 
    6871extern char axisdev_filename[256], reflectivity_filename[256];
    6972

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

    r1535 r1614  
    118118  r[0] = ph->u;
    119119  r[1] = ph->v;
    120   r[2] = (float) sqrt(1.0 - r[0]*r[0] - r[1]*r[1]);
     120
     121// AM 11/2002: fixed line below: u v are the direction cosines of the
     122// *downgoing* photon. Hence, third component must be negative!
     123// This was a serious bug affecting all versions before 0.6 (see TDAS
     124// note on Reflector program 0.6).
     125
     126  r[2] = (float) -sqrt(1.0 - r[0]*r[0] - r[1]*r[1]);
    121127
    122128  /* get photon time and production height */   
    … …  
    494500  /*
    495501   if we still have the photon, we continue with the reflexion;       
    496    we calculate normal vector in this point
    497    (and normalize, with the sign changed)
     502   we calculate normal vector in this point and normalize:
    498503  */
    499504
    500505  rnor[0] = 2.0f*xcut[0];
    501506  rnor[1] = 2.0f*xcut[1];
    502   rnor[2] = (float) (2.0*(xcut[2] - 2.0*ct_Focal[i_mirror]));
     507  rnor[2] = (float) (2.0*(xcut[2] - 2.0*ct_data[i_mirror].f));
    503508
    504509  /* CBC */
    … …  
    506511  /* CBC */
    507512
    508   rnorm    = -NORM( rnor );
     513// Changed AM, 11/2002: now we use the normal vector going "outwards"
     514// from inside the sphere (=removed minus sign in normalization below).
     515// It is easier to do so, since now the vector rm indicating the
     516// photon direction also goes from the front to the back of the mirror.
     517
     518  rnorm    = NORM( rnor );
    509519  rnor[0] /= rnorm;
    510520  rnor[1] /= rnorm;
    … …  
    525535  */
    526536
    527   calpha = (float) fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
     537// AM 11/2002: removed absolute value in scalar
     538// product below (it is now unnecessary):
     539
     540  calpha = (float) (rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
    528541
    529542  /* CBC */
    … …  
    532545       
    533546  /*  finally!!! we have the reflected trajectory of the photon */
     547
    534548
    535549  rrefl[0] = (float) (2.0*rnor[0]*calpha - rm[0]);
    … …  
    633647   calculate angle of incidence between tray. and camera plane
    634648   the camera plane is
    635    0 y + 0 y + z - ct_Focal = 0 => (A,B,C,D) = (0,0,1,-ct_Focal)
     649   0 x + 0 y + z - ct_Focal_mean = 0 => (A,B,C,D) = (0,0,1,-ct_Focal_mean)
    636650   from Table 3.20 "Tasch. der Math."
    637651  */
    … …  
    679693  /* Output */
    680694
    681 /*  cph->w   = wl;  */
    682   cph->x   = xcam[0];
    683   cph->y   = xcam[1];
     695  /* AM Nov 2002: Added one further change of coordinates so that the camera
     696   * images have the "right" orientation: they will appear as seen by an
     697   * observer on ground, standing behind the mirror dish and looking towards
     698   * the camera. Formerly cph->x and cph->y were simply xcam[0] and xcam[1].
     699   */
     700
     701  cph->x   = -xcam[1];
     702  cph->y   = -xcam[0];
     703
    684704  cph->u   = r[0];
    685705  cph->v   = r[1];

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

    r1535 r1614  
    120120    extern void makeOmegaI(float theta, float phi);
    121121
     122    extern int wobble_position;
    122123
    123124    /*  Various counters: phs = absphs + refphs[0..3] + cphs  */
    … …  
    127128    cphs;                       /*  Number of cphotons          */
    128129
     130    /* AM, Nov 2002 */
     131    /* Counters of the number of cphs produced by different particle types: */
     132    long
     133    mu_cphs,                    /* by muons  */
     134    el_cphs,                    /* by electrons  */
     135    other_cphs;                 /* by other particles  */
     136
     137
     138    int parent_id;              /* Id of the particle producing the C-photon */
     139
    129140    FILE *tmpf=NULL;            /*  Temp fp to handle o/f       */
    130141    size_t read;                /*  items read: != 1 on error   */
    … …  
    140151    theta;                      /*  photon zenith angle */
    141152
     153    float tel_theta, tel_phi;   /* Indicate telescope orientation, following
     154                                 * Corsika's criterium: theta is the zenith
     155                                 * angle, but phi is the azimuth of the vector
     156                                 * along the telescope axis going from the
     157                                 * camera to the dish, measured from north,
     158                                 * anticlockwise (0-2pi).
     159                                 */
     160
    142161    /*  Reset counters and set fileptr  */
    143162    phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0;
    … …  
    147166    TranslateHeader(rheadp, cheadp);
    148167
    149     /*  Calculate OmegaCT matrices  */
    150168    if (is_Fixed_Target)
    151       { makeOmega(fixed_Theta, fixed_Phi);
    152         makeOmegaI(fixed_Theta, fixed_Phi);   }
    153     else
    154       { makeOmega(cheadp->Theta, cheadp->Phi);
    155         makeOmegaI(cheadp->Theta, cheadp->Phi);   }
     169      {
     170        tel_theta = fixed_Theta;
     171        tel_phi   = fixed_Phi;
     172      }
     173    else  /* Use the direction of the primary if no telescope
     174           * orientation is given: */
     175      {
     176        tel_theta = cheadp->Theta;
     177        tel_phi   = cheadp->Phi;
     178      }
     179
     180    /* AM, Nov 2002: added wobble mode option. It refers
     181     * to wobble position +-1, 3 defined in TDAS  01-05 */
     182
     183    switch (wobble_position)
     184      {
     185      case 1:
     186        tel_theta = acos(cos(tel_theta)/1.000024369);
     187        tel_phi  -= atan(0.006981317/sin(tel_theta));
     188        break;
     189      case -1:
     190        tel_theta = acos(cos(tel_theta)/1.000024369);
     191        tel_phi  += atan(0.006981317/sin(tel_theta));
     192        break;
     193      case 3:
     194        tel_theta += 0.006981317;  /* 0.4 degrees */
     195        break;
     196      default:
     197        break;
     198      }
     199
     200    /* Calculate OmegaCT matrices  */
     201    /* AM 11/2002 : added PI to 2nd argument in makeOmega, to take
     202     * into account that theta and phi from Corsika indicate the
     203     * direction of the momentum (hence, downgoing) of the primary
     204     * particle, phi measured from positive x axis, anticlockwise, and
     205     * theta measured from negative z axis (see Corsika manual). In the
     206     * call to makeOmega, the second argument must be the azimuth of the
     207     * telescope up-going pointing direction.
     208     */
     209
     210    makeOmega(tel_theta, tel_phi+M_PI);
     211    makeOmegaI(tel_theta, tel_phi+M_PI);
     212
    156213    memcpy( OmegaCT, Omega, 9*sizeof(float) );
    157214    memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
    … …  
    175232    rheadp->CorePos[1][0] += Telescope_y;
    176233
     234    /* Initialize cphoton counters: */
     235    el_cphs = mu_cphs = other_cphs = 0;
     236
    177237    /*  Loop on data blocks  */
    178238
    … …  
    214274
    215275            CPhotons[cphs].w = Photons[ph].w;
     276
     277            /* AM Nov 2002: now we read the type of particle which produced
     278             * the Cherenkov photon :
     279             */
     280            parent_id = (int)Photons[ph].w/100000;
    216281            Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
    217 
     282 
    218283            /* AM Nov 2001: we found that sometimes the value
    219284               stored in Photons[ph].w is not correct, and can result in a
    … …  
    244309            /*  Increment number of photons     */
    245310            phs++;
    246            
     311
    247312            /*  Calculate some quantities  */
    248313            theta = (float) acos(sqrt(
    … …  
    269334                      CPhotons[cphs].y);
    270335
     336                /* AM Nov 2002: We now count for every event how many
     337                 * photons were produced by each type.
     338                 */
     339
     340                switch (parent_id)
     341                  {
     342                  case 2:
     343                    el_cphs ++;
     344                    break;
     345                  case 3:
     346                    el_cphs ++;
     347                    break;
     348                  case 5:
     349                    mu_cphs ++;
     350                    break;
     351                  case 6:
     352                    mu_cphs ++;
     353                    break;
     354                  default:
     355                    other_cphs++;
     356                  }
    271357
    272358                /*  Update first/last arrival times  */
    … …  
    314400        rheadp->OutOfChamPhs = refphs[3];
    315401        rheadp->CPhotons     = (long) overflow * NR_OF_CPHOTONS + cphs;
     402
     403        if (rheadp->CPhotons > 0.)
     404          {
     405            /* Fill in relative fraction of each particle type producing
     406             * Cherenkov light in this even:
     407             */
     408            rheadp->elec_cph_fraction = (float)el_cphs/(float)rheadp->CPhotons;
     409            rheadp->muon_cph_fraction = (float)mu_cphs/(float)rheadp->CPhotons;
     410            rheadp->other_cph_fraction = (float)other_cphs/(float)rheadp->CPhotons;
     411          }
     412        else
     413          {
     414            rheadp->elec_cph_fraction = 0.;
     415            rheadp->muon_cph_fraction = 0.;
     416            rheadp->other_cph_fraction = 0.;
     417          }
     418
    316419        fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile);
    317420
    … …  
    363466  static int newFile = TRUE;  /*  if TRUE, check if cerfile is valid    */
    364467
    365   CerRunHeader cRunHead;
    366  
     468  RflRunHeader RflRunHead;
     469
     470  extern int atmModel;      /* current atm. model */
     471
    367472  do
    368473    {
    … …  
    384489            {
    385490              /* Write Reflector "run header" (one per cer file!): */
    386               memcpy(&cRunHead, cheadp, sizeof(CerEventHeader));
    387               fwrite(&cRunHead, sizeof(CerRunHeader), 1, rflfile);
     491              memcpy(&RflRunHead, cheadp, sizeof(RflRunHead));
     492              RflRunHead.wobble_mode = wobble_position;
     493              RflRunHead.atmospheric_model = atmModel;
     494              fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile);
    388495            }
    389496
    … …  
    503610    {   free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]");  }
    504611
    505     if (ct_Focal)
    506     {   free(ct_Focal);         Log(VECT_FREE__LOG, "ct_Focal");  }
    507 
    508612    if (CPhotons)
    509613    {   free(CPhotons);         Log(VECT_FREE__LOG, "CPhotons");  }
Note: See TracChangeset for help on using the changeset viewer.