]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDclusterizerV1.cxx
Merge with TRD-develop
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
index afb19bba02bee660ff7eb93f2ef1872970bccc0f..e76f11e6b709ceca4c2c871d21de5116c3493ef5 100644 (file)
 
 /*
 $Log$
+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
+
 */
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -24,13 +60,20 @@ $Log$
 ///////////////////////////////////////////////////////////////////////////////
 
 #include <TF1.h>
+#include <TTree.h>
+#include <TH1.h>
 
+#include "AliRun.h"
+
+#include "AliTRD.h"
 #include "AliTRDclusterizerV1.h"
 #include "AliTRDmatrix.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDdigitizer.h"
 #include "AliTRDrecPoint.h"
-#include "AliTRDdataArray.h"
+#include "AliTRDdataArrayF.h"
+#include "AliTRDdataArrayI.h"
+#include "AliTRDdigitsManager.h"
 
 ClassImp(AliTRDclusterizerV1)
 
@@ -41,7 +84,7 @@ AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
   // AliTRDclusterizerV1 default constructor
   //
 
-  fDigitsArray = NULL;
+  fDigitsManager = NULL;
 
 }
 
@@ -53,23 +96,64 @@ AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title
   // AliTRDclusterizerV1 default constructor
   //
 
-  fDigitsArray = NULL;
+  fDigitsManager = new AliTRDdigitsManager();
 
   Init();
 
 }
 
