TWIST MOFIA MANUAL
. Analysis Branch
The procedure dplot is the TWIST event analysis subroutine. From dplot
calls are made to analyze the event starting with the TDC unpacking,
filling the histograms, filtering the event, etc.
.1 TDC Unpacking
From dplot a call is made to the PUBLIC function tdcunp in module
tdc_mod, in order to unpack the TDCs. This module also contains the
declarations and initializations (and filling) of the data structures
associated with the TDC hit structures which are shown in Figure
9. Table 5 shows a brief description of the tdc_type structure. This
type has two instantiations, DCtdc(ihit) and PCtdc(ihit). The first
two components of the structure are the time and width of the TDC
signal. The time stored in these structures is in nanoseconds, and is
measured relative to the delayed (by ~10 ?s) trigger signal.The third
component, flag, is an integer that is assigned a zero value for
normal TDC signals, and a non-zero value if the TDC signal has a
peculiar characteristic (such as a leading edge but no trailing edge,
etc). The last two components in this structure are two pointers. The
first, wireP, points back to the wire geometry structure and,
therefore, provides the link between the hit structure and the
geometry structures. To point to the wire and plane numbers for a
particular DC TDC hit, for example, we have
INTEGER (i4), POINTER :: iwP, ipP
iwP => DCtdc(ihit)%wireP%iwire
ipP => DCtdc(ihit)%wireP%planeP%iplane
and so on. In the above example we could have defined integer
variables iw and ip (instead of pointers iwP and ipP); in which case
we have
INTEGER (i4) :: iw , ip
iw = DCtdc(ihit)%wireP%iwire
ip = DCtdc(ihit)%wireP%planeP%iplane
However, since copying data from a structure to a new variable is
generally less efficient than using a pointer, it is recommended that
the above style be used.
The contents of tdc_type are well documented in mainf90/tdc_mod.f90.
A description of the "flag" values and their meanings is in a table
here:
http://twist.triumf.ca/private/offline/mofia/general/tdc_type_flags.html
The second pointer in the tdc_type structure, whitsP, provides the
link to the whits_type structure. This structure is also defined and
filled in the module tdc_mod. The first element in this structure,
nhits, is the number of TDC hits on that wire and the second element,
hits(index), contains the hit number as it appears in the hit list in
the tdc_type structure. This is best understood through an
example. The number of hits on wire 32 in DC plane 26 is then
DCwhits(26,32)%nhits
To assign two pointers, ihitP and timeP, to the hit index and TDC time
for the second hit on this wire we have
INTEGER (i4), POINTER:: ihitP
REAL (r4):: timeP
ihitP => DCwhits(26,32)%hits(2)
timeP => DCtdc(ihitP)%time
The above example demonstrates another case of assigning pointers to a
component of a structure. While the last two lines may have been
combined in one, so that
timeP => DCtdc(DCwhits(26,32)%hits(2))%time
the style of assigning a pointer makes the code more readable. It is
also more efficient, since typically such statements would be inside
do loops (looping over planes, wires, hits, etc), and hence the hit
number is extracted from the structures once by assigning it to a
pointer(ihitP, in this case), and the pointer is then used within that
loop from that point on. Another useful technique is to assign a
pointer to a structure (as opposed to a component of a structure in
the examples above). This also provides faster execution time; and
improves the readability of the code. For example, one can declare the
pointer whitP and assign it to a structure
DO iPlane = 1, Ndplanes
DO iWire = 1, Ndwires(iPlane)
whitP => DCWhits(iPlane,iWire)
nhitsP => whitP%nhits
END DO
END DO
Scintillator hit times are stored in the SCtdc structure which is
similar to the TDC structure above, and differs only in that the
scintP pointer to the scintillator geometry replaces wireP.
The two structures SCadc and PCadc hold ADC data from the PACT
modules, unpacked and interpreted to give the energy deposited in the
detector. These structures are similar to the TDC structures above
with the two components time and width replaced by the component
e_lost which holds the energy deposited in KeV.
Additionally the raw TDC data is unpacked into DCtdc_raw, PCtdc_raw,
SCtdc_raw, PUtdc_raw, PCadc_raw and SCadc_raw whenever the namelist
variable rawout > 0. This is intended for use in debugging and
possibly determination of time calibration. The structure PUtdc_raw
holds pulser information.
Procedures to add and remove entries from these structures, as well as
the type definitions, are found in tdc_mod.f90 in the directory
mainf90. Tables relating scintillator numbers, plane and wire numbers,
as well as pulsar information to the TDC slot and address locations
are handled in tdcmap_mod.f90 in the directory mainf90. This code
reads in the mapping information through CFM. The CFM types
corresponding to these map files are FBC1_MAP, FBC2_MAP, and FBC3_MAP.
.2 Filtering
Following the call to tdcunp, a call is made to FillRawHists (and its
user counterpart, uFillRawHists) where some histograms are filled
before any event filtering is done. The event is then filtered by
calling the PUBLIC subroutine Filters in module
mainf90/filters_mod.f90, followed by a call to FillHists and
uFillHists to fill some histograms after event filtering.
The public subroutine Filters makes calls to three PRIVATE
subroutines: FiltersInit initializes the various counters used in this
module, FiltersCounters calculates these counters, and FiltersApply is
the subroutine where the event filters are applied. These filters
include: scintillator filters, drift chamber filters, proportional
chamber filters and RF filters. Several tests are performed on each
event to determine whether it passes or fails a filter/cut. Values for
these cuts are determined by the user through the namelist variables
which will be discussed below. Counters are maintained for events
failing a certain cut and statistics of these failures may be
displayed by typing show fail at the MOFIA command line, as in the
example below
MOFIA> show fail
Event Statistics:
ICFAIL= 0 ( GOOD EVENT ) 2072 2072
ICFAIL= 2 ( NO HITS IN CHAMBERS ) 5 5
ICFAIL=11 ( BAD EVENT ) 684 684
ICFAIL=12 ( BAD EVENT ) 13 13
ICFAIL=14 ( BAD EVENT ) 225 225
ICFAIL=15 ( BAD EVENT ) 1 1
The first column of numbers shows the number of events failing the
filters since the last analyze command was entered while the second
column shows the total since the MOFIA session was started. In
addition, a histogram (ID=5000) is incremented every time an event
fails a filter.
The best source for documentation is in mainf90/filters_mod.f90
itself; search for "TEST NN" where NN is the failure code (e.g. "TEST
12").
.3 Cross Talk
Following event filtering a call is made to the PUBLIC subroutine
Xtalk in the module mainf90/xtalk_mod.f90. Xtalk calls two PRIVATE
subroutines, XtalkInit which is used to initialize some counters, and
XtalkAnalyze which performs an analysis of each hit to determine if it
is a cross talk hit, increments the cross talk counters for each plane
and wire in the DC, and removes the hit from the data structures if it
is determined to be a cross talk hit. Three criteria are used to
determine whether a hit is a likely cross talk hit. First, the hit is
required to have a TDC width shorter than a user imposed value
determined by the variable DC_XTALK_WCUT in the namelist
DCCUTS. Second, the hit is required to occur in a cell adjacent to one
with a hit that has a TDC width longer than DC_XTALK_WCUT. Third, the
suspected cross talk hit is required to coincide in time with the hit
in the adjacent cell. Note, however, that the user must provide a
value for DC_XTALK_WCUT, otherwise no cross talk analysis will be
performed.
The user can initialize or print out the cross talk percentages at any
time during a MOFIA session through one of the "func" commands on the
MOFIA command line (type "func" alone in mofia for the list of
func's). Output files will be created and a message will be displayed
on the screen notifying the user of the file names.
.4 Calibrations
The calibrations branches of MOFIA are controlled by a set of flags in
the namelist DCCFLAGS (DC Calibration FLAGS). These branches are not
normally executed when running MOFIA, and the user has to turn a
specific flag on to execute a given calibration branch. The
calibrations flags control MOFIA branches that are used to compute
efficiency, plane positions, wire positions, plane rotations, DC
chamber resolution and time zero.
MOFIA> show name dccflags
.4.1 Efficiency
The efficiency code relies on the tracking to compute plane and wire
efficiencies for both the DCs and the PCs. After a track is
successfully reconstructed a call is made to the PUBLIC subroutines
EffDC and EffPC in the module user/efficiency_mod.f90. These
subroutines make calls to EffDCinit and EffPCinit to initialize the
efficiency counters, and to EffDCcalc and EffPCcalc to calculate the
efficiencies. The call to EffDC and EffPC, however, is controlled by
the namelist flag FindEff in namelist efficiency, so that if this
namelist variable is set to false the efficiency code will not get
executed. In this case if the user attempts to output the efficiency
counter (using func 12) a message will appear on the screen informing
the user that the FindEff flag has to be turned on before the data is
analyzed. The show command may be used to display the contents of
namelist efficiency
MOFIA> show name efficiency
The efficiency code uses the reconstructed track parameters to
traverse through the detector and find the cells intersected by the
track. It then checks whether these cells have hits and increments the
appropriate counters shown in the structures diagram of Figure 11. In
order to avoid edge effects (i.e. cases where the track has actually
exited the active area of the chamber, but the track parameters point
to the first/last cell due to tracking errors) two parameters are
provided to impose an edge cut. These are RadiusCutDense which allows
the user to set a radius cut on planes in the dense stack that have
only 48 instrumented wires, and RadiusCutSparse to place a radius cut
on the sparse stack in which all 80 wires are instrumented. The
variable CellCut allows the user to determine how many adjacent cells
should be investigated if no hit is found where the reconstructed
track parameters point.
.4.2 Time Zero
.4.3 Alignments
Both rotational and translational alignments are currently implemented
in MOFIA. When any of the variables FindPlanePos, FindWirePos or
FindPlaneRot in the namelist DCCFLAGS is turned on a call is made to
the PUBLIC subroutine Align in the module mainf90/align_mod.f90. Align
calls AlignInit, a PRIVATE initialization subroutine in the module,
followed by a call to one or more of the PRIVATE subroutines
AlignPlaneShifts, AlignWireShifts, and AlignPlaneRotations depending
on whether the corresponding namelist flags are turned on. These are
the subroutines that perform calculations of plane and wire
translational position corrections (in the U and V directions) and the
plane rotational corrections (around the z-axis). Once the plane
translational and rotational alignments are determined, the results
are stored in calibration files that are read into MOFIA through
CFM. The CFM types corresponding to these calibrations are DC_PPC and
DC_PRC, for the DC plane position corrections and plane rotation
corrections, respectively.
The command show name alignment shows the contents of the alignment
namelist.
MOFIA> show name alignment
.4.3.1 Translational Alignments
After analyzing a number of events determined by the user through the
variable nEventsPlane in the namelist alignment the subroutine
AlignPlaneShifts is called to determine an average residual for each
plane from tracks reconstructed successfully. The subroutine proceeds
by installing the average residuals as a correction to the plane and
wire positions. The tracking then proceeds as normal until another
nEventsPlane are analyzed and the call is made again to
AlignPlaneShifts to determine and install the new corrections. This
iteration continues until the specified runs (or events) are
analyzed. The user can also iterate on the same file by simply using
the MOFIA commands to rewind and reanalyze the file as many times as
the user desires. The variable FixPlanes in namelist alignment also
allows the user to fix two U and two V planes if the user so desires,
and to specify which planes should be fixed through the variables
FixedPlane1, FixedPlane2, FixedPlane3 and FixedPlane4 in the same
namelist. The corrections for the fixed planes will not be installed,
instead their nominal positions will be used.
.4.3.2 Rotational Alignments
For rotational alignments the subroutine AlignPlaneRotations is called
once nEventsPlane are analyzed. In this case an average residual is
determined for each plane, binned in 7 bins along the length of a
wire. For each plane the average residuals in each length bin can be
used to determine the plane rotational correction.
.4.4 Resolution
Resolution for DME is strongly dependent on drift distance being
mainly influenced by two processes: ionization statistics and
diffusion. At distances close to the wire ionization statistics
dominate, while diffusion becomes dominant near the edge of the cell
where the electric field is weak. The strong variation in resolution
across the cell requires the binning of the residuals along the
distance from the wire. The resolution code currently uses a bin width
of 100 microns.
When the calibration flag FindRes in the namelist DCCFLAGS is turned
on, a call is made to the PUBLIC subroutine Resolution in module
mainf90/resolution_mod.f90. This subroutine makes a call to
ResolutionDC which proceeds to calculate the resolution by examining
the residuals in each of the distance bins. One of three methods can
be chosen by the user to analyze the residuals: fitting the binned
residuals to a gaussian, performing a squared sum calculation, or
determining the FWHM of the binned distributions. Each of these
methods has its advantages and disadvantages. For example, the
gaussian method would be appropriate if the resolution is constant
within a certain bin (or the bin width is infinitesmal). Details on
the resolution calculations will be provided in a technical note.
The method used to analyze the residuals is determined by the variable
calcResidualMethod in namelist rezs. Also in this namelist the user
can choose the number of events required to be analyzed before the
residuals are calculated. This is determined by the variable
nEventsMax.
Once the resolution (defined by sigma of a guassian fit, sigma from a
squared sum, or sigma from FWHM) is determined, the tracking proceeds
again using the new resolutions in the track fitting until nEventsMax
are analyzed at which point the subroutine Resolution is called again
to determine and install the new resolution. The iterations continue
until all the specified events/runs are analyzed with each iteration
using the resolutions determined from the previous iteration in the
track fit.
The command show name rezs shows the contents of the rezs namelist
MOFIA> show name rezs
.5 Pattern Recognition
The purpose of the pattern recognition is to assign hits to tracks and
provide starting track parameters to the fitting routine. More
information on Pattern Recognition is coming 'soon' from Jim Musser;
see the main Mofia Documentation page.
.6 Tracking
Both a chi squared fit and a Kalman filter are used for tracking
purposes. Currently the Kalman filter is used for straight track
fitting and the chi squared fit is used for helix track fitting. The
intention is to modify the Kalman filter in the future to handle helix
tracks as well.
.6.1 Chi Squared Fit
Ongoing development tests of the Chi Squared fitter by Konstantin
Olchanski can be found on the web at
http://twist.triumf.ca/~olchansk/twist/
.6.2 Kalman Filter
Details on the Kalman filter are coming 'soon' from Maher Quraan.