[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzerKStandard.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 "AliCaloRawAnalyzerKStandard.h"
#include "AliCaloBunchInfo.h"
#include "AliCaloFitResults.h"
#include "AliLog.h"
#include "TMath.h"
#include <stdexcept>
#include <iostream>
#include "TF1.h"
#include "TGraph.h"
#include "TRandom.h"


using namespace std;


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


AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzer("Chi Square ( kStandard )", "KStandard"),
						 fkEulerSquared(7.389056098930650227),
						 fTf1(0),
						 fTau(2.35),
						 fFixTau(kTRUE)
{
  
  fAlgo = Algo::kStandard;
  //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);
    }
}


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


AliCaloFitResults
AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo>  &bunchlist, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
{
  
  Float_t pedEstimate  = 0;
  short maxADC = 0;
  Int_t first = 0;
  Int_t last = 0;
  Int_t bunchIndex = 0;
  Float_t ampEstimate = 0;
  short timeEstimate  = 0;
  Float_t time = 0;
  Float_t amp=0;
  Float_t chi2 = 0;
  Int_t ndf = 0;
  Bool_t fitDone = kFALSE;

  
  int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate, 
					maxADC, timeEstimate, pedEstimate, first, last,   fAmpCut ); 
  
  
  if (ampEstimate >= fAmpCut  ) 
    { 
      time = timeEstimate; 
      Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1); 
      amp = ampEstimate; 
      
      if ( nsamples > 1 && maxADC< OVERFLOWCUT ) 
	{ 
	  FitRaw(first, last, amp, time, chi2, fitDone);
	  time += timebinOffset;
	  timeEstimate += timebinOffset;
	  ndf = nsamples - 2;
	}
    } 
  if ( fitDone ) 
    { 
      Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
      Float_t timeDiff = time - timeEstimate;
      
      if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) 
	{
	  amp     = ampEstimate;
	  time    = timeEstimate; 
	  fitDone = kFALSE;
	} 
    }  
  if (amp >= fAmpCut ) 
    { 
      if ( ! fitDone) 
	{ 
	  amp += (0.5 - gRandom->Rndm()); 
	}
      //Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
      // lowGain  = in.IsLowGain();
     
      time = time * TIMEBINWITH; 
      
      /////////////!!!!!!!!!time -= in.GetL1Phase();

      time -= fL1Phase;

      //    AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
      //  AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf); 
      

      return AliCaloFitResults( -99, -99, fAlgo , amp, time,
				time, chi2, ndf, Ret::kDummy );
      
      
      //	AliCaloFitSubarray(index, maxrev, first, last));
    
    }
  
  
  return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
}


/*  
   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 );
*/


