the MIXT geometry (IHEP+GPS2) has been introduced
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
index 397f5e9..156b57c 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
+/* $Log:
+   1 October 2000. Yuri Kharlov:
+     AreNeighbours()
+     PPSD upper layer is considered if number of layers>1
+
+   18 October 2000. Yuri Kharlov:
+     AliPHOSClusterizerv1()
+     CPV clusterizing parameters added
+
+     MakeClusters()
+     After first PPSD digit remove EMC digits only once
+*/
+
 //_________________________________________________________________________
-// A brief description of the class
-//*-- Author : Yves Schutz  SUBATECH 
+//  Implementation version 1 of the clusterization algorithm 
+// 
+//*-- Author: Yves Schutz (SUBATECH) 
 //////////////////////////////////////////////////////////////////////////////
 
+#include <assert.h>
+
 // --- ROOT system ---
 
 #include "TMath.h" 
 
 // --- Standard library ---
 
-#include <iostream>
+#include <iostream.h>
 
 // --- AliRoot header files ---
 
@@ -32,6 +50,7 @@
 #include "AliPHOSDigit.h"
 #include "AliPHOSEmcRecPoint.h"
 #include "AliPHOSPpsdRecPoint.h"
+#include "AliPHOSCpvRecPoint.h"
 #include "AliPHOSv0.h" 
 #include "AliRun.h" 
 
@@ -40,33 +59,42 @@ ClassImp(AliPHOSClusterizerv1)
 //____________________________________________________________________________
 AliPHOSClusterizerv1::AliPHOSClusterizerv1()
 {
+  // default ctor (to be used)
+
   fA                       = 0.;
   fB                       = 0.01 ;
+  fGeom                    = AliPHOSGeometry::GetInstance();
   fNumberOfEmcClusters     = 0 ; 
   fNumberOfPpsdClusters    = 0 ; 
   fEmcClusteringThreshold  = 0.1;   
   fEmcEnergyThreshold      = 0.01;    
-  fPpsdClusteringThreshold = 0.00000015; 
-  fPpsdEnergyThreshold     = 0.0000001;  
+  fPpsdClusteringThreshold = 0.0;
+  fPpsdEnergyThreshold     = 0.1;  
+  fCpvClusteringThreshold  = 0.0;
+  fCpvEnergyThreshold      = 0.1;  
   fW0                      = 4.5 ;
   fLocMaxCut               = 0.06 ;
+  fW0CPV                   = 4.5 ;
+  fLocMaxCutCPV            = 0.06 ;
 }
 
 //____________________________________________________________________________
 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)
 {
-  // neigbours are defined as digits having at least common vertex
-  // The order of A and B in AreNeighbours(A,B) is important: first (A) should be digit 
-  // in cluster, which compared with digits not clusterized yet  
-  Int_t rv = 0 ; 
+  // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
+  //                                       = 1 are neighbour
+  //                                       = 2 are not neighbour but do not continue searching
+  // neighbours are defined as digits having at least common vertex
+  // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster 
+  //                                      which is compared to a digit (d2)  not yet in a cluster  
 
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
+  Int_t rv = 0 ; 
 
   Int_t relid1[4] ; 
-  geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
+  fGeom->AbsToRelNumbering(d1->GetId(), relid1) ; 
 
   Int_t relid2[4] ; 
-  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
+  fGeom->AbsToRelNumbering(d2->GetId(), relid2) ; 
  
   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module 
     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
@@ -87,25 +115,26 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)
       rv=2 ;
 
   }
