Changeset 8680 for trunk/MagicSoft

Timestamp:

08/20/07 11:04:58 (17 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• .rootrc (modified) (1 diff)
• Changelog (modified) (1 diff)
• NEWS (modified) (3 diffs)
• mbase/MStatusDisplay.cc (modified) (2 diffs)
• mbase/MStatusDisplay.h (modified) (1 diff)
• mhflux/MMcSpectrumWeight.cc (modified) (18 diffs)
• mhflux/MMcSpectrumWeight.h (modified) (3 diffs)
• mjobs/MJSpectrum.cc (modified) (38 diffs)
• mjobs/MJSpectrum.h (modified) (2 diffs)
• sponde.rc (modified) (1 diff)

trunk/MagicSoft/Mars/.rootrc

@@ -45 +45 @@
 Root.Html.SearchEngine:      http://magic.astro.uni-wuerzburg.de/mars/search.html
 Root.Html.ViewCVS:           http://www.astro.uni-wuerzburg.de/cgi-bin/viewcvs.cgi/
+
+#############################################################################
+#                                                                           #
+# Some default changes                                                      #
+#                                                                           #
+#############################################################################
+
+#Canvas.HighLightColor:  3
+#Canvas.MoveOpaque:      0
+#Canvas.ResizeOpaque:    0
+Canvas.ShowEventStatus: Yes
+#Canvas.ShowToolBar:     Yes
+#Canvas.ShowEditor:      Yes
+#Canvas.AutoExec:        No

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2007/08/20 Thomas Bretz
+
+   * .rootrc:
+     - added some comments about defaults
+     - set the ShowEventStatus for the canvases to yes as default
+
+   * sponde.cc:
+     - removed the refill option (it was just a dummy)
+     - removed the accurate option. It didn't give more accurate 
+       results at all
+
+   * mbase/MStatusDisplay.[h,cc]:
+     - added an update option to SetProgressBarPosition
+
+   * mhflux/MMcSpectrumWeight.[h,cc]:
+     - allow to give a container name to GetFormula*
+     - changed default for spectral slope from -9 to -99
+     - allow to give integral range to GetSpec*Integral
+     - added a new member function CompeleteEnergySpectrum which completes
+       a simulated spectrum starting at an energy fEnergyMin down to
+       an energy emin.
+     - added two options ("new" and "old") to Print function
+     - do not stop anymore if lower energy boundary changes
+
+   * mjobs/MJSpectrum.[h,cc]:
+     - cleaned the code in general
+     - removed fRefill (was not used in the code at all)
+     - added MJSpectrum to global ListOfCleanups to handle
+       the display more properly
+     - removed reading of the first (it was the second!) 
+       MMcCorsikaRunHeader. It is now read for each file individually
+     - The read monte carlo events are now weighted with the mc
+       production area (events per area)
+     - incomplete (to lower energies) spectra are completed
+     - removed accurate mode, it was not more accurate
+     - we fit the spectrum now from the first to the last bin
+     - for comparison crab and 1553 are plotted
+     - changed the processing such that first the MCs are processed
+       and then the spectrum is refilled
+     - now the MC distribution from OriginalMC is read only once
+     - added new tab showing the basic event distribution
+
+
 
  2007/08/19 Thomas Bretz

trunk/MagicSoft/Mars/NEWS

@@ -223 +223 @@
      give any improvement in accuracy, only decreased execution speed.
 
+   - sponde: the so called "simple"-mode has been removed. It didn't
+     give any improvement in simple.
+
+   - sponde: the so called "refill"-mode has been removed. It was anyhow
+     not implemented.
+
    - sponde: now checks whether the theta distribution of your on-data
      and the theta-distribution of your Monte Carlo sample (after
@@ -236 +242 @@
      different maximum impact parameters.
 
+   - sponde: If MC files with different lower energy limits are used
+     the primary MC spectrum is artificially completed down to the
+     lowest energy used at all. WARNING: that this gives correct 
+     collection areas ONLY if none of the events in this region would
+     survive your cuts at all.
+
    - sponde: the output file now contains more information about
      the spectrum (eg. the full 2D collection area histogram).
@@ -248 +260 @@
      The same information you get from the error bars of the weighted
      histograms, but this is less intuitive.
+
+   - sponde: The spectrum plots now show the crab- and 1553-spectrum
+     for comparison. It is not meant to show these curves in
+     publications, they are only for production.
+
+   - sponde: The OriginalMC tree with the events produced by corsika 
+     is now processed only once
 
 

trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc

@@ -602 +602 @@
 // is assumed to be (0,1)
 //
-void MStatusDisplay::SetProgressBarPosition(Float_t p)
+void MStatusDisplay::SetProgressBarPosition(Float_t p, Bool_t upd)
 {
     fBar->SetPosition(p);
+    if (upd)
+        gClient->ProcessEventsFor(fBar);
 }
 
@@ -685 +687 @@
     // p==NULL means: Take gClient->GetRoot() if not in batch mode
     // see TGWindow::TGWindow()
+
+    // Make sure that the display is removed via RecursiveRemove
+    // from whereever possible.
+    SetBit(kMustCleanup);
 
     //

trunk/MagicSoft/Mars/mbase/MStatusDisplay.h

@@ -168 +168 @@
      void SetUpdateTime(Long_t t);
 
-     void SetProgressBarPosition(Float_t p);
+     void SetProgressBarPosition(Float_t p, Bool_t upd=kFALSE);
      TGProgressBar *GetBar() const { return (TGProgressBar*)fBar; }
 

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

@@ -19 +19 @@
 !   Author(s): Marcos Lopez 10/2003 <mailto:marcos@gae.ucm.es>
 !
-!   Copyright: MAGIC Software Development, 2000-2006
+!   Copyright: MAGIC Software Development, 2000-2007
 !
 !
@@ -40 +40 @@
 //  Method:
 //  -------
-//                                 
+//
 //   - Corsika spectrun: dN/dE = A * E^(a)
 //     with a = fOldSlope, and A = N/integral{E*de} from ELowLim to EUppLim
@@ -72 +72 @@
 //    MMcEvt
 //    MMcCorsikaRunHeader
+//    [MPointingPos]
+//    [MHillas]
 //
 //  Output Container:
@@ -81 +83 @@
 #include <TF1.h>
 #include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
 #include <TSpline.h>
 
@@ -109 +113 @@
     fNameMcEvt   = "MMcEvt";
 
-    fNewSlope    = -9;
-    fOldSlope    = -9;
+    fNewSlope    = -99;
+    fOldSlope    = -99;
 
     fEnergyMin   = -1;
@@ -210 +214 @@
 // The slope is returned as %.3f
 //
-TString MMcSpectrumWeight::GetFormulaSpecOld() const
-{
-    return Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fOldSlope);
+TString MMcSpectrumWeight::GetFormulaSpecOld(const char *name) const
+{
+    return Form("pow(%s.fEnergy, %.3f)", name, fOldSlope);
 }
 
@@ -225 +229 @@
 // unchanged.
 //
-TString MMcSpectrumWeight::GetFormulaSpecNew() const
-{
-    TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula.Data();
+TString MMcSpectrumWeight::GetFormulaSpecNew(const char *name) const
+{
+    TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", name, fNewSlope) : fFormula.Data();
     if (!fFormula.IsNull())
-        str.ReplaceAll("X", Form("(%s.fEnergy)", fNameMcEvt.Data()));
+        str.ReplaceAll("X", Form("(%s.fEnergy)", name));
 
     return str;
@@ -258 +262 @@
 // Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600))
 //
-TString MMcSpectrumWeight::GetFormulaWeights() const
+TString MMcSpectrumWeight::GetFormulaWeights(const char *name) const
 {
     if (GetFormulaSpecOld()==GetFormulaSpecNew())
@@ -268 +272 @@
     const Double_t norm = fNorm*iold/inew;
 
-    return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data());
+    return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew(name).Data(), GetFormulaSpecOld(name).Data());
 }
 
@@ -276 +280 @@
 // GetFormulaSpecNewX() describing the destination spectrum
 //
-Double_t MMcSpectrumWeight::GetSpecNewIntegral() const
+Double_t MMcSpectrumWeight::GetSpecNewIntegral(Double_t emin, Double_t emax) const
 {
     TF1 funcnew("Dummy", GetFormulaSpecNewX());
-    return funcnew.Integral(fEnergyMin, fEnergyMax);
+    return funcnew.Integral(emin, emax);
 }
 
@@ -287 +291 @@
 // GetFormulaSpecOldX() describing the simulated spectrum
 //
-Double_t MMcSpectrumWeight::GetSpecOldIntegral() const
+Double_t MMcSpectrumWeight::GetSpecOldIntegral(Double_t emin, Double_t emax) const
 {
     TF1 funcold("Dummy", GetFormulaSpecOldX());
-    return funcold.Integral(fEnergyMin, fEnergyMax);
+    return funcold.Integral(emin, emax);
 }
 
@@ -357 +361 @@
     if (fEnergyMax>fEnergyMin && !fAllowChange)
     {
+        if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10)
+        {
+            *fLog << err;
+            *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change (";
+            *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl;
+            return kFALSE;
+        }
+
         if (TMath::Abs(fOldSlope-rh.GetSlopeSpec())>1e-10)
         {
@@ -364 +376 @@
             return kFALSE;
         }
-        if (TMath::Abs(fEnergyMin-rh.GetELowLim())>1e-10)
-        {
-            *fLog << err;
-            *fLog << "ERROR - The minimum simulated Monte Carlo energy is not allowed to change ";
-            *fLog << "(" << fEnergyMin << " --> " << rh.GetELowLim() << ")... abort." << endl;
-            return kFALSE;
-        }
-        if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10)
-        {
-            *fLog << err;
-            *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change (";
-            *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl;
-            return kFALSE;
-        }
-        return kTRUE;
+
+        // No change happened
+        if (TMath::Abs(fEnergyMin-rh.GetELowLim())<=1e-10)
+            return kTRUE;
+
+        // The lower energy limit has changed
+        *fLog << warn;
+        *fLog << "The minimum simulated Monte Carlo energy has changed from ";
+        *fLog << fEnergyMin << "GeV to " << rh.GetELowLim() << "GeV." << endl;
+        fEnergyMin = rh.GetELowLim();
     }
 
@@ -385 +392 @@
     fEnergyMax = rh.GetEUppLim();
 
-    if (fNewSlope==-9)
-    {
-        *fLog << inf << "The new slope of the power law is undefined (-9)... no weighting applied." << endl;
+    if (fNewSlope==-99)
+    {
+        *fLog << inf << "A new slope for the power law has not been defined... no weighting applied." << endl;
         fNewSlope = fOldSlope;
     }
@@ -403 +410 @@
 // ---------------------------------------------------------------------------
 //
+// completes a simulated spectrum starting at an energy fEnergyMin down to
+// an energy emin.
+//
+// It is assumed that the contents of MMcSpectrumWeight for the new spectrum
+// correctly describe the spectrum within the histogram, and fEnergyMin
+// and fEnergyMax correctly describe the range.
+//
+// In the 1D case it is assumed that the x-axis is a zenith angle binning.
+// In the 2D case the x-axis is assumed to be zenith angle, the y-axis
+// to be energy.
+//
+void MMcSpectrumWeight::CompleteEnergySpectrum(TH1 &h, Double_t emin) const
+{
+    if (h.InheritsFrom(TH3::Class()))
+    {
+        *fLog << "ERROR - MMcSpctrumWeight::CompleteEnergySpectrum doesn't support TH3." << endl;
+        return;
+    }
+
+    if (fEnergyMin <= emin)
+        return;
+
+    const Double_t norm = GetSpecNewIntegral();
+
+    if (!h.InheritsFrom(TH2::Class()))
+    {
+        h.Scale(GetSpecNewIntegral(emin, fEnergyMax)/norm);
+        return;
+    }
+
+    const TAxis &axey = *h.GetYaxis();
+
+    const Int_t first = axey.FindFixBin(emin);
+    const Int_t last  = axey.FindFixBin(fEnergyMin);
+
+    for (int x=1; x<=h.GetNbinsX(); x++)
+    {
+         const Double_t f = h.Integral(x, x, -1, 9999)/norm;
+
+        for (int y=first; y<=last; y++)
+        {
+            const Double_t lo = axey.GetBinLowEdge(y)  <emin       ? emin       : axey.GetBinLowEdge(y);
+            const Double_t hi = axey.GetBinLowEdge(y+1)>fEnergyMin ? fEnergyMin : axey.GetBinLowEdge(y+1);
+
+            h.AddBinContent(h.GetBin(x, y), f*GetSpecNewIntegral(lo, hi));
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+//
 // The current contants are printed
 //
 void MMcSpectrumWeight::Print(Option_t *o) const
 {
+    const TString opt(o);
+
+    const Bool_t hasnew = opt.Contains("new") || opt.IsNull();
+    const Bool_t hasold = opt.Contains("old") || opt.IsNull();
+
     *fLog << all << GetDescriptor() << endl;
-    *fLog << " Simulated energy range:   " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
-    *fLog << " Simulated spectral slope: " << fOldSlope << endl;
-    *fLog << " New spectral slope:       " << fNewSlope << endl;
+
+    if (hasold)
+    {
+        *fLog << " Simulated energy range:   " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
+        *fLog << " Simulated spectral slope: ";
+        if (fOldSlope==-99)
+            *fLog << "undefined" << endl;
+        else
+            *fLog << fOldSlope << endl;
+    }
+    if (hasnew)
+    {
+        *fLog << " New spectral slope:       ";
+        if (fNewSlope==-99)
+            *fLog << "undefined" << endl;
+        else
+            *fLog << fNewSlope << endl;
+    }
     *fLog << " Additional user norm.:    " << fNorm << endl;
     *fLog << " Spectra are normalized:   " << (fNormEnergy<0?"by integral":Form("at %.1fGeV", fNormEnergy)) << endl;
-    *fLog << " Old Spectrum:     " << GetFormulaSpecOldX() << "   (I=" << GetSpecOldIntegral() << ")" << endl;
-    *fLog << " New Spectrum:     " << GetFormulaSpecNewX() << "   (I=" << GetSpecNewIntegral() << ")" << endl;
+    if (hasold)
+    {
+        *fLog << " Old Spectrum: " << GetFormulaSpecOldX();
+        if (fEnergyMin>=0 && fEnergyMax>0)
+            *fLog << "   (I=" << GetSpecOldIntegral() << ")";
+        *fLog << endl;
+    }
+    if (hasnew)
+    {
+        *fLog << " New Spectrum: " << GetFormulaSpecNewX();
+        if (fEnergyMin>=0 && fEnergyMax>0)
+            *fLog << "   (I=" << GetSpecNewIntegral() << ")";
+        *fLog << endl;
+    }
     if (fFunc)
-        *fLog << " Weight function:  " << fFunc->GetTitle()   << endl;
+        *fLog << " Weight func:  " << fFunc->GetTitle()   << endl;
 }
 
@@ -436 +526 @@
 }
 
+/*
+ * This could be used to improve the Zd-weighting within a bin.
+ * Another option is to use more bins, or both.
+ * Note that it seems unnecessary, because the shape within the
+ * theta-bins should be similar in data and Monte Carlo... hopefully.
+ *
+void MMcSpectrumWeight::InitZdWeights()
+{
+    TH2D w(*fWeightsZd);
+
+    for (int i=1; i<=w.GetNbinsX(); i++)
+    {
+         const Double_t tmin = w.GetBinLowEdge(i)  *TMath::DegToRad();
+         const Double_t tmax = w.GetBinLowEdge(i+1)*TMath::DegToRad();
+
+         const Double_t wdth  = tmax-tmin;
+         const Double_t integ = cos(tmin)-cos(tmax);
+
+         w.SetBinContent(i, w.GetBinContent(i)*wdth/integ);
+    }
+
+    //  const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
+    //  const Double_t theta = fPointing->GetZd()*TMath::DegToRad();
+    //  w = sin(theta)*w.GetBinContent(i);
+}
+*/
+
 // ----------------------------------------------------------------------------
 //
@@ -447 +564 @@
     if (fWeightsZd)
     {
-        /*
-         const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
-
-         Double_t tmin = fWeightsZd->GetXaxis()->GetBinLowEdge(i)  *TMath::DegToRad();
-         Double_t tmax = fWeightsZd->GetXaxis()->GetBinLowEdge(i+1)*TMath::DegToRad();
-         Double_t wdth = fWeightsZd->GetXaxis()->GetBinWidth(i)    *TMath::DegToRad();
-
-         Double_t cont = fWeightsZd->GetBinContent(i);
-
-         Double_t integ = cos(tmin)-cos(tmax);
-
-         w = sin(tmax)*cont/integ*wdth;
-         */
-
         const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
         w = fWeightsZd->GetBinContent(i);
@@ -472 +575 @@
 
     const Double_t e = fMcEvt->GetEnergy();
+
     fWeight->SetVal(fFunc->Eval(e)*w);
 

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h

@@ -77 +77 @@
 
     // Getter
-    TString GetFormulaSpecOld() const;
-    TString GetFormulaSpecNew() const;
-    TString GetFormulaWeights() const;
+    TString GetFormulaSpecOld(const char *name) const;
+    TString GetFormulaSpecNew(const char *name) const;
+    TString GetFormulaWeights(const char *name) const;
+    TString GetFormulaSpecOld() const { return GetFormulaSpecOld(fNameMcEvt); }
+    TString GetFormulaSpecNew() const { return GetFormulaSpecNew(fNameMcEvt); }
+    TString GetFormulaWeights() const { return GetFormulaWeights(fNameMcEvt); }
 
     TString GetFormulaSpecOldX() const { return ReplaceX(GetFormulaSpecOld()); }
@@ -85 +88 @@
     TString GetFormulaWeightsX() const { return ReplaceX(GetFormulaWeights()); }
 
-    Double_t GetSpecNewIntegral() const;
-    Double_t GetSpecOldIntegral() const;
+    Double_t GetSpecNewIntegral(Double_t emin, Double_t emax) const;
+    Double_t GetSpecOldIntegral(Double_t emin, Double_t emax) const;
+
+    Double_t GetSpecNewIntegral() const { return GetSpecNewIntegral(fEnergyMin, fEnergyMax); }
+    Double_t GetSpecOldIntegral() const { return GetSpecOldIntegral(fEnergyMin, fEnergyMax); }
 
     Double_t CalcSpecNew(Double_t e) const;
@@ -93 +99 @@
     Double_t GetEnergyMin() const { return fEnergyMin; }
     Double_t GetEnergyMax() const { return fEnergyMax; }
+
+    // Functions
+    void CompleteEnergySpectrum(TH1 &h, Double_t emin) const;
 
     // TObject

trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

@@ -18 +18 @@
 !   Author(s): Thomas Bretz, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2005
+!   Copyright: MAGIC Software Development, 2000-2007
 !
 !
@@ -62 +62 @@
 #include "../mhflux/MHAlpha.h"
 #include "../mhflux/MHCollectionArea.h"
-//#include "../mhflux/MHThreshold.h"
 #include "../mhflux/MHEnergyEst.h"
 #include "../mhflux/MMcSpectrumWeight.h"
@@ -80 +79 @@
 #include "MFillH.h"
 #include "MHillasCalc.h"
-//#include "MSrcPosCalc.h"
 #include "MContinue.h"
 
@@ -89 +87 @@
 MJSpectrum::MJSpectrum(const char *name, const char *title)
     : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0),
-    fEstimateEnergy(0), fCalcHadronness(0), fRefill(kFALSE),
-    /*fSimpleMode(kTRUE),*/ fForceTheta(kFALSE),
-    fRawMc(kFALSE), fNoThetaWeights(kFALSE)
+    fEstimateEnergy(0), fCalcHadronness(0),  fForceTheta(kFALSE)
 {
     fName  = name  ? name  : "MJSpectrum";
     fTitle = title ? title : "Standard program to calculate spectrum";
-}
-
+
+    // Make sure that fDisplay is maintained properly
+    // (i.e. removed via RecursiveRemove if closed)
+    gROOT->GetListOfCleanups()->Add(this);
+}
+
+// --------------------------------------------------------------------------
+//
+// Delete all the fCut* data members and fCalcHadronness/fEstimateEnergy
+//
 MJSpectrum::~MJSpectrum()
 {
@@ -128 +132 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Read a MTask named name into task from the open file. If this task is
+// required mustexist can be set. Depending on success kTRUE or kFALSE is
+// returned. If the task is a MContinue SetAllowEmpty is called.
+// The name of the task is set to name.
+//
 Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const
 {
+    // If a task is set delete task
     if (task)
     {
@@ -136 +148 @@
     }
 
+    // Try to get task from file
     task = (MTask*)gFile->Get(name);
     if (!task)
@@ -144 +157 @@
         return kFALSE;
     }
+
+    // Check if task inherits from MTask
     if (!task->InheritsFrom(MTask::Class()))
     {
@@ -151 +166 @@
     }
 
+    // Set name of task
     task->SetName(name);
 
+    // SetAllowEmpty for MContinue tasks
     if (dynamic_cast<MContinue*>(task))
         dynamic_cast<MContinue*>(task)->SetAllowEmpty();
@@ -159 +176 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Print setup used for the MC processing, namely MAlphaFitter ans all fCut*.
+//
 void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
 {
@@ -194 +215 @@
     }
 
-    // Read the first corsika RunHeader from the MC files
-    TChain chain("RunHeaders");
-    if (!set.AddFilesOn(chain))
-        return kFALSE;
-
-    MMcCorsikaRunHeader *h=0;
-    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
-    chain.GetEntry(1);
-
-    if (!h)
-    {
-        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
-        return kFALSE;
-    }
-
-    // Propagate the run header to MMcSpectrumWeight
-    if (!w.Set(*h))
-    {
-        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
-        return kFALSE;
-    }
-
-    // Print new setup of MMcSpectrumWeight
-    w.Print();
+    w.Print("new");
+
     return kTRUE;
 }
@@ -281 +280 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// return maximum value of MMcRunHeader.fImpactMax stored in the RunHeaders
+// of the files from the MC dataset
+//
+Bool_t MJSpectrum::AnalyzeMC(const MDataSet &set, Float_t &impactmax, Float_t &emin/*, Float_t emax*/) const
+{
+    if (fDisplay)
+        fDisplay->SetStatusLine1("Analyzing Monte Carlo headers...");
+
+    // ----- Create chain -----
+    *fLog << inf << "Getting files... " << flush;
+    TChain chain("RunHeaders");
+    if (!set.AddFilesOn(chain))
+        return kFALSE;
+    *fLog << "done. " << endl;
+
+    *fLog << all;
+    *fLog << "Searching for maximum impact... " << flush;
+    impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax");
+    *fLog << "found " << impactmax/100 << "m" << endl;
+
+    *fLog << "Searching for minimum lower energy limit... " << flush;
+    emin = chain.GetMinimum("MMcCorsikaRunHeader.fELowLim");
+    *fLog << "found " << emin << "GeV" << endl;
+
+    *fLog << "Searching for maximum lower energy limit... " << flush;
+    const Float_t emin2 = chain.GetMaximum("MMcCorsikaRunHeader.fELowLim");
+    *fLog << "found " << emin2 << "GeV" << endl;
+
+    if (emin2>emin)
+    {
+        *fLog << warn;
+        *fLog << "WARNING - You are using files with different lower limits for the simulated" << endl;
+        *fLog << "          energy, i.e. that the collection area and your correction" << endl;
+        *fLog << "          coefficients between " << emin << "GeV and ";
+        *fLog << emin2 << "GeV might be wrong." << endl;
+    }
+
+/*
+    *fLog << "Getting maximum energy... " << flush;
+    emax = chain.GetMaximum("MMcCorsikaRunHeader.fEUppLim");
+    *fLog << "found " << emax << "GeV" << endl;
+ */
+    return kTRUE;
+}
+
 Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
 {
@@ -286 +332 @@
     fLog->Separator("Compiling original MC distribution");
 
-    // The name of the input container is MMcEvtBasic
-    weight.SetNameMcEvt("MMcEvtBasic");
-    // Get the corresponding rule for the weights
-    const TString w(weight.GetFormulaWeights());
-    // Set back to MMcEvt
-    weight.SetNameMcEvt();
-
-    *fLog << inf << "Using weights: " << w << endl;
-    *fLog << "Please stand by, this may take a while..." << endl;
-
-    if (fDisplay)
-        fDisplay->SetStatusLine1("Getting maximum impact...");
-
-    // ----- Create chain -----
-    *fLog << "Getting files... " << flush;
-    TChain chain("RunHeaders");
-    if (!set.AddFilesOn(chain))
-        return kFALSE;
-    *fLog << "done. " << endl;
-
-    *fLog << "Getting maximum impact... " << flush;
-    const Double_t impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax");
-    *fLog << "found " << impactmax/100 << "m" << endl;
-
-    *fLog << "Getting files... " << flush;
+    Float_t impactmax=0, Emin=0;
+    if (!AnalyzeMC(set, impactmax, Emin))
+        return kFALSE;
+
+    *fLog << inf << "Getting files... " << flush;
     MDirIter iter;
     if (!set.AddFilesOn(iter))
         return kFALSE;
     *fLog << "done. " << endl;
+
+    const Int_t tot = iter.GetNumEntries();
 
     // Prepare histogram
@@ -339 +367 @@
 
     if (fDisplay)
-        fDisplay->SetStatusLine1("Reading OriginalMC");
+        fDisplay->SetStatusLine1("Reading OriginalMC distribution...");
 
     TH1 *hfill = (TH1*)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC");
     hfill->SetDirectory(0);
 
+    *fLog << all << "Compiling simulated source spectrum... please stand by, this may take a while." << endl;
+
+    Double_t evts = 0;
+    Int_t    num  = 0;
+
+    // Reading this with a eventloop is five times slower :(
     TString fname;
-    while (1)
-    {
+    while (fDisplay)
+    {
+        if (fDisplay)
+            fDisplay->SetProgressBarPosition(Float_t(num++)/tot);
+
+        // Get next filename
         fname = iter.Next();
         if (fname.IsNull())
             break;
 
+        if (fDisplay)
+            fDisplay->SetStatusLine2(fname);
+
+        // open file
         TFile file(fname);
         if (file.IsZombie())
@@ -358 +400 @@
         }
 
+        // Get tree OriginalMC
+        TTree *t = dynamic_cast<TTree*>(file.Get("OriginalMC"));
+        if (!t)
+        {
+            *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl;
+            return kFALSE;
+        }
+        if (t->GetEntries()==0)
+        {
+            *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl;
+            continue;
+        }
+
+        // Get tree RunHeaders
         TTree *rh = dynamic_cast<TTree*>(file.Get("RunHeaders"));
         if (!rh)
@@ -364 +420 @@
             return kFALSE;
         }
-
         if (rh->GetEntries()!=1)
         {
@@ -371 +426 @@
         }
 
-        TTree *t = dynamic_cast<TTree*>(file.Get("OriginalMC"));
-        if (!t)
+        // Get corsika run header
+        MMcCorsikaRunHeader *head=0;
+        rh->SetBranchAddress("MMcCorsikaRunHeader.", &head);
+        rh->GetEntry(0);
+        if (!head)
         {
-            *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl;
+            *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from " << fname << "." << endl;
             return kFALSE;
         }
 
+        // Get the maximum impact parameter of this file. Due to different
+        // production areas an additional scale-factor is applied.
+        //
+        // Because it is assumed that the efficiency outside the production
+        // area is nearly zero no additional weight has to be applied to the
+        // events after cuts.
+        const Double_t impact = rh->GetMaximum("MMcRunHeader.fImpactMax");
+        const Double_t scale  = impactmax/impact;
+ 
+        // Propagate the run header to MMcSpectrumWeight
+        if (!weight.Set(*head))
+        {
+            *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
+            return kFALSE;
+        }
+
+        // Get the corresponding rule for the weights
+        const TString weights(weight.GetFormulaWeights("MMcEvtBasic"));
+
+        // No we found everything... go on reading contents
         *fLog << inf << "Reading OriginalMC of " << fname << endl;
-
-        if (t->GetEntries()==0)
-        {
-            *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl;
-            continue;
-        }
-
-        // Why does it crash if closed within loop (no update?)
-        //if (fDisplay)
-        //    fDisplay->SetStatusLine2(fname);
-
-        // Normalize by deviding through production area
-        const Double_t impact  = rh->GetMaximum("MMcRunHeader.fImpactMax");
-        const Double_t scale   = impactmax/impact;
-        const TString  weights = Form("%.16e*(%s)", scale*scale, w.Data());
 
         // Fill histogram from tree
@@ -403 +466 @@
         }
         hfill->SetDirectory(0);
+
+        evts += hfill->GetEntries();
+
+        // ----- Complete incomplete energy ranges -----
+        // ----- and apply production area weights -----
+        weight.CompleteEnergySpectrum(*hfill, Emin);
+
+        hfill->Scale(scale*scale);
+
+        // Add weighted data from file to final histogram
         h.Add(hfill);
     }
     delete hfill;
 
-    *fLog << "Total number of events from file: " << h.GetEntries() << endl;
+    *fLog << all << "Total number of events from files in Monte Carlo dataset: " << evts << endl;
     if (fDisplay)
         fDisplay->SetStatusLine2("done.");
 
-    if (h.GetEntries()>0)
+    if (!fDisplay || h.GetEntries()>0)
         return kTRUE;
 
-    *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
-
+    *fLog << err << "ERROR - Histogram with distribution from OriginalMC empty..." << endl;
     return kFALSE;
 }
 
-Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
-{
+void MJSpectrum::GetThetaDistribution(TH1D &temp1, TH2D &hist2) const
+{
+    TH1D &temp2 = *hist2.ProjectionX();
+
     // Display some stuff
     if (fDisplay)
@@ -430 +504 @@
         c.cd(1);
         gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
         temp1.DrawCopy();
 
@@ -435 +511 @@
         c.cd(2);
         gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
         temp2.SetName("NVsTheta");
         temp2.DrawCopy("hist");
@@ -440 +518 @@
         c.cd(4);
         gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
 
         c.cd(3);
         gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
     }
 
     // Calculate the Probability
     temp1.Divide(&temp2);
-    temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());
+    temp1.Scale(1./temp1.Integral());
 
     // Some cosmetics: Name, Axis, etc.
@@ -456 +538 @@
         temp1.DrawCopy("hist");
 
-    return kTRUE;
+    delete &temp2;
 }
 
@@ -492 +574 @@
     }
 
