[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"
#include "AliHLTPHOSConstants.h"

using namespace PhosHLTConst;


ClassImp(AliHLTPHOSClusterizer);

/**
* Main constructor
**/
AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():fPHOSModule(0), fThreshold(0), fClusterThreshold(0), 
					       fHighGainFactor(0.005), fLowGainFactor(0.08),
					       fArraySize(3), fMultiplicity(fArraySize*fArraySize)
{
  //See header file for documentation
  
}//end

AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &):fPHOSModule(0), 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
AliHLTPHOSClusterizer::BuildCellEnergyArray(AliHLTPHOSRcuCellEnergyDataStruct* cellData, 
					    AliHLTPHOSRecPointListDataStruct* recPointList)
{
  //BuildCellEnergyArray

  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_ZROWS_RCU;
      xMod = x + (cellData->fRcuX)*N_XCOLUMNS_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 
AliHLTPHOSClusterizer::CreateRecPointStructArray(AliHLTPHOSRecPointDataStruct* recPointStructArrayPtr, 
						 AliHLTPHOSRecPointListDataStruct* list, 
						 int nPoints) 

{
  //CreateRecPointStructArray

  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;
  Int_t i = 0;
  Int_t j = 0;
  Int_t a = 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_ZROWS_MOD/2*/ || (z+fArraySize/2) >= N_ZROWS_MOD/*/2*/) 
	{
	  edgeFlagRows = 1;
	  continue;
	}
      if((x-fArraySize/2) < 0/*= - N_COLUMNS_MOD/2*/ || (x+fArraySize/2) >= N_XCOLUMNS_MOD) 
	{
	  edgeFlagCols = 1;
	  continue;
	}
      

      if(!flag) energiesList = new Float_t[fMultiplicity];

      flag = 0;
      k = 0;	      
	  
      for(i = -fArraySize/2; i <= fArraySize/2; i++)
	{
	  if(flag) break;
	  for(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)
	{
	  for(a = 0; a < k; a++) 
	    {
	      (recPointStructArrayPtr[nRecPoints].fEnergiesListPtr)[a] = energiesList[a];
	    }
	  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++;
	}
    }
  if(energiesList)
    {
      delete [] energiesList;
      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
AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
{
  //CalculateCenterOfGravity
  //Copied from offline code, and modified

  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]) ) ) ;
	  //printf("log in cog: %f\n", log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]));
	  xt += xi * w ;
	  zt += zj * w ;
	  wtot += w ;	  
	}
    }
  xt /= wtot ;
  zt /= wtot ;
  
  //printf("wtot cog: %f\n", wtot);

  recPointPtr->fX = xt;
  recPointPtr->fZ = zt;
  
  return 0;

}//end CalculateCenterOfGravity
  
