Index: /trunk/Mars/msimreflector/MLens.cc =================================================================== --- /trunk/Mars/msimreflector/MLens.cc (revision 19591) +++ /trunk/Mars/msimreflector/MLens.cc (revision 19592) @@ -33,18 +33,9 @@ #include -#include // atof (Ubuntu 8.10) - -#include -#include -#include - #include "MQuaternion.h" -#include "MMirror.h" #include "MLog.h" #include "MLogManip.h" -#include "MH.h" - ClassImp(MLens); @@ -60,35 +51,5 @@ fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)"; - fMirrors.SetOwner(); -} - -// -------------------------------------------------------------------------- -// -// Set the SigmaPSF of all mirrors currently stored. -// -void MLens::SetSigmaPSF(Double_t psf) -{ - fMirrors.R__FOR_EACH(MMirror, SetSigmaPSF)(psf); -} - -// -------------------------------------------------------------------------- -// -// 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 MLens::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; - } + fMaxR = 25; } @@ -100,32 +61,5 @@ Double_t MLens::GetA() const { - Double_t A = 0; - - TIter Next(&fMirrors); - MMirror *m = 0; - while ((m=static_cast(Next()))) - A += m->GetA(); - - return A; -} - -// -------------------------------------------------------------------------- -// -// 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 **MLens::GetFirstPtr() const -{ - return (const MMirror**)fMirrors.GetObjectRef(0); -} - -// -------------------------------------------------------------------------- -// -// Get number of mirrors. There should be no holes in the array! -// -const UInt_t MLens::GetNumMirrors() const -{ - return fMirrors.GetEntriesFast(); + return fMaxR*fMaxR*TMath::Pi(); } @@ -142,390 +76,198 @@ } -// 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 MLens::ExecuteOptics(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; + fTanTheta = tan(theta); + + fTanTheta2[0] = fTanTheta*fTanTheta; + fTanTheta2[1] = fTanTheta*fTanTheta * fPeakZ; + fTanTheta2[2] = fTanTheta*fTanTheta * fPeakZ*fPeakZ; } - - return -1; -} - -// -------------------------------------------------------------------------- -// -// I/O helper for ReadFile to avoid calling "delete arr" before every return -// -MMirror *MLens::EvalTokens(TObjArray &arr, Double_t defpsf) const -{ - if (arr.GetEntries()<9) +}; + +float CalcIntersection(const MQuaternion &p, const MQuaternion &u, const Groove &g, const double &fThickness) +{ + // p is the position in the plane of the lens (px, py, 0, t) + // u is the direction of the photon (ux, uy, dz, dt) + + // Assuming a cone with peak at z and opening angle theta + + // Radius at distance z+dz from the tip and at distance r from the axis of the cone + // r = (z-dz) * tan(theta) + + // Defining + // zc := z+dz + + // Surface of the cone + // (X) ( tan(theta) * sin(phi) ) + // (Y) = zc * ( tan(theta) * cos(phi) ) + // (Z) ( 1 ) + + // Line + // (X) (u.x) (p.x) + // (Y) = dz * (u.y) + (p.y) + // (Z) (u.z) ( 0 ) + + // Equality + // + // X^2+Y^2 = r^2 + // + // (dz*u.x+p.x)^2 + (dz*u.y+p.y)^2 = (z-dz)^2 * tan(theta)^2 + // + // dz^2*ux^2 + px^2 + 2*dz*ux*px + dz^2*uy^2 + py^2 + 2*dz*uy*py = (z^2*tan(theta)^2 + dz^2*tan(theta)^2 - 2*z*dz*tan(theta)^2 + // + // dz^2*(ux^2 + uy^2 - tan(theta)^2) + 2*dz*(ux*px + uy*py + z*tan(theta)^2) + (px^2+py^2 - z^2*tan(theta)^2) = 0 + + // t := tan(theta) + // a := ux^2 + uy^2 - t^2 + // b/2 := ux*px + uy*py - z *t^2 := B + // c := px^2 + py^2 - z^2*t^2 + + // dz = [ -b +- sqrt(b^2 - 4*a*c) ] / [2*a] + // dz = [ -B +- sqrt(B^2 - a*c) ] / a + + const double a = u.R2() - g.fTanTheta2[0]; + const double B = u.XYvector()*p.XYvector() + g.fTanTheta2[1]; + const double c = p.R2() - g.fTanTheta2[2]; + + const double B2 = B*B; + const double ac = a*c; + + if (B2GetName()), - atof(arr[1]->GetName()), - atof(arr[2]->GetName())); - - const TVector3 norm(atof(arr[3]->GetName()), - atof(arr[4]->GetName()), - atof(arr[5]->GetName())); - - UInt_t n = 6; - - TString six = arr[n]->GetName(); - - Char_t shape = 'S'; - if (!six.IsFloat()) - { - shape = six[0]; - n++; - } - - const Double_t F = atof(arr[n++]->GetName()); - - const TString val = arr[n++]->GetName(); - - const Double_t psf = val.IsFloat() ? val.Atof() : -1; - - n += val.IsFloat() ? 1 : 0; - - TString type = arr[n-1]->GetName(); - type.Prepend("MMirror"); - - for (UInt_t i=0; iInheritsFrom(MMirror::Class())) - { - *fLog << err << "ERROR - Cannot create new instance of class " << type << ": " << endl; - *fLog << "Class doesn't inherit from MMirror." << endl; - return 0; - } - - MMirror *m = static_cast(cls->New()); - if (!m) - { - *fLog << err << "ERROR - 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 0; - } - - m->SetPosition(pos); - m->SetNorm(norm); - m->SetShape(shape); - m->SetFocalLength(F); - m->SetSigmaPSF(psf>=0 ? psf : defpsf); - - if (m->ReadM(arr)>=0) - return m; - - *fLog << err << "ERROR - " << type << "::ReadM failed." << endl; - - delete m; + (sqrtBac-B) / a, + (sqrtBac+B) / a + }; + + // Only one dz[i] can be within the lens (due to geometrical reasons + + if (dz[0]<0 && dz[0]>fThickness) + return dz[0]; + + if (dz[1]<0 && dz[1]>fThickness) + return dz[1]; + return 0; } -// -------------------------------------------------------------------------- -// -// Read a reflector setup from a file. This needs improvemtn. -// -// The file structur is like: -// -// x y z nx ny nz F [psf] Type ... -// -// x: x-coordinate of the mirror's center -// y: y-coordinate of the mirror's center -// z: z-coordinate of the mirror's center -// nx: x-component of the normal vecor in the mirror center -// ny: y-component of the normal vecor in the mirror center -// nz: z-component of the normal vecor in the mirror center -// [shape]: S for spherical , P for parabolic -// F: Focal distance of a spherical mirror -// [psf]: This number is the psf given in the units of x,y,z and -// defined at the focal distance F. It can be used to overwrite -// the second argument given in ReadFile for individual mirrors. -// Type: A instance of a mirrot of the class Type MMirrorType is created -// (Type can be, for example, Hex for for MMirrorHex). -// ...: Additional arguments as defined in MMirrorType::ReadM -// -// -// Coordinate System: -// The coordinate system is local in the reflectors frame. -// It is defined viewing from the back side of the reflector -// towards the camera. (x "right", y "up", z from reflector to camera) -// Note, that it is left-handed! -// -Bool_t MLens::ReadFile(TString fname, Double_t defpsf) -{ - SetTitle(fname); - fMirrors.Delete(); - - gSystem->ExpandPathName(fname); - - ifstream fin(fname); - if (!fin) - { - *fLog << err << "Cannot open file " << fname << ": "; - *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl; - return kFALSE; - } - -/* - 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); - */ - - Int_t cnt = 0; - - while (1) - { - TString line; - line.ReadLine(fin); - if (!fin) - break; - - // Count lines - cnt++; - - // Skip comments - if (line.BeginsWith("#")) - continue; - - // Remove leading and trailing whitespaces - line=line.Strip(TString::kBoth); - - // Skip empty lines - if (line.IsNull()) - continue; - - // Tokenize line - TObjArray *arr = line.Tokenize(' '); - - // Evaluate tokens - MMirror *m = EvalTokens(*arr, defpsf); - - // Delete now obsolete array - delete arr; - - // Check if a new mirror could be created successfully - if (!m) - { - *fLog << err << "Error in line " << cnt << ": " << line << endl; - return kFALSE; - } - - // Add new mirror to array - fMirrors.Add(m); - } - - 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]]); - } - }*/ -} - -// ------------------------------------------------------------------------ -// -Bool_t MLens::WriteFile(TString fname) const -{ - gSystem->ExpandPathName(fname); - - ofstream fout(fname); - if (!fout) - { - *fLog << err << "Cannot open file " << fname << ": "; - *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl; - return kFALSE; - } - - TIter Next(&fMirrors); - MMirror *m = 0; - while ((m=static_cast(Next()))) - { - // x: x-coordinate of the mirror's center - // y: y-coordinate of the mirror's center - // z: z-coordinate of the mirror's center - // nx: x-component of the normal vecor in the mirror center - // ny: y-component of the normal vecor in the mirror center - // nz: z-component of the normal vecor in the mirror center - // [shape]: S for spherical , P for parabolic - // F: Focal distance of a spherical mirror - // [psf]: This number is the psf given in the units of x,y,z and - // defined at the focal distance F. It can be used to overwrite - // the second argument given in ReadFile for individual mirrors. - // Type: A instance of a mirrot of the class Type MMirrorType is created - // (Type can be, for example, Hex for for MMirrorHex). - // ...: Additional arguments as defined in MMirrorType::ReadM - - fout << setw(12) << m->X() << " "; - fout << setw(12) << m->Y() << " "; - fout << setw(12) << m->Z() << " "; - fout << setw(12) << m->Nx() << " "; - fout << setw(12) << m->Ny() << " "; - fout << setw(12) << m->Nz() << " "; - if (m->GetShape()==1) - fout << "P "; - fout << m->GetFocalLength() << " "; - // cout << m->GetSigmaPSF() << " "; - fout << m->ClassName()+7 << " "; - m->WriteM(fout); - fout << endl; - } - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Print the collection of mirrors -// -void MLens::Print(Option_t *o) const -{ - *fLog << all << GetDescriptor() << " (" << GetNumMirrors() << "): " << endl; - - fMirrors.Print(o); -} - -// -------------------------------------------------------------------------- -// -// Paint the collection of mirrors -// -void MLens::Paint(Option_t *o) -{ - if (!TString(o).Contains("same", TString::kIgnoreCase)) - MH::SetPadRange(-fMaxR*1.01, -fMaxR*1.01, fMaxR*1.01, fMaxR*1.01); - - fMirrors.Paint(o); - - TEllipse e; - e.SetFillStyle(0); - e.SetLineColor(17); - e.PaintEllipse(0, 0, fMaxR, fMaxR, 0, 360, 0); -} - -// -------------------------------------------------------------------------- -// -// SigmaPSF: -1 -// FileName: reflector.txt -// -// SigmaPSF can be used to set a default for the psf of the mirrors -// read from the file. Note, that this can be overwritten for individual -// mirrors in the file. The SigmaPSF is the sigma of a 1D-Gauss fitted -// to the radial profile of the light distribution. -// -// For details on the file structure see MLens::ReadFile -// -Int_t MLens::ReadEnv(const TEnv &env, TString prefix, Bool_t print) -{ - Bool_t rc = kFALSE; - - Double_t psf = -1; - if (IsEnvDefined(env, prefix, "SetSigmaPSF", print)) - { - rc = kTRUE; - psf = GetEnvValue(env, prefix, "SetSigmaPSF", -1); - } - - if (IsEnvDefined(env, prefix, "FileName", print)) - { - rc = kTRUE; - if (!ReadFile(GetEnvValue(env, prefix, "FileName", ""), psf)) - return kERROR; - } - - return rc; -} +double CalcRefractiveIndex(const double &lambda) +{ + // https://refractiveindex.info/?shelf=organic&book=poly(methyl_methacrylate)&page=Szczurowski + + // n^2-1=\frac{0.99654λ^2}{λ^2-0.00787}+\frac{0.18964λ^2}{λ^2-0.02191}+\frac{0.00411λ^2}{λ^2-3.85727} + + const double l2 = lambda*lambda; + + const double c0 = 0.99654/(1-0.00787/l2); + const double c1 = 0.18964/(1-0.02191/l2); + const double c2 = 0.00411/(1-3.85727/l2); + + return sqrt(1+c0+c1+c2); +} + +void ApplyRefraction(MQuaternion &u, const TVector3 &n, const double &n1, const double &n2) +{ + // u: incoming direction + // n: normal vector of surface + // n2: refractive index of new(?) medium + // n1: refractive index of old(?) medium + + const double r = n2/n1; + const double c = n*u.fVectorPart; + + const double R = 1 - r*r*(1-c*c); + + const auto v = r*c + sqrt(R<0 ? 0 : R); + + u = r*u - n*v; +} + +Int_t MLens::ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &wavelength) const +{ + const double F = 50; // [mm] focal length + const double R = 25; // [mm] radius of lens + const int N = 25; // [cnt] Nuber of grooves within radius + + const double w = R/N; // [mm] Width of a single groove + + const double n = CalcRefractiveIndex(wavelength==0 ? 450 : wavelength); + + const int nth = TMath::FloorNint(p.R()/w); // Ray will hit the nth groove + + const double r0 = nth*w; // Lower radius of nth groove + const double r1 = (nth+1)*w; // Upper radius of nth groove + const double cen = (nth+0.5)*w; // Center of nth groove + + // FIXME: Do we have to check the nth-1 and nth+1 groove as well? + + // Angle to the center of the groove + const double phi = acos(cen / (F/1.125)); + + // theta peak_z + const Groove g( TMath::Pi()/2 - phi, r0*tan(phi)); + + // Calculate the insersection between the ray and the cone + float dz = CalcIntersection(p, u, g, -3); + + // No groove was hit + if (dz>=0) + return -1; + + // Propagate to the hit point + p.PropagateZ(u, dz); + + // Check where the ray has hit + const double pr = p.R(); + + // If the ray has not hit within the right radius.. reject + if (pr<=r0 || pr>r1) + return -1; + + // Normal vector of the surface at the hit point + // FIXME: Surface roughness? + TVector3 norm; + norm.SetMagThetaPhi(1, g.fTheta+TMath::Pi()/2, p.XYvector().Phi()); + + // Apply refraction at lens entrance (change directional vector) + // FIXME: Surace roughness + ApplyRefraction(u, norm, 1, n); + + // Propagate ray until bottom of lens + const double v = u.fRealPart; // c changes in the medium + u.fRealPart = n*v; + p.PropagateZ(u, -3); + u.fRealPart = v; // Set back to c + + // New normal surface (bottom of lens) + norm.SetMagThetaPhi(1, 0, 0); + + // Apply refraction at lens exit (change directional vector) + // FIXME: Surace roughness + ApplyRefraction(u, norm, n, 1); + + // To make this consistent with a mirror system, we now change our coordinate system + // Rays from the lens to the camera are up-going (positive sign) + u.fVectorPart.SetZ(-u.Z()); + + // This is only correct if the focal distance is calculated from the inner lens surface! + p.fVectorPart.SetZ(0); + + return nth; // Keep track of groove index +} Index: /trunk/Mars/msimreflector/MLens.h =================================================================== --- /trunk/Mars/msimreflector/MLens.h (revision 19591) +++ /trunk/Mars/msimreflector/MLens.h (revision 19592) @@ -6,42 +6,15 @@ #endif -#ifndef ROOT_TObjArray -#include -#endif - class MQuaternion; -class MMirror; class MLens : public MOptics { 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(); - // Helper for I/O - MMirror *EvalTokens(TObjArray &arr, Double_t defpsf) const; - - // MParContainer - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); - public: MLens(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 defpsf=-1); - Bool_t WriteFile(TString fname) const; Double_t GetMaxR() const { return fMaxR; } @@ -50,14 +23,7 @@ virtual Bool_t CanHit(const MQuaternion &p) const; - Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u) const; + Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u, const Short_t &) const; - void SetSigmaPSF(Double_t psf); - - virtual Bool_t IsValid() const { return fMirrors.GetEntries(); } - - - // TObject - void Paint(Option_t *o); - void Print(Option_t *o) const; + Bool_t IsValid() const { return kTRUE; } ClassDef(MLens, 1) // Parameter container storing the description of a lens