+    // Check whether histogram is empty
     if (proj.GetMaximum()==0)
     {
@@ -497 +580 @@
         *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap" << endl;
         *fLog << "        with the Zenith Angle distribution of your observation." << endl;
-        *fLog << "        Maybe the energy binning is not defined or wrong." << endl;
+        *fLog << "        Maybe the energy binning is undefined or wrong (from ";
+        *fLog << h2.GetYaxis()->GetXmin() << "GeV to " << h2.GetYaxis()->GetXmax() << "GeV)" << endl;
         theta->SetLineColor(kRed);
         return kFALSE;;
     }
 
+    // scale histogram and set new maximum for display
     proj.Scale(theta->GetMaximum()/proj.GetMaximum());
     theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
 
+    // draw project
     proj.Draw("same");
 
@@ -688 +774 @@
     fill0.SetFilter(&sel1);
 
-    if (!fRawMc)
+    //if (!fRawMc)
         tlist1.AddToList(&sel1);
     tlist1.AddToList(&fill0);
@@ -767 +853 @@
 TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
 {
-    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
-    f.SetParameter(0, -2.5);
-    f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000)));
+    Axis_t lo, hi;
+    MH::GetRangeUser(spectrum, lo, hi);
+
+    TF1 f("f", "[1]*(x/500)^[0]", lo, hi);
+    f.SetParameter(0, -3.0);
+    f.SetParameter(1, spectrum.GetMaximum());
     f.SetLineColor(kBlue);
