Changeset 6978

Timestamp:

04/25/05 15:37:35 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (1 diff)
• mbase/MStatusDisplay.cc (modified) (2 diffs)
• mbase/MStatusDisplay.h (modified) (2 diffs)
• mhbase/MFillH.cc (modified) (1 diff)
• mhbase/MH.cc (modified) (1 diff)
• mhbase/MH.h (modified) (1 diff)
• mhbase/MH3.h (modified) (1 diff)
• mhflux/MHFalseSource.h (modified) (1 diff)
• mjobs/MJCut.cc (modified) (13 diffs)
• mjobs/MJCut.h (modified) (2 diffs)
• mjobs/MJOptimize.cc (modified) (1 diff)
• mjobs/MJSpectrum.cc (modified) (18 diffs)
• mjobs/MJSpectrum.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -21 +21 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2005/04/25 Thomas Bretz
+
+   * ganymed.cc:
+     - changed policy of writing the resulting events to the result file
+
+   * sponde.cc:
+     - added commandline option to use all monte carlos
+     - added command line option to read the MCs more accurate
+
+   * sponde.rc:
+     - added
+
+   * mbase/MStatusDisplay.[h,cc]:
+     - added some code to get Tab by name
+     - fixed a typo in a status line output
+
+   * mhbase/MH.[h,cc], mhbase/MH3.[h,cc], mhflux/MHFalseSource.h,
+     mhist/MHCamEvent.[h,cc], mhist/MHCamEventRot.h,
+     mhist/MHEvent.h, mhist/MHStarMap.h, mhist/MHTriggerLvl0.[h,cc],
+     mhistmc/MHMcTriggerLvl2.[h,cc], mhvstime/MHPixVsTime.[h,cc],
+     mhvstime/MHSectorVsTime.[h,cc], mimage/MHHillas.[h,cc],
+     mimage/MHHillasExt.[h,cc], mimage/MHHillasSrc.[h,cc],
+     mimage/MHImagePar.[h,cc], mimage/MHNewImagePar.[h,cc]:
+     - changed GetHistByName to be const-qualified to be compatible
+       with FindObject
+     - added some FindObject function to call GetHistByName
+
+   * mhflux/MHAlpha.[h,cc]:
+     - changed such, that it can be forced to display the excess
+       events versus size
+
+   * mjobs/MJCut.[h,cc]:
+     - display number of excess events versus Size per default
+     - removed energy estimator
+
+   * mjobs/MJOptimize.cc:
+     - display number of excess events verss size after optimization
+
+   * mjobs/MJSpectrum.[h,cc]:
+     - implemented setting up energy estimator
+     - replaced some gLog by fLog
+     - display comparison of image parameters
+
+
 
  2005/04/22 Thomas Bretz

trunk/MagicSoft/Mars/NEWS

