EMCAL/AliCaloRawAnalyzerLMS.cxx

   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 "AliCaloRawAnalyzerLMS.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
  37 using namespace std;
  38
  39 ClassImp( AliCaloRawAnalyzerLMS )
  40
  41 AliCaloRawAnalyzerLMS::AliCaloRawAnalyzerLMS() : AliCaloRawAnalyzerFitter("Chi Square Fit", "LMS")
  42 {
  43   fAlgo = Algo::kLMS;
  44 }
  45
  46
  47 AliCaloRawAnalyzerLMS::~AliCaloRawAnalyzerLMS()
  48 {
  49   // delete fTf1;
  50 }
  51
  52
  53 AliCaloFitResults
  54 AliCaloRawAnalyzerLMS::Evaluate( const vector<AliCaloBunchInfo>  &bunchvector, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  55 {
  56   // Extracting signal parameters using fitting
  57   short maxampindex; //index of maximum amplitude
  58   short maxamp; //Maximum amplitude
  59   int index = SelectBunch( bunchvector,  &maxampindex,  &maxamp );
  60
  61   if( index >= 0)
  62     {
  63       Float_t  ped  = ReverseAndSubtractPed( &(bunchvector.at(index))  ,  altrocfg1, altrocfg2, fReversed  );
  64       Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(),  fReversed );
  65       short maxrev = maxampindex  -  bunchvector.at(index).GetStartBin();
  66       // timebinOffset is timebin value at maximum (maxrev)
  67       short timebinOffset = maxampindex - (bunchvector.at(index).GetLength()-1);
  68       if(  maxf < fAmpCut  ||  ( maxamp - ped) > fOverflowCut  ) // (maxamp - ped) > fOverflowCut = Close to saturation (use low gain then)
  69         {
  70           return  AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset);
  71         }
  72       else if ( maxf >= fAmpCut )
  73         {
  74           int first = 0;
  75           int last = 0;
  76           SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last, fFitArrayCut);
  77           int nsamples =  last - first + 1;
  78
  79           if( ( nsamples  )  >= fNsampleCut )
  80             {
  81               Float_t tmax = (maxrev - first); // local tmax estimate
  82               TGraph *graph =  new TGraph(  nsamples, fXaxis,  &fReversed[first] );
  83               fTf1->SetParameter(0, maxf*fkEulerSquared );
  84               fTf1->SetParameter(1, tmax - fTau);
  85               // set rather loose parameter limits
  86               fTf1->SetParLimits(0, 0.5*maxf*fkEulerSquared, 2*maxf*fkEulerSquared );
  87               fTf1->SetParLimits(1, tmax - fTau - 4, tmax - fTau + 4);
  88
  89               if (fFixTau) {
  90                 fTf1->FixParameter(2, fTau);
  91               }
  92               else {
  93                 fTf1->ReleaseParameter(2); // allow par. to vary
  94                 fTf1->SetParameter(2, fTau);
  95               }
  96
  97               Short_t tmpStatus = 0;
  98               try {
  99                 tmpStatus =  graph->Fit(fTf1, "Q0RW");
 100               }
 101               catch (const std::exception & e) {
 102                 AliError( Form("TGraph Fit exception %s", e.what()) );
 103                 return AliCaloFitResults( maxamp, ped, Ret::kNoFit, maxf, timebinOffset,
 104                                           timebinOffset, Ret::kDummy, Ret::kDummy, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 105               }
 106
 107               if( fVerbose == true )
 108                 {
 109                   AliCaloRawAnalyzer::PrintBunch( bunchvector.at(index) );
 110                   PrintFitResult( fTf1 ) ;
 111                 }
 112               // global tmax
 113               tmax = fTf1->GetParameter(1) + timebinOffset - (maxrev - first) // abs. t0
 114                 + fTf1->GetParameter(2); // +tau, makes sum tmax
 115
 116                 delete graph;
 117                 return AliCaloFitResults( maxamp, ped , Ret::kFitPar,
 118                                           fTf1->GetParameter(0)/fkEulerSquared,
 119                                           tmax,
 120                                           timebinOffset,
 121                                           fTf1->GetChisquare(),
 122                                           fTf1->GetNDF(),
 123                                           Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 124                 //     delete graph;
 125
 126             }
 127           else
 128             {
 129
 130               Float_t chi2 = CalculateChi2(maxf, maxrev, first, last);
 131               Int_t ndf = last - first - 1; // nsamples - 2
 132               return AliCaloFitResults( maxamp, ped, Ret::kCrude, maxf, timebinOffset,
 133                                         timebinOffset, chi2, ndf, Ret::kDummy, AliCaloFitSubarray(index, maxrev, first, last) );
 134             }
 135         } // ampcut
 136     }
 137   return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
 138
 139 }
 140