[u/mrichter/AliRoot.git] / HLT / PHOS / AliHLTPHOSClusterizer.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Authors: Øystein Djuvsland <oysteind@ift.uib.no>                       *
 *                                                                        *
 * 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.                  *
 **************************************************************************/


/** @file   AliHLTPHOSClusterizer.cxx
    @author Øystein Djuvsland
    @date   
    @brief  A temporary clusterizer for PHOS
*/


#include "AliHLTPHOSPhysicsDefinitions.h"
#include "AliHLTPHOSClusterizer.h"
#include "AliHLTPHOSCommonDefs.h"
#include "TVector3.h"
#include "TMath.h"
#include "AliHLTPHOSRcuCellEnergyDataStruct.h"
#include "AliHLTPHOSRecPointListDataStruct.h"
#include "AliHLTPHOSValidCellDataStruct.h"
#include "AliHLTPHOSRecPointDataStruct.h"
#include "AliHLTPHOSClusterDataStruct.h"


ClassImp(AliHLTPHOSClusterizer);

/**
* Main constructor
**/
AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():fPHOSModule(-1), fThreshold(0), fClusterThreshold(0), 
					       fHighGainFactor(0.005), fLowGainFactor(0.08),
					       fArraySize(3), fMultiplicity(fArraySize*fArraySize)
{
  //Main constructor
  
}//end

AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &):fPHOSModule(-1), fThreshold(0), fClusterThreshold(0), 
									    fHighGainFactor(0.005), fLowGainFactor(0.08),
									    fArraySize(3), fMultiplicity(fArraySize*fArraySize)
{
  //Copy constructor, not implemented
}//end

AliHLTPHOSClusterizer:: ~AliHLTPHOSClusterizer()  
{
  //Destructor
}

/**
* Building a 2D array of cell energies of the PHOS detector
* @param cellData object containing the cell energies from one event
* @param recPointList list to be filled with coordinates of local maxima in the detector
**/

Int_t
AliHLTPHOSClusterizer::BuildCellEnergyArray(AliHLTPHOSRcuCellEnergyDataStruct* cellData, 
					    AliHLTPHOSRecPointListDataStruct* recPointList)
{
  //Build the cell energy array of the detector

  Int_t x = 0;
  Int_t z = 0;
  Int_t gain = 0;
  Int_t xMod = 0;
  Int_t zMod = 0;
  Int_t index = 0;
  Double_t energyCount = 0;

  for(Int_t cell = 0; cell < cellData->fCnt; cell++)
    {
     
      z = (cellData->fValidData[cell]).fZ;
      x = (cellData->fValidData[cell]).fX;
      gain = (cellData->fValidData[cell]).fGain;

      zMod = z + (cellData->fRcuZ)*N_ROWS_RCU;
      xMod = x + (cellData->fRcuX)*N_COLUMNS_RCU;
      
      energyCount = (cellData->fValidData[cell]).fEnergy;

      if(gain == 0 && energyCount < 1023) 
	{
	  fEnergyArray[xMod][zMod] = fHighGainFactor * energyCount;
	 
	  if(fEnergyArray[xMod][zMod] > fClusterThreshold)
	    {	      
	      recPointList[index].fZ = zMod;
	      recPointList[index].fX = xMod;
	      
	      for(Int_t j = 0; j < index; j++)
		{
		  if(recPointList[j].fZ == zMod && recPointList[j].fX == xMod)
		    recPointList[j].fZ = -1;
		}
	      index++;
	    }
	  
	}
      else if(gain == 1 && fEnergyArray[xMod][zMod] == 0) 
	{
	  fEnergyArray[xMod][zMod] = fLowGainFactor * energyCount;
	   if(fEnergyArray[xMod][zMod] > fClusterThreshold)
	    {	      
	      recPointList[index].fZ = zMod;
	      recPointList[index].fX = xMod;
	      recPointList[index].fModule = cellData->fModuleID;
	      index++;
	    }
	}

    }  
  
  fPHOSModule = cellData->fModuleID;

  return index;
}//end BuildCellEnergyArray


