Index: /trunk/Mars/fact/plots/quality.C
===================================================================
--- /trunk/Mars/fact/plots/quality.C	(revision 17651)
+++ /trunk/Mars/fact/plots/quality.C	(revision 17652)
@@ -1,2 +1,5 @@
+#include <algorithm>
+#include <functional>
+
 Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0)
 {
@@ -87,59 +90,8 @@
     return vecp[vecp.size()-1].second;
 }
-/*
+
 Float_t Prediction(Double_t time)
 {
     Double_t jd = time + 40587 + 2400000.5;
-
-    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
-
-    // Nova::EquPosn pos = FindPointing(time);
-
-    // get local position of moon
-    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
-
-    // Distance between source and moon
-    //const double angle = Nova::GetAngularSeparation(moon, hrzm)*TMath::DegToRad();
-
-    // Current prediction
-    //double cang = sin(angle);
-    double calt = sin(hrzm.alt*TMath::DegToRad());
-
-    double disk = Nova::GetLunarDisk(jd);
-
-    //                                                                           semi-major axis
-    double lc = calt*pow(disk, 2.2)*pow(Nova::GetLunarEarthDist(jd)/384400, -2);///sqrt(cang);
-
-    return lc>0 ? 4+103*lc : 4;
-}*/
-
-Float_t Prediction(Double_t time)
-{
-    Double_t jd = time + 40587 + 2400000.5;
-    /*
-    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
-
-    Nova::EquPosn pos = FindPointing(time);
-
-    // get local position of moon
-    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
-
-    // Distance between source and moon
-    const double angle = Nova::GetAngularSeparation(moon, pos);
-
-    // Distance between earth and moon relative to major semi-axis
-    const double dist  = Nova::GetLunarEarthDist(jd)/384400;
-
-    // Current prediction
-    double cang = 1-sin(angle*TMath::DegToRad());
-    double calt = sin(hrzm.alt*TMath::DegToRad());
-
-    double disk = Nova::GetLunarDisk(jd);
-
-    // light condition
-    double lc = sqrt(calt)*pow(disk, 2.3)*pow(cang+1, 1)*pow(dist, -2);
-
-    return lc>0 ? 7.2 + 69*lc : 7.2;
-    */
 
     // Sun properties
@@ -161,19 +113,16 @@
     // Current prediction
     double sin_malt  = hrzm.alt<0 ? 0 : sin(hrzm.alt*TMath::DegToRad());
-    double sin_mdist = sin(angle*TMath::DegToRad());
-    double cos_salt  = cos(hrzs.zd*TMath::DegToRad());
-
-    double c0 = pow(disk,     2.52);
-    double c1 = pow(sin_malt, 0.72);
-    double c2 = pow(edist,   -2.00);
-    double c3 = exp(1.46*(1-sin_mdist)*(1-sin_mdist));
-    double c4 = exp(33.0)*exp(-20.1*(1-cos_salt)*(1-cos_salt));
-
-    double cur = 6.4 + 96.9*c0*c1*c2*c3 + c4;
-
-   // cout << cur << " " << hrzm.alt << " " << c1 << endl;
+    double cos_mdist = cos(angle*TMath::DegToRad());
+    double sin_szd   = sin(hrzs.zd*TMath::DegToRad());
+
+    double c0 = pow(disk,      2.63);
+    double c1 = pow(sin_malt,  0.60);
+    double c2 = pow(edist,    -2.00);
+    double c3 = exp(0.67*cos_mdist*cos_mdist*cos_mdist*cos_mdist);
+    double c4 = exp(-97.8+105.8*sin_szd*sin_szd);
+
+    double cur = 6.2 + 95.7*c0*c1*c2*c3 + c4;
 
     return cur;
-
 }
 
@@ -228,4 +177,6 @@
     Double_t time;
     Float_t  Imed;
+    Float_t  Idev;
+    Float_t  I[416];
 
     if (!file.SetPtrAddress("Time",  &time))
@@ -233,20 +184,42 @@
 
     if (!file.SetPtrAddress("I_med", &Imed))
+        return -1;
+
+    if (!file.SetPtrAddress("I_dev", &Idev))
+        return -1;
+
+    if (!file.SetPtrAddress("I", I))
         return -1;
 
     TGraph g1;
     TGraph g2;
-    g1.SetName("Currents");
+    TGraph g3;
+    TGraph g4;
+    TGraph g5;
+    g1.SetName("CurrentsMed");
     g2.SetName("Prediction");
+    g3.SetName("CurrentsDev");
+    g4.SetName("CurrentsMax-4");
+    g5.SetName("CurrentsMax");
 
     while (file.GetNextRow())
         if (Contains(vec, time))
         {
+            // crazy pixels
+            I[66]  = 0;
+            I[191] = 0;
+            I[193] = 0;
+
+            sort(I, I+320);
+
             g1.SetPoint(g1.GetN(), time*24*3600, Imed);
             g2.SetPoint(g2.GetN(), time*24*3600, Prediction(time));
+            g3.SetPoint(g3.GetN(), time*24*3600, Imed+Idev);
+            g4.SetPoint(g4.GetN(), time*24*3600, I[315]);
+            g5.SetPoint(g5.GetN(), time*24*3600, I[319]);
         }
 
     g1.SetMinimum(0);
-    g1.SetMaximum(99);
+    g1.SetMaximum(149);
 
     g1.SetMarkerStyle(kFullDotMedium);
@@ -260,13 +233,15 @@
     g1.DrawClone("AP");
 
+    g5.SetMarkerColor(kGray);
+    g5.DrawClone("P");
+
+    g4.SetMarkerColor(kGray+1);
+    g4.DrawClone("P");
+
+    g3.SetMarkerColor(kGray+2);
+    g3.DrawClone("P");
+
     g2.SetMarkerColor(kBlue);
     g2.SetMarkerStyle(kFullDotMedium);
-    g2.GetXaxis()->SetTimeDisplay(true);
-    g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
-    g2.GetXaxis()->SetLabelSize(0.12);
-    g2.GetYaxis()->SetLabelSize(0.1);
-    g2.GetYaxis()->SetTitle("CURRENT");
-    g2.GetYaxis()->SetTitleOffset(0.2);
-    g2.GetYaxis()->SetTitleSize(0.1);
     g2.DrawClone("P");
 
@@ -647,5 +622,5 @@
 }
 
-void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, TString outpath="./")
+void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, const char *outpath="quality")
 {
     // To get correct dates in the histogram you have to add
@@ -772,6 +747,4 @@
     cout << PlotTemperature4(runs, fname) << endl;
 
-    c->SaveAs(Form("%s/%04d%02d%02d-quality.png", outpath.Data(), y, m, d));
-}
-
-//20130314_141
+    c->SaveAs(Form("%s/%04d%02d%02d.png", outpath, y, m, d));
+}