int
AliHLTPHOSClusterizer::CalculateMoments(AliHLTPHOSRecPointDataStruct* recPointPtr, Bool_t axisOnly)
{
  //Copied from offline code, and modified
  
  
  // Calculate the shower moments in the eigen reference system
  // M2x, M2z, M3x, M4z
  // Calculate the angle between the shower position vector and the eigen vector
  
  Double_t wtot = 0. ;
  Double_t x    = 0.;
  Double_t z    = 0.;
  Double_t dxx  = 0.;
  Double_t dzz  = 0.;
  Double_t dxz  = 0.;
  Double_t lambda0=0, lambda1=0;
  Float_t logWeight = 4.5;
  
  //Int_t iDigit;
  
  // 1) Find covariance matrix elements:
  //    || dxx dxz ||
  //    || dxz dzz ||
  
  Int_t xi;
  Int_t zj;
  Double_t w;
  
  Int_t xc = recPointPtr->fCoordinatesPtr[0];
  Int_t zc = recPointPtr->fCoordinatesPtr[1];
  
  for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
    {
      xi = xc + i;
      for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	{
	  zj = zc + j;
	  if (fEnergyArray[xi][zj] > 0) 
	    {
	      w     = TMath::Max(0.,logWeight + log( fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
	      //printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
	      x    += w * xi ;
	      z    += w * zj ; 
	      dxx  += w * xi * xi ;
	      dzz  += w * zj * zj ;
	      dxz  += w * xi * zj ; 
	      wtot += w ;
	    }
	}
    }
  //printf("wtot: %f\n", wtot);
  if (wtot>0) 
    {
      x   /= wtot ;
      z   /= wtot ;
      dxx /= wtot ;
      dzz /= wtot ;
      dxz /= wtot ;
      dxx -= x * x ;
      dzz -= z * z ;
      dxz -= x * z ;
      
  // 2) Find covariance matrix eigen values lambda0 and lambda1

      lambda0 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
      lambda1 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
    }
    
  recPointPtr->fM2x   = lambda0;
  recPointPtr->fM2z   = lambda1;
  //  printf("lambda0 = %f -- lambda1 = %f\n", lambda0, lambda1);
  // 3) Find covariance matrix eigen vectors e0 and e1
  if(!axisOnly)
    {
      TVector2 e0,e1;
      if (dxz != 0)
	e0.Set(1.,(lambda0-dxx)/dxz);
      else
	e0.Set(0.,1.);
      
      e0 = e0.Unit();
      e1.Set(-e0.Y(),e0.X());
      
      // 4) Rotate cluster tensor from (x,z) to (e0,e1) system
      //    and calculate moments M3x and M4z
      
      Float_t cosPhi = e0.X();
      Float_t sinPhi = e0.Y();
      
      Float_t xiPHOS ;
      Float_t zjPHOS ;
      Double_t dx3, dz3, dz4;
      wtot = 0.;
      x    = 0.;
      z    = 0.;
      dxx  = 0.;
      dzz  = 0.;
      dxz  = 0.;
      dx3  = 0.;
      dz3  = 0.;
      dz4  = 0.;
      for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	{
	  for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	    {
	      xi = xc + i;
	      zj = zc + j;
	      xiPHOS = (Float_t)xi*cosPhi + (Float_t)zj*sinPhi;
	      zjPHOS = (Float_t)zj*cosPhi - (Float_t)xi*sinPhi;
	      //	      xi = (Float_t)xi*cosPhi + (Float_t)zj*sinPhi;
	      //	      zj = (Float_t)zj*cosPhi - (Float_t)xi*sinPhi;
	      if (fEnergyArray[xi][zj] > 0) 
		{
		  w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
		  //  printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
			  
		  x    += w * xi ;
		  z    += w * zj ; 
		  dxx  += w * xi * xi ;
		  dzz  += w * zj * zj ;
		  dxz  += w * xi * zj ; 
		  dx3  += w * xi * xi * xi;
		  dz3  += w * zj * zj * zj ;
		  dz4  += w * zj * zj * zj * zj ;
		  wtot += w ;
		}
	    }
	}
      if (wtot>0) 
	{
	  x   /= wtot ;
	  z   /= wtot ;
	  dxx /= wtot ;
	  dzz /= wtot ;
	  dxz /= wtot ;
	  dx3 /= wtot ;
	  dz3 /= wtot ;
	  dz4 /= wtot ;
	  dx3 += -3*dxx*x + 2*x*x*x;
	  dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z;
	  dxx -= x * x ;
	  dzz -= z * z ;
	  dxz -= x * z ;
	}
  
      // 5) Find an angle between cluster center vector and eigen vector e0
      
      Float_t phi = 0;
      if ( (x*x+z*z) > 0 ) 
	phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z));

      recPointPtr->fM3x   = dx3;
      recPointPtr->fM4z   = dz4;
      recPointPtr->fPhixe = phi;
      //      printf("%f\n",phi);
    }
  return 0;
}


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

  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 
AliHLTPHOSClusterizer::ResetCellEnergyArray()