/**
* Creating an array of rec points
* @param recPointStructArrayPtr array to store the rec points
* @param list list of rec points
* @param nPoints number of points
**/
Int_t 
AliHLTPHOSClusterizer::CreateRecPointStructArray(AliHLTPHOSRecPointDataStruct* recPointStructArrayPtr, 
						 AliHLTPHOSRecPointListDataStruct* list, 
						 Int_t nPoints) 

{
  //Create the rec point struct array

  Int_t flag = 0;
  Int_t edgeFlagRows = 0;
  Int_t edgeFlagCols = 0;
  Int_t k = 0;
  Int_t nRecPoints = 0;
  Int_t z = 0;
  Int_t x = 0;

  Float_t* energiesList = NULL;
  
  for(Int_t point = 0; point < nPoints; point++)
    {
      z = list[point].fZ;
      x = list[point].fX;
      if(z == -1) continue;

      if((z-fArraySize/2) < 0/*= - N_ROWS_MOD/2*/ || (z+fArraySize/2) >= N_ROWS_MOD/*/2*/) 
	{
	  edgeFlagRows = 1;
	  continue;
	}
      if((x-fArraySize/2) < 0/*= - N_COLUMNS_MOD/2*/ || (x+fArraySize/2) >= N_COLUMNS_MOD) 
	{
	  edgeFlagCols = 1;
	  continue;
	}
      

      if(!flag) energiesList = new Float_t[fMultiplicity];
      flag = 0;
      k = 0;	      
	  
      for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	{
	  if(flag) break;
	  for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	    {			  
	      
	      if(fEnergyArray[x+i][z+j] > fEnergyArray[x][z] && abs(i) < 2 && abs(j) < 2)
		{
		  flag = 1;
		  break;
		}
	      energiesList[k] = fEnergyArray[x+i][z+j];
	      k++;
	    }
	}
	  
      
      if(!flag && k)
	{
	  recPointStructArrayPtr[nRecPoints].fEnergiesListPtr = energiesList;
	  recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[0] = x;
	  recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[1] = z;
	  recPointStructArrayPtr[nRecPoints].fX = x;
	  recPointStructArrayPtr[nRecPoints].fZ = z;
	  recPointStructArrayPtr[nRecPoints].fMultiplicity = fMultiplicity;
	  recPointStructArrayPtr[nRecPoints].fPHOSModule = list[point].fModule;
	  nRecPoints++;
	}
    }


  energiesList = NULL;

  return nRecPoints;
  
}//end CreateRecPointStructArray


/**
* Calculating the center of gravity of a rec point
* Not working well at this point!
* @param recPointPtr pointer to the rec point
**/
Int_t
AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
{
  //Calculate the center of gravity

  Float_t xt = 0;
  Float_t zt = 0;
  Float_t xi = 0;
  Float_t zj = 0;
  Float_t w = 0;
  Float_t w0 = 4.5;
  Float_t wtot = 0;

  Int_t x = recPointPtr->fCoordinatesPtr[0];
  Int_t z = recPointPtr->fCoordinatesPtr[1];
  
  for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
    {
      for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	{			  
	  xi = x + i;
	  zj = z + j;
	  w = TMath::Max( (Float_t)0., (Float_t)(w0 + log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]) ) ) ;
	  xt += xi * w ;
	  zt += zj * w ;
	  wtot += w ;
	  
	}
    }
  xt /= wtot ;
  zt /= wtot ;
  
  recPointPtr->fX = xt;
  recPointPtr->fZ = zt;
  
  return 0;

}//end CalculateCenterOfGravity
  

/**
* Create a simpler data struct for a rec point, containing
* only coordinates, energy and timing
* @param recPointPtr pointer to the rec point
* @param clusterStructPtr pointer to the cluster struct
**/
Int_t
AliHLTPHOSClusterizer::ClusterizeStruct(AliHLTPHOSRecPointDataStruct* recPointPtr, AliHLTPHOSClusterDataStruct* clusterStructPtr)
{
  //Simplify the rec points

  Float_t clusterEnergy = 0;
  Float_t* energiesListPtr = recPointPtr->fEnergiesListPtr;
  
  for(Int_t i = 0; i < recPointPtr->fMultiplicity; i++)
    {
      clusterEnergy += energiesListPtr[i];
    }

  clusterStructPtr->fLocalPositionPtr[0] = recPointPtr->fX;
  clusterStructPtr->fLocalPositionPtr[1] = recPointPtr->fZ;
  clusterStructPtr->fClusterEnergy = clusterEnergy;
  clusterStructPtr->fPHOSModule = recPointPtr->fPHOSModule;

  return 0;

}//end ClusterizeStruct


