Changeset 8680 for trunk/MagicSoft/Mars/mjobs

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mjobs

Files:

• MJSpectrum.cc (modified) (38 diffs)
• MJSpectrum.h (modified) (2 diffs)

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);