Changeset 9953 for trunk/Cosy


Ignore:
Timestamp:
09/30/10 09:57:16 (14 years ago)
Author:
tbretz
Message:
Simplified and unified plotting macros.
Location:
trunk/Cosy/tpoint
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosy/tpoint/plot.C

    r8547 r9953  
     1#include <iostream>
    12#include <iomanip>
    2 
    3 Double_t GetResidual(Double_t fRawEl,  Double_t fRawAz,
    4                      Double_t fStarEl, Double_t fStarAz)
    5 {
    6     fRawEl *= TMath::DegToRad();
    7     fRawAz *= TMath::DegToRad();
    8     fStarEl *= TMath::DegToRad();
    9     fStarAz *= TMath::DegToRad();
    10     Double_t del = fRawEl-fStarEl;
    11     Double_t daz = fRawAz-fStarAz;
    12     Double_t dphi2 = daz/2.;
    13     Double_t cos2  = cos(dphi2)*cos(dphi2);
    14     Double_t sin2  = sin(dphi2)*sin(dphi2);
    15     Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
    16 
    17     Double_t dist = acos(d);
    18 
    19     return dist * TMath::RadToDeg();
    20 }
    21 
    22 void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
    23     {
    24         TView *view = pad->GetView();
    25 
    26         if (!view)
    27         {
    28             cout << "No View!" << endl;
    29             return;
    30         }
    31 
    32         TMarker mark0;
    33         TMarker mark1;
    34         mark0.SetMarkerStyle(kStar);
    35         mark1.SetMarkerStyle(kStar);
    36         mark1.SetMarkerColor(kRed);
    37    
    38         r0 /= 90;
    39         r1 /= 90;
    40         phi0 *= TMath::DegToRad();
    41         phi1 *= TMath::DegToRad();
    42    
    43         Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
    44         Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
    45 
    46         mark0.DrawMarker(x0[0], x0[1]);
    47         mark1.DrawMarker(x1[0], x1[1]);
    48 
    49         return;
    50         Double_t y0[3], y1[3];
    51    
    52         view->WCtoNDC(x0, y0);
    53         view->WCtoNDC(x1, y1);
    54    
    55         mark0.DrawMarker(y0[0], y0[1]);
    56         mark1.DrawMarker(y1[0], y1[1]);
    57     }
    58 
    59 int fill(const char *fname, TGraph *g, TH1 *h)
    60 {
    61 /*
    62     TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360,  10, 0, 90);
    63     h2res1.SetBit(TH1::kNoStats);
    64     h2res1.DrawCopy("surf1pol");
    65     gPad->Modified();
    66     gPad->Update();
    67     gPad->SetTheta(90);
    68     gPad->SetPhi(-90);
    69 
    70     DrawMarker(gPad, 45, 0, 0, 0);
    71     gPad->Modified();
    72     gPad->Update();
    73 
    74     return;
    75 
    76   */
     3#include <fstream>
     4
     5#include <TGraph.h>
     6#include <TCanvas.h>
     7#include <TSystem.h>
     8
     9#include "mars/MTime.h"
     10#include "mars/MMath.h"
     11#include "mars/MAstro.h"
     12#include "mars/MDirIter.h"
     13#include "mars/MStatusDisplay.h"
     14
     15using namespace std;
     16
     17int ReadFile(const char *fname, TGraph &g1, TGraph &g2, TGraph &g0)
     18{
    7719    ifstream fin(fname);
    7820
    79     cout << "Reading " << setw(23) << fname << "..." << flush;
    80 
     21    cout << "Reading " << setw(23) << gSystem->BaseName(fname) << "..." << flush;
     22
     23    MTime t;
    8124    while (1)
    8225    {
    … …  
    9033
    9134        Float_t alt, az, dalt, daz, mjd;
    92         sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
    93                &alt, &az, &dalt, &daz, &mjd);
    94 
    95         if (dalt==0/* ||  GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
    96             continue;
    97 
    98         //cout << dalt << " " << daz << " " << GetResidual(alt, az, alt+dalt, az+daz) << endl;
    99 
    100         mjd -=  53140.097505;
    101 
    102         Double_t res = GetResidual(alt, az, alt-dalt, az-daz);
    103 
    104         g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
    105         g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
    106         g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
    107 
    108         h->Fill(res);
    109     }
    110 
    111     cout << "done (" << setw(3) << (Int_t)h->GetEntries() << "/";
    112     cout << setw(4) << g[0].GetN() << ") " << flush;
    113 
    114     return g[0].GetN();
    115 }
    116 
    117 struct Description_t
    118 {
    119     const char *fName;
    120     const char *fTitle;
    121     const char *fFile;
    122 };
    123 
    124 const Int_t counts = 22;
    125 Description_t desc[counts] =
    126 {
    127     // No good pointing model has been applied yet
    128   //  {"0404",  "TPoints Residuals 4/2004" ,   "tpoint/tpoint0404.txt"},
    129   //  {"0405",  "TPoints Residuals 5/2004" ,   "tpoint/tpoint0405.txt"},
    130   //  {"04081", "TPoints Residuals 8/2004-1" , "tpoint/tpoint0408-1.txt"},
    131     {"04082", "TPoints Residuals 8/2004-2" , "tpoint/tpoint0408-2.txt"},
    132     {"0409",  "TPoints Residuals 9/2004" ,   "tpoint/tpoint0409.txt"},
    133     // Culmination tests
    134   //  {"0410",  "TPoints Residuals 10/2004" ,  "tpoint/tpoint0410.txt"},
    135     {"0411",  "TPoints Residuals 11/2004" ,  "tpoint/tpoint0411.txt"},
    136     {"0412",  "TPoints Residuals 12/2004" ,  "tpoint/tpoint0412.txt"},
    137     // Worse pointing due to realignment of the mirror
    138     {"0503",  "TPoints Residuals 3/2005" ,   "tpoint/tpoint0503.txt"},
    139     {"0504",  "TPoints Residuals 4/2005" ,   "tpoint/tpoint0504.txt"},
    140     {"05051", "TPoints Residuals 5/2005-1" , "tpoint/tpoint0505-1.txt"},
    141     // Mirror alignment has been fixed
    142     {"05052", "TPoints Residuals 5/2005-2" , "tpoint/tpoint0505-2.txt"},
    143     // Fixes to pointing model due to fixing a screw
    144     {"0506", "TPoints Residuals 6/2005" , "tpoint/tpoint0506.txt"},
    145 
    146     {"0508", "TPoints Residuals 8/2005" , "tpoint/tpoint0508.txt"},
    147     {"0509", "TPoints Residuals 9/2005" , "tpoint/tpoint0509.txt"},
    148     // Quick-and-dirty mirror alignment
    149     {"0510", "TPoints Residuals 10/2005" , "tpoint/tpoint0510.txt"},
    150     // New mirror alignment after Tenerife meeting
    151     {"0511", "TPoints Residuals 11/2005" , "tpoint/tpoint0511.txt"},
    152 
    153     {"0512", "TPoints Residuals 12/2005" , "tpoint/tpoint0512.txt"},
    154     {"0601", "TPoints Residuals 1/2006" ,  "tpoint/tpoint0601.txt"},
    155     {"0603", "TPoints Residuals 3/2006" ,  "tpoint/tpoint0603.txt"},
    156     {"0604", "TPoints Residuals 4/2006" ,  "tpoint/tpoint0604.txt"},
    157     //{"0605", "TPoints Residuals 5/2006" ,  "tpoint/tpoint0605.txt"},
    158     //{"0606", "TPoints Residuals 6/2006" ,  "tpoint/tpoint0606.txt"},
    159     {"0701", "TPoints Residuals 1/2007" ,  "tpoint/tpoint0701.txt"},
    160     {"0702", "TPoints Residuals 2/2007" ,  "tpoint/tpoint0702.txt"},
    161     {"0703", "TPoints Residuals 3/2007" ,  "tpoint/tpoint0703.txt"},
    162     {"0704", "TPoints Residuals 4/2007" ,  "tpoint/tpoint0704.txt"},
    163     {"0705", "TPoints Residuals 5/2007" ,  "tpoint/tpoint0705.txt"},
    164 };
    165 
    166 void plot()
    167 {
    168     TGraph g[3];
    169 
    170     MBinning bins(10, 0, 0.2);
    171 
    172     TArrayI n(counts);
    173 
    174     TH1F hx[counts];
    175     for (int i=0; i<counts; i++)
    176     {
    177         hx[i].SetNameTitle(desc[i].fName, desc[i].fTitle);
    178         hx[i].SetDirectory(0);
    179         bins.Apply(hx[i]);
    180 
    181         cout << setw(2) << i << ": " << flush;
    182         n[i] = fill(desc[i].fFile, g, &hx[i]);
    183         cout << " Mean: " << setw(5) << setprecision(2) << hx[i].GetMean() << " deg +/- " <<  hx[i].GetRMS()<<endl;
    184     }
    185 
    186     g[0].SetMarkerColor(kGreen);
    187     g[1].SetMarkerColor(kMagenta);
    188     g[2].SetMarkerColor(kBlack);
    189     //g[2].SetLineColor(kBlack);
    190     g[0].SetMarkerStyle(kFullDotMedium);
    191     g[1].SetMarkerStyle(kFullDotMedium);
    192     g[2].SetMarkerStyle(kFullDotLarge);
    193 
    194     // --------- First Canvas ----------
    195 
    196     new TCanvas("Vs Time", "");
    197 
    198     TObject *obj[4];
     35        if (sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
     36                   &az, &alt, &dalt, &daz, &mjd)!=5)
     37            continue;
     38
     39        if (dalt==0)
     40            continue;
     41
     42        Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz)*60;
     43
     44        if (res>12)
     45            continue;
     46
     47        t.SetMjd(mjd);
     48
     49        g0.SetPoint(g0.GetN(), t.GetAxisTime(), g0.GetN());
     50        g1.SetPoint(g1.GetN(), t.GetAxisTime(), res);
     51        g2.SetPoint(g2.GetN(), g2.GetN(),       res);
     52    }
     53
     54    cout << g1.GetN() << endl;
     55
     56    return g1.GetN();
     57}
     58
     59void DrawGrid(TH1 &h)
     60{
     61    TGraph l;
     62    l.SetLineStyle(kDashed);
     63    l.SetLineWidth(2);
     64    l.SetLineColor(kBlue);
     65
     66    Int_t style[] = { kSolid, 9, 7 };
     67
     68    Double_t x[2] = { h.GetXaxis()->GetXmin(), h.GetXaxis()->GetXmax() };
     69    Double_t y[2];
    19970
    20071    for (int i=0; i<3; i++)
    20172    {
    202         g[i].SetFillColor(kWhite);
    203         g[i].SetLineColor(kWhite);
    204         obj[i] = g[i].Clone();
    205         obj[i]->SetBit(kCanDelete);
    206         obj[i]->Draw(i==0?"AP":"P");
    207     }
    208 
    209     TLegend leg(0.905, 0.86, 0.99, 0.99);
    210     leg.AddEntry(obj[0], "  \\Delta\\theta");
    211     leg.AddEntry(obj[1], "  \\Delta\\phi");
    212     leg.AddEntry(obj[2], "  \\Delta");
    213     leg.DrawClone()->SetBit(kCanDelete);
    214 
    215     TLine l;
     73        y[0] = y[1] = i+1;
     74        l.SetLineStyle(style[i]);
     75        l.DrawGraph(2, x, y, "L");
     76    }
     77
     78    if (x[0]<1e6)
     79        return;
     80
     81    Double_t X[2];
     82    Double_t Y[2] = { h.GetMinimum(), h.GetMaximum() };
     83
     84    l.SetLineWidth(1);
     85
     86    for (int i=2004; i<2020; i++)
     87    {
     88        for (int j=0; j<12; j++)
     89        {
     90            X[0] = X[1] = MTime(i, 1+j, 1).GetAxisTime();
     91
     92            if (X[0]<x[0] || X[0]>x[1])
     93                continue;
     94
     95            l.SetLineStyle(j==0 ? kDashed : kDotted);
     96            l.DrawGraph(2, X, Y, "L");
     97        }
     98    }
     99}
     100
     101void DrawTimes(TH1 &h)
     102{
     103    TGraph l;
    216104    l.SetLineColor(kBlue);
    217     for (int i=0; i<n.GetSize(); i++)
    218         l.DrawLine(n[i], 0, n[i], 0.2);
    219 
    220     // --------- Second Canvas ----------
    221 
    222     new TCanvas("Distrib", "");
    223 
    224     Double_t max=0;
    225     for (int i=0; i<n.GetSize(); i++)
    226     {
    227         if (hx[i].GetEntries()==0)
    228         {
    229             cout << "Skip #" << i << endl;
    230             continue;
    231         }
    232 
    233         hx[i].Scale(1./hx[i].GetEntries());
    234         max = TMath::Max(max, hx[i].GetMaximum());
    235     }
    236     for (int i=0; i<n.GetSize(); i++)
    237     {
    238         hx[i].SetMaximum(max*1.05);
    239         if (i<6)
    240             hx[i].SetLineColor(kRed+i);
    241         else
    242             hx[i].SetMarkerStyle(kPlus+i-6);
    243 
    244 
    245     }
    246 
    247     for (int i=0; i<counts; i++)
    248         hx[i].DrawCopy(i==0?"LP":"LPsame");
    249 
    250     return;
    251 
    252     for (int i=0; i<n.GetSize(); i++)
    253     {
    254         cout << "Mean:  " << Form("%.3f +- %.3f", hx[i].GetMean(), hx[i].GetRMS());
    255         cout << "   (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")" << endl;
    256     }
    257 }
     105
     106    Double_t x[2];
     107    Double_t y[2] = { h.GetMinimum(), h.GetMaximum() };
     108
     109    for (int i=0; dates[i+1]>0; i++)
     110    {
     111        x[0] = x[1] = dates[i];
     112        l.DrawGraph(2, x, y, "L");
     113    }
     114
     115    DrawGrid(h);
     116}
     117
     118Int_t GetBin(TGraph &times, Int_t i)
     119{
     120    return TMath::BinarySearch(times.GetN(), times.GetX(), dates[i]);;
     121}
     122
     123void DrawIndices(TGraph &times, TH1 &h)
     124{
     125    TGraph l;
     126    l.SetLineStyle(kDashed);
     127    l.SetLineColor(kBlue);
     128
     129    Double_t x[2];
     130    Double_t y[2] = { h.GetMinimum(), h.GetMaximum() };
     131
     132    for (int i=0; dates[i+1]>0; i++)
     133    {
     134        Int_t bin = GetBin(times, i);
     135        if (bin<0)
     136            continue;
     137
     138        x[0] = x[1] = bin+0.5;
     139        l.DrawGraph(2, x, y, "L");
     140    }
     141
     142    DrawGrid(h);
     143}
     144
     145void plot(MDirIter &Next)
     146{
     147    Next.Sort();
     148
     149    TGraph g1;
     150    TGraph g2;
     151    TGraph g0;
     152
     153    while (1)
     154    {
     155        TString name=Next();
     156        if (name.IsNull())
     157            break;
     158
     159        ReadFile(name, g1, g2, g0);
     160    }
     161
     162    MStatusDisplay *d = new MStatusDisplay;
     163
     164    gROOT->SetSelectedPad(0);
     165    TCanvas &c0 = d->AddTab("TimeLine");
     166    c0.SetGrid();
     167
     168    g0.SetTitle("TimeLine");
     169    g0.GetXaxis()->SetTimeDisplay(1);
     170    g0.SetMarkerStyle(kFullDotMedium);
     171    g0.SetLineWidth(1);
     172    g0.SetLineColor(kGreen);
     173    g0.DrawClone("ALP");
     174
     175    DrawTimes(*g0.GetHistogram());
     176
     177
     178    gROOT->SetSelectedPad(0);
     179    TCanvas &c1 = d->AddTab("ResDate");
     180    //c1.SetGrid();
     181
     182    g1.SetTitle("Measured Residual vs. Date");
     183    g1.GetXaxis()->SetTimeDisplay(1);
     184//    g1.SetMaximum(10);
     185    g1.SetMarkerStyle(kFullDotMedium);
     186    g1.DrawClone("AP");
     187
     188    DrawTimes(*g1.GetHistogram());
     189
     190
     191    gROOT->SetSelectedPad(0);
     192    TCanvas &c2 = d->AddTab("ResIdx");
     193    c2.SetGrid();
     194
     195    g2.SetTitle("Measured Residual vs. Index");
     196//    g2.SetMaximum(0.1);
     197    g2.SetMarkerStyle(kFullDotMedium);
     198    g2.DrawClone("AP");
     199
     200    DrawIndices(g0, *g2.GetHistogram());
     201
     202    Double_t prob[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
     203    Double_t quant[4];
     204
     205    int n;
     206    for (n=0; dates[n]>0; n++);
     207
     208    TH1F h0("Q98", "", n-1, dates);
     209    TH1F h1("Q95", "", n-1, dates);
     210    TH1F h2("Q68", "", n-1, dates);
     211    TH1F h3("Q50", "", n-1, dates);
     212
     213    Int_t first = 0;
     214
     215    cout << "n=" << n << endl;
     216
     217    for (int i=0; dates[i+1]>0; i++)
     218    {
     219        Int_t bin0 = GetBin(g0, i);
     220        Int_t bin1 = GetBin(g0, i+1);
     221
     222        cout << i << " " << bin0 << " " << bin1 << endl;
     223        if (bin0<0)
     224            continue;
     225
     226        if (bin0==bin1)
     227            continue;
     228
     229        TMath::Quantiles(bin1-bin0, 4, g2.GetY()+bin0, quant, prob, kFALSE, NULL, 7);
     230
     231        h0.SetBinContent(i+1, quant[0]);
     232        h1.SetBinContent(i+1, quant[1]);
     233        h2.SetBinContent(i+1, quant[2]);
     234        h3.SetBinContent(i+1, quant[3]);
     235
     236        if (first==0)
     237            first = i;
     238    }
     239
     240
     241    gROOT->SetSelectedPad(0);
     242    TCanvas &c4 = d->AddTab("Quantiles");
     243    //c4.SetGrid();
     244
     245    h2.SetFillColor(19);
     246    h1.SetFillColor(16);
     247    h0.SetFillColor(13);
     248
     249    h2.SetTitle("Measured Residual-Distribution vs. Date");
     250    h2.SetYTitle("Residual [arcmin]");
     251    h2.SetStats(kFALSE);
     252    h2.GetXaxis()->SetTimeDisplay(1);
     253    h2.GetXaxis()->SetTimeFormat("%y/%m");
     254    h2.GetYaxis()->CenterTitle();
     255    h2.SetMinimum(0);
     256    //.SetMaximum(h2.GetMaximum()*1.2);
     257    h2.SetMarkerStyle(kFullDotMedium);
     258    h2.GetXaxis()->SetRange(first, h3.GetNbinsX());
     259    h2.DrawCopy();
     260
     261    //.DrawCopy("same");
     262    h1.DrawCopy("same");
     263    h0.DrawCopy("same");
     264
     265    TGraph g;
     266    for (int i=0; i<h2.GetNbinsX(); i++)
     267    {
     268        Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };
     269        Double_t y[2] = { 0, h2.GetBinContent(i+1) };
     270
     271        g.DrawGraph(2, x, y, "L");
     272    }
     273
     274    DrawGrid(h2);
     275
     276    gROOT->SetSelectedPad(0);
     277    TCanvas &c5 = d->AddTab("Combined");
     278
     279    h2.DrawCopy();
     280    h1.DrawCopy("same");
     281    h0.DrawCopy("same");
     282
     283    for (int i=0; i<h2.GetNbinsX(); i++)
     284    {
     285        Double_t x[2] = { h2.GetBinLowEdge(i+1), h2.GetBinLowEdge(i+1) };
     286        Double_t y[2] = { 0, h2.GetBinContent(i+1) };
     287
     288        g.DrawGraph(2, x, y, "L");
     289    }
     290
     291    g1.DrawClone("P");
     292
     293    DrawGrid(h2);
     294
     295
     296}

  • trunk/Cosy/tpoint/plot_m1.C

    r9512 r9953  
    1 #include <iomanip>
     1#include "mars/MTime.h"
    22
    3 void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
    4     {
    5         TView *view = pad->GetView();
     3using namespace std;
    64
    7         if (!view)
    8         {
    9             cout << "No View!" << endl;
    10             return;
    11         }
    12 
    13         TMarker mark0;
    14         TMarker mark1;
    15         mark0.SetMarkerStyle(kStar);
    16         mark1.SetMarkerStyle(kStar);
    17         mark1.SetMarkerColor(kRed);
    18    
    19         r0 /= 90;
    20         r1 /= 90;
    21         phi0 *= TMath::DegToRad();
    22         phi1 *= TMath::DegToRad();
    23    
    24         Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
    25         Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
    26 
    27         mark0.DrawMarker(x0[0], x0[1]);
    28         mark1.DrawMarker(x1[0], x1[1]);
    29 
    30         return;
    31         Double_t y0[3], y1[3];
    32    
    33         view->WCtoNDC(x0, y0);
    34         view->WCtoNDC(x1, y1);
    35    
    36         mark0.DrawMarker(y0[0], y0[1]);
    37         mark1.DrawMarker(y1[0], y1[1]);
    38     }
    39 
    40 int fill(const char *fname, TGraph *g, TH1 *h)
    41 {
    425/*
    43     TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360,  10, 0, 90);
    44     h2res1.SetBit(TH1::kNoStats);
    45     h2res1.DrawCopy("surf1pol");
    46     gPad->Modified();
    47     gPad->Update();
    48     gPad->SetTheta(90);
    49     gPad->SetPhi(-90);
    50 
    51     DrawMarker(gPad, 45, 0, 0, 0);
    52     gPad->Modified();
    53     gPad->Update();
    54 
    55     return;
    56 
    57   */
    58     ifstream fin(fname);
    59 
    60     cout << "Reading " << setw(23) << gSystem->BaseName(fname) << "..." << flush;
    61 
    62     while (1)
    63     {
    64         TString str;
    65         str.ReadLine(fin);
    66         if (!fin)
    67             break;
    68 
    69         if (str.Contains("#"))
    70             continue;
    71 
    72         Float_t alt, az, dalt, daz, mjd;
    73         sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
    74                &az, &alt, &dalt, &daz, &mjd);
    75 
    76         if (dalt==0/* ||  GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
    77             continue;
    78 
    79         mjd -=  53140.097505;
    80 
    81         Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz);
    82 
    83         g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
    84         g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
    85         g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
    86 
    87         h->Fill(res);
    88     }
    89 
    90     cout << "done (" << setw(3) << (Int_t)h->GetEntries() << "/";
    91     cout << setw(4) << g[0].GetN() << ") " << flush;
    92 
    93     return g[0].GetN();
    94 }
    95 
    96 struct Description_t
    97 {
    98     const char *fName;
    99     const char *fTitle;
    100     const char *fFile;
    101 };
    102 
    103 const Int_t counts = 29+10+18+1+13+11+31-27+3+5+18;
    104 Description_t desc[counts] =
    105 {
    106     /*
    107      //   29. Apr. 2004    ~25800
    108      //    5. Aug. 2004    ~32000
    109      //   19. Aug. 2004    ~33530
    110      //    7. Jun. 2005    ~57650
    111      //    8. Jun. 2005
    112      //    9. Jun. 2005    ~57860
    113      //   12. Sep. 2005    ~68338
    114      //   24. Nov. 2005    ~75562
    115      //   17. Oct. 2006   ~103130
    116      //   17. Jun. 2007   ~248193
    117      */
    118 
    1196    // Culmination tests
    1207    //{"0411",  "TPoints Residuals 11/2004" ,  "tpoint/tpoint0411.txt"},
    … …  
    14633    //{"05111", "TPoints Residuals 11/2005-1" , "tpoint/tpoint0511-1.txt"},
    14734
    148     // New pointing model installed (24.11.2005)
    149     {"05112", "TPoints Residuals 11/2005-2" , "tpoint/tpoint0511-2.txt"},
    150     {"+0512", "TPoints Residuals 12/2005" , "tpoint/tpoint0512.txt"},
    151     {"+0601", "TPoints Residuals 1/2006" ,  "tpoint/tpoint0601.txt"},
    152     {"+0603-1", "TPoints Residuals 3/2006-1" ,  "tpoint/tpoint0603-1.txt"},
     35    // [...]
     36*/
    15337
    154     // Changes to the mirror
    155     {"0603-2", "TPoints Residuals 3/2006-2" ,  "tpoint/tpoint0603-2.txt"},
    156     {"+0604", "TPoints Residuals 4/2006" ,  "tpoint/tpoint0604.txt"},
    157     {"+0607", "TPoints Residuals  7/2006" ,  "tpoint/tpoint0607.txt"},
    158     {"+0608", "TPoints Residuals  8/2006" ,  "tpoint/tpoint0608.txt"},
    159     {"+0609", "TPoints Residuals  9/2006" ,  "tpoint/tpoint0609.txt"},
    160     {"+0610", "TPoints Residuals 10/2006" ,  "tpoint/tpoint0610.txt"},
     38Double_t dates[] = {
     39    MTime(2005,  3, 20).GetAxisTime(),
     40    MTime(2005,  4, 29).GetAxisTime(),
     41    MTime(2005,  5, 25).GetAxisTime(),
     42    MTime(2005,  6,  8).GetAxisTime(), // New pointing model
     43    MTime(2005,  8, 15).GetAxisTime(),
     44    //        MTime(2005,  9, 12).GetAxisTime(), // New pointing model
     45    MTime(2005, 11, 10).GetAxisTime(), // New mirror alignment after Tenerife meeting
     46    MTime(2005, 11, 24).GetAxisTime(), // New pointing model
     47    MTime(2006,  3, 19).GetAxisTime(), // Changes to the mirror
     48    //    2006,  4, 23                 // Mirror refocussing
     49    MTime(2006, 10, 17).GetAxisTime(), // New pointing model
     50    MTime(2007,  6, 17).GetAxisTime(), // New pointing model
     51    MTime(2007,  8,  4).GetAxisTime(), // Mirror refocussing
     52    MTime(2007, 10, 18).GetAxisTime(), // New pointing model
     53    MTime(2008,  1, 14).GetAxisTime(), // New pointing model
     54    MTime(2008,  6, 11).GetAxisTime(), // New pointing model
     55    MTime(2008,  6, 18).GetAxisTime(), // New pointing model
     56    //    2009,  3,  7                 // New pointing model
     57    //    2009,  5, 14                 // New pointing model
     58    // Are we missing TPoints between 080618 and 090501??
     59    MTime(2009,  5,  1).GetAxisTime(), // Drive upgrade started
     60    MTime(2009,  5, 11).GetAxisTime(), // Upgrade finished
     61    MTime(2009,  5, 12).GetAxisTime(), // First new pointing model
     62    MTime(2009,  5, 13).GetAxisTime(), // Second new pointing model
     63    MTime(2009,  6, 11).GetAxisTime(), // Jump (reason unknown)
     64    MTime(2009,  7, 23).GetAxisTime(), // New LUTs
     65    MTime(2009,  8, 17).GetAxisTime(), // New pointing model
     66    MTime(2010, 02, 01).GetAxisTime(), // New pointing model
     67    //    2010, 02, 03                 // New LUTs for M2(!)
     68    MTime(2010, 02, 27).GetAxisTime(), // New pointing model
     69    MTime(2010, 03, 31).GetAxisTime(), // New pointing model
     70    //    2010, 06, 14                 // New LUTs
     71    MTime(2010, 07, 03).GetAxisTime(), // New LUTs
     72    MTime(2010,  9, 29).GetAxisTime(), // New pointing model
     73    MTime(-1).GetAxisTime(),
     74    -1
     75};
    16176
    162     // New pointing model: 6/10/17
    163     {"0611", "TPoints Residuals 11/2006" ,  "tpoint/tpoint0611.txt"},
    164     {"+0612", "TPoints Residuals 12/2006" ,  "tpoint/tpoint0612.txt"},
    165     {"+0701", "TPoints Residuals  1/2007" ,  "tpoint/tpoint0701.txt"},
    166     {"+0702", "TPoints Residuals  2/2007" ,  "tpoint/tpoint0702.txt"},
    167     {"+0703", "TPoints Residuals  3/2007" ,  "tpoint/tpoint0703.txt"},
    168 
    169     {"+0704", "TPoints Residuals  4/2007" ,  "tpoint/tpoint0704.txt"},
    170     {"+0705", "TPoints Residuals  5/2007" ,  "tpoint/tpoint0705.txt"},
    171     {"+0706", "TPoints Residuals  6/2007" ,  "tpoint/tpoint0706.txt"},
    172     {"+07071", "TPoints Residuals  7/2007-1" ,  "tpoint/tpoint0707-1.txt"},
    173     {"+07072", "TPoints Residuals  7/2007-2" ,  "tpoint/tpoint0707-2.txt"},
    174 
    175     // New pointing model: 07/06/17
    176     {"",      "",                         ""},
    177 
    178     // AMC adjust:         07/08/04
    179     {"0708",  "TPoints Residuals 7/8/23", "tpoint/m1/tpoints20070823.txt"},
    180 
    181     // New pointing model: 07/10/18
    182     {"",      "",                         ""},
    183 
    184     // New pointing model: 08/01/14
    185     {"0801",  "TPoints Residuals 7/8/23", "tpoint/m1/tpoints20080115.txt"},
    186 
    187     // New pointing model 080611
    188     {"0806", "TPoints Residuals 7/8/23", "tpoint/tpoint0806.txt"},
    189 
    190     // New pointing model 080618
    191     {"0807",  "",                         "tpoint/tpoint0807.txt"},
    192     {"+0807", "",                         "tpoint/tpoint0808.txt"},
    193     {"+0807", "",                         "tpoint/tpoint0811.txt"},
    194     {"+0807", "",                         "tpoint/tpoint0812.txt"},
    195     {"+0901", "",                         "tpoint/tpoint0901.txt"},
    196 
    197     // Drive upgrade started
    198     {"", "",                         ""},
    199 
    200     // Drive Upgrade finished 090511
    201     {"090512", "TPoints Residuals 7/8/23",  "tpoint/m1/2009_05_12/tpoint_20090511_220118.txt"},
    202     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090511_221333.txt"},
    203     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090511_222434.txt"},
    204     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090511_223152.txt"},
    205     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090511_223706.txt"},
    206     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090511_224056.txt"},
    207     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090511_231933.txt"},
    208     {"+090512", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_12/tpoint_20090512_011131.txt"},
    209 
    210     // First new pointing model 090512
    211     {"090513", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_13/tpoint_20090512_210644.txt"},
    212     {"+090513", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_13/tpoint_20090513_002757.txt"},
    213     {"+090513", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_13/tpoint_20090513_025124.txt"},
    214     {"+090513", "TPoints Residuals 7/8/23", "tpoint/m1/2009_05_13/tpoint_20090513_033505.txt"},
    215 
    216     // Second new pointing model 090513
    217     {"090514", "TPoints 09/05/14", "tpoint/m1/2009_05_14/tpoint_20090514_011903.txt"},
    218 
    219     {"090517", "TPoints 09/05/17" ,"tpoint/m1/2009_05_17/tpoint_20090516_234825.txt"},
    220     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_05_17/tpoint_20090517_023340.txt"},
    221     // New period
    222     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_01/tpoint_20090531_215148.txt"},
    223     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_01/tpoint_20090531_222549.txt"},
    224     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_02/tpoint_20090601_223009.txt"},
    225     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_03/tpoint_20090602_213509.txt"},
    226     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_03/tpoint_20090603_011936.txt"},
    227     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_04/tpoint_20090603_215840.txt"},
    228     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_04/tpoint_20090603_230510.txt"},
    229     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_05/tpoint_20090604_215943.txt"},
    230     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_05/tpoint_20090604_232320.txt"},
    231     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_10/tpoint_20090609_231948.txt"},
    232     {"+090517","TPoints 09/05/17", "tpoint/m1/2009_06_11/tpoint_20090611_011155.txt"},
    233 
    234     // ------- Something happened (reason unknown) --------
    235 
    236     {"090613","TPoints 09/05/17", "tpoint/m1/2009_06_11/tpoint_20090611_023625.txt"},
    237 
    238     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_06_13/tpoint_20090613_033838.txt"},
    239     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_06_14/tpoint_20090614_021257.txt"},
    240     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_06_18/tpoint_20090618_041433.txt"},
    241 
    242     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_02/tpoint_20090701_215304.txt"},
    243     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_02/tpoint_20090701_222059.txt"},
    244     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_02/tpoint_20090701_224051.txt"},
    245     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_02/tpoint_20090701_225615.txt"},
    246     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_02/tpoint_20090701_230946.txt"},
    247 
    248     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_03/tpoint_20090702_225940.txt"},
    249     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_04/tpoint_20090703_224721.txt"},
    250     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_05/tpoint_20090705_012638.txt"},
    251 
    252     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_12/tpoint_20090712_014300.txt"},
    253     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_12/tpoint_20090712_024710.txt"},
    254     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_13/tpoint_20090713_004241.txt"},
    255     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_13/tpoint_20090713_025934.txt"},
    256 
    257     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_14/tpoint_20090714_024729.txt"},
    258     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_15/tpoint_20090715_021320.txt"},
    259     {"+090613","TPoints 09/05/17", "tpoint/m1/2009_07_15/tpoint_20090715_025237.txt"},
    260 
    261     // ----- AMC adjust (23.7.) -----
    262 
    263     {"090723","TPoints 09/05/17", "tpoint/m1/2009_07_24/tpoint_20090724_051142.txt"},
    264     {"+090723", "TPoints 08/2009" , "tpoint/m1/2009_08_02/tpoint_20090801_223848.txt"},
    265     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_03/tpoint_20090802_224434.txt"},
    266     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_04/tpoint_20090803_230214.txt"},
    267     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_10/tpoint_20090809_233524.txt"},
    268     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_10/tpoint_20090810_014642.txt"},
    269     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_10/tpoint_20090810_020616.txt"},
    270     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_11/tpoint_20090810_232543.txt"},
    271     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_12/tpoint_20090812_001846.txt"},
    272     {"+090723","TPoints 08/2009" , "tpoint/m1/2009_08_13/tpoint_20090813_030807.txt"},
    273 
    274     // ------ 09/08/17 new pointing model -------
    275 
    276     {"090817","TPoints", "tpoint/m1/2009_09_01/tpoint_20090831_233034.txt"},
    277     {"+090817","TPoints", "tpoint/m1/2009_09_02/tpoint_20090901_231918.txt"},
    278     {"+090817","TPoints", "tpoint/m1/2009_09_02/tpoint_20090902_005204.txt"},
    279     {"+090817","TPoints", "tpoint/m1/2009_09_08/tpoint_20090907_213404.txt"},
    280     {"+090817","TPoints", "tpoint/m1/2009_09_09/tpoint_20090908_225027.txt"},
    281     {"+090817","TPoints", "tpoint/m1/2009_09_09/tpoint_20090908_231720.txt"},
    282     {"+090817","TPoints", "tpoint/m1/2009_09_09/tpoint_20090909_010649.txt"},
    283     {"+090817","TPoints", "tpoint/m1/2009_09_10/tpoint_20090909_235222.txt"},
    284     {"+090817","TPoints", "tpoint/m1/2009_09_10/tpoint_20090910_025054.txt"},
    285     {"+090817","TPoints", "tpoint/m1/2009_09_12/tpoint_20090912_004826.txt"},
    286     {"+090817","TPoints", "tpoint/m1/2009_09_17/tpoint_20090917_015404.txt"},
    287     {"+090817","TPoints", "tpoint/m1/2009_09_18/tpoint_20090917_202449.txt"},
    288     {"+090817","TPoints", "tpoint/m1/2009_09_23/tpoint_20090922_223445.txt"},
    289     {"+090817","TPoints", "tpoint/m1/2009_09_24/tpoint_20090923_204444.txt"},
    290     {"+090817","TPoints", "tpoint/m1/2009_09_24/tpoint_20090924_041824.txt"},
    291     {"+090817","TPoints", "tpoint/m1/2009_09_25/tpoint_20090924_225725.txt"},
    292     {"+090817","TPoints", "tpoint/m1/2009_09_26/tpoint_20090925_212537.txt"},
    293     {"+090817","TPoints", "tpoint/m1/2009_09_27/tpoint_20090926_212332.txt"},
    294 };
     77#include "plot.C"
    29578
    29679void plot_m1()
    29780{
    298     TGraph g[3];
    29981
    300     MBinning bins(100, 0, 0.2);
     82    MDirIter Next;
     83    Next.AddDirectory("tpoint/m1", "tpoint*.txt", -1);
    30184
    302     TArrayI n(counts);
    303 
    304     TH1F hx[counts];
    305     Int_t num = -1;
    306     for (int i=0; i<counts; i++)
    307     {
    308         if (desc[i].fName[0]!='+')
    309         {
    310             num++;
    311 
    312             hx[num].SetNameTitle(desc[i].fName, desc[i].fTitle);
    313             hx[num].SetDirectory(0);
    314             bins.Apply(hx[num]);
    315         }
    316 
    317         cout << setw(2) << num << ": " << flush;
    318         n[num] = fill(desc[i].fFile, g, &hx[num]);
    319         cout << " Mean: " << setw(5) << setprecision(3) << hx[num].GetMean()*3600 << " sec +/- " <<  hx[num].GetRMS()*3600 << endl;
    320     }
    321 
    322     n.Set(++num);
    323 
    324     g[0].SetMarkerColor(kGreen);
    325     g[1].SetMarkerColor(kMagenta);
    326     g[2].SetMarkerColor(kBlack);
    327     //g[2].SetLineColor(kBlack);
    328     g[0].SetMarkerStyle(kFullDotMedium);
    329     g[1].SetMarkerStyle(kFullDotMedium);
    330     g[2].SetMarkerStyle(kFullDotLarge);
    331 
    332     // --------- First Canvas ----------
    333 
    334     new TCanvas("Vs Time", "");
    335 
    336     TObject *obj[4];
    337 
    338     for (int i=0; i<3; i++)
    339     {
    340         g[i].SetFillColor(kWhite);
    341         g[i].SetLineColor(kWhite);
    342         obj[i] = g[i].Clone();
    343         obj[i]->SetBit(kCanDelete);
    344         obj[i]->Draw(i==0?"AP":"P");
    345     }
    346 
    347     TLegend leg(0.905, 0.86, 0.99, 0.99);
    348     leg.AddEntry(obj[0], "  \\Delta\\theta");
    349     leg.AddEntry(obj[1], "  \\Delta\\phi");
    350     leg.AddEntry(obj[2], "  \\Delta");
    351     leg.DrawClone()->SetBit(kCanDelete);
    352 
    353     TLine l;
    354     l.SetLineColor(kBlue);
    355     for (int i=0; i<n.GetSize(); i++)
    356         l.DrawLine(n[i], 0, n[i], 0.2);
    357 
    358     // --------- Second Canvas ----------
    359 
    360     new TCanvas("Distrib", "");
    361 
    362     Double_t max=0;
    363     for (int i=0; i<n.GetSize(); i++)
    364     {
    365         if (hx[i].GetEntries()==0)
    366         {
    367             cout << "Skip #" << i << endl;
    368             continue;
    369         }
    370 
    371         hx[i].Scale(1./hx[i].GetEntries());
    372         max = TMath::Max(max, hx[i].GetMaximum());
    373     }
    374     for (int i=0; i<n.GetSize(); i++)
    375     {
    376         hx[i].SetMaximum(max*1.05);
    377         if (i<6)
    378             hx[i].SetLineColor(kRed+i);
    379         else
    380             hx[i].SetMarkerStyle(kPlus+i-6);
    381 
    382 
    383     }
    384 
    385     for (int i=0; i<counts; i++)
    386         hx[i].DrawCopy(i==0?"LP":"LPsame");
    387 
    388     Double_t time[] = {
    389         MTime(2005,  3, 20).GetAxisTime(),
    390         MTime(2005,  4, 29).GetAxisTime(),
    391         MTime(2005,  5, 25).GetAxisTime(),
    392         MTime(2005,  6,  8).GetAxisTime(),
    393         MTime(2005,  8, 15).GetAxisTime(),
    394 //        MTime(2005,  9, 12).GetAxisTime(),
    395         MTime(2005, 11, 10).GetAxisTime(),
    396         MTime(2005, 11, 24).GetAxisTime(),
    397         MTime(2006,  3, 19).GetAxisTime(),
    398         MTime(2006, 10, 17).GetAxisTime(),
    399         MTime(2007,  6, 17).GetAxisTime(),
    400         MTime(2007,  8,  4).GetAxisTime(),
    401         MTime(2007, 10, 18).GetAxisTime(),
    402         MTime(2008,  1, 14).GetAxisTime(),
    403         MTime(2008,  6, 11).GetAxisTime(),
    404         MTime(2008,  6, 18).GetAxisTime(),
    405         // Are we missing TPoints between 080618 and 090501??
    406         MTime(2009,  5, 1).GetAxisTime(),
    407         MTime(2009,  5, 11).GetAxisTime(),
    408         MTime(2009,  5, 12).GetAxisTime(),
    409         MTime(2009,  5, 13).GetAxisTime(),
    410         MTime(2009,  5, 16).GetAxisTime(),
    411         MTime(2009,  6, 11).GetAxisTime(),
    412         MTime(2009,  7, 23).GetAxisTime(),
    413         MTime(2009,  8, 17).GetAxisTime(),
    414         MTime(2009, 12, 31).GetAxisTime(),
    415     };
    416 
    417     TH1D histres[4];
    418 
    419     MBinning bins;
    420     bins.SetEdges(TArrayD(24, time));
    421 
    422     bins.Apply(histres[0]);
    423     bins.Apply(histres[1]);
    424     bins.Apply(histres[2]);
    425     bins.Apply(histres[3]);
    426 
    427     TGraphAsymmErrors result[4];
    428     TGraph resultm;
    429     for (int i=0; i<n.GetSize(); i++)
    430     {
    431         cout << i+1 << " - Mean:  " << Form("%5.1f +- %5.1f", hx[i].GetMean()*3600, hx[i].GetRMS()*3600);
    432         cout << "   (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")   " << (int)hx[i].GetEntries() << endl;
    433 
    434         if (hx[i].GetEntries()<1)
    435             continue;
    436 
    437         /*
    438          TF1 fg("fg", "gaus", 0, 1);
    439          hx[i].Fit(&fg);
    440          result2.SetPoint(result.GetN(), result.GetN()+7, hx[i].GetMean());
    441          result2.SetPointError(result.GetN()-1, 0, hx[i].GetRMS()/sqrt(hx[i].GetEntries()));
    442          result.SetPoint(result.GetN(), result.GetN()+7, fg.GetParameter(1));
    443          result.SetPointError(result.GetN()-1, 0, fg.GetParameter(2));
    444          */
    445 
    446         //Double_t q[4] = { MMath::GaussProb(0.5), MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
    447         Double_t q[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
    448         Double_t rc[4];
    449         hx[i].GetQuantiles(4, rc, q);
    450 
    451         for (int j=0; j<4; j++)
    452         {
    453             histres[j].SetBinContent(i+1, 60*rc[j]);
    454 
    455             result[j].SetPoint(i, time[i], 60*rc[j]);
    456             result[j].SetPointError(i, 0, 0, 60*rc[j], 0);
    457 
    458             //result[j].SetPoint(result[j].GetN(), result[j].GetN()+1, 60*rc[j]);
    459             //result[j].SetPointError(result[j].GetN()-1, 0, 0, 60*rc[j], 0);
    460         }
    461 
    462 //        result.SetPoint(result.GetN(), result.GetN()+1, rc[1]);
    463 //        result2.SetPoint(result2.GetN(), result2.GetN()+1, rc[2]);
    464 
    465         //result.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
    466         //result2.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
    467 
    468 //        result.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
    469 //        result2.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
    470 
    471 //        result2.SetPointError(result.GetN()-1, 0, 0, rc[2], 0);
    472 
    473 //        resultm.SetPoint(resultm.GetN(), resultm.GetN()+1, 60*hx[i].GetMean());
    474         resultm.SetPoint(resultm.GetN(), (time[i]+time[i+1])/2, 60*hx[i].GetMean());
    475 
    476     }
    477 
    478     new TCanvas;
    479 
    480     gPad->SetBorderMode(0);
    481     gPad->SetFrameBorderMode(0);
    482     gPad->SetFillColor(kWhite);
    483     gPad->SetRightMargin(0.01);
    484     gPad->SetTopMargin(0.01);
    485     gPad->SetLeftMargin(0.06);
    486     gPad->SetGridy();
    487 
    488     //Int_t col[] = { 12, 15, 17, 19 };
    489     //Int_t col[] = { 12, 16, 18, 0 };
    490     Int_t col[] = { 13, 16, 19, 0 };
    491 
    492     TH1 *h = &histres[2];//result[3].GetHistogram();
    493     h->SetXTitle("");
    494     h->SetYTitle("Residual / arcmin");
    495     h->SetBit(TH1::kNoStats);
    496     h->GetXaxis()->CenterTitle();
    497     h->GetYaxis()->CenterTitle();
    498     h->GetYaxis()->SetTitleOffset(0.75);
    499 //    h->GetXaxis()->SetTimeFormat("%m/%y %F1995-01-01 00:00:00 GMT");
    500 //    h->GetXaxis()->SetTimeDisplay(1);
    501 
    502     h->GetXaxis()->SetLabelColor(kWhite);
    503 
    504     TLine line;
    505 
    506     for (int j=2; j>=0; j--)
    507     {
    508         histres[j].SetMinimum(0);
    509         histres[j].SetFillColor(col[j]);//12+2*j);
    510 
    511         histres[j].DrawCopy(j==2?"":"same");
    512 
    513         /*
    514         //result[j].SetLineColor(kBlue);
    515         //result[j].SetLineWidth(2);
    516         //result[j].SetMarkerColor(kBlue);
    517         result[j].SetMinimum(0);
    518         //result2.SetMarkerStyle(kFullDotMedium);
    519         //result[j].SetMarkerStyle(23);
    520         result[j].SetFillColor(col[j]);//12+2*j);
    521         result[j].DrawClone(j==3 ? "ABX" : "B"); // E3 B
    522         */
    523     }
    524 
    525 
    526     resultm.SetMarkerStyle(20);
    527     resultm.SetMarkerSize(0.8);
    528     resultm.DrawClone("P");
    529 
    530     line.DrawLine(time[0], 0, time[0], histres[2].GetMaximum());
    531     for (int i=0; i<bins.GetNumBins()-1; i++)
    532         line.DrawLine(time[i+1], 0, time[i+1], histres[2].GetBinContent(i+1));
    533 
    534     TText txt;
    535     txt.SetTextSize(0.037);
    536 //    txt.SetTextAngle(-45);
    537     txt.SetTextAlign(23);
    538     for (int m=4; m<13; m++)
    539     {
    540         Double_t monl = MTime(2005, m,   1,0).GetAxisTime();
    541         Double_t mont = MTime(2005, m,   15,0).GetAxisTime();
    542 //        txt.DrawText(mon, -0.12, Form("%02d/05", m));
    543         txt.DrawText(mont, -0.10, Form("%d", m));
    544         line.DrawLine(monl, -0.12, monl, 0.12);
    545     }
    546     for (int m=1; m<13; m++)
    547     {
    548         Double_t monl = MTime(2006, m,   1,0).GetAxisTime();
    549         Double_t mont = MTime(2006, m,   15,0).GetAxisTime();
    550 //        txt.DrawText(mon, -0.12, Form("%02d/06", m));
    551         txt.DrawText(mont, -0.10, Form("%d", m));
    552         line.DrawLine(monl, -0.12, monl, 0.12);
    553     }
    554     for (int m=1; m<13; m++)
    555     {
    556         Double_t monl = MTime(2007, m,   1,0).GetAxisTime();
    557         Double_t mont = MTime(2007, m,   15,0).GetAxisTime();
    558 //        txt.DrawText(mon, -0.12, Form("%02d/07", m));
    559         txt.DrawText(mont, -0.10, Form("%d", m));
    560         line.DrawLine(monl, -0.12, monl, 0.12);
    561     }
    562     for (int m=1; m<13; m++)
    563     {
    564         Double_t monl = MTime(2008, m,   1,0).GetAxisTime();
    565         Double_t mont = MTime(2008, m,   15,0).GetAxisTime();
    566 //        txt.DrawText(mon, -0.12, Form("%02d/07", m));
    567         txt.DrawText(mont, -0.10, Form("%d", m));
    568         line.DrawLine(monl, -0.12, monl, 0.12);
    569     }
    570 
    571     for (int m=1; m<13; m++)
    572     {
    573         Double_t monl = MTime(2009, m,   1,0).GetAxisTime();
    574         Double_t mont = MTime(2009, m,   15,0).GetAxisTime();
    575 //        txt.DrawText(mon, -0.12, Form("%02d/07", m));
    576         txt.DrawText(mont, -0.10, Form("%d", m));
    577         line.DrawLine(monl, -0.12, monl, 0.12);
    578     }
    579 
    580     Double_t y6 = MTime(2006,1,1,0).GetAxisTime();
    581     Double_t y7 = MTime(2007,1,1,0).GetAxisTime();
    582     Double_t y8 = MTime(2008,1,1,0).GetAxisTime();
    583     Double_t y9 = MTime(2009,1,1,0).GetAxisTime();
    584     Double_t y0 = MTime(2010,1,1,0).GetAxisTime();
    585 //    Double_t y0 = MTime(2009,5,1,0).GetAxisTime();
    586 
    587     txt.SetTextSize(0.042);
    588     txt.DrawText(y6-(y7-y6)/2, -0.6, "2005");
    589     txt.DrawText((y6+y7)/2,    -0.6, "2006");
    590     txt.DrawText((y7+y8)/2,    -0.6, "2007");
    591     txt.DrawText((y8+y9)/2,    -0.6, "2008");
    592     txt.DrawText((y9+y0)/2,    -0.6, "2009");
    593 
    594     line.DrawLine(y6, -0.7, y6, 0.26);
    595     line.DrawLine(y7, -0.7, y7, 0.26);
    596     line.DrawLine(y8, -0.7, y8, 0.26);
    597     line.DrawLine(y9, -0.7, y9, 0.26);
    598     line.SetLineStyle(3);
    599     line.DrawLine(y6, 0.26, y6, 1.05*histres[2].GetMaximum());
    600     line.DrawLine(y7, 0.26, y7, 1.05*histres[2].GetMaximum());
    601     line.DrawLine(y8, 0.26, y8, 1.05*histres[2].GetMaximum());
    602     line.DrawLine(y9, 0.26, y9, 1.05*histres[2].GetMaximum());
    603 
    604     line.SetLineColor(kBlue);
    605     line.SetLineWidth(2);
    606     line.SetLineStyle(kSolid);
    607     line.DrawLine(time[0], 1*360/16384.*60, time[bins.GetNumBins()], 1*360/16384.*60);
    608     line.SetLineStyle(9);
    609     line.DrawLine(time[0], 2*360/16384.*60, time[bins.GetNumBins()], 2*360/16384.*60);
    610     line.SetLineStyle(7);
    611     line.DrawLine(time[0], 3*360/16384.*60, time[bins.GetNumBins()], 3*360/16384.*60);
    612 
    613     /*    result.SetMinimum(-0.06);
    614     result.SetMarkerStyle(kFullDotMedium);
    615     result.DrawClone("LP");*/
     85    plot(Next);
    61686}

  • trunk/Cosy/tpoint/plot_m2.C

    r9512 r9953  
    1 #include <iomanip>
     1#include "mars/MTime.h"
    22
    3 void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
    4     {
    5         TView *view = pad->GetView();
    6 
    7         if (!view)
    8         {
    9             cout << "No View!" << endl;
    10             return;
    11         }
    12 
    13         TMarker mark0;
    14         TMarker mark1;
    15         mark0.SetMarkerStyle(kStar);
    16         mark1.SetMarkerStyle(kStar);
    17         mark1.SetMarkerColor(kRed);
    18    
    19         r0 /= 90;
    20         r1 /= 90;
    21         phi0 *= TMath::DegToRad();
    22         phi1 *= TMath::DegToRad();
    23    
    24         Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
    25         Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};
    26 
    27         mark0.DrawMarker(x0[0], x0[1]);
    28         mark1.DrawMarker(x1[0], x1[1]);
    29 
    30         return;
    31         Double_t y0[3], y1[3];
    32    
    33         view->WCtoNDC(x0, y0);
    34         view->WCtoNDC(x1, y1);
    35    
    36         mark0.DrawMarker(y0[0], y0[1]);
    37         mark1.DrawMarker(y1[0], y1[1]);
    38     }
    39 
    40 int fill(const char *fname, TGraph *g, TH1 *h)
    41 {
    42 /*
    43     TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360,  10, 0, 90);
    44     h2res1.SetBit(TH1::kNoStats);
    45     h2res1.DrawCopy("surf1pol");
    46     gPad->Modified();
    47     gPad->Update();
    48     gPad->SetTheta(90);
    49     gPad->SetPhi(-90);
    50 
    51     DrawMarker(gPad, 45, 0, 0, 0);
    52     gPad->Modified();
    53     gPad->Update();
    54 
    55     return;
    56 
    57   */
    58     ifstream fin(fname);
    59 
    60     cout << "Reading " << setw(23) << fname << "..." << flush;
    61 
    62     while (1)
    63     {
    64         TString str;
    65         str.ReadLine(fin);
    66         if (!fin)
    67             break;
    68 
    69         if (str.Contains("#"))
    70             continue;
    71 
    72         Float_t alt, az, dalt, daz, mjd;
    73         sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
    74                &az, &alt, &dalt, &daz, &mjd);
    75 
    76         if (dalt==0/* ||  GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
    77             continue;
    78 
    79         mjd -=  53140.097505;
    80 
    81         Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz);
    82 
    83         g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
    84         g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
    85         g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);
    86 
    87         h->Fill(res);
    88     }
    89 
    90     cout << "done (" << setw(3) << (Int_t)h->GetEntries() << "/";
    91     cout << setw(4) << g[0].GetN() << ") " << flush;
    92 
    93     return g[0].GetN();
    94 }
    95 
    96 struct Description_t
    97 {
    98     const char *fName;
    99     const char *fTitle;
    100     const char *fFile;
     3Double_t dates[] = {
     4    MTime(2009,  3, 30).GetAxisTime(), // Start of M2 observations
     5    MTime(2009,  4, 11).GetAxisTime(), // New pointing model
     6    MTime(2009,  5, 11).GetAxisTime(), // New pointing model
     7    MTime(2009,  6, 11).GetAxisTime(), // Something happened (readon unknown)
     8    MTime(2009,  8, 17).GetAxisTime(), // New pointing model
     9    MTime(2009,  9,  9).GetAxisTime(), // Something happened (reason unknown)
     10    MTime(2010,  2,  1).GetAxisTime(), // New pointing model
     11    // MTime(2010,  2,  3).GetAxisTime(), // New LUTs
     12    MTime(2010,  2, 26).GetAxisTime(), // New pointing model
     13    MTime(2010,  3, 31).GetAxisTime(), // New pointing model
     14    MTime(-1).GetAxisTime(),
     15    -1
    10116};
    10217
    103 const Int_t counts = 46+15;
    104 Description_t desc[counts] =
    105 {
    106     {"090401",  "TPoints Residuals 8/2004-2" , "tpoint/m2/first/tpoints_m2_1.txt"},
    107     {"+090401",  "TPoints Residuals 9/2004" ,  "tpoint/m2/first/tpoints_m2_2.txt"},
    108     {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090403_025841.txt"},
    109     {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090403_214619.txt"},
    110     {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090404_235955.txt"},
    111     {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090405_015920.txt"},
    112     // New pointing model
    113     {"090402",   "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090411_231731.txt"},
    114     {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090412_232919.txt"},
    115     {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090414_001746.txt"},
    116     {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090418_012742.txt"},
    117     {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090418_024412.txt"},
    118     // New pointing model
    119     {"090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_12/tpoint_20090512_011933.txt"},
    120     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_12/tpoint_20090512_013112.txt"},
    121     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090512_210644.txt"},
    122     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_002758.txt"},
    123     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_025123.txt"},
    124     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_033504.txt"},
    125     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_14/tpoint_20090514_013332.txt"},
    126     // Just start of a new observation period
    127     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_01/tpoint_20090531_215139.txt"},
    128     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_01/tpoint_20090531_222548.txt"},
    129     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_02/tpoint_20090601_223001.txt"},
    130     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090602_213508.txt"},
    131     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090602_230944.txt"},
    132     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090603_000616.txt"},
    133     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090603_011935.txt"},
    134     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_04/tpoint_20090603_215841.txt"},
    135     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_05/tpoint_20090604_215941.txt"},
    136 
    137     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_10/tpoint_20090609_232320.txt"},
    138     {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_11/tpoint_20090611_011148.txt"},
    139 
    140     // Something happened (reason unknown)
    141 
    142     {"090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_11/tpoint_20090611_023623.txt"},
    143     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_13/tpoint_20090613_033839.txt"},
    144     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_14/tpoint_20090614_021304.txt"},
    145     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_18/tpoint_20090618_041436.txt"},
    146     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_12/tpoint_20090712_014259.txt"},
    147     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_12/tpoint_20090712_024706.txt"},
    148     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_13/tpoint_20090713_004244.txt"},
    149     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_13/tpoint_20090713_025925.txt"},
    150     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_14/tpoint_20090714_022322.txt"},
    151     {"+090611",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_07_15/tpoint_20090715_021316.txt"},
    152     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_02/tpoint_20090801_225305.txt"},
    153     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_03/tpoint_20090802_224437.txt"},
    154     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_04/tpoint_20090803_231853.txt"},
    155     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_11/tpoint_20090810_234604.txt"},
    156     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_11/tpoint_20090811_020844.txt"},
    157     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_12/tpoint_20090812_001850.txt"},
    158     {"+090611",  "TPoints Residuals 08/2009" , "tpoint/m2/2009_08_13/tpoint_20090813_030805.txt"},
    159 
    160     // ------ 09/08/17 new pointing model -------
    161 
    162     {"090817",  "TPoints" , "tpoint/m2/2009_09_01/tpoint_20090831_233035.txt"},
    163     {"+090817",  "TPoints" , "tpoint/m2/2009_09_02/tpoint_20090901_231921.txt"},
    164     {"+090817",  "TPoints" , "tpoint/m2/2009_09_08/tpoint_20090907_213405.txt"},
    165     {"+090817",  "TPoints" , "tpoint/m2/2009_09_09/tpoint_20090908_225029.txt"},
    166     {"+090817",  "TPoints" , "tpoint/m2/2009_09_09/tpoint_20090908_231713.txt"},
    167 
    168     // ------ 09/08/17 What happened? -------
    169 
    170     {"090910",  "TPoints" , "tpoint/m2/2009_09_10/tpoint_20090909_235225.txt"},
    171     {"+090910",  "TPoints" , "tpoint/m2/2009_09_10/tpoint_20090910_025048.txt"},
    172     {"+090910",  "TPoints" , "tpoint/m2/2009_09_12/tpoint_20090912_014815.txt"},
    173     {"+090910",  "TPoints" , "tpoint/m2/2009_09_17/tpoint_20090917_015406.txt"},
    174     {"+090910",  "TPoints" , "tpoint/m2/2009_09_18/tpoint_20090917_202447.txt"},
    175     {"+090910",  "TPoints" , "tpoint/m2/2009_09_23/tpoint_20090922_223447.txt"},
    176     {"+090910",  "TPoints" , "tpoint/m2/2009_09_24/tpoint_20090923_204442.txt"},
    177     {"+090910",  "TPoints" , "tpoint/m2/2009_09_25/tpoint_20090925_002702.txt"},
    178     {"+090910",  "TPoints" , "tpoint/m2/2009_09_26/tpoint_20090925_212533.txt"},
    179     {"+090910",  "TPoints" , "tpoint/m2/2009_09_27/tpoint_20090926_212331.txt"},
    180 };
     18#include "plot.C"
    18119
    18220void plot_m2()
    18321{
    184     TGraph g[3];
     22    MDirIter Next;
     23    Next.AddDirectory("tpoint/m2", "tpoint*.txt", -1);
    18524
    186     MBinning bins(100, 0, 0.2);
    187 
    188     TArrayI n(counts);
    189 
    190     TH1F hx[counts];
    191     Int_t num = -1;
    192     for (int i=0; i<counts; i++)
    193     {
    194         if (desc[i].fName[0]!='+')
    195         {
    196             num++;
    197 
    198             hx[num].SetNameTitle(desc[i].fName, desc[i].fTitle);
    199             hx[num].SetDirectory(0);
    200             bins.Apply(hx[num]);
    201         }
    202 
    203         cout << setw(2) << num << ": " << flush;
    204         n[num] = fill(desc[i].fFile, g, &hx[num]);
    205         cout << " Mean: " << setw(5) << setprecision(2) << hx[num].GetMean() << " deg +/- " <<  hx[num].GetRMS()<< endl;
    206     }
    207 
    208     n.Set(++num);
    209 
    210     g[0].SetMarkerColor(kGreen);
    211     g[1].SetMarkerColor(kMagenta);
    212     g[2].SetMarkerColor(kBlack);
    213     //g[2].SetLineColor(kBlack);
    214     g[0].SetMarkerStyle(kFullDotMedium);
    215     g[1].SetMarkerStyle(kFullDotMedium);
    216     g[2].SetMarkerStyle(kFullDotLarge);
    217 
    218     // --------- First Canvas ----------
    219 
    220     new TCanvas("Vs Time", "");
    221 
    222     TObject *obj[4];
    223 
    224     for (int i=0; i<3; i++)
    225     {
    226         g[i].SetFillColor(kWhite);
    227         g[i].SetLineColor(kWhite);
    228         obj[i] = g[i].Clone();
    229         obj[i]->SetBit(kCanDelete);
    230         obj[i]->Draw(i==0?"AP":"P");
    231     }
    232 
    233     TLegend leg(0.905, 0.86, 0.99, 0.99);
    234     leg.AddEntry(obj[0], "  \\Delta\\theta");
    235     leg.AddEntry(obj[1], "  \\Delta\\phi");
    236     leg.AddEntry(obj[2], "  \\Delta");
    237     leg.DrawClone()->SetBit(kCanDelete);
    238 
    239     TLine l;
    240     l.SetLineColor(kBlue);
    241     for (int i=0; i<n.GetSize(); i++)
    242         l.DrawLine(n[i], 0, n[i], 0.2);
    243 
    244     // --------- Second Canvas ----------
    245 
    246     new TCanvas("Distrib", "");
    247 
    248     Double_t max=0;
    249     for (int i=0; i<n.GetSize(); i++)
    250     {
    251         if (hx[i].GetEntries()==0)
    252         {
    253             cout << "Skip #" << i << endl;
    254             continue;
    255         }
    256 
    257         hx[i].Scale(1./hx[i].GetEntries());
    258         max = TMath::Max(max, hx[i].GetMaximum());
    259     }
    260     for (int i=0; i<n.GetSize(); i++)
    261     {
    262         hx[i].SetMaximum(max*1.05);
    263         if (i<6)
    264             hx[i].SetLineColor(kRed+i);
    265         else
    266             hx[i].SetMarkerStyle(kPlus+i-6);
    267 
    268 
    269     }
    270 
    271     for (int i=0; i<counts; i++)
    272         hx[i].DrawCopy(i==0?"LP":"LPsame");
    273 
    274     Double_t time[] = {
    275         MTime(2009,  3, 30).GetAxisTime(),
    276         MTime(2009,  4, 11).GetAxisTime(),
    277         MTime(2009,  5, 11).GetAxisTime(),
    278         MTime(2009,  6, 11).GetAxisTime(),
    279         MTime(2009,  8, 17).GetAxisTime(),
    280         MTime(2009,  9,  9).GetAxisTime(),
    281         MTime(2009, 12, 31).GetAxisTime(),
    282     };
    283 
    284     TH1D histres[4];
    285 
    286     MBinning bins;
    287     bins.SetEdges(TArrayD(7, time));
    288 
    289     bins.Apply(histres[0]);
    290     bins.Apply(histres[1]);
    291     bins.Apply(histres[2]);
    292     bins.Apply(histres[3]);
    293 
    294     TGraphAsymmErrors result[4];
    295     TGraph resultm;
    296     for (int i=0; i<n.GetSize(); i++)
    297     {
    298         cout << i+1 << " - Mean:  " << Form("%.4f +- %.4f", hx[i].GetMean(), hx[i].GetRMS());
    299         cout << "   (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")   " << (int)hx[i].GetEntries() << endl;
    300 
    301         /*
    302          TF1 fg("fg", "gaus", 0, 1);
    303          hx[i].Fit(&fg);
    304          result2.SetPoint(result.GetN(), result.GetN()+7, hx[i].GetMean());
    305          result2.SetPointError(result.GetN()-1, 0, hx[i].GetRMS()/sqrt(hx[i].GetEntries()));
    306          result.SetPoint(result.GetN(), result.GetN()+7, fg.GetParameter(1));
    307          result.SetPointError(result.GetN()-1, 0, fg.GetParameter(2));
    308          */
    309 
    310         //Double_t q[4] = { MMath::GaussProb(0.5), MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
    311         Double_t q[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
    312         Double_t rc[4];
    313         hx[i].GetQuantiles(4, rc, q);
    314 
    315         for (int j=0; j<4; j++)
    316         {
    317             histres[j].SetBinContent(i+1, 60*rc[j]);
    318 
    319             result[j].SetPoint(i, time[i], 60*rc[j]);
    320             result[j].SetPointError(i, 0, 0, 60*rc[j], 0);
    321 
    322             //result[j].SetPoint(result[j].GetN(), result[j].GetN()+1, 60*rc[j]);
    323             //result[j].SetPointError(result[j].GetN()-1, 0, 0, 60*rc[j], 0);
    324         }
    325 
    326 //        result.SetPoint(result.GetN(), result.GetN()+1, rc[1]);
    327 //        result2.SetPoint(result2.GetN(), result2.GetN()+1, rc[2]);
    328 
    329         //result.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
    330         //result2.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
    331 
    332 //        result.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
    333 //        result2.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
    334 
    335 //        result2.SetPointError(result.GetN()-1, 0, 0, rc[2], 0);
    336 
    337 //        resultm.SetPoint(resultm.GetN(), resultm.GetN()+1, 60*hx[i].GetMean());
    338         resultm.SetPoint(resultm.GetN(), (time[i]+time[i+1])/2, 60*hx[i].GetMean());
    339 
    340     }
    341 
    342     new TCanvas;
    343 
    344     gPad->SetBorderMode(0);
    345     gPad->SetFrameBorderMode(0);
    346     gPad->SetFillColor(kWhite);
    347     gPad->SetRightMargin(0.01);
    348     gPad->SetTopMargin(0.01);
    349     gPad->SetLeftMargin(0.06);
    350     gPad->SetGridy();
    351 
    352     //Int_t col[] = { 12, 15, 17, 19 };
    353     //Int_t col[] = { 12, 16, 18, 0 };
    354     Int_t col[] = { 13, 16, 19, 0 };
    355 
    356     TH1 *h = &histres[2];//result[3].GetHistogram();
    357     h->SetXTitle("");
    358     h->SetYTitle("Residual / arcmin");
    359     h->SetBit(TH1::kNoStats);
    360     h->GetXaxis()->CenterTitle();
    361     h->GetYaxis()->CenterTitle();
    362     h->GetYaxis()->SetTitleOffset(0.75);
    363 //    h->GetXaxis()->SetTimeFormat("%m/%y %F1995-01-01 00:00:00 GMT");
    364 //    h->GetXaxis()->SetTimeDisplay(1);
    365 
    366     h->GetXaxis()->SetLabelColor(kWhite);
    367 
    368     TLine line;
    369 
    370     for (int j=2; j>=0; j--)
    371     {
    372         histres[j].SetMinimum(0);
    373         histres[j].SetFillColor(col[j]);//12+2*j);
    374 
    375         histres[j].DrawCopy(j==2?"":"same");
    376 
    377         /*
    378         //result[j].SetLineColor(kBlue);
    379         //result[j].SetLineWidth(2);
    380         //result[j].SetMarkerColor(kBlue);
    381         result[j].SetMinimum(0);
    382         //result2.SetMarkerStyle(kFullDotMedium);
    383         //result[j].SetMarkerStyle(23);
    384         result[j].SetFillColor(col[j]);//12+2*j);
    385         result[j].DrawClone(j==3 ? "ABX" : "B"); // E3 B
    386         */
    387     }
    388 
    389 
    390     resultm.SetMarkerStyle(20);
    391     resultm.SetMarkerSize(0.8);
    392     resultm.DrawClone("P");
    393 
    394     line.DrawLine(time[0], 0, time[0], histres[2].GetMaximum());
    395     for (int i=0; i<bins.GetNumBins()-1; i++)
    396         line.DrawLine(time[i+1], 0, time[i+1], histres[2].GetBinContent(i+1));
    397 
    398     TText txt;
    399     txt.SetTextSize(0.037);
    400 //    txt.SetTextAngle(-45);
    401     txt.SetTextAlign(23);
    402     /*
    403     for (int m=4; m<13; m++)
    404     {
    405         Double_t monl = MTime(2005, m,   1,0).GetAxisTime();
    406         Double_t mont = MTime(2005, m,   15,0).GetAxisTime();
    407 //        txt.DrawText(mon, -0.12, Form("%02d/05", m));
    408         txt.DrawText(mont, -0.10, Form("%d", m));
    409         line.DrawLine(monl, -0.12, monl, 0.12);
    410     }
    411     for (int m=1; m<13; m++)
    412     {
    413         Double_t monl = MTime(2006, m,   1,0).GetAxisTime();
    414         Double_t mont = MTime(2006, m,   15,0).GetAxisTime();
    415 //        txt.DrawText(mon, -0.12, Form("%02d/06", m));
    416         txt.DrawText(mont, -0.10, Form("%d", m));
    417         line.DrawLine(monl, -0.12, monl, 0.12);
    418     }
    419     for (int m=1; m<13; m++)
    420     {
    421         Double_t monl = MTime(2007, m,   1,0).GetAxisTime();
    422         Double_t mont = MTime(2007, m,   15,0).GetAxisTime();
    423 //        txt.DrawText(mon, -0.12, Form("%02d/07", m));
    424         txt.DrawText(mont, -0.10, Form("%d", m));
    425         line.DrawLine(monl, -0.12, monl, 0.12);
    426     }*/
    427     for (int m=4; m<13; m++)
    428     {
    429         Double_t monl = MTime(2009, m,   1,0).GetAxisTime();
    430         Double_t mont = MTime(2009, m,   15,0).GetAxisTime();
    431 //        txt.DrawText(mon, -0.12, Form("%02d/07", m));
    432         txt.DrawText(mont, -0.10, Form("%d", m));
    433         line.DrawLine(monl, -0.12, monl, 0.12);
    434     }
    435 
    436 //    Double_t y6 = MTime(2006,1,1,0).GetAxisTime();
    437 //    Double_t y7 = MTime(2007,1,1,0).GetAxisTime();
    438 //    Double_t y8 = MTime(2008,1,1,0).GetAxisTime();
    439     Double_t y9 = MTime(2009,1,1,0).GetAxisTime();
    440     Double_t y0 = MTime(2010,1,1,0).GetAxisTime();
    441 
    442     txt.SetTextSize(0.042);
    443 //    txt.DrawText(y6-(y7-y6)/2, -0.6, "2005");
    444 //    txt.DrawText((y6+y7)/2,    -0.6, "2006");
    445 //    txt.DrawText((y7+y8)/2,    -0.6, "2007");
    446     txt.DrawText((y9+y0)/2,    -0.6, "2009");
    447 
    448 //    line.DrawLine(y6, -0.7, y6, 0.26);
    449 //    line.DrawLine(y7, -0.7, y7, 0.26);
    450 //    line.DrawLine(y8, -0.7, y8, 0.26);
    451     line.DrawLine(y9, -0.7, y9, 0.26);
    452     line.SetLineStyle(3);
    453 //    line.DrawLine(y6, 0.26, y6, 1.05*histres[2].GetMaximum());
    454 //    line.DrawLine(y7, 0.26, y7, 1.05*histres[2].GetMaximum());
    455 //    line.DrawLine(y8, 0.26, y8, 1.05*histres[2].GetMaximum());
    456     line.DrawLine(y9, 0.26, y9, 1.05*histres[2].GetMaximum());
    457 
    458     line.SetLineColor(kBlue);
    459     line.SetLineWidth(2);
    460     line.SetLineStyle(kSolid);
    461     line.DrawLine(time[0], 1*360/16384.*60, time[bins.GetNumBins()], 1*360/16384.*60);
    462     line.SetLineStyle(9);
    463     line.DrawLine(time[0], 2*360/16384.*60, time[bins.GetNumBins()], 2*360/16384.*60);
    464     line.SetLineStyle(7);
    465     line.DrawLine(time[0], 3*360/16384.*60, time[bins.GetNumBins()], 3*360/16384.*60);
    466 
    467     /*    result.SetMinimum(-0.06);
    468     result.SetMarkerStyle(kFullDotMedium);
    469     result.DrawClone("LP");*/
     25    plot(Next);
    47026}
Note: See TracChangeset for help on using the changeset viewer.