New Clusterizer that makes cluster grouping neighbour cells around highest energy...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Aug 2010 07:31:42 +0000 (07:31 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Aug 2010 07:31:42 +0000 (07:31 +0000)
14 files changed:
EMCAL/AliEMCALClusterizer.cxx
EMCAL/AliEMCALClusterizer.h
EMCAL/AliEMCALClusterizerNxN.cxx [new file with mode: 0644]
EMCAL/AliEMCALClusterizerNxN.h [new file with mode: 0644]
EMCAL/AliEMCALClusterizerv1.cxx
EMCAL/AliEMCALClusterizerv1.h
EMCAL/AliEMCALRecParam.cxx
EMCAL/AliEMCALRecParam.h
EMCAL/AliEMCALReconstructor.cxx
EMCAL/AliEMCALReconstructor.h
EMCAL/EMCALrecLinkDef.h
EMCAL/libEMCALrec.pkg
EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C
EMCAL/macros/RecParamDB/PrintEMCALRecParam.C

index bf4bda5..90b0cb4 100644 (file)
 //  Base class for the clusterization algorithm (pure abstract)
 //*--
 //*-- Author: Yves Schutz  SUBATECH 
-// Modif: 
-//  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
-//                           of new  IO (à la PHOS)
+// 
+//   Clusterization mother class. Contains common methods/data members of different 
+//   clusterizers. GCB 2010
 //////////////////////////////////////////////////////////////////////////////
 
 // --- ROOT system ---
 #include "TClonesArray.h"
 #include "TTree.h"
+#include <TFile.h> 
+class TFolder;
+#include <TMath.h> 
+#include <TMinuit.h>
+#include <TTree.h> 
+class TSystem; 
+#include <TBenchmark.h>
+#include <TBrowser.h>
+#include <TROOT.h>
 
 // --- Standard library ---
-
+#include <cassert>
 
 // --- AliRoot header files ---
 #include "AliEMCALClusterizer.h"
+#include "AliEMCALReconstructor.h"
+#include "AliRunLoader.h"
+#include "AliRun.h"
 #include "AliLog.h"
+#include "AliEMCAL.h"
+#include "AliEMCALRecPoint.h"
+#include "AliEMCALRecParam.h"
+#include "AliEMCALGeometry.h"
+#include "AliEMCALRecParam.h"
+#include "AliCDBManager.h"
+#include "AliCaloCalibPedestal.h"
+#include "AliEMCALCalibData.h"
+class AliCDBStorage;
+#include "AliCDBEntry.h"
 
 ClassImp(AliEMCALClusterizer)
 
@@ -41,12 +63,75 @@ ClassImp(AliEMCALClusterizer)
 AliEMCALClusterizer::AliEMCALClusterizer():
   fDigitsArr(NULL),
   fTreeR(NULL),
-  fRecPoints(NULL)
+  fRecPoints(NULL),
+  fGeom(NULL),
+  fCalibData(NULL), 
+  fCaloPed(NULL),
+  fADCchannelECA(0.),fADCpedestalECA(0.),
+  fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
+  fDefaultInit(kFALSE),fToUnfold(kFALSE),
+  fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
+  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
 {
   // ctor
+  
+  Init();
+  
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry): 
+  fDigitsArr(NULL),
+  fTreeR(NULL),
+  fRecPoints(NULL),
+  fGeom(geometry),
+  fCalibData(NULL), 
+  fCaloPed(NULL),
+  fADCchannelECA(0.),fADCpedestalECA(0.),
+  fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
+  fDefaultInit(kFALSE),fToUnfold(kFALSE),
+  fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
+  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
+{
+  // 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
+  // change needed by HLT - MP
+  
+  // Note for the future: the use on runloader should be avoided or optional at least
+  // another way is to make Init virtual and protected at least such that the deriving classes can overload
+  // Init() ;
+  //
+  
+  if (!fGeom)
+  {
+    AliFatal("Geometry not initialized.");
+  }
+  
 }
 
 //____________________________________________________________________________
+AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped): 
+  fDigitsArr(NULL),
+  fTreeR(NULL),
+  fRecPoints(NULL),
+  fGeom(geometry),
+  fCalibData(calib),
+  fCaloPed(caloped),
+  fADCchannelECA(0.),fADCpedestalECA(0.),
+  fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
+  fDefaultInit(kFALSE),fToUnfold(kFALSE),
+  fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
+  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
+{
+       // ctor, geometry and calibration are initialized elsewhere.
+       
+       if (!fGeom)
+               AliFatal("Geometry not initialized.");
+  
+}
+
+
+//____________________________________________________________________________
 AliEMCALClusterizer::~AliEMCALClusterizer()
 {
   // dtor
@@ -61,6 +146,240 @@ AliEMCALClusterizer::~AliEMCALClusterizer()
 }
 
 //____________________________________________________________________________
