Changeset 19904 for trunk/FACT++

Timestamp:

12/15/19 16:23:56 (5 years ago)

Author:

tbretz

Message:

Implemented the Flux Sigma and control plots to check if a point should be replaced by an UL.

Location:

trunk/FACT++

Files:

• spectrum/display.C (modified) (4 diffs)
• spectrum/spectrum.sql (modified) (1 diff)
• src/spectrum.cc (modified) (4 diffs)

trunk/FACT++/spectrum/display.C

@@ -32 +32 @@
     fx.SetLineColor(kBlack);
 
-    TF1 f0("f0", "0", 0, 90);
+    TF1 f0("f0", "[0]", 0, 90);
+    f0.SetParameter(0, 0);
     f0.SetLineStyle(kDashed);
     f0.SetLineColor(kBlack);
@@ -211 +212 @@
 
     // --------------------------------
+    f0.SetLineColor(kRed);
+    // --------------------------------
 
     c = new TCanvas("Integral Spectrum", "Integral Spectrum");
-    c->SetLogy();
+    c->Divide(2,2);
+    c->cd(4);
+    gPad->SetLogy();
     file.GetObject("Data/Energy/Differential/IntegratedSpectrum", h1);
     h1->SetLineColor(kGray);
     h1->SetMarkerColor(kGray);
+    if (h1->GetMinimum()<=0)
+    {
+        h1->SetMinimum(h1->GetMinimum(0)*0.1);
+        cout << "WARNING: Integral Spectrum contains negative flux points!\n" << endl;
+    }
     h1->DrawCopy("P");
     file.GetObject("Data/Energy/Integral/Spectrum", h1);
@@ -238 +248 @@
     //h1->DrawCopy("P same");
 
+    c->cd(1);
+    gPad->SetGridx();
+    file.GetObject("Data/Energy/Integral/SigmaFlux", h1);
+    h1->SetMaximum(5);
+    h1->SetMarkerStyle(kFullDotLarge);
+    h1->DrawCopy("P");
+    f0.SetParameter(0, 1);
+    f0.DrawCopy("same");
+
+    c->cd(3);
+    gPad->SetGridx();
+    file.GetObject("Data/Energy/Integral/BackgroundI", h1);
+    h1->Scale(5);
+    h1->SetMaximum(30);
+    h1->SetLineColor(kRed);
+    h1->DrawCopy("P");
+    file.GetObject("Data/Energy/Integral/SignalI", h1);
+    h1->DrawCopy("P same");
+    file.GetObject("Data/Energy/Integral/ExcessI", h1);
+    h1->SetLineColor(kBlue);
+    h1->DrawCopy("P same");
+
+    f0.SetParameter(0, 10);
+    f0.DrawCopy("same");
+
     // --------------------------------
 
     c = new TCanvas("Differential Spectrum", "Differential Spectrum");
-    c->SetLogy();
+    c->Divide(2,2);
+    c->cd(4);
+    gPad->SetLogy();
     file.GetObject("Data/Energy/Differential/Spectrum", h1);
+    if (h1->GetMinimum()<=0)
+    {
+        h1->SetMinimum(h1->GetMinimum(0)*0.1);
+        cout << "WARNING: Differential Spectrum contains negative flux points!\n" << endl;
+    }
     h1->DrawCopy("P");
 
@@ -261 +303 @@
     //h1->SetMarkerStyle(22);
     //h1->DrawCopy("P same");
+
+    c->cd(1);
+    gPad->SetGridx();
+    file.GetObject("Data/Energy/Differential/SigmaFlux", h1);
+    h1->SetMaximum(5);
+    h1->SetMarkerStyle(kFullDotLarge);
+    h1->DrawCopy("P");
+    f0.SetParameter(0, 1);
+    f0.DrawCopy("same");
+
+    c->cd(3);
+    gPad->SetGridx();
+    file.GetObject("Data/Energy/Differential/Background", h1);
+    h1->Scale(5);
+    h1->SetMaximum(30);
+    h1->SetLineColor(kRed);
+    h1->DrawCopy("P");
+    file.GetObject("Data/Energy/Differential/Signal", h1);
+    h1->DrawCopy("P same");
+    file.GetObject("Data/Energy/Differential/Excess", h1);
+    h1->SetLineColor(kBlue);
+    h1->DrawCopy("P same");
+
+    f0.SetParameter(0, 10);
+    f0.DrawCopy("same");
 }

trunk/FACT++/spectrum/spectrum.sql

@@ -74 +74 @@
       %106:join2
 ),
