Ciao a tutti,
As reported yesterday at the Rome neutrino working group meeting,
a new version of the ERP calibration code has been completed and will be
put to work in order to produce the new CALMOD calibration constants
for runs 16543 (14-OCT-1998) - today. Please note that the numbers
given below are taken from the 1st week already analyzed with the new code
(CALMOD reference run # 16543).
I would like to list some of the features of this - hopefully
better and improved - new code. In addition, I would like (once again)
to go into some detail on how the code works.
(1) The code is still written in FORTRAN - yet, it is no longer a single
program but instead is split in 4 files:
DISK$MACRO:[MACROCAL.CALCODE.NEW.AXP]muon.for (main program, 253 lines)
DISK$MACRO:[MACROCAL.CALCODE.NEW.AXP]muon_lib.for (55 subroutines and
functions, 6012 lines)
DISK$MACRO:[MACROCAL.CALCODE.NEW.AXP]muon.inc (common blocks shared
among subroutines, functions and main program, 102 lines)
DISK$MACRO:[MACROCAL.CALCODE.NEW.AXP]calmod.inc (common blocks providing
the interface with CALMOD libraries, 230 lines)
(2) Great effort has been put in documenting the code. Structured
programming is used - in fact there is not a single "GOTO"
command in this new code.
(3) The data used for all calibrations are taken from the .TRK files that
correspond to 1 week's worth of data (~30 runs).
As a reminder, the .TRK files are produced at the end of each run from
the BMON procedure by applying standard tracking to all events where
ERP hits are present and storing the streamer-tube position of the
track for each scintillator counter with an ERP hit.
These .TRK files are first processed in order to keep ONLY 2-face and
3-face hits with exactly 1 ERP hit per face, with consistent streamer
tube track and a minimum pathlength in MACRO of 200cm.
(4) For all calibrations there are fundumental "sanity" checks: that the
hit is a real muon (and not some calibration event), that the ADC and
TDC values are within their dynamic range etc.
(5) The code is split in 2 major parts:
- Energy calibration (ADC-A gains, ADC-D gains & pedestals)
- Timing calibrations (refractive index, TDC-H offsets, TDC-L offsets and
gains, InterERP TDC offsets and gains)
a) Energy calibration is performed in 3 phases:
(E1) L/R ADC-A gain ratio (and optionally, pedestal) calibration.
This is achieved by using muon hits in the center 600cm. of the
scintillator counters. The position in the counter, as reconstructed
by the streamer-tube track, is combined with the counter response
function in order to obtain the "expected p.e. ratio" for each hit.
Then, we try to find the relative ADC-A gain ratio between side-0
and side-1 (and, optionally the corresponding pedestals) such that
the difference between the "ERP-calculated p.e. ratio" and the
"measured p.e. ratio" is minimum (in the LMS sense).
MINUIT is used, where 2 error functions are tested:
( (ADC0 - ped0)*g/(ADC1 - ped1) - ratio ) ^2
( log((ADC0-ped0)*g/(ADC1-ped1) - log(ratio) ) ^2
where ped0,ped1 are the ADC-A pedestals for the 2 sides,
g is the relative gain of the 2 sides (gain0/gain1) and
ratio is the "expected p.e. ratio" mentioned above.
The two error functions reach practically the same minimum.
By fixing the 2 pedestals, minization is achieved only w.r.t. g.
Note that the MINUIT minization is applied twice: once with and
the second time without limits to guarantee convergence and to
obtain better error estimates (see MINUIT manual for details).
In addition, adaptive 3-sigma cuts are applied to the data points
before the final fit result is accepted.
The resulting L/R gain ratios are nicely peaked aroung 1.0, with
some boxes giving ratios outside the [0.5,2.0] range: these counters
are either Attico verticals, where we couldn't set their gain for
the last 2 years (they have no clear s.p.e. peak), or horizontal
counters where, for some reason, one of the 2 PMTs is not working
properly.
Note in particular the very good match (<5% difference) of these
ratios with the ones calculated by PHRASE using a similar procedure.
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i) The old code was only using hits in the center 100 cm of the
counter, requiring a (mean) p.e. ratio = 1.0. Thus, it was
too sensitive to low statistic fluctuations.
(ii) A single fit was applied (no sigma-cuts).
(E2) Absolute energy calibration - done through a fit to the Landau curve.
Only tracks with at least 1cm in the scintillator counter are used.
Then, a histogram of energy depositions (MeV/cm) is filled for
each counter. A fit with respect to the Landau distribution is
then performed - actually, 2 fits; one with and another without
bounds. Note that 3 free parameters are used for the
Landau curve: total area (proportional to the number of hits),
distribution "width" and distribution peak.
The last parameter (peak) is then required to be = 1.8 MeV/cm.
If different, all ERP ADC gains are normalized to give this value.
Then, a new histogram is filled to verify that this is indeed the
case; if not, the procedure is applied iteratively.
Final fit results are very good, with values of chi^2 around 0.7
and no box with chi^2 more than 1.32 (2B04).
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i)There was a single fit to the Landau curve, where the limits
of the histogram to be filled were fixed (0 - 10 MeV/cm).
The new code is adaptively determining the limits of the
histogram, by calculating the sample mean of the points
and then using the range (0 - 5*mean) to build the histogram.
In this way, boxes 1E05 and 1N05 that had "escaped" to
unphysical regions (50 MeV/cm) were never calibrated because
the corresponding histograms resulted empty.
(ii) Fitting is done iteratively until it converges to the
desired value (old code was performing a single fit.)
New code is using 200 bins for the histogram - old was
using 100 bins.
(E3) ADC-D / ADC-A equalization - done through a line fit of the
linear range of the two ADC values for each side of a scintillator
counter. Note that the resulting slope should be (and indeed, is)
very close to the value 0.091 (11:1), since the "dynode" signal
is obtained from the "anode" signal through a resistor divider.
The trickiest part is to choose the hits to use for calibration.
For this reason, adaptive 2-sigma cut is applied (the result is
only accepted when the data points within 2-sigma result in the
same fit parameters between 2 successive iterations).
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i) Multiple fits guarantee convergence.
(ii) Fitting is performing in the ADC-D vs. ADC-A domain instead
of the p.e. domain (in theory equivalent - in practice it is
easier to control "good" hits based on their raw ADC values)
b) Timing calibration is performed in 4 phases
(T1) TDC-H L/R offset and refr. index calibration - done through a
line fit of the streamer-tube position (STZ) in the counter and
that reconstructed from the TDC-H's (TDCZ).
We only use hits in the center 800cm of the counter (to avoid
problems related to the end-of-counter light propagation) with
consistent track (i.e. crossing the scintillator counter) and
TDC-H and TDC-L values (in the range [200,3000] ticks).
The fit is iteratively performed as follows:
6 "passes" are applied, where on each pass stricter conditions
on the TDCZ and STZ aggreement is required: for the first 2
passes all hits are accepted as follows:
1st, 2nd pass : infinity
3rd pass : 150cm
4th pass : 100cm
5th, 6th pass : 50cm
For each pass, adaptive 3-sigma cut is applied to the data points
to remove outliers. This actually turns out to be a much stricter
condition than the ones mentioned above, since the average sigma
for these fits is about 9cm, with NO box exceeding 20cm.
The resulting refr. indices are very nicely peaked around the
value 1.56, with a spread of 0.02. There are only 3 boxes that
exceed the value 1.65: 2B04 (1.75), 3B07 (1.65) and 3B10 (1.65),
3 of the 4 "former body leakers".
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i) Addition of 3-sigma cuts, application of all corrections
to refr. index and L/R TDC offset during various iterations
and only check for very small changes at the end (previous
code was prohibiting small changes DURING iterations).
(T2) TDC-H intra-SM box offset calibration.
Maybe the "heart" of the whole ERP muon calibration, since it is
using the "fundumental" condition of beta = 1.0 for through-going
muons.
Only 2-face or 3-face hits within the same SM (actually, the same
ERP supervisor, i.e. North and South face boxes are excluded)
are used, where all scintillator counters are crossed by the
reconstructed streamer-tube track, with a tolerance between STZ
and TDCZ of 50cm.
Note that 3-face hits are decomposed in 3 distinct pairs.
Only pairs with at least 200cm pathlength are kept.
At first, all such hits are stored in a big buffer, containing
all the necessary information for every hit: the 2 boxes (a and b)
involved, ordered from top to bottom, assuming all such hits are
downward going muons, time-of flight as reconstructed by TDC-H's
tof = (tb0+tb1)/2.0 - (ta0+ta1)/2.0 and pathlength.
Note that this buffer contains 100,000 - 150,000 entries for a
week's worth of muon data, with ~ 1000 pairs for horizontal counters
and 400 hits for vertical ones.
Then, an iterative algorithm is applied in 5 passes.
Scanning through the buffer mentioned above, for each box the
difference in time, dt, for all muons that reconstruct a beta in
the range [1.0-pass*0.1,1.0+pass*0.1] is calculated as follows:
dt = tof(nsec) - pathlength(cm)/29.98
The mean and sigma for this dt quantity is then calculated for
all boxes, and the corresponding mean is then subtracted from
each box.
This procedure is iteratively applied until ALL boxes have
an average dt < 0.00005 nsecs.
In addition, 3-sigma cuts are applied to the hits to remove
outliers.
Note that the resulting offsets have an arbitrary "SM offset"
to them - at the end it is calculated and subtracted from
all boxes such that the sum of all TDC-H offsets for each
SM is equal to 0.0
The typical sigma at the end of this calibration is 0.6 nsecs,
with NO box exceeding 1.0 nsec (the worst one is 2B04 with ~ 1.0 nsec).
It's worth mentioning that the resulting box offsets have a clear
systematic, with the ATTICO boxes having an offset roughly 40 nsecs
before the LOWER boxes, for all 6 SMs.
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i) Addition of 3-sigma cuts, application of all corrections
to box TDC offset during various iterations and only
check for very small changes at the end (previous code
was prohibiting small changes DURING iterations).
(ii) All boxes are corrected for each iteration - previous
code was only adjusting the box with the worst offset
at each iteration. As a result, previous code required
~15,000 iterations - new one ~ 100 iterations
(iii) The old code requested that the sum of ALL MACRO box
offsets is equal to zero - the new one that the sum
of ALL BOXES OF EACH SM is equal to zero. As a result,
on the previous box offsets, there was a strange 80 nsec
offset for all SM3 boxes, with some negative offset (-15nsec)
for the other 5 SMs.
(T3) InterERP constants and N/S TDC-H box offset calibration.
Very similar to intra-SM box offset calibration.
Only hits that cross two neighboring SMs are considered.
The same accumulation buffer mentioned above is filled with
all necessary information. Note that the 14 N/S boxes are
treated as 14 different combinations.
MINUIT minimization is then performed w.r.t. 29 parameters:
10 interERP TDC slopes (N face, SM1, SM2 on uVAX#1, SM2 or uVAX#2,
SM3, SM4, SM5 on uVAX#2, SM5 on uVAX#3, SM6 and S face)
5 SM-pair offsets (SM1/SM2, SM2/SM3, SM3/SM4, SM4/SM5 and SM5/SM6).
14 N/S box offsets (1N01/SM1, 1N02/SM1, ... SM6/6S06, SM6/6S07).
The error function to minimize is again the sum of dt^2, where
dt is calculated as above, adding the interERP terms.
Again, MINUIT is applied iteratively, discarding at each iteration
hits with dt outside the 3-sigma window, until there is no
improvement in the residual error.
For each iteration, MINUIT is run twice: with and without bounds
to ensure convergence.
It takes roughly 10 iterations to converge, with a typical
final sigma of about 0.8 nsecs, with the South face having the
worst performance of ~ 0.9 nsecs.
From the 14 N/S box offsets, a global "face" offset is extracted
and then these boxes are normalized to hace a 0.0 mean offset.
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i) Previous code was doing one 3-sigma cut - new one
iteratively. Old code was performing bounded minimization,
new code bounded + unbounded.
(T4) TDC-L offset and gain calibration.
Using "consistent" muon hits (with |TDC-H - TDC-L| < 80 ticks),
a line fit between time-walk-corrected TDC-L and TDC-H values
is performed for each side of a scintillator counter separately.
Then, the slope and intercept is interpreted in terms of TDC gain
and offset. Note that the slope of the fit, which should be equal
to the ratio of TDC-L / TDC-H gains, should be the same with the
one calculated from the LED data. If not, then there is either
some bad time-walk constant (for small differences, of the order
of 3-5%) or some hardware problem.
This line fit results in remarkable consistency of the two methods
(muon data / LED data), with NO counter having discrepancy more
than 3%.
The typical sigma of this fit is of the order of 1.7 ticks (~0.25
nsecs) with only 13 channels (sides of counters) exceeding 5 ticks
(1 nsec).
The worst counter is 2B04-0, with 11 ticks (~2nsecs).
Just like the previous calibrations, the line-fit is performed
iteratively, applying 3-sigma cuts to discard outliers between
iterations. Less than 5 iterations usually suffice to converge.
DIFFERENCES WITH RESPECT TO THE OLD VERSION OF THE CODE:
(i) Previous code was performing one line fit in the time
domain. New one is iteratively applying 3-sigma cuts in
the TDC (tick) domain that is easier to understand in
terms of bad hits and outliers.
(6) Please note that the new code treats ALL 476 MACRO boxes, without
excluding "bad" ones (N/S, Attico verticals etc.) By the same token,
people who do analysis should use all scintillator counters. If
by doing so there are some "bad" boxes emerging that need to be
excluded or re-calibrated to give satisfactory results, they
should immediately communicate their "micro-cuts" to us, the
calibration crew at Gran Sasso: using the same data, one should get
the same quantities; unless there are problems in the software.
(6) The new code is producing plots of all the quantities mentioned above,
together with a detailed printout of all fit quantities: number of
total hits, number of hits passing the various selection criteria
and sigma cuts, fit results and quality (error/sigma/chi^2) values.
All these quantities are available to everybody and should shortly
be automatically posted on the web ("MACRO ON LINE").
A short messages, giving the highlight of every week will also be
produced automatically and mailed to the interested parties.
(7) The idea, once we catch up with the current data, is to have
all calibration code run automatically every week and CALMOD ('work' DB)
be automatically updated with the resulting constants. In this way,
people who do analyses can immediately run their software and let
un know about the quality of calibration constants.
(7) Since this is a rather major revision of the calibration code,
any feedback is more than welcome. Any comments, questions and
such should be directed to
katsavou@lngs.infn.it & zaccheo@lngs.infn.it
Ioannis Katsavounidis