[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 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] AliEMCALClusterizerFixedWindow * cl = new AliEMCALClusterizerFixedWindow("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

#include "AliEMCALClusterizerFixedWindow.h"

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

// --- Standard library ---
#include <cassert>
//#include <iostream>
//#include <fstream>

// --- AliRoot header files ---
#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"

ClassImp(AliEMCALClusterizerFixedWindow)

//____________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow()
: AliEMCALClusterizer(), 
nPhi(4), 
nEta(4), 
shiftPhi(2), 
shiftEta(2),
fTRUshift(0), 
clusters_array(0),
fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
{
	// ctor with the indication of the file where header Tree and digits Tree are stored

}

//____________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
: AliEMCALClusterizer(geometry), 
nPhi(4), 
nEta(4), 
shiftPhi(2), 
shiftEta(2), 
fTRUshift(0), 
clusters_array(0),
fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
{
	// 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
}

//____________________________________________________________________________
AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
: AliEMCALClusterizer(geometry, calib, caloped),
nPhi(4), 
nEta(4), 
shiftPhi(2), 
shiftEta(2), 
fTRUshift(0), 
clusters_array(0),
fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
{
	// ctor, geometry and calibration are initialized elsewhere.
}

//____________________________________________________________________________
AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
{
	// dtor
  delete fClustersInfo;
}

//____________________________________________________________________________
void AliEMCALClusterizerFixedWindow::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"));
	}    
}

//____________________________________________________________________________
void AliEMCALClusterizerFixedWindow::MakeClusters()
{
	// Make clusters
	
	if (fGeom == 0) 
		AliFatal("Did not get geometry from EMCALLoader");
	
	fNumberOfECAClusters = 0;
	fRecPoints->Delete();
  
  if (fClustersInfo->GetLastElementId() > 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 < nEta)
	{
		AliFatal(Form("Error: nEta = %d is greater than nEtaDigits = %d.",nEta,nEtaDigits));
		return;
	}
	if (nPhiDigits < nPhi)
	{
		AliFatal(Form("Error: nPhi = %d is greater than nPhiDigits = %d.",nPhi,nPhiDigits));
		return;
	}
	if (nEtaDigits % shiftEta != 0)
	{
		AliFatal(Form("Error: shiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",shiftEta,nEtaDigits));
		return;
	}
	if (nPhiDigits % shiftPhi != 0)
	{
		AliFatal(Form("Error: shiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",shiftPhi,nPhiDigits));
		return;
	}
	if (nEta % shiftEta != 0)
	{
		AliFatal(Form("Error: shiftEta = %d is not divisor of nEta = %d.",shiftEta,nEta));
		return;
	}
	if (nPhi % shiftPhi != 0)
	{
		AliFatal(Form("Error: shiftPhi = %d is not divisor of nPhi = %d).",shiftPhi,nPhi));
		return;
	}
  
  Int_t maxiShiftPhi = nPhi / shiftPhi;
  Int_t maxiShiftEta = nEta / shiftEta;
	
	Int_t nDigitsCluster = nPhi * nEta;
  
  Int_t nClusEtaNoShift = nEtaDigits / nEta;
  Int_t nClusPhiNoShift = nPhiDigits / nPhi;
  
  Int_t nClusters =  nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
  
  Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
  
  if (!clusters_array)
  {
    clusters_array = new AliEMCALDigit**[nTotalClus];
    for (Int_t i = 0; i < nTotalClus; i++)
    {
      clusters_array[i] = NULL;
    }
  }
  
  AliEMCALDigit *digit = 0;
  
  for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++)
  {
    Int_t nClusPhi = (nPhiDigits - shiftPhi * ishiftPhi) / nPhi;
    
    for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++)
    {
      
      Int_t nClusEta = (nEtaDigits - shiftEta * ishiftEta) / nEta; 
      
      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 - shiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
        
        Int_t iTRUphi = iphi_eff / nPhiDigits;
        
        iphi_eff -= iTRUphi * nPhiDigits;
        
        Int_t iClusPhi = iphi_eff / nPhi; 
        
        if (iphi_eff < 0 || iClusPhi >= nClusPhi) 
          continue;
        
        Int_t ieta_eff = ieta - shiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
        
        Int_t iTRUeta = ieta_eff / nEtaDigits;
        
        ieta_eff -= iTRUeta * nEtaDigits;
        
        Int_t iClusEta = ieta_eff / nEta; 
        
        if (ieta_eff < 0 || iClusEta >= nClusEta) 
          continue;
        
        iphi_eff += iTRUphi * nPhiDigits;
        iClusPhi = iphi_eff / nPhi; 
        
        ieta_eff += iTRUeta * nEtaDigits;
        iClusEta = ieta_eff / nEta; 
        
        Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi; 
        Int_t iDigit = iphi_eff % nPhi + (ieta_eff % nEta) * nPhi;

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