-    spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped
+    f.SetLineWidth(2);
+    spectrum.Fit(&f, "NIR"); // M skipped
     f.DrawCopy("same");
 
@@ -881 +971 @@
     spec = (TH1D*)spectrum.DrawCopy();
 
-    // Calculate Spectrum * E^2
-    for (int i=0; i<spectrum.GetNbinsX(); i++)
-    {
-        const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)*spectrum.GetBinLowEdge(i+2))*1e-3;
-
-        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
-        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *e*e);
-    }
+    TF1 fc("Crab", "7.0e-6*(x/300)^(-2.31-0.26*log10(x/300))", 100, 6000);
+    fc.SetLineStyle(7);
+    fc.SetLineWidth(2);
+    fc.SetLineColor(14);
+    fc.DrawCopy("same");
+
+    TF1 fa("PG1553", "1.8e-6*(x/200)^-4.21", 90, 600);
+    static_cast<const TAttLine&>(fc).Copy(fa);
+    fa.DrawCopy("same");
 
     // Display dN/dE*E^2 for conveinience
@@ -896 +987 @@
     gPad->SetGrid();
 
+    // Calculate Spectrum * E^2
+    for (int i=0; i<spectrum.GetNbinsX(); i++)
+    {
+        const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)*spectrum.GetBinLowEdge(i+2))*1e-3;
+
+        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
+        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *e*e);
+    }
+
     spectrum.SetName("FluxStd");
     spectrum.SetMarkerStyle(kFullDotMedium);
     spectrum.SetTitle("Differential flux times E^{2}");
