Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7538) +++ trunk/MagicSoft/Mars/Changelog (revision 7539) @@ -63,4 +63,7 @@ - added GetForest + * mjtrain/MJTrainSeparation.[h,cc]: + - addedsome code for upcomming automatic event selection + Index: trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc (revision 7538) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc (revision 7539) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz 11/2005 ! -! Copyright: MAGIC Software Development, 2005 +! Copyright: MAGIC Software Development, 2006 ! ! @@ -30,5 +30,7 @@ #include "MJTrainSeparation.h" +#include #include +#include #include #include @@ -136,4 +138,261 @@ } +/* +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,4 +407,10 @@ return kFALSE; } + + // ----------------------- Auto Train? ---------------------- + + if (fAutoTrain) + if (!AutoTrain()) + return kFALSE; // --------------------- Setup files -------------------- Index: trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h =================================================================== --- trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h (revision 7538) +++ trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.h (revision 7539) @@ -24,5 +24,11 @@ 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,4 +50,6 @@ } + void EnableAutoTrain(Bool_t b=kTRUE) { fAutoTrain = kTRUE; } + Bool_t Train(const char *out);