Changeset 7539 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
02/27/06 14:06:48 (19 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r7538 r7539  
    6363     - added GetForest
    6464
     65   * mjtrain/MJTrainSeparation.[h,cc]:
     66     - addedsome code for upcomming automatic event selection
     67
    6568
    6669
  • trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc

    r7420 r7539  
    1818!   Author(s): Thomas Bretz 11/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2005
     20!   Copyright: MAGIC Software Development, 2006
    2121!
    2222!
     
    3030#include "MJTrainSeparation.h"
    3131
     32#include <TF1.h>
    3233#include <TH2.h>
     34#include <TChain.h>
    3335#include <TGraph.h>
    3436#include <TVirtualPad.h>
     
    136138}
    137139
     140/*
     141Bool_t  MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
     142{
     143    fLog->Separator("Initialize energy weighting");
     144
     145    if (!CheckEnv(w))
     146    {
     147        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
     148        return kFALSE;
     149    }
     150
     151    TChain chain("RunHeaders");
     152    set.AddFilesOn(chain);
     153
     154    MMcCorsikaRunHeader *h=0;
     155    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
     156    chain.GetEntry(1);
     157
     158    if (!h)
     159    {
     160        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
     161        return kFALSE;
     162    }
     163
     164    if (!w.Set(*h))
     165    {
     166        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
     167        return kFALSE;
     168    }
     169
     170    w.Print();
     171    return kTRUE;
     172}
     173
     174Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
     175{
     176    // Some debug output
     177    fLog->Separator("Compiling original MC distribution");
     178
     179    weight.SetNameMcEvt("MMcEvtBasic");
     180    const TString w(weight.GetFormulaWeights());
     181    weight.SetNameMcEvt();
     182
     183    *fLog << inf << "Using weights: " << w << endl;
     184    *fLog << "Please stand by, this may take a while..." << flush;
     185
     186    if (fDisplay)
     187        fDisplay->SetStatusLine1("Compiling MC distribution...");
     188
     189    // Create chain
     190    TChain chain("OriginalMC");
     191    set.AddFilesOn(chain);
     192
     193    // Prepare histogram
     194    h.Reset();
     195
     196    // Fill histogram from chain
     197    h.SetDirectory(gROOT);
     198    if (h.InheritsFrom(TH2::Class()))
     199    {
     200        h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
     201        h.SetXTitle("\\Theta [\\circ]");
     202        h.SetYTitle("E [GeV]");
     203        h.SetZTitle("Counts");
     204        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
     205    }
     206    else
     207    {
     208        h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
     209        h.SetXTitle("\\Theta [\\circ]");
     210        h.SetYTitle("Counts");
     211        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
     212    }
     213    h.SetDirectory(0);
     214
     215    *fLog << "done." << endl;
     216    if (fDisplay)
     217        fDisplay->SetStatusLine2("done.");
     218
     219    if (h.GetEntries()>0)
     220        return kTRUE;
     221
     222    *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
     223
     224    return h.GetEntries()>0;
     225}
     226*/
     227
     228Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
     229{
     230    TChain chain("OriginalMC");
     231    set.AddFilesOn(chain);
     232
     233    min = chain.GetMinimum("MMcEvtBasic.fEnergy");
     234    max = chain.GetMaximum("MMcEvtBasic.fEnergy");
     235
     236    num = chain.GetEntries();
     237
     238    if (num<100)
     239        *fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;
     240
     241    return num>=100;
     242}
     243
     244Double_t MJTrainSeparation::GetDataRate(MDataSet &set) const
     245{
     246    TChain chain1("Events");
     247    set.AddFilesOff(chain1);
     248
     249    const Double_t num = chain1.GetEntries();
     250    if (num<100)
     251    {
     252        *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
     253        return -1;
     254    }
     255
     256    TChain chain("EffectiveOnTime");
     257    set.AddFilesOff(chain);
     258
     259    TH1F h;
     260    h.SetDirectory(gROOT);
     261    h.SetNameTitle("OnTime", "Effective on-time");
     262    chain.Draw("MEffectiveOnTime.fVal>>OnTime", "", "goff");
     263    h.SetDirectory(0);
     264
     265    if (h.Integral()<1)
     266    {
     267        *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
     268        return -1;
     269    }
     270
     271    *fLog << inf << "Found " << num << " events in " << h.Integral();
     272    *fLog << "s (" << num/h.Integral() << "Hz)" << endl;
     273
     274    return num/h.Integral();
     275}
     276
     277/*
     278 Scale:
     279
     280
     281 TF1 fold("old", "x^(-2.6)", emin, emax);
     282 TF1 fnew("new", "x^(-4.0)", emin, emax);
     283
     284 TF1 q("q", "new/old", emin, emax);
     285
     286 Double_t scale = 1./q.GetMaximum(emin, emax);
     287
     288 // Anzahl produzierter Events vor MFEnergySlope:
     289 Double_t nold = fold.Integral(emin, emax);
     290
     291 // Anzahl produzierter Events nach MFEnergySlope:
     292 Double_t nnew = fnew.Integral(emin, emax)*scale;
     293
     294 class MFSpectrum : MMcSpectrumWeight
     295 {
     296 Double_t fScale;
     297 Bool_t   fResult;
     298
     299 MFSpectrum::MFSpectrum(const char *name, const char *title)
     300 {
     301    fName  = name  ? name  : "MMcSpectrumWeight";
     302    fTitle = title ? title : "Task to calculate weights to change the energy spectrum";
     303
     304    Init(fName, fTitle);
     305
     306 }
     307
     308 Int_t PreProcess(MParList *pList)
     309 {
     310     Int_t rc = MFSpectrumWeight::PreProcess(pList);
     311     if (rc!=kTRUE)
     312        return rc;
     313
     314     fScale = fEval->GetMaximum(fEnergyMin, fEnergyMax);
     315
     316     return kTRUE;
     317 }
     318
     319 Int_t Process()
     320 {
     321    const Double_t e = fMcEvt->GetEnergy();
     322
     323    Double_t prob = fFunc->Eval(e)/fScale;
     324
     325    const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope);
     326    const Float_t Nrnd = ;
     327
     328    fResult = Nexp >= gRandom->Uniform();
     329 }
     330
     331 }
     332
     333
     334 */
     335
     336Bool_t MJTrainSeparation::AutoTrain()
     337{
     338    Double_t num, min, max;
     339    if (!GetEventsProduced(fDataSetTrain, num, min, max))
     340        return kFALSE;
     341
     342    *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;
     343
     344    // Target spectrum
     345    TF1 flx("Flux", "[0]*(x/1000)^(-2.6)", min, max);
     346    flx.SetParameter(0, 1e-5);
     347
     348    // Number n0 of events this spectrum would produce per s and m^2
     349    const Double_t n0 = flx.Integral(min, max);    //[#]
     350
     351    // Area produced in MC
     352    const Double_t A = TMath::Pi()*300*300;        //[m²]
     353
     354    // Rate R of events this spectrum would produce per s
     355    const Double_t R = n0*A;                       //[Hz]
     356
     357    // Number N of events produced (in trainings sample)
     358    const Double_t N = num;                        //[#]
     359
     360    // This correponds to an observation time T [s]
     361    const Double_t T = N/R;                        //[s]
     362
     363    // With an average data rate after star of
     364    const Double_t r = GetDataRate(fDataSetTrain); //[Hz]
     365
     366    // this yields a number of n events to be read for training
     367    const Double_t n = r*T;                        //[#]
     368
     369    if (r<0)
     370        return kFALSE;
     371
     372    *fLog << "Calculated a total Monte Carlo observation time of " << T << "s" << endl;
     373    *fLog << "For a data rate of " << r << "Hz this corresponds to " << n << " data events." << endl;
     374
     375    fNumTrainOn  = (UInt_t)-1;
     376    fNumTrainOff = TMath::Nint(n);
     377
     378    /*
     379     An event rate dependent selection?
     380     ----------------------------------
     381     Total average data rate:      R
     382     Goal number of events:        N
     383     Number of data events:        N0
     384     Rate assigned to single evt:  r
     385
     386     Selection probability: N/N0 * r/R
     387
     388     f := N/N0 * r
     389
     390     MF f("f * MEventRate.fRate < rand");
     391     */
     392
     393
     394    return kTRUE;
     395}
     396
    138397Bool_t MJTrainSeparation::Train(const char *out)
    139398{
     
    148407        return kFALSE;
    149408    }
     409
     410    // ----------------------- Auto Train? ----------------------
     411
     412    if (fAutoTrain)
     413        if (!AutoTrain())
     414            return kFALSE;
    150415
    151416    // --------------------- Setup files --------------------
  • trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h

    r7420 r7539  
    2424    UInt_t fNumTestOff;
    2525
     26    Bool_t fAutoTrain;
     27
    2628    void DisplayResult(MH3 &h31, MH3 &h32);
     29
     30    Bool_t   GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const;
     31    Double_t GetDataRate(MDataSet &set) const;
     32    Bool_t   AutoTrain();
    2733
    2834public:
     
    4450    }
    4551
     52    void EnableAutoTrain(Bool_t b=kTRUE) { fAutoTrain = kTRUE; }
     53
    4654    Bool_t Train(const char *out);
    4755
Note: See TracChangeset for help on using the changeset viewer.