-    spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
+    spectrum.SetYTitle("E^{2}#cdot dN/dE [N#cdot TeV/sm^{2}]");
     spectrum.SetDirectory(0);
     spectrum.DrawCopy();
+
+    TF1 fc2("Crab*E^2", "x^2*Crab*1e-6", 100, 6000);
+    static_cast<const TAttLine&>(fc).Copy(fc2);
+    fc2.DrawCopy("same");
+
+    TF1 fa2("PG1553*E^2", "x^2*PG1553*1e-6", 100, 6000);
+    static_cast<const TAttLine&>(fc).Copy(fa2);
+    fa2.DrawCopy("same");
 
     // Fit Spectrum
@@ -943 +1051 @@
     return res;
   */
-/*
-    // Plot other spectra  from Whipple
-    f.SetParameter(0, -2.45);
-    f.SetParameter(1, 3.3e-7);
-    f.SetRange(300, 8000);
-    f.SetLineColor(kBlack);
-    f.SetLineStyle(kDashed);
-    f.DrawCopy("same");
-
-    // Plot other spectra  from Cangaroo
-    f.SetParameter(0, -2.53);
-    f.SetParameter(1, 2.0e-7);
-    f.SetRange(7000, 50000);
-    f.SetLineColor(kBlack);
-    f.SetLineStyle(kDashed);
-    f.DrawCopy("same");
-
-    // Plot other spectra  from Robert
-    f.SetParameter(0, -2.59);
-    f.SetParameter(1, 2.58e-7);
-    f.SetRange(150, 1500);
-    f.SetLineColor(kBlack);
-    f.SetLineStyle(kDashed);
-    f.DrawCopy("same");
-
-    // Plot other spectra  from HEGRA
-    f.SetParameter(0, -2.61);
-    f.SetParameter(1, 2.7e-7);
-    f.SetRange(1000, 20000);
-    f.SetLineColor(kBlack);
-    f.SetLineStyle(kDashed);
-    f.DrawCopy("same");
- */
 }
 
