several changes: 1) check for overflow before fitting, not using fit results 2) keep...
authordsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 May 2010 12:30:29 +0000 (12:30 +0000)
committerdsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 26 May 2010 12:30:29 +0000 (12:30 +0000)
EMCAL/AliCaloRawAnalyzer.cxx
EMCAL/AliCaloRawAnalyzer.h
EMCAL/AliCaloRawAnalyzerCrude.cxx
EMCAL/AliCaloRawAnalyzerFastFit.cxx
EMCAL/AliCaloRawAnalyzerLMS.cxx
EMCAL/AliCaloRawAnalyzerNN.cxx
EMCAL/AliCaloRawAnalyzerPeakFinder.cxx
EMCAL/AliEMCALDigit.cxx
EMCAL/AliEMCALDigit.h
EMCAL/AliEMCALRawUtils.cxx
EMCAL/AliEMCALRawUtils.h

index 2c68461..cae7b21 100644 (file)
@@ -40,6 +40,7 @@ AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort)
                                                                                   fFitArrayCut(5),
                                                                                   fAmpCut(4),
                                                                                   fNsampleCut(5),
+                                                                                  fOverflowCut(950),
                                                                                   fNsamplePed(3),
                                                                                   fIsZerosupressed( false ),
                                                                                   fVerbose( false )
@@ -335,7 +336,7 @@ AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
   }
 
   int nsamples =  last - first + 1;
-  //  printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i \n", first, last, nsamples); 
+  // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time); 
 
   Int_t x = 0;
   Double_t chi2 = 0;
@@ -350,7 +351,7 @@ AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
     }
     dy    = fReversed[x] - f; 
     chi2 += dy*dy;
-    //    printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); 
+    // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); 
   }
 
   if (adcErr>0.0) { // weight chi2
index a9df4e6..7791c59 100644 (file)
@@ -54,12 +54,14 @@ class  AliCaloRawAnalyzer : public TObject
   void SetAmpCut(const Float_t cut) { fAmpCut = cut ; } ;
   void SetFitArrayCut(const Int_t cut) { fFitArrayCut = cut ; } ;
   void SetNsampleCut(const Int_t cut) { fNsampleCut = cut ; } ;
+  void SetOverflowCut(const Int_t cut) { fOverflowCut = cut ; } ;
   void SetNsamplePed(const Int_t i) { fNsamplePed = i ; } ;
 
   bool GetIsZeroSuppressed() const { return fIsZerosupressed;} ;
   Float_t GetAmpCut() const { return fAmpCut; } ;
   Int_t GetFitArrayCut() const { return fFitArrayCut; } ;
   Int_t GetNsampleCut() const { return fNsampleCut; } ;
+  Int_t GetOverflowCut() const { return fOverflowCut; } ;
   Int_t GetNsamplePed() const { return fNsamplePed; } ;
 
   // access to array info
@@ -93,6 +95,7 @@ class  AliCaloRawAnalyzer : public TObject
   int fFitArrayCut;  //Cut on ADC value (after ped. subtraction) for signals used for fit
   Float_t fAmpCut;   //Max ADC - pedestal must be higher than this befor attemting to extract the amplitude 
   int fNsampleCut;   //Minimum number of sample require before attemting to extract signal parameters 
+  int fOverflowCut; // value when ADC starts to saturate
   int fNsamplePed;   //Number of samples used for pedestal calculation (first in bunch) 
   bool fIsZerosupressed; //Wether or not the data is zeros supressed, by default its assumed that the baseline is also subtracted if set to true
   bool fVerbose;     //Print debug information to std out if set to true
index 01ba110..979b2e2 100644 (file)
@@ -25,6 +25,7 @@
 #include "AliCaloRawAnalyzerCrude.h"
 #include "AliCaloFitResults.h"
 #include "AliCaloBunchInfo.h"
+#include "TMath.h"
 
 using namespace std;
 
@@ -44,37 +45,39 @@ AliCaloRawAnalyzerCrude::~AliCaloRawAnalyzerCrude()
 
 
 AliCaloFitResults
-AliCaloRawAnalyzerCrude::Evaluate(const vector<AliCaloBunchInfo> &bunchvector, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/)
+AliCaloRawAnalyzerCrude::Evaluate(const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2)
 {
   // Evaluation of signal parameters
-  if( bunchvector.size()  <=  0 )
+  short maxampindex; //index of maximum amplitude
+  short maxamp; //Maximum amplitude
+  int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
+  if( index >= 0)
     {
-      return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
-    }
+      Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
+      Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
+      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, AliCaloFitResults::kCrude, maxf, timebinOffset);
+       }
+      else if ( maxf >= fAmpCut ) // no if statement needed really; keep for readability
+       {
+         int first = 0;
+         int last = 0;
+         int maxrev =  maxampindex -  bunchvector.at(index).GetStartBin();
+         SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev , &first, &last);
+
+         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) ); 
+       } // ampcut
+    } // bunch index    
+
+  return AliCaloFitResults(AliCaloFitResults::kInvalid , AliCaloFitResults::kInvalid);
 
