Changeset 8674

Timestamp:

08/18/07 12:04:44 (17 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (1 diff)
• mbase/MEnv.cc (modified) (2 diffs)
• mbase/MEnv.h (modified) (2 diffs)
• mfileio/MReadTree.cc (modified) (2 diffs)
• mfileio/MWriteRootFile.cc (modified) (1 diff)
• mhflux/MHCollectionArea.cc (modified) (3 diffs)
• mhflux/MMcSpectrumWeight.cc (modified) (4 diffs)
• mjobs/MDataSet.cc (modified) (2 diffs)
• mjobs/MDataSet.h (modified) (2 diffs)
• mjobs/MJSpectrum.cc (modified) (14 diffs)
• mjobs/MSequence.cc (modified) (2 diffs)
• mjobs/MSequence.h (modified) (1 diff)
• sponde.cc (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2007/08/18 Thomas Bretz
+
+   * sponde.cc:
+     - added new option --force-theta
+
+   * mbase/MEnv.[h,cc], mjobs/MSequence.[h,cc], mjobs/MDataSet.[h,cc]:
+     - GetName noe returns only the filename not the whole path. The
+       old behaviour made it impossible to access the container from
+       the file.
+     - GetRcName now returns the whole path/name.
+     - Print now outputs also path and file-name
+
+   * mfileio/MReadTree.cc, mfileio/MWriteRootFile.cc:
+     - fixed typos in comments
+
+   * mhflux/MHCollectionArea.cc:
+     - reset fCorsikaVersion to 0 in PreProcess
+     - print old and new Cosika version if mismatch is found
+
+   * mhflux/MMcSpectrumWeight.cc:
+     - replaced the %.16f by %.16e. This is more accurate in cases
+       with high exponents
+     - added some sample/test code to weight the Zenith Angle
+       according to the sin-distribution produced by Corsika.
+       Currently not in use
+
+   * mjobs/MJSpectrum.cc:
+     - removed the simple/accurate mode. There is now reason why
+       the previous "accurate"-mode should be more accurate at all.
+       It is only slower
+     - Reading the OriginalMC tree now is done such that the
+       events are properly weighted by the production area. This
+       allowes to use different impact paramters from dfferent files.
+     - a check has been implemented which compared the zenith angle
+       distribution of the data and the resulting monte carlo data.
+       Execution of the program can be forced with a new option.
+     - write more information to output file.
+
+
 
  2007/08/17 Thomas Bretz

trunk/MagicSoft/Mars/NEWS

@@ -221 +221 @@
      contents where exchanged.... fixed.
 
+   - sponde: the so called "accurate"-mode has been removed. It didn't
+     give any improvement in accuracy, only decreased execution speed.
+
+   - sponde: now checks whether the theta distribution of your on-data
+     and the theta-distribution of your Monte Carlo sample (after
+     weighting) fits. If it doesn't fit properly (eg. the Monte Carlo
+     sample is incomplete) execution is stopped. Execution can be forced
+     using the new option --force-theta. Use this option with care!
+
+   - sponde: Proper collection areas can now be constructed also from
+     Monte Carlo samples generated with different maximum impact
+     parameters. Note that in previous version you neither got 
+     a warning or failure, nor was there any obvious sign that the
+     collection area was overestimated due to usage of files with 
+     different maximum impact parameters.
+
+   - sponde: the output file now contains more information about
+     the spectrum (eg. the full 2D collection area histogram).
+     Note, that this information can only be written to the file
+     if it is stored automatically via command line argument.
+     If you only store the status display from within the display
+     the information is lost.
+
 
 

trunk/MagicSoft/Mars/mbase/MEnv.cc

@@ -71 +71 @@
     if (GetEntries()<=0 && !fname.IsNull() && fname!=name)
         ReadFile(fname, kEnvLocal);;
+}
+
+//---------------------------------------------------------------------------
+//
+// Make sure that the name used for writing doesn't contain a full path
+//
+const char *MEnv::GetName() const
+{
+    const char *pos = strrchr(GetRcName(), '/');
+    return pos>0 ? pos+1 : GetRcName();
 }
 
@@ -847 +857 @@
 //---------------------------------------------------------------------------
 //
