[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 "AliEMCALClusterizerFixedWindow.h"

ClassImp(AliEMCALClusterizerFixedWindow)

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

//__________________________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
  : AliEMCALClusterizer(geometry), 
    fNphi(4), 
    fNeta(4), 
    fShiftPhi(2), 
    fShiftEta(2),
    fTRUshift(0),
    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),
    fClustersArray(0)
{
  // Constructor
}

//__________________________________________________________________________________________
AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
{
  // Destructor
  
  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);
    }
  }
  
  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"));
  }    
}

//__________________________________________________________________________________________
void AliEMCALClusterizerFixedWindow::MakeClusters()
{
  // Make clusters
	
  if (fGeom == 0) 
    AliFatal("Did not get geometry from EMCALLoader");
	
  fNumberOfECAClusters = 0;
  fRecPoints->Delete();
	
  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;
    }
  }
  
  // Set up TObjArray with pointers to digits to work on calibrated digits 
  TObjArray *digitsC = new TObjArray();
  AliEMCALDigit *digit;
  Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
  TIter nextdigit(fDigitsArr);
  while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigit()))) { // calibrate and clean up digits
    dEnergyCalibrated =  digit->GetAmplitude();
    time              =  digit->GetTime();
    Calibrate(dEnergyCalibrated, time, digit->GetId());
    digit->SetCalibAmp(dEnergyCalibrated);
    digit->SetTime(time);
    if (dEnergyCalibrated < fMinECut) {
      continue;
    }
    else if (!fGeom->CheckAbsCellId(digit->GetId())) {
      continue;
    }
    else {
      ehs += dEnergyCalibrated;
      digitsC->AddLast(digit);
    }
  } 
  
  AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
                  fDigitsArr->GetEntries(),fMinECut,ehs));
   
  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 nextdigitC(digitsC);
      while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
        
        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;
          }
        }
        
        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
	
  if(fNumberOfECAClusters >= fRecPoints->GetSize()) 
    fRecPoints->Expand(fNumberOfECAClusters+1);
  
  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) {
      iRecPoint++;       
      recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
      recPoint->SetUniqueID(iCluster);

      for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
        if (fClustersArray[iCluster][iDigit] == NULL) continue;
        digit = fClustersArray[iCluster][iDigit];
        recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
        fClustersArray[iCluster][iDigit] = NULL;
      }
    }
    
    delete[] fClustersArray[iCluster];
    fClustersArray[iCluster] = NULL;
  }

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