| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "na64util/numerical/langaus.hh" | ||
| 2 | |||
| 3 | #if defined(ROOT_FOUND) && ROOT_FOUND | ||
| 4 | |||
| 5 | #include <TROOT.h> | ||
| 6 | #include <TFitResult.h> | ||
| 7 | |||
| 8 | namespace na64dp { | ||
| 9 | namespace util { | ||
| 10 | |||
| 11 | const Double_t gLandauFMax = -0.22278298; | ||
| 12 | |||
| 13 | Double_t | ||
| 14 | ✗ | langaus( Double_t x, Double_t width | |
| 15 | , Double_t mp | ||
| 16 | , Double_t area | ||
| 17 | , Double_t sigma | ||
| 18 | , Int_t np | ||
| 19 | , Double_t sc | ||
| 20 | ) { | ||
| 21 | // Numeric constants | ||
| 22 | ✗ | const Double_t invsq2pi = 1/sqrt(2*M_PI); | |
| 23 | // Variables | ||
| 24 | Double_t xx; | ||
| 25 | Double_t mpc; | ||
| 26 | Double_t fland; | ||
| 27 | ✗ | Double_t sum = 0.0; | |
| 28 | Double_t xlow,xupp; | ||
| 29 | Double_t step; | ||
| 30 | Double_t i; | ||
| 31 | // MP shift correction | ||
| 32 | ✗ | mpc = mp - gLandauFMax * width; | |
| 33 | // Range of convolution integral | ||
| 34 | ✗ | xlow = x - sc * sigma; | |
| 35 | ✗ | xupp = x + sc * sigma; | |
| 36 | ✗ | step = (xupp-xlow) / np; | |
| 37 | // Convolution integral of Landau and Gaussian by sum | ||
| 38 | ✗ | for(i=1.0; i<=np/2.; i++) { | |
| 39 | ✗ | xx = xlow + (i-.5) * step; | |
| 40 | ✗ | fland = TMath::Landau(xx, mpc, width) / width; | |
| 41 | ✗ | sum += fland * TMath::Gaus(x, xx, sigma); | |
| 42 | ✗ | xx = xupp - (i-.5) * step; | |
| 43 | ✗ | fland = TMath::Landau(xx, mpc, width) / width; | |
| 44 | ✗ | sum += fland * TMath::Gaus(x, xx, sigma); | |
| 45 | } | ||
| 46 | ✗ | return (area * step * sum * invsq2pi / sigma); | |
| 47 | } | ||
| 48 | |||
| 49 | Double_t | ||
| 50 | ✗ | langausfun(Double_t *x, Double_t *par) { | |
| 51 | ✗ | return langaus( *x, par[0], par[1], par[2], par[3] ); | |
| 52 | } | ||
| 53 | |||
| 54 | std::pair<TF1 *, bool> | ||
| 55 | ✗ | langausfit( TH1F * his | |
| 56 | , Double_t * fitrange | ||
| 57 | , Double_t * startvalues | ||
| 58 | , Double_t * parlimitslo | ||
| 59 | , Double_t * parlimitshi | ||
| 60 | , const char * fitOpts | ||
| 61 | , Double_t * fitparams | ||
| 62 | , Double_t * fiterrors | ||
| 63 | , Double_t & chiSqr | ||
| 64 | , Int_t & ndf ) { | ||
| 65 | Char_t FunName[128]; | ||
| 66 | ✗ | snprintf( FunName, sizeof(FunName), "Fitfcn_%s", his->GetName() ); | |
| 67 | |||
| 68 | TF1 *ffitold = (TF1*) ROOT::GetROOT() | ||
| 69 | ✗ | ->GetListOfFunctions() | |
| 70 | ✗ | ->FindObject(FunName); | |
| 71 | ✗ | if( ffitold ) | |
| 72 | ✗ | delete ffitold; | |
| 73 | ✗ | TF1 *ffit = new TF1(FunName, langausfun, fitrange[0], fitrange[1], 4); | |
| 74 | |||
| 75 | ✗ | ffit->SetParameters(startvalues); | |
| 76 | ✗ | ffit->SetParNames( "Width", "MP", "Area", "GSigma" ); | |
| 77 | ✗ | for(Int_t i=0; i<4; i++) { | |
| 78 | ✗ | ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); | |
| 79 | } | ||
| 80 | |||
| 81 | // fit within specified range, use ParLimits, do not plot | ||
| 82 | ✗ | auto result = his->Fit(FunName, fitOpts ? fitOpts : "SLRQ0"); | |
| 83 | |||
| 84 | // Obtain output | ||
| 85 | ✗ | ffit->GetParameters(fitparams); | |
| 86 | ✗ | for(Int_t i=0; i<4; i++) | |
| 87 | ✗ | fiterrors[i] = ffit->GetParError(i); | |
| 88 | ✗ | chiSqr = ffit->GetChisquare(); | |
| 89 | ✗ | ndf = ffit->GetNDF(); | |
| 90 | ✗ | return {ffit, result->IsValid()}; | |
| 91 | } | ||
| 92 | |||
| 93 | } | ||
| 94 | } | ||
| 95 | |||
| 96 | #endif | ||
| 97 | |||
| 98 |