more standardizing of result return codes + always return tmax info + loose parameter...
authordsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 19 Mar 2010 00:34:32 +0000 (00:34 +0000)
committerdsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 19 Mar 2010 00:34:32 +0000 (00:34 +0000)
EMCAL/AliCaloFitResults.cxx
EMCAL/AliCaloFitResults.h
EMCAL/AliCaloRawAnalyzer.cxx
EMCAL/AliCaloRawAnalyzer.h
EMCAL/AliCaloRawAnalyzerCrude.cxx
EMCAL/AliCaloRawAnalyzerFastFit.cxx
EMCAL/AliCaloRawAnalyzerLMS.cxx
EMCAL/AliCaloRawAnalyzerNN.cxx
EMCAL/AliCaloRawAnalyzerPeakFinder.cxx
EMCAL/AliEMCALRawUtils.cxx
EMCAL/AliEMCALRawUtils.h

index 2d0348b..37b149e 100644 (file)
 // Applied. fStatus holds information on wether or not 
 // The signal was fitted sucessfully. fStatus might have a different meaning If other 
 // procedures than  A different meaning Fitting is applied 
+
 AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, 
                                     const Short_t fitstatus, const Float_t  amp,  
-                                    const Float_t t0,  const Float_t chi,  
+                                    const Float_t time,  const Int_t maxTimebin, const Float_t chi,  
                                     const Int_t ndf, Int_t minSig,
                                     const AliCaloFitSubarray fitSubarray) : 
   fMaxSig(maxSig),
   fPed(ped), 
   fStatus(fitstatus),
   fAmpSig(amp),
-  fT0(t0),
+  fTime(time),
+  fMaxTimebin(maxTimebin),
   fChi2Sig(chi),
   fNdfSig(ndf),
   fMinSig(minSig),
@@ -47,34 +49,50 @@ AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped,
 
 AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, 
                                     const Short_t fitstatus, const Float_t  amp,  
-                                    const Float_t t0,  const Float_t chi,  
+                                    const Float_t time, const Int_t maxTimebin, const Float_t chi,  
                                     const Int_t ndf, Int_t minSig ) : 
   fMaxSig(maxSig),
   fPed(ped), 
   fStatus(fitstatus),
   fAmpSig(amp),
-  fT0(t0),
+  fTime(time),
+  fMaxTimebin(maxTimebin),
   fChi2Sig(chi),
   fNdfSig(ndf),
   fMinSig(minSig),
   fFitSubarray(kDummy)  
 {
-
 }
 
 
+AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Float_t ped, 
+                                    const Short_t fitstatus, const Float_t  amp,  
+                                    const Int_t maxTimebin) : 
+  fMaxSig(maxSig),
+  fPed(ped), 
+  fStatus(fitstatus),
+  fAmpSig(amp),
+  fTime(maxTimebin),
+  fMaxTimebin(maxTimebin),
+  fChi2Sig( kNoFit ),
+  fNdfSig( kNoFit ),
+  fMinSig( kNoFit ),
+  fFitSubarray( kNoFit ) 
+{
+}
 
 
 AliCaloFitResults::AliCaloFitResults(const Int_t maxSig, const Int_t minSig) : 
   fMaxSig(maxSig),
-  fPed(kNoFit),
-  fStatus( kNoFit ),
-  fAmpSig( kNoFit ), 
-  fT0(kNoFit),  
-  fChi2Sig( kNoFit ),
-  fNdfSig( kNoFit),
+  fPed( kInvalid ),
+  fStatus( kInvalid ),
+  fAmpSig( kInvalid ), 
+  fTime( kInvalid ),
+  fMaxTimebin( kInvalid ),
+  fChi2Sig( kInvalid ),
+  fNdfSig( kInvalid),
   fMinSig (minSig),
-  fFitSubarray(kNoFit)  
+  fFitSubarray(kInvalid)  
 {
 
 }