@@ -1179 +1254 @@
     *fLog << endl;
 
+    // Setup everything which is read from the ganymed file
     MBinning bins1("BinningAlpha");
     MBinning bins2("BinningEnergyEst");
@@ -1191 +1267 @@
     MBinning binsa("BinningAsym");
     MBinning binsb("BinningConc1");
-
 
     MAlphaFitter fit;
@@ -1209 +1284 @@
     plist.AddToList(&fit);
 
-    TH1D temp1, size;
-    const Float_t ontime = ReadInput(plist, temp1, size);
+    // Read from the ganymed file
+    TH1D htheta, size;
+    const Float_t ontime = ReadInput(plist, htheta, size);
     if (ontime<0)
     {
@@ -1217 +1293 @@
     }
 
-    plist.AddToList(&bins2);
-
-    // Initialize weighting to a new spectrum
-    // as defined in the resource file
+    // Set Zenith angle binning to binning from the ganymed-file
+    bins3.SetEdges(htheta, 'x');
+
+    // Read energy binning from resource file
+    if (!CheckEnv(bins2))
+    {
+        *fLog << err << "ERROR - Reading energy binning from resources failed." << endl;
+        return kFALSE;
+    }
+    plist.AddToList(&bins2); // For later use in MC processing
+
+    // Initialize weighting to a new spectrum as defined in the resource file
     MMcSpectrumWeight weight;
     if (!InitWeighting(set, weight))
         return kFALSE;
 
