f7336fa3 |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
3 | * * |
4 | * Author: The ALICE Off-line Project. * |
5 | * Contributors are mentioned in the code where appropriate. * |
6 | * * |
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. The authors make no claims * |
12 | * about the suitability of this software for any purpose. It is * |
13 | * provided "as is" without express or implied warranty. * |
14 | **************************************************************************/ |
15 | |
16 | /* |
17 | $Log$ |
3e1a3ad8 |
18 | Revision 1.9 2000/11/01 14:53:20 cblume |
19 | Merge with TRD-develop |
20 | |
793ff80c |
21 | Revision 1.1.4.5 2000/10/15 23:40:01 cblume |
22 | Remove AliTRDconst |
23 | |
24 | Revision 1.1.4.4 2000/10/06 16:49:46 cblume |
25 | Made Getters const |
26 | |
27 | Revision 1.1.4.3 2000/10/04 16:34:58 cblume |
28 | Replace include files by forward declarations |
29 | |
30 | Revision 1.1.4.2 2000/09/22 14:49:49 cblume |
31 | Adapted to tracking code |
32 | |
33 | Revision 1.8 2000/10/02 21:28:19 fca |
34 | Removal of useless dependecies via forward declarations |
35 | |
94de3818 |
36 | Revision 1.7 2000/06/27 13:08:50 cblume |
37 | Changed to Copy(TObject &A) to appease the HP-compiler |
38 | |
43da34c0 |
39 | Revision 1.6 2000/06/09 11:10:07 cblume |
40 | Compiler warnings and coding conventions, next round |
41 | |
dd9a6ee3 |
42 | Revision 1.5 2000/06/08 18:32:58 cblume |
43 | Make code compliant to coding conventions |
44 | |
8230f242 |
45 | Revision 1.4 2000/06/07 16:27:01 cblume |
46 | Try to remove compiler warnings on Sun and HP |
47 | |
9d0b222b |
48 | Revision 1.3 2000/05/08 16:17:27 cblume |
49 | Merge TRD-develop |
50 | |
6f1e466d |
51 | Revision 1.1.4.1 2000/05/08 15:09:01 cblume |
52 | Introduce AliTRDdigitsManager |
53 | |
c0dd96c3 |
54 | Revision 1.1 2000/02/28 18:58:54 cblume |
55 | Add new TRD classes |
56 | |
f7336fa3 |
57 | */ |
58 | |
59 | /////////////////////////////////////////////////////////////////////////////// |
60 | // // |
61 | // TRD cluster finder for the slow simulator. |
62 | // // |
63 | /////////////////////////////////////////////////////////////////////////////// |
64 | |
65 | #include <TF1.h> |
94de3818 |
66 | #include <TTree.h> |
793ff80c |
67 | #include <TH1.h> |
f7336fa3 |
68 | |
793ff80c |
69 | #include "AliRun.h" |
70 | |
71 | #include "AliTRD.h" |
f7336fa3 |
72 | #include "AliTRDclusterizerV1.h" |
73 | #include "AliTRDmatrix.h" |
74 | #include "AliTRDgeometry.h" |
75 | #include "AliTRDdigitizer.h" |
6f1e466d |
76 | #include "AliTRDdataArrayF.h" |
793ff80c |
77 | #include "AliTRDdataArrayI.h" |
78 | #include "AliTRDdigitsManager.h" |
f7336fa3 |
79 | |
80 | ClassImp(AliTRDclusterizerV1) |
81 | |
82 | //_____________________________________________________________________________ |
83 | AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer() |
84 | { |
85 | // |
86 | // AliTRDclusterizerV1 default constructor |
87 | // |
88 | |
6f1e466d |
89 | fDigitsManager = NULL; |
f7336fa3 |
90 | |
3e1a3ad8 |
91 | fClusMaxThresh = 0; |
92 | fClusSigThresh = 0; |
93 | |
f7336fa3 |
94 | } |
95 | |
96 | //_____________________________________________________________________________ |
97 | AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title) |
98 | :AliTRDclusterizer(name,title) |
99 | { |
100 | // |
101 | // AliTRDclusterizerV1 default constructor |
102 | // |
103 | |
6f1e466d |
104 | fDigitsManager = new AliTRDdigitsManager(); |
f7336fa3 |
105 | |
106 | Init(); |
107 | |
108 | } |
109 | |
8230f242 |
110 | //_____________________________________________________________________________ |
dd9a6ee3 |
111 | AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c) |
8230f242 |
112 | { |
113 | // |
114 | // AliTRDclusterizerV1 copy constructor |
115 | // |
116 | |
dd9a6ee3 |
117 | ((AliTRDclusterizerV1 &) c).Copy(*this); |
8230f242 |
118 | |
119 | } |
120 | |
f7336fa3 |
121 | //_____________________________________________________________________________ |
122 | AliTRDclusterizerV1::~AliTRDclusterizerV1() |
123 | { |
8230f242 |
124 | // |
125 | // AliTRDclusterizerV1 destructor |
126 | // |
f7336fa3 |
127 | |
6f1e466d |
128 | if (fDigitsManager) { |
129 | delete fDigitsManager; |
f7336fa3 |
130 | } |
131 | |
132 | } |
133 | |
dd9a6ee3 |
134 | //_____________________________________________________________________________ |
135 | AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c) |
136 | { |
137 | // |
138 | // Assignment operator |
139 | // |
140 | |
141 | if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this); |
142 | return *this; |
143 | |
144 | } |
145 | |
8230f242 |
146 | //_____________________________________________________________________________ |
43da34c0 |
147 | void AliTRDclusterizerV1::Copy(TObject &c) |
8230f242 |
148 | { |
149 | // |
150 | // Copy function |
151 | // |
152 | |
43da34c0 |
153 | ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh; |
154 | ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh; |
43da34c0 |
155 | ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL; |
8230f242 |
156 | |
157 | AliTRDclusterizer::Copy(c); |
158 | |
159 | } |
160 | |
f7336fa3 |
161 | //_____________________________________________________________________________ |
162 | void AliTRDclusterizerV1::Init() |
163 | { |
164 | // |
165 | // Initializes the cluster finder |
166 | // |
167 | |
168 | // The default parameter for the clustering |
3e1a3ad8 |
169 | fClusMaxThresh = 3; |
170 | fClusSigThresh = 1; |
f7336fa3 |
171 | |
172 | } |
173 | |
174 | //_____________________________________________________________________________ |
175 | Bool_t AliTRDclusterizerV1::ReadDigits() |
176 | { |
177 | // |
178 | // Reads the digits arrays from the input aliroot file |
179 | // |
180 | |
181 | if (!fInputFile) { |
182 | printf("AliTRDclusterizerV1::ReadDigits -- "); |
183 | printf("No input file open\n"); |
184 | return kFALSE; |
185 | } |
186 | |
f7336fa3 |
187 | // Read in the digit arrays |
6f1e466d |
188 | return (fDigitsManager->ReadDigits()); |
f7336fa3 |
189 | |
190 | } |
191 | |
192 | //_____________________________________________________________________________ |
793ff80c |
193 | Bool_t AliTRDclusterizerV1::MakeClusters() |
f7336fa3 |
194 | { |
195 | // |
196 | // Generates the cluster. |
197 | // |
198 | |
199 | Int_t row, col, time; |
200 | |
3e1a3ad8 |
201 | if (fTRD->IsVersion() != 1) { |
f7336fa3 |
202 | printf("AliTRDclusterizerV1::MakeCluster -- "); |
203 | printf("TRD must be version 1 (slow simulator).\n"); |
204 | return kFALSE; |
205 | } |
206 | |
207 | // Get the geometry |
3e1a3ad8 |
208 | AliTRDgeometry *geo = fTRD->GetGeometry(); |
f7336fa3 |
209 | |
210 | printf("AliTRDclusterizerV1::MakeCluster -- "); |
211 | printf("Start creating clusters.\n"); |
212 | |
8230f242 |
213 | AliTRDdataArrayI *digits; |
793ff80c |
214 | AliTRDdataArrayI *track0; |
215 | AliTRDdataArrayI *track1; |
216 | AliTRDdataArrayI *track2; |
f7336fa3 |
217 | |
3e1a3ad8 |
218 | // Threshold value for the maximum |
219 | Int_t maxThresh = fClusMaxThresh; |
220 | // Threshold value for the digit signal |
221 | Int_t sigThresh = fClusSigThresh; |
f7336fa3 |
222 | |
223 | // Iteration limit for unfolding procedure |
8230f242 |
224 | const Float_t kEpsilon = 0.01; |
f7336fa3 |
225 | |
8230f242 |
226 | const Int_t kNclus = 3; |
227 | const Int_t kNsig = 5; |
3e1a3ad8 |
228 | const Int_t kNtrack = 3 * kNclus; |
229 | |
230 | Float_t padSignal[kNsig]; |
231 | Float_t clusterSignal[kNclus]; |
232 | Float_t clusterPads[kNclus]; |
233 | Int_t clusterDigit[kNclus]; |
234 | Int_t clusterTracks[kNtrack]; |
f7336fa3 |
235 | |
236 | Int_t chamBeg = 0; |
793ff80c |
237 | Int_t chamEnd = AliTRDgeometry::Ncham(); |
3e1a3ad8 |
238 | if (fTRD->GetSensChamber() >= 0) { |
239 | chamBeg = fTRD->GetSensChamber(); |
6f1e466d |
240 | chamEnd = chamBeg + 1; |
f7336fa3 |
241 | } |
242 | Int_t planBeg = 0; |
793ff80c |
243 | Int_t planEnd = AliTRDgeometry::Nplan(); |
3e1a3ad8 |
244 | if (fTRD->GetSensPlane() >= 0) { |
245 | planBeg = fTRD->GetSensPlane(); |
f7336fa3 |
246 | planEnd = planBeg + 1; |
247 | } |
248 | Int_t sectBeg = 0; |
793ff80c |
249 | Int_t sectEnd = AliTRDgeometry::Nsect(); |
f7336fa3 |
250 | |
3e1a3ad8 |
251 | // Start clustering in every chamber |
f7336fa3 |
252 | for (Int_t icham = chamBeg; icham < chamEnd; icham++) { |
253 | for (Int_t iplan = planBeg; iplan < planEnd; iplan++) { |
254 | for (Int_t isect = sectBeg; isect < sectEnd; isect++) { |
255 | |
3e1a3ad8 |
256 | if (fTRD->GetSensSector() >= 0) { |
257 | Int_t sens1 = fTRD->GetSensSector(); |
258 | Int_t sens2 = sens1 + fTRD->GetSensSectorRange(); |
793ff80c |
259 | sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) |
260 | * AliTRDgeometry::Nsect(); |
dd9a6ee3 |
261 | if (sens1 < sens2) { |
9d0b222b |
262 | if ((isect < sens1) || (isect >= sens2)) continue; |
dd9a6ee3 |
263 | } |
264 | else { |
9d0b222b |
265 | if ((isect < sens1) && (isect >= sens2)) continue; |
dd9a6ee3 |
266 | } |
9d0b222b |
267 | } |
268 | |
8230f242 |
269 | Int_t idet = geo->GetDetector(iplan,icham,isect); |
f7336fa3 |
270 | |
3e1a3ad8 |
271 | Int_t nClusters = 0; |
272 | Int_t nClustersUnfold = 0; |
273 | |
f7336fa3 |
274 | printf("AliTRDclusterizerV1::MakeCluster -- "); |
275 | printf("Analyzing chamber %d, plane %d, sector %d.\n" |
3e1a3ad8 |
276 | ,icham,iplan,isect); |
f7336fa3 |
277 | |
3e1a3ad8 |
278 | Int_t nRowMax = geo->GetRowMax(iplan,icham,isect); |
279 | Int_t nColMax = geo->GetColMax(iplan); |
280 | Int_t nTimeBefore = geo->GetTimeBefore(); |
281 | Int_t nTimeTotal = geo->GetTimeTotal(); |
f7336fa3 |
282 | |
3e1a3ad8 |
283 | // Get the digits |
8230f242 |
284 | digits = fDigitsManager->GetDigits(idet); |
3e1a3ad8 |
285 | digits->Expand(); |
793ff80c |
286 | track0 = fDigitsManager->GetDictionary(idet,0); |
3e1a3ad8 |
287 | track0->Expand(); |
793ff80c |
288 | track1 = fDigitsManager->GetDictionary(idet,1); |
3e1a3ad8 |
289 | track1->Expand(); |
793ff80c |
290 | track2 = fDigitsManager->GetDictionary(idet,2); |
3e1a3ad8 |
291 | track2->Expand(); |
292 | |
293 | // Loop through the chamber and find the maxima |
294 | for ( row = 0; row < nRowMax; row++) { |
295 | for ( col = 2; col < nColMax; col++) { |
296 | for (time = 0; time < nTimeTotal; time++) { |
297 | |
298 | Int_t signalL = digits->GetDataUnchecked(row,col ,time); |
299 | Int_t signalM = digits->GetDataUnchecked(row,col-1,time); |
300 | Int_t signalR = digits->GetDataUnchecked(row,col-2,time); |
301 | |
302 | // Look for the maximum |
303 | if ((signalM >= maxThresh) && |
304 | (signalL >= sigThresh) && |
305 | (signalR >= sigThresh)) { |
306 | if ((signalL < signalM) && (signalR < signalM)) { |
307 | // Maximum found, mark the position by a negative signal |
308 | digits->SetDataUnchecked(row,col-1,time,-signalM); |
309 | } |
310 | } |
311 | |
312 | } |
313 | } |
314 | } |
315 | |
316 | // Now check the maxima and calculate the cluster position |
317 | for ( row = 0; row < nRowMax ; row++) { |
318 | for ( col = 1; col < nColMax-1; col++) { |
319 | for (time = 0; time < nTimeTotal; time++) { |
320 | |
321 | // Maximum found ? |
322 | if (digits->GetDataUnchecked(row,col,time) < 0) { |
f7336fa3 |
323 | |
9d0b222b |
324 | Int_t iPad; |
8230f242 |
325 | for (iPad = 0; iPad < kNclus; iPad++) { |
3e1a3ad8 |
326 | Int_t iPadCol = col - 1 + iPad; |
327 | clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row |
328 | ,iPadCol |
329 | ,time)); |
330 | clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time); |
331 | clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1; |
332 | clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1; |
333 | clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1; |
f7336fa3 |
334 | } |
335 | |
3e1a3ad8 |
336 | // Neighbouring maximum on right side? |
337 | Int_t iType = 0; |
338 | if (col < nColMax-3) { |
339 | if (digits->GetDataUnchecked(row,col+2,time) < 0) { |
340 | for (iPad = 0; iPad < kNsig; iPad++) { |
341 | padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row |
342 | ,col-1+iPad |
343 | ,time)); |
f7336fa3 |
344 | } |
3e1a3ad8 |
345 | // Unfold the two maxima and set the signal on |
346 | // the overlapping pad to the ratio |
347 | clusterSignal[2] *= Unfold(kEpsilon,padSignal); |
348 | iType = 1; |
f7336fa3 |
349 | } |
350 | } |
f7336fa3 |
351 | |
3e1a3ad8 |
352 | Float_t clusterCharge = clusterSignal[0] |
353 | + clusterSignal[1] |
354 | + clusterSignal[2]; |
355 | |
356 | // Calculate the position of the cluster by using the |
357 | // center of gravity method |
358 | clusterPads[0] = row + 0.5; |
359 | clusterPads[1] = col + 0.5 |
360 | + (clusterSignal[2] - clusterSignal[0]) |
361 | / clusterCharge; |
362 | // Take the shift of the additional time bins into account |
363 | clusterPads[2] = time - nTimeBefore + 0.5; |
364 | |
365 | if (fVerbose) { |
366 | printf("-----------------------------------------------------------\n"); |
367 | printf("Create cluster no. %d\n",nClusters); |
368 | printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0] |
369 | ,clusterPads[1] |
370 | ,clusterPads[2]); |
371 | printf("Indices: %d, %d, %d\n",clusterDigit[0] |
372 | ,clusterDigit[1] |
373 | ,clusterDigit[2]); |
374 | printf("Total charge = %f\n",clusterCharge); |
375 | printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0] |
376 | ,clusterTracks[1] |
377 | ,clusterTracks[2]); |
378 | printf(" pad1 %d, %d, %d\n",clusterTracks[3] |
379 | ,clusterTracks[4] |
380 | ,clusterTracks[5]); |
381 | printf(" pad2 %d, %d, %d\n",clusterTracks[6] |
382 | ,clusterTracks[7] |
383 | ,clusterTracks[8]); |
f7336fa3 |
384 | } |
385 | |
3e1a3ad8 |
386 | nClusters++; |
387 | if (iType == 1) nClustersUnfold++; |
f7336fa3 |
388 | |
389 | // Add the cluster to the output array |
3e1a3ad8 |
390 | fTRD->AddCluster(clusterPads |
391 | ,clusterDigit |
392 | ,idet |
393 | ,clusterCharge |
394 | ,clusterTracks |
395 | ,iType); |
f7336fa3 |
396 | |
397 | } |
3e1a3ad8 |
398 | } |
399 | } |
400 | } |
f7336fa3 |
401 | |
3e1a3ad8 |
402 | // Compress the arrays |
403 | digits->Compress(1,0); |
404 | track0->Compress(1,0); |
405 | track1->Compress(1,0); |
406 | track2->Compress(1,0); |
f7336fa3 |
407 | |
3e1a3ad8 |
408 | // Write the cluster and reset the array |
793ff80c |
409 | WriteClusters(idet); |
3e1a3ad8 |
410 | fTRD->ResetRecPoints(); |
793ff80c |
411 | |
3e1a3ad8 |
412 | printf("AliTRDclusterizerV1::MakeCluster -- "); |
413 | printf("Found %d (%d unfolded) clusters.\n" |
414 | ,nClusters,nClustersUnfold); |
f7336fa3 |
415 | |
3e1a3ad8 |
416 | } |
417 | } |
418 | } |
f7336fa3 |
419 | |
f7336fa3 |
420 | printf("AliTRDclusterizerV1::MakeCluster -- "); |
421 | printf("Done.\n"); |
422 | |
423 | return kTRUE; |
424 | |
425 | } |
426 | |
427 | //_____________________________________________________________________________ |
428 | Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal) |
429 | { |
430 | // |
431 | // Method to unfold neighbouring maxima. |
432 | // The charge ratio on the overlapping pad is calculated |
433 | // until there is no more change within the range given by eps. |
434 | // The resulting ratio is then returned to the calling method. |
435 | // |
436 | |
3e1a3ad8 |
437 | Int_t itStep = 0; // Count iteration steps |
f7336fa3 |
438 | |
3e1a3ad8 |
439 | Float_t ratio = 0.5; // Start value for ratio |
440 | Float_t prevRatio = 0; // Store previous ratio |
f7336fa3 |
441 | |
3e1a3ad8 |
442 | Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal |
443 | Float_t newRightSignal[3] = {0}; // Array to store right cluster signal |
f7336fa3 |
444 | |
3e1a3ad8 |
445 | // Start the iteration |
f7336fa3 |
446 | while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { |
447 | |
448 | itStep++; |
449 | prevRatio = ratio; |
450 | |
3e1a3ad8 |
451 | // Cluster position according to charge ratio |
452 | Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) |
453 | / (padSignal[0] + padSignal[1] + ratio*padSignal[2]); |
454 | Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) |
455 | / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]); |
f7336fa3 |
456 | |
3e1a3ad8 |
457 | // Set cluster charge ratio |
f7336fa3 |
458 | Float_t ampLeft = padSignal[1]; |
459 | Float_t ampRight = padSignal[3]; |
460 | |
3e1a3ad8 |
461 | // Apply pad response to parameters |
462 | newLeftSignal[0] = ampLeft * PadResponse(-1 - maxLeft); |
463 | newLeftSignal[1] = ampLeft * PadResponse( 0 - maxLeft); |
464 | newLeftSignal[2] = ampLeft * PadResponse( 1 - maxLeft); |
f7336fa3 |
465 | |
3e1a3ad8 |
466 | newRightSignal[0] = ampRight * PadResponse(-1 - maxRight); |
467 | newRightSignal[1] = ampRight * PadResponse( 0 - maxRight); |
468 | newRightSignal[2] = ampRight * PadResponse( 1 - maxRight); |
f7336fa3 |
469 | |
3e1a3ad8 |
470 | // Calculate new overlapping ratio |
471 | ratio = newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]); |
f7336fa3 |
472 | |
473 | } |
474 | |
475 | return ratio; |
476 | |
477 | } |
478 | |
479 | //_____________________________________________________________________________ |
480 | Float_t AliTRDclusterizerV1::PadResponse(Float_t x) |
481 | { |
482 | // |
483 | // The pad response for the chevron pads. |
484 | // We use a simple Gaussian approximation which should be good |
485 | // enough for our purpose. |
3e1a3ad8 |
486 | // Updated for new PRF 1/5/01. |
f7336fa3 |
487 | // |
488 | |
489 | // The parameters for the response function |
3e1a3ad8 |
490 | const Float_t kA = 0.8303; |
491 | const Float_t kB = -0.00392; |
492 | const Float_t kC = 0.472 * 0.472; |
493 | const Float_t kD = 2.19; |
f7336fa3 |
494 | |
3e1a3ad8 |
495 | Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC))); |
f7336fa3 |
496 | |
497 | return (pr); |
498 | |
499 | } |