[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerFixedWindow.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

// --- Root ---
#include <TMath.h> 
#include <TMinuit.h>
#include <TTree.h> 
#include <TBenchmark.h>
#include <TBrowser.h>
#include <TROOT.h>
#include <TClonesArray.h>
#include <TH1I.h>

// --- AliRoot ---
#include "AliLog.h"
#include "AliEMCALRecPoint.h"
#include "AliEMCALDigit.h"
#include "AliEMCALGeometry.h"
#include "AliCaloCalibPedestal.h"
#include "AliEMCALCalibData.h"
#include "AliESDCaloCluster.h"
#include "AliEMCALUnfolding.h"
#include "AliEMCALFixedWindowClusterInfo.h"

#include "AliEMCALClusterizerFixedWindow.h"

ClassImp(AliEMCALClusterizerFixedWindow)

//__________________________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow()
  : AliEMCALClusterizer(), 
    fNphi(4), 
    fNeta(4), 
    fShiftPhi(2), 
    fShiftEta(2),
    fTRUshift(0),
    fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo")),
    fClustersArray(0)
{
	// Constructor
}

//__________________________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
  : AliEMCALClusterizer(geometry), 
    fNphi(4), 
    fNeta(4), 
    fShiftPhi(2), 
    fShiftEta(2),
    fTRUshift(0),
    fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo")),
    fClustersArray(0)
{
  // Constructor
}

//__________________________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
  : AliEMCALClusterizer(geometry, calib, caloped),
    fNphi(4), 
    fNeta(4), 
    fShiftPhi(2), 
    fShiftEta(2),
    fTRUshift(0),
    fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo")),
    fClustersArray(0)
{
  // Constructor
}

//__________________________________________________________________________________________
AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
{
	// Destructor
  
  delete fClustersInfo;
  delete fClustersArray;
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n) 
{
  // Set fNphi; if clusterizer already initialized gives a warning and does nothing
  
  if (fClustersArray)
    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
  else
    fNphi = n;
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n) 
{
  // Set fNeta; if clusterizer already initialized gives a warning and does nothing
  
  if (fClustersArray)
    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
  else
    fNeta = n;
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s) 
{
  // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
  
  if (fClustersArray)
    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
  else
    fShiftPhi = s;
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s) 
{
  // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
  
  if (fClustersArray)
    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
  else
    fShiftEta = s;
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b) 
{
  // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
  
  if (fClustersArray)
    AliWarning("Clusterizer already initialized. Unable to change the parameters.");
  else
    fTRUshift = b;
}


//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
{
	// Steering method to perform clusterization for the current event 
	
	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);
		}
	}
	
	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"));
	}    
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::MakeClusters()
{
	// Make clusters
	
	if (fGeom == 0) 
		AliFatal("Did not get geometry from EMCALLoader");
	
	fNumberOfECAClusters = 0;
	fRecPoints->Delete();
  
  if (fClustersInfo->GetLastElementPosition() > 0)
    fClustersInfo->Clear();
	
	Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
	
	// Defining geometry and clusterization parameter
	Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
	Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
  
  Int_t nTRUPhi = 1;
  Int_t nTRUEta = 1;
  
  Int_t nEtaDigits = nEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
  Int_t nPhiDigits = nPhiDigitsSupMod * fGeom->GetNPhiSuperModule();    
  
  if (fTRUshift){
    nTRUPhi = fGeom->GetNPhiSuperModule() * 3;
    nTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
    nEtaDigits /= nTRUEta;
    nPhiDigits /= nTRUPhi;
  }

  // Check if clusterizer parameter are compatible with calorimeter geometry
  if (nEtaDigits < fNeta){
		AliFatal(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,nEtaDigits));
		return;
	}
	if (nPhiDigits < fNphi){
		AliFatal(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,nPhiDigits));
		return;
	}
	if (nEtaDigits % fShiftEta != 0){
		AliFatal(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,nEtaDigits));
		return;
	}
	if (nPhiDigits % fShiftPhi != 0){
		AliFatal(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,nPhiDigits));
		return;
	}
	if (fNeta % fShiftEta != 0){
		AliFatal(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
		return;
	}
	if (fNphi % fShiftPhi != 0){
		AliFatal(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
		return;
	}
  
  Int_t maxiShiftPhi = fNphi / fShiftPhi;
  Int_t maxiShiftEta = fNeta / fShiftEta;
	
	Int_t nDigitsCluster = fNphi * fNeta;
  
  Int_t nClusEtaNoShift = nEtaDigits / fNeta;
  Int_t nClusPhiNoShift = nPhiDigits / fNphi;
  
  Int_t nClusters =  nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
  
  Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
  
  if (!fClustersArray){
    fClustersArray = new AliEMCALDigit**[nTotalClus];
    for (Int_t i = 0; i < nTotalClus; i++)
    {
      fClustersArray[i] = NULL;
    }
  }
  
  AliEMCALDigit *digit = 0;
  
  for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++){
    Int_t nClusPhi = (nPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
    
    for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++){
      
      Int_t nClusEta = (nEtaDigits - fShiftEta * ishiftEta) / fNeta; 
      
      Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
      
      TIter nextdigit(fDigitsArr);
      
      nextdigit.Reset();
      
      while ((digit = static_cast<AliEMCALDigit*>(nextdigit()))){
        fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
        fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
        
        Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
        
        Int_t iTRUphi = iphi_eff / nPhiDigits;
        
        iphi_eff -= iTRUphi * nPhiDigits;
        
        Int_t iClusPhi = iphi_eff / fNphi; 
        
        if (iphi_eff < 0 || iClusPhi >= nClusPhi) 
          continue;
        
        Int_t ieta_eff = ieta - fShiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
        
        Int_t iTRUeta = ieta_eff / nEtaDigits;
        
        ieta_eff -= iTRUeta * nEtaDigits;
        
        Int_t iClusEta = ieta_eff / fNeta; 
        
        if (ieta_eff < 0 || iClusEta >= nClusEta) 
          continue;
        
        iphi_eff += iTRUphi * nPhiDigits;
        iClusPhi = iphi_eff / fNphi; 
        
        ieta_eff += iTRUeta * nEtaDigits;
        iClusEta = ieta_eff / fNeta; 
        
        Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi; 
        Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;

        if (iCluster >= nClusters){
          AliWarning(Form("iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
          return;
        }
        
        iCluster += iTotalClus;
        
        if (fClustersArray[iCluster] == NULL){
          fNumberOfECAClusters++;
          fClustersArray[iCluster] = new AliEMCALDigit*[nDigitsCluster];
          for (Int_t i = 0; i < nDigitsCluster; i++){
            fClustersArray[iCluster][i] = NULL;
          }
          
          fClustersInfo->Add(iCluster, -1, iClusEta, iClusPhi);
        }
        
        if (fClustersArray[iCluster][iDigit] != NULL){
          AliWarning("Digit already added!");
          return;
        }
        
        fClustersArray[iCluster][iDigit] = digit;
        
      } // loop on digit
      
    } // loop on eta shift
    
	} // loop on phi shift
	
  Int_t iRecPoint = 0;
	for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++){
    
		if (fClustersArray[iCluster] == NULL) continue;
		
		(*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
		AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
		
		if (recPoint) {
      if (fClustersInfo->ContainsIndex(iRecPoint))
        AliFatal(Form("ERROR: index present already, %d", iRecPoint));
      
      fClustersInfo->SetIndexFromId(iCluster, iRecPoint);
      
			iRecPoint++;       
			recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);

			for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
				if (fClustersArray[iCluster][iDigit] == NULL) continue;
				digit = fClustersArray[iCluster][iDigit];
        Float_t dEnergyCalibrated = digit->GetAmplitude();
        Float_t time              = digit->GetTime();
        Calibrate(dEnergyCalibrated,time,digit->GetId());
				digit->SetCalibAmp(dEnergyCalibrated);
				recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE); //Time or TimeR?
        fClustersArray[iCluster][iDigit] = NULL;
			}
		}
    
    delete[] fClustersArray[iCluster];
    fClustersArray[iCluster] = NULL;
	}
  
	AliDebug(1, Form("MakeClusters: Number of digits %d  -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
	
	AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast())); 
}
Commit	Line	Data
89504c57	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
27bad386	16	// This class derives from AliEMCALClustrerizer
89504c57	17
27bad386	18	// --- Root ---
89504c57	19	#include <TMath.h>
	20	#include <TMinuit.h>
	21	#include <TTree.h>
	22	#include <TBenchmark.h>
	23	#include <TBrowser.h>
	24	#include <TROOT.h>
	25	#include <TClonesArray.h>
	26	#include <TH1I.h>
	27