+// Add name and full path to output
+//
+void MEnv::PrintEnv(EEnvLevel level) const
+{
+    cout << "# Path: " << GetRcName() << endl;
+    cout << "# Name: " << GetName() << endl;
+    TEnv::PrintEnv(level);
+}
+
+//---------------------------------------------------------------------------
+//
 // Print resources which have never been touched (for debugging)
 //

trunk/MagicSoft/Mars/mbase/MEnv.h

@@ -38 +38 @@
     const char *GetValue(const char *name, const char *dflt);
 
-    const char *GetName() const { return GetRcName(); }
+    const char *GetName() const { return strrchr(GetRcName(), '/')?strrchr(GetRcName(), '/')+1:GetRcName(); }
 
     Int_t       GetColor(const char *name, Int_t dftl);
@@ -67 +67 @@
     void        AddEnv(const TEnv &env, Bool_t overwrite=kTRUE);
 
+    void        PrintEnv(EEnvLevel level = kEnvAll) const;
     void        Print(Option_t *option) const { TEnv::Print(option); }
-    void        Print() const { TEnv::PrintEnv(kEnvLocal); } //*MENU*
+    void        Print() const { PrintEnv(kEnvLocal); } //*MENU*
 
     void PrintUntouched() const;

trunk/MagicSoft/Mars/mfileio/MReadTree.cc

@@ -267 +267 @@
 //  list gets a new address. Therefor we reset all the autodel flags.
 //
-//  MAKE SURE THAT ALL YOPUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY
+//  MAKE SURE THAT ALL YOUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY
 //
 Bool_t MReadTree::Notify()
@@ -295 +295 @@
     // list gets a new address. Therefor we reset all the autodel flags.
     //
-    // MAKE SURE THAT ALL YOPUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY
+    // MAKE SURE THAT ALL YOUR CUSTOM STREAMERS TAKE CARE OF ALL MEMORY
     //
     TIter NextB(fTree->GetListOfBranches());

trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc

@@ -35 +35 @@
 // There is a special mode of operation which opens a new file for each new
 // file read by the reading task (opening the new file is initiated by
-// ReInit()) For more details se the corresponding constructor.
+// ReInit()) For more details see the corresponding constructor.
 //
 // Memory based trees

trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc

@@ -183 +183 @@
     MH::SetBinning(&fHistAll, &binst, &binse);
 
-    fMcAreaRadius = -1;
+    fMcAreaRadius   = -1;
+    fCorsikaVersion =  0;
 
     return kTRUE;
@@ -205 +206 @@
 
     if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
-        *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
+    {
+        *fLog << warn;
+        *fLog << "Warning - Read files have different Corsika versions..." << endl;
+        *fLog << "          Last file=" << fCorsikaVersion << "   New file=" << runheader->GetCorsikaVersion() << endl;
+    }
     fCorsikaVersion = runheader->GetCorsikaVersion();
 
@@ -285 +290 @@
         //                         GetCollectionAreaAbs(), fMcAreaRadius);
         const TString txt = Form("r_{max}=%.0fm --> A_{max}=%.0fm^{2}",
-                                 GetCollectionAreaAbs(), fMcAreaRadius);
+                                 fMcAreaRadius, GetCollectionAreaAbs());
 
         TLatex text(0.31, 0.95, txt);

trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

@@ -254 +254 @@
 // are equal only fNorm is returned.
 //
-// The normalization constant is returned as %.16f
+// The normalization constant is returned as %.16e
 //
 // Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600))
