Changeset 8882 for trunk/MagicSoft/Mars/mjobs

Timestamp:

03/18/08 17:41:53 (17 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mjobs

Files:

• MJSpectrum.cc (modified) (20 diffs)
• MJSpectrum.h (modified) (4 diffs)
• MJob.cc (modified) (1 diff)

trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

@@ -38 +38 @@
 #include <TLine.h>
 #include <TFile.h>
+#include <TGraph.h>
 #include <TChain.h>
 #include <TLatex.h>
@@ -56 +57 @@
 #include "MHn.h"
 #include "MBinning.h"
+#include "MParameters.h"
 #include "MDataSet.h"
 #include "MMcCorsikaRunHeader.h"
@@ -79 +81 @@
 #include "MFDataMember.h"
 #include "MContinue.h"
+#include "MWriteRootFile.h"
 
 ClassImp(MJSpectrum);
@@ -239 +242 @@
     }
 
-    TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",  "TH1D", "OnTime");
-    TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
-    if (!vstime || !size)
+    TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",      "TH1D", "OnTime");
+    TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess",     "TH1D", "Hist");
+    TH1D *time   = (TH1D*)arr.FindObjectInCanvas("ExcessTime", "TH1D", "Hist");
+    if (!vstime || !size || !time)
         return -1;
 
@@ -280 +284 @@
         return -1;
 
-    return vstime->Integral();
+    if (vstime->GetBinContent(0)>0)
+    {
+        *fLog << err << "ERROR - Undeflow bin of OnTime histogram not empty as it ought to be." << endl;
+        return -1;
+    }
+
+    const Double_t ofl = vstime->GetBinContent(vstime->GetNbinsX());
+    const Double_t eff = vstime->Integral()+ofl;
+    if (ofl>0)
+        *fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100*ofl/eff, ofl) << " of the eff. observation time is out of histogram range." << endl;
+
+    const Double_t obs = time->GetXaxis()->GetXmax()-time->GetXaxis()->GetXmin();
+    if (fForceRunTime)
+    {
+        *fLog << inf;
+        *fLog << "Total run time: " << obs/60 << "min" << endl;
+        *fLog << "Eff. obs. time: " << eff/60 << "min  (" << Form("%.1f%%", 100*eff/obs) << ")" << endl;
+    }
+
+    return fForceRunTime ? obs : eff;
 }
 
@@ -479 +502 @@
         // ----- Complete incomplete energy ranges -----
         // ----- and apply production area weights -----
-        weight.CompleteEnergySpectrum(*hfill, Emin);
-
-        hfill->Scale(scale*scale);
+        weight.CompleteEnergySpectrum(*hfill, Emin, scale*scale);
 
         // Add weighted data from file to final histogram
@@ -537 +558 @@
     // Calculate the Probability
     temp1.Divide(&temp2);
-    temp1.Scale(1./temp1.Integral());
+    temp1.Scale(1./temp1.Integral(1, temp1.GetNbinsX()));
 
     // Some cosmetics: Name, Axis, etc.
