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

Timestamp:

11/14/02 21:39:04 (22 years ago)

Author:

bigongia

Message:

*** empty log message ***

File:

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

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

@@ -120 +120 @@
     extern void makeOmegaI(float theta, float phi);
 
+    extern int wobble_position;
 
     /*  Various counters: phs = absphs + refphs[0..3] + cphs  */
@@ -127 +128 @@
     cphs;                       /*  Number of cphotons          */
 
+    /* AM, Nov 2002 */
+    /* Counters of the number of cphs produced by different particle types: */
+    long 
+    mu_cphs,                    /* by muons  */
+    el_cphs,                    /* by electrons  */
+    other_cphs;                 /* by other particles  */
+
+
+    int parent_id;              /* Id of the particle producing the C-photon */
+
     FILE *tmpf=NULL;            /*  Temp fp to handle o/f       */
     size_t read;                /*  items read: != 1 on error   */
@@ -140 +151 @@
     theta;                      /*  photon zenith angle */
 
+    float tel_theta, tel_phi;   /* Indicate telescope orientation, following
+                                 * Corsika's criterium: theta is the zenith
+                                 * angle, but phi is the azimuth of the vector
+                                 * along the telescope axis going from the
+                                 * camera to the dish, measured from north,
+                                 * anticlockwise (0-2pi).
+                                 */
+
     /*  Reset counters and set fileptr  */
     phs= absphs= refphs[0]= refphs[1]= refphs[2]= refphs[3]= cphs= 0;
@@ -147 +166 @@
     TranslateHeader(rheadp, cheadp);
 
-    /*  Calculate OmegaCT matrices  */
     if (is_Fixed_Target)
-      { makeOmega(fixed_Theta, fixed_Phi);
-        makeOmegaI(fixed_Theta, fixed_Phi);   }
-    else
-      { makeOmega(cheadp->Theta, cheadp->Phi);
-        makeOmegaI(cheadp->Theta, cheadp->Phi);   }
+      { 
+        tel_theta = fixed_Theta;
+        tel_phi   = fixed_Phi;
+      }
+    else  /* Use the direction of the primary if no telescope 
+           * orientation is given: */
+      { 
+        tel_theta = cheadp->Theta;
+        tel_phi   = cheadp->Phi;
+      }
+
+    /* AM, Nov 2002: added wobble mode option. It refers 
+     * to wobble position +-1, 3 defined in TDAS  01-05 */
+
+    switch (wobble_position)
+      {
+      case 1:
+        tel_theta = acos(cos(tel_theta)/1.000024369);
+        tel_phi  -= atan(0.006981317/sin(tel_theta));
+        break;
+      case -1:
+        tel_theta = acos(cos(tel_theta)/1.000024369);
+        tel_phi  += atan(0.006981317/sin(tel_theta));
+        break;
+      case 3:
+        tel_theta += 0.006981317;  /* 0.4 degrees */
+        break;
+      default:
+        break;
+      } 
+
+    /* Calculate OmegaCT matrices  */
+    /* AM 11/2002 : added PI to 2nd argument in makeOmega, to take 
+     * into account that theta and phi from Corsika indicate the
+     * direction of the momentum (hence, downgoing) of the primary 
+     * particle, phi measured from positive x axis, anticlockwise, and
+     * theta measured from negative z axis (see Corsika manual). In the
+     * call to makeOmega, the second argument must be the azimuth of the 
+     * telescope up-going pointing direction.
+     */
+
+    makeOmega(tel_theta, tel_phi+M_PI);
+    makeOmegaI(tel_theta, tel_phi+M_PI);
+
     memcpy( OmegaCT, Omega, 9*sizeof(float) );
     memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
@@ -175 +232 @@
     rheadp->CorePos[1][0] += Telescope_y;
 
+    /* Initialize cphoton counters: */
+    el_cphs = mu_cphs = other_cphs = 0;
+
     /*  Loop on data blocks  */
 
@@ -214 +274 @@
 
             CPhotons[cphs].w = Photons[ph].w;
+
+            /* AM Nov 2002: now we read the type of particle which produced
+             * the Cherenkov photon :
+             */
+            parent_id = (int)Photons[ph].w/100000;
             Photons[ph].w = wlen = (float) fmod(Photons[ph].w, 1000.);
-
+ 
             /* AM Nov 2001: we found that sometimes the value
                stored in Photons[ph].w is not correct, and can result in a 
@@ -244 +309 @@
             /*  Increment number of photons     */
             phs++;
-            
+
             /*  Calculate some quantities  */
             theta = (float) acos(sqrt(
@@ -269 +334 @@
                       CPhotons[cphs].y);
 
+                /* AM Nov 2002: We now count for every event how many
+                 * photons were produced by each type.
+                 */
+
+                switch (parent_id)
+                  {
+                  case 2:
+                    el_cphs ++;
+                    break;
+                  case 3:
+                    el_cphs ++;
+                    break;
+                  case 5:
+                    mu_cphs ++;
+                    break;
+                  case 6:
+                    mu_cphs ++;
+                    break;
+                  default:
+                    other_cphs++;
+                  }
 
                 /*  Update first/last arrival times  */
@@ -314 +400 @@
         rheadp->OutOfChamPhs = refphs[3];
         rheadp->CPhotons     = (long) overflow * NR_OF_CPHOTONS + cphs;
+
+        if (rheadp->CPhotons > 0.)
+          {
+            /* Fill in relative fraction of each particle type producing 
+             * Cherenkov light in this even:
+             */
+            rheadp->elec_cph_fraction = (float)el_cphs/(float)rheadp->CPhotons;
+            rheadp->muon_cph_fraction = (float)mu_cphs/(float)rheadp->CPhotons;
+            rheadp->other_cph_fraction = (float)other_cphs/(float)rheadp->CPhotons;
+          }
+        else
+          {
+            rheadp->elec_cph_fraction = 0.;
+            rheadp->muon_cph_fraction = 0.;
+            rheadp->other_cph_fraction = 0.;
+          }
+
         fwrite(rheadp, sizeof(RflEventHeader), 1, rflfile);
 
@@ -363 +466 @@
   static int newFile = TRUE;  /*  if TRUE, check if cerfile is valid    */
 
-  CerRunHeader cRunHead;
-  
+  RflRunHeader RflRunHead;
+
+  extern int atmModel;      /* current atm. model */
+
   do
     {
@@ -384 +489 @@
             {
               /* Write Reflector "run header" (one per cer file!): */
-              memcpy(&cRunHead, cheadp, sizeof(CerEventHeader));
-              fwrite(&cRunHead, sizeof(CerRunHeader), 1, rflfile);
+              memcpy(&RflRunHead, cheadp, sizeof(RflRunHead));
+              RflRunHead.wobble_mode = wobble_position;
+              RflRunHead.atmospheric_model = atmModel;
+              fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile);
             }
 
@@ -503 +610 @@
     {   free(AxisDeviation[1]); Log(VECT_FREE__LOG, "AxisDeviation[1]");  }
 
-    if (ct_Focal)
-    {   free(ct_Focal);         Log(VECT_FREE__LOG, "ct_Focal");  }
-
     if (CPhotons)
     {   free(CPhotons);         Log(VECT_FREE__LOG, "CPhotons");  }