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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * 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.                  *
 **************************************************************************/

/*
 
 
Author: R. GUERNANE LPSC Grenoble CNRS/IN2P3
*/

#include "AliCaloRawAnalyzerFakeALTRO.h"
#include "AliCaloBunchInfo.h"
#include "AliCaloFitResults.h"
#include "AliLog.h"
#include "TMath.h"
#include <stdexcept>
#include <iostream>
#include "TF1.h"
#include "TGraph.h"
#include "AliCaloConstants.h"

using namespace std;


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


AliCaloRawAnalyzerFakeALTRO::AliCaloRawAnalyzerFakeALTRO() : AliCaloRawAnalyzerFitter("Chi Square Fit", "LMS")
// fkEulerSquared(7.389056098930650227),
//						 fTf1(0),
//						 fTau(2.35),
//						 fFixTau(kTRUE)
{
  //comment

  /*
  for(int i=0; i < ALTROMAXSAMPLES; 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);
  }
  */
  
}


AliCaloRawAnalyzerFakeALTRO::~AliCaloRawAnalyzerFakeALTRO()
{
  //delete fTf1;
}


/*
void 
AliCaloRawAnalyzerFakeALTRO::InitFormula( TF1* f)
{
  f = new TF1( "myformula", "[0]*((x - [1])/[2])^2*exp(-2*(x -[1])/[2])",  0, 30 ); 
}
*/

AliCaloFitResults
AliCaloRawAnalyzerFakeALTRO::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, Ret::kCrude, maxf, timebinOffset);
 	}            
      else if ( maxf >= fAmpCut )
	{
	  int first = 0;
	  int last = 0;
	  SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, fFitArrayCut );
	  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, Ret::kNoFit, maxf, timebinOffset,
					  timebinOffset, Ret::kDummy, Ret::kDummy,  Ret::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 , Ret::kFitPar,
					  fTf1->GetParameter(0)/fkEulerSquared, 
					  tmax,
					  timebinOffset,  
					  fTf1->GetChisquare(), 
					  fTf1->GetNDF(),
					  Ret::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, Ret::kCrude, maxf, timebinOffset,
					timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
	    }
        } // ampcut
    }
  return AliCaloFitResults(  Ret::kInvalid,  Ret::kInvalid );
  
}


/*
void 
AliCaloRawAnalyzerFakeALTRO::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
de39a0ff	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17
	18
	19	Author: R. GUERNANE LPSC Grenoble CNRS/IN2P3
	20	*/
	21
	22	#include "AliCaloRawAnalyzerFakeALTRO.h"
	23	#include "AliCaloBunchInfo.h"
	24	#include "AliCaloFitResults.h"
	25	#include "AliLog.h"
	26	#include "TMath.h"
	27	#include <stdexcept>
	28	#include <iostream>
	29	#include "TF1.h"
	30	#include "TGraph.h"
168c7b3c	31	#include "AliCaloConstants.h"
de39a0ff	32
	33	using namespace std;
	34
	35
396baaf6	36	//#define BAD 4 //CRAP PTH
de39a0ff	37
	38	ClassImp( AliCaloRawAnalyzerFakeALTRO )
	39
	40
396baaf6	41	AliCaloRawAnalyzerFakeALTRO::AliCaloRawAnalyzerFakeALTRO() : AliCaloRawAnalyzerFitter("Chi Square Fit", "LMS")
	42	// fkEulerSquared(7.389056098930650227),
	43	// fTf1(0),
	44	// fTau(2.35),
	45	// fFixTau(kTRUE)
de39a0ff	46	{
de39a0ff	47	//comment
396baaf6	48
396baaf6	49	/*
92d9f317	50	for(int i=0; i < ALTROMAXSAMPLES; i++)
de39a0ff	51	{
	52	fXaxis[i] = i;
	53	}
	54
	55	fTf1 = new TF1( "myformula", "[0]((x - [1])/[2])^2exp(-2*(x -[1])/[2])", 0, 30 );
	56	if (fFixTau) {
	57	fTf1->FixParameter(2, fTau);
	58	}
	59	else {
	60	fTf1->ReleaseParameter(2); // allow par. to vary
	61	fTf1->SetParameter(2, fTau);
	62	}
396baaf6	63	*/
396baaf6	64
de39a0ff	65	}
	66
	67
	68	AliCaloRawAnalyzerFakeALTRO::~AliCaloRawAnalyzerFakeALTRO()
	69	{
396baaf6	70	//delete fTf1;
de39a0ff	71	}
	72
	73