+Float_t  AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId) 
+{
+  
+  // Convert digitized amplitude into energy.
+  // Calibration parameters are taken from calibration data base for raw data,
+  // or from digitizer parameters for simulated data.
+  if(fCalibData){
+    
+    if (fGeom==0)
+      AliFatal("Did not get geometry from EMCALLoader") ;
+    
+    Int_t iSupMod = -1;
+    Int_t nModule = -1;
+    Int_t nIphi   = -1;
+    Int_t nIeta   = -1;
+    Int_t iphi    = -1;
+    Int_t ieta    = -1;
+    
+    Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
+    if(!bCell) {
+      fGeom->PrintGeometry();
+      Error("Calibrate()"," Wrong cell id number : %i", absId);
+      assert(0);
+    }
+    
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
+         
+    // Check if channel is bad (dead or hot), in this case return 0.   
+    // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
+    // for the moment keep it here but remember to do the selection at the sdigitizer level 
+    // and remove it from here
+    Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
+    if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
+                 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
+                 return 0;
+    }
+    //Check if time is too large or too small, indication of a noisy channel, remove in this case
+    if(time > fTimeMax || time < fTimeMin) return 0;
+         
+    fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
+    fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
+    
+    return -fADCpedestalECA + amp * fADCchannelECA ;        
+    
+  }
+  else //Return energy with default parameters if calibration is not available
+    return -fADCpedestalECA + amp * fADCchannelECA ; 
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizer::GetCalibrationParameters() 
+{
+  // Set calibration parameters:
+  // if calibration database exists, they are read from database,
+  // otherwise, they are taken from digitizer.
+  //
+  // It is a user responsilibity to open CDB before reconstruction, 
+  // for example: 
+  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+  
+  //Check if calibration is stored in data base
+  
+  if(!fCalibData)
+  {
+    AliCDBEntry *entry = (AliCDBEntry*) 
+    AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+    if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
+  }
+  
+  if(!fCalibData)
+    AliFatal("Calibration parameters not found in CDB!");
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizer::GetCaloCalibPedestal() 
+{
+       // Set calibration parameters:
+       // if calibration database exists, they are read from database,
+       // otherwise, they are taken from digitizer.
+       //
+       // It is a user responsilibity to open CDB before reconstruction, 
+       // for example: 
+       // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
+       
+       //Check if calibration is stored in data base
+       
+       if(!fCaloPed)
+  {
+               AliCDBEntry *entry = (AliCDBEntry*) 
+               AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+               if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
+  }
+       
+       if(!fCaloPed)
+               AliFatal("Pedestal info not found in CDB!");
+       
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizer::Init()
+{
+  // Make all memory allocations which can not be done in default constructor.
+  // Attach the Clusterizer task to the list of EMCAL tasks
+  
+  AliRunLoader *rl = AliRunLoader::Instance();
+  if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
+    fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  else 
+    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
+  
+  AliDebug(1,Form("geom %p",fGeom));
+  
+  if(!gMinuit) 
+    gMinuit = new TMinuit(100) ;
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizer::InitParameters()
+{ 
+  // Initializes the parameters for the Clusterizer
+  fNumberOfECAClusters = 0 ;
+  fCalibData           = 0 ;
+  fCaloPed             = 0 ;
+       
+  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
+  if(!recParam) {
+    AliFatal("Reconstruction parameters for EMCAL not set!");
+  } else {
+    fECAClusteringThreshold = recParam->GetClusteringThreshold();
+    fECAW0                  = recParam->GetW0();
+    fMinECut                = recParam->GetMinECut();    
+    fToUnfold               = recParam->GetUnfold();
+    if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!"); 
+    fECALocMaxCut           = recParam->GetLocMaxCut();
+    fTimeCut                = recParam->GetTimeCut();
+    fTimeMin                = recParam->GetTimeMin();
+    fTimeMax                = recParam->GetTimeMax();
+    
+    AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
+                    fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
+  }
+  
+}
+
+
+//____________________________________________________________________________
+void AliEMCALClusterizer::Print(Option_t * /*option*/)const
+{
+  // Print clusterizer parameters
+  
+  TString message("\n") ; 
+  
+  if( strcmp(GetName(), "") !=0 ){
+    
+    // Print parameters
+    
+    TString taskName(Version()) ;
+    
+    printf("--------------- "); 
+    printf("%s",taskName.Data()) ; 
+    printf(" "); 
+    printf("Clusterizing digits: "); 
+    printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
+    printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
+    if(fToUnfold)
+      printf("\nUnfolding on\n");
+    else
+      printf("\nUnfolding off\n");
+    
+    printf("------------------------------------------------------------------"); 
+  }
+  else
+    printf("AliEMCALClusterizer not initialized ") ;
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
+{
+  // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
+  if(strstr(option,"deb")) {
+    printf("PrintRecPoints: Clusterization result:") ; 
+    
+    printf("           Found %d ECA Rec Points\n ", 
+           fRecPoints->GetEntriesFast()) ; 
+  }
+  
+  if(strstr(option,"all")) {
+    if(strstr(option,"deb")) {
+      printf("\n-----------------------------------------------------------------------\n") ;
+      printf("Clusters in ECAL section\n") ;
+      printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
+    }
+    Int_t index; 
+    for (index =  0 ; index < fRecPoints->GetEntries() ; index++) {
+      AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ; 
+      TVector3  globalpos;  
+      //rp->GetGlobalPosition(globalpos);
+      TVector3  localpos;  
+      rp->GetLocalPosition(localpos);
+      Float_t lambda[2]; 
+      rp->GetElipsAxis(lambda);
+      
+      Int_t nprimaries=0;
+      Int_t * primaries = rp->GetPrimaries(nprimaries);
+      
+      if(strstr(option,"deb")) 
+        printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
+               rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
+               globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
+               rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
+      if(strstr(option,"deb")){ 
+        for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
+          printf("%d ", primaries[iprimary] ) ; 
+        }
+      }
+    }
+    
+    if(strstr(option,"deb"))
+      printf("\n-----------------------------------------------------------------------\n");
+  }
+}
+
+//___________________________________________________________________
+void  AliEMCALClusterizer::PrintRecoInfo()
+{
+  //print reco version
+  printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
+  
+}
+
+//____________________________________________________________________________
 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
 {
   // Read the digits from the input tree
@@ -87,3 +406,6 @@ void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
   Int_t bufsize = 32000;
   fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
 }
+
+
+
index 388213f..33584e6 100644 (file)
@@ -19,6 +19,9 @@ class TTree;
 // --- Standard library ---
 
 // --- AliRoot header files ---
+class AliEMCALGeometry ;
+class AliEMCALCalibData ;
+class AliCaloCalibPedestal ;
 
 class AliEMCALClusterizer : public TObject {
 
@@ -26,22 +29,47 @@ public:
 
   AliEMCALClusterizer() ;        // default ctor
   virtual ~AliEMCALClusterizer() ; // dtorEM
+  AliEMCALClusterizer(AliEMCALGeometry* geometry);
+  AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * pedestal);
 
   virtual void    Digits2Clusters(Option_t *option) = 0;
 
-  virtual Float_t GetTimeCut() const = 0;
-
-  virtual void SetECAClusteringThreshold(Float_t) = 0;
-  virtual void SetECALocalMaxCut(Float_t)         = 0;
-  virtual void SetECALogWeight(Float_t)           = 0;
-  virtual void SetTimeCut(Float_t)                = 0;
-  virtual void SetUnfolding(Bool_t)               = 0;
-
-  virtual const char * Version() const {Warning("Version", "Not Defined") ; return 0 ; } 
+  virtual Float_t Calibrate(const Float_t amp, const Float_t time, const Int_t cellId) ;  // Tranforms Amp to energy 
+  virtual void    Init() ;
+  virtual void    InitParameters() ; //{ AliInfo("Overload this method."); }
+
+  //Get/Set reconstruction parameters
+  virtual void  GetCalibrationParameters(void) ;
+  virtual void  GetCaloCalibPedestal(void) ;
+  virtual void  SetCalibrationParameters(AliEMCALCalibData * calib)   { fCalibData = calib ; }
+  virtual void  SetCaloCalibPedestal(AliCaloCalibPedestal  * caloped) { fCaloPed   = caloped ; }
+  
+  virtual Float_t GetTimeMin()           const { return fTimeMin ; }
+  virtual Float_t GetTimeMax()           const { return fTimeMax ; }
+  virtual Float_t GetTimeCut()           const { return fTimeCut ; }
+  virtual void    GetNumberOfClustersFound(int numb )const { numb = fNumberOfECAClusters ;} 
+  virtual Float_t GetECAClusteringThreshold()        const { return fECAClusteringThreshold;}  
+  virtual Float_t GetECALocalMaxCut()                const { return fECALocMaxCut;} 
+  virtual Float_t GetECALogWeight()                  const { return fECAW0;}
+  virtual Float_t GetMinECut()                       const { return fMinECut;}
+
+  virtual void SetTimeMin(Float_t t)                                           { fTimeMin = t ;}
+  virtual void SetTimeMax(Float_t t)                                           { fTimeMax = t ;}
+  virtual void SetTimeCut(Float_t t)                                           { fTimeCut = t ;}
+  virtual void SetECAClusteringThreshold(Float_t th)  { fECAClusteringThreshold = th ; }
+  virtual void SetMinECut(Float_t mine)               { fMinECut = mine; }
+  virtual void SetECALocalMaxCut(Float_t cut)         { fECALocMaxCut = cut ; }
+  virtual void SetECALogWeight(Float_t w)             { fECAW0 = w ; }
+  virtual void SetUnfolding(Bool_t toUnfold = kTRUE ) {fToUnfold = toUnfold ;}  
 
   virtual void SetInput(TTree *digitsTree);
   virtual void SetOutput(TTree *clustersTree);
-  virtual void InitParameters() { AliInfo("Overload this method."); }
+  
+  virtual void Print(Option_t * option)const ;
+  virtual void PrintRecPoints(Option_t * option);
+  virtual void PrintRecoInfo();                        //*MENU*
+  
+  virtual const char * Version() const {Warning("Version", "Not Defined") ; return 0 ; } 
 
 protected:
 
@@ -50,12 +78,34 @@ protected:
   TClonesArray *fDigitsArr; // Array with EMCAL digits
   TTree *fTreeR;            // Tree with output clusters
   TObjArray    *fRecPoints; // Array with EMCAL clusters
-
+  
+  AliEMCALGeometry     * fGeom;       //! pointer to geometry for utilities
+  AliEMCALCalibData    * fCalibData ; //! Calibration database if aval
+  AliCaloCalibPedestal * fCaloPed   ; //! Tower status map if aval
+  
+  Float_t fADCchannelECA ;           // width of one ADC channel for EC section (GeV)
+  Float_t fADCpedestalECA ;          // pedestal of ADC for EC section (GeV) 
+
+  Float_t fTimeMin ;                 // Minimum time of physical signal in a cell/digit
+  Float_t fTimeMax ;                 // Maximum time of physical signal in a cell/digit
+  Float_t fTimeCut ;                 // Maximum time difference between the digits inside EMC cluster
+
+  Bool_t  fDefaultInit;              //! Says if the task was created by defaut ctor (only parameters are initialized)
+  Bool_t  fToUnfold ;                // To perform unfolding 
+  Int_t   fNumberOfECAClusters ;     // number of clusters found in EC section
+  
+  Float_t fECAClusteringThreshold ;  // minimum energy to seed a EC digit in a cluster
+  Float_t fECALocMaxCut ;            // minimum energy difference to distinguish local maxima in a cluster
+  Float_t fECAW0 ;                   // logarithmic weight for the cluster center of gravity calculation
+  Float_t fMinECut;                  // Minimum energy for a digit to be a member of a cluster
+  
 private:
   AliEMCALClusterizer(const AliEMCALClusterizer &); //copy ctor
   AliEMCALClusterizer & operator = (const AliEMCALClusterizer &);
-
-  ClassDef(AliEMCALClusterizer,1)  // Clusterization algorithm class 
+  
+  
+  
+  ClassDef(AliEMCALClusterizer,2)  // Clusterization algorithm class 
 } ;
 
 #endif // AliEMCALCLUSTERIZER_H
diff --git a/EMCAL/AliEMCALClusterizerNxN.cxx b/EMCAL/AliEMCALClusterizerNxN.cxx
new file mode 100644 (file)
index 0000000..98eebec
--- /dev/null
@@ -0,0 +1,432 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
+// Algorithm:
+// 1. peek the most energetic cell
+// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
+// 3. remove the cells contributing to the cluster
+// 4. start from 1 for the remaining clusters
+// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
+// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy 
+// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
+// Use Case:
+//  root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")  
+//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
+//               //reads gAlice from header file "..."                      
+//  root [1] cl->ExecuteTask()  
+//               //finds RecPoints in all events stored in galice.root
+//  root [2] cl->SetDigitsBranch("digits2") 
+//               //sets another title for Digitis (input) branch
+//  root [3] cl->SetRecPointsBranch("recp2")  
+//               //sets another title four output branches
+//  root [4] cl->SetTowerLocalMaxCut(0.03)  
+//               //set clusterization parameters
+//  root [5] cl->ExecuteTask("deb all time")  
+//               //once more finds RecPoints options are 
+//               // deb - print number of found rec points
+//               // deb all - print number of found RecPoints and some their characteristics 
+//               // time - print benchmarking results
+
+// --- ROOT system ---
+#include <TMath.h> 
+#include <TMinuit.h>
+#include <TTree.h> 
+#include <TBenchmark.h>
+#include <TBrowser.h>
+#include <TROOT.h>
+#include <TClonesArray.h>
+
+// --- Standard library ---
+#include <cassert>
+
+// --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliEMCALClusterizerNxN.h"
+#include "AliEMCALRecPoint.h"
+#include "AliEMCALDigit.h"
+#include "AliEMCALGeometry.h"
+#include "AliCaloCalibPedestal.h"
+#include "AliEMCALCalibData.h"
+#include "AliESDCaloCluster.h"
+
+ClassImp(AliEMCALClusterizerNxN)
+
+Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;
+
+//____________________________________________________________________________
+AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
+  : AliEMCALClusterizer()
+{
+  // ctor with the indication of the file where header Tree and digits Tree are stored
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
+  : AliEMCALClusterizer(geometry)
+{
+  // 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
+  // change needed by HLT - MP
+
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
+: AliEMCALClusterizer(geometry, calib, caloped)
+{
+       // ctor, geometry and calibration are initialized elsewhere.
+                               
+}
+
+
+//____________________________________________________________________________
+  AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
+{
+  // dtor
+}
+
+
+//____________________________________________________________________________
+void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
+{
+  // Steering method to perform clusterization for the current event 
+  // in AliEMCALLoader
+
+  if(strstr(option,"tim"))
+    gBenchmark->Start("EMCALClusterizer"); 
+  
+  if(strstr(option,"print"))
+    Print("") ; 
+  //Get calibration parameters from file or digitizer default values.
+  GetCalibrationParameters() ;
+
+  //Get dead channel map from file or digitizer default values.
+  GetCaloCalibPedestal() ;
+       
+  fNumberOfECAClusters = 0;
+
+  MakeClusters() ;  //only the real clusters
+
+  if(fToUnfold)
+    MakeUnfolding() ;
+
+  Int_t index ;
+
+  //Evaluate position, dispersion and other RecPoint properties for EC section                      
+  for(index = 0; index < fRecPoints->GetEntries(); index++) 
+    {
+      dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
+      AliDebug(5, Form("MAX INDEX %d ", dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->GetMaximalEnergyIndex()));
+      //For each rec.point set the distance to the nearest bad crystal
+      dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
+    }
+  
+  fRecPoints->Sort() ;
+
+  for(index = 0; index < fRecPoints->GetEntries(); index++) 
+    {
+      (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
+      (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
+    }
+  
+  fTreeR->Fill();
+  
+  if(strstr(option,"deb") || strstr(option,"all"))  
+    PrintRecPoints(option) ;
+
+  AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
+
+  fRecPoints->Delete();
+
+  if(strstr(option,"tim")){
+    gBenchmark->Stop("EMCALClusterizer");
+    printf("Exec took %f seconds for Clusterizing", 
+          gBenchmark->GetCpuTime("EMCALClusterizer"));
+  }    
+}
+
+//____________________________________________________________________________
+Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
+{
+       // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
+       //                                       = 1 are neighbour
+       //                                       = 2 is in different SM; continue searching 
+       // In case it is in different SM, but same phi rack, check if neigbours at eta=0
+       // neighbours are defined as digits having at least a common side 
+       // 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  
+       
+       static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
+       static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
+       static Int_t rowdiff=0, coldiff=0;
+
+       shared = kFALSE;
+       
+       fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
+       fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
+       fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
+       fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
+       
+       //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
+       if(nSupMod1 != nSupMod2 ) 
+         {
+           //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
+           Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
+           Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
+           
+           if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
+           
+           // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
+           // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
+           if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
+           else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
+           
+           shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
+           
+       }//Different SM, same phi
+       
+       rowdiff = TMath::Abs(iphi1 - iphi2);  
+       coldiff = TMath::Abs(ieta1 - ieta2) ;  
+
+       // neighbours +-1 in col and row
+       if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
+         {
+           
+           AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
+                            d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
+           
+           return 1;
+         }//Neighbours
+       else 
+         {
+           AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
+                            d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
+           shared = kFALSE;
+           return 2 ; 
+         }//Not neighbours
+}
+
+//____________________________________________________________________________
+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
+
+  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())) ) {
+    digitsC.AddLast(digit);
+  }
+
+  TIter nextdigitC(&digitsC);
+
+  AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f)\n",
+                 fDigitsArr->GetEntries(),fMinECut));
+
+  Bool_t bDone = kFALSE;
+  while ( bDone != kTRUE )
+    {
+      //first sort the digits:
+      Int_t iMaxEnergyDigit = -1;
+      Float_t dMaxEnergyDigit = -1;
+      AliEMCALDigit *pMaxEnergyDigit = 0;
+      nextdigitC.Reset();
+      while ( (digit = dynamic_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  )
+         if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
+           {
+             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));
+
+      // keep the candidate digits in a list
+      TList clusterDigitList;
+      clusterDigitList.SetOwner(kFALSE);
+      clusterDigitList.AddLast(pMaxEnergyDigit);        
+
+      Double_t clusterCandidateEnergy = dMaxEnergyDigit;
+
+      // now loop over the resto 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());
+         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 !!!!! 
+               {      
+                 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) ;
+      
+         AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ; 
+         fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
+         recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ; 
+         fNumberOfECAClusters++ ;       
+         recPoint->SetClusterType(AliESDCaloCluster::kEMCALClusterv1);
+
+         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());
+             AliDebug(5, Form(" Adding digit %d", digit->GetId()));
+             // note: this way the sharing info is lost!
+             recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
+             digitsC.Remove(digit);              
+           }
+       }
+      else
+       {
+         // we do not want to start clustering in the same spot!
+         // but in this case we may NOT reuse this seed for another cluster!
+         // need a better bookeeping?
+         digitsC.Remove(pMaxEnergyDigit);
+       }
+
+      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())); 
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerNxN::MakeUnfolding()
+{
+  // Unfolds clusters using the shape of an ElectroMagnetic shower
+  // Performs unfolding of all clusters
+               
+  if(fNumberOfECAClusters > 0){
+    if (fGeom==0)
+      AliFatal("Did not get geometry from EMCALLoader") ;
+    Int_t nModulesToUnfold = fGeom->GetNCells();
+
+    Int_t numberofNotUnfolded = fNumberOfECAClusters ;
+    Int_t index ;
+    for(index = 0 ; index < numberofNotUnfolded ; index++){
+
+      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
+
+      TVector3 gpos;
+      Int_t absId;
+      recPoint->GetGlobalPosition(gpos);
+      fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
+      if(absId > nModulesToUnfold)
+        break ;
+
+      Int_t nMultipl = recPoint->GetMultiplicity() ;
+      AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
+      Float_t * maxAtEnergy = new Float_t[nMultipl] ;
+      Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
+
+      if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
+        //UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
+        fRecPoints->Remove(recPoint);
+        fRecPoints->Compress() ;
+        index-- ;
+        fNumberOfECAClusters-- ;
+        numberofNotUnfolded-- ;
+      }
+      else{
+        recPoint->SetNExMax(1) ; //Only one local maximum
+      }
+
+      delete[] maxAt ;
+      delete[] maxAtEnergy ;
+    }
+  }
+  // End of Unfolding of clusters
+}
+
+//____________________________________________________________________________
+Double_t  AliEMCALClusterizerNxN::ShowerShape(Double_t x, Double_t y)
+{ 
+  // Shape of the shower
+  // If you change this function, change also the gradient evaluation in ChiSquare()
+  //
+  Double_t r = sqrt(x*x+y*y);
+  Double_t r133  = TMath::Power(r, 1.33) ;
+  Double_t r669  = TMath::Power(r, 6.69) ;
+  Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
+  return shape ;
+}
+
+//____________________________________________________________________________
+void  AliEMCALClusterizerNxN::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, 
+                                           Int_t /*nMax*/, 
+                                           AliEMCALDigit ** /*maxAt*/, 
+                                           Float_t * /*maxAtEnergy*/)
+{
+  //
+  // Performs the unfolding of a cluster with nMax overlapping showers 
+  //
+  AliWarning("Not implemented. To be.");
+}
+
+//___________________________________________________________________
+void   AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
+{
+  //
+  // input is calibrated - the case when we run already on ESD
+  //
+  AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;
+}
diff --git a/EMCAL/AliEMCALClusterizerNxN.h b/EMCAL/AliEMCALClusterizerNxN.h
new file mode 100644 (file)
index 0000000..1d4c893
--- /dev/null
@@ -0,0 +1,68 @@
+#ifndef ALIEMCALCLUSTERIZERNXN_H
+#define ALIEMCALCLUSTERIZERNXN_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliEMCALClusterizerNxN.h 41181 2010-05-12 13:58:06Z gconesab $ */
+
+//_________________________________________________________________________
+// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
+// Algorithm:
+// 1. peek the most energetic cell
+// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
+// 3. remove the cells contributing to the cluster
+// 4. start from 1 for the remaining clusters
+// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
+// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy 
+// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
+
+#include "AliEMCALClusterizer.h"
+class AliEMCALRecPoint ; 
+class AliEMCALDigit ;
+
+class AliEMCALClusterizerNxN : public AliEMCALClusterizer {
+  
+public:
+  
+  AliEMCALClusterizerNxN() ;         
+  AliEMCALClusterizerNxN(AliEMCALGeometry* geometry);
+  AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * pedestal);
+       
+  virtual ~AliEMCALClusterizerNxN()  ;
+
+  virtual Int_t   AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared)const ; 
+                               // Checks if digits are in neighbour cells 
+
+   virtual void   Digits2Clusters(Option_t *option);                // Does the job
+
+
+  static Double_t ShowerShape(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; 
+                                            //class member function (not object member function)
+  static void     SetInputCalibrated(Bool_t val);
+
+  virtual const char * Version() const { return "clu-NxN" ; }  
+  
+protected:
+
+  virtual void   MakeClusters();            
+
+  static Bool_t fgkIsInputCalibrated; // to enable reclusterization from ESD cells
+
+private:
+  AliEMCALClusterizerNxN(const AliEMCALClusterizerNxN &); //copy ctor
+  AliEMCALClusterizerNxN & operator = (const AliEMCALClusterizerNxN &);
+
+
+  virtual void   MakeUnfolding();
+  void           UnfoldCluster(AliEMCALRecPoint * iniEmc, Int_t Nmax, 
+                               AliEMCALDigit ** maxAt,
+                               Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package
+
+private:
+
+   ClassDef(AliEMCALClusterizerNxN,1)   // Clusterizer implementation version 1
+
+};
+
+#endif // AliEMCALCLUSTERIZERNXN_H
index 2a321eb..c834a99 100644 (file)
 /* $Id$ */
 
 //-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
-//  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
-//                           of new  IO (à la PHOS)
-//  Mar 2007, Aleksei Pavlinov - new algoritmh of pseudo clusters
+//--         Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
 //////////////////////////////////////////////////////////////////////////////
 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and 
 //  unfolds the clusters having several local maxima.  
 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all 
 //  parameters including input digits branch title, thresholds etc.)
-//  This TTask is normally called from Reconstructioner, but can as well be used in 
-//  standalone mode.
-// Use Case:
-//  root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")  
-//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-//               //reads gAlice from header file "..."                      
-//  root [1] cl->ExecuteTask()  
-//               //finds RecPoints in all events stored in galice.root
-//  root [2] cl->SetDigitsBranch("digits2") 
-//               //sets another title for Digitis (input) branch
-//  root [3] cl->SetRecPointsBranch("recp2")  
-//               //sets another title four output branches
-//  root [4] cl->SetTowerLocalMaxCut(0.03)  
-//               //set clusterization parameters
-//  root [5] cl->ExecuteTask("deb all time")  
-//               //once more finds RecPoints options are 
-//               // deb - print number of found rec points
-//               // deb all - print number of found RecPoints and some their characteristics 
-//               // time - print benchmarking results
+//
 
 // --- ROOT system ---
-#include <cassert>
 
-class TROOT;
-#include <TH1.h>
 #include <TFile.h> 
-class TFolder;
 #include <TMath.h> 
 #include <TMinuit.h>
 #include <TTree.h> 
-class TSystem; 
 #include <TBenchmark.h>
 #include <TBrowser.h>
 #include <TROOT.h>
+#include <TList.h>
+#include <TClonesArray.h>
 
 // --- Standard library ---
-
+#include <cassert>
 
 // --- AliRoot header files ---
-#include "AliRunLoader.h"
-#include "AliRun.h"
-#include "AliESD.h"
+#include "AliLog.h"
 #include "AliEMCALClusterizerv1.h"
 #include "AliEMCALRecPoint.h"
 #include "AliEMCALDigit.h"
-#include "AliEMCALDigitizer.h"
-#include "AliEMCAL.h"
 #include "AliEMCALGeometry.h"
-#include "AliEMCALRecParam.h"
-#include "AliEMCALReconstructor.h"
-#include "AliCDBManager.h"
 #include "AliCaloCalibPedestal.h"
 #include "AliEMCALCalibData.h"
-class AliCDBStorage;
-#include "AliCDBEntry.h"
+#include "AliESDCaloCluster.h"
 
 ClassImp(AliEMCALClusterizerv1)
 
 //____________________________________________________________________________
-AliEMCALClusterizerv1::AliEMCALClusterizerv1()
-  : AliEMCALClusterizer(),
-    fGeom(0),
-    fDefaultInit(kFALSE),
-    fToUnfold(kFALSE),
-    fNumberOfECAClusters(0),fCalibData(0),fCaloPed(0),
-    fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
-    fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.)
+AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   
@@ -100,45 +62,20 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1()
 
 //____________________________________________________________________________
 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