{
  //ResetCellEnergyArray

  for(Int_t x = 0; x < N_ZROWS_MOD; x++)
    {
      for(Int_t z = 0; z < N_XCOLUMNS_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
6e709a0d	15
aac22523	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"
2ef3c547	34	#include "AliHLTPHOSConstants.h"
	35
	36	using namespace PhosHLTConst;
aac22523	37
	38
	39	ClassImp(AliHLTPHOSClusterizer);
	40
	41	/**
	42	* Main constructor
	43	**/
6e709a0d	44	AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():fPHOSModule(0), fThreshold(0), fClusterThreshold(0),
91b95d47	45	fHighGainFactor(0.005), fLowGainFactor(0.08),
aac22523	46	fArraySize(3), fMultiplicity(fArraySize*fArraySize)
aac22523	47	{
6e709a0d	48	//See header file for documentation
aac22523	49
	50	}//end
	51
6e709a0d	52	AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &):fPHOSModule(0), fThreshold(0), fClusterThreshold(0),
91b95d47	53	fHighGainFactor(0.005), fLowGainFactor(0.08),
aac22523	54	fArraySize(3), fMultiplicity(fArraySize*fArraySize)
aac22523	55	{
91b95d47	56	//Copy constructor, not implemented
aac22523	57	}//end
	58
	59	AliHLTPHOSClusterizer:: ~AliHLTPHOSClusterizer()
	60	{
91b95d47	61	//Destructor
aac22523	62	}
	63
	64	/**
	65	* Building a 2D array of cell energies of the PHOS detector
	66	* @param cellData object containing the cell energies from one event
	67	* @param recPointList list to be filled with coordinates of local maxima in the detector
	68	**/
	69
6e709a0d	70	int
aac22523	71	AliHLTPHOSClusterizer::BuildCellEnergyArray(AliHLTPHOSRcuCellEnergyDataStruct* cellData,
	72	AliHLTPHOSRecPointListDataStruct* recPointList)
	73	{
6e709a0d	74	//BuildCellEnergyArray
aac22523	75
	76	Int_t x = 0;
	77	Int_t z = 0;
	78	Int_t gain = 0;
	79	Int_t xMod = 0;
	80	Int_t zMod = 0;
	81	Int_t index = 0;
	82	Double_t energyCount = 0;
	83
	84	for(Int_t cell = 0; cell < cellData->fCnt; cell++)
	85	{
	86
	87	z = (cellData->fValidData[cell]).fZ;
	88	x = (cellData->fValidData[cell]).fX;
	89	gain = (cellData->fValidData[cell]).fGain;
	90
442af5b7	91	zMod = z + (cellData->fRcuZ)*N_ZROWS_RCU;
442af5b7	92	xMod = x + (cellData->fRcuX)*N_XCOLUMNS_RCU;
aac22523	93
	94	energyCount = (cellData->fValidData[cell]).fEnergy;
	95
	96	if(gain == 0 && energyCount < 1023)
	97	{
	98	fEnergyArray[xMod][zMod] = fHighGainFactor * energyCount;
	99
	100	if(fEnergyArray[xMod][zMod] > fClusterThreshold)
	101	{
	102	recPointList[index].fZ = zMod;
	103	recPointList[index].fX = xMod;
	104
	105	for(Int_t j = 0; j < index; j++)
	106	{
	107	if(recPointList[j].fZ == zMod && recPointList[j].fX == xMod)
	108	recPointList[j].fZ = -1;
	109	}
	110	index++;
	111	}
	112
	113	}
	114	else if(gain == 1 && fEnergyArray[xMod][zMod] == 0)
	115	{
	116	fEnergyArray[xMod][zMod] = fLowGainFactor * energyCount;
	117	if(fEnergyArray[xMod][zMod] > fClusterThreshold)
	118	{
	119	recPointList[index].fZ = zMod;
	120	recPointList[index].fX = xMod;
	121	recPointList[index].fModule = cellData->fModuleID;
	122	index++;
	123	}
	124	}
	125
	126	}
	127
	128	fPHOSModule = cellData->fModuleID;
	129
	130	return index;
	131	}//end BuildCellEnergyArray
	132
	133
	134
	135	/**
	136	* Creating an array of rec points
	137	* @param recPointStructArrayPtr array to store the rec points
	138	* @param list list of rec points
	139	* @param nPoints number of points
	140	**/
6e709a0d	141	int
aac22523	142	AliHLTPHOSClusterizer::CreateRecPointStructArray(AliHLTPHOSRecPointDataStruct* recPointStructArrayPtr,
aac22523	143	AliHLTPHOSRecPointListDataStruct* list,
6e709a0d	144	int nPoints)
aac22523	145
aac22523	146	{
6e709a0d	147	//CreateRecPointStructArray
91b95d47	148
aac22523	149	Int_t flag = 0;
	150	Int_t edgeFlagRows = 0;
	151	Int_t edgeFlagCols = 0;
	152	Int_t k = 0;
	153	Int_t nRecPoints = 0;
	154	Int_t z = 0;
	155	Int_t x = 0;
6e709a0d	156	Int_t i = 0;
	157	Int_t j = 0;
	158	Int_t a = 0;
aac22523	159
aac22523	160	Float_t* energiesList = NULL;
6e709a0d	161
aac22523	162	for(Int_t point = 0; point < nPoints; point++)
	163	{
	164	z = list[point].fZ;
	165	x = list[point].fX;
	166	if(z == -1) continue;
	167
442af5b7	168	if((z-fArraySize/2) < 0/= - N_ZROWS_MOD/2/ \|\| (z+fArraySize/2) >= N_ZROWS_MOD//2/)
aac22523	169	{
	170	edgeFlagRows = 1;
	171	continue;
	172	}
442af5b7	173	if((x-fArraySize/2) < 0/= - N_COLUMNS_MOD/2/ \|\| (x+fArraySize/2) >= N_XCOLUMNS_MOD)
aac22523	174	{
	175	edgeFlagCols = 1;
	176	continue;
	177	}
	178
	179
	180	if(!flag) energiesList = new Float_t[fMultiplicity];
6e709a0d	181
aac22523	182	flag = 0;
	183	k = 0;
	184
6e709a0d	185	for(i = -fArraySize/2; i <= fArraySize/2; i++)
aac22523	186	{
aac22523	187	if(flag) break;
6e709a0d	188	for(j = -fArraySize/2; j <= fArraySize/2; j++)
aac22523	189	{
	190
	191	if(fEnergyArray[x+i][z+j] > fEnergyArray[x][z] && abs(i) < 2 && abs(j) < 2)
	192	{
	193	flag = 1;
	194	break;
	195	}
6e709a0d	196
aac22523	197	energiesList[k] = fEnergyArray[x+i][z+j];
	198	k++;
	199	}
	200	}
	201
	202
	203	if(!flag && k)
	204	{
6e709a0d	205	for(a = 0; a < k; a++)
	206	{
	207	(recPointStructArrayPtr[nRecPoints].fEnergiesListPtr)[a] = energiesList[a];
	208	}
aac22523	209	recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[0] = x;
	210	recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[1] = z;
	211	recPointStructArrayPtr[nRecPoints].fX = x;
	212	recPointStructArrayPtr[nRecPoints].fZ = z;
6e709a0d	213	// recPointStructArrayPtr[nRecPoints].fMultiplicity = fMultiplicity;
aac22523	214	recPointStructArrayPtr[nRecPoints].fPHOSModule = list[point].fModule;
	215	nRecPoints++;
	216	}
	217	}
6e709a0d	218	if(energiesList)
	219	{
	220	delete [] energiesList;
	221	energiesList = NULL;
	222	}
aac22523	223
	224	return nRecPoints;
	225
	226	}//end CreateRecPointStructArray
	227
	228
	229	/**
	230	* Calculating the center of gravity of a rec point
	231	* Not working well at this point!
	232	* @param recPointPtr pointer to the rec point
	233	**/
6e709a0d	234	int
aac22523	235	AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
aac22523	236	{
6e709a0d	237	//CalculateCenterOfGravity
6e709a0d	238	//Copied from offline code, and modified
91b95d47	239
aac22523	240	Float_t xt = 0;
	241	Float_t zt = 0;
	242	Float_t xi = 0;
	243	Float_t zj = 0;
	244	Float_t w = 0;
	245	Float_t w0 = 4.5;
	246	Float_t wtot = 0;
	247
	248	Int_t x = recPointPtr->fCoordinatesPtr[0];
	249	Int_t z = recPointPtr->fCoordinatesPtr[1];
	250
	251	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	252	{
	253	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	254	{
	255	xi = x + i;
	256	zj = z + j;
309b8fca	257	w = TMath::Max( (Float_t)0., (Float_t)(w0 + log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]) ) ) ;
6e709a0d	258	//printf("log in cog: %f\n", log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]));
aac22523	259	xt += xi * w ;
aac22523	260	zt += zj * w ;
6e709a0d	261	wtot += w ;
aac22523	262	}
	263	}
	264	xt /= wtot ;
	265	zt /= wtot ;
	266
6e709a0d	267	//printf("wtot cog: %f\n", wtot);
6e709a0d	268
aac22523	269	recPointPtr->fX = xt;
	270	recPointPtr->fZ = zt;
	271
	272	return 0;
	273
	274	}//end CalculateCenterOfGravity
	275
6e709a0d	276	int
	277	AliHLTPHOSClusterizer::CalculateMoments(AliHLTPHOSRecPointDataStruct* recPointPtr, Bool_t axisOnly)
	278	{
	279	//Copied from offline code, and modified
	280
	281
	282	// Calculate the shower moments in the eigen reference system
	283	// M2x, M2z, M3x, M4z
	284	// Calculate the angle between the shower position vector and the eigen vector
	285
	286	Double_t wtot = 0. ;
	287	Double_t x = 0.;
	288	Double_t z = 0.;
	289	Double_t dxx = 0.;
	290	Double_t dzz = 0.;
	291	Double_t dxz = 0.;
	292	Double_t lambda0=0, lambda1=0;
	293	Float_t logWeight = 4.5;
	294
	295	//Int_t iDigit;
	296
	297	// 1) Find covariance matrix elements:
	298	// \|\| dxx dxz \|\|
	299	// \|\| dxz dzz \|\|
	300
	301	Int_t xi;
	302	Int_t zj;
	303	Double_t w;
	304
	305	Int_t xc = recPointPtr->fCoordinatesPtr[0];
	306	Int_t zc = recPointPtr->fCoordinatesPtr[1];
	307
	308	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	309	{
	310	xi = xc + i;
	311	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	312	{
	313	zj = zc + j;
	314	if (fEnergyArray[xi][zj] > 0)
	315	{
	316	w = TMath::Max(0.,logWeight + log( fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
	317	//printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
	318	x += w * xi ;
	319	z += w * zj ;
	320	dxx += w * xi * xi ;
	321	dzz += w * zj * zj ;
	322	dxz += w * xi * zj ;
	323	wtot += w ;
	324	}
	325	}
	326	}
	327	//printf("wtot: %f\n", wtot);
	328	if (wtot>0)
	329	{
	330	x /= wtot ;
	331	z /= wtot ;
	332	dxx /= wtot ;
	333	dzz /= wtot ;
	334	dxz /= wtot ;
	335	dxx -= x * x ;
	336	dzz -= z * z ;
	337	dxz -= x * z ;
	338
	339	// 2) Find covariance matrix eigen values lambda0 and lambda1
340
341	lambda0 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
342	lambda1 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
343	}
344
345	recPointPtr->fM2x = lambda0;
346	recPointPtr->fM2z = lambda1;
347	// printf("lambda0 = %f -- lambda1 = %f\n", lambda0, lambda1);
348	// 3) Find covariance matrix eigen vectors e0 and e1
349	if(!axisOnly)
350	{
351	TVector2 e0,e1;
352	if (dxz != 0)
353	e0.Set(1.,(lambda0-dxx)/dxz);
354	else
355	e0.Set(0.,1.);
356
357	e0 = e0.Unit();
358	e1.Set(-e0.Y(),e0.X());
359
360	// 4) Rotate cluster tensor from (x,z) to (e0,e1) system
361	// and calculate moments M3x and M4z
362
363	Float_t cosPhi = e0.X();
364	Float_t sinPhi = e0.Y();
365
366	Float_t xiPHOS ;
367	Float_t zjPHOS ;
368	Double_t dx3, dz3, dz4;
369	wtot = 0.;
370	x = 0.;
371	z = 0.;
372	dxx = 0.;
373	dzz = 0.;
374	dxz = 0.;
375	dx3 = 0.;
376	dz3 = 0.;
377	dz4 = 0.;
378	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
379	{
380	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
381	{
382	xi = xc + i;
383	zj = zc + j;
384	xiPHOS = (Float_t)xicosPhi + (Float_t)zjsinPhi;
385	zjPHOS = (Float_t)zjcosPhi - (Float_t)xisinPhi;
386	// xi = (Float_t)xicosPhi + (Float_t)zjsinPhi;
387	// zj = (Float_t)zjcosPhi - (Float_t)xisinPhi;
388	if (fEnergyArray[xi][zj] > 0)
389	{
390	w = TMath::Max(0.,logWeight+TMath::Log(fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
391	// printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
392
393	x += w * xi ;
394	z += w * zj ;
395	dxx += w * xi * xi ;
396	dzz += w * zj * zj ;
397	dxz += w * xi * zj ;
398	dx3 += w * xi * xi * xi;
399	dz3 += w * zj * zj * zj ;
400	dz4 += w * zj * zj * zj * zj ;
401	wtot += w ;
402	}
403	}
404	}
405	if (wtot>0)
406	{
407	x /= wtot ;
408	z /= wtot ;
409	dxx /= wtot ;
410	dzz /= wtot ;
411	dxz /= wtot ;
412	dx3 /= wtot ;
413	dz3 /= wtot ;
414	dz4 /= wtot ;
415	dx3 += -3dxxx + 2xx*x;
416	dz4 += -4dz3z + 6dzzzz -3zzz*z;
417	dxx -= x * x ;
418	dzz -= z * z ;
419	dxz -= x * z ;
420	}
421
422	// 5) Find an angle between cluster center vector and eigen vector e0
423
424	Float_t phi = 0;
425	if ( (xx+zz) > 0 )
426	phi = TMath::ACos ((xe0.X() + ze0.Y()) / sqrt(xx + zz));
427
428	recPointPtr->fM3x = dx3;
429	recPointPtr->fM4z = dz4;
430	recPointPtr->fPhixe = phi;
431	// printf("%f\n",phi);
432	}
433	return 0;
434	}
aac22523	435
	436
	437	/**
	438	* Create a simpler data struct for a rec point, containing
	439	* only coordinates, energy and timing
	440	* @param recPointPtr pointer to the rec point
6e709a0d	441	* @clusterStructPtr pointer to the cluster struct
aac22523	442	**/
6e709a0d	443	int
aac22523	444	AliHLTPHOSClusterizer::ClusterizeStruct(AliHLTPHOSRecPointDataStruct* recPointPtr, AliHLTPHOSClusterDataStruct* clusterStructPtr)
aac22523	445	{
6e709a0d	446	//ClusterizeStruct
aac22523	447
	448	Float_t clusterEnergy = 0;
	449	Float_t* energiesListPtr = recPointPtr->fEnergiesListPtr;
	450
	451	for(Int_t i = 0; i < recPointPtr->fMultiplicity; i++)
	452	{
	453	clusterEnergy += energiesListPtr[i];
	454	}
	455
	456	clusterStructPtr->fLocalPositionPtr[0] = recPointPtr->fX;
	457	clusterStructPtr->fLocalPositionPtr[1] = recPointPtr->fZ;
	458	clusterStructPtr->fClusterEnergy = clusterEnergy;
	459	clusterStructPtr->fPHOSModule = recPointPtr->fPHOSModule;
	460
	461	return 0;
	462
	463	}//end ClusterizeStruct
	464
	465
	466
	467
	468	/**
	469	* Resets the cell energy array
	470	**/
6e709a0d	471	int
aac22523	472	AliHLTPHOSClusterizer::ResetCellEnergyArray()
	473
	474	{
6e709a0d	475	//ResetCellEnergyArray
91b95d47	476
442af5b7	477	for(Int_t x = 0; x < N_ZROWS_MOD; x++)
aac22523	478	{
442af5b7	479	for(Int_t z = 0; z < N_XCOLUMNS_MOD; z++)
aac22523	480	{
	481	fEnergyArray[x][z] = 0;
	482	}
	483	}
	484
	485	return 0;
	486
	487	}//end ResetCellEnergyArray
	488