Changeset 8047 for trunk/MagicSoft/Mars/mjobs

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (modified) (16 diffs)

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