@@ -261 +261 @@
 {
     if (GetFormulaSpecOld()==GetFormulaSpecNew())
-        return Form("%.16f", fNorm);
+        return Form("%.16e", fNorm);
 
     const Double_t iold = fNormEnergy<0 ? GetSpecOldIntegral() : CalcSpecOld(fNormEnergy);
@@ -268 +268 @@
     const Double_t norm = fNorm*iold/inew;
 
-    return Form("%.16f*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data());
+    return Form("%.16e*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data());
 }
 
@@ -447 +447 @@
     if (fWeightsZd)
     {
+        /*
+         const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
+
+         Double_t tmin = fWeightsZd->GetXaxis()->GetBinLowEdge(i)  *TMath::DegToRad();
+         Double_t tmax = fWeightsZd->GetXaxis()->GetBinLowEdge(i+1)*TMath::DegToRad();
+         Double_t wdth = fWeightsZd->GetXaxis()->GetBinWidth(i)    *TMath::DegToRad();
+
+         Double_t cont = fWeightsZd->GetBinContent(i);
+
+         Double_t integ = cos(tmin)-cos(tmax);
+
+         w = sin(tmax)*cont/integ*wdth;
+         */
+
         const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
         w = fWeightsZd->GetBinContent(i);

trunk/MagicSoft/Mars/mjobs/MDataSet.cc

@@ -274 +274 @@
 }
 
+//---------------------------------------------------------------------------
+//
+// Make sure that the name used for writing doesn't contain a full path
+//
+const char *MDataSet::GetName() const
+{
+    const char *pos = strrchr(GetRcName(), '/');
+    return pos>0 ? pos+1 : GetRcName();
+}
+
+
 // --------------------------------------------------------------------------
 //
@@ -299 +310 @@
         return;
     }
+    gLog << "# Path: " << GetRcName() << endl;
+    gLog << "# Name: " << GetName() << endl;
+    gLog << endl;
     gLog << "AnalysisNumber: " << fNumAnalysis << endl << endl;
 

trunk/MagicSoft/Mars/mjobs/MDataSet.h

@@ -67 +67 @@
     MDataSet(Int_t num=(UInt_t)-1) : fNumAnalysis(num) { }
     MDataSet(const char *fname, TString sequences="", TString data="");
+
+    const char *GetName() const;
+    const char *GetRcName() const { return fName; }
 
     void Copy(TObject &obj) const
@@ -153 +156 @@
 
     // TObject
-    void Print(Option_t *o="") const; //*MENU*
+    void Print(Option_t *o) const;
+    void Print() const { Print(); } //*MENU*
 
     ClassDef(MDataSet, 2)

trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

@@ -47 +47 @@
 #include "MLogManip.h"
 
+#include "MDirIter.h"
+
 #include "MStatusArray.h"
 #include "MStatusDisplay.h"
@@ -88 +90 @@
     : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0),
     fEstimateEnergy(0), fCalcHadronness(0), fRefill(kFALSE),
-    fSimpleMode(kTRUE), fRawMc(kFALSE), fNoThetaWeights(kFALSE)
+    /*fSimpleMode(kTRUE),*/ fForceTheta(kFALSE),
+    fRawMc(kFALSE), fNoThetaWeights(kFALSE)
 {
     fName  = name  ? name  : "MJSpectrum";
@@ -291 +294 @@
 
     *fLog << inf << "Using weights: " << w << endl;
-    *fLog << "Please stand by, this may take a while..." << flush;
+    *fLog << "Please stand by, this may take a while..." << endl;
 
     if (fDisplay)
-        fDisplay->SetStatusLine1("Compiling MC distribution...");
+        fDisplay->SetStatusLine1("Getting maximum impact...");
 
     // Create chain
-    TChain chain("OriginalMC");
+    *fLog << "Getting files... " << flush;
+    TChain chain("RunHeaders");
     if (!set.AddFilesOn(chain))
         return kFALSE;
+    *fLog << "done. " << endl;
+
+    *fLog << "Getting maximum impact... " << flush;
+    const Double_t impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax");
+    *fLog << "found " << impactmax/100 << "m" << endl;
+
+    *fLog << "Getting files... " << flush;
+    MDirIter iter;
+    if (!set.AddFilesOn(iter))
+        return kFALSE;
+    *fLog << "done. " << endl;
 
     // Prepare histogram
     h.Reset();
     h.Sumw2();
-
-    // Fill histogram from chain
-    h.SetDirectory(gROOT);
     if (h.InheritsFrom(TH2::Class()))
     {
@@ -312 +324 @@
         h.SetXTitle("\\Theta [\\circ]");
         h.SetYTitle("E [GeV]");
-        h.SetZTitle("Counts");
-        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
+        h.SetZTitle(Form("Counts normalized to r=%.0fm", impactmax/100));
     }
     else
@@ -319 +330 @@
         h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
         h.SetXTitle("\\Theta [\\circ]");
-        h.SetYTitle("Counts");
-        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
+        h.SetYTitle(Form("Counts normalized to r=%.0fm", impactmax/100));
     }
     h.SetDirectory(0);
 
