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