Implemented NxM cluster finder

 - One can specify the neighbors to be accepted
   - Default is 3x3
   - Use SetNRowDiff(2), SetNColDiff(2) to set 5x5
 - To achieve 15-20% speed up store calibrated energy in AliEMCALDigit (not streamed!)
 - Tested that the 3x3 before and after this optimization yield exactly the same results

diff --git a/EMCAL/AliEMCALClusterizerNxN.cxx b/EMCAL/AliEMCALClusterizerNxN.cxx
index b3d0b122b29105cb43af7aaafbcdc21562dd4344..63e9ca8c0d68de3126e82246fcacf696bc098197 100644 (file)
--- a/EMCAL/AliEMCALClusterizerNxN.cxx
+++ b/EMCAL/AliEMCALClusterizerNxN.cxx
@@ -67,14 +67,14 @@ ClassImp(AliEMCALClusterizerNxN)
 
 //____________________________________________________________________________
 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
-  : AliEMCALClusterizer()
+: AliEMCALClusterizer(), fNRowDiff(1), fNColDiff(1)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
 }
 
 //____________________________________________________________________________
 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
-  : AliEMCALClusterizer(geometry)
+  : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   // use this contructor to avoid usage of Init() which uses runloader
@@ -84,20 +84,17 @@ AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
 
 //____________________________________________________________________________
 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
-: AliEMCALClusterizer(geometry, calib, caloped)
+: AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1)
 {
-       // ctor, geometry and calibration are initialized elsewhere.
-                               
+  // ctor, geometry and calibration are initialized elsewhere.
 }
 
-
 //____________________________________________________________________________
-  AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
+AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
 {
   // dtor
 }
 
-
 //____________________________________________________________________________
 void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
 {
@@ -210,7 +207,7 @@ Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit *
   coldiff = TMath::Abs(ieta1 - ieta2) ;  
   
   // neighbours +-1 in col and row
-  if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
+  if ( TMath::Abs(coldiff) <= fNColDiff && TMath::Abs(rowdiff) <= fNRowDiff)
     {
       
       AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
@@ -230,20 +227,20 @@ Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit *
 //____________________________________________________________________________
 void AliEMCALClusterizerNxN::MakeClusters()
 {
-  // Steering method to construct the clusters stored in a list of Reconstructed Points
-  // A cluster is defined as a list of neighbour digits
-  // Mar 03, 2007 by PAI
+  // Make clusters
   
-  if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
+  if (fGeom==0) 
+    AliFatal("Did not get geometry from EMCALLoader");
   
   fRecPoints->Clear();
   
   // Set up TObjArray with pointers to digits to work on 
-  //TObjArray *digitsC = new TObjArray();
   TObjArray digitsC;
   TIter nextdigit(fDigitsArr);
   AliEMCALDigit *digit = 0;
   while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
+    Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
+    digit->SetCalibAmp(dEnergyCalibrated);
     digitsC.AddLast(digit);
   }
   
@@ -260,31 +257,27 @@ void AliEMCALClusterizerNxN::MakeClusters()
     Float_t dMaxEnergyDigit = -1;
     AliEMCALDigit *pMaxEnergyDigit = 0;
     nextdigitC.Reset();
-    while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) 
+    while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) ) 
     { // scan over the list of digitsC
-      Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
-      //AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
-      
-      //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold  )
+      Float_t dEnergyCalibrated = digit->GetCalibAmp();
+
       if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
-           {
-             if (dEnergyCalibrated > dMaxEnergyDigit)
+      {
+        if (dEnergyCalibrated > dMaxEnergyDigit) 
         {
           dMaxEnergyDigit = dEnergyCalibrated;
           iMaxEnergyDigit = digit->GetId();
           pMaxEnergyDigit = digit;
         }
-           }
+      }
     }
-    
     if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0) 
     {
       bDone = kTRUE;
       continue;
     }
     
-    AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
-    AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));
+    AliDebug (2, Form("Max digit found: %1.5f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
     
     // keep the candidate digits in a list
     TList clusterDigitList;
@@ -293,34 +286,33 @@ void AliEMCALClusterizerNxN::MakeClusters()
     
     Double_t clusterCandidateEnergy = dMaxEnergyDigit;
     
-    // now loop over the resto of the digits and cluster into NxN cluster 
+    // now loop over the rest of the digits and cluster into NxN cluster 
     // we do not actually cluster yet: we keep them in the list clusterDigitList
     nextdigitC.Reset();
     while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) 
     { // scan over the list of digitsC
       if (digit == pMaxEnergyDigit) continue;
-      Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
+      Float_t dEnergyCalibrated = digit->GetCalibAmp();
       AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
-      //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold  )
       if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0  )
-           {
-             Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
-             if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
-             Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
-             if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!! 
+      {
+        Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
+        if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
+        Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
+        if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!! 
         {      
           clusterDigitList.AddLast(digit) ;
           clusterCandidateEnergy += dEnergyCalibrated;
         }
-           }
+      }
     }// loop over the next digits
     
     // start a cluster here only if a cluster energy is larger than clustering threshold
-    //if (clusterCandidateEnergy > 0.1)
     AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
     if (clusterCandidateEnergy > fECAClusteringThreshold)
     {
-      if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
+      if(fNumberOfECAClusters >= fRecPoints->GetSize()) 
+        fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
       
       AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ; 
       fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
@@ -332,9 +324,8 @@ void AliEMCALClusterizerNxN::MakeClusters()
         AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
         for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
         {
-          
           digit = (AliEMCALDigit*)clusterDigitList.At(idig);
-          Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
+          Float_t dEnergyCalibrated = digit->GetCalibAmp();
           AliDebug(5, Form(" Adding digit %d", digit->GetId()));
           // note: this way the sharing info is lost!
           recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
@@ -353,8 +344,6 @@ void AliEMCALClusterizerNxN::MakeClusters()
     AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));      
   } // while ! done 
   