-    *fLog << "done." << endl;
+    const TString cont = h.InheritsFrom(TH2::Class()) ?
+        "MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC" :
+        "MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC";
+
+    if (fDisplay)
+        fDisplay->SetStatusLine1("Reading OriginalMC");
+
+    TH1 *hfill = (TH1*)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC");
+    hfill->SetDirectory(0);
+
+    TString fname;
+    while (1)
+    {
+        fname = iter.Next();
+        if (fname.IsNull())
+            break;
+
+        TFile file(fname);
+        if (file.IsZombie())
+        {
+            *fLog << err << "ERROR - Couldn't open file " << fname << endl;
+            return kFALSE;
+        }
+
+        TTree *rh = dynamic_cast<TTree*>(file.Get("RunHeaders"));
+        if (!rh)
+        {
+            *fLog << err << "ERROR - File " << fname << " doesn't contain tree RunHeaders." << endl;
+            return kFALSE;
+        }
+
+        if (rh->GetEntries()!=1)
+        {
+            *fLog << err << "ERROR - RunHeaders of " << fname << " doesn't contain exactly one entry." << endl;
+            return kFALSE;
+        }
+
+        TTree *t = dynamic_cast<TTree*>(file.Get("OriginalMC"));
+        if (!t)
+        {
+            *fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl;
+            return kFALSE;
+        }
+
+        *fLog << inf << "Reading OriginalMC of " << fname << endl;
+
+        if (t->GetEntries()==0)
+        {
+            *fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl;
+            continue;
+        }
+
+        // Why does it crash if closed within loop (no update?)
+        //if (fDisplay)
+        //    fDisplay->SetStatusLine2(fname);
+
+        // Normalize by deviding through production area
+        const Double_t impact  = rh->GetMaximum("MMcRunHeader.fImpactMax");
+        const Double_t scale   = impactmax/impact;
+        const TString  weights = Form("%.16e*(%s)", scale*scale, w.Data());
+
+        // Fill histogram from tree
+        hfill->SetDirectory(&file);
+        t->Draw(cont, weights, "goff");
+        hfill->SetDirectory(0);
+        h.Add(hfill);
+    }
+    delete hfill;
+
+    *fLog << "Total number of events from file: " << h.GetEntries() << endl;
     if (fDisplay)
         fDisplay->SetStatusLine2("done.");
@@ -333 +412 @@
     *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
 
-    return h.GetEntries()>0;
+    return kFALSE;
 }
 
@@ -412 +491 @@
     {
         *fLog << err;
-        *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap";
-        *fLog << "        with the the Zenith Angle distribution of your observation.";
+        *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap" << endl;
+        *fLog << "        with the Zenith Angle distribution of your observation." << endl;
+        *fLog << "        Maybe the energy binning is not defined or wrong." << endl;
         theta->SetLineColor(kRed);
         return kFALSE;;
@@ -424 +504 @@
 
     // Compare both histograms
-    const Double_t prob = proj.Chi2Test(theta, "");
+    *fLog << inf << "Comparing theta distributions for data and MCs." << endl;
+    const Double_t prob = proj.Chi2Test(theta, "P");
     if (prob==1)
         return kTRUE;
@@ -431 +512 @@
     {
         *fLog << inf;
-        *fLog << "The Zenith Angle distribution of your Monte Carlos fits well";
-        *fLog << "with the the Zenith Angle distribution of your observation.";
+        *fLog << "The Zenith Angle distribution of your Monte Carlos fits well" << endl;
+        *fLog << "with the Zenith Angle distribution of your observation." << endl;
+        *fLog << "The probability for identical Theta distributions is " << prob << endl;
         return kTRUE;
     }
 
-    // <0.01: passen gar nicht
-    // ==1: passen perfekt
-    // >0.99: passer sehr gut
-
     if (prob<0.01)
     {
         *fLog << err;
-        *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit";
-        *fLog << "        with the the Zenith Angle distribution of your observation.";
+        *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit" << endl;
+        *fLog << "        with the Zenith Angle distribution of your observation." << endl;
+        *fLog << "        The probability for identical Theta distributions is " << prob << endl;
+        if (!fForceTheta)
+            *fLog << "        To force processing use --force-theta (with care!)" << endl;
         theta->SetLineColor(kRed);
-        return kFALSE;
+        return fForceTheta;
     }
 
     *fLog << warn;
-    *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well";
-    *fLog << "          with the the Zenith Angle distribution of your observation.";
+    *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well" << endl;
+    *fLog << "          with the Zenith Angle distribution of your observation." << endl;
+    *fLog << "          The probability for identical Theta distributions is " << prob << endl;
     return kTRUE;
 }
