Changeset 8047

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (1 diff)
• mhbase/MH.cc (modified) (1 diff)
• mhbase/MH.h (modified) (1 diff)
• mhflux/MHCollectionArea.cc (modified) (9 diffs)
• mhflux/MHCollectionArea.h (modified) (1 diff)
• mhflux/MHEnergyEst.cc (modified) (3 diffs)
• mhflux/MHThreshold.cc (modified) (1 diff)
• mhflux/MMcSpectrumWeight.cc (modified) (2 diffs)
• mjobs/MJSpectrum.cc (modified) (16 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2006/10/11 Thomas Bretz
+
+   * mhbase/MH.[h,cc]:
+     - added a function to calculate binomial errors including weights
+       (this was added in root 5.13/04, but necessary for older versions)
+
+   * mhflux/MHCollectionArea.[h,cc]:
+     - added Sumw2() to the constructor so that the weights array gets
+       correctly initialize
+     - replaced the calculation of the binomial errors by the
+       corresponding root-function and the new MH function
+     - made sure that in all histogram operations the errors are
+       properly propagated
+     - let ReInit determine fMcRadius from MMcRunHeader
+
+   * mhflux/MHEnergyEst.cc, mhflux/MHThreshold.cc
+     - fixed the order in the constructor such that the Sumw2() does
+       correctly initialize the weights array
+
+   * mhflux/MMcSpectrumWeight.cc:
+     - a minor code reordering
+
+   * mjobs/MJSpectrum.cc:
+     - made sure that the histogram with the corsika spectrum has
+       the errors initialized and thus takes the weights correctly
+       into account
+     - corresponding to this changed some draw option to get the
+       same plots (hist) as before
+     - added a lot of comments to the code
+     - when the zenith angle weights are applied to the MC distribution
+       make sure that also the errors are correctly treated.
+
+
+
  2006/10/10 Thomas Bretz
 

trunk/MagicSoft/Mars/NEWS

@@ -3 +3 @@
  *** Version  <cvs>
 
+   - sponde: In the calculation of the collection area(s) and the 
+     distribution for MOnte Carlo and estimated energy the error
+     calculation was wrong because root didn't take the errors
+     properly into account... fixed.
 
 

trunk/MagicSoft/Mars/mhbase/MH.cc

@@ -1316 +1316 @@
 // --------------------------------------------------------------------------
 //
+// This is a workaround for rrot-version <5.13/04 to get correct
+// binomial errors even if weights have been used, do:
+//    h->Divide(h1, h2, 5, 2, "b");
+//    MH::SetBinomialErrors(*h, *h1, *h2, 5, 2);
+//
+// see http://root.cern.ch/phpBB2/viewtopic.php?p=14818
+//
+void MH::SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1, Double_t c2)
+{
+    for (Int_t binx=0; binx<=hres.GetNbinsX()+1; binx++)
+    {
+        const Double_t b1 = h1.GetBinContent(binx);
+        const Double_t b2 = h2.GetBinContent(binx);
+        const Double_t w  = c2*b2 ? (c1*b1)/(c2*b2) : 0;
+        const Double_t e1 = h2.GetBinError(binx);
+        const Double_t e2 = h1.GetBinError(binx);
+
+        const Double_t rc = ((1-2*w)*e1*e1+w*w*e2*e2)/(b2*b2);
+
+        hres.SetBinError(binx, TMath::Sqrt(TMath::Abs(rc)));
+    }
+}
+
+// --------------------------------------------------------------------------
+//
 // See MTask::PrintSkipped
 //

trunk/MagicSoft/Mars/mhbase/MH.h

@@ -76 +76 @@
     static void SetBinning(TH1 *h, const TH1 *x);
 
+    static void SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1=1, Double_t c2=1);
+
     static void RemoveFirstBin(TH1 &h);
 

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

trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

@@ -176 +176 @@
     fLog->Separator("Initialize energy weighting");
 