-  : AliEMCALClusterizer(),
-    fGeom(geometry),
-    fDefaultInit(kFALSE),
-    fToUnfold(kFALSE),
-    fNumberOfECAClusters(0),fCalibData(0), fCaloPed(0),
-    fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
-    fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.)
+  : AliEMCALClusterizer(geometry)
 {
   // 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
   // change needed by HLT - MP
 
-  // Note for the future: the use on runloader should be avoided or optional at least
-  // another way is to make Init virtual and protected at least such that the deriving classes can overload
-  // Init() ;
-  //
-
-  if (!fGeom)
-    {
-      AliFatal("Geometry not initialized.");
-    }
-
 }
 
 //____________________________________________________________________________
 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
-: AliEMCALClusterizer(),
-fGeom(geometry),
-fDefaultInit(kFALSE),
-fToUnfold(kFALSE),
-fNumberOfECAClusters(0),fCalibData(calib), fCaloPed(caloped),
-fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
-fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.)
+: AliEMCALClusterizer(geometry, calib, caloped)
 {
        // ctor, geometry and calibration are initialized elsewhere.
-       
-       if (!fGeom)
-               AliFatal("Geometry not initialized.");
-                       
+                               
 }
 
 
@@ -149,58 +86,6 @@ fECAW0(0.),fTimeCut(1.),fTimeMin(-1.),fTimeMax(1.),fMinECut(0.)
 }
 
 //____________________________________________________________________________