+//_____________________________________________________________________________
+AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
+{
+  //
+  // AliTRDclusterizerV1 copy constructor
+  //
+
+  ((AliTRDclusterizerV1 &) c).Copy(*this);
+
+}
+
 //_____________________________________________________________________________
 AliTRDclusterizerV1::~AliTRDclusterizerV1()
 {
+  //
+  // AliTRDclusterizerV1 destructor
+  //
 
-  if (fDigitsArray) {
-    fDigitsArray->Delete();
-    delete fDigitsArray;
+  if (fDigitsManager) {
+    delete fDigitsManager;
   }
 
 }
 
+//_____________________________________________________________________________
+AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
+  return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDclusterizerV1::Copy(TObject &c)
+{
+  //
+  // Copy function
+  //
+
+  ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
+  ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
+  ((AliTRDclusterizerV1 &) c).fClusMethod    = fClusMethod;
+  ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
+
+  AliTRDclusterizer::Copy(c);
+
+}
+
 //_____________________________________________________________________________
 void AliTRDclusterizerV1::Init()
 {
@@ -97,16 +181,13 @@ Bool_t AliTRDclusterizerV1::ReadDigits()
     return kFALSE;
   }
 
-  // Create a new segment array for the digits 
-  fDigitsArray = new AliTRDsegmentArray(kNsect*kNplan*kNcham);
-
   // Read in the digit arrays
-  return (fDigitsArray->LoadArray("TRDdigits"));  
+  return (fDigitsManager->ReadDigits());  
 
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDclusterizerV1::MakeCluster()
+Bool_t AliTRDclusterizerV1::MakeClusters()
 {
   //
   // Generates the cluster.
@@ -115,20 +196,23 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
   Int_t row, col, time;
 
   // Get the pointer to the detector class and check for version 1
-  AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
-  if (TRD->IsVersion() != 1) {
+  AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
+  if (trd->IsVersion() != 1) {
     printf("AliTRDclusterizerV1::MakeCluster -- ");
     printf("TRD must be version 1 (slow simulator).\n");
     return kFALSE; 
   }
 
   // Get the geometry
-  AliTRDgeometry *Geo = TRD->GetGeometry();
+  AliTRDgeometry *geo = trd->GetGeometry();
 
   printf("AliTRDclusterizerV1::MakeCluster -- ");
   printf("Start creating clusters.\n");
 
-  AliTRDdataArray *Digits;
+  AliTRDdataArrayI *digits;
+  AliTRDdataArrayI *track0;
+  AliTRDdataArrayI *track1;
+  AliTRDdataArrayI *track2; 
 
   // Parameters
   Float_t maxThresh        = fClusMaxThresh;   // threshold value for maximum
@@ -136,45 +220,54 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
   Int_t   clusteringMethod = fClusMethod;      // clustering method option (for testing)
 
   // Iteration limit for unfolding procedure
-  const Float_t epsilon = 0.01;             
+  const Float_t kEpsilon = 0.01;             
 
-  const Int_t   nClus   = 3;  
-  const Int_t   nSig    = 5;
+  const Int_t   kNclus   = 3;  
+  const Int_t   kNsig    = 5;
 
   Int_t chamBeg = 0;
-  Int_t chamEnd = kNcham;
-  if (TRD->GetSensChamber() >= 0) {
-    chamBeg = TRD->GetSensChamber();
-    chamEnd = chamEnd + 1;
+  Int_t chamEnd = AliTRDgeometry::Ncham();
+  if (trd->GetSensChamber()  >= 0) {
+    chamBeg = trd->GetSensChamber();
+    chamEnd = chamBeg + 1;
   }
   Int_t planBeg = 0;
-  Int_t planEnd = kNplan;
-  if (TRD->GetSensPlane()   >= 0) {
-    planBeg = TRD->GetSensPlane();
+  Int_t planEnd = AliTRDgeometry::Nplan();
+  if (trd->GetSensPlane()    >= 0) {
+    planBeg = trd->GetSensPlane();
     planEnd = planBeg + 1;
   }
   Int_t sectBeg = 0;
-  Int_t sectEnd = kNsect;
-  if (TRD->GetSensSector()  >= 0) {
-    sectBeg = TRD->GetSensSector();
-    sectEnd = sectBeg + 1;
-  }
+  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++) {
 
-        Int_t idet = Geo->GetDetector(iplan,icham,isect);
+        if (trd->GetSensSector() >= 0) {
+          Int_t sens1 = trd->GetSensSector();
+          Int_t sens2 = sens1 + trd->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;
         printf("AliTRDclusterizerV1::MakeCluster -- ");
         printf("Analyzing chamber %d, plane %d, sector %d.\n"
                ,icham,iplan,isect);
 
-        Int_t   nRowMax  = Geo->GetRowMax(iplan,icham,isect);
-        Int_t   nColMax  = Geo->GetColMax(iplan);
-        Int_t   nTimeMax = Geo->GetTimeMax();
+        Int_t   nRowMax  = geo->GetRowMax(iplan,icham,isect);
+        Int_t   nColMax  = geo->GetColMax(iplan);
+        Int_t   nTimeMax = geo->GetTimeMax();
 
         // Create a detector matrix to keep maxima
         AliTRDmatrix *digitMatrix  = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
@@ -183,16 +276,27 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
         AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
                                                      ,isect,icham,iplan);
 
+        // Create a matrix for track indexes
+        AliTRDmatrix *trackMatrix  = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
+                                                     ,isect,icham,iplan);
+
         // Read in the digits
-        Digits = (AliTRDdataArray *) fDigitsArray->At(idet);
+        digits = fDigitsManager->GetDigits(idet);
+        track0 = fDigitsManager->GetDictionary(idet,0);
+        track1 = fDigitsManager->GetDictionary(idet,1);
+        track2 = fDigitsManager->GetDictionary(idet,2); 
 
         // Loop through the detector pixel
         for (time = 0; time < nTimeMax; time++) {
           for ( col = 0;  col <  nColMax;  col++) {
             for ( row = 0;  row <  nRowMax;  row++) {
 
-              Int_t signal = Digits->GetData(row,col,time);
-              Int_t index  = Digits->GetIndex(row,col,time);
+              Int_t signal = digits->GetData(row,col,time);
+              Int_t index  = digits->GetIndex(row,col,time);
+              Int_t t[3] = {-1};
+              t[0] = track0->GetData(row,col,time) - 1;
+              t[1] = track1->GetData(row,col,time) - 1;
+              t[2] = track2->GetData(row,col,time) - 1;  
 
               // Fill the detector matrix
               if (signal > signalThresh) {
@@ -200,8 +304,10 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
                 digitMatrix->SetSignal(row,col,time,signal);
                // Store the digits number
                 digitMatrix->AddTrack(row,col,time,index);
+                for(Int_t i = 0; i < 3; i++) {
+                  trackMatrix->AddTrack(row,col,time,t[i]);
+                }      
               }
-
            }
          }
        }
@@ -237,31 +343,38 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
                   && (digitMatrix->GetSignal(row,col,time) > maxThresh)) {
 
                 // Ratio resulting from unfolding
-                Float_t ratio                =  0;    
+                Float_t ratio                 =  0;    
                 // Signals on max and neighbouring pads
-                Float_t padSignal[nSig]      = {0};   
+                Float_t padSignal[kNsig]      = {0};   
                 // Signals from cluster
-                Float_t clusterSignal[nClus] = {0};
+                Float_t clusterSignal[kNclus] = {0};
                 // Cluster pad info
-                Float_t clusterPads[nClus]   = {0};   
+                Float_t clusterPads[kNclus]   = {0};   
                 // Cluster digit info
-                Int_t   clusterDigit[nClus]  = {0};
+                Int_t   clusterDigit[kNclus]  = {0};
+                // Cluster MC tracks info
+                const Int_t nt = kNclus*3;
+                Int_t clusterTracks[nt] = {-1};   
 
-                for (Int_t iPad = 0; iPad < nClus; iPad++) {
+                Int_t iPad;
+                for (iPad = 0; iPad < kNclus; iPad++) {
                   clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
                   clusterDigit[iPad]  = digitMatrix->GetTrack(row,col-1+iPad,time,0);
+                  for (Int_t j = 0; j < 3; j++) {
+                    clusterTracks[iPad*3+j] = trackMatrix->GetTrack(row,col-1+iPad,time,j);
+                  }  
                 }
 
                 // neighbouring maximum on right side?
                 if (col < nColMax - 2) {
                   if (maximaMatrix->GetSignal(row,col + 2,time) > 0) {
 
-                    for (Int_t iPad = 0; iPad < 5; iPad++) {
+                    for (iPad = 0; iPad < 5; iPad++) {
                       padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
                     }
 
                     // unfold:
-                    ratio = Unfold(epsilon, padSignal);
+                    ratio = Unfold(kEpsilon, padSignal);
 
                     // set signal on overlapping pad to ratio
                     clusterSignal[2] *= ratio;
@@ -274,10 +387,27 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
                 case 1:
                   // method 1: simply center of mass
                   clusterPads[0] = row + 0.5;
-                  clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) /
-                                   (clusterSignal[1] + clusterSignal[2] + clusterSignal[3]);
+                  clusterPads[1] = col + 0.5 + (clusterSignal[2] - clusterSignal[0]) /
+                                   (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]);
                   clusterPads[2] = time + 0.5;
 
+
+/*               printf("col = %d, left = %f, center = %f, right = %f, 
+                         final =%f \n", col, 
+                         digitMatrix->GetSignal(row,col-1,time),
+                         digitMatrix->GetSignal(row,col,time),
+                         digitMatrix->GetSignal(row,col+1,time),
+                         clusterPads[1]);
+
+                 printf("col = %d, sig(0) = %f, sig(1) = %f, sig(2) = %f, 
+                         final =%f \n", col, 
+                         clusterSignal[0],
+                         clusterSignal[1],
+                         clusterSignal[2],
+                         clusterPads[1]);
+
+*/
+
                   nClusters++;
                   break;
                 case 2:
@@ -307,7 +437,7 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
                                         + clusterSignal[2];
 
                 // Add the cluster to the output array 
-                TRD->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);
+                trd->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge,clusterTracks);
 
               }
             }  // time
@@ -317,24 +447,28 @@ Bool_t AliTRDclusterizerV1::MakeCluster()
         printf("AliTRDclusterizerV1::MakeCluster -- ");
         printf("Number of clusters found: %d\n",nClusters);
 
+       WriteClusters(idet);
+       trd->ResetRecPoints();
+
         delete digitMatrix;
         delete maximaMatrix;
+        delete trackMatrix;
 
       }          // isect
     }            // iplan
   }              // icham
 
