PWG2/EBYE/LRC/AliLRCFit.cxx

   1 /**************************************************************************
   2  * Author: Panos Christakoglou.                                           *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //-------------------------------------------------------------------------
  15 //    Description:
  16 //    This class include into LRC library for Long-Range Correlation analysis
  17 //    it makes fit of the 1d histogramm
  18 //    calculates ax+b coefficients with error and hi square
  19 //    Origin: Petr Naumenko, SPbSU-CERN, Petr.Naoumenko@cern.ch
  20 //-------------------------------------------------------------------------
  21
  22 /* $Id$ */
  23
  24 #include "AliLRCFit.h"
  25
  26
  27 ClassImp(AliLRCFit)
  28
  29
  30 AliLRCFit::AliLRCFit():fN (0), fTrueN(0), fNmin (0), fNmax(0), fS1(.0), fSz(.0), fSfz(.0), fSf(.0), fSf2(.0), fhi2(.0), fw(.0), fz(.0), f(.0), fnum(.0), fdf(.0), fdelta(.0), fa(.0), fb(.0), fda(.0), fdb(.0), fda1(.0), fdb1(.0), fxmin(.0), fxmax(.0){
  31 //Empty constructor
  32 }
  33
  34 AliLRCFit::AliLRCFit(TH1D* h):fN (0), fTrueN(0), fNmin (0), fNmax(0), fS1(.0), fSz(.0), fSfz(.0), fSf(.0), fSf2(.0), fhi2(.0), fw(.0), fz(.0), f(.0), fnum(.0), fdf(.0), fdelta(.0), fa(.0), fb(.0), fda(.0), fdb(.0), fda1(.0), fdb1(.0), fxmin(.0), fxmax(.0){
  35     //Constructor make fit of 1d histogramm
  36     fxmin = h->GetXaxis()->GetXmin();
  37     fxmax = h->GetXaxis()->GetXmax();
  38     fN=h->GetNbinsX();
  39     fNmin=1;
  40     fNmax=fN;
  41     fnum = h->GetXaxis()->GetXmin();
  42     fdf = (h->GetXaxis()->GetXmax()-fnum)/fN;
  43     for(int i=1; i<=fN; i++){
  44         f=i*fdf-fdf/2+fnum;
  45         fw = h->GetBinError(i);
  46         if(fw!=0){
  47             fz = h->GetBinContent(i);
  48             fS1 = fS1 + 1/(fw*fw);
  49             fSz = fSz + fz/(fw*fw);
  50             fSfz = fSfz + f*fz/(fw*fw);
  51             fSf = fSf + f/(fw*fw);
  52             fSf2 = fSf2 + (f*f)/(fw*fw);
  53         }
  54     }
  55     fdelta = fS1*fSf2 - fSf*fSf;
  56     fb = (fS1*fSfz - fSf*fSz)/fdelta;
  57     fa = ((fSf2-fSf)*fSz - (fSf-fS1)*fSfz)/fdelta;
  58     fda = sqrt((fS1+fSf2-2*fSf)/fdelta);
  59     fdb = sqrt(fS1/fdelta);
  60     fdb1 = 0;
  61     fda1 = 0;
  62     f = h->GetXaxis()->GetXmin();
  63     for(int i=1; i<=fN; i++){
  64         f=i*fdf-fdf/2+fnum;
  65         fw = h->GetBinError(i);
  66         if(fw!=0){
  67             fz = h->GetBinContent(i);
  68             fdb1 = fdb1 + ((fS1*f - fSf)*(fS1*f - fSf)/(fw*fw)) * ((fz-fa-fb*(f-1))*(fz-fa-fb*(f-1))/(fw*fw));
  69             fda1 = fda1 + (((fSf2-fSf)-(fSf-fS1)*f)*((fSf2-fSf)-(fSf-fS1)*f)/(fw*fw)) * ((fz-fa+fb-fb*f)*(fz-fa+fb-fb*f)/(fw*fw));
  70             fhi2 = fhi2 + ((fz-(fa-fb)-fb*f) * (fz-(fa-fb)-fb*f)) / (fw*fw);
  71             fTrueN++;
  72         }
  73     }
  74     fdb1 = sqrt(fdb1/(fdelta*fdelta));
  75     fda1 = sqrt(fda1/(fdelta*fdelta));
  76     fhi2 = fhi2 / (fTrueN-2);
  77 }
  78
  79 AliLRCFit::AliLRCFit(TH1D *h, double xmin, double xmax):fN (0), fTrueN(0), fNmin (0), fNmax(0), fS1(.0), fSz(.0), fSfz(.0), fSf(.0), fSf2(.0), fhi2(.0), fw(.0), fz(.0), f(.0), fnum(.0), fdf(.0), fdelta(.0), fa(.0), fb(.0), fda(.0), fdb(.0), fda1(.0), fdb1(.0), fxmin(.0), fxmax(.0){
  80      //Constructor make fit of 1d histogramm between xmin and xmax
  81     fxmin = xmin;
  82     fxmax = xmax;
  83     fN=h->GetNbinsX();
  84     fnum = h->GetXaxis()->GetXmin();
  85     fdf = (h->GetXaxis()->GetXmax()-fnum)/fN;
  86     fNmin=int((xmin-f)/fdf)+1;
  87     fNmax=int((xmax-f)/fdf)+1;
  88     for(int i=fNmin; i<=fNmax; i++){
  89         f=i*fdf-fdf/2+fnum;
  90         fw = h->GetBinError(i);
  91         if(fw!=0){
  92             fz = h->GetBinContent(i);
  93             fS1 = fS1 + 1/(fw*fw);
  94             fSz = fSz + fz/(fw*fw);
  95             fSfz = fSfz + f*fz/(fw*fw);
  96             fSf = fSf + f/(fw*fw);
  97             fSf2 = fSf2 + (f*f)/(fw*fw);
  98         }
  99     }
 100     fdelta = fS1*fSf2 - fSf*fSf;
 101     fb = (fS1*fSfz - fSf*fSz)/fdelta;
 102     fa = ((fSf2-fSf)*fSz - (fSf-fS1)*fSfz)/fdelta;
 103     fda = sqrt((fS1+fSf2-2*fSf)/fdelta);
 104     fdb = sqrt(fS1/fdelta);
 105     fdb1 = 0;
 106     fda1 = 0;
 107     for(int i=fNmin; i<=fNmax; i++){
 108         f=i*fdf-fdf/2+fnum;
 109         fw = h->GetBinError(i);
 110         if(fw!=0){
 111             fz = h->GetBinContent(i);
 112             fdb1 = fdb1 + ((fS1*f - fSf)*(fS1*f - fSf)/(fw*fw)) * ((fz-fa-fb*(f-1))*(fz-fa-fb*(f-1))/(fw*fw));
 113             fda1 = fda1 + (((fSf2-fSf)-(fSf-fS1)*f)*((fSf2-fSf)-(fSf-fS1)*f)/(fw*fw)) * ((fz-fa+fb-fb*f)*(fz-fa+fb-fb*f)/(fw*fw));
 114             fhi2 = fhi2 + ((fz-(fa-fb)-fb*f) * (fz-(fa-fb)-fb*f)) / (fw*fw);
 115             fTrueN++;
 116         }
 117     }
 118     fdb1 = sqrt(fdb1/(fdelta*fdelta));
 119     fda1 = sqrt(fda1/(fdelta*fdelta));
 120     fhi2 = fhi2 / (fTrueN-2);
 121 }
 122
 123 AliLRCFit::~AliLRCFit() {
 124 }
 125
 126 double AliLRCFit::Geta() const {return fa;}
 127 double AliLRCFit::Getb() const {return fb;}
 128 double AliLRCFit::Getda() const {return fda;}
 129 double AliLRCFit::Getdb() const {return fdb;}
 130 double AliLRCFit::Getda1() const {return fda1;}
 131 double AliLRCFit::Getdb1() const {return fdb1;}
 132 double AliLRCFit::Gethi2() const {return fhi2;}
 133 double AliLRCFit::Getf() const {return f;}
 134 double AliLRCFit::Getxmin() const {return fxmin;}
 135 double AliLRCFit::Getxmax() const {return fxmax;}
 136 int AliLRCFit::GetN() const {return fN;}
 137
 138
 139
 140