Clusterizer introduced
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jan 2002 07:43:19 +0000 (07:43 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jan 2002 07:43:19 +0000 (07:43 +0000)
EMCAL/AliEMCALClusterizer.cxx [new file with mode: 0644]
EMCAL/AliEMCALClusterizer.h [new file with mode: 0644]
EMCAL/AliEMCALClusterizerv1.cxx [new file with mode: 0644]
EMCAL/AliEMCALClusterizerv1.h [new file with mode: 0644]
EMCAL/EMCALLinkDef.h
EMCAL/Makefile

diff --git a/EMCAL/AliEMCALClusterizer.cxx b/EMCAL/AliEMCALClusterizer.cxx
new file mode 100644 (file)
index 0000000..325aef6
--- /dev/null
@@ -0,0 +1,54 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//  Base class for the clusterization algorithm (pure abstract)
+//*--
+//*-- Author: Yves Schutz  SUBATECH 
+//////////////////////////////////////////////////////////////////////////////
+
+// --- ROOT system ---
+
+
+// --- Standard library ---
+
+
+
+// --- AliRoot header files ---
+
+#include "AliEMCALClusterizer.h"
+
+ClassImp(AliEMCALClusterizer)
+
+//____________________________________________________________________________
+  AliEMCALClusterizer::AliEMCALClusterizer():TTask("","")
+{
+  // ctor
+}
+//____________________________________________________________________________
+AliEMCALClusterizer::AliEMCALClusterizer(const char* headerFile, const char* name):
+TTask(name, headerFile)
+{
+  // ctor
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizer::~AliEMCALClusterizer()
+{
+  // dtor
+}
diff --git a/EMCAL/AliEMCALClusterizer.h b/EMCAL/AliEMCALClusterizer.h
new file mode 100644 (file)
index 0000000..04583eb
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef ALIEMCALCLUSTERIZER_H
+#define ALIEMCALCLUSTERIZER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+                            
+/* $Id$ */
+
+//_________________________________________________________________________
+//  Base class for the clusterization algorithm (pure abstract)
+//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
+// --- ROOT system ---
+
+#include "TTask.h" 
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+//#include "AliEMCALDigit.h"
+
+class AliEMCALClusterizer : public TTask {
+
+public:
+
+  AliEMCALClusterizer() ;        // default ctor
+  AliEMCALClusterizer(const char * headerFile, const char * name) ;
+  virtual ~AliEMCALClusterizer() ; // dtor
+
+  virtual Float_t GetEmcClusteringThreshold()const = 0 ; 
+  virtual Float_t GetEmcLocalMaxCut()const = 0 ; 
+  virtual Float_t GetEmcLogWeight()const = 0 ; 
+  virtual Float_t GetTimeGate() const = 0 ;
+  virtual Float_t GetCpvClusteringThreshold()const = 0 ; 
+  virtual Float_t GetCpvLocalMaxCut()const = 0 ; 
+  virtual Float_t GetCpvLogWeight()const = 0 ; 
+  virtual char *  GetRecPointsBranch() const = 0 ;
+  virtual const Int_t GetRecPointsInRun()  const = 0 ; 
+  virtual char *  GetDigitsBranch() const = 0 ;
+
+  virtual void MakeClusters() = 0 ; 
+  virtual void Print(Option_t * option)const = 0;
+
+  virtual void SetTowerClusteringThreshold(Float_t cluth) = 0 ; 
+  virtual void SetTowerLocalMaxCut(Float_t cut) = 0 ; 
+  virtual void SetTowerLogWeight(Float_t w) = 0 ; 
+  virtual void SetTimeGate(Float_t gate) = 0 ;
+  virtual void SetPreShoClusteringThreshold(Float_t cluth) = 0 ; 
+  virtual void SetPreShoLocalMaxCut(Float_t cut) = 0 ; 
+  virtual void SetPreShoLogWeight(Float_t w) = 0 ; 
+  virtual void SetDigitsBranch(const char * title) = 0 ;
+  virtual void SetRecPointsBranch(const char *title) = 0 ;
+  virtual void SetUnfolding(Bool_t toUnfold ) = 0 ;
+  virtual const char * Version() const = 0 ;  
+
+  ClassDef(AliEMCALClusterizer,1)  // Clusterization algorithm class 
+
+} ;
+
+#endif // AliEMCALCLUSTERIZER_H
diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx
new file mode 100644 (file)
index 0000000..fbfd5aa
--- /dev/null
@@ -0,0 +1,1094 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+/* $Log:
+   1 October 2000. Yuri Kharlov:
+     AreNeighbours()
+     PPSD upper layer is considered if number of layers>1
+
+   18 October 2000. Yuri Kharlov:
+     AliEMCALClusterizerv1()
+     CPV clusterizing parameters added
+
+     MakeClusters()
+     After first PPSD digit remove EMC digits only once
+*/
+//*-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
+//////////////////////////////////////////////////////////////////////////////
+//  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 "TROOT.h" 
+#include "TFile.h" 
+#include "TFolder.h" 
+#include "TMath.h" 
+#include "TMinuit.h"
+#include "TTree.h" 
+#include "TSystem.h" 
+#include "TBenchmark.h"
+
+// --- Standard library ---
+
+#include <iostream.h>
+#include <iomanip.h>
+
+// --- AliRoot header files ---
+
+#include "AliEMCALClusterizerv1.h"
+#include "AliEMCALDigit.h"
+#include "AliEMCALDigitizer.h"
+#include "AliEMCALTowerRecPoint.h"
+#include "AliEMCAL.h"
+#include "AliEMCALGetter.h"
+#include "AliRun.h"
+
+ClassImp(AliEMCALClusterizerv1)
+  
+//____________________________________________________________________________
+  AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
+{
+  // default ctor (to be used mainly by Streamer)
+  
+  fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ; 
+  
+  fPreShoClusteringThreshold = 0.0001;
+  fTowerClusteringThreshold  = 0.2;   
+  
+  fTowerLocMaxCut  = 0.03 ;
+  fPreShoLocMaxCut = 0.03 ;
+  
+  fW0     = 4.5 ;
+  fW0CPV  = 4.0 ;
+
+  fTimeGate = 1.e-8 ; 
+
+  fToUnfold = kFALSE ;
+
+  fHeaderFileName = fDigitsBranchTitle  = "" ;
+  fRecPointsInRun          = 0 ;   
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile,const char* name)
+:AliEMCALClusterizer(headerFile, name)
+{
+  // ctor with the indication of the file where header Tree and digits Tree are stored
+  
+
+  fNumberOfPreShoClusters = fNumberOfTowerClusters = 0 ; 
+
+  
+  fPreShoClusteringThreshold  = 0.0001;
+  fTowerClusteringThreshold   = 0.2;   
+  
+  fTowerLocMaxCut  = 0.03 ;
+  fPreShoLocMaxCut = 0.03 ;
+  
+  fW0     = 4.5 ;
+  fW0CPV  = 4.0 ;
+
+  fTimeGate = 1.e-8 ; 
+  
+  fToUnfold = kFALSE ;
+  
+  fHeaderFileName     = GetTitle() ; 
+  fDigitsBranchTitle  = GetName() ;
+  
+  TString clusterizerName( GetName()) ; 
+  clusterizerName.Append(":") ; 
+  clusterizerName.Append(Version()) ; 
+  SetName(clusterizerName) ;
+  fRecPointsInRun          = 0 ; 
+
+  Init() ;
+
+}
+//____________________________________________________________________________
+  AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
+{
+}
+//____________________________________________________________________________
+Float_t  AliEMCALClusterizerv1::Calibrate(Int_t amp, Bool_t inpresho) const
+{
+  if ( inpresho ) // calibrate as pre shower
+     return -fADCpedestalPreSho + amp * fADCchannelPreSho ; 
+
+  else //calibrate as tower 
+    return -fADCpedestalTower + amp * fADCchannelTower ;                
+}
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::Exec(Option_t * option)
+{
+  // Steering method
+
+  if( strcmp(GetName(), "")== 0 ) 
+    Init() ;
+
+  if(strstr(option,"tim"))
+    gBenchmark->Start("EMCALClusterizer"); 
+  
+  if(strstr(option,"print"))
+    Print("") ; 
+
+  gAlice->GetEvent(0) ;
+  
+ //check, if the branch with name of this" already exits?
+  TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
+  TIter next(lob) ; 
+  TBranch * branch = 0 ;  
+  Bool_t emcaltowerfound = kFALSE, emcalpreshofound = kFALSE, clusterizerfound = kFALSE ; 
+
+  TString branchname = GetName() ;
+  branchname.Remove(branchname.Index(Version())-1) ;
+  
+  while ( (branch = (TBranch*)next()) && (!emcaltowerfound || !emcalpreshofound || !clusterizerfound) ) {
+    if ( (strcmp(branch->GetName(), "EMCALTowerRP")==0) && (strcmp(branch->GetTitle(), branchname.Data())==0) ) 
+      emcaltowerfound = kTRUE ;
+    
+    else if ( (strcmp(branch->GetName(), "EMCALPreShoRP")==0) && (strcmp(branch->GetTitle(), branchname.Data())==0) ) 
+      emcalpreshofound = kTRUE ;
+   
+    else if ((strcmp(branch->GetName(), "AliEMCALClusterizer")==0) && (strcmp(branch->GetTitle(), GetName())==0) ) 
+      clusterizerfound = kTRUE ; 
+  }
+
+  if ( emcalpreshofound || emcaltowerfound || clusterizerfound ) {
+    cerr << "WARNING: AliEMCALClusterizer::Exec -> Tower(PreSho)RecPoints and/or Clusterizer branch with name " 
+        << branchname.Data() << " already exits" << endl ;
+    return ; 
+  }       
+
+  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+  Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
+  Int_t ievent ;
+
+  for(ievent = 0; ievent < nevents; ievent++){
+
+    if(ievent == 0)
+      GetCalibrationParameters() ;
+
+    fNumberOfTowerClusters = fNumberOfPreShoClusters = 0 ;
+   
+    gime->Event(ievent,"D") ;
+        
+    MakeClusters() ;
+    
+    if(fToUnfold)
+      MakeUnfolding() ;
+
+    WriteRecPoints(ievent) ;
+
+    if(strstr(option,"deb"))  
+      PrintRecPoints(option) ;
+
+    //increment the total number of digits per run 
+    fRecPointsInRun += gime->TowerRecPoints()->GetEntriesFast() ;  
+    fRecPointsInRun += gime->PreShowerRecPoints()->GetEntriesFast() ;  
+ }
+  
+  if(strstr(option,"tim")){
+    gBenchmark->Stop("EMCALClusterizer");
+    cout << "AliEMCALClusterizer:" << endl ;
+    cout << "  took " << gBenchmark->GetCpuTime("EMCALClusterizer") << " seconds for Clusterizing " 
+        <<  gBenchmark->GetCpuTime("EMCALClusterizer")/nevents << " seconds per event " << endl ;
+    cout << endl ;
+  }
+  
+}
+
+//____________________________________________________________________________
+Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
+                                   Int_t nPar, Float_t * fitparameters) const
+{ 
+  // Calls TMinuit to fit the energy distribution of a cluster with several maxima 
+  // The initial values for fitting procedure are set equal to the positions of local maxima.
+  // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
+
+  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
+  TClonesArray * digits = gime->Digits() ; 
+  
+
+  gMinuit->mncler();                     // Reset Minuit's list of paramters
+  gMinuit->SetPrintLevel(-1) ;           // No Printout
+  gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;  
+                                         // To set the address of the minimization function 
+
+  TList * toMinuit = new TList();
+  toMinuit->AddAt(emcRP,0) ;
+  toMinuit->AddAt(digits,1) ;
+  
+  gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
+
+  // filling initial values for fit parameters
+  AliEMCALDigit * digit ;
+
+  Int_t ierflg  = 0; 
+  Int_t index   = 0 ;
+  Int_t nDigits = (Int_t) nPar / 3 ;
+
+  Int_t iDigit ;
+
+  AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
+
+  for(iDigit = 0; iDigit < nDigits; iDigit++){
+    digit = (AliEMCALDigit *) maxAt[iDigit]; 
+
+    Int_t relid[4] ;
+    Float_t x = 0.;
+    Float_t z = 0.;
+    geom->AbsToRelNumbering(digit->GetId(), relid) ;
+    geom->PosInAlice(relid, x, z) ;
+
+    Float_t energy = maxAtEnergy[iDigit] ;
+
+    gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
+    index++ ;   
+    if(ierflg != 0){ 
+      cout << "EMCAL Unfolding>  Unable to set initial value for fit procedure : x = " << x << endl ;
+      return kFALSE;
+    }
+    gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
+    index++ ;   
+    if(ierflg != 0){
+      cout << "EMCAL Unfolding>  Unable to set initial value for fit procedure : z = " << z << endl ;
+      return kFALSE;
+    }
+    gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
+    index++ ;   
+    if(ierflg != 0){
+      cout << "EMCAL Unfolding>  Unable to set initial value for fit procedure : energy = " << energy << endl ;      
+      return kFALSE;
+    }
+  }
+
+  Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
+                      //  depends on it. 
+  Double_t p1 = 1.0 ;
+  Double_t p2 = 0.0 ;
+
+  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
+  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
+  gMinuit->SetMaxIterations(5);
+  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
+
+  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
+
+  if(ierflg == 4){  // Minimum not found   
+    cout << "EMCAL Unfolding>  Fit not converged, cluster abandoned "<< endl ;      
+    return kFALSE ;
+  }            
+  for(index = 0; index < nPar; index++){
+    Double_t err ;
+    Double_t val ;
+    gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
+    fitparameters[index] = val ;
+   }
+
+  delete toMinuit ;
+  return kTRUE;
+
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::GetCalibrationParameters() 
+{
+  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
+  TString branchname = GetName() ;
+  branchname.Remove(branchname.Index(Version())-1) ;
+
+  AliEMCALDigitizer * dig = gime->Digitizer(branchname) ;
+
+  fADCchannelTower   = dig->GetTowerchannel() ;
+  fADCpedestalTower  = dig->GetTowerpedestal();
+
+  fADCchannelPreSho  = dig->GetPreShochannel() ;
+  fADCpedestalPreSho = dig->GetPreShopedestal() ; 
+
+}
+//____________________________________________________________________________
+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
+  
+  if ( strcmp(GetTitle(), "") == 0 )
+    SetTitle("galice.root") ;
+
+  TString branchname = GetName() ;
+  branchname.Remove(branchname.Index(Version())-1) ;
+
+  AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname, "update") ; 
+  if ( gime == 0 ) {
+    cerr << "ERROR: AliEMCALClusterizerv1::Init -> Could not obtain the Getter object !" << endl ; 
+    return ;
+  } 
+    
+  const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
+  fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
+
+  if(!gMinuit) 
+    gMinuit = new TMinuit(100) ;
+
+  gime->PostClusterizer(this) ;
+  // create a folder on the white board 
+  gime->PostRecPoints(branchname ) ;
+
+  gime->PostDigits(branchname) ;
+  gime->PostDigitizer(branchname) ;
+  
+}
+
+//____________________________________________________________________________
+Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
+{
+  // 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 a 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  
+
+   AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ;
+
+  Int_t rv = 0 ; 
+
+  Int_t relid1[4] ; 
+  geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
+
+  Int_t relid2[4] ; 
+  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
+  if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same EMCAL Arm 
+    Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
+    Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;  
+    
+    if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
+      if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
+      rv = 1 ; 
+    }
+    else {
+      if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) 
+       rv = 2; //  Difference in row numbers is too large to look further 
+    }
+
+  } 
+  else {
+    
+    if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )  
+      rv=2 ;
+
+  }
+
+  return rv ; 
+}
+
+
+//____________________________________________________________________________
+Bool_t AliEMCALClusterizerv1::IsInTower(AliEMCALDigit * digit) const
+{
+  // Tells if (true) or not (false) the digit is in a EMCAL-Tower 
+  Bool_t rv = kFALSE ; 
+  if (!digit->IsInPreShower()) 
+    rv = kTRUE; 
+  return rv ; 
+}
+
+//____________________________________________________________________________
+Bool_t AliEMCALClusterizerv1::IsInPreShower(AliEMCALDigit * digit) const
+{
+  // Tells if (true) or not (false) the digit is in a EMCAL-PreShower
+  Bool_t rv = kFALSE ; 
+  if (digit->IsInPreShower()) 
+    rv = kTRUE; 
+  return rv ; 
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
+{
+
+  // Creates new branches with given title
+  // fills and writes into TreeR.
+
+  TString branchName(GetName() ) ;
+  branchName.Remove(branchName.Index(Version())-1) ;
+
+  AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; 
+  TObjArray * towerRecPoints = gime->TowerRecPoints(branchName) ; 
+  TObjArray * preshoRecPoints = gime->PreShowerRecPoints(branchName) ; 
+  TClonesArray * digits = gime->Digits(branchName) ; 
+
+  Int_t index ;
+  //Evaluate position, dispersion and other RecPoint properties...
+  for(index = 0; index < towerRecPoints->GetEntries(); index++)
+    (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->EvalAll(fW0,digits) ;
+
+  towerRecPoints->Sort() ;
+
+  for(index = 0; index < towerRecPoints->GetEntries(); index++)
+    (dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(index)))->SetIndexInList(index) ;
+
+  towerRecPoints->Expand(towerRecPoints->GetEntriesFast()) ; 
+
+  //Now the same for CPV
+  for(index = 0; index < preshoRecPoints->GetEntries(); index++)
+    (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->EvalAll(fW0CPV,digits)  ;
+
+  preshoRecPoints->Sort() ;
+
+  for(index = 0; index < preshoRecPoints->GetEntries(); index++)
+    (dynamic_cast<AliEMCALRecPoint *>(preshoRecPoints->At(index)))->SetIndexInList(index) ;
+
+  preshoRecPoints->Expand(preshoRecPoints->GetEntriesFast()) ;
+  
+  //Make branches in TreeR for RecPoints and Clusterizer
+  char * filename = 0;
+  if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){   //generating file name
+    filename = new char[strlen(gAlice->GetBaseFile())+20] ;
+    sprintf(filename,"%s/EMCAL.Reco.root",gAlice->GetBaseFile()) ;
+  }
+  
+  //Make new branches
+  TDirectory *cwd = gDirectory;
+  
+  Int_t bufferSize = 32000 ;    
+  Int_t splitlevel = 0 ;
+
+  //First EMC
+  TBranch * emcBranch = gAlice->TreeR()->Branch("EMCALTowerRP","TObjArray",&towerRecPoints,bufferSize,splitlevel);
+  emcBranch->SetTitle(branchName);
+  if (filename) {
+    emcBranch->SetFile(filename);
+    TIter next( emcBranch->GetListOfBranches());
+    TBranch * sb ;
+    while ((sb=(TBranch*)next())) {
+      sb->SetFile(filename);
+    }   
+    
+    cwd->cd();
+  }
+    
+  //Now CPV branch
+  TBranch * cpvBranch = gAlice->TreeR()->Branch("EMCALPreShoRP","TObjArray",&preshoRecPoints,bufferSize,splitlevel);
+  cpvBranch->SetTitle(branchName);
+  if (filename) {
+    cpvBranch->SetFile(filename);
+    TIter next( cpvBranch->GetListOfBranches());
+    TBranch * sb;
+    while ((sb=(TBranch*)next())) {
+      sb->SetFile(filename);
+    }   
+    cwd->cd();
+  } 
+    
+  //And Finally  clusterizer branch
+  AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(branchName) ;
+  TBranch * clusterizerBranch = gAlice->TreeR()->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1",
+                                             &cl,bufferSize,splitlevel);
+  clusterizerBranch->SetTitle(branchName);
+  if (filename) {
+    clusterizerBranch->SetFile(filename);
+    TIter next( clusterizerBranch->GetListOfBranches());
+    TBranch * sb ;
+    while ((sb=(TBranch*)next())) {
+      sb->SetFile(filename);
+    }   
+    cwd->cd();
+  }
+  emcBranch        ->Fill() ;
+  cpvBranch        ->Fill() ;
+  clusterizerBranch->Fill() ;
+
+  gAlice->TreeR()->Write(0,kOverwrite) ;  
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::MakeClusters()
+{
+  // Steering method to construct the clusters stored in a list of Reconstructed Points
+  // A cluster is defined as a list of neighbour digits
+  
+  TString branchName(GetName()) ; 
+  branchName.Remove(branchName.Index(Version())-1) ; 
+  
+  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
+  
+  TObjArray * towerRecPoints  = gime->TowerRecPoints(branchName) ; 
+  TObjArray * preshoRecPoints = gime->PreShowerRecPoints(branchName) ; 
+  towerRecPoints->Delete() ;
+  preshoRecPoints->Delete() ;
+  
+  TClonesArray * digits = gime->Digits(branchName) ; 
+  TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
+  
+  
+  // Clusterization starts  
+  
+  TIter nextdigit(digitsC) ; 
+  AliEMCALDigit * digit ; 
+  Bool_t notremoved = kTRUE ;
+  
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
+    AliEMCALRecPoint * clu = 0 ; 
+    
+    TArrayI clusterdigitslist(1500) ;   
+    Int_t index ;
+    if (( IsInTower (digit)  && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fTowerClusteringThreshold  ) || 
+        ( IsInPreShower (digit) && Calibrate(digit->GetAmp(),digit->IsInPreShower()) > fPreShoClusteringThreshold  ) ) {
+      
+      Int_t iDigitInCluster = 0 ; 
+      
+      if  ( IsInTower(digit) ) {   
+       // start a new Tower RecPoint
+       if(fNumberOfTowerClusters >= towerRecPoints->GetSize()) 
+         towerRecPoints->Expand(2*fNumberOfTowerClusters+1) ;
+       
+       towerRecPoints->AddAt(new  AliEMCALTowerRecPoint(""), fNumberOfTowerClusters) ;
+       clu = dynamic_cast<AliEMCALTowerRecPoint *>(towerRecPoints->At(fNumberOfTowerClusters)) ; 
+       fNumberOfTowerClusters++ ; 
+       clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower())) ; 
+       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;  
+       iDigitInCluster++ ; 
+       digitsC->Remove(digit) ; 
+       
+      } else { 
+       
+       // start a new Pre Shower cluster
+       if(fNumberOfPreShoClusters >= preshoRecPoints->GetSize()) 
+         preshoRecPoints->Expand(2*fNumberOfPreShoClusters+1);
+       
+       preshoRecPoints->AddAt(new AliEMCALTowerRecPoint(""), fNumberOfPreShoClusters) ;
+       
+       clu =  dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(fNumberOfPreShoClusters))  ;  
+       fNumberOfPreShoClusters++ ; 
+       clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->IsInPreShower() ) );     
+       clusterdigitslist[iDigitInCluster] = digit->GetIndexInList()  ; 
+       iDigitInCluster++ ; 
+       digitsC->Remove(digit) ; 
+       nextdigit.Reset() ;
+       
+       // Here we remove remaining Tower digits, which cannot make a cluster
+       
+       if( notremoved ) { 
+         while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
+           if( IsInTower(digit) )
+             digitsC->Remove(digit) ;
+           else 
+             break ; 
+         }
+         notremoved = kFALSE ;
+       }
+       
+      } // else        
+      
+      nextdigit.Reset() ;
+      
+      AliEMCALDigit * digitN ; 
+      index = 0 ;
+      while (index < iDigitInCluster){ // scan over digits already in cluster 
+       digit =  (AliEMCALDigit*)digits->At(clusterdigitslist[index])  ;      
+       index++ ; 
+        while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits 
+         Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!!
+         switch (ineb ) {
+          case 0 :   // not a neighbour
+           break ;
+         case 1 :   // are neighbours 
+           clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->IsInPreShower() ) ) ;
+           clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; 
+           iDigitInCluster++ ; 
+           digitsC->Remove(digitN) ;
+           break ;
+          case 2 :   // too far from each other
+           goto endofloop;   
+         } // switch
+         
+       } // while digitN
+       
+      endofloop: ;
+       nextdigit.Reset() ; 
+       
+      } // loop over cluster     
+      
+    } // energy theshold  
+    
+    
+  } // while digit
+  
+  delete digitsC ;
+  
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::MakeUnfolding()
+{
+  Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
+  
+//   // Unfolds clusters using the shape of an ElectroMagnetic shower
+//   // Performs unfolding of all EMC/CPV clusters
+
+//   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
+  
+//   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
+//   TObjArray * emcRecPoints = gime->TowerRecPoints() ; 
+//   TObjArray * cpvRecPoints = gime->PreShoRecPoints() ; 
+//   TClonesArray * digits = gime->Digits() ; 
+  
+//   // Unfold first EMC clusters 
+//   if(fNumberOfTowerClusters > 0){
+
+//     Int_t nModulesToUnfold = geom->GetNModules() ; 
+
+//     Int_t numberofNotUnfolded = fNumberOfTowerClusters ; 
+//     Int_t index ;   
+//     for(index = 0 ; index < numberofNotUnfolded ; index++){
+      
+//       AliEMCALTowerRecPoint * emcRecPoint = (AliEMCALTowerRecPoint *) emcRecPoints->At(index) ;
+//       if(emcRecPoint->GetEMCALMod()> nModulesToUnfold)
+//     break ;
+      
+//       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
+//       Int_t * maxAt = new Int_t[nMultipl] ;
+//       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
+//       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fTowerLocMaxCut,digits) ;
+      
+//       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
+//     UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+//     emcRecPoints->Remove(emcRecPoint); 
+//     emcRecPoints->Compress() ;
+//     index-- ;
+//     fNumberOfTowerClusters -- ;
+//     numberofNotUnfolded-- ;
+//       }
+      
+//       delete[] maxAt ; 
+//       delete[] maxAtEnergy ; 
+//     }
+//   } 
+//   // Unfolding of EMC clusters finished
+
+
+//   // Unfold now CPV clusters
+//   if(fNumberOfPreShoClusters > 0){
+    
+//     Int_t nModulesToUnfold = geom->GetNModules() ;
+
+//     Int_t numberofPreShoNotUnfolded = fNumberOfPreShoClusters ;     
+//     Int_t index ;   
+//     for(index = 0 ; index < numberofPreShoNotUnfolded ; index++){
+      
+//       AliEMCALRecPoint * recPoint = (AliEMCALRecPoint *) cpvRecPoints->At(index) ;
+
+//       if(recPoint->GetEMCALMod()> nModulesToUnfold)
+//     break ;
+      
+//       AliEMCALTowerRecPoint * emcRecPoint = (AliEMCALTowerRecPoint*) recPoint ; 
+      
+//       Int_t nMultipl = emcRecPoint->GetMultiplicity() ; 
+//       Int_t * maxAt = new Int_t[nMultipl] ;
+//       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
+//       Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fPreShoLocMaxCut,digits) ;
+      
+//       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0       
+//     UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+//     cpvRecPoints->Remove(emcRecPoint); 
+//     cpvRecPoints->Compress() ;
+//     index-- ;
+//     numberofPreShoNotUnfolded-- ;
+//     fNumberOfPreShoClusters-- ;
+//       }
+      
+//       delete[] maxAt ; 
+//       delete[] maxAtEnergy ; 
+//     } 
+//   }
+//   //Unfolding of PreSho clusters finished
+  
+}
+
+//____________________________________________________________________________
+Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t r)
+{ 
+  // Shape of the shower (see EMCAL TDR)
+  // If you change this function, change also the gradient evaluation in ChiSquare()
+
+  Double_t r4    = r*r*r*r ;
+  Double_t r295  = TMath::Power(r, 2.95) ;
+  Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
+  return shape ;
+}
+
+//____________________________________________________________________________
+void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower, 
+                                                Int_t nMax, 
+                                                int * maxAt, 
+                                                Float_t * maxAtEnergy)
+{
+  // Performs the unfolding of a cluster with nMax overlapping showers 
+  
+  Fatal("AliEMCALClusterizerv1::UnfoldCluster", "--> Unfolding not implemented") ;
+
+ //  AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
+//   const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
+//   const TClonesArray * digits = gime->Digits() ; 
+//   TObjArray * emcRecPoints = gime->TowerRecPoints() ; 
+//   TObjArray * cpvRecPoints = gime->PreShoRecPoints() ; 
+
+//   Int_t nPar = 3 * nMax ;
+//   Float_t * fitparameters = new Float_t[nPar] ;
+
+//   Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
+//   if( !rv ) {
+//     // Fit failed, return and remove cluster
+//     delete[] fitparameters ; 
+//     return ;
+//   }
+
+//   // create ufolded rec points and fill them with new energy lists
+//   // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
+//   // and later correct this number in acordance with actual energy deposition
+
+//   Int_t nDigits = iniTower->GetMultiplicity() ;  
+//   Float_t * efit = new Float_t[nDigits] ;
+//   Float_t xDigit=0.,zDigit=0.,distance=0. ;
+//   Float_t xpar=0.,zpar=0.,epar=0.  ;
+//   Int_t relid[4] ;
+//   AliEMCALDigit * digit = 0 ;
+//   Int_t * emcDigits = iniTower->GetDigitsList() ;
+
+//   Int_t iparam ;
+//   Int_t iDigit ;
+//   for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+//     digit = (AliEMCALDigit*) digits->At(emcDigits[iDigit] ) ;   
+//     geom->AbsToRelNumbering(digit->GetId(), relid) ;
+//     geom->RelPosInModule(relid, xDigit, zDigit) ;
+//     efit[iDigit] = 0;
+
+//     iparam = 0 ;    
+//     while(iparam < nPar ){
+//       xpar = fitparameters[iparam] ;
+//       zpar = fitparameters[iparam+1] ;
+//       epar = fitparameters[iparam+2] ;
+//       iparam += 3 ;
+//       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
+//       distance =  TMath::Sqrt(distance) ;
+//       efit[iDigit] += epar * ShowerShape(distance) ;
+//     }
+//   }
+  
+
+//   // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
+//   // so that energy deposited in each cell is distributed betwin new clusters proportionally
+//   // to its contribution to efit
+
+//   Float_t * emcEnergies = iniTower->GetEnergiesList() ;
+//   Float_t ratio ;
+
+//   iparam = 0 ;
+//   while(iparam < nPar ){
+//     xpar = fitparameters[iparam] ;
+//     zpar = fitparameters[iparam+1] ;
+//     epar = fitparameters[iparam+2] ;
+//     iparam += 3 ;    
+    
+//     AliEMCALTowerRecPoint * emcRP = 0 ;  
+
+//     if(iniTower->IsTower()){ //create new entries in fTowerRecPoints...
+      
+//       if(fNumberOfTowerClusters >= emcRecPoints->GetSize())
+//     emcRecPoints->Expand(2*fNumberOfTowerClusters) ;
+      
+//       (*emcRecPoints)[fNumberOfTowerClusters] = new AliEMCALTowerRecPoint("") ;
+//       emcRP = (AliEMCALTowerRecPoint *) emcRecPoints->At(fNumberOfTowerClusters);
+//       fNumberOfTowerClusters++ ;
+//     }
+//     else{//create new entries in fPreShoRecPoints
+//       if(fNumberOfPreShoClusters >= cpvRecPoints->GetSize())
+//     cpvRecPoints->Expand(2*fNumberOfPreShoClusters) ;
+      
+//       (*cpvRecPoints)[fNumberOfPreShoClusters] = new AliEMCALPreShoRecPoint("") ;
+//       emcRP = (AliEMCALTowerRecPoint *) cpvRecPoints->At(fNumberOfPreShoClusters);
+//       fNumberOfPreShoClusters++ ;
+//     }
+    
+//     Float_t eDigit ;
+//     for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
+//       digit = (AliEMCALDigit*) digits->At( emcDigits[iDigit] ) ; 
+//       geom->AbsToRelNumbering(digit->GetId(), relid) ;
+//       geom->RelPosInModule(relid, xDigit, zDigit) ;
+//       distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
+//       distance =  TMath::Sqrt(distance) ;
+//       ratio = epar * ShowerShape(distance) / efit[iDigit] ; 
+//       eDigit = emcEnergies[iDigit] * ratio ;
+//       emcRP->AddDigit( *digit, eDigit ) ;
+//     }       
+//   }
+//   delete[] fitparameters ; 
+//   delete[] efit ; 
+  
+}
+
+//_____________________________________________________________________________
+void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
+{
+  // Calculates the Chi square for the cluster unfolding minimization
+  // Number of parameters, Gradient, Chi squared, parameters, what to do
+
+  abort() ; 
+ //  Fatal("AliEMCALClusterizerv1::UnfoldingChiSquare","-->Unfolding not implemented") ;
+
+//   TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
+
+//   AliEMCALTowerRecPoint * emcRP = (AliEMCALTowerRecPoint*) toMinuit->At(0)  ;
+//   TClonesArray * digits = (TClonesArray*)toMinuit->At(1)  ;
+
+
+  
+//   //  AliEMCALTowerRecPoint * emcRP = (AliEMCALTowerRecPoint *) gMinuit->GetObjectFit() ; // TowerRecPoint to fit
+
+//   Int_t * emcDigits     = emcRP->GetDigitsList() ;
+
+//   Int_t nOdigits = emcRP->GetDigitsMultiplicity() ; 
+
+//   Float_t * emcEnergies = emcRP->GetEnergiesList() ;
+
+//   const AliEMCALGeometry * geom = AliEMCALGetter::GetInstance()->EMCALGeometry() ; 
+//   fret = 0. ;     
+//   Int_t iparam ;
+
+//   if(iflag == 2)
+//     for(iparam = 0 ; iparam < nPar ; iparam++)    
+//       Grad[iparam] = 0 ; // Will evaluate gradient
+  
+//   Double_t efit ;    
+
+//   AliEMCALDigit * digit ;
+//   Int_t iDigit ;
+
+//   for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
+
+//     digit = (AliEMCALDigit*) digits->At( emcDigits[iDigit] ) ; 
+
+//     Int_t relid[4] ;
+//     Float_t xDigit ;
+//     Float_t zDigit ;
+
+//     geom->AbsToRelNumbering(digit->GetId(), relid) ;
+
+//     geom->RelPosInModule(relid, xDigit, zDigit) ;
+
+//      if(iflag == 2){  // calculate gradient
+//        Int_t iParam = 0 ;
+//        efit = 0 ;
+//        while(iParam < nPar ){
+//      Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
+//      iParam++ ; 
+//      distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ; 
+//      distance = TMath::Sqrt( distance ) ; 
+//      iParam++ ;      
+//      efit += x[iParam] * ShowerShape(distance) ;
+//      iParam++ ;
+//        }
+//        Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E) 
+//        iParam = 0 ;
+//        while(iParam < nPar ){
+//      Double_t xpar = x[iParam] ;
+//      Double_t zpar = x[iParam+1] ;
+//      Double_t epar = x[iParam+2] ;
+//      Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
+//      Double_t shape = sum * ShowerShape(dr) ;
+//      Double_t r4 = dr*dr*dr*dr ;
+//      Double_t r295 = TMath::Power(dr,2.95) ;
+//      Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
+//                                      0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
+        
+//      Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x    
+//      iParam++ ; 
+//      Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z         
+//      iParam++ ; 
+//      Grad[iParam] += shape ;                                  // Derivative over energy             
+//      iParam++ ; 
+//        }
+//      }
+//      efit = 0;
+//      iparam = 0 ;
+
+//      while(iparam < nPar ){
+//        Double_t xpar = x[iparam] ;
+//        Double_t zpar = x[iparam+1] ;
+//        Double_t epar = x[iparam+2] ;
+//        iparam += 3 ;
+//        Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar)  ;
+//        distance =  TMath::Sqrt(distance) ;
+//        efit += epar * ShowerShape(distance) ;
+//      }
+
+//      fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ; 
+//      // Here we assume, that sigma = sqrt(E)
+//   }
+
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::Print(Option_t * option)const
+{
+  // Print clusterizer parameters
+
+  if( strcmp(GetName(), "") !=0 ){
+    
+    // Print parameters
+    TString taskName(GetName()) ; 
+    taskName.ReplaceAll(Version(), "") ;
+
+    cout << "---------------"<< taskName.Data() << " " << GetTitle()<< "-----------" << endl 
+        << "Clusterizing digits from the file: " << fHeaderFileName.Data() << endl 
+        << "                           Branch: " << fDigitsBranchTitle.Data() << endl 
+        << endl 
+        << "                       EMC Clustering threshold = " << fTowerClusteringThreshold << endl
+        << "                       EMC Local Maximum cut    = " << fTowerLocMaxCut << endl
+        << "                       EMC Logarothmic weight   = " << fW0 << endl
+        << endl
+        << "                       CPV Clustering threshold = " << fPreShoClusteringThreshold << endl
+        << "                       CPV Local Maximum cut    = " << fPreShoLocMaxCut << endl
+       << "                       CPV Logarothmic weight   = " << fW0CPV << endl
+        << endl ;
+    if(fToUnfold)
+      cout << " Unfolding on " << endl ;
+    else
+      cout << " Unfolding off " << endl ;
+    
+    cout << "------------------------------------------------------------------" <<endl ;
+  }
+  else
+    cout << " AliEMCALClusterizerv1 not initialized " << endl ;
+}
+//____________________________________________________________________________
+void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
+{
+  // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
+
+  TObjArray * towerRecPoints = AliEMCALGetter::GetInstance()->TowerRecPoints() ; 
+  TObjArray * preshoRecPoints = AliEMCALGetter::GetInstance()->PreShowerRecPoints() ; 
+
+  cout << "AliEMCALClusterizerv1: : event "<<gAlice->GetEvNumber() << endl ;
+  cout << "       Found "<< towerRecPoints->GetEntriesFast() << " TOWER Rec Points and " 
+          << preshoRecPoints->GetEntriesFast() << " PRE SHOWER RecPoints" << endl ;
+
+  fRecPointsInRun +=  towerRecPoints->GetEntriesFast() ; 
+  fRecPointsInRun +=  preshoRecPoints->GetEntriesFast() ; 
+
+  if(strstr(option,"all")) {
+
+    cout << "Tower clusters " << endl ;
+    cout << " Index  Ene(MeV)   Multi  Module     phi     r  theta    Lambda 1   Lambda 2  # of prim  Primaries list "      <<  endl;      
+    
+    Int_t index ;
+    for (index = 0 ; index < towerRecPoints->GetEntries() ; index++) {
+      AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(towerRecPoints->At(index)) ; 
+      TVector3  globalpos;  
+      rp->GetGlobalPosition(globalpos);
+      Float_t lambda[2]; 
+      rp->GetElipsAxis(lambda);
+      Int_t * primaries; 
+      Int_t nprimaries;
+      primaries = rp->GetPrimaries(nprimaries);
+
+      cout << setw(4) << rp->GetIndexInList() << "   " 
+          << setw(7) << setprecision(3) << rp->GetEnergy() << "      "
+          << setw(3) << rp->GetMultiplicity() << "      " 
+          << setw(1) << rp->GetEMCALArm() << "     " 
+          << setw(5) << setprecision(4) << globalpos.X() << "  " 
+          << setw(5) << setprecision(4) << globalpos.Y() << "  " 
+          << setw(5) << setprecision(4) << globalpos.Z() << "     "
+          << setw(4) << setprecision(2) << lambda[0]  << "  "
+          << setw(4) << setprecision(2) << lambda[1]  << "  "
+          << setw(2) << nprimaries << "  " ;
+     
+      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
+       cout << setw(4) <<   primaries[iprimary] << "  "  ;
+      cout << endl ;    
+    }
+
+    //Now plot Pre shower recPoints
+
+    cout << "-----------------------------------------------------------------------"<<endl ;
+
+    cout << "PreShower clusters " << endl ;
+    cout << " Index  Ene(MeV)   Multi  Module     phi     r  theta    Lambda 1   Lambda 2  # of prim  Primaries list "      <<  endl;      
+    
+    for (index = 0 ; index < preshoRecPoints->GetEntries() ; index++) {
+      AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(preshoRecPoints->At(index)) ; 
+      TVector3  globalpos;  
+      rp->GetGlobalPosition(globalpos);
+      Float_t lambda[2]; 
+      rp->GetElipsAxis(lambda);
+      Int_t * primaries; 
+      Int_t nprimaries;
+      primaries = rp->GetPrimaries(nprimaries);
+
+      cout << setw(4) << rp->GetIndexInList() << "   " 
+          << setw(7) << setprecision(3) << rp->GetEnergy() << "      "
+          << setw(3) << rp->GetMultiplicity() << "      " 
+          << setw(1) << rp->GetEMCALArm() << "     " 
+          << setw(5) << setprecision(4) << globalpos.X() << "  " 
+          << setw(5) << setprecision(4) << globalpos.Y() << "  " 
+          << setw(5) << setprecision(4) << globalpos.Z() << "     "
+          << setw(4) << setprecision(2) << lambda[0]  << "  "
+          << setw(4) << setprecision(2) << lambda[1]  << "  "
+          << setw(2) << nprimaries << "  " ;
+     
+      for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
+       cout << setw(4) <<   primaries[iprimary] << "  "  ;
+      cout << endl ;    
+    }
+
+    cout << "-----------------------------------------------------------------------"<<endl ;
+  }
+}
+
diff --git a/EMCAL/AliEMCALClusterizerv1.h b/EMCAL/AliEMCALClusterizerv1.h
new file mode 100644 (file)
index 0000000..1e15dbd
--- /dev/null
@@ -0,0 +1,129 @@
+#ifndef ALIEMCALCLUSTERIZERV1_H
+#define ALIEMCALCLUSTERIZERV1_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//  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
+//
+//*-- Author: Yves Schutz (SUBATECH)
+
+// --- ROOT system ---
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+
+#include "AliEMCALClusterizer.h"
+class AliEMCALTowerRecPoint ; 
+class AliEMCALDigit ;
+class AliEMCALDigitizer ;
+class AliEMCALGeometry ;
+
+
+class AliEMCALClusterizerv1 : public AliEMCALClusterizer {
+  
+public:
+  
+  AliEMCALClusterizerv1() ;         
+  AliEMCALClusterizerv1(const char * headerFile, const char * name = "Default");
+  virtual ~AliEMCALClusterizerv1()  ;
+  
+  virtual Int_t   AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const ; 
+                               // Checks if digits are in neighbour cells 
+
+  virtual Float_t Calibrate(Int_t amp, Bool_t inpresho)const ;  // Tranforms Amp to energy 
+
+  virtual void    GetNumberOfClustersFound(int * numb )const{  numb[0] = fNumberOfTowerClusters ; 
+                                                               numb[1] = fNumberOfPreShoClusters ; }
+
+  virtual Float_t GetEmcClusteringThreshold()const{ return fTowerClusteringThreshold;}
+  virtual Float_t GetEmcLocalMaxCut()const        { return fTowerLocMaxCut;} 
+  virtual Float_t GetEmcLogWeight()const          { return fW0;}  
+  virtual Float_t GetTimeGate() const             { return fTimeGate ; }
+  virtual Float_t GetCpvClusteringThreshold()const{ return fPreShoClusteringThreshold;  } 
+  virtual Float_t GetCpvLocalMaxCut()const        { return fPreShoLocMaxCut;} 
+  virtual Float_t GetCpvLogWeight()const          { return fW0CPV;}  
+  virtual char *  GetRecPointsBranch() const      { return (char*) fRecPointsBranchTitle.Data() ;}
+  virtual const Int_t GetRecPointsInRun() const  {return fRecPointsInRun ;} 
+  virtual char *  GetDigitsBranch() const         { return (char*) fDigitsBranchTitle.Data() ;}
+
+  void    Exec(Option_t *option);                // Does the job
+
+  virtual void Print(Option_t * option)const ;
+
+  virtual void SetTowerClusteringThreshold(Float_t cluth)  { fTowerClusteringThreshold = cluth ; }
+  virtual void SetTowerLocalMaxCut(Float_t cut)            { fTowerLocMaxCut = cut ; }
+  virtual void SetTowerLogWeight(Float_t w)                { fW0 = w ; }
+  virtual void SetTimeGate(Float_t gate)                   { fTimeGate = gate ;}
+  virtual void SetPreShoClusteringThreshold(Float_t cluth) { fPreShoClusteringThreshold = cluth ; }
+  virtual void SetPreShoLocalMaxCut(Float_t cut)           { fPreShoLocMaxCut = cut ; }
+  virtual void SetPreShoLogWeight(Float_t w)               { fW0CPV = w ; }
+  virtual void SetDigitsBranch(const char * title) { fDigitsBranchTitle = title  ;}
+  virtual void SetRecPointsBranch(const char *title){fRecPointsBranchTitle = title; }
+  virtual void SetUnfolding(Bool_t toUnfold = kTRUE ) {fToUnfold = toUnfold ;}  
+  static Double_t ShowerShape(Double_t r) ; // 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" ; }  
+
+protected:
+
+  void           WriteRecPoints(Int_t event) ;
+  virtual void   MakeClusters( ) ;            
+  virtual Bool_t IsInTower (AliEMCALDigit * digit)const ;     // Tells if id digit is in Tower
+  virtual Bool_t IsInPreShower (AliEMCALDigit * digit)const ;     // Tells if id digit is in PreShower
+
+  
+private:
+
+  void    GetCalibrationParameters(void) ;
+  
+  Bool_t  FindFit(AliEMCALTowerRecPoint * emcRP, int * MaxAt, Float_t * maxAtEnergy, 
+                 Int_t NPar, Float_t * FitParametres) const; //Used in UnfoldClusters, calls TMinuit
+  void Init() ;
+
+  virtual void   MakeUnfolding() ;
+  void           UnfoldCluster(AliEMCALTowerRecPoint * iniEmc,Int_t Nmax, 
+                      int * maxAt,Float_t * maxAtEnergy ) ; //Unfolds cluster using TMinuit package
+  void           PrintRecPoints(Option_t * option) ;
+
+private:
+
+  TString fHeaderFileName ;          // name of the file which contains gAlice, Tree headers etc.
+  TString fDigitsBranchTitle ;       // name of the file, where digits branch is stored
+  TString fRecPointsBranchTitle ;    // name of the file, where RecPoints branchs are stored
+
+  Int_t   fNTowers ;                 // number of Towers in EMCAL
+
+  Bool_t  fToUnfold ;                // To perform unfolding 
+
+  Int_t   fNumberOfTowerClusters ;     // number of Tower clusters found 
+  Int_t   fNumberOfPreShoClusters ;    // number of PreShower clusters found
+  Float_t fADCchannelTower ;           // width of one ADC channel for Tower (GeV)
+  Float_t fADCpedestalTower ;          // pedestal of ADC for Tower (GeV) 
+  Float_t fADCchannelPreSho ;          // width of one ADC channel for Pre Shower (GeV)
+  Float_t fADCpedestalPreSho ;         // pedestal of ADC for PreShower (GeV)
+
+  Float_t fTowerClusteringThreshold ;  // minimum energy to include a EMC digit in a cluster
+  Float_t fPreShoClusteringThreshold ; // minimum energy to include a CPV digit in a cluster
+  Float_t fTowerLocMaxCut ;            // minimum energy difference to distinguish local maxima in a cluster
+  Float_t fW0 ;                        // logarithmic weight for the cluster center of gravity calculation
+  Float_t fPreShoLocMaxCut ;           //  minimum energy difference to distinguish local maxima in a CPV cluster
+  Float_t fW0CPV ;                   // logarithmic weight for the CPV cluster center of gravity calculation
+  Int_t fRecPointsInRun ;            //! Total number of recpoints in one run
+  Float_t fTimeGate ;                // Maximum time difference between the digits in ont EMC cluster
+    
+  ClassDef(AliEMCALClusterizerv1,1)   // Clusterizer implementation version 1
+
+};
+
+#endif // AliEMCALCLUSTERIZERV1_H
index 27f4400..2c115d9 100644 (file)
 #pragma link C++ class AliEMCALHit+;
 #pragma link C++ class AliEMCALDigit+;
 #pragma link C++ class AliEMCALTick+;
-//#pragma link C++ class AliEMCALRecPoint+;
-//#pragma link C++ class AliEMCALEmcRecPoint+;
+#pragma link C++ class AliEMCALRecPoint+;
+#pragma link C++ class AliEMCALTowerRecPoint+;
 #pragma link C++ class AliEMCALSDigitizer+;
 #pragma link C++ class AliEMCALDigitizer+;
-//#pragma link C++ class AliEMCALClusterizer+;
+#pragma link C++ class AliEMCALClusterizer+;
+#pragma link C++ class AliEMCALClusterizerv1+;
 #pragma link C++ class AliEMCALGetter+;
 #pragma link C++ class AliEMCALJetFinder+;
 #pragma link C++ class AliEMCALJet+;
index c1d7d95..2e25ff0 100644 (file)
@@ -15,8 +15,8 @@ SRCS          =  AliEMCAL.cxx  AliEMCALGeometry.cxx AliEMCALv0.cxx \
                 AliEMCALGetter.cxx \
                 AliEMCALJetFinder.cxx AliEMCALJet.cxx AliEMCALFast.cxx \
                 AliEMCALGetter.cxx AliEMCALTick.cxx \
-                AliEMCALHadronCorrection.cxx AliEMCALHadronCorrectionv0.cxx 
-#                AliEMCALRecPoint.cxx AliEMCALEmcRecPoint.cxx AliEMCALClusterizer.cxx
+                AliEMCALRecPoint.cxx AliEMCALTowerRecPoint.cxx \
+                AliEMCALClusterizer.cxx AliEMCALClusterizerv1.cxx
 
 # FORTRAN sources