Added many comments and some documentation, fixed up some coding violations
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
CommitLineData
b0f5e3fc 1#include <stdio.h>
8a00af9a 2#include <stdlib.h>
fd61217e 3#include <iostream.h>
b0f5e3fc 4#include <TObjArray.h>
1ca7869b 5#include <TRandom.h>
fd61217e 6#include <TMath.h>
b0f5e3fc 7
fd61217e 8#include "AliITSmodule.h"
9#include "AliITSMapA2.h"
9e8d6423 10#include "AliITSsegmentationSSD.h"
11#include "AliITSresponseSSD.h"
b0f5e3fc 12#include "AliITSsimulationSSD.h"
fd61217e 13//#include "AliITSdictSSD.h"
b0f5e3fc 14#include "AliITSdcsSSD.h"
15#include "AliITS.h"
16#include "AliRun.h"
0315d466 17#include "AliITSgeom.h"
b0f5e3fc 18
19ClassImp(AliITSsimulationSSD);
fd61217e 20////////////////////////////////////////////////////////////////////////
21// Version: 0
22// Written by Enrico Fragiacomo
23// July 2000
24//
25// AliITSsimulationSSD is the simulation of SSDs.
26
0315d466 27//----------------------------------------------------------------------
b0f5e3fc 28AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
29 AliITSresponse *resp){
0315d466 30 // Constructor
b0f5e3fc 31
32 fSegmentation = seg;
33 fResponse = resp;
0315d466 34 Float_t noise[2] = {0.,0.};
35 fResponse->GetNoiseParam(noise[0],noise[1]); // retrieves noise parameters
36 fDCS = new AliITSdcsSSD(seg,resp);
b0f5e3fc 37
38 fNstrips = fSegmentation->Npx();
39 fPitch = fSegmentation->Dpx(0);
fd61217e 40 fMapA2 = new AliITSMapA2(fSegmentation);
41
fd61217e 42 fSteps = 100; // still hard-wired - set in SetDetParam and get it via
b0f5e3fc 43 // fDCS together with the others eventually
b0f5e3fc 44}
0315d466 45//______________________________________________________________________
b0f5e3fc 46AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD
47 &source){
0315d466 48 // Operator =
fd61217e 49
b0f5e3fc 50 if(this==&source) return *this;
51
0315d466 52 this->fDCS = new AliITSdcsSSD(*(source.fDCS));
53 this->fMapA2 = source.fMapA2;
b0f5e3fc 54 this->fNstrips = source.fNstrips;
55 this->fPitch = source.fPitch;
56 this->fSteps = source.fSteps;
57 return *this;
58}
0315d466 59//_____________________________________________________________---------
b0f5e3fc 60AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
0315d466 61 // copy constructor
fd61217e 62
0315d466 63 *this = source;
b0f5e3fc 64}
0315d466 65//______________________________________________________________________
b0f5e3fc 66AliITSsimulationSSD::~AliITSsimulationSSD() {
0315d466 67 // destructor
fd61217e 68 delete fMapA2;
b0f5e3fc 69 delete fDCS;
0315d466 70}
71//_______________________________________________________________-------
b0f5e3fc 72void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
73 Int_t dummy) {
0315d466 74 // Digitizes hits for one SSD module
75
76 Int_t lay, lad, detect;
77 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
78 AliITSgeom *geom = aliITS->GetITSgeom();
79 geom->GetModuleId(module,lay, lad, detect);
80 if ( lay == 6 )((AliITSsegmentationSSD*)fSegmentation)->SetLayer(6);
81 if ( lay == 5 )((AliITSsegmentationSSD*)fSegmentation)->SetLayer(5);
82
83 TObjArray *hits = mod->GetHits();
84 Int_t nhits = hits->GetEntriesFast();
85 if (!nhits) return;
86 //cout<<"!! module, nhits ="<<module<<","<<nhits<<endl;
fd61217e 87
0315d466 88 Double_t x0=0.0, y0=0.0, z0=0.0;
89 Double_t x1=0.0, y1=0.0, z1=0.0;
90 Double_t de=0.0;
91 Int_t maxNdigits = 2*fNstrips;
92 Float_t **pList = new Float_t* [maxNdigits];
93 memset(pList,0,sizeof(Float_t*)*maxNdigits);
94 Int_t indexRange[4] = {0,0,0,0};
95 static Bool_t first=kTRUE;
96 Int_t lasttrack = -2;
97 Int_t idtrack = -2;
b0f5e3fc 98
0315d466 99 for(Int_t i=0; i<nhits; i++) {
100 // LineSegmentL returns 0 if the hit is entering
101 // If hits is exiting returns positions of entering and exiting hits
102 // Returns also energy loss
103
104 if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
105 HitToDigit(module, x0, y0, z0, x1, y1, z1, de, indexRange, first);
fd61217e 106
0315d466 107 if (lasttrack != idtrack || i==(nhits-1)) {
108 GetList(idtrack,pList,indexRange);
109 first=kTRUE;
110 } // end if
111 lasttrack=idtrack;
112 } // end if
113 } // end loop over hits
fd61217e 114
0315d466 115 ApplyNoise();
116 ApplyCoupling();
e8189707 117
0315d466 118 ChargeToSignal(pList);
b0f5e3fc 119
0315d466 120 fMapA2->ClearMap();
b0f5e3fc 121}
0315d466 122//----------------------------------------------------------------------
fd61217e 123void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0,
124 Double_t z0, Double_t x1, Double_t y1,
125 Double_t z1, Double_t de,
126 Int_t *indexRange, Bool_t first) {
0315d466 127 // Turns hits in SSD module into one or more digits.
128
129 Float_t tang[2] = {0.0,0.0};
130 fSegmentation->Angles(tang[0], tang[1]);// stereo<< -> tan(stereo)~=stereo
131 Double_t x, y, z;
132 Double_t dex=0.0, dey=0.0, dez=0.0;
133 Double_t pairs;
134 Double_t ionE = 3.62E-9; // ionization energy of Si (GeV)
b0f5e3fc 135
0315d466 136 Double_t sigma[2] = {0.,0.};// standard deviation of the diffusion gaussian
137
138 Double_t D[2] = {11.,30.}; // diffusion constant {h,e} (cm**2/sec)
139 Double_t tdrift[2] = {0.,0.}; // time of drift
140 Double_t vdrift[2] = {0.86E6,2.28E6}; // drift velocity (cm/sec)
141 Double_t w;
142 Double_t inf[2], sup[2], par0[2];
143 // Steps in the module are determined "manually" (i.e. No Geant)
144 // NumOfSteps divide path between entering and exiting hits in steps
145 Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez);
fd61217e 146
0315d466 147 // Enery loss is equally distributed among steps
148 de = de/numOfSteps;
149 pairs = de/ionE; // e-h pairs generated
150
151 for(Int_t j=0; j<numOfSteps; j++) { // stepping
152 // cout<<"step number ="<<j<<endl;
153 x = x0 + (j+0.5)*dex;
154 y = y0 + (j+0.5)*dey;
155 //if ( y > (seg->Dy()/2 +10)*1.0E-4 ) {
156 if ( y > (fSegmentation->Dy()/2+10)*1.0E-4 ) {
157 // check if particle is within the detector
158 cout<<"AliITSsimulationSSD::HitToDigit: Warning: hit "
159 "out of detector y0,y,dey,j ="
160 <<y0<<","<<y<<","<<dey<<","<<j<<endl;
161 return;
162 } // end if
163 z = z0 + (j+0.5)*dez;
164
165 // calculate drift time
166 // y is the minimum path
167 tdrift[0] = (y+(fSegmentation->Dy()*1.0E-4)/2) / vdrift[0];
168 tdrift[1] = ((fSegmentation->Dy()*1.0E-4)/2-y) / vdrift[1];
169
170 for(Int_t k=0; k<2; k++) { // both sides remember: 0=Pside 1=Nside
171
172 tang[k]=TMath::Tan(tang[k]);
173
174 // w is the coord. perpendicular to the strips
175 if(k==0) {
176 //w=(x+(seg->Dx()*1.0E-4)/2)-(z+(seg->Dz()*1.0E-4)/2)*tang[k];
177 w = (x+(fSegmentation->Dx()*1.0E-4)/2) -
178 (z+(fSegmentation->Dz()*1.0E-4)/2)*tang[k];
179 }else{
180 //w =(x+(seg->Dx()*1.0E-4)/2)+(z-(seg->Dz()*1.0E-4)/2)*tang[k];
181 w = (x+(fSegmentation->Dx()*1.0E-4)/2) +
182 (z-(fSegmentation->Dz()*1.0E-4)/2)*tang[k];
183 //cout<<"k,x,z,w ="<<k<<","<<x<<","<<z<<","<<w<<endl;
184 } // end if
185 w = w / (fPitch*1.0E-4); // w is converted in units of pitch
186
187 if((w<(-0.5)) || (w>(fNstrips-0.5))) {
188 // this check rejects hits in regions not covered by strips
189 // 0.5 takes into account boundaries
190 if(k==0) cout<<"AliITSsimulationSSD::HitToDigit: "
191 "Warning: no strip in this region of P side"
192 <<endl;
193 else cout<<"AliITSsimulationSSD::HitToDigit: "
194 "Warning: no strip in this region of N side"<<endl;
195 return;
196 } // end if
197
198 // sigma is the standard deviation of the diffusion gaussian
199
200 if(tdrift[k]<0) return;
201
202 sigma[k] = TMath::Sqrt(2*D[k]*tdrift[k]);
203 sigma[k] = sigma[k] /(fPitch*1.0E-4); //units of Pitch
204 if(sigma[k]==0.0) {
205 cout<<"AliITSsimulationSSD::DigitiseModule: Error: sigma=0"
206 <<endl;
207 exit(0);
208 } // end if
209
210 par0[k] = pairs;
211 // we integrate the diffusion gaussian from -3sigma to 3sigma
212 inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average
213 sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
214 // IntegrateGaussian does the actual
215 // integration of diffusion gaussian
216 IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k],
217 indexRange, first);
218 } // end for loop over side (0=Pside, 1=Nside)
219 } // end stepping
220 //delete seg;
b0f5e3fc 221}
0315d466 222//____________________________________________________________________--
b0f5e3fc 223void AliITSsimulationSSD::ApplyNoise() {
0315d466 224 // Apply Noise.
fd61217e 225
0315d466 226 Float_t signal;
227 Float_t noise[2] = {0.,0.};
228 fResponse->GetNoiseParam(noise[0],noise[1]); // retrieves noise parameters
fd61217e 229
0315d466 230 for(Int_t k=0;k<2;k++){ // both sides (0=Pside, 1=Nside)
231 for(Int_t ix=0;ix<fNstrips;ix++){ // loop over strips
232 signal = (Float_t) fMapA2->GetSignal(k,ix);// retrieves signal
233 // from map
fd61217e 234
0315d466 235 signal += gRandom->Gaus(0,noise[k]);// add noise to signal
236 if(signal<0.) signal=0.0; // in case noise is negative...
fd61217e 237
0315d466 238 fMapA2->SetHit(k,ix,(Double_t)signal); // give back signal to map
239 } // loop over strip
240 } // loop over k (P or N side)
b0f5e3fc 241}
0315d466 242//______________________________________________________________________
b0f5e3fc 243void AliITSsimulationSSD::ApplyCoupling() {
0315d466 244 // Apply the effect of electronic coupling between channels
245 Float_t signal, signalLeft=0, signalRight=0;
246
247 for(Int_t ix=0;ix<fNstrips;ix++){
248 if(ix>0.) signalLeft = (Float_t) fMapA2->GetSignal(0,ix-1)*
249 fDCS->GetCouplingPL();
250 else signalLeft = 0.0;
251 if(ix<(fNstrips-1)) signalRight = (Float_t) fMapA2->GetSignal(0,ix+1)*
252 fDCS->GetCouplingPR();
253 else signalRight = 0.0;
254 signal = (Float_t) fMapA2->GetSignal(0,ix);
255 signal += signalLeft + signalRight;
256 fMapA2->SetHit(0,ix,(Double_t)signal);
fd61217e 257
0315d466 258 if(ix>0.) signalLeft = (Float_t) fMapA2->GetSignal(1,ix-1)*
259 fDCS->GetCouplingNL();
260 else signalLeft = 0.0;
261 if(ix<(fNstrips-1)) signalRight = (Float_t) fMapA2->GetSignal(1,ix+1)*
262 fDCS->GetCouplingNR();
263 else signalRight = 0.0;
264 signal = (Float_t) fMapA2->GetSignal(1,ix);
265 signal += signalLeft + signalRight;
266 fMapA2->SetHit(1,ix,(Double_t)signal);
267 } // loop over strips
b0f5e3fc 268}
0315d466 269//______________________________________________________________________
fd61217e 270Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
0315d466 271 // Computes the integral of a gaussian using Error Function
272 Float_t sqrt2 = TMath::Sqrt(2.0);
273 Float_t sigm2 = sqrt2*s;
274 Float_t integral;
fd61217e 275
0315d466 276 integral = 0.5 * TMath::Erf( (x - av) / sigm2);
277 return integral;
278}
279//______________________________________________________________________
fd61217e 280void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w,
281 Double_t sigma,
282 Double_t inf, Double_t sup,
283 Int_t *indexRange, Bool_t first) {
0315d466 284 // integrate the diffusion gaussian
285 // remind: inf and sup are w-3sigma and w+3sigma
286 // we could define them here instead of passing them
287 // this way we are free to introduce asimmetry
288
289 Double_t a=0.0, b=0.0;
290 Double_t signal = 0.0, dXCharge1 = 0.0, dXCharge2 = 0.0;
291 // dXCharge1 and 2 are the charge to two neighbouring strips
292 // Watch that we only involve at least two strips
293 // Numbers greater than 2 of strips in a cluster depend on
294 // geometry of the track and delta rays, not charge diffusion!
fd61217e 295
0315d466 296 Double_t strip = TMath::Floor(w); // clostest strip on the left
297
298 if ( TMath::Abs((strip - w)) < 0.5) {
299 // gaussian mean is closer to strip on the left
300 a = inf; // integration starting point
301 if((strip+0.5)<=sup) {
302 // this means that the tail of the gaussian goes beyond
303 // the middle point between strips ---> part of the signal
304 // is given to the strip on the right
305 b = strip + 0.5; // integration stopping point
306 dXCharge1 = F( w, b, sigma) - F(w, a, sigma);
307 dXCharge2 = F( w, sup, sigma) - F(w ,b, sigma);
308 }else {
309 // this means that all the charge is given to the strip on the left
310 b = sup;
311 dXCharge1 = 0.9973; // gaussian integral at 3 sigmas
312 dXCharge2 = 0.0;
313 } // end if
314
315 dXCharge1 = par * dXCharge1;// normalize by mean of number of carriers
316 dXCharge2 = par * dXCharge2;
317
318 // for the time being, signal is the charge
319 // in ChargeToSignal signal is converted in ADC channel
320 signal = fMapA2->GetSignal(k,strip);
321 signal += dXCharge1;
322
323 fMapA2->SetHit(k,strip,(Double_t)signal);
324 if(((Int_t) strip) < (fNstrips-1)) {
325 // strip doesn't have to be the last (remind: last=fNstrips-1)
326 // otherwise part of the charge is lost
327 signal = fMapA2->GetSignal(k,(strip+1));
328 signal += dXCharge2;
329 fMapA2->SetHit(k,(strip+1),(Double_t)signal);
330 } // end if
fd61217e 331
0315d466 332 if(dXCharge1 > 1.) {
333 if (first) {
334 indexRange[k*2+0]=indexRange[k*2+1]=(Int_t) strip;
335 first=kFALSE;
336 } // end if first
337
338 indexRange[k*2+0]=TMath::Min(indexRange[k*2+0],(Int_t) strip);
339 indexRange[k*2+1]=TMath::Max(indexRange[k*2+1],(Int_t) strip);
340 } // dXCharge > 1 e-
341
342 }else{
343 // gaussian mean is closer to strip on the right
344 strip++; // move to strip on the rigth
345 b = sup; // now you know where to stop integrating
346 if((strip-0.5)>=inf) {
347 // tail of diffusion gaussian on the left goes left of
348 // middle point between strips
349 a = strip - 0.5; // integration starting point
350 dXCharge1 = F(w, b, sigma) - F(w, a, sigma);
351 dXCharge2 = F(w, a, sigma) - F(w, inf, sigma);
352 }else {
353 a = inf;
354 dXCharge1 = 0.9973; // gaussian integral at 3 sigmas
355 dXCharge2 = 0.0;
356 } // end if
fd61217e 357
0315d466 358 dXCharge1 = par * dXCharge1; // normalize by means of carriers
359 dXCharge2 = par * dXCharge2;
360
361 // for the time being, signal is the charge
362 // in ChargeToSignal signal is converted in ADC channel
363 signal = fMapA2->GetSignal(k,strip);
364 signal += dXCharge1;
365 fMapA2->SetHit(k,strip,(Double_t)signal);
366 if(((Int_t) strip) > 0) {
367 // strip doesn't have to be the first
368 // otherwise part of the charge is lost
369 signal = fMapA2->GetSignal(k,(strip-1));
370 signal += dXCharge2;
371 fMapA2->SetHit(k,(strip-1),(Double_t)signal);
372 } // end if
fd61217e 373
0315d466 374 if(dXCharge1 > 1.) {
375 if (first) {
376 indexRange[k*2+0]=indexRange[k*2+1]=(Int_t) strip;
377 first=kFALSE;
378 } // end if first
379
380 indexRange[k*2+0]=TMath::Min(indexRange[k*2+0],(Int_t) strip);
381 indexRange[k*2+1]=TMath::Max(indexRange[k*2+1],(Int_t) strip);
382 } // dXCharge > 1 e-
383 } // end if
b0f5e3fc 384}
0315d466 385//______________________________________________________________________
fd61217e 386Int_t AliITSsimulationSSD::NumOfSteps(Double_t x, Double_t y, Double_t z,
0315d466 387 Double_t & dex,Double_t & dey,Double_t & dez){
388 // number of steps
389 // it also returns steps for each coord
390 //AliITSsegmentationSSD *seg = new AliITSsegmentationSSD();
fd61217e 391
0315d466 392 Double_t step = 25E-4;
393 //step = (Double_t) seg->GetStepSize(); // step size (cm)
394 Int_t numOfSteps = (Int_t) (TMath::Sqrt(x*x+y*y+z*z)/step);
fd61217e 395
0315d466 396 if (numOfSteps < 1) numOfSteps = 1; // one step, at least
fd61217e 397
0315d466 398 // we could condition the stepping depending on the incident angle
399 // of the track
400 dex = x/numOfSteps;
401 dey = y/numOfSteps;
402 dez = z/numOfSteps;
b0f5e3fc 403
0315d466 404 return numOfSteps;
fd61217e 405}
0315d466 406//----------------------------------------------------------------------
407void AliITSsimulationSSD::GetList(Int_t label,Float_t **pList,
408 Int_t *indexRange) {
409 // loop over nonzero digits
410 Int_t ix,globalIndex;
411 Float_t signal=0.;
412 Float_t highest,middle,lowest;
fd61217e 413 // printf("SPD-GetList: indexRange[0] indexRange[1] indexRange[2] indexRange[3] %d %d %d %d\n",indexRange[0], indexRange[1], indexRange[2], indexRange[3]);
414
0315d466 415 for(Int_t k=0; k<2; k++) {
416 for(ix=indexRange[k*2+0];ix<indexRange[k*2+1]+1;ix++){
417 if(indexRange[k*2+0]<indexRange[k*2+1])
418 signal=fMapA2->GetSignal(k,ix);
fd61217e 419
0315d466 420 globalIndex = k*fNstrips+ix; // globalIndex starts from 0!
421 if(!pList[globalIndex]){
422 //
423 //Create new list (6 elements-3 signals and 3 tracks+total sig)
424 //
425 pList[globalIndex] = new Float_t [6];
426 // set list to -1
427 *pList[globalIndex] = -2.;
428 *(pList[globalIndex]+1) = -2.;
429 *(pList[globalIndex]+2) = -2.;
430 *(pList[globalIndex]+3) = 0.;
431 *(pList[globalIndex]+4) = 0.;
432 *(pList[globalIndex]+5) = 0.;
433 *pList[globalIndex] = (float)label;
434 *(pList[globalIndex]+3) = signal;
435 }else{
436 // check the signal magnitude
437 highest = *(pList[globalIndex]+3);
438 middle = *(pList[globalIndex]+4);
439 lowest = *(pList[globalIndex]+5);
440 signal -= (highest+middle+lowest);
441 //
442 // compare the new signal with already existing list
443 //
444 if(signal<lowest) continue; // neglect this track
445 if (signal>highest){
446 *(pList[globalIndex]+5) = middle;
447 *(pList[globalIndex]+4) = highest;
448 *(pList[globalIndex]+3) = signal;
449 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
450 *(pList[globalIndex]+1) = *pList[globalIndex];
451 *pList[globalIndex] = label;
452 }else if (signal>middle){
453 *(pList[globalIndex]+5) = middle;
454 *(pList[globalIndex]+4) = signal;
455 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
456 *(pList[globalIndex]+1) = label;
457 }else{
458 *(pList[globalIndex]+5) = signal;
459 *(pList[globalIndex]+2) = label;
460 } // end if
461 } // end if
462 } // end of loop pixels in x
463 } // end of loop over pixels in z
fd61217e 464}
0315d466 465//----------------------------------------------------------------------
fd61217e 466void AliITSsimulationSSD::ChargeToSignal(Float_t **pList) {
0315d466 467 // charge to signal
468 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
469 Float_t threshold = 0.;
470 Int_t digits[3], tracks[3],hits[3],gi,j1;
471 Float_t charges[3];
472 Float_t signal,phys;
473 Float_t noise[2] = {0.,0.};
474
475 fResponse->GetNoiseParam(noise[0],noise[1]);
fd61217e 476
0315d466 477 for(Int_t k=0;k<2;k++){ // both sides (0=Pside, 1=Nside)
478
479 // Threshold for zero-suppression
480 // It can be defined in AliITSresponseSSD
481 // threshold = (Float_t)fResponse->MinVal(k);
482 // I prefer to think adjusting the threshold "manually", looking
483 // at the scope, and considering noise standard deviation
484 threshold = 4.0*noise[k]; // 4 times noise is a choice
485 for(Int_t ix=0;ix<fNstrips;ix++){ // loop over strips
486
487 signal = (Float_t) fMapA2->GetSignal(k,ix);
488
489 gi =k*fNstrips+ix; // global index
490 if (signal > threshold) {
491 digits[0]=k;
492 digits[1]=ix;
493
494 // convert to ADC signal
495 // conversion factor are rather arbitrary (need tuning)
496 // minimum ionizing particle--> ~30000 pairs--> ADC channel 50
497 signal = signal*50.0/30000.0;
498 if(signal>1000.) signal = 1000.0;//if exceeding, accumulate
499 // last one
500 digits[2]=(Int_t) signal;
501
502 //gi =k*fNstrips+ix; // global index
503 for(j1=0;j1<3;j1++){
504 if (pList[gi]) {
505 tracks[j1] = (Int_t)(*(pList[gi]+j1));
506 } else {
507 tracks[j1]=-2; //noise
508 } // end if pList
509 charges[j1] = 0;
510 } // end for j1
511
512 phys=0;
513
514 hits[0]=0;
515 hits[1]=0;
516 hits[2]=0;
517 // finally add digit
518 aliITS->AddSimDigit(2,phys,digits,tracks,hits,charges);
519
520 //if(pList[gi]) delete [] pList[gi];
521 } // end if signal > threshold
522 if(pList[gi]) delete [] pList[gi];
523 } // end for ix
524 } // end for k
525 delete [] pList;
b0f5e3fc 526}