@@ -21 +21 @@
    - Changed the default paths for calibrated data and image files.
      (The implemented access to these files doesn't yet exist)
+
+   - ganymed (MJCut and MJOptimize) now displayes the number of 
+     excess events versus size. The energy estimation is done in
+     MJSpectrum (sponde)
 
 

trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc

@@ -178 +178 @@
 // --------------------------------------------------------------------------
 
+TGCompositeFrame *MStatusDisplay::GetTabContainer(const char *name) const
+{
+#if ROOT_VERSION_CODE < ROOT_VERSION(4,03,05)
+    if (!fTab)
+        return 0;
+
+   TGFrameElement *el;
+   TGTabElement *tab = 0;
+   TGCompositeFrame *comp = 0;
+
+   TIter next(fTab->GetList());
+   next();           // skip first container
+
+   while ((el = (TGFrameElement *) next())) {
+      el = (TGFrameElement *) next();
+      comp = (TGCompositeFrame *) el->fFrame;
+      next();
+      tab = (TGTabElement *)el->fFrame;
+      if (name == tab->GetText()->GetString()) {
+         return comp;
+      }
+   }
+
+   return 0;
+#else
+   return fTab ? fTab->GetTabContainer(name) : 0;
+#endif
+}
+
+TGTabElement *MStatusDisplay::GetTabTab(const char *name) const
+{
+#if ROOT_VERSION_CODE < ROOT_VERSION(4,03,05)
+    if (!fTab)
+        return 0;
+
+   TGFrameElement *el;
+   TGTabElement *tab = 0;
+
+   TIter next(fTab->GetList());
+   next();           // skip first container
+
+   while ((el = (TGFrameElement *) next())) {
+      next();
+      tab = (TGTabElement *)el->fFrame;
+      if (name == tab->GetText()->GetString()) {
+         return tab;
+      }
+   }
+
+   return 0;
+#else
+   return fTab ? fTab->GetTabTab(name) : 0;
+#endif
+}
+// --------------------------------------------------------------------------
+
+
 // --------------------------------------------------------------------------
 //
@@ -1341 +1398 @@
 
     case kLogAppend:
-        SetStatusLine1("Appending logg...");
+        SetStatusLine1("Appending log...");
         SetStatusLine2("");
         *fLog << inf << "Appending log... " << flush;

trunk/MagicSoft/Mars/mbase/MStatusDisplay.h

@@ -30 +30 @@
 class TGTextView;
 class TGStatusBar;
+class TGTabElement;
 class TGProgressBar;
 class TGHProgressBar;
@@ -122 +123 @@
     Bool_t HandleEvent(Event_t *event);
 
+    TGCompositeFrame *GetTabContainer(const char *name) const;
+    TGTabElement     *GetTabTab(const char *name) const;
+
     Bool_t HandleTimer(TTimer *timer=NULL);
     void UpdateTab(TGCompositeFrame *f);

trunk/MagicSoft/Mars/mhbase/MFillH.cc

@@ -543 +543 @@
      */
 
-    TVirtualPad *save = gPad;
-    if (fCanvas)
-        fCanvas->cd();
+//    TVirtualPad *save = gPad;
+//    if (fCanvas)
+//        fCanvas->cd();
 
     const Bool_t rc = fH->Fill(fParContainer, fWeight?fWeight->GetVal():1);
     fH->SetNumExecutions(GetNumExecutions()+1);
 
-    if (save && fCanvas)
-        save->cd();
+//    if (save && fCanvas)
+//        save->cd();
     return rc;
 } 

trunk/MagicSoft/Mars/mhbase/MH.cc

@@ -117 +117 @@
 // a pointer to a root histogram from the MH-derived class.
 //
