Utility routines

File Rec_Util.f

Subroutine Ilin_Fit_3d (rs,xs,ys,zs,cx,cy,cz,nnpt,ind,rpt,xpt,ypt,zpt,wpt,icos,chisq,ierror)

User routine to find parameters for 3-D line through the array of 3-D points with energy as a weight - cluster axes.

There are several models can be applied to find the main cluster axes. One of them used 3-D fit to find the line with the smallest $\chi^2$ of distances from that line to all cluster points. Another one used something similar to the trust definition in PYTHIA. (The programs are available optionally) Here it was chosen the tensor of inertia model of cluster as solid body. Actually program finds the main tensor axis.

********************************************************************
*     Called old CERNLIB routines to find matrix eigenvalues       *
********************************************************************
c  rs,xs,ys,zs -- weighted cluster center INPUT  = OUTPUT 
c                                - radius and coordinates
c  cx,cy,cz    -- direct cosines  INPUT  = zero approximation
c                                 OUTPUT = result of fit
c   nnpt        -- number of points in fit array
c   ind         -- number of sampling coordinate
c   rpt,xpt,ypt,zpt -- arraies of points to fit - radius and coordinates
c   wpt         -- array of weights of poits
c   icos        -- number of axis the closest to cosines zero approximation
c   chisq       -- chi-sqare OUTPUT
c   ierror      -- .NOT. zero if something wrong during the fit procedure
c----------------------------------------------------------------------

---------------

SUBROUTINE GKEY_PACK_HBAR

(IMOD,IX,IY,IZ,KEYG)

c-----------------------------------------------------------------------
c              Valid for Barrel part of Hadron Calorimeter 
c-----------------------------------------------------------------------
c  Input -->
c         IMOD  --> module number 1 -- 32
c         IX    --> X-coordinate cell number in module coordinate system
c         IY    --> Y-coordinate cell number in module coordinate system
c         IZ    --> Z-coordinate cell number in module coordinate system
C  Output  -->  KEYG -- PACKED GEOMETRICAL KEY FOR CALORIMETER HITS
c-----------------------------------------------------------------------

---------------

SUBROUTINE GKEY_UNPACK_HBAR

(KEYG,IMOD,IX,IY,IZ)

c-----------------------------------------------------------------------
c              Valid for Barrel part of Hadron Calorimeter 
c-----------------------------------------------------------------------
C  Input  -->  KEYG -- PACKED GEOMETRICAL KEY FOR CALORIMETER HITS
c
c  Output -->
c         IMOD  --> module number 1 -- 32
c         IX    --> X-coordinate cell number in module coordinate system
c         IY    --> Y-coordinate cell number in module coordinate system
c         IZ    --> Z-coordinate cell number in module coordinate system
c-----------------------------------------------------------------------

---------------

SUBROUTINE GKEY_PACK_HCAP

(K_quad,K_xreg,K_yreg,K_zreg,IX,IY,IZ,KEYG)

c-----------------------------------------------------------------------
c              Valid for EndCap part of Hadron Calorimeter 
c-----------------------------------------------------------------------
c  Input -->
C         K\_quad--> Quadrant number
C         K\_xreg--> Number of X region of quadrant
C         K\_yreg--> Number of Y region of quadrant 
C         K\_zreg--> Number of Z region of quadrant 
c         IX    --> X-coordinate cell number in module coordinate system
c         IY    --> Y-coordinate cell number in module coordinate system
c         IZ    --> Z-coordinate cell number in module coordinate system
C  Output  -->  KEYG -- PACKED GEOMETRICAL KEY FOR CALORIMETER HITS
c-----------------------------------------------------------------------

---------------

SUBROUTINE GKEY_UNPACK_HCAP

(KEYG,K_quad,K_xreg,K_yreg,K_zreg,IX,IY,IZ)

