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 }