396baaf6	74	/*
	75	void
	76	AliCaloRawAnalyzerFakeALTRO::InitFormula( TF1* f)
	77	{
	78	f = new TF1( "myformula", "[0]((x - [1])/[2])^2exp(-2*(x -[1])/[2])", 0, 30 );
	79	}
	80	*/
de39a0ff	81
	82	AliCaloFitResults
	83	AliCaloRawAnalyzerFakeALTRO::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 )
	84	{
	85	// Extracting signal parameters using fitting
	86	short maxampindex; //index of maximum amplitude
	87	short maxamp; //Maximum amplitude
	88	int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
	89
	90	if( index >= 0)
	91	{
	92	Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
	93	Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
	94	short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
	95	// timebinOffset is timebin value at maximum (maxrev)
	96	short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
	97	if( maxf < fAmpCut \|\| ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
	98	{
168c7b3c	99	return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
de39a0ff	100	}
	101	else if ( maxf >= fAmpCut )
	102	{
	103	int first = 0;
	104	int last = 0;
92d9f317	105	SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut );
de39a0ff	106	int nsamples = last - first + 1;
	107
	108	if( ( nsamples ) >= fNsampleCut )
	109	{
	110	Float_t tmax = (maxrev - first); // local tmax estimate
	111	TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
	112	fTf1->SetParameter(0, maxf*fkEulerSquared );
	113	fTf1->SetParameter(1, tmax - fTau);
	114	// set rather loose parameter limits
	115	fTf1->SetParLimits(0, 0.5maxffkEulerSquared, 2maxffkEulerSquared );
	116	fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
	117
	118	if (fFixTau) {
	119	fTf1->FixParameter(2, fTau);
	120	}
	121	else {
	122	fTf1->ReleaseParameter(2); // allow par. to vary
	123	fTf1->SetParameter(2, fTau);
	124	}
	125
	126	Short_t tmpStatus = 0;
	127	try {
	128	tmpStatus = graph->Fit(fTf1, "Q0RW");
	129	}
	130	catch (const std::exception & e) {
	131	AliError( Form("TGraph Fit exception %s", e.what()) );
168c7b3c	132	return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
168c7b3c	133	timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
de39a0ff	134	}
	135
	136	if( fVerbose == true )
	137	{
	138	AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
	139	PrintFitResult( fTf1 ) ;
	140	}
	141	// global tmax
	142	tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
	143	+ fTf1->GetParameter(2); // +tau, makes sum tmax
	144
	145	delete graph;
168c7b3c	146	return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
de39a0ff	147	fTf1->GetParameter(0)/fkEulerSquared,
	148	tmax,
	149	timebinOffset,
	150	fTf1->GetChisquare(),
	151	fTf1->GetNDF(),
168c7b3c	152	Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
de39a0ff	153
	154	// delete graph;
	155
	156	}
	157	else
	158	{
	159	Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
	160	Int_t ndf = last - first - 1; // nsamples - 2
168c7b3c	161	return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
168c7b3c	162	timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
de39a0ff	163	}
	164	} // ampcut
	165	}
168c7b3c	166	return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
de39a0ff	167
	168	}
	169
	170
396baaf6	171
396baaf6	172	/*
de39a0ff	173	void
	174	AliCaloRawAnalyzerFakeALTRO::PrintFitResult(const TF1 *f) const
	175	{
	176	//comment
	177	cout << endl;
	178	cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
	179	cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
	180	cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
	181	cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
	182	// cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
	183	cout << endl << endl;
	184	}
396baaf6	185	*/
de39a0ff	186