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