]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDclusterizerV1.cxx
Reconstruction of RAW data. Introduction of cluster finder (A. de Caro)
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
index 6c793795353f3becbd2457512bc77c0b334f55cc..c5fff1f5e931a23ddb8536e8a19077d93ee39ac7 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.12  2001/05/21 17:42:58  hristov
-Constant casted to avoid the ambiguity
-
-Revision 1.11  2001/05/21 16:45:47  hristov
-Last minute changes (C.Blume)
-
-Revision 1.10  2001/05/07 08:06:44  cblume
-Speedup of the code. Create only AliTRDcluster
-
-Revision 1.9  2000/11/01 14:53:20  cblume
-Merge with TRD-develop
-
-Revision 1.1.4.5  2000/10/15 23:40:01  cblume
-Remove AliTRDconst
-
-Revision 1.1.4.4  2000/10/06 16:49:46  cblume
-Made Getters const
-
-Revision 1.1.4.3  2000/10/04 16:34:58  cblume
-Replace include files by forward declarations
-
-Revision 1.1.4.2  2000/09/22 14:49:49  cblume
-Adapted to tracking code
-
-Revision 1.8  2000/10/02 21:28:19  fca
-Removal of useless dependecies via forward declarations
-
-Revision 1.7  2000/06/27 13:08:50  cblume
-Changed to Copy(TObject &A) to appease the HP-compiler
-
-Revision 1.6  2000/06/09 11:10:07  cblume
-Compiler warnings and coding conventions, next round
-
-Revision 1.5  2000/06/08 18:32:58  cblume
-Make code compliant to coding conventions
-
-Revision 1.4  2000/06/07 16:27:01  cblume
-Try to remove compiler warnings on Sun and HP
-
-Revision 1.3  2000/05/08 16:17:27  cblume
-Merge TRD-develop
-
-Revision 1.1.4.1  2000/05/08 15:09:01  cblume
-Introduce AliTRDdigitsManager
-
-Revision 1.1  2000/02/28 18:58:54  cblume
-Add new TRD classes
-
-*/
+/* $Id$ */
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
@@ -77,15 +27,17 @@ Add new TRD classes
 #include <TFile.h>
 
 #include "AliRun.h"
+#include "AliRunLoader.h"
+#include "AliLoader.h"
 
-#include "AliTRD.h"
 #include "AliTRDclusterizerV1.h"
 #include "AliTRDmatrix.h"
 #include "AliTRDgeometry.h"
-#include "AliTRDdigitizer.h"
 #include "AliTRDdataArrayF.h"
 #include "AliTRDdataArrayI.h"
 #include "AliTRDdigitsManager.h"
+#include "AliTRDparameter.h"
+#include "AliTRDpadPlane.h"
 
 ClassImp(AliTRDclusterizerV1)
 
@@ -96,12 +48,7 @@ AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
   // AliTRDclusterizerV1 default constructor
   //
 
-  fDigitsManager = NULL;
-
-  fClusMaxThresh = 0;
-  fClusSigThresh = 0;
-
-  fUseLUT        = kFALSE;
+  fDigitsManager = 0;
 
 }
 
@@ -114,13 +61,13 @@ AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title
   //
 
   fDigitsManager = new AliTRDdigitsManager();
-
-  Init();
+  fDigitsManager->CreateArrays();
 
 }
 
 //_____________________________________________________________________________
 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
+:AliTRDclusterizer(c)
 {
   //
   // AliTRDclusterizerV1 copy constructor
@@ -139,6 +86,7 @@ AliTRDclusterizerV1::~AliTRDclusterizerV1()
 
   if (fDigitsManager) {
     delete fDigitsManager;
+    fDigitsManager = NULL;
   }
 
 }
@@ -156,63 +104,18 @@ AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c
 }
 
 //_____________________________________________________________________________
