EMCAL/AliCaloRawAnalyzerKStandard.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 "AliCaloRawAnalyzerKStandard.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 #include "TRandom.h"
  37
  38 #include "AliEMCALRawResponse.h"
  39
  40
  41 using namespace std;
  42
  43 ClassImp( AliCaloRawAnalyzerKStandard )
  44
  45
  46 AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard")
  47 {
  48
  49   fAlgo = Algo::kStandard;
  50 }
  51
  52
  53 AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
  54 {
  55   //  delete fTf1;
  56 }
  57
  58
  59 AliCaloFitResults
  60 AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo>  &bunchlist, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  61 {
  62   //Evaluation Amplitude and TOF
  63   Float_t pedEstimate  = 0;
  64   short maxADC = 0;
  65   Int_t first = 0;
  66   Int_t last = 0;
  67   Int_t bunchIndex = 0;
  68   Float_t ampEstimate = 0;
  69   short timeEstimate  = 0;
  70   Float_t time = 0;
  71   Float_t amp=0;
  72   Float_t chi2 = 0;
  73   Int_t ndf = 0;
  74   Bool_t fitDone = kFALSE;
  75   int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
  76                                         maxADC, timeEstimate, pedEstimate, first, last,   (int)fAmpCut );
  77
  78
  79   if (ampEstimate >= fAmpCut  )
  80     {
  81       time = timeEstimate;
  82       Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
  83       amp = ampEstimate;
  84
  85       if ( nsamples > 1 && maxADC< OVERFLOWCUT )
  86         {
  87           FitRaw(first, last, amp, time, chi2, fitDone);
  88           time += timebinOffset;
  89           timeEstimate += timebinOffset;
  90           ndf = nsamples - 2;
  91         }
  92     }
  93   if ( fitDone )
  94     {
  95       Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
  96       Float_t timeDiff = time - timeEstimate;
  97
  98       if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) )
  99         {
 100           amp     = ampEstimate;
 101           time    = timeEstimate;
 102           fitDone = kFALSE;
 103         }
 104     }
 105   if (amp >= fAmpCut )
 106     {
 107       if ( ! fitDone)
 108         {
 109           amp += (0.5 - gRandom->Rndm());
 110         }
 111       time = time * TIMEBINWITH;
 112       time -= fL1Phase;
 113
 114       return AliCaloFitResults( -99, -99, fAlgo , amp, time,
 115                                 (int)time, chi2, ndf, Ret::kDummy );
 116      }
 117   return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
 118 }
 119
 120
 121 void
 122  AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
 123 {
 124   // Fits the raw signal time distribution
 125   int nsamples = lastTimeBin - firstTimeBin + 1;
 126   fitDone = kFALSE;
 127   if (nsamples < 3) { return; }
 128
 129   TGraph *gSig =  new TGraph( nsamples);
 130
 131   for (int i=0; i<nsamples; i++)
 132     {
 133       Int_t timebin = firstTimeBin + i;
 134       gSig->SetPoint(i, timebin, GetReversed(timebin));
 135     }
 136
 137   TF1 * signalF = new TF1("signal", AliEMCALRawResponse::RawResponseFunction, 0, TIMEBINS , 5);
 138
 139   signalF->SetParameters(10.,5., TAU  ,ORDER,0.); //set all defaults once, just to be safe
 140   signalF->SetParNames("amp","t0","tau","N","ped");
 141   signalF->FixParameter(2,TAU);
 142   signalF->FixParameter(3,ORDER);
 143   signalF->FixParameter(4, 0);
 144   signalF->SetParameter(1, time);
 145   signalF->SetParameter(0, amp);
 146   signalF->SetParLimits(0, 0.5*amp, 2*amp );
 147   signalF->SetParLimits(1, time - 4, time + 4);
 148
 149   try {
 150     gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
 151     amp  = signalF->GetParameter(0);
 152     time = signalF->GetParameter(1);
 153     chi2 = signalF->GetChisquare();
 154     fitDone = kTRUE;
 155   }
 156   catch (const std::exception & e)
 157     {
 158       AliError( Form("TGraph Fit exception %s", e.what()) );
 159       // stay with default amp and time in case of exception, i.e. no special action required
 160       fitDone = kFALSE;
 161   }
 162
 163   delete signalF;
 164   delete gSig; // delete TGraph
 165   return;
 166 }
 167
 168