/*
  // 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 
AliCaloRawAnalyzerKStandard::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;
}


//____________________________________________________________________________ 
void
 AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const 
{ // Fits the raw signal time distribution
  
  //--------------------------------------------------
  //Do the fit, different fitting algorithms available
  //--------------------------------------------------
  
  // fprintf(fp, "%s:%d:%s\n", __FILE__, __LINE__, __FUNCTION__ );

  int nsamples = lastTimeBin - firstTimeBin + 1;
  fitDone = kFALSE;
  
  // switch(fFittingAlgorithm) 
  //   {
  //  case Algo::kStandard:
  //    {
  if (nsamples < 3) { return; } // nothing much to fit
  //printf("Standard fitter \n");

  // Create Graph to hold data we will fit 
  
  TGraph *gSig =  new TGraph( nsamples); 
 
  for (int i=0; i<nsamples; i++) 
    {
      Int_t timebin = firstTimeBin + i;    
      gSig->SetPoint(i, timebin, GetReversed(timebin)); 
    }
      
  TF1 * signalF = new TF1("signal", RawResponseFunction, 0, TIMEBINS , 5);
  signalF->SetParameters(10.,5., TAU  ,ORDER,0.); //set all defaults once, just to be safe
  signalF->SetParNames("amp","t0","tau","N","ped");
  signalF->FixParameter(2,TAU); // tau in units of time bin
  signalF->FixParameter(3,ORDER); // order
  signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here 
  signalF->SetParameter(1, time);
  signalF->SetParameter(0, amp);
  // set rather loose parameter limits
  signalF->SetParLimits(0, 0.5*amp, 2*amp );
  signalF->SetParLimits(1, time - 4, time + 4); 
      
  try {			
    gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
    // assign fit results
    amp  = signalF->GetParameter(0); 
    time = signalF->GetParameter(1);
    chi2 = signalF->GetChisquare();
    fitDone = kTRUE;
  }
  catch (const std::exception & e) {
    AliError( Form("TGraph Fit exception %s", e.what()) ); 
    // stay with default amp and time in case of exception, i.e. no special action required
    fitDone = kFALSE;
  }
  delete signalF;
      
  //printf("Std   : Amp %f, time %g\n",amp, time);
  delete gSig; // delete TGraph
      
  //	break;
  //    }//kStandard Fitter
  //----------------------------
 
  /*
    case Algo::kLogFit:
    {
    if (nsamples < 3) { return; } // nothing much to fit
    //printf("LogFit \n");
      
    // Create Graph to hold data we will fit 
    TGraph *gSigLog =  new TGraph( nsamples); 
    for (int i=0; i<nsamples; i++) {
    Int_t timebin = firstTimeBin + i;    
    gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) ); 
    }
      
    TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0,  TIMEBINS , 5);
    signalFLog->SetParameters(2.3, 5.,TAU,ORDER,0.); //set all defaults once, just to be safe
    signalFLog->SetParNames("amplog","t0","tau","N","ped");
    signalFLog->FixParameter(2,TAU); // tau in units of time bin
    signalFLog->FixParameter(3, ORDER); // order
    signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here 
    signalFLog->SetParameter(1, time);
    if (amp>=1) {
    signalFLog->SetParameter(0, TMath::Log(amp));
    }
      
    gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points
      
    // assign fit results
    Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp
    amp = TMath::Exp(amplog);
    time = signalFLog->GetParameter(1);
    fitDone = kTRUE;
      
    delete signalFLog;
    //printf("LogFit: Amp %f, time %g\n",amp, time);
    delete gSigLog; 
    break;
    } //kLogFit 
      //----------------------------	
      //----------------------------
      }//switch fitting algorithms
  */
  return;
}


//__________________________________________________________________
void 
AliCaloRawAnalyzerKStandard::FitParabola(const TGraph *gSig, Float_t & amp) const 
{
  //BEG YS alternative methods to calculate the amplitude
  Double_t * ymx = gSig->GetX() ; 
  Double_t * ymy = gSig->GetY() ; 
  const Int_t kN = 3 ; 
  Double_t ymMaxX[kN] = {0., 0., 0.} ; 
  Double_t ymMaxY[kN] = {0., 0., 0.} ; 
  Double_t ymax = 0. ; 
  // find the maximum amplitude
  Int_t ymiMax = 0 ;  
  for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
    if (ymy[ymi] > ymMaxY[0] ) {
      ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
      ymMaxX[0] = ymx[ymi] ;
      ymiMax = ymi ; 
    }
  }
  // find the maximum by fitting a parabola through the max and the two adjacent samples
  if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
    ymMaxY[1] = ymy[ymiMax+1] ;
    ymMaxY[2] = ymy[ymiMax-1] ; 
    ymMaxX[1] = ymx[ymiMax+1] ;
    ymMaxX[2] = ymx[ymiMax-1] ; 
    if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
      //fit a parabola through the 3 points y= a+bx+x*x*x
      Double_t sy = 0 ; 
      Double_t sx = 0 ; 
      Double_t sx2 = 0 ; 
      Double_t sx3 = 0 ; 
      Double_t sx4 = 0 ; 
      Double_t sxy = 0 ; 
      Double_t sx2y = 0 ; 
      for (Int_t i = 0; i < kN ; i++) {
        sy += ymMaxY[i] ; 
        sx += ymMaxX[i] ; 		
        sx2 += ymMaxX[i]*ymMaxX[i] ; 
        sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; 
        sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ; 
        sxy += ymMaxX[i]*ymMaxY[i] ; 
        sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ; 
      }
      Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2); 
      Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
      Double_t c  = cN / cD ; 
      Double_t b  = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
      Double_t a  = (sy-b*sx-c*sx2)/kN  ;
      Double_t xmax = -b/(2*c) ; 
      ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
      amp = ymax;
    }
  }
  
  Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ; 
  if (diff > 0.1) 
    amp = ymMaxY[0] ; 
  //printf("Yves   : Amp %f, time %g\n",amp, time);
  //END YS
  return;
}


