Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of...
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
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$
18 Revision 1.12  2001/05/21 17:42:58  hristov
19 Constant casted to avoid the ambiguity
20
21 Revision 1.11  2001/05/21 16:45:47  hristov
22 Last minute changes (C.Blume)
23
24 Revision 1.10  2001/05/07 08:06:44  cblume
25 Speedup of the code. Create only AliTRDcluster
26
27 Revision 1.9  2000/11/01 14:53:20  cblume
28 Merge with TRD-develop
29
30 Revision 1.1.4.5  2000/10/15 23:40:01  cblume
31 Remove AliTRDconst
32
33 Revision 1.1.4.4  2000/10/06 16:49:46  cblume
34 Made Getters const
35
36 Revision 1.1.4.3  2000/10/04 16:34:58  cblume
37 Replace include files by forward declarations
38
39 Revision 1.1.4.2  2000/09/22 14:49:49  cblume
40 Adapted to tracking code
41
42 Revision 1.8  2000/10/02 21:28:19  fca
43 Removal of useless dependecies via forward declarations
44
45 Revision 1.7  2000/06/27 13:08:50  cblume
46 Changed to Copy(TObject &A) to appease the HP-compiler
47
48 Revision 1.6  2000/06/09 11:10:07  cblume
49 Compiler warnings and coding conventions, next round
50
51 Revision 1.5  2000/06/08 18:32:58  cblume
52 Make code compliant to coding conventions
53
54 Revision 1.4  2000/06/07 16:27:01  cblume
55 Try to remove compiler warnings on Sun and HP
56
57 Revision 1.3  2000/05/08 16:17:27  cblume
58 Merge TRD-develop
59
60 Revision 1.1.4.1  2000/05/08 15:09:01  cblume
61 Introduce AliTRDdigitsManager
62
63 Revision 1.1  2000/02/28 18:58:54  cblume
64 Add new TRD classes
65
66 */
67
68 ///////////////////////////////////////////////////////////////////////////////
69 //                                                                           //
70 // TRD cluster finder for the slow simulator. 
71 //                                                                           //
72 ///////////////////////////////////////////////////////////////////////////////
73
74 #include <TF1.h>
75 #include <TTree.h>
76 #include <TH1.h>
77 #include <TFile.h>
78
79 #include "AliRun.h"
80
81 #include "AliTRD.h"
82 #include "AliTRDclusterizerV1.h"
83 #include "AliTRDmatrix.h"
84 #include "AliTRDgeometry.h"
85 #include "AliTRDdigitizer.h"
86 #include "AliTRDdataArrayF.h"
87 #include "AliTRDdataArrayI.h"
88 #include "AliTRDdigitsManager.h"
89
90 ClassImp(AliTRDclusterizerV1)
91
92 //_____________________________________________________________________________
93 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
94 {
95   //
96   // AliTRDclusterizerV1 default constructor
97   //
98
99   fDigitsManager = NULL;
100
101   fClusMaxThresh = 0;
102   fClusSigThresh = 0;
103
104   fUseLUT        = kFALSE;
105
106 }
107
108 //_____________________________________________________________________________
109 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
110                     :AliTRDclusterizer(name,title)
111 {
112   //
113   // AliTRDclusterizerV1 default constructor
114   //
115
116   fDigitsManager = new AliTRDdigitsManager();
117
118   Init();
119
120 }
121
122 //_____________________________________________________________________________
123 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
124 {
125   //
126   // AliTRDclusterizerV1 copy constructor
127   //
128
129   ((AliTRDclusterizerV1 &) c).Copy(*this);
130
131 }
132
133 //_____________________________________________________________________________
134 AliTRDclusterizerV1::~AliTRDclusterizerV1()
135 {
136   //
137   // AliTRDclusterizerV1 destructor
138   //
139
140   if (fDigitsManager) {
141     delete fDigitsManager;
142   }
143
144 }
145
146 //_____________________________________________________________________________
147 AliTRDclusterizerV1 &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
158 //_____________________________________________________________________________
159 void AliTRDclusterizerV1::Copy(TObject &c)
160 {
161   //
162   // Copy function
163   //
164
165   ((AliTRDclusterizerV1 &) c).fUseLUT        = fUseLUT;
166   ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
167   ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
168   ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
169   for (Int_t ilut = 0; ilut < kNlut; ilut++) {
170     ((AliTRDclusterizerV1 &) c).fLUT[ilut] = fLUT[ilut];
171   }
172
173   AliTRDclusterizer::Copy(c);
174
175 }
176
177 //_____________________________________________________________________________
178 void AliTRDclusterizerV1::Init()
179 {
180   //
181   // Initializes the cluster finder
182   //
183
184   // The default parameter for the clustering
185   fClusMaxThresh = 3;
186   fClusSigThresh = 1;
187
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
214 }
215
216 //_____________________________________________________________________________
217 Bool_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
229   // Read in the digit arrays
230   return (fDigitsManager->ReadDigits());  
231
232 }
233
234 //_____________________________________________________________________________
235 Bool_t AliTRDclusterizerV1::MakeClusters()
236 {
237   //
238   // Generates the cluster.
239   //
240
241   Int_t row, col, time;
242
243   if (fTRD->IsVersion() != 1) {
244     printf("AliTRDclusterizerV1::MakeCluster -- ");
245     printf("TRD must be version 1 (slow simulator).\n");
246     return kFALSE; 
247   }
248
249   // Get the geometry
250   AliTRDgeometry *geo = fTRD->GetGeometry();
251
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  
263   printf("AliTRDclusterizerV1::MakeCluster -- ");
264   printf("Start creating clusters.\n");
265
266   AliTRDdataArrayI *digits;
267   AliTRDdataArrayI *track0;
268   AliTRDdataArrayI *track1;
269   AliTRDdataArrayI *track2; 
270
271   // Threshold value for the maximum
272   Int_t maxThresh = fClusMaxThresh;   
273   // Threshold value for the digit signal
274   Int_t sigThresh = fClusSigThresh;   
275
276   // Iteration limit for unfolding procedure
277   const Float_t kEpsilon = 0.01;             
278
279   const Int_t   kNclus   = 3;  
280   const Int_t   kNsig    = 5;
281   const Int_t   kNtrack  = 3 * kNclus;
282
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
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];   
298
299   Int_t chamBeg = 0;
300   Int_t chamEnd = AliTRDgeometry::Ncham();
301   if (fTRD->GetSensChamber()  >= 0) {
302     chamBeg = fTRD->GetSensChamber();
303     chamEnd = chamBeg + 1;
304   }
305   Int_t planBeg = 0;
306   Int_t planEnd = AliTRDgeometry::Nplan();
307   if (fTRD->GetSensPlane()    >= 0) {
308     planBeg = fTRD->GetSensPlane();
309     planEnd = planBeg + 1;
310   }
311   Int_t sectBeg = 0;
312   Int_t sectEnd = AliTRDgeometry::Nsect();
313
314   // Start clustering in every chamber
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
319         if (fTRD->GetSensSector() >= 0) {
320           Int_t sens1 = fTRD->GetSensSector();
321           Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
322           sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
323                  * AliTRDgeometry::Nsect();
324           if (sens1 < sens2) {
325             if ((isect < sens1) || (isect >= sens2)) continue;
326           }
327           else {
328             if ((isect < sens1) && (isect >= sens2)) continue;
329           }
330         }
331
332         Int_t idet = geo->GetDetector(iplan,icham,isect);
333
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;
340
341         printf("AliTRDclusterizerV1::MakeCluster -- ");
342         printf("Analyzing chamber %d, plane %d, sector %d.\n"
343               ,icham,iplan,isect);
344
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();  
349
350         // Get the digits
351         digits = fDigitsManager->GetDigits(idet);
352         digits->Expand();
353         track0 = fDigitsManager->GetDictionary(idet,0);
354         track0->Expand();
355         track1 = fDigitsManager->GetDictionary(idet,1);
356         track1->Expand();
357         track2 = fDigitsManager->GetDictionary(idet,2); 
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
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));
368  
369               // Look for the maximum
370               if (signalM >= maxThresh) {
371                 if (((signalL >= sigThresh) &&
372                      (signalL <  signalM))  ||
373                     ((signalR >= sigThresh) &&
374                      (signalR <  signalM))) {
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++) {
386           for (time = 0; time < nTimeTotal; time++) {
387             for ( col = 1;  col <  nColMax-1;  col++) {
388
389               // Maximum found ?             
390               if (digits->GetDataUnchecked(row,col,time) < 0) {
391
392                 Int_t iPad;
393                 for (iPad = 0; iPad < kNclus; iPad++) {
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;
402                 }
403
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;
450                 if (col < nColMax-3) {
451                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
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));
481                   }
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;
489                 }
490
491                 Float_t clusterCharge = clusterSignal[0]
492                                       + clusterSignal[1]
493                                       + clusterSignal[2];
494                 
495                 // The position of the cluster
496                 clusterPads[0] = row + 0.5;
497                 // Take the shift of the additional time bins into account
498                 clusterPads[2] = time - nTimeBefore + 0.5;
499
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
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                                        
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]);
570                   printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
571                 }
572
573                 // Add the cluster to the output array 
574                 fTRD->AddCluster(clusterPads
575                                 ,clusterDigit
576                                 ,idet
577                                 ,clusterCharge
578                                 ,clusterTracks
579                                 ,clusterSigmaY2
580                                 ,iType);
581
582               }
583             } 
584           }   
585         }     
586
587         // Compress the arrays
588         digits->Compress(1,0);
589         track0->Compress(1,0);
590         track1->Compress(1,0);
591         track2->Compress(1,0);
592
593         // Write the cluster and reset the array
594         WriteClusters(idet);
595         fTRD->ResetRecPoints();
596
597         printf("AliTRDclusterizerV1::MakeCluster -- ");
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);
605
606       }    
607     }      
608   }        
609
610   printf("AliTRDclusterizerV1::MakeCluster -- ");
611   printf("Done.\n");
612
613   return kTRUE;
614
615 }
616
617 //_____________________________________________________________________________
618 Float_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
627   Int_t   itStep            = 0;      // Count iteration steps
628
629   Float_t ratio             = 0.5;    // Start value for ratio
630   Float_t prevRatio         = 0;      // Store previous ratio
631
632   Float_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
633   Float_t newRightSignal[3] = {0};    // Array to store right cluster signal
634
635   // Start the iteration
636   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
637
638     itStep++;
639     prevRatio = ratio;
640
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]);
646
647     // Set cluster charge ratio
648     Float_t ampLeft  = padSignal[1] / PadResponse(0 - maxLeft );
649     Float_t ampRight = padSignal[3] / PadResponse(0 - maxRight);
650
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);
655
656     newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
657     newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
658     newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
659
660     // Calculate new overlapping ratio
661     ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] / 
662                           (newLeftSignal[2] + newRightSignal[0]));
663
664   }
665
666   return ratio;
667
668 }
669
670 //_____________________________________________________________________________
671 Float_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.
677   // Updated for new PRF 1/5/01.
678   //
679
680   // The parameters for the response function
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;
685
686   Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));
687
688   return (pr);
689
690 }