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