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 using namespace std;
  39
  40 ClassImp( AliCaloRawAnalyzerKStandard )
  41
  42
  43 AliCaloRawAnalyzerKStandard::AliCaloRawAnalyzerKStandard() : AliCaloRawAnalyzerFitter("Chi Square ( kStandard )", "KStandard")
  44 {
  45
  46   fAlgo = Algo::kStandard;
  47 }
  48
  49
  50 AliCaloRawAnalyzerKStandard::~AliCaloRawAnalyzerKStandard()
  51 {
  52   //  delete fTf1;
  53 }
  54
  55
  56 AliCaloFitResults
  57 AliCaloRawAnalyzerKStandard::Evaluate( const vector<AliCaloBunchInfo>  &bunchlist, const UInt_t altrocfg1,  const UInt_t altrocfg2 )
  58 {
  59   Float_t pedEstimate  = 0;
  60   short maxADC = 0;
  61   Int_t first = 0;
  62   Int_t last = 0;
  63   Int_t bunchIndex = 0;
  64   Float_t ampEstimate = 0;
  65   short timeEstimate  = 0;
  66   Float_t time = 0;
  67   Float_t amp=0;
  68   Float_t chi2 = 0;
  69   Int_t ndf = 0;
  70   Bool_t fitDone = kFALSE;
  71   int nsamples = PreFitEvaluateSamples( bunchlist, altrocfg1, altrocfg2, bunchIndex, ampEstimate,
  72                                         maxADC, timeEstimate, pedEstimate, first, last,   (int)fAmpCut );
  73
  74
  75   if (ampEstimate >= fAmpCut  )
  76     {
  77       time = timeEstimate;
  78       Int_t timebinOffset = bunchlist.at(bunchIndex).GetStartBin() - (bunchlist.at(bunchIndex).GetLength()-1);
  79       amp = ampEstimate;
  80
  81       if ( nsamples > 1 && maxADC< OVERFLOWCUT )
  82         {
  83           FitRaw(first, last, amp, time, chi2, fitDone);
  84           time += timebinOffset;
  85           timeEstimate += timebinOffset;
  86           ndf = nsamples - 2;
  87         }
  88     }
  89   if ( fitDone )
  90     {
  91       Float_t ampAsymm = (amp - ampEstimate)/(amp + ampEstimate);
  92       Float_t timeDiff = time - timeEstimate;
  93
  94       if ( (TMath::Abs(ampAsymm) > 0.1) || (TMath::Abs(timeDiff) > 2) )
  95         {
  96           amp     = ampEstimate;
  97           time    = timeEstimate;
  98           fitDone = kFALSE;
  99         }
 100     }
 101   if (amp >= fAmpCut )
 102     {
 103       if ( ! fitDone)
 104         {
 105           amp += (0.5 - gRandom->Rndm());
 106         }
 107       time = time * TIMEBINWITH;
 108       time -= fL1Phase;
 109
 110       return AliCaloFitResults( -99, -99, fAlgo , amp, time,
 111                                 (int)time, chi2, ndf, Ret::kDummy );
 112      }
 113   return AliCaloFitResults( Ret::kInvalid, Ret::kInvalid );
 114 }
 115
 116
 117 //____________________________________________________________________________
 118 void
 119  AliCaloRawAnalyzerKStandard::FitRaw(const Int_t firstTimeBin, const Int_t lastTimeBin, Float_t & amp, Float_t & time, Float_t & chi2, Bool_t & fitDone) const
 120 {
 121   // Fits the raw signal time distribution
 122   int nsamples = lastTimeBin - firstTimeBin + 1;
 123   fitDone = kFALSE;
 124   if (nsamples < 3) { return; }
 125
 126   TGraph *gSig =  new TGraph( nsamples);
 127
 128   for (int i=0; i<nsamples; i++)
 129     {
 130       Int_t timebin = firstTimeBin + i;
 131       gSig->SetPoint(i, timebin, GetReversed(timebin));
 132     }
 133
 134   TF1 * signalF = new TF1("signal", RawResponseFunction, 0, TIMEBINS , 5);
 135   signalF->SetParameters(10.,5., TAU  ,ORDER,0.); //set all defaults once, just to be safe
 136   signalF->SetParNames("amp","t0","tau","N","ped");
 137   signalF->FixParameter(2,TAU);
 138   signalF->FixParameter(3,ORDER);
 139   signalF->FixParameter(4, 0);
 140   signalF->SetParameter(1, time);
 141   signalF->SetParameter(0, amp);
 142   signalF->SetParLimits(0, 0.5*amp, 2*amp );
 143   signalF->SetParLimits(1, time - 4, time + 4);
 144
 145   try {
 146     gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points
 147     amp  = signalF->GetParameter(0);
 148     time = signalF->GetParameter(1);
 149     chi2 = signalF->GetChisquare();
 150     fitDone = kTRUE;
 151   }
 152   catch (const std::exception & e) {
 153     AliError( Form("TGraph Fit exception %s", e.what()) );
 154     // stay with default amp and time in case of exception, i.e. no special action required
 155     fitDone = kFALSE;
 156   }
 157
 158   delete signalF;
 159   delete gSig; // delete TGraph
 160   return;
 161 }
 162
 163
 164 //__________________________________________________________________
 165 void
 166 AliCaloRawAnalyzerKStandard::FitParabola(const TGraph *gSig, Float_t & amp) const
 167 {
 168   //BEG YS alternative methods to calculate the amplitude
 169   Double_t * ymx = gSig->GetX() ;
 170   Double_t * ymy = gSig->GetY() ;
 171   const Int_t kN = 3 ;
 172   Double_t ymMaxX[kN] = {0., 0., 0.} ;
 173   Double_t ymMaxY[kN] = {0., 0., 0.} ;
 174   Double_t ymax = 0. ;
 175   // find the maximum amplitude
 176   Int_t ymiMax = 0 ;
 177   for (Int_t ymi = 0; ymi < gSig->GetN(); ymi++) {
 178     if (ymy[ymi] > ymMaxY[0] ) {
 179       ymMaxY[0] = ymy[ymi] ; //<========== This is the maximum amplitude
 180       ymMaxX[0] = ymx[ymi] ;
 181       ymiMax = ymi ;
 182     }
 183   }
 184   // find the maximum by fitting a parabola through the max and the two adjacent samples
 185   if ( ymiMax < gSig->GetN()-1 && ymiMax > 0) {
 186     ymMaxY[1] = ymy[ymiMax+1] ;
 187     ymMaxY[2] = ymy[ymiMax-1] ;
 188     ymMaxX[1] = ymx[ymiMax+1] ;
 189     ymMaxX[2] = ymx[ymiMax-1] ;
 190     if (ymMaxY[0]*ymMaxY[1]*ymMaxY[2] > 0) {
 191       //fit a parabola through the 3 points y= a+bx+x*x*x
 192       Double_t sy = 0 ;
 193       Double_t sx = 0 ;
 194       Double_t sx2 = 0 ;
 195       Double_t sx3 = 0 ;
 196       Double_t sx4 = 0 ;
 197       Double_t sxy = 0 ;
 198       Double_t sx2y = 0 ;
 199       for (Int_t i = 0; i < kN ; i++) {
 200         sy += ymMaxY[i] ;
 201         sx += ymMaxX[i] ;
 202         sx2 += ymMaxX[i]*ymMaxX[i] ;
 203         sx3 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
 204         sx4 += ymMaxX[i]*ymMaxX[i]*ymMaxX[i]*ymMaxX[i] ;
 205         sxy += ymMaxX[i]*ymMaxY[i] ;
 206         sx2y += ymMaxX[i]*ymMaxX[i]*ymMaxY[i] ;
 207       }
 208       Double_t cN = (sx2y*kN-sy*sx2)*(sx3*sx-sx2*sx2)-(sx2y*sx-sxy*sx2)*(sx3*kN-sx*sx2);
 209       Double_t cD = (sx4*kN-sx2*sx2)*(sx3*sx-sx2*sx2)-(sx4*sx-sx3*sx2)*(sx3*kN-sx*sx2) ;
 210       Double_t c  = cN / cD ;
 211       Double_t b  = ((sx2y*kN-sy*sx2)-c*(sx4*kN-sx2*sx2))/(sx3*kN-sx*sx2) ;
 212       Double_t a  = (sy-b*sx-c*sx2)/kN  ;
 213       Double_t xmax = -b/(2*c) ;
 214       ymax = a + b*xmax + c*xmax*xmax ;//<========== This is the maximum amplitude
 215       amp = ymax;
 216     }
 217   }
 218
 219   Double_t diff = TMath::Abs(1-ymMaxY[0]/amp) ;
 220   if (diff > 0.1)
 221     {
 222       amp = ymMaxY[0] ;
 223     }
 224
 225   return;
 226 }
 227
 228
 229 //__________________________________________________________________
 230 Double_t
 231 AliCaloRawAnalyzerKStandard::RawResponseFunction(Double_t *x, Double_t *par)
 232 {
 233   Double_t signal = 0.;
 234   Double_t tau    = par[2];
 235   Double_t n      = par[3];
 236   Double_t ped    = par[4];
 237   Double_t xx     = ( x[0] - par[1] + tau ) / tau ;
 238
 239   if (xx <= 0)
 240     signal = ped ;
 241   else {
 242     signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
 243   }
 244   return signal ;
 245 }
 246