index bd141bf..3a6cfce 100644 (file)
 class  AliCaloFitResults
 {
  public:
-  enum kReturnCode {kDummy=-1, kCrude=-9, kNoFit=-99, kInvalid=-9999};// possible return values
+  enum kReturnCode {kFitPar=1, kDummy=-1, kCrude=-9, kNoFit=-99, kInvalid=-9999};// possible return values
+  // kFitPar: method fit or parametrization was used
+  // kDummy: just a filler parameter, if e.g. chi2 not available
+  // kCrude: maximum was used
+  // kNoFit: maximum was used, exception handling for fit invoked
+  // kInvalid: could not even look for maximum
 
   explicit AliCaloFitResults( const Int_t maxSig, 
                              const Float_t ped, 
                              const Short_t fitStatus, 
                              const Float_t  amp, 
-                             const Float_t t0,
+                             const Float_t time,
+                             const Int_t maxTimebin,
                              const Float_t chi, 
                              const Int_t ndf, 
                              const Int_t minSig, 
@@ -46,13 +52,21 @@ class  AliCaloFitResults
                              const Float_t ped, 
                              const Short_t fitStatus, 
                              const Float_t  amp, 
-                             const Float_t t0,
+                             const Float_t time,
+                             const Int_t maxTimebin,
                              const Float_t chi, 
                              const Int_t ndf, 
                              const Int_t minSig = kDummy); 
 
+  // shorter interface when no fit is done
+  explicit AliCaloFitResults( const Int_t maxSig, 
+                             const Float_t ped, 
+                             const Short_t fitStatus, 
+                             const Float_t  amp, 
+                             const Int_t maxTimebin); 
+
+  // minimum interface
   explicit AliCaloFitResults( const Int_t maxSig, const Int_t minSig );
-  //AliCaloFitResults( const Int_t maxSig, const Int_t minSig );
 
 
   virtual  ~AliCaloFitResults();
@@ -61,7 +75,9 @@ class  AliCaloFitResults
   Int_t  GetMinSig() const { return fMinSig;};
   Int_t  GetStatus() const  { return fStatus;};
   Float_t   GetAmp() const {  return fAmpSig; };
-  Float_t   GetTof() const {  return fT0; }; 
+  Float_t   GetTof() const {  return fTime; }; 
+  Float_t   GetTime() const {  return fTime; };
+  Int_t   GetMaxTimebin() const {  return fMaxTimebin; };
   Float_t   GetChisSquare() const { return fChi2Sig;};
   Int_t  GetNdf() const { return fNdfSig; };
   AliCaloFitSubarray  GetFitSubarray() const { return fFitSubarray; };
@@ -72,7 +88,8 @@ class  AliCaloFitResults
   Float_t    fPed;      //Pedestal 
   Int_t   fStatus;   //Sucess or failure of fitting pocedure
   Float_t    fAmpSig;   //Amplitude in entities of ADC counts
-  Float_t    fT0;       //Start time of signal in entities of sample intervals 
+  Float_t    fTime;       //peak/max time of signal in entities of sample intervals 
+  Int_t    fMaxTimebin;  // timebin with maximum ADC value
   Float_t    fChi2Sig;  //Chi Square of fit 
   Int_t   fNdfSig;   //Number of degrees of freedom of fit
   Int_t   fMinSig;   //Pedestal 
index 4dbd2c6..38b317b 100644 (file)
@@ -192,6 +192,33 @@ AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxind
 }
 
 
+bool  
+AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
+{
+  // a bunch is considered invalid if the maximum is in the first or last time-bin
+  short tmpmax = -1;
+  int tmpindex = -1;
+  const UShort_t *sig = bunch->GetData();
+
+  for(int i=0; i < bunch->GetLength(); i++ )
+    {
+      if( sig[i] > tmpmax )
+       {
+         tmpmax   =  sig[i];
+         tmpindex =  i; 
+       }
+    }
+  
+  bool bunchOK = true;
+  if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
+    {
+      bunchOK = false;
+    }
+  
+  return  bunchOK;
+}
+
+
 int 
 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) const
 {
@@ -216,10 +243,13 @@ AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,sho
        }
     }
  
-  
-  //  *maxampbin     = indx;
-  //  *maxamplitude  = max;
+  if (bunchindex >= 0) {
+    bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
+    if (! bunchOK) { 
+      bunchindex = -1; 
+    }
+  }
+
   return  bunchindex;
 }
 
@@ -272,7 +302,7 @@ AliCaloFitResults
 AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo>  &/*bunchvector*/, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/)
 { // method to do the selection of what should possibly be fitted
   // not implemented for base class
-  return AliCaloFitResults( 0, 0, 0, 0, 0, 0, 0 );
+  return AliCaloFitResults( 0, 0 );
 }
 
 
index 2575d3b..7fab35e 100644 (file)
@@ -82,6 +82,7 @@ class  AliCaloRawAnalyzer : public TObject
  protected:
   short Max( const AliCaloBunchInfo *const bunch, int *const maxindex) const;
   UShort_t Max(const UShort_t *data, const int length ) const;