27bad386	28	// --- AliRoot ---
89504c57	29	#include "AliLog.h"
	30	#include "AliEMCALRecPoint.h"
	31	#include "AliEMCALDigit.h"
	32	#include "AliEMCALGeometry.h"
	33	#include "AliCaloCalibPedestal.h"
	34	#include "AliEMCALCalibData.h"
	35	#include "AliESDCaloCluster.h"
	36	#include "AliEMCALUnfolding.h"
	37	#include "AliEMCALFixedWindowClusterInfo.h"
	38
27bad386	39	#include "AliEMCALClusterizerFixedWindow.h"
27bad386	40
89504c57	41	ClassImp(AliEMCALClusterizerFixedWindow)
89504c57	42
27bad386	43	//__________________________________________________________________________________________
89504c57	44	AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow()
27bad386	45	: AliEMCALClusterizer(),
	46	fNphi(4),
	47	fNeta(4),
	48	fShiftPhi(2),
	49	fShiftEta(2),
	50	fTRUshift(0),
	51	fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo")),
	52	fClustersArray(0)
89504c57	53	{
27bad386	54	// Constructor
89504c57	55	}
89504c57	56
27bad386	57	//__________________________________________________________________________________________
89504c57	58	AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
27bad386	59	: AliEMCALClusterizer(geometry),
	60	fNphi(4),
	61	fNeta(4),
	62	fShiftPhi(2),
	63	fShiftEta(2),
	64	fTRUshift(0),
	65	fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo")),
	66	fClustersArray(0)
