[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloFourier.cxx

// $Id: AliHLTPHOSFourier.cxx 34951 2009-09-23 14:35:38Z phille $

/**************************************************************************
 * This file is property of and copyright by the Experimental Nuclear     *
 * Physics Group, Dep. of Physics                                         *
 * University of Oslo, Norway, 2007                                       *
 *                                                                        *
 * Author: Per Thomas Hille <perthi@fys.uio.no> for the ALICE HLT Project.*
 * Contributors are mentioned in the code where appropriate.              *
 * Please report bugs to perthi@fys.uio.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.                  *
 **************************************************************************/
#include "AliHLTCaloFourier.h"

//#include "AliHLTCaloRcuFFTDataStruct.h"

#include  "AliHLTCaloConstants.h"

ClassImp(AliHLTCaloFourier);


AliHLTCaloFourier::AliHLTCaloFourier() :  fFFT_own(0),
					  fFFTInputArray(0),
					  fFFTOutputArray(0),
					  fIsFirstChannel(true),
					  fFixedDataSize(0), 
					  fFFTOupuStruct(),
					  fCurrentEvent(0)
{
  

}


AliHLTCaloFourier::AliHLTCaloFourier(const AliHLTCaloFourier&) : fFFT_own(0),
								 fFFTInputArray(0),
								 fFFTOutputArray(0),
								 fIsFirstChannel(true),
								 fFixedDataSize(0),
								 fFFTOupuStruct(), 
								 fCurrentEvent(0)

{
  
}


AliHLTCaloFourier::~AliHLTCaloFourier()
{

}


AliHLTCaloRcuFFTDataStruct 
AliHLTCaloFourier::GetPSD()
{
  return fFFTOupuStruct;
}


void 
AliHLTCaloFourier::ProcessFourier(const Int_t *data, const int length, const int /*z*/, const int /*x*/, const int gain, const int event)
{
  Double_t  re = 0;
  Double_t  im = 0;

  if( (event > 0 ) && (event != fCurrentEvent ))
    {
      fCurrentEvent = event;
      ResetEventPSD(gain);
    }

  if(fIsFirstChannel == true)
    {
      fCurrentEvent = event;
      fIsFirstChannel = false;
      fFixedDataSize = length;
      Int_t n_size = fFixedDataSize +1;
      fFFT_own = TVirtualFFT::FFT(1, &n_size, "R2C ES K"); 
      Init();
  }

  if( CheckSignal(data, length) == true)
    {
      Int2Double(data, fFFTInputArray,  fFixedDataSize );
      fFFT_own->SetPoints( fFFTInputArray );	 
      fFFT_own->Transform();

      for(int j=0; j < length; j++)
	{
	  fFFT_own->GetPointComplex(j,  re,  im);
	  fFFTOupuStruct.fGlobalAccumulatedPSD[gain][j] +=  EvaluateMagnitude(re, im);
	  fFFTOupuStruct.fGlobalLastPSD[gain][j] +=  EvaluateMagnitude(re, im);

	}
    }
  //  printf("AliHLTCaloFourier::ProcessFourier;  (z, x, gain)  =  (%d, %d, %d),  length = %d",  z,  x,  gain, length);
}

void 
AliHLTCaloFourier::ResetEventPSD(const int gain)
{
  cout << " AliHLTCaloFourier::ResetEventPS, resetting event PSD "<< endl;
  for(int i = 0;  i < fFixedDataSize; i++ )
    {
      fFFTOupuStruct.fGlobalLastPSD[gain][i] = 0;
    }
}


bool 
//AliHLTCaloFourier::CheckSignal(const UInt_t *data, const int length)
AliHLTCaloFourier::CheckSignal(const Int_t *data, const int length)
{
  //  UInt_t tmpMax =  Max(  const_cast< UInt_t *>(data), length);
  //  UInt_t tmpMin =  Min(  const_cast< UInt_t *>(data), length);


  Int_t tmpMax =  Max(  const_cast< Int_t *>(data), length);
  Int_t tmpMin =  Min(  const_cast< Int_t *>(data), length);


  // if( (tmpMax -tmpMin) > 200)
  if( (tmpMax -tmpMin) > 100)
   {
      cout << "FourierAna::CheckSignal min = "<< tmpMin << "  max =  " << tmpMax << endl;
    }
  
  if( (tmpMax >= AliHLTCaloConstants::GetMAXBINVALUE() ) || tmpMin < 1 )
    {
      cout << "ERROR, FourierAna::CheckSignal failed, signal out of range, min= "<< tmpMin << "max = " << tmpMax << endl;
      return false;
    }
  else
    {
      return true;
    }
  
  //return true;
}


void 
AliHLTCaloFourier::Init()
{
  fFFTInputArray = new double[fFixedDataSize];
  fFFTOutputArray = new double[fFixedDataSize];
  
  for(int gain = 0; gain <  AliHLTCaloConstants::GetNGAINS(); gain ++)
    {
      fFFTOupuStruct.fDataLength = fFixedDataSize;
 
      for(int k= 0; k <fFixedDataSize; k++ )
	{
	  fFFTInputArray[k] = 0;
	  fFFTOutputArray[k] = 0;

	}
      for(int i=0; i <  AliHLTCaloConstants::GetALTROMAXSAMPLES()  ; i++)
	{
	  fFFTOupuStruct.fGlobalAccumulatedPSD[gain][i] = 0;
	  fFFTOupuStruct.fGlobalLastPSD[gain][i] = 0;
	}
    }
}


double 
AliHLTCaloFourier::EvaluateMagnitude(const double re, const double im)
{
  return re*re + im*im;
}


void 
AliHLTCaloFourier::Int2Double(const Int_t *inputarray, double *outputarray, const int size)
{
  for(int i=0; i< size; i++)
    {
      outputarray[i] = (double)inputarray[i];
    }
}
Commit	Line	Data
3800a654	1	// $Id: AliHLTPHOSFourier.cxx 34951 2009-09-23 14:35:38Z phille $
	2
	3	/**************************************************************************
	4	* This file is property of and copyright by the Experimental Nuclear *
	5	* Physics Group, Dep. of Physics *
	6	* University of Oslo, Norway, 2007 *
	7	* *
	8	* Author: Per Thomas Hille <perthi@fys.uio.no> for the ALICE HLT Project.*
	9	* Contributors are mentioned in the code where appropriate. *
	10	* Please report bugs to perthi@fys.uio.no *
	11	* *
	12	* Permission to use, copy, modify and distribute this software and its *
	13	* documentation strictly for non-commercial purposes is hereby granted *
	14	* without fee, provided that the above copyright notice appears in all *
	15	* copies and that both the copyright notice and this permission notice *
	16	* appear in the supporting documentation. The authors make no claims *
	17	* about the suitability of this software for any purpose. It is *
	18	* provided "as is" without express or implied warranty. *
	19	**************************************************************************/
	20	#include "AliHLTCaloFourier.h"
	21
	22	//#include "AliHLTCaloRcuFFTDataStruct.h"
	23
76ffc0ce	24	#include "AliHLTCaloConstants.h"
76ffc0ce	25
4f4b7ba4	26	ClassImp(AliHLTCaloFourier);
4f4b7ba4	27
76ffc0ce	28
76ffc0ce	29
3800a654	30	AliHLTCaloFourier::AliHLTCaloFourier() : fFFT_own(0),
	31	fFFTInputArray(0),
	32	fFFTOutputArray(0),
	33	fIsFirstChannel(true),
	34	fFixedDataSize(0),
	35	fFFTOupuStruct(),
	36	fCurrentEvent(0)
	37	{
	38
	39
	40	}
	41
	42
76ffc0ce	43
76ffc0ce	44
3800a654	45	AliHLTCaloFourier::AliHLTCaloFourier(const AliHLTCaloFourier&) : fFFT_own(0),
	46	fFFTInputArray(0),
	47	fFFTOutputArray(0),
	48	fIsFirstChannel(true),
	49	fFixedDataSize(0),
	50	fFFTOupuStruct(),
	51	fCurrentEvent(0)
	52
	53	{
	54
	55	}
	56
	57
76ffc0ce	58
3800a654	59	AliHLTCaloFourier::~AliHLTCaloFourier()
	60	{
	61
	62	}
	63
	64
	65
	66	AliHLTCaloRcuFFTDataStruct
	67	AliHLTCaloFourier::GetPSD()
	68	{
	69	return fFFTOupuStruct;
	70	}
	71
	72
	73
	74	void
	75	AliHLTCaloFourier::ProcessFourier(const Int_t data, const int length, const int /z/, const int /x*/, const int gain, const int event)
	76	{
	77	Double_t re = 0;
	78	Double_t im = 0;
	79
	80	if( (event > 0 ) && (event != fCurrentEvent ))
	81	{
	82	fCurrentEvent = event;
	83	ResetEventPSD(gain);
	84	}
	85
	86	if(fIsFirstChannel == true)
	87	{
	88	fCurrentEvent = event;
	89	fIsFirstChannel = false;
	90	fFixedDataSize = length;
	91	Int_t n_size = fFixedDataSize +1;
	92	fFFT_own = TVirtualFFT::FFT(1, &n_size, "R2C ES K");
	93	Init();
	94	}
	95
	96	if( CheckSignal(data, length) == true)
	97	{
	98	Int2Double(data, fFFTInputArray, fFixedDataSize );
	99	fFFT_own->SetPoints( fFFTInputArray );
	100	fFFT_own->Transform();
	101
	102	for(int j=0; j < length; j++)
	103	{
	104	fFFT_own->GetPointComplex(j, re, im);
	105	fFFTOupuStruct.fGlobalAccumulatedPSD[gain][j] += EvaluateMagnitude(re, im);
	106	fFFTOupuStruct.fGlobalLastPSD[gain][j] += EvaluateMagnitude(re, im);
	107
	108	}
	109	}
	110	// printf("AliHLTCaloFourier::ProcessFourier; (z, x, gain) = (%d, %d, %d), length = %d", z, x, gain, length);
	111	}
	112
	113	void
	114	AliHLTCaloFourier::ResetEventPSD(const int gain)
	115	{
	116	cout << " AliHLTCaloFourier::ResetEventPS, resetting event PSD "<< endl;
	117	for(int i = 0; i < fFixedDataSize; i++ )
	118	{
	119	fFFTOupuStruct.fGlobalLastPSD[gain][i] = 0;
	120	}
	121	}
	122
123
124	bool
125	//AliHLTCaloFourier::CheckSignal(const UInt_t *data, const int length)
126	AliHLTCaloFourier::CheckSignal(const Int_t *data, const int length)
127	{
128	// UInt_t tmpMax = Max( const_cast< UInt_t *>(data), length);
129	// UInt_t tmpMin = Min( const_cast< UInt_t *>(data), length);
130
131
132	Int_t tmpMax = Max( const_cast< Int_t *>(data), length);
133	Int_t tmpMin = Min( const_cast< Int_t *>(data), length);
134
135
136	// if( (tmpMax -tmpMin) > 200)
137	if( (tmpMax -tmpMin) > 100)
138	{
139	cout << "FourierAna::CheckSignal min = "<< tmpMin << " max = " << tmpMax << endl;
140	}
141
76ffc0ce	142	if( (tmpMax >= AliHLTCaloConstants::GetMAXBINVALUE() ) \|\| tmpMin < 1 )
3800a654	143	{
	144	cout << "ERROR, FourierAna::CheckSignal failed, signal out of range, min= "<< tmpMin << "max = " << tmpMax << endl;
	145	return false;
	146	}
	147	else
	148	{
	149	return true;
	150	}
	151
	152	//return true;
	153	}
	154
	155
	156	void
	157	AliHLTCaloFourier::Init()
	158	{
	159	fFFTInputArray = new double[fFixedDataSize];
	160	fFFTOutputArray = new double[fFixedDataSize];
	161
76ffc0ce	162	for(int gain = 0; gain < AliHLTCaloConstants::GetNGAINS(); gain ++)
3800a654	163	{
	164	fFFTOupuStruct.fDataLength = fFixedDataSize;
	165
	166	for(int k= 0; k <fFixedDataSize; k++ )
	167	{
	168	fFFTInputArray[k] = 0;
	169	fFFTOutputArray[k] = 0;
	170
	171	}
76ffc0ce	172	for(int i=0; i < AliHLTCaloConstants::GetALTROMAXSAMPLES() ; i++)
3800a654	173	{
	174	fFFTOupuStruct.fGlobalAccumulatedPSD[gain][i] = 0;
	175	fFFTOupuStruct.fGlobalLastPSD[gain][i] = 0;
	176	}
	177	}
	178	}
	179
	180
	181	double
	182	AliHLTCaloFourier::EvaluateMagnitude(const double re, const double im)
	183	{
	184	return rere + imim;
	185	}
	186
	187
	188	void
	189	AliHLTCaloFourier::Int2Double(const Int_t inputarray, double outputarray, const int size)
	190	{
	191	for(int i=0; i< size; i++)
	192	{
	193	outputarray[i] = (double)inputarray[i];
	194	}
	195	}