-  //delete digitsC ; //nope we use an object
-  
   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
 }
 

diff --git a/EMCAL/AliEMCALClusterizerNxN.h b/EMCAL/AliEMCALClusterizerNxN.h
index 007844d497a9a2f4e1bf3ba52894e90f02e61e37..d5b2491a827982d2fd13da8907a3a8a194731eba 100644 (file)
--- a/EMCAL/AliEMCALClusterizerNxN.h
+++ b/EMCAL/AliEMCALClusterizerNxN.h
@@ -37,6 +37,11 @@ public:
 
   virtual const char * Version() const { return "clu-NxN" ; }  
   
+  void SetNRowDiff(Int_t nd) { fNRowDiff = nd; }
+  void SetNColDiff(Int_t nd) { fNColDiff = nd; }
+  Int_t GetNRowDiff() const { return fNRowDiff; } 
+  Int_t GetNColDiff() const { return fNColDiff; } 
+
 protected:
 
   virtual void   MakeClusters();            
@@ -45,8 +50,10 @@ private:
   AliEMCALClusterizerNxN(const AliEMCALClusterizerNxN &); //copy ctor
   AliEMCALClusterizerNxN & operator = (const AliEMCALClusterizerNxN &);
 
-   ClassDef(AliEMCALClusterizerNxN,2)   // Clusterizer implementation version 1
+  Int_t  fNRowDiff;  //how many neighbors to consider along row (phi)
+  Int_t  fNColDiff;  //how many neighbors to consider along col (eta)
 
+  ClassDef(AliEMCALClusterizerNxN,3)   // Clusterizer implementation version 1
 };
 
 #endif // AliEMCALCLUSTERIZERNXN_H

diff --git a/EMCAL/AliEMCALDigit.cxx b/EMCAL/AliEMCALDigit.cxx
index fce7b7b0205b29658ca5ab6c28f4b6278fb8a091..55723d157dab8eb85bd3712ed8abe4590765809d 100644 (file)
--- a/EMCAL/AliEMCALDigit.cxx
+++ b/EMCAL/AliEMCALDigit.cxx
@@ -60,8 +60,8 @@ AliDigitNew(),
   fTimeR(0.),
   fChi2(0.),
   fNDF(0),
-  fDigitType(kUnknown)
-
+  fDigitType(kUnknown),
+  fAmpCalib(-1)
 {
   // default ctor 
 
@@ -102,7 +102,8 @@ AliEMCALDigit::AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t dig
     fTimeR(time),
     fChi2(chi2),
     fNDF(ndf),
-    fDigitType(type)
+    fDigitType(type),
+    fAmpCalib(-1)
 {  
   // ctor with all data 
 
@@ -165,7 +166,8 @@ AliEMCALDigit::AliEMCALDigit(const AliEMCALDigit & digit)
     fTimeR(digit.fTimeR), 
     fChi2(digit.fChi2), 
     fNDF(digit.fNDF),
-    fDigitType(digit.fDigitType)
+    fDigitType(digit.fDigitType),
+    fAmpCalib(digit.fAmpCalib)
 {
   // copy ctor
   // data memebrs of the base class (AliNewDigit)

diff --git a/EMCAL/AliEMCALDigit.h b/EMCAL/AliEMCALDigit.h
index c5771b756626b293d8f11708e9835ba90f8b4b73..e672f8b06d66dbcc41e767d08d8aa9d12fe79a51 100644 (file)
--- a/EMCAL/AliEMCALDigit.h
+++ b/EMCAL/AliEMCALDigit.h
@@ -68,7 +68,7 @@ class AliEMCALDigit : public AliDigitNew {
   void    SetNDF(Int_t ndf)         { fNDF       = ndf  ;}
   void    SetType(Int_t t)          { fDigitType = t    ;}
   void    ShiftPrimary(Int_t shift); // shift to separate different TreeK in merging
-       
+
   //Raw time sample
   //ALTRO
   Int_t   GetNALTROSamplesLG() const {if(fDigitType==kLG)return fNSamples; else return 0;}
@@ -84,6 +84,9 @@ class AliEMCALDigit : public AliDigitNew {
   void SetALTROSamplesLG (const Int_t nSamplesLG, Int_t *samplesLG);
   void SetFALTROSamples  (const Int_t nSamples,   Int_t *samples) { if(fDigitType==kTrigger) SetALTROSamplesLG(nSamples, samples);} 
 
+  void    SetCalibAmp(Float_t amp)  { fAmpCalib = amp; }
+  Double_t GetCalibAmp()   const    { return fAmpCalib; }
+
   void Print(const Option_t* /*opt*/) const;
        
  private: 
@@ -111,9 +114,9 @@ class AliEMCALDigit : public AliDigitNew {
   Int_t   fNDF;         // Fit Number of Degrees of Freedom
        
   Int_t fDigitType;     // This is a trigger digit(0), HG (1) or LG (3)
-       
-  ClassDef(AliEMCALDigit,5)   // Digit in EMCAL 
+  Float_t fAmpCalib;    //!Calibrated energy
 
+  ClassDef(AliEMCALDigit,6)   // Digit in EMCAL 
 } ;
 
 #endif //  ALIEMCALDIGIT_H