[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)

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){
    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) ;
      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())); 
}

//___________________________________________________________________
void   AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
{
  //
  // input is calibrated - the case when we run already on ESD
  //
  AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;
}
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
	68	Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;
	69
	70	//____________________________________________________________________________
	71	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
	72	: AliEMCALClusterizer()
	73	{
	74	// ctor with the indication of the file where header Tree and digits Tree are stored
	75	}
	76
	77	//____________________________________________________________________________
	78	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
	79	: AliEMCALClusterizer(geometry)
	80	{
	81	// ctor with the indication of the file where header Tree and digits Tree are stored
	82	// use this contructor to avoid usage of Init() which uses runloader
	83	// change needed by HLT - MP
	84
	85	}
	86
	87	//____________________________________________________________________________
	88	AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
	89	: AliEMCALClusterizer(geometry, calib, caloped)
	90	{
	91	// ctor, geometry and calibration are initialized elsewhere.
	92
	93	}
	94
	95
	96	//____________________________________________________________________________
	97	AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
	98	{
	99	// dtor
	100	}
	101
	102
	103	//____________________________________________________________________________
	104	void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
	105	{
	106	// Steering method to perform clusterization for the current event
	107	// in AliEMCALLoader
51e4198e	108
ee08edde	109	if(strstr(option,"tim"))
	110	gBenchmark->Start("EMCALClusterizer");
	111
	112	if(strstr(option,"print"))
	113	Print("") ;
51e4198e	114
ee08edde	115	//Get calibration parameters from file or digitizer default values.
ee08edde	116	GetCalibrationParameters() ;
51e4198e	117
ee08edde	118	//Get dead channel map from file or digitizer default values.
	119	GetCaloCalibPedestal() ;
	120
	121	fNumberOfECAClusters = 0;
51e4198e	122
ee08edde	123	MakeClusters() ; //only the real clusters
51e4198e	124
65bec413	125	if(fToUnfold){
	126	fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
	127	fClusterUnfolding->MakeUnfolding();
	128	}
51e4198e	129
65bec413	130	//Evaluate position, dispersion and other RecPoint properties for EC section
ee08edde	131	Int_t index ;
ee08edde	132	for(index = 0; index < fRecPoints->GetEntries(); index++)
51e4198e	133	{
	134	AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
	135	if(rp){
	136	rp->EvalAll(fECAW0,fDigitsArr) ;
	137	AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
ee08edde	138	//For each rec.point set the distance to the nearest bad crystal
51e4198e	139	rp->EvalDistanceToBadChannels(fCaloPed);
ee08edde	140	}
51e4198e	141	}
ee08edde	142
ee08edde	143	fRecPoints->Sort() ;
51e4198e	144
ee08edde	145	for(index = 0; index < fRecPoints->GetEntries(); index++)
51e4198e	146	{
7e1d9a9b	147	AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
	148	if(rp){
	149	rp->SetIndexInList(index) ;
	150	rp->Print();
	151	}
	152	else AliFatal("RecPoint NULL!!");
51e4198e	153	}
ee08edde	154
	155	fTreeR->Fill();
	156
	157	if(strstr(option,"deb") \|\| strstr(option,"all"))
	158	PrintRecPoints(option) ;
51e4198e	159
ee08edde	160	AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
51e4198e	161
ee08edde	162	fRecPoints->Delete();
51e4198e	163
ee08edde	164	if(strstr(option,"tim")){
	165	gBenchmark->Stop("EMCALClusterizer");
	166	printf("Exec took %f seconds for Clusterizing",
51e4198e	167	gBenchmark->GetCpuTime("EMCALClusterizer"));
ee08edde	168	}
	169	}
	170
	171	//____________________________________________________________________________
	172	Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
	173	{
c526514b	174	// Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
	175	// = 1 are neighbour
	176	// = 2 is in different SM; continue searching
	177	// In case it is in different SM, but same phi rack, check if neigbours at eta=0
	178	// neighbours are defined as digits having at least a common side
	179	// The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
	180	// which is compared to a digit (d2) not yet in a cluster
	181
	182	static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
	183	static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
	184	static Int_t rowdiff=0, coldiff=0;
	185
	186	shared = kFALSE;
	187
	188	fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
	189	fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
	190	fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
	191	fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
	192
	193	//If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
	194	if(nSupMod1 != nSupMod2 )
	195	{
	196	//Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
	197	Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
	198	Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
	199
	200	if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
	201
	202	// In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
	203	// C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
	204	if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
	205	else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
	206
	207	shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
	208
	209	}//Different SM, same phi
	210
	211	rowdiff = TMath::Abs(iphi1 - iphi2);
	212	coldiff = TMath::Abs(ieta1 - ieta2) ;
	213
	214	// neighbours +-1 in col and row
	215	if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
	216	{
	217
	218	AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
	219	d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
	220
	221	return 1;
	222	}//Neighbours
	223	else
	224	{
	225	AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
	226	d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
	227	shared = kFALSE;
	228	return 2 ;
	229	}//Not neighbours
ee08edde	230	}
	231
	232	//____________________________________________________________________________
	233	void AliEMCALClusterizerNxN::MakeClusters()
	234	{
	235	// Steering method to construct the clusters stored in a list of Reconstructed Points
	236	// A cluster is defined as a list of neighbour digits
	237	// Mar 03, 2007 by PAI
51e4198e	238
ee08edde	239	if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
51e4198e	240
ee08edde	241	fRecPoints->Clear();
51e4198e	242
ee08edde	243	// Set up TObjArray with pointers to digits to work on
	244	//TObjArray *digitsC = new TObjArray();
	245	TObjArray digitsC;
	246	TIter nextdigit(fDigitsArr);
	247	AliEMCALDigit *digit = 0;
	248	while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
	249	digitsC.AddLast(digit);
	250	}