-  
+
+  //Do NOT clusterize upper PPSD  
+  if( IsInPpsd(d1) && IsInPpsd(d2) &&
+     relid1[1] > 0                 &&
+     relid1[1] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ;
+
   return rv ; 
 }
 
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::FillandSort(const DigitsList * dl, TObjArray * tl) 
 {
-  // copies the digits with energy above thershold and sorts the list
+  // Copies the digits with energy above thershold and sorts the list
   // according to increasing Id number
 
-  cout << "HOLA 1" << endl;
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
   Int_t relid[4] ;  
   
   TIter next(dl) ; 
   AliPHOSDigit * digit ;
   
-
   while ( (digit = (AliPHOSDigit *)next()) ) { 
 
 //     cout << " clusterizerv1 " << endl ;
@@ -116,7 +145,7 @@ void AliPHOSClusterizerv1::FillandSort(const DigitsList * dl, TObjArray * tl)
 
     Int_t id    = digit->GetId() ; 
     Float_t ene = Calibrate(digit->GetAmp()) ; 
-    geom->AbsToRelNumbering(id, relid) ;
+    fGeom->AbsToRelNumbering(id, relid) ;
     if(relid[1]==0){ // EMC
       if ( ene > fEmcEnergyThreshold )
        tl->Add(digit) ;
@@ -128,45 +157,77 @@ void AliPHOSClusterizerv1::FillandSort(const DigitsList * dl, TObjArray * tl)
     }
 
   }
-  cout << tl << endl;
-
-  TIter next2(tl) ; Int_t iii=0;
-  cout << "Entries is " << tl->GetEntries() << tl->GetLast() << endl;
-  while ( (digit = (AliPHOSDigit *)next2()) ) { 
-    cout << iii++ << "   "  << digit << " Amplitud is  " << digit->GetAmp() << endl;
-  }
-  cout << "HOLA de nuevo " << tl << endl;
-  tl->Print();
-  cout << "sort" << endl ;
   tl->Sort() ; 
 }
 
 //____________________________________________________________________________
 void AliPHOSClusterizerv1:: GetNumberOfClustersFound(Int_t * numb) 
 {
- numb[0] = fNumberOfEmcClusters ; 
- numb[1] = fNumberOfPpsdClusters ; 
+  // Fills numb with the number of EMC  (numb[0]) clusters found
+  //                               PPSD (numb[1]) clusters found
+
+  numb[0] = fNumberOfEmcClusters ; 
+  numb[1] = fNumberOfPpsdClusters ; 
 }
 
 //____________________________________________________________________________
 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) 
 {
+  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
   Bool_t rv = kFALSE ; 
 
-  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;  
+  Int_t relid[4] ; 
+  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+
+  if ( relid[1] == 0  ) rv = kTRUE; 
+
+  return rv ; 
+}
+
+//____________________________________________________________________________
+Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) 
+{
+  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
+  Bool_t rv = kFALSE ; 
 
   Int_t relid[4] ; 
-  geom->AbsToRelNumbering(digit->GetId(), relid) ; 
+  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
 
-  if ( relid[1] == 0  )
-    rv = kTRUE; 
+  if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE; 
 
   return rv ; 
 }
 
 //____________________________________________________________________________
