mEmcDefGeom.c

Go to the documentation of this file.

#include "mEmcDefGeom.h"
#include "emlLib.h"

#include "gsl/gsl_math.h"

#include <math.h>

long type_of_call mEmcDefGeom_(
  TABLE_HEAD_ST  *dEmcGeaParams_h,  DEMCGEAPARAMS_ST    *dEmcGeaParams ,
  TABLE_HEAD_ST   *dEmcGeometry_h,   DEMCGEOMETRY_ST     *dEmcGeometry )
{
/*:>--------------------------------------------------------------------
**: ROUTINE:    mEmcDefGeom_
**: DESCRIPTION: Physics Analysis Module ANSI C template.
**:             This is an ANSI C Physics Analysis Module template
**:             automatically generated by stic from mEmcDefGeom.idl.
**:             Please edit comments and code.
**: ARGUMENTS:
**:       IN:
**:    INOUT:
**:      OUT:
**:      dEmcGeaParams    - PLEASE FILL IN DESCRIPTION HERE
**:     dEmcGeaParams_h   - header Structure for dEmcGeaParams
**:       dEmcGeometry    - PLEASE FILL IN DESCRIPTION HERE
**:      dEmcGeometry_h   - header Structure for dEmcGeometry
**: RETURNS:    STAF Condition Value
**:>------------------------------------------------------------------*/
  /*              */

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        subroutine EMC_GINI_NEW()
        
        Authour :
                 This routine has been converted to C on Mon Mar  2  1998 by:
                 
                 Phool Chand            phool@phenix.barc.ernet.in
                 Dipanwita Dutta        dipa@phenix.barc.ernet.in

        Modifications:
                 Gabor David, May 16 1998
        
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/*     Major revision: November 19, 1998, G. David

       All through the chain we adopt much of the PHENIX numbering scheme
       West arm is arm=0, East arm is arm=1
       West arm consists of former sectors 1-4
       East arm consists of former sectors 5-8, BUT IN REVERSE ORDER
       Former    Current
       5         3  (PbSc)
       6         2  (PbSc)
       7         1  (PbGl)
       8         0  (PbGl)

       Numbering of towers within a sector:
       Columns (z direction) are numbered 0-71 (0-95 for PbGl)
       Rows (y, or phi direction) are numbered 0-35 (0-47 for PbGl)

       HOWEVER:
       In the West arm column 0 is the column at the highest (positive) z,

       whereas
       in the East arm column 0 is the column at the lowest (negative) z
       (as before).

       This is done in order to be (more) conform with hardware numbering,
       where ALL sectors have their 0,0 element at the lower left corner
       WHEN VIEWED FROM BEHIND

       Note, that this new numbering is different from the way geometry
       is defined is PISA, where sectors are numbered 1-8, continuously
       in phi, and tower numbers start at the lowest negative z.
       This change effects mEmcGeaMakeRaw, too.

*/

/********************************************************/

#define max1(x,y)       ( ( int ) ( (x) > (y)  ) ? (x) : (y) )
#define amax1(x,y)      ( ( (x) > (y)  ) ? (x) : (y) )
#define max(x,y)        ( ( (x) > (y)  ) ? (x) : (y) )
#define min(x,y)        ( ( (x) < (y)  ) ? (x) : (y) )

   int size;

   int pcons1, pcons2;

/****************************************************************/

#define max_chanz       96
#define max_chany       48
#define max_subdet      8
#define max_fields      10


   float        ra_det[120][max_subdet];

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        c       The number of walls (and detector types) defined can vary;
        c       the array I_TRANS translates the "wall,subdetector" pair
        c       where wall can be 1 to 8, subdetector is 1 for PbSc,
        c       2 for PbGl
        c       to a single number (1 to 8), referenced in emc_geom and
        c       the hit arrays
        c       
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   int  i_trans[8][2]   =
                {
                1, 0,
                2, 0,
                3, 0,
                4, 0,
                5, 0,
                6, 0,
                0, 7,
                0, 8
                };

   int  i_limits[16]    = { 72,36,72,36,72,36,72,36,72,36,72,36, 96,48,96,48 };

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        c
        c       fields in emc_geom:
        c       1,z,y,det --> x position of the center of the front face
        c       2,z,y,det --> y position of the center of the front face
        c       3,z,y,det --> z position of the center of the front face
        c       4,z,y,det --> theta at the center of the front face
        c       5,z,y,det --> phi at the center of the front face
        c       6,z,y,det --> distance (r) at the center of the front face
        c       7,z,y,det --> e(x) -> unit vector for the axis of the tower
        c       8,z,y,det --> e(y) -> unit vector for the axis of the tower
        c       9,z,y,det --> e(z) -> unit vector for the axis of the tower
        c       
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   float        emc_geom[max_fields][max_chanz][max_chany][max_subdet];
   int  i,j,iz,iy,i1;
   int ll;
   float          r_work1;
   float          r_appx,r_appy,r_appz;
   float r_measthe = 0.0;
   float r_measphi = 0.0;
   
   float   r_sectheta,r_secphi,r_dist,r_distz,r_disty;


/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        cgd***      EXTERNAL    GFDETU
        c
        c
        c     MAIN GEOMETRY parameter file (phnx.par) segment
        c       
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   int  emc_walls       =       8;       /*      number of sectors, lately (4 in each arm)      */

   float        emc_cutgam      =       0.001;
   float        emc_cutele      =       0.001;
   float        emc_cutneu      =       0.001;
   float        emc_cuthad      =       0.001;
   float        emc_cutmuo      =       0.001;
   float        ra_sampfrac[23][2]      =
        {
        .0001, .2443, .0002, .2432, .0005, .2313, .0008, .2204, 
        .0010, .2054, .0012, .1981, .0015, .1889, .0018, .1809,
        .0020, .1760, .0022, .1717, .0025, .1650, .0028, .1607, 
        .0030, .1565, .0040, .1429, .0050, .1307, .0060, .1214,
        .0070, .1130, .0080, .1058, .0090, .0993, .0100, .0936, 
        .0120, .0824, .0150, .0710, .0200, .0576
        };

   float        r_cutgam,r_sampfrac;
   int  l_automatic_threshold   =       TRUE;   
   int  ind1,ind2;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     New calibration: three cvolu_opt(5,8)  --  AUAU, PPLO, PPHI
     currently four GCUTS thresholds 0.001, 0.005, 0.02, 0.1    
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   /* +++++++++++++++++++++++++++++++++++++++++++++++++++

      New threshold setting 'PPME' introduced to accomodate
      high pt photons in Pythia events - G. David, Feb. 11, 2000

    +++++++++++++++++++++++++++++++++++++++++++++++++++++  */
   float        ra_newcalib[4][4]       =
                        {
                                0.2285, 0.1891, 0.1715, 0.1428,
                                0.2285, 0.1891, 0.1715, 0.1428,
                                0.1872, 0.1658, 0.1715, 0.1387,
                                0.1382, 0.1593, 0.1715, 0.1295
                };      


/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        c
        c
        c
        c
        c     End of Lead Scintillator (PbSc) geometry parameters
        c
        c
        c*******************************************************************
        c
        c     When calculating the impact position you have to add the
        c     depth if the insensitive material in front (but already
        c     inside the detector volume)
        c       
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   float        rpos_plus_pbsc  =       0.0;
   float        rpos_plus_pbgl  =       0.0;

   int  istaf;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        c
        c
        cgd***  DATA emc_iwall_max/8/,emc_itype_max/3/,emc_i1_max/8/
        c
        c       Indices within supermodule
        c
        cgd***  DATA emc_ind1_max_sc/18/,emc_ind2_max_sc/144/
        cgd***  DATA emc_ind1_max_gl/192/,emc_ind2_max_gl/24/
        c
        c       Misc
        c
        cgd***  DATA emc_itrack_max/32767/,emc_spart_max/62/,
        cgd***     1       emc_ncycle_max/31/
        c
        c
        c     THIS IS HOW PARAMETERS WERE WRITTEN OUT IN emc.f
        c
        c     Defining parameters for cells and supercells for each sector
        c     and calorimeter type (max 8 sectors, types: Shish-Kebab, PbGlass)
        c
        c
        c          udetpar(1) = float(emc_walls)     ! number of walls from PHNX.PAR
        c          udetpar(2) = float(emc_opt)       ! EMCal option from PHNX.PAR
        c          udetpar(3) = float(iwall)        ! wall number (1-4)
        c          udetpar(4) = float(itype)        ! detector type (Sh-K,PbGl)
        c          udetpar(5) = angle               ! phi angle of wall center
        c          udetpar(6) = rpos                ! radial position of wall center
        c          udetpar(7) = zc_start             ! center of first cell, z coor.
        c          udetpar(8) = yc_start             ! center of first cell, y coor.
        c          udetpar(9) = lsiz                ! long. size of a cell
        c          udetpar(10) = tsiz               ! transverse size of a cell
        c          udetpar(11) = FLOAT(no_modz)     ! No. of cells in z in a supermod.
        c          udetpar(12) = FLOAT(no_mody)     ! No. of cells in y in a supermod.
        c          udetpar(13) = FLOAT(no_smodz)    ! No. of supermods. in z / wall
        c          udetpar(14) = FLOAT(no_smody)    ! No. of supermods. in y / wall
        c          udetpar(15) = 0.0                ! 
        c          udetpar(16) = 0.0                !
        c          udetpar(17) = 0.0                !
        c          udetpar(18) = FLOAT(scint_emc_med)  ! Shish-Kebab scint. medium
        c          udetpar(19) = 0.0                !
        c          udetpar(20) = 0.0                !
        c          udetpar(21) = 0.0
        c          udetpar(22) = FLOAT(emc_debug)
        c          udetpar(23) = gcuts(1)
        c          udetpar(24) = gcuts(2)
        c          udetpar(25) = gcuts(3)
        c          udetpar(26) = gcuts(4)
        c          udetpar(27) = gcuts(5)
        c          udetpar(28) = 0.0     !
        c          udetpar(29) = 0.0     !
        c          udetpar(30) = emc_r_min_sc ! bitp lower limit, PbSc
        c          udetpar(31) = emc_r_max_sc ! bitp upper limit, PbSc
        c          udetpar(32) = emc_r_step ! bitp stepsize, PbSc
        c          udetpar(33) = emc_z_min ! bitp lower limit
        c          udetpar(34) = emc_z_max ! bitp upper limit
        c          udetpar(35) = emc_z_step ! bitp stepsize
        c          udetpar(36) = emc_x_min_sc ! bitp lower limit, PbSc
        c          udetpar(37) = emc_x_max_sc ! bitp upper limit, PbSc
        c          udetpar(38) = emc_x_step ! bitp stepsize, PbSc
        c          udetpar(39) = 0.0     !
        c          udetpar(40) = emc_dele_max_sc ! bitp dE upper limit, PbSc
        c          udetpar(41) = emc_dele_step_sc ! bitp dE upper limit, PbSc
        c          udetpar(42) = emc_tof_min ! bitp lower limit
        c          udetpar(43) = emc_tof_max ! bitp upper limit
        c          udetpar(44) = emc_tof_step ! bitp stepsize
        c          udetpar(45) = 0.0     !
        c          udetpar(46) = 0.0     !
        c          udetpar(47) = 0.0     !
        c          udetpar(48) = 0.0     !
        c          udetpar(49) = 0.0     !
        c          udetpar(50) = FLOAT(emc_ind1_max_sc) ! bitp tower ind. 
        c          udetpar(51) = FLOAT(emc_ind2_max_sc) ! bitp tower ind. 
        c          udetpar(52) = FLOAT(emc_iwall_max) ! 
        c          udetpar(53) = FLOAT(emc_itype_max) ! 
        c          udetpar(54) = FLOAT(emc_i1_max) ! 
        c          udetpar(55) = 0.0     !
        c          udetpar(56) = 0.0     !
        c          udetpar(57) = 0.0     !
        c          udetpar(58) = 0.0     !
        c          udetpar(59) = 0.0     !
        c          udetpar(60) = FLOAT(emc_itrack_max) ! 
        c          udetpar(61) = FLOAT(emc_spart_max) ! 
        c          udetpar(62) = FLOAT(emc_ncycle_max) ! 
        c          udetpar(63) = 0.0     !
        c          udetpar(64) = 0.0     !
        c          udetpar(65) = emc_cutgam ! 
        c          udetpar(66) = emc_cutele ! 
        c          udetpar(67) = emc_cutneu ! 
        c          udetpar(68) = emc_cuthad ! 
        c          udetpar(69) = emc_cutmuo ! 
        c          udetpar(70) = 0.0     !
        c          udetpar(71) = 0.0     !
        c          udetpar(72) = 0.0     !
        c          udetpar(73) = 0.0     !
        c          udetpar(74) = 0.0     !
        c          udetpar(75) = 0.0     !
        c          udetpar(76) = 0.0     !
        c          udetpar(77) = 0.0     !
        c          udetpar(78) = 0.0     !
        c          udetpar(79) = 0.0     !
        c          udetpar(80) = 0.0     !
        c
        c     Field added to ra_det
        c     (eventually to be overwritten by user parameters
        c
        c          ra_det(81,j) = COSD(phi angle)
        c          ra_det(82,j) = SIND(phi angle)
        c          ra_det(83,j) = Attenuation length
        c          ra_det(84,j) = Speed of light (cm/ns)
        c
        c     Get the very first detector parameter set
        c
        cgd***          CALL GFDETU('EMC ','EC11',nudetpar,nw,udetpar(1))
        cgd***          IF(nw <= 0) THEN
        cgd***           CALL EMC_GINI_FILL()
        cgd***           GOTO 1117
        cgd***          ENDIF
        c       
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

/*     Executable */


   double d_theta,d_phi,d_work,d_work0,d_work1,d_work2;
   double todeg,torad;



   torad = M_PI / 180.0;
   todeg = 180.0 / M_PI;

   /* This routine produces the default dEmcGeometry and dEmcGeaParams
      tables if no PISA input is present - at the request of Jeff 
      Mitchell, Jule 16-18 1999, G. David
      Eventually it should become obsolete when the same data will
      be fetched from the database
   */


   if(dEmcGeaParams_h->nok > 0)
   {
      return(STAFCV_OK);
   }

   /* This has been put in because of initialization problem */
   for (ll = 0; ll < sizeof(ra_det)/sizeof(ra_det[0][0]); ll++)
     {
       *((float *)ra_det + ll) = 0.0;
     }
   
   /* --------------------------------------------- */
   
   /* Put in values by hand */

   /* First wall */

   for(j=0;j<8;j++)
     {
       /* Number of sectors */
       ra_det[0][j] = 8.0;
       /* emc-opt */
       ra_det[1][j] = 3.0;
       /* Sector number */
       ra_det[2][j] = (float) j+1;
     }
   /* Sector type */


   for(j=0;j<6;j++){ ra_det[3][j]=1.0;}
   for(j=6;j<8;j++){ ra_det[3][j]=2.0;}
   /* Sector phi angle - 180 degrees!!! */
   ra_det[4][0] = 202.632;
   ra_det[4][1] = 180.000;
   ra_det[4][2] = 157.368;
   ra_det[4][3] = 134.738;
   ra_det[4][4] =  45.264;
   ra_det[4][5] =  22.632;
   ra_det[4][6] =   0.000;
   ra_det[4][7] = -22.000;
   /* Sector rpos */
   for(j=0;j<6;j++){ ra_det[5][j]= 510.0;}
   for(j=6;j<8;j++){ ra_det[5][j]= 552.64;}
   /* Sector zc_start */
   for(j=0;j<6;j++){ ra_det[6][j]= -196.74;}
   for(j=6;j<8;j++){ ra_det[6][j]= -193.23;}
      /* Correction following Maxim Volkov's suggestion
         G. David, Nov. 21, 2000 */
      if(j>5) {ra_det[5][j] = 543.2;}

   /* Sector yc_start */
   ra_det[7][0] = -285.771;
   ra_det[7][1] =  -96.985;
   ra_det[7][2] =  106.738;
   ra_det[7][3] =  294.022;
   ra_det[7][4] =  294.022;
   ra_det[7][5] =  106.738;
   ra_det[7][6] =  -95.598;
   ra_det[7][7] =  -291.674;

   /* Longitudinal size */
   for(j=0;j<6;j++){ ra_det[8][j]= 38.934;}
   for(j=6;j<8;j++){ ra_det[8][j]= 51.200;}
   /* Transverse size */
   for(j=0;j<6;j++){ ra_det[9][j]=  5.542;}
   /* Maxim Volkov discovered a discrepancy between actual size in
      PbGl tower size in PISA and the size passed in GeaParams.  
      True sizes are in z direction 4.104 cm and in y direction 4.1105 cm.
      Since the towers are everywhere assumed to be a square and
      using the z dimension in y only leaves max. 3 mm discrepancy,
      we use the z true z dimension for both z and y now.
      This has to disappear anyway a.s.a.p. when we introduce the
      true (survey) geometry and propagate it back to PISA.
      G. David, Nov. 21, 2000
   for(j=6;j<8;j++){ ra_det[9][j]=  4.068;}
   */
   for(j=6;j<8;j++){ ra_det[9][j]=  4.104;}


   /* Modules in z in supermodule */
   for(j=0;j<6;j++){ ra_det[10][j]=  12.0;}
   for(j=6;j<8;j++){ ra_det[10][j]=   6.0;}
   /* Modules in y in supermodule */
   for(j=0;j<6;j++){ ra_det[11][j]=  12.0;}
   for(j=6;j<8;j++){ ra_det[11][j]=   4.0;}
   /* Supermodules in z */
   for(j=0;j<6;j++){ ra_det[12][j]=   6.0;}
   for(j=6;j<8;j++){ ra_det[12][j]=  16.0;}
   /* Supermodules in y */
   for(j=0;j<6;j++){ ra_det[13][j]=   3.0;}
   for(j=6;j<8;j++){ ra_det[13][j]=  12.0;}
   /* Reserved fields , zero now */
   for(j=0;j<8;j++){ ra_det[14][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[15][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[16][j]=   0.0;}
   /* Scint active medium */
   for(j=0;j<6;j++){ ra_det[17][j]= 881.0;}
   for(j=6;j<8;j++){ ra_det[17][j]= 810.0;}
   /* Reserved fields , zero now */
   for(j=0;j<8;j++){ ra_det[18][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[19][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[20][j]=   0.0;}
   /* Emc debug option */
   for(j=0;j<8;j++){ ra_det[21][j]=   1.0;}
   /* GCUTS default values */
   for ( i = 0; i < 5; i++ )
     {
       ra_det[22+i][j] = 0.001;
     }
   /* Reserved fields , zero now */
   for(j=0;j<8;j++){ ra_det[27][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[28][j]=   0.0;}
   /* Min. radius for bitpacking */
   for(j=0;j<6;j++){ ra_det[29][j]= 510.0;}
   for(j=6;j<8;j++){ ra_det[29][j]= 535.0;}
   /* Max. radius for bitpacking */
   for(j=0;j<6;j++){ ra_det[30][j]= 580.0;}
   for(j=6;j<8;j++){ ra_det[30][j]= 605.0;}
   /* r stepsize for bitpacking */
   for(j=0;j<8;j++){ ra_det[31][j]=   0.07;}
   /* Min. z for bitpacking */
   for(j=0;j<8;j++){ ra_det[32][j]= -200.0;}
   /* Max. z for bitpacking */
   for(j=0;j<8;j++){ ra_det[33][j]=  200.0;}
   /* z stepsize for bitpacking */
   for(j=0;j<8;j++){ ra_det[34][j]=   0.1;}
   /* Min. x for bitpacking */
   for(j=0;j<6;j++){ ra_det[35][j]=  280.0;}
   for(j=6;j<8;j++){ ra_det[35][j]=  460.0;}
   /* Max. x for bitpacking */
   for(j=0;j<6;j++){ ra_det[36][j]=  560.0;}
   for(j=6;j<8;j++){ ra_det[36][j]=  740.0;}
   /* x stepsize for bitpacking */
   for(j=0;j<6;j++){ ra_det[37][j]=    0.07;}
   for(j=6;j<8;j++){ ra_det[37][j]=    0.07;}
   /* Reserved field , zero now */
   for(j=0;j<8;j++){ ra_det[38][j]=   0.0;}
   /* Max. E loss for bitpacking */
   for(j=0;j<6;j++){ ra_det[39][j]=    0.0327;}
   for(j=6;j<8;j++){ ra_det[39][j]=    0.256;}
   /* E loss stepsize for bitpacking */
   for(j=0;j<6;j++){ ra_det[40][j]=    1.0e-6;}
   for(j=6;j<8;j++){ ra_det[40][j]=    8.0e-6;}
   /* Tofmin for bitpacking */
   for(j=0;j<8;j++){ ra_det[41][j]=   17.0;}
   /* Tofmax for bitpacking */
   for(j=0;j<8;j++){ ra_det[42][j]=  345.0;}
   /* Tofstep for bitpacking */
   for(j=0;j<8;j++){ ra_det[43][j]=    0.01;}
   /* Reserved fields , zero now */
   for(j=0;j<8;j++){ ra_det[44][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[45][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[46][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[47][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[48][j]=   0.0;}
   /* Max. number of supermodules in sector */
   for(j=0;j<6;j++){ ra_det[49][j]=  18.0;}
   for(j=6;j<8;j++){ ra_det[49][j]= 192.0;}
   /* Max. number of towers in supermodule */
   for(j=0;j<6;j++){ ra_det[50][j]= 144.0;}
   for(j=6;j<8;j++){ ra_det[50][j]=  24.0;}
   /* Max. number of sectors for bitpacking */
   for(j=0;j<8;j++){ ra_det[51][j]=   8.0;}
   /* Max. number of sector types for bitpacking */
   for(j=0;j<8;j++){ ra_det[52][j]=   3.0;}
   /* Max. number ... for bitpacking */
   for(j=0;j<8;j++){ ra_det[53][j]=   8.0;}
   /* Reserved fields , zero now */
   for(j=0;j<8;j++){ ra_det[54][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[55][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[56][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[57][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[58][j]=   0.0;}
   /* Max. number of tracks for bitpacking */
   for(j=0;j<8;j++){ ra_det[59][j]= 32767.0;}
   /* Max. number of spart for bitpacking */
   for(j=0;j<8;j++){ ra_det[60][j]=    62.0;}
   /* Max. number of ncycle for bitpacking */
   for(j=0;j<8;j++){ ra_det[61][j]=    31.0;}
   /* Reserved fields , zero now */
   for(j=0;j<8;j++){ ra_det[62][j]=   0.0;}
   for(j=0;j<8;j++){ ra_det[63][j]=   0.0;}
   /* AuAu special cuts */
   for(j=0;j<8;j++){ ra_det[64][j]=   0.001;}
   for(j=0;j<8;j++){ ra_det[65][j]=   0.001;}
   for(j=0;j<8;j++){ ra_det[66][j]=   0.001;}
   for(j=0;j<8;j++){ ra_det[67][j]=   0.001;}
   for(j=0;j<8;j++){ ra_det[68][j]=   0.001;}
   for ( i = 0; i < 11; i++ )    /*      Last elements of array */
     {
       for(j=0;j<6;j++) 
         if(i<8)
           {
             ra_det[69+i][j] = 1.0e-4;
           }
       else
         {
         ra_det[69+i][j] = 0.0;
         }
       for(j=6;j<8;j++) ra_det[69+i][j] = 0.0;

     }  

   /* End putting in values by hand.  Now calculate derived quantities */

   for ( j = 0; j < emc_walls; j++ )
   {


     if( ra_det[70][0] == 0.0)
       {
         l_automatic_threshold =  FALSE;
       }
     if(l_automatic_threshold)
       {
         r_sampfrac = 0.2285;        /*  To make sure you always get something  */
         ind1 = 0;
         ind2 = 0 ;
         if(j < 6)
           {
             if(ra_det[70][j] >= 0.0029) ind2 = 1;
             if(ra_det[70][j] >= 0.0049) ind2 = 2;
             if(ra_det[70][j] >= 0.0099) ind2 = 3;

             if(ra_det[22][j] <= 0.0001) ind1 = 0;
             if(ra_det[22][j] > 0.0001 && ra_det[22][j] <= 0.005) ind1 = 1;
             if(ra_det[22][j] > 0.005 && ra_det[22][j] <= 0.02) ind1 = 2;
             if(ra_det[22][j] > 0.02 && ra_det[22][j] <= 0.1) ind1 = 3;

             r_sampfrac = ra_newcalib[ind1][ind2];
           }
       }
     else
       {
         r_sampfrac = .2443;
         r_cutgam = ra_det[22][j];
         if(r_cutgam <= ra_sampfrac[0][0] )
           {
             r_sampfrac = ra_sampfrac[0][1];
           }
         else
           {
             if(r_cutgam >= ra_sampfrac[22][0] )
               {
                 r_sampfrac = ra_sampfrac[22][1];
               }
             else
               {
                 
                 for ( i = 1; i < 22; i++)
                 {
                   if(r_cutgam == ra_sampfrac[i][0] )
                     {
                       r_sampfrac = ra_sampfrac[i][1];
                     }
                   else
                     {
                       if ( r_cutgam > ra_sampfrac[i][0] && r_cutgam < ra_sampfrac[i+1][0] )
                         {
                           r_sampfrac = ra_sampfrac[i+1][1] + 
                             ( ra_sampfrac[i][1] - ra_sampfrac[i+1][1] ) *
                             ( ra_sampfrac[i+1][0]-r_cutgam ) /
                             ( ra_sampfrac[i+1][0] - ra_sampfrac[i][0] );
                         }
                     }
                 }
               }
           }
       }
     if(r_sampfrac == 0.0) r_sampfrac = 0.2054;
     
     d_work = ra_det[4][j];
     d_work = torad * d_work;
     d_work1 = cos(d_work);
     d_work2 = sin(d_work);
     ra_det[80][j] = d_work1;
     ra_det[81][j] = d_work2;
     
     ra_det[82][j] = 200.0;                /*    Attenuation length     */
     ra_det[83][j] = 17.0;                 /*    Speed of light in medium       */
     ra_det[86][j] = 0.0;                  /*    Noise  */
     ra_det[87][j] = 1.0;                  /*    MeV/ADC (FAKE VALUE!   !       !       )       */
     ra_det[89][j] = 10.0;                /*     Pulse decay time (ns)  */
     ra_det[90][j] = 10.0;                /*     Tune for energy dependence     */
     ra_det[91][j] = 3.0;                 /*     Delay (?)      */
     ra_det[92][j] = 3.0;                 /*     tfac (?)       */
     
     ra_det[93][j] = 0.04;                 /*    TOF thresh     */
     
     ra_det[94][j] = 3.0;                 /*     Time to reach max at 1 GeV     */
     if(j < 6)
       {
         ra_det[84][j] = r_sampfrac;                /*   Sampling fraction      */
         
         /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            c ! This is for 1 MeV tracking
            c ! threshold
            c ! Correct it later with
            c ! proper functional form
            c
            ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
         
         ra_det[84][j] = ra_det[84][j] * 0.87; /*        To correct for other losses    */
         ra_det[85][j] = 1500.0;              /*         Photoelectrons er GeV  */
       }
     else
       {
         ra_det[84][j] = 1.0;
         ra_det[85][j] = 600.0;              /*  Photoelectrons er GeV  */
       }
     if( ra_det[90][j] > 0.0)
       {
         r_work1 = ra_det[90][j] * log(2.0);
         ra_det[92][j] = ra_det[94][j] * r_work1 / (1.0 + r_work1);
       }
     
     /* Put this is STAF table  */
     
     for ( i = 0; i < 120; i++ )
       {
         dEmcGeaParams[j].detarray[i] = ra_det[i][j];
       }
   }        /*   End of loop over j = 1, emc_walls, filling ra_det      */
   
   dEmcGeaParams_h->nok = emc_walls;
   
   /*   Now fill emc_geom       */
   
   for (ll = 0; ll < sizeof(emc_geom)/sizeof(emc_geom[0][0][0][0]); ll++)
     {
       *((float *)emc_geom + ll) = 0.0;
     }

   
   
   for ( i1 = 0; i1 < emc_walls; i1++ )
   {
      for ( iz = 0; iz <  ra_det[10][i1]*ra_det[12][i1]; iz++ )  /*      Loop in z, all tow.    */
      {
         for ( iy = 0; iy <  ra_det[11][i1]*ra_det[13][i1]; iy++ ) /*    Loop in y, all tow.    */
         {
            
/*      Apparent z position is the center of the first cell + ...       */

            r_appz = ra_det[6][i1] +   iz  * ra_det[9][i1];
            
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   c Apparent y position is the center of the first cell + ...
   c times the COSine of the wall inclination
   c
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

            d_work1 = ra_det[80][i1];
            d_work = fabs(d_work1);

            r_appy = ra_det[7][i1] +   iy  * ra_det[9][i1] * d_work;
            
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   c Apparent x position is
   c RPOS * COSine of wall inclination
   c - YPOSisiton in the coo.system of the wall * SINe of
   c wall inclination
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

            r_appx = ra_det[5][i1] * d_work - 
                    (- (ra_det[11][i1]*ra_det[13][i1] - 1.0) 
                     * ra_det[9][i1] / 2.0
                     +   iy  * ra_det[9][i1]) * ra_det[81][i1];

            if(i1 > 3) r_appx = - r_appx;
            d_work1 = r_appx*r_appx + r_appy*r_appy;
            d_work1 = sqrt(d_work1);
            d_work2 = r_appz;
            d_work0 = r_appx;
            
            if(d_work1 > 0.0)
            {

               d_work = atan2(d_work1,d_work2);
               r_measthe = d_work;
               r_measthe = todeg * r_measthe;
               
               d_work = acos(d_work0/d_work1);
               r_measphi = d_work;
               r_measphi = todeg * r_measphi;
               
               if(r_appy < 0.0) r_measphi = 360.0 - r_measphi;
            }


            emc_geom[0][iz][iy][i1] = r_appx;
            emc_geom[1][iz][iy][i1] = r_appy;
            emc_geom[2][iz][iy][i1] = r_appz;
            emc_geom[3][iz][iy][i1] = r_measthe;
            emc_geom[4][iz][iy][i1] = r_measphi;
            d_work1 = r_appx*r_appx + r_appy*r_appy + r_appz*r_appz;
            d_work = sqrt(d_work1);
            r_work1 = d_work;
            
            emc_geom[5][iz][iy][i1] = r_work1;
            emc_geom[6][iz][iy][i1] = r_work1 / 30.0;    /*      Photon flash   */

            /* Write unit vector corresponding to the acis of the tower */

            emc_geom[7][iz][iy][i1] = ra_det[80][i1];
            emc_geom[8][iz][iy][i1] = ra_det[81][i1];
            emc_geom[9][iz][iy][i1] = 0.0;
            
            
                 
         }/*     Loop over iy   */
      }/*        Loop over iz   */
   }/*   Loop over i1   */

/*      Write dEmcGeometry STAF table  
        the array emc_geom is indexed by sectors 1-8 at this point,
        but dEmcGeometry is written out according to a different
        numbering scheme
*/

   istaf = 0;
   for ( i1 = 0; i1 < emc_walls; i1++)
   {
      for ( iy = 0; iy < max_chany; iy++ )
      {
         for ( iz = 0; iz < max_chanz; iz++ )
         {
            if(emc_geom[5][iz][iy][i1] > 0.0)
            {
               dEmcGeometry[istaf].id = istaf + 1;
               dEmcGeometry[istaf].hwkey = 0;
               if(i1 < 4)
               {
                  dEmcGeometry[istaf].arm = 0;
                  dEmcGeometry[istaf].sector    = i1;
               }
               else
               {
                 /* Sector numbering from bottom up in East Arm */
                  dEmcGeometry[istaf].arm = 1;
                  dEmcGeometry[istaf].sector    = 7 - i1;
               }

               /*  Numbering reversed in West Arm */
               if(dEmcGeometry[istaf].arm == 0  && iz <= 71)
                 {
                   dEmcGeometry[istaf].ind[0]    = 71 - iz;
                 }
               else
                 {
                   dEmcGeometry[istaf].ind[0]    = iz;
                 }

               dEmcGeometry[istaf].ind[1]    = iy;

               dEmcGeometry[istaf].swkey = 
                 100000 * dEmcGeometry[istaf].arm +
                 10000 *  dEmcGeometry[istaf].sector +
                 100 *    iy+ 
                          iz;


               for (i = 0; i < 3; i++)
                 {
                   dEmcGeometry[istaf].nomxyz[i] = emc_geom[i][iz][iy][i1];
                   dEmcGeometry[istaf].actxyz[i] = emc_geom[i][iz][iy][i1];
                   dEmcGeometry[istaf].nomunitv[i] = emc_geom[7+i][iz][iy][i1];
                   dEmcGeometry[istaf].actunitv[i] = emc_geom[7+i][iz][iy][i1];
                 }

               dEmcGeometry[istaf].nomtheta  = emc_geom[3][iz][iy][i1];
               dEmcGeometry[istaf].nomphi    = emc_geom[4][iz][iy][i1];
               dEmcGeometry[istaf].acttheta  = emc_geom[3][iz][iy][i1];
               dEmcGeometry[istaf].actphi    = emc_geom[4][iz][iy][i1];
               dEmcGeometry[istaf].nomdist   = emc_geom[5][iz][iy][i1];
               dEmcGeometry[istaf].actdist   = emc_geom[5][iz][iy][i1];
               dEmcGeometry[istaf].nomflash  = emc_geom[6][iz][iy][i1];
               dEmcGeometry[istaf].actflash  = emc_geom[6][iz][iy][i1];

/*      Get theta, phi within the sector (local) coordinate system) */
/*      But buddy, this is c, indices are shifted by one */

/*      Yes, and at some point one should check if iz,iy,i1 is correct... */
/*      or it should be iz+1, etc. */

               r_dist = ra_det[5][i1];
               r_distz = (  iz
                          - ra_det[10][i1] * ra_det[12][i1] / 2.0)
                  * ra_det[9][i1];
               r_disty = (  iy - ra_det[11][i1] * ra_det[15][i1]
                          - ra_det[11][i1] * ra_det[13][i1] / 2.0)
                  * ra_det[9][i1];

               d_work1 = r_distz;
               d_work2 = r_disty;
               d_work = r_dist;
               d_theta = atan2(d_work1,d_work);
               d_phi = atan2(d_work2,d_work);
               d_theta = todeg * d_theta;
               d_phi = todeg * d_phi;
               r_sectheta = d_theta;
               r_secphi = d_phi;
               dEmcGeometry[istaf].sectheta = r_sectheta;
               dEmcGeometry[istaf].secphi = r_secphi;
               
            
               
               istaf = istaf + 1;
               
            }/*  Photon flash > 0.0, tower defined      */
         }/*     Loop over iz   */
      }/*        Loop over iy   */
   }/*   Loop over iz   */
   

   dEmcGeometry_h->nok = istaf;
   
   return (STAFCV_OK);
}