PWGLF/FORWARD/analysis2/corrs/Reference.C

   1 #include <TH1.h>
   2 #include <TF1.h>
   3 #include <TMath.h>
   4 #include <TGraph.h>
   5 #include <TStyle.h>
   6 #include <TLatex.h>
   7
   8 struct Reference
   9 {
  10   static Double_t function(Double_t* x, Double_t* par)
  11   {
  12     const Double_t invSq2Pi = 1. / TMath::Sqrt(TMath::TwoPi());
  13     const Double_t mpShift  = -0.22278298;
  14     const Int_t    np       = 100;
  15     const Double_t sc       = 5;
  16
  17     Double_t       sum      = 0;
  18     Double_t       x0       = x[0];
  19     Double_t       xi       = par[0];
  20     Double_t       mpc      = par[1] - mpShift * xi;
  21     Double_t       area     = par[2];
  22     Double_t       sigma    = par[3];
  23     Double_t       xLow     = x0 - sc * sigma;
  24     Double_t       xUpp     = x0 + sc * sigma;
  25     Double_t       step     = (xUpp - xLow) / np;
  26
  27     for (Int_t i = 1; i <= np/2; i++) {
  28       Double_t x1 = xLow + (i - .5) * step;
  29       Double_t x2 = xUpp - (i - .5) * step;
  30
  31       sum += TMath::Landau(x1, mpc, xi, true) * TMath::Gaus(x0,x1,sigma);
  32       sum += TMath::Landau(x2, mpc, xi, true) * TMath::Gaus(x0,x2,sigma);
  33     }
  34
  35     return area * step * sum * invSq2Pi / sigma;
  36   }
  37
  38   static TF1* fit(TH1*      hist,
  39                   Double_t* range,
  40                   Double_t* guess,
  41                   Double_t* lowLimits,
  42                   Double_t* uppLimits,
  43                   Double_t* params,
  44                   Double_t* errors,
  45                   Double_t& chi2,
  46                   Int_t&    ndf)
  47   {
  48     TF1* func = new TF1(Form("func%s", hist->GetName()), &function,
  49                         range[0], range[1], 4);
  50     func->SetParameters(guess);
  51     func->SetParNames("#xi", "#Delta_{p}", "C", "#sigma");
  52     for (Int_t i = 0; i < 4; i++)
  53       func->SetParLimits(i, lowLimits[i], uppLimits[i]);
  54
  55     hist->Fit(func, "RB0");
  56
  57     for (Int_t i = 0; i < 4; i++) {
  58       params[i] = func->GetParameter(i);
  59       errors[i] = func->GetParError(i);
  60     }
  61     chi2 = func->GetChisquare();
  62     ndf  = func->GetNDF();
  63
  64     return func;
  65   }
  66   static Bool_t search(TF1*   f,    Double_t start, Double_t step,
  67                        Bool_t peak, Double_t& res)
  68   {
  69     const Int_t maxCalls = 10000;
  70     Double_t    lOld     = -2;
  71     Double_t    l        = peak ? -1 : -1e300;
  72     Int_t       i        = 0;
  73     Double_t    x        = 0;
  74     Double_t    cur      = start;
  75
  76     while ((i < maxCalls) &&
  77            TMath::Abs(l - lOld) > 1e-6 &&
  78            TMath::Abs(step) > 1e-8) {
  79       i++;
  80       lOld = l;
  81       x    = cur + step;
  82       l    = f->Eval(x);
  83       if (!peak) l = TMath::Abs(l-res);
  84       // Printf("mode=%d x=%g step=%g l=%g lOld=%g", mode, x, step, l, lOld);
  85
  86       if ((peak && l < lOld) || (!peak && l > lOld))
  87         step = -step/10; // Go the other way in smaller steps
  88
  89       cur += step;
  90     }
  91     if (i >= maxCalls) return false;
  92     res = x;
  93     return true;
  94   }
  95   static Int_t find(TF1* f, Int_t iXi, Int_t iMpc,
  96                     Double_t& maxX, Double_t& fwhm)
  97   {
  98
  99     Double_t    xi       = f->GetParameter(iXi); // params[0];
 100     Double_t    mpc      = f->GetParameter(iMpc); // params[1];
 101     if (!search(f, mpc-0.1*xi,  0.05 * xi, true, maxX))  return false;
 102     Double_t    left     = f->Eval(maxX) / 2;
 103     Double_t    right    = left;
 104     if (!search(f, maxX-.5*xi, xi,         false, left))  return false;
 105     if (!search(f, maxX+xi,    -xi,        false, right)) return false;
 106     fwhm = right - left;
 107
 108     return true;
 109   }
 110   static TH1* histo()
 111   {
 112     Int_t data[100] = {0,    0,  0,  0,  0,  0,  2,  6, 11, 18,
 113                        18,  55, 90,141,255,323,454,563,681,737,
 114                        821,796,832,720,637,558,519,460,357,291,
 115                        279,241,212,153,164,139,106, 95, 91, 76,
 116                        80,  80, 59, 58, 51, 30, 49, 23, 35, 28,
 117                        23,  22, 27, 27, 24, 20, 16, 17, 14, 20,
 118                        12,  12, 13, 10, 17,  7,  6, 12,  6, 12,
 119                        4,    9,  9, 10,  3,  4,  5,  2,  4,  1,
 120                        5,    5,  1,  7,  1,  6,  3,  3,  3,  4,
 121                        5,    4,  4,  2,  2,  7,  2,  4};
 122     TH1D *ret = new TH1D("snr","Signal-to-noise",100,0,100);
 123     for (Int_t i = 0; i < 100; i++) ret->Fill(i, data[i]);
 124
 125     ret->SetFillColor(kRed+1);
 126     ret->SetFillStyle(3001);
 127     ret->SetXTitle("#Delta");
 128     ret->SetDirectory(0);
 129
 130     return ret;
 131   }
 132   static void maxFwhm(TF1* f, Int_t iXi, Int_t iMpc)
 133   {
 134     printf("Find max and FWHM ...");
 135     Double_t maxX = 0;
 136     Double_t fwhm = 0;
 137     if (!find(f, iXi, iMpc, maxX, fwhm))
 138       Printf(" failed");
 139     else {
 140       Printf(" Max: %f FWHM: %f", maxX, fwhm);
 141       Double_t y = f->Eval(maxX);
 142       TGraph* maxG = new TGraph(1);
 143       maxG->SetPoint(0, maxX, y);
 144       maxG->SetMarkerStyle(20);
 145       maxG->SetMarkerColor(kBlue+2);
 146       maxG->Draw("p");
 147
 148       TLatex* maxT = new TLatex(1.2*maxX, y, Form("Max: %f", maxX));
 149       maxT->SetTextAlign(11);
 150       maxT->Draw();
 151
 152       TGraph* fwhmG = new TGraph(2);
 153       fwhmG->SetPoint(0, maxX-fwhm/2, y/2);
 154       fwhmG->SetPoint(1, maxX+fwhm/2, y/2);
 155       fwhmG->SetMarkerStyle(21);
 156       fwhmG->SetMarkerColor(kMagenta+2);
 157       fwhmG->Draw("pl");
 158
 159       TLatex* fwhmT = new TLatex(1.2*(maxX+fwhm/2), y/2,
 160                                  Form("FWHM: %f", fwhm));
 161       fwhmT->SetTextAlign(11);
 162       fwhmT->Draw();
 163     }
 164
 165   }
 166   static void test()
 167   {
 168     TH1* hist = histo();
 169
 170     Double_t range[] = { .3 * hist->GetMean(), 3 * hist->GetMean() };
 171     Double_t guess[] = { 1.8, 20., 5e4, 3.0 };
 172     Double_t low[]   = { 0.5,  5.,  1., 0.4 };
 173     Double_t upp[]   = { 5.0, 50., 1e6, 5.0 };
 174     Double_t chi2    = 0;
 175     Int_t    ndf     = 0;
 176     Double_t params[4];
 177     Double_t errors[4];
 178
 179     printf("Fitting ...");
 180     TF1* f = fit(hist, range, guess, low, upp, params, errors, chi2, ndf);
 181     Printf(" done");
 182     f->Print();
 183
 184     printf("Drawing ...");
 185     gStyle->SetOptStat(1111);
 186     gStyle->SetOptFit(111);
 187     hist->Draw();
 188     f->Draw("same");
 189     Printf(" done");
 190
 191     maxFwhm(f, 0, 1);
 192   }
 193 };