mEmcGeaParams.c

Go to the documentation of this file.

#include "mEmcGeaParams.h"
#include "emlLib.h"
#include <math.h>

long 
mEmcGeaParams_(
  TABLE_HEAD_ST         *emcpar_h,         EMCPAR_ST           *emcpar ,
  TABLE_HEAD_ST  *dEmcGeaParams_h,  DEMCGEAPARAMS_ST    *dEmcGeaParams ,
  TABLE_HEAD_ST   *dEmcGeometry_h,   DEMCGEOMETRY_ST     *dEmcGeometry )
{
/*:>--------------------------------------------------------------------
**: ROUTINE:    mEmcGeaParams_
**: DESCRIPTION: Physics Analysis Module ANSI C template.
**:             This is an ANSI C Physics Analysis Module template
**:             automatically generated by stic from mEmcGeaParams.idl.
**:             Please edit comments and code.
**: ARGUMENTS:
**:       IN:
**:             emcpar    - PLEASE FILL IN DESCRIPTION HERE
**:            emcpar_h   - header Structure for emcpar
**:    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.

*/

   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;


        /*      CHARACTER*4 cudet */

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        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;

   /*  Modified sampling fractions according to recalibration
       August 1999, G. David (modifications entered Jan. 17, 2000)
       Modified values: .0010  (AUAU)
                        .0030  (PPLO)
                        .0100  (PPHI)

   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        ra_sampfrac[23][2]      =
        {
        .0001, .2443, .0002, .2432, .0005, .2313, .0008, .2204, 
        .0010, .2251, .0012, .1981, .0015, .1889, .0018, .1809,
        .0020, .1760, .0022, .1717, .0025, .1650, .0028, .1607, 
        .0030, .1872, .0040, .1429, .0050, .1307, .0060, .1214,
        .0070, .1130, .0080, .1058, .0090, .0993, .0100, .1426, 
        .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;
   

   float ra_translate[3][8];
   

   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;

   


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


   if(emcpar_h->nok < 1 || emcpar_h->nok > 8)
   {
      return ( STAFCV_BAD);
   }    

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        c     If nothing is passed in fields 71-76, you are working with
        c     the old calibration logic
        c
        cgd***          IF(udetpar(71) == 0.0)
        cgd***           l_automatic_threshold =  FALSE 
        cgd***          ENDIF
        c
        cgd***          emc_walls = INT(udetpar(1))
        c
        cgd***          DO i = 1,80
        cgd***           ra_det(i,1) = udetpar(i)
        cgd***          ENDDO
        c
        c
        cgd***          IF(emc_walls > 1)
        cgd***           DO j = 2,MIN(6,emc_walls)
        cgd***            WRITE(cudet,'(a2,i1,a1)')
        cgd***     1       'EC',j,'1'
        cgd***            CALL GFDETU('EMC ',cudet,nudetpar,nw,udetpar(1))
        cgd***            DO i = 1,80
        cgd***             ra_det(i,j) = udetpar(i)
        cgd***            ENDDO
        cgd***           ENDDO
        cgd***          ENDIF
        c
        c     If there is PbGl
        c
        cgd***          IF(emc_walls > 6)
        cgd***           DO j = 7,MIN(8,emc_walls)
        cgd***            WRITE(cudet,'(a2,i1,a1)')
        cgd***     1       'EC',j,'2'
        cgd***            CALL GFDETU('EMC ',cudet,nudetpar,nw,udetpar(1))
        cgd***            DO i = 1,80
        cgd***             ra_det(i,j) = udetpar(i)
        cgd***            ENDDO
        cgd***           ENDDO
        cgd***          ENDIF
        c
        cgd*** 1117     CONTINUE
        c
        c
        c-----    Needed for STAF
        c
        c       
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */

   emc_walls = emcpar[1].emc_walls;

        
   for (ll = 0; ll < sizeof(ra_det)/sizeof(ra_det[0][0]); ll++)
     {
       *((float *)ra_det + ll) = 0.0;
     }
   /* This has been put bacause of initialization problem       */


   for ( j = 0; j < emc_walls; j++ )
   {
      ra_det[0][j] = emcpar[j].emc_walls;
      ra_det[1][j] = emcpar[j].emc_opt;
      ra_det[2][j] = emcpar[j].iwall;
      ra_det[3][j] = emcpar[j].itype;
      ra_det[4][j] = emcpar[j].angle;
      ra_det[5][j] = emcpar[j].rpos;
      /* Correction following Maxim Volkov's suggestion
         G. David, Nov. 21, 2000 */
      if(j>5) {ra_det[5][j] = 543.2;}
        
      ra_det[6][j] = emcpar[j].zc_start;
      ra_det[7][j] = emcpar[j].yc_start;
      ra_det[8][j] = emcpar[j].lsiz;
      ra_det[9][j] = emcpar[j].tsiz;
   /* 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
   */
      if(j>5) {ra_det[9][j]=  4.104;}
   

      ra_det[10][j] = emcpar[j].no_modz;
      ra_det[11][j] = emcpar[j].no_mody;
      ra_det[12][j] = emcpar[j].no_smodz;
      ra_det[13][j] = emcpar[j].no_smody;
      ra_det[14][j] = 0.0;
      ra_det[15][j] = 0.0;
      ra_det[16][j] = 0.0;
      ra_det[17][j] = emcpar[j].scint_emc_med;
      ra_det[18][j] = 0.0;
      ra_det[19][j] = 0.0;
      ra_det[20][j] = 0.0;
      ra_det[21][j] = emcpar[j].emc_debug;

      /* These are the CUTGAM, CUTELE, etc. values as defined
         in the CUTS data card */

      for ( i = 0; i < 5; i++ )
      {
         ra_det[22+i][j] = emcpar[j].gcuts[i];
      }
      ra_det[27][j] = 0.0;
      ra_det[28][j] = 0.0;
      ra_det[29][j] = emcpar[j].emc_r_min_sc;
      ra_det[30][j] = emcpar[j].emc_r_max_sc;
      ra_det[31][j] = emcpar[j].emc_r_step;
      ra_det[32][j] = emcpar[j].emc_z_min;
      ra_det[33][j] = emcpar[j].emc_z_max;
      ra_det[34][j] = emcpar[j].emc_z_step;
      ra_det[35][j] = emcpar[j].emc_x_min_sc;
      ra_det[36][j] = emcpar[j].emc_x_max_sc;
      ra_det[37][j] = emcpar[j].emc_x_step;
      ra_det[38][j] = 0.0;
      ra_det[39][j] = emcpar[j].emc_dele_max_sc;
      ra_det[40][j] = emcpar[j].emc_dele_step_sc;
      ra_det[41][j] = emcpar[j].emc_tof_min;
      ra_det[42][j] = emcpar[j].emc_tof_max;
      ra_det[43][j] = emcpar[j].emc_tof_step;
      ra_det[44][j] = 0.0;
      ra_det[45][j] = 0.0;
      ra_det[46][j] = 0.0;
      ra_det[47][j] = 0.0;
      ra_det[48][j] = 0.0;
      ra_det[49][j] = emcpar[j].emc_ind1_max_sc;
      ra_det[50][j] = emcpar[j].emc_ind2_max_sc;
      ra_det[51][j] = emcpar[j].emc_iwall_max;
      ra_det[52][j] = emcpar[j].emc_itype_max;
      ra_det[53][j] = emcpar[j].emc_i1_max;
      ra_det[54][j] = 0.0;
      ra_det[55][j] = 0.0;
      ra_det[56][j] = 0.0;
      ra_det[57][j] = 0.0;
      ra_det[58][j] = 0.0;
      ra_det[59][j] = emcpar[j].emc_itrack_max;
      ra_det[60][j] = emcpar[j].emc_spart_max;
      ra_det[61][j] = emcpar[j].emc_ncycle_max;
      ra_det[62][j] = 0.0;
      ra_det[63][j] = 0.0;
      ra_det[64][j] = emcpar[j].emc_cutgam;
      ra_det[65][j] = emcpar[j].emc_cutele;
      ra_det[66][j] = emcpar[j].emc_cutneu;
      ra_det[67][j] = emcpar[j].emc_cuthad;
      ra_det[68][j] = emcpar[j].emc_cutmuo;
      for ( i = 0; i < 11; i++ )    /*   Last elements of array */
      {
         ra_det[69+i][j] = emcpar[j].array[21+i];
      } 
      /* The following data could also be put in ra_det[14+i][j] */
      for ( i = 0; i < 3; i++ )    /*    Carriage translation   */
      {
         ra_translate[i][j] = emcpar[j].array[i];
         ra_det[14+i][j] = emcpar[j].array[i];
      } 
      /*      
              ra_det[79] is "emc_response_option"
      */

      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)
         {
           /* Screw this damn Linux that doesn't know what == means
              i.e. the comparison precision is too high,
              0.01 is not equal to 0.00999999978, but then again,
              it is unable to write out 0.01...  Jan 17, 2000 G. David */
           /*
            if(ra_det[70][j] == 0.003) ind2 = 1;
            if(ra_det[70][j] == 0.01) ind2 = 2;
           */
           /* Modified to accomodate 'PPME' setting
              G. David, February 11, 2000

            if(ra_det[70][j] >= 0.0029) ind2 = 1;
            if(ra_det[70][j] >= 0.0099) ind2 = 2;
           */

            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.2256;
      
/* --------------------------------------
      ra_det[80][j] = cosd(ra_det[4][j]);
      ra_det[81][j] = sind(ra_det[4][j]);
------------------------------------------*/
      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;
      
         

      /* Changed at Sasha Bazilevsky's request, Dec. 5,2000 G. David*/
      ra_det[82][j] = 120.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 (?)       */
      
      /*        cgd*** ra_det(94,j) = 0.05  ! TOF thresh        */
           
      ra_det[93][j] = 0.04;                 /*   TOF thresh     */
      
      /*        cgd*** ra_det(94,j) = 0.08 ! 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 < ( int ) (ra_det[10][i1]*ra_det[12][i1]); iz++ )  /*     Loop in z, all tow.    */
      {
         for ( iy = 0; iy < ( int ) (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] + ( float ) ( 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
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */


/*----------
            r_appy = ra_det[7][i1] + ( float ) ( iy ) * ra_det[9][i1] * fabsf(ra_det[80][i1]);
------------*/
            d_work1 = ra_det[80][i1];
            d_work = fabs(d_work1);

            r_appy = ra_det[7][i1] + ( float ) ( 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] * fabsf(ra_det[80][i1]) - 
                    (- (ra_det[11][i1]*ra_det[13][i1] - 1.0) 
                     * ra_det[9][i1] / 2.0
                     + ( float ) ( iy ) * ra_det[9][i1]) * ra_det[81][i1];
------------*/
            r_appx = ra_det[5][i1] * d_work - 
                    (- (ra_det[11][i1]*ra_det[13][i1] - 1.0) 
                     * ra_det[9][i1] / 2.0
                     + ( float ) ( iy ) * ra_det[9][i1]) * ra_det[81][i1];

/*          if(i1 < 4) r_appx = - r_appx;    */

            if(i1 > 3) r_appx = - r_appx;

            /* Add carriage translation */

            r_appx = r_appx + ra_translate[0][i1];
            

/*          r_work1 = sqrtf(r_appx*r_appx + r_appy*r_appy);  */
            d_work1 = r_appx*r_appx + r_appy*r_appy;
            d_work1 = sqrt(d_work1);
/*          r_work1 = d_work1;  */
            d_work2 = r_appz;
            d_work0 = r_appx;
            
/*----------
        c    if(r_work1 > 0.0)
        c    {
        c       r_measthe = atan2d(r_work1,r_appz);
        c       r_measphi = acosd(r_appx/r_work1);
        c       if(r_appy < 0.0) r_measphi = 360.0 - r_measphi;
        c    }
------------*/
            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;
/*          r_work1 = sqrtf(r_appx*r_appx + r_appy*r_appy + r_appz*r_appz); */
            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 = ( (float) (iz)
                          - ra_det[10][i1] * ra_det[12][i1] / 2.0)
                  * ra_det[9][i1];
               r_disty = ( (float) (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);
}