1 /*******************************************************************************
2 * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
4 * Author: The IceCube RALICE-based Offline Project.
5 * Contributors are mentioned in the code where appropriate.
7 * Permission to use, copy, modify and distribute this software and its
8 * documentation strictly for non-commercial purposes is hereby granted
9 * without fee, provided that the above copyright notice appears in all
10 * copies and that both the copyright notice and this permission notice
11 * appear in the supporting documentation.
12 * The authors make no claims about the suitability of this software for
13 * any purpose. It is provided "as is" without express or implied warranty.
14 *******************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////
20 // TTask derived class to perform hit extraction from waveforms.
22 // The code in this processor is based on the algorithms as developed by
23 // Nick van Eijndhoven and Garmt de Vries-Uiterweerd (Utrecht University, The Netherlands).
25 // Procedure applied for Amanda TWR data :
26 // ---------------------------------------
28 // 1) The waveform is fed to a TSpectrum object, and the peak locations
29 // are determined with the TSpectrum::Search() function.
31 // 2) The waveform is divided into regions corresponding to the peaks found by
32 // TSpectrum. The region boundary between two peaks is at the location of
33 // the minimum between the two peaks.
35 // 3) For each region the "effective baseline" (used in the
36 // evaluation of the leading edge value) is determined as :
37 // effective baseline = fBasefracXXX * value at lower region boundary.
38 // This takes into account the effect from the previous pulse.
39 // For the first pulse, the effective baseline is equal to the overall
42 // 4) For each region, the point of steepest rise between the lower region
43 // boundary and the peak location is determined. The tangent at this point
44 // is extrapolated to the effective baseline. The point of intersection yields the
47 // 5) For each region the range of charge integration is determined as :
48 // - Start of integration at the lower region boundary or at the leading edge,
49 // whichever comes last;
50 // - End of integration at the upper region boundary or at the point where the
51 // signal drops below the overall baseline, whichever comes first.
53 // 6) For each region the integrated charge is determined as :
54 // Sum over bins in integration range of (value in bin - overall baseline).
56 // 7) For each pulse the quality is evaluated by requiring that :
57 // peak location - lower region boundary > lower region boundary - leading edge.
58 // For a too shallow steepest rise, the leading edge value is unreliable, in
59 // which case the pulse is merged with the previous pulse.
61 // 8) Each pulse is checked for saturation and discarded if necessary.
63 // 9) TSpectrum needs a minimum number of bins for its Search function, otherwise
64 // the clipping window is too large, which causes an error. If a waveform does
65 // not contain enough bins, the following alternative approach is used :
66 // - A loop over all bins is performed.
67 // - As soon as the signal exceeds a given threshold, a pulse is started.
68 // - The pulse ends when the signal drops below the threshold again.
69 // - While looping, the charge is integrated for each pulse.
71 // The defaults of the various parameters can be changed by the corresponding
72 // Set memberfunctions.
74 // Information about the actual parameter settings can be found in the event
75 // structure itself via the device named "IceMakeHits".
77 //--- Author: Nick van Eijndhoven and Garmt de Vries-Uiterweerd 15-jan-2007 Utrecht University
78 //- Modified: GdV $Date$ Utrecht University
79 ///////////////////////////////////////////////////////////////////////////
81 #include "IceMakeHits.h"
82 #include "Riostream.h"
84 ClassImp(IceMakeHits) // Class implementation to enable ROOT I/O
86 IceMakeHits::IceMakeHits(const char* name,const char* title) : TTask(name,title)
88 // Default constructor.
97 ///////////////////////////////////////////////////////////////////////////
98 IceMakeHits::~IceMakeHits()
100 // Default destructor.
102 ///////////////////////////////////////////////////////////////////////////
103 void IceMakeHits::SetBasefracA(Float_t val)
105 // Set baseline fractional update for Amanda TWR extraction.
106 // The default as set in the constructor of this class is 0.5.
109 ///////////////////////////////////////////////////////////////////////////
110 void IceMakeHits::SetSigmaA(Float_t val)
112 // Set clipping window width for Amanda TWR extraction.
113 // The default as set in the constructor of this class is 1.5.
116 ///////////////////////////////////////////////////////////////////////////
117 void IceMakeHits::SetMaxPeaksA(Int_t val)
119 // Set maximum number of peaks in a waveform for Amanda TWR extraction.
120 // The default as set in the constructor of this class is 10.
123 ///////////////////////////////////////////////////////////////////////////
124 void IceMakeHits::SetMinPulseHeightA(Float_t val)
126 // Set minimum required pulse height for Amanda TWR extraction.
127 // This is used only for narrow pulses that cannot be handled with TSpectrum.
128 // The default as set in the constructor of this class is 50.
129 fMinPulseHeightA=val;
131 ///////////////////////////////////////////////////////////////////////////
132 void IceMakeHits::SetThresholdA(Float_t val)
134 // Set threshold for use in analysis of narrow pulses for Amanda TWR extraction.
135 // A peak is assumed to start when the signal rises above threshold*maxval,
136 // where maxval is the maximum value found in the waveform.
137 // The default as set in the constructor of this class is 0.2.
140 ///////////////////////////////////////////////////////////////////////////
141 void IceMakeHits::Exec(Option_t* opt)
143 // Implementation of the hit cleaning procedures.
146 AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
150 fEvt=(IceEvent*)parent->GetObject("IceEvent");
153 fDaq=(AliDevice*)fEvt->GetDevice("Daq");
156 // Storage of the used parameters in the IceMakeHits device
158 params.SetNameTitle("IceMakeHits","IceMakeHits processor parameters");
159 params.SetSlotName("BasefracA",1);
160 params.SetSlotName("SigmaA",2);
161 params.SetSlotName("MaxPeaksA",3);
162 params.SetSignal(fBasefracA,1);
163 params.SetSignal(fSigmaA,2);
164 params.SetSignal(fMaxPeaksA,3);
166 fEvt->AddDevice(params);
173 ///////////////////////////////////////////////////////////////////////////
174 void IceMakeHits::Amanda()
176 // Hit extraction from the Amanda TWR data.
178 // Arrays for storing info
179 Float_t* baseline=new Float_t[fMaxPeaksA];
180 Int_t* lowend=new Int_t[fMaxPeaksA];
181 Int_t* upend=new Int_t[fMaxPeaksA];
182 Int_t* startcharge=new Int_t[fMaxPeaksA];
183 Int_t* stopcharge=new Int_t[fMaxPeaksA];
184 Int_t* status=new Int_t[fMaxPeaksA]; // 0=OK, 1=rejected, 2=saturation
185 Float_t* leadingedge=new Float_t[fMaxPeaksA];
186 Float_t* charge=new Float_t[fMaxPeaksA];
187 Float_t* tot=new Float_t[fMaxPeaksA];
189 // Some objects and variables we will need
191 TSpectrum spec(fMaxPeaksA);
192 Int_t nrIterations=(Int_t)(7*fSigmaA+0.5); // Number of iterations used in TSpectrum::SearchHighRes()
193 Int_t npeaks=0, ibin=0, lookforsteepestuntilbin=0, steep=0;
194 Float_t maxval=0, rise=0, rc=0, yyy=0;
195 Bool_t pulsegoingon=false;
196 Int_t* index=new Int_t[fMaxPeaksA];
198 // Update the DAQ device data for the nature of the hit info
199 fDaq->AddNamedSlot("TWR");
200 fDaq->SetSignal(1,"TWR");
201 Int_t idx=fDaq->GetSlotIndex("Muon");
202 if (idx) fDaq->SetSignal(0,idx);
204 // All Amanda OMs with a signal
205 TObjArray* aoms=fEvt->GetDevices("IceAOM");
208 // OM, waveform and hit
212 hit.SetSlotName("ADC",1);
213 hit.SetSlotName("LE",2);
214 hit.SetSlotName("TOT",3);
216 // Loop over all fired OMs and extract the hit info
217 for (Int_t iom=0; iom<aoms->GetEntries(); iom++)
219 omx=(IceAOM*)aoms->At(iom);
221 // Remove all existing hits of this OM
223 // Should we skip OMs that we know from the dbase to have problems ?
224 //// if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
226 // Investigate all waveforms for this OM
227 for (Int_t iwf=1; iwf<=omx->GetNwaveforms(); iwf++)
229 wf=omx->GetWaveform(iwf);
231 maxval=wf->GetMaximum();
232 // Check if clipping window is not too large
233 if(wf->GetNbinsX() > 2*nrIterations+1)
235 // Find peaks with TSpectrum
236 npeaks=spec.Search(wf,fSigmaA,"goff");
237 // Discard waveform if no or too many peaks found
238 if(npeaks<1 || npeaks>fMaxPeaksA) continue;
240 // Get differential of WF
242 for(ibin=2;ibin<diff.GetNbinsX();ibin++)
244 diff.SetBinContent(ibin,wf->GetBinContent(ibin)-wf->GetBinContent(ibin-1));
246 diff.SetBinContent(1,0);
247 // Set baseline and lower end for first peak,
251 // Sort peaks in time
252 TMath::Sort(npeaks,spec.GetPositionX(),index,false);
253 // For each of the peaks,
254 for(Int_t ipeak=0; ipeak<npeaks; ipeak++)
256 // Find baseline and region around peak
258 // (Second and later peaks: lower edge = upper edge previous peak,
259 // baseline is average of previous baseline and minimum value between two
263 lowend[ipeak]=upend[ipeak-1]+1;
264 baseline[ipeak]=fBasefracA*foo.GetBinContent(lowend[ipeak]);
266 // (Upper edge range is minimum between this and next peak)
269 foo.SetAxisRange(spec.GetPositionX()[index[ipeak]],spec.GetPositionX()[index[ipeak+1]]);
270 upend[ipeak]=foo.GetMinimumBin()-1;
272 // (Last peak: upper edge is end of histo)
275 upend[ipeak]=wf->GetNbinsX();
277 // Find steepest rise
278 lookforsteepestuntilbin=wf->FindBin(spec.GetPositionX()[index[ipeak]]);
280 // Look for steepest rise between lower edge and peak position
281 foo.SetAxisRange(wf->GetBinCenter(lowend[ipeak]),wf->GetBinCenter(lookforsteepestuntilbin));
282 steep=foo.GetMaximumBin();
283 rise=foo.GetBinContent(steep);
285 // Extrapolate tangent to find leading edge
286 yyy=wf->GetBinContent(steep)-baseline[ipeak];
287 rc=rise/foo.GetBinWidth(steep);
288 if(rc>0) leadingedge[ipeak]=wf->GetBinCenter(steep)-yyy/rc; else leadingedge[ipeak]=0;
290 // Determine peak status
292 // Check for saturation
293 if(rc<0.1 && wf->GetBinContent(wf->FindBin(spec.GetPositionX()[index[ipeak]])) == maxval)
297 // Check quality: LE should not be too far below lower edge
298 // Otherwise, ignore this peak and set baseline back to what it was
299 else if(wf->GetBinLowEdge(lowend[ipeak]) - leadingedge[ipeak] > spec.GetPositionX()[index[ipeak]] - wf->GetBinLowEdge(lowend[ipeak]))
302 if(ipeak>0) baseline[ipeak]=baseline[ipeak-1];
305 // Start charge integration at LE, or at lower edge of range
306 startcharge[ipeak]=wf->FindBin(leadingedge[ipeak]);
307 if(lowend[ipeak]>startcharge[ipeak]) startcharge[ipeak]=lowend[ipeak];
309 // Integrate charge until pulse drop below baseline, or else until edge of range
310 stopcharge[ipeak]=upend[ipeak];
311 for(ibin=wf->FindBin(spec.GetPositionX()[index[ipeak]]); ibin<=upend[ipeak]; ibin++)
313 if(wf->GetBinContent(ibin)<0)
315 stopcharge[ipeak]=ibin-1;
320 // Determine time over threshold
321 tot[ipeak]=wf->GetBinLowEdge(stopcharge[ipeak]+1)-wf->GetBinLowEdge(startcharge[ipeak]);
325 for(ibin=startcharge[ipeak]; ibin<=stopcharge[ipeak]; ibin++)
327 charge[ipeak]+=wf->GetBinContent(ibin);
330 } // end loop over peaks
332 // Check all peaks, from latest to earliest
333 for(int ipeak=npeaks-1; ipeak>=0; ipeak--)
336 // If this peak was rejected, add charge and TOT to previous peak (if there is one)
337 if(status[ipeak]==1 && ipeak>0)
339 charge[ipeak-1]+=charge[ipeak];
341 tot[ipeak-1]+=tot[ipeak];
345 // If this peak is OK, add hit info
349 hit.SetSignal(charge[ipeak],"ADC");
350 hit.SetSignal(leadingedge[ipeak],"LE");
351 hit.SetSignal(tot[ipeak],"TOT");
355 } // end loop over peaks
356 // If number of bins too small, use different method
360 // If maximum value high enough to suspect presence of peak,
361 if(maxval>fMinPulseHeightA)
366 for(ibin=1; ibin<=wf->GetNbinsX(); ibin++)
368 // If bin content above threshold, start pulse
369 if(wf->GetBinContent(ibin)>fThresholdA*maxval){
374 leadingedge[npeaks]=wf->GetBinLowEdge(ibin);
375 charge[npeaks]=wf->GetBinContent(ibin);
380 charge[npeaks]+=wf->GetBinContent(ibin);
388 tot[npeaks]=wf->GetBinLowEdge(ibin)-leadingedge[npeaks];
390 // Store pulse information
392 hit.SetSignal(charge[npeaks],"ADC");
393 hit.SetSignal(leadingedge[npeaks],"LE");
394 hit.SetSignal(tot[npeaks],"TOT");
397 // Get ready for next pulse
403 } // End of loop over bins
406 } // End of alternative method for narrow pulses
415 delete[] startcharge;
418 delete[] leadingedge;
423 ///////////////////////////////////////////////////////////////////////////
424 void IceMakeHits::InIce()
426 // Hit extraction from IceCube InIce ATWD data.
428 ///////////////////////////////////////////////////////////////////////////
429 void IceMakeHits::IceTop()
431 // Hit extraction from IceTop ATWD data.
433 ///////////////////////////////////////////////////////////////////////////