Changeset 8882 for trunk


Ignore:
Timestamp:
03/18/08 17:41:53 (17 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r8881 r8882  
    1919                                                 -*-*- END OF LINE -*-*-
    2020
     21 2008/03/18 Thomas Bretz
     22
     23   * sponde.cc:
     24     - added new option "--force-runtime"
     25
     26   * mbase/MEnv.h:
     27     - added WriteFile to context menu
     28
     29   * mjobs/MJSpectrum.[h,cc]:
     30     - added a new option to force using the runtime instead of the
     31       effective observation time (this might bw wrong for very
     32       short datasets)
     33     - added a check if the effective observation time is out of
     34       the histogram range... print a warning if so and include
     35       the overflow bins into the eff. obs time
     36     - added an estimated sensitivity curve for high and low za
     37       to the spectrum plots
     38     - added description text for 1553 and crab spectrum
     39     - write out the MC events after cuts including their weights
     40     - do not fit at 1TeV but 500GeV instead
     41
     42   * mjobs/MJob.cc:
     43     - check in WriteContainer whether the file is already open
     44
     45   * mpointing/MPointingDevCalc.cc:
     46     - added some more comments
     47
     48
     49
    2150 2008/03/14 Daniel Hoehne
    2251
    2352   * datacenter/macros/filldotrun.C:
    2453     - inserted new arehucas version
     54
     55
     56
     57 2008/03/04 Thomas Bretz
     58
     59   * condor/program.submit, condor/macro.submit, condor/script.submit:
     60     - added
    2561
    2662

  • trunk/MagicSoft/Mars/mbase/MEnv.h

    r8741 r8882  
    7373    Int_t       ReadFile(const char *fname, EEnvLevel level);
    7474
     75    Int_t       WriteFile(const char *filename, EEnvLevel level) { return TEnv::WriteFile(filename, level); }
     76    Int_t       WriteFile(const char *filename) { return WriteFile(filename, kEnvLocal); } //*MENU*
     77
    7578    void        PrintEnv(EEnvLevel level = kEnvAll) const;
    7679    void        Print(Option_t *option) const { TEnv::Print(option); }

  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

    r8750 r8882  
    3838#include <TLine.h>
    3939#include <TFile.h>
     40#include <TGraph.h>
    4041#include <TChain.h>
    4142#include <TLatex.h>
    … …  
    5657#include "MHn.h"
    5758#include "MBinning.h"
     59#include "MParameters.h"
    5860#include "MDataSet.h"
    5961#include "MMcCorsikaRunHeader.h"
    … …  
    7981#include "MFDataMember.h"
    8082#include "MContinue.h"
     83#include "MWriteRootFile.h"
    8184
    8285ClassImp(MJSpectrum);
    … …  
    239242    }
    240243
    241     TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",  "TH1D", "OnTime");
    242     TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
    243     if (!vstime || !size)
     244    TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",      "TH1D", "OnTime");
     245    TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess",     "TH1D", "Hist");
     246    TH1D *time   = (TH1D*)arr.FindObjectInCanvas("ExcessTime", "TH1D", "Hist");
     247    if (!vstime || !size || !time)
    244248        return -1;
    245249
    … …  
    280284        return -1;
    281285
    282     return vstime->Integral();
     286    if (vstime->GetBinContent(0)>0)
     287    {
     288        *fLog << err << "ERROR - Undeflow bin of OnTime histogram not empty as it ought to be." << endl;
     289        return -1;
     290    }
     291
     292    const Double_t ofl = vstime->GetBinContent(vstime->GetNbinsX());
     293    const Double_t eff = vstime->Integral()+ofl;
     294    if (ofl>0)
     295        *fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100*ofl/eff, ofl) << " of the eff. observation time is out of histogram range." << endl;
     296
     297    const Double_t obs = time->GetXaxis()->GetXmax()-time->GetXaxis()->GetXmin();
     298    if (fForceRunTime)
     299    {
     300        *fLog << inf;
     301        *fLog << "Total run time: " << obs/60 << "min" << endl;
     302        *fLog << "Eff. obs. time: " << eff/60 << "min  (" << Form("%.1f%%", 100*eff/obs) << ")" << endl;
     303    }
     304
     305    return fForceRunTime ? obs : eff;
    283306}
    284307
    … …  
    479502        // ----- Complete incomplete energy ranges -----
    480503        // ----- and apply production area weights -----
    481         weight.CompleteEnergySpectrum(*hfill, Emin);
    482 
    483         hfill->Scale(scale*scale);
     504        weight.CompleteEnergySpectrum(*hfill, Emin, scale*scale);
    484505
    485506        // Add weighted data from file to final histogram
    … …  
    537558    // Calculate the Probability
    538559    temp1.Divide(&temp2);
    539     temp1.Scale(1./temp1.Integral());
     560    temp1.Scale(1./temp1.Integral(1, temp1.GetNbinsX()));
    540561
    541562    // Some cosmetics: Name, Axis, etc.
    … …  
    849870        {
    850871            const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",  f1, f1);
    851             const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
     872            const TString fmt1 = Form("(\\frac{E}{500GeV})^{%%%sf#pm%%%sf}", f0, f0);
    852873
    853874            str  = Form(fmt0.Data(), p1/exp, e1/exp, np);
    … …  
    966987    spectrum.SetBit(TH1::kNoStats);
    967988
    968     for (int i=0; i<excess.GetNbinsX(); i++)
    969     {
    970         spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
    971         spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
    972 
    973         spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
    974         spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
    975     }
     989    // Minimum number of excessevents to get 3sigma in 1h
     990    TF1 sensl("SensLZA", "85*(x/200)^(-0.55)", 100, 6000);
     991    TF1 sensh("SensHZA", "35*(x/200)^(-0.70)", 100, 1000);
     992    //sens.SetLineColor(kBlue);
     993    //sens.DrawClone("Csame");
    976994
    977995    c1.cd(4);
    … …  
    981999    gPad->SetGridx();
    9821000    gPad->SetGridy();
     1001
     1002    TGraph gsensl;//("Sensitivity LZA", "");
     1003    TGraph gsensh;//("Sensitivity HZA", "");
     1004
     1005    const Double_t sqrton = TMath::Sqrt(ontime/3600.);
     1006
     1007    for (int i=0; i<excess.GetNbinsX(); i++)
     1008    {
     1009        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
     1010        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
     1011
     1012        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
     1013        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
     1014
     1015        if (collarea.GetBinContent(i+1)<=0)
     1016            continue;
     1017
     1018        const Double_t cen = spectrum.GetBinCenter(i+1);
     1019        gsensl.SetPoint(gsensl.GetN(), cen, sensl.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime);
     1020        gsensh.SetPoint(gsensh.GetN(), cen, sensh.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime);
     1021
     1022        cout << cen << "   " << sensl.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime;
     1023        cout <<  "   " << sensh.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime;
     1024        cout << endl;
     1025    }
     1026
    9831027    spectrum.SetMinimum(1e-12);
    9841028    spectrum.SetXTitle("E [GeV]");
    … …  
    9871031
    9881032    TLatex tex;
     1033    tex.SetTextColor(13);
    9891034
    9901035    TF1 fc("Crab", "7.0e-6*(x/300)^(-2.31-0.26*log10(x/300))", 100, 6000);
    991     fc.SetLineStyle(7);
     1036    fc.SetLineStyle(9);
    9921037    fc.SetLineWidth(2);
    9931038    fc.SetLineColor(14);
    9941039    fc.DrawCopy("same");
    9951040
    996     tex.DrawText(90, fc.Eval(100), "Crab (\\Gamma=-2.31)");
     1041    tex.DrawLatex(1250, fc.Eval(1250), "Crab/\\Gamma=-2.3");
    9971042
    9981043    TF1 fa("PG1553", "1.8e-6*(x/200)^-4.21", 90, 600);
    … …  
    10001045    fa.DrawCopy("same");
    10011046
    1002     tex.DrawText(90, fa.Eval(90), "PG1553 (\\Gamma=-4.21)");
     1047    tex.SetTextAlign(23);
     1048    tex.DrawLatex(600, fa.Eval(600), "PG1553/\\Gamma=-4.2");
     1049
     1050    gROOT->SetSelectedPad(0);
     1051
     1052    gsensl.SetLineStyle(5);
     1053    gsensl.SetLineColor(14);
     1054    gsensl.SetLineWidth(2);
     1055    gsensl.DrawClone("C")->SetBit(kCanDelete);
     1056
     1057    gsensh.SetLineStyle(5);
     1058    gsensh.SetLineColor(14);
     1059    gsensh.SetLineWidth(2);
     1060    gsensh.DrawClone("C")->SetBit(kCanDelete);
    10031061
    10041062    // Display dN/dE*E^2 for conveinience
    … …  
    10151073        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
    10161074        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *e*e);
     1075    }
     1076
     1077    for (int i=0; i<gsensl.GetN(); i++)
     1078    {
     1079        const Double_t e = gsensl.GetX()[i]*1e-3;
     1080
     1081        gsensl.GetY()[i] *= e*e;
     1082        gsensh.GetY()[i] *= e*e;
    10171083    }
    10181084
    … …  
    10311097    static_cast<const TAttLine&>(fc).Copy(fa2);
    10321098    fa2.DrawCopy("same");
     1099
     1100    gsensl.DrawClone("C")->SetBit(kCanDelete);
     1101    gsensh.DrawClone("C")->SetBit(kCanDelete);
    10331102
    10341103    // Fit Spectrum
    … …  
    13171386    hist.InitTitle(";Slope;Disp-Dist [\\circ];");
    13181387    hist.SetDrawOption("colz profx");
     1388}
     1389
     1390// --------------------------------------------------------------------------
     1391//
     1392// Setup write to write:
     1393//     container         tree       optional?
     1394//  --------------     ----------  -----------
     1395//   "MHillas"      to  "Events"
     1396//   "MHillasSrc"   to  "Events"
     1397//   "Hadronness"   to  "Events"       yes
     1398//   "MEnergyEst"   to  "Events"       yes
     1399//   "DataType"     to  "Events"
     1400//
     1401void MJSpectrum::SetupWriter(MWriteRootFile *write/*, const char *name*/) const
     1402{
     1403    if (!write)
     1404        return;
     1405
     1406    //write->SetName(name);
     1407    write->AddContainer("MHillas",        "Events");
     1408    write->AddContainer("MHillasSrc",     "Events");
     1409    write->AddContainer("MHillasExt",     "Events");
     1410    //write->AddContainer("MPointingPos",   "Events");
     1411    write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
     1412    write->AddContainer("MImagePar",      "Events", kFALSE);
     1413    write->AddContainer("MNewImagePar",   "Events", kFALSE);
     1414    write->AddContainer("MNewImagePar2",  "Events", kFALSE);
     1415    write->AddContainer("Hadronness",     "Events", kFALSE);
     1416    write->AddContainer("MSrcPosCam",     "Events", kFALSE);
     1417    write->AddContainer("MSrcPosAnti",    "Events", kFALSE);
     1418    write->AddContainer("ThetaSquared",   "Events", kFALSE);
     1419    write->AddContainer("OpticalAxis",    "Events", kFALSE);
     1420    write->AddContainer("Disp",           "Events", kFALSE);
     1421    write->AddContainer("Ghostbuster",    "Events", kFALSE);
     1422    write->AddContainer("MEnergyEst",     "Events", kFALSE);
     1423    write->AddContainer("MTime",          "Events", kFALSE);
     1424    write->AddContainer("MMcEvt",         "Events", kFALSE);
     1425    write->AddContainer("MWeight",        "Events");
     1426    write->AddContainer("DataType",       "Events");
     1427    write->AddContainer("RunNumber",      "Events");
     1428    write->AddContainer("EvtNumber",      "Events");
    13191429}
    13201430
    … …  
    15981708    taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
    15991709
     1710    MWriteRootFile write(GetPathOut());
     1711    SetupWriter(&write);
     1712
     1713    MParameterI par("DataType");
     1714    plist.AddToList(&par);
     1715    par.SetVal(2);
     1716
    16001717    tlist2.AddToList(&read);
    16011718    // If no weighting should be applied but the events should
    … …  
    16211738    tlist2.AddToList(fCut3);
    16221739    tlist2.AddToList(&taskenv2);
     1740    if (!GetPathOut().IsNull())
     1741        tlist2.AddToList(&write);
    16231742    tlist2.AddToList(&fill31);
    16241743    tlist2.AddToList(&fill4);
    … …  
    16701789
    16711790    // Spectrum fitted (convert res[1] from TeV to GeV)
    1672     TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
     1791    TF1 flx("flx", Form("%e*pow(x/500, %f)", res[1]/500, res[0]));
    16731792
    16741793    // Number of events this spectrum would produce per s and m^2
    … …  
    17011820    cont.Add((TObject*)GetEnv()); // const_cast
    17021821    cont.Add((TObject*)&set);     // const_cast
     1822    cont.Add(plist.FindObject("MAlphaFitter"));
    17031823    cont.Add(&area0);
    17041824    cont.Add(&area1);
    … …  
    17081828        cont.Add(fDisplay);
    17091829
    1710     if (!WriteContainer(cont, "", "RECREATE"))
     1830    if (!WriteContainer(cont, "", GetPathOut().IsNull()?"RECREATE":"UPDATE"))
    17111831    {
    17121832        *fLog << err << GetDescriptor() << ": Writing result failed." << endl;

  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.h

    r8709 r8882  
    1919class MAlphaFitter;
    2020class MStatusArray;
     21class MWriteRootFile;
    2122class MHCollectionArea;
    2223class MMcSpectrumWeight;
    … …  
    3637
    3738    Bool_t fForceTheta;
     39    Bool_t fForceRunTime;
    3840
    3941    // Setup Histograms
    … …  
    5052    Bool_t   Refill(MParList &plist, TH1D &h) /*const*/;
    5153    Bool_t   InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
     54
     55    // Write output
     56    void     SetupWriter(MWriteRootFile *write/*, const char *name*/) const;
    5257
    5358    // Display Output
    … …  
    6772    Bool_t Process(const MDataSet &set);
    6873
    69     void ForceTheta(Bool_t b=kTRUE) { fForceTheta=b; }
     74    void ForceTheta(Bool_t b=kTRUE)   { fForceTheta=b; }
     75    void ForceRunTime(Bool_t b=kTRUE) { fForceRunTime=b; }
    7076
    7177    void SetEnergyEstimator(const MTask *task);

  • trunk/MagicSoft/Mars/mjobs/MJob.cc

    r8719 r8882  
    431431    title += fName;
    432432
     433    // In case the update-option is selected check whether
     434    // the file is already open
     435    if (TString(option).Contains("update", TString::kIgnoreCase))
     436    {
     437        TFile *file = dynamic_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(oname));
     438        if (file && file->IsOpen() && !file->IsZombie())
     439        {
     440            *fLog << inf << "Open file found." << endl;
     441            file->cd();
     442            return WriteContainer(cont);
     443        }
     444    }
     445
     446    // Open a new file with the defined option for writing
    433447    TFile file(oname, option, title, compr);
    434     if (!file.IsOpen())
     448    if (!file.IsOpen() || file.IsZombie())
    435449    {
    436450        *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;

  • trunk/MagicSoft/Mars/mpointing/MPointingDevCalc.cc

    r8767 r8882  
    144144// New pointing models have been installed (if the pointing model
    145145// is different, than the previous one it might mean, that also
    146 // the starguider calibration is different.)
     146// the starguider calibration is different.) The date only means
     147// day-time (noon) at which the model has been changed.
    147148//
    148149//   29. Apr. 2004    ~25800
    … …  
    156157//   17. Oct. 2006   ~103130
    157158//   17. Jun. 2007   ~248193
    158 //   18. Oct. 2007
     159//   18. Oct. 2007   ~291039(?)
    159160//
    160161// From 2.2.2006 beginnig of the night (21:05h, run >=81855) to 24.2.2006

  • trunk/MagicSoft/Mars/sponde.cc

    r8767 r8882  
    5151    gLog << " Operation Mode:" << endl;
    5252    gLog << "   --force-theta             Force execution even with non-fitting theta distributions." << endl;
     53    gLog << "   --force-runtime           Force usage of runtime instead of eff. observation time." << endl;
    5354    gLog << endl;
    5455    gLog << " Options:" << endl;
    … …  
    120121
    121122    const Bool_t  kForceTheta    =  arg.HasOnlyAndRemove("--force-theta");
     123    const Bool_t  kForceRunTime  =  arg.HasOnlyAndRemove("--force-runtime");
    122124
    123125    //
    … …  
    244246
    245247        job.ForceTheta(kForceTheta);
     248        job.ForceRunTime(kForceRunTime);
    246249
    247250        if (!job.Process(seq))
Note: See TracChangeset for help on using the changeset viewer.