mEmcGeaTowerEval.c

Go to the documentation of this file.

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

long mEmcGeaTowerEval_(
  TABLE_HEAD_ST   *dEmcGeometry_h,   DEMCGEOMETRY_ST     *dEmcGeometry ,
  TABLE_HEAD_ST      *dEmcEvent_h,      DEMCEVENT_ST        *dEmcEvent ,
  TABLE_HEAD_ST *dEmcCalibTower_h, DEMCCALIBTOWER_ST   *dEmcCalibTower ,
  TABLE_HEAD_ST   *dEmcGeaTrack_h,   DEMCGEATRACK_ST     *dEmcGeaTrack ,
  TABLE_HEAD_ST *dEmcGeaTowerTrack_h, DEMCGEATOWERTRACK_ST *dEmcGeaTowerTrack ,
  TABLE_HEAD_ST *dEmcGeaTowerEval_h, DEMCGEATOWEREVAL_ST *dEmcGeaTowerEval )
{
/*:>--------------------------------------------------------------------
**: ROUTINE:    mEmcGeaTowerEval_
**: DESCRIPTION: Physics Analysis Module ANSI C template.
**:             This is an ANSI C Physics Analysis Module template
**:             automatically generated by stic from mEmcGeaTowerEval.idl.
**:             Please edit comments and code.
**: AUTHOR:     hpl - H.P. Lovecraft, hplovecraft@cthulhu.void
**: ARGUMENTS:
**:       IN:
**:       dEmcGeometry    - PLEASE FILL IN DESCRIPTION HERE
**:      dEmcGeometry_h   - header Structure for dEmcGeometry
**:          dEmcEvent    - PLEASE FILL IN DESCRIPTION HERE
**:         dEmcEvent_h   - header Structure for dEmcEvent
**:     dEmcCalibTower    - PLEASE FILL IN DESCRIPTION HERE
**:    dEmcCalibTower_h   - header Structure for dEmcCalibTower
**:       dEmcGeaTrack    - PLEASE FILL IN DESCRIPTION HERE
**:      dEmcGeaTrack_h   - header Structure for dEmcGeaTrack
**:  dEmcGeaTowerTrack    - PLEASE FILL IN DESCRIPTION HERE
**:  dEmcGeaTowerTrack_h   - header Structure for dEmcGeaTowerTrack
**:    INOUT:
**:      OUT:
**:   dEmcGeaTowerEval    - PLEASE FILL IN DESCRIPTION HERE
**:  dEmcGeaTowerEval_h   - header Structure for dEmcGeaTowerEval
**: RETURNS:    STAF Condition Value
**:>------------------------------------------------------------------*/

  int i_twrev_id;
  int i_calib_id,i_gtr_id;
  int ll;
  int i,j,k;
  int twrkey,arm,sector,indy,indz;

  /* ra_track[j][0] -> PID
     ra_track[j][1] -> ptot
     ra_track[j][2] -> vertex x
     ra_track[j][3] -> vertex y
     ra_track[j][4] -> vertex z
     ra_track[j][5] -> ancestry level
     ra_track[j][6] ->               */
  float ra_track[3][7];

  /* ia_track[j][0] -> true track number
     ia_track[j][1] -> true track number of parent
     ia_track[j][2] -> PID of parent */
  int   ia_track[3][3];

  static int l_first = 1;
  static float emc_geom[2][4][48][96][3];
  float event_vertex[3];
  
  double d_work1,d_work2,d_work3,d_work4,d_work;
  float r_work1;
  static float r_speedoflight = 30.0;

  
  int size;

  /* Sanity check */

  if(dEmcGeometry_h->nok <=0) return (STAFCV_BAD);
  if(dEmcEvent_h->nok <=0) return (STAFCV_BAD);
  if(dEmcGeaTrack_h->nok <=0) return (STAFCV_BAD);
  if(dEmcGeaTowerTrack_h->nok <=0) return (STAFCV_BAD);

  /* Initalization */

  if(l_first == 1)
    {
      l_first = 0;
      for (ll = 0; ll < sizeof(emc_geom)/sizeof(emc_geom[0][0][0][0][0]); ll++)
        {
          *((float *)emc_geom + ll) = 0.0;
        }
      for( i = 0; i < dEmcGeometry_h->nok; i++)
        {
          arm = dEmcGeometry[i].arm;
          sector = dEmcGeometry[i].sector;
          indy = dEmcGeometry[i].ind[1];
          indz = dEmcGeometry[i].ind[0];
          for( j = 0; j < 3; j++)
            {
              emc_geom[arm][sector][indy][indz][j] = 
                dEmcGeometry[i].actxyz[j];
            }
        }
    }

  i_twrev_id = 0;
  for( i = 0; i < 3; i++)
    {
      event_vertex[i] = dEmcEvent[0].xyz[i];
    }
  


  /* Loop over all entries in TowerTrack */
  
  for( i = 0; i < dEmcGeaTowerTrack_h->nok; i++)
    {
      for (ll = 0; ll < sizeof(ra_track)/sizeof(ra_track[0][0]); ll++)
        {
          *((float *)ra_track + ll) = 0.0;
        }
      for (ll = 0; ll < sizeof(ia_track)/sizeof(ia_track[0][0]); ll++)
        {
          *((int *)ia_track + ll) = 0;
        }

      /* Translate software key to arm, sector, y and z index */
      
      twrkey = dEmcGeaTowerTrack[i].twrkey;
      arm = twrkey / 100000;
      sector = (twrkey - arm * 100000) / 10000;
      indy = (twrkey - arm * 100000 - sector * 10000) / 100;
      indz = twrkey - arm * 100000 - sector * 10000 - indy * 100;


      /* Get the row index of the corresponding tower in CalibTower
         to write tof, measured e */

      i_calib_id = -1;
      j = 0;
      while( i_calib_id < 0 && j < dEmcCalibTower_h->nok)
        {
          if(dEmcCalibTower[j].swkey == twrkey) i_calib_id = j;
          j++;
        }


      /* Loop over tracks in TowerTrack
         get the row index of the corresponding track in GeaTrack
         to retrieve ptot, ancestry, pid, etc. */
      j = 0;
      while( j <= 2 && dEmcGeaTowerTrack[i].trkno[j] > 0)
        {
          ia_track[j][0] = dEmcGeaTowerTrack[i].trkno[j];
          
          /* Search track and get its vertex, ancestry, etc. */
          i_gtr_id = -1;
          k = 0;
          while( i_gtr_id < 0 && k < dEmcGeaTrack_h->nok)
            {
              if( ia_track[j][0] == dEmcGeaTrack[k].trkno) i_gtr_id = k;
              k++;
            }
          if(i_gtr_id > -1)
            {
              ia_track[j][1] = dEmcGeaTrack[i_gtr_id].itparent;
              ia_track[j][2] = dEmcGeaTrack[i_gtr_id].idparent;
              ra_track[j][0] = dEmcGeaTrack[i_gtr_id].pid;
              ra_track[j][1] = dEmcGeaTrack[i_gtr_id].ptot;
              ra_track[j][2] = dEmcGeaTrack[i_gtr_id].xyz[0];
              ra_track[j][3] = dEmcGeaTrack[i_gtr_id].xyz[1];
              ra_track[j][4] = dEmcGeaTrack[i_gtr_id].xyz[2];
              ra_track[j][5] = dEmcGeaTrack[i_gtr_id].anclvl;
            } /* Entry in dEmcGeaTrack found */
          j++;
        }  /* Loop over j, trkno[j] in dEmcGeaTowerTrack */
      
      /* Write dEmcGeaTowerEval */

      dEmcGeaTowerEval[i_twrev_id].id = i_twrev_id;
      dEmcGeaTowerEval[i_twrev_id].input = dEmcEvent[0].serialno;
      dEmcGeaTowerEval[i_twrev_id].twrkey = twrkey;
      dEmcGeaTowerEval[i_twrev_id].arm = arm;
      dEmcGeaTowerEval[i_twrev_id].sector = sector;
      dEmcGeaTowerEval[i_twrev_id].indy = indy;
      dEmcGeaTowerEval[i_twrev_id].indz = indz;

      for( j = 0; j < 3; j++)
        {
          dEmcGeaTowerEval[i_twrev_id].trkno[j] = ia_track[j][0];
          dEmcGeaTowerEval[i_twrev_id].itparent[j] = ia_track[j][1];
          dEmcGeaTowerEval[i_twrev_id].idparent[j] = ia_track[j][2];
          dEmcGeaTowerEval[i_twrev_id].pid[j] = ra_track[j][0];
          dEmcGeaTowerEval[i_twrev_id].ptot[j] = ra_track[j][1];
          dEmcGeaTowerEval[i_twrev_id].vertex[j][0] = ra_track[j][2];
          dEmcGeaTowerEval[i_twrev_id].vertex[j][1] = ra_track[j][3];
          dEmcGeaTowerEval[i_twrev_id].vertex[j][2] = ra_track[j][4];
          dEmcGeaTowerEval[i_twrev_id].ancestry[j] = ra_track[j][5];
        }

      /* Calculate nominal (event vertex) and actual (track vertex)
         distances (using the highest Edep particle
         and the sin theta values */

      d_work1 = 0.0;
      d_work2 = 0.0;
      d_work3 = 0.0;
      d_work4 = 0.0;
      for( j = 0; j < 3; j++)
        {
          d_work = emc_geom[arm][sector][indy][indz][j] - event_vertex[j];
          d_work1 = d_work1 + d_work * d_work;
          if( j < 2 ) d_work3 = d_work3 + d_work * d_work;

          /* Use coordinates of the first (highest Edep) track */
          d_work = emc_geom[arm][sector][indy][indz][j] - ra_track[0][j+2];
          d_work2 = d_work2 + d_work * d_work;
          if( j < 2 ) d_work4 = d_work4 + d_work * d_work;
        }
      d_work1 = sqrt(d_work1);  /* Distance from event vertex */
      d_work2 = sqrt(d_work2);  /* Distance from track vertex */
      d_work3 = sqrt(d_work3);  /* xy plane distance from event vertex */
      d_work4 = sqrt(d_work4);  /* xy plane distance from track vertex */

      if (d_work1 > 0.0 && d_work2 > 0.0)
        {
          dEmcGeaTowerEval[i_twrev_id].dist_nom = d_work1;
          dEmcGeaTowerEval[i_twrev_id].dist_act = d_work2;
          dEmcGeaTowerEval[i_twrev_id].sinthe_nom = d_work3 / d_work1;
          dEmcGeaTowerEval[i_twrev_id].sinthe_act = d_work4 / d_work2;
        }
      else      
        {
          dEmcGeaTowerEval[i_twrev_id].dist_nom = 0.0;
          dEmcGeaTowerEval[i_twrev_id].dist_act = 0.0;
          dEmcGeaTowerEval[i_twrev_id].sinthe_nom = 0.0;
          dEmcGeaTowerEval[i_twrev_id].sinthe_act = 0.0;
        }
      
      /* Transfer info from dEmcGeaTowerTrack */

      for( j = 0; j < 3; j++)
        {
          dEmcGeaTowerEval[i_twrev_id].trkno[j] = 
            dEmcGeaTowerTrack[i].trkno[j];
          dEmcGeaTowerEval[i_twrev_id].edep[j] = 
            dEmcGeaTowerTrack[i].edep[j];
          r_work1 = 
            dEmcGeaTowerTrack[i].toffirst[j] - d_work1 / r_speedoflight;
          if(r_work1 > 0.0)
            {
              dEmcGeaTowerEval[i_twrev_id].toffirst[j] = r_work1;
            }
          else
            {
              dEmcGeaTowerEval[i_twrev_id].toffirst[j] = 0.0;
            }
        }
      

      if(i_calib_id >= 0)
        {
          dEmcGeaTowerEval[i_twrev_id].mease = 
            dEmcCalibTower[i_calib_id].ecal;
          r_work1 = dEmcCalibTower[i_calib_id].tof - d_work1 / r_speedoflight;
          if(r_work1 > 0.0)
            {
              dEmcGeaTowerEval[i_twrev_id].tof = r_work1;
            }
          else
            {
              dEmcGeaTowerEval[i_twrev_id].tof = 0.0;
            }
        }
      else      
        {
          dEmcGeaTowerEval[i_twrev_id].mease = 0.0;
          dEmcGeaTowerEval[i_twrev_id].tof = 0.0;
        }
      
      i_twrev_id = i_twrev_id + 1;
      


    }  /* Loop over i, all rows of dEmcGeaTowerTrack */
  

  dEmcGeaTowerEval_h->nok = i_twrev_id;
  

  return STAFCV_OK;
}