-Float_t  AliEMCALClusterizerv1::Calibrate(const Float_t amp, const Float_t time, const Int_t absId) 
-{
-  // Convert digitized amplitude into energy.
-  // Calibration parameters are taken from calibration data base for raw data,
-  // or from digitizer parameters for simulated data.
-
-  if(fCalibData){
-    
-    if (fGeom==0)
-      AliFatal("Did not get geometry from EMCALLoader") ;
-    
-    Int_t iSupMod = -1;
-    Int_t nModule = -1;
-    Int_t nIphi   = -1;
-    Int_t nIeta   = -1;
-    Int_t iphi    = -1;
-    Int_t ieta    = -1;
-    
-    Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
-    if(!bCell) {
-      fGeom->PrintGeometry();
-      Error("Calibrate()"," Wrong cell id number : %i", absId);
-      assert(0);
-    }
-
-    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
-         
-       // Check if channel is bad (dead or hot), in this case return 0.        
-       // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
-       // for the moment keep it here but remember to do the selection at the sdigitizer level 
-       // and remove it from here
-       Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
-       if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
-                 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
-                 return 0;
-       }
-       //Check if time is too large or too small, indication of a noisy channel, remove in this case
-       if(time > fTimeMax || time < fTimeMin) return 0;
-         
-    fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
-    fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
-  
-   return -fADCpedestalECA + amp * fADCchannelECA ;        
-  }
-  else //Return energy with default parameters if calibration is not available
-    return -fADCpedestalECA + amp * fADCchannelECA ; 
-  
-}
-
-//____________________________________________________________________________
 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
 {
   // Steering method to perform clusterization for the current event 
@@ -339,8 +224,8 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit
     return kFALSE ;
   }
   for(index = 0; index < nPar; index++){
-    Double_t err ;
-    Double_t val ;
+    Double_t err = 0. ;
+    Double_t val = 0. ;
     gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
     fitparameters[index] = val ;
   }
