mEmcGeaClusterEval2.c

Go to the documentation of this file.

#include "mEmcGeaClusterEval2.h"
#include "emlLib.h"

long type_of_call mEmcGeaClusterEval2_(
  TABLE_HEAD_ST      *dEmcEvent_h,      DEMCEVENT_ST        *dEmcEvent ,
  TABLE_HEAD_ST   *dEmcGeaTrack_h,   DEMCGEATRACK_ST     *dEmcGeaTrack ,
  TABLE_HEAD_ST *dEmcGeaTowerTrack_h, DEMCGEATOWERTRACK_ST *dEmcGeaTowerTrack ,
  TABLE_HEAD_ST *dEmcClusterExt_h, DEMCCLUSTEREXT_ST   *dEmcClusterExt ,
  TABLE_HEAD_ST *dEmcGeaTrackCluster_h, DEMCGEATRACKCLUSTER_ST *dEmcGeaTrackCluster ,
  TABLE_HEAD_ST *dEmcGeaClusterTrack_h, DEMCGEACLUSTERTRACK_ST *dEmcGeaClusterTrack )
{
/*:>--------------------------------------------------------------------
**: ROUTINE:    mEmcGeaClusterEval2_
**: DESCRIPTION: Physics Analysis Module ANSI C template.
**:             This is an ANSI C Physics Analysis Module template
**:             automatically generated by stic from mEmcGeaClusterEval2.idl.
**:             Please edit comments and code.
**: AUTHOR:     hpl - H.P. Lovecraft, hplovecraft@cthulhu.void
**: ARGUMENTS:
**:       IN:
**:          dEmcEvent    - PLEASE FILL IN DESCRIPTION HERE
**:         dEmcEvent_h   - header Structure for dEmcEvent
**:       dEmcGeaTrack    - PLEASE FILL IN DESCRIPTION HERE
**:      dEmcGeaTrack_h   - header Structure for dEmcGeaTrack
**:  dEmcGeaTowerTrack    - PLEASE FILL IN DESCRIPTION HERE
**:  dEmcGeaTowerTrack_h   - header Structure for dEmcGeaTowerTrack
**:     dEmcClusterExt    - PLEASE FILL IN DESCRIPTION HERE
**:    dEmcClusterExt_h   - header Structure for dEmcClusterExt
**:    INOUT:
**:      OUT:
**:  dEmcGeaTrackCluster    - PLEASE FILL IN DESCRIPTION HERE
**:  dEmcGeaTrackCluster_h   - header Structure for dEmcGeaTrackCluster
**:  dEmcGeaClusterTrack    - PLEASE FILL IN DESCRIPTION HERE
**:  dEmcGeaClusterTrack_h   - header Structure for dEmcGeaClusterTrack
**: RETURNS:    STAF Condition Value
**:>------------------------------------------------------------------*/



  int i_cltr_id,i_trcl_id;
  int i,j,j1,j2,k,k1,k2;
  int ll;
  int l_found,l_found1,l_found2;
  int i_track;
  
  
   int size;

   /* ia_track will store the track numbers of the tracks contributing
      to a specific cluster 
   */
   int   ia_track[50];

   /* ra_track will store data of tracks listed in ia_track.
      Field: 0 - energy (cumulates "edep" field of TowerTrack)
             1 - pid
             1 - ptot
             3-5 - impact xyz
             6 - twrhit
             7 - edep ("edep" field of GeaTrack)
             8 - anclvl
             9-11 - vertex xyz
   */
   float ra_track[12][50];
   
   int i_size;
   int i_size1;
   int i_size2;
   int i_size3;
   
   int itemp,i_twrkey;
   float r_e;
   float r_work;
   float r_edep;
   
   float ra_work[12];
   int   i_work1;
   float nom_edep;
   
   int ia_tower[4][16];  
   /* GEANT track number 1,2,3 and towerindex for max. 16 contributors */
   float ra_tower[5][16];
   /* GEANT deposited energy 1,2,3, attributed energy, total GEANT energy */


   double ytemp1,ytemp2,ztemp1,ztemp2,ydist,zdist;

   int isec,iindy,iindz;
   


   /*     Executable  */



  if(dEmcGeaTrack_h->nok <=0) return (STAFCV_BAD);
  if(dEmcGeaTowerTrack_h->nok <=0) return (STAFCV_BAD);
  if(dEmcClusterExt_h->nok <=0) return (STAFCV_BAD);
  
  i_cltr_id = 0;
  i_size = 600;
  i_size1 = 50;
  i_size2 = 64;
  i_size3 = 80;
  
  
  /* First part.  Figure out cluster overlaps, i.e. which track deposited */
  /* energy in a specific cluster */

  /* Loop over all clusters found */

  for( i = 0; i < dEmcClusterExt_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]); ll++)
        {
          *(ia_track + ll) = 0;
        }
      for (ll = 0; ll < sizeof(ia_tower)/sizeof(ia_tower[0][0]); ll++)
        {
          *((int *)ia_tower + ll) = 0;
        }
      for (ll = 0; ll < sizeof(ra_tower)/sizeof(ra_tower[0][0]); ll++)
        {
          *((float *)ra_tower + ll) = 0.0;
        }

      /* Loop over all towers included in a cluster */


      for( j = 0; j < 16 && dEmcClusterExt[i].partesum[j] > 0.0; j++)
        {
          i_twrkey = dEmcClusterExt[i].twrlist[j];
          if(j==0)
            {
              ra_tower[3][j] = dEmcClusterExt[i].partesum[j];
            }
          else
            {
              ra_tower[3][j] = dEmcClusterExt[i].partesum[j]
                - dEmcClusterExt[i].partesum[j-1];
            }
          /* Energy in the tower attributed to this cluster */
          
          /* Loop over all tracks contributing to the tower */

          l_found = -1;
          k = 0;
          while( l_found < 0 && k < dEmcGeaTowerTrack_h->nok )
            {
              if( i_twrkey == dEmcGeaTowerTrack[k].twrkey)
                {
                  l_found = k;
                }
              k = k + 1;
            } /* while Loop over GeaTowerTrack to find contributor */


          if( l_found >= 0)  /* Contributor found */
            {
              ia_tower[3][j] = i_twrkey;
              
              for( k = 0; k <= 2; k++)
                {
                  i_track = dEmcGeaTowerTrack[l_found].trkno[k];
                  ia_tower[k][j] = i_track;
                  ra_tower[k][j] = dEmcGeaTowerTrack[l_found].edep[k];
                  ra_tower[4][j] = ra_tower[4][j] + ra_tower[k][j];
                  
                  if( i_track > 0)
                    {
                      l_found1 = -1;
                      k1 = 0;
                      
                      while( l_found1 < 0 && k1 < i_size1)
                        {
                          if( (ia_track[k1] == i_track) || (ia_track[k1] == 0))
                              l_found1 = k1;
                          k1 = k1 + 1;

                        }  /* while Loop over ia_track (already registered?) */

                      /* Let's hope you either found the track or there is
                         still free space in the array */

                      if( l_found1 >= 0)
                        {
                          ia_track[l_found1] = i_track;
                          ra_track[0][l_found1] = ra_track[0][l_found1] +
                            dEmcGeaTowerTrack[l_found].edep[k];

                          /* Total energy in this tower */
                        }
                      
                    } /* trkno[k] > 0 */
                  
                } /* for Loop over k within a specific GeaTowerTrack record */
              
              
            }  /* l_found > 0, found tower */

        } /* for Loop over j, non-zero towers in a specific cluster */

      /*  You are done with all tower contributors of a cluster */
      /* Sort them decreasing */

      for( j = 0; j < i_size1 - 1; j++)
        {
          for( k = j + 1; k < i_size1; k++)

            {
              if(ra_track[0][k-1] < ra_track[0][k])
                {
                  i_work1 = ia_track[k-1];
                  ia_track[k-1] = ia_track[k];
                  ia_track[k] = i_work1;
                  for ( k1 = 0; k1 <= 11; k1++) ra_work[k1] = ra_track[k1][k-1];
                  for ( k1 = 0; k1 <= 11; k1++) ra_track[k1][k-1] 
                                                 = ra_track[k1][k];
                  for ( k1 = 0; k1 <= 11; k1++) ra_track[k1][k] = ra_work[k1];
                }  /* End interchange of elements */
            }      /* End sort inner loop */
        }          /* End sort outer loop */

      /* Search pid and ptot for the top 3 tracks */

      for( j = 0; j <= 2; j++)
        {
          l_found2 = -1;
          k = 0;
          
          while( k < dEmcGeaTrack_h->nok && l_found2 < 0 && ia_track[j] > 0)
            {
              if(dEmcGeaTrack[k].trkno == ia_track[j])
                {
                  l_found2 = k;
                }
              k = k + 1;
              
            } /* for Loop over k, entries in dEmcGeaTrack */

          if( l_found2 >= 0)
            {
              ra_track[1][j] = dEmcGeaTrack[l_found2].pid;

              /*              ra_track[2][j] = dEmcGeaTrack[l_found2].ekin; */
              ra_track[2][j] = dEmcGeaTrack[l_found2].ptot;
              ra_track[3][j] = dEmcGeaTrack[l_found2].impxyz[0];
              ra_track[4][j] = dEmcGeaTrack[l_found2].impxyz[1];
              ra_track[5][j] = dEmcGeaTrack[l_found2].impxyz[2];
              ra_track[6][j] = dEmcGeaTrack[l_found2].twrhit;
              ra_track[7][j] = dEmcGeaTrack[l_found2].edep;
              ra_track[8][j] = dEmcGeaTrack[l_found2].anclvl;
              ra_track[9][j] = dEmcGeaTrack[l_found2].xyz[0];
              ra_track[10][j] = dEmcGeaTrack[l_found2].xyz[1];
              ra_track[11][j] = dEmcGeaTrack[l_found2].xyz[2];
              
              
            }
          
        }    /* for Loop over j to search for pid at ptot of top 3 tracks */
      
      
      /* Write output record for this cluster */
      
      dEmcGeaClusterTrack[i_cltr_id].id = i_cltr_id;
      dEmcGeaClusterTrack[i_cltr_id].clusid = i;
      dEmcGeaClusterTrack[i_cltr_id].input = 0;
      dEmcGeaClusterTrack[i_cltr_id].type =
        dEmcClusterExt[i].type;
      dEmcGeaClusterTrack[i_cltr_id].arm =
        dEmcClusterExt[i].arm;
      dEmcGeaClusterTrack[i_cltr_id].sector =
        dEmcClusterExt[i].sector;
      
      
      dEmcGeaClusterTrack[i_cltr_id].charged =
            dEmcClusterExt[i].charged;
      for ( k=0; k<3; k++)
        {
          dEmcGeaClusterTrack[i_cltr_id].pc3proj[k] =
            dEmcClusterExt[i].pc3proj[k];
          
        }
      dEmcGeaClusterTrack[i_cltr_id].chi2_sh = 
            dEmcClusterExt[i].chi2_sh;
      /*
      dEmcGeaClusterTrack[i_cltr_id].prob_photon_sh = 
            dEmcClusterExt[i].prob_photon_sh;
      */
      dEmcGeaClusterTrack[i_cltr_id].prob_photon_sh = 
            dEmcClusterExt[i].prob_photon;
      for ( k=0; k<2; k++)
        {
                dEmcGeaClusterTrack[i_cltr_id].e_sh[k] = 
                  dEmcClusterExt[i].e_sh[k];
        }

      for( k = 0; k<= 2; k++)
        {
          dEmcGeaClusterTrack[i_cltr_id].trkno[k] = ia_track[k];
          /* Get the sum of ......
           */

          i_track = ia_track[k];
          j = 0;
          r_edep = 0.0;
          
          while(j<16 && dEmcClusterExt[i].partesum[j] > 0.0)
            {
              i_twrkey = dEmcClusterExt[i].twrlist[j];
              for( j1 = 0; j1 < 16; j1++ )
                
                {
                  if(ia_tower[3][j1] == i_twrkey)
                    {
                      for ( j2 = 0; j2 < 3; j2++)
                        {
                          if (ia_tower[j2][j1] == i_track)
                            {
                              r_edep = r_edep +
                                ra_tower[j2][j1] * ra_tower[3][j1] 
                                / ra_tower[4][j1];
                            }
                        }
                    }
                }
              j++;
            }
          dEmcGeaClusterTrack[i_cltr_id].edep[k] = r_edep;

          /*
          dEmcGeaClusterTrack[i_cltr_id].edep[k] = ra_track[0][k];
          */

          dEmcGeaClusterTrack[i_cltr_id].pid[k] = ra_track[1][k];

          
          dEmcGeaClusterTrack[i_cltr_id].ptot[k] = ra_track[2][k];
          dEmcGeaClusterTrack[i_cltr_id].tracktwrhit[k] = ra_track[6][k];
          dEmcGeaClusterTrack[i_cltr_id].edep_nom[k] = ra_track[7][k];
          dEmcGeaClusterTrack[i_cltr_id].ancestry[k] = ra_track[8][k];
          for( j = 0; j <= 2; j++)
            {
              dEmcGeaClusterTrack[i_cltr_id].xyz[j][k] = ra_track[3+j][k];
              dEmcGeaClusterTrack[i_cltr_id].vertex[j][k] = ra_track[9+j][k];
            }

          dEmcGeaClusterTrack[i_cltr_id].measxyz[k] = 
            dEmcClusterExt[i].xyz[k];
        }

      dEmcGeaClusterTrack[i_cltr_id].mease = 
            dEmcClusterExt[i].e;
      if(dEmcClusterExt[i].e > 0.0)
        {
          for( k =0; k <= 2; k++)
            {
              dEmcGeaClusterTrack[i_cltr_id].efrac[k] =
                dEmcGeaClusterTrack[i_cltr_id].edep[k] /
                dEmcGeaClusterTrack[i_cltr_id].mease;
            }
        }
      dEmcGeaClusterTrack[i_cltr_id].tof = dEmcClusterExt[i].tof;
      dEmcGeaClusterTrack[i_cltr_id].tofmin = dEmcClusterExt[i].tofmin;
      dEmcGeaClusterTrack[i_cltr_id].tofmax = dEmcClusterExt[i].tofmax;
      dEmcGeaClusterTrack[i_cltr_id].etof = dEmcClusterExt[i].ecent;
      dEmcGeaClusterTrack[i_cltr_id].twrhit = dEmcClusterExt[i].twrhit;
      for( k=0; k<2; k++)
        {
          dEmcGeaClusterTrack[i_cltr_id].disp[k] = 
            dEmcClusterExt[i].disp[k];
          dEmcGeaClusterTrack[i_cltr_id].padisp[k] = 
            dEmcClusterExt[i].padisp[k];
        }
      for( k=0; k<8; k++)
        {
          dEmcGeaClusterTrack[i_cltr_id].partesum[k] = 
            dEmcClusterExt[i].partesum[k];
          dEmcGeaClusterTrack[i_cltr_id].chglist[k] = 
            dEmcClusterExt[i].chglist[k];
        }
      

      /* Increment STAF row counter */
      i_cltr_id = i_cltr_id + 1;

      /* Silly fix to avoid all 0 outputs - something may be fishy
         if Lead Glass sector 1 ??? */

      /* This correction seems to have gone wrong !!!
      if (dEmcGeaClusterTrack[i_cltr_id - 1].pid[k] == 0)
        {
          i_cltr_id = i_cltr_id - 1;
        }

      */
      
    }    /* Loop over i, all clusters recognized in dEmcClusterExt */
  
  dEmcGeaClusterTrack_h->nok = i_cltr_id;
  

  /* -------------------------------------------------------------   */

  /* Second part.  Loop over all tracks that deposited anything in the 
     calorimeter.  Point to clusters where a specific track contributed */

  i_trcl_id = 0;

  for ( i = 0; i < dEmcGeaTrack_h->nok; i++)
    {
      
      l_found = -1;   /* Found at least one cluster where it contributed */
      for (ll = 0; ll < sizeof(ra_track)/sizeof(ra_track[0][0]); ll++)
        {
          *((float *)ra_track + ll) = 0.0;
        }

      /* Since 0 is a valid pointer now, "not found" is indicated
         by -1 */
      for(j = 0; j <= 2; j++) ia_track[j] = -1;
      i_track = dEmcGeaTrack[i].trkno;
      nom_edep = dEmcGeaTrack[i].edep;
      
      for ( j = 0; j < dEmcGeaClusterTrack_h->nok; j++)
        {
          for ( k = 0; k <= 2; k++)
            {
              if ( dEmcGeaClusterTrack[j].trkno[k] == i_track
                   && l_found < i_size1 - 1 )
                {
                  l_found = l_found + 1;
                  ia_track[l_found] = dEmcGeaClusterTrack[j].clusid;
                  ra_track[0][l_found] = dEmcGeaClusterTrack[j].pid[k];
                  ra_track[1][l_found] = dEmcGeaClusterTrack[j].ptot[k];
                  ra_track[2][l_found] = dEmcGeaClusterTrack[j].edep[k];
                  /*
                  ra_track[3][l_found] = dEmcGeaClusterTrack[j].efrac[k];
                  */
                  if(nom_edep > 0.1)
                    {
                      
                      ra_track[3][l_found] = 
                        dEmcGeaClusterTrack[j].edep[k] / nom_edep;
                    }
                  else
                    {
                      ra_track[3][l_found] = 0.0;
                    }
                }  /* New cluster found where it contributed */
              
            }  /* Loop over k = 0, 2, within one dEmcGeaClusterTrack record */
          
        }   /* Loop over j = 0,dEmcGeaClusterTrack_h.nok   */

      /* Sort descending the clusters where the track contributed */

      if (l_found > - 1 )
        {
          for (j = 0; j < l_found; j++)
            {
              for ( k = 1; k <= l_found; k++) 
                {
                  if (ra_track[3][k-1] < ra_track[3][k])
                    {
                      i_work1 = ia_track[k-1];
                      for ( k1 = 0; k1 <= 11; k1++) ra_work[k1] = 
                                                     ra_track[k1][k-1];
                      for ( k1 = 0; k1 <= 11; k1++) ra_track[k1][k-1] 
                                                     = ra_track[k1][k];
                      for ( k1 = 0; k1 <= 11; k1++) ra_track[k1][k] = 
                                                     ra_work[k1];
                      ia_track[k-1] = ia_track[k];
                      ia_track[k] = i_work1;
                    }
                }  /* End loop over k, inner sort loop */
            }      /* End loop over j, outer sort loop */
      
          /* Write out a new record of dEmcGeaTrackCluster */
          
          dEmcGeaTrackCluster[i_trcl_id].id = i_trcl_id + 1;
          dEmcGeaTrackCluster[i_trcl_id].trkno = i_track;
          dEmcGeaTrackCluster[i_trcl_id].track_ptr = i;
          dEmcGeaTrackCluster[i_trcl_id].nom_edep = nom_edep;
          dEmcGeaTrackCluster[i_trcl_id].input = 0;
          for (k = 0; k <= 2; k++)
            {
              dEmcGeaTrackCluster[i_trcl_id].clusid[k] = ia_track[k];
              dEmcGeaTrackCluster[i_trcl_id].pid = ra_track[0][0];
              dEmcGeaTrackCluster[i_trcl_id].ptot = ra_track[1][0];
              dEmcGeaTrackCluster[i_trcl_id].edep[k] = ra_track[2][k];
              dEmcGeaTrackCluster[i_trcl_id].efrac[k] = ra_track[3][k];
            }

          /*
          printf("true = %f %f \n",dEmcGeaTrackCluster[i_trcl_id].ptot,
                 dEmcGeaTrackCluster[i_trcl_id].edep[0]);
          */      
          i_trcl_id = i_trcl_id + 1;
          

        }   /* End if(l_found) - found at least one cluster where contrib. */
      

      
    }   /*  Loop over i = 0,dEmcGeaTrack_h.nok , all tracks depositing E */
  

  dEmcGeaTrackCluster_h->nok = i_trcl_id;
  

    return STAFCV_OK;
}