]> 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
Adaption to new fluka common blocks (E. Futo)
[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$
ad1fb852 18Revision 1.18 2002/10/14 14:57:43 hristov
19Merging the VirtualMC branch to the main development branch (HEAD)
20
b9d0a01d 21Revision 1.15.6.2 2002/07/24 10:09:30 alibrary
22Updating VirtualMC
23
24Revision 1.17 2002/06/12 09:54:35 cblume
25Update of tracking code provided by Sergei
26
5443e65e 27Revision 1.16 2002/03/25 20:01:30 cblume
28Introduce parameter class
29
17b26de4 30Revision 1.15 2001/11/14 12:09:11 cblume
31Use correct name for digitizer
32
47517f42 33Revision 1.14 2001/11/14 10:50:45 cblume
34Changes in digits IO. Add merging of summable digits
35
abaf1f1d 36Revision 1.13 2001/05/28 17:07:58 hristov
37Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
38
a819a5f7 39Revision 1.12 2001/05/21 17:42:58 hristov
40Constant casted to avoid the ambiguity
41
26edf6a4 42Revision 1.11 2001/05/21 16:45:47 hristov
43Last minute changes (C.Blume)
44
db30bf0f 45Revision 1.10 2001/05/07 08:06:44 cblume
46Speedup of the code. Create only AliTRDcluster
47
3e1a3ad8 48Revision 1.9 2000/11/01 14:53:20 cblume
49Merge with TRD-develop
50
793ff80c 51Revision 1.1.4.5 2000/10/15 23:40:01 cblume
52Remove AliTRDconst
53
54Revision 1.1.4.4 2000/10/06 16:49:46 cblume
55Made Getters const
56
57Revision 1.1.4.3 2000/10/04 16:34:58 cblume
58Replace include files by forward declarations
59
60Revision 1.1.4.2 2000/09/22 14:49:49 cblume
61Adapted to tracking code
62
63Revision 1.8 2000/10/02 21:28:19 fca
64Removal of useless dependecies via forward declarations
65
94de3818 66Revision 1.7 2000/06/27 13:08:50 cblume
67Changed to Copy(TObject &A) to appease the HP-compiler
68
43da34c0 69Revision 1.6 2000/06/09 11:10:07 cblume
70Compiler warnings and coding conventions, next round
71
dd9a6ee3 72Revision 1.5 2000/06/08 18:32:58 cblume
73Make code compliant to coding conventions
74
8230f242 75Revision 1.4 2000/06/07 16:27:01 cblume
76Try to remove compiler warnings on Sun and HP
77
9d0b222b 78Revision 1.3 2000/05/08 16:17:27 cblume
79Merge TRD-develop
80
6f1e466d 81Revision 1.1.4.1 2000/05/08 15:09:01 cblume
82Introduce AliTRDdigitsManager
83
c0dd96c3 84Revision 1.1 2000/02/28 18:58:54 cblume
85Add new TRD classes
86
f7336fa3 87*/
88
89///////////////////////////////////////////////////////////////////////////////
90// //
91// TRD cluster finder for the slow simulator.
92// //
93///////////////////////////////////////////////////////////////////////////////
94
95#include <TF1.h>
94de3818 96#include <TTree.h>
793ff80c 97#include <TH1.h>
a819a5f7 98#include <TFile.h>
f7336fa3 99
793ff80c 100#include "AliRun.h"
101
102#include "AliTRD.h"
f7336fa3 103#include "AliTRDclusterizerV1.h"
104#include "AliTRDmatrix.h"
105#include "AliTRDgeometry.h"
106#include "AliTRDdigitizer.h"
6f1e466d 107#include "AliTRDdataArrayF.h"
793ff80c 108#include "AliTRDdataArrayI.h"
109#include "AliTRDdigitsManager.h"
17b26de4 110#include "AliTRDparameter.h"
f7336fa3 111
112ClassImp(AliTRDclusterizerV1)
113
114//_____________________________________________________________________________
115AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
116{
117 //
118 // AliTRDclusterizerV1 default constructor
119 //
120
17b26de4 121 fDigitsManager = 0;
db30bf0f 122
f7336fa3 123}
124
125//_____________________________________________________________________________
126AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
127 :AliTRDclusterizer(name,title)
128{
129 //
130 // AliTRDclusterizerV1 default constructor
131 //
132
6f1e466d 133 fDigitsManager = new AliTRDdigitsManager();
17b26de4 134 fDigitsManager->CreateArrays();
f7336fa3 135
136}
137
8230f242 138//_____________________________________________________________________________
dd9a6ee3 139AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
8230f242 140{
141 //
142 // AliTRDclusterizerV1 copy constructor
143 //
144
dd9a6ee3 145 ((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242 146
147}
148
f7336fa3 149//_____________________________________________________________________________
150AliTRDclusterizerV1::~AliTRDclusterizerV1()
151{
8230f242 152 //
153 // AliTRDclusterizerV1 destructor
154 //
f7336fa3 155
6f1e466d 156 if (fDigitsManager) {
157 delete fDigitsManager;
abaf1f1d 158 fDigitsManager = NULL;
f7336fa3 159 }
160
161}
162
dd9a6ee3 163//_____________________________________________________________________________
164AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
165{
166 //
167 // Assignment operator
168 //
169
170 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
171 return *this;
172
173}
174
8230f242 175//_____________________________________________________________________________
43da34c0 176void AliTRDclusterizerV1::Copy(TObject &c)
8230f242 177{
178 //
179 // Copy function
180 //
181
17b26de4 182 ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
8230f242 183
184 AliTRDclusterizer::Copy(c);
185
186}
187
f7336fa3 188//_____________________________________________________________________________
189Bool_t AliTRDclusterizerV1::ReadDigits()
190{
191 //
192 // Reads the digits arrays from the input aliroot file
193 //
194
195 if (!fInputFile) {
17b26de4 196 printf("<AliTRDclusterizerV1::ReadDigits> ");
f7336fa3 197 printf("No input file open\n");
198 return kFALSE;
199 }
200
abaf1f1d 201 fDigitsManager->Open(fInputFile->GetName());
ad1fb852 202 fDigitsManager->SetEvent(fEvent);
abaf1f1d 203
f7336fa3 204 // Read in the digit arrays
6f1e466d 205 return (fDigitsManager->ReadDigits());
f7336fa3 206
207}
208
209//_____________________________________________________________________________
793ff80c 210Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3 211{
212 //
213 // Generates the cluster.
214 //
215
216 Int_t row, col, time;
217
3e1a3ad8 218 if (fTRD->IsVersion() != 1) {
17b26de4 219 printf("<AliTRDclusterizerV1::MakeCluster> ");
f7336fa3 220 printf("TRD must be version 1 (slow simulator).\n");
221 return kFALSE;
222 }
223
224 // Get the geometry
3e1a3ad8 225 AliTRDgeometry *geo = fTRD->GetGeometry();
f7336fa3 226
17b26de4 227 // Create a default parameter class if none is defined
228 if (!fPar) {
229 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
5443e65e 230 printf("<AliTRDclusterizerV1::MakeCluster> ");
231 printf("Create the default parameter object.\n");
17b26de4 232 }
233
234 Float_t timeBinSize = fPar->GetTimeBinSize();
a819a5f7 235 // Half of ampl.region
236 const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
237
17b26de4 238 Float_t omegaTau = fPar->GetOmegaTau();
47517f42 239 if (fVerbose > 0) {
17b26de4 240 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 241 printf("OmegaTau = %f \n",omegaTau);
17b26de4 242 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 243 printf("Start creating clusters.\n");
244 }
f7336fa3 245
8230f242 246 AliTRDdataArrayI *digits;
793ff80c 247 AliTRDdataArrayI *track0;
248 AliTRDdataArrayI *track1;
249 AliTRDdataArrayI *track2;
f7336fa3 250
3e1a3ad8 251 // Threshold value for the maximum
17b26de4 252 Int_t maxThresh = fPar->GetClusMaxThresh();
3e1a3ad8 253 // Threshold value for the digit signal
17b26de4 254 Int_t sigThresh = fPar->GetClusSigThresh();
f7336fa3 255
256 // Iteration limit for unfolding procedure
8230f242 257 const Float_t kEpsilon = 0.01;
f7336fa3 258
8230f242 259 const Int_t kNclus = 3;
260 const Int_t kNsig = 5;
3e1a3ad8 261 const Int_t kNtrack = 3 * kNclus;
262
db30bf0f 263 Int_t iType = 0;
264 Int_t iUnfold = 0;
265
266 Float_t ratioLeft = 1.0;
267 Float_t ratioRight = 1.0;
268
3e1a3ad8 269 Float_t padSignal[kNsig];
270 Float_t clusterSignal[kNclus];
271 Float_t clusterPads[kNclus];
272 Int_t clusterDigit[kNclus];
273 Int_t clusterTracks[kNtrack];
f7336fa3 274
275 Int_t chamBeg = 0;
793ff80c 276 Int_t chamEnd = AliTRDgeometry::Ncham();
3e1a3ad8 277 if (fTRD->GetSensChamber() >= 0) {
278 chamBeg = fTRD->GetSensChamber();
6f1e466d 279 chamEnd = chamBeg + 1;
f7336fa3 280 }
281 Int_t planBeg = 0;
793ff80c 282 Int_t planEnd = AliTRDgeometry::Nplan();
3e1a3ad8 283 if (fTRD->GetSensPlane() >= 0) {
284 planBeg = fTRD->GetSensPlane();
f7336fa3 285 planEnd = planBeg + 1;
286 }
287 Int_t sectBeg = 0;
793ff80c 288 Int_t sectEnd = AliTRDgeometry::Nsect();
f7336fa3 289
3e1a3ad8 290 // Start clustering in every chamber
f7336fa3 291 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
292 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
293 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
294
3e1a3ad8 295 if (fTRD->GetSensSector() >= 0) {
296 Int_t sens1 = fTRD->GetSensSector();
297 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
793ff80c 298 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
299 * AliTRDgeometry::Nsect();
dd9a6ee3 300 if (sens1 < sens2) {
9d0b222b 301 if ((isect < sens1) || (isect >= sens2)) continue;
dd9a6ee3 302 }
303 else {
9d0b222b 304 if ((isect < sens1) && (isect >= sens2)) continue;
dd9a6ee3 305 }
9d0b222b 306 }
307
8230f242 308 Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3 309
db30bf0f 310 Int_t nClusters = 0;
311 Int_t nClusters2pad = 0;
312 Int_t nClusters3pad = 0;
313 Int_t nClusters4pad = 0;
314 Int_t nClusters5pad = 0;
315 Int_t nClustersLarge = 0;
3e1a3ad8 316
47517f42 317 if (fVerbose > 0) {
17b26de4 318 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 319 printf("Analyzing chamber %d, plane %d, sector %d.\n"
320 ,icham,iplan,isect);
321 }
f7336fa3 322
5443e65e 323 Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect);
324 Int_t nColMax = fPar->GetColMax(iplan);
325 Int_t nTimeBefore = fPar->GetTimeBefore();
326 Int_t nTimeTotal = fPar->GetTimeTotal();
327
328 Float_t row0 = fPar->GetRow0(iplan,icham,isect);
329 Float_t col0 = fPar->GetCol0(iplan);
330 Float_t rowSize = fPar->GetRowPadSize(iplan,icham,isect);
331 Float_t colSize = fPar->GetColPadSize(iplan);
f7336fa3 332
3e1a3ad8 333 // Get the digits
8230f242 334 digits = fDigitsManager->GetDigits(idet);
3e1a3ad8 335 digits->Expand();
793ff80c 336 track0 = fDigitsManager->GetDictionary(idet,0);
3e1a3ad8 337 track0->Expand();
793ff80c 338 track1 = fDigitsManager->GetDictionary(idet,1);
3e1a3ad8 339 track1->Expand();
793ff80c 340 track2 = fDigitsManager->GetDictionary(idet,2);
3e1a3ad8 341 track2->Expand();
342
343 // Loop through the chamber and find the maxima
344 for ( row = 0; row < nRowMax; row++) {
345 for ( col = 2; col < nColMax; col++) {
346 for (time = 0; time < nTimeTotal; time++) {
347
a819a5f7 348 Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
349 Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
350 Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
3e1a3ad8 351
352 // Look for the maximum
db30bf0f 353 if (signalM >= maxThresh) {
354 if (((signalL >= sigThresh) &&
355 (signalL < signalM)) ||
356 ((signalR >= sigThresh) &&
357 (signalR < signalM))) {
3e1a3ad8 358 // Maximum found, mark the position by a negative signal
359 digits->SetDataUnchecked(row,col-1,time,-signalM);
360 }
361 }
362
363 }
364 }
365 }
366
367 // Now check the maxima and calculate the cluster position
368 for ( row = 0; row < nRowMax ; row++) {
db30bf0f 369 for (time = 0; time < nTimeTotal; time++) {
370 for ( col = 1; col < nColMax-1; col++) {
3e1a3ad8 371
372 // Maximum found ?
373 if (digits->GetDataUnchecked(row,col,time) < 0) {
f7336fa3 374
9d0b222b 375 Int_t iPad;
8230f242 376 for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8 377 Int_t iPadCol = col - 1 + iPad;
378 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
379 ,iPadCol
380 ,time));
381 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
382 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
383 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
384 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
f7336fa3 385 }
386
db30bf0f 387 // Count the number of pads in the cluster
388 Int_t nPadCount = 0;
389 Int_t ii = 0;
390 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
391 >= sigThresh) {
392 nPadCount++;
393 ii++;
394 if (col-ii < 0) break;
395 }
396 ii = 0;
397 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
398 >= sigThresh) {
399 nPadCount++;
400 ii++;
401 if (col+ii+1 >= nColMax) break;
402 }
403
404 nClusters++;
405 switch (nPadCount) {
406 case 2:
407 iType = 0;
408 nClusters2pad++;
409 break;
410 case 3:
411 iType = 1;
412 nClusters3pad++;
413 break;
414 case 4:
415 iType = 2;
416 nClusters4pad++;
417 break;
418 case 5:
419 iType = 3;
420 nClusters5pad++;
421 break;
422 default:
423 iType = 4;
424 nClustersLarge++;
425 break;
426 };
427
428 // Don't analyze large clusters
429 //if (iType == 4) continue;
430
431 // Look for 5 pad cluster with minimum in the middle
432 Bool_t fivePadCluster = kFALSE;
3e1a3ad8 433 if (col < nColMax-3) {
434 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
db30bf0f 435 fivePadCluster = kTRUE;
436 }
437 if ((fivePadCluster) && (col < nColMax-5)) {
438 if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
439 fivePadCluster = kFALSE;
440 }
441 }
442 if ((fivePadCluster) && (col > 1)) {
443 if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
444 fivePadCluster = kFALSE;
445 }
446 }
447 }
448
449 // 5 pad cluster
450 // Modify the signal of the overlapping pad for the left part
451 // of the cluster which remains from a previous unfolding
452 if (iUnfold) {
453 clusterSignal[0] *= ratioLeft;
454 iType = 3;
455 iUnfold = 0;
456 }
457
458 // Unfold the 5 pad cluster
459 if (fivePadCluster) {
460 for (iPad = 0; iPad < kNsig; iPad++) {
461 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
462 ,col-1+iPad
463 ,time));
f7336fa3 464 }
db30bf0f 465 // Unfold the two maxima and set the signal on
466 // the overlapping pad to the ratio
17b26de4 467 ratioRight = Unfold(kEpsilon,iplan,padSignal);
db30bf0f 468 ratioLeft = 1.0 - ratioRight;
469 clusterSignal[2] *= ratioRight;
470 iType = 3;
471 iUnfold = 1;
f7336fa3 472 }
f7336fa3 473
3e1a3ad8 474 Float_t clusterCharge = clusterSignal[0]
475 + clusterSignal[1]
476 + clusterSignal[2];
477
db30bf0f 478 // The position of the cluster
3e1a3ad8 479 clusterPads[0] = row + 0.5;
3e1a3ad8 480 // Take the shift of the additional time bins into account
481 clusterPads[2] = time - nTimeBefore + 0.5;
482
17b26de4 483 if (fPar->LUTOn()) {
db30bf0f 484
485 // Calculate the position of the cluster by using the
486 // lookup table method
5443e65e 487 clusterPads[1] = col + 0.5
488 + fPar->LUTposition(iplan,clusterSignal[0]
17b26de4 489 ,clusterSignal[1]
490 ,clusterSignal[2]);
db30bf0f 491
492 }
493 else {
494
495 // Calculate the position of the cluster by using the
496 // center of gravity method
497 clusterPads[1] = col + 0.5
498 + (clusterSignal[2] - clusterSignal[0])
499 / clusterCharge;
500
501 }
502
5443e65e 503 Float_t q0 = clusterSignal[0];
504 Float_t q1 = clusterSignal[1];
505 Float_t q2 = clusterSignal[2];
506 Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
507 (clusterCharge*clusterCharge);
a819a5f7 508
509 // Correct for ExB displacement
17b26de4 510 if (fPar->ExBOn()) {
a819a5f7 511 Int_t local_time_bin = (Int_t) clusterPads[2];
512 Float_t driftLength = local_time_bin * timeBinSize + kAmWidth;
17b26de4 513 Float_t colSize = fPar->GetColPadSize(iplan);
a819a5f7 514 Float_t deltaY = omegaTau*driftLength/colSize;
515 clusterPads[1] = clusterPads[1] - deltaY;
516 }
517
47517f42 518 if (fVerbose > 1) {
3e1a3ad8 519 printf("-----------------------------------------------------------\n");
520 printf("Create cluster no. %d\n",nClusters);
521 printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
522 ,clusterPads[1]
523 ,clusterPads[2]);
524 printf("Indices: %d, %d, %d\n",clusterDigit[0]
525 ,clusterDigit[1]
526 ,clusterDigit[2]);
527 printf("Total charge = %f\n",clusterCharge);
528 printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
529 ,clusterTracks[1]
530 ,clusterTracks[2]);
531 printf(" pad1 %d, %d, %d\n",clusterTracks[3]
532 ,clusterTracks[4]
533 ,clusterTracks[5]);
534 printf(" pad2 %d, %d, %d\n",clusterTracks[6]
535 ,clusterTracks[7]
536 ,clusterTracks[8]);
db30bf0f 537 printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
f7336fa3 538 }
539
5443e65e 540 // Calculate the position and the error
541 Float_t clusterPos[3];
542 clusterPos[0] = clusterPads[1] * colSize + col0;
543 clusterPos[1] = clusterPads[0] * rowSize + row0;
544 clusterPos[2] = clusterPads[2];
545 Float_t clusterSig[2];
546 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
547 clusterSig[1] = rowSize * rowSize / 12.;
548
f7336fa3 549 // Add the cluster to the output array
5443e65e 550 fTRD->AddCluster(clusterPos
3e1a3ad8 551 ,idet
552 ,clusterCharge
553 ,clusterTracks
5443e65e 554 ,clusterSig
3e1a3ad8 555 ,iType);
f7336fa3 556
557 }
3e1a3ad8 558 }
559 }
560 }
f7336fa3 561
3e1a3ad8 562 // Compress the arrays
563 digits->Compress(1,0);
564 track0->Compress(1,0);
565 track1->Compress(1,0);
566 track2->Compress(1,0);
f7336fa3 567
3e1a3ad8 568 // Write the cluster and reset the array
793ff80c 569 WriteClusters(idet);
3e1a3ad8 570 fTRD->ResetRecPoints();
793ff80c 571
47517f42 572 if (fVerbose > 0) {
17b26de4 573 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 574 printf("Found %d clusters in total.\n"
575 ,nClusters);
576 printf(" 2pad: %d\n",nClusters2pad);
577 printf(" 3pad: %d\n",nClusters3pad);
578 printf(" 4pad: %d\n",nClusters4pad);
579 printf(" 5pad: %d\n",nClusters5pad);
580 printf(" Large: %d\n",nClustersLarge);
581 }
f7336fa3 582
3e1a3ad8 583 }
584 }
585 }
f7336fa3 586
47517f42 587 if (fVerbose > 0) {
17b26de4 588 printf("<AliTRDclusterizerV1::MakeCluster> ");
47517f42 589 printf("Done.\n");
590 }
f7336fa3 591
592 return kTRUE;
593
594}
595
596//_____________________________________________________________________________
17b26de4 597Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal)
f7336fa3 598{
599 //
600 // Method to unfold neighbouring maxima.
601 // The charge ratio on the overlapping pad is calculated
602 // until there is no more change within the range given by eps.
603 // The resulting ratio is then returned to the calling method.
604 //
605
17b26de4 606 Int_t irc = 0;
3e1a3ad8 607 Int_t itStep = 0; // Count iteration steps
f7336fa3 608
3e1a3ad8 609 Float_t ratio = 0.5; // Start value for ratio
610 Float_t prevRatio = 0; // Store previous ratio
f7336fa3 611
3e1a3ad8 612 Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
613 Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
17b26de4 614 Float_t newSignal[3] = {0};
f7336fa3 615
3e1a3ad8 616 // Start the iteration
f7336fa3 617 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
618
619 itStep++;
620 prevRatio = ratio;
621
3e1a3ad8 622 // Cluster position according to charge ratio
623 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
624 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
625 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
626 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3 627
3e1a3ad8 628 // Set cluster charge ratio
17b26de4 629 irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
630 Float_t ampLeft = padSignal[1] / newSignal[1];
631 irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
632 Float_t ampRight = padSignal[3] / newSignal[1];
f7336fa3 633
3e1a3ad8 634 // Apply pad response to parameters
17b26de4 635 irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
636 irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
f7336fa3 637
3e1a3ad8 638 // Calculate new overlapping ratio
26edf6a4 639 ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] /
db30bf0f 640 (newLeftSignal[2] + newRightSignal[0]));
f7336fa3 641
642 }
643
644 return ratio;
645
646}
647