@@ -849 +870 @@
         {
             const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",  f1, f1);
-            const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
+            const TString fmt1 = Form("(\\frac{E}{500GeV})^{%%%sf#pm%%%sf}", f0, f0);
 
             str  = Form(fmt0.Data(), p1/exp, e1/exp, np);
@@ -966 +987 @@
     spectrum.SetBit(TH1::kNoStats);
 
-    for (int i=0; i<excess.GetNbinsX(); i++)
-    {
-        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
-        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
-
-        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
-        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
-    }
+    // Minimum number of excessevents to get 3sigma in 1h
+    TF1 sensl("SensLZA", "85*(x/200)^(-0.55)", 100, 6000);
+    TF1 sensh("SensHZA", "35*(x/200)^(-0.70)", 100, 1000);
+    //sens.SetLineColor(kBlue);
+    //sens.DrawClone("Csame");
 
     c1.cd(4);
@@ -981 +999 @@
     gPad->SetGridx();
     gPad->SetGridy();
+
+    TGraph gsensl;//("Sensitivity LZA", "");
+    TGraph gsensh;//("Sensitivity HZA", "");
+
+    const Double_t sqrton = TMath::Sqrt(ontime/3600.);
+
+    for (int i=0; i<excess.GetNbinsX(); i++)
+    {
+        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
+        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
+
+        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
+        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
+
+        if (collarea.GetBinContent(i+1)<=0)
+            continue;
+
+        const Double_t cen = spectrum.GetBinCenter(i+1);
+        gsensl.SetPoint(gsensl.GetN(), cen, sensl.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime);
+        gsensh.SetPoint(gsensh.GetN(), cen, sensh.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime);
+
+        cout << cen << "   " << sensl.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime;
+        cout <<  "   " << sensh.Eval(cen)*sqrton/spectrum.GetBinWidth(i+1)*1000/collarea.GetBinContent(i+1)/ontime;
+        cout << endl;
+    }
+
     spectrum.SetMinimum(1e-12);
     spectrum.SetXTitle("E [GeV]");
@@ -987 +1031 @@
 
     TLatex tex;
+    tex.SetTextColor(13);
 
     TF1 fc("Crab", "7.0e-6*(x/300)^(-2.31-0.26*log10(x/300))", 100, 6000);
-    fc.SetLineStyle(7);
+    fc.SetLineStyle(9);
     fc.SetLineWidth(2);
     fc.SetLineColor(14);
     fc.DrawCopy("same");
 
-    tex.DrawText(90, fc.Eval(100), "Crab (\\Gamma=-2.31)");
+    tex.DrawLatex(1250, fc.Eval(1250), "Crab/\\Gamma=-2.3");
 
     TF1 fa("PG1553", "1.8e-6*(x/200)^-4.21", 90, 600);
@@ -1000 +1045 @@
     fa.DrawCopy("same");
 
-    tex.DrawText(90, fa.Eval(90), "PG1553 (\\Gamma=-4.21)");
+    tex.SetTextAlign(23);
+    tex.DrawLatex(600, fa.Eval(600), "PG1553/\\Gamma=-4.2");
+
+    gROOT->SetSelectedPad(0);
+
+    gsensl.SetLineStyle(5);
+    gsensl.SetLineColor(14);
+    gsensl.SetLineWidth(2);
+    gsensl.DrawClone("C")->SetBit(kCanDelete);
+
+    gsensh.SetLineStyle(5);
+    gsensh.SetLineColor(14);
+    gsensh.SetLineWidth(2);
+    gsensh.DrawClone("C")->SetBit(kCanDelete);
 
     // Display dN/dE*E^2 for conveinience
@@ -1015 +1073 @@
         spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
         spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *e*e);
+    }
+
+    for (int i=0; i<gsensl.GetN(); i++)
+    {
+        const Double_t e = gsensl.GetX()[i]*1e-3;
+
+        gsensl.GetY()[i] *= e*e;
+        gsensh.GetY()[i] *= e*e;
     }
 
@@ -1031 +1097 @@
     static_cast<const TAttLine&>(fc).Copy(fa2);
     fa2.DrawCopy("same");
+
+    gsensl.DrawClone("C")->SetBit(kCanDelete);
+    gsensh.DrawClone("C")->SetBit(kCanDelete);
 
     // Fit Spectrum
@@ -1317 +1386 @@
     hist.InitTitle(";Slope;Disp-Dist [\\circ];");
     hist.SetDrawOption("colz profx");
+}
+
+// --------------------------------------------------------------------------
+//
+// Setup write to write:
+//     container         tree       optional?
+//  --------------     ----------  -----------
+//   "MHillas"      to  "Events"
+//   "MHillasSrc"   to  "Events"
+//   "Hadronness"   to  "Events"       yes
+//   "MEnergyEst"   to  "Events"       yes
+//   "DataType"     to  "Events"
+//
+void MJSpectrum::SetupWriter(MWriteRootFile *write/*, const char *name*/) const
+{
+    if (!write)
+        return;
+
+    //write->SetName(name);
+    write->AddContainer("MHillas",        "Events");
+    write->AddContainer("MHillasSrc",     "Events");
+    write->AddContainer("MHillasExt",     "Events");
+    //write->AddContainer("MPointingPos",   "Events");
+    write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
+    write->AddContainer("MImagePar",      "Events", kFALSE);
+    write->AddContainer("MNewImagePar",   "Events", kFALSE);
+    write->AddContainer("MNewImagePar2",  "Events", kFALSE);
+    write->AddContainer("Hadronness",     "Events", kFALSE);
+    write->AddContainer("MSrcPosCam",     "Events", kFALSE);
+    write->AddContainer("MSrcPosAnti",    "Events", kFALSE);
+    write->AddContainer("ThetaSquared",   "Events", kFALSE);
+    write->AddContainer("OpticalAxis",    "Events", kFALSE);
+    write->AddContainer("Disp",           "Events", kFALSE);
+    write->AddContainer("Ghostbuster",    "Events", kFALSE);
+    write->AddContainer("MEnergyEst",     "Events", kFALSE);
+    write->AddContainer("MTime",          "Events", kFALSE);
+    write->AddContainer("MMcEvt",         "Events", kFALSE);
+    write->AddContainer("MWeight",        "Events");
+    write->AddContainer("DataType",       "Events");
+    write->AddContainer("RunNumber",      "Events");
+    write->AddContainer("EvtNumber",      "Events");
 }
 
