[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(), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
{
  // ctor with the indication of the file where header Tree and digits Tree are stored
}

//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
  : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
{
  // ctor with the indication of the file where header Tree and digits Tree are stored
  // use this contructor to avoid usage of Init() which uses runloader
  // change needed by HLT - MP
}

//____________________________________________________________________________
AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
: AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
{
  // 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();
	
  MakeClusters();  //only the real clusters
  
  if (fToUnfold) {
    fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
    fClusterUnfolding->MakeUnfolding();
  }
  
  //Evaluate position, dispersion and other RecPoint properties for EC section 
  for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { 
    AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
    if (rp) {
      rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
      AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
      //For each rec.point set the distance to the nearest bad crystal
      if (fCaloPed)
        rp->EvalDistanceToBadChannels(fCaloPed);
    }
  }
  
  fRecPoints->Sort();
  
  for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
    AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
    if (rp) {
      rp->SetIndexInList(index);
    }
    else AliFatal("RecPoint NULL!!");
  }
  
  if (fTreeR)
    fTreeR->Fill();
  
  if (strstr(option,"deb") || strstr(option,"all"))  
    PrintRecPoints(option);
  
  AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
  
  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  
  
  if (fEnergyGrad) { //false by default
    if (d2->GetCalibAmp()>d1->GetCalibAmp())
      return 3; // energy of neighboring cell should be smaller in order to become a neighbor
  }

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