-void AliTRDclusterizerV1::Copy(TObject &c)
+void AliTRDclusterizerV1::Copy(TObject &c) const
 {
   //
   // Copy function
   //
 
-  ((AliTRDclusterizerV1 &) c).fUseLUT        = fUseLUT;
-  ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
-  ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
-  ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
-  for (Int_t ilut = 0; ilut < kNlut; ilut++) {
-    ((AliTRDclusterizerV1 &) c).fLUT[ilut] = fLUT[ilut];
-  }
+  ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
 
   AliTRDclusterizer::Copy(c);
 
 }
 
-//_____________________________________________________________________________
-void AliTRDclusterizerV1::Init()
-{
-  //
-  // Initializes the cluster finder
-  //
-
-  // The default parameter for the clustering
-  fClusMaxThresh = 3;
-  fClusSigThresh = 1;
-
-  // Use the lookup table for the position determination
-  fUseLUT        = kTRUE;
-
-  // The lookup table from Bogdan
-  Float_t lut[128] = {  
-    0.0068,  0.0198,  0.0318,  0.0432,  0.0538,  0.0642,  0.0742,  0.0838,
-    0.0932,  0.1023,  0.1107,  0.1187,  0.1268,  0.1347,  0.1423,  0.1493,  
-    0.1562,  0.1632,  0.1698,  0.1762,  0.1828,  0.1887,  0.1947,  0.2002,  
-    0.2062,  0.2118,  0.2173,  0.2222,  0.2278,  0.2327,  0.2377,  0.2428,  
-    0.2473,  0.2522,  0.2567,  0.2612,  0.2657,  0.2697,  0.2743,  0.2783,  
-    0.2822,  0.2862,  0.2903,  0.2943,  0.2982,  0.3018,  0.3058,  0.3092,  
-    0.3128,  0.3167,  0.3203,  0.3237,  0.3268,  0.3302,  0.3338,  0.3368,  
-    0.3402,  0.3433,  0.3462,  0.3492,  0.3528,  0.3557,  0.3587,  0.3613,  
-    0.3643,  0.3672,  0.3702,  0.3728,  0.3758,  0.3783,  0.3812,  0.3837,  
-    0.3862,  0.3887,  0.3918,  0.3943,  0.3968,  0.3993,  0.4017,  0.4042,  
-    0.4067,  0.4087,  0.4112,  0.4137,  0.4157,  0.4182,  0.4207,  0.4227,  
-    0.4252,  0.4272,  0.4293,  0.4317,  0.4338,  0.4358,  0.4383,  0.4403,  
-    0.4423,  0.4442,  0.4462,  0.4482,  0.4502,  0.4523,  0.4543,  0.4563,  
-    0.4582,  0.4602,  0.4622,  0.4638,  0.4658,  0.4678,  0.4697,  0.4712,  
-    0.4733,  0.4753,  0.4767,  0.4787,  0.4803,  0.4823,  0.4837,  0.4857,  
-    0.4873,  0.4888,  0.4908,  0.4922,  0.4942,  0.4958,  0.4972,  0.4988  
-  }; 
-  for (Int_t ilut = 0; ilut < kNlut; ilut++) {
-    fLUT[ilut] = lut[ilut];
-  }
-
-}
-
 //_____________________________________________________________________________
 Bool_t AliTRDclusterizerV1::ReadDigits()
 {
@@ -220,14 +123,16 @@ Bool_t AliTRDclusterizerV1::ReadDigits()
   // Reads the digits arrays from the input aliroot file
   //
 
-  if (!fInputFile) {
-    printf("AliTRDclusterizerV1::ReadDigits -- ");
+  if (!fRunLoader) {
+    printf("<AliTRDclusterizerV1::ReadDigits> ");
     printf("No input file open\n");
     return kFALSE;
   }
+  AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
+  if (!loader->TreeD()) loader->LoadDigits();
 
   // Read in the digit arrays
-  return (fDigitsManager->ReadDigits());  
+  return (fDigitsManager->ReadDigits(loader->TreeD()));
 
 }
 
@@ -240,28 +145,37 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
 
   Int_t row, col, time;
 
+  /*
   if (fTRD->IsVersion() != 1) {
-    printf("AliTRDclusterizerV1::MakeCluster -- ");
+    printf("<AliTRDclusterizerV1::MakeCluster> ");
     printf("TRD must be version 1 (slow simulator).\n");
     return kFALSE; 
   }
+  */
 
   // Get the geometry
-  AliTRDgeometry *geo = fTRD->GetGeometry();
+  AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
 
-  Float_t timeBinSize = geo->GetTimeBinSize();
+  // Create a default parameter class if none is defined
+  if (!fPar) {
+    fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
+    printf("<AliTRDclusterizerV1::MakeCluster> ");
+    printf("Create the default parameter object.\n");
+  }
+  fPar->Init();
+
+  //Float_t timeBinSize = fPar->GetDriftVelocity()
+  //                    / fPar->GetSamplingFrequency();
   // Half of ampl.region
-  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
-
-  AliTRDdigitizer *digitizer = (AliTRDdigitizer*) fInputFile->Get("digitizer");
-  printf("AliTRDclusterizerV1::MakeCluster -- ");
-  printf("Got digitizer\n");
-  Float_t omegaTau = digitizer->GetOmegaTau();
-  printf("AliTRDclusterizerV1::MakeCluster -- ");
-  printf("OmegaTau = %f \n",omegaTau);
-  printf("AliTRDclusterizerV1::MakeCluster -- ");
-  printf("Start creating clusters.\n");
+  //  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
+
+  Float_t omegaTau = fPar->GetOmegaTau();
+  if (fVerbose > 0) {
+    printf("<AliTRDclusterizerV1::MakeCluster> ");
+    printf("OmegaTau = %f \n",omegaTau);
+    printf("<AliTRDclusterizerV1::MakeCluster> ");
+    printf("Start creating clusters.\n");
+  } 
 
   AliTRDdataArrayI *digits;
   AliTRDdataArrayI *track0;
@@ -269,10 +183,9 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   AliTRDdataArrayI *track2; 
 
   // Threshold value for the maximum
-  Int_t maxThresh = fClusMaxThresh;   
+  Int_t maxThresh = fPar->GetClusMaxThresh();   
   // Threshold value for the digit signal
-  Int_t sigThresh = fClusSigThresh;   
-
+  Int_t sigThresh = fPar->GetClusSigThresh();   
   // Iteration limit for unfolding procedure
   const Float_t kEpsilon = 0.01;             
 
@@ -280,55 +193,30 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   const Int_t   kNsig    = 5;
   const Int_t   kNtrack  = 3 * kNclus;
 
-  // For the LUT
-  const Float_t kLUTmin  = 0.106113;
-  const Float_t kLUTmax  = 0.995415;
-
-  Int_t   iType          = 0;
-  Int_t   iUnfold        = 0;
+  Int_t    iType         = 0;
+  Int_t    iUnfold       = 0;  
+  Double_t ratioLeft     = 1.0;
+  Double_t ratioRight    = 1.0;
 
-  Float_t ratioLeft      = 1.0;
-  Float_t ratioRight     = 1.0;
-
-  Float_t padSignal[kNsig];   
-  Float_t clusterSignal[kNclus];
-  Float_t clusterPads[kNclus];   
-  Int_t   clusterDigit[kNclus];
-  Int_t   clusterTracks[kNtrack];   
-
-  Int_t chamBeg = 0;
-  Int_t chamEnd = AliTRDgeometry::Ncham();
-  if (fTRD->GetSensChamber()  >= 0) {
-    chamBeg = fTRD->GetSensChamber();
-    chamEnd = chamBeg + 1;
-  }
-  Int_t planBeg = 0;
-  Int_t planEnd = AliTRDgeometry::Nplan();
-  if (fTRD->GetSensPlane()    >= 0) {
-    planBeg = fTRD->GetSensPlane();
-    planEnd = planBeg + 1;
-  }
-  Int_t sectBeg = 0;
-  Int_t sectEnd = AliTRDgeometry::Nsect();
+  //
+  Double_t padSignal[kNsig];   
+  Double_t clusterSignal[kNclus];
+  Double_t clusterPads[kNclus];   
+  Int_t    clusterDigit[kNclus];
+  Int_t    clusterTracks[kNtrack];   
+
+  Int_t    chamBeg = 0;
+  Int_t    chamEnd = AliTRDgeometry::Ncham();
+  Int_t    planBeg = 0;
+  Int_t    planEnd = AliTRDgeometry::Nplan();
+  Int_t    sectBeg = 0;
+  Int_t    sectEnd = AliTRDgeometry::Nsect();
 
   // Start clustering in every chamber
   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
 
-        if (fTRD->GetSensSector() >= 0) {
-          Int_t sens1 = fTRD->GetSensSector();
-          Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
-          sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) 
-                 * AliTRDgeometry::Nsect();
-          if (sens1 < sens2) {
-            if ((isect < sens1) || (isect >= sens2)) continue;
-         }
-          else {
-            if ((isect < sens1) && (isect >= sens2)) continue;
-         }
-       }
-
         Int_t idet = geo->GetDetector(iplan,icham,isect);
 
         Int_t nClusters      = 0;
