[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

code for estimating momentum resolution

To: halld-offline@jlab.org
Subject: code for estimating momentum resolution
From: "Mark M. Ito" <marki@jlab.org>
Date: Tue, 12 Feb 2008 14:46:53 -0500
Organization: Jefferson Lab
Reply-To: "Mark M. Ito" <marki@jlab.org>
Sender: owner-halld-offline@jlab.org
User-Agent: Thunderbird 2.0.0.9 (X11/20071115)

I've checked in a set of FORTRAN routines to estimate momentum and angle 
resolution in the CDC and FDC (separately for now). You can get the code 
from subversion:

  svn co https://halldsvn.jlab.org/repos/trunk/home/marki/gluex/rezest

Look for rezest.F in the "rezest" directory created by the 
aforementioned command. Or you can get it from  my personal web area:

  http://www.jlab.org/~marki/halld/rezest.F

or extract the attachment to this message.

Documentation on use is in comments at the beginning of the source code.

C These routines calculate resolutions in various track parameters due
C to both position resolution and multiple scattering individually as
C well as the total (i. e., their sum in quadrature).
C
C The subroutines are:
C
C      SUBROUTINE REZEST_CURVE(DET, LAMBDA, L, EPSILON, N_M, P, M, N_RL,
C     X     B, DK_MS, DK_RES, DK)
C
C      SUBROUTINE REZEST_AZIMUTH(DET, LAMBDA, L, EPSILON, N_M, P, M,
C     X     N_RL, N_RL_FRONT, B, R_M, DK, DPHI_MS, DPHI_RES, DPHI_K,
C     x     DPHI)
C
C      SUBROUTINE REZEST_POLAR(DET, LAMBDA, L, EPSILON, STEREO,
C     X     N_M, P, M, N_RL, N_RL_FRONT, DTHETA_MS, DTHETA_RES,
C     X     DTHETA)
C
C All have the following REAL*4 input arguments in common:
C
C     DET      Detector, CHARACTER*3, either 'CDC' or 'FDC'
C     LAMBDA   Dip angle, difference in polar angle in lab between track
C              and pi/2 (i. e., 90 degrees) (radians)
C     L        Length of detector
C              CDC: the distance from the inner radius of the chamber to the
C                   outer radius,
C              FDC: the distance from the front of the chamber to the back of
C                   the chamber
C              (meters)
C     EPSILON  Position resolution (meters)
C     N_M      Number of measurements (real number)
C     P        Magnitude of total momentum (GeV/c)
C     M        Mass of the particle (GeV)
C     N_RL     Number of radiation lengths in the detector (real number)
C
C In addition REZEST_CURVE and REZEST_AZIMUTH require as input the REAL*4
C argument:
C
C     B        Magnetic field (Tesla)
C
C REZEST_AZIMUTH and REZEST_POLAR require a REAL*4 argument:
C
C     N_RL_FRONT Number of radiations lengths in front of the detector
C
C Finally, REZEST_AZIMUTH requires a REAL*4 argument:
C
C     DK       Resolution in curvature, as output from REZEST_CURVE (see below)
C              (inverse meters)
C
C Output arguments are all REAL*4:
C
C REZEST_CURVE: DK_MS, DK_RES, DK
C               curvature resolution due to multiple scattering, position
C               resolution and their sum in quadrature respectively
C               (inverse meters)
C
C REZEST_AZIMUTH: DPHI_MS, DPHI_RES, DPHI
C               resolution in azimuthal angle due to multiple scattering,
C               position resolution and their sum in quadrature respectively
C               (radians)
C
C REZEST_POLAR: DTHETA_MS, DTHETA_RES, DTHETA
C               resolution in polar angle due to multiple scattering,
C               position resolution and their sum in quadrature respectively
C               (radians)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE REZEST_CURVE(DET, LAMBDA, L, EPSILON, N_M, P, M, N_RL,
     X     B, DK_MS, DK_RES, DK)
C     Calculates resolution in curvature
      IMPLICIT NONE
      INCLUDE "lun.par"
      CHARACTER*3 DET
      REAL*4 LAMBDA, L, EPSILON, N_M, P, M, N_RL, B
      REAL*4 E, GAMMA, BETA
      REAL*4 L_PRIME, DK_RES, DK_MS, DK, N_RL_TOTAL
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      L_prime = L*cos(lambda)
      dk_res = epsilon/L_prime**2*sqrt(720/(n_m + 4))
      IF (DET .EQ. 'CDC') THEN
         N_RL_TOTAL = N_RL/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         N_RL_TOTAL = N_RL/COS(3.14159/2.0 - LAMBDA)
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP   
      ENDIF
      dk_ms = 0.016/(L*p*beta*cos(lambda)**2)*sqrt(n_rl_TOTAL)
      dk = sqrt(dk_res**2 + dk_ms**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE REZEST_AZIMUTH(DET, LAMBDA, L, EPSILON, N_M, P, M,
     X     N_RL, N_RL_FRONT, B, R_M, DK, DPHI_MS, DPHI_RES, DPHI_K,
     x     DPHI)
      IMPLICIT NONE
      INCLUDE "lun.par"
      CHARACTER*3 DET
      REAL*4 LAMBDA, L, EPSILON, N_M, P, M, N_RL, N_RL_FRONT, B, R_M,
     X     DK
      REAL*4 E, GAMMA, BETA
      REAL*4 THETA0, C, DPHI_MS, DPHI_RES, DPHI_K, DPHI,
     X     N_RL_TOTAL, L_PRIME
      DATA C/1.0/
      L_PRIME = L*COS(LAMBDA)
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      IF (DET .EQ. 'CDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(3.14159/2.0 - LAMBDA)
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP   
      ENDIF
      N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      THETA0 = 0.0136/(BETA*C*P)*SQRT(N_RL_TOTAL)
     X     *(1.0 + 0.038*LOG(N_RL_TOTAL))
      DPHI_MS = THETA0/SQRT(3.0)
      DPHI_RES = EPSILON/L_PRIME*SQRT(12*(N_M - 1)/(N_M*(N_M + 1)))
      DPHI_K = R_M*DK
      DPHI = SQRT(DPHI_RES**2 + DPHI_MS**2 + DPHI_K**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE REZEST_POLAR(DET, LAMBDA, L, EPSILON, STEREO,
     X     N_M, P, M, N_RL, N_RL_FRONT, DTHETA_MS, DTHETA_RES,
     X     DTHETA)
      IMPLICIT NONE
      INCLUDE "lun.par"
      CHARACTER*3 DET
      REAL*4 LAMBDA, L, EPSILON, STEREO, N_M, P, M, N_RL,
     X     N_RL_FRONT
      REAL*4 E, GAMMA, BETA
      REAL*4 C, N_RL_TOTAL, THETA0, L_EFF, DTHETA_MS, DTHETA_RES, DTHETA
      DATA C/1.0/
      CALL KINEMATICS(P, M, E, GAMMA, BETA)
      IF (DET .EQ. 'CDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(LAMBDA)
      ELSE IF (DET .EQ. 'FDC') THEN
         N_RL_TOTAL = (N_RL + N_RL_FRONT)/COS(3.14159/2.0 - LAMBDA)
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP   
      ENDIF
      THETA0 = 0.0136/(BETA*C*P)*SQRT(N_RL_TOTAL)
     X     *(1.0 + 0.038*LOG(N_RL_TOTAL))
      DTHETA_MS = THETA0/SQRT(3.0)
      IF (DET .EQ. 'CDC') THEN
         L_EFF = SIN(STEREO)*L
      ELSE IF (DET .EQ. 'FDC') THEN
         L_EFF = L
      ELSE
         WRITE (6,*) 'BAD DETECTOR TYPE'
         STOP
      ENDIF
      DTHETA_RES = EPSILON/L_EFF
     X     *SQRT(12*(N_M - 1)/(N_M*(N_M + 1)))
      DTHETA = SQRT(DTHETA_RES**2 + DTHETA_MS**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE KINEMATICS(P, M, E, GAMMA, BETA)
      REAL*4 P, M
      REAL*4 E, GAMMA, BETA
      REAL*4 C
      DATA C/1.0/
      E = sqrt(p**2 + m**2)
      gamma = E/m
      beta = sqrt(1.0 - 1/gamma**2)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

begin:vcard
fn:Mark M. Ito
n:Ito;Mark
org:Jefferson Lab
adr:Suite 8;;12000 Jefferson Avenue;Newport News;VA;23601;USA
email;internet:marki@jlab.org
title:Staff Scientist
tel;work:757-269-5295
tel;fax:757-269-6331
tel;pager:757-584-5295
url:http://www.jlab.org/~marki
version:2.1
end:vcard

Prev by Date: Software meeting today
Next by Date: web interface for browsing subversion repository
Prev by thread: Re: Software meeting today
Next by thread: web interface for browsing subversion repository
Index(es):
- Date
- Thread