X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliCaloRawAnalyzerLMS.cxx;h=e22017085dd5f66e8408bfdc2628ba5c0b558e8c;hb=03fdac6e164e394d567245fb75d1afc8489795af;hp=f6fa73b7a088d92d5ea8d48de000aa62afb091b7;hpb=3b8fd9fe1c7b2d56d712b2f1086ddea59a4c4d27;p=u%2Fmrichter%2FAliRoot.git

diff --git a/EMCAL/AliCaloRawAnalyzerLMS.cxx b/EMCAL/AliCaloRawAnalyzerLMS.cxx
index f6fa73b7a08..e22017085dd 100644
--- a/EMCAL/AliCaloRawAnalyzerLMS.cxx
+++ b/EMCAL/AliCaloRawAnalyzerLMS.cxx
@@ -29,6 +29,7 @@
 #include "AliCaloFitResults.h"
 #include "AliLog.h"
 #include "TMath.h"
+#include <stdexcept>
 #include <iostream>
 #include "TF1.h"
 #include "TGraph.h"
@@ -41,10 +42,11 @@ using namespace std;
 ClassImp( AliCaloRawAnalyzerLMS )
 
 
-AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit"),
+AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square Fit", "LMS"),
 						 fkEulerSquared(7.389056098930650227),
-						 fSig(0),
-						 fTf1(0)
+						 fTf1(0),
+						 fTau(2.35),
+						 fFixTau(kTRUE)
 {
   //comment
   for(int i=0; i < MAXSAMPLES; i++)
@@ -52,9 +54,15 @@ AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzer("Chi Square
       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 ); 
-
+  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);
+  }
+ 
 }
 
 
@@ -76,65 +84,80 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
   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();
+      // 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)
 	{
-	  short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
-	  short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
+	  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 )
+	  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, 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); 
- 
+	      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 ,    tmpStatus,  
-					 fTf1->GetParameter(0)/fkEulerSquared, 
-					 fTf1->GetParameter(1) + timebinOffset,  
-					 fTf1->GetChisquare(), 
-					 fTf1->GetNDF());
-		
-		
+		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
 	    {
-	      return AliCaloFitResults( maxamp, ped, -1, maxf, -1, -1, -1 );
+	      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) ); 
 	    }
-	}
-      else
-	{
-	  return AliCaloFitResults( maxamp , ped, -1, maxf, -1, -1, -1 );
-	}       
-    }
-  else
-    {
-      return AliCaloFitResults( -99, -99 );
+        } // ampcut
     }
-
-  return AliCaloFitResults( -99, -99 );
+  return AliCaloFitResults( AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
   
 }
 
@@ -151,3 +174,4 @@ AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
   //  cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus()  << ",.. !!!!" << endl << endl;
   cout << endl << endl;
 }
+