+  bool CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch) const;
   bool IsInTimeRange( const int maxindex ) const;
   Float_t  ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t altrocfg1,  const UInt_t altrocfg2, double *outarray ) const;
   int  SelectBunch( const vector<AliCaloBunchInfo> &bunchvector, short *const maxampbin, short *const maxamplitude ) const;
index 2cf79c8..191a60d 100644 (file)
@@ -48,11 +48,11 @@ AliCaloRawAnalyzerCrude::Evaluate(const vector<AliCaloBunchInfo> &bunchvector, c
   // Evaluation of signal parameters
   if( bunchvector.size()  <=  0 )
     {
-      return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid , AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
+      return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
     }
 
   Int_t amp = 0;
-  Float_t tof = -99;
+  Int_t tof = -99;
   const UShort_t *sig;
   
   double ped = EvaluatePedestal( bunchvector.at(0).GetData(), bunchvector.at(0).GetLength() ) ;
@@ -72,7 +72,7 @@ AliCaloRawAnalyzerCrude::Evaluate(const vector<AliCaloBunchInfo> &bunchvector, c
 
   //:EvaluatePedestal(const UShort_t * const data, const int length )
   //  double ped = EvaluatePedestal(sig, length) ;
-  return  AliCaloFitResults(amp, ped, AliCaloFitResults::kNoFit, amp - ped, tof, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit );
+  return  AliCaloFitResults(amp, ped, AliCaloFitResults::kCrude, amp - ped, tof);
   
 } //end Crude
 
index ae3a20d..48e2b5e 100644 (file)
@@ -69,7 +69,6 @@ AliCaloRawAnalyzerFastFit::Evaluate( const vector<AliCaloBunchInfo> &bunchvector
       Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
       int maxrev =  maxampindex -  bunchvector.at(index).GetStartBin();
       short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
-
       if ( maxf > fAmpCut )
        {
          SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev , &first, &last);
@@ -87,26 +86,27 @@ AliCaloRawAnalyzerFastFit::Evaluate( const vector<AliCaloBunchInfo> &bunchvector
              Double_t eSignal = 1; // nominal 1 ADC error
              Double_t dAmp = maxf; 
              Double_t eAmp = 0;
-             Double_t dTime = 0;
+             Double_t dTime0 = 0;
              Double_t eTime = 0;
              Double_t chi2 = 0;
              Double_t dTau = 2.35; // time-bin units
              
              AliCaloFastAltroFitv0::FastFit(fXAxis, ordered , nsamples,
-                                            eSignal, dTau, dAmp, eAmp, dTime, eTime, chi2);
+                                            eSignal, dTau, dAmp, eAmp, dTime0, eTime, chi2);
           
-             return AliCaloFitResults(maxamp, ped, AliCaloFitResults::kDummy,  dAmp, dTime+timebinOffset,  chi2,  AliCaloFitResults::kDummy,
+             Double_t dTimeMax = dTime0 + timebinOffset - (maxrev - first) // abs. t0
+               + dTau; // +tau, makes sum tmax
+             return AliCaloFitResults(maxamp, ped, AliCaloFitResults::kFitPar,  dAmp, dTimeMax, timebinOffset, chi2,  AliCaloFitResults::kDummy,
                                       AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
            } // samplecut
          else 
            {
-             return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit,
-                                       AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); 
+             return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
            }
        } // ampcut
       else 
        {
-         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit);
+         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
        }
     } // bunch index    
 
index 0be5fd0..ddeadb9 100644 (file)
@@ -29,6 +29,7 @@
 #include "AliCaloFitResults.h"
 #include "AliLog.h"
 #include "TMath.h"
+#include <stdexcept>
 #include <iostream>
 #include "TF1.h"
 #include "TGraph.h"
@@ -87,8 +88,8 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
       int last;
       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 )
        {
          SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);
@@ -96,10 +97,13 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
          
          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(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);
@@ -109,18 +113,29 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
                fTf1->SetParameter(2, fTau);
              }
 
-             Short_t tmpStatus =  graph->Fit(fTf1, "Q0RW");
-            
+             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); 
+             }
+
              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,  
+               return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar,
                                          fTf1->GetParameter(0)/fkEulerSquared, 
-                                         fTf1->GetParameter(1) + timebinOffset,  
+                                         tmax,
+                                         timebinOffset,  
                                          fTf1->GetChisquare(), 
                                          fTf1->GetNDF(),
                                          AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
