Changeset 1300 for trunk/MagicSoft

Timestamp:

04/25/02 11:21:46 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhist/MH.cc (modified) (8 diffs)
• mhist/MH.h (modified) (2 diffs)
• mhist/MH3.cc (modified) (1 diff)
• mhist/MHMcCollectionArea.cc (modified) (2 diffs)
• mhist/MHMcCollectionArea.h (modified) (1 diff)
• mmontecarlo/MMcCollectionAreaCalc.cc (modified) (4 diffs)
• mmontecarlo/MMcCollectionAreaCalc.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+ 2002/04/25: Thomas Bretz
+
+   * mmontecarlo/MMcCollectionAreaCalc.[h,cc]:
+     - counts now the number of simulated showers
+     - implemented some sanity checks (corsika version, etc)
+
+   * mhist/MMcCollectionArea.[h,cc]:
+     - added a first implementation of a calculation using only triggered
+       events
+
+   * mhist/MH.[h,cc]:
+     - changed the first argument in SetBinning (according to the number
+       of axis) to TH2 or TH3
+
+   * mhist/MH2.cc:
+     - changed the first argument in SetBinning (according to the number
+       of axis) to TH2 or TH3
+
+
+
  2002/04/24: Thomas Bretz
 

trunk/MagicSoft/Mars/mhist/MH.cc

@@ -50 +50 @@
 
 #include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
 #include <TCanvas.h>
 
@@ -136 +138 @@
 }
 
-void MH::SetBinning(TH1 *h, const MBinning *binsx, const MBinning *binsy)
+void MH::SetBinning(TH2 *h, const MBinning *binsx, const MBinning *binsy)
 {
     TAxis &x = *h->GetXaxis();
@@ -169 +171 @@
 }
 
-void MH::SetBinning(TH1 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz)
+void MH::SetBinning(TH3 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz)
 {
     //
@@ -214 +216 @@
 }
 
-void MH::SetBinning(TH1 *h, const TArrayD *binsx, const TArrayD *binsy)
+void MH::SetBinning(TH2 *h, const TArrayD *binsx, const TArrayD *binsy)
 {
     MBinning bx;
@@ -223 +225 @@
 }
 
-void MH::SetBinning(TH1 *h, const TArrayD *binsx, const TArrayD *binsy, const TArrayD *binsz)
+void MH::SetBinning(TH3 *h, const TArrayD *binsx, const TArrayD *binsy, const TArrayD *binsz)
 {
     MBinning bx;
@@ -245 +247 @@
 }
 
-void MH::SetBinning(TH1 *h, const TAxis *binsx, const TAxis *binsy)
+void MH::SetBinning(TH2 *h, const TAxis *binsx, const TAxis *binsy)
 {
     const Int_t nx = binsx->GetNbins();
@@ -260 +262 @@
 }
 
-void MH::SetBinning(TH1 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz)
+void MH::SetBinning(TH3 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz)
 {
     const Int_t nx = binsx->GetNbins();
@@ -281 +283 @@
 void MH::SetBinning(TH1 *h, TH1 *x)
 {
-    SetBinning(h, x->GetXaxis(), x->GetYaxis(), x->GetZaxis());
-}
+    if (h->InheritsFrom(TH3::Class()) && x->InheritsFrom(TH3::Class()))
+    {
+        SetBinning((TH3*)h, x->GetXaxis(), x->GetYaxis(), x->GetZaxis());
+        return;
+    }
+    if (h->InheritsFrom(TH2::Class()) && x->InheritsFrom(TH2::Class()))
+    {
+        SetBinning((TH2*)h, x->GetXaxis(), x->GetYaxis());
+        return;
+    }
+    if (h->InheritsFrom(TH1::Class()) && x->InheritsFrom(TH1::Class()))
+    {
+        SetBinning(h, x->GetXaxis());
+        return;
+    }
+}

trunk/MagicSoft/Mars/mhist/MH.h

@@ -7 +7 @@
 
 class TH1;
+class TH2;
+class TH3;
 class TAxis;
 class TArrayD;
@@ -28 +30 @@
 
     static void SetBinning(TH1 *h, const MBinning *binsx);
