Changeset 18282


Ignore:
Timestamp:
08/18/15 10:28:44 (9 years ago)
Author:
Jens Buss
Message:
merged Trunk into Branch
Location:
branches/MarsGapdTimeJitter
Files:
34 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/MarsGapdTimeJitter

  • branches/MarsGapdTimeJitter/Makefile.rules

    r17980 r18282  
    3939        @echo " - Generating dictionary $(CINT)Cint.cc"
    4040        rootcint -f $(CINT)Cint.cc \
    41         -c $(INCLUDES) $(LIBNOVA_INCLUDE_PATH) $(DEFINES) $(HEADERS) $(CINT)Incl.h $(CINT)LinkDef.h
     41        -c -p $(INCLUDES) $(LIBNOVA_INCLUDE_PATH) $(DEFINES) $(HEADERS) $(CINT)Incl.h $(CINT)LinkDef.h
    4242
    4343%.d:

  • branches/MarsGapdTimeJitter/fact/analysis/callisto.C

    r17894 r18282  
    11#include "MLogManip.h"
    22
    3 int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output", bool use_delays=false)
     3int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output",
     4//             TString drstime = "/gpfs0/fact/processing/drs_time_calib/20111115.time.drs.fits",
     5             const char *drstime = "/gpfs0/fact/processing/drs_time_calib/20111115.time.drs.fits",
     6             const char *delays="resources/delays-20150217.txt")
    47{
    58    // ======================================================
    … …  
    107110    }
    108111
    109     TString timfile = drs.GetFileName(0, MSequence::kFitsDat);
    110112    TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs);
    111113    TString pedfile = seq.GetFileName(0, MSequence::kFitsPed);
    … …  
    115117    gLog << "DRS calib     300: " << drsfile << '\n';
    116118    gLog << "DRS calib    1024: " << drs1024 << "\n\n";
     119    gLog << "DRS calib    Time: " << drstime << "\n\n";
     120
     121    MDrsCalibrationTime drscalibtime;
     122    if (!drscalibtime.ReadFits(drstime)))
     123        return 5;
    117124
    118125    MDrsCalibration drscalib300;
    119126    if (!drscalib300.ReadFits(drsfile.Data()))
    120         return 5;
     127        return 6;
    121128
    122129    MDrsCalibration drscalib1024;
    123130    if (!drscalib1024.ReadFits(drs1024.Data()))
    124         return 6;
     131        return 7;
    125132
    126133    gLog << all;
    127     gLog << "Time calibration : " << timfile << '\n';
    128134    gLog << "Pedestal     file: " << pedfile << '\n';
    129135    gLog << "Light Pulser file: " << calfile << '\n' << endl;
    … …  
    135141    {
    136142        gLog << err << "ERROR - Sequence valid but without files." << endl;
    137         return 7;
     143        return 8;
    138144    }
    139145    iter.Print();
    … …  
    195201    // hrate.DefaultLabelY("ERROR");
    196202
    197     MDrsCalibrationTime timecam;
     203//    MDrsCalibrationTime timecam;
    198204
    199205    gStyle->SetOptFit(kTRUE);
    200206
    201207    // ======================================================
    202 
     208/*
    203209    gLog << endl;
    204210    gLog.Separator("Processing DRS timing calibration run");
    … …  
    240246    if (!loop0.GetDisplay())
    241247        return 9;
    242 
     248*/
    243249    /*
    244250     MHDrsCalibrationT *t = (MHDrsCalibrationT*)plist4.FindObject("MHDrsCalibrationT");
    … …  
    258264    plist1.AddToList(&drscalib300);
    259265    plist1.AddToList(&badpixels);
    260     plist1.AddToList(&timecam);
     266    plist1.AddToList(&drscalibtime);
    261267
    262268    MEvtLoop loop1("DetermineCalConst");
    … …  
    348354    tlist1.AddToList(&fill1d);
    349355
    350     if (!loop1.Eventloop(max1))
    351         return 10;
    352 
    353     if (!loop1.GetDisplay())
    354         return 11;
    355 
    356     if (use_delays)
    357         timecam.SetDelays(*evt1h.GetHist());
     356    //if (!loop1.Eventloop(max1))
     357    //    return 10;
     358
     359    //if (!loop1.GetDisplay())
     360    //    return 11;
     361
     362    if (delays)
     363    {
     364        TGraph g(delays);
     365        if (g.GetN()!=1440)
     366        {
     367            gLog << err << "Error reading file " << delays << endl;
     368            return 41;
     369        }
     370
     371        drscalibtime.SetDelays(g);
     372    }
    358373
    359374    // ======================================================
    … …  
    368383    plist3.AddToList(&drscalib300);
    369384    plist3.AddToList(&badpixels);
    370     plist3.AddToList(&timecam);
     385    plist3.AddToList(&drscalibtime);
    371386
    372387    MEvtLoop loop3("DetermineRndmPed");
    … …  
    442457    plist4.AddToList(&drscalib300);
    443458    plist4.AddToList(&badpixels);
    444     plist4.AddToList(&timecam);
     459    plist4.AddToList(&drscalibtime);
    445460
    446461    MEvtLoop loop4("DetermineExtractedPed");
    … …  
    510525    plist5.AddToList(&drscalib300);
    511526    plist5.AddToList(&badpixels);
    512     plist5.AddToList(&timecam);
     527    plist5.AddToList(&drscalibtime);
    513528
    514529    MEvtLoop loop5("CalibratingData");

  • branches/MarsGapdTimeJitter/fact/analysis/callisto_data.C

    r18078 r18282  
    9090    // -------------------------------------------------------
    9191
    92     TFile file(drstime);
    93     if (file.IsZombie())
     92    MDrsCalibrationTime timecam;
     93    if (!timecam.ReadFits(drstime))
    9494    {
    95         gLog << err << "ERROR - Could not open " << drstime << endl;
     95        gLog << err << "ERROR - Could not get MDrsCalibrationTime from " << drstime << endl;
    9696        return 21;
    9797    }
    9898
    99     MDrsCalibrationTime *timecam = 0;
    100     file.GetObject("MDrsCalibrationTime", timecam);
    101     if (!timecam)
     99    if (delays)
    102100    {
    103         gLog << err << "ERROR - Could not get MDrsCalibrationTime from " << drstime << endl;
    104         return 22;
     101        TGraph g(delays);
     102        if (g.GetN()!=1440)
     103        {
     104            gLog << err << "Error reading file " << delays << endl;
     105            return 22;
     106        }
     107
     108        timecam.SetDelays(g);
    105109    }
    106110
    … …  
    113117    // -------------------------------------------------------
    114118
    115     //if (use_delays)
    116     //    timecam.SetDelays(*evt1h.GetHist());
     119    if (delays)
     120    {
     121        TGraph g(delays);
     122        if (g.GetN()!=1440)
     123        {
     124            gLog << err << "Error reading file " << delays << endl;
     125            return 41;
     126        }
     127
     128        timecam.SetDelays(g);
     129    }
    117130
    118131    // -------------------------------------------------------
    … …  
    194207    plist5.AddToList(&drscalib300);
    195208    plist5.AddToList(&badpixels);
    196     plist5.AddToList(timecam);
     209    plist5.AddToList(&timecam);
    197210
    198211    MEvtLoop loop5("CalibratingData");

  • branches/MarsGapdTimeJitter/fact/plots/plotratescan.C

    r15243 r18282  
    129129{
    130130    TString query;
    131     query += "SELECT fTimeBegin, fTimeEnd, fVoltageIsOn, fOvervoltage, fCurrentMedMean, fNight ";
     131    query += "SELECT fTimeBegin, fTimeEnd, fVoltageIsOn, fOvervoltage, fCurrentMedMean, fNight, fAzMin, fAzMax, fZdMin, fZdMax ";
    132132    query += "FROM Ratescan WHERE fRatescanID=";
    133133    query += search_id;
    … …  
    149149    const char *time_end = (*row)[1];
    150150    const char *night    = (*row)[5];
     151    const char *az_beg = (*row)[6];
     152    const char *az_end = (*row)[7];
     153    const char *zd_beg = (*row)[8];
     154    const char *zd_end = (*row)[9];
    151155
    152156    int   voltage_on  = (*row)[2] ? atoi((*row)[2]) :   -1;
    … …  
    167171    leg.SetBorderSize(1);
    168172    leg.SetFillColor(kWhite);
    169     leg.AddText("Ratescan");
    170     leg.AddText("");
    171     leg.AddText(Form("Begin    %s", time_beg));
    172     leg.AddText(Form("End       %s", time_end));
    173     leg.AddText("");
     173    leg.SetTextAlign(12);
     174    leg.AddText(Form("Ratescan %d ", search_id));
     175    //leg.AddText("");
     176    leg.AddText(Form("Begin       %s", time_beg));
     177    leg.AddText(Form("End          %s", time_end));
     178    leg.AddText(Form("Az            %s#circ to %s#circ", az_beg, az_end));
     179    leg.AddText(Form("Zd             %s#circ to %s#circ", zd_beg, zd_end));
     180    //leg.AddText("");
    174181    if (voltage_on==0)
    175182        leg.AddText("Voltage off");
    … …  
    179186            leg.AddText(Form("Current     <I_{med}>  =  %.1f #muA", current));
    180187        if (overvoltage>-70)
    181             leg.AddText(Form("Voltage         #DeltaU  =  %+.2f V", overvoltage));
     188            leg.AddText(Form("Voltage   #DeltaU  =  %+.2f V", overvoltage));
    182189    }
    183190
    … …  
    209216    g.DrawClone("PL");
    210217
     218    TGraph gr("good_ratescan_edit.txt", "%lg %lg");
     219    gr.SetLineColor(12);
     220    gr.DrawClone("L");
     221
    211222    c->Write();
    212223
    … …  
    216227    name += time_beg;
    217228
    218     c->SaveAs(name+".pdf");
    219     c->SaveAs(name+".eps");
    220     c->SaveAs(name+".png");
     229    //c->SaveAs(name+".pdf");
     230    //c->SaveAs(name+".eps");
     231    //c->SaveAs("/loc_data/analysis/"+name+".png");
     232    c->SaveAs(Form("/loc_data/analysis/ratescans/%04d/%02d/%02d/%06d_%d.png", atoi(night)/10000, (atoi(night)/100)%100, atoi(night)%100, atoi(night), search_id));
    221233
    222234    delete c;
    … …  
    305317
    306318    TString oname = Form("%06d-ratescan.root", night);
    307 
    308     cout << "   " << oname << '\n' << endl;
    309 
     319    //cout << "   " << oname << '\n' << endl;
    310320    TFile rootfile(oname.Data(), "recreate");
    311321

  • branches/MarsGapdTimeJitter/fact/plots/quality.C

    r17956 r18282  
    9595    Double_t jd = time + 40587 + 2400000.5;
    9696
     97    // Sun properties           
     98    Nova::EquPosn  sun  = Nova::GetSolarEquCoords(jd);                   
     99    Nova::ZdAzPosn hrzs = Nova::GetHrzFromEqu(sun, jd);                 
     100
    97101    // Get source position
    98102    Nova::EquPosn pos = FindPointing(time);
    99103
    100     return FACT::PredictI(jd, pos);
     104
     105    // Moon properties
     106    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
     107    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
     108    double        disk = Nova::GetLunarDisk(jd);
     109
     110    // Derived moon properties
     111    double angle = Nova::GetAngularSeparation(moon, pos);
     112    double edist = Nova::GetLunarEarthDist(jd)/384400;
     113
     114    // Current prediction
     115    double sin_malt  = hrzm.alt<0 ? 0 : sin(hrzm.alt*TMath::DegToRad());
     116    double cos_mdist = cos(angle*TMath::DegToRad());
     117    double sin_szd   = sin(hrzs.zd*TMath::DegToRad());
     118
     119    double c0 = pow(disk,      2.63);
     120    double c1 = pow(sin_malt,  0.60);
     121    double c2 = pow(edist,    -2.00);
     122    double c3 = exp(0.67*cos_mdist*cos_mdist*cos_mdist*cos_mdist);
     123    double c4 = exp(-97.8+105.8*sin_szd*sin_szd);
     124
     125    double cur = 6.2 + 95.7*c0*c1*c2*c3 + c4;
     126
     127    return cur;
    101128}
    102129
    … …  
    340367}
    341368
     369Int_t PlotSqm(TArrayD **vec, TString fname)
     370{
     371    fname += ".SQM_CONTROL_DATA.fits";
     372
     373    fits file(fname.Data());
     374    if (!file)
     375    {
     376        cerr << fname << ": " << gSystem->GetError() << endl;
     377        return -2;
     378    }
     379
     380    //cout << fname << endl;
     381
     382    Double_t time;
     383    Float_t mag; // magnitude
     384
     385    if (!file.SetPtrAddress("Time",  &time))
     386        return -1;
     387    if (!file.SetPtrAddress("Mag",  &mag))
     388        return -1;
     389
     390    //const int marker_style = kFullDotMedium;
     391    const int marker_style = kDot;
     392
     393    TGraph g1;
     394    g1.SetName("SQM");
     395
     396
     397    bool found_first_time = false;
     398    while (file.GetNextRow())
     399        if (Contains(vec, time))
     400        {
     401            g1.SetPoint(g1.GetN(), time*24*3600, mag);
     402        }
     403
     404
     405    g1.SetMinimum(19.0);
     406    g1.SetMaximum(21.5);
     407    g1.SetMarkerColor(kBlack);
     408    g1.SetMarkerStyle(marker_style);
     409    g1.GetXaxis()->SetTimeDisplay(true);
     410    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
     411    g1.GetXaxis()->SetLabelSize(0.12);
     412    g1.GetYaxis()->SetLabelSize(0.1);
     413    g1.GetYaxis()->SetTitle("SQM [mag]");
     414    g1.GetYaxis()->SetTitleOffset(0.2);
     415    g1.GetYaxis()->SetTitleSize(0.1);
     416    g1.DrawClone("AP");
     417
     418    return 0;
     419}
     420
     421Int_t PlotSqmLinear(TArrayD **vec, TString fname)
     422{
     423    fname += ".SQM_CONTROL_DATA.fits";
     424
     425    fits file(fname.Data());
     426    if (!file)
     427    {
     428        cerr << fname << ": " << gSystem->GetError() << endl;
     429        return -2;
     430    }
     431
     432    //cout << fname << endl;
     433
     434    Double_t time;
     435    Float_t mag; // magnitude
     436
     437    if (!file.SetPtrAddress("Time",  &time))
     438        return -1;
     439    if (!file.SetPtrAddress("Mag",  &mag))
     440        return -1;
     441
     442    //const int marker_style = kFullDotMedium;
     443    const int marker_style = kDot;
     444
     445    TGraph g1;
     446    g1.SetName("SQM");
     447
     448
     449    bool found_first_time = false;
     450    while (file.GetNextRow())
     451        if (Contains(vec, time))
     452        {
     453            Double_t mag_double = 4.4 * pow(10, 20) * pow( 10, mag*(-0.4));
     454            g1.SetPoint(g1.GetN(), time*24*3600, mag_double);
     455        }
     456
     457
     458    g1.SetMinimum(1.e12);
     459    g1.SetMaximum(5.e12);
     460    g1.SetMarkerColor(kBlack);
     461    g1.SetMarkerStyle(marker_style);
     462    g1.GetXaxis()->SetTimeDisplay(true);
     463    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
     464    g1.GetXaxis()->SetLabelSize(0.12);
     465    g1.GetYaxis()->SetLabelSize(0.1);
     466    g1.GetYaxis()->SetTitle("SQM lin [phot./(s sr m^2)]");
     467    g1.GetYaxis()->SetTitleOffset(0.2);
     468    g1.GetYaxis()->SetTitleSize(0.1);
     469    g1.DrawClone("AP");
     470
     471    return 0;
     472}
     473
     474Int_t PlotHumidity(TArrayD **vec, TString fname)
     475{
     476    fname += ".FSC_CONTROL_HUMIDITY.fits";
     477
     478    fits file(fname.Data());
     479    if (!file)
     480    {
     481        cerr << fname << ": " << gSystem->GetError() << endl;
     482        return -2;
     483    }
     484
     485    //cout << fname << endl;
     486
     487    Double_t time;
     488    Float_t H[4];
     489
     490    if (!file.SetPtrAddress("Time",  &time))
     491        return -1;
     492    if (!file.SetPtrAddress("H",  H))
     493        return -1;
     494
     495    //const int marker_style = kFullDotMedium;
     496    const int marker_style = kDot;
     497
     498    TGraph g1;
     499    TGraph g2;
     500    TGraph g3;
     501    TGraph g4;
     502    //TGraph g5;
     503    g1.SetName("H0");
     504    g2.SetName("H1");
     505    g3.SetName("H2");
     506    g4.SetName("H3");
     507    //g5.SetName("PFmini");
     508
     509
     510    Double_t first_time, last_time;
     511    bool found_first_time = false;
     512    while (file.GetNextRow())
     513        if (Contains(vec, time))
     514        {
     515            if (!found_first_time){
     516              first_time = time*24*3600;
     517              found_first_time = true;
     518            }           
     519            g1.SetPoint(g1.GetN(), time*24*3600, H[0]);
     520            g2.SetPoint(g2.GetN(), time*24*3600, H[1]);
     521            g3.SetPoint(g3.GetN(), time*24*3600, H[2]);
     522            g4.SetPoint(g4.GetN(), time*24*3600, H[3]);
     523            //g5.SetPoint(g5.GetN(), time*24*3600, I[319]);
     524            last_time = time*24*3600;
     525        }
     526
     527
     528    g1.SetMinimum(10);
     529    g1.SetMaximum(80);
     530    g1.SetMarkerColor(kAzure);
     531    g1.SetMarkerStyle(marker_style);
     532    g1.GetXaxis()->SetTimeDisplay(true);
     533    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
     534    g1.GetXaxis()->SetLabelSize(0.12);
     535    g1.GetYaxis()->SetLabelSize(0.1);
     536    g1.GetYaxis()->SetTitle("HUMITIDY");
     537    g1.GetYaxis()->SetTitleOffset(0.2);
     538    g1.GetYaxis()->SetTitleSize(0.1);
     539    g1.DrawClone("AP");
     540
     541    g2.SetMarkerColor(kAzure+1);
     542    g2.SetMarkerStyle(marker_style);
     543    g2.DrawClone("P");
     544   
     545    g3.SetMarkerColor(kAzure+3);
     546    g3.SetMarkerStyle(marker_style);
     547    g3.DrawClone("P");
     548   
     549    g4.SetMarkerColor(kAzure+6);
     550    g4.SetMarkerStyle(marker_style);
     551    g4.DrawClone("P");
     552   
     553    //g5.SetMarkerColor(kAzure+1);
     554    //g5.SetMarkerStyle(kFullDotMedium);
     555    //g5.DrawClone("P");
     556   
     557    g1.DrawClone("P");
     558
     559    TLine l1(first_time-600, 40, last_time+600, 40);
     560    l1.SetLineColor(kOrange);
     561    l1.DrawClone();
     562    TText t1(first_time-600, 41, "Please, note in logbook");
     563    t1.SetTextSize(0.1);
     564    t1.DrawClone();
     565
     566   
     567    TLine l2(first_time-600, 55, last_time+600, 55);
     568    l2.SetLineColor(kRed);
     569    l2.DrawClone();
     570    TText t2(first_time-600, 56, "Please, report to fact-online");
     571    t2.SetTextSize(0.1);
     572    t2.DrawClone();
     573
     574
     575    return 0;
     576}
     577
     578Int_t PlotHumidity2(TArrayD **vec, TString fname)
     579{
     580    fname += ".PFMINI_CONTROL_DATA.fits";
     581
     582    fits file(fname.Data());
     583    if (!file)
     584    {
     585        cerr << fname << ": " << gSystem->GetError() << endl;
     586        return -2;
     587    }
     588
     589    //cout << fname << endl;
     590
     591    Double_t time;
     592    Float_t H;
     593
     594    if (!file.SetPtrAddress("Time",  &time))
     595        return -1;
     596    if (!file.SetPtrAddress("Humidity", &H))
     597        return -1;
     598
     599    const int marker_style = kFullDotMedium;
     600    //const int marker_style = kDot;
     601
     602    TGraph g1;
     603    g1.SetName("PFmini");
     604
     605
     606    while (file.GetNextRow())
     607        if (Contains(vec, time))
     608        {
     609            g1.SetPoint(g1.GetN(), time*24*3600, H);
     610        }
     611   
     612    g1.SetMarkerStyle(marker_style);
     613    g1.SetMarkerColor(kGreen);
     614    g1.DrawClone("P");
     615    return 0;
     616}
    342617
    343618Int_t PlotPointing(TArrayD **vec, TString fname)
    … …  
    420695    g.GetXaxis()->SetLabelSize(0.1);
    421696    g.GetYaxis()->SetLabelSize(0.1);
    422     g.GetYaxis()->SetTitle("TEMP");
     697    g.GetYaxis()->SetTitle("TEMPERATURE");
    423698    g.GetYaxis()->SetTitleOffset(0.2);
    424699    g.GetYaxis()->SetTitleSize(0.1);
    … …  
    542817    if (!file.SetPtrAddress("Time",  &time))
    543818        return -1;
    544     if (!file.SetPtrAddress("Data1", temp) &&
    545         !file.SetPtrAddress("temp", temp))
     819//    if (!file.SetPtrAddress("Data1", temp) &&
     820//        !file.SetPtrAddress("temp", temp))
     821    if (!file.SetPtrAddress("temp", temp))
    546822        return -1;
    547823
    … …  
    669945    }
    670946
    671     TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960);
    672     c->Divide(1, 6, 1e-5, 1e-5);
     947    //check if the sqm was already installed on the telescope                                                                                                       
     948    TCanvas *c = NULL;
     949    TString datestring = Form("%04d%02d%02d", y, m, d);
     950    if( datestring.Atoi() > 20140723 ) {
     951      TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 1280);
     952      c->Divide(1, 8, 1e-5, 1e-5);
     953    }else{
     954      TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 1120);
     955      c->Divide(1, 7, 1e-5, 1e-5);
     956    }
    673957
    674958    gROOT->SetSelectedPad(0);
    … …  
    7161000    gPad->SetRightMargin(0.001);
    7171001    gPad->SetLeftMargin(0.04);
     1002    gPad->SetBottomMargin(0);
    7181003    cout << PlotTemperature1(runs, fname) << endl;
    7191004    cout << PlotTemperature2(runs, fname) << endl;
    … …  
    7211006    cout << PlotTemperature4(runs, fname) << endl;
    7221007
     1008    gROOT->SetSelectedPad(0);
     1009    c->cd(7);
     1010    gPad->SetGrid();
     1011    gPad->SetTopMargin(0);
     1012    gPad->SetBottomMargin(0);
     1013    gPad->SetRightMargin(0.001);
     1014    gPad->SetLeftMargin(0.04);
     1015    cout << PlotHumidity(runs, fname) << endl;
     1016    cout << PlotHumidity2(runs, fname) << endl;
     1017
     1018    //check if the sqm was already installed
     1019    if( datestring.Atoi() > 20140723 ) {
     1020      gROOT->SetSelectedPad(0);
     1021      c->cd(7);
     1022      gPad->SetGrid();
     1023      gPad->SetTopMargin(0);
     1024      gPad->SetBottomMargin(0);
     1025      gPad->SetRightMargin(0.001);
     1026      gPad->SetLeftMargin(0.04);
     1027      cout << PlotHumidity(runs, fname) << endl;
     1028      cout << PlotHumidity2(runs, fname) << endl;
     1029
     1030      gROOT->SetSelectedPad(0);
     1031      c->cd(8);
     1032      gPad->SetGrid();
     1033      gPad->SetTopMargin(0);
     1034      gPad->SetRightMargin(0.001);
     1035      gPad->SetLeftMargin(0.04);
     1036      cout << PlotSqm(runs, fname) << endl;
     1037    }else{
     1038
     1039      gROOT->SetSelectedPad(0);
     1040      c->cd(7);
     1041      gPad->SetGrid();
     1042      gPad->SetTopMargin(0);
     1043      gPad->SetRightMargin(0.001);
     1044      gPad->SetLeftMargin(0.04);
     1045      cout << PlotHumidity(runs, fname) << endl;
     1046      cout << PlotHumidity2(runs, fname) << endl;
     1047    }
     1048
    7231049    c->SaveAs(Form("%s/%04d%02d%02d.png", outpath, y, m, d));
    7241050

  • branches/MarsGapdTimeJitter/fact/processing/buildseqentries.C

    r15237 r18282  
    584584
    585585        TString query;
    586         query  = "SELECT fRunID FROM runinfo WHERE fExcludedFDAKEY=1 AND ";
     586        query  = "SELECT fRunID FROM RunInfo WHERE fExcludedFDAKEY=1 AND ";
    587587        query += cond;
    588588        query += " AND fRunTypeKEY=4 AND fRunID BETWEEN ";
    … …  
    599599
    600600            // Get DRS file with pedestal (roi<1024)
    601             query  = "SELECT MAX(fRunID) FROM runinfo WHERE ";
     601            query  = "SELECT MAX(fRunID) FROM RunInfo WHERE ";
    602602            query += cond;
    603603            query += " AND fHasDrsFile=1 AND (fDrsStep IS NULL OR fDrsStep=2) AND fRunTypeKEY=2";

  • branches/MarsGapdTimeJitter/fact/processing/drstemp.C

    r17145 r18282  
    163163
    164164    float temp[160];
    165     if (file.SetPtrAddress("temp", temp, 82))
    166     {
    167         drstemp82(file, beg, end);
    168         return;
    169     }
     165//    if (file.SetPtrAddress("temp", temp, 82))
     166//    {
     167//        drstemp82(file, beg, end);
     168//        return;
     169//    }
    170170
    171171    file.SetPtrAddress("temp", temp, 160);

  • branches/MarsGapdTimeJitter/fact/processing/fillratescan.C

    r17102 r18282  
    4848Float_t GetOffset(TString fname, Double_t beg, Double_t end)
    4949{
    50     fname.ReplaceAll("RATE_SCAN_DATA", "FEEDBACK_DEVIATION");
     50    if (end < 15937)
     51        fname.ReplaceAll("RATE_SCAN_DATA", "FEEDBACK_DEVIATION");
     52    else
     53        fname.ReplaceAll("RATE_SCAN_DATA", "FEEDBACK_CALIBRATED_CURRENTS");
    5154
    5255    fits file(fname.Data());
    … …  
    6164        return -100;
    6265
    63     if (!file.SetPtrAddress("DeltaUser", &delta))
    64         return -100;
     66    if (end < 15937)
     67    {
     68        if (!file.SetPtrAddress("DeltaUser", &delta))
     69            return -100;
     70    }
     71    else
     72    {
     73        if (!file.SetPtrAddress("U_nom", &delta))
     74            return -100;
     75    }
    6576
    6677    //cout << "Search for: " << beg-15773 << " " << end-15773 << endl;
    … …  
    8091Float_t GetCurrent(TString fname, Double_t beg, Double_t end)
    8192{
    82     fname.ReplaceAll("RATE_SCAN_DATA", "CALIBRATED_CURRENTS");
    83     fname = gSystem->BaseName(fname.Data());
    84 
    85     fname.Prepend("/scratch_nfs/calibrated_currents/");
     93    fname.ReplaceAll("RATE_SCAN_DATA", "FEEDBACK_CALIBRATED_CURRENTS");
     94
     95    //the next two lines are needed for ISDC and data before 11.3.2013
     96    //fname.ReplaceAll("RATE_SCAN_DATA", "CALIBRATED_CURRENTS");
     97    //fname = gSystem->BaseName(fname.Data());
     98    //fname.Prepend("/scratch_nfs/calibrated_currents/");
    8699
    87100    fits file(fname.Data());
    … …  
    326339        !file.SetPtrAddress(string(old ? "Data3" : "RelOnTime"), &ontime))
    327340        return -1;
    328     */
     341        */
    329342    Double_t  *ptime   = file.SetPtrAddress("Time");
    330343    ULong64_t *pid     = file.SetPtrAddress(old ? "Data0" : "Id");

  • branches/MarsGapdTimeJitter/mbase/MStatusDisplay.cc

    r17786 r18282  
    26222622
    26232623    // FIXME: Use mkstemp instead
    2624     if (!mktemp(const_cast<char*>(tmp.Data())))
     2624    if (!mkstemp(const_cast<char*>(tmp.Data())))
    26252625    {
    26262626        *fLog << err << "ERROR - MStatusDisplay::UpdatePSHeader: mktemp failed." << endl;

  • branches/MarsGapdTimeJitter/mcore/DrsCalib.h

    r17757 r18282  
    964964
    965965        fStat.resize(samples*channels);
     966    }
     967
     968    void Reset()
     969    {
     970        for (auto it=fStat.begin(); it!=fStat.end(); it++)
     971        {
     972            it->first = 0;
     973            it->second = 0;
     974        }
    966975    }
    967976
    … …  
    11441153                const double valw = it->second;
    11451154
    1146                 it->first = sumv>0 ? n*(1-s*sumw/sumv) : 0;
     1155                it->first  = sumv>0 ? n*(1-s*sumw/sumv) : 0;
    11471156
    11481157                sumv += valv;
    … …  
    11631172                sumw += valw;
    11641173
    1165                 it->first = sumv>0 ? n*(s*sumw/sumv-1) : 0;
     1174                it->first  = sumv>0 ? n*(s*sumw/sumv-1) : 0;
    11661175            }
    11671176
    … …  
    11971206        return pos-Offset(ch, pos);
    11981207    }
     1208
     1209    /*
     1210    // See MDrsCalibrationTime
     1211    std::string ReadFitsImp(const std::string &str)
     1212    {
     1213        fits file(str);
     1214        if (!file)
     1215        {
     1216            std::ostringstream msg;
     1217            msg << "Could not open file '" << str << "': " << strerror(errno);
     1218            return msg.str();
     1219        }
     1220
     1221        if (file.GetStr("TELESCOP")!="FACT")
     1222        {
     1223            std::ostringstream msg;
     1224            msg << "Reading '" << str << "' failed: Not a valid FACT file (TELESCOP not FACT in header)";
     1225            return msg.str();
     1226        }
     1227
     1228        if (file.GetNumRows()!=1)
     1229        {
     1230            std::ostringstream msg;
     1231            msg << "Reading '" << str << "' failed: Number of rows in table is not 1.";
     1232            return msg.str();
     1233        }
     1234
     1235        fNumSamples  = file.GetUInt("NROI");
     1236        fNumChannels = file.GetUInt("NCH");
     1237        fNumEntries  = file.GetUInt("NBTIME");
     1238
     1239        fStat.resize(fNumSamples*fNumChannels);
     1240
     1241        double *f = reinterpret_cast<double*>(file.SetPtrAddress("CellWidthMean"));
     1242        double *s = reinterpret_cast<double*>(file.SetPtrAddress("CellWidthRms"));
     1243
     1244        if (!file.GetNextRow())
     1245        {
     1246            std::ostringstream msg;
     1247            msg << "Reading data from " << str << " failed.";
     1248            return msg.str();
     1249        }
     1250
     1251        for (uint32_t i=0; i<fNumSamples*fNumChannels; i++)
     1252        {
     1253            fStat[i].first  = f[i];
     1254            fStat[i].second = s[i];
     1255        }
     1256
     1257        return std::string();
     1258    }
     1259
     1260    std::string WriteFitsImp(const std::string &filename, uint32_t night=0) const
     1261    {
     1262        ofits file(filename.c_str());
     1263        if (!file)
     1264        {
     1265            std::ostringstream msg;
     1266            msg << "Could not open file '" << filename << "': " << strerror(errno);
     1267            return msg.str();
     1268        }
     1269
     1270        file.SetDefaultKeys();
     1271        file.AddColumnDouble(fStat.size(), "CellWidthMean", "ratio", "Relative cell width mean");
     1272        file.AddColumnDouble(fStat.size(), "CellWidthRms",  "ratio", "Relative cell width rms");
     1273
     1274        file.SetInt("ADCRANGE", 2000,    "Dynamic range of the ADC in mV");
     1275        file.SetInt("DACRANGE", 2500,    "Dynamic range of the DAC in mV");
     1276        file.SetInt("ADC",      12,      "Resolution of ADC in bits");
     1277        file.SetInt("DAC",      16,      "Resolution of DAC in bits");
     1278        file.SetInt("NPIX",     1440,    "Number of channels in the camera");
     1279        file.SetInt("NTM",      0,       "Number of time marker channels");
     1280        file.SetInt("NROI",     fNumSamples,  "Region of interest");
     1281        file.SetInt("NCH",      fNumChannels, "Number of chips");
     1282        file.SetInt("NBTIME",   fNumEntries,  "Num of entries for time calibration");
     1283
     1284        file.WriteTableHeader("DrsCellTimes");
     1285        if (night>0)
     1286            file.SetInt("NIGHT", night, "Night as int");
     1287
     1288        std::vector<double> data(fNumSamples*fNumChannels*2);
     1289
     1290        for (uint32_t i=0; i<fNumSamples*fNumChannels; i++)
     1291        {
     1292            data[i] = fStat[i].first;
     1293            data[fNumSamples*fNumChannels+i] = fStat[i].second;
     1294        }
     1295
     1296        if (!file.WriteRow(data.data(), data.size()*sizeof(double)))
     1297        {
     1298            std::ostringstream msg;
     1299            msg << "Writing data to " << filename << " failed.";
     1300            return msg.str();
     1301        }
     1302
     1303        return std::string();
     1304    }*/
    11991305};
    12001306

  • branches/MarsGapdTimeJitter/mcore/factfits.h

    r17284 r18282  
    3737            readDrsCalib(fname);
    3838    }
     39
     40        const std::vector<int16_t> &GetOffsetCalibration() const { return fOffsetCalibration; }
    3941
    4042private:

  • branches/MarsGapdTimeJitter/mdrs/MCalibrateDrsTimes.cc

    r17759 r18282  
    152152            continue;
    153153
    154         const Float_t signal = (*fSignals)[sw].GetArrivalTime();
     154        const Float_t signal = (*fSignals)[sw].GetArrivalTimeHiGain();
     155        const Float_t slope  = (*fSignals)[sw].GetArrivalTimeHiGainError();
    155156        const Float_t offset = fCalib ? fCalib->GetOffset(hw, start[hw], signal) : 0;
    156         const Float_t delay  = fCalib ? fCalib->GetDelay(sw) : 0;
     157        const Float_t offset2 = (fCalib && (signal-slope)>=0) ? fCalib->GetOffset(hw, start[hw], signal-slope) : 0;
     158        const Float_t delay  = fCalib ? fCalib->GetDelay(hw) : 0;
    157159
    158160        //if (fIsTimeMarker)
    … …  
    160162
    161163        // convert from slices to ns
    162         const Float_t utime = 1000*(signal       )/fFreq-delay; // [ns]
    163         const Float_t time  = 1000*(signal-offset)/fFreq-delay; // [ns]
     164        const Float_t utime      = 1000*(signal       )/fFreq-delay;  // [ns]
     165        const Float_t time       = 1000*(signal-offset)/fFreq-delay;  // [ns]
     166        const Float_t slopecal   = (slope-offset+offset2)<0 ? -1 : 1000*(slope-offset+offset2)/fFreq; // [ns]
     167        const Float_t uslope     = slope<0 ? -1 : 1000*(slope)/fFreq;                // [ns]
    164168
    165169        /*
    … …  
    172176        {
    173177            (*fArrivalTime)[idx[j]].SetArrivalTime(time);
     178            (*fArrivalTime)[idx[j]].SetTimeSlope(slopecal);
    174179            if (fArrivalTimeU)
     180            {
    175181                (*fArrivalTimeU)[idx[j]].SetArrivalTime(utime);
     182                (*fArrivalTimeU)[idx[j]].SetTimeSlope(uslope);
     183            }
    176184        }
    177185    }

  • branches/MarsGapdTimeJitter/mdrs/MDrsCalibration.h

    r17758 r18282  
    3030    }
    3131
    32     bool ReadFits(TString str)
     32    bool ReadFits(TString fname)
    3333    {
    34         gSystem->ExpandPathName(str);
     34        gSystem->ExpandPathName(fname);
    3535
    36         const std::string msg = ReadFitsImp(str.Data());
     36        std::string msg;
     37        try
     38        {
     39            msg = ReadFitsImp(fname.Data());
     40        }
     41        catch (const std::exception &e)
     42        {
     43            msg = e.what();
     44        }
     45
    3746        if (msg.empty())
     47        {
     48            *fLog << inf << "Read DRS calibration file " << fname << std::endl;
    3849            return true;
     50        }
    3951
    40         *fLog << err << msg << std::endl;
     52        *fLog << err << "Error reading from " << fname << ": " << msg << std::endl;
    4153        return false;
    4254    }

  • branches/MarsGapdTimeJitter/mdrs/MDrsCalibrationTime.cc

    r14922 r18282  
    1616!
    1717!
    18 !   Author(s): Thomas Bretz  2013 <mailto:thomas.bretz@epfl.ch>
     18!   Author(s): Thomas Bretz  2013 <mailto:tbretz@physik.rwth-aachen.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2013
     20!   Copyright: MAGIC Software Development, 2000-2015
    2121!
    2222!
    … …  
    3131
    3232#include <TH1.h>
     33#include <TGraph.h>
    3334
    3435ClassImp(MDrsCalibrationTime);
    … …  
    4849    return true;
    4950}
     51
     52void MDrsCalibrationTime::SetDelays(const TGraph &g)
     53{
     54    fDelays.assign(g.GetY(), g.GetY()+g.GetN());
     55}
     56
     57bool MDrsCalibrationTime::ReadFits(TString fname)
     58{
     59    gSystem->ExpandPathName(fname);
     60
     61    try
     62    {
     63        fits file(fname.Data());
     64        if (!file)
     65            throw runtime_error(strerror(errno));
     66
     67        if (file.GetStr("TELESCOP")!="FACT")
     68        {
     69            std::ostringstream msg;
     70            msg << "Not a valid FACT file (TELESCOP not FACT in header)";
     71            throw runtime_error(msg.str());
     72        }
     73
     74        if (file.GetNumRows()!=1)
     75            throw runtime_error("Number of rows in table is not 1.");
     76
     77        fNumSamples  = file.GetUInt("NROI");
     78        fNumChannels = file.GetUInt("NCH");
     79        fNumEntries  = file.GetUInt("NBTIME");
     80
     81        const double *d = reinterpret_cast<double*>(file.SetPtrAddress("CellOffset"));
     82        if (!file.GetNextRow())
     83            throw runtime_error("Read error.");
     84
     85        fOffsets.assign(d, d+fNumSamples*fNumChannels),
     86        fDelays.resize(0);
     87    }
     88    catch (const exception &e)
     89    {
     90        *fLog << err << "Error reading from " << fname << ": " << e.what() << endl;
     91        return false;
     92    }
     93
     94    *fLog << inf << "Read DRS calibration file " << fname << endl;
     95    return true;
     96}
     97
     98bool MDrsCalibrationTime::WriteFits(const string &fname) const
     99{
     100    const Bool_t exists = !gSystem->AccessPathName(fname.c_str(), kFileExists);
     101    if (exists)
     102    {
     103        *fLog << err << "File '" << fname << "' already exists." << endl;
     104        return false;
     105    }
     106
     107    try
     108    {
     109        ofits file(fname.c_str());
     110        if (!file)
     111            throw runtime_error(strerror(errno));
     112
     113        file.SetDefaultKeys();
     114        file.AddColumnDouble(fNumSamples*fNumChannels, "CellOffset", "samples", "Integral cell offset");
     115
     116        file.SetInt("ADCRANGE", 2000,    "Dynamic range of the ADC in mV");
     117        file.SetInt("DACRANGE", 2500,    "Dynamic range of the DAC in mV");
     118        file.SetInt("ADC",      12,      "Resolution of ADC in bits");
     119        file.SetInt("DAC",      16,      "Resolution of DAC in bits");
     120        file.SetInt("NPIX",     1440,    "Number of channels in the camera");
     121        file.SetInt("NTM",      0,       "Number of time marker channels");
     122        file.SetInt("NROI",     fNumSamples,  "Region of interest");
     123        file.SetInt("NCH",      fNumChannels, "Number of chips");
     124        file.SetInt("NBTIME",   fNumEntries,  "Num of entries for time calibration");
     125
     126        file.WriteTableHeader("DrsCellTimes");
     127        //file.SetInt("NIGHT", night, "Night as int");
     128
     129        /*
     130        file.SetStr("DATE-OBS", fDateObs, "First event of whole DRS calibration");
     131        file.SetStr("DATE-END", fDateEnd, "Last event of whole DRS calibration");
     132        file.SetStr("RUN-BEG", fDateRunBeg[0], "First event of run 0");
     133        file.SetStr("RUN-END", fDateRunEnd[0], "Last event of run 0");
     134        */
     135
     136        if (!file.WriteRow(fOffsets.data(), fOffsets.size()*sizeof(double)))
     137            throw runtime_error("Write error.");
     138    }
     139    catch (const exception &e)
     140    {
     141        *fLog << err << "Error writing to " << fname << ": " << e.what() << endl;
     142        return false;
     143    }
     144
     145    *fLog << inf << "Wrote DRS calibration file " << fname << endl;
     146    return true;
     147}

  • branches/MarsGapdTimeJitter/mdrs/MDrsCalibrationTime.h

    r17697 r18282  
    1010
    1111class TH1;
     12class TGraph;
    1213
    13 class MDrsCalibrationTime : public MParContainer, public DrsCalibrateTime
     14class MDrsCalibrationTime : public MParContainer//, public DrsCalibrateTime
    1415{
     16    int64_t fNumEntries;
     17
     18    size_t fNumSamples;
     19    size_t fNumChannels;
     20
     21    std::vector<double> fOffsets;
    1522    std::vector<double> fDelays;
    1623
    … …  
    1825    MDrsCalibrationTime(const char *name=0, const char *title=0)
    1926    {
    20         fName  = name ? name : "MDrsCalibrationTime";
     27        fName  = name  ? name : "MDrsCalibrationTime";
    2128        fTitle = title ? title : "";
    2229    }
    … …  
    2431    void InitSize(uint16_t channels, uint16_t samples)
    2532    {
    26         //fDelays.clear();
    27         //fDelays.resize(channels);
    28 
    29         DrsCalibrateTime::InitSize(channels, samples);
     33        fNumSamples  = samples;
     34        fNumChannels = channels;
    3035    }
    3136
    3237    void SetCalibration(const DrsCalibrateTime &cal)
    3338    {
    34         *static_cast<DrsCalibrateTime*>(this) = cal;
     39        fNumEntries  = cal.fNumEntries,
     40        fNumSamples  = cal.fNumSamples,
     41        fNumChannels = cal.fNumChannels,
     42
     43        fOffsets.resize(fNumSamples*fNumChannels);
     44
     45        for (size_t c=0; c<fNumChannels; c++)
     46            for (size_t s=0; s<fNumSamples; s++)
     47                fOffsets[c*fNumSamples+s] = cal.Offset(c, s);
    3548    }
    3649
    3750    bool SetDelays(const TH1 &cam);
     51    void SetDelays(const TGraph &g);
    3852
    39     double GetOffset(int hw, int spos, float tm) const
     53    double GetOffset(uint32_t hw, uint32_t spos, float tm) const
    4054    {
    41         return Offset(hw/9, fmod(tm+spos, 1024)) - Offset(hw/9, spos);
     55        const uint32_t ch = (hw/9)*fNumSamples;
     56        return fOffsets[ch + fmod(tm+spos, fNumSamples)] - fOffsets[ch + spos];
    4257    }
    4358
    44     double GetDelay(int sw) const
     59    double GetDelay(int hw) const
    4560    {
    46         return fDelays.size()==0 ? 0 : fDelays[sw];
     61        return fDelays.size()==0 ? 0 : fDelays[hw];
    4762    }
    4863
    49     ClassDef(MDrsCalibrationTime, 2) // A list of histograms storing the Fadc spektrum of one pixel
     64    bool ReadFits(TString fname);
     65    bool WriteFits(const std::string &fname) const;
     66
     67    ClassDef(MDrsCalibrationTime, 3) // A list of histograms storing the Fadc spektrum of one pixel
    5068};
    5169

  • branches/MarsGapdTimeJitter/mdrs/MHDrsCalibrationTime.cc

    r17760 r18282  
    1616!
    1717!
    18 !   Author(s): Thomas Bretz 2013 <mailto:tbretz@phys.ethz.ch>
    19 !
    20 !   Copyright: MAGIC Software Development, 2000-2014
     18!   Author(s): Thomas Bretz 2013 <mailto:tbretz@physik.rwth-aachen.de>
     19!
     20!   Copyright: MAGIC Software Development, 2000-2015
    2121!
    2222!
    … …  
    3737#include "MStatusDisplay.h"
    3838
     39#include "MDrsCalibrationTime.h"
    3940#include "MPedestalSubtractedEvt.h"
    4041

  • branches/MarsGapdTimeJitter/mdrs/MHDrsCalibrationTime.h

    r14922 r18282  
    22#define MARS_MHDrsCalibrationTime
    33
    4 #ifndef MARS_DrsCalibrationTime
    5 #include "MDrsCalibrationTime.h"
     4#ifndef MARS_DrsCalib
     5#include "DrsCalib.h"
    66#endif
    77
    … …  
    2323    MDrsCalibrationTime    *fCal;     //!
    2424
    25     MDrsCalibrationTime fData; //
     25    DrsCalibrateTime fData; //
    2626
    2727    void InitHistogram();

  • branches/MarsGapdTimeJitter/mfileio/FileIOLinkDef.h

    r18009 r18282  
    1919#pragma link C++ class MWriteRootFile+;
    2020#pragma link C++ class MWriteFitsFile+;
     21
    2122#pragma link C++ class MMatrix+;
    2223

  • branches/MarsGapdTimeJitter/mfileio/MFitsArray.h

    r14935 r18282  
    3030
    3131public:
     32    virtual ~MArrayHelperBase() { };
     33
    3234   ofits* GetFitsTable()     {return fFitsTable;}
    3335   UInt_t *      GetArraySizePtr()  {return &fArraySize;}

  • branches/MarsGapdTimeJitter/mfileio/MWriteFitsFile.cc

    r17868 r18282  
    519519
    520520         // initialize all columns of the sub-table, defined by the current container
     521         TString containerName = i_subTable->second.GetContainer()->GetName();
    521522         TClass * cl = i_subTable->second.GetContainer()->IsA();
    522          if (!InitColumns(i_table->first, TString(cl->GetName()) + ".",  fitsTable,
    523                           i_subTable->second.GetContainer(), cl)     ) 
     523         if (!InitColumns(i_table->first, containerName + ".",  fitsTable,
     524                          i_subTable->second.GetContainer(), cl)     )
    524525            return kFALSE;
    525526

  • branches/MarsGapdTimeJitter/mfilter/MFMagicCuts.cc

    r18104 r18282  
    316316
    317317    //fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");
    318     fMap[kELeakage] = fMatrix->AddColumn("NewImagePar.fLeakage1");
     318    fMap[kELeakage] = fMatrix->AddColumn("MNewImagePar.fLeakage1");
    319319
    320320    fMap[kEAlpha]   = fMatrix->AddColumn("MHillasSrc.fAlpha");

  • branches/MarsGapdTimeJitter/mhist/MHEvent.cc

    r13365 r18282  
    194194        fHist->SetName("Island Index");
    195195        fHist->SetYTitle("Index");
     196        fHist->SetPrettyPalette();
     197        break;
     198     case kEvtTimeSlope:
     199     case kEvtTimeSlopeCleaned:
     200        fHist->SetName("Time Slope");
     201        fHist->SetYTitle("delta_t [ns]");
    196202        fHist->SetPrettyPalette();
    197203        break;
    … …  
    295301        fHist->SetCamContent(*event, 5);
    296302        break;
     303    case kEvtTimeSlope:
     304        fHist->SetCamContent(*event, 13);
     305        break;
     306    case kEvtTimeSlopeCleaned:
     307        fHist->SetCamContent(*event, 14);
     308        break;
    297309    default:
    298310        *fLog << "ERROR - Case " << (int)fType << " not implemented..." << endl;

  • branches/MarsGapdTimeJitter/mhist/MHEvent.h

    r13365 r18282  
    2727        kEvtCleaningLevels, kEvtCleaningData,
    2828        kEvtIdxMax, kEvtArrTime, kEvtArrTimeCleaned,
    29         kEvtTrigPix, kEvtIslandIndex
     29        kEvtTrigPix, kEvtIslandIndex, kEvtTimeSlope,
     30        kEvtTimeSlopeCleaned
    3031    };
    3132

  • branches/MarsGapdTimeJitter/mimage/MNewImagePar.cc

    r9374 r18282  
    220220        const Double_t dzx   =  hillas.GetCosDelta()*dx + hillas.GetSinDelta()*dy; // [mm]
    221221        const Double_t dzy   = -hillas.GetSinDelta()*dx + hillas.GetCosDelta()*dy; // [mm]
    222         const Double_t dz    =  gpix.GetT()*gpix.GetT()/4;
    223         const Double_t tana  =  dzy*dzy/(dzx*dzx);
    224         const Double_t distr =  (1+tana)/(rl + tana*rw);
    225         if (distr>dist0-dz || dzx==0)
     222        const Double_t dz    =  gpix.GetT()/2;
     223        const Double_t distr =  (dzy*dzy+dzx*dzx)/(dzx*dzx*rl + dzy*dzy*rw);
     224        if ((dzx==0 && dzy==0) || sqrt(distr)>sqrt(dist0)-dz)
    226225            fConcCore += nphot;
    227226

  • branches/MarsGapdTimeJitter/mjobs/MJSimulation.cc

    r18107 r18282  
    244244    write.AddContainer("MRawEvtHeader",       "Events");
    245245    write.AddContainer("MMcEvt",              "Events", kFALSE);
     246    write.AddContainer("MMcEvtBasic",         "Events", kFALSE);
    246247    write.AddContainer("IncidentAngle",       "Events", kFALSE);
    247248    write.AddContainer("MPointingPos",        "Events", kFALSE);
     249    write.AddContainer("MSimSourcePos",       "Events", kFALSE);
    248250}
    249251
    … …  
    439441    MParSpline splinecones("ConesAngularAcceptance");
    440442    MParSpline splinecones2("ConesTransmission");
     443    MParSpline splineAdditionalPhotonAcceptance("AdditionalPhotonAcceptance");
    441444    plist.AddToList(&splinepde);
    442445    plist.AddToList(&splinemirror);
    443446    plist.AddToList(&splinecones);
    444447    plist.AddToList(&splinecones2);
     448    plist.AddToList(&splineAdditionalPhotonAcceptance);
    445449
    446450    const TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
    … …  
    458462    MSimAbsorption cones("SimConesAngularAcceptance");
    459463    MSimAbsorption cones2("SimConesTransmission");
     464    MSimAbsorption additionalPhotonAcceptance("SimAdditionalPhotonAcceptance");
    460465    absapd.SetParName("PhotonDetectionEfficiency");
    461466    absmir.SetParName("MirrorReflectivity");
    … …  
    463468    cones.SetUseTheta();
    464469    cones2.SetParName("ConesTransmission");
     470    additionalPhotonAcceptance.SetParName("AdditionalPhotonAcceptance");
    465471 
    466472    MSimPointingPos pointing;
    … …  
    654660    write3af.AddContainer("MRawEvtData",      "Events");
    655661    write3af.AddContainer("MTruePhotonsPerPixelCont", "Events");
    656     write3af.AddContainer("MPhotonEvent","Events");
    657662
    658663    write3ar.AddContainer("ElectronicNoise",  "RunHeaders", kTRUE, 1);
    … …  
    804809            tasks.AddToList(&absapd);
    805810            tasks.AddToList(&absmir);
     811            tasks.AddToList(&additionalPhotonAcceptance);
    806812            if (0)
    807813                tasks.AddToList(&simatm); // FASTER?

  • branches/MarsGapdTimeJitter/mjoptim/MJOptimizeCuts.cc

    r8907 r18282  
    9595
    9696// Parameter container
    97 #include "MGeomCamMagic.h"
     97#include "MGeomCamFACT.h"
    9898#include "MParameters.h"
    9999#include "MFilterList.h"
    … …  
    150150    MParList parlist;
    151151
    152     MGeomCamMagic geom; // For GetConvMm2Deg
     152    MGeomCamFACT geom; // For GetConvMm2Deg
    153153    parlist.AddToList(&geom);
    154154
    … …  
    274274    MParList parlist;
    275275
    276     MGeomCamMagic geom; // For GetConvMm2Deg
     276    MGeomCamFACT geom; // For GetConvMm2Deg
    277277    parlist.AddToList(&geom);
    278278

  • branches/MarsGapdTimeJitter/msignal/MExtractFACT.cc

    r17835 r18282  
    171171        //
    172172        Float_t max  = -1;
     173        Float_t tmax  = -1;
    173174        if (pmax>pbeg && pmax<pend-1)
    174175        {
    … …  
    194195                        max = exp(a + b*dx + c*dx*dx);
    195196
    196                         // Time is position of maximum
    197                         //time = dx;
    198                         //time += Int_t(pmax-ptr);
     197                        // tmax is position of maximum
     198                        tmax = dx;
     199                        tmax += Int_t(pmax-ptr);
    199200                    }
    200201                }
    … …  
    203204
    204205        Float_t time = -1;
     206        Float_t slope = -1;
    205207        if (max>=0)
    206208        {
    … …  
    208210
    209211            // Time is position of half hight leading edge
    210             pend = std::max(pbeg, pmax-10);
     212            pend = std::max(pbeg, pmax-15);
    211213            for (;pmax>=pend; pmax--)
    212214                if (*pmax<max/2)
    … …  
    217219                time = (max/2-pmax[0])/(pmax[1]-pmax[0]);
    218220                time += Int_t(pmax-ptr);
     221                slope = tmax - time;
    219222            }
    220223        }
    … …  
    223226        {
    224227            time = gRandom->Uniform(rangehi)+fHiGainFirst+1;
    225             max  = pbeg[uint16_t(time)];
     228            max  = ptr[uint16_t(time)];
    226229        }
    227230
    … …  
    232235        pix.SetGainSaturation(0);
    233236
    234         tix.SetArrivalTime(time);
     237        tix.SetArrivalTime(time, slope);
    235238        tix.SetGainSaturation(0);
    236239    }

  • branches/MarsGapdTimeJitter/msignal/MSignalCam.cc

    r9573 r18282  
    584584// 10: as 0, but returns kFALSE if signal <=0
    585585// 11: as 8, but returns kFALSE if signal <=0
     586// 12: time slope
    586587//
    587588Bool_t MSignalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
    … …  
    593594
    594595    // Used inlcudes status unampped
    595     if (!pix.IsPixelUsed() && (type<6 || type==8))
     596    if (!pix.IsPixelUsed() && (type<6 || type==8 || type==14))
    596597        return kFALSE;
    597598
    … …  
    668669        return pix.GetNumPhotons()>0;
    669670
     671    case 13: // time slope
     672        val = pix.GetTimeSlope();
     673        break;
     674
     675    case 14: // time slope
     676        if (pix.IsPixelUnmapped())
     677            return kFALSE;
     678        val = pix.GetTimeSlope();
     679        break;
     680
    670681    case 9:
    671682    default:

  • branches/MarsGapdTimeJitter/msignal/MSignalPix.cc

    r12938 r18282  
    6464//    size of the calibrated data by roughly 0.5%
    6565//
     66// Version 8:
     67// ----------
     68//  - added new time variable fTimeSlope describing the width of the rise time
     69//
    6670////////////////////////////////////////////////////////////////////////////
    6771#include "MSignalPix.h"
    … …  
    8084MSignalPix::MSignalPix(Float_t phot, Float_t errphot) :
    8185    fIsCore(kFALSE), fRing(1), fIdxIsland(-1),
    82     fPhot(phot), fErrPhot(errphot), fArrivalTime(-1)
     86    fPhot(phot), fErrPhot(errphot), fArrivalTime(-1),
     87    fTimeSlope(-1)
    8388{
    8489    MMath::ReducePrecision(fPhot);
    … …  
    9499    fErrPhot     =  0;
    95100    fArrivalTime = -1;
     101    fTimeSlope = -1;
    96102}
    97103

  • branches/MarsGapdTimeJitter/msignal/MSignalPix.h

    r8528 r18282  
    1919    Float_t  fErrPhot;       // the error of fPhot
    2020    Float_t  fArrivalTime;   // Calibrated Arrival Time
     21    Float_t  fTimeSlope;     // Time between half rise time and position of maximum
    2122
    2223public:
    … …  
    2425    MSignalPix(const MSignalPix &pix)
    2526        : fIsCore(pix.fIsCore), fRing(pix.fRing), fIdxIsland(pix.fIdxIsland),
    26         fPhot(pix.fPhot), fErrPhot(pix.fErrPhot), fArrivalTime(pix.fArrivalTime)
     27        fPhot(pix.fPhot), fErrPhot(pix.fErrPhot), fArrivalTime(pix.fArrivalTime),
     28        fTimeSlope(pix.fTimeSlope)
    2729    {
    2830    }
    … …  
    3941        pix.fErrPhot     = fErrPhot;
    4042        pix.fArrivalTime = fArrivalTime;
     43        pix.fTimeSlope   = fTimeSlope;
    4144    }
    4245    void    Print(Option_t *opt = NULL) const;
    … …  
    4649    Float_t GetErrorPhot() const          { return fErrPhot; }
    4750    Float_t GetArrivalTime() const        { return fArrivalTime; }
     51    Float_t GetTimeSlope() const          { return fTimeSlope; }
    4852
    4953    Bool_t  IsPixelUsed() const           { return fRing>0; }
    … …  
    6569    void    Set(Float_t np, Float_t ep)   { MMath::ReducePrecision(np); MMath::ReducePrecision(ep);  fPhot = np; fErrPhot = ep; }
    6670    void    SetArrivalTime(Float_t tm)    { fArrivalTime = tm; }
     71    void    SetTimeSlope(Float_t ts)      { fTimeSlope = ts; }
    6772
    6873    //void    AddNumPhotons(Float_t f)      { fPhot += f; }
    6974    //void    Scale(Float_t f)              { fPhot /= f; }
    7075
    71     ClassDef(MSignalPix, 7)  // class containing information about the Cerenkov Photons in a pixel
     76    ClassDef(MSignalPix, 8)  // class containing information about the Cerenkov Photons in a pixel
    7277};
    7378

  • branches/MarsGapdTimeJitter/msim/MSimPointingPos.cc

    r17846 r18282  
    3939//
    4040// If no view cone option was given and off-target observations are switched
    41 // on by setting fOffTargetDistance!=0 the poitnting position is calculated:
     41// on by setting fOffTargetDistance!=0 the pointing position is calculated:
    4242//
    4343//  1) fOffTargetDistance < 0:
    … …  
    5151//     fOffTargetDistance. (phi==0 is the direction of positive theta)
    5252//
     53// The original zenith and azimuth coordinates of the shower axis are stored in
     54// the MSimSourcePos container.
     55//
     56// If the view cone option was given and off-target observations are switched on
     57// the orientation is fixed to the main direction around the view cone was
     58// produced.
     59// In addition a 'quasi'-simulated source position is calculated,
     60// depending on fOffTargetDistance and fOffTargetPhi (see 1) and 2) above).
     61// The corresponding zenith and azimuth coordinates are stored in the
     62// MSimSourcePos container. This is of course not a physical source position,
     63// but it can be used to determine the performance of wobble analysis on
     64// background events (which are homogenous distributed).
     65//
     66//
    5367//
    5468//  Input Containers:
    … …  
    5872//  Output Containers:
    5973//   MPointingPos
     74//   MSimSourcePos
    6075//
    6176//////////////////////////////////////////////////////////////////////////////
    … …  
    8499//
    85100MSimPointingPos::MSimPointingPos(const char* name, const char *title)
    86     : fRunHeader(0), fEvtHeader(0), fPointing(0),
     101    : fRunHeader(0), fEvtHeader(0), fPointing(0), fSimSourcePosition(0),
    87102    fOffTargetDistance(0), fOffTargetPhi(-1)
    88103
    … …  
    140155    fPointing = (MPointingPos*)pList->FindCreateObj("MPointingPos");
    141156    if (!fPointing)
     157        return kFALSE;
     158
     159    fSimSourcePosition = (MPointingPos*)pList->FindCreateObj("MPointingPos","MSimSourcePos");
     160    if (!fSimSourcePosition)
    142161        return kFALSE;
    143162
    … …  
    180199    {
    181200        *fLog << warn;
    182         *fLog << "WARNING - Combining the view cone option with off-target observations doesn't make sense." << endl;
    183         *fLog << "          Option for off-target observations will be ignored." << endl;
     201        *fLog << "WARNING - Combining the view cone option with off-target pointing can lead to not homogenous events." << endl;
     202        *fLog << "          A simulated source position according to the off-target parameters will be created instead." << endl;
    184203    }
    185204    // FIXME: Check also the enlightened region on the ground!
    … …  
    209228
    210229    // Local sky coordinates (direction of telescope axis)
    211     Double_t zd = viewcone ? fRunHeader->GetZdMin() : fEvtHeader->GetZd()*TMath::RadToDeg();  // x==north
    212     Double_t az = viewcone ? fRunHeader->GetAzMin() : fEvtHeader->GetAz()*TMath::RadToDeg();  // y==west
    213 
    214     if (!viewcone)
    215     {
    216         Double_t dtheta, dphi;
    217         GetDelta(dtheta, dphi);
    218 
    219         const Double_t theta = zd*TMath::DegToRad();
    220         const Double_t phi   = az*TMath::DegToRad();
    221 
    222         TVector3 src, pnt;
    223         src.SetMagThetaPhi(1, theta,        phi);
    224         pnt.SetMagThetaPhi(1, theta+dtheta, phi);
    225 
    226         pnt.Rotate(dphi, src);
    227 
    228         zd = pnt.Theta()*TMath::RadToDeg();
    229         az = pnt.Phi()  *TMath::RadToDeg();
    230     }
     230    Double_t zdCorsika = viewcone ? fRunHeader->GetZdMin() : fEvtHeader->GetZd()*TMath::RadToDeg();  // x==north
     231    Double_t azCorsika = viewcone ? fRunHeader->GetAzMin() : fEvtHeader->GetAz()*TMath::RadToDeg();  // y==west
     232
     233    // Calculate off target position in local sky coordinates
     234    Double_t dtheta, dphi;
     235    GetDelta(dtheta, dphi);
     236
     237    const Double_t theta = zdCorsika*TMath::DegToRad();
     238    const Double_t phi   = azCorsika*TMath::DegToRad();
     239
     240    TVector3 originalVector, wobbleVector;
     241    originalVector.SetMagThetaPhi(1, theta,        phi);
     242    wobbleVector.SetMagThetaPhi(1, theta+dtheta, phi);
     243
     244    wobbleVector.Rotate(dphi, originalVector);
     245
     246    Double_t zdWobble, azWobble;
     247    zdWobble = wobbleVector.Theta()*TMath::RadToDeg();
     248    azWobble = wobbleVector.Phi()  *TMath::RadToDeg();
    231249
    232250    // Transform the corsika coordinate system (north is magnetic north)
    233251    // into the telescopes local coordinate system. Note, that all vectors
    234252    // are already correctly oriented.
    235     az += fRunHeader->GetMagneticFieldAz()*TMath::RadToDeg();
    236 
    237     // Setup the pointing position
    238     fPointing->SetLocalPosition(zd, az);
    239 
    240     // Calculate incident angle between magnetic field direction
    241     // and pointing direction ( phi and theta? )
     253    azCorsika += fRunHeader->GetMagneticFieldAz()*TMath::RadToDeg();
     254    azWobble += fRunHeader->GetMagneticFieldAz()*TMath::RadToDeg();
     255
     256    // Setup the pointing and the simulated source position
     257    if (!viewcone)
     258    {
     259        fPointing->SetLocalPosition(zdWobble, azWobble);
     260        fSimSourcePosition->SetLocalPosition(zdCorsika, azCorsika);
     261    }
     262    else
     263    {
     264        fPointing->SetLocalPosition(zdCorsika, azCorsika);
     265        fSimSourcePosition->SetLocalPosition(zdWobble, azWobble);
     266    }
    242267
    243268    return kTRUE;

  • branches/MarsGapdTimeJitter/msim/MSimPointingPos.h

    r9367 r18282  
    1616    MCorsikaRunHeader *fRunHeader;  //! Header storing event information
    1717    MCorsikaEvtHeader *fEvtHeader;  //! Header storing event information
    18     MPointingPos      *fPointing;   //! Output storing telescope poiting position in local (telescope) coordinate system
     18    MPointingPos      *fPointing;   //! Output storing telescope pointing position in local (telescope) coordinate system
     19    MPointingPos      *fSimSourcePosition;   //! Output storing simulated source pointing position in local (telescope) coordinate system
    1920
    2021    Double_t fOffTargetDistance;    // [rad] Distance of the observed off-target position from the source

  • branches/MarsGapdTimeJitter/msimcamera/MSimCamera.cc

    r18159 r18282  
    155155    }
    156156    // Check all entries for inf and nan. Those are not accepted here.
    157     for( std::vector< double > row : fFixTimeOffsetsBetweenPixelsInNs->fM ){
    158         for( double number : row){
    159 
    160             if( std::isnan(number) || std::isinf(number) ){
    161 
     157    for( int row_index=0; row_index<fFixTimeOffsetsBetweenPixelsInNs->fM.size(); row_index++){
     158        std::vector<double> row = fFixTimeOffsetsBetweenPixelsInNs->fM.at(row_index);
     159        for( int col_index=0; col_index<row.size(); col_index++){
     160            double specific_delay = row.at(col_index);
     161            if( std::isnan(specific_delay) || std::isinf(specific_delay) ){
    162162                *fLog << err << "In Source: "<< __FILE__ <<" in line: ";
    163163                *fLog << __LINE__;
    164164                *fLog << " in function: "<< __func__ <<"\n";
    165                 *fLog << "There is a non normal number in the fix temporal ";
    166                 *fLog << "pixel offsets. This is at least one number is ";
    167                 *fLog << "NaN or Inf. This here is >"<< number;
     165                *fLog << "There is a non normal specific_delay in the fix temporal ";
     166                *fLog << "pixel offsets. This is that at least one specific_delay is ";
     167                *fLog << "NaN or Inf. This here is >"<< specific_delay;
    168168                *fLog << "<... aborting." << endl;               
    169169                return kFALSE;
    170170            }
    171171        }
     172
    172173    }
    173174    // -------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.