-void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, RecPointsList * emcl, RecPointsList * ppsdl)
+Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) 
 {
+  // Tells if (true) or not (false) the digit is in a PHOS-EMC module
+  Bool_t rv = kFALSE ; 
+
+  Int_t relid[4] ; 
+  fGeom->AbsToRelNumbering(digit->GetId(), relid) ; 
+
+  if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE; 
+
+  return rv ; 
+}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, 
+                                       AliPHOSRecPoint::RecPointsList * emcl, 
+                                       AliPHOSRecPoint::RecPointsList * ppsdl, 
+                                       AliPHOSRecPoint::RecPointsList * cpvl)
+{
+  // Steering method to construct the clusters stored in a list of Reconstructed Points
+  // A cluster is defined as a list of neighbour digits
+
+  fNumberOfEmcClusters  = 0 ;
+  fNumberOfPpsdClusters = 0 ;
+  fNumberOfCpvClusters  = 0 ;
+
   // Fill and sort the working digits list
   TObjArray tempodigitslist( dl->GetEntries() ) ;
   FillandSort(dl, &tempodigitslist) ; 
@@ -178,31 +239,47 @@ void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, RecPointsList * e
 
   while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digits
     AliPHOSRecPoint * clu ; 
-   
-    int * clusterdigitslist[dl->GetEntries()] ;   
+
+    AliPHOSDigit ** clusterdigitslist = new AliPHOSDigit*[dl->GetEntries()] ;   
     Int_t index ;
-    if (( ( IsInEmc(digit) ) && ( Calibrate(digit->GetAmp() ) > fEmcClusteringThreshold ) ) || 
-        ( ( !IsInEmc(digit) ) && ( Calibrate(digit->GetAmp() ) > fPpsdClusteringThreshold ) ) ) {
+    if (( IsInEmc (digit) && Calibrate(digit->GetAmp()) > fEmcClusteringThreshold  ) || 
+        ( IsInPpsd(digit) && Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) ||
+        ( IsInCpv (digit) && Calibrate(digit->GetAmp()) > fCpvClusteringThreshold  ) ) {
   
       Int_t iDigitInCluster = 0 ; 
 
       if  ( IsInEmc(digit) ) {   
-        new ((*emcl)[fNumberOfEmcClusters]) AliPHOSEmcRecPoint(fW0, fLocMaxCut) ; // start a new EMC RecPoint
-       clu = (AliPHOSEmcRecPoint *) (*emcl)[fNumberOfEmcClusters] ; 
+       // start a new EMC RecPoint
+       if(fNumberOfEmcClusters >= emcl->GetSize()) emcl->Expand(2*fNumberOfEmcClusters+1) ;
+       (*emcl)[fNumberOfEmcClusters] = new  AliPHOSEmcRecPoint(fW0, fLocMaxCut) ;
+       clu = (AliPHOSEmcRecPoint *) emcl->At(fNumberOfEmcClusters) ; 
        fNumberOfEmcClusters++ ; 
        clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ; 
-
-       clusterdigitslist[iDigitInCluster] = (int* ) digit ;    
+       clusterdigitslist[iDigitInCluster] = digit ;    
        iDigitInCluster++ ; 
        tempodigitslist.Remove(digit) ; 
-      }
 
-      else { 
-       new ((*ppsdl)[fNumberOfPpsdClusters]) AliPHOSPpsdRecPoint() ; // start a new PPSD cluster
-       clu =  (AliPHOSPpsdRecPoint *) ppsdl->At(fNumberOfPpsdClusters)  ;  
-       fNumberOfPpsdClusters++ ; 
+      } else { 
+       
+       // start a new PPSD cluster
+       if(fNumberOfPpsdClusters >= ppsdl->GetSize()) ppsdl->Expand(2*fNumberOfPpsdClusters+1);
+       if(fNumberOfCpvClusters  >= cpvl ->GetSize()) cpvl ->Expand(2*fNumberOfCpvClusters +1);
+       if      (IsInPpsd(digit)) {
+         (*ppsdl)[fNumberOfPpsdClusters] = new AliPHOSPpsdRecPoint() ;
+         clu =  (AliPHOSPpsdRecPoint *) ppsdl->At(fNumberOfPpsdClusters)  ;  
+         fNumberOfPpsdClusters++ ; 
+       }
+       else if (IsInCpv(digit) ) {
+         (*cpvl) [fNumberOfCpvClusters]  = new AliPHOSCpvRecPoint(fW0CPV, fLocMaxCutCPV) ;
+         clu =  (AliPHOSCpvRecPoint  *) cpvl ->At(fNumberOfCpvClusters)  ;  
+         fNumberOfCpvClusters++ ; 
+       }
+       else {
+         cout << "AliPHOSClusterizerv1::MakeClusters: unknown configuration " << fGeom->GetName() << endl;
+         assert(0==1);
+       }
        clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;    
-       clusterdigitslist[iDigitInCluster] = (int* ) digit ;    
+       clusterdigitslist[iDigitInCluster] = digit  ;   
        iDigitInCluster++ ; 
        tempodigitslist.Remove(digit) ; 
         nextdigit.Reset() ;
@@ -210,17 +287,14 @@ void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, RecPointsList * e
        // Here we remove resting EMC digits, which cannot make cluster
 
         if( notremoved ) { 
-         
          while( ( digit = (AliPHOSDigit *)nextdigit() ) ) {
-           
             if( IsInEmc(digit) ) 
              tempodigitslist.Remove(digit) ;
             else 
              break ;
-         
-         } // while digit  
-         
-       } // if notremoved 
+         }
+         notremoved = kFALSE ;
+       }
        
       } // else        
       
@@ -229,16 +303,16 @@ void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, RecPointsList * e
       AliPHOSDigit * digitN ; 
       index = 0 ;
       while (index < iDigitInCluster){ // scan over digits already in cluster 
-       digit = (AliPHOSDigit *) clusterdigitslist[index]  ;      
+       digit =  clusterdigitslist[index]  ;      
        index++ ; 
         while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits 
-         Int_t ineb = AreNeighbours(digit, digitN);   //  call (digit,digitN) in THAT oder !!!!!
+         Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
           switch (ineb ) {
           case 0 :   // not a neighbour
-           break ;      
+           break ;
          case 1 :   // are neighbours 
-           clu->AddDigit( *digitN, Calibrate( digitN->GetAmp() ) ) ;
-           clusterdigitslist[iDigitInCluster] =(int*) digitN ; 
+           clu->AddDigit(*digitN, Calibrate( digitN->GetAmp() ) ) ;
+           clusterdigitslist[iDigitInCluster] = digitN ; 
            iDigitInCluster++ ; 
            tempodigitslist.Remove(digitN) ;
            break ;
@@ -252,20 +326,30 @@ void AliPHOSClusterizerv1::MakeClusters(const DigitsList * dl, RecPointsList * e
        nextdigit.Reset() ; 
        
       } // loop over cluster     
-   }  //below energy theshold  
-  
+    }  //below energy theshold  
+    
+    delete[] clusterdigitslist ; 
+    
   } // while digit
 
   tempodigitslist.Clear() ; 
+
+  ppsdl->Sort() ;
+  Int_t index ;
+  for(index = 0; index < ppsdl->GetEntries(); index++)
+    ((AliPHOSPpsdRecPoint *)ppsdl->At(index))->SetIndexInList(index) ;
 }
 
 //____________________________________________________________________________
 void AliPHOSClusterizerv1::PrintParameters() 
 {
+  // Print the energy thresholds 
+
   cout << "PHOS Clusterizer version 1 :" << endl 
        << "                       EMC  Clustering threshold = " << fEmcClusteringThreshold << endl
        << "                       EMC  Energy threshold     = " << fEmcEnergyThreshold << endl                  
        << "                      PPSD  Clustering threshold = " << fPpsdClusteringThreshold << endl
-       << "                      PPSD  Energy threshold     = " << fPpsdEnergyThreshold << endl ;                
+       << "                      PPSD  Energy threshold     = " << fPpsdEnergyThreshold << endl
+       << "                       CPV  Clustering threshold = " << fCpvClusteringThreshold << endl
+       << "                       CPV  Energy threshold     = " << fCpvEnergyThreshold << endl ;                
 }