PWG2/RESONANCES/AliRsnFitResult.cxx

   1 //
   2 //
   3 //
   4
   5 #include <TF1.h>
   6 #include <TH1.h>
   7 #include <TMath.h>
   8
   9 #include "AliLog.h"
  10 #include <AliRsnFitResult.h>
  11
  12 ClassImp(AliRsnFitResult)
  13
  14 //__________________________________________________________________________________________________
  15 AliRsnFitResult::AliRsnFitResult(const char *name, ESignalType sType, EBackgroundType bType) :
  16   TNamed(name, ""),
  17   fSignalType(sType),
  18   fBackgroundType(bType)
  19 {
  20 //
  21 // Default constructor
  22 //
  23
  24   InitFunctions();
  25 }
  26
  27 //__________________________________________________________________________________________________
  28 Double_t AliRsnFitResult::NormSquareRoot(Double_t *x, Double_t *par)
  29 {
  30 //
  31 // Computes a square-root like function normalized
  32 // to the integral in the background range specified
  33 // for this object.
  34 // This is obtained by dividing the function value by
  35 // its integral analytically computed in that range.
  36 // Returns 0 in all cases where a floating point exception
  37 // could be raised by computing square root of a negative number.
  38 //
  39
  40   if (x[0] < par[1]) return 0.0;
  41
  42   Double_t fcnVal = TMath::Sqrt(x[0] - par[1]);
  43
  44   Double_t fcnMax = fPeakRange[1] - par[1];
  45   Double_t fcnMin = fPeakRange[0] - par[1];
  46   Double_t fcnInt = (2.0 / 3.0) * ((fcnMax*TMath::Sqrt(fcnMax) - fcnMin*TMath::Sqrt(fcnMin)));
  47
  48   return par[0] * fcnVal / fcnInt;
  49 }
  50
  51 //__________________________________________________________________________________________________
  52 Double_t AliRsnFitResult::NormLinear(Double_t *x, Double_t *par)
  53 {
  54 //
  55 // Computes a 1st order polynomial function normalized
  56 // to the integral in the background range specified
  57 // for this object.
  58 // This is obtained by dividing the function value by
  59 // its integral analytically computed in that range.
  60 // Returns 0 in all cases where a floating point exception
  61 // could be raised by computing square root of a negative number.
  62 //
  63
  64   Double_t fcnVal = 1.0 + x[0] * par[1];
  65   Double_t fcnInt = (fPeakRange[1] - fPeakRange[0]) + 0.5*par[1]*(fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0]);
  66
  67   return par[0] * fcnVal / fcnInt;
  68 }
  69
  70 //__________________________________________________________________________________________________
  71 Double_t AliRsnFitResult::NormPoly2(Double_t *x, Double_t *par)
  72 {
  73 //
  74 // Computes a 2nd order polynomial function normalized
  75 // to the integral in the background range specified
  76 // for this object.
  77 // This is obtained by dividing the function value by
  78 // its integral analytically computed in that range.
  79 // Returns 0 in all cases where a floating point exception
  80 // could be raised by computing square root of a negative number.
  81 //
  82
  83   Double_t fcnVal = 1.0 + x[0] * par[1] + x[0]*x[0] * par[2];
  84   Double_t fcnInt = (fPeakRange[1] - fPeakRange[0]) + par[1]*(fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0])/2.0 + par[2]*(fPeakRange[1]*fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0]*fPeakRange[0])/3.0;
  85
  86   return par[0] * fcnVal / fcnInt;
  87 }
  88
  89 //__________________________________________________________________________________________________
  90 Double_t AliRsnFitResult::NormBreitWigner(Double_t *x, Double_t *par)
  91 {
  92 //
  93 // Computes a Breit-Wigner function normalized
  94 // to the integral in the background range specified
  95 // for this object.
  96 // This is obtained by dividing the function value by
  97 // its integral analytically computed in that range.
  98 // Returns 0 in all cases where a floating point exception
  99 // could be raised by computing square root of a negative number.
 100 //
 101
 102   Double_t fcnVal = 1.0 / ((x[0] - par[1])*(x[0] - par[1]) + 0.25*par[2]*par[2]);
 103   Double_t fcnInt = 2.0 / par[2] * (TMath::ATan((fPeakRange[1] - par[1]) / (0.5 * par[2])) - TMath::ATan((fPeakRange[0] - par[1]) / (0.5 * par[2])));
 104
 105   return par[0] * fcnVal / fcnInt;
 106 }
 107
 108 //__________________________________________________________________________________________________
 109 Double_t AliRsnFitResult::NormGaus(Double_t *x, Double_t *par)
 110 {
 111 //
 112 // Computes a Gaussian function normalized
 113 // to the integral in the background range specified
 114 // for this object.
 115 // This is obtained by dividing the function value by
 116 // its integral analytically computed in that range.
 117 // Returns 0 in all cases where a floating point exception
 118 // could be raised by computing square root of a negative number.
 119 //
 120
 121   Double_t fcnVal      = TMath::Gaus(x[0], par[1], par[2], kTRUE);
 122   Double_t fcnIntLeft  = 0.5 * TMath::Erf(((par[1] - fPeakRange[0]) / par[2]) / TMath::Sqrt(2.0));
 123   Double_t fcnIntRight = 0.5 * TMath::Erf(((fPeakRange[1] - par[1]) / par[2]) / TMath::Sqrt(2.0));
 124   Double_t fcnInt      = fcnIntLeft + fcnIntRight;
 125
 126   return par[0] * fcnVal / fcnInt;
 127 }
 128
 129 //__________________________________________________________________________________________________
 130 Double_t AliRsnFitResult::Signal(Double_t *x, Double_t *par)
 131 {
 132 //
 133 // Computes the signal function according to chosen option
 134 //
 135
 136   switch (fSignalType)
 137   {
 138     case kBreitWigner:
 139       return NormBreitWigner(x, par);
 140     case kGaus:
 141       return NormGaus(x, par);
 142     default:
 143       AliError("Invalid signal function");
 144       return 0.0;
 145   }
 146 }
 147
 148 //__________________________________________________________________________________________________
 149 Double_t AliRsnFitResult::Background(Double_t *x, Double_t *par)
 150 {
 151 //
 152 // Computes the background function according to chosen option
 153 //
 154
 155   switch (fBackgroundType)
 156   {
 157     case kSquareRoot:
 158       return NormSquareRoot(x, par);
 159     case kLinear:
 160       return NormLinear(x, par);
 161     case kPoly2:
 162       return NormPoly2(x, par);
 163     default:
 164       AliError("Invalid background function");
 165       return 0.0;
 166   }
 167 }
 168
 169 //__________________________________________________________________________________________________
 170 Double_t AliRsnFitResult::Sum(Double_t *x, Double_t *par)
 171 {
 172 //
 173 // Computes the sum of the signal and the background chosen.
 174 // First 3 parameters are for the signal (they are always 3),
 175 // and other parameters, up to the necessary amount, are for
 176 // the background.
 177 //
 178
 179   Double_t signal     = Signal(x, par);
 180   Double_t background = Background(x, &par[3]);
 181
 182   return signal + background;
 183 }
 184
 185 //__________________________________________________________________________________________________
 186 Int_t AliRsnFitResult::GetNParBackground()
 187 {
 188 //
 189 // Tells how many parameters has the background function
 190 //
 191
 192   switch (fBackgroundType)
 193   {
 194     case kSquareRoot: return 2;
 195     case kLinear    : return 2;
 196     case kPoly2     : return 3;
 197     default         : return 0;
 198   }
 199 }
 200
 201 //__________________________________________________________________________________________________
 202 Bool_t AliRsnFitResult::InitFunctions()
 203 {
 204 //
 205 // Initialize functions
 206 //
 207
 208   //if (!fSignalFcn) fSignalFcn = new TF1(Form("sg%s", GetName()), Signal, fFullRange[0], fFullRange[1], 3);
 209   //if (!fBackgroundFcn) fBackgroundFcn = new TF1(Form("bg%s", GetName()), Background, fFullRange[0], fFullRange[1], GetNParBackground());
 210   //if (!fSumFcn) fSumFcn = new TF1(Form("sum%s", GetName()), Sum, fFullRange[0], fFullRange[1], 3+GetNParBackground());
 211 }
 212
 213 //__________________________________________________________________________________________________
 214 void AliRsnFitResult::SetHistogram(TH1F *histogram)
 215 {
 216 //
 217 // Clones the passed histogram to proceed with fits
 218 //
 219
 220   if (fHistogram) delete fHistogram;
 221
 222   fHistogram = (TH1D*)histogram->Clone();
 223   fHistogram->SetName(Form("h_%s_%s", GetName(), histogram->GetName()));
 224   fHistogram->SetTitle(histogram->GetTitle());
 225 }
 226
 227 //__________________________________________________________________________________________________
 228 Bool_t AliRsnFitResult::SingleFit(const char *opt, Double_t mass, Double_t width)
 229 {
 230 //
 231 // Executes a single fit of the function
 232 //
 233 /*
 234   // require histogram
 235   if (!fHistogram)
 236   {
 237     AliError("Require an initialized histogram!");
 238     return kFALSE;
 239   }
 240
 241   // if necessary, initialize functions
 242   if (!fSignalFcn || !fBackgroundFcn || !fSumFcn) InitFunctions();
 243
 244   // step #0: know how many parameter we have
 245   Int_t npar = GetNParBackground();
 246
 247   // step #1:
 248   // fit outside peak to roughly initialize the background
 249   for (Int_t i = 0; i < npar; i++) fBackgroundFcn->SetParameter(i, 0.5);
 250   status = (Int_t)h->Fit(fcnOut, opt);
 251   if (status) return kFALSE;
 252
 253   // step #2:
 254   // estimate signal/background in the peak
 255   Int_t    imax = h->GetMaximumBin();
 256   Double_t xmax = h->GetBinCenter(imax);
 257   Double_t ymax = h->GetMaximum() - fBackgroundFcn->Eval(xmax);
 258
 259   // step #3:
 260   // fit whole range with signal + background
 261   // which is initialized to the reference values for mass and width
 262   fSignalFcn->SetParameter(0, ymax);
 263   fSignalFcn->SetParameter(1, mass);
 264   fSignalFcn->SetParameter(2, width);
 265   for (Int_t i = 0; i < npar; i++) fcnSum->SetParameter(i + 3, fBackgroundFcn->GetParameter(i));
 266   status = (Int_t)h->Fit(fcnSum, opt);
 267   if (status) return kFALSE;
 268
 269   // get integral and its error here
 270   fValue[kSumIntegralError] = fSumFcn->IntegralError(fitPeakRangeMin, fitPeakRangeMax);
 271   fValue[kSumIntegral]      = fSumFcn->Integral     (fitPeakRangeMin, fitPeakRangeMax);
 272   fValue[kChi2]             = fSumFcn->GetChisquare();
 273   fValue[kNDF]              = fSumFcn->GetNDF();
 274 */
 275   return kTRUE;
 276 }