[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 <stdexcept>
#include <iostream>
#include "TF1.h"
#include "TGraph.h"

using namespace std;


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


AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit", "LMS"),
						 fkEulerSquared(7.389056098930650227),
						 fTf1(0),
						 fTau(2.35),
						 fFixTau(kTRUE)
{
  //comment
  for(int i=0; i < MAXSAMPLES; i++)
    {
      fXaxis[i] = i;
    }
  
  fTf1 = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])",  0, 30 ); 
  if (fFixTau) {
    fTf1->FixParameter(2, fTau);
  }
  else {
    fTf1->ReleaseParameter(2); // allow par. to vary
    fTf1->SetParameter(2, fTau);
  }
 
}


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  );
      Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
      short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
      // timebinOffset is timebin value at maximum (maxrev)
      short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
      if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
	{
	  return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
 	}            
      else if ( maxf >= fAmpCut )
	{
	  int first = 0;
	  int last = 0;
	  SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);
	  int nsamples =  last - first + 1;
	  
	  if( ( nsamples  )  >= fNsampleCut )
	    {
	      Float_t tmax = (maxrev - first); // local tmax estimate
	      TGraph *graph =  new TGraph(  nsamples, fXaxis,  &fReversed[first] );
	      fTf1->SetParameter(0, maxf*fkEulerSquared );
	      fTf1->SetParameter(1, tmax - fTau); 
	      // set rather loose parameter limits
	      fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared );
	      fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4); 

	      if (fFixTau) {
		fTf1->FixParameter(2, fTau);
	      }
	      else {
		fTf1->ReleaseParameter(2); // allow par. to vary
		fTf1->SetParameter(2, fTau);
	      }

	      Short_t tmpStatus = 0;
	      try {
		tmpStatus =  graph->Fit(fTf1, "Q0RW");
	      }
	      catch (const std::exception & e) {
		AliError( Form("TGraph Fit exception %s", e.what()) ); 
		return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, timebinOffset,
					  timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
	      }

	      if( fVerbose == true )
		{
		  AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) ); 
		  PrintFitResult( fTf1 ) ;
		}  
	      // global tmax
	      tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
		+ fTf1->GetParameter(2); // +tau, makes sum tmax
	      
	        delete graph;
		return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar,
					  fTf1->GetParameter(0)/fkEulerSquared, 
					  tmax,
					  timebinOffset,  
					  fTf1->GetChisquare(), 
					  fTf1->GetNDF(),
					  AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
				
		//     delete graph;
	
	    }
	  else
	    {
	      Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
	      Int_t ndf = last - first - 1; // nsamples - 2
	      return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset,
					timebinOffset, chi2, ndf, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
	    }
        } // ampcut
    }
  return AliCaloFitResults( AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
  
}


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"
d655d7dd	31	#include "TMath.h"
507751ce	32	#include <stdexcept>
d655d7dd	33	#include <iostream>
	34	#include "TF1.h"
	35	#include "TGraph.h"
	36
	37	using namespace std;
	38
	39
	40	#define BAD 4 //CRAP PTH
e37e3c84	41
e37e3c84	42	ClassImp( AliCaloRawAnalyzerLMS )
d655d7dd	43
d655d7dd	44
48a2e3eb	45	AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit", "LMS"),
57839add	46	fkEulerSquared(7.389056098930650227),
f57baa2d	47	fTf1(0),
	48	fTau(2.35),
	49	fFixTau(kTRUE)
d655d7dd	50	{
	51	//comment
	52	for(int i=0; i < MAXSAMPLES; i++)
	53	{
	54	fXaxis[i] = i;
	55	}
	56
f57baa2d	57	fTf1 = new TF1( "myformula", "[0]((x - [1])/[2])^2exp(-2*(x -[1])/[2])", 0, 30 );
	58	if (fFixTau) {
	59	fTf1->FixParameter(2, fTau);
	60	}
	61	else {
	62	fTf1->ReleaseParameter(2); // allow par. to vary
	63	fTf1->SetParameter(2, fTau);
	64	}
	65
d655d7dd	66	}
	67
	68
57839add	69	AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
d655d7dd	70	{
	71	delete fTf1;
	72	}
	73
	74
	75
57839add	76	AliCaloFitResults
57839add	77	AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
d655d7dd	78	{
	79	// Extracting signal parameters using fitting
	80	short maxampindex; //index of maximum amplitude
	81	short maxamp; //Maximum amplitude
	82	int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
	83
	84	if( index >= 0)
	85	{
	86	Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
d655d7dd	87	Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
f57baa2d	88	short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
507751ce	89	// timebinOffset is timebin value at maximum (maxrev)
f57baa2d	90	short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
2cd0ffda	91	if( maxf < fAmpCut \|\| ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
d655d7dd	92	{
2cd0ffda	93	return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
	94	}
	95	else if ( maxf >= fAmpCut )
	96	{
	97	int first = 0;
	98	int last = 0;
3b8fd9fe	99	SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
3b8fd9fe	100	int nsamples = last - first + 1;
d655d7dd	101
f57baa2d	102	if( ( nsamples ) >= fNsampleCut )
d655d7dd	103	{
507751ce	104	Float_t tmax = (maxrev - first); // local tmax estimate
3b8fd9fe	105	TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
d655d7dd	106	fTf1->SetParameter(0, maxf*fkEulerSquared );
507751ce	107	fTf1->SetParameter(1, tmax - fTau);
	108	// set rather loose parameter limits
	109	fTf1->SetParLimits(0, 0.5maxffkEulerSquared, 2maxffkEulerSquared );
	110	fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
f57baa2d	111
	112	if (fFixTau) {
	113	fTf1->FixParameter(2, fTau);
	114	}
	115	else {
	116	fTf1->ReleaseParameter(2); // allow par. to vary
	117	fTf1->SetParameter(2, fTau);
	118	}
d655d7dd	119
507751ce	120	Short_t tmpStatus = 0;
	121	try {
	122	tmpStatus = graph->Fit(fTf1, "Q0RW");
	123	}
	124	catch (const std::exception & e) {
	125	AliError( Form("TGraph Fit exception %s", e.what()) );
2cd0ffda	126	return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, timebinOffset,
2cd0ffda	127	timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
507751ce	128	}
507751ce	129
d655d7dd	130	if( fVerbose == true )
d655d7dd	131	{
57839add	132	AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
d655d7dd	133	PrintFitResult( fTf1 ) ;
d655d7dd	134	}
507751ce	135	// global tmax
	136	tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
	137	+ fTf1->GetParameter(2); // +tau, makes sum tmax
d655d7dd	138
d655d7dd	139	delete graph;
507751ce	140	return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar,
f57baa2d	141	fTf1->GetParameter(0)/fkEulerSquared,
507751ce	142	tmax,
507751ce	143	timebinOffset,
f57baa2d	144	fTf1->GetChisquare(),
	145	fTf1->GetNDF(),
	146	AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
	147
d655d7dd	148	// delete graph;
	149
	150	}
	151	else
	152	{
2cd0ffda	153	Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
	154	Int_t ndf = last - first - 1; // nsamples - 2
	155	return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset,
	156	timebinOffset, chi2, ndf, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
d655d7dd	157	}
2cd0ffda	158	} // ampcut
d655d7dd	159	}
f57baa2d	160	return AliCaloFitResults( AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
d655d7dd	161
	162	}
	163
	164
	165	void
57839add	166	AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
d655d7dd	167	{
	168	//comment
	169	cout << endl;
	170	cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
	171	cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
	172	cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
	173	cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
	174	// cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
	175	cout << endl << endl;
	176	}
f57baa2d	177