Index: trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h
===================================================================
--- trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h	(revision 6932)
+++ trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h	(revision 6938)
@@ -11,4 +11,5 @@
 #pragma link C++ class MHFalseSource+;
 #pragma link C++ class MHEffectiveOnTime+;
+#pragma link C++ class MHCollectionArea+;
 
 #endif
Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc
===================================================================
--- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc	(revision 6938)
+++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc	(revision 6938)
@@ -0,0 +1,396 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   Author(s): Harald Kornmayer 1/2001
+!
+!   Copyright: MAGIC Software Development, 2000-2002
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+//  MHCollectionArea
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MHCollectionArea.h" 
+
+#include <TCanvas.h>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MH.h"
+#include "MBinning.h"
+
+#include "MMcEvt.hxx"
+#include "MMcRunHeader.hxx"
+#include "MMcCorsikaRunHeader.h"
+
+#include "MParList.h"
+#include "MParameters.h"
+
+ClassImp(MHCollectionArea);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+//  Creates the three necessary histograms:
+//   - selected showers (input)
+//   - all showers (input)
+//   - collection area (result)
+//
+MHCollectionArea::MHCollectionArea(const char *name, const char *title)
+  : fMcEvt(0), fEnergy(0), fMcAreaRadius(300.), fIsExtern(kFALSE)
+{ 
+    //   initialize the histogram for the distribution r vs E
+    //
+    //   we set the energy range from 2 Gev to 20000 GeV (in log 4 orders
+    //   of magnitude) and for each order we take 25 subdivision --> 100 xbins
+    //
+    //   we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins
+    //
+    fName  = name  ? name  : "MHCollectionArea";
+    fTitle = title ? title : "Collection Area vs. Energy/Theta";
+
+    fHistSel.SetName("SelEvts");
+    fHistSel.SetTitle("Number of Events after cuts");
+    fHistSel.SetXTitle("\\Theta [deg]");
+    fHistSel.SetYTitle("E [GeV]");
+    fHistSel.SetDirectory(NULL);
+    fHistSel.UseCurrentStyle();
+
+    fHistAll.SetName("AllEvts");
+    fHistAll.SetTitle("Number of events produced");
+    fHistAll.SetXTitle("\\Theta [deg]");
+    fHistAll.SetYTitle("E [GeV]");
+    fHistAll.SetDirectory(NULL);
+    fHistAll.UseCurrentStyle();
+
+    fHEnergy.SetName("CollEnergy");
+    fHEnergy.SetTitle("Collection Area vs Energy");
+    fHEnergy.SetXTitle("E [GeV]");
+    fHEnergy.SetYTitle("A [m^{2}]");
+    fHEnergy.SetDirectory(NULL);
+    fHEnergy.UseCurrentStyle();
+
+    MBinning binsa, binse, binst;
+    binse.SetEdgesLog(15, 10, 1000000);
+    binst.SetEdgesCos(50, 0, 60);
+
+    binse.Apply(fHEnergy);
+
+    MH::SetBinning(&fHistSel, &binst, &binse);
+    MH::SetBinning(&fHistAll, &binst, &binse);
+}
+
+// --------------------------------------------------------------------------
+//
+//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
+//  flag
+//
+void MHCollectionArea::CalcEfficiency()
+{
+    TH1D *hsel = fHistSel.ProjectionY();
+    TH1D *hall = fHistAll.ProjectionY();
+
+    const Int_t nbinx = hsel->GetNbinsX();
+
+    // -----------------------------------------------------------
+    //
+    // Impact parameter range:  TO BE FIXED! Impact parameter shoud be
+    // read from run header, but it is not yet in!!
+    //
+    const Float_t r1 = 0;
+    const Float_t r2 = fMcAreaRadius;
+
+    //*fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <<endl<<endl;
+    const Float_t total_area = TMath::Pi() * (r2*r2 - r1*r1);
+
+    for (Int_t ix=1; ix<=nbinx; ix++)
+    {
+        const Float_t Na = hall->GetBinContent(ix);
+
+        if (Na <= 0)
+            continue;
+
+        const Float_t Ns = hsel->GetBinContent(ix);
+
+        // Since Na is an estimate of the total number of showers generated
+        // in the energy bin, it may happen that Ns (triggered showers) is
+        // larger than Na. In that case, the bin is skipped:
+
+        if (Na < Ns)
+            continue;
+
+        const Double_t eff = Ns/Na;
+
+        const Double_t efferr = TMath::Sqrt((1.-eff)*Ns)/Na;
+	
+	const Float_t col_area  =  eff * total_area;
+	const Float_t col_area_error =  efferr * total_area;
+
+        fHEnergy.SetBinContent(ix, col_area);
+        fHEnergy.SetBinError(ix, col_area_error);
+    }
+
+    delete hsel;
+    delete hall;
+}
+
+Bool_t MHCollectionArea::SetupFill(const MParList *pl)
+{
+    fHistSel.Reset();
+    if (!fIsExtern)
+        fHistAll.Reset();
+
+    fMcEvt = (MMcEvt*)pl->FindObject("MMcEvt");
+    if (!fMcEvt)
+    {
+        *fLog << err << "MMcEvt not found... abort." << endl;
+        return kFALSE;
+    }
+
+    fEnergy = (MParameterD*)pl->FindObject("MEnergyEst", "MParameterD");
+    if (!fEnergy)
+    {
+        *fLog << err << "MEnergyEst [MParameterD] not found... abort." << endl;
+        return kFALSE;
+    }
+
+    MBinning binst, binse;
+    binst.SetEdges(fHistAll, 'x');
+    binse.SetEdges(fHistAll, 'y');
+
+    //if (!fIsExtern)
+    {
+        MBinning *bins = (MBinning*)pl->FindObject("BinningTheta", "MBinning");
+        if (bins)
+            binst.SetEdges(*bins);
+
+        bins = (MBinning*)pl->FindObject("BinningEnergy", "MBinning");
+        if (bins)
+            binse.SetEdges(*bins);
+    }
+
+    binse.Apply(fHEnergy);
+
+    MH::SetBinning(&fHistSel, &binst, &binse);
+    MH::SetBinning(&fHistAll, &binst, &binse);
+
+    return kTRUE;
+}
+
+Bool_t MHCollectionArea::ReInit(MParList *plist)
+{
+    if (fIsExtern)
+        return kTRUE;
+
+    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
+    if (!runheader)
+    {
+        *fLog << err << "MMcRunHeader not found... abort." << endl;
+        return kFALSE;
+    }
+
+    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
+    *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
+
+ /*
+    if (fTheta>=0 && fTheta!=runheader->GetTelesTheta())
+    {
+        *fLog << warn << "Warning - Read files have different TelesTheta (";
+        *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl;
+    }
+    fTheta = runheader->GetTelesTheta();
+  */
+    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
+        *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
+    fCorsikaVersion = runheader->GetCorsikaVersion();
+
+    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
+    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
+
+    MMcCorsikaRunHeader *crh = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
+    if (!crh)
+    {
+        *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
+        return kFALSE;
+    }
+
+    //
+    // Calculate approximately the original number of events in each
+    // energy bin:
+    //
+    const Float_t emin = crh->GetELowLim();
+    const Float_t emax = crh->GetEUppLim();
+    const Float_t expo = 1 + crh->GetSlopeSpec();
+    const Float_t k = runheader->GetNumSimulatedShowers() /
+        (pow(emax,expo) - pow(emin,expo));
+
+    const Int_t nbiny = fHistAll.GetNbinsY();
+
+    TAxis &axe = *fHistAll.GetYaxis();
+    for (Int_t i = 1; i <= nbiny; i++)
+    {
+        const Float_t e1 = axe.GetBinLowEdge(i);
+        const Float_t e2 = axe.GetBinLowEdge(i+1);
+
+        if (e1 < emin || e2 > emax)
+            continue;
+
+        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
+        //
+        // We fill the i-th energy bin, with the total number of events
+        // Second argument of Fill would be impact parameter of each
+        // event, but we don't really need it for the collection area,
+        // so we just put a dummy value (1.)
+        //
+
+        const Float_t energy = (e1+e2)/2.;
+        fHistAll.Fill(20, energy, events);
+        // you have  MMcRunHeader.fShowerThetaMin and MMcRunHeader.fShowerThetaMax
+    }
+
+    return kTRUE;
+}
+
+void MHCollectionArea::Paint(Option_t *option)
+{
+    TVirtualPad *pad = gPad;
+
+    TH1 *h=0;
+
+    pad->cd(1);
+    //if (gPad->FindObject("ProjAllX"))
+    //    fHistAll.ProjectionX("ProjAllX", -1, 9999, "E");
+    if (gPad->FindObject("ProjSelX"))
+        fHistSel.ProjectionX("ProjSelX", -1, 9999, "E");
+
+    pad->cd(2);
+    if (gPad->FindObject("ProjAllY"))
+        h=fHistAll.ProjectionY("ProjAllY", -1, 9999, "E");
+    if (gPad->FindObject("ProjSelY"))
+        fHistSel.ProjectionY("ProjSelY", -1, 9999, "E");
+
+    if (h->GetMaximum()>0)
+    {
+        gPad->SetLogx();
+        gPad->SetLogy();
+    }
+    pad->cd(4);
+    CalcEfficiency();
+    if (fHEnergy.GetMaximum()>0)
+    {
+        gPad->SetLogx();
+        gPad->SetLogy();
+    }
+}
+
+void MHCollectionArea::Draw(Option_t *option)
+{
+    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
+
+    // Do the projection before painting the histograms into
+    // the individual pads
+    AppendPad();
+
+    pad->SetBorderMode(0);
+    pad->Divide(2,2);
+
+    TH1D *h=0;
+
+    if (fHistSel.GetNbinsX()>1)
+    {
+        pad->cd(1);
+        gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
+        /*
+        h = fHistAll.ProjectionX("ProjAllX", -1, 9999, "E");
+        h->SetXTitle("\\Theta [\\circ]");
+        h->SetDirectory(NULL);
+        h->SetLineColor(kGreen);
+        h->SetBit(kCanDelete);
+        h->Draw();
+        */
+        h = fHistSel.ProjectionX("ProjSelX", -1, 9999, "E");
+        h->SetXTitle("\\Theta [\\circ]");
+        h->SetDirectory(NULL);
+        h->SetLineColor(kRed);
+        h->SetBit(kCanDelete);
+        h->Draw("hist"/*"same"*/);
+    }
+    else
+        delete pad->GetPad(1);
+
+    if (fHistSel.GetNbinsY()>1)
+    {
+        pad->cd(2);
+        gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
+
+        h = fHistAll.ProjectionY("ProjAllY", -1, 9999, "E");
+        h->SetDirectory(NULL);
+        h->SetLineColor(kGreen);
+        h->SetXTitle("E [GeV]");
+        h->SetBit(kCanDelete);
+        h->Draw();
+
+        h = fHistSel.ProjectionY("ProjSelY", -1, 9999, "E");
+        h->SetDirectory(NULL);
+        h->SetLineColor(kRed);
+        h->SetBit(kCanDelete);
+        h->Draw("same");
+    }
+    else
+        delete pad->GetPad(2);
+
+    if (fHEnergy.GetNbinsX()>1)
+    {
+        pad->cd(4);
+        gPad->SetBorderMode(0);
+        gPad->SetGridx();
+        gPad->SetGridy();
+        fHEnergy.Draw();
+    }
+    else
+        delete pad->GetPad(4);
+/*
+    if (fHistAll.GetNbinsY()>1)
+    {
+        pad->cd(4);
+    }
+    else
+        delete pad->GetPad(4);
+  */
+}
+
+Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
+{
+    const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/;
+    const Double_t theta  = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
+
+    //*fLog << energy << " " << theta << endl;
+
+    fHistSel.Fill(theta, energy);
+
+    return kTRUE;
+}
Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h
===================================================================
--- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h	(revision 6938)
+++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h	(revision 6938)
@@ -0,0 +1,76 @@
+#ifndef MARS_MHCollectionArea
+#define MARS_MHCollectionArea
+
+#ifndef MARS_MH
+#include "MH.h"
+#endif
+
+#ifndef ROOT_TH1
+#include <TH1.h>
+#endif
+
+#ifndef ROOT_TH2
+#include <TH2.h>
+#endif
+
+class MMcEvt;
+class MParameterD;
+
+class MHCollectionArea : public MH
+{
+private:
+    MMcEvt *fMcEvt;
+    MParameterD *fEnergy;
+
+    TH2D fHistSel;
+    TH2D fHistAll;
+
+    TH1D fHEnergy;
+
+    Float_t fMcAreaRadius; // [m] Radius of circle within which the cores of Corsika events have been uniformly distributed.
+
+    UInt_t fTotalNumSimulatedShowers;
+    UInt_t fCorsikaVersion;
+    Bool_t fAllEvtsTriggered;
+    Bool_t fIsExtern;
+
+    void Calc(TH2D &hsel, TH2D &hall);
+
+public:
+    MHCollectionArea(const char *name=NULL, const char *title=NULL);
+
+    Bool_t SetupFill(const MParList *pList);
+    Bool_t ReInit(MParList *pList);
+    Bool_t Fill(const MParContainer *par, const Stat_t weight=1);
+
+    void Draw(Option_t *option="");
+    void Paint(Option_t *option="");
+
+    void CalcEfficiency();
+
+    void SetHistAll(const TH2D &h) { h.Copy(fHistAll); fIsExtern=kTRUE; }
+
+    TH2D &GetHistAll()             { return fHistAll; }
+    const TH2D &GetHistAll() const { return fHistAll; }
+
+/*
+    void DrawAll(Option_t *option="");
+    void DrawSel(Option_t *option="");
+
+    const TH1D *GetHist()       { return fHistCol; }
+    const TH1D *GetHist() const { return fHistCol; }
+
+
+
+    void CalcEfficiency();
+    void CalcEfficiency2();
+
+    void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall);
+    void Calc(const MHMcEfficiency &heff);
+  */
+    void SetMcAreaRadius(Float_t x) { fMcAreaRadius = x; }
+
+    ClassDef(MHCollectionArea, 1)  // Data Container to calculate Collection Area
+};
+
+#endif
Index: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h
===================================================================
--- trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h	(revision 6932)
+++ trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h	(revision 6938)
@@ -48,4 +48,6 @@
     MHEnergyEst(const char *name=NULL, const char *title=NULL);
 
+    const TH3D &GetHEnergy() const { return fHEnergy; }
+
     void InitMapping(MHMatrix *mat);
     void StopMapping();
Index: trunk/MagicSoft/Mars/mhflux/Makefile
===================================================================
--- trunk/MagicSoft/Mars/mhflux/Makefile	(revision 6932)
+++ trunk/MagicSoft/Mars/mhflux/Makefile	(revision 6938)
@@ -29,4 +29,5 @@
            MHEnergyEst.cc \
 	   MHEffectiveOnTime.cc \
+           MHCollectionArea.cc \
            MHFalseSource.cc
 
