Index: /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 2117) +++ /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 2118) @@ -10,4 +10,5 @@ #pragma link C++ class MFAlpha+; #pragma link C++ class MFEventSelector+; +#pragma link C++ class MFEventSelector2+; #pragma link C++ class MFTriggerLvl1+; #pragma link C++ class MFTriggerLvl2+; Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 2117) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 2118) @@ -16,7 +16,7 @@ ! ! -! Author(s): Thomas Bretz 01/2002 -! -! Copyright: MAGIC Software Development, 2000-2002 +! Author(s): Thomas Bretz, 01/2002 +! +! Copyright: MAGIC Software Development, 2000-2003 ! ! @@ -103,12 +103,4 @@ // -------------------------------------------------------------------------- // -// Destructor. Delete filters. -// -MFEventSelector::~MFEventSelector() -{ -} - -// -------------------------------------------------------------------------- -// // PreProcess all filters. // Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 2117) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 2118) @@ -32,14 +32,9 @@ enum { kNumTotalFromFile = BIT(14) }; - /* - enum { kUseFixedRatio=BIT(14) } - */ Int_t fErrors[2]; public: - // MFEventSelector(); MFEventSelector(const char *name=NULL, const char *title=NULL); - ~MFEventSelector(); Bool_t IsExpressionTrue() const { return fResult; } @@ -49,8 +44,4 @@ void SetSelectionRatio(Float_t f) { fSelRatio = f; } - /* - void SetUseFixedRatio(Bool_t b=kTRUE) { b ? SetBit(kUseFixedRatio) : ResetBit(kUseFixedRatio); } - */ - ClassDef(MFEventSelector, 0) // A Filter to select events from files }; Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc (revision 2118) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector2.cc (revision 2118) @@ -0,0 +1,422 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 01/2002 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MFEventSelector2 +// +// This is a filter to make a selection of events from a file. +// +// see Construtor for more instructions and the example below: +// +// -------------------------------------------------------------------- +// +// void select() +// { +// MParList plist; +// MTaskList tlist; +// +// MStatusDisplay *d=new MStatusDisplay; +// +// plist.AddToList(&tlist); +// +// MReadTree read("Events", "myinputfile.root"); +// read.DisableAutoScheme(); +// // Accelerate execution... +// // read.EnableBranch("MMcEvt.fTelescopeTheta"); +// +// // create nominal distribution (theta converted into degrees) +// MH3 nomdist("r2d(MMcEvt.fTelescopeTheta)"); +// MBinning binsx; +// binsx.SetEdges(5, 0, 45); // five bins from 0deg to 45deg +// MH::SetBinning(&nomdist.GetHist(), &binsx); +// +// // use this to create a nominal distribution in 2D +// // MH3 nomdist("r2d(MMcEvt.fTelescopeTheta)", "MMcEvt.fEnergy"); +// // MBinning binsy; +// // binsy.SetEdgesLog(5, 10, 10000); +// // MH::SetBinning((TH2*)&nomdist.GetHist(), &binsx, &binsy); +// +// // Fill the nominal distribution with whatever you want: +// for (int i=0; iRndm(). You may +// control this procedure using the global object gRandom. +// +// Because of the random numbers this works best for huge samples... +// +// Don't try to use this filter for the reading task or as a selector +// in the reading task: This won't work! +// +// Remark: You can also use the filter together with MContinue +// +///////////////////////////////////////////////////////////////////////////// +#include "MFEventSelector2.h" + +#include // INT_MAX + +#include // gRandom +#include // TCanvas + +#include "MH3.h" // MH3 +#include "MRead.h" // MRead +#include "MEvtLoop.h" // MEvtLoop +#include "MTaskList.h" // MTaskList +#include "MBinning.h" // MBinning +#include "MFillH.h" // MFillH +#include "MParList.h" // MParList +#include "MStatusDisplay.h" // MStatusDisplay + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MFEventSelector2); + +const TString MFEventSelector2::gsDefName = "MFEventSelector2"; +const TString MFEventSelector2::gsDefTitle = "Filter to select events with a given distribution"; + +// -------------------------------------------------------------------------- +// +// Constructor. Takes a reference to an MH3 which gives you +// 1) The nominal distribution. The distribution is renormalized, so +// that the absolute values doesn't matter. To crop the distribution +// to a nominal value of total events use SetNumMax +// 2) The dimension of the MH3 determins the dimension in which the +// event selector will work, eg +// MH3 hist("MMcEvt.fTelescopeTheta", "MMcEvt.fEnergy"); +// Would result in a redistribution of Theta and Energy. +// 3) The given rules are the variables which are used for the +// redistribution, eg: +// MH3 hist("MMcEvt.fTelescopeTheta"); +// would result in redistributing Theta. +// MH3 hist("cos(MMcEvt.fTelescopeTheta)"); +// would result in redistributing cos(Theta). +// +MFEventSelector2::MFEventSelector2(MH3 &hist, const char *name=NULL, const char *title=NULL) +: fHistOrig(NULL), fHistNom(&hist), fHistRes(NULL), + fDataX(hist.GetRule('x')), fDataY(hist.GetRule('y')), + fDataZ(hist.GetRule('z')), fNumMax(-1) +{ + fName = name ? (TString)name : gsDefName; + fTitle = title ? (TString)title : gsDefTitle; +} + +// -------------------------------------------------------------------------- +// +// Destructor. Deletes fHistRes if allocated. +// +MFEventSelector2::~MFEventSelector2() +{ + if (fHistRes) + delete fHistRes; +} + +// -------------------------------------------------------------------------- +// +// Recreate a MH3 from fHistNom used as a template. Copy the Binning +// from fHistNom to the new histogram, and return a pointer to the TH1 +// base class of the MH3. +// +TH1 &MFEventSelector2::InitHistogram(MH3* &hist) +{ + // if fHistRes is already allocated delete it first + if (hist) + delete hist; + + // duplicate the fHistNom histogram + hist = (MH3*)fHistNom->New(); + + // copy binning from one histogram to the other one + MH::SetBinning(&hist->GetHist(), &fHistNom->GetHist()); + + return hist->GetHist(); +} + +// -------------------------------------------------------------------------- +// +// Try to read the present distribution from the file. Therefore the +// Reading task of the present loop is used in a new eventloop. +// +Bool_t MFEventSelector2::Read(MRead &read) +{ + if (read.GetEntries() > INT_MAX) // FIXME: LONG_MAX ??? + { + *fLog << err << "INT_MAX exceeded." << endl; + return kFALSE; + } + + MEvtLoop run(GetName()); + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + run.SetParList(&plist); + + MBinning binsx("BinningMH3X"); + MBinning binsy("BinningMH3Y"); + MBinning binsz("BinningMH3Z"); + binsx.SetEdges(fHistNom->GetHist(), 'x'); + binsy.SetEdges(fHistNom->GetHist(), 'y'); + binsz.SetEdges(fHistNom->GetHist(), 'z'); + plist.AddToList(&binsx); + plist.AddToList(&binsy); + plist.AddToList(&binsz); + + MFillH fill(fHistOrig); + fill.SetBit(MFillH::kDoNotDisplay); + tlist.AddToList(&read); + tlist.AddToList(&fill); + run.SetDisplay(fDisplay); + if (!run.Eventloop()) + { + *fLog << err << dbginf << "Evtloop failed." << endl; + return kFALSE; + } + + return read.Rewind(); +} + +// -------------------------------------------------------------------------- +// +// After reading the histograms and arrays used for the random event +// selction are created. If a MStatusDisplay is set the histograms are +// displayed there. +// +void MFEventSelector2::PrepareHistograms() +{ + TH1 &ho = fHistOrig->GetHist(); + TH1 &hn = fHistNom->GetHist(); + + fHistNom->SetTitle("Users Nominal Distribution"); + fHistRes->SetTitle("Resulting Distribution"); + + // normalize to number of counts in primary distribution + hn.Scale(1./hn.Integral()); + + MH3 *h3 = NULL; + TH1 &hist = InitHistogram(h3); + + hist.Divide(&hn, &ho); + hist.Scale(1./hist.GetMaximum()); + + if (fDisplay) + { + fCanvas->Clear(); + fCanvas->Divide(2,2); + + fCanvas->cd(1); + gPad->SetBorderMode(0); + hn.DrawCopy(); + + fCanvas->cd(2); + gPad->SetBorderMode(0); + ho.DrawCopy(); + } + hn.Multiply(&ho, &hist); + hn.SetTitle("Resulting Nominal Distribution"); + + if (fNumMax>0) + { + cout << "SCALE: " << fNumMax/hn.Integral() << endl; + cout << "SCALE: " << fNumMax << endl; + cout << "SCALE: " << hn.Integral() << endl; + hn.Scale(fNumMax/hn.Integral()); + } + + hn.SetEntries(hn.Integral()+0.5); + if (fCanvas) + { + fCanvas->cd(3); + gPad->SetBorderMode(0); + hn.DrawCopy(); + + fCanvas->cd(4); + gPad->SetBorderMode(0); + fHistRes->Draw(); + } + delete h3; + + const Int_t num = fHistRes->GetNbins(); + fIs.Set(num); + fNom.Set(num); + for (int i=0; iGetDimension()) + { + case 3: + if (!fDataZ.PreProcess(parlist)) + { + *fLog << err << "Preprocessing of rule for z-axis failed... abort." << endl; + return kFALSE; + } + case 2: + if (!fDataY.PreProcess(parlist)) + { + *fLog << err << "Preprocessing of rule for y-axis failed... abort." << endl; + return kFALSE; + } + case 1: + if (!fDataX.PreProcess(parlist)) + { + *fLog << err << "Preprocessing of rule for x-axis failed... abort." << endl; + return kFALSE; + } + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// PreProcess the filter. Means: +// 1) Preprocess the rules +// 2) Read The present distribution from the file. +// 3) Initialize the histogram for the resulting distribution +// 4) Prepare the random selection +// 5) Repreprocess the reading task. +// +Bool_t MFEventSelector2::PreProcess(MParList *parlist) +{ + MTaskList *tasklist = (MTaskList*)parlist->FindObject("MTaskList"); + if (!tasklist) + { + *fLog << err << "MTaskList not found... abort." << endl; + return kFALSE; + } + + MRead *read = (MRead*)tasklist->FindObject("MRead"); + if (!read) + { + *fLog << err << "MRead not found... abort." << endl; + return kFALSE; + } + + if (!PreProcessData(parlist)) + return kFALSE; + + InitHistogram(fHistOrig); + InitHistogram(fHistRes); + + fHistOrig->SetTitle("Primary Distribution"); + + // Initialize online display if requested + fCanvas = fDisplay ? &fDisplay->AddTab(GetName()) : NULL; + if (fCanvas) + fHistOrig->Draw(); + + // Read primary distribution + if (!Read(*read)) + return kFALSE; + + // Prepare histograms and arrays for selection + PrepareHistograms(); + + return read->CallPreProcess(parlist); +} + +// -------------------------------------------------------------------------- +// +// fIs[i] contains the distribution of the events still to be read from +// the file. fNom[i] contains the number of events in each bin which +// are requested. +// The events are selected by: +// gRandom->Rndm()*fIs[bin]<=fNom[bin] +// +Bool_t MFEventSelector2::Process() +{ + fResult = kFALSE; + + // get x,y and z (0 if fData not valid) + const Double_t valx=fDataX.GetValue(); + const Double_t valy=fDataY.GetValue(); + const Double_t valz=fDataZ.GetValue(); + + // get corresponding bin number + const Int_t bin = fHistNom->FindFixBin(valx, valy, valz)-1; + + // under- and overflow bins are not counted + if (bin<0) + return kTRUE; + + if (gRandom->Rndm()*fIs[bin]<=fNom[bin]) + { + // how many events do we still want to read in this bin + fNom[bin]-=1; + fResult = kTRUE; + + // fill bin (same as Fill(valx, valy, valz)) + TH1 &h = fHistRes->GetHist(); + h.AddBinContent(bin+1); + h.SetEntries(h.GetEntries()+1); + } + // how many events are still pending to be read + fIs[bin]-=1; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Update online display if set. +// +Bool_t MFEventSelector2::PostProcess() +{ + if (!fCanvas || !fDisplay) + return kTRUE; + + fCanvas->cd(4); + fHistRes->DrawClone("nonew"); + fCanvas->Modified(); + fCanvas->Update(); + + return kTRUE; +} Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h (revision 2118) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector2.h (revision 2118) @@ -0,0 +1,64 @@ +#ifndef MARS_MFEventSelector2 +#define MARS_MFEventSelector2 + +///////////////////////////////////////////////////////////////////////////// +// // +// MFEventSelector2 // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif +#ifndef ROOT_TArrayL +#include +#endif +#ifndef MARS_MDataChain +#include "MDataChain.h" +#endif + +class TH1; +class TCanvas; + +class MH3; +class MRead; + +class MFEventSelector2 : public MFilter +{ +private: + static const TString gsDefName; + static const TString gsDefTitle; + + MH3 *fHistOrig; // original distribution of the events + MH3 *fHistNom; // nominal distribution + MH3 *fHistRes; // resulting distribution + MDataChain fDataX; // data rule for the x-axis + MDataChain fDataY; // data rule for the y-axis + MDataChain fDataZ; // data rule for the z-axis + Long_t fNumMax; // Maximum number of selected events + TArrayL fIs; //! array storing the original distribution + TArrayL fNom; //! array storing the nominal distribution + TCanvas *fCanvas; //! canvas for online display + + Bool_t fResult; + + TH1 &InitHistogram(MH3* &hist); + Bool_t Read(MRead &read); + void PrepareHistograms(); + Bool_t PreProcessData(MParList *parlist); + + Bool_t PreProcess(MParList *parlist); + Bool_t Process(); + Bool_t PostProcess(); + +public: + MFEventSelector2(MH3 &hist, const char *name=NULL, const char *title=NULL); + ~MFEventSelector2(); + + void SetNumMax(Long_t max=-1) { fNumMax = max; } + Bool_t IsExpressionTrue() const { return fResult; } + + ClassDef(MFEventSelector2, 0) +}; + +#endif Index: /trunk/MagicSoft/Mars/mfilter/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mfilter/Makefile (revision 2117) +++ /trunk/MagicSoft/Mars/mfilter/Makefile (revision 2118) @@ -21,5 +21,5 @@ INCLUDES = -I. -I../mbase -I../mraw -I../mmc -I../mdata -I../manalysis \ - -I../mfileio -I../mgeom -I../mimage + -I../mfileio -I../mgeom -I../mimage -I../mhist -I../mmain # @code @@ -37,4 +37,5 @@ MFilterList.cc \ MFEventSelector.cc \ + MFEventSelector2.cc \ MFDataChain.cc \ MFDataMember.cc \