Changeset 19875
- Timestamp:
- 11/18/19 11:35:07 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/hawc/fit_singles.C
r19861 r19875 23 23 24 24 // -------------------------------------------------------------------------- 25 // after fit_spectrum_bt2b.C 25 26 26 27 // Fit function for a single pe spectrum … … 35 36 const Double_t expo = par[6]; 36 37 38 const Double_t P = cross*TMath::Exp(-cross); 39 37 40 Double_t y = 0; 38 41 for (int N=1; N<14; N++) … … 41 44 const Double_t sigN = TMath::Sqrt(N*sigma*sigma + noise*noise); 42 45 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); 44 47 45 48 y += TMath::Gaus(xx[0], muN, sigN) * p / sigN; … … 56 59 Double_t expo = f.GetParameter(6); 57 60 61 const Double_t P = cross*TMath::Exp(-cross); 62 58 63 Double_t y = 0; 59 64 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); 61 66 62 67 return y / (y + 1); … … 71 76 const Double_t expo = func.GetParameter(6); 72 77 78 const Double_t P = cross*TMath::Exp(-cross); 79 73 80 Double_t sum = 0; 74 81 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); 76 83 77 84 const Double_t scale = hist.GetBinWidth(1); … … 123 130 // (Currently the integral is switched off for the 1440 individual 124 131 // spectra in all cases) 125 bool fast = true; // Switch off using integral132 bool fast = false; // Switch off using integral 126 133 127 134 // Values which should be read from the file but not available atm … … 357 364 Double_t sigma_est = fwhmSum/2.3548; // FWHM = 2*sqrt(2*ln(2))*sigma 358 365 359 Double_t fitmin = maxpos- 3*sigma_est; // was 3366 Double_t fitmin = maxpos-sigma_est; // was 3*sigma_est 360 367 Double_t fitmax = hSum.GetXaxis()->GetXmax(); 361 368 … … 369 376 if (!fixednoise) 370 377 func.SetParLimits(5, 0, 150); // Noise 378 func.SetParLimits(6, 0, 2); // Expo 371 379 372 380 func.SetParameter(0, ampl); // Amplitude 373 func.SetParameter(1, maxpos *1.04);// Gain381 func.SetParameter(1, maxpos); // Gain 374 382 func.SetParameter(2, 0.10); // Sigma 375 383 func.SetParameter(3, 0.25); // Crosstalk … … 378 386 func.FixParameter(5, -1); // Noise 379 387 else 380 func.SetParameter(5, 11.5);// Noise381 382 func.SetParameter(6, 0. 6);// Expo383 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); 385 393 hSum.Fit(&func, fast?"N0QSR":"IN0QSR"); 386 394
Note:
See TracChangeset
for help on using the changeset viewer.