Changeset 19875


Ignore:
Timestamp:
Nov 18, 2019, 11:35:07 AM (10 months ago)
Author:
tbretz
Message:
This was the wrong version of the function. This updates it to be identical with the paper.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/hawc/fit_singles.C

    r19861 r19875  
    2323
    2424// --------------------------------------------------------------------------
     25// after fit_spectrum_bt2b.C
    2526
    2627// Fit function for a single pe spectrum
     
    3536    const Double_t expo  = par[6];
    3637
     38    const Double_t P = cross*TMath::Exp(-cross);
     39
    3740    Double_t y = 0;
    3841    for (int N=1; N<14; N++)
     
    4144        const Double_t sigN = TMath::Sqrt(N*sigma*sigma + noise*noise);
    4245
    43         const Double_t p = TMath::Power(cross, N-1) * TMath::Power(N, -expo);
     46        const Double_t p = TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
    4447
    4548        y += TMath::Gaus(xx[0], muN, sigN) * p / sigN;
     
    5659    Double_t expo  = f.GetParameter(6);
    5760
     61    const Double_t P = cross*TMath::Exp(-cross);
     62
    5863    Double_t y = 0;
    5964    for (int N=2; N<14; N++)
    60         y +=  TMath::Power(cross,  N-1) * TMath::Power(N, -expo);
     65        y += TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
    6166
    6267    return y / (y + 1);
     
    7176    const Double_t expo  = func.GetParameter(6);
    7277
     78    const Double_t P = cross*TMath::Exp(-cross);
     79
    7380    Double_t sum = 0;
    7481    for (int N=1; N<14; N++)
    75         sum += TMath::Power(cross, N-1) * TMath::Power(N, -expo);
     82        sum += TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
    7683
    7784    const Double_t scale = hist.GetBinWidth(1);
     
    123130    // (Currently the integral is switched off for the 1440 individual
    124131    // spectra in all cases)
    125     bool fast = true; // Switch off using integral
     132    bool fast = false; // Switch off using integral
    126133
    127134    // Values which should be read from the file but not available atm
     
    357364    Double_t sigma_est = fwhmSum/2.3548; // FWHM = 2*sqrt(2*ln(2))*sigma
    358365
    359     Double_t fitmin = maxpos-3*sigma_est; // was 3
     366    Double_t fitmin = maxpos-sigma_est; // was 3*sigma_est
    360367    Double_t fitmax = hSum.GetXaxis()->GetXmax();
    361368
     
    369376    if (!fixednoise)
    370377        func.SetParLimits(5, 0, 150);      // Noise
     378    func.SetParLimits(6, 0, 2);            // Expo
    371379
    372380    func.SetParameter(0, ampl);                         // Amplitude
    373     func.SetParameter(1, maxpos*1.04);                  // Gain
     381    func.SetParameter(1, maxpos);                       // Gain
    374382    func.SetParameter(2, 0.10);                         // Sigma
    375383    func.SetParameter(3, 0.25);                         // Crosstalk
     
    378386        func.FixParameter(5, -1);                       // Noise
    379387    else
    380         func.SetParameter(5, 11.5);                     // Noise
    381 
    382     func.SetParameter(6, 0.6);                          // Expo
    383 
    384     func.SetRange(maxpos-20, fitmax);
     388        func.SetParameter(5,  0.1*maxpos);              // Noise
     389
     390    func.SetParameter(6, 0.95);                         // Expo
     391
     392    func.SetRange(fitmin, fitmax);
    385393    hSum.Fit(&func, fast?"N0QSR":"IN0QSR");
    386394
Note: See TracChangeset for help on using the changeset viewer.