@@ -338,14 +226,18 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
         Int_t nClusters5pad  = 0;
         Int_t nClustersLarge = 0;
 
-        printf("AliTRDclusterizerV1::MakeCluster -- ");
-        printf("Analyzing chamber %d, plane %d, sector %d.\n"
-              ,icham,iplan,isect);
+        if (fVerbose > 0) {
+          printf("<AliTRDclusterizerV1::MakeCluster> ");
+          printf("Analyzing chamber %d, plane %d, sector %d.\n"
+                ,icham,iplan,isect);
+       }
+
+        Int_t    nRowMax     = fPar->GetRowMax(iplan,icham,isect);
+        Int_t    nColMax     = fPar->GetColMax(iplan);
+        Int_t    nTimeBefore = fPar->GetTimeBefore();
+        Int_t    nTimeTotal  = fPar->GetTimeTotal();  
 
-        Int_t nRowMax     = geo->GetRowMax(iplan,icham,isect);
-        Int_t nColMax     = geo->GetColMax(iplan);
-        Int_t nTimeBefore = geo->GetTimeBefore();
-        Int_t nTimeTotal  = geo->GetTimeTotal();  
+        AliTRDpadPlane *padPlane = fPar->GetPadPlane(iplan,icham);
 
         // Get the digits
         digits = fDigitsManager->GetDigits(idet);