-  Int_t amp = 0;
-  Int_t tof = -99;
-  const UShort_t *sig;
-  
-  double ped = EvaluatePedestal( bunchvector.at(0).GetData(), bunchvector.at(0).GetLength() ) ;
-
-  for( unsigned int i= 0; i < bunchvector.size(); ++i)
-    {
-      sig = bunchvector.at(i).GetData();
-      int length = bunchvector.at(i).GetLength(); 
-      
-      for(int j = 0; j < length; j ++)
-       if( sig[j] > amp  )
-         {
-           amp   = sig[j];
-           tof   = bunchvector.at(i).GetStartBin() - j;                     
-         }
-    }
-
-  //:EvaluatePedestal(const UShort_t * const data, const int length )
-  //  double ped = EvaluatePedestal(sig, length) ;
-  return  AliCaloFitResults(amp, ped, AliCaloFitResults::kCrude, amp - ped, tof);
-  
 } //end Crude
 
 
index 4a7d1d8..86d68d7 100644 (file)
@@ -64,13 +64,19 @@ AliCaloRawAnalyzerFastFit::Evaluate( const vector<AliCaloBunchInfo> &bunchvector
   if( index >= 0)
     {
       Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
-      int first = 0;
-      int last = 0;
       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 )
+
+      if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
+       {
+         return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
+       }
+      else if ( maxf >= fAmpCut ) // no if statement needed really; keep for readability
        {
+         int first = 0;
+         int last = 0;
+         int maxrev =  maxampindex -  bunchvector.at(index).GetStartBin();
+
          SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev , &first, &last);
          int nsamples =  last - first + 1;
 
@@ -101,13 +107,12 @@ AliCaloRawAnalyzerFastFit::Evaluate( const vector<AliCaloBunchInfo> &bunchvector
            } // samplecut
          else 
            {
-             return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
+             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) ); 
            }
        } // ampcut
-      else 
-       {
-         return AliCaloFitResults( maxamp , ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
-       }
     } // bunch index    
 
   return AliCaloFitResults(AliCaloFitResults::kInvalid , AliCaloFitResults::kInvalid);
index 455ec98..e220170 100644 (file)
@@ -84,14 +84,18 @@ 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 );
       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 )
+      if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
        {
+         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;
          
@@ -119,7 +123,8 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
              }
              catch (const std::exception & e) {
                AliError( Form("TGraph Fit exception %s", e.what()) ); 
-               return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, timebinOffset); 
+               return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kNoFit, maxf, timebinOffset,
+                                         timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
              }
 
              if( fVerbose == true )
@@ -145,15 +150,13 @@ AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, c
            }
          else
            {
-             return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
+             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, AliCaloFitResults::kCrude, maxf, timebinOffset);
-       }       
+        } // ampcut
     }