-TH1 *MH::GetHistByName(const TString name)
+TH1 *MH::GetHistByName(const TString name) const
 {
     *fLog << warn << GetDescriptor() << ": GetHistByName not overloaded! Can't be used!" << endl;

trunk/MagicSoft/Mars/mhbase/MH.h

@@ -49 +49 @@
     virtual TString GetDataMember() const { return ""; }
 
-    virtual TH1 *GetHistByName(const TString name);
+    virtual TH1 *GetHistByName(const TString name) const;
 
     static TCanvas *MakeDefCanvas(TString name="", const char *title="",

trunk/MagicSoft/Mars/mhbase/MH3.h

@@ -72 +72 @@
     const TH1 &GetHist() const { return *fHist; }
 
-    TH1 *GetHistByName(const TString name="") { return fHist; }
+    TH1 *GetHistByName(const TString name="") const { return fHist; }
+    TObject *FindObject(const TObject *obj) const { return 0; }
+    TObject *FindObject(const char *name) const
+    {
+        return (TObject*)GetHistByName(name);
+    }
 
     void SetColors() const;

trunk/MagicSoft/Mars/mhflux/MHFalseSource.h

@@ -59 +59 @@
     MHFalseSource(const char *name=NULL, const char *title=NULL);
 
-    TH1 *GetHistByName(const TString name) { return &fHist; }
-
     void FitSignificance(Float_t sigint=10, Float_t sigmax=75, Float_t bgmin=45, Float_t bgmax=85, Byte_t polynom=2); //*MENU*
     void FitSignificanceStd() { FitSignificance(); } //*MENU*

trunk/MagicSoft/Mars/mjobs/MJCut.cc

@@ -48 +48 @@
 #include "MPrint.h"
 #include "MContinue.h"
-#include "MEnergyEstimate.h"
+//#include "MEnergyEstimate.h"
 #include "MTaskEnv.h"
 #include "MSrcPosCalc.h"
@@ -56 +56 @@
 
 #include "../mhflux/MAlphaFitter.h"
+#include "../mhflux/MHAlpha.h"
+
 #include "MH3.h"
 #include "MBinning.h"
@@ -71 +73 @@
 // Default constructor.  Set defaults for fStoreSummary, fStoreresult,
 // fWriteOnly, fIsWobble and fFullDisplay to kFALSE and initialize
-// fEstimateEnergy and fCalcHadronness with NULL.
+// /*fEstimateEnergy and*/ fCalcHadronness with NULL.
 //
 MJCut::MJCut(const char *name, const char *title)
-    : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE),
+    : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
     fIsWobble(kFALSE), fIsMonteCarlo(kFALSE),  fFullDisplay(kFALSE), /*fSubstraction(kFALSE),*/
-    fEstimateEnergy(0), fCalcHadronness(0)
+    /*fEstimateEnergy(0),*/ fCalcHadronness(0)
 {
     fName  = name  ? name  : "MJCut";
@@ -88 +90 @@
 MJCut::~MJCut()
 {
-    if (fEstimateEnergy)
-        delete fEstimateEnergy;
+    //if (fEstimateEnergy)
+    //    delete fEstimateEnergy;
     if (fCalcHadronness)
         delete fCalcHadronness;
@@ -132 +134 @@
 // Setup a task estimating the energy. The given task is cloned.
 //
+/*
 void MJCut::SetEnergyEstimator(const MTask *task)
 {
@@ -138 +141 @@
     fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
 }
+*/
 
 // --------------------------------------------------------------------------
@@ -246 +250 @@
     TObjArray arr;
 
+    // Save all MBinnings
     TIter Next(plist);
     TObject *o=0;
@@ -252 +257 @@
             arr.Add(o);
 
+    // Save also the result, not only the setup
+    const MHAlpha *halpha = (MHAlpha*)plist.FindObject("MHAlpha");
+    if (halpha)
+        arr.Add((TObject*)(&halpha->GetAlphaFitter()));
+
     const TString fname(GetOutputFile(num));
 
+    // If requested, write to already open output file
     if (fNameResult.IsNull() && fStoreResult)
     {
@@ -413 +424 @@
         fit.SetScaleMode(MAlphaFitter::kNone);
 
+    MHAlpha halphaoff("MHAlphaOff");
+    halphaoff.ForceUsingSize();
+    plist.AddToList(&halphaoff);
+
     MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
     MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
@@ -452 +467 @@
     SetupWriter(write1, "WriteAfterCut3");
 
-
+/*
     MEnergyEstimate est;
 
     MTaskEnv taskenv1("EstimateEnergy");
     taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
-
+  */
     MTaskEnv taskenv2("CalcHadronness");
     taskenv2.SetDefault(fCalcHadronness);
@@ -495 +510 @@
     if (fIsWobble)
         tlist2.AddToList(&hcalc2);
-    tlist2.AddToList(&taskenv1);
+    //tlist2.AddToList(&taskenv1);
     tlist2.AddToList(&taskenv2);
     tlist2.AddToList(&cont0);
@@ -539 +554 @@
 
     TObjArray cont;
-    cont.Add(&fit);
     cont.Add(&cont0);
     cont.Add(&cont1);
     cont.Add(&cont2);
     cont.Add(&cont3);
-    if (taskenv1.GetTask())
-        cont.Add(taskenv1.GetTask());
+    //if (taskenv1.GetTask())
+    //    cont.Add(taskenv1.GetTask());
     if (taskenv2.GetTask())
         cont.Add(taskenv2.GetTask());
@@ -641 +655 @@
     }
     */
+    MHAlpha halphaon("MHAlpha");
+    halphaon.ForceUsingSize();
+    plist.AddToList(&halphaon);
+
     MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
     MFillH ffs2("MHFalseSource", "MHillas", "FillFS");

trunk/MagicSoft/Mars/mjobs/MJCut.h

@@ -26 +26 @@
     TString fNameOutput;
 
-    MTask *fEstimateEnergy;
+    //MTask *fEstimateEnergy;
     MTask *fCalcHadronness;
 
@@ -57 +57 @@
     void SetNameOutFile(const char *name="")      { fNameOutput=name; }
 
-    void SetEnergyEstimator(const MTask *task=0);
+    //void SetEnergyEstimator(const MTask *task=0);
     void SetHadronnessCalculator(const MTask *task=0);
 

trunk/MagicSoft/Mars/mjobs/MJOptimize.cc

@@ -914 +914 @@
     histof.SkipHistTheta();
     //histof.SkipHistEnergy();
-    histon.InitMapping(&m, 0);
-    histof.InitMapping(&m, 0);
+    histon.ForceUsingSize();
+    histof.ForceUsingSize();
+    histon.InitMapping(&m, 1);
+    histof.InitMapping(&m, 1);
 
     if (filter && filter->InheritsFrom(MFMagicCuts::Class()))

trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

@@ -68 +68 @@
 #include "MFEventSelector2.h"
 #include "MFDataMember.h"
+#include "MEnergyEstimate.h"
+#include "MTaskEnv.h"
 #include "MFillH.h"
 #include "MHillasCalc.h"
@@ -79 +81 @@
 MJSpectrum::MJSpectrum(const char *name, const char *title)
     : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
-    fRefill(kFALSE), fSimpleMode(kTRUE)
+    fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
 {
     fName  = name  ? name  : "MJSpectrum";
@@ -99 +101 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Setup a task estimating the energy. The given task is cloned.
+//
+void MJSpectrum::SetEnergyEstimator(const MTask *task)
+{
+    if (fEstimateEnergy)
+        delete fEstimateEnergy;
+    fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
+}
+
 Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
 {
@@ -125 +138 @@
 }
 
-Bool_t MJSpectrum::ReadInput(const MParList &plist)
-{
-    const TString fname = fPathIn;
-
-    *fLog << inf << "Reading from file: " << fname << endl;
-
-    TFile file(fname, "READ");
+void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
+{
+    fLog->Separator("Alpha Fitter");
+    *fLog << all;
+    fit.Print();
+
+    fLog->Separator("Used Cuts");
+    fCut0->Print();
+    fCut1->Print();
+    fCut2->Print();
+    fCut3->Print();
+
+    //gLog.Separator("Energy Estimator");
+    //fEstimateEnergy->Print();
+}
+
+Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
+{
+    *fLog << inf << "Reading from file: " << fPathIn << endl;
+
+    TFile file(fPathIn, "READ");
     if (!file.IsOpen())
     {
-        *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
-        return kFALSE;
-    }
+        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
+        return -1;
+    }
+
+    MStatusArray arr;
+    if (arr.Read()<=0)
+    {
+        *fLog << "MStatusDisplay not found in file... abort." << endl;
+        return -1;
+    }
+
+    TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",  "TH1D", "OnTime");
+    TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "MHAlpha");
+    if (!vstime || !size)
+        return -1;
+
+    vstime->Copy(h1);
+    size->Copy(h2);
+    h1.SetDirectory(0);
+    h2.SetDirectory(0);
+
+    if (fDisplay)
+        arr.DisplayIn(*fDisplay, "MHAlpha");
 
     if (!ReadTask(fCut0, "Cut0"))
-        return kFALSE;
+        return -1;
     if (!ReadTask(fCut1, "Cut1"))
-        return kFALSE;
+        return -1;
     if (!ReadTask(fCut2, "Cut2"))
-        return kFALSE;
+        return -1;
     if (!ReadTask(fCut3, "Cut3"))
-        return kFALSE;
-
-    if (!ReadTask(fEstimateEnergy, "EstimateEnergy"))
-        return kFALSE;
-
-    TObjArray arr;
+        return -1;
+
+    TObjArray arrread;
 
     TIter Next(plist);
@@ -156 +200 @@
     while ((o=Next()))
         if (o->InheritsFrom(MBinning::Class()))
-            arr.Add(o);
-
-    arr.Add(plist.FindObject("MAlphaFitter"));
-
-    return ReadContainer(arr);
-}
-
-void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
-{
-    gLog.Separator("Alpha Fitter");
-    *fLog << all;
-    fit.Print();
-
-    gLog.Separator("Used Cuts");
-    fCut0->Print();
-    fCut1->Print();
-    fCut2->Print();
-    fCut3->Print();
-
-    gLog.Separator("Energy Estimator");
-    fEstimateEnergy->Print();
-}
-
-Float_t MJSpectrum::ReadHistograms(TH1D &h1, TH1D &h2) const
-{
-    TFile file(fPathIn, "READ");
-    if (!file.IsOpen())
-    {
-        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
+            arrread.Add(o);
+
+    arrread.Add(plist.FindObject("MAlphaFitter"));
+
+    if (!ReadContainer(arrread))
         return -1;
-    }
-
-    MStatusArray *arr = (MStatusArray*)file.Get("MStatusDisplay");
-    if (!arr)
-    {
-        gLog << "MStatusDisplay not found in file... abort." << endl;
-        return -1;
-    }
-
-    TH1D *vstime = (TH1D*)arr->FindObjectInCanvas("Theta",        "TH1D", "OnTime");
-    TH1D *excen  = (TH1D*)arr->FindObjectInCanvas("ExcessEnergy", "TH1D", "MHAlpha");
-    if (!vstime || !excen)
-        return -1;
-
-    vstime->Copy(h1);
-    excen->Copy(h2);
-    h1.SetDirectory(0);
-    h2.SetDirectory(0);
-
-    if (fDisplay)
-        arr->DisplayIn(*fDisplay, "MHAlpha");
-
-    delete arr;
 
     return vstime->Integral();
@@ -303 +300 @@
 void MJSpectrum::DisplayResult(const TH2D &h2) const
 {
-    if (!fDisplay->CdCanvas("ZdDist"))
+    if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
         return;
 
@@ -342 +339 @@
     read.AddFile(fPathIn);
 
-    MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
-    MFillH fill2("MHAlpha", "MHillasSrc", "FillAlpha");
+    MEnergyEstimate est;
+    MTaskEnv taskenv1("EstimateEnergy");
+    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
+
+    MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlphaOff");
+    MFillH fill2("MHAlpha",              "MHillasSrc", "FillAlphaOn");
 
     MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
@@ -352 +353 @@
 
     tlist.AddToList(&read);
-    tlist.AddToList(fEstimateEnergy);
+    tlist.AddToList(&taskenv1);
     tlist.AddToList(&f0);
     tlist.AddToList(&f1);
@@ -389 +390 @@
     halpha->GetHEnergy().Copy(h2);
     h2.SetDirectory(0);
+
+    return kTRUE;
+}
+
+Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
+{
+    MTaskList tlist1;
+    plist.Replace(&tlist1);
+
+    MReadMarsFile readmc("OriginalMC");
+    //readmc.DisableAutoScheme();
+    set.AddFilesOn(readmc);
+    readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
+    readmc.EnableBranch("MMcEvtBasic.fEnergy");
+
+    mh1.SetLogy();
+    mh1.SetLogz();
+    mh1.SetName("ThetaE");
+
+    MFillH fill0(&mh1);
+    //fill0.SetDrawOption("projx only");
+
+    MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
+    MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
+    if (bins2 && bins3)
+    {
+        bins2->SetName("BinningThetaEY");
+        bins3->SetName("BinningThetaEX");
+    }
+    tlist1.AddToList(&readmc);
+
+    temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
+    MH3 mh3mc(temp1);
+
+    MFEventSelector2 sel1(mh3mc);
+    sel1.SetHistIsProbability();
+
+    fill0.SetFilter(&sel1);
+
+    if (!fRawMc)
+        tlist1.AddToList(&sel1);
+    tlist1.AddToList(&fill0);
+
+    MEvtLoop loop1(fName);
+    loop1.SetParList(&plist);
+    loop1.SetLogStream(fLog);
+    loop1.SetDisplay(fDisplay);
+
+    if (!SetupEnv(loop1))
+        return kFALSE;
+
+    if (!loop1.Eventloop(fMaxEvents))
+    {
+        *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
+        return kFALSE;
+    }
+
+    tlist1.PrintStatistics();
+
+    if (!loop1.GetDisplay())
+    {
+        *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
+        return kFALSE;
+    }
+
+    if (bins2 && bins3)
+    {
+        bins2->SetName("BinningEnergyEst");
+        bins3->SetName("BinningTheta");
+    }
+    return kTRUE;
+}
+
+void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
+{
+    TH1D collarea(area.GetHEnergy());
+    TH1D spectrum(excess);
+    TH1D weights;
+    hest.GetWeights(weights);
+
+    cout << "Effective On time: " << ontime << "s" << endl;
+
+    spectrum.SetDirectory(NULL);
+    spectrum.SetBit(kCanDelete);
+    spectrum.Scale(1./ontime);
+    spectrum.Divide(&collarea);
+    spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
+    spectrum.SetYTitle("N/sm^{2}");
+
+    TCanvas &c1 = fDisplay->AddTab("Spectrum");
+    c1.Divide(2,2);
+    c1.cd(1);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGridx();
+    gPad->SetGridy();
+    collarea.DrawCopy();
+
+    c1.cd(2);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGridx();
+    gPad->SetGridy();
+    spectrum.DrawCopy();
+
+    c1.cd(3);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGridx();
+    gPad->SetGridy();
+    weights.DrawCopy();
+
+    //spectrum.Divide(&weights);
+    spectrum.Multiply(&weights);
+    spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
+
+    for (int i=0; i<excess.GetNbinsX(); i++)
+    {
+        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);
+    }
+
+    c1.cd(4);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGridx();
+    gPad->SetGridy();
+    spectrum.SetXTitle("E [GeV]");
+    spectrum.SetYTitle("N/TeVsm^{2}");
+    spectrum.DrawCopy();
+
+    TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
+    f.SetParameter(0, -2.87);
+    f.SetParameter(1, 1.9e-6);
+    f.SetLineColor(kGreen);
+    spectrum.Fit(&f, "NI", "", 55, 2e4);
+    f.DrawCopy("same");
+
+    /*
+     TString str;
+     str += "(1.68#pm0.15)10^{-7}";
+     str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
+     str += "\\frac{ph}{TeVm^{2}s}";
+
+     TLatex tex;
+     tex.DrawLatex(2e2, 7e-5, str);
+     */
+}
+
+Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
+{
+    cout << name << endl;
+    TString same(name);
+    same += "Same";
+
+    TH1 *h1  = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
+    TH1 *h2  = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
+    if (!h1 || !h2)
+        return kFALSE;
+
+    TObject *obj = plist.FindObject(plot);
+    if (!obj)
+    {
+        *fLog << warn << plot << " not in parameter list... skipping." << endl;
+        return kFALSE;
+    }
+
+    TH1 *h3  = (TH1*)obj->FindObject(name);
+    if (!h3)
+    {
+        *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
+        return kFALSE;
+    }
+
+
+    const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
+    const Double_t scale = fit ? fit->GetScaleFactor() : 1;
+
+    gPad->SetBorderMode(0);
+    h2->SetLineColor(kBlack);
+    h3->SetLineColor(kBlue);
+    h2->Add(h1, -scale);
+
+    h2->Scale(1./h2->Integral());
+    h3->Scale(1./h3->Integral());
+
+    h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
+
+    h2 = h2->DrawCopy();
+    h3 = h3->DrawCopy("same");
+
+    // Don't do this on the original object!
+    h2->SetStats(kFALSE);
+    h3->SetStats(kFALSE);
+
+    return kTRUE;
+}
+
+Bool_t MJSpectrum::DisplaySize(MParList &plist) const
+{
+    *fLog << inf << "Reading from file: " << fPathIn << endl;
+
+    cout << "Opening..." << endl;
+    TFile file(fPathIn, "READ");
+    if (!file.IsOpen())
+    {
+        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
+        return kFALSE;
+    }
+
+    cout << "Reading..." << endl;
+    file.cd();
+    MStatusArray arr;
+    if (arr.Read()<=0)
+    {
+        *fLog << "MStatusDisplay not found in file... abort." << endl;
+        return kFALSE;
+    }
+
+    cout << "Searching..." << endl;
+
+    TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "MHAlpha");
+    if (!excess)
+        return kFALSE;
+
+    cout << "Displaying..." << endl;
+
+    // ------------------- Plot excess versus size ------------------- 
+
+    TCanvas &c = fDisplay->AddTab("Excess");
+    c.Divide(3,2);
+    c.cd(1);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGridx();
+    gPad->SetGridy();
+
+    excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
+    excess->SetStats(kFALSE);
+    excess->Scale(1./excess->Integral());
+    excess->DrawCopy();
+
+    TObject *o=0;
+    if ((o=plist.FindObject("ExcessSize")))
+    {
+        TH1F *histsel = (TH1F*)o->FindObject("");
+        histsel->SetStats(kFALSE);
+        histsel->Scale(1./histsel->Integral());
+        histsel->SetLineColor(kBlue);
+        histsel->SetBit(kCanDelete);
+        histsel->DrawCopy("same");
+    }
+
+    // -------------- Comparison of Image Parameters --------------
+    c.cd(2);
+    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost");
+
+    c.cd(3);
+    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
+
+    c.cd(4);
+    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost");
+
+    c.cd(5);
+    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost");
+
+    c.cd(6);
+    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost");
 
     return kTRUE;
@@ -432 +706 @@
     plist.AddToList(&fit);
 
-    if (!ReadInput(plist))
-        return kFALSE;
+    TH1D temp1, size;
+    const Float_t ontime = ReadInput(plist, temp1, size);
+    if (ontime<0)
+    {
+        *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
+        return kFALSE;
+    }
 
     PrintSetup(fit);
-
-    TH1D temp1, excess;
-    const Float_t ontime = ReadHistograms(temp1, excess);
-    if (ontime<0)
-    {
-        *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
-        return kFALSE;
-    }
+    bins3.SetEdges(temp1, 'x');
 
     TH1D temp2(temp1);
@@ -452 +724 @@
         return kFALSE;
 
-    if (fRefill && !Refill(plist, temp2))
+    TH1D excess;
+    if (!Refill(plist, excess))
         return kFALSE;
 
@@ -460 +733 @@
     {
         hist.UseCurrentStyle();
-        MH::SetBinning(&hist, temp1.GetXaxis(), excess.GetXaxis());
+        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
         if (!ReadOrigMCDistribution(set, hist))
             return kFALSE;
 
-        for (int y=0; y<hist.GetNbinsY(); y++)
-            for (int x=0; x<hist.GetNbinsX(); x++)
-                hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
+        if (!fRawMc)
+        {
+            for (int y=0; y<hist.GetNbinsY(); y++)
+                for (int x=0; x<hist.GetNbinsX(); x++)
+                    hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
+        }
     }
     else
-    {
-        MTaskList tlist1;
-        plist.Replace(&tlist1);
-
-        MReadMarsFile readmc("OriginalMC");
-        //readmc.DisableAutoScheme();
-        set.AddFilesOn(readmc);
-        readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
-        readmc.EnableBranch("MMcEvtBasic.fEnergy");
-
-        mh1.SetLogy();
-        mh1.SetLogz();
-        mh1.SetName("ThetaE");
-
-        MFillH fill0(&mh1);
-        //fill0.SetDrawOption("projx only");
-
-        MBinning binsx(bins3, "BinningThetaEX");
-        MBinning binsy(bins2, "BinningThetaEY");
-        plist.AddToList(&binsx);
-        plist.AddToList(&binsy);
-        tlist1.AddToList(&readmc);
-
-        temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
-        MH3 mh3mc(temp1);
-
-        MFEventSelector2 sel1(mh3mc);
-        sel1.SetHistIsProbability();
-
-        fill0.SetFilter(&sel1);
-
-        tlist1.AddToList(&sel1);
-        tlist1.AddToList(&fill0);
-
-        MEvtLoop loop1(fName);
-        loop1.SetParList(&plist);
-        loop1.SetLogStream(fLog);
-        loop1.SetDisplay(fDisplay);
-
-        if (!SetupEnv(loop1))
+        if (!IntermediateLoop(plist, mh1, temp1, set))
             return kFALSE;
-
-        if (!loop1.Eventloop(fMaxEvents))
-        {
-            *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
-            return kFALSE;
-        }
-
-        tlist1.PrintStatistics();
-
-        if (!loop1.GetDisplay())
-        {
-            *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
-            return kFALSE;
-        }
-    }
 
     DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
@@ -572 +794 @@
     MFillH fill4(&hest, "", "FillEnergyEst");
 
+    MH3 hsize("MHillas.fSize");
+    //MH3 henergy("MEnergyEst.fVal");
+    hsize.SetName("ExcessSize");
+    //henergy.SetName("EnergyEst");
+    MBinning bins(size, "BinningExcessSize");
+    plist.AddToList(&hsize);
+    //plist.AddToList(&henergy);
+    plist.AddToList(&bins);
+
     MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
     MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
@@ -579 +810 @@
     MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
     MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
+    MFillH fill8a("ExcessSize     [MH3]",           "",             "FillExcessSize");
+    //MFillH fill9a("EnergyEst      [MH3]",           "",             "FillExcessEEst");
     fill1a.SetNameTab("PreCut");
     fill2a.SetNameTab("PostCut");
@@ -586 +819 @@
     fill6a.SetNameTab("ImgPar");
     fill7a.SetNameTab("NewPar");
+    fill8a.SetBit(MFillH::kDoNotDisplay);
+    //fill9a.SetBit(MFillH::kDoNotDisplay);
+
+    MEnergyEstimate est;
+    MTaskEnv taskenv1("EstimateEnergy");
+    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
 
     tlist2.AddToList(&read);
-    tlist2.AddToList(&contsel);
+    if (!fRawMc)
+        tlist2.AddToList(&contsel);
     tlist2.AddToList(&calc);
     tlist2.AddToList(&hcalc1);
@@ -597 +837 @@
     tlist2.AddToList(fCut2);
     tlist2.AddToList(fCut3);
-    tlist2.AddToList(fEstimateEnergy);
+    tlist2.AddToList(&taskenv1);
     tlist2.AddToList(&fill3);
     tlist2.AddToList(&fill4);
@@ -606 +846 @@
     tlist2.AddToList(&fill6a);
     tlist2.AddToList(&fill7a);
+    tlist2.AddToList(&fill8a);
+    //tlist2.AddToList(&fill9a);
 
     MEvtLoop loop2(fName);
@@ -631 +873 @@
     // -------------------------- Spectrum ----------------------------
 
-    TH1D collarea(area.GetHEnergy());
-    TH1D weights;
-    hest.GetWeights(weights);
-
-    cout << "Effective On time: " << ontime << "s" << endl;
-
-    excess.SetDirectory(NULL);
-    excess.SetBit(kCanDelete);
-    excess.Scale(1./ontime);
-    excess.Divide(&collarea);
-    excess.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
-    excess.SetYTitle("N/sm^{2}");
-
-    TCanvas &c = fDisplay->AddTab("Spectrum");
-    c.Divide(2,2);
-    c.cd(1);
-    gPad->SetBorderMode(0);
-    gPad->SetLogx();
-    gPad->SetLogy();
-    gPad->SetGridx();
-    gPad->SetGridy();
-    collarea.DrawCopy();
-
-    c.cd(2);
-    gPad->SetBorderMode(0);
-    gPad->SetLogx();
-    gPad->SetLogy();
-    gPad->SetGridx();
-    gPad->SetGridy();
-    excess.DrawCopy();
-
-    c.cd(3);
-    gPad->SetBorderMode(0);
-    gPad->SetLogx();
-    gPad->SetLogy();
-    gPad->SetGridx();
-    gPad->SetGridy();
-    weights.DrawCopy();
-
-    excess.Divide(&weights);
-    //excess.Multiply(&weights);
-    excess.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
-
-    for (int i=0; i<excess.GetNbinsX(); i++)
-    {
-        excess.SetBinContent(i+1, excess.GetBinContent(i+1)/excess.GetBinWidth(i+1)*1000);
-        excess.SetBinError(i+1,   excess.GetBinError(i+1)/  excess.GetBinWidth(i+1)*1000);
-    }
-
-    excess.SetYTitle("N/TeVsm^{2}");
-
-    c.cd(4);
-    gPad->SetBorderMode(0);
-    gPad->SetLogx();
-    gPad->SetLogy();
-    gPad->SetGridx();
-    gPad->SetGridy();
-    excess.DrawCopy();
-
-    TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
-    f.SetParameter(0, -2.87);
-    f.SetParameter(1, 1.9e-6);
-    f.SetLineColor(kGreen);
-    excess.Fit(&f, "NI", "", 55, 2e4);
-    f.DrawCopy("same");
+    DisplaySpectrum(area, excess, hest, ontime);
+    DisplaySize(plist);
 
     if (!fPathOut.IsNull())
         fDisplay->SaveAsRoot(fPathOut);
 
-    /*
-     TString str;
-     str += "(1.68#pm0.15)10^{-7}";
-     str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
-     str += "\\frac{ph}{TeVm^{2}s}";
-
-     TLatex tex;
-     tex.DrawLatex(2e2, 7e-5, str);
-     */
-
     return kTRUE;
 }

trunk/MagicSoft/Mars/mjobs/MJSpectrum.h

@@ -14 +14 @@
 class MParList;
 class MDataSet;
+class MHEnergyEst;
 class MAlphaFitter;
+class MStatusArray;
+class MHCollectionArea;
 
 class MJSpectrum : public MJob
@@ -27 +30 @@
     Bool_t fRefill;
     Bool_t fSimpleMode;
+    Bool_t fRawMc;
 
+    // Read Input
     Bool_t  ReadTask(MTask* &task, const char *name) const;
-    Bool_t  ReadInput(const MParList &plist);
-    void    PrintSetup(const MAlphaFitter &fit) const;
-    Float_t ReadHistograms(TH1D &h1, TH1D &h2) const;
+    Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size);
     Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const;
     Bool_t  GetThetaDistribution(TH1D &temp1, TH1D &temp2) const;
+    Bool_t  Refill(MParList &plist, TH1D &h) const;
+
+    // Display Output
+    void    PrintSetup(const MAlphaFitter &fit) const;
     void    DisplayResult(const TH2D &mh1) const;
-    Bool_t  Refill(MParList &plist, TH1D &h) const;
+    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set) const;
+    void    DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
+    Bool_t  DisplaySize(MParList &plist) const;
+    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const;
 
 public:
@@ -45 +55 @@
     void EnableRefilling(Bool_t b=kTRUE)  { fRefill=b; }
     void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; }
+    void EnableRawMc(Bool_t b=kTRUE)      { fRawMc=b; }
+
+    void SetEnergyEstimator(const MTask *task);
 
     ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum