Changeset 14705

Timestamp:

11/27/12 10:33:18 (12 years ago)

Author:

tbretz

Message:

Added more command-line options; added a plot for the estimated relative energy threshold; apply zd and current limits on either plot; small improvements to the plot layout

File:

• trunk/FACT++/src/makeplots.cc (modified) (8 diffs)

trunk/FACT++/src/makeplots.cc

@@ -54 +54 @@
         ("date-time", var<string>(), "SQL time (UTC)")
         ("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
+        ("max-current", var<double>(75), "Maximum current to display in other plots.")
+        ("max-zd", var<double>(75), "Maximum zenith distance to display in other plots")
+        ("no-limits", po_switch(), "Switch off limits in plots")
         ;
 
@@ -102 +105 @@
         time.SetFromStr(conf.Get<string>("date-time"));
 
+    const double max_current = conf.Get<double>("max-current");
+    const double max_zd      = conf.Get<double>("max-zd");
+    const double no_limits   = conf.Get<bool>("no-limits");
+
     ln_rst_time sun_astronomical;
     ln_get_solar_rst_horizon(time.JD(), &observer, -12, &sun_astronomical);
@@ -140 +147 @@
     hframe.SetStats(kFALSE);
     hframe.GetXaxis()->SetTimeFormat("%Hh%M'");
-    hframe.GetXaxis()->SetTitle("Time");
+    hframe.GetXaxis()->SetTitle("Time [UTC]");
     hframe.GetXaxis()->CenterTitle();
     hframe.GetYaxis()->CenterTitle();
     hframe.GetXaxis()->SetTimeDisplay(true);
-    hframe.GetXaxis()->SetLabelSize(0.025);
+    hframe.GetYaxis()->SetTitleSize(0.040);
+    hframe.GetXaxis()->SetTitleSize(0.040);
+    hframe.GetXaxis()->SetTitleOffset(1.1);
+    hframe.GetYaxis()->SetLabelSize(0.040);
+    hframe.GetXaxis()->SetLabelSize(0.040);
 
     TCanvas c1;
+    gPad->SetLeftMargin(0.085);
     gPad->SetRightMargin(0.01);
     gPad->SetTopMargin(0.03);
@@ -156 +168 @@
 
     TCanvas c2;
+    gPad->SetLeftMargin(0.085);
     gPad->SetRightMargin(0.01);
     gPad->SetTopMargin(0.03);
@@ -164 +177 @@
     hframe.DrawCopy();
 
+    TCanvas c3;
+    gPad->SetLeftMargin(0.085);
+    gPad->SetRightMargin(0.01);
+    gPad->SetTopMargin(0.03);
+    gPad->SetGrid();
+    gPad->SetLogy();
+    hframe.GetYaxis()->SetTitle("Estimated relative threshold");
+    hframe.SetMinimum(0.9);
+    hframe.SetMaximum(180);
+    hframe.DrawCopy();
+
     Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
     Int_t style[] = { kSolid, kDashed, kDotted };
@@ -182 +206 @@
 
         // Create graphs
-        TGraph g1, g2;
+        TGraph g1, g2, g3;
         g1.SetName(name.data());
         g2.SetName(name.data());
+        g3.SetName(name.data());
         g1.SetLineWidth(2);
         g2.SetLineWidth(2);
+        g3.SetLineWidth(2);
         g1.SetLineStyle(style[cnt/6]);
         g2.SetLineStyle(style[cnt/6]);
+        g3.SetLineStyle(style[cnt/6]);
         g1.SetLineColor(color[cnt%6]);
         g2.SetLineColor(color[cnt%6]);
+        g3.SetLineColor(color[cnt%6]);
 
         // Loop over 24 hours
@@ -204 +232 @@
             ln_get_hrz_from_equ(&pos, &observer, jd+h, &hrz);
 
-            // Check if source is visible
-            if (hrz.alt<15)
-                continue;
-
-            // Add point to curve
-            const double axis = (mjd+h)*24*3600;
-            g1.SetPoint(g1.GetN(), axis, hrz.alt);
-
-            // Get moon properties and estimate current
+            // Get moon properties and
             ln_equ_posn moon  = fMoonCoords[i].first;
             const double disk = fMoonCoords[i].second;
 
-            ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrz);
-
+            ln_hrz_posn hrzm;
+            ln_get_hrz_from_equ(&moon, &observer, jd+h, &hrzm);
+
+            // Distance between source and moon
             const double angle = Angle(moon.ra, moon.dec, pos.ra, pos.dec);
 
-            const double lc = angle*hrz.alt*pow(disk, 6)/360/360;
-
-            // Add point to curve
-            g2.SetPoint(g2.GetN(), axis, lc>0 ? 7.7+4942*lc : 7.7);
+            // Current prediction
+            const double lc = angle*hrzm.alt*pow(disk, 6)/360/360;
+            const double cur = lc>0 ? 7.7+4942*lc : 7.7;
+
+            // Relative  energy threshold prediction
+            const double cs = cos((90+hrz.alt)*M_PI/180);
+            const double ratio = (10.*sqrt(409600.*cs*cs+9009.) + 6400.*cs - 60.)/10.;
+
+            // Add points to curve
+            const double axis = (mjd+h)*24*3600;
+
+            // If there is a gap of more than one bin, start a new curve
+            if (g1.GetN()>0 && axis-g1.GetX()[g1.GetN()-1]>450)
+            {
+                c1.cd();
+                ((TGraph*)g1.DrawClone("C"))->SetBit(kCanDelete);
+                while (g1.GetN())
+                    g1.RemovePoint(0);
+            }
+
+            if (g2.GetN()>0 && axis-g2.GetX()[g2.GetN()-1]>450)
+            {
+                c2.cd();
+                ((TGraph*)g2.DrawClone("C"))->SetBit(kCanDelete);
+                while (g2.GetN())
+                    g2.RemovePoint(0);
+            }
+
+            if (g3.GetN()>0 && axis-g3.GetX()[g3.GetN()-1]>450)
+            {
+                c3.cd();
+                ((TGraph*)g3.DrawClone("C"))->SetBit(kCanDelete);
+                while (g3.GetN())
+                    g3.RemovePoint(0);
+            }
+
+            // Add data
+            if (no_limits || cur<max_current)
+                g1.SetPoint(g1.GetN(), axis, hrz.alt);
+
+            if (no_limits || 90-hrz.alt<max_zd)
+                g2.SetPoint(g2.GetN(), axis, cur);
+
+            if (no_limits || (cur<max_current && 90-hrz.alt<max_zd))
+                g3.SetPoint(g3.GetN(), axis, ratio*cur/7.7);
         }
 
@@ -237 +300 @@
         g->SetBit(kCanDelete);
 
+        c3.cd();
+        g = (TGraph*)g3.DrawClone("C");
+        g->SetBit(kCanDelete);
+
         leg.AddEntry(g, name.data(), "l");
     }
 
     // Save three plots
-    TCanvas c3;
+    TCanvas c4;
     leg.Draw();
 
     c1.SaveAs("test1.eps");
     c2.SaveAs("test2.eps");
-    c3.SaveAs("legend.eps");
+    c3.SaveAs("test3.eps");
+    c4.SaveAs("legend.eps");
 
     c1.SaveAs("test1.root");
     c2.SaveAs("test2.root");
+    c3.SaveAs("test3.root");
 
     return 0;