-    bins3.SetEdges(temp1, 'x');
-
-    // Read the MC distribution as produced by corsika into
-    // temp2 (1D) and apply the weights previously determined
-    TH1D temp2(temp1);
-    if (!ReadOrigMCDistribution(set, temp2, weight))
-        return kFALSE;
-
-    // Calculate the weights according to the zenith angle distribution
-    // of the raw-data which have to be applied to the MC events
-    if (!GetThetaDistribution(temp1, temp2))
-        return kFALSE;
-
-    // Tell the weighting task about the ZA-weights
-    if (!fNoThetaWeights)
-        weight.SetWeightsZd(&temp1);
+    // Print Theta and energy binning for cross-checks
+    *fLog << all << endl;
+    bins2.Print();
+    bins3.Print();
+
+    // Now we read the MC distribution as produced by corsika
+    // vs zenith angle and energy.
+    // Weight for the new energy spectrum defined in MMcSpectumWeight
+    // are applied.
+    // Also correction for different lower energy bounds and
+    // different production areas (impact parameters) are applied.
+    TH2D hist;
+    hist.UseCurrentStyle();
+    MH::SetBinning(&hist, &bins3, &bins2);
+    if (!ReadOrigMCDistribution(set, hist, weight))
+        return kFALSE;
+
+    // Check if user has closed the display
+    if (!fDisplay)
+        return kTRUE;
+
+    // Display histograms and calculate za-weights into htheta
+    GetThetaDistribution(htheta, hist);
+
+    // Give the zenoith angle weights to the weighting task
+    weight.SetWeightsZd(&htheta);
+
+    // No we apply the the zenith-angle-weights to the corsika produced
+    // MC distribution. Unfortunately this must be done manually
+    // because we are multiplying column by column
+    for (int y=0; y<hist.GetNbinsY(); y++)
+        for (int x=0; x<hist.GetNbinsX(); x++)
+        {
+            hist.SetBinContent(x, y, hist.GetBinContent(x, y)*htheta.GetBinContent(x));
+            hist.SetBinError(x, y,   hist.GetBinError(x, y)  *htheta.GetBinContent(x));
+        }
+
+    // Display the resulting distribution and check it matches
+    // the observation time distribution (this could be false
+    // for example if you miss MCs of some zenith angles, which you have
+    // data for)
+    if (!DisplayResult(hist))
+        return kFALSE;
+
+    // -------------- Fill excess events versus energy ---------------
 
     // Refill excess histogram to determin the excess events
