Changeset 7539

Timestamp:

02/27/06 14:06:48 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mjtrain/MJTrainSeparation.cc (modified) (4 diffs)
• mjtrain/MJTrainSeparation.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -63 +63 @@
      - added GetForest
 
+   * mjtrain/MJTrainSeparation.[h,cc]:
+     - addedsome code for upcomming automatic event selection
+
 
 

trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc

@@ -18 +18 @@
 !   Author(s): Thomas Bretz 11/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2005
+!   Copyright: MAGIC Software Development, 2006
 !
 !
@@ -30 +30 @@
 #include "MJTrainSeparation.h"
 
+#include <TF1.h>
 #include <TH2.h>
+#include <TChain.h>
 #include <TGraph.h>
 #include <TVirtualPad.h>
@@ -136 +138 @@
 }
 
+/*
+Bool_t  MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
+{
+    fLog->Separator("Initialize energy weighting");
+
+    if (!CheckEnv(w))
+    {
+        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
+        return kFALSE;
+    }
+
+    TChain chain("RunHeaders");
+    set.AddFilesOn(chain);
+
+    MMcCorsikaRunHeader *h=0;
+    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
+    chain.GetEntry(1);
+
+    if (!h)
+    {
+        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
+        return kFALSE;
+    }
+
+    if (!w.Set(*h))
+    {
+        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
+        return kFALSE;
+    }
+
+    w.Print();
+    return kTRUE;
+}
+
+Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
+{
+    // Some debug output
+    fLog->Separator("Compiling original MC distribution");
+
+    weight.SetNameMcEvt("MMcEvtBasic");
+    const TString w(weight.GetFormulaWeights());
+    weight.SetNameMcEvt();
+
+    *fLog << inf << "Using weights: " << w << endl;
+    *fLog << "Please stand by, this may take a while..." << flush;
+
+    if (fDisplay)
+        fDisplay->SetStatusLine1("Compiling MC distribution...");
+
+    // Create chain
+    TChain chain("OriginalMC");
+    set.AddFilesOn(chain);
+
+    // Prepare histogram
+    h.Reset();
+
+    // Fill histogram from chain
+    h.SetDirectory(gROOT);
+    if (h.InheritsFrom(TH2::Class()))
+    {
+        h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
+        h.SetXTitle("\\Theta [\\circ]");
+        h.SetYTitle("E [GeV]");
+        h.SetZTitle("Counts");
+        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
+    }
+    else
+    {
+        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.SetDirectory(0);
+
+    *fLog << "done." << endl;
+    if (fDisplay)
+        fDisplay->SetStatusLine2("done.");
+
+    if (h.GetEntries()>0)
+        return kTRUE;
+
+    *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
+
+    return h.GetEntries()>0;
+}
+*/
+
+Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
+{
+    TChain chain("OriginalMC");
+    set.AddFilesOn(chain);
+
+    min = chain.GetMinimum("MMcEvtBasic.fEnergy");
+    max = chain.GetMaximum("MMcEvtBasic.fEnergy");
+
+    num = chain.GetEntries();
+
+    if (num<100)
+        *fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;
+
+    return num>=100;
+}
+
+Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const
+{
+    TChain chain1("Events");
+    set.AddFilesOff(chain1);
+
+    const Double_t num = chain1.GetEntries();
+    if (num<100)
+    {
+        *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
+        return -1;
+    }
+
+    TChain chain("EffectiveOnTime");
+    set.AddFilesOff(chain);
+
+    TH1F h;
+    h.SetDirectory(gROOT);
+    h.SetNameTitle("OnTime", "Effective on-time");
+    chain.Draw("MEffectiveOnTime.fVal>>OnTime", "", "goff");
+    h.SetDirectory(0);
+
+    if (h.Integral()<1)
+    {
+        *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
+        return -1;
+    }
+
+    *fLog << inf << "Found " << num << " events in " << h.Integral();
+    *fLog << "s (" << num/h.Integral() << "Hz)" << endl;
+
+    return num/h.Integral();
+}
+
+/*
+ Scale:
+
+
+ TF1 fold("old", "x^(-2.6)", emin, emax);
+ TF1 fnew("new", "x^(-4.0)", emin, emax);
+
+ TF1 q("q", "new/old", emin, emax);
+
+ Double_t scale = 1./q.GetMaximum(emin, emax);
+
+ // Anzahl produzierter Events vor MFEnergySlope:
+ Double_t nold = fold.Integral(emin, emax);
+
+ // Anzahl produzierter Events nach MFEnergySlope:
+ Double_t nnew = fnew.Integral(emin, emax)*scale;
+
+ class MFSpectrum : MMcSpectrumWeight
+ {
+ Double_t fScale;
+ Bool_t   fResult;
+
+ MFSpectrum::MFSpectrum(const char *name, const char *title)
+ {
+    fName  = name  ? name  : "MMcSpectrumWeight";
+    fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 
+
+    Init(fName, fTitle);
+
+ }
+
+ Int_t PreProcess(MParList *pList)
+ {
+     Int_t rc = MFSpectrumWeight::PreProcess(pList);
+     if (rc!=kTRUE)
+        return rc;
+
+     fScale = fEval->GetMaximum(fEnergyMin, fEnergyMax);
+
+     return kTRUE;
+ }
+
+ Int_t Process()
+ {
+    const Double_t e = fMcEvt->GetEnergy();
+
+    Double_t prob = fFunc->Eval(e)/fScale;
+
+    const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope);
+    const Float_t Nrnd = ;
+
+    fResult = Nexp >= gRandom->Uniform();
+ }
+
+ }
+
+
+ */
+
+Bool_t MJTrainSeparation::AutoTrain()
+{
+    Double_t num, min, max;
+    if (!GetEventsProduced(fDataSetTrain, num, min, max))
+        return kFALSE;
+
+    *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;
+
+    // Target spectrum
+    TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max);
+    flx.SetParameter(0, 1e-5);
+
+    // Number n0 of events this spectrum would produce per s and m^2
+    const Double_t n0 = flx.Integral(min, max);    //[#]
+
+    // Area produced in MC
+    const Double_t A = TMath::Pi()*300*300;        //[m²]
+
+    // Rate R of events this spectrum would produce per s
+    const Double_t R = n0*A;                       //[Hz]
+
+    // Number N of events produced (in trainings sample)
+    const Double_t N = num;                        //[#]
+
+    // This correponds to an observation time T [s]
+    const Double_t T = N/R;                        //[s]
+
+    // With an average data rate after star of
+    const Double_t r = GetDataRate(fDataSetTrain); //[Hz]
+
+    // this yields a number of n events to be read for training
+    const Double_t n = r*T;                        //[#]
+
+    if (r<0)
+        return kFALSE;
+
+    *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl;
+    *fLog << "For a data rate of " << r << "Hz this corresponds to " << n << " data events." << endl;
+
+    fNumTrainOn  = (UInt_t)-1;
+    fNumTrainOff = TMath::Nint(n);
+
+    /*
+     An event rate dependent selection?
+     ----------------------------------
+     Total average data rate:      R
+     Goal number of events:        N
+     Number of data events:        N0
+     Rate assigned to single evt:  r
+
+     Selection probability: N/N0 * r/R
+
+     f := N/N0 * r
+
+     MF f("f * MEventRate.fRate < rand");
+     */
+
+
+    return kTRUE;
+}
+
 Bool_t MJTrainSeparation::Train(const char *out)
 {
@@ -148 +407 @@
         return kFALSE;
     }
+
+    // ----------------------- Auto Train? ----------------------
+
+    if (fAutoTrain)
+        if (!AutoTrain())
+            return kFALSE;
 
     // --------------------- Setup files --------------------

trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h

@@ -24 +24 @@
     UInt_t fNumTestOff;
 
+    Bool_t fAutoTrain;
+
     void DisplayResult(MH3 &h31, MH3 &h32);
+
+    Bool_t   GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const;
+    Double_t GetDataRate(MDataSet &set) const;
+    Bool_t   AutoTrain();
 
 public:
@@ -44 +50 @@
     }
 
+    void EnableAutoTrain(Bool_t b=kTRUE) { fAutoTrain = kTRUE; }
+
     Bool_t Train(const char *out);