+    // Read Resources
     if (!CheckEnv(w))
     {
@@ -182 +183 @@
     }
 
+    // Read the first corsika RunHeader from the MC files
     TChain chain("RunHeaders");
     set.AddFilesOn(chain);
@@ -195 +197 @@
     }
 
+    // Propagate the run header to MMcSpectrumWeight
     if (!w.Set(*h))
     {
@@ -201 +204 @@
     }
 
+    // Print new setup of MMcSpectrumWeight
     w.Print();
     return kTRUE;
@@ -266 +270 @@
     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();
 
@@ -282 +289 @@
     // Prepare histogram
     h.Reset();
+    h.Sumw2();
 
     // Fill histogram from chain
@@ -331 +339 @@
         gPad->SetBorderMode(0);
         temp2.SetName("NVsTheta");
-        temp2.DrawCopy();
+        temp2.DrawCopy("hist");
 
         c.cd(4);
@@ -349 +357 @@
     temp1.SetYTitle("Probability");
     if (fDisplay)
-        temp1.DrawCopy();
+        temp1.DrawCopy("hist");
 
     return kTRUE;
@@ -580 +588 @@
     case 0:
         {
-            const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",   f1, f1);
+            const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",  f1, f1);
             const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
 
@@ -638 +646 @@
 TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
 {
+    // Create copies of the histograms
     TH1D collarea(area.GetHEnergy());
     TH1D spectrum(excess);
     TH1D weights;
+
+    // Get spill-over corrections from energy estimation
     hest.GetWeights(weights);
 
+    // Print effective on-time
     cout << "Effective On time: " << ontime << "s" << endl;
 
+    // Setup spectrum plot
+    spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
+    spectrum.SetYTitle("N/sm^{2}");
     spectrum.SetDirectory(NULL);
     spectrum.SetBit(kCanDelete);
+
+    // Divide by collection are and on-time
     spectrum.Scale(1./ontime);
     spectrum.Divide(&collarea);
-    spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
-    spectrum.SetYTitle("N/sm^{2}");
-
+
+    // Draw spectrum before applying spill-over corrections
     TCanvas &c1 = fDisplay->AddTab("Spectrum");
     c1.Divide(2,2);
@@ -679 +695 @@
     weights.DrawCopy();
 
+    // Apply spill-over correction (done't take the errors into account)
+    // They are supposed to be identical with the errors of the
+    // collection area and cancel out.
     //spectrum.Multiply(&weights);
     spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
@@ -840 +859 @@
 
     //h2->Scale(1./ontime);   //h2->Integral());
-    h3->Scale(scale);    //h3->Integral());
+    h3->Scale(scale);         //h3->Integral());
 
     h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
@@ -1001 +1020 @@
     plist.AddToList(&bins2);
 
+    // Initialize weighting to a new spectru
+    // m as defined in the resource file
     MMcSpectrumWeight weight;
     if (!InitWeighting(set, weight))
         return kFALSE;
 
+    // Print the setup of the MAlphaFitter
     PrintSetup(fit);
     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);
 
+    // Refill excess histogram to determin the excess events
     TH1D excess;
     if (!Refill(plist, excess))
@@ -1026 +1054 @@
     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()*/);
@@ -1031 +1062 @@
             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.SetEntries(hist.Integral());
+                    hist.SetBinError(x, y,   hist.GetBinError(x, y)  *temp1.GetBinContent(x));
+                }
         }
     }
     else
     {
+        // This rereads the original MC spectrumand aaplies both
+        // weights spectral weights and ZA-weights.
         weight.SetNameMcEvt("MMcEvtBasic");
         if (!IntermediateLoop(plist, mh1, temp1, set, weight))
@@ -1146 +1184 @@
 
     tlist2.AddToList(&read);
+    // If no weighting should be applied but the events should
+    // be thrown away according to the theta distribution
+    // it is enabled here
     if (!fRawMc && fNoThetaWeights)
         tlist2.AddToList(&contsel);