@@ -130,13 +145,12 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
            }
          else
            {
-             return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit,
-                                       AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); 
+             return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
            }
        }
       else
        {
-         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit);
+         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
        }       
     }
 
index ac4565a..5d5f9fb 100644 (file)
@@ -62,7 +62,7 @@ AliCaloRawAnalyzerNN::Evaluate( const vector<AliCaloBunchInfo> &bunchvector,
   // The eveluation of  Peak position and amplitude using the Neural Network
   if( bunchvector.size()  <=  0 )
     {
-      return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid , AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
+      return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
     } 
  
   short maxampindex;
@@ -90,8 +90,7 @@ AliCaloRawAnalyzerNN::Evaluate( const vector<AliCaloBunchInfo> &bunchvector,
     {
       if (  ( maxrev   - first) < 2  &&  (last -   maxrev ) < 2)
        {
-         return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit,
-                                   AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); 
+         return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
        }
       else
        {
@@ -105,13 +104,12 @@ AliCaloRawAnalyzerNN::Evaluate( const vector<AliCaloBunchInfo> &bunchvector,
          double amp = (maxamp - ped)*fNeuralNet->Value( 0,  fNNInput[0],  fNNInput[1], fNNInput[2], fNNInput[3], fNNInput[4]);
          double tof = (fNeuralNet->Value( 1,  fNNInput[0],  fNNInput[1], fNNInput[2], fNNInput[3], fNNInput[4]) + timebinOffset ) ;
 
-         return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kDummy, amp , tof, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
+         return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, amp , tof, timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
                                    AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 
        }
     }
-  return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit,
-                           AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); 
+  return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
 }
 
 
index 7891011..be45052 100644 (file)
@@ -133,7 +133,7 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvec
       if(  maxf < fAmpCut  ||  ( maxamp - ped) > 900  )         // (maxamp - ped) > 900 = Close to saturation (use low gain then)
        {
          //      cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped  << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex  << endl;
-         return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit );
+         return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
        }
       
       int first;
@@ -221,18 +221,17 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvec
              //      tof = tof/fAmpA[tmpindex];
 
              
-             return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kDummy, fAmpA[tmpindex], tof, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
+             return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, fAmpA[tmpindex], tof, timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
                                        AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );  
            }
          else
            {
-             return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit,
-                                       AliCaloFitResults::kNoFit, AliCaloFitSubarray(index, maxrev, first, last) ); 
+             return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kCrude, maxf, timebinOffset); 
            }
        }
       else 
        {
-         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kNoFit, maxf, maxrev+timebinOffset, AliCaloFitResults::kNoFit, AliCaloFitResults::kNoFit);
+         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
        }
     }
    //  cout << __FILE__ << __LINE__ <<  "WARNING, returning amp = -1 "  <<  endl;
index 09b92c7..4367eb7 100644 (file)
@@ -27,6 +27,7 @@
 //*-- Author: Marco van Leeuwen (LBL)
 
 #include "AliEMCALRawUtils.h"
+#include <stdexcept>
   
 #include "TF1.h"
 #include "TGraph.h"
@@ -376,27 +377,33 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
       } // loop over bunches
 
    