/**
* Resets the cell energy array
**/ 
Int_t 
AliHLTPHOSClusterizer::ResetCellEnergyArray()

{
  //Reset the cell energy array

  for(Int_t x = 0; x < N_ROWS_MOD; x++)
    {
      for(Int_t z = 0; z < N_COLUMNS_MOD; z++)
	{
	  fEnergyArray[x][z] = 0;
	}
    }

  return 0;

}//end ResetCellEnergyArray
Commit	Line	Data
aac22523	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Authors: Øystein Djuvsland <oysteind@ift.uib.no> *
	5	* *
	6	* Permission to use, copy, modify and distribute this software and its *
	7	* documentation strictly for non-commercial purposes is hereby granted *
	8	* without fee, provided that the above copyright notice appears in all *
	9	* copies and that both the copyright notice and this permission notice *
	10	* appear in the supporting documentation. The authors make no claims *
	11	* about the suitability of this software for any purpose. It is *
	12	* provided "as is" without express or implied warranty. *
	13	**************************************************************************/
	14
	15
	16	/** @file AliHLTPHOSClusterizer.cxx
	17	@author Øystein Djuvsland
	18	@date
	19	@brief A temporary clusterizer for PHOS
	20	*/
	21
	22
	23
	24	#include "AliHLTPHOSPhysicsDefinitions.h"
	25	#include "AliHLTPHOSClusterizer.h"
	26	#include "AliHLTPHOSCommonDefs.h"
	27	#include "TVector3.h"
aac22523	28	#include "TMath.h"
91b95d47	29	#include "AliHLTPHOSRcuCellEnergyDataStruct.h"
	30	#include "AliHLTPHOSRecPointListDataStruct.h"
	31	#include "AliHLTPHOSValidCellDataStruct.h"
	32	#include "AliHLTPHOSRecPointDataStruct.h"
	33	#include "AliHLTPHOSClusterDataStruct.h"
aac22523	34
	35
	36	ClassImp(AliHLTPHOSClusterizer);
	37
	38	/**
	39	* Main constructor
	40	**/
91b95d47	41	AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():fPHOSModule(-1), fThreshold(0), fClusterThreshold(0),
91b95d47	42	fHighGainFactor(0.005), fLowGainFactor(0.08),
aac22523	43	fArraySize(3), fMultiplicity(fArraySize*fArraySize)
aac22523	44	{
91b95d47	45	//Main constructor
aac22523	46
	47	}//end
	48
