Introduce parameter class and regions of interest for merging
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
CommitLineData
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$
47517f42 18Revision 1.14 2001/11/14 10:50:45 cblume
19Changes in digits IO. Add merging of summable digits
20
abaf1f1d 21Revision 1.13 2001/05/28 17:07:58 hristov
22Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
23
a819a5f7 24Revision 1.12 2001/05/21 17:42:58 hristov
25Constant casted to avoid the ambiguity
26
26edf6a4 27Revision 1.11 2001/05/21 16:45:47 hristov
28Last minute changes (C.Blume)
29
db30bf0f 30Revision 1.10 2001/05/07 08:06:44 cblume
31Speedup of the code. Create only AliTRDcluster
32
3e1a3ad8 33Revision 1.9 2000/11/01 14:53:20 cblume
34Merge with TRD-develop
35
793ff80c 36Revision 1.1.4.5 2000/10/15 23:40:01 cblume
37Remove AliTRDconst
38
39Revision 1.1.4.4 2000/10/06 16:49:46 cblume
40Made Getters const
41
42Revision 1.1.4.3 2000/10/04 16:34:58 cblume
43Replace include files by forward declarations
44
45Revision 1.1.4.2 2000/09/22 14:49:49 cblume
46Adapted to tracking code
47
48Revision 1.8 2000/10/02 21:28:19 fca
49Removal of useless dependecies via forward declarations
50
94de3818 51Revision 1.7 2000/06/27 13:08:50 cblume
52Changed to Copy(TObject &A) to appease the HP-compiler
53
43da34c0 54Revision 1.6 2000/06/09 11:10:07 cblume
55Compiler warnings and coding conventions, next round
56
dd9a6ee3 57Revision 1.5 2000/06/08 18:32:58 cblume
58Make code compliant to coding conventions
59
8230f242 60Revision 1.4 2000/06/07 16:27:01 cblume
61Try to remove compiler warnings on Sun and HP
62
9d0b222b 63Revision 1.3 2000/05/08 16:17:27 cblume
64Merge TRD-develop
65
6f1e466d 66Revision 1.1.4.1 2000/05/08 15:09:01 cblume
67Introduce AliTRDdigitsManager
68
c0dd96c3 69Revision 1.1 2000/02/28 18:58:54 cblume
70Add new TRD classes
71
f7336fa3 72*/
73
74///////////////////////////////////////////////////////////////////////////////
75// //
76// TRD cluster finder for the slow simulator.
77// //
78///////////////////////////////////////////////////////////////////////////////
79
80#include <TF1.h>
94de3818 81#include <TTree.h>
793ff80c 82#include <TH1.h>
a819a5f7 83#include <TFile.h>
f7336fa3 84
793ff80c 85#include "AliRun.h"
86
87#include "AliTRD.h"
f7336fa3 88#include "AliTRDclusterizerV1.h"
89#include "AliTRDmatrix.h"
90#include "AliTRDgeometry.h"
91#include "AliTRDdigitizer.h"
6f1e466d 92#include "AliTRDdataArrayF.h"
793ff80c 93#include "AliTRDdataArrayI.h"
94#include "AliTRDdigitsManager.h"
f7336fa3 95
96ClassImp(AliTRDclusterizerV1)
97
98//_____________________________________________________________________________
99AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
100{
101 //
102 // AliTRDclusterizerV1 default constructor
103 //
104
6f1e466d 105 fDigitsManager = NULL;
f7336fa3 106
3e1a3ad8 107 fClusMaxThresh = 0;
108 fClusSigThresh = 0;
109
db30bf0f 110 fUseLUT = kFALSE;
111
f7336fa3 112}
113
114//_____________________________________________________________________________
115AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
116 :AliTRDclusterizer(name,title)
117{
118 //
119 // AliTRDclusterizerV1 default constructor
120 //
121
6f1e466d 122 fDigitsManager = new AliTRDdigitsManager();
f7336fa3 123
124 Init();
125
126}
127
8230f242 128//_____________________________________________________________________________
dd9a6ee3 129AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
8230f242 130{
131 //
132 // AliTRDclusterizerV1 copy constructor
133 //
134
dd9a6ee3 135 ((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242 136
137}
138
f7336fa3 139//_____________________________________________________________________________
140AliTRDclusterizerV1::~AliTRDclusterizerV1()
141{
8230f242 142 //
143 // AliTRDclusterizerV1 destructor
144 //
f7336fa3 145
6f1e466d 146 if (fDigitsManager) {
147 delete fDigitsManager;
abaf1f1d 148 fDigitsManager = NULL;
f7336fa3 149 }
150
151}
152
dd9a6ee3 153//_____________________________________________________________________________
154AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
155{
156 //
157 // Assignment operator
158 //
159
160 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
161 return *this;
162
163}
164
8230f242 165//_____________________________________________________________________________
43da34c0 166void AliTRDclusterizerV1::Copy(TObject &c)
8230f242 167{
168 //
169 // Copy function
170 //
171
db30bf0f 172 ((AliTRDclusterizerV1 &) c).fUseLUT = fUseLUT;
43da34c0 173 ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
174 ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
43da34c0 175 ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
db30bf0f 176 for (Int_t ilut = 0; ilut < kNlut; ilut++) {
177 ((AliTRDclusterizerV1 &) c).fLUT[ilut] = fLUT[ilut];
178 }
8230f242 179
180 AliTRDclusterizer::Copy(c);
181
182}
183
f7336fa3 184//_____________________________________________________________________________
185void AliTRDclusterizerV1::Init()
186{
187 //
188 // Initializes the cluster finder
189 //
190
191 // The default parameter for the clustering
3e1a3ad8 192 fClusMaxThresh = 3;
193 fClusSigThresh = 1;
f7336fa3 194
db30bf0f 195 // Use the lookup table for the position determination
196 fUseLUT = kTRUE;
197
198 // The lookup table from Bogdan
199 Float_t lut[128] = {
200 0.0068, 0.0198, 0.0318, 0.0432, 0.0538, 0.0642, 0.0742, 0.0838,
201 0.0932, 0.1023, 0.1107, 0.1187, 0.1268, 0.1347, 0.1423, 0.1493,
202 0.1562, 0.1632, 0.1698, 0.1762, 0.1828, 0.1887, 0.1947, 0.2002,
203 0.2062, 0.2118, 0.2173, 0.2222, 0.2278, 0.2327, 0.2377, 0.2428,
204 0.2473, 0.2522, 0.2567, 0.2612, 0.2657, 0.2697, 0.2743, 0.2783,
205 0.2822, 0.2862, 0.2903, 0.2943, 0.2982, 0.3018, 0.3058, 0.3092,
206 0.3128, 0.3167, 0.3203, 0.3237, 0.3268, 0.3302, 0.3338, 0.3368,
207 0.3402, 0.3433, 0.3462, 0.3492, 0.3528, 0.3557, 0.3587, 0.3613,
208 0.3643, 0.3672, 0.3702, 0.3728, 0.3758, 0.3783, 0.3812, 0.3837,
209 0.3862, 0.3887, 0.3918, 0.3943, 0.3968, 0.3993, 0.4017, 0.4042,
210 0.4067, 0.4087, 0.4112, 0.4137, 0.4157, 0.4182, 0.4207, 0.4227,
211 0.4252, 0.4272, 0.4293, 0.4317, 0.4338, 0.4358, 0.4383, 0.4403,
212 0.4423, 0.4442, 0.4462, 0.4482, 0.4502, 0.4523, 0.4543, 0.4563,
213 0.4582, 0.4602, 0.4622, 0.4638, 0.4658, 0.4678, 0.4697, 0.4712,
214 0.4733, 0.4753, 0.4767, 0.4787, 0.4803, 0.4823, 0.4837, 0.4857,
215 0.4873, 0.4888, 0.4908, 0.4922, 0.4942, 0.4958, 0.4972, 0.4988
216 };
217 for (Int_t ilut = 0; ilut < kNlut; ilut++) {
218 fLUT[ilut] = lut[ilut];
219 }
220
abaf1f1d 221 fDigitsManager->CreateArrays();
222
f7336fa3 223}
224
225//_____________________________________________________________________________
226Bool_t AliTRDclusterizerV1::ReadDigits()
227{
228 //
229 // Reads the digits arrays from the input aliroot file
230 //
231
232 if (!fInputFile) {
233 printf("AliTRDclusterizerV1::ReadDigits -- ");
234 printf("No input file open\n");
235 return kFALSE;
236 }
237
abaf1f1d 238 fDigitsManager->Open(fInputFile->GetName());
239
f7336fa3 240 // Read in the digit arrays
6f1e466d 241 return (fDigitsManager->ReadDigits());
f7336fa3 242
243}
244
245//_____________________________________________________________________________
793ff80c 246Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3 247{
248 //
249 // Generates the cluster.
250 //
251
252 Int_t row, col, time;
253
3e1a3ad8 254 if (fTRD->IsVersion() != 1) {
f7336fa3 255 printf("AliTRDclusterizerV1::MakeCluster -- ");
256 printf("TRD must be version 1 (slow simulator).\n");
257 return kFALSE;
258 }
259
260 // Get the geometry
3e1a3ad8 261 AliTRDgeometry *geo = fTRD->GetGeometry();
f7336fa3 262
a819a5f7 263 Float_t timeBinSize = geo->GetTimeBinSize();
264 // Half of ampl.region
265 const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
266
47517f42 267 AliTRDdigitizer *digitizer = (AliTRDdigitizer*) fInputFile->Get("TRDdigitizer");
a819a5f7 268 Float_t omegaTau = digitizer->GetOmegaTau();
47517f42 269 if (fVerbose > 0) {
270 printf("AliTRDclusterizerV1::MakeCluster -- ");
271 printf("OmegaTau = %f \n",omegaTau);
272 printf("AliTRDclusterizerV1::MakeCluster -- ");
273 printf("Start creating clusters.\n");
274 }
f7336fa3 275
8230f242 276 AliTRDdataArrayI *digits;
793ff80c 277 AliTRDdataArrayI *track0;
278 AliTRDdataArrayI *track1;
279 AliTRDdataArrayI *track2;
f7336fa3 280
3e1a3ad8 281 // Threshold value for the maximum
282 Int_t maxThresh = fClusMaxThresh;
283 // Threshold value for the digit signal
284 Int_t sigThresh = fClusSigThresh;
f7336fa3 285
286 // Iteration limit for unfolding procedure
8230f242 287 const Float_t kEpsilon = 0.01;
f7336fa3 288
8230f242 289 const Int_t kNclus = 3;
290 const Int_t kNsig = 5;
3e1a3ad8 291 const Int_t kNtrack = 3 * kNclus;
292
db30bf0f 293 // For the LUT
294 const Float_t kLUTmin = 0.106113;
295 const Float_t kLUTmax = 0.995415;
296
297 Int_t iType = 0;
298 Int_t iUnfold = 0;
299
300 Float_t ratioLeft = 1.0;
301 Float_t ratioRight = 1.0;
302
3e1a3ad8 303 Float_t padSignal[kNsig];
304 Float_t clusterSignal[kNclus];
305 Float_t clusterPads[kNclus];
306 Int_t clusterDigit[kNclus];
307 Int_t clusterTracks[kNtrack];
f7336fa3 308
309 Int_t chamBeg = 0;
793ff80c 310 Int_t chamEnd = AliTRDgeometry::Ncham();
3e1a3ad8 311 if (fTRD->GetSensChamber() >= 0) {
312 chamBeg = fTRD->GetSensChamber();
6f1e466d 313 chamEnd = chamBeg + 1;
f7336fa3 314 }
315 Int_t planBeg = 0;
793ff80c 316 Int_t planEnd = AliTRDgeometry::Nplan();
3e1a3ad8 317 if (fTRD->GetSensPlane() >= 0) {
318 planBeg = fTRD->GetSensPlane();
f7336fa3 319 planEnd = planBeg + 1;
320 }
321 Int_t sectBeg = 0;
793ff80c 322 Int_t sectEnd = AliTRDgeometry::Nsect();
f7336fa3 323
3e1a3ad8 324 // Start clustering in every chamber
f7336fa3 325 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
326 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
327 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
328
3e1a3ad8 329 if (fTRD->GetSensSector() >= 0) {
330 Int_t sens1 = fTRD->GetSensSector();
331 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
793ff80c 332 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
333 * AliTRDgeometry::Nsect();
dd9a6ee3 334 if (sens1 < sens2) {
9d0b222b 335 if ((isect < sens1) || (isect >= sens2)) continue;
dd9a6ee3 336 }
337 else {
9d0b222b 338 if ((isect < sens1) && (isect >= sens2)) continue;
dd9a6ee3 339 }
9d0b222b 340 }
341
8230f242 342 Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3 343
db30bf0f 344 Int_t nClusters = 0;
345 Int_t nClusters2pad = 0;
346 Int_t nClusters3pad = 0;
347 Int_t nClusters4pad = 0;
348 Int_t nClusters5pad = 0;
349 Int_t nClustersLarge = 0;
3e1a3ad8 350
47517f42 351 if (fVerbose > 0) {
352 printf("AliTRDclusterizerV1::MakeCluster -- ");
353 printf("Analyzing chamber %d, plane %d, sector %d.\n"
354 ,icham,iplan,isect);
355 }
f7336fa3 356
3e1a3ad8 357 Int_t nRowMax = geo->GetRowMax(iplan,icham,isect);
358 Int_t nColMax = geo->GetColMax(iplan);
359 Int_t nTimeBefore = geo->GetTimeBefore();
360 Int_t nTimeTotal = geo->GetTimeTotal();
f7336fa3 361
3e1a3ad8 362 // Get the digits
8230f242 363 digits = fDigitsManager->GetDigits(idet);
3e1a3ad8 364 digits->Expand();
793ff80c 365 track0 = fDigitsManager->GetDictionary(idet,0);
3e1a3ad8 366 track0->Expand();
793ff80c 367 track1 = fDigitsManager->GetDictionary(idet,1);
3e1a3ad8 368 track1->Expand();
793ff80c 369 track2 = fDigitsManager->GetDictionary(idet,2);
3e1a3ad8 370 track2->Expand();
371
372 // Loop through the chamber and find the maxima
373 for ( row = 0; row < nRowMax; row++) {
374 for ( col = 2; col < nColMax; col++) {
375 for (time = 0; time < nTimeTotal; time++) {
376
a819a5f7 377 Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
378 Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
379 Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
3e1a3ad8 380
381 // Look for the maximum
db30bf0f 382 if (signalM >= maxThresh) {
383 if (((signalL >= sigThresh) &&
384 (signalL < signalM)) ||
385 ((signalR >= sigThresh) &&
386 (signalR < signalM))) {
3e1a3ad8 387 // Maximum found, mark the position by a negative signal
388 digits->SetDataUnchecked(row,col-1,time,-signalM);
389 }
390 }
391
392 }
393 }
394 }
395
396 // Now check the maxima and calculate the cluster position
397 for ( row = 0; row < nRowMax ; row++) {
db30bf0f 398 for (time = 0; time < nTimeTotal; time++) {
399 for ( col = 1; col < nColMax-1; col++) {
3e1a3ad8 400
401 // Maximum found ?
402 if (digits->GetDataUnchecked(row,col,time) < 0) {
f7336fa3 403
9d0b222b 404 Int_t iPad;
8230f242 405 for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8 406 Int_t iPadCol = col - 1 + iPad;
407 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
408 ,iPadCol
409 ,time));
410 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
411 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
412 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
413 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
f7336fa3 414 }
415
db30bf0f 416 // Count the number of pads in the cluster
417 Int_t nPadCount = 0;
418 Int_t ii = 0;
419 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
420 >= sigThresh) {
421 nPadCount++;
422 ii++;
423 if (col-ii < 0) break;
424 }
425 ii = 0;
426 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
427 >= sigThresh) {
428 nPadCount++;
429 ii++;
430 if (col+ii+1 >= nColMax) break;
431 }
432
433 nClusters++;
434 switch (nPadCount) {
435 case 2:
436 iType = 0;
437 nClusters2pad++;
438 break;
439 case 3:
440 iType = 1;
441 nClusters3pad++;
442 break;
443 case 4:
444 iType = 2;
445 nClusters4pad++;
446 break;
447 case 5:
448 iType = 3;
449 nClusters5pad++;
450 break;
451 default:
452 iType = 4;
453 nClustersLarge++;
454 break;
455 };
456
457 // Don't analyze large clusters
458 //if (iType == 4) continue;
459
460 // Look for 5 pad cluster with minimum in the middle
461 Bool_t fivePadCluster = kFALSE;
3e1a3ad8 462 if (col < nColMax-3) {
463 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
db30bf0f 464 fivePadCluster = kTRUE;
465 }
466 if ((fivePadCluster) && (col < nColMax-5)) {
467 if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
468 fivePadCluster = kFALSE;
469 }
470 }
471 if ((fivePadCluster) && (col > 1)) {
472 if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
473 fivePadCluster = kFALSE;
474 }
475 }
476 }
477
478 // 5 pad cluster
479 // Modify the signal of the overlapping pad for the left part
480 // of the cluster which remains from a previous unfolding
481 if (iUnfold) {
482 clusterSignal[0] *= ratioLeft;
483 iType = 3;
484 iUnfold = 0;
485 }
486
487 // Unfold the 5 pad cluster
488 if (fivePadCluster) {
489 for (iPad = 0; iPad < kNsig; iPad++) {
490 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
491 ,col-1+iPad
492 ,time));
f7336fa3 493 }
db30bf0f 494 // Unfold the two maxima and set the signal on
495 // the overlapping pad to the ratio
496 ratioRight = Unfold(kEpsilon,padSignal);
497 ratioLeft = 1.0 - ratioRight;
498 clusterSignal[2] *= ratioRight;
499 iType = 3;
500 iUnfold = 1;
f7336fa3 501 }
f7336fa3 502
3e1a3ad8 503 Float_t clusterCharge = clusterSignal[0]
504 + clusterSignal[1]
505 + clusterSignal[2];
506
db30bf0f 507 // The position of the cluster
3e1a3ad8 508 clusterPads[0] = row + 0.5;
3e1a3ad8 509 // Take the shift of the additional time bins into account
510 clusterPads[2] = time - nTimeBefore + 0.5;
511
db30bf0f 512 if (fUseLUT) {
513
514 // Calculate the position of the cluster by using the
515 // lookup table method
516 Float_t ratioLUT;
517 Float_t signLUT;
518 Float_t lut = 0.0;
519 if (clusterSignal[0] > clusterSignal[2]) {
520 ratioLUT = clusterSignal[0] / clusterSignal[1];
521 signLUT = -1.0;
522 }
523 else {
524 ratioLUT = clusterSignal[2] / clusterSignal[1];
525 signLUT = 1.0;
526 }
527 if (ratioLUT < kLUTmin) {
528 lut = 0.0;
529 }
530 else if (ratioLUT > kLUTmax) {
531 lut = 0.5;
532 }
533 else {
534 Int_t indexLUT = TMath::Nint ((kNlut-1) * (ratioLUT - kLUTmin)
535 / (kLUTmax - kLUTmin));
536 lut = fLUT[indexLUT];
537 }
538 clusterPads[1] = col + 0.5 + signLUT * lut;
539
540 }
541 else {
542
543 // Calculate the position of the cluster by using the
544 // center of gravity method
545 clusterPads[1] = col + 0.5
546 + (clusterSignal[2] - clusterSignal[0])
547 / clusterCharge;
548
549 }
550
a819a5f7 551 Float_t clusterSigmaY2 = (clusterSignal[2] + clusterSignal[0]) / clusterCharge
552 - (clusterPads[1]-col-0.5) * (clusterPads[1]-col-0.5);
553
554 // Correct for ExB displacement
555 if (digitizer->GetExB()) {
556 Int_t local_time_bin = (Int_t) clusterPads[2];
557 Float_t driftLength = local_time_bin * timeBinSize + kAmWidth;
558 Float_t colSize = geo->GetColPadSize(iplan);
559 Float_t deltaY = omegaTau*driftLength/colSize;
560 clusterPads[1] = clusterPads[1] - deltaY;
561 }
562
47517f42 563 if (fVerbose > 1) {
3e1a3ad8 564 printf("-----------------------------------------------------------\n");
565 printf("Create cluster no. %d\n",nClusters);
566 printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
567 ,clusterPads[1]
568 ,clusterPads[2]);
569 printf("Indices: %d, %d, %d\n",clusterDigit[0]
570 ,clusterDigit[1]
571 ,clusterDigit[2]);
572 printf("Total charge = %f\n",clusterCharge);
573 printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
574 ,clusterTracks[1]
575 ,clusterTracks[2]);
576 printf(" pad1 %d, %d, %d\n",clusterTracks[3]
577 ,clusterTracks[4]
578 ,clusterTracks[5]);
579 printf(" pad2 %d, %d, %d\n",clusterTracks[6]
580 ,clusterTracks[7]
581 ,clusterTracks[8]);
db30bf0f 582 printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
f7336fa3 583 }
584
f7336fa3 585 // Add the cluster to the output array
3e1a3ad8 586 fTRD->AddCluster(clusterPads
587 ,clusterDigit
588 ,idet
589 ,clusterCharge
590 ,clusterTracks
a819a5f7 591 ,clusterSigmaY2
3e1a3ad8 592 ,iType);
f7336fa3 593
594 }
3e1a3ad8 595 }
596 }
597 }
f7336fa3 598
3e1a3ad8 599 // Compress the arrays
600 digits->Compress(1,0);
601 track0->Compress(1,0);
602 track1->Compress(1,0);
603 track2->Compress(1,0);
f7336fa3 604
3e1a3ad8 605 // Write the cluster and reset the array
793ff80c 606 WriteClusters(idet);
3e1a3ad8 607 fTRD->ResetRecPoints();
793ff80c 608
47517f42 609 if (fVerbose > 0) {
610 printf("AliTRDclusterizerV1::MakeCluster -- ");
611 printf("Found %d clusters in total.\n"
612 ,nClusters);
613 printf(" 2pad: %d\n",nClusters2pad);
614 printf(" 3pad: %d\n",nClusters3pad);
615 printf(" 4pad: %d\n",nClusters4pad);
616 printf(" 5pad: %d\n",nClusters5pad);
617 printf(" Large: %d\n",nClustersLarge);
618 }
f7336fa3 619
3e1a3ad8 620 }
621 }
622 }
f7336fa3 623
47517f42 624 if (fVerbose > 0) {
625 printf("AliTRDclusterizerV1::MakeCluster -- ");
626 printf("Done.\n");
627 }
f7336fa3 628
629 return kTRUE;
630
631}
632
633//_____________________________________________________________________________
634Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
635{
636 //
637 // Method to unfold neighbouring maxima.
638 // The charge ratio on the overlapping pad is calculated
639 // until there is no more change within the range given by eps.
640 // The resulting ratio is then returned to the calling method.
641 //
642
3e1a3ad8 643 Int_t itStep = 0; // Count iteration steps
f7336fa3 644
3e1a3ad8 645 Float_t ratio = 0.5; // Start value for ratio
646 Float_t prevRatio = 0; // Store previous ratio
f7336fa3 647
3e1a3ad8 648 Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
649 Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
f7336fa3 650
3e1a3ad8 651 // Start the iteration
f7336fa3 652 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
653
654 itStep++;
655 prevRatio = ratio;
656
3e1a3ad8 657 // Cluster position according to charge ratio
658 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
659 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
660 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
661 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3 662
3e1a3ad8 663 // Set cluster charge ratio
a819a5f7 664 Float_t ampLeft = padSignal[1] / PadResponse(0 - maxLeft );
665 Float_t ampRight = padSignal[3] / PadResponse(0 - maxRight);
f7336fa3 666
3e1a3ad8 667 // Apply pad response to parameters
668 newLeftSignal[0] = ampLeft * PadResponse(-1 - maxLeft);
669 newLeftSignal[1] = ampLeft * PadResponse( 0 - maxLeft);
670 newLeftSignal[2] = ampLeft * PadResponse( 1 - maxLeft);
f7336fa3 671
3e1a3ad8 672 newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
673 newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
674 newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
f7336fa3 675
3e1a3ad8 676 // Calculate new overlapping ratio
26edf6a4 677 ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] /
db30bf0f 678 (newLeftSignal[2] + newRightSignal[0]));
f7336fa3 679
680 }
681
682 return ratio;
683
684}
685
686//_____________________________________________________________________________
687Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
688{
689 //
690 // The pad response for the chevron pads.
691 // We use a simple Gaussian approximation which should be good
692 // enough for our purpose.
3e1a3ad8 693 // Updated for new PRF 1/5/01.
f7336fa3 694 //
695
696 // The parameters for the response function
3e1a3ad8 697 const Float_t kA = 0.8303;
698 const Float_t kB = -0.00392;
699 const Float_t kC = 0.472 * 0.472;
700 const Float_t kD = 2.19;
f7336fa3 701
3e1a3ad8 702 Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));
f7336fa3 703
704 return (pr);
705
706}