89504c57	67	{
27bad386	68	// Constructor
89504c57	69	}
89504c57	70
27bad386	71	//__________________________________________________________________________________________
89504c57	72	AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
27bad386	73	: AliEMCALClusterizer(geometry, calib, caloped),
	74	fNphi(4),
	75	fNeta(4),
	76	fShiftPhi(2),
	77	fShiftEta(2),
	78	fTRUshift(0),
	79	fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo")),
	80	fClustersArray(0)
89504c57	81	{
27bad386	82	// Constructor
89504c57	83	}
89504c57	84
27bad386	85	//__________________________________________________________________________________________
89504c57	86	AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
89504c57	87	{
27bad386	88	// Destructor
27bad386	89
89504c57	90	delete fClustersInfo;
27bad386	91	delete fClustersArray;
	92	}
	93
	94	//__________________________________________________________________________________________
	95	void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n)
	96	{
	97	// Set fNphi; if clusterizer already initialized gives a warning and does nothing
	98
	99	if (fClustersArray)
	100	AliWarning("Clusterizer already initialized. Unable to change the parameters.");
	101	else
	102	fNphi = n;
	103	}
	104
	105	//__________________________________________________________________________________________
	106	void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n)
	107	{
	108	// Set fNeta; if clusterizer already initialized gives a warning and does nothing
	109
	110	if (fClustersArray)
	111	AliWarning("Clusterizer already initialized. Unable to change the parameters.");
	112	else
	113	fNeta = n;
	114	}
	115
	116	//__________________________________________________________________________________________
	117	void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s)
	118	{
	119	// Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
	120
	121	if (fClustersArray)
	122	AliWarning("Clusterizer already initialized. Unable to change the parameters.");
	123	else
	124	fShiftPhi = s;
	125	}
	126
	127	//__________________________________________________________________________________________
	128	void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s)
	129	{
	130	// Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
	131
	132	if (fClustersArray)
	133	AliWarning("Clusterizer already initialized. Unable to change the parameters.");
	134	else
	135	fShiftEta = s;
	136	}
	137
	138	//__________________________________________________________________________________________
	139	void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b)
	140	{
	141	// Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
	142
	143	if (fClustersArray)
	144	AliWarning("Clusterizer already initialized. Unable to change the parameters.");
	145	else
	146	fTRUshift = b;
89504c57	147	}
89504c57	148
27bad386	149
27bad386	150	//__________________________________________________________________________________________
89504c57	151	void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
	152	{
	153	// Steering method to perform clusterization for the current event
89504c57	154
	155	if (strstr(option,"tim"))
	156	gBenchmark->Start("EMCALClusterizer");
	157
	158	if (strstr(option,"print"))
	159	Print("");
	160
	161	//Get calibration parameters from file or digitizer default values.
	162	GetCalibrationParameters();
	163
	164	//Get dead channel map from file or digitizer default values.
	165	GetCaloCalibPedestal();
	166
	167	MakeClusters(); //only the real clusters
	168
	169	if (fToUnfold) {
	170	fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
	171	fClusterUnfolding->MakeUnfolding();
	172	}
	173
	174	//Evaluate position, dispersion and other RecPoint properties for EC section
	175	for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
	176	AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
	177	if (rp) {
	178	rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
	179	AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
	180	//For each rec.point set the distance to the nearest bad crystal
	181	if (fCaloPed)
	182	rp->EvalDistanceToBadChannels(fCaloPed);
	183	}
	184	}
	185
