Changeset 1300 for trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc (modified) (2 diffs)

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)
+{
+}
+*/