-
   return AliCaloFitResults( AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid );
   
 }
index 5d5f9fb..606fe36 100644 (file)
@@ -75,22 +75,30 @@ AliCaloRawAnalyzerNN::Evaluate( const vector<AliCaloBunchInfo> &bunchvector,
       return AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
     }
   
-  int first = 0;
-  int last = 0;
-  Float_t ped = ReverseAndSubtractPed( &(bunchvector.at( index ) )  ,  altrocfg1, altrocfg2, fReversed  );
-  
-  short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
+  Float_t ped = ReverseAndSubtractPed( &(bunchvector.at( index ) )  ,  altrocfg1, altrocfg2, fReversed  );  
   short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
   double maxf =  maxamp - ped;
 
+  if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
+    {
+      return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
+    }
+
+  int first = 0;
+  int last = 0; 
+  short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
   SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev , &first, &last);
 
+  Float_t chi2 = 0;
+  Int_t ndf = 0;
   if(maxrev  < 1000 )
     {
       if (  ( maxrev   - first) < 2  &&  (last -   maxrev ) < 2)
        {
-         return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
+         chi2 = CalculateChi2(maxf, maxrev, first, last);
+         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
        {
@@ -104,12 +112,19 @@ 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::kFitPar, amp , tof, timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
+         // use local-array time for chi2 estimate
+         chi2 = CalculateChi2(amp, tof-timebinOffset+maxrev, first, last);
+         ndf = last - first - 1; // nsamples - 2
+         return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, amp , tof, timebinOffset, chi2, ndf,
                                    AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 
        }
     }
-  return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset); 
+  chi2 = CalculateChi2(maxf, maxrev, first, last);
+  ndf = last - first - 1; // nsamples - 2
+  return AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset,
+                           timebinOffset, chi2, ndf, AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) ); 
+
 }
 
 
index beb6b3a..e40d505 100644 (file)
@@ -138,22 +138,21 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvec
     {
       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();
       short timebinOffset = maxampindex - (bunchvector.at( index ).GetLength()-1); 
             
-      if(  maxf <= fAmpCut  ||  ( maxamp - ped) > 900  )        // (maxamp - ped) > 900 = Close to saturation (use low gain then)
+      if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
        {
          return  AliCaloFitResults( maxamp, ped, AliCaloFitResults::kCrude, maxf, timebinOffset);
-       }
-      
-      int first;
-      int last;
-      
-      if ( maxf >= fAmpCut )
-       {         
+       }            
+      else if ( maxf >= fAmpCut )
+       {
+         int first = 0;
+         int last = 0;
+         short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();     
          SelectSubarray( fReversed,  bunchvector.at(index).GetLength(), maxrev, &first, &last);
          int nsamples =  last - first;
-         if( ( nsamples  )  >= fNsampleCut )
+
+         if( ( nsamples  )  >= fNsampleCut ) // no if statement needed really; keep for readability
            {
              int startbin = bunchvector.at(index).GetStartBin();  
              int n = last - first;  
@@ -199,19 +198,21 @@ AliCaloRawAnalyzerPeakFinder::Evaluate( const vector<AliCaloBunchInfo> &bunchvec
              
              tof = timebinOffset - 0.01*tof/fAmp; // clock ticks
              
+             // use local-array time for chi2 estimate
+             Float_t chi2 = CalculateChi2(fAmp, tof-timebinOffset+maxrev, first, last);
+             Int_t ndf = last - first - 1; // nsamples - 2
              return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kFitPar, fAmp, tof, 
-                                       timebinOffset, AliCaloFitResults::kDummy, AliCaloFitResults::kDummy,
+                                       timebinOffset, chi2, ndf,
                                        AliCaloFitResults::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );  
            }
          else
            {
-             return AliCaloFitResults( maxamp, ped , AliCaloFitResults::kCrude, maxf, timebinOffset); 
+             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, AliCaloFitResults::kCrude, maxf, timebinOffset);
-       }
+       } // ampcut
     }
   return  AliCaloFitResults(AliCaloFitResults::kInvalid, AliCaloFitResults::kInvalid);
 }
