Timing Reconstruction of ZDC Pre-Processor Data and its Calibration
S.White-Aug.2,2010

This notebook runs as a time calibration application on ATLAS ZDC data taken with LED calibration pulses in a series of 1 nsec delay steps using PHOS4 (see N. Massol, ATLAS Calibration Delay Chip Study, ATL-ELEC-2004-001 and references therein). PHOS4 steps are assumed to be calibrated but, if calibration data are available, they could be used to correct the vector "delay". Some studies of A. Poblaguev suggest small differences in timing behavior-such as waveform shape- between physics and LED data but, for the moment, the best calibration we have is from the LED delay scans. It is difficult to do timing calibration directly with physics data since they are concentrated in a narrow time region. This may be possible by varying the clock phase to the PPM.
	The algorithm we propose for time reconstruction was chosen to give a consistent method over the full available time range. The ZDC data can be used to obtain 80 MHz sampling of the ZDC waveform, which is approximately 100 nsec wide due to the 320m long ethernet ("DRAKA") cables we used to bring signals from the tunnel. At the present time we are only working with 40 MHz samples but, as Andrei showed in the ZDC performance plots ATL-COM-LUM-2010-022 , this already gives ~210 psec resolution/PMT with an algorithm designed for t~0. The purpose of the present algorithm is to extend this performance to all t.
	For a bandlimited analog signal with sparse sampling there are standard practices in electronics based on the Nyquist-Shannon sampling theorem. This was originally discussed  in Claude Shannon's paper "Communication in the Presence of Noise",Proc. Institute of Radio Engineers, vol. 37, no. 1, pp. 10–21, Jan. 1949. Reprinted as classic paper in: Proc. IEEE, vol. 86, no. 2, (Feb. 1998), which is a classic in information theory. Shannon derives an optimum reconstruction algorithm based on the Sinc function, where Sinc[x]=Sin[x]/x and->1 as x->0,MouseOver. This function is scaled to give roots at all other sampled points (assuming equal spacing) and each sampled point contributes a Sinc function with an amplitude equal to the sample and centered on that point to the series representation of the reconstructed waveform. Tektronix digital scopes use this algorithm on-chip to produce high resolution waveforms. Shannon's sampling theorem states that once the step size reaches 1/(2f), where f is the highest frequency component no further information is gained by finer sampling. ZDC data, even at 80 MHz are not quite at this limit. The following timing algorithm is based on Shannon's method.

Needs["PlotLegends`"];
SetDirectory[NotebookDirectory[]];
slices=ConstantArray[0,{22,6}];
timeslices=ConstantArray[0,{22,6,2}];

time=Range[0,125,25];
delay={0,1,2,3,4,5,6,7,8,9,11,13,14,15,16,17,18,19,20,21,22,23};

(*slices[[1]]={190,610,375,200,125,80};
slices[[2]]={160,620,380,205,130,95};
slices[[3]]={140,615,390,210,125,80};
slices[[4]]={120,615,395,210,130,85};
slices[[5]]={97,620,405,220,130,80};
slices[[6]]={80,612,420,225,140,90};
slices[[7]]={62,610,425,235,140,95};
slices[[8]]={50,605,435,235,145,95};
slices[[9]]={37,590,450,240,150,97};
slices[[10]]={30,575,460,245,150,97};
slices[[11]]={15,550,485,260,155,100};
(*slices[[12]]={12,530,590,265,160,100};*)
slices[[12]]={4,495,495,275,160,100};
slices[[13]]={2,495,515,275,165,105};
slices[[14]]={2,465,520,275,165,110};
slices[[15]]={2,445,525,290,170,110};
slices[[16]]={2,420,570,315,180,120};
slices[[17]]={2,385,550,310,175,115};
slices[[18]]={2,365,565,320,180,115};
slices[[19]]={2,335,575,325,185,120};
slices[[20]]={2,300,590,330,185,120};
slices[[21]]={2,280,595,340,195,125};
slices[[22]]={2,245,600,350,200,125};*)
(*Export["DelayScan.dat",slices]*)