@@ -351,105 +236,6 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit
 }
 
 //____________________________________________________________________________
-void AliEMCALClusterizerv1::GetCalibrationParameters() 
-{
-  // Set calibration parameters:
-  // if calibration database exists, they are read from database,
-  // otherwise, they are taken from digitizer.
-  //
-  // It is a user responsilibity to open CDB before reconstruction, 
-  // for example: 
-  // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
-
-  //Check if calibration is stored in data base
-
-  if(!fCalibData)
-    {
-      AliCDBEntry *entry = (AliCDBEntry*) 
-       AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
-      if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
-    }
-  
-  if(!fCalibData)
-    AliFatal("Calibration parameters not found in CDB!");
-}
-
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::GetCaloCalibPedestal() 
-{
-       // Set calibration parameters:
-       // if calibration database exists, they are read from database,
-       // otherwise, they are taken from digitizer.
-       //
-       // It is a user responsilibity to open CDB before reconstruction, 
-       // for example: 
-       // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
-       
-       //Check if calibration is stored in data base
-       
-       if(!fCaloPed)
-    {
-               AliCDBEntry *entry = (AliCDBEntry*) 
-               AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
-               if (entry) fCaloPed =  (AliCaloCalibPedestal*) entry->GetObject();
-    }
-       
-       if(!fCaloPed)
-               AliFatal("Pedestal info not found in CDB!");
-       
-}
-
-
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::Init()
-{
-  // Make all memory allocations which can not be done in default constructor.
-  // Attach the Clusterizer task to the list of EMCAL tasks
-  
-  AliRunLoader *rl = AliRunLoader::Instance();
-  if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
-    fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
-  else 
-    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
-
-  AliDebug(1,Form("geom %p",fGeom));
-
-  if(!gMinuit) 
-    gMinuit = new TMinuit(100) ;
-
-}
-
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::InitParameters()
-{ 
-  // Initializes the parameters for the Clusterizer
-  fNumberOfECAClusters = 0;
-
-  fCalibData               = 0 ;
-  fCaloPed                 = 0 ;
-       
-  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
-  if(!recParam) {
-    AliFatal("Reconstruction parameters for EMCAL not set!");
-  } else {
-    fECAClusteringThreshold = recParam->GetClusteringThreshold();
-    fECAW0                  = recParam->GetW0();
-    fMinECut                = recParam->GetMinECut();    
-    fToUnfold               = recParam->GetUnfold();
-    if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!"); 
-    fECALocMaxCut           = recParam->GetLocMaxCut();
-    fTimeCut                = recParam->GetTimeCut();
-    fTimeMin                = recParam->GetTimeMin();
-    fTimeMax                = recParam->GetTimeMax();
-
-    AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
-                   fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
-  }
-
-}
-
-//____________________________________________________________________________
 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
 {
        // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
@@ -462,7 +248,6 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d
        
        static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
        static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
-       static Int_t rowdiff, coldiff;
 
        shared = kFALSE;
        
@@ -488,8 +273,8 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d
                
        }//Different SM, same phi
        
