remove dependency to aliroot libraries, access of ESDEvent object through abstract...

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

diff --git a/EMCAL/AliCaloRawAnalyzerLMS.cxx b/EMCAL/AliCaloRawAnalyzerLMS.cxx
index f6fa73b7a088d92d5ea8d48de000aa62afb091b7..412630147609817b83a0db5d5e6639dc42eaadcb 100644 (file)
--- 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,20 +42,30 @@ 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)
 {
+  
+  fAlgo = Algo::kLMS;
   //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 ); 
-
+  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 +87,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)
+       {
+         return  AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
+       }            
+      else if ( maxf >= fAmpCut )
        {
-         short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
-         short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
+         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, 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 ,    tmpStatus,  
-                                        fTf1->GetParameter(0)/fkEulerSquared, 
-                                        fTf1->GetParameter(1) + timebinOffset,  
-                                        fTf1->GetChisquare(), 
-                                        fTf1->GetNDF());
-               
-               
+               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
            {
-             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, Ret::kCrude, maxf, timebinOffset,
+                                       timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
            }
-       }
-      else
-       {
-         return AliCaloFitResults( maxamp , ped, -1, maxf, -1, -1, -1 );
-       }       
+        } // ampcut
     }
-  else
-    {
-      return AliCaloFitResults( -99, -99 );
-    }
-
-  return AliCaloFitResults( -99, -99 );
+  return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
   
 }
 
@@ -151,3 +177,4 @@ AliCaloRawAnalyzerLMS::PrintFitResult(const TF1 *f) const
   //  cout << __FILE__ << __LINE__ << "STATUS = " << f->GetStatus()  << ",.. !!!!" << endl << endl;
   cout << endl << endl;
 }
+