@@ -1247 +1360 @@
         return kFALSE;
 
-    // Print the setup and result of the MAlphaFitter
-    // Print used cuts
+    // Print the setup and result of the MAlphaFitter, print used cuts
     PrintSetup(fit);
-
-    TH2D hist;
-    MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
-    //if (fSimpleMode)
-    {
-        // Now we read the MC distribution as produced by corsika again
-        // with the spectral weights applied into a histogram vs.
-        // zenith angle and energy
-        hist.UseCurrentStyle();
-        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
-        if (!ReadOrigMCDistribution(set, hist, weight))
-            return kFALSE;
-
-        // No we apply the the zenith-angle-weights to the corsika produced
-        // MC distribution. Unfortunately this must be done manually
-        // because we are multiplying column by column
-        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));
-                    hist.SetBinError(x, y,   hist.GetBinError(x, y)  *temp1.GetBinContent(x));
-                }
-        }
-    }/*
-    else
-    {
-        // This rereads the original MC spectrum and aplies both
-        // weights, spectral weights and ZA-weights.
-        weight.SetNameMcEvt("MMcEvtBasic");
-        if (!IntermediateLoop(plist, mh1, temp1, set, weight))
-            return kFALSE;
-        weight.SetNameMcEvt();
-    }*/
-
-    if (!DisplayResult(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/))
-        return kFALSE;
 
     // ------------------------- Final loop --------------------------