-       rowdiff = TMath::Abs(iphi1 - iphi2);  
-       coldiff = TMath::Abs(ieta1 - ieta2) ;  
+       Int_t rowdiff = TMath::Abs(iphi1 - iphi2);  
+       Int_t coldiff = TMath::Abs(ieta1 - ieta2) ;  
 
        // neighbours with at least common side; May 11, 2007
        if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {  
@@ -609,7 +394,7 @@ void AliEMCALClusterizerv1::MakeUnfolding()
       AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
 
       TVector3 gpos;
-      Int_t absId;
+      Int_t absId = -1;
       recPoint->GetGlobalPosition(gpos);
       fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
       if(absId > nModulesToUnfold)
@@ -686,15 +471,14 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
 
   AliEMCALDigit * digit = 0 ;
   Int_t * digitsList = iniTower->GetDigitsList() ;
-
-  Int_t iparam ;
+  
+  Int_t iparam = 0 ;
   Int_t iDigit ;
   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
     digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
     fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
     efit[iDigit] = 0;
 
-    iparam = 0 ;
     while(iparam < nPar ){
       xpar = fitparameters[iparam] ;
       zpar = fitparameters[iparam+1] ;
@@ -710,7 +494,7 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
   // to its contribution to efit
 
   Float_t * energiesList = iniTower->GetEnergiesList() ;
-  Float_t ratio ;
+  Float_t ratio = 0 ;
 
   iparam = 0 ;
   while(iparam < nPar ){
@@ -729,7 +513,7 @@ void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower,
     fNumberOfECAClusters++ ;
     recPoint->SetNExMax((Int_t)nPar/3) ;
 
-    Float_t eDigit ;
+    Float_t eDigit = 0. ;
     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
       digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
       fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
@@ -774,7 +558,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
     for(iparam = 0 ; iparam < nPar ; iparam++)
       Grad[iparam] = 0 ; // Will evaluate gradient
 
-  Double_t efit ;
+  Double_t efit = 0. ;
 
   AliEMCALDigit * digit ;
   Int_t iDigit ;
@@ -791,7 +575,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
 
     if(iflag == 2){  // calculate gradient
       Int_t iParam = 0 ;
-      efit = 0 ;
+      efit = 0. ;
       while(iParam < nPar ){
         Double_t dx = (xDigit - x[iParam]) ;
         iParam++ ;
@@ -837,86 +621,3 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
     // Here we assume, that sigma = sqrt(E) 
   }
 }
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
-{
-  // Print clusterizer parameters
-
-  TString message("\n") ; 
-  
-  if( strcmp(GetName(), "") !=0 ){
-    
-    // Print parameters
-    TString taskName(Version()) ;
-    
-    printf("--------------- "); 
-    printf("%s",taskName.Data()) ; 
-    printf(" "); 
-    printf("Clusterizing digits: "); 
-    printf("\n                       ECA Local Maximum cut    = %f", fECALocMaxCut); 
-    printf("\n                       ECA Logarithmic weight   = %f", fECAW0); 
-    if(fToUnfold)
-      printf("\nUnfolding on\n");
-    else
-      printf("\nUnfolding off\n");
-    
-    printf("------------------------------------------------------------------"); 
-  }
-  else
-    printf("AliEMCALClusterizerv1 not initialized ") ;
-}
-
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
-{
-  // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
-  if(strstr(option,"deb")) {
-    printf("PrintRecPoints: Clusterization result:") ; 
-  
-    printf("           Found %d ECA Rec Points\n ", 
-        fRecPoints->GetEntriesFast()) ; 
-  }
-
-  if(strstr(option,"all")) {
-    if(strstr(option,"deb")) {
-      printf("\n-----------------------------------------------------------------------\n") ;
-      printf("Clusters in ECAL section\n") ;
-      printf("Index    Ene(GeV) Multi Module     GX    GY   GZ  lX    lY   lZ   Dispersion Lambda 1   Lambda 2  # of prim  Primaries list\n") ;
-    }
-   Int_t index =0;
-
-    for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
-      AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ; 
-      TVector3  globalpos;  
-      //rp->GetGlobalPosition(globalpos);
-      TVector3  localpos;  
-      rp->GetLocalPosition(localpos);
-      Float_t lambda[2]; 
-      rp->GetElipsAxis(lambda);
-      Int_t * primaries; 
-      Int_t nprimaries;
-      primaries = rp->GetPrimaries(nprimaries);
-      if(strstr(option,"deb")) 
-      printf("\n%6d  %8.4f  %3d     %4.1f    %4.1f %4.1f  %4.1f %4.1f %4.1f    %4.1f   %4f  %4f    %2d     : ", 
-            rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
-            globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
-            rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
-      if(strstr(option,"deb")){ 
-        for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
-         printf("%d ", primaries[iprimary] ) ; 
-        }
-      }
-    }
-
-    if(strstr(option,"deb"))
-    printf("\n-----------------------------------------------------------------------\n");
-  }
-}
-
-//___________________________________________________________________
-void  AliEMCALClusterizerv1::PrintRecoInfo()
-{
-  printf(" AliEMCALClusterizerv1::PrintRecoInfo() : version %s \n", Version() );
-
-}
index 3efdf53..f7b94eb 100644 (file)
@@ -9,28 +9,20 @@
 //  Implementation version 1 of the clusterization algorithm                     
 //  Performs clusterization (collects neighbouring active cells) and 
 //  unfolding of the clusters with several local maxima.  
-//  results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
-//  PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer
+//  results are stored in TreeR
 //
 //*-- Author: Yves Schutz (SUBATECH)
-// Modif: 
-//  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
-//                           of new  IO (à la PHOS)
+//--          Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
+                        
 
 // --- ROOT system ---
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
-
 #include "AliEMCALClusterizer.h"
-class TH1F;
 class AliEMCALRecPoint ; 
 class AliEMCALDigit ;
-class AliEMCALDigitizer ;
-class AliEMCALGeometry ;
-class AliEMCALCalibData ;
-class AliCaloCalibPedestal ;
 
 class AliEMCALClusterizerv1 : public AliEMCALClusterizer {
   
@@ -43,43 +35,15 @@ public:
   virtual ~AliEMCALClusterizerv1()  ;
 
   virtual Int_t   AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared)const ; 
-                               // Checks if digits are in neighbour cells 
-
-  virtual Float_t Calibrate(const Float_t amp, const Float_t time, const Int_t cellId) ;  // Tranforms Amp to energy 
-
-  virtual void    GetNumberOfClustersFound(int numb )const { numb = fNumberOfECAClusters ;} 
-  virtual Float_t GetECAClusteringThreshold()        const { return fECAClusteringThreshold;}  
-  virtual Float_t GetECALocalMaxCut()                const { return fECALocMaxCut;} 
-  virtual Float_t GetECALogWeight()                  const { return fECAW0;}
-  virtual Float_t GetMinECut()                       const { return fMinECut;}
-
-  virtual Float_t GetTimeCut()                       const { return fTimeCut ; }
-  virtual Float_t GetTimeMin()                       const { return fTimeMin ; }
-  virtual Float_t GetTimeMax()                       const { return fTimeMax ; }
-
+                               // Checks if digits are in neighbour cells
   virtual void    Digits2Clusters(Option_t *option);                // Does the job
 
-  virtual void    Print(Option_t * option)const ;
-
-  virtual void SetECAClusteringThreshold(Float_t cluth)  { fECAClusteringThreshold = cluth ; }
-  virtual void SetMinECut(Float_t mine)                  { fMinECut = mine; }
-  virtual void SetECALocalMaxCut(Float_t cut)            { fECALocMaxCut = cut ; }
-  virtual void SetECALogWeight(Float_t w)                { fECAW0 = w ; }
-  virtual void SetTimeCut(Float_t t)                                    { fTimeCut = t ;}
-  virtual void SetTimeMin(Float_t t)                                    { fTimeMin = t ;}
-  virtual void SetTimeMax(Float_t t)                                    { fTimeMax = t ;}
-
-  virtual void SetUnfolding(Bool_t toUnfold = kTRUE )    {fToUnfold = toUnfold ;}  
   static Double_t ShowerShape(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; 
                                             //class member function (not object member function)
   static void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)  ;
                                             // Chi^2 of the fit. Should be static to be passes to MINUIT
   virtual const char * Version() const { return "clu-v1" ; }  
 
-  void   PrintRecoInfo();                        //*MENU*
-  void   SetCalibrationParameters(AliEMCALCalibData * calib)   { fCalibData = calib ; }
-  void   SetCaloCalibPedestal(AliCaloCalibPedestal  * caloped) { fCaloPed   = caloped ; }
-
 protected:
 
   virtual void   MakeClusters();            
@@ -88,44 +52,19 @@ private:
   AliEMCALClusterizerv1(const AliEMCALClusterizerv1 &); //copy ctor
   AliEMCALClusterizerv1 & operator = (const AliEMCALClusterizerv1 &);
 
-  void    GetCalibrationParameters(void) ;
-  void    GetCaloCalibPedestal(void) ;
 
   Bool_t  FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** MaxAt, const Float_t * maxAtEnergy, 
                  Int_t NPar, Float_t * FitParametres) const; //Used in UnfoldClusters, calls TMinuit
-  void Init() ;
-  void InitParameters();
 
   virtual void   MakeUnfolding();
   void           UnfoldCluster(AliEMCALRecPoint * iniEmc, Int_t Nmax, 
-                              AliEMCALDigit ** maxAt,
-                              Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package
-  void           PrintRecPoints(Option_t * option) ;
+                               AliEMCALDigit ** maxAt,
+                               Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package
 
 private:
-  AliEMCALGeometry* fGeom;           //! pointer to geometry for utilities
-
-  Bool_t  fDefaultInit;              //! Says if the task was created by defaut ctor (only parameters are initialized)
-  Bool_t  fToUnfold ;                // To perform unfolding 
-  Int_t   fNumberOfECAClusters ;     // number of clusters found in EC section
-
-  //Calibration parameters... to be replaced by database 
-
-  AliEMCALCalibData    * fCalibData ;   //! Calibration database if aval
-  AliCaloCalibPedestal * fCaloPed   ;   //! Tower status map if aval
 
-  Float_t fADCchannelECA ;          // width of one ADC channel for EC section (GeV)
-  Float_t fADCpedestalECA ;         // pedestal of ADC for EC section (GeV) 
-  Float_t fECAClusteringThreshold ;  // minimum energy to seed a EC digit in a cluster
-  Float_t fECALocMaxCut ;            // minimum energy difference to distinguish local maxima in a cluster
-  Float_t fECAW0 ;                   // logarithmic weight for the cluster center of gravity calculation
-  Float_t fTimeCut ;                 // Maximum time difference between the digits inside EMC cluster
-  Float_t fTimeMin ;                 // Minimum time of physical signal in a cell/digiy
-  Float_t fTimeMax ;                 // Maximum time of physical signal in a cell/digit
-  Float_t fMinECut;                  // Minimum energy for a digit to be a member of a cluster
 
-  ClassDef(AliEMCALClusterizerv1,9)   // Clusterizer implementation version 1
+  ClassDef(AliEMCALClusterizerv1,10)   // Clusterizer implementation version 1
 
 };
 
index 20ec7cb..ca73559 100644 (file)
@@ -34,7 +34,7 @@
 
 ClassImp(AliEMCALRecParam)
   
-  TObjArray* AliEMCALRecParam::fgkMaps =0; //ALTRO mappings 
+TObjArray* AliEMCALRecParam::fgkMaps =0; //ALTRO mappings 
 
 AliEMCALRecParam::AliEMCALRecParam() :
   AliDetectorRecoParam(),
@@ -46,6 +46,7 @@ AliEMCALRecParam::AliEMCALRecParam() :
   fTimeCut(1.),// high value, accept all
   fTimeMin(-1.),// small value, accept all
   fTimeMax(1.),// high value, accept all//clustering
+  fClusterizerFlag(AliEMCALRecParam::kClusterizerv1),
   fTrkCutX(6.0), 
   fTrkCutY(6.0), 
   fTrkCutZ(6.0),  
@@ -86,7 +87,7 @@ AliEMCALRecParam::AliEMCALRecParam() :
   // Pb Pb
   
   // as a first step, all array elements are initialized to 0.0
-  Int_t i, j;
+  Int_t i=0, j=0;
   for (i = 0; i < 6; i++) {
     for (j = 0; j < 6; j++) {      
       fGamma[i][j] =  fPiZero[i][j] = fHadron[i][j] = 0.; 
@@ -234,6 +235,7 @@ AliEMCALRecParam::AliEMCALRecParam(const AliEMCALRecParam& rp) :
   fTimeCut(rp.fTimeCut), 
   fTimeMin(rp.fTimeMin),
   fTimeMax(rp.fTimeMax),//clustering
+  fClusterizerFlag(rp.fClusterizerFlag),
   fTrkCutX(rp.fTrkCutX), 
   fTrkCutY(rp.fTrkCutY), 
   fTrkCutZ(rp.fTrkCutZ),  
@@ -256,7 +258,7 @@ AliEMCALRecParam::AliEMCALRecParam(const AliEMCALRecParam& rp) :
   //copy constructor
   
   //PID values
-  Int_t i, j;
+  Int_t i=0, j=0;
   for (i = 0; i < 6; i++) {
     for (j = 0; j < 6; j++) {
       fGamma[i][j]       = rp.fGamma[i][j];
@@ -287,6 +289,7 @@ AliEMCALRecParam& AliEMCALRecParam::operator = (const AliEMCALRecParam& rp)
     fTimeCut   = rp.fTimeCut;
     fTimeMax   = rp.fTimeMax;
     fTimeMin   = rp.fTimeMin;//clustering
+    fClusterizerFlag = rp.fClusterizerFlag;
     fTrkCutX   = rp.fTrkCutX;
     fTrkCutY   = rp.fTrkCutY;
     fTrkCutZ   = rp.fTrkCutZ;
@@ -307,7 +310,7 @@ AliEMCALRecParam& AliEMCALRecParam::operator = (const AliEMCALRecParam& rp)
     fFitLEDEvents      = rp.fFitLEDEvents;//raw signal
          
     //PID values
-    Int_t i, j;
+    Int_t i=0, j=0;
     for (i = 0; i < 6; i++) {
       for (j = 0; j < 6; j++) {
        fGamma[i][j]       = rp.fGamma[i][j];
@@ -383,7 +386,7 @@ AliEMCALRecParam* AliEMCALRecParam::GetLowFluxParam()
   
   //PID parameters for pp  implemented 
   // as a first step, all array elements are initialized to 0.0
-  Int_t i, j;
+  Int_t i=0, j=0;
   for (i = 0; i < 6; i++) {
     for (j = 0; j < 6; j++) {
       params->SetGamma(i,j,0.);
@@ -554,6 +557,7 @@ void AliEMCALRecParam::Print(Option_t * opt) const
   // if nothing is specified print all, if "reco" just clusterization/track matching
   // if "pid", just PID parameters, if "raw", just raw utils parameters.
   if(!strcmp("",opt) || !strcmp("reco",opt)){
+    AliInfo(Form("Clusterizer selected: %d", fClusterizerFlag));
     AliInfo(Form("Clusterization parameters :\n fClusteringThreshold=%.3f,\n fW0=%.3f,\n fMinECut=%.3f,\n fUnfold=%d,\n fLocMaxCut=%.3f,\n fTimeCut=%2.1f ns\n fTimeMin=%2.1f ns\n fTimeMax=%2.1f ns\n",
                 fClusteringThreshold,fW0,fMinECut,fUnfold,fLocMaxCut,fTimeCut*1.e9,fTimeMin*1e9,fTimeMax*1e9));
     
index 48af057..3e09c71 100644 (file)
 class AliEMCALRecParam : public AliDetectorRecoParam
 {
  public:
+
+  enum AliEMCALClusterizerFlag
+  {
+    kClusterizerv1  = 0,
+    kClusterizerNxN = 1   
+  };
   
   AliEMCALRecParam() ;
   AliEMCALRecParam(const AliEMCALRecParam& recParam);
@@ -113,16 +119,19 @@ class AliEMCALRecParam : public AliDetectorRecoParam
   Bool_t   UseFALTRO()            const {return fUseFALTRO; }
   Bool_t   FitLEDEvents()         const {return fFitLEDEvents; }
 
-       virtual void Print(Option_t * option="") const ;
+  virtual void Print(Option_t * option="") const ;
   
   static AliEMCALRecParam* GetDefaultParameters();
   static AliEMCALRecParam* GetLowFluxParam();
   static AliEMCALRecParam* GetHighFluxParam();
   static AliEMCALRecParam* GetCalibParam();
   static AliEMCALRecParam* GetCosmicParam();
-  
+
   static const  TObjArray* GetMappings();
   
+  void    SetClusterizerFlag(Short_t val) { fClusterizerFlag = val;  }
+  Short_t GetClusterizerFlag() const      { return fClusterizerFlag; }
+  
  private:
   //Clustering
   Float_t fClusteringThreshold ; // Minimum energy to seed a EC digit in a cluster
@@ -133,6 +142,7 @@ class AliEMCALRecParam : public AliDetectorRecoParam
   Float_t fTimeCut ;             // Maximum time of digits with respect to EMC cluster max.
   Float_t fTimeMin ;             // Minimum time of digits
   Float_t fTimeMax ;             // Maximum time of digits
+  Short_t fClusterizerFlag ;     // Choice of the clusterizer; Default selection (v1) is zero
 
   //PID (Guenole)
   Double_t fGamma[6][6];         // Parameter to Compute PID for photons     
@@ -169,7 +179,7 @@ class AliEMCALRecParam : public AliDetectorRecoParam
        
   static TObjArray* fgkMaps;       // ALTRO mappings for RCU0..RCUX
   
-  ClassDef(AliEMCALRecParam,12)     // Reconstruction parameters
+  ClassDef(AliEMCALRecParam,13)     // Reconstruction parameters
     
     } ;
 
index 63094c3..048721c 100644 (file)
@@ -45,6 +45,7 @@
 #include "AliEMCALRawUtils.h"
 #include "AliEMCALDigit.h"
 #include "AliEMCALClusterizerv1.h"
+#include "AliEMCALClusterizerNxN.h"
 #include "AliEMCALRecPoint.h"
 #include "AliEMCALPID.h"
 #include "AliEMCALTrigger.h"
@@ -106,13 +107,11 @@ AliEMCALReconstructor::AliEMCALReconstructor()
                if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
     }
        
-       if(!fPedestalData)
-               AliFatal("Dead map not found in CDB!");
-       
-       
-  //Init the clusterizer with geometry and calibration pointers, avoid doing it twice.
-  fgClusterizer = new AliEMCALClusterizerv1(fGeom, fCalibData,fPedestalData); 
+  if(!fPedestalData)
+    AliFatal("Dead map not found in CDB!");
        
+  InitClusterizer();
+               
   if(!fGeom) AliFatal(Form("Could not get geometry!"));
 
   AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance();
@@ -143,6 +142,32 @@ AliEMCALReconstructor::~AliEMCALReconstructor()
 // }
 
 //____________________________________________________________________________
+void AliEMCALReconstructor::InitClusterizer() 
+{
+  //Init the clusterizer with geometry and calibration pointers, avoid doing it twice.
+  
+  AliEMCALRecParam *recParam = NULL;
+  AliCDBEntry *entry = (AliCDBEntry*) 
+  AliCDBManager::Instance()->Get("EMCAL/Calib/RecoParam");
+  //Get The reco param for the default event specie
+  if (entry) 
+   recParam = (AliEMCALRecParam*)((TObjArray *) entry->GetObject())->At(0);
+    
+  if(!recParam)  
+    AliFatal("RecoParam not found in CDB!");
+    
+  if (recParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
+  {
+    fgClusterizer = new AliEMCALClusterizerv1(fGeom, fCalibData,fPedestalData); 
+  }
+  else
+  {
+    fgClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData,fPedestalData); 
+  }
+  
+}
+
+//____________________________________________________________________________
 void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
 {
   // method called by AliReconstruction; 
@@ -154,6 +179,7 @@ void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree)
   AliCodeTimerAuto("",0)
 
   ReadDigitsArrayFromTree(digitsTree);
+  
   fgClusterizer->InitParameters();
   fgClusterizer->SetOutput(clustersTree);
  
@@ -287,7 +313,7 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   for (Int_t idig = 0 ; idig < nDigits ; idig++) {
     const AliEMCALDigit * dig = (const AliEMCALDigit*)digits->At(idig);
     if(dig->GetAmplitude() > 0 ){
-         energy = (static_cast<AliEMCALClusterizerv1*> (fgClusterizer))->Calibrate(dig->GetAmplitude(),dig->GetTime(),dig->GetId()); //TimeR or Time?
+         energy = fgClusterizer->Calibrate(dig->GetAmplitude(),dig->GetTime(),dig->GetId()); //TimeR or Time?
          if(energy > 0){ //Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them  
                  emcCells.SetCell(idignew,dig->GetId(),energy, dig->GetTime());   
                  idignew++;
@@ -381,13 +407,13 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       ec->SetPosition(xyz);
       ec->SetE(clust->GetEnergy());
                
-         //Distance to the nearest bad crystal
-         ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower()); 
+      //Distance to the nearest bad crystal
+      ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower()); 
 
       ec->SetNCells(newCellMult);
       //Change type of list from short to ushort
       UShort_t *newAbsIdList  = new UShort_t[newCellMult];
-      Double_t *newFracList  = new Double_t[newCellMult];
+      Double_t *newFracList   = new Double_t[newCellMult];
       for(Int_t i = 0; i < newCellMult ; i++) {
         newAbsIdList[i]=absIdList[i];
         newFracList[i]=fracList[i];
index 9b60075..73b6f72 100644 (file)
@@ -54,7 +54,8 @@ public:
   virtual ~AliEMCALReconstructor() ; //dtor
 
   virtual  void Init() {;}
-
+  virtual  void InitClusterizer();
+  
   Bool_t       Debug() const { return fDebug ; }
 
   using AliReconstructor::FillESD;
@@ -63,7 +64,7 @@ public:
   AliTracker*  CreateTracker () const 
   {return new AliEMCALTracker;} 
   using AliReconstructor::Reconstruct;
-  virtual void Reconstruct(TTree* digitsTree, TTree* clustersTree) const;
+  virtual void Reconstruct(TTree* digitsTree, TTree* clustersTree) const ;
 
   virtual Bool_t             HasDigitConversion() const {return kTRUE;};
   virtual void               ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const;
index 6979aaf..2053121 100644 (file)
@@ -6,6 +6,7 @@
 #pragma link C++ class AliEMCALReconstructor+;
 #pragma link C++ class AliEMCALClusterizer+;
 #pragma link C++ class AliEMCALClusterizerv1+;
+#pragma link C++ class AliEMCALClusterizerNxN+;
 #pragma link C++ class AliEMCALTrack+;
 #pragma link C++ class AliEMCALTracker+;
 #pragma link C++ class AliEMCALPID+;
index 6ce12bf..f1e9258 100644 (file)
@@ -5,6 +5,7 @@ SRCS =\
 AliEMCALReconstructor.cxx \
 AliEMCALClusterizer.cxx \
 AliEMCALClusterizerv1.cxx \
+AliEMCALClusterizerNxN.cxx \
 AliEMCALTrack.cxx \
 AliEMCALTracker.cxx \
 AliEMCALPID.cxx \
@@ -20,8 +21,6 @@ DHDR:=EMCALrecLinkDef.h
 # have to tune
 EINCLUDE:=PYTHIA6 RAW EVGEN THijing VZERO
 
-EXPORT:=AliEMCALTrack.h
-
 ifeq (win32gcc,$(ALICE_TARGET))
 PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) \
                          -lEMCALUtils -lEMCALbase -lEMCALsim -lSTEER -lESD -lCDB -lSTEERBase \
index c9229c6..d1f3b54 100644 (file)
@@ -127,6 +127,8 @@ AliEMCALRecParam* GetHighMultiplicityParameters()
   //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
   
   AliEMCALRecParam* params =  AliEMCALRecParam::GetDefaultParameters();
+  //params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN);
+  params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1);
   
   //Clusterization
   params->SetClusteringThreshold(0.5);
@@ -306,6 +308,9 @@ AliEMCALRecParam* GetLowMultiplicityParameters()
   //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetLowFluxParam();
   
   AliEMCALRecParam* params =  AliEMCALRecParam::GetDefaultParameters();
+  //params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerNxN);
+  params->SetClusterizerFlag(AliEMCALRecParam::kClusterizerv1);
+  
   params->SetClusteringThreshold(0.1); // 100 MeV                                             
   params->SetMinECut(0.01);  //10 MeV          
   params->SetTimeCut(250e-9);//250 ns
index 00176fe..fe225fe 100644 (file)
@@ -4,7 +4,7 @@
 // Author: Gustavo Conesa (INFN-LNF)
 
 
-void PrintEMCALRecParam(char * file = "Run0_999999999_v0_s0.root"){
+void PrintEMCALRecParam(char * file = "$ALICE_ROOT/OCDB/EMCAL/Calib/RecoParam/Run0_999999999_v0_s0.root"){
 
 
 TFile * f = new TFile(file,"READ");
@@ -24,9 +24,10 @@ cout<<"============== "<<rparam->GetName()<<" ==============="<<endl;
 
 cout<<"================================================"<<endl;
 
-
-rparam->Print("");
-
+//rparam->Print("reco");//Print only clusterizer parameters
+//rparam->Print("pid");//Print only pid parameters
+//rparam->Print("raw");//Print only raw digitization parameters
+rparam->Print("");// Print all
 
 
 }