slices={{190,610,375,200,125,80},{160,620,380,205,130,95},{140,615,390,210,125,80},
{120,615,395,210,130,85},{97,620,405,220,130,80},{80,612,420,225,140,90},{62,610,425,235,140,95},
{50,605,435,235,145,95},{37,590,450,240,150,97},{30,575,460,245,150,97},{15,550,485,260,155,100},
{4,495,495,275,160,100},{2,495,515,275,165,105},{2,465,520,275,165,110},{2,445,525,290,170,110},
{2,420,570,315,180,120},{2,385,550,310,175,115},{2,365,565,320,180,115},{2,335,575,325,185,120},
{2,300,590,330,185,120},{2,280,595,340,195,125},{2,245,600,350,200,125}};

Do[
Do[timeslices[[i,j,1]]=time[[j]];timeslices[[i,j,2]]=slices[[i,j]],{j,6}],{i,22}]

From the raw PPM data and even approximations up to interpolation order 3 it is hard to see a clear trend in the delay steps. An animation of this can be found at:

The Sinc function should be scaled to give roots at the next sample.

FindRoot[Sinc[t],{t,1}]

Use sinc interpolation with roots at neighboring sample times. An animation can be found at:

tpeak=ConstantArray[0,22];peakcts=ConstantArray[0,22];peakint=ConstantArray[0,22];a2=ConstantArray[0,22];