@@ -359,19 +251,28 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
 
         // Loop through the chamber and find the maxima 
         for ( row = 0;  row <  nRowMax;    row++) {
-          for ( col = 2;  col <  nColMax;    col++) {
+         for ( col = 2;  col <  nColMax;    col++) {
+           //for ( col = 4;  col <  nColMax-2;    col++) {
             for (time = 0; time < nTimeTotal; time++) {
 
               Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col  ,time));
               Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
               Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
  
+//           // Look for the maximum
+//               if (signalM >= maxThresh) {
+//                 if (((signalL >= sigThresh) &&
+//                      (signalL <  signalM))  ||
+//                     ((signalR >= sigThresh) &&
+//                      (signalR <  signalM))) {
+//                   // Maximum found, mark the position by a negative signal
+//                   digits->SetDataUnchecked(row,col-1,time,-signalM);
+//             }
+//           }
              // Look for the maximum
               if (signalM >= maxThresh) {
-                if (((signalL >= sigThresh) &&
-                     (signalL <  signalM))  ||
-                    ((signalR >= sigThresh) &&
-                     (signalR <  signalM))) {
+                if ( (TMath::Abs(signalL)<=signalM) && (TMath::Abs(signalR)<=signalM) && 
+                    (TMath::Abs(signalL)+TMath::Abs(signalR))>sigThresh ) {
                   // Maximum found, mark the position by a negative signal
                   digits->SetDataUnchecked(row,col-1,time,-signalM);
                }
@@ -442,10 +343,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                   break;
                };
 
-               // Don't analyze large clusters
-                //if (iType == 4) continue;
-
-                // Look for 5 pad cluster with minimum in the middle
+                // Look for 5 pad cluster with minimum in the middle
                 Bool_t fivePadCluster = kFALSE;
                 if (col < nColMax-3) {
                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
@@ -468,7 +366,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                // of the cluster which remains from a previous unfolding
                 if (iUnfold) {
                   clusterSignal[0] *= ratioLeft;
-                  iType   = 3;
+                  iType   = 5;
                   iUnfold = 0;
                }
 
@@ -481,103 +379,75 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                   }
                   // Unfold the two maxima and set the signal on 
                   // the overlapping pad to the ratio
-                  ratioRight        = Unfold(kEpsilon,padSignal);
+                  ratioRight        = Unfold(kEpsilon,iplan,padSignal);
                   ratioLeft         = 1.0 - ratioRight; 
                   clusterSignal[2] *= ratioRight;
-                  iType   = 3;
+                  iType   = 5;
                   iUnfold = 1;
                 }
 
-                Float_t clusterCharge = clusterSignal[0]
-                                      + clusterSignal[1]
-                                      + clusterSignal[2];
+
+                Double_t clusterCharge = clusterSignal[0]
+                                       + clusterSignal[1]
+                                       + clusterSignal[2];
                 
                // The position of the cluster
                 clusterPads[0] = row + 0.5;
                // Take the shift of the additional time bins into account
                 clusterPads[2] = time - nTimeBefore + 0.5;
 
-                if (fUseLUT) {
-
+                if (fPar->LUTOn()) {
                  // Calculate the position of the cluster by using the
                  // lookup table method
-                  Float_t ratioLUT;
-                  Float_t signLUT;
-                  Float_t lut = 0.0;
-                  if (clusterSignal[0] > clusterSignal[2]) {
-                    ratioLUT = clusterSignal[0] / clusterSignal[1];
-                    signLUT  = -1.0;
-                 }
-                  else {
-                    ratioLUT = clusterSignal[2] / clusterSignal[1];
-                    signLUT  =  1.0;
-                 }
-                  if      (ratioLUT < kLUTmin) {
-                    lut = 0.0;
-                 }
-                  else if (ratioLUT > kLUTmax) {
-                    lut = 0.5;
-                 }
-                  else {
-                    Int_t indexLUT = TMath::Nint ((kNlut-1) * (ratioLUT - kLUTmin)  
-                                                           / (kLUTmax  - kLUTmin)); 
-                    lut = fLUT[indexLUT];
-                 }
-                  clusterPads[1] = col + 0.5 + signLUT * lut;
-
+                  clusterPads[1] =
+                                  fPar->LUTposition(iplan,clusterSignal[0]
+                                                          ,clusterSignal[1]
+                                                         ,clusterSignal[2]);
                }
                else {
-
                  // Calculate the position of the cluster by using the
                  // center of gravity method
-                  clusterPads[1] = col + 0.5 
-                                 + (clusterSignal[2] - clusterSignal[0]) 
-                                / clusterCharge;
+                 for (Int_t i=0;i<5;i++) padSignal[i]=0;
+                 padSignal[2] = TMath::Abs(digits->GetDataUnchecked(row,col,time));   // central  pad
+                 padSignal[1] = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); // left     pad
+                 padSignal[3] = TMath::Abs(digits->GetDataUnchecked(row,col+1,time)); // right    pad
+                 if (col>2 &&TMath::Abs(digits->GetDataUnchecked(row,col-2,time)<padSignal[1])){
+                   padSignal[0] = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
+                 }
+                 if (col<nColMax-3 &&TMath::Abs(digits->GetDataUnchecked(row,col+2,time)<padSignal[3])){
+                   padSignal[4] = TMath::Abs(digits->GetDataUnchecked(row,col+2,time));
+                 }               
+                 clusterPads[1] =  GetCOG(padSignal);
+                 Double_t check = fPar->LUTposition(iplan,clusterSignal[0]
+                                                          ,clusterSignal[1]
+                                                         ,clusterSignal[2]);
+                 //              Float_t diff = clusterPads[1] -  check;
 
                }
 
-                Float_t clusterSigmaY2 = (clusterSignal[2] + clusterSignal[0]) / clusterCharge 
-                                       - (clusterPads[1]-col-0.5) * (clusterPads[1]-col-0.5);
-
-                // Correct for ExB displacement
-                if (digitizer->GetExB()) { 
-                  Int_t   local_time_bin = (Int_t) clusterPads[2];
-                  Float_t driftLength    = local_time_bin * timeBinSize + kAmWidth;
-                  Float_t colSize        = geo->GetColPadSize(iplan);
-                  Float_t deltaY         = omegaTau*driftLength/colSize;
-                  clusterPads[1]         = clusterPads[1] - deltaY;
-                }
-                                       
-                if (fVerbose) {
-                  printf("-----------------------------------------------------------\n");
-                  printf("Create cluster no. %d\n",nClusters);
-                  printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
-                                                                   ,clusterPads[1]
-                                                                    ,clusterPads[2]);
-                  printf("Indices: %d, %d, %d\n",clusterDigit[0]
-                                               ,clusterDigit[1]
-                                                ,clusterDigit[2]);
-                  printf("Total charge = %f\n",clusterCharge);
-                  printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
-                                                   ,clusterTracks[1]
-                                                    ,clusterTracks[2]);
-                  printf("        pad1 %d, %d, %d\n",clusterTracks[3]
-                                                   ,clusterTracks[4]
-                                                    ,clusterTracks[5]);
-                  printf("        pad2 %d, %d, %d\n",clusterTracks[6]
-                                                   ,clusterTracks[7]
-                                                    ,clusterTracks[8]);
-                  printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
-                }
-
+                Double_t q0 = clusterSignal[0];
+               Double_t q1 = clusterSignal[1];
+                Double_t q2 = clusterSignal[2];
+                Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
+                                          (clusterCharge*clusterCharge);
+
+               // Calculate the position and the error
+                Double_t colSize = padPlane->GetColSize(col);
+                Double_t rowSize = padPlane->GetRowSize(row);
+                Double_t clusterPos[3];
+               clusterPos[0] = padPlane->GetColPos(col) + (clusterPads[1]-0.5)*colSize;  // MI change
+               clusterPos[1] = padPlane->GetRowPos(row) -0.5*rowSize; //MI change
+                clusterPos[2] = clusterPads[2];
+                Double_t clusterSig[2];
+                clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
+                clusterSig[1] = rowSize * rowSize / 12.;                                       
                 // Add the cluster to the output array 
