Changeset 8674 for trunk/MagicSoft/Mars/mjobs

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mjobs

Files:

• MDataSet.cc (modified) (2 diffs)
• MDataSet.h (modified) (2 diffs)
• MJSpectrum.cc (modified) (14 diffs)
• MSequence.cc (modified) (2 diffs)
• MSequence.h (modified) (1 diff)

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