-    static void SetBinning(TH1 *h, const MBinning *binsx, const MBinning *binsy);
-    static void SetBinning(TH1 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz);
+    static void SetBinning(TH2 *h, const MBinning *binsx, const MBinning *binsy);
+    static void SetBinning(TH3 *h, const MBinning *binsx, const MBinning *binsy, const MBinning *binsz);
 
     static void SetBinning(TH1 *h, const TArrayD *binsx);
-    static void SetBinning(TH1 *h, const TArrayD *binsx, const TArrayD *binsy);
-    static void SetBinning(TH1 *h, const TArrayD *binsx, const TArrayD *binsy, const TArrayD *binsz);
+    static void SetBinning(TH2 *h, const TArrayD *binsx, const TArrayD *binsy);
+    static void SetBinning(TH3 *h, const TArrayD *binsx, const TArrayD *binsy, const TArrayD *binsz);
 
     static void SetBinning(TH1 *h, const TAxis *binsx);
-    static void SetBinning(TH1 *h, const TAxis *binsx, const TAxis *binsy);
-    static void SetBinning(TH1 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz);
+    static void SetBinning(TH2 *h, const TAxis *binsx, const TAxis *binsy);
+    static void SetBinning(TH3 *h, const TAxis *binsx, const TAxis *binsy, const TAxis *binsz);
 
     static void SetBinning(TH1 *h, TH1 *x);

trunk/MagicSoft/Mars/mhist/MH3.cc

@@ -229 +229 @@
     case 2:
         fHist->SetTitle(title+" (2D)");
-        SetBinning(fHist, binsx, binsy);
+        SetBinning((TH2*)fHist, binsx, binsy);
         return kTRUE;
     case 3:
         fHist->SetTitle(title+" (3D)");
-        SetBinning(fHist, binsx, binsy, binsz);
+        SetBinning((TH3*)fHist, binsx, binsy, binsz);
         return kTRUE;
     }

trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc

@@ -196 +196 @@
 }
 
+// --------------------------------------------------------------------------
+//
+//  Calculate the Efficiency (collection area) and set the 'ReadyToSave'
+//  flag
+//
+void MHMcCollectionArea::CalcEfficiency()
+{
+    MHMcEfficiency heff;
+    heff.Calc(*fHistSel, *fHistAll);
+
+    Calc(heff);
+}
+
+void MHMcCollectionArea::CalcEfficiency(UInt_t numevts, Float_t angle)
+{
+    // Here we estimate the total number of showers in each energy bin
+    // known the energy range and spectral index of the generated showers
+    // (set now by hand since the info is not available in the header!)
+
+    TH1D &histsel = *fHistSel->ProjectionX();
+
+    TH1D histall;
+
+    TAxis &xaxis = *histsel.GetXaxis();
+
+    MH::SetBinning(&histall, &xaxis);
+    MH::SetBinning(fHistCol, &xaxis);
+
+    const Float_t emin  = 10.;
+    const Float_t emax  = 30000.;  // Energies in GeV.
+    const Float_t index = 2.6;     // Differential spectral Index.
+
+    const Float_t expo = 1.-index;
+
+    const Float_t k = (Float_t)numevts / (pow(emax,expo) - pow(emin,expo));
+
+    const Int_t nbinx = xaxis.GetNbins();
+
+    for (Int_t i=1; i<=nbinx; i++)
+    {
+        const Float_t e1 = histall.GetBinLowEdge(i);
+        const Float_t e2 = histall.GetBinLowEdge(i+1);
+
+        if (e1 < emin || e2 > emax)
+            continue;
+
+        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
+
+        histall.SetBinContent(i, events);
+        histall.SetBinError(i, sqrt(events));
+
+    }
+
+    // -----------------------------------------------------------
+
+    // Impact parameter range:
+    const Float_t r1 = 0;
+    const Float_t r2 = 400;
+
+    const Float_t dr = TMath::Pi() * (r2*r2 - r1*r1);
+
+    angle *= TMath::Pi()/180;
+
+    for (Int_t ix=1; ix<=nbinx; ix++)
+    {
+        const Float_t Na = histall.GetBinContent(ix);
+
+        if (Na <= 0)
+            continue;
+
+        const Float_t Ns = histsel.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 err = sqrt((1.-eff)*Ns)/Na;
+
+        const Float_t area = dr * cos(angle);
+
+        fHistCol->SetBinContent(ix, eff*area);
+        fHistCol->SetBinError(ix, err*area);
+    }
+
+    delete &histsel;
+
+    SetReadyToSave();
+}
+
+void MHMcCollectionArea::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall)
+{
+    MHMcEfficiency heff;
+    heff.Calc(*mcsel.GetHist(), *mcall.GetHist());
+
+    Calc(heff);
+}
+
 void MHMcCollectionArea::Calc(const MHMcEfficiency &heff)
 {
@@ -242 +344 @@
 //  flag
 //
-void MHMcCollectionArea::CalcEfficiency()
-{
-    MHMcEfficiency heff;
-    heff.Calc(*fHistSel, *fHistAll);
-
-    Calc(heff);
-}
-
-void MHMcCollectionArea::Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall)
-{
-    MHMcEfficiency heff;
-    heff.Calc(*mcsel.GetHist(), *mcall.GetHist());
-
-    Calc(heff);
-}
+/*
+void MHMcCollectionArea::Calc(MHMcEnergyImpact &mcsel, UInt_t numevents, Float_t angle)
+{
+}
+*/

trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h

@@ -38 +38 @@
 
     void CalcEfficiency();
+    void CalcEfficiency(UInt_t allevts, Float_t theta);
 
     void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall);

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc

@@ -33 +33 @@
 #include "MMcEvt.hxx"
 #include "MMcTrig.hxx"
+#include "MMcRunHeader.hxx"
 
 #include "MHMcCollectionArea.h"
@@ -73 +74 @@
         return kFALSE;
 
+
+    fTotalNumSimulatedShowers = 0;
+    fCorsikaVersion           = 0;
+    fAllEvtsTriggered         = kFALSE;
+
+    return kTRUE;
+}
+
+Bool_t MMcCollectionAreaCalc::ReInit(MParList *plist)
+{
+    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
+    if (!runheader)
+    {
+        *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
+        return kFALSE;
+    }
+
+    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
+
+    *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
+
+    if (fTheta>=0 && fTheta!=runheader->GetTelesTheta())
+        *fLog << warn << dbginf << "Warning - Read files have different TelesTheta... exit." << endl;
+
+    fTheta = runheader->GetTelesTheta();
+
+    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
+        *fLog << warn << dbginf << "Warning - Read files have different Corsika versions... exit." << endl;
+
+    fCorsikaVersion = runheader->GetCorsikaVersion();
+
+    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
+
+    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
+
     return kTRUE;
 }
@@ -81 +117 @@
     const Float_t impact = fMcEvt->GetImpact()/100.;
 
-    fCollArea->FillAll(energy, impact);
+    if (!fAllEvtsTriggered)
+        fCollArea->FillAll(energy, impact);
 
     if (fMcTrig->GetFirstLevel() <= 0)
@@ -96 +133 @@
     //   do the calculation of the effectiv area
     //
-    fCollArea->CalcEfficiency();
+    *fLog << inf << "Calculation Collection Area..." << endl;
+
+    if (!fAllEvtsTriggered)
+        fCollArea->CalcEfficiency();
+    else
+    {
+        *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl;
+        fCollArea->CalcEfficiency(fTotalNumSimulatedShowers,
+                                  fCorsikaVersion == 5200 ? fTheta : 0);
+    }
 
     return kTRUE;
 }
+

trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.h

@@ -21 +21 @@
     TString fObjName;
 
+    UInt_t fTotalNumSimulatedShowers;
+    UInt_t fCorsikaVersion;
+    Float_t fTheta;
+
+    Bool_t fAllEvtsTriggered;
+
 public:
     MMcCollectionAreaCalc(const char *input=NULL,
                           const char *name=NULL, const char *title=NULL);
+
+    Bool_t ReInit(MParList *plist);
 
     Bool_t PreProcess(MParList *pList);