-                fTRD->AddCluster(clusterPads
-                                ,clusterDigit
-                                ,idet
-                                ,clusterCharge
-                                ,clusterTracks
-                               ,clusterSigmaY2
-                                ,iType);
+                AddCluster(clusterPos
+                         ,idet
+                         ,clusterCharge
+                         ,clusterTracks
+                         ,clusterSig
+                          ,iType,clusterPads[1]);
 
               }
             } 
@@ -587,35 +457,39 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
        // Compress the arrays
         digits->Compress(1,0);
         track0->Compress(1,0);
-        track1->Compress(1,0);
+       track1->Compress(1,0);
         track2->Compress(1,0);
 
         // Write the cluster and reset the array
        WriteClusters(idet);
-       fTRD->ResetRecPoints();
-
-        printf("AliTRDclusterizerV1::MakeCluster -- ");
-        printf("Found %d clusters in total.\n"
-              ,nClusters);
-        printf("                                    2pad:  %d\n",nClusters2pad);
-        printf("                                    3pad:  %d\n",nClusters3pad);
-        printf("                                    4pad:  %d\n",nClusters4pad);
-        printf("                                    5pad:  %d\n",nClusters5pad);
-        printf("                                    Large: %d\n",nClustersLarge);
-
+       ResetRecPoints();
       }    
     }      
   }        
 
