Index: /trunk/Mars/hawc/fit_singles.C
===================================================================
--- /trunk/Mars/hawc/fit_singles.C	(revision 19874)
+++ /trunk/Mars/hawc/fit_singles.C	(revision 19875)
@@ -23,4 +23,5 @@
 
 // --------------------------------------------------------------------------
+// after fit_spectrum_bt2b.C
 
 // Fit function for a single pe spectrum
@@ -35,4 +36,6 @@
     const Double_t expo  = par[6];
 
+    const Double_t P = cross*TMath::Exp(-cross);
+
     Double_t y = 0;
     for (int N=1; N<14; N++)
@@ -41,5 +44,5 @@
         const Double_t sigN = TMath::Sqrt(N*sigma*sigma + noise*noise);
 
-        const Double_t p = TMath::Power(cross, N-1) * TMath::Power(N, -expo);
+        const Double_t p = TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
 
         y += TMath::Gaus(xx[0], muN, sigN) * p / sigN;
@@ -56,7 +59,9 @@
     Double_t expo  = f.GetParameter(6);
 
+    const Double_t P = cross*TMath::Exp(-cross);
+
     Double_t y = 0;
     for (int N=2; N<14; N++)
-        y +=  TMath::Power(cross,  N-1) * TMath::Power(N, -expo);
+        y += TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
 
     return y / (y + 1);
@@ -71,7 +76,9 @@
     const Double_t expo  = func.GetParameter(6);
 
+    const Double_t P = cross*TMath::Exp(-cross);
+
     Double_t sum = 0;
     for (int N=1; N<14; N++)
-        sum += TMath::Power(cross, N-1) * TMath::Power(N, -expo);
+        sum += TMath::Power(N*P, N-1)/TMath::Power(TMath::Factorial(N-1), expo);
 
     const Double_t scale = hist.GetBinWidth(1);
@@ -123,5 +130,5 @@
     // (Currently the integral is switched off for the 1440 individual
     // spectra in all cases)
-    bool fast = true; // Switch off using integral
+    bool fast = false; // Switch off using integral
 
     // Values which should be read from the file but not available atm
@@ -357,5 +364,5 @@
     Double_t sigma_est = fwhmSum/2.3548; // FWHM = 2*sqrt(2*ln(2))*sigma
 
-    Double_t fitmin = maxpos-3*sigma_est; // was 3
+    Double_t fitmin = maxpos-sigma_est; // was 3*sigma_est
     Double_t fitmax = hSum.GetXaxis()->GetXmax();
 
@@ -369,7 +376,8 @@
     if (!fixednoise)
         func.SetParLimits(5, 0, 150);      // Noise
+    func.SetParLimits(6, 0, 2);            // Expo
 
     func.SetParameter(0, ampl);                         // Amplitude
-    func.SetParameter(1, maxpos*1.04);                  // Gain
+    func.SetParameter(1, maxpos);                       // Gain
     func.SetParameter(2, 0.10);                         // Sigma
     func.SetParameter(3, 0.25);                         // Crosstalk
@@ -378,9 +386,9 @@
         func.FixParameter(5, -1);                       // Noise
     else
-        func.SetParameter(5, 11.5);                     // Noise
-
-    func.SetParameter(6, 0.6);                          // Expo
-
-    func.SetRange(maxpos-20, fitmax);
+        func.SetParameter(5,  0.1*maxpos);              // Noise
+
+    func.SetParameter(6, 0.95);                         // Expo
+
+    func.SetRange(fitmin, fitmax);
     hSum.Fit(&func, fast?"N0QSR":"IN0QSR");