Changeset 14181 for fact

Timestamp:

06/15/12 14:59:24 (12 years ago)

Author:

Jens Buss

Message:

fitt to sum spectrum added

File:

• fact/tools/marsmacros/singlepe.C (modified) (5 diffs)

fact/tools/marsmacros/singlepe.C

@@ -593 +593 @@
     func.SetLineColor(kBlue);
 
+    // define fit function for Amplitudespectrum
+    TF1 func2("sum_spektrum", fcn, 0, 2000, 5);
+    func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
+    func2.SetLineColor(kBlue);
+
+    // define fit function for Amplitudespectrum
+    TF1 func3("gain_sum_spektrum", fcn, 0, 10, 5);
+    func3.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
+    func3.SetLineColor(kBlue);
+
+
+
     // Map for which pixel shall be plotted and which not
     TArrayC usePixel(1440);
@@ -695 +707 @@
             hist->SetName(Form("Pix%d", pixel));
             hist->DrawCopy()->SetDirectory(0);
-            //func.DrawCopy("same");
-        }
+            func.DrawCopy("SAME");
+        }
+
+
 
         count_ok++;
@@ -702 +716 @@
         delete hist;
     }
+
+    //------------------fill histograms-----------------------
 
     hcam.SetUsed(usePixel);
@@ -710 +726 @@
     hcam.DrawCopy();
 
+    gROOT->SetSelectedPad(0);
     TCanvas &c11 = d->AddTab("SumHist");
     gPad->SetLogy();
     hSum1.DrawCopy();
 
+    //------------------------fit sum spectrum-------------------------------
+    const Int_t    maxbin   = hSum1.GetMaximumBin();
+    const Double_t binwidth = hSum1.GetBinWidth(maxbin);
+    const Double_t ampl     = hSum1.GetBinContent(maxbin);
+
+
+    double fwhm2 = 0;
+    for (int i=1; i<maxbin; i++)
+        if (hSum1.GetBinContent(maxbin-i)+hSum1.GetBinContent(maxbin+i)<ampl*2)
+        {
+            fwhm2 = (2*i+1)*binwidth;
+            break;
+        }
+
+    if (fwhm2==0)
+    {
+        cout << "Could not determine start value for sigma." << endl;
+    }
+
+
+    Double_t par[5] =
+    {
+        ampl, hSum1.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10
+    };
+
+    func2.SetParameters(par);
+
+    const double min = par[1]-fwhm2*1.4;
+    const double max = par[1]*5;
+
+    hSum1.Fit(&func2, "NOQS", "", min, max);
+    func2.DrawCopy("SAME");
+    //-------------------END OF: fit sum spectrum-----------------------------
+
+
     TCanvas &c12 = d->AddTab("GainHist");
     gPad->SetLogy();
     hSum2.DrawCopy();
+    const Double_t fMaxPos  = hSum2.GetBinCenter(maxbin);
+
+    //------------------fit gausses to peaks -----------------------
+
+    UInt_t nfunc = 5;
+    TF1 **gaus = new TF1*[nfunc];
+
+    for (UInt_t nr = 0; nr<nfunc; nr++)
+    {
+      char fname[20];
+      sprintf(fname,"gaus%d",nr);
+      gaus[nr] = new TF1(fname,"gaus",  fMaxPos*(nr+1)-0.6, fMaxPos*(nr+1)+0.6);
+        hSum2.Fit(gaus[nr], "N0QS", "", fMaxPos*(nr+1)-0.6, fMaxPos*(nr+1)+0.6);
+        gaus[nr]->DrawCopy("SAME");
+//        delete gaus[nr];
+    }
+
+
+    //------------------fit gain corrected sum spectrum-----------------------
+
+    const Int_t    maxbin2   = hSum2.GetMaximumBin();
+    const Double_t binwidth2 = hSum2.GetBinWidth(maxbin);
+    const Double_t ampl2     = hSum2.GetBinContent(maxbin);
+    cout << "AMPL: " << ampl2 << endl;
+
+    fwhm2 = 0;
+    for (int i=1; i<maxbin; i++)
+        if (hSum2.GetBinContent(maxbin2-i)+hSum2.GetBinContent(maxbin2+i)<ampl2*2)
+        {
+            fwhm2 = (2*i+1)*binwidth2;
+            break;
+        }
+
+    if (fwhm2==0)
+    {
+        cout << "Could not determine start value for sigma." << endl;
+    }
+
+
+    Double_t par2[5] =
+    {
+        ampl2, hSum2.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, 0
+    };
+
+        cout << "AMPL: " << ampl2 << endl;
+            cout << "maxpos: " << hSum2.GetBinCenter(maxbin2) << endl;
+                cout << "fwhm: " << fwhm2*1.4 << endl;
+
+    func3.SetParameters(par2);
+
+    const double min2 = par2[1]-fwhm2*1.4;
+    const double max2 = par2[1]*5;
+
+    hSum2.Fit(&func3, "N0QS", "", min2, max2);
+    func3.DrawCopy("SAME");
+    //-------------------END OF: fit sum spectrum-----------------------------
 
     TCanvas &c2 = d->AddTab("Result");
@@ -765 +873 @@
 
     Double_t y = 0;
-    for (int order = 1; order < 5; order++)
+    for (int order = 1; order < 14; order++)
     {
         Double_t arg   = (x - order*gain - shift)/sigma;