[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_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 
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_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(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_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
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
	91	zMod = z + (cellData->fRcuZ)*N_ROWS_RCU;
	92	xMod = x + (cellData->fRcuX)*N_COLUMNS_RCU;
	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
	168	if((z-fArraySize/2) < 0/= - N_ROWS_MOD/2/ \|\| (z+fArraySize/2) >= N_ROWS_MOD//2/)
	169	{
	170	edgeFlagRows = 1;
	171	continue;
	172	}
	173	if((x-fArraySize/2) < 0/= - N_COLUMNS_MOD/2/ \|\| (x+fArraySize/2) >= N_COLUMNS_MOD)
	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	}
	218
6e709a0d	219	if(energiesList)
	220	{
	221	delete [] energiesList;
	222	energiesList = NULL;
	223	}
aac22523	224
	225	return nRecPoints;
	226
	227	}//end CreateRecPointStructArray
	228
	229
	230	/**
	231	* Calculating the center of gravity of a rec point
	232	* Not working well at this point!
	233	* @param recPointPtr pointer to the rec point
	234	**/
6e709a0d	235	int
aac22523	236	AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
aac22523	237	{
6e709a0d	238	//CalculateCenterOfGravity
6e709a0d	239	//Copied from offline code, and modified
91b95d47	240
aac22523	241	Float_t xt = 0;
	242	Float_t zt = 0;
	243	Float_t xi = 0;
	244	Float_t zj = 0;
	245	Float_t w = 0;
	246	Float_t w0 = 4.5;
	247	Float_t wtot = 0;
	248
	249	Int_t x = recPointPtr->fCoordinatesPtr[0];
	250	Int_t z = recPointPtr->fCoordinatesPtr[1];
	251
	252	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	253	{
	254	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	255	{
	256	xi = x + i;
	257	zj = z + j;
309b8fca	258	w = TMath::Max( (Float_t)0., (Float_t)(w0 + log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]) ) ) ;
6e709a0d	259	//printf("log in cog: %f\n", log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]));
aac22523	260	xt += xi * w ;
aac22523	261	zt += zj * w ;
6e709a0d	262	wtot += w ;
aac22523	263	}
	264	}
	265	xt /= wtot ;
	266	zt /= wtot ;
	267
6e709a0d	268	//printf("wtot cog: %f\n", wtot);
6e709a0d	269
aac22523	270	recPointPtr->fX = xt;
	271	recPointPtr->fZ = zt;
	272
	273	return 0;
	274
	275	}//end CalculateCenterOfGravity
	276
