[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerPeakFinderV2.cxx

/**************************************************************************
 * This file is property of and copyright by                              *
 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
 *                                                                        *
 * Primary Author: Per Thomas Hille  <perthomas.hille@yale.edu>           *
 *                                                                        *
 * Contributors are mentioned in the code where appropriate.              *
 * Please report bugs to p.t.hille@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.                  *
 **************************************************************************/

// The Peak-Finder algorithm
// The amplitude is extracted  as a
// weighted sum of the samples using the 
// best possible weights.
// The wights is calculated only once and the
// Actual extraction of amplitude and peak position
// Is done with a simple vector multiplication, allowing for
// Extreemely fast computations. 

#include "AliCaloRawAnalyzerPeakFinderV2.h"
#include "AliCaloBunchInfo.h"
#include "AliCaloFitResults.h"
#include <iostream>
#include "unistd.h"
#include "TMath.h"
#include "AliLog.h"

using namespace std;

ClassImp( AliCaloRawAnalyzerPeakFinderV2 )

AliCaloRawAnalyzerPeakFinderV2::AliCaloRawAnalyzerPeakFinderV2() :AliCaloRawAnalyzer("Peak-FinderV2", "PF")
//    fTof(0), 
//							      fAmp(0)
{
  //comment

  fNsampleCut = 5;

  for(int i=0; i < MAXSTART; i++)
    {
      for(int j=0; j < SAMPLERANGE; j++ )
	{
	  fPFAmpVectors[i][j] = new double[100];
	  fPFTofVectors[i][j] = new double[100];
	  fPFAmpVectorsCoarse[i][j] = new double[100];
	  fPFTofVectorsCoarse[i][j] = new double[100];

	  
	  for(int k=0; k < 100; k++ )
	    {
	      fPFAmpVectors[i][j][k] = 0; 
	      fPFTofVectors[i][j][k] = 0;
	      fPFAmpVectorsCoarse[i][j][k] = 0;
	      fPFTofVectorsCoarse[i][j][k] = 0; 
	    }
	}
    }

  LoadVectors();

}


AliCaloRawAnalyzerPeakFinderV2::~AliCaloRawAnalyzerPeakFinderV2()
{
  //comment
  for(int i=0; i < MAXSTART; i++)
    {
      for(int j=0; j < SAMPLERANGE; j++ )
	{
	  delete[] fPFAmpVectors[i][j];
	  delete[] fPFTofVectors[i][j];
	  delete[] fPFAmpVectorsCoarse[i][j];
	  delete[] fPFTofVectorsCoarse[i][j];
	}
    }
}


Double_t  
AliCaloRawAnalyzerPeakFinderV2::ScanCoarse(const Double_t *const array, const int length ) const
{
  Double_t tmpTof = 0;
  Double_t tmpAmp= 0;

  for(int i=0; i < length; i++)
    {
      tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i]; 
      tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i]; 
    }
  
  tmpTof = tmpTof / tmpAmp ;
  return tmpTof;
}


AliCaloFitResults 
AliCaloRawAnalyzerPeakFinderV2::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
{
  // Extracting the amplitude using the Peak-FinderV2 algorithm
  // The amplitude is a weighted sum of the samples using 
  // optimum weights.

  short maxampindex; //index of maximum amplitude
  short maxamp; //Maximum amplitude
  //  fAmp = 0;
  fAmpA[0] = 0;
  fAmpA[1] = 0;
  fAmpA[2] = 0;

  int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
 
  if( index >= 0)
    {
      Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
      Float_t maxf = TMath::MaxElement(   bunchvector.at(index).GetLength(),  fReversed );
      
      if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )	 // (maxamp - ped) > 900 = Close to saturation (use low gain then)
	{
	  //	  cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
	  return  AliCaloFitResults( maxamp, ped,  -1, maxf,   maxampindex, -1, -1 );
 	}
      
      int first;
      int last;
      
      if ( maxf > fAmpCut )
	{
	  SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxampindex -  bunchvector.at(index).GetStartBin(), &first, &last);
	  int nsamples =  last - first;
	  if( ( nsamples  )  >= fNsampleCut )
	    {
	      int startbin = bunchvector.at(index).GetStartBin();  
	      int n = last - first;  
	      int pfindex = n - fNsampleCut; 
	      pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
	     
	      int dt =  maxampindex - startbin -1; 
	   
	      
	      //    cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
	    
	 
	      //  Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
	      //  cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
	      

	      //	      cout << __FILE__ << __LINE__ << ",  dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;

	      for( int i=0; i < SAMPLERANGE; i++ )
		{
		  for( int j = 0; j < 3; j++ )
		    {
		      //    fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j]; 
		      fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt  +i +j -1 ];
		    }
		}
	      
	      double diff = 9999;
	      int tmpindex = 0;

	      for(int k=0; k < 3; k ++)
		{
		  //		  cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
		  if(  TMath::Abs(fAmpA[k] - ( maxamp - ped) )  < diff)
		    {
		      diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
		      tmpindex = k; 
		    }
		}
	      
	      Float_t tmptof = ScanCoarse( &fReversed[dt] , n );

	      //      Double_t tofnew = PolTof(tmptof) + ( dt + startbin  )*100  ;
	      //	     int dt =  maxampindex - startbin -1;    
	    
	      //	      Double_t tofnew = PolTof(tmptof) + startbin*100 ;

	      //      Double_t tofnew = PolTof(tmptof)  +  maxampindex*100 ;

	      
	      //	      Double_t tofnew = PolTof(tmptof)  + (dt + startbin + tmpindex )*100 ; 
	      Double_t tofnew = PolTof(tmptof)  ; 

	      //	      tmptof= tofnew;

	      //	      Double_t tofnew = PolTof(tmptof) + maxampindex  ;


	      if(tmptof < 0 )
		{
		  cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;  
		}

  
	      if( tmptof < -1 )
		{
		  tmpindex = 0;
		}
	      else
		if( tmptof  > -1 && tmptof < 100 )
		  {
		    tmpindex =1;
		  }
		else
		  {
		    tmpindex = 2;
		  }
	      

	      double tof = 0;
	      
	      for(int k=0; k < SAMPLERANGE; k++   )
		{
		  tof +=  fPFTofVectors[0][pfindex][k]*fReversed[ dt  +k + tmpindex -1 ];   
		}
	      
	      //	      cout << __FILE__ << __LINE__ <<  "tofRaw =   "<< tof /  fAmpA[tmpindex]  << endl;
	      
	      // tof = tof /  fAmpA[tmpindex] +  (dt + startbin)*100;
	      
	      if( TMath::Abs(  (maxf - fAmpA[tmpindex])/maxf )  >   0.1 )
		{
		  fAmpA[tmpindex] = maxf;
		}

	      // tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex]; // ns
	      tof = dt + startbin + tmpindex - 0.01*tof/fAmpA[tmpindex]; // clock ticks
	  
	      
	      return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof  , -2, -3 ); 
	    }
	  else
	    {
	      return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 ); 
	    }
	}
    }
   //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
  return  AliCaloFitResults(-1, -1);
}


Double_t
AliCaloRawAnalyzerPeakFinderV2::PolTof( const double fx1 ) const
{
  
  //Newtons method
  Double_t tolerance = 0.01;
  //  cout << "************************"  << endl;
  Double_t fx2 = PolValue( fx1 );
  Double_t tmpfx1 = fx1; 

  while(  TMath::Abs( fx2 - fx1  ) >  tolerance )
    {
      Double_t der = PolDerivative( tmpfx1 ); 
      //     tmpx = der*( x - val) +x;
      tmpfx1 = ( fx1 - fx2)/der  +tmpfx1;
      
      //    tmpx = der*( val - tmpx ) +tmpx;
      fx2 = PolValue(  tmpfx1 );
      //     cout << __FILE__ << __LINE__ <<  "Der =\t" << der << "  x=\t"<<x<<"\tval="<<val <<  endl;
      //     cout << __FILE__ << __LINE__ <<  "Der =\t" << der << "  tmpx=\t"<< tmpfx1 <<"\tval="<< fx2 <<  endl;
    }
  
  //  cout << __FILE__ << __LINE__ << "CONVERGED !! fx1 = "<< fx1 <<"  tmpfx1 = "<< tmpfx1  <<"  f(tmpfx1) = "<< fx2 << endl;  

  //  cout << "************************"  << endl;
  
  return tmpfx1;

}


Double_t 
AliCaloRawAnalyzerPeakFinderV2::PolValue(const Double_t x) const
{
  static Double_t p0 = -55.69;
  static Double_t p1 = 4.718;
  static Double_t p2 = -0.05587;
  static Double_t p3 = 0.0003185;
  static Double_t p4 = -7.91E-7;
  static Double_t p5 = 7.576E-10;
  
  //  return  p0 + p1*x + p2*TMath::Power(x, 2) + p3*TMath::Power(x, 3) + p4*TMath::Power(x, 4)  + p5*TMath::Power(x, 5);
  
  return  p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x  + p5*x*x*x*x*x;
  
}


Double_t 
AliCaloRawAnalyzerPeakFinderV2::PolDerivative(const Double_t x) const
{
  static Double_t dp0 = 0;
  static Double_t dp1 = 4.718;
  static Double_t dp2 = -0.11174;
  static Double_t dp3 = 0.0009555;
  static Double_t dp4 = -3.164E-6;
  static Double_t dp5 = 3.788E-9;

  //  return  dp0 + dp1 + dp2*x + dp3*TMath::Power(x, 2) + dp4*TMath::Power(x, 3)  + dp5*TMath::Power(x, 4);
  
  
  return  dp0 + dp1 + dp2*x + dp3*x*x + dp4*x*x*x  + dp5*x*x*x*x;
  
}


void 
AliCaloRawAnalyzerPeakFinderV2::LoadVectors()
{
  //Read in the Peak finder vecors from file
  for(int i = 0;  i < MAXSTART ; i++)
    {
      for( int j=0; j < SAMPLERANGE; j++)
	{
	  char filenameCoarse[256];
	  char filename[256];
	  
	  int n = j+fNsampleCut;

	  //	  double start = (double)i+0.5;
	  double start = (double)i+0;

	  sprintf(filename,        "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
	  sprintf(filenameCoarse,  "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);

	  FILE *fp  =  fopen(filename, "r");
	  FILE *fpc =  fopen(filenameCoarse, "r");

	  if( fp == 0 )
	    {
	      AliFatal( Form( "could not open file: %s", filename ) );
	    }
	
	  if(fpc == 0)
	    {
	      AliFatal( Form( "could not open file: %s", filenameCoarse ) );
	    }
	  else
	    {
	      for(int m = 0; m < n ; m++ )
		{
		  cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
 
		  fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
		  //		  fPFAmpVectorsCoarse[i][j][m] = 1;
		  fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
		}
	      
	      fscanf(fp,   "\n" );
	      fscanf(fpc,  "\n" );
	      
	      for(int m = 0; m < n ; m++ )
		{
		  //  fPFTofVectors[i][j][m] = 1;

		  fscanf(fp, "%lf\t",   &fPFTofVectors[i][j][m]  );
		  fscanf(fpc, "%lf\t",  &fPFTofVectorsCoarse[i][j][m]  );  
		  //  fPFTofVectorsCoarse[i][j][m] = 1;  
		}
	      fclose (fp);
	      fclose (fpc);
	    }
	}
    }
}
Commit	Line	Data
48a2e3eb	1	/**************************************************************************
	2	* This file is property of and copyright by *
	3	* the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
	4	* *
	5	* Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
	6	* *
	7	* Contributors are mentioned in the code where appropriate. *
	8	* Please report bugs to p.t.hille@fys.uio.no *
	9	* *
	10	* Permission to use, copy, modify and distribute this software and its *
	11	* documentation strictly for non-commercial purposes is hereby granted *
	12	* without fee, provided that the above copyright notice appears in all *
	13	* copies and that both the copyright notice and this permission notice *
	14	* appear in the supporting documentation. The authors make no claims *
	15	* about the suitability of this software for any purpose. It is *
	16	* provided "as is" without express or implied warranty. *
	17	**************************************************************************/
	18
	19	// The Peak-Finder algorithm
	20	// The amplitude is extracted as a
	21	// weighted sum of the samples using the
	22	// best possible weights.
	23	// The wights is calculated only once and the
	24	// Actual extraction of amplitude and peak position
	25	// Is done with a simple vector multiplication, allowing for
	26	// Extreemely fast computations.
	27
	28	#include "AliCaloRawAnalyzerPeakFinderV2.h"
	29	#include "AliCaloBunchInfo.h"
	30	#include "AliCaloFitResults.h"
	31	#include <iostream>
	32	#include "unistd.h"
	33	#include "TMath.h"
	34	#include "AliLog.h"
	35
	36	using namespace std;
	37
	38	ClassImp( AliCaloRawAnalyzerPeakFinderV2 )
	39
	40	AliCaloRawAnalyzerPeakFinderV2::AliCaloRawAnalyzerPeakFinderV2() :AliCaloRawAnalyzer("Peak-FinderV2", "PF")
	41	// fTof(0),
	42	// fAmp(0)
	43	{
	44	//comment
	45
	46	fNsampleCut = 5;
	47
	48	for(int i=0; i < MAXSTART; i++)
	49	{
	50	for(int j=0; j < SAMPLERANGE; j++ )
	51	{
	52	fPFAmpVectors[i][j] = new double[100];
	53	fPFTofVectors[i][j] = new double[100];
	54	fPFAmpVectorsCoarse[i][j] = new double[100];
	55	fPFTofVectorsCoarse[i][j] = new double[100];
	56
	57
	58	for(int k=0; k < 100; k++ )
	59	{
	60	fPFAmpVectors[i][j][k] = 0;
	61	fPFTofVectors[i][j][k] = 0;
	62	fPFAmpVectorsCoarse[i][j][k] = 0;
	63	fPFTofVectorsCoarse[i][j][k] = 0;
	64	}
65	}
66	}
67
68	LoadVectors();
69
70	}
71
72
73	AliCaloRawAnalyzerPeakFinderV2::~AliCaloRawAnalyzerPeakFinderV2()
74	{
75	//comment
76	for(int i=0; i < MAXSTART; i++)
77	{
78	for(int j=0; j < SAMPLERANGE; j++ )
79	{
80	delete[] fPFAmpVectors[i][j];
81	delete[] fPFTofVectors[i][j];
82	delete[] fPFAmpVectorsCoarse[i][j];
83	delete[] fPFTofVectorsCoarse[i][j];
84	}
85	}
86	}
87
88
89	Double_t
90	AliCaloRawAnalyzerPeakFinderV2::ScanCoarse(const Double_t *const array, const int length ) const
91	{
92	Double_t tmpTof = 0;
93	Double_t tmpAmp= 0;
94
95	for(int i=0; i < length; i++)
96	{
97	tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i];
98	tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i];
99	}
100
101	tmpTof = tmpTof / tmpAmp ;
102	return tmpTof;
103	}
104
105
106	AliCaloFitResults
107	AliCaloRawAnalyzerPeakFinderV2::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
108	{
109	// Extracting the amplitude using the Peak-FinderV2 algorithm
110	// The amplitude is a weighted sum of the samples using
111	// optimum weights.
112
113	short maxampindex; //index of maximum amplitude
114	short maxamp; //Maximum amplitude
115	// fAmp = 0;
116	fAmpA[0] = 0;
117	fAmpA[1] = 0;
118	fAmpA[2] = 0;
119
120	int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
121
122	if( index >= 0)
123	{
124	Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
125	Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
126
127	if( maxf < fAmpCut \|\| ( maxamp - ped) > 900 ) // (maxamp - ped) > 900 = Close to saturation (use low gain then)
128	{
129	// cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex << endl;
130	return AliCaloFitResults( maxamp, ped, -1, maxf, maxampindex, -1, -1 );
131	}
132
133	int first;
134	int last;
135
136	if ( maxf > fAmpCut )
137	{
138	SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex - bunchvector.at(index).GetStartBin(), &first, &last);
139	int nsamples = last - first;
140	if( ( nsamples ) >= fNsampleCut )
141	{
142	int startbin = bunchvector.at(index).GetStartBin();
143	int n = last - first;
144	int pfindex = n - fNsampleCut;
145	pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex;
146
147	int dt = maxampindex - startbin -1;
148
149
150	// cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl;
151
152
153
154	// Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
155	// cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
156
157
158	// cout << __FILE__ << __LINE__ << ", dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl;
159
160	for( int i=0; i < SAMPLERANGE; i++ )
161	{
162	for( int j = 0; j < 3; j++ )
163	{
164	// fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j];
165	fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt +i +j -1 ];
166	}
167	}
168
169	double diff = 9999;
170	int tmpindex = 0;
171
172	for(int k=0; k < 3; k ++)
173	{
174	// cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl;
175	if( TMath::Abs(fAmpA[k] - ( maxamp - ped) ) < diff)
176	{
177	diff = TMath::Abs(fAmpA[k] - ( maxamp - ped));
178	tmpindex = k;
179	}
180	}
181
182	Float_t tmptof = ScanCoarse( &fReversed[dt] , n );
183
184	// Double_t tofnew = PolTof(tmptof) + ( dt + startbin )*100 ;
185	// int dt = maxampindex - startbin -1;
186
187	// Double_t tofnew = PolTof(tmptof) + startbin*100 ;
188
189	// Double_t tofnew = PolTof(tmptof) + maxampindex*100 ;
190
191
192	// Double_t tofnew = PolTof(tmptof) + (dt + startbin + tmpindex )*100 ;
193	Double_t tofnew = PolTof(tmptof) ;
194
195	// tmptof= tofnew;
196
197	// Double_t tofnew = PolTof(tmptof) + maxampindex ;
198
199
200	if(tmptof < 0 )
201	{
202	cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl;
203	}
204
205
206	if( tmptof < -1 )
207	{
208	tmpindex = 0;
209	}
210	else
211	if( tmptof > -1 && tmptof < 100 )
212	{
213	tmpindex =1;
214	}
215	else
216	{
217	tmpindex = 2;
218	}
219
220
221
222	double tof = 0;
223
224	for(int k=0; k < SAMPLERANGE; k++ )
225	{
226	tof += fPFTofVectors[0][pfindex][k]*fReversed[ dt +k + tmpindex -1 ];
227	}
228
229	// cout << __FILE__ << __LINE__ << "tofRaw = "<< tof / fAmpA[tmpindex] << endl;
230
231	// tof = tof / fAmpA[tmpindex] + (dt + startbin)*100;
232
233	if( TMath::Abs( (maxf - fAmpA[tmpindex])/maxf ) > 0.1 )
234	{
235	fAmpA[tmpindex] = maxf;
236	}
237
238	// tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex]; // ns
239	tof = dt + startbin + tmpindex - 0.01*tof/fAmpA[tmpindex]; // clock ticks
240
241
242	return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof , -2, -3 );
243	}
244	else
245	{
246	return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 );
247	}
248	}
249	}
250	// cout << __FILE__ << __LINE__ << "WARNING, returning amp = -1 " << endl;
251	return AliCaloFitResults(-1, -1);
252	}
253
254
255	Double_t
256	AliCaloRawAnalyzerPeakFinderV2::PolTof( const double fx1 ) const
257	{
258
259	//Newtons method
260	Double_t tolerance = 0.01;
261	// cout << "************************" << endl;
262	Double_t fx2 = PolValue( fx1 );
263	Double_t tmpfx1 = fx1;
264
265	while( TMath::Abs( fx2 - fx1 ) > tolerance )
266	{
267	Double_t der = PolDerivative( tmpfx1 );
268	// tmpx = der*( x - val) +x;
269	tmpfx1 = ( fx1 - fx2)/der +tmpfx1;
270
271	// tmpx = der*( val - tmpx ) +tmpx;
272	fx2 = PolValue( tmpfx1 );
273	// cout << __FILE__ << __LINE__ << "Der =\t" << der << " x=\t"<<x<<"\tval="<<val << endl;
274	// cout << __FILE__ << __LINE__ << "Der =\t" << der << " tmpx=\t"<< tmpfx1 <<"\tval="<< fx2 << endl;
275	}
276
277	// cout << __FILE__ << __LINE__ << "CONVERGED !! fx1 = "<< fx1 <<" tmpfx1 = "<< tmpfx1 <<" f(tmpfx1) = "<< fx2 << endl;
278
279	// cout << "************************" << endl;
280
281	return tmpfx1;
282
283	}
284
285
286	Double_t
287	AliCaloRawAnalyzerPeakFinderV2::PolValue(const Double_t x) const
288	{
289	static Double_t p0 = -55.69;
290	static Double_t p1 = 4.718;
291	static Double_t p2 = -0.05587;
292	static Double_t p3 = 0.0003185;
293	static Double_t p4 = -7.91E-7;
294	static Double_t p5 = 7.576E-10;
295
296	// return p0 + p1x + p2TMath::Power(x, 2) + p3TMath::Power(x, 3) + p4TMath::Power(x, 4) + p5*TMath::Power(x, 5);
297
298	return p0 + p1x + p2xx + p3xxx + p4xxxx + p5xxxx*x;
299
300	}
301
302
303	Double_t
304	AliCaloRawAnalyzerPeakFinderV2::PolDerivative(const Double_t x) const
305	{
306	static Double_t dp0 = 0;
307	static Double_t dp1 = 4.718;
308	static Double_t dp2 = -0.11174;
309	static Double_t dp3 = 0.0009555;
310	static Double_t dp4 = -3.164E-6;
311	static Double_t dp5 = 3.788E-9;
312
313	// return dp0 + dp1 + dp2x + dp3TMath::Power(x, 2) + dp4TMath::Power(x, 3) + dp5TMath::Power(x, 4);
314
315
316
317	return dp0 + dp1 + dp2x + dp3xx + dp4xxx + dp5xxxx;
318
319	}
320
321
322	void
323	AliCaloRawAnalyzerPeakFinderV2::LoadVectors()
324	{
325	//Read in the Peak finder vecors from file
326	for(int i = 0; i < MAXSTART ; i++)
327	{
328	for( int j=0; j < SAMPLERANGE; j++)
329	{
330	char filenameCoarse[256];
331	char filename[256];
332
333	int n = j+fNsampleCut;
334
335	// double start = (double)i+0.5;
336	double start = (double)i+0;
337
338	sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n);
339	sprintf(filenameCoarse, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n);
340
341	FILE *fp = fopen(filename, "r");
342	FILE *fpc = fopen(filenameCoarse, "r");
343
344	if( fp == 0 )
345	{
346	AliFatal( Form( "could not open file: %s", filename ) );
347	}
348
349	if(fpc == 0)
350	{
351	AliFatal( Form( "could not open file: %s", filenameCoarse ) );
352	}
353	else
354	{
355	for(int m = 0; m < n ; m++ )
356	{
357	cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl;
358
359	fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] );
360	// fPFAmpVectorsCoarse[i][j][m] = 1;
361	fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] );
362	}
363
364	fscanf(fp, "\n" );
365	fscanf(fpc, "\n" );
366
367	for(int m = 0; m < n ; m++ )
368	{
369	// fPFTofVectors[i][j][m] = 1;
370
371	fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] );
372	fscanf(fpc, "%lf\t", &fPFTofVectorsCoarse[i][j][m] );
373	// fPFTofVectorsCoarse[i][j][m] = 1;
374	}
375	fclose (fp);
376	fclose (fpc);
377	}
378	}
379	}
380	}
381
382