91b95d47	49	AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &):fPHOSModule(-1), fThreshold(0), fClusterThreshold(0),
91b95d47	50	fHighGainFactor(0.005), fLowGainFactor(0.08),
aac22523	51	fArraySize(3), fMultiplicity(fArraySize*fArraySize)
aac22523	52	{
91b95d47	53	//Copy constructor, not implemented
aac22523	54	}//end
	55
	56	AliHLTPHOSClusterizer:: ~AliHLTPHOSClusterizer()
	57	{
91b95d47	58	//Destructor
aac22523	59	}
	60
	61	/**
	62	* Building a 2D array of cell energies of the PHOS detector
	63	* @param cellData object containing the cell energies from one event
	64	* @param recPointList list to be filled with coordinates of local maxima in the detector
	65	**/
	66
	67	Int_t
	68	AliHLTPHOSClusterizer::BuildCellEnergyArray(AliHLTPHOSRcuCellEnergyDataStruct* cellData,
	69	AliHLTPHOSRecPointListDataStruct* recPointList)
	70	{
91b95d47	71	//Build the cell energy array of the detector
aac22523	72
	73	Int_t x = 0;
	74	Int_t z = 0;
	75	Int_t gain = 0;
	76	Int_t xMod = 0;
	77	Int_t zMod = 0;
	78	Int_t index = 0;
	79	Double_t energyCount = 0;
	80
	81	for(Int_t cell = 0; cell < cellData->fCnt; cell++)
	82	{
	83
	84	z = (cellData->fValidData[cell]).fZ;
	85	x = (cellData->fValidData[cell]).fX;
	86	gain = (cellData->fValidData[cell]).fGain;
	87
	88	zMod = z + (cellData->fRcuZ)*N_ROWS_RCU;
	89	xMod = x + (cellData->fRcuX)*N_COLUMNS_RCU;
	90
	91	energyCount = (cellData->fValidData[cell]).fEnergy;
	92
	93	if(gain == 0 && energyCount < 1023)
	94	{
	95	fEnergyArray[xMod][zMod] = fHighGainFactor * energyCount;
	96
	97	if(fEnergyArray[xMod][zMod] > fClusterThreshold)
	98	{
	99	recPointList[index].fZ = zMod;
	100	recPointList[index].fX = xMod;
	101
	102	for(Int_t j = 0; j < index; j++)
	103	{
	104	if(recPointList[j].fZ == zMod && recPointList[j].fX == xMod)
	105	recPointList[j].fZ = -1;
	106	}
	107	index++;
	108	}
	109
	110	}
	111	else if(gain == 1 && fEnergyArray[xMod][zMod] == 0)
	112	{
	113	fEnergyArray[xMod][zMod] = fLowGainFactor * energyCount;
	114	if(fEnergyArray[xMod][zMod] > fClusterThreshold)
	115	{
	116	recPointList[index].fZ = zMod;
	117	recPointList[index].fX = xMod;
	118	recPointList[index].fModule = cellData->fModuleID;
	119	index++;
	120	}
	121	}
	122
	123	}
	124
	125	fPHOSModule = cellData->fModuleID;
	126
	127	return index;
	128	}//end BuildCellEnergyArray
	129
	130
	131
	132	/**
	133	* Creating an array of rec points
	134	* @param recPointStructArrayPtr array to store the rec points
	135	* @param list list of rec points
136	* @param nPoints number of points
137	**/
138	Int_t
139	AliHLTPHOSClusterizer::CreateRecPointStructArray(AliHLTPHOSRecPointDataStruct* recPointStructArrayPtr,
140	AliHLTPHOSRecPointListDataStruct* list,
141	Int_t nPoints)
142
143	{
91b95d47	144	//Create the rec point struct array
91b95d47	145
aac22523	146	Int_t flag = 0;
	147	Int_t edgeFlagRows = 0;
	148	Int_t edgeFlagCols = 0;
	149	Int_t k = 0;
	150	Int_t nRecPoints = 0;
	151	Int_t z = 0;
	152	Int_t x = 0;
	153
	154	Float_t* energiesList = NULL;
	155
	156	for(Int_t point = 0; point < nPoints; point++)
	157	{
	158	z = list[point].fZ;
	159	x = list[point].fX;
	160	if(z == -1) continue;
	161
	162	if((z-fArraySize/2) < 0/= - N_ROWS_MOD/2/ \|\| (z+fArraySize/2) >= N_ROWS_MOD//2/)
	163	{
	164	edgeFlagRows = 1;
	165	continue;
	166	}
	167	if((x-fArraySize/2) < 0/= - N_COLUMNS_MOD/2/ \|\| (x+fArraySize/2) >= N_COLUMNS_MOD)
	168	{
	169	edgeFlagCols = 1;
	170	continue;
	171	}
	172
	173
	174	if(!flag) energiesList = new Float_t[fMultiplicity];
	175	flag = 0;
	176	k = 0;
	177
	178	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	179	{
	180	if(flag) break;
	181	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	182	{
	183
	184	if(fEnergyArray[x+i][z+j] > fEnergyArray[x][z] && abs(i) < 2 && abs(j) < 2)
	185	{
	186	flag = 1;
	187	break;
	188	}
	189	energiesList[k] = fEnergyArray[x+i][z+j];
	190	k++;
	191	}
	192	}
	193
	194
	195	if(!flag && k)
	196	{
	197	recPointStructArrayPtr[nRecPoints].fEnergiesListPtr = energiesList;
	198	recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[0] = x;
	199	recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[1] = z;
	200	recPointStructArrayPtr[nRecPoints].fX = x;
	201	recPointStructArrayPtr[nRecPoints].fZ = z;
	202	recPointStructArrayPtr[nRecPoints].fMultiplicity = fMultiplicity;
	203	recPointStructArrayPtr[nRecPoints].fPHOSModule = list[point].fModule;
	204	nRecPoints++;
	205	}
	206	}
	207
	208
	209
210	energiesList = NULL;
211
212	return nRecPoints;
213
214	}//end CreateRecPointStructArray
215
216
217	/**
218	* Calculating the center of gravity of a rec point
219	* Not working well at this point!
220	* @param recPointPtr pointer to the rec point
221	**/
222	Int_t
223	AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
224	{
91b95d47	225	//Calculate the center of gravity
91b95d47	226
aac22523	227	Float_t xt = 0;
	228	Float_t zt = 0;
	229	Float_t xi = 0;
	230	Float_t zj = 0;
	231	Float_t w = 0;
	232	Float_t w0 = 4.5;
	233	Float_t wtot = 0;
	234
	235	Int_t x = recPointPtr->fCoordinatesPtr[0];
	236	Int_t z = recPointPtr->fCoordinatesPtr[1];
	237
	238	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	239	{
	240	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	241	{
	242	xi = x + i;
	243	zj = z + j;
309b8fca	244	w = TMath::Max( (Float_t)0., (Float_t)(w0 + log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]) ) ) ;
aac22523	245	xt += xi * w ;
	246	zt += zj * w ;
	247	wtot += w ;
	248
	249	}
	250	}
	251	xt /= wtot ;
	252	zt /= wtot ;
	253
	254	recPointPtr->fX = xt;
	255	recPointPtr->fZ = zt;
	256
	257	return 0;
	258
	259	}//end CalculateCenterOfGravity
	260
	261
	262
	263	/**
	264	* Create a simpler data struct for a rec point, containing
	265	* only coordinates, energy and timing
	266	* @param recPointPtr pointer to the rec point
9be2600f	267	* @param clusterStructPtr pointer to the cluster struct
aac22523	268	**/
	269	Int_t
	270	AliHLTPHOSClusterizer::ClusterizeStruct(AliHLTPHOSRecPointDataStruct* recPointPtr, AliHLTPHOSClusterDataStruct* clusterStructPtr)
	271	{
91b95d47	272	//Simplify the rec points
aac22523	273
	274	Float_t clusterEnergy = 0;
	275	Float_t* energiesListPtr = recPointPtr->fEnergiesListPtr;
	276
	277	for(Int_t i = 0; i < recPointPtr->fMultiplicity; i++)
	278	{
	279	clusterEnergy += energiesListPtr[i];
	280	}
	281
	282	clusterStructPtr->fLocalPositionPtr[0] = recPointPtr->fX;
	283	clusterStructPtr->fLocalPositionPtr[1] = recPointPtr->fZ;
	284	clusterStructPtr->fClusterEnergy = clusterEnergy;
	285	clusterStructPtr->fPHOSModule = recPointPtr->fPHOSModule;
	286
	287	return 0;
	288
	289	}//end ClusterizeStruct
	290
	291
	292
	293
	294	/**
	295	* Resets the cell energy array
	296	**/
	297	Int_t
	298	AliHLTPHOSClusterizer::ResetCellEnergyArray()
	299
	300	{
91b95d47	301	//Reset the cell energy array
91b95d47	302
aac22523	303	for(Int_t x = 0; x < N_ROWS_MOD; x++)
	304	{
	305	for(Int_t z = 0; z < N_COLUMNS_MOD; z++)
	306	{
	307	fEnergyArray[x][z] = 0;
	308	}
	309	}
	310
	311	return 0;
	312
	313	}//end ResetCellEnergyArray
	314