Index: /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 1983) +++ /trunk/MagicSoft/Mars/mfilter/FilterLinkDef.h (revision 1984) @@ -20,3 +20,5 @@ #pragma link C++ class MFCT1SelFinal+; +#pragma link C++ class MFEnergySlope+; + #endif Index: /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc (revision 1984) +++ /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc (revision 1984) @@ -0,0 +1,228 @@ +/* ======================================================================== *\ +! +! * +! * 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): Antonio Stamerra 02/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! Filter to select events with a energy slope. 02/2003 +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MFEnergySlope // +// // +// A filter to select MC events (generated with a energy slope MCSlope) // +// with a different energy slope NewSlope set by the user. // +// // +// The new slope is set by the user with SetSlope(). // +// Only negative slopes are admitted; positive ones are internally // +// converted. // +// Events are selected following the new energy power slope, and the // +// sample is normalized to the number of events at: // +// 1. McMaxEnergy, if abs(NewSlope) < abs(McSlope); // +// 2. McMinEnergy, if abs(NewSlope) > abs(McSlope); // +// Mc{Max,Min}Energy are set with SetMcMinEnergy() and SetMcMaxEnergy(); // +// with GeV values. // +// Default values are the min/max energy of the MC sample. // +// // +// With this filter the statistics of the MC sample is reduced. // +// camera ver.0.6 and reflector ver.0.6 are required to fetch // +// the correct MC information // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MFEnergySlope.h" + +//#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + + +#include "MMcEvt.hxx" +#include "MMcCorsikaRunHeader.h" + +ClassImp(MFEnergySlope); + +// -------------------------------------------------------------------------- +// +// Constructor +// +MFEnergySlope::MFEnergySlope(const char *name, const char *title): + fNumSelectedEvts(0), fNewSlope(-1), fMcMinEnergy(-1.), fMcMaxEnergy(-1.) +{ + // fContName = cname; + fName = name ? name : "MFEnergySlope"; + fTitle = title ? title : "Filter to select energy with given slope"; +} + +// -------------------------------------------------------------------------- +// +// Destructor +// +MFEnergySlope::~MFEnergySlope() +{ +} + +// -------------------------------------------------------------------------- +// +// Preprocess +// +// MC slope and min/max energy are read +// Normalization factor is computed +// +Bool_t MFEnergySlope::PreProcess(MParList *pList) +{ + + MMcCorsikaRunHeader *runheader = (MMcCorsikaRunHeader*)pList->FindObject("MMcCorsikaRunHeader"); + + if (!runheader) + { + *fLog << err << dbginf << fName << " [MMcCorsikaRunHeader] not found... aborting." << endl; + return kFALSE; + } + // + // Read info from the MC sample (it must be generated with + // reflector ver.0.6 and camera ver. 0.6) + // + fMcSlope = runheader->GetSlopeSpec(); + if (fMcMinEnergy<0) + fMcMinEnergy = runheader->GetELowLim(); + if (fMcMinEnergy<0) + fMcMaxEnergy = runheader->GetEUppLim(); + + *fLog << "MFEnergySlope::PreProcess; fetched MC info:" << endl << + " E Slope: "<< fMcSlope << endl << + " E Min: "<< fMcMinEnergy << endl << + " E Max: "<< fMcMaxEnergy << endl << + " New E Slope: "<< fNewSlope << endl; + + // Slope is used with positive values in the code + if (fNewSlope < 0) + fNewSlope *= -1; + if (fMcSlope < 0) + fMcSlope *= -1; + + + // Set normalization on energy + fN0 = pow(fNewSlope>fMcSlope?fMcMinEnergy:fMcMaxEnergy,fNewSlope-fMcSlope); + + *fLog << "Normalization factor:" <FindObject("MMcEvt"); + if (!fEvt) + { + *fLog << err << dbginf << fName << " [MMcEvt] not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Select events randomly according to the MC ("old") and required ("new") +// energy slope. +// Old E slope is fMcSlope +// New E slope is set by the user (fval; fNewSlope) +// If old and new energy slope are the same skip the selection. +// The MC energy slope and lower and upper limits are taken from the +// run header (requires reflector ver.>0.6 and camera ver.>0.6) +// +Bool_t MFEnergySlope::Process() +{ + fResult = kTRUE; + + // Energy slopes are the same: skip it + if (fNewSlope == fMcSlope) + return kTRUE; + + // The value of the normalized spectrum is compared with a + // random value in [0,1]; + // if the latter is higher the event is accepted + const Float_t energy = fEvt->GetEnergy(); + + /* + // + // If energy bounds different from MC ones are selected, then + // events outside these bounds are rejected, as their slope has + // not been changed. + // + if (energy > fMcMaxEnergy || energy < fMcMinEnergy) + { + fResult = kFALSE; + return kTRUE; + } + */ + + const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope); + const Float_t Nrnd = gRandom->Uniform(); + + fResult = Nexp > Nrnd; + + if (fResult) + { + fNumSelectedEvts++; + return kTRUE; + } + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Postprocess all filters. +// +Bool_t MFEnergySlope::PostProcess() +{ + //--------------------------------- + + return kTRUE; +} + +void MFEnergySlope::StreamPrimitive(ofstream &out) const +{ + /* + if (fEvt) + fEvt->SavePrimitive(out); + + out << " MFEnergySlope " << GetUniqueName() << "("; + + if (fEvt) + out << "&" << fEvt->GetUniqueName(); + else + out << "\"" << fContName << "\""; + + out << ", New Energy Slope= " << fValue << " Normalization factor: " + << fN0 <<" );" << endl; + */ +} + + + + + + + + + Index: /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc~ =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc~ (revision 1984) +++ /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc~ (revision 1984) @@ -0,0 +1,225 @@ +/* ======================================================================== *\ +! +! * +! * 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): Antonio Stamerra 02/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! Filter to select events with a energy slope. 02/2003 +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MFEnergySlope // +// // +// A filter to select MC events (generated with a energy slope MCSlope) // +// with a different energy slope NewSlope set by the user. // +// // +// The new slope is set by the user with SetSlope(). // +// Only negative slopes are admitted; positive ones are internally // +// converted. // +// Events are selected following the new energy power slope, and the // +// sample is normalized to the number of events at: // +// 1. McMaxEnergy, if abs(NewSlope) < abs(McSlope); // +// 2. McMinEnergy, if abs(NewSlope) > abs(McSlope); // +// Mc{Max,Min}Energy are set with SetMcMinEnergy() and SetMcMaxEnergy(); // +// with GeV values. // +// Default values are the min/max energy of the MC sample. // +// // +// With this filter the statistics of the MC sample is reduced. // +// camera ver.0.6 and reflector ver.0.6 are required to fetch // +// the correct MC information // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MFEnergySlope.h" + +//#include +#include + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + + +#include "MMcEvt.hxx" +#include "MMcCorsikaRunHeader.h" + +ClassImp(MFEnergySlope); + +// -------------------------------------------------------------------------- +// +// Constructor +// +MFEnergySlope::MFEnergySlope(const char *name, const char *title): + fNumSelectedEvts(0), fNewSlope(-1), fMcMinEnergy(-1.), fMcMaxEnergy(-1.) +{ + // fContName = cname; + fName = name ? name : "MFEnergySlope"; + fTitle = title ? title : "Filter to select energy with given slope"; +} + +// -------------------------------------------------------------------------- +// +// Destructor +// +MFEnergySlope::~MFEnergySlope() +{ +} + +// -------------------------------------------------------------------------- +// +// Preprocess +// +// MC slope and min/max energy are read +// Normalization factor is computed +// +Bool_t MFEnergySlope::PreProcess(MParList *pList) +{ + + MMcCorsikaRunHeader *runheader = (MMcCorsikaRunHeader*)pList->FindObject("MMcCorsikaRunHeader"); + + if (!runheader) + { + *fLog << err << dbginf << fName << " [MMcCorsikaRunHeader] not found... aborting." << endl; + return kFALSE; + } + // + // Read info from the MC sample (it must be generated with + // reflector ver.0.6 and camera ver. 0.6) + // + fMcSlope = runheader->GetSlopeSpec(); + if (fMcMinEnergy<0) + fMcMinEnergy = runheader->GetELowLim(); + if (fMcMinEnergy<0) + fMcMaxEnergy = runheader->GetEUppLim(); + + *fLog << "MFEnergySlope::PreProcess; fetched MC info:" << endl << + " E Slope: "<< fMcSlope << endl << + " E Min: "<< fMcMinEnergy << endl << + " E Max: "<< fMcMaxEnergy << endl << + " New E Slope: "<< fNewSlope << endl; + + // Slope is used with positive values in the code + if (fNewSlope < 0) + fNewSlope *= -1; + if (fMcSlope < 0) + fMcSlope *= -1; + + + // Set normalization on energy + fN0 = pow(fNewSlope>fMcSlope?fMcMinEnergy:fMcMaxEnergy,fNewSlope-fMcSlope); + + *fLog << "Normalization factor:" <FindObject("MMcEvt"); + if (!fEvt) + { + *fLog << err << dbginf << fName << " [MMcEvt] not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Select events randomly according to the MC ("old") and required ("new") +// energy slope. +// Old E slope is fMcSlope +// New E slope is set by the user (fval; fNewSlope) +// If old and new energy slope are the same skip the selection. +// The MC energy slope and lower and upper limits are taken from the +// run header (requires reflector ver.>0.6 and camera ver.>0.6) +// +Bool_t MFEnergySlope::Process() +{ + fResult = kTRUE; + + // Energy slopes are the same: skip it + if (fNewSlope == fMcSlope) + return kTRUE; + + // The value of the normalized spectrum is compared with a + // random value in [0,1]; + // if the latter is higher the event is accepted + const Float_t energy = fEvt->GetEnergy(); + + // + // If energy bounds different from MC ones are selected, then + // events outside these bounds are rejected, as their slope has + // not been changed. + // + if (energy > fMcMaxEnergy || energy < fMcMinEnergy) + { + fResult = kFALSE; + return kTRUE; + } + const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope); + const Float_t Nrnd = gRandom->Uniform(); + + fResult = Nexp > Nrnd; + + if (fResult) + { + fNumSelectedEvts++; + return kTRUE; + } + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Postprocess all filters. +// +Bool_t MFEnergySlope::PostProcess() +{ + //--------------------------------- + + return kTRUE; +} + +void MFEnergySlope::StreamPrimitive(ofstream &out) const +{ + /* + if (fEvt) + fEvt->SavePrimitive(out); + + out << " MFEnergySlope " << GetUniqueName() << "("; + + if (fEvt) + out << "&" << fEvt->GetUniqueName(); + else + out << "\"" << fContName << "\""; + + out << ", New Energy Slope= " << fValue << " Normalization factor: " + << fN0 <<" );" << endl; + */ +} + + + + + + + + + Index: /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.h (revision 1984) +++ /trunk/MagicSoft/Mars/mfilter/MFEnergySlope.h (revision 1984) @@ -0,0 +1,74 @@ +#ifndef MARS_MFEnergySlope +#define MARS_MFEnergySlope +///////////////////////////////////////////////////////////////////////////// +// // +// MFEnergySlope // +// // +// auth. A.stamerra // +// created 30.01.03 // +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +class MParList; +class MMcEvt; +class MMcCorsikaRunHeader; + +class MFEnergySlope : public MFilter +{ +private: + + Int_t fNumSelectedEvts; + Int_t fErrors[2]; + + MMcEvt *fEvt; + // MMcCorsikaRunHeader *fMcRunHeader; + + Bool_t fResult; + Float_t fNewSlope; // New slope set by user + + Float_t fMcSlope; // Original energy slope from MC data + Float_t fMcMinEnergy; // Starting energy of MC data + Float_t fMcMaxEnergy; // Ending energy of MC data + + Float_t fN0; // Normalization factor + + // void Init(const Float_t val, + // const char *name, const char *title); + + void StreamPrimitive(ofstream &out) const; + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + +public: + + // MFEnergySlope(const char *name=NULL, const char *title=NULL); + // MFEnergySlope(); + MFEnergySlope(const char *name=NULL, const char *title=NULL); + ~MFEnergySlope(); + + + Bool_t IsExpressionTrue() const { return fResult; } + + // Slope is used with positive values in the code + void SetNewSlope(Float_t f) {fNewSlope = TMath::Abs(f);} + void SetMcSlope(Float_t f) {fMcSlope = TMath::Abs(f);} + + void SetMcMinEnergy(Float_t f) {fMcMinEnergy = f;} + void SetMcMaxEnergy(Float_t f) {fMcMaxEnergy = f;} + + ClassDef(MFEnergySlope, 0) // A Filter to select events with a given energy slope +}; + +#endif + + + + + + Index: /trunk/MagicSoft/Mars/mfilter/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mfilter/Makefile (revision 1983) +++ /trunk/MagicSoft/Mars/mfilter/Makefile (revision 1984) @@ -43,6 +43,6 @@ MFCT1SelBasic.cc \ MFCT1SelStandard.cc \ - MFCT1SelFinal.cc -# MFDataChain.cc \ + MFCT1SelFinal.cc \ + MFEnergySlope.cc SRCS = $(SRCFILES)