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