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

/**************************************************************************
 * This file is property of and copyright by                              *
 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
 *                                                                        *
 * Primary Author: Per Thomas Hille <p.t.hille@fys.uio.no>                *
 *                                                                        *
 * 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.                  *
 **************************************************************************/


// Extraction of amplitude and peak position
// FRom CALO raw data using
// least square fit for the
// Moment assuming identical and 
// independent errors (equivalent with chi square)
// 

#include "AliCaloRawAnalyzerLMS.h"
#include "AliCaloBunchInfo.h"
#include "AliCaloFitResults.h"
#include "AliLog.h"
#include "TMath.h"
#include <iostream>
#include "TF1.h"
#include "TGraph.h"

using namespace std;


#define BAD 4  //CRAP PTH
 
ClassImp( AliCaloRawAnalyzerLMS )


AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit"),
						 fkEulerSquared(7.389056098930650227),
						 fSig(0),
						 fTf1(0)
{
  //comment
  for(int i=0; i < MAXSAMPLES; i++)
    {
      fXaxis[i] = i;
    }
  
  fTf1 = new TF1( "myformula", "[0]*((x - [1])/2.35)^2*exp(-2*(x -[1])/2.35)",  0, 30 ); 
  //  fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^[3]*exp(-[3]*(x -[1])/[2])",  0, 30 ); 

}


AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
{
  delete fTf1;
}


AliCaloFitResults
AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
{
  // Extracting signal parameters using fitting
  short maxampindex; //index of maximum amplitude
  short maxamp; //Maximum amplitude
  int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
  
  if( index >= 0)
    {
      Float_t  ped  = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
      int first;
      int last;
      Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );

      if ( maxf > fAmpCut )
	{
	  short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
	  short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
	  SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);
	  int nsamples =  last - first + 1;
	  
	  if( ( nsamples  )  > fNsampleCut )
	    {
	      
	      TGraph *graph =  new TGraph(  nsamples, fXaxis,  &fReversed[first] );
	      fTf1->SetParameter(0, maxf*fkEulerSquared );
	      fTf1->SetParameter(1, 0.2);
	      //	      fTf1->SetParameter(2,  2);
	      //	      fTf1->SetParameter(2,  3); 
	    
	      //     return   AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1); 

	      Short_t tmpStatus =  graph->Fit(fTf1, "Q0RW");
	     
	      //       return   AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1); 
 
	      if( fVerbose == true )
		{
		  AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) ); 
		  PrintFitResult( fTf1 ) ;
		}  
	      
	        delete graph;
		return AliCaloFitResults( maxamp, ped ,    tmpStatus,  
					 fTf1->GetParameter(0)/fkEulerSquared, 
					 fTf1->GetParameter(1) + timebinOffset,  
					 fTf1->GetChisquare(), 
					 fTf1->GetNDF());
		
		
		//     delete graph;
	
	    }
	  else
	    {
	      return AliCaloFitResults( maxamp, ped, -1, maxf, -1, -1, -1 );
	    }
	}
      else
	{
	  return AliCaloFitResults( maxamp , ped, -1, maxf, -1, -1, -1 );
	}       
    }
  else
    {
      return AliCaloFitResults( -99, -99 );
    }

  return AliCaloFitResults( -99, -99 );
  
}


void 
AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
{
  //comment
  cout << endl;
  cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
  cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
  cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
  cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
  //  cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus()  << ",.. !!!!" << endl << endl;
  cout << endl << endl;
}
Commit	Line	Data
d655d7dd	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 <p.t.hille@fys.uio.no> *
	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
	20	// Extraction of amplitude and peak position
57839add	21	// FRom CALO raw data using
d655d7dd	22	// least square fit for the
	23	// Moment assuming identical and
	24	// independent errors (equivalent with chi square)
	25	//
	26
57839add	27	#include "AliCaloRawAnalyzerLMS.h"
	28	#include "AliCaloBunchInfo.h"
	29	#include "AliCaloFitResults.h"