ShannonFilm=Animate[slice=slices[[i]];
shannon[t_]:=Sum[slice[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
Plot[shannon[t],{t,0,120},ImageSize->{500,300},PlotRange->{{0,125},{0,650}},Frame->True,GridLines->{time,None},PlotStyle->Thick,FrameLabel->{Style["time(nsec)",18],Style["reconstructed amplitude",18],Style["Waveform Reconstruction, Shannon method",18]},LabelStyle->Directive[Black,Bold,FontSize->18]],{{i,1,"delay step(nsec)"},1,22,1}]

tpeak=ConstantArray[0,22];peakcts=ConstantArray[0,22];peakint=ConstantArray[0,22];a2=ConstantArray[0,22];
ShFilm=Table[Manipulate[slice=slices[[i]];
shannon[t_]:=Sum[slice[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
Plot[shannon[t],{t,0,120},ImageSize->{500,300},Frame->True,GridLines->{time,None},PlotStyle->Thick,FrameLabel->{Style["time(nsec)",18],Style["rconstructed amplitude",18],Style["Waveform Reconstruction, Shannon method",18]},LabelStyle->Directive[Black,Bold,FontSize->18]],{{i,ii,""},1,22,1}],{ii,1,21,1}];

Export["ShannonFilm.gif",ShFilm,"DisplayDurations"->3/21]

Find time and amplitude at the Maximum  of the interpolation function. Both the peak and the integral of the function give stable measurement of the energy deposit in a ZDC module - better than using 1 or 2 slices. This is the ZDC signal reconstruction method.

Do[slice=slices[[i]];
shannon[t_]:=Sum[slice[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->20];
peakcts[[i]]=a[[1]];peakint[[i]]=NIntegrate[shannon[t],{t,0,120}]/60;a2[[i]]=slices[[i,2]];
tpeak[[i]]=t/.a[[2]];,{i,1,22}];

ListPlot[{a2,peakcts,peakint},Joined->True,PlotRange->{{0,22},{0,650}},ImageSize->{500,300},Frame->True,PlotStyle->Thick,FrameLabel->{Style["delay step(nsec)",18],Style["rconstructed Energy(a.u.)",18],Style["Energy Reconstruction, Shannon method vs. A2",18]},LabelStyle->Directive[Black,Bold,FontSize->18],PlotLegend->{"a2","reconstructed peak","Integral,Shannon"},LegendPosition->{-.6,-.2},LegendTextSpace->2.5,LegendLabel->Style["Energy Method",14],LegendLabelSpace->.5,LegendOrientation->Vertical,LegendBackground->LightPurple,LegendShadow->{.02,-.02},LegendSize->{0.5,0.5}]

peakfind=ConstantArray[0,{22,2}];

Do[peakfind[[i,2]]=delay[[i]];peakfind[[i,1]]=tpeak[[i]]-tpeak[[1]];,{i,1,22}];

We then plot actual PHOS4 delay time vs. the time obtained from this signal reconstruction method.There is a clear non - linearity in the interval 0 -> 25 nsec.Without a correction one makes an error in time scale of ~2 in the region where most physics data are concentrated.

ListPlot[peakfind,PlotRange->{{0,24},{0,25}},ImageSize->{500,300},Frame->True,(*Joined->True,*)PlotStyle->PointSize[.02],FrameLabel->{Style["reconstructed time( nsec)",18],Style["PHOS4 delay(nsec)",18],Style["uncorrected, Shannon method, peak",18]},LabelStyle->Directive[Black,Bold,FontSize->18]]

In the following we find a piecewise function which maps PPM reconstructed time to "true" time based on the PHOS4 delays.

peaksub=ConstantArray[0,{10,2}];
Do[peaksub[[i,2]]=delay[[i]];peaksub[[i,1]]=peakfind[[i,1]];,{i,1,10}];
Pfit1[t_]=Fit[peaksub,{1,t,t^2,t^3},t];
Plot[Pfit1[t],{t,0,6},PlotStyle->Thick,Epilog->{Directive[PointSize[Large],Red],Point[peaksub]}]

The derivative at t=0 is 1.5. Therefore an uncalibrated algorithm, like the one used in Fig.8 of our ATLAS note, will have the wrong time scale by this factor.

Pfit1'[0]

peaksub=ConstantArray[0,{8,2}];
Do[peaksub[[i,2]]=delay[[i+9]];peaksub[[i,1]]=peakfind[[i+9,1]];,{i,1,8}];
Pfit2[t_]=Fit[peaksub,{1,t,t^2,t^3},t];
Plot[Pfit2[t],{t,5,17},PlotStyle->Thick,Epilog->{Directive[PointSize[Large],Red],Point[peaksub]}]

peaksub=ConstantArray[0,{7,2}];
Do[peaksub[[i,2]]=delay[[i+15]];peaksub[[i,1]]=peakfind[[i+15,1]];,{i,1,7}];
Pfit3[t_]=Fit[peaksub,{1,t,t^2,t^3},t];
Plot[Pfit3[t],{t,15,25},PlotStyle->Thick,Epilog->{Directive[PointSize[Large],Red],Point[peaksub]}]

The output of this program is the following function which should be used to correct PPM time, reconstructed by the above method, for physics data. It gives small residuals which are partly due to the low statistics used in this study. It represents the best knowledge we have of ZDC time calibration. The calibration could be run separately on all channels.

CalibrationFunction[t_]:=Piecewise[{{Pfit1[t],t<5},{Pfit2[t],t>5&&t<16},{Pfit3[t],t>16}}]

CalibrationFunction[t]

Plot[CalibrationFunction[t],{t,-1,25},ImageSize->{500,300},Frame->True,FrameLabel->{Style["PPM Reconstructed time",18],Style["PHOS4 delay",18],Style["Derivation of Calibration Function",18]},LabelStyle->Directive[Black,Bold,FontSize->18],PlotStyle->Thick,Epilog->{Directive[PointSize[Large],Red],Point[peakfind]}]
Residual=ConstantArray[0,22];

Do[Residual[[i]]=peakfind[[i,2]]-CalibrationFunction[peakfind[[i,1]]],{i,1,22}]
{Histogram[Residual,10],ListPlot[Residual,PlotRange->{{-1,24},{-1,1}},PlotStyle->PointSize[.03]]}

tp=ConstantArray[0,10];tpraw=ConstantArray[0,10];Resol=ConstantArray[0,22];
Do[
i=1;
slice=slices[[i]];
Do[
slicenoise=Table[slice[[jj]]+RandomReal[2],{jj,6}];
shannon[t_]:=Sum[slicenoise[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->9];
tpraw[[n]]=t/.a[[2]];
,{n,1,10}];
tp=Table[tpraw[[k]]-Mod[tpraw[[k]],25]+CalibrationFunction[Mod[tpraw[[k]],25]],{k,10}];
If[i==1,ListPlot[Transpose[{tpraw,tp}]],];
tpeak[[i]]=Mean[tp];
tp;
Resol[[i]]=StandardDeviation[tp];
(*Do[*)
i=3;
slice=slices[[i]];
Do[
slicenoise=Table[slice[[jj]]+RandomReal[2],{jj,6}];
shannon[t_]:=Sum[slicenoise[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->9];
tpraw[[n]]=t/.a[[2]];
,{n,1,10}];
tp=Table[tpraw[[k]]-Mod[tpraw[[k]],25]+CalibrationFunction[Mod[tpraw[[k]],25]],{k,10}];
If[i==1,ListPlot[Transpose[{tpraw,tp}]],];
tpeak[[i]]=Mean[tp];
tp;
Resol[[i]]=StandardDeviation[tp];
i=5;
slice=slices[[i]];
Do[
slicenoise=Table[slice[[jj]]+RandomReal[2],{jj,6}];
shannon[t_]:=Sum[slicenoise[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->9];
tpraw[[n]]=t/.a[[2]];
,{n,1,10}];
tp=Table[tpraw[[k]]-Mod[tpraw[[k]],25]+CalibrationFunction[Mod[tpraw[[k]],25]],{k,10}];
If[i==1,ListPlot[Transpose[{tpraw,tp}]],];
tpeak[[i]]=Mean[tp];
tp;
Resol[[i]]=StandardDeviation[tp];
i=7;
slice=slices[[i]];
Do[
slicenoise=Table[slice[[jj]]+RandomReal[2],{jj,6}];
shannon[t_]:=Sum[slicenoise[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->9];
tpraw[[n]]=t/.a[[2]];
,{n,1,10}];
tp=Table[tpraw[[k]]-Mod[tpraw[[k]],25]+CalibrationFunction[Mod[tpraw[[k]],25]],{k,10}];
If[i==1,ListPlot[Transpose[{tpraw,tp}]],];
tpeak[[i]]=Mean[tp];
tp;
Resol[[i]]=StandardDeviation[tp];
i=9;
slice=slices[[i]];
Do[
slicenoise=Table[slice[[jj]]+RandomReal[2],{jj,6}];
shannon[t_]:=Sum[slicenoise[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->9];
tpraw[[n]]=t/.a[[2]];
,{n,1,10}];
tp=Table[tpraw[[k]]-Mod[tpraw[[k]],25]+CalibrationFunction[Mod[tpraw[[k]],25]],{k,10}];
If[i==1,ListPlot[Transpose[{tpraw,tp}]],];
tpeak[[i]]=Mean[tp];
tp;
Resol[[i]]=StandardDeviation[tp];
i=11;
slice=slices[[i]];
Do[
slicenoise=Table[slice[[jj]]+RandomReal[2],{jj,6}];
shannon[t_]:=Sum[slicenoise[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
a=NMaximize[{shannon[t],10<t<60},t,WorkingPrecision->9];
tpraw[[n]]=t/.a[[2]];
,{n,1,10}];
tp=Table[tpraw[[k]]-Mod[tpraw[[k]],25]+CalibrationFunction[Mod[tpraw[[k]],25]],{k,10}];
If[i==1,ListPlot[Transpose[{tpraw,tp}]],];
tpeak[[i]]=Mean[tp];
tp;
Resol[[i]]=StandardDeviation[tp]
,{i,3,3}];
Resol

ListPlot[Resol]

Addendum: Leading Edge Timing
The same methods were tried for leading edge constant fraction algorithms. All require similar non-linearity corrections, as seen below. On a Macintosh G5 running Mathematica 7.0, with working precision of 20 the leading edge algorithms (using FindRoot) run 10 times faster (9 milliseconds) than the peak algorithm (which uses NMaximize). This is a topic for further study.

50% Constant Fraction

thalf=ConstantArray[0,22];th=ConstantArray[0,{22,2}];
Timing[Do[slice=slices[[i]];
shannon[t_]:=Sum[slice[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
shsub[t_]:=(shannon[t]-peakcts[[i]]/2);
b=FindRoot[shsub[t],{t,10},WorkingPrecision->20];thalf[[i]]=t/.b[[1]];
th[[i,1]]=delay[[i]];th[[i,2]]=thalf[[i]]-5.758;,{i,1,22}];]

75% Constant Fraction

tqua=ConstantArray[0,22];tq=ConstantArray[0,{22,2}];
Do[slice=slices[[i]];
shannon[t_]:=Sum[slice[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
shsub[t_]:=(shannon[t]-peakcts[[i]]*.75);
b=FindRoot[shsub[t],{t,20}];tqua[[i]]=t/.b[[1]];
tq[[i,1]]=delay[[i]];tq[[i,2]]=tqua[[i]]-13.1446;,{i,1,22}];

20% Constant Fraction

teight=ConstantArray[0,22];te=ConstantArray[0,{22,2}];
Do[slice=slices[[i]];
shannon[t_]:=Sum[slice[[j]]*Sinc[Pi*(t-time[[j]])/25],{j,1,6}];
shsub[t_]:=(shannon[t]-peakcts[[i]]*.2);
b=FindRoot[shsub[t],{t,10}];teight[[i]]=t/.b[[1]];
te[[i,1]]=delay[[i]];te[[i,2]]=teight[[i]]+3.84;,{i,1,22}];
{ListPlot[th],ListPlot[tq],ListPlot[te]}

(*ListPlot[{tq,th,te},Joined->True,PlotRange->{{0,22},{0,22}},ImageSize->{500,300},Frame->True,PlotStyle->Thick,FrameLabel->{Style["delay step(nsec)",18],Style["uncorrected reconstructed time(nsec)",18],Style["edge method, uncorrected",18]},LabelStyle->Directive[Black,Bold,FontSize->18],PlotLegend->{"25%","50%","80%"}]*)

The following is a template for non-Mathematica implementation:

c	Reconstruction Routine
	Subroutine ZDC_Rec(Slice,nSlice,time,Energy,method,status_code)
c	Slice=address of 1st slice
c	nSlice= number of slices(5,6, or 7)
c	time=reconstructed time depending on method
c	Energy=Integral of waveform in units of energy or PPM counts
c	method=Reconstruction method:
c				1=peak, not linearized
c				2=peak, linearized
c				3=leading edge from 50% constant fraction (not linearized)
c				4=leading edge from 75% constant fraction (not linearized)
c				5=leading edge from 20% constant fraction (not linearized)
c	status_code:
c				0= normal return
c				1=time out of range? Probably not needed.
c				2=amplitude in a region likely to yield poorer timing resolution
	shannon[τ] = Σ (n = 1, nSlice)]slice(n)*Sinc[Pi*(τ - (n-1*25))/25]
	Energy=(Subscript[Σ, (n=1, nSlice)])slice(n)
	t=FindMaximum[shannon[τ], starting at time=50 or maximum slice]
	peak=shannon[t]
	status_code=0
	Go To(1,2,3,4,5) method
1	time=t
	return
2	time=
(Piecewise[{{0.05683627487211110724860754038130085448`14.159053595649281+1.53276368648127925515764204895952799188`15.440255802587249 t+0.22980834373505649499444710063000431752`14.990277351584679 (t^2)-0.03656366471191917165312482196647781764`15.017450350511368 (t^3), t<5}, {-1.28974955223497072897926985556163683083`15.045345725799459+2.65357869933460418084110241462447658569`15.887895476808172 t-0.13712401461402091041472481255492873603`15.61383882578186 (t^2)+0.00281211806384422081427932826912270437`15.407090932179447 (t^3), t>5&&t<16}, {-42.18688740322649664921002481082898043578`15.642502116011583+8.61752055701945705349939015890034206935`15.770529113051126 t-0.42592398060653285221504730872878903511`15.751550579986562 (t^2)+0.00753849507486616550061783241461220913`15.753249992664164 (t^3), t>16}, {0, True}}])
	return
3	tempfunction[τ]=
	time=Findroot[shannon[τ] -peak/2, (in interval from 0-> 50)]
	return
4	tempfunction[τ]=
	time=Findroot[shannon[τ] -peak*.75, (in interval from 0-> 50)]
	return
5	tempfunction[τ]=
	time=Findroot[shannon[τ] -peak*.2, (in interval from 0-> 50)]
	return
	end