index 9856a69..74b4f42 100644 (file)
@@ -82,7 +82,7 @@ AliDigitNew(),
 }
 
 //____________________________________________________________________________
-AliEMCALDigit::AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t digEnergy, Float_t time, Int_t type, Int_t index, Float_t dE) 
+AliEMCALDigit::AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t digEnergy, Float_t time, Int_t type, Int_t index, Float_t chi2, Int_t ndf, Float_t dE) 
   : AliDigitNew(),
     fAmpFloat(digEnergy),
     fNSamples(0),
@@ -100,8 +100,8 @@ AliEMCALDigit::AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t dig
     fMaxIter(5),
     fTime(time),
     fTimeR(time),
-    fChi2(0.),
-    fNDF(0),
+    fChi2(chi2),
+    fNDF(ndf),
     fDigitType(type)
 {  
   // ctor with all data 
index 00effe8..37a52e3 100644 (file)
@@ -31,7 +31,7 @@ class AliEMCALDigit : public AliDigitNew {
  public:
   
   AliEMCALDigit() ;
-  AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t digEnergy, Float_t time, Int_t type,Int_t index = -1, Float_t dE = 0) ;
+  AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t digEnergy, Float_t time, Int_t type,Int_t index = -1, Float_t chi2=0, Int_t ndf=0, Float_t dE = 0) ;
   AliEMCALDigit(const AliEMCALDigit & digit) ;
   virtual ~AliEMCALDigit() ;
 
index 6d12115..77d2ee7 100644 (file)
@@ -340,7 +340,8 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
   reader->Select("EMCAL",0,43); // 43 = AliEMCALGeoParams::fgkLastAltroDDL
 
   // fRawAnalyzer setup
-  fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done
+  fRawAnalyzer->SetNsampleCut(5); // requirement for fits to be done, for the new methods
+  fRawAnalyzer->SetOverflowCut(fgkOverflowCut);
   fRawAnalyzer->SetAmpCut(fNoiseThreshold);
   fRawAnalyzer->SetFitArrayCut(fNoiseThreshold);
   fRawAnalyzer->SetIsZeroSuppressed(true); // TMP - should use stream->IsZeroSuppressed(), or altro cfg registers later
@@ -356,16 +357,6 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
 
     while (in.NextChannel()) {
 
-/*
-         Int_t    hhwAdd    = in.GetHWAddress();
-         UShort_t iiBranch  = ( hhwAdd >> 11 ) & 0x1; // 0/1
-         UShort_t iiFEC     = ( hhwAdd >>  7 ) & 0xF;
-         UShort_t iiChip    = ( hhwAdd >>  4 ) & 0x7;
-         UShort_t iiChannel =   hhwAdd         & 0xF;
-                
-         if ( !( iiBranch == 0 && iiFEC == 1 && iiChip == 3 && ( iiChannel >= 8 && iiChannel <= 15 ) ) && !( iiBranch == 1 && iiFEC == 0 && in.GetColumn() == 0 ) ) continue;
-*/
-               
       //Check if the signal  is high or low gain and then do the fit, 
       //if it  is from TRU or LEDMon do not fit
       caloFlag = in.GetCaloFlag();
@@ -391,7 +382,9 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
        short timeEstimate  = 0;
        Float_t ampEstimate = 0;
        Bool_t fitDone = kFALSE;
-               
+       Float_t chi2 = 0;
+       Int_t ndf = 0;
+
       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()); 
@@ -400,9 +393,11 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
        time         = fitResults.GetTime();
        timeEstimate = fitResults.GetMaxTimebin();
        ampEstimate  = fitResults.GetMaxSig();
+       chi2 = fitResults.GetChi2();
+       ndf = fitResults.GetNdf();
        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, 
@@ -428,10 +423,11 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
          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, fitDone);
+         if ( nsamples > 1 && maxADC<fgkOverflowCut ) { // possibly something to fit
+           FitRaw(first, last, amp, time, chi2, fitDone);
            time += timebinOffset;
            timeEstimate += timebinOffset;
+           ndf = nsamples - 2;
          }
          
        } // ampEstimate check
@@ -450,7 +446,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
        } 
       } // fitDone
     
-      if (amp >= fNoiseThreshold  && amp<fgkRawSignalOverflow) { // something to be stored
+      if (amp >= fNoiseThreshold) { // 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]
        }
@@ -465,7 +461,7 @@ void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr,
 
        AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp));
        // printf("Added tower: SM %d, row %d, column %d, amp %3.2f\n",in.GetModule(), in.GetRow(), in.GetColumn(),amp);
-       AddDigit(digitsArr, id, lowGain, amp, time); 
+       AddDigit(digitsArr, id, lowGain, amp, time, chi2, ndf); 
       }
       
        }//ALTRO
@@ -525,7 +521,7 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t timeSam
 }
 
 //____________________________________________________________________________ 
-void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time) {
+void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf) {
   //
   // Add a new digit. 
   // This routine checks whether a digit exists already for this tower 
@@ -546,7 +542,7 @@ void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain
                        type = AliEMCALDigit::kLGnoHG;
                } 
                Int_t idigit = digitsArr->GetEntries();
-               new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit) ; 
+               new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, type, idigit, chi2, ndf); 
                AliDebug(2,Form("Add digit Id %d for the first time, type %d", id, type));
   }//digit added first time
   else { // a digit already exists, check range 
@@ -591,11 +587,16 @@ void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
       AliDebug(1,Form("Remove digit with id %d, LGnoHG",digit->GetId()));
       digitsArr->Remove(digit);
     }
-    //Check if time if too large or too small, remove if so
+    //Check if time is too large or too small, remove if so
     else if(fTimeMin > digit->GetTime() || fTimeMax < digit->GetTime()) {
       digitsArr->Remove(digit);
       AliDebug(1,Form("Remove digit with id %d, Bad Time %e",digit->GetId(), digit->GetTime()));
     }
+    // Check if Chi2 is undefined
+    else if (0 > digit->GetChi2()) {
+      digitsArr->Remove(digit);
+      AliDebug(1,Form("Remove digit with id %d, Bad Chi2 %e",digit->GetId(), digit->GetChi2()));
+    }
     //Good digit, just reassign the index of the digit in case there was a previous removal
     else {
       digit->SetIndexInList(n);        
@@ -609,7 +610,7 @@ void AliEMCALRawUtils::TrimDigits(TClonesArray *digitsArr)
 }
        
 //____________________________________________________________________________ 
-void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const 
+void AliEMCALRawUtils::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
   
   //--------------------------------------------------
@@ -648,9 +649,7 @@ void AliEMCALRawUtils::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin,
        // 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
+       chi2 = signalF->GetChisquare();
        fitDone = kTRUE;
       }
       catch (const std::exception & e) {
index ee92e6d..c6940d6 100644 (file)
@@ -43,7 +43,7 @@ class AliEMCALRawUtils : public TObject {
   void Raw2Digits(AliRawReader *reader, TClonesArray *digitsArr, const AliCaloCalibPedestal* pedbadmap,
                                  TClonesArray *digitsTRG=0x0);
 
-  void AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time);
+  void AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Float_t amp, Float_t time, Float_t chi2, Int_t ndf);
   void AddDigit(TClonesArray *digitsArr, Int_t id, Int_t timeSamples[], Int_t nSamples);
   void TrimDigits(TClonesArray *digitsArr);
 
@@ -95,7 +95,7 @@ class AliEMCALRawUtils : public TObject {
 
   // Signal shape functions
        
-  void FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Bool_t & fitDone) const ;
+  void FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, 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);