@@ -1598 +1708 @@
     taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
 
+    MWriteRootFile write(GetPathOut());
+    SetupWriter(&write);
+
+    MParameterI par("DataType");
+    plist.AddToList(&par);
+    par.SetVal(2);
+
     tlist2.AddToList(&read);
     // If no weighting should be applied but the events should
@@ -1621 +1738 @@
     tlist2.AddToList(fCut3);
     tlist2.AddToList(&taskenv2);
+    if (!GetPathOut().IsNull())
+        tlist2.AddToList(&write);
     tlist2.AddToList(&fill31);
     tlist2.AddToList(&fill4);
@@ -1670 +1789 @@
 
     // Spectrum fitted (convert res[1] from TeV to GeV)
-    TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
+    TF1 flx("flx", Form("%e*pow(x/500, %f)", res[1]/500, res[0]));
 
     // Number of events this spectrum would produce per s and m^2
@@ -1701 +1820 @@
     cont.Add((TObject*)GetEnv()); // const_cast
     cont.Add((TObject*)&set);     // const_cast
+    cont.Add(plist.FindObject("MAlphaFitter"));
     cont.Add(&area0);
     cont.Add(&area1);
@@ -1708 +1828 @@
         cont.Add(fDisplay);
 
-    if (!WriteContainer(cont, "", "RECREATE"))
+    if (!WriteContainer(cont, "", GetPathOut().IsNull()?"RECREATE":"UPDATE"))
     {
         *fLog << err << GetDescriptor() << ": Writing result failed." << endl;

trunk/MagicSoft/Mars/mjobs/MJSpectrum.h

@@ -19 +19 @@
 class MAlphaFitter;
 class MStatusArray;
+class MWriteRootFile;
 class MHCollectionArea;
 class MMcSpectrumWeight;
@@ -36 +37 @@
 
     Bool_t fForceTheta;
+    Bool_t fForceRunTime;
 
     // Setup Histograms
@@ -50 +52 @@
     Bool_t   Refill(MParList &plist, TH1D &h) /*const*/;
     Bool_t   InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
+
+    // Write output
+    void     SetupWriter(MWriteRootFile *write/*, const char *name*/) const;
 
     // Display Output
@@ -67 +72 @@
     Bool_t Process(const MDataSet &set);
 
-    void ForceTheta(Bool_t b=kTRUE) { fForceTheta=b; }
+    void ForceTheta(Bool_t b=kTRUE)   { fForceTheta=b; }
+    void ForceRunTime(Bool_t b=kTRUE) { fForceRunTime=b; }
 
     void SetEnergyEstimator(const MTask *task);

trunk/MagicSoft/Mars/mjobs/MJob.cc

@@ -431 +431 @@
     title += fName;
 
+    // In case the update-option is selected check whether
+    // the file is already open
+    if (TString(option).Contains("update", TString::kIgnoreCase))
+    {
+        TFile *file = dynamic_cast<TFile*>(gROOT->GetListOfFiles()->FindObject(oname));
+        if (file && file->IsOpen() && !file->IsZombie())
+        {
+            *fLog << inf << "Open file found." << endl;
+            file->cd();
+            return WriteContainer(cont);
+        }
+    }
+
+    // Open a new file with the defined option for writing
     TFile file(oname, option, title, compr);
-    if (!file.IsOpen())
+    if (!file.IsOpen() || file.IsZombie())
     {
         *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;