51e4198e	251
ee08edde	252	TIter nextdigitC(&digitsC);
51e4198e	253
ee08edde	254	AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f)\n",
51e4198e	255	fDigitsArr->GetEntries(),fMinECut));
51e4198e	256
ee08edde	257	Bool_t bDone = kFALSE;
ee08edde	258	while ( bDone != kTRUE )
51e4198e	259	{
	260	//first sort the digits:
	261	Int_t iMaxEnergyDigit = -1;
	262	Float_t dMaxEnergyDigit = -1;
	263	AliEMCALDigit *pMaxEnergyDigit = 0;
	264	nextdigitC.Reset();
	265	while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
	266	{ // scan over the list of digitsC
	267	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	268	//AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
	269
	270	//if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
	271	if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
ee08edde	272	{
ee08edde	273	if (dEnergyCalibrated > dMaxEnergyDigit)
51e4198e	274	{
	275	dMaxEnergyDigit = dEnergyCalibrated;
	276	iMaxEnergyDigit = digit->GetId();
	277	pMaxEnergyDigit = digit;
	278	}
ee08edde	279	}
51e4198e	280	}
	281
	282	if (iMaxEnergyDigit < 0 \|\| digitsC.GetEntries() <= 0)
	283	{
	284	bDone = kTRUE;
	285	continue;
	286	}
	287
	288	AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
	289	AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));
	290
	291	// keep the candidate digits in a list
	292	TList clusterDigitList;
	293	clusterDigitList.SetOwner(kFALSE);
	294	clusterDigitList.AddLast(pMaxEnergyDigit);
	295
	296	Double_t clusterCandidateEnergy = dMaxEnergyDigit;
	297
	298	// now loop over the resto of the digits and cluster into NxN cluster
	299	// we do not actually cluster yet: we keep them in the list clusterDigitList
	300	nextdigitC.Reset();
	301	while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
	302	{ // scan over the list of digitsC
	303	if (digit == pMaxEnergyDigit) continue;
	304	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	305	AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
	306	//if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold )
	307	if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 )
ee08edde	308	{
	309	Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
	310	if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
	311	Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
	312	if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
51e4198e	313	{
	314	clusterDigitList.AddLast(digit) ;
	315	clusterCandidateEnergy += dEnergyCalibrated;
	316	}
ee08edde	317	}
51e4198e	318	}// loop over the next digits
	319
	320	// start a cluster here only if a cluster energy is larger than clustering threshold
	321	//if (clusterCandidateEnergy > 0.1)
	322	AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
	323	if (clusterCandidateEnergy > fECAClusteringThreshold)
	324	{
	325	if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
ee08edde	326
51e4198e	327	AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
	328	fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
	329	recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
	330	if(recPoint){
	331	fNumberOfECAClusters++ ;
	332	recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
	333
	334	AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
	335	for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
	336	{
	337
	338	digit = (AliEMCALDigit*)clusterDigitList.At(idig);
	339	Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
	340	AliDebug(5, Form(" Adding digit %d", digit->GetId()));
	341	// note: this way the sharing info is lost!
	342	recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
	343	digitsC.Remove(digit);
	344	}
	345	}// recpoint
	346	}
	347	else
	348	{
	349	// we do not want to start clustering in the same spot!
	350	// but in this case we may NOT reuse this seed for another cluster!
	351	// need a better bookeeping?
	352	digitsC.Remove(pMaxEnergyDigit);
	353	}
	354
	355	AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
	356	} // while ! done
	357
ee08edde	358	//delete digitsC ; //nope we use an object
	359
	360	AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
	361	}
	362
ee08edde	363	//___________________________________________________________________
	364	void AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
	365	{
	366	//
	367	// input is calibrated - the case when we run already on ESD
	368	//
	369	AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;
	370	}