+Flux AS
+(
+   SELECT -- Return final result
+      *,
+
+      -- Differetial Spectrum
+
+      SimExcess/SimFluxW  AS  Efficiency,
+
+      Excess/SimExcess/Width/AreaTime  AS  ExcessRatio,
+         1/SQRT(
+             + POW(ErrExcess    / Excess,    2)
+             + POW(ErrSimExcess / SimExcess, 2)
+         )  AS  SigmaExcessRatio,
+
+
+      Excess/SimExcess*SimFluxW/Width/AreaTime  AS  Flux,
+         1/SQRT(
+             + POW(ErrExcess    / Excess,    2)
+             + POW(ErrSimExcess / SimExcess, 2)
+             + POW(ErrSimFluxW  / SimFluxW,  2)
+         )  AS  SigmaFlux,
+
+      -- Integral Spectrum
+
+      SimExcessI/SimFluxI  AS  EfficiencyI,
+
+      ExcessI/SimExcessI/AreaTime  AS  ExcessRatioI,
+         1/SQRT(
+             + POW(ErrExcessI    / ExcessI,    2)
+             + POW(ErrSimExcessI / SimExcessI, 2)
+         )  AS  SigmaExcessRatioI,
+
+
+      ExcessI/SimExcessI*SimFluxI/AreaTime  AS  FluxI,
+         1/SQRT(
+             + POW(ErrExcessI    / ExcessI,    2)
+             + POW(ErrSimExcessI / SimExcessI, 2)
+             + POW(ErrSimFluxI   / SimFluxI,   2)
+         )  AS  SigmaFluxI
+
+   FROM
+      CombinedData
+),
 Spectrum AS
 (
-SELECT -- Return final result
-   *,
+   SELECT -- Return final result
+      *,
 
-   -- Differetial Spectrum
-
-   SimExcess/SimFluxW  AS  Efficiency,
-
-   Excess/SimExcess/Width/AreaTime  AS  ExcessRatio,
-   Excess/SimExcess/Width/AreaTime
-      * SQRT(
-          + POW(ErrExcess    / Excess,    2)
-          + POW(ErrSimExcess / SimExcess, 2)
-        )  AS  ErrExcessRatio,
-
-
-   Excess/SimExcess*SimFluxW/Width/AreaTime  AS  Flux,
-   Excess/SimExcess*SimFluxW/Width/AreaTime
-      * SQRT(
-          + POW(ErrExcess    / Excess,    2)
-          + POW(ErrSimExcess / SimExcess, 2)
-          + POW(ErrSimFluxW  / SimFluxW,  2)
-        )  AS  ErrFlux,
-
-   -- Integral Spectrum
-
-   SimExcessI/SimFluxI  AS  EfficiencyI,
-
-   ExcessI/SimExcessI/AreaTime  AS  ExcessRatioI,
-   ExcessI/SimExcessI/AreaTime
-      * SQRT(
-          + POW(ErrExcessI    / ExcessI,    2)
-          + POW(ErrSimExcessI / SimExcessI, 2)
-        )  AS  ErrExcessRatioI,
-
-
-   ExcessI/SimExcessI*SimFluxI/AreaTime  AS  FluxI,
-   ExcessI/SimExcessI*SimFluxI/AreaTime
-      * SQRT(
-          + POW(ErrExcessI    / ExcessI,    2)
-          + POW(ErrSimExcessI / SimExcessI, 2)
-          + POW(ErrSimFluxI   / SimFluxI,   2)
-        )  AS  ErrFluxI
-
-FROM
-   CombinedData
-WHERE
-   Excess>0
+      ExcessRatio /SigmaExcessRatio   AS  ErrExcessRatio,
+      ExcessRatioI/SigmaExcessRatioI  AS  ErrExcessRatioI,
+      ABS(Flux)   /SigmaFlux          AS  ErrFlux,
+      ABS(FluxI)  /SigmaFluxI         AS  ErrFluxI
+   FROM
+      Flux
 )
 SELECT

trunk/FACT++/src/spectrum.cc

@@ -2060 +2060 @@
             map<size_t, double> rolke_int;
 
+            if (verbose>0)
+                cout << "Bin     Excess     Significance          Flux ErrFlux" << endl;
+
             for (auto ir=res13.cbegin(); ir!=res13.cend(); ir++)
             {
@@ -2094 +2097 @@
                 if (verbose>0)
                 {
-                    cout << setw(5)  << row["center"]  << ": ";
-                    cout << setw(10) << row["Excess"]  << " ";
-                    cout << setw(10) << row["Flux"]    << " ";
-                    cout << setw(10) << row["ErrFlux"] << " ";
-                    cout << setw(10) << row["Significance"] << '\n';
+                    cout << setw(5)  << row["center"] << ":";
+                    cout << " " << setw(10) << row["Excess"];
+                    cout << " " << setw(10) << row["Significance"];
+                    cout << " " << setw(10) << row["Flux"];
+                    cout << " " << setw(10) << row["ErrFlux"];
+                    cout << endl;
                 }
             }
@@ -2169 +2173 @@
             WriteHistogram(connection, hist);
 
+            hist.name = "SigmaFlux";
+            hist.v    = "SigmaFlux";
+            hist.err  = "";
+            WriteHistogram(connection, hist);
+
             if (*ib=="Energy")
             {
@@ -2177 +2186 @@
                 hist.v     = "FluxI";
                 hist.err   = "ErrFluxI";
+                WriteHistogram(connection, hist);
+
+                hist.dir   = "Data/Energy/Integral";
+                hist.name  = "SigmaFlux";
+                hist.v     = "SigmaFluxI";
+                hist.err   = "";
                 WriteHistogram(connection, hist);