c-----------------------------------------------------------------------
c              Valid for EndCap part of Hadron Calorimeter 
c-----------------------------------------------------------------------
C  Input  -->  KEYG -- PACKED GEOMETRICAL KEY FOR CALORIMETER HITS
c
c  Output -->
C         K\_quad--> Quadrant number
C         K\_xreg--> Number of X region of quadrant
C         K\_yreg--> Number of Y region of quadrant 
C         K\_zreg--> Number of Z region of quadrant 
c         IX    --> X-coordinate cell number in module coordinate system
c         IY    --> Y-coordinate cell number in module coordinate system
c         IZ    --> Z-coordinate cell number in module coordinate system
c-----------------------------------------------------------------------

---------------

Subroutine Mom_Helix_Dist (xp,yp,zp,cx,cy,cz,p,qq,xhi,yhi,zhi,dist)

3-D distance from hit to helix track model with momentum p(cx,cy,cz) at point r(xp,yp,zp)

C.    ******************************************************************
C.    *  input: xp,yp,zp particle 3-d position  [cm]                   *
C.    *         cx,cy,cz momentum direction cosines                    *
C.    *         qq,p  - charge and momentum  [GeV/c]                   *
C.    *         xhi,yhi,zhi - HIT  coordinates                         *
C.    *  output:  dist - 3-d distance from helix to hit-point          *
C.    *          in field Bz30 = 3 Tesla                               *
C.    *          in field Bz40 = 4 Tesla                               *
C.    ******************************************************************

Method:

Let us write some the simplest definitions:
the distance along $Z$ is $z = \Theta b$ and distance along 2-dim circle is $a = \Theta R$ from these two equation and taken into account that the same point has the same time to reach it one can derive in some approximation:
$t = D \Theta R m/p_{xy} = D \Theta b m/p_z$; where: $m$ is particle mass:
then we have $b = R (p_z/p_{xy})$

Some definitions:

$$\mathit{p_{xy}} = \sqrt{\mathit{p_x}^{2} + \mathit{p_y}^{2}}... ...c {\mathit{p}_{z}}{ \left\vert \! \,q\,B\, \! \right\vert }} ;$$

