Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 5334) +++ trunk/MagicSoft/Mars/Changelog (revision 5335) @@ -20,4 +20,15 @@ -*-*- END OF LINE -*-*- + 2004/11/02: Abelardo Moralejo + + * mfilter/MFSize.[h,cc] + - added. A filter to modify the Size spectrum of a sample through + cutting out events as a function of their Size. + + * mfilter/Makefile, FilterLinkDef.h + - added new class MFSize. + + * mtemp/mpadova/macros/resize.C + - Added. Example of how to use the new filter. 2004/11/01: Markus Gaug Index: trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 5334) +++ trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 5335) @@ -21,4 +21,6 @@ #pragma link C++ class MFCosmics+; +#pragma link C++ class MFSize+; + #pragma link C++ class MFEnergySlope+; Index: trunk/MagicSoft/Mars/mfilter/MFSize.cc =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSize.cc (revision 5335) +++ trunk/MagicSoft/Mars/mfilter/MFSize.cc (revision 5335) @@ -0,0 +1,158 @@ +/* ======================================================================== *\ +! +! * +! * 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): Abelardo Moralejo 10/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MFSize // +// // +// A filter to select events as a function of Size. Depending on the Size // +// value, a certain fraction of events is rejected. A histogram (fProb) // +// must be provided, containing the acceptance probability of events as a // +// function of log10(Size) (which must be the x axis of the histogram). // +// There is no restriction on the histogram binning or ranges (events // +// with Size outside histogram ranges will simply be accepted). Histogram // +// bin contents should be between 0. and 1. If they are not, negative // +// will be equivalent to 0 (events always rejected), and values above 1 // +// will be equivalent to 1 (events always accepted). A warning will be // +// displayed by MFSize::SetProb in that case. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MFSize.h" + +#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHillas.h" + +ClassImp(MFSize); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Constructor +// +MFSize::MFSize(const char *name, const char *title): fNumSelectedEvts(0) +{ + fName = name ? name : "MFSize"; + fTitle = title ? title : "Filter to select events according to Size"; + + // + // Default fProb histogram: 20 bins in the log10(Size) range [2,4] + // Values are calculated to make the Size distribution of MC gammas + // in the standard MC sample to look like the distribution for MC + // hadrons (p+He) + // + fProb = new TH1F("ProbVsLog10Size", "", 20, 2., 4.); + + Float_t frac[20] = {0.144928, 0.171911, 0.18763, 0.209461, 0.253219, + 0.305425, 0.384593, 0.485204, 0.612452, 0.708138, + 0.754522, 0.728028, 0.774046, 0.791422, 0.808775, + 0.896842, 1., 1., 1., 1.}; + + for (Int_t i = 0; i < 20; i++) + fProb->SetBinContent(i+1, frac[i]); +} + +// -------------------------------------------------------------------------- +// +// SetProb: Sets the histogram containing the probability acceptance for +// events as a function of log10(Size). Checks that the values are between +// 0 and 1. +// + +void MFSize::SetProb(TH1F* h) +{ + fProb = h; + + for (Int_t i = 1; i <= fProb->GetNbinsX(); i++) + { + if (fProb->GetBinContent(i) > 1.) + *fLog << warn << dbginf << fName << + "Probability value above 1 found in histogram. " << + "Value will be interpreted as 1." << endl; + + else if(fProb->GetBinContent(i) < 0. ) + *fLog << warn << dbginf << fName << + "Probability value below 0 found in histogram. " << + "Value will be interpreted as 0." << endl; + } +} + +// -------------------------------------------------------------------------- +// +// Preprocess +// +// Look for the MHillas container +// +Int_t MFSize::PreProcess(MParList *pList) +{ + fHillas = (MHillas*)pList->FindObject(AddSerialNumber("MHillas")); + if (!fHillas) + { + *fLog << err << dbginf << fName << AddSerialNumber("MHillas") << + " not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Process +// +Int_t MFSize::Process() +{ + fResult = kTRUE; + + Int_t ibin = fProb->FindBin(log10(fHillas->GetSize())); + + // + // If value is outside histogram range, accept event + // + if (ibin > fProb->GetNbinsX() || ibin < 1) + { + fNumSelectedEvts++; + return kTRUE; + } + + Float_t SelFrac = fProb->GetBinContent(ibin); + + const Float_t Nrnd = gRandom->Uniform(); + + fResult = SelFrac > Nrnd; + + if (!fResult) + return kTRUE; + + fNumSelectedEvts++; + return kTRUE; +} + Index: trunk/MagicSoft/Mars/mfilter/MFSize.h =================================================================== --- trunk/MagicSoft/Mars/mfilter/MFSize.h (revision 5335) +++ trunk/MagicSoft/Mars/mfilter/MFSize.h (revision 5335) @@ -0,0 +1,47 @@ +#ifndef MARS_MFSize +#define MARS_MFSize +///////////////////////////////////////////////////////////////////////////// +// // +// MFSize // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +#include + +class MParList; +class MHillas; + +class MFSize : public MFilter +{ +private: + Int_t fNumSelectedEvts; // counter for number of selected events + + MHillas *fHillas; // Hillas parameters container + TH1F* fProb; // Acceptance probability vs log10 size + + Bool_t fResult; // Result returned by IsExpressionTrue + + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MFSize(const char *name=NULL, const char *title=NULL); + + void SetProb(TH1F* h); + + Bool_t IsExpressionTrue() const { return fResult; } + + ClassDef(MFSize, 0) // A Filter to select events as a function of Size +}; + +#endif + + + + + + Index: trunk/MagicSoft/Mars/mfilter/Makefile =================================================================== --- trunk/MagicSoft/Mars/mfilter/Makefile (revision 5334) +++ trunk/MagicSoft/Mars/mfilter/Makefile (revision 5335) @@ -32,5 +32,6 @@ MFSoftwareTrigger.cc \ MFCosmics.cc \ - MFEnergySlope.cc + MFEnergySlope.cc \ + MFSize.cc ############################################################ Index: trunk/MagicSoft/Mars/mtemp/mpadova/macros/resize.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mpadova/macros/resize.C (revision 5335) +++ trunk/MagicSoft/Mars/mtemp/mpadova/macros/resize.C (revision 5335) @@ -0,0 +1,96 @@ +/* ======================================================================== *\ +! +! * +! * 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. +! * +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + + +void resize() +{ + MParList plist; + MTaskList tlist; + + plist.AddToList(&tlist); + + MReadMarsFile read("Events", "star_gamma_train.root"); + + read.DisableAutoScheme(); + tlist.AddToList(&read); + + MWriteRootFile write("star_gamma_train_new.root", "recreate"); + + write.AddContainer("MMcEvt", "Events"); + write.AddContainer("MHillas", "Events"); + write.AddContainer("MHillasExt", "Events"); + write.AddContainer("MImagePar", "Events"); + write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MHillasSrc", "Events"); + + write.AddContainer("MGeomCam", "RunHeaders"); + write.AddContainer("MMcConfigRunHeader", "RunHeaders"); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders"); + write.AddContainer("MMcFadcHeader", "RunHeaders"); + write.AddContainer("MMcTrigHeader", "RunHeaders"); + + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MSrcPosCam", "RunHeaders"); + write.AddContainer("MMcRunHeader", "RunHeaders"); + + + // + // Histogram: event acceptance probability vs log10(Size). Binning + // and range are completely free. Beyond the histogram limits all + // events are accepted + // + + TH1F* h = new TH1F("h", "", 20, 2., 4.); + + Float_t frac[20] = {0.144928, 0.171911, 0.18763, 0.209461, 0.253219, + 0.305425, 0.384593, 0.485204, 0.612452, 0.708138, + 0.754522, 0.728028, 0.774046, 0.791422, 0.808775, + 0.896842, 1., 1., 1., 1.}; + + for (Int_t i = 0; i < 20; i++) + h->SetBinContent(i+1, frac[i]); + + MFSize fsize; + fsize.SetProb(h); + + write.SetFilter(&fsize); + + tlist.AddToList(&fsize); + tlist.AddToList(&write); + + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + MProgressBar *bar = new MProgressBar; + evtloop.SetProgressBar(bar); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + bar->DestroyWindow(); + + return; +}