[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerNxN.cxx

/**************************************************************************
 * 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"
#include "AliEMCALUnfolding.h"

ClassImp(AliEMCALClusterizerNxN)

//____________________________________________________________________________
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){
    fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
    fClusterUnfolding->MakeUnfolding();
  }
  
  //Evaluate position, dispersion and other RecPoint properties for EC section 
  Int_t index ;
  for(index = 0; index < fRecPoints->GetEntries(); index++) 
  { 
    AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
    if(rp){
      rp->EvalAll(fECAW0,fDigitsArr,fgkIsInputCalibrated) ;
      AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
      //For each rec.point set the distance to the nearest bad crystal
      rp->EvalDistanceToBadChannels(fCaloPed);
    }
  }
  
  fRecPoints->Sort() ;
  
  for(index = 0; index < fRecPoints->GetEntries(); index++) 
  {
    AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
    if(rp){
      rp->SetIndexInList(index) ;
      rp->Print();
    }
    else AliFatal("RecPoint NULL!!");
  }
  
  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)) ; 
      if(recPoint){
        fNumberOfECAClusters++ ;       
        recPoint->SetClusterType(AliVCluster::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); 		  
        }
      }// recpoint
    }
    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())); 
}
Commit	Line	Data
ee08edde	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
	17	// Algorithm:
	18	// 1. peek the most energetic cell
	19	// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
	20	// 3. remove the cells contributing to the cluster
	21	// 4. start from 1 for the remaining clusters
	22	// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
	23	// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy
	24	// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
	25	// Use Case:
	26	// root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")
	27	// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
	28	// //reads gAlice from header file "..."
	29	// root [1] cl->ExecuteTask()
	30	// //finds RecPoints in all events stored in galice.root
	31	// root [2] cl->SetDigitsBranch("digits2")
	32	// //sets another title for Digitis (input) branch
	33	// root [3] cl->SetRecPointsBranch("recp2")
	34	// //sets another title four output branches
	35	// root [4] cl->SetTowerLocalMaxCut(0.03)
	36	// //set clusterization parameters
	37	// root [5] cl->ExecuteTask("deb all time")
	38	// //once more finds RecPoints options are
	39	// // deb - print number of found rec points
	40	// // deb all - print number of found RecPoints and some their characteristics
	41	// // time - print benchmarking results
	42
	43	// --- ROOT system ---
	44	#include <TMath.h>
	45	#include <TMinuit.h>
	46	#include <TTree.h>
	47	#include <TBenchmark.h>
	48	#include <TBrowser.h>
	49	#include <TROOT.h>
	50	#include <TClonesArray.h>
	51
	52	// --- Standard library ---
	53	#include <cassert>
	54
	55	// --- AliRoot header files ---
	56	#include "AliLog.h"
	57	#include "AliEMCALClusterizerNxN.h"
	58	#include "AliEMCALRecPoint.h"
	59	#include "AliEMCALDigit.h"
	60	#include "AliEMCALGeometry.h"
	61	#include "AliCaloCalibPedestal.h"
	62	#include "AliEMCALCalibData.h"
	63	#include "AliESDCaloCluster.h"
65bec413	64	#include "AliEMCALUnfolding.h"
ee08edde	65
	66	ClassImp(AliEMCALClusterizerNxN)
	67
ee08edde	68	//____________________________________________________________________________
	69	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
	70	: AliEMCALClusterizer()
	71	{
	72	// ctor with the indication of the file where header Tree and digits Tree are stored
	73	}
	74
	75	//____________________________________________________________________________
	76	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
	77	: AliEMCALClusterizer(geometry)
	78	{
	79	// ctor with the indication of the file where header Tree and digits Tree are stored
	80	// use this contructor to avoid usage of Init() which uses runloader
	81	// change needed by HLT - MP
	82
	83	}
	84
	85	//____________________________________________________________________________
	86	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
	87	: AliEMCALClusterizer(geometry, calib, caloped)
	88	{
	89	// ctor, geometry and calibration are initialized elsewhere.
	90
	91	}
	92
	93
	94	//____________________________________________________________________________
	95	AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
	96	{
	97	// dtor
	98	}
	99
	100
	101	//____________________________________________________________________________
	102	void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
	103	{
	104	// Steering method to perform clusterization for the current event
	105	// in AliEMCALLoader
51e4198e	106
ee08edde	107	if(strstr(option,"tim"))
	108	gBenchmark->Start("EMCALClusterizer");
	109
	110	if(strstr(option,"print"))
	111	Print("") ;
51e4198e	112
ee08edde	113	//Get calibration parameters from file or digitizer default values.
ee08edde	114	GetCalibrationParameters() ;
51e4198e	115
ee08edde	116	//Get dead channel map from file or digitizer default values.
	117	GetCaloCalibPedestal() ;
	118
	119	fNumberOfECAClusters = 0;
51e4198e	120
ee08edde	121	MakeClusters() ; //only the real clusters
51e4198e	122
65bec413	123	if(fToUnfold){
	124	fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
	125	fClusterUnfolding->MakeUnfolding();
	126	}
51e4198e	127
65bec413	128	//Evaluate position, dispersion and other RecPoint properties for EC section
ee08edde	129	Int_t index ;
ee08edde	130	for(index = 0; index < fRecPoints->GetEntries(); index++)
51e4198e	131	{
	132	AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
	133	if(rp){
0d0d6b98	134	rp->EvalAll(fECAW0,fDigitsArr,fgkIsInputCalibrated) ;
51e4198e	135	AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
ee08edde	136	//For each rec.point set the distance to the nearest bad crystal
51e4198e	137	rp->EvalDistanceToBadChannels(fCaloPed);
ee08edde	138	}
51e4198e	139	}
ee08edde	140
ee08edde	141	fRecPoints->Sort() ;
51e4198e	142
ee08edde	143	for(index = 0; index < fRecPoints->GetEntries(); index++)
51e4198e	144	{
7e1d9a9b	145	AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
	146	if(rp){
	147	rp->SetIndexInList(index) ;
	148	rp->Print();
	149	}
	150	else AliFatal("RecPoint NULL!!");
51e4198e	151	}
ee08edde	152
	153	fTreeR->Fill();
	154
	155	if(strstr(option,"deb") \|\| strstr(option,"all"))
	156	PrintRecPoints(option) ;
51e4198e	157
ee08edde	158	AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
51e4198e	159
ee08edde	160	fRecPoints->Delete();
51e4198e	161
ee08edde	162	if(strstr(option,"tim")){
	163	gBenchmark->Stop("EMCALClusterizer");
	164	printf("Exec took %f seconds for Clusterizing",
51e4198e	165	gBenchmark->GetCpuTime("EMCALClusterizer"));
ee08edde	166	}
	167	}
	168
	169	//____________________________________________________________________________
	170	Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
	171	{
c526514b	172	// Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
	173	// = 1 are neighbour
	174	// = 2 is in different SM; continue searching
	175	// In case it is in different SM, but same phi rack, check if neigbours at eta=0
	176	// neighbours are defined as digits having at least a common side
	177	// The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
	178	// which is compared to a digit (d2) not yet in a cluster
	179
	180	static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
	181	static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
	182	static Int_t rowdiff=0, coldiff=0;
	183
	184	shared = kFALSE;
	185
	186	fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
	187	fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
	188	fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
	189	fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
	190
	191	//If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
	192	if(nSupMod1 != nSupMod2 )
	193	{
	194	//Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
	195	Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
	196	Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
	197
	198	if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
	199
	200	// In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
	201	// C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
	202	if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
	203	else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
	204
	205	shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
	206
	207	}//Different SM, same phi
	208
	209	rowdiff = TMath::Abs(iphi1 - iphi2);
	210	coldiff = TMath::Abs(ieta1 - ieta2) ;
	211
	212	// neighbours +-1 in col and row
	213	if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
	214	{
	215
	216	AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
	217	d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
	218
	219	return 1;
	220	}//Neighbours
	221	else
	222	{
	223	AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
	224	d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
	225	shared = kFALSE;
	226	return 2 ;
	227	}//Not neighbours
ee08edde	228	}
	229
	230	//____________________________________________________________________________
	231	void AliEMCALClusterizerNxN::MakeClusters()
	232	{
	233	// Steering method to construct the clusters stored in a list of Reconstructed Points
	234	// A cluster is defined as a list of neighbour digits
	235	// Mar 03, 2007 by PAI
51e4198e	236
ee08edde	237	if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
51e4198e	238
ee08edde	239	fRecPoints->Clear();
51e4198e	240
ee08edde	241	// Set up TObjArray with pointers to digits to work on
	242	//TObjArray *digitsC = new TObjArray();
	243	TObjArray digitsC;
	244	TIter nextdigit(fDigitsArr);
	245	AliEMCALDigit *digit = 0;
	246	while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
	247	digitsC.AddLast(digit);
	248	}
51e4198e	249
ee08edde	250	TIter nextdigitC(&digitsC);
51e4198e	251
ee08edde	252	AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f)\n",
51e4198e	253	fDigitsArr->GetEntries(),fMinECut));
51e4198e	254
ee08edde	255	Bool_t bDone = kFALSE;
ee08edde	256	while ( bDone != kTRUE )
51e4198e	257	{
	258	//first sort the digits:
	259	Int_t iMaxEnergyDigit = -1;
	260	Float_t dMaxEnergyDigit = -1;
	261	AliEMCALDigit *pMaxEnergyDigit = 0;
	262	nextdigitC.Reset();
	263	while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
	264	{ // scan over the list of digitsC
	265	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	266	//AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
	267
	268	//if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
	269	if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
ee08edde	270	{
ee08edde	271	if (dEnergyCalibrated > dMaxEnergyDigit)
51e4198e	272	{
	273	dMaxEnergyDigit = dEnergyCalibrated;
	274	iMaxEnergyDigit = digit->GetId();
	275	pMaxEnergyDigit = digit;
	276	}
ee08edde	277	}
51e4198e	278	}
	279
	280	if (iMaxEnergyDigit < 0 \|\| digitsC.GetEntries() <= 0)
	281	{
	282	bDone = kTRUE;
	283	continue;
	284	}
	285
	286	AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
	287	AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));
	288
	289	// keep the candidate digits in a list
	290	TList clusterDigitList;
	291	clusterDigitList.SetOwner(kFALSE);
	292	clusterDigitList.AddLast(pMaxEnergyDigit);
	293
	294	Double_t clusterCandidateEnergy = dMaxEnergyDigit;
	295
	296	// now loop over the resto of the digits and cluster into NxN cluster
	297	// we do not actually cluster yet: we keep them in the list clusterDigitList
	298	nextdigitC.Reset();
	299	while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
	300	{ // scan over the list of digitsC
	301	if (digit == pMaxEnergyDigit) continue;
	302	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	303	AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
	304	//if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
	305	if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 )
ee08edde	306	{
	307	Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
	308	if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
	309	Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
	310	if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
51e4198e	311	{
	312	clusterDigitList.AddLast(digit) ;
	313	clusterCandidateEnergy += dEnergyCalibrated;
	314	}
ee08edde	315	}
51e4198e	316	}// loop over the next digits
	317
	318	// start a cluster here only if a cluster energy is larger than clustering threshold
	319	//if (clusterCandidateEnergy > 0.1)
	320	AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
	321	if (clusterCandidateEnergy > fECAClusteringThreshold)
	322	{
	323	if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
ee08edde	324
51e4198e	325	AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
	326	fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
	327	recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
	328	if(recPoint){
	329	fNumberOfECAClusters++ ;
	330	recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
	331
	332	AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
	333	for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
	334	{
	335
	336	digit = (AliEMCALDigit*)clusterDigitList.At(idig);
	337	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	338	AliDebug(5, Form(" Adding digit %d", digit->GetId()));
	339	// note: this way the sharing info is lost!
	340	recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
	341	digitsC.Remove(digit);
	342	}
	343	}// recpoint
	344	}
	345	else
	346	{
	347	// we do not want to start clustering in the same spot!
	348	// but in this case we may NOT reuse this seed for another cluster!
	349	// need a better bookeeping?
	350	digitsC.Remove(pMaxEnergyDigit);
	351	}
	352
	353	AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
	354	} // while ! done
	355
ee08edde	356	//delete digitsC ; //nope we use an object
	357
	358	AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
	359	}
	360