@@ -1304 +1378 @@
 
     // Selector to get correct (final) theta-distribution
-    temp1.SetXTitle("MPointingPos.fZd");
-
-    MH3 mh3(temp1);
-
-    MFEventSelector2 sel2(mh3);
-    sel2.SetHistIsProbability();
-
-    MContinue contsel(&sel2);
-    contsel.SetInverted();
+    //temp1.SetXTitle("MPointingPos.fZd");
+    //
+    //MH3 mh3(temp1);
+    //
+    //MFEventSelector2 sel2(mh3);
+    //sel2.SetHistIsProbability();
+    //
+    //MContinue contsel(&sel2);
+    //contsel.SetInverted();
 
     // Get correct source position
@@ -1318 +1392 @@
 
     // Calculate corresponding Hillas parameters
-    /*
-     MHillasCalc hcalc1;
-     MHillasCalc hcalc2("MHillasCalcAnti");
-     hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
-     hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
-     hcalc2.SetNameHillasSrc("MHillasSrcAnti");
-     hcalc2.SetNameSrcPosCam("MSrcPosAnti");
-     */
+    //MHillasCalc hcalc1;
+    //MHillasCalc hcalc2("MHillasCalcAnti");
+    //hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
+    //hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
+    //hcalc2.SetNameHillasSrc("MHillasSrcAnti");
+    //hcalc2.SetNameSrcPosCam("MSrcPosAnti");
 
     // Fill collection area and energy estimator (unfolding)
@@ -1331 +1403 @@
     MHCollectionArea area0("TriggerArea");
     MHCollectionArea area1;
-    area0.SetHistAll(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/);
-    area1.SetHistAll(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/);
+    area0.SetHistAll(hist);
+    area1.SetHistAll(hist);
 
     MHEnergyEst      hest;
@@ -1371 +1443 @@
     MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
     MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
+    fill0a.SetNameTab("EvtDist");
     fill1a.SetNameTab("PreCut");
     fill2a.SetNameTab("PostCut");
@@ -1399 +1472 @@
     // be thrown away according to the theta distribution
     // it is enabled here
-    if (!fRawMc && fNoThetaWeights)
-        tlist2.AddToList(&contsel);
+    //if (!fRawMc && fNoThetaWeights)
+    //    tlist2.AddToList(&contsel);
     //tlist2.AddToList(&calc);
     //tlist2.AddToList(&hcalc1);

trunk/MagicSoft/Mars/mjobs/MJSpectrum.h

@@ -33 +33 @@
     MTask *fCalcHadronness;
 
-    Bool_t fRefill;
-    //Bool_t fSimpleMode;
     Bool_t fForceTheta;
-    Bool_t fRawMc;
-    Bool_t fNoThetaWeights;
 
     // Read Input
-    Bool_t  ReadTask(MTask* &task, const char *name, Bool_t mustexist=kTRUE) const;
-    Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size);
-    Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &w) const;
-    Bool_t  GetThetaDistribution(TH1D &temp1, TH1D &temp2) const;
-    Bool_t  Refill(MParList &plist, TH1D &h) /*const*/;
-    Bool_t  InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
+    Bool_t   ReadTask(MTask* &task, const char *name, Bool_t mustexist=kTRUE) const;
+    Float_t  ReadInput(MParList &plist, TH1D &h1, TH1D &size);
+    Bool_t   AnalyzeMC(const MDataSet &set, Float_t &impactmax, Float_t &emin/*, Float_t emax*/) const;
+    Bool_t   ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &w) const;
+    void     GetThetaDistribution(TH1D &temp1, TH2D &temp2) const;
+    Bool_t   Refill(MParList &plist, TH1D &h) /*const*/;
+    Bool_t   InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
 
     // Display Output
@@ -63 +60 @@
     Bool_t Process(const MDataSet &set);
 
-    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 ForceTheta(Bool_t b=kTRUE)       { fForceTheta=b; }
+    void ForceTheta(Bool_t b=kTRUE) { fForceTheta=b; }
 
     void SetEnergyEstimator(const MTask *task);

trunk/MagicSoft/Mars/sponde.rc

@@ -1 +1 @@
 EstimateEnergy.Rule: (0.380075+(MPointingPos.fZd*MPointingPos.fZd*0.00109028))*pow(MHillas.fSize,0.892462)
 
+BinningSize.Raw:        18   49   63600 log
+BinningEnergyEst.Raw:   18   53   35800 log
+
 #MMcSpectrumWeight.NewSlope: -2.26