d655d7dd	30	#include "AliLog.h"
	31	#include "TMath.h"
	32	#include <iostream>
	33	#include "TF1.h"
	34	#include "TGraph.h"
	35
	36	using namespace std;
	37
	38
	39	#define BAD 4 //CRAP PTH
e37e3c84	40
e37e3c84	41	ClassImp( AliCaloRawAnalyzerLMS )
d655d7dd	42
d655d7dd	43
e37e3c84	44	AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit"),
57839add	45	fkEulerSquared(7.389056098930650227),
	46	fSig(0),
	47	fTf1(0)
d655d7dd	48	{
	49	//comment
	50	for(int i=0; i < MAXSAMPLES; i++)
	51	{
	52	fXaxis[i] = i;
	53	}
	54
	55	fTf1 = new TF1( "myformula", "[0]((x - [1])/2.35)^2exp(-2*(x -[1])/2.35)", 0, 30 );
	56	// fTf1 = new TF1( "myformula", "[0]((x - [1])/[2])^[3]exp(-[3]*(x -[1])/[2])", 0, 30 );
	57
	58	}
	59
	60
57839add	61	AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
d655d7dd	62	{
	63	delete fTf1;
	64	}
	65
	66
	67
57839add	68	AliCaloFitResults
57839add	69	AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
d655d7dd	70	{
	71	// Extracting signal parameters using fitting
	72	short maxampindex; //index of maximum amplitude
	73	short maxamp; //Maximum amplitude
	74	int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
	75
	76	if( index >= 0)
	77	{
	78	Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
	79	int first;
	80	int last;
	81	Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
	82
	83	if ( maxf > fAmpCut )
	84	{
3b8fd9fe	85	short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
	86	short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
	87	SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
	88	int nsamples = last - first + 1;
d655d7dd	89
	90	if( ( nsamples ) > fNsampleCut )
	91	{
	92
3b8fd9fe	93	TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
d655d7dd	94	fTf1->SetParameter(0, maxf*fkEulerSquared );
57839add	95	fTf1->SetParameter(1, 0.2);
d655d7dd	96	// fTf1->SetParameter(2, 2);
	97	// fTf1->SetParameter(2, 3);
	98
57839add	99	// return AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1);
d655d7dd	100
	101	Short_t tmpStatus = graph->Fit(fTf1, "Q0RW");
	102
57839add	103	// return AliCaloFitResults( -1 , -1 , -1, -1, -1, -1 , -1);
d655d7dd	104
	105	if( fVerbose == true )
	106	{
57839add	107	AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
d655d7dd	108	PrintFitResult( fTf1 ) ;
	109	}
	110
	111	delete graph;
57839add	112	return AliCaloFitResults( maxamp, ped , tmpStatus,
d655d7dd	113	fTf1->GetParameter(0)/fkEulerSquared,
3b8fd9fe	114	fTf1->GetParameter(1) + timebinOffset,
d655d7dd	115	fTf1->GetChisquare(),
	116	fTf1->GetNDF());
	117
	118
	119	// delete graph;
	120
	121	}
	122	else
	123	{
57839add	124	return AliCaloFitResults( maxamp, ped, -1, maxf, -1, -1, -1 );
d655d7dd	125	}
	126	}
	127	else
	128	{
57839add	129	return AliCaloFitResults( maxamp , ped, -1, maxf, -1, -1, -1 );
d655d7dd	130	}
	131	}
	132	else
	133	{
57839add	134	return AliCaloFitResults( -99, -99 );
d655d7dd	135	}
d655d7dd	136
57839add	137	return AliCaloFitResults( -99, -99 );
d655d7dd	138
	139	}
	140
	141
	142	void
57839add	143	AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
d655d7dd	144	{
	145	//comment
	146	cout << endl;
	147	cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
	148	cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
	149	cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
	150	cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
	151	// cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
	152	cout << endl << endl;
	153	}