Changeset 17521 for trunk/Mars/fact


Ignore:
Timestamp:
01/20/14 17:42:35 (11 years ago)
Author:
Daniela Dorner
Message:
new version to run automatically every morning
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/fact/plots/quality.C

    r15244 r17521  
    4747
    4848    TGraph g;
     49    g.SetName("Threshold");
    4950
    5051    while (file.GetNextRow())
     
    5455    }
    5556
     57    g.SetMinimum(281);
    5658    g.SetMarkerStyle(kFullDotMedium);
    5759    g.GetXaxis()->SetTimeDisplay(true);
     
    7072#include <pair>
    7173
    72 vector<pair<double, ln_equ_posn>> vec;
    73 
    74 ln_equ_posn FindPointing(Double_t time)
    75 {
    76     for (int i=0; i<vec.size(); i++)
    77         if (time<vec[i].first)
    78             return i==0 ? ln_equ_posn() : vec[i-1].second;
    79 
    80     return vec[vec.size()-1].second;
    81 }
    82 
     74vector<pair<double, Nova::EquPosn>> vecp;
     75
     76Nova::EquPosn FindPointing(Double_t time)
     77{
     78    for (int i=0; i<vecp.size(); i++)
     79        if (time<vecp[i].first)
     80        {
     81            if (i==0)
     82                return Nova::EquPosn();
     83            else
     84                return vecp[i-1].second;
     85        }
     86
     87    return vecp[vecp.size()-1].second;
     88}
     89/*
    8390Float_t Prediction(Double_t time)
    8491{
     
    93100
    94101    // Distance between source and moon
    95     //const double angle =
    96     //    MAstro::AngularDistance(TMath::DegToRad()*(90-moon.dec),
    97     //                            TMath::DegToRad()*moon.ra,
    98     //                            TMath::DegToRad()*(90-pos.dec),
    99     //                            TMath::DegToRad()*pos.ra);
     102    //const double angle = Nova::GetAngularSeparation(moon, hrzm)*TMath::DegToRad();
    100103
    101104    // Current prediction
     
    105108    double disk = Nova::GetLunarDisk(jd);
    106109
    107     double lc = /*cang==0 ? 0 :*/ calt*pow(disk, 2.5);///sqrt(cang);
    108 
    109     return 5+103*lc>4.5 ? 5+103*lc : 4.5;
     110    //                                                                           semi-major axis
     111    double lc = calt*pow(disk, 2.2)*pow(Nova::GetLunarEarthDist(jd)/384400, -2);///sqrt(cang);
     112
     113    return lc>0 ? 4+103*lc : 4;
     114}*/
     115
     116Float_t Prediction(Double_t time)
     117{
     118    Double_t jd = time + 40587 + 2400000.5;
     119    /*
     120    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
     121
     122    Nova::EquPosn pos = FindPointing(time);
     123
     124    // get local position of moon
     125    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
     126
     127    // Distance between source and moon
     128    const double angle = Nova::GetAngularSeparation(moon, pos);
     129
     130    // Distance between earth and moon relative to major semi-axis
     131    const double dist  = Nova::GetLunarEarthDist(jd)/384400;
     132
     133    // Current prediction
     134    double cang = 1-sin(angle*TMath::DegToRad());
     135    double calt = sin(hrzm.alt*TMath::DegToRad());
     136
     137    double disk = Nova::GetLunarDisk(jd);
     138
     139    // light condition
     140    double lc = sqrt(calt)*pow(disk, 2.3)*pow(cang+1, 1)*pow(dist, -2);
     141
     142    return lc>0 ? 7.2 + 69*lc : 7.2;
     143    */
     144
     145    // Sun properties
     146    Nova::EquPosn  sun  = Nova::GetSolarEquCoords(jd);
     147    Nova::ZdAzPosn hrzs = Nova::GetHrzFromEqu(sun, jd);
     148
     149    // Get source position
     150    Nova::EquPosn pos = FindPointing(time);
     151
     152    // Moon properties
     153    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
     154    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
     155    double        disk = Nova::GetLunarDisk(jd);
     156
     157    // Derived moon properties
     158    double angle = Nova::GetAngularSeparation(moon, pos);
     159    double edist = Nova::GetLunarEarthDist(jd)/384400;
     160
     161    // Current prediction
     162    double sin_malt  = hrzm.alt<0 ? 0 : sin(hrzm.alt*TMath::DegToRad());
     163    double sin_mdist = sin(angle*TMath::DegToRad());
     164    double cos_salt  = cos(hrzs.zd*TMath::DegToRad());
     165
     166    double c0 = pow(disk,     2.52);
     167    double c1 = pow(sin_malt, 0.72);
     168    double c2 = pow(edist,   -2.00);
     169    double c3 = exp(1.46*(1-sin_mdist)*(1-sin_mdist));
     170    double c4 = exp(33.0)*exp(-20.1*(1-cos_salt)*(1-cos_salt));
     171
     172    double cur = 6.4 + 96.9*c0*c1*c2*c3 + c4;
     173
     174   // cout << cur << " " << hrzm.alt << " " << c1 << endl;
     175
     176    return cur;
     177
    110178}
    111179
     
    131199    while (file.GetNextRow())
    132200    {
    133         ln_equ_posn p;
     201        Nova::EquPosn p;
    134202        p.ra  = ra*15;
    135203        p.dec = dec;
    136204
    137         vec.push_back(make_pair(time, p));
     205        vecp.push_back(make_pair(time, p));
    138206    }
    139207
     
    169237    TGraph g1;
    170238    TGraph g2;
     239    g1.SetName("Currents");
     240    g2.SetName("Prediction");
    171241
    172242    while (file.GetNextRow())
     
    234304
    235305    TGraph g1, g2;
     306    g1.SetName("TriggerRate");
     307    g2.SetName("RelOnTime");
    236308
    237309    while (file.GetNextRow())
     
    246318
    247319    g1.SetMinimum(0);
    248     g1.SetMaximum(200);
     320    g1.SetMaximum(269);
    249321    g1.SetMarkerStyle(kFullDotMedium);
    250322    g1.GetXaxis()->SetTimeDisplay(true);
     
    260332    gPad->GetCanvas()->cd(4);
    261333
    262     gPad->SetGridx();
     334    gPad->SetGrid();
    263335    gPad->SetTopMargin(0);
    264336    gPad->SetBottomMargin(0);
     
    278350    g2.DrawClone("AP");
    279351
    280 
    281     return 0;
    282 }
     352    return 0;
     353}
     354
     355void PlotRateQC(UInt_t night, MSQLServer &serv)
     356{
     357    TString query =
     358        "LEFT JOIN AnalysisResultsRunLP USING(fNight, fRunID) "
     359        "WHERE fRunTypeKey=1 AND NOT ISNULL (AnalysisResultsRunLP.fNumEvtsAfterCleaning) AND fNight=";
     360    query += night;
     361
     362    TTree *t = serv.GetTree("RunInfo", query);
     363    if (!t)
     364        return;
     365
     366    int save = gErrorIgnoreLevel;
     367    gErrorIgnoreLevel = kFatal;
     368
     369    gROOT->SetSelectedPad(0);
     370    gPad->GetCanvas()->cd(3);
     371
     372    t->Draw("AnalysisResultsRunLP.fNumEvtsAfterCleaning/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same");
     373    TGraph *g = (TGraph*)gPad->GetPrimitive("Graph");
     374    if (g)
     375    {
     376        g->SetName("CleaningRate");
     377        g->SetMarkerColor(kRed);
     378        g->SetMarkerStyle(29);//kFullDotMedium);
     379    }
     380
     381    t->Draw("AnalysisResultsRunLP.fNumEvtsAfterQualCuts/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same");
     382    g = (TGraph*)gPad->GetPrimitive("Graph");
     383    if (g)
     384    {
     385        g->SetName("RateAfterQC");
     386        g->SetMarkerColor(kBlue);
     387        g->SetMarkerStyle(29);//kFullDotMedium);
     388    }
     389
     390    gErrorIgnoreLevel = save;
     391}
     392
    283393
    284394Int_t PlotPointing(TArrayD **vec, TString fname)
     
    304414
    305415    TGraph g;
     416    g.SetName("Zd");
    306417
    307418    while (file.GetNextRow())
     
    309420            g.SetPoint(g.GetN(), time*24*3600, 90-zd);
    310421
    311     g.SetMinimum(0);
     422    g.SetMinimum(1);
    312423    g.SetMaximum(90);
    313424    g.SetMarkerStyle(kFullDotMedium);
    314425    g.GetXaxis()->SetTimeDisplay(true);
    315426    g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    316     g.GetXaxis()->SetLabelSize(0.1);
     427    g.GetXaxis()->SetLabelSize(0.12);
    317428    g.GetYaxis()->SetLabelSize(0.1);
    318429    g.GetYaxis()->SetTitle("ELEVATION");
     
    323434    return 0;
    324435}
    325 /*
    326 Int_t PlotMem(TArrayD **vec, TString fname)
    327 {
    328     fname += ".FAD_CONTROL_STATISTICS1.fits";
     436
     437Int_t PlotTemperature1(TArrayD **vec, TString fname)
     438{
     439    fname += ".TEMPERATURE_DATA.fits";
    329440
    330441    fits file(fname.Data());
     
    338449
    339450    Double_t time;
    340     UInt_t buf[5];
    341     ULong64_t mem[4];
    342     UInt_t rate[2];
     451    Float_t  temp;
    343452
    344453    if (!file.SetPtrAddress("Time",  &time))
    345454        return -1;
    346     if (!file.SetPtrAddress("bufferInfo", buf) ||
    347         !file.SetPtrAddress("memInfo",    mem) ||
    348         !file.SetPtrAddress("rateNew",    rate))
     455    if (!file.SetPtrAddress("T", &temp))
    349456        return -1;
    350457
    351458    TGraph g;
     459    g.SetName("ContainerTemp");
    352460
    353461    while (file.GetNextRow())
    354462        if (Contains(vec, time))
    355             g.SetPoint(g.GetN(), time*24*3600, mem[1]);
    356 
     463            g.SetPoint(g.GetN(), time*24*3600, temp);
     464
     465    g.SetMinimum(-5);
     466    g.SetMaximum(49);
    357467    g.SetMarkerStyle(kFullDotMedium);
     468    g.SetMarkerColor(kRed);
    358469    g.GetXaxis()->SetTimeDisplay(true);
    359     g.GetXaxis()->SetTimeFormat("%H:%M");
    360     g.GetXaxis()->SetTimeOffset(0, "gmt");
    361     g.GetXaxis()->SetLabelSize(0.12);
     470    g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
     471    g.GetXaxis()->SetLabelSize(0.1);
    362472    g.GetYaxis()->SetLabelSize(0.1);
    363     //g.GetYaxis()->SetTitle("ELEVATION");
     473    g.GetYaxis()->SetTitle("TEMP");
    364474    g.GetYaxis()->SetTitleOffset(0.2);
    365475    g.GetYaxis()->SetTitleSize(0.1);
     
    368478    return 0;
    369479}
    370 */
    371 void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0)
     480
     481Int_t PlotTemperature2(TArrayD **vec, TString fname)
     482{
     483    fname += ".MAGIC_WEATHER_DATA.fits";
     484
     485    fits file(fname.Data());
     486    if (!file)
     487    {
     488        cerr << fname << ": " << gSystem->GetError() << endl;
     489        return -2;
     490    }
     491
     492    //cout << fname << endl;
     493
     494    Double_t time;
     495    Float_t  temp;
     496
     497    if (!file.SetPtrAddress("Time",  &time))
     498        return -1;
     499    if (!file.SetPtrAddress("T", &temp))
     500        return -1;
     501
     502    TGraph g;
     503    g.SetName("OutsideTemp");
     504
     505    while (file.GetNextRow())
     506        if (Contains(vec, time))
     507            g.SetPoint(g.GetN(), time*24*3600, temp);
     508
     509    g.SetMarkerStyle(kFullDotMedium);
     510    g.SetMarkerColor(kBlue);
     511    g.DrawClone("P");
     512
     513    return 0;
     514}
     515
     516Int_t PlotTemperature3(TArrayD **vec, TString fname)
     517{
     518    fname += ".FSC_CONTROL_TEMPERATURE.fits";
     519
     520    fits file(fname.Data());
     521    if (!file)
     522    {
     523        cerr << fname << ": " << gSystem->GetError() << endl;
     524        return -2;
     525    }
     526
     527    //cout << fname << endl;
     528
     529    Double_t time;
     530    Float_t  temp[31];
     531
     532    if (!file.SetPtrAddress("Time",  &time))
     533        return -1;
     534    if (!file.SetPtrAddress("T_sens", temp))
     535        return -1;
     536
     537    TGraph g, g1, g2;
     538    g.SetName("SensorTempAvg");
     539    g1.SetName("SensorTempMin");
     540    g2.SetName("SensorTempMax");
     541
     542    while (file.GetNextRow())
     543        if (Contains(vec, time))
     544        {
     545            float min =  100;
     546            float max = -100;
     547            double avg = 0;
     548            int num = 0;
     549            for (int i=0; i<31; i++)
     550                if (temp[i]!=0)
     551                {
     552                    avg += temp[i];
     553                    num++;
     554
     555                    min = TMath::Min(min, temp[i]);
     556                    max = TMath::Max(max, temp[i]);
     557                }
     558
     559            g.SetPoint(g.GetN(), time*24*3600, avg/num);
     560            g1.SetPoint(g1.GetN(), time*24*3600, min);
     561            g2.SetPoint(g2.GetN(), time*24*3600, max);
     562        }
     563
     564    g.SetMarkerStyle(kFullDotMedium);
     565    g.DrawClone("P");
     566
     567    /*
     568     g1.SetLineWidth(1);
     569     g1.DrawClone("L");
     570
     571     g2.SetLineWidth(1);
     572     g2.DrawClone("L");
     573    */
     574    return 0;
     575}
     576
     577Int_t PlotTemperature4(TArrayD **vec, TString fname)
     578{
     579    fname += ".FAD_CONTROL_TEMPERATURE.fits";
     580
     581    fits file(fname.Data());
     582    if (!file)
     583    {
     584        cerr << fname << ": " << gSystem->GetError() << endl;
     585        return -2;
     586    }
     587
     588    //cout << fname << endl;
     589
     590    Double_t time;
     591    Float_t  temp[160];
     592
     593    if (!file.SetPtrAddress("Time",  &time))
     594        return -1;
     595    if (!file.SetPtrAddress("Data1", temp) &&
     596        !file.SetPtrAddress("temp", temp))
     597        return -1;
     598
     599    Int_t num = file.GetN("temp")==0 ? file.GetN("Data1") : file.GetN("temp");
     600    Int_t beg = num==82 ? 2 : 0;
     601
     602    TGraphErrors g1;
     603    TGraph g2,g3;
     604
     605    g1.SetName("FadTempAvg");
     606    g2.SetName("FadTempMin");
     607    g3.SetName("FadTempMax");
     608
     609    while (file.GetNextRow())
     610        if (Contains(vec, time))
     611        {
     612            double avg = 0;
     613            double rms = 0;
     614            float min =  100;
     615            float max = -100;
     616            for (int i=beg; i<num; i++)
     617            {
     618                avg += temp[i];
     619                rms += temp[i]*temp[i];
     620
     621                min = TMath::Min(min, temp[i]);
     622                max = TMath::Max(max, temp[i]);
     623            }
     624
     625            avg /= num-beg;
     626            rms /= num-beg;
     627
     628            g1.SetPoint(g1.GetN(), time*24*3600, avg);
     629            g1.SetPointError(g1.GetN()-1, 0, sqrt(rms-avg*avg));
     630
     631            g2.SetPoint(g2.GetN(), time*24*3600, min);
     632            g3.SetPoint(g3.GetN(), time*24*3600, max);
     633        }
     634
     635    g1.SetLineColor(kGreen);
     636    g1.DrawClone("[]");
     637
     638    g2.SetLineColor(kGreen);
     639    g2.SetLineWidth(1);
     640    g2.DrawClone("L");
     641
     642    g3.SetLineColor(kGreen);
     643    g3.SetLineWidth(1);
     644    g3.DrawClone("L");
     645
     646    return 0;
     647}
     648
     649void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, TString outpath="./")
    372650{
    373651    // To get correct dates in the histogram you have to add
     
    377655    {
    378656        UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt();
    379         y = nt/10000;
     657        y =  nt/10000;
    380658        m = (nt/100)%100;
    381         d = nt%100;
     659        d =  nt%100;
    382660
    383661        cout << y << "/" << m << "/" << d << endl;
     
    409687    {
    410688        TString query;
    411         query += "SELECT fRunID, fRunStart, fRunStop FROM runinfo";
     689        query += "SELECT fRunID, fRunStart, fRunStop FROM RunInfo";
    412690        query += " WHERE fNight=";
    413691        query += night;
    414         query += " AND fRunTypeKEY=1 ORDER BY fRunStart";
     692        query += " AND fRunTypeKey=1 ORDER BY fRunStart";
    415693
    416694        TSQLResult *res = serv.Query(query);
     
    438716        if (n==0)
    439717            cout << "WARNING - No data runs in db, displaying all data." << endl;
     718        else
     719            cout << "Num: " << n << "\n" << endl;
    440720    }
    441721
    442722    TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960);
    443     c->Divide(1, 5, 1e-5, 1e-5);
     723    c->Divide(1, 6, 1e-5, 1e-5);
    444724
    445725    gROOT->SetSelectedPad(0);
     726
    446727    c->cd(1);
    447     gPad->SetGridx();
     728    gPad->SetGrid();
    448729    gPad->SetTopMargin(0);
    449730    gPad->SetRightMargin(0.001);
     
    454735    gROOT->SetSelectedPad(0);
    455736    c->cd(2);
    456     gPad->SetGridx();
     737    gPad->SetGrid();
    457738    gPad->SetTopMargin(0);
    458739    gPad->SetRightMargin(0.001);
     
    463744    gROOT->SetSelectedPad(0);
    464745    c->cd(3);
    465     gPad->SetGridx();
     746    gPad->SetGrid();
    466747    gPad->SetTopMargin(0);
    467748    gPad->SetBottomMargin(0);
     
    469750    gPad->SetLeftMargin(0.04);
    470751    cout << PlotRate(runs, fname) << endl;
     752    cout << PlotRateQC(night, serv) << endl;
    471753
    472754    gROOT->SetSelectedPad(0);
    473755    c->cd(5);
    474     gPad->SetGridx();
     756    gPad->SetGrid();
     757    gPad->SetTopMargin(0);
     758    gPad->SetBottomMargin(0);
     759    gPad->SetRightMargin(0.001);
     760    gPad->SetLeftMargin(0.04);
     761    cout << PlotPointing(runs, fname) << endl;
     762
     763    gROOT->SetSelectedPad(0);
     764    c->cd(6);
     765    gPad->SetGrid();
    475766    gPad->SetTopMargin(0);
    476767    gPad->SetRightMargin(0.001);
    477768    gPad->SetLeftMargin(0.04);
    478     cout << PlotPointing(runs, fname) << endl;
    479 /*
    480     gROOT->SetSelectedPad(0);
    481     c->cd(6);
    482     gPad->SetGridx();
    483     gPad->SetTopMargin(0);
    484     gPad->SetRightMargin(0.001);
    485     gPad->SetLeftMargin(0.04);
    486     cout << PlotMem(runs, fname) << endl;
    487 */
    488     c->SaveAs(Form("%04d%02d%02d-quality.pdf", y, m, d));
     769    cout << PlotTemperature1(runs, fname) << endl;
     770    cout << PlotTemperature2(runs, fname) << endl;
     771    cout << PlotTemperature3(runs, fname) << endl;
     772    cout << PlotTemperature4(runs, fname) << endl;
     773
     774    c->SaveAs(Form("%s/%04d%02d%02d-quality.png", outpath.Data(), y, m, d));
    489775}
    490776
Note: See TracChangeset for help on using the changeset viewer.