-  printf("AliTRDclusterizerV1::MakeCluster -- ");
-  printf("Done.\n");
+  if (fVerbose > 0) {
+    printf("<AliTRDclusterizerV1::MakeCluster> ");
+    printf("Done.\n");
+  }
 
   return kTRUE;
 
 }
 
+Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
+{
+  //
+  // get COG position
+  // used for clusters with more than 3 pads - where LUT not applicable
+  Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4];
+  Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum;
+  return res;            
+}
+
+
+
 //_____________________________________________________________________________
-Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
+Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
 {
   //
   // Method to unfold neighbouring maxima.
@@ -624,13 +498,15 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
   // The resulting ratio is then returned to the calling method.
   //
 
-  Int_t   itStep            = 0;      // Count iteration steps
+  Int_t   irc                = 0;
+  Int_t   itStep             = 0;      // Count iteration steps
 
-  Float_t ratio             = 0.5;    // Start value for ratio
-  Float_t prevRatio         = 0;      // Store previous ratio
+  Double_t ratio             = 0.5;    // Start value for ratio
+  Double_t prevRatio         = 0;      // Store previous ratio
 
-  Float_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
-  Float_t newRightSignal[3] = {0};    // Array to store right cluster signal
+  Double_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
+  Double_t newRightSignal[3] = {0};    // Array to store right cluster signal
+  Double_t newSignal[3]      = {0};
 
   // Start the iteration
   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
@@ -639,26 +515,23 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
     prevRatio = ratio;
 
     // Cluster position according to charge ratio
-    Float_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
-                     / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
-    Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
-                     / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
+    Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
+                      / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
+    Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
+                      / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 
     // Set cluster charge ratio
-    Float_t ampLeft  = padSignal[1] / PadResponse(0 - maxLeft );
-    Float_t ampRight = padSignal[3] / PadResponse(0 - maxRight);
+    irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
+    Double_t ampLeft  = padSignal[1] / newSignal[1];
+    irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
+    Double_t ampRight = padSignal[3] / newSignal[1];
 
     // Apply pad response to parameters
-    newLeftSignal[0]  = ampLeft  * PadResponse(-1 - maxLeft);
-    newLeftSignal[1]  = ampLeft  * PadResponse( 0 - maxLeft);
-    newLeftSignal[2]  = ampLeft  * PadResponse( 1 - maxLeft);
-
-    newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
-    newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
-    newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
+    irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
+    irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
 
     // Calculate new overlapping ratio
-    ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] / 
+    ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / 
                           (newLeftSignal[2] + newRightSignal[0]));
 
   }
@@ -667,24 +540,3 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
 
 }
 
-//_____________________________________________________________________________
-Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
-{
-  //
-  // The pad response for the chevron pads. 
-  // We use a simple Gaussian approximation which should be good
-  // enough for our purpose.
-  // Updated for new PRF 1/5/01.
-  //
-
-  // The parameters for the response function
-  const Float_t kA  =  0.8303; 
-  const Float_t kB  = -0.00392;
-  const Float_t kC  =  0.472 * 0.472;
-  const Float_t kD  =  2.19;
-
-  Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));
-
-  return (pr);
-
-}