6e709a0d	277	int
	278	AliHLTPHOSClusterizer::CalculateMoments(AliHLTPHOSRecPointDataStruct* recPointPtr, Bool_t axisOnly)
	279	{
	280	//Copied from offline code, and modified
	281
	282
	283	// Calculate the shower moments in the eigen reference system
	284	// M2x, M2z, M3x, M4z
	285	// Calculate the angle between the shower position vector and the eigen vector
	286
	287	Double_t wtot = 0. ;
	288	Double_t x = 0.;
	289	Double_t z = 0.;
	290	Double_t dxx = 0.;
	291	Double_t dzz = 0.;
	292	Double_t dxz = 0.;
	293	Double_t lambda0=0, lambda1=0;
	294	Float_t logWeight = 4.5;
	295
	296	//Int_t iDigit;
	297
	298	// 1) Find covariance matrix elements:
	299	// \|\| dxx dxz \|\|
	300	// \|\| dxz dzz \|\|
	301
	302	Int_t xi;
	303	Int_t zj;
	304	Double_t w;
	305
	306	Int_t xc = recPointPtr->fCoordinatesPtr[0];
	307	Int_t zc = recPointPtr->fCoordinatesPtr[1];
	308
	309	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
	310	{
	311	xi = xc + i;
	312	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
	313	{
	314	zj = zc + j;
	315	if (fEnergyArray[xi][zj] > 0)
	316	{
	317	w = TMath::Max(0.,logWeight + log( fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
	318	//printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
	319	x += w * xi ;
	320	z += w * zj ;
	321	dxx += w * xi * xi ;
	322	dzz += w * zj * zj ;
	323	dxz += w * xi * zj ;
	324	wtot += w ;
	325	}
	326	}
	327	}
	328	//printf("wtot: %f\n", wtot);
	329	if (wtot>0)
	330	{
	331	x /= wtot ;
	332	z /= wtot ;
	333	dxx /= wtot ;
	334	dzz /= wtot ;
	335	dxz /= wtot ;
	336	dxx -= x * x ;
	337	dzz -= z * z ;
	338	dxz -= x * z ;
	339
	340	// 2) Find covariance matrix eigen values lambda0 and lambda1
341
342	lambda0 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
343	lambda1 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ;
344	}
345
346	recPointPtr->fM2x = lambda0;
347	recPointPtr->fM2z = lambda1;
348	// printf("lambda0 = %f -- lambda1 = %f\n", lambda0, lambda1);
349	// 3) Find covariance matrix eigen vectors e0 and e1
350	if(!axisOnly)
351	{
352	TVector2 e0,e1;
353	if (dxz != 0)
354	e0.Set(1.,(lambda0-dxx)/dxz);
355	else
356	e0.Set(0.,1.);
357
358	e0 = e0.Unit();
359	e1.Set(-e0.Y(),e0.X());
360
361	// 4) Rotate cluster tensor from (x,z) to (e0,e1) system
362	// and calculate moments M3x and M4z
363
364	Float_t cosPhi = e0.X();
365	Float_t sinPhi = e0.Y();
366
367	Float_t xiPHOS ;
368	Float_t zjPHOS ;
369	Double_t dx3, dz3, dz4;
370	wtot = 0.;
371	x = 0.;
372	z = 0.;
373	dxx = 0.;
374	dzz = 0.;
375	dxz = 0.;
376	dx3 = 0.;
377	dz3 = 0.;
378	dz4 = 0.;
379	for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
380	{
381	for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
382	{
383	xi = xc + i;
384	zj = zc + j;
385	xiPHOS = (Float_t)xicosPhi + (Float_t)zjsinPhi;
386	zjPHOS = (Float_t)zjcosPhi - (Float_t)xisinPhi;
387	// xi = (Float_t)xicosPhi + (Float_t)zjsinPhi;
388	// zj = (Float_t)zjcosPhi - (Float_t)xisinPhi;
389	if (fEnergyArray[xi][zj] > 0)
390	{
391	w = TMath::Max(0.,logWeight+TMath::Log(fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
392	// printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
393
394	x += w * xi ;
395	z += w * zj ;
396	dxx += w * xi * xi ;
397	dzz += w * zj * zj ;
398	dxz += w * xi * zj ;
399	dx3 += w * xi * xi * xi;
400	dz3 += w * zj * zj * zj ;
401	dz4 += w * zj * zj * zj * zj ;
402	wtot += w ;
403	}
404	}
405	}
406	if (wtot>0)
407	{
408	x /= wtot ;
409	z /= wtot ;
410	dxx /= wtot ;
411	dzz /= wtot ;
412	dxz /= wtot ;
413	dx3 /= wtot ;
414	dz3 /= wtot ;
415	dz4 /= wtot ;
416	dx3 += -3dxxx + 2xx*x;
417	dz4 += -4dz3z + 6dzzzz -3zzz*z;
418	dxx -= x * x ;
419	dzz -= z * z ;
420	dxz -= x * z ;
421	}
422
423	// 5) Find an angle between cluster center vector and eigen vector e0
424
425	Float_t phi = 0;
426	if ( (xx+zz) > 0 )
427	phi = TMath::ACos ((xe0.X() + ze0.Y()) / sqrt(xx + zz));
428
429	recPointPtr->fM3x = dx3;
430	recPointPtr->fM4z = dz4;
431	recPointPtr->fPhixe = phi;
432	// printf("%f\n",phi);
433	}
434	return 0;
435	}
aac22523	436
	437
	438	/**
	439	* Create a simpler data struct for a rec point, containing
	440	* only coordinates, energy and timing
	441	* @param recPointPtr pointer to the rec point
6e709a0d	442	* @clusterStructPtr pointer to the cluster struct
aac22523	443	**/
6e709a0d	444	int
aac22523	445	AliHLTPHOSClusterizer::ClusterizeStruct(AliHLTPHOSRecPointDataStruct* recPointPtr, AliHLTPHOSClusterDataStruct* clusterStructPtr)
aac22523	446	{
6e709a0d	447	//ClusterizeStruct
aac22523	448
	449	Float_t clusterEnergy = 0;
	450	Float_t* energiesListPtr = recPointPtr->fEnergiesListPtr;
	451
	452	for(Int_t i = 0; i < recPointPtr->fMultiplicity; i++)
	453	{
	454	clusterEnergy += energiesListPtr[i];
	455	}
	456
	457	clusterStructPtr->fLocalPositionPtr[0] = recPointPtr->fX;
	458	clusterStructPtr->fLocalPositionPtr[1] = recPointPtr->fZ;
	459	clusterStructPtr->fClusterEnergy = clusterEnergy;
	460	clusterStructPtr->fPHOSModule = recPointPtr->fPHOSModule;
	461
	462	return 0;
	463
	464	}//end ClusterizeStruct
	465
	466
	467
	468
	469	/**
	470	* Resets the cell energy array
	471	**/
6e709a0d	472	int
aac22523	473	AliHLTPHOSClusterizer::ResetCellEnergyArray()
	474
	475	{
6e709a0d	476	//ResetCellEnergyArray
91b95d47	477
aac22523	478	for(Int_t x = 0; x < N_ROWS_MOD; x++)
	479	{
	480	for(Int_t z = 0; z < N_COLUMNS_MOD; z++)
	481	{
	482	fEnergyArray[x][z] = 0;
	483	}
	484	}
	485
	486	return 0;
	487
	488	}//end ResetCellEnergyArray
	489