@@ -1167 +1249 @@
     TH2D hist;
     MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
-    if (fSimpleMode)
+    //if (fSimpleMode)
     {
         // Now we read the MC distribution as produced by corsika again
@@ -1189 +1271 @@
                 }
         }
-    }
+    }/*
     else
     {
@@ -1198 +1280 @@
             return kFALSE;
         weight.SetNameMcEvt();
-    }
-
-    if (!DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist()))
+    }*/
+
+    if (!DisplayResult(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/))
         return kFALSE;
 
@@ -1245 +1327 @@
     MHCollectionArea area0("TriggerArea");
     MHCollectionArea area1;
-    area0.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
-    area1.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
+    area0.SetHistAll(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/);
+    area1.SetHistAll(/*fSimpleMode ?*/ hist /*: (TH2D&)mh1.GetHist()*/);
 
     MHEnergyEst      hest;
@@ -1400 +1482 @@
     TObjArray cont;
     cont.Add((TObject*)GetEnv());
+    cont.Add(&area0);
+    cont.Add(&area1);
+    cont.Add(&hest);
+    cont.Add(const_cast<TEnv*>(GetEnv()));
 
     if (fDisplay)

trunk/MagicSoft/Mars/mjobs/MSequence.cc

@@ -478 +478 @@
 }
 
+//---------------------------------------------------------------------------
+//
+// Make sure that the name used for writing doesn't contain a full path
+//
+const char *MSequence::GetName() const
+{
+    const char *pos = strrchr(GetRcName(), '/');
+    return pos>0 ? pos+1 : GetRcName();
+}
+
+
 // --------------------------------------------------------------------------
 //
@@ -490 +501 @@
         return;
     }
+    gLog << "# Path: " << GetRcName() << endl;
+    gLog << "# Name: " << GetName() << endl;
+    gLog << endl;
     gLog << "Sequence:       " << fSequence << endl;
     if (fMonteCarlo)

trunk/MagicSoft/Mars/mjobs/MSequence.h

@@ -78 +78 @@
 
     // TObject
-    void Print(Option_t *o="") const; //*MENU*
+    void Print(Option_t *o) const;
+    void Print() const { Print(); } //*MENU*
+
+    const char *GetName() const;
+    const char *GetRcName() const { return fName; }
 
     // Genaral interface

trunk/MagicSoft/Mars/sponde.cc

@@ -68 +68 @@
     gLog << "   --print-files             Print Files taken from Sequences ('+' found, '-' missing)" << endl;
     gLog << "   --config=sponde.rc        Resource file [default=sponde.rc]" << endl;
+    gLog << "   --force-theta             Force execution even with non-fitting theta distributions." << endl;
     gLog << endl;
     gLog << "   --version, -V             Show startup message with version number" << endl;
@@ -122 +123 @@
     const Bool_t  kRefill        =  arg.HasOnlyAndRemove("--refill");
 //    const Bool_t  kSimple        = !arg.HasOnlyAndRemove("--accurate");
+    const Bool_t  kForceTheta    =  arg.HasOnlyAndRemove("--force-theta");
     const Bool_t  kRawMc         =  arg.HasOnlyAndRemove("--raw-mc");
 
@@ -223 +225 @@
     job.SetPathIn(kInfile);
 
+    job.ForceTheta(kForceTheta);
     job.EnableRefilling(kRefill);
 //    job.EnableSimpleMode(kSimple);