$$\mathit{x}_{c} = \mathit{x}_{p} + {\displaystyle \frac {\mat... ...\mathit{y}_{p} - {\displaystyle \frac {\mathit{p}_{x}}{q\,B}}$$

$$\mathit{R^2+b^2} = {\displaystyle \frac {\mathit{px}^{2} + ... ...thit{pz}^{2}}{ \left\vert \! \,q\,B\, \! \right\vert ^{2}}}$$

Helix equation with axes along Z and with zero $\Theta$ at the point $x_p, y_p, z_p$ exactly:

$$\left\{ \begin{array}{l} x = \cos{(q\; \Theta)}(x_p-x_c) - \... ...p-x_c) + y_c; \\ z = b\; \Theta + z_p; \\ \end{array}\right.$$

$$dx = \frac{\partial x}{\partial \Theta}; \; \; dy = \frac{\p... ... \Theta}; \; \; dz = \frac{\partial z}{\partial \Theta}; \; \;$$

Equation of normal plane to helix through the hit point:

$$(x_h-x)dx/(R^2+b^2) + (y_h-y)dy/(R^2+b^2) + (z_h-z)dz/(R^2+b^2)=0$$

After substitutions:

$$\begin{array}{l} - p_y\; q\; B x_h \sin(q\; \Theta)\; \vert ... ...\Theta + p_z\; q\; B^2 z_p\; \vert q\; B\vert = 0 \end{array}$$

The second order approximation (Tailor raw) of equation of normal plane to helix passes through the hit point near by $\Theta=0$ is

$$\begin{array}{l} \left( {\displaystyle \frac {( - {\displays... ..._x}^{2} + \mathit{p_y} ^{2} + \mathit{p_z}^{2}}} =0 \end{array}$$

The solution of this equation gives the 3-D nearest to the hit point at the helix with rather good accuracy.

---------------

Subroutine Mom_Helix_Dist_Array (xp,yp,zp,cx,cy,cz,p,qq,nh,xhi,yhi,zhi,dist)

The same as Mom_Helix_Dist but for array of hits.

---------------

Subroutine Helix_dist

(xs,ys,zs,xh,yh,zh,x0,y0,rad,b,phi0,dphi,xm,ym,zm,dist)

C.
C.    ******************************************************************
C.    *                                                                *
C.    *     Distance from arbitrary point to  Helix                    *
C.    *    Helix path through the center of coordinate system          *
C.    *                                                                *
C.    *  input: xs,ys,zs - cluster center coordinates                  *
C.    *         cx,cy,cz - direction cosines from cluster fit          *
C.    *         xh,yh,zh - arbitrary poit coordinates                  *
C.    *         x0, y0, rad - circle helix parameter                   *
C.    *         b helix slope                                          *
C.    *         phi0 angle at the helix at 0,0,0 point                 *
C.    *         dphi angle between 0,0,0 and center of cluster         *
C.    *               along circle                                     *
C.    *                                                                *
C.    *   output: xm,ym,zm - point at the helix closests to hit        *
C.    *           dist - distance from helix to hit                    *
C.    *                                                                *
C.    ******************************************************************

---------------

SUBROUTINE Mom_Helix

(xp,yp,zp,cx,cy,cz,p,q,x0,y0,d0,rad,b,phi0,dphi,z0)

C.
C.    ******************************************************************
C.    *                                                                *
C.    *           Helix track model calculation                        *
C.    *  for momentum p(cx,cy,cz) at point r(xp,yp,zp)                 *
C.    *          in field Bz30 = 3 Tesla                               *
C.    *          in field Bz40 = 4 Tesla                               *
C.    *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*
C.    *      This program has a mistake in the d0 definition           *
C.    *       it was used one more approach d0 lies at Z=0             *
C.    *!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*
C.    *                                                                *
C.    *                                                                *
C.    *  input: xp,yp,zp particle 3-d position  [cm]                   *
C.    *         cx,cy,cz momentum direction cosines                    *
C.    *         q,p  - charge and momentum  [GeV/c]                    *
C.    *  output: x0, y0, rad - circle helix parameter                  *
C.    *          d0 - P.C.A. distance                                  *
C.    *          b helix slope                                         *
C.    *          phi0 angle at the helix circle at P.C.A.              *
C.    *          dphi angle between P.C.A. and xp,yp,zp                *
C.    *               along circle                                     *
C.    *                                                                *
C.    ******************************************************************

---------------

SUBROUTINE Clust_Helix

(xs,ys,zs,CX,CY,CZ,x0,y0,rad,b,phi0,dphi,pxy,p)

C.
C.    ******************************************************************
C.    *                                                                *
C.    *           Helix track model calculation                        *
C.    *       goes through the xs,ys,zs and cx,cy,cz  and 0,0,0        *
C.    *                                                                *
C.    *       Author : V.L. Morgunov                                   *
C.    *       created  25-Oct-2000  see helix.mws                      *
C.    *                                                                *
C.    *    Helix path through the center of coordinate system          *
C.    *                                                                *
C.    *  input: xs,ys,zs cluster center coordinates                    *
C.    *         cx,cy,cz direction cosines from cluster fit            *
C.    *  output: x0, y0, rad - circle helix parameter                  *
C.    *          b helix slope                                         *
C.    *          phi0 angle at the helix circle at 0,0,0 point         *
C.    *          dphi angle between 0,0,0 and center of cluster        *
C.    *               along circle                                     *
C.    *          pxy - transversal to field momentum                   *
C.    *          p   - derived momentum                                *
C.    *                                                                *
C.    ******************************************************************

---------------

SUBROUTINE helix_geant_draw

(charge,xp,yp,zp,cx,cy,cz,p,icol)

C.    ******************************************************************
C.    *   Calling sequence CALL GHELX3(FIELDM*CHARGE,STEP,VECT,VOUT)   *
C.    *                                                                *
C.    *  Tracking routine in a constant field oriented along axis 3    *
C.    *  Tracking is performed with a conventional helix step method   *
C.    *   input                                                        *
C.    *     STEP = arc length of the step asked                        *
C.    *     VECT = input vector (position,direction cos and momentum)  *
C.    *     CHARGE=  electric charge of the particle                   *
C.    *   output                                                       *
C.    *     VOUT = same as VECT after completion of the step           *
C.    *  made from   GHELX3                                            *
c.    *  Authors    R.Brun, M.Hansroul,  Rewriten  V.Perevoztchikov    *
C.    ******************************************************************
C.    ------------------------------------------------------------------
C       units are kgauss,centimeters,gev/c
C.    ------------------------------------------------------------------

---------------

SUBROUTINE Mhelix_draw

(x0,y0,r,zs,phi0,dphi,icol)

C.
C.    ******************************************************************
C.    *                                                                *
C.    *       Draw a helix track model                                 *
C.    *       goes through the xs,ys,zs and cx,cy,cz                   *
C.    *                                                                *
C.    *       Rewriten from   Author : P.Zanarini                      *
c     *   Input parameters:                                            *
c     *         x0,y0,r  Helix center coordinates and its radius       *
c     *         Zs    -- shift from zero point to cluster center       *
c     *         Phi   -- initial angle to draw                         *
c     *         dPhi  -- angle to draw                                 *
c     *   All of them are calculated by Clust\_Helix subroutine        *
C.    ******************************************************************
C.

---------------

SUBROUTINE Mcircle

(R,x0,y0,z0,icol)

C.
C.    ******************************************************************
C.    *       Draw a circle of radius R centered on x0,y0,z0           *
C.    *      Rewriten from (GEANT 321) Author : P.Zanarini             *
C.    ******************************************************************
C.

---------------

SUBROUTINE draw_cross

(X,Y,Z,SIZSYM,icol)

C.    ******************************************************************
C.    *       Draw an HIT point at (X,Y,Z)                             *
c.    *    Was made from GDRHIT (GEANT 321)                            *
C.    ******************************************************************

---------------

SUBROUTINE draw_hit

(X,Y,Z,SIZSYM,icol)

C.    ******************************************************************
C.    *       Draw an HIT point at (X,Y,Z)                             *
c.    *    Was made from GDRHIT (GEANT 321)                            *
C.    ******************************************************************

---------------

subroutine Point_proj

(key,x0,y0,z0,anx,any,anz,n,x,y,z,xx,yy)

c---------------------------------------------------------
c   Project of set of n points {x,y,z} to the reference
c  plate and convert its coordinates to the frame at that plate
c---------------------------------------------------------
c   INPUT 
c          key = 0 -- calculate transfer matrix
c          key = 1 -- transfer points to the reference plate
c          x0, y0, z0 -- point at reference plate
c          anx,any,anz -- direct vecrtor of reference plate
c          n - array of x,y,z point to convert
c   OUTPUT
c         -- xx,yy,zz resulting array of new coordinates
c---------------------------------------------------------

---------------

subroutine PAIR_TO_GROUP

(npair,ia,ib,ngr,mgrup,kgrup,lgrup)

C.******************************************************************
c   INPUT
C     Npair       -- number of track pairs
c     Ia, Ib      -- track addresses lists of close-by pairs
c-------------------------------------------------------------------
c   OUTPUT
c     Ngrup       -- number of groups found
c     Mgrup(i)    -- number of members in group
c     Kgrup(i)    -- number of pairs belong to group
c     Lgrup(i,j)  -- Track adresses of group
C.******************************************************************
c   Lgrup has to be exactly equal of the same variable in tpctrack.inc
C.******************************************************************
c     No more then 50 groups and 20 members in any group
c------------------------------------------------------------

---------------

Subroutine Point_line

(xp,yp,zp,cx,cy,cz,xv,yv,zv,xm,ym,zm)

c----------------------------------------------------------------------
c       Point at line that closest to point near by line  
c----------------------------------------------------------------------
c       xp, yp, zp -- point is at the line cx,cy,cz
c       xv, yv, zv -- point that is out of line
c       xm, ym, zm -- closest point at the line
c----------------------------------------------------------------------

---------------

Function sdist

(xp,yp,zp,cx,cy,cz,xv,yv,zv)

c----------------------------------------------------------------------
c        Distance from line to point
c----------------------------------------------------------------------
c       xp, yp, zp -- point is at the line
c       xv, yv, zv -- point is out of line

---------------

subroutine Matrix_sort

(a,n,m,k,i,j)

********************************************************************
*     Utility routine to find K smallest values                    *
*    from N x M values of array A and their addresses              *
*     INPUT  A -- Array with tested numbers                        *
*            N -- first index size in A                            *
*            M -- second indexs size in A                          *
*            K -- number of needed minimal numbers                 *
*     OUTPUT I -- list of first indexes in increasing order of value*
*            J -- list of second indexes in increasing order of val.*
********************************************************************

---------------

subroutine Min_amp

(n,a,m,b,ia)

********************************************************************
*     Utility routine to find M smallest values                    *
*    from N amplitudes array and thier addresses                   *
*          Can be used if M << N                                   *
*     in opposite case better use SORTZV (CERN-lib)                *
*     INPUT                                                        *
*            N -- initial array size                               *
*            A -- Array with tested numbers                        *
*            M -- Number of needed min numbers                     *
*    OUTPUT                                                        *
*            B -- Sorted array                                     *
*            Ia -- addresses of sorted numbers in initial array    *
********************************************************************

---------------

subroutine Max_amp

(n,a,m,b,ia)

********************************************************************
*     Utility routine to find M biggest values                     *
*    from N amplitudes array and thier addresses                   *
*          Can be used if M << N                                   *
*     in opposite case better use SORTZV (CERN-lib)                *
*     INPUT  N -- initial array size                               *
*            A -- Array with tested numbers                        *
*            M -- Number of needed min numbers                     *
*            B -- Sorted array                                     *
*            Ia -- addresses of sorted numbers in initial array    *
********************************************************************

---------------

function angle_dist

(t1,p1,t2,p2)

********************************************************************
*    Calculate angle between two vectors at sphere                 *
* Input                                                            *
*   T(i) and P(i) --  spherical coordinates of vectors [degree]    *
* Output                                                           *
*    angle\_dist -- Angle between vectors  [degree]                *
********************************************************************

---------------

subroutine hit_move_rot

(rcyl,rs,xs,ys,zs,cx,cy,cz,xhi,yhi,zhi,xmr,ymr,zmr)
********************************************************************
*     user routine called from Find\_HB\_clust                     *
*  to find points at the cylinder surface around  Super-Clusters   *
*  in the coordinate system with axes of Super-Cluster as Z-axes   *
c      See MAPLE file cyleq2.mws                                   *
********************************************************************

---------------

subroutine hit_move

(rcyl,rs,xs,ys,zs,cx,cy,cz,xhi,yhi,zhi,xm,ym,zm)

********************************************************************
*     user routine called from Find\_HB\_clust                     *
*  to find points at the cylinder surface around  Super-Clusters   *
c==========================================================
c      See MAPLE file cyleq2.mws
c  Input parameters
c   rcyl -- Radius of cylinder around Super-Cluster 
c   xs,ys,zs -- Center of Super-Cluster  
c   rs,cx,cy,cz -- it's radius and Direct COSines
c   xhi,yhi,zhi -- Hit coordinates
c  Output parameters
c   xm,ym,zm -- Projection of hit to the cylinder surface 
c               around Super-Cluster
c==========================================================
      rr = 1.0/rs
c  Coordinate of point at the Super-Cluster axes 
c     that is perpendicular to Hit point

---------------

subroutine hit_rot

(cx,cy,cz,xm,ym,zm,xmr,ymr,zmr)

********************************************************************
*     user routine called from Find\_HB\_clust                       *
*  to find points in the coordinate system                         *
*     with axes of Super-Cluster as Z-axes                         *
c      See MAPLE file cyleq2.mws
c  Input parameters
c   Rotation matrix to the coordinate system
c     with axes of Super-Cluster as Z-axes
c==========================================================
c  Direction of S-Cl axes become to be +Z' axes
c    and X' axes has a + signum in old coordinate system
c   if is(cy>0.0) then e11:=cy/sqrt(cx^2+cy^2); 
c                else e11:=-cy/sqrt(cx^2+cy^2); fi;
c   if is(cy>0.0) then e12:=-cx/sqrt(cx^2+cy^2); 
c                 else e12:=cx/sqrt(cx^2+cy^2); fi;
c   e13 := 0;
c   e21 := e13*e32-e12*e33;
c   e22 :=  e11*e33-e13*e31;
c   e23 := e12*e31-e11*e32;
c   e31 := cx;
c   e32 := cy;
c   e33 := cz;
c    Rotation matrix is inverse to defined above
c==========================================================

---------------

function coosh

(x,step)

********************************************************************
*  Shift coordinate to the center of cell with width = STEP        *
********************************************************************

---------------

SUBROUTINE FIT_LONG_SHOWER

(n,x,y,a,b,c,chi)

c=========================================================================
c               Utility routine
c       Chi--square method for calculation of
c       Longitudilal Electromagnetic Shower model
c       parameters A,B,C in E(x) = C*(X^A)*exp(-B*X)
c           See MAPLE file shower.mws
c       N    --- number of points
c       X(i) --- distance along shower
c       Y(i) --- amplitudes along shower
c       W(i) --- weight
c       Chi  --- Chi-square in the histogram space but not the logarithmic one
c        it is some model dependent chisqare with distribution around unit 
c=========================================================================

---------------

SUBROUTINE FIT_LINEAR_S

(n,x,y,s,a,b)

c=========================================================================
c               Utility routine
c       Chi--square method for calculation of
c       parameters A,B in linear equation Y(x) = A*X + B
c           See MAPLE file linear.mws
c       N    --- number of points
c       X(i) --- distance along shower
c       Y(i) --- amplitudes along shower
c       S(i) --- error of amplitude
c=========================================================================
*   Here have to be added -- correct zero amplitudes account
*      and   chi-square calculations for fitted curve
c=========================================================================

---------------

SUBROUTINE FIT_LINEAR_W

(n,x,y,w,a,b)

c=========================================================================
c               Utility routine
c       Chi--square method for calculation of
c       parameters A,B in linear equation Y(x) = A*X + B
c           See MAPLE file linear\_w.mws
c       N          ---  number     of points
c       X(i), Y(i) ---  coordinate of point
c       W(i)       ---  weight     of point 
c=========================================================================
*   Here have to be added -- correct zero amplitudes account
*      and   chi-square calculations for fitted curve
c=========================================================================

---------------

SUBROUTINE FIT_EXP

(n,x,y,s,a,b)

c=========================================================================
c               Utility routine
c       Chi--square method for calculation of
c       parameters A,B in equation Y(x) = A*exp(-B*X)
c       Or Transversal Electromagnetic Shower model
c           See MAPLE file linear.mws
c       N    --- number of points
c       X(i) --- distance along shower
c       Y(i) --- amplitudes along shower
c       S(i) --- error of amplitude
c=========================================================================
*   Here have to be added -- correct zero amplitudes account
*      and   chi-square calculations for fitted curve
c=========================================================================

---------------

SUBROUTINE SDIRCO

(X1,Y1,Z1,X2,Y2,Z2,CX,CY,CZ)

c----------------------------------------------------------------------
C          Line equation THROUGH THE TWO POINTS 
c----------------------------------------------------------------------