Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 9427) +++ /trunk/MagicSoft/Mars/Changelog (revision 9428) @@ -74,4 +74,6 @@ - added histogram for MNewImagePar - added tab "Info2" + - added all MHPhotonEvent to the paremeter list to allow access + of the ReadEnv * msim/MSimAbsorption.[h,cc]: @@ -92,4 +94,10 @@ * msimcamera/MSimAPD.[h,cc]: - take the acidental rate from a container + + * mbase/Makefile, mbase/BaseLinkDef.h: + - added MParSpline + + * mbase/MParSpline.[h,cc]: + - added Index: /trunk/MagicSoft/Mars/mbase/BaseLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/BaseLinkDef.h (revision 9427) +++ /trunk/MagicSoft/Mars/mbase/BaseLinkDef.h (revision 9428) @@ -93,4 +93,6 @@ #pragma link C++ class MParEmulated+; +#pragma link C++ class MParSpline+; + // Tool tasks #pragma link C++ class MClone+; Index: /trunk/MagicSoft/Mars/mbase/MParSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MParSpline.cc (revision 9428) +++ /trunk/MagicSoft/Mars/mbase/MParSpline.cc (revision 9428) @@ -0,0 +1,323 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of CheObs, the Modular 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 appears 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, 4/2009 +! +! Copyright: CheObs Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MParSpline +// +// Parameter container to store a MParSpline +// +////////////////////////////////////////////////////////////////////////////// +#include "MParSpline.h" + +#include // MParSpline + +#include "MLog.h" +#include "MLogManip.h" + +#include "MSpline3.h" + +ClassImp(MParSpline); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor +// +MParSpline::MParSpline(const char *name, const char *title) : fSpline(0) +{ + fName = name ? name : "MParSpline"; + fTitle = title ? title : "Parameter container to store a MSpline3"; +} + +// -------------------------------------------------------------------------- +// +// Return Xmin of the spline. 0 if fSpline==NULL +// +Double_t MParSpline::GetXmin() const +{ + return fSpline ? fSpline->GetXmin() : 0; +} + +// -------------------------------------------------------------------------- +// +// Return Xmax of the spline. 0 if fSpline==NULL +// +Double_t MParSpline::GetXmax() const +{ + return fSpline ? fSpline->GetXmax() : 0; +} + +// -------------------------------------------------------------------------- +// +// Return the width of the range in which the spline is defined. +// +Double_t MParSpline::GetWidth() const +{ + return fSpline ? fSpline->GetXmax()-fSpline->GetXmin() : 0; +} + +// -------------------------------------------------------------------------- +// +// Return the result of the spline at val. +// The user is reposible to make sure that fSpline is valid (!=NULL) +// +Double_t MParSpline::Eval(Double_t val) const +{ + return fSpline->Eval(val); +} + +// -------------------------------------------------------------------------- +// +// Delete fSpline and set to 0. +// +void MParSpline::Clear(Option_t *) +{ + if (fSpline) + delete fSpline; + fSpline=0; +} + +// -------------------------------------------------------------------------- +// +// Read a TGraph from a file and return a newly allocated MSpline3. +// +MSpline3 *MParSpline::ReadSpline(const char *fname) const +{ + *fLog << inf << "Reading spline from " << fname << endl; + + TGraph g(fname); + if (g.GetN()==0) + { + *fLog << err << "ERROR - No data points from " << fname << "." << endl; + return 0; + } + + // option: b1/e1 b2/e2 (first second derivative?) + // option: valbeg/valend (first second derivative?) + + MSpline3 *s = new MSpline3(g); + s->SetTitle(fname); + + return s; +} + +// -------------------------------------------------------------------------- +// +// Initializes a spline from min to max with n points with 1. +// +/* +void MParSpline::InitUnity(UInt_t n, Float_t min, Float_t max) +{ + // Delete the existing spline + Clear(); + + // We need at least two points (the edges) + if (n<2) + return; + + // Initialize an array with the x-values + const MBinning bins(n-1, min, max); + + // Initialize an array with all one + MArrayD y(n); + y.Reset(1); + + // Set the spline + fSpline = new MSpline3(bins.GetEdges(), y.GetArray(), n); +} +*/ + +// -------------------------------------------------------------------------- +// +// Takes the existing fSpline and multiplies it with the given spline. +// As x-points the values from fSpline are used. +// +void MParSpline::Multiply(const TSpline3 &spline) +{ + if (!fSpline) + { + Clear(); + // WARNING WARNING WARNING: This is a very dangerous cast! + fSpline = (MSpline3*)spline.Clone(); + return; + } + + // Initialize a TGraph with the number of knots from a TSpline + TGraph g(fSpline->GetNp()); + + // Loop over all knots + for (int i=0; iGetNp(); i++) + { + // Get th i-th knot + Double_t x, y; + fSpline->GetKnot(i, x, y); + + // Multiply y by the value from the given spline + y *= spline.Eval(x); + + // push the point "back" + g.SetPoint(i, x, y); + } + + // Clear the spline and initialize a new one from the updated TGraph + Clear(); + fSpline = new MSpline3(g); +} + +// -------------------------------------------------------------------------- +// +// Initializes a TSpline3 with n knots x and y. Call Multiply for it. +// +void MParSpline::Multiply(UInt_t n, const Double_t *x, const Double_t *y) +{ + const TSpline3 spline("Spline", const_cast(x), const_cast(y), n); + Multiply(spline); +} + +// -------------------------------------------------------------------------- +// +// Read a Spline from a file using ReadSpline. +// Call Multiply for it. +// +void MParSpline::Multiply(const char *fname) +{ + const MSpline3 *spline = ReadSpline(fname); + if (!spline) + return; + + Multiply(*spline); + + delete spline; +} + +// -------------------------------------------------------------------------- +// +// Read a spline from a file and set fSpline accfordingly. +// +Bool_t MParSpline::ReadFile(const char *fname) +{ + MSpline3 *spline = ReadSpline(fname); + if (!spline) + return kFALSE; + + + // option: b1/e1 b2/e2 (first second derivative?) + // option: valbeg/valend (first second derivative?) + + Clear(); + fSpline = spline; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the spline to a function. For details see MSpline3 +// +void MParSpline::SetFunction(const TF1 &f) +{ + // FIXME: Use TF1 directly? (In most cases this seems to be slower) + + // option: b1/e1 b2/e2 (first second derivative?) + // option: valbeg/valend (first second derivative?) + + // if (f.GetNpar()==0) + // No SUPPORT + + Clear(); + fSpline = new MSpline3(f);//, fRunHeader->GetFreqSampling()/1000.); + fSpline->SetTitle(f.GetTitle()); + + // FIXME: Make a boundary condition e1b1,0,0 (First der, at Xmin and Xmax==0) + // FIXME: Force the spline to be 0 at Xmin and Xmax? +} + +Bool_t MParSpline::SetFunction(const char *func, Int_t n, Double_t xmin, Double_t xmax) +{ + // FIXME: Use TF1 directly? (In most cases this seems to be slower) + TF1 f("f", func, xmin, xmax); + f.SetNpx(n); + + SetFunction(f); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// FileName: reflectivity.txt +// Function.Name: exp(-(x/2)^2/2) +// Function.Npx: 25 +// Function.Xmin: -8 +// Function.Xmax: 8 +// +Int_t MParSpline::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "FileName", print)) + { + rc = kTRUE; + if (!ReadFile(GetEnvValue(env, prefix, "FileName", ""))) + return kERROR; + } + + if (IsEnvDefined(env, prefix, "Function.Name", print)) + { + rc = kTRUE; + + Int_t npx = 25; + Float_t xmin = -8; + Float_t xmax = 8; + TString func = "exp(-(x/2)^2/2)"; + + if (IsEnvDefined(env, prefix, "Function.Npx", print)) + npx = GetEnvValue(env, prefix, "Function.Npx", npx); + if (IsEnvDefined(env, prefix, "Function.Xmin", print)) + xmin = GetEnvValue(env, prefix, "Function.Xmin", xmin); + if (IsEnvDefined(env, prefix, "Function.Xmax", print)) + xmax = GetEnvValue(env, prefix, "Function.Xmax", xmax); + + SetFunction(GetEnvValue(env, prefix, "Function.Name", func), npx, xmin, xmax); + } + return rc; +} + +void MParSpline::Paint(Option_t *) +{ + fSpline->SetMarkerStyle(kFullDotMedium); + + if (!fSpline->GetHistogram()) + fSpline->Paint(); + + TH1 *h = fSpline->GetHistogram(); + if (!h) + return; + + //h->SetXTitle("t [ns]"); + //h->SetYTitle("a.u."); + + fSpline->Paint("PC"); +} Index: /trunk/MagicSoft/Mars/mbase/MParSpline.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MParSpline.h (revision 9428) +++ /trunk/MagicSoft/Mars/mbase/MParSpline.h (revision 9428) @@ -0,0 +1,55 @@ +#ifndef MARS_MParSpline +#define MARS_MParSpline + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class TF1; // MParSpline +class TSpline3; +class MSpline3; + +class MParSpline : public MParContainer +{ +private: + MSpline3 *fSpline; // Spline hold by this container + + // MParSpline + MSpline3 *ReadSpline(const char *fname) const; + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + +public: + MParSpline(const char *name=0, const char *title=0); + + // MParSpline + const MSpline3 *GetSpline() const { return fSpline; } + + Double_t GetXmin() const; + Double_t GetXmax() const; + Double_t GetWidth() const; + + Bool_t IsValid() const { return fSpline ? kTRUE : kFALSE; } + + Double_t Eval(Double_t val) const; + + void SetFunction(const TF1 &f); + Bool_t SetFunction(const char *func, Int_t n, Double_t xmin, Double_t xmax); + + //void InitUnity(UInt_t n, Float_t min, Float_t max); + + void Multiply(const TSpline3 &spline); + void Multiply(UInt_t n, const Double_t *x, const Double_t *y); + void Multiply(const char *fname); + + Bool_t ReadFile(const char *fname); + + // TObject + void Clear(Option_t *o=""); + void Paint(Option_t *o=""); + + ClassDef(MParSpline, 1) // Container to hold a MSpline3 +}; + +#endif Index: /trunk/MagicSoft/Mars/mbase/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mbase/Makefile (revision 9427) +++ /trunk/MagicSoft/Mars/mbase/Makefile (revision 9428) @@ -31,4 +31,5 @@ MLogPlugin.cc \ MLogHtml.cc \ + MParSpline.cc \ MParList.cc \ MParContainer.cc \ Index: /trunk/MagicSoft/Mars/mjobs/MJSimulation.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJSimulation.cc (revision 9427) +++ /trunk/MagicSoft/Mars/mjobs/MJSimulation.cc (revision 9428) @@ -430,12 +430,20 @@ filled.SetNameTab("MedDev", "Time between first and last photon hitting a detector"); - MHPhotonEvent planeG(1); // Get from MaxImpact - MHPhotonEvent plane0(2); // Get from MReflector - MHPhotonEvent plane1(2); - MHPhotonEvent plane2(2); - MHPhotonEvent plane3(2); - MHPhotonEvent plane4(2); - MHPhotonEvent planeF1(3); // Get from MGeomCam - MHPhotonEvent planeF2(header.IsPointRun()?4:3); // Get from MGeomCam + MHPhotonEvent planeG(1, "HPhotonEventGround"); // Get from MaxImpact + MHPhotonEvent plane0(2, "HMirrorPlane0"); // Get from MReflector + //MHPhotonEvent plane1(2, "HMirrorPlane1"); + MHPhotonEvent plane2(2, "HMirrorPlane2"); + MHPhotonEvent plane3(2, "HMirrorPlane3"); + MHPhotonEvent plane4(2, "HMirrorPlane4"); + MHPhotonEvent planeF1(3, "HPhotonEventCamera"); // Get from MGeomCam + MHPhotonEvent planeF2(header.IsPointRun()?4:3, "HPhotonEventCones"); // Get from MGeomCam + + plist.AddToList(&planeG); + plist.AddToList(&plane0); + plist.AddToList(&plane2); + plist.AddToList(&plane3); + plist.AddToList(&plane4); + plist.AddToList(&planeF1); + plist.AddToList(&planeF2);; //MHPSF psf;