89504c57	186	for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
	187	AliEMCALRecPoint rp = dynamic_cast<AliEMCALRecPoint >(fRecPoints->At(index));
	188	if (rp) {
	189	rp->SetIndexInList(index);
	190	}
	191	else AliFatal("RecPoint NULL!!");
	192	}
	193
	194	if (fTreeR)
	195	fTreeR->Fill();
	196
	197	if (strstr(option,"deb") \|\| strstr(option,"all"))
	198	PrintRecPoints(option);
	199
	200	AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
	201
	202	if (strstr(option,"tim")) {
	203	gBenchmark->Stop("EMCALClusterizer");
	204	printf("Exec took %f seconds for Clusterizing",
	205	gBenchmark->GetCpuTime("EMCALClusterizer"));
	206	}
	207	}
	208
27bad386	209	//__________________________________________________________________________________________
89504c57	210	void AliEMCALClusterizerFixedWindow::MakeClusters()
	211	{
	212	// Make clusters
	213
	214	if (fGeom == 0)
	215	AliFatal("Did not get geometry from EMCALLoader");
	216
	217	fNumberOfECAClusters = 0;
	218	fRecPoints->Delete();
	219
27bad386	220	if (fClustersInfo->GetLastElementPosition() > 0)
89504c57	221	fClustersInfo->Clear();
	222
	223	Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
	224
	225	// Defining geometry and clusterization parameter
	226	Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
	227	Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
	228
	229	Int_t nTRUPhi = 1;
	230	Int_t nTRUEta = 1;
	231
	232	Int_t nEtaDigits = nEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
	233	Int_t nPhiDigits = nPhiDigitsSupMod * fGeom->GetNPhiSuperModule();
	234
27bad386	235	if (fTRUshift){
89504c57	236	nTRUPhi = fGeom->GetNPhiSuperModule() * 3;
	237	nTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
	238	nEtaDigits /= nTRUEta;
	239	nPhiDigits /= nTRUPhi;
	240	}
	241
	242	// Check if clusterizer parameter are compatible with calorimeter geometry
27bad386	243	if (nEtaDigits < fNeta){
27bad386	244	AliFatal(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,nEtaDigits));
89504c57	245	return;
89504c57	246	}
27bad386	247	if (nPhiDigits < fNphi){
27bad386	248	AliFatal(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,nPhiDigits));
89504c57	249	return;
89504c57	250	}
27bad386	251	if (nEtaDigits % fShiftEta != 0){
27bad386	252	AliFatal(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,nEtaDigits));
89504c57	253	return;
89504c57	254	}
27bad386	255	if (nPhiDigits % fShiftPhi != 0){
27bad386	256	AliFatal(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,nPhiDigits));
89504c57	257	return;
89504c57	258	}
27bad386	259	if (fNeta % fShiftEta != 0){
27bad386	260	AliFatal(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
89504c57	261	return;
89504c57	262	}
27bad386	263	if (fNphi % fShiftPhi != 0){
27bad386	264	AliFatal(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
89504c57	265	return;
	266	}
	267
27bad386	268	Int_t maxiShiftPhi = fNphi / fShiftPhi;
27bad386	269	Int_t maxiShiftEta = fNeta / fShiftEta;
89504c57	270
27bad386	271	Int_t nDigitsCluster = fNphi * fNeta;
89504c57	272
27bad386	273	Int_t nClusEtaNoShift = nEtaDigits / fNeta;
27bad386	274	Int_t nClusPhiNoShift = nPhiDigits / fNphi;
89504c57	275
	276	Int_t nClusters = nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
	277
	278	Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
	279
27bad386	280	if (!fClustersArray){
27bad386	281	fClustersArray = new AliEMCALDigit**[nTotalClus];
89504c57	282	for (Int_t i = 0; i < nTotalClus; i++)
89504c57	283	{
27bad386	284	fClustersArray[i] = NULL;
89504c57	285	}
	286	}
	287
	288	AliEMCALDigit *digit = 0;
	289
27bad386	290	for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++){
27bad386	291	Int_t nClusPhi = (nPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
89504c57	292
27bad386	293	for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++){
89504c57	294
27bad386	295	Int_t nClusEta = (nEtaDigits - fShiftEta * ishiftEta) / fNeta;
89504c57	296
	297	Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
	298
	299	TIter nextdigit(fDigitsArr);
	300
	301	nextdigit.Reset();
	302
27bad386	303	while ((digit = static_cast<AliEMCALDigit*>(nextdigit()))){
89504c57	304	fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
	305	fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
	306
27bad386	307	Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
89504c57	308
	309	Int_t iTRUphi = iphi_eff / nPhiDigits;
	310
	311	iphi_eff -= iTRUphi * nPhiDigits;
	312
27bad386	313	Int_t iClusPhi = iphi_eff / fNphi;
89504c57	314
	315	if (iphi_eff < 0 \|\| iClusPhi >= nClusPhi)
	316	continue;
	317
27bad386	318	Int_t ieta_eff = ieta - fShiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
89504c57	319
	320	Int_t iTRUeta = ieta_eff / nEtaDigits;
	321
	322	ieta_eff -= iTRUeta * nEtaDigits;
	323
27bad386	324	Int_t iClusEta = ieta_eff / fNeta;
89504c57	325
	326	if (ieta_eff < 0 \|\| iClusEta >= nClusEta)
	327	continue;
	328
	329	iphi_eff += iTRUphi * nPhiDigits;
27bad386	330	iClusPhi = iphi_eff / fNphi;
89504c57	331
89504c57	332	ieta_eff += iTRUeta * nEtaDigits;
27bad386	333	iClusEta = ieta_eff / fNeta;
89504c57	334
89504c57	335	Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi;
27bad386	336	Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
89504c57	337
27bad386	338	if (iCluster >= nClusters){
27bad386	339	AliWarning(Form("iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
89504c57	340	return;
	341	}
	342
	343	iCluster += iTotalClus;
	344
27bad386	345	if (fClustersArray[iCluster] == NULL){
89504c57	346	fNumberOfECAClusters++;
27bad386	347	fClustersArray[iCluster] = new AliEMCALDigit*[nDigitsCluster];
	348	for (Int_t i = 0; i < nDigitsCluster; i++){
	349	fClustersArray[iCluster][i] = NULL;
89504c57	350	}
	351
	352	fClustersInfo->Add(iCluster, -1, iClusEta, iClusPhi);
	353	}
	354
27bad386	355	if (fClustersArray[iCluster][iDigit] != NULL){
27bad386	356	AliWarning("Digit already added!");
89504c57	357	return;
	358	}
	359
27bad386	360	fClustersArray[iCluster][iDigit] = digit;
89504c57	361
	362	} // loop on digit
	363
	364	} // loop on eta shift
	365
	366	} // loop on phi shift
	367
	368	Int_t iRecPoint = 0;
27bad386	369	for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++){
89504c57	370
27bad386	371	if (fClustersArray[iCluster] == NULL) continue;
89504c57	372
	373	(*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
	374	AliEMCALRecPoint recPoint = dynamic_cast<AliEMCALRecPoint> (fRecPoints->At(iRecPoint));
	375
27bad386	376	if (recPoint) {
89504c57	377	if (fClustersInfo->ContainsIndex(iRecPoint))
	378	AliFatal(Form("ERROR: index present already, %d", iRecPoint));
	379
	380	fClustersInfo->SetIndexFromId(iCluster, iRecPoint);
	381
	382	iRecPoint++;
	383	recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
27bad386	384
	385	for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
	386	if (fClustersArray[iCluster][iDigit] == NULL) continue;
	387	digit = fClustersArray[iCluster][iDigit];
89504c57	388	Float_t dEnergyCalibrated = digit->GetAmplitude();
	389	Float_t time = digit->GetTime();
	390	Calibrate(dEnergyCalibrated,time,digit->GetId());
	391	digit->SetCalibAmp(dEnergyCalibrated);
	392	recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE); //Time or TimeR?
27bad386	393	fClustersArray[iCluster][iDigit] = NULL;
89504c57	394	}
	395	}
	396
27bad386	397	delete[] fClustersArray[iCluster];
27bad386	398	fClustersArray[iCluster] = NULL;
89504c57	399	}
89504c57	400
27bad386	401	AliDebug(1, Form("MakeClusters: Number of digits %d -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
89504c57	402
27bad386	403	AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast()));
89504c57	404	}