-  printf("AliTRDclusterizerV1::MakeCluster -- ");
-  printf("Total number of points found: %d\n"
-        ,TRD->RecPoints()->GetEntries());
+//   printf("AliTRDclusterizerV1::MakeCluster -- ");
+//   printf("Total number of points found: %d\n"
+//         ,trd->RecPoints()->GetEntries());
 
-  // Get the pointer to the cluster branch
-  TTree *ClusterTree = gAlice->TreeR(); 
+//   // Get the pointer to the cluster branch
+//   TTree *clusterTree = gAlice->TreeR(); 
 
-  // Fill the cluster-branch
-  printf("AliTRDclusterizerV1::MakeCluster -- ");
-  printf("Fill the cluster tree.\n");
-  ClusterTree->Fill();
+//   // Fill the cluster-branch
+//   printf("AliTRDclusterizerV1::MakeCluster -- ");
+//   printf("Fill the cluster tree.\n");
+//   clusterTree->Fill();
   printf("AliTRDclusterizerV1::MakeCluster -- ");
   printf("Done.\n");
 
@@ -404,12 +538,12 @@ Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
   //
 
   // The parameters for the response function
-  const Float_t aa  =  0.8872;
-  const Float_t bb  = -0.00573;
-  const Float_t cc  =  0.454;
-  const Float_t cc2 =  cc*cc;
+  const Float_t kA  =  0.8872;
+  const Float_t kB  = -0.00573;
+  const Float_t kC  =  0.454;
+  const Float_t kC2 =  kC*kC;
 
-  Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2)));
+  Float_t pr = kA * (kB + TMath::Exp(-x*x / (2. * kC2)));
 
   return (pr);