Changeset 8047 for trunk/MagicSoft/Mars/mhflux

Timestamp:

10/11/06 08:55:36 (18 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• MHCollectionArea.cc (modified) (9 diffs)
• MHCollectionArea.h (modified) (1 diff)
• MHEnergyEst.cc (modified) (3 diffs)
• MHThreshold.cc (modified) (1 diff)
• MMcSpectrumWeight.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc

@@ -64 +64 @@
 //
 MHCollectionArea::MHCollectionArea(const char *name, const char *title)
-  : fMcEvt(0), /*fEnergy(0),*/ fMcAreaRadius(300.), fIsExtern(kFALSE)
+  : fMcEvt(0), fMcAreaRadius(-1), fIsExtern(kFALSE)
 { 
     //   initialize the histogram for the distribution r vs E
@@ -82 +82 @@
     fHistSel.SetDirectory(NULL);
     fHistSel.UseCurrentStyle();
+    fHistSel.SetLineColor(kBlue);
 
     fHistAll.SetName("AllEvts");
@@ -105 +106 @@
     MH::SetBinning(&fHistSel, &binst, &binse);
     MH::SetBinning(&fHistAll, &binst, &binse);
+
+    // For some unknown reasons this must be called after
+    // the binning has been initialized at least once
+    fHistSel.Sumw2();
+    fHistAll.Sumw2();
+    fHEnergy.Sumw2();
 }
 
@@ -114 +121 @@
 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!!
+    TH1D *hsel = fHistSel.ProjectionY("Spy", -1, 9999, "E");;
+    TH1D *hall = fHistAll.ProjectionY("Apy", -1, 9999, "E");
+
+    //
+    // Impact parameter range.
     //
     const Float_t totalarea = GetCollectionAreaAbs();//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 colarea      =  eff    * totalarea;
-        const Float_t colareaerror =  efferr * totalarea;
-
-        fHEnergy.SetBinContent(ix, colarea);
-        fHEnergy.SetBinError(ix,   colareaerror);
-    }
+    // "b" option: calculate binomial errors
+    fHEnergy.Divide(hsel, hall, totalarea, 1, "b");
+#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
+    MH::SetBinomialErrors(fHEnergy, *hsel, *hall, totalarea, 1);
+#endif
 
     delete hsel;
     delete hall;
 }
+
 
 Bool_t MHCollectionArea::SetupFill(const MParList *pl)
@@ -196 +180 @@
     MH::SetBinning(&fHistAll, &binst, &binse);
 
+    fMcAreaRadius = -1;
+
     return kTRUE;
 }
@@ -201 +187 @@
 Bool_t MHCollectionArea::ReInit(MParList *plist)
 {
+    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
+    if (!runheader)
+    {
+        *fLog << err << "MMcRunHeader not found... abort." << endl;
+        return kFALSE;
+    }
+
+    // FIXME: Does this need some weighting with the number of produced events?
+    if (runheader->GetImpactMax()>fMcAreaRadius*100)
+    {
+        fMcAreaRadius = 0.01*runheader->GetImpactMax(); // cm->m
+        *fLog << inf << "Maximum simulated impact: " << fMcAreaRadius << "m" << endl;
+    }
+
     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;
@@ -352 +337 @@
     if (h1 && h2 && h)
     {
-        h->Divide(h2, h1);
-        h->SetMinimum(0);
+            h->Divide(h2, h1, 1, 1, "b");
+#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
+            MH::SetBinomialErrors(*h, *h2, *h1);
+#endif
+            h->SetMinimum(0);
     }
 
@@ -433 +421 @@
         gPad->SetLogx();
         h = h2->DrawCopy();
-        h->Divide(h1);
+        h->Divide(h2, h1, 1, 1, "b");
+#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
+        MH::SetBinomialErrors(*h, *h2, *h1);
+#endif
         h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency");
         h->SetDirectory(NULL);
@@ -456 +447 @@
 Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
 {
-    const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/;
+    const Double_t energy = fMcEvt->GetEnergy();
     const Double_t theta  = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
 
-    //*fLog << energy << " " << theta << endl;
-
     fHistSel.Fill(theta, energy, weight);
 

trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h

@@ -60 +60 @@
     Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; }
 
-/*
-    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; }
 

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -94 +94 @@
     fHImpact.SetYTitle("E_{est}/E_{mc}-1");
 
-    fHEnergy.Sumw2();
-    fHImpact.Sumw2();
-    fHResolution.Sumw2();
-
     MBinning binsi, binse, binst, binsr;
     binse.SetEdgesLog(25, 10, 1000000);
@@ -107 +103 @@
     SetBinning(&fHImpact,     &binsi, &binsr);
     SetBinning(&fHResolution, &binse, &binse, &binsr);
+
+    // For some unknown reasons this must be called after
+    // the binning has been initialized at least once
+    fHEnergy.Sumw2();
+    fHResolution.Sumw2();
+    fHImpact.Sumw2();
 }
 
@@ -243 +245 @@
 {
     // Project into EnergyEst_ey
+    // the "e" ensures that errors are calculated
     TH1D *h1 = (TH1D*)fHEnergy.Project3D("rtlprmft_ex"); // E_Est
     TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey"); // E_mc

trunk/MagicSoft/Mars/mhflux/MHThreshold.cc

@@ -79 +79 @@
     fHEnergy.SetDirectory(NULL);
     fHEnergy.UseCurrentStyle();
-    fHEnergy.Sumw2();
 
     MBinning binse(50, 20, 50000, "", "log");
     binse.Apply(fHEnergy);
+
+    fHEnergy.Sumw2();
 }
 

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

@@ -443 +443 @@
 Int_t MMcSpectrumWeight::Process()
 {
-    const Double_t e = fMcEvt->GetEnergy();
-
     Double_t w = 1;
 
@@ -459 +457 @@
     }
 
+    const Double_t e = fMcEvt->GetEnergy();
     fWeight->SetVal(fFunc->Eval(e)*w);