-    if ( caloFlag < 2 ){ // ALTRO
+      if ( caloFlag < 2 ){ // ALTRO
                
-         Float_t time = 0; 
-         Float_t amp = 0; 
+       Float_t time = 0; 
+       Float_t amp = 0; 
+       short timeEstimate = 0;
+       Float_t ampEstimate = 0;
+       Bool_t fitDone = kFALSE;
                
       if ( fFittingAlgorithm == kFastFit || fFittingAlgorithm == kNeuralNet || fFittingAlgorithm == kLMS || fFittingAlgorithm == kPeakFinder || fFittingAlgorithm == kCrude) {
        // all functionality to determine amp and time etc is encapsulated inside the Evaluate call for these methods 
        AliCaloFitResults fitResults = fRawAnalyzer->Evaluate( bunchlist, in.GetAltroCFG1(), in.GetAltroCFG2()); 
 
        amp = fitResults.GetAmp();
-       time = fitResults.GetTof();
+       time = fitResults.GetTime();
+       timeEstimate = fitResults.GetMaxTimebin();
+       ampEstimate = fitResults.GetMaxSig();
+       if (fitResults.GetStatus() == AliCaloFitResults::kFitPar) {
+         fitDone = kTRUE;
+       } 
       }
       else { // for the other methods we for now use the functionality of 
        // AliCaloRawAnalyzer as well, to select samples and prepare for fits, 
        // if it looks like there is something to fit
 
        // parameters init.
-       Float_t ampEstimate  = 0;
+       Float_t pedEstimate  = 0;
        short maxADC = 0;
-       short timeEstimate = 0;
-       Float_t pedEstimate = 0;
        Int_t first = 0;
        Int_t last = 0;
        Int_t bunchIndex = 0;
@@ -412,35 +419,35 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
          
          time = timeEstimate; // maxrev in AliCaloRawAnalyzer speak; comes with an offset w.r.t. real timebin
          Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1); 
-
          amp = ampEstimate; 
          
          if ( nsamples > 1 ) { // possibly something to fit
-           FitRaw(first, last, amp, time);
+           FitRaw(first, last, amp, time, fitDone);
            time += timebinOffset;
+           timeEstimate += timebinOffset;
          }
          
-         if ( amp>0 && time>0 ) { // brief sanity check of fit results
-           
-           // check fit results: should be consistent with initial estimates
-           // more magic numbers, but very loose cuts, for now..
-           // We have checked that amp and ampEstimate values are positive so division for assymmetry
-           // calculation should be OK/safe
-           Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
-           if ( (TMath::Abs(ampAsymm) > 0.1) ) {
-             AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations ped %f max-ped %f time %d",
-                             amp, time, pedEstimate, ampEstimate, timeEstimate));
-             
-             // what should do we do then? skip this channel or assign the simple estimate? 
-             // for now just overwrite the fit results with the simple estimate
-             amp = ampEstimate;
-             time = timeEstimate; 
-           } // asymm check
-         } // amp & time check
        } // ampEstimate check
       } // method selection
+
+      if ( fitDone ) { // brief sanity check of fit results        
+       Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
+       Float_t timeDiff = time - timeEstimate;
+       if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) ) {
+         // AliDebug(2,Form("Fit results amp %f time %f not consistent with expectations amp %f time %d", amp, time, ampEstimate, timeEstimate));
+         
+         // for now just overwrite the fit results with the simple/initial estimate
+         amp = ampEstimate;
+         time = timeEstimate; 
+         fitDone = kFALSE;
+       } 
+      } // fitDone
     
       if (amp > fNoiseThreshold  && amp<fgkRawSignalOverflow) { // something to be stored
+       if ( ! fitDone) { // smear ADC with +- 0.5 uniform (avoid discrete effects)
+         amp += (0.5 - gRandom->Rndm()); // Rndm generates a number in ]0,1]
+       }
+
        Int_t id =  fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ;
        lowGain = in.IsLowGain();
 
@@ -544,13 +551,14 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain
 }
 
 //____________________________________________________________________________ 
-void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time) const 
+void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const 
 { // Fits the raw signal time distribution
   
   //--------------------------------------------------
   //Do the fit, different fitting algorithms available
   //--------------------------------------------------
   int nsamples = lastTimeBin - firstTimeBin + 1;
+  fitDone = kFALSE;
 
   switch(fFittingAlgorithm) {
   case kStandard:
@@ -573,18 +581,27 @@ void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin,
       signalF->FixParameter(4, 0); // pedestal should be subtracted when we get here 
       signalF->SetParameter(1, time);
       signalF->SetParameter(0, amp);
-                               
-      gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
-                               
-      // assign fit results
-      amp = signalF->GetParameter(0); 
-      time = signalF->GetParameter(1);
-
+      // 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);
+
+       // cross-check with ParabolaFit to see if the results make sense
+       FitParabola(gSig, amp); // amp is possibly updated
+       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;
 
-      // cross-check with ParabolaFit to see if the results make sense
-      FitParabola(gSig, amp); // amp is possibly updated
-
       //printf("Std   : Amp %f, time %g\n",amp, time);
       delete gSig; // delete TGraph
                                
@@ -620,6 +637,7 @@ void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin,
       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);
index 6d45676..371c56c 100644 (file)
@@ -88,7 +88,7 @@ class AliEMCALRawUtils : public TObject {
 
   // Signal shape functions
        
-  void FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time) const ;
+  void FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const ;
   void FitParabola(const TGraph *gSig, Float_t & amp) const ; 
   static Double_t RawResponseFunction(Double_t *x, Double_t *par); 
   static Double_t RawResponseFunctionLog(Double_t *x, Double_t *par);