Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 9235) +++ trunk/MagicSoft/Mars/Changelog (revision 9236) @@ -23,6 +23,7 @@ * mbase/BaseLinkDef.h, mbase/Makefile: - added MQuaternion - - * mbase/MQuaternion.[h,cc]: + - added MReflection + + * mbase/MQuaternion.[h,cc], mbase/MReflection.[h,cc]: - added @@ -56,4 +57,9 @@ - we don't have to use a MHexagon anymore caluclating DistanceToPrimitive + + * msimreflector/MMirror.[h,cc], msimreflector/MMirrorHex.[h,cc], + msimreflector/MMirrorSquare.[h,cc], msimreflector/MMirrorDisk.[h,cc], + msimreflector/MReflector.[h,cc], msimreflector/MSimReflector.[h,cc]: + - added Index: trunk/MagicSoft/Mars/msimreflector/MMirror.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirror.cc (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirror.cc (revision 9236) @@ -0,0 +1,134 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2008 +! +! Copyright: Software Development, 2000-2008 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMirror +// +// Note, that we could use basic geometry classes instead, but especially +// CanHit is time critical. So this class is (should be) optimized for +// execution speed. +// +// This base class provides the code to calculate a spherical mirror +// (ExecuteMirror) and to scatter in a way to get a proper PSF. +// Furthermore it stored the geometry of a mirror. +// +// ------------------------------------------------------------------------ +// +// Bool_t CanHit(const MQuaternion &p) const; +// +// This is a very rough estimate of whether a photon at a position p +// can hit a mirror. The position might be off in z and the photon +// still has to follow its trajectory. Nevertheless we can fairly assume +// the the way to travel in x/y is pretty small so we can give a rather +// good estimate of whether the photon can hit. +// +// Never throw away a photon whihc can hit the mirror! +// +// ------------------------------------------------------------------------ +// +// Bool_t HasHit(const MQuaternion &p) const; +// +// Check if the given position coincides with the mirror. The position +// is assumed to be the incident point on the mirror's surface. +// +// The coordinates are in the mirrors coordinate frame. +// +// The action should coincide with what is painted in Paint() +// +// ------------------------------------------------------------------------ +// +// void Paint(Option_t *opt) +// +// Paint the mirror in x/y. +// +// The graphic should coincide with the action in HasHit +// +// ------------------------------------------------------------------------ +// +// Int_t ReadM(const TObjArray &tok); +// +// Read the mirror's setup from a file. The first eight tokens should be +// ignored. (This could be fixed!) +// +////////////////////////////////////////////////////////////////////////////// +#include "MMirror.h" + +#include + +#include "MQuaternion.h" + +ClassImp(MMirror); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Return the TVector2 which is the x/y position of the mirror minus +// q.XYvector/(; +// +TVector2 MMirror::operator-(const MQuaternion &q) const +{ + return TVector2(X()-q.X(), Y()-q.Y()); +} + +// -------------------------------------------------------------------------- +// +// Return the TVector2 which is the difference of this mirror and the +// given mirror +// +TVector2 MMirror::operator-(const MMirror &m) const +{ + return TVector2(X()-m.X(), Y()-m.Y()); +} + +// -------------------------------------------------------------------------- +// +// Simulate the PSF. Therefor we smear out the given normal vector +// with a gaussian. +// +// Returns a vector which can be added to the normal vector. +// +// FIXME: What is the correct focal distance to be given here? +// Can the smearing be imporved? +// +TVector3 MMirror::SimPSF(const TVector3 &n, Double_t F, Double_t psf) const +{ + //const TVector3 n( x, y, -d) // Normal vector of the mirror + const TVector3 xy(-n.Y(), n.X(), 0); // Normal vector in x/y plane + + Double_t gx, gy; + gRandom->Rannor(gx, gy); // 2D random Gauss distribution + + psf /= 2; // The factor two because of the doubleing of the angle in the reflection + psf /= F; // Scale the Gauss to the size of the PSF + psf *= n.Z(); // Normalize the addon vector to the normal vector + //psf *= n.Mag(); // Alternative! (Gaussian projected on the surface of a sphere) + + TVector3 dn(gx*psf, gy*psf, 0); // Instead of psf/F also atan(psf/F) might make sense + + dn.Rotate(-n.Theta(), xy); // Tilt the gauss-vector to the normal vector + + return dn; // Return the vector to be added to the normal vector +} Index: trunk/MagicSoft/Mars/msimreflector/MMirror.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirror.h (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirror.h (revision 9236) @@ -0,0 +1,91 @@ +#ifndef MARS_MMirror +#define MARS_MMirror + +#ifndef ROOT_TRotation +#include +#endif + +class MQuaternion; + +class MMirror : public TObject +{ + friend void operator/=(MQuaternion &, const MMirror &); + friend void operator*=(MQuaternion &, const MMirror &); + friend void operator-=(MQuaternion &, const MMirror &); + friend void operator+=(MQuaternion &, const MMirror &); + +private: + TVector3 fPos; + TVector3 fNorm; // faster without + + TRotation fTilt; + + // ----- Spherical mirror data members ----- + Double_t fFocalLength; + Double_t fSigmaPSF; + + // MMirror *fNeighbors[964]; + +public: + MMirror() : fSigmaPSF(-1) + { + } + + // ----- Mirror basic functions ----- + TVector2 operator-(const MQuaternion &q) const;// { return TVector2(X()-q.X(), Y()-q.Y()); } + TVector2 operator-(const MMirror &m) const;// { return TVector2(X()-m.X(), Y()-m.Y()); } + + void SetPosition(const TVector3 &v) { fPos = v; } + void SetNorm(const TVector3 &n) { + fNorm = n; + + fTilt = TRotation(); + // Can be simplified?? rotate the mirror along + // perpendicular to its normal projected to x/y and z + // by its "zenith angle" + fTilt.Rotate(-n.Theta(), TVector3(-n.Y(), n.X(), 0)); + } + + Double_t X() const { return fPos.X(); } + Double_t Y() const { return fPos.Y(); } + + TVector2 GetPosXY() const { return fPos.XYvector(); } + const TVector3 &GetPos() const { return fPos; } + const TVector3 &GetNorm() const { return fNorm; } + + Double_t GetDist() const { return fPos.Perp(); } + + virtual Double_t GetMaxR() const=0;// { return TMath::Max(fMaxRX, fMaxRY); } + + TVector3 SimPSF(const TVector3 &n, Double_t F, Double_t psf) const; + + Bool_t ExecuteMirror(MQuaternion &p, MQuaternion &u) const; + + // ----- Basic function for parabolic mirror ----- + void SetFocalLength(Double_t f) { fFocalLength = f; } + Double_t GetFocalLength() const { return fFocalLength; } + + Bool_t ExecuteReflection(MQuaternion &p, MQuaternion &u) const; + + // ----- Mirror specialized functions ----- + + // This should fit with Paint() + virtual Bool_t HasHit(const MQuaternion &p) const=0; + // This should fit with HasHit() + //void Paint(Option_t *)=0; + virtual Bool_t CanHit(const MQuaternion &p) const=0; + + virtual Int_t ReadM(const TObjArray &tok)=0; + + /* + Bool_t IsSortable() const { return kTRUE; } + Int_t Compare(const TObject *obj) const + { + MMirror &m = (MMirror&)*obj; + return m.fPos.Mag2() +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMirrorDisk +// +// A disk like spherical mirror. +// +////////////////////////////////////////////////////////////////////////////// +#include "MMirrorDisk.h" + +#include +#include + +#include "MQuaternion.h" + +ClassImp(MMirrorDisk); + +using namespace std; + +// ------------------------------------------------------------------------ +// +// This is a very rough estimate of whether a photon at a position p +// can hit a mirror. The position might be off in z and the photon +// still has to follow its trajectory. Nevertheless we can fairly assume +// the the way to travel in x/y is pretty small so we can give a rather +// good estimate of whether the photon can hit. +// +// never throw away a photon whihc can hit the mirror! +// +Bool_t MMirrorDisk::CanHit(const MQuaternion &p) const +{ + // p is given in the reflectors coordinate frame. This is meant + // to be a fast check to sort out all mirrors which we can omit + // without time consuming calculations. + + // Check if this mirror can be hit at all + const Double_t dx = TMath::Abs(p.X()-X()); + if (dx>fR*1.05) + return kFALSE; + + const Double_t dy = TMath::Abs(p.Y()-Y()); + if (dy>fR*1.05) + return kFALSE; + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +// Check if the given position coincides with the mirror. The position +// is assumed to be the incident point on the mirror's surface. +// +// The coordinates are in the mirrors coordinate frame. +// +// The action should coincide with what is painted in Paint() +// +Bool_t MMirrorDisk::HasHit(const MQuaternion &p) const +{ + // p is the incident point in the mirror in the mirror's + // coordinate frame + + // Black spot in the mirror center (here we can fairly ignore + // the distance from the plane to the mirror surface, as long + // as the black spot does not become too large) + if (p.R2()<0.5*0.5) + return kFALSE; + + // Check if the photon has really hit the square mirror + return p.R()GetFillColor()); + + if (!TString(opt).Contains("border", TString::kIgnoreCase)) + e.PaintEllipse(X(), Y(), 0.5, 0.5, 0, 360, 0); +} + +// ------------------------------------------------------------------------ +// +// Read the mirror's setup from a file. The first eight tokens should be +// ignored. (This could be fixed!) +// +// Here we read: fR +// +Int_t MMirrorDisk::ReadM(const TObjArray &tok) +{ + if (tok.GetSize()<9) + return -1; + + Double_t r = atof(tok[8]->GetName()); + + if (r<=0) + return -1; + + fR = r; + + return 1; +} Index: trunk/MagicSoft/Mars/msimreflector/MMirrorDisk.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirrorDisk.h (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirrorDisk.h (revision 9236) @@ -0,0 +1,28 @@ +#ifndef MARS_MMirrorDisk +#define MARS_MMirrorDisk + +#ifndef MARS_MMirror +#include "MMirror.h" +#endif + +class MMirrorDisk : public MMirror +{ +private: + Double_t fR; + +public: + MMirrorDisk() : fR(24.75) { } + + Double_t GetMaxR() const { return fR; } + + // This should fit with Paint() + Bool_t HasHit(const MQuaternion &p) const; + // This should fit with HasHit() + void Paint(Option_t *); + Bool_t CanHit(const MQuaternion &p) const; + Int_t ReadM(const TObjArray &tok); + + ClassDef(MMirrorDisk, 0) // A spherical disk type mirror +}; + +#endif Index: trunk/MagicSoft/Mars/msimreflector/MMirrorHex.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirrorHex.cc (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirrorHex.cc (revision 9236) @@ -0,0 +1,157 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2008 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMirrorHex +// +// A hexagonal spherical mirror +// +////////////////////////////////////////////////////////////////////////////// +#include "MMirrorHex.h" + +#include +#include + +#include "MHexagon.h" +#include "MQuaternion.h" + +ClassImp(MMirrorHex); + +using namespace std; + +const Double_t MMirrorHex::fgCos30 = TMath::Cos(30*TMath::DegToRad()); +const Double_t MMirrorHex::fgCos60 = TMath::Cos(60*TMath::DegToRad()); +const Double_t MMirrorHex::fgSin60 = TMath::Sin(60*TMath::DegToRad()); + +// ------------------------------------------------------------------------ +// +// This is a very rough estimate of whether a photon at a position p +// can hit a mirror. The position might be off in z and the photon +// still has to follow its trajectory. Nevertheless we can fairly assume +// the the way to travel in x/y is pretty small so we can give a rather +// good estimate of whether the photon can hit. +// +// never throw away a photon whihc can hit the mirror! +// +Bool_t MMirrorHex::CanHit(const MQuaternion &p) const +{ + // p is given in the reflectors coordinate frame. This is meant + // to be a fast check to sort out all mirrors which we can omit + // without time consuming calculations. + + return TMath::Hypot(p.X()-X(), p.Y()-X())<1.05*fMaxR; +} + +// ------------------------------------------------------------------------ +// +// Check if the given position coincides with the mirror. The position +// is assumed to be the incident point on the mirror's surface. +// +// The coordinates are in the mirrors coordinate frame. +// +// The action should coincide with what is painted in Paint() +// +Bool_t MMirrorHex::HasHit(const MQuaternion &p) const +{ + // p is the incident point in the mirror in the mirror's + // coordinate frame + + // Black spot in the mirror center (here we can fairly ignore + // the distance from the plane to the mirror surface, as long + // as the black spot does not become too large) + if (p.R2()<0.5*0.5) + return kFALSE; + + // + // Now check if point is outside of hexagon; just check x coordinate + // in three coordinate systems: the default one, in which two sides of + // the hexagon are paralel to the y axis (see camera displays) and two + // more, rotated with respect to that one by +- 60 degrees. + // + if (TMath::Abs(X())>fD/*/2*/) + return kFALSE; + + const Double_t dxc = p.X()*fgCos60; + const Double_t dys = p.Y()*fgSin60; + + if (TMath::Abs(dxc+dys)>fD/*/2*/) + return kFALSE; + + if (TMath::Abs(dxc-dys)>fD/*/2*/) + return kFALSE; + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +// Paint the mirror in x/y. +// +// The graphic should coincide with the action in HasHit +// +void MMirrorHex::Paint(Option_t *opt) +{ + MHexagon h; + TEllipse e; + if (!TString(opt).Contains("line", TString::kIgnoreCase)) + { + h.SetLineStyle(0); + h.SetFillColor(17); + e.SetLineStyle(0); + e.SetFillColor(gPad->GetFillColor()); + } + else + { + h.SetFillStyle(0); + e.SetFillStyle(0); + } + + h.PaintHexagon(X(), Y(), fD*2); + + if (!TString(opt).Contains("border", TString::kIgnoreCase)) + e.PaintEllipse(X(), Y(), 0.5, 0.5, 0, 360, 0); +} + +// ------------------------------------------------------------------------ +// +// Read the mirror's setup from a file. The first eight tokens should be +// ignored. (This could be fixed!) +// +// Here we read: D +// +Int_t MMirrorHex::ReadM(const TObjArray &tok) +{ + if (tok.GetSize()<9) + return -1; + + const Double_t d = atof(tok[8]->GetName()); + + if (d<=0) + return -1; + + SetD(d); + + return 1; +} Index: trunk/MagicSoft/Mars/msimreflector/MMirrorHex.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirrorHex.h (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirrorHex.h (revision 9236) @@ -0,0 +1,35 @@ +#ifndef MARS_MMirrorHex +#define MARS_MMirrorHex + +#ifndef MARS_MMirror +#include "MMirror.h" +#endif + +class MMirrorHex : public MMirror +{ +private: + const static Double_t fgCos30; + const static Double_t fgCos60; + const static Double_t fgSin60; + + Double_t fD; // half diameter D or better distance between opposite sides + Double_t fMaxR; + + void SetD(Double_t d) { fD=d/2; fMaxR=fD/fgCos30; } + +public: + MMirrorHex() { SetD(24.75); } + + Double_t GetMaxR() const { return fMaxR; } + + // This should fit with Paint() + Bool_t HasHit(const MQuaternion &p) const; + // This should fit with HasHit() + void Paint(Option_t *); + Bool_t CanHit(const MQuaternion &p) const; + Int_t ReadM(const TObjArray &tok); + + ClassDef(MMirrorHex, 0) // A spherical hexagon type mirror +}; + +#endif Index: trunk/MagicSoft/Mars/msimreflector/MMirrorSquare.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirrorSquare.cc (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirrorSquare.cc (revision 9236) @@ -0,0 +1,149 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2008 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMirror +// +// A square type spherical mirror +// +////////////////////////////////////////////////////////////////////////////// +#include "MMirrorSquare.h" + +#include +#include +#include + +#include "MQuaternion.h" + +ClassImp(MMirrorSquare); + +using namespace std; + +// ------------------------------------------------------------------------ +// +// This is a very rough estimate of whether a photon at a position p +// can hit a mirror. The position might be off in z and the photon +// still has to follow its trajectory. Nevertheless we can fairly assume +// the the way to travel in x/y is pretty small so we can give a rather +// good estimate of whether the photon can hit. +// +// never throw away a photon whihc can hit the mirror! +// +Bool_t MMirrorSquare::CanHit(const MQuaternion &p) const +{ + // p is given in the reflectors coordinate frame. This is meant + // to be a fast check to sort out all mirrors which we can omit + // without time consuming calculations. + + // Check if this mirror can be hit at all + const Double_t dx = TMath::Abs(p.X()-X()); + if (dx>fSideLength*1.05) + return kFALSE; + + const Double_t dy = TMath::Abs(p.Y()-Y()); + if (dy>fSideLength*1.05) + return kFALSE; + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +// Check if the given position coincides with the mirror. The position +// is assumed to be the incident point on the mirror's surface. +// +// The coordinates are in the mirrors coordinate frame. +// +// The action should coincide with what is painted in Paint() +// +Bool_t MMirrorSquare::HasHit(const MQuaternion &p) const +{ + // p is the incident point in the mirror in the mirror's + // coordinate frame + + // Black spot in the mirror center (here we can fairly ignore + // the distance from the plane to the mirror surface, as long + // as the black spot does not become too large) + if (p.R2()<0.5*0.5) + return kFALSE; + + // Check if the photon has really hit the square mirror + if (TMath::Max(fabs(p.X()), fabs(p.Y()))>fSideLength) + return kFALSE; + + return kTRUE; +} + +// ------------------------------------------------------------------------ +// +// Paint the mirror in x/y. +// +// The graphic should coincide with the action in HasHit +// +void MMirrorSquare::Paint(Option_t *opt) +{ + TBox b; + TEllipse e; + + if (!TString(opt).Contains("line", TString::kIgnoreCase)) + { + b.SetFillColor(17); + b.SetLineStyle(0); + e.SetLineStyle(0); + e.SetFillColor(gPad->GetFillColor()); + } + else + { + b.SetFillStyle(0); + e.SetFillStyle(0); + } + + b.PaintBox(X()-fSideLength, Y()-fSideLength, X()+fSideLength, Y()+fSideLength); + + if (!TString(opt).Contains("border", TString::kIgnoreCase)) + e.PaintEllipse(X(), Y(), 0.5, 0.5, 0, 360, 0); +} + +// ------------------------------------------------------------------------ +// +// Read the mirror's setup from a file. The first eight tokens should be +// ignored. (This could be fixed!) +// +// Here we read: L +// +Int_t MMirrorSquare::ReadM(const TObjArray &tok) +{ + if (tok.GetSize()<9) + return -1; + + Double_t l = atof(tok[8]->GetName()); + + if (l<=0) + return -1; + + fSideLength = l/2; + + return 1; +} Index: trunk/MagicSoft/Mars/msimreflector/MMirrorSquare.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MMirrorSquare.h (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MMirrorSquare.h (revision 9236) @@ -0,0 +1,27 @@ +#ifndef MARS_MMirrorSquare +#define MARS_MMirrorSquare + +#ifndef MARS_MMirror +#include "MMirror.h" +#endif + +class MMirrorSquare : public MMirror +{ +private: + Double_t fSideLength; // HALF of the side length! +public: + MMirrorSquare() : fSideLength(24.75) { } + + Double_t GetMaxR() const { return TMath::Sqrt(2.)*fSideLength; } + + // This should fit with Paint() + Bool_t HasHit(const MQuaternion &p) const; + // This should fit with HasHit() + void Paint(Option_t *); + Bool_t CanHit(const MQuaternion &p) const; + Int_t ReadM(const TObjArray &tok); + + ClassDef(MMirrorSquare, 0) // A spherical square type mirror +}; + +#endif Index: trunk/MagicSoft/Mars/msimreflector/MReflector.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MReflector.cc (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MReflector.cc (revision 9236) @@ -0,0 +1,263 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2008 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MReflector +// +////////////////////////////////////////////////////////////////////////////// +#include "MReflector.h" + +#include +#include + +#include + +#include "MQuaternion.h" +#include "MMirror.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MReflector); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor +// +MReflector::MReflector(const char *name, const char *title) +{ + fName = name ? name : "MReflector"; + fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)"; + + fMirrors.SetOwner(); +} + +// -------------------------------------------------------------------------- +// +// Calculate the maximum radius of th ereflector. This is not meant as +// a precise number but as a rough estimate e.g. to bin a histogram. +// +void MReflector::InitMaxR() +{ + fMaxR = 0; + + TIter Next(&fMirrors); + MMirror *m = 0; + while ((m=static_cast(Next()))) + { + // Take into account the maximum incident angle 8eg 10deg) and + // the theta-angle of the mirror and the z-distance. + const Double_t r = m->GetDist()+1.5*m->GetMaxR(); + if (r > fMaxR) + fMaxR = r; + } +} + +// -------------------------------------------------------------------------- +// +// Get the pointer to the first mirror. This is a very dangerous way of +// access, but the fastest possible. because it is the most often called +// function in ExecuteReflector we have to have a very fast access. +// +const MMirror **MReflector::GetFirstPtr() const +{ + return (const MMirror**)fMirrors.GetObjectRef(0); +} + +// -------------------------------------------------------------------------- +// +// Get number of mirrors. There should be no holes in the array! +// +const UInt_t MReflector::GetNumMirrors() const +{ + return fMirrors.GetEntriesFast(); +} + +// -------------------------------------------------------------------------- +// +// Check with a rough estimate whether a photon can hit the reflector. +// +Bool_t MReflector::CanHit(const MQuaternion &p) const +{ + // p is given in the reflectory coordinate frame. This is meant as a + // fast check without lengthy calculations to omit all photons which + // cannot hit the reflector at all + return p.R2()ExpandPathName(fname); + + ifstream fin(fname); + if (!fin) + { + *fLog << err << "Cannot open file " << fname << ": "; + *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl; + } + +/* + Int_t idx[964]; + Int_t srt[964]; + for (int i=0; i<964; i++) + srt[i] = gRandom->Integer(964); + + TMath::Sort(964, srt, idx); + */ + + fMirrors.Delete(); + + while (1) + { + TString line; + line.ReadLine(fin); + if (!fin) + break; + + if (line.BeginsWith("#")) + { + //cout << line << endl; + continue; + } + + line=line.Strip(TString::kBoth); + + if (line.IsNull()) + continue; + + TObjArray *arr = line.Tokenize(' '); + + if (arr->GetSize()<8) + { + cout << "Skip3: " <GetName()), + atof((*arr)[1]->GetName()), + atof((*arr)[2]->GetName())); + + const TVector3 norm(atof((*arr)[3]->GetName()), + atof((*arr)[4]->GetName()), + atof((*arr)[5]->GetName())); + + const Double_t F = atof((*arr)[6]->GetName()); + + TString type = (*arr)[7]->GetName(); + type.Prepend("MMirror"); + + TString msg; + TClass *cls = MParContainer::GetClass(type); + if (!cls) + { + *fLog << err << dbginf << "ERROR - Class " << type << " not in dictionary." << endl; + return kFALSE; + } + + if (!cls->InheritsFrom(MMirror::Class())) + { + *fLog << err << dbginf << "Cannot create new instance of class " << type << ": " << endl; + *fLog << "Class doesn't inherit from MMirror." << endl; + return kFALSE; + } + + MMirror *m = (MMirror*)cls->New(); + if (!m) + { + *fLog << err << dbginf << "Cannot create new instance of class " << type << ": " << endl; + *fLog << " - Class has no default constructor." << endl; + *fLog << " - An abstract member functions of a base class is not overwritten." << endl; + return kFALSE; + } + + m->SetFocalLength(F); + m->SetPosition(pos); + m->SetNorm(norm); + + Int_t n = m->ReadM(*arr); + if (n<=0) + { + *fLog << err << dbginf << "ERROR - ReadM failed." << endl; + return kFALSE; + } + + fMirrors.Add(m); + + //maxr = TMath::Max(maxr, TMath::Hypot(pos[i].X()+24.75, pos[i].Y()+24.75)); + //maxr = TMath::Max(maxr, TMath::Hypot(pos.X()+24.75, pos.Y()+24.75)); + + delete arr; + } + + InitMaxR(); + + return kTRUE; + +// fMirrors.Sort(); +/* + for (int i=0; i<964; i++) + { + MMirror &ref = (MMirror&)*fMirrors[i]; + + TArrayD dist(964); + for (int j=0; j<964; j++) + { + const MMirror &mir = (MMirror&)*fMirrors[j]; + dist[j] = (ref-mir).Mod(); + } + + TArrayI idx(964); + TMath::Sort(964, dist.GetArray(), idx.GetArray(), kFALSE); + + for (int j=0; j<964; j++) + { + ref.fNeighbors.Add(fMirrors[idx[j]]); + } + }*/ +} + +// -------------------------------------------------------------------------- +// +// Paint the collection of mirrors +// +void MReflector::Paint(Option_t *o) +{ + fMirrors.Paint(o); + /* + TIter Next(&fMirrors); + MMirror *m = 0; + + while ((m=static_cast(Next()))) + m->Paint(o);*/ +} + Index: trunk/MagicSoft/Mars/msimreflector/MReflector.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MReflector.h (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MReflector.h (revision 9236) @@ -0,0 +1,51 @@ +#ifndef MARS_MReflector +#define MARS_MReflector + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TObjArray +#include +#endif + +class MQuaternion; +class MMirror; + +class MReflector : public MParContainer +{ +private: + //Simple Array // 5.1s + TObjArray fMirrors; // 6.1s (Pointer) + //TObjArray fMirrors; // 6.1s (GetObjectRef) + //TObjArray fMirrors; // 8.3s (Next) + //TObjArray fMirrors; // 10.1s (UncheckedAt) + //TList fMirrors; // 10.7s + //TOrdCollection fMirrors; // 23.4s + + Double_t fMaxR; + + void InitMaxR(); + +public: + MReflector(const char *name=NULL, const char *title=NULL); + + const MMirror **GetFirstPtr() const; + const UInt_t GetNumMirrors() const; + + const MMirror *GetMirror(UInt_t idx) const { return idx>=GetNumMirrors()?0:*(GetFirstPtr()+idx); } + + Bool_t ReadFile(TString fname); + + Double_t GetMaxR() const { return fMaxR; } + + virtual Bool_t CanHit(const MQuaternion &p) const; + + Int_t ExecuteReflector(MQuaternion &p, MQuaternion &u) const; + + void Paint(Option_t *o); + + ClassDef(MReflector, 0) // Parameter container storing a collection of several mirrors (reflector) +}; + +#endif Index: trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc (revision 9236) @@ -0,0 +1,547 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2008 +! +! Copyright: Software Development, 2000-2008 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSimReflector +// +////////////////////////////////////////////////////////////////////////////// +#include "MSimReflector.h" + +#include + +#include "MGeomCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MQuaternion.h" +#include "MMirror.h" +#include "MReflector.h" +#include "MReflection.h" + +#include "MCorsikaEvtHeader.h" +//#include "MCorsikaRunHeader.h" + +#include "MPhotonEvent.h" +#include "MPhotonData.h" + +#include "MPointingPos.h" + +ClassImp(MSimReflector); + +using namespace std; + +// USEFUL CORSIKA OPTIONS: +// NOCLONG + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MSimReflector::MSimReflector(const char* name, const char *title) + : fEvt(0), fMirror0(0), fMirror1(0), fMirror2(0), fMirror3(0), + fMirror4(0), /*fRunHeader(0),*/ fEvtHeader(0), fReflector(0), + fGeomCam(0), fPointing(0) +{ + fName = name ? name : "MSimReflector"; + fTitle = title ? title : "Task to calculate reflection os a mirror"; +} + +// -------------------------------------------------------------------------- +// +// Search for the necessary parameter containers. +// +Int_t MSimReflector::PreProcess(MParList *pList) +{ + fMirror0 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane0"); + if (!fMirror0) + return kFALSE; + fMirror1 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane1"); + if (!fMirror1) + return kFALSE; + fMirror2 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane2"); + if (!fMirror2) + return kFALSE; + fMirror3 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane3"); + if (!fMirror3) + return kFALSE; + fMirror4 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane4"); + if (!fMirror4) + return kFALSE; + + fReflector = (MReflector*)pList->FindObject("MReflector"); + if (!fReflector) + { + *fLog << err << "MReflector not found... aborting." << endl; + return kFALSE; + } + + fGeomCam = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam"); + if (!fGeomCam) + { + *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl; + + fGeomCam = (MGeomCam*)pList->FindCreateObj(fNameGeomCam); + if (!fGeomCam) + return kFALSE; + } + + fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent"); + if (!fEvt) + { + *fLog << err << "MPhotonEvent not found... aborting." << endl; + return kFALSE; + } + /* + fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); + if (!fRunHeader) + { + *fLog << err << "MCorsikaRunHeader not found... aborting." << endl; + return kFALSE; + } + */ + fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader"); + if (!fEvtHeader) + { + *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl; + return kFALSE; + } + + fPointing = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointing) + { + *fLog << err << "MPointingPos not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// The main point of calculating the reflection is to determine the +// coincidence point of the particle trajectory on the mirror surface. +// +// If the position and the trajectory of a particle is known it is enough +// to calculate the z-value of coincidence. x and y are then well defined. +// +// Since the problem (mirror) has a rotational symmetry we only have to care +// about the distance from the z-axis. +// +// Given: +// +// p: position vector of particle (z=0) +// u: direction vector of particle +// F: Focal distance of the mirror +// +// We define: +// +// q := (px, py ) +// v := (ux/uz, uy/uz) +// r^2 := x^2 + y^2 +// +// +// Distance from z-axis: +// --------------------- +// +// q' = q - z*v (z>0) +// +// Calculate distance r (|q|) +// +// r^2 = (px-z*ux)^2 + (py-z*uy)^2 +// r^2 = px^2+py^2 + z^2*(ux^2+uy^2) - 2*z*(px*ux+py*uy) +// r^2 = |q|^2 + z^2*|v|^2 - 2*z* q*v +// +// +// Spherical Mirror Surface: (distance of surface point from 0/0/0) +// ------------------------- +// +// Sphere: r^2 + z^2 = R^2 | Parabola: z = p*r^2 +// Mirror: r^2 + (z-R)^2 = R^2 | Mirror: z = p*r^2 +// | +// Focal length: F=R/2 | Focal length: F = 1/4p +// | +// r^2 + (z-2*F)^2 = (2*F)^2 | z = F/4*r^2 +// | +// z = -sqrt(4*F*F - r*r) + 2*F | +// z-2*F = -sqrt(4*F*F - r*r) | +// (z-2*F)^2 = 4*F*F - r*r | +// z^2-4*F*z+4*F^2 = 4*F*F - r*r (4F^2-r^2>0) | z - F/4*r^2 = 0 +// z^2-4*F*z+r^2 = 0 +// +// Find the z for which our particle has the same distance from the z-axis +// as the mirror surface. +// +// substitute r^2 +// +// +// Equation to solve: +// ------------------ +// +// z^2*(1+|v|^2) - 2*z*(2*F+q*v) + |q|^2 = 0 | z^2*|v|^2 - 2*(2/F+q*v)*z + |q|^2 = 0 +// +// z = (-b +- sqrt(b*b - 4ac))/(2*a) +// +// a = 1+|v|^2 | a = |v|^2 +// b = - 2*(2*F+q*v) | b = - 2*(2/F+q*v) +// c = |q|^2 | c = |q|^2 +// | +// +// substitute b := 2*b' +// +// z = (-2*b' +- 2*sqrt(b'*b' - ac))/(2*a) +// z = (- b' +- sqrt(b'*b' - ac))/a +// z = (-b'/a +- sqrt(b'*b' - ac))/a +// +// substitute f := b'/a +// +// z = (-f +- sqrt(f^2 - c/a) +// +// ======================================================================================= +// +// After z of the incident point has been determined the position p is +// propagated along u to the plane with z=z. Now it is checked if the +// mirror was really hit (this is implemented in HasHit). +// From the position on the surface and the mirrors curvature we can +// now calculate the normal vector at the incident point. +// This normal vector is smeared out with MMirror::PSF (basically a +// random gaussian) and then the trajectory is reflected on the +// resulting normal vector. +// +Bool_t MMirror::ExecuteReflection(MQuaternion &p, MQuaternion &u) const +{ + // If the z-componenet of the direction vector is normalized to 1 + // the calculation of the incident points becomes very simple and + // the resulting z is just the z-coordinate of the incident point. + const TVector2 v(u.XYvector()/u.Z()); + const TVector2 q(p.XYvector()); + + // Radius of curvature + const Double_t G = 2*fFocalLength; + + // Find the incident point of the vector to the mirror + // u corresponds to downwaqrd going particles, thus we use -u here + const Double_t b = G - q*v; + const Double_t a = v.Mod2(); + const Double_t c = q.Mod2(); + + // Solution for q spherical (a+1) (parabolic mirror (a) instead of (a+1)) + const Double_t f = b/(a+1); + const Double_t g = c/(a+1); + + // Solution of second order polynomial (transformed: a>0) + // (The second solution can be omitted, it is the intersection + // with the upper part of the sphere) + // const Double_t dz = a==0 ? c/(2*b) : f - TMath::Sqrt(f*f - g); + const Double_t z = f - TMath::Sqrt(f*f - g); + + // Move the photon along its trajectory to the x/y plane of the + // mirror's coordinate frame. Therefor stretch the vector + // until its z-component is the distance from the vector origin + // until the vector hits the mirror surface. + // p += z/u.Z()*u; + // p is at the mirror plane and we want to propagate back to the mirror surface + p.PropagateZ(u, z); + + // MirrorShape: Now check if the photon really hit the mirror or just missed it + if (!HasHit(p)) + return kFALSE; + + // Get normal vector for reflection by calculating the derivatives + // of a spherical mirror along x and y + const Double_t d = TMath::Sqrt(G*G - p.R2()); + + // This is a normal vector at the incident point + TVector3 n(p.X(), p.Y(), -d); + // This is the obvious solution for the normal vector + // TVector3 n(-p.X()/d, -p.Y()/d, 1)); + + if (fSigmaPSF>0) + n += SimPSF(n, fFocalLength, fSigmaPSF); + + // Changes also the sign of the z-direction of flight + // This is faster giving identical results + u *= MReflection(n); + //u *= MReflection(p.X(), p.Y(), -d); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Converts the coordinates into the coordinate frame of the mirror. +// Executes the reflection calling ExecuteReflection and converts +// the coordinates back. +// Depending on whether the mirror was hit kTRUE or kFALSE is returned. +// It the mirror was not hit the result coordinates are wrong. +// +Bool_t MMirror::ExecuteMirror(MQuaternion &p, MQuaternion &u) const +{ + // Move the mirror to the point of origin and rotate the position into + // the individual mirrors coordinate frame. + // Rotate the direction vector into the mirror's coordinate frame + p -= fPos; + p *= fTilt; + u *= fTilt; + + // Move the photon along its trajectory to the x/y plane of the + // mirror's coordinate frame. Therefor stretch the vector + // until its z-component vanishes. + //p -= p.Z()/u.Z()*u; + + // p is at the reflector plane and we want to propagate back to the mirror plane + p.PropagateZ0(u); + + // Now try to propagate the photon from the plane to the mirror + // and reflect its direction vector on the mirror. + if (!ExecuteReflection(p, u)) + return kFALSE; + + // Derotate from mirror coordinates and shift the photon back to + // reflector coordinates. + // Derotate the direction vector + u *= fTilt.Inverse(); + p *= fTilt.Inverse(); + p += fPos; + + return kTRUE; +} + +// Jeder Spiegel sollte eine Liste aller andern Spiegel in der +// reihenfolge Ihrer Entfernung enthalten. Wir starten mit der Suche +// immer beim zuletzt getroffenen Spiegel! +// +// -------------------------------------------------------------------------- +// +// Loops over all mirrors of the reflector. After doing a rough check +// whether the mirror can be hit at all the reflection is executed +// calling the ExecuteMirror function of the mirrors. +// +// If a mirror was hit its index is retuened, -1 otherwise. +// +// FIXME: Do to lopping over all mirrors for all photons this is the +// most time consuming function in teh reflector simulation. By a more +// intelligent way of finding the right mirror then just testing all +// this could be accelerated a lot. +// +Int_t MReflector::ExecuteReflector(MQuaternion &p, MQuaternion &u) const +{ + //static const TObjArray *arr = &((MMirror*)fMirrors[0])->fNeighbors; + + // This way of access is somuch faster than the program is + // a few percent slower if accessed by UncheckedAt + const MMirror **s = GetFirstPtr(); + const MMirror **e = s+GetNumMirrors(); + //const MMirror **s = (const MMirror**)fMirrors.GetObjectRef(0); + //const MMirror **e = s+fMirrors.GetEntriesFast(); + //const MMirror **s = (const MMirror**)arr->GetObjectRef(0); + //const MMirror **e = s+arr->GetEntriesFast(); + + // Loop over all mirrors + for (const MMirror **m=s; mfNeighbors; + + return m-s; + } + + return -1; +} + +// -------------------------------------------------------------------------- +// +// Converts the photons into the telscope coordinate frame using the +// pointing position from MPointingPos. +// +// Reflects all photons on all mirrors and stores the final photons on +// the focal plane. Also intermediate photons are stored for debugging. +// +Int_t MSimReflector::Process() +{ + // Get arrays from event container + TClonesArray &arr = fEvt->GetArray(); + + TClonesArray &cpy0 = fMirror0->GetArray(); + TClonesArray &cpy1 = fMirror1->GetArray(); + TClonesArray &cpy2 = fMirror2->GetArray(); + TClonesArray &cpy3 = fMirror3->GetArray(); + TClonesArray &cpy4 = fMirror4->GetArray(); + cpy0.ExpandCreateFast(arr.GetEntriesFast()); + cpy1.ExpandCreateFast(arr.GetEntriesFast()); + cpy2.ExpandCreateFast(arr.GetEntriesFast()); + cpy3.ExpandCreateFast(arr.GetEntriesFast()); + cpy4.ExpandCreateFast(arr.GetEntriesFast()); + + // Initialize mirror properties + const Double_t F = fGeomCam->GetCameraDist()*100; // Focal length [cm] + + // Local sky coordinates (direction of telescope axis) + const Double_t zd = fPointing->GetZdRad(); // x==north + const Double_t az = fPointing->GetAzRad(); + + // Rotation matrix to derotate sky + TRotation rot; // The signs are positive because we align the incident point on ground to the telescope axis + rot.RotateZ( az); // Rotate point on ground to align it with the telescope axis + rot.RotateX( zd); // tilt the point from ground to make it parallel to the mirror plane + + // Now: viewed from the backside of the mirror + // x is pointing downwards + // y is pointing left + + // Rotate around z counterclockwise + // rot.RotateZ(TMath::Pi()/2); // Rotate x to point right / y downwards / z from mirror to camera + // rot.RotateX(TMath::Pi()); // Flip -y and y (also flips Z :( ) + + // Now: viewed from the backside of the mirror + // x is pointing upwards + // y is pointing right + + const TVector3 impact(fEvtHeader->GetX(), fEvtHeader->GetY(), 0); + + // Counter for number of total and final events + UInt_t cnt[6] = { 0, 0, 0, 0, 0, 0 }; + + const Int_t num = arr.GetEntriesFast(); + for (Int_t idx=0; idx(arr.UncheckedAt(idx)); + + // w is pointing away from the direction the photon comes from + // CORSIKA-orig: x(north), y(west), z(up), t(time) + // NOW: x(east), y(north), z(up), t(time) + MQuaternion p(dat->GetPosQ()); // z=0 + MQuaternion w(dat->GetDirQ()); // z<0 + + // Shift the coordinate system to the telescope. Corsika's + // coordinate system is always w.r.t. to the particle axis + p -= impact; + + // Rotate the coordinates into the reflector's coordinate system. + // It is assumed that the z-plane is parallel to the focal plane. + // (The reflector coordinate system is defined by the telescope orientation) + p *= rot; + w *= rot; + + // ---> Simulate star-light! + // w.fVectorPart.SetXYZ(0.2/17, 0.2/17, -(1-TMath::Hypot(0.3, 0.2)/17)); + + // Now propagate the photon to the z-plane in the new coordinate system + p.PropagateZ0(w); + + // Store new position and direction in the reflector's coordinate frame + dat->SetPosition(p); + dat->SetDirection(w); + + *static_cast(cpy0.UncheckedAt(cnt[0]++)) = *dat; + + // Check if the photon has hit the camera housing and holding + if (fGeomCam->HitFrame(p, w)) + continue; + + // FIXME: Do we really need this one?? + *static_cast(cpy1.UncheckedAt(cnt[1]++)) = *dat; + + // Check if the reflector can be hit at all + if (!fReflector->CanHit(p)) + continue; + + *static_cast(cpy2.UncheckedAt(cnt[2]++)) = *dat; + + // Now execute the reflection of the photon on the mirrors' surfaces + const Int_t num = fReflector->ExecuteReflector(p, w); + if (num<0) + continue; + + // Set new position and direction (w.r.t. to the reflector's coordinate system) + // Set also the index of the mirror which was hit as tag. + dat->SetTag(num); + dat->SetPosition(p); + dat->SetDirection(w); + + *static_cast(cpy3.UncheckedAt(cnt[3]++)) = *dat; + + // Propagate the photon along its trajectory to the focal plane z=F + p.PropagateZ(w, F); + + // Store new position + dat->SetPosition(p); + + *static_cast(cpy4.UncheckedAt(cnt[4]++)) = *dat; + + // Discard all photons which definitly can not hit the detector surface + if (!fGeomCam->HitDetector(p)) + continue; + + // Copy this event to the next 'new' in the list + *static_cast(arr.UncheckedAt(cnt[5]++)) = *dat; + } + + // Now we shrink the array to a storable size (for details see + // MPhotonEvent::Shrink). + fMirror0->Shrink(cnt[0]); + fMirror1->Shrink(cnt[1]); + fMirror2->Shrink(cnt[2]); + fMirror3->Shrink(cnt[3]); + fMirror4->Shrink(cnt[4]); + fEvt->Shrink(cnt[5]); + + // Doesn't seem to be too time consuming. But we could also sort later! + // (after cones, inside the camera) + fEvt->Sort(); + + // FIXME FIXME FIXME: Set maxindex, first and last time. + + return kTRUE; +} + Index: trunk/MagicSoft/Mars/msimreflector/MSimReflector.h =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MSimReflector.h (revision 9236) +++ trunk/MagicSoft/Mars/msimreflector/MSimReflector.h (revision 9236) @@ -0,0 +1,47 @@ +#ifndef MARS_MSimReflector +#define MARS_MSimReflector + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MParList; +class MGeomCam; +class MPointingPos; +class MPhotonEvent; +class MCorsikaEvtHeader; + +class MReflector; + +class MSimReflector : public MTask +{ +private: + MPhotonEvent *fEvt; //! Event storing the photons + MPhotonEvent *fMirror0; //! Event storing the photons in the mirror plane (w/o camera shadow) + MPhotonEvent *fMirror1; //! Event storing the photons in the mirror plane (w/ camera shadow) + MPhotonEvent *fMirror2; //! Event storing the photons in the mirror plane (w/ camera shadow) + MPhotonEvent *fMirror3; //! Event storing the photons in the mirror plane (w/ camera shadow) + MPhotonEvent *fMirror4; //! Event storing the photons in the mirror plane (w/ camera shadow) + //MCorsikaRunHeader *fRunHeader; //! Header storing event information + MCorsikaEvtHeader *fEvtHeader; //! Header storing event information + + MReflector *fReflector; //! + MGeomCam *fGeomCam; //! + MPointingPos *fPointing; //! + + TString fNameGeomCam; // Name of the geometry container storing the APD gemeotry + + // MTask + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MSimReflector(const char *name=NULL, const char *title=NULL); + + // MSimReflector + void SetNameGeomCam(const char *name="MGeomCam") { fNameGeomCam = name; } + + ClassDef(MSimReflector, 0) // Task to calculate reflection on a mirror +}; + +#endif