//__________________________________________________________________
Double_t 
AliCaloRawAnalyzerKStandard::RawResponseFunction(Double_t *x, Double_t *par)
{
  // Matches version used in 2007 beam test
  //
  // Shape of the electronics raw reponse:
  // It is a semi-gaussian, 2nd order Gamma function of the general form
  //
  // xx = (t - t0 + tau) / tau  [xx is just a convenient help variable]
  // F = A * (xx**N * exp( N * ( 1 - xx) )   for xx >= 0
  // F = 0                                   for xx < 0 
  //
  // parameters:
  // A:   par[0]   // Amplitude = peak value
  // t0:  par[1]
  // tau: par[2]
  // N:   par[3]
  // ped: par[4]
  //
  Double_t signal = 0.;
  Double_t tau    = par[2];
  Double_t n      = par[3];
  Double_t ped    = par[4];
  Double_t xx     = ( x[0] - par[1] + tau ) / tau ;

  if (xx <= 0) 
    signal = ped ;  
  else {  
    signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; 
  }
  return signal ;  
}
Commit	Line	Data
92d9f317	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
	21	// FRom CALO raw data using
	22	// least square fit for the
	23	// Moment assuming identical and
	24	// independent errors (equivalent with chi square)
	25	//
	26
	27	#include "AliCaloRawAnalyzerKStandard.h"
	28	#include "AliCaloBunchInfo.h"
	29	#include "AliCaloFitResults.h"
	30	#include "AliLog.h"
	31	#include "TMath.h"
	32	#include <stdexcept>
	33	#include <iostream>
	34	#include "TF1.h"
	35	#include "TGraph.h"
	36	#include "TRandom.h"
	37
	38
	39	using namespace std;
	40
	41
	42	#define BAD 4 //CRAP PTH
	43
	44	ClassImp( AliCaloRawAnalyzerKStandard )
	45
	46
	47	AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzer("Chi Square ( kStandard )", "KStandard"),
	48	fkEulerSquared(7.389056098930650227),
	49	fTf1(0),
	50	fTau(2.35),
	51	fFixTau(kTRUE)
	52	{
	53
	54	fAlgo = Algo::kStandard;
	55	//comment
	56	for(int i=0; i < ALTROMAXSAMPLES; i++)
	57	{
	58	fXaxis[i] = i;
	59	}
	60
	61	fTf1 = new TF1( "myformula", "[0]((x - [1])/[2])^2exp(-2*(x -[1])/[2])", 0, 30 );
	62	if (fFixTau)
	63	{
	64	fTf1->FixParameter(2, fTau);
65	}
66	else
67	{
68	fTf1->ReleaseParameter(2); // allow par. to vary
69	fTf1->SetParameter(2, fTau);
70	}
71	}
72
73
74	AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
75	{
76	delete fTf1;
77	}
78
79
80
81	AliCaloFitResults
82	AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo> &bunchlist, const UInt_t altrocfg1, const UInt_t altrocfg2 )
83	{
84
85	Float_t pedEstimate = 0;
86	short maxADC = 0;
87	Int_t first = 0;
88	Int_t last = 0;
89	Int_t bunchIndex = 0;
90	Float_t ampEstimate = 0;
91	short timeEstimate = 0;
92	Float_t time = 0;
93	Float_t amp=0;
94	Float_t chi2 = 0;
95	Int_t ndf = 0;
96	Bool_t fitDone = kFALSE;
97
98
99	int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
100	maxADC, timeEstimate, pedEstimate, first, last, fAmpCut );
101
102
103	if (ampEstimate >= fAmpCut )
104	{
105	time = timeEstimate;
106	Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
107	amp = ampEstimate;
108
109	if ( nsamples > 1 && maxADC< OVERFLOWCUT )
110	{
111	FitRaw(first, last, amp, time, chi2, fitDone);
112	time += timebinOffset;
113	timeEstimate += timebinOffset;
114	ndf = nsamples - 2;
115	}
116	}
117	if ( fitDone )
118	{
119	Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
120	Float_t timeDiff = time - timeEstimate;
121
122	if ( (TMath::Abs(ampAsymm) > 0.1) \|\| (TMath::Abs(timeDiff) > 2) )
123	{
124	amp = ampEstimate;
125	time = timeEstimate;
126	fitDone = kFALSE;
127	}
128	}
129	if (amp >= fAmpCut )
130	{
131	if ( ! fitDone)
132	{
133	amp += (0.5 - gRandom->Rndm());
134	}
135	//Int_t id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
136	// lowGain = in.IsLowGain();
137
138	time = time * TIMEBINWITH;
139
140	/////////////!!!!!!!!!time -= in.GetL1Phase();
141
142	time -= fL1Phase;
143
144	// AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
145	// AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf);
146
147
148	return AliCaloFitResults( -99, -99, fAlgo , amp, time,
149	time, chi2, ndf, Ret::kDummy );
150
151
152	// AliCaloFitSubarray(index, maxrev, first, last));
153
154	}
155
156
157	return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
158	}
159
160
161
162
163	/*
164	return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
165	timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
166	}
167	} // ampcut
168	}
169	return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
170	*/
171
172
173	/*
174	// Extracting signal parameters using fitting
175	short maxampindex; //index of maximum amplitude
176	short maxamp; //Maximum amplitude
177	int index = SelectBunch( bunchvector, &maxampindex, &maxamp );
178
179	if( index >= 0)
180	{
181	Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed );
182	Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
183	short maxrev = maxampindex - bunchvector.at(index).GetStartBin();
184	// timebinOffset is timebin value at maximum (maxrev)
185	short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
186	if( maxf < fAmpCut \|\| ( maxamp - ped) > fOverflowCut ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
187	{
188	return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
189	}
190	else if ( maxf >= fAmpCut )
191	{
192	int first = 0;
193	int last = 0;
194	SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, fFitArrayCut);
195	int nsamples = last - first + 1;
196
197	if( ( nsamples ) >= fNsampleCut )
198	{
199	Float_t tmax = (maxrev - first); // local tmax estimate
200	TGraph *graph = new TGraph( nsamples, fXaxis, &fReversed[first] );
201	fTf1->SetParameter(0, maxf*fkEulerSquared );
202	fTf1->SetParameter(1, tmax - fTau);
203	// set rather loose parameter limits
204	fTf1->SetParLimits(0, 0.5maxffkEulerSquared, 2maxffkEulerSquared );
205	fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
206
207	if (fFixTau) {
208	fTf1->FixParameter(2, fTau);
209	}
210	else {
211	fTf1->ReleaseParameter(2); // allow par. to vary
212	fTf1->SetParameter(2, fTau);
213	}
214
215	Short_t tmpStatus = 0;
216	try {
217	tmpStatus = graph->Fit(fTf1, "Q0RW");
218	}
219	catch (const std::exception & e) {
220	AliError( Form("TGraph Fit exception %s", e.what()) );
221	return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
222	timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
223	}
224
225	if( fVerbose == true )
226	{
227	AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
228	PrintFitResult( fTf1 ) ;
229	}
230	// global tmax
231	tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
232	+ fTf1->GetParameter(2); // +tau, makes sum tmax
233
234	delete graph;
235	return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
236	fTf1->GetParameter(0)/fkEulerSquared,
237	tmax,
238	timebinOffset,
239	fTf1->GetChisquare(),
240	fTf1->GetNDF(),
241	Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
242
243	// delete graph;
244
245	}
246	else
247	{
248
249	Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
250	Int_t ndf = last - first - 1; // nsamples - 2
251	return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
252	timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
253	}
254	} // ampcut
255	}
256	return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
257
258	}
259	*/
260
261
262
263	void
264	AliCaloRawAnalyzerKStandard::PrintFitResult(const TF1 *f) const
265	{
266	//comment
267	cout << endl;
268	cout << __FILE__ << __LINE__ << "Using this samplerange we get" << endl;
269	cout << __FILE__ << __LINE__ << "AMPLITUDE = " << f->GetParameter(0)/fkEulerSquared << ",.. !!!!" << endl;
270	cout << __FILE__ << __LINE__ << "TOF = " << f->GetParameter(1) << ",.. !!!!" << endl;
271	cout << __FILE__ << __LINE__ << "NDF = " << f->GetNDF() << ",.. !!!!" << endl;
272	// cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus() << ",.. !!!!" << endl << endl;
273	cout << endl << endl;
274	}
275
276
277
278
279
280	//____________________________________________________________________________
281	void
282	AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
283	{ // Fits the raw signal time distribution
284
285	//--------------------------------------------------
286	//Do the fit, different fitting algorithms available
287	//--------------------------------------------------
288
289	// fprintf(fp, "%s:%d:%s\n", __FILE__, __LINE__, __FUNCTION__ );
290
291	int nsamples = lastTimeBin - firstTimeBin + 1;
292	fitDone = kFALSE;
293
294	// switch(fFittingAlgorithm)
295	// {
296	// case Algo::kStandard:
297	// {
298	if (nsamples < 3) { return; } // nothing much to fit
299	//printf("Standard fitter \n");
300
301	// Create Graph to hold data we will fit
302
303	TGraph *gSig = new TGraph( nsamples);
304
305	for (int i=0; i<nsamples; i++)
306	{
307	Int_t timebin = firstTimeBin + i;
308	gSig->SetPoint(i, timebin, GetReversed(timebin));
309	}
310
311	TF1 * signalF = new TF1("signal", RawResponseFunction, 0, TIMEBINS , 5);
312	signalF->SetParameters(10.,5., TAU ,ORDER,0.); //set all defaults once, just to be safe
313	signalF->SetParNames("amp","t0","tau","N","ped");
314	signalF->FixParameter(2,TAU); // tau in units of time bin
315	signalF->FixParameter(3,ORDER); // order
316	signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here
317	signalF->SetParameter(1, time);
318	signalF->SetParameter(0, amp);
319	// set rather loose parameter limits
320	signalF->SetParLimits(0, 0.5amp, 2amp );
321	signalF->SetParLimits(1, time - 4, time + 4);
322
323	try {
324	gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
325	// assign fit results
326	amp = signalF->GetParameter(0);
327	time = signalF->GetParameter(1);
328	chi2 = signalF->GetChisquare();
329	fitDone = kTRUE;
330	}
331	catch (const std::exception & e) {
332	AliError( Form("TGraph Fit exception %s", e.what()) );
333	// stay with default amp and time in case of exception, i.e. no special action required
334	fitDone = kFALSE;
335	}
336	delete signalF;
337
338	//printf("Std : Amp %f, time %g\n",amp, time);
339	delete gSig; // delete TGraph
340
341	// break;
342	// }//kStandard Fitter
343	//----------------------------
344
345	/*
346	case Algo::kLogFit:
347	{
348	if (nsamples < 3) { return; } // nothing much to fit
349	//printf("LogFit \n");
350
351	// Create Graph to hold data we will fit
352	TGraph *gSigLog = new TGraph( nsamples);
353	for (int i=0; i<nsamples; i++) {
354	Int_t timebin = firstTimeBin + i;
355	gSigLog->SetPoint(timebin, timebin, TMath::Log(fRawAnalyzer->GetReversed(timebin) ) );
356	}
357
358	TF1 * signalFLog = new TF1("signalLog", RawResponseFunctionLog, 0, TIMEBINS , 5);
359	signalFLog->SetParameters(2.3, 5.,TAU,ORDER,0.); //set all defaults once, just to be safe
360	signalFLog->SetParNames("amplog","t0","tau","N","ped");
361	signalFLog->FixParameter(2,TAU); // tau in units of time bin
362	signalFLog->FixParameter(3, ORDER); // order
363	signalFLog->FixParameter(4, 0); // pedestal should be subtracted when we get here
364	signalFLog->SetParameter(1, time);
365	if (amp>=1) {
366	signalFLog->SetParameter(0, TMath::Log(amp));
367	}
368
369	gSigLog->Fit(signalFLog, "QROW"); // Note option 'W': equal errors on all points
370
371	// assign fit results
372	Double_t amplog = signalFLog->GetParameter(0); //Not Amp, but Log of Amp
373	amp = TMath::Exp(amplog);
374	time = signalFLog->GetParameter(1);
375	fitDone = kTRUE;
376
377	delete signalFLog;
378	//printf("LogFit: Amp %f, time %g\n",amp, time);
379	delete gSigLog;
380	break;
381	} //kLogFit
382	//----------------------------
383	//----------------------------
384	}//switch fitting algorithms
385	*/
386	return;
387	}
388
389
390	//__________________________________________________________________
391	void
392	AliCaloRawAnalyzerKStandard::FitParabola(const TGraph *gSig, Float_t & amp) const
393	{
394	//BEG YS alternative methods to calculate the amplitude
395	Double_t * ymx = gSig->GetX() ;
396	Double_t * ymy = gSig->GetY() ;
397	const Int_t kN = 3 ;
398	Double_t ymMaxX[kN] = {0., 0., 0.} ;
399	Double_t ymMaxY[kN] = {0., 0., 0.} ;
400	Double_t ymax = 0. ;
401	// find the maximum amplitude
402	Int_t ymiMax = 0 ;
403	for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
404	if (ymy[ymi] > ymMaxY[0] ) {
405	ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
406	ymMaxX[0] = ymx[ymi] ;
407	ymiMax = ymi ;
408	}
409	}
410	// find the maximum by fitting a parabola through the max and the two adjacent samples
411	if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
412	ymMaxY[1] = ymy[ymiMax+1] ;
413	ymMaxY[2] = ymy[ymiMax-1] ;
414	ymMaxX[1] = ymx[ymiMax+1] ;
415	ymMaxX[2] = ymx[ymiMax-1] ;
416	if (ymMaxY[0]ymMaxY[1]ymMaxY[2] > 0) {
417	//fit a parabola through the 3 points y= a+bx+xxx
418	Double_t sy = 0 ;
419	Double_t sx = 0 ;
420	Double_t sx2 = 0 ;
421	Double_t sx3 = 0 ;
422	Double_t sx4 = 0 ;
423	Double_t sxy = 0 ;
424	Double_t sx2y = 0 ;
425	for (Int_t i = 0; i < kN ; i++) {
426	sy += ymMaxY[i] ;
427	sx += ymMaxX[i] ;
428	sx2 += ymMaxX[i]*ymMaxX[i] ;
429	sx3 += ymMaxX[i]ymMaxX[i]ymMaxX[i] ;
430	sx4 += ymMaxX[i]ymMaxX[i]ymMaxX[i]*ymMaxX[i] ;
431	sxy += ymMaxX[i]*ymMaxY[i] ;
432	sx2y += ymMaxX[i]ymMaxX[i]ymMaxY[i] ;
433	}
434	Double_t cN = (sx2ykN-sysx2)(sx3sx-sx2sx2)-(sx2ysx-sxysx2)(sx3kN-sxsx2);
435	Double_t cD = (sx4kN-sx2sx2)(sx3sx-sx2sx2)-(sx4sx-sx3sx2)(sx3kN-sxsx2) ;
436	Double_t c = cN / cD ;
437	Double_t b = ((sx2ykN-sysx2)-c(sx4kN-sx2sx2))/(sx3kN-sx*sx2) ;
438	Double_t a = (sy-bsx-csx2)/kN ;
439	Double_t xmax = -b/(2*c) ;
440	ymax = a + bxmax + cxmax*xmax ;//<========== This is the maximum amplitude
441	amp = ymax;
442	}
443	}
444
445	Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
446	if (diff > 0.1)
447	amp = ymMaxY[0] ;
448	//printf("Yves : Amp %f, time %g\n",amp, time);
449	//END YS
450	return;
451	}
452
453
454
455	//__________________________________________________________________
456	Double_t
457	AliCaloRawAnalyzerKStandard::RawResponseFunction(Double_t x, Double_t par)
458	{
459	// Matches version used in 2007 beam test
460	//
461	// Shape of the electronics raw reponse:
462	// It is a semi-gaussian, 2nd order Gamma function of the general form
463	//
464	// xx = (t - t0 + tau) / tau [xx is just a convenient help variable]
465	// F = A * (xx*N exp( N * ( 1 - xx) ) for xx >= 0
466	// F = 0 for xx < 0
467	//
468	// parameters:
469	// A: par[0] // Amplitude = peak value
470	// t0: par[1]
471	// tau: par[2]
472	// N: par[3]
473	// ped: par[4]
474	//
475	Double_t signal = 0.;
476	Double_t tau = par[2];
477	Double_t n = par[3];
478	Double_t ped = par[4];
479	Double_t xx = ( x[0] - par[1] + tau ) / tau ;
480
481	if (xx <= 0)
482	signal = ped ;
483	else {
484	signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
485	}
486	return signal ;
487	}
488