mEmcGeaClusterEval.c

Go to the documentation of this file.

/*:>--------------------------------------------------------------------
**: FILE:       mEmcGeaClusterEval.c.template
**: HISTORY:
**:             00jan93-v000a-hpl- Created by stic Version
**:  Id: idl.y,v 1.3 1998/04/20 19:56:35 leitch Exp  
**:<------------------------------------------------------------------*/
#include "mEmcGeaClusterEval.h"
#include "emlLib.h"

long mEmcGeaClusterEval_(
  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 *dEmcClusterLocalExt_h, DEMCCLUSTERLOCALEXT_ST *dEmcClusterLocalExt ,
  TABLE_HEAD_ST *dEmcGeaTrackCluster_h, DEMCGEATRACKCLUSTER_ST *dEmcGeaTrackCluster ,
  TABLE_HEAD_ST *dEmcGeaClusterTrack_h, DEMCGEACLUSTERTRACK_ST *dEmcGeaClusterTrack )
{
/*:>--------------------------------------------------------------------
**: ROUTINE:    mEmcGeaClusterEval_
**: DESCRIPTION: Physics Analysis Module ANSI C template.
**:             This is an ANSI C Physics Analysis Module template
**:             automatically generated by stic from mEmcGeaClusterEval.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
**:  dEmcClusterLocalExt    - PLEASE FILL IN DESCRIPTION HERE
**:  dEmcClusterLocalExt_h   - header Structure for dEmcClusterLocalExt
**:    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,k,k1;
  int l_found,l_found1,l_found2;
  int i_track;
  int ll;   
  int i_twrkey;
  float ra_work[12];
  int   i_work1;
  float nom_edep;

  // Global boolean if one contributor is ever found in the cluster
  // needed for the real+simul merging !
  int AtLeastOneContributorFound = 0 ; 

  int maxNumOfTracksPerTower = 50;

  // We consider a realistic maximum of *3* tracks contributing to the same tower 
  // The most common case is, of course, 1 ...
  int maxNumOfTracksKeptPerTower = 3;  

  /* ia_track will store the track numbers of the tracks contributing
     to a specific cluster 
  */
  int   ia_track[maxNumOfTracksPerTower];
  
  /* ra_track will store data of tracks listed in ia_track.
     Field: 0 - energy (cumulates "edep" field of TowerTrack)
            1 - pid
            2 - ptot
            3-5 - impact xyz
            6 - twrhit
            7 - edep ("edep" field of GeaTrack)
            8 - anclvl
            9-11 - vertex xyz
  */
  float ra_track[12][maxNumOfTracksPerTower];
  
  if(dEmcGeaTrack_h->nok <=0) return (STAFCV_BAD);
  if(dEmcGeaTowerTrack_h->nok <=0) return (STAFCV_BAD);
  if(dEmcClusterLocalExt_h->nok <=0) return (STAFCV_BAD);
  
  i_cltr_id = 0;

  /* 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 < dEmcClusterLocalExt_h->nok; i++)
    {

      // reset "AtLeastOneContributorFound" for each cluster
      AtLeastOneContributorFound = 0 ; 

      /* Zeroing temporary arrays */

      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;
        }

      /* Loop over all (16) towers included in a cluster */

      for( j = 0; j < 16 && dEmcClusterLocalExt[i].partesum[j] > 0.0; j++)
        {
          i_twrkey = dEmcClusterLocalExt[i].twrlist[j];
          
          /* Loop over all tracks contributing to a given tower */

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

          if( l_found >= 0)  /* Contributor found */
            {

              /* We consider a maximum of *3* tracks contributing to the same tower */

              for( k = 0; k <= maxNumOfTracksKeptPerTower-1; k++)
                {
                  i_track = dEmcGeaTowerTrack[l_found].trkno[k];
                  if( i_track > 0)
                    {
                      l_found1 = -1;
                      k1 = 0;
                      
                      while( l_found1 < 0 && k1 < maxNumOfTracksPerTower)
                        {
                          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];
                        }

                    } /* 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 */

      if( AtLeastOneContributorFound )      
        {
          /* Sort them decreasing */
          for( j = 0; j < maxNumOfTracksPerTower - 1; j++)
            {
              /*          for( k = j + 1; k < maxNumOfTracksPerTower; k++) */
              for( k = 1; k < maxNumOfTracksPerTower; k++)
                {
                  if(ra_track[0][k-1] < ra_track[0][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 interchange of elements */
                }      /* End sort inner loop */
            }          /* End sort outer loop */
          
          /* Search pid and ptot for the top 3 tracks */
          
          for( j = 0; j <= maxNumOfTracksKeptPerTower-1; 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 */

        }  /* Contributor found */

      /* Write output record for this cluster */

      // In the case of embedded sim+real data, you often DO NOT find GEANT 
      // contributors to a reconstructed cluster. At this point we have, 
      // thus, 2 options:
      //
      // 1) Do not apply the AtLeastOneContributorFound condition:
      //    Fill dEmcGeaClusterTrack in all cases (though the *pure* dEmcGeaClusterTrack 
      //    fields below would be empty in most of the cases ...).
      //    In the analysis afterwards, the way to distinguish a "pure" ClusterTrack 
      //    (with some GEANT track in it) and a ClusterLocal without contribution
      //    will be, e.g., asking that the first track not be there (trkno_1 ==0).
      
      // 2) Apply the AtLeastOneContributorFound condition: Fill only dEmcGeaClusterTrack
      //    when a GEANT contributor is found. The resulting dEmcGeaClusterTrack
      //    tables are much less heavy, but we have to propagate the (real) event
      //    cluster multiplicity information in an adhoc manner ...
      //
      // We choose now option 2) and propagate the clus. mult. info in the 
      // unused "Charged" field of dEmcGeaClusterTrack D.d'E. (2002)
      // 
      
      if( AtLeastOneContributorFound )      
        {
          dEmcGeaClusterTrack[i_cltr_id].id = i_cltr_id;
          dEmcGeaClusterTrack[i_cltr_id].clusid = i;
          dEmcGeaClusterTrack[i_cltr_id].input = 0;

          // Fields common to dEmcClusterLocalExt
          dEmcGeaClusterTrack[i_cltr_id].evno = dEmcClusterLocalExt[i].evno;
          dEmcGeaClusterTrack[i_cltr_id].keycent = dEmcClusterLocalExt[i].twrlist[0];
          dEmcGeaClusterTrack[i_cltr_id].type = dEmcClusterLocalExt[i].type;
          dEmcGeaClusterTrack[i_cltr_id].arm  = dEmcClusterLocalExt[i].arm;
          dEmcGeaClusterTrack[i_cltr_id].sector = dEmcClusterLocalExt[i].sector;  
          dEmcGeaClusterTrack[i_cltr_id].mease = dEmcClusterLocalExt[i].e;
          
          // The ecore field of dEmcGeaClusterTrack stores now ecore for PbSc and ecorr for PbGl !
          if (dEmcGeaClusterTrack[i_cltr_id].type == 1 ) // PbSc
            dEmcGeaClusterTrack[i_cltr_id].ecore = dEmcClusterLocalExt[i].ecore;
          if (dEmcGeaClusterTrack[i_cltr_id].type == 2)  // PbGl
            dEmcGeaClusterTrack[i_cltr_id].ecore = dEmcClusterLocalExt[i].ecorr;
          
          dEmcGeaClusterTrack[i_cltr_id].tof = dEmcClusterLocalExt[i].tof;
          dEmcGeaClusterTrack[i_cltr_id].tofmin = dEmcClusterLocalExt[i].tofmin;
          dEmcGeaClusterTrack[i_cltr_id].tofmax = dEmcClusterLocalExt[i].tofmax;
          // the energy of the first-TOF tower is the energy of the highest tower:
          dEmcGeaClusterTrack[i_cltr_id].etof = dEmcClusterLocalExt[i].ecent;
          dEmcGeaClusterTrack[i_cltr_id].etofmin = dEmcClusterLocalExt[i].etofmin;
          dEmcGeaClusterTrack[i_cltr_id].etofmax = dEmcClusterLocalExt[i].etofmax;
          dEmcGeaClusterTrack[i_cltr_id].twrhit = dEmcClusterLocalExt[i].twrhit;
          
          for ( k=0 ; k<= 2 ; k++)
            {
              dEmcGeaClusterTrack[i_cltr_id].e_sh[k] = dEmcClusterLocalExt[i].e_sh[k];
              dEmcGeaClusterTrack[i_cltr_id].disp[k] = dEmcClusterLocalExt[i].disp[k];
              dEmcGeaClusterTrack[i_cltr_id].padisp[k] = dEmcClusterLocalExt[i].padisp[k];
              dEmcGeaClusterTrack[i_cltr_id].measxyz[k] = dEmcClusterLocalExt[i].xyz[k];
            } 
          for( k=0; k<=8; k++)
            {
              dEmcGeaClusterTrack[i_cltr_id].partesum[k] = dEmcClusterLocalExt[i].partesum[k];
            }
          
          // Fields common to dEmcClusterLocalExt but somehow different ...
          dEmcGeaClusterTrack[i_cltr_id].chi2_sh = dEmcClusterLocalExt[i].chi2 ;
          dEmcGeaClusterTrack[i_cltr_id].prob_photon_sh = dEmcClusterLocalExt[i].prob_photon ;
          
          // "Pure" dEmcGeaClusterTrack fields
          
          for( k = 0 ; k<= maxNumOfTracksKeptPerTower-1 ; k++)
            {
              dEmcGeaClusterTrack[i_cltr_id].trkno[k] = ia_track[k];
              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];
                }
            }
          
          if(dEmcGeaClusterTrack[i_cltr_id].mease > 0.0)
            {
              for( k =0; k <= maxNumOfTracksKeptPerTower-1 ; k++)
                {
                  dEmcGeaClusterTrack[i_cltr_id].efrac[k] =
                    dEmcGeaClusterTrack[i_cltr_id].edep[k] /
                    dEmcGeaClusterTrack[i_cltr_id].mease;
                }
            }
          
          // "Pure" dEmcGeaClusterTrack fields not used

          // WATCH-OUT: we (ab)use of this field to store the real evt. cluster info
          // dEmcGeaClusterTrack[i_cltr_id].charged = 0;
          dEmcGeaClusterTrack[i_cltr_id].charged = (long)dEmcClusterLocalExt_h->nok;

          for( k=0; k<8; k++)
            {
              dEmcGeaClusterTrack[i_cltr_id].chglist[k] = 0.;
            }

          // WATCH-OUT: we (ab)use of this field to store the dead/warn map info
          //for( k=0; k<3; k++)
          //  {
          //   dEmcGeaClusterTrack[i_cltr_id].pc3proj[k] = 0.;
          //  }

          dEmcGeaClusterTrack[i_cltr_id].pc3proj[0] = (float)dEmcClusterLocalExt[i].deadmap; // (21+9 bits int stored in float)
          dEmcGeaClusterTrack[i_cltr_id].pc3proj[1] = (float)dEmcClusterLocalExt[i].warnmap; // (21+9 bits int stored in float)


          dEmcGeaClusterTrack[i_cltr_id].pc3proj[2] = 0;

          /* Increment STAF row counter */
          i_cltr_id = i_cltr_id + 1;
      
        } /* End of  if( AtLeastOneContributorFound ) , i.e. you found at least one GEANT
             contributor to this cluster */
      
    }    /* Loop over i, all clusters recognized in dEmcClusterLocalExt */
  
  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 */
      AtLeastOneContributorFound = 0 ; 

      // Zeroing arrays
      for (ll = 0; ll < sizeof(ra_track)/sizeof(ra_track[0][0]); ll++)
        {
          *((float *)ra_track + ll) = 0;
        }

      /* Since 0 is a valid pointer now, "not found" is indicated
         by -1 */
      for(j = 0; j <= maxNumOfTracksKeptPerTower-1 ; 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 <= maxNumOfTracksKeptPerTower-1 ; k++)
            {
              if ( dEmcGeaClusterTrack[j].trkno[k] == i_track
                   && l_found < maxNumOfTracksPerTower - 1 )
                {
                  l_found = l_found + 1; 
                  AtLeastOneContributorFound = 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 ( AtLeastOneContributorFound )
        {
          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 <= maxNumOfTracksKeptPerTower-1 ; 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];
            }

          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;
}