]> 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 39990a4eb2b2332dee38dfc89334e6c1753d5834..c5fff1f5e931a23ddb8536e8a19077d93ee39ac7 100644 (file)
@@ -37,6 +37,7 @@
 #include "AliTRDdataArrayI.h"
 #include "AliTRDdigitsManager.h"
 #include "AliTRDparameter.h"
+#include "AliTRDpadPlane.h"
 
 ClassImp(AliTRDclusterizerV1)
 
@@ -161,10 +162,12 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
     printf("<AliTRDclusterizerV1::MakeCluster> ");
     printf("Create the default parameter object.\n");
   }
+  fPar->Init();
 
-  Float_t timeBinSize = fPar->GetTimeBinSize();
+  //Float_t timeBinSize = fPar->GetDriftVelocity()
+  //                    / fPar->GetSamplingFrequency();
   // Half of ampl.region
-  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
+  //  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
 
   Float_t omegaTau = fPar->GetOmegaTau();
   if (fVerbose > 0) {
@@ -183,7 +186,6 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   Int_t maxThresh = fPar->GetClusMaxThresh();   
   // Threshold value for the digit signal
   Int_t sigThresh = fPar->GetClusSigThresh();   
-
   // Iteration limit for unfolding procedure
   const Float_t kEpsilon = 0.01;             
 
@@ -191,24 +193,24 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   const Int_t   kNsig    = 5;
   const Int_t   kNtrack  = 3 * kNclus;
 
-  Int_t   iType          = 0;
-  Int_t   iUnfold        = 0;
-
-  Float_t ratioLeft      = 1.0;
-  Float_t ratioRight     = 1.0;
+  Int_t    iType         = 0;
+  Int_t    iUnfold       = 0;  
+  Double_t ratioLeft     = 1.0;
+  Double_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();
-  Int_t planBeg = 0;
-  Int_t planEnd = AliTRDgeometry::Nplan();
-  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++) {
@@ -230,15 +232,12 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                 ,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     = fPar->GetRowMax(iplan,icham,isect);
+        Int_t    nColMax     = fPar->GetColMax(iplan);
+        Int_t    nTimeBefore = fPar->GetTimeBefore();
+        Int_t    nTimeTotal  = fPar->GetTimeTotal();  
 
-        Float_t row0        = fPar->GetRow0(iplan,icham,isect);
-        Float_t col0        = fPar->GetCol0(iplan);
-        Float_t rowSize     = fPar->GetRowPadSize(iplan,icham,isect);
-        Float_t colSize     = fPar->GetColPadSize(iplan);
+        AliTRDpadPlane *padPlane = fPar->GetPadPlane(iplan,icham);
 
         // Get the digits
         digits = fDigitsManager->GetDigits(idet);
@@ -252,20 +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 = 4;  col <  nColMax-2;    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);
                }
@@ -336,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) {
@@ -362,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;
                }
 
@@ -378,13 +382,14 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                   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;
@@ -392,78 +397,57 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                 clusterPads[2] = time - nTimeBefore + 0.5;
 
                 if (fPar->LUTOn()) {
-
                  // Calculate the position of the cluster by using the
                  // lookup table method
-                  clusterPads[1] = col + 0.5
-                                 + fPar->LUTposition(iplan,clusterSignal[0]
+                  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 q0 = clusterSignal[0];
-               Float_t q1 = clusterSignal[1];
-                Float_t q2 = clusterSignal[2];
-                Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
-                                         (clusterCharge*clusterCharge);
-
-                // Correct for ExB displacement
-                if (fPar->ExBOn()) { 
-                  Int_t   local_time_bin = (Int_t) clusterPads[2];
-                  Float_t driftLength    = local_time_bin * timeBinSize + kAmWidth;
-                  Float_t colSize        = fPar->GetColPadSize(iplan);
-                  Float_t deltaY         = omegaTau*driftLength/colSize;
-                  clusterPads[1]         = clusterPads[1] - deltaY;
-                }
-                                       
-                if (fVerbose > 1) {
-                  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
-                Float_t clusterPos[3];
-                clusterPos[0] = clusterPads[1] * colSize + col0;
-                clusterPos[1] = clusterPads[0] * rowSize + row0;
+                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];
-                Float_t clusterSig[2];
+                Double_t clusterSig[2];
                 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
-                clusterSig[1] = rowSize * rowSize / 12.;
-
+                clusterSig[1] = rowSize * rowSize / 12.;                                       
                 // Add the cluster to the output array 
                 AddCluster(clusterPos
-                          ,idet
-                          ,clusterCharge
-                          ,clusterTracks
-                          ,clusterSig
-                          ,iType);
+                         ,idet
+                         ,clusterCharge
+                         ,clusterTracks
+                         ,clusterSig
+                          ,iType,clusterPads[1]);
 
               }
             } 
@@ -473,24 +457,12 @@ 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);
        ResetRecPoints();
-
-        if (fVerbose > 0) {
-          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);
-       }
-
       }    
     }      
   }        
@@ -504,8 +476,20 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
 
 }
 
+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, Int_t plane, Float_t* padSignal)
+Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
 {
   //
   // Method to unfold neighbouring maxima.
@@ -514,15 +498,15 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal
   // The resulting ratio is then returned to the calling method.
   //
 
-  Int_t   irc               = 0;
-  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
-  Float_t newSignal[3]      = {0};
+  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)) {
@@ -531,23 +515,23 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, 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
     irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
-    Float_t ampLeft  = padSignal[1] / newSignal[1];
+    Double_t ampLeft  = padSignal[1] / newSignal[1];
     irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
-    Float_t ampRight = padSignal[3] / newSignal[1];
+    Double_t ampRight = padSignal[3] / newSignal[1];
 
     // Apply pad response to parameters
     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]));
 
   }