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