1  /* ======================================================================== *\


2  !


3  ! *


4  ! * This file is part of CheObs, the Modular Analysis and Reconstruction


5  ! * Software. It is distributed to you in the hope that it can be a useful


6  ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.


7  ! * It is distributed WITHOUT ANY WARRANTY.


8  ! *


9  ! * Permission to use, copy, modify and distribute this software and its


10  ! * documentation for any purpose is hereby granted without fee,


11  ! * provided that the above copyright notice appears in all copies and


12  ! * that both that copyright notice and this permission notice appear


13  ! * in supporting documentation. It is provided "as is" without express


14  ! * or implied warranty.


15  ! *


16  !


17  !


18  ! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uniwuerzburg.de>


19  !


20  ! Copyright: CheObs Software Development, 20002009


21  !


22  !


23  \* ======================================================================== */


24 


25  //////////////////////////////////////////////////////////////////////////////


26  //


27  // MSimReflector


28  //


29  // fDetectorMargin is a margin (in mm) which is given to the


30  // MGeomCam::HitDetector. It should define a margin around the area


31  // defined in HitDetector on the focal plane in which photons are kept.


32  // Usually this can be 0 because photons not hitting the detector are


33  // obsolete except they can later be "moved" inside the detector, e.g.


34  // if you use MSimPSF to emulate a PSF by moving photons randomly


35  // on the focal plane.


36  //


37  //////////////////////////////////////////////////////////////////////////////


38  #include "MSimReflector.h"


39 


40  #include <TMath.h>


41  #include <TRandom.h>


42 


43  #include "MGeomCam.h"


44 


45  #include "MLog.h"


46  #include "MLogManip.h"


47 


48  #include "MParList.h"


49 


50  #include "MQuaternion.h"


51  #include "MMirror.h"


52  #include "MReflector.h"


53  #include "MReflection.h"


54 


55  #include "MCorsikaEvtHeader.h"


56  //#include "MCorsikaRunHeader.h"


57 


58  #include "MPhotonEvent.h"


59  #include "MPhotonData.h"


60 


61  #include "MPointingPos.h"


62 


63  ClassImp(MSimReflector);


64 


65  using namespace std;


66 


67  // USEFUL CORSIKA OPTIONS:


68  // NOCLONG


69 


70  // 


71  //


72  // Default Constructor.


73  //


74  MSimReflector::MSimReflector(const char* name, const char *title)


75  : fEvt(0), fMirror0(0), fMirror1(0), fMirror2(0), fMirror3(0),


76  fMirror4(0), /*fRunHeader(0),*/ fEvtHeader(0), fReflector(0),


77  fGeomCam(0), fPointing(0), fNameReflector("MReflector"),


78  fDetectorMargin(0)


79  {


80  fName = name ? name : "MSimReflector";


81  fTitle = title ? title : "Task to calculate reflection os a mirror";


82  }


83 


84  // 


85  //


86  // Search for the necessary parameter containers.


87  //


88  Int_t MSimReflector::PreProcess(MParList *pList)


89  {


90  fMirror0 = (MPhotonEvent*)pList>FindCreateObj("MPhotonEvent", "MirrorPlane0");


91  if (!fMirror0)


92  return kFALSE;


93  fMirror1 = (MPhotonEvent*)pList>FindCreateObj("MPhotonEvent", "MirrorPlane1");


94  if (!fMirror1)


95  return kFALSE;


96  fMirror2 = (MPhotonEvent*)pList>FindCreateObj("MPhotonEvent", "MirrorPlane2");


97  if (!fMirror2)


98  return kFALSE;


99  fMirror3 = (MPhotonEvent*)pList>FindCreateObj("MPhotonEvent", "MirrorPlane3");


100  if (!fMirror3)


101  return kFALSE;


102  fMirror4 = (MPhotonEvent*)pList>FindCreateObj("MPhotonEvent", "MirrorPlane4");


103  if (!fMirror4)


104  return kFALSE;


105 


106  fReflector = (MReflector*)pList>FindObject(fNameReflector, "MReflector");


107  if (!fReflector)


108  {


109  *fLog << err << fNameReflector << " [MReflector] not found..." << endl;


110  return kFALSE;


111  }


112 


113  if (fReflector>GetNumMirrors()==0)


114  {


115  *fLog << err << "ERROR  Reflector '" << fNameReflector << "' doesn't contain a single mirror." << endl;


116  return kFALSE;


117  }


118 


119  fGeomCam = (MGeomCam*)pList>FindObject(fNameGeomCam, "MGeomCam");


120  if (!fGeomCam)


121  {


122  if (!fNameGeomCam.IsNull())


123  *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;


124 


125  fGeomCam = (MGeomCam*)pList>FindObject("MGeomCam");


126  if (!fGeomCam)


127  {


128  *fLog << err << "MGeomCam not found... aborting." << endl;


129  return kFALSE;


130  }


131  }


132 


133  fEvt = (MPhotonEvent*)pList>FindObject("MPhotonEvent");


134  if (!fEvt)


135  {


136  *fLog << err << "MPhotonEvent not found... aborting." << endl;


137  return kFALSE;


138  }


139  /*


140  fRunHeader = (MCorsikaRunHeader*)pList>FindObject("MCorsikaRunHeader");


141  if (!fRunHeader)


142  {


143  *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;


144  return kFALSE;


145  }


146  */


147  fEvtHeader = (MCorsikaEvtHeader*)pList>FindObject("MCorsikaEvtHeader");


148  if (!fEvtHeader)


149  {


150  *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl;


151  return kFALSE;


152  }


153 


154  fPointing = (MPointingPos*)pList>FindObject(/*"PointingCorsika",*/ "MPointingPos");


155  if (!fPointing)


156  {


157  *fLog << err << "MPointingPos not found... aborting." << endl;


158  return kFALSE;


159  }


160 


161  return kTRUE;


162  }


163 


164  // 


165  //


166  // The main point of calculating the reflection is to determine the


167  // coincidence point of the particle trajectory on the mirror surface.


168  //


169  // If the position and the trajectory of a particle is known it is enough


170  // to calculate the zvalue of coincidence. x and y are then well defined.


171  //


172  // Since the problem (mirror) has a rotational symmetry we only have to care


173  // about the distance from the zaxis.


174  //


175  // Given:


176  //


177  // p: position vector of particle (z=0)


178  // u: direction vector of particle


179  // F: Focal distance of the mirror


180  //


181  // We define:


182  //


183  // q := (px, py )


184  // v := (ux/uz, uy/uz)


185  // r^2 := x^2 + y^2


186  //


187  //


188  // Distance from zaxis:


189  // 


190  //


191  // q' = q  z*v (z>0)


192  //


193  // Calculate distance r (q)


194  //


195  // r^2 = (pxz*ux)^2 + (pyz*uy)^2


196  // r^2 = px^2+py^2 + z^2*(ux^2+uy^2)  2*z*(px*ux+py*uy)


197  // r^2 = q^2 + z^2*v^2  2*z* q*v


198  //


199  //


200  // Spherical Mirror Surface: (distance of surface point from 0/0/0)


201  // 


202  //


203  // Sphere: r^2 + z^2 = R^2  Parabola: z = p*r^2


204  // Mirror: r^2 + (zR)^2 = R^2  Mirror: z = p*r^2


205  // 


206  // Focal length: F=R/2  Focal length: F = 1/4p


207  // 


208  // r^2 + (z2*F)^2 = (2*F)^2  z = F/4*r^2


209  // 


210  // z = sqrt(4*F*F  r*r) + 2*F 


211  // z2*F = sqrt(4*F*F  r*r) 


212  // (z2*F)^2 = 4*F*F  r*r 


213  // z^24*F*z+4*F^2 = 4*F*F  r*r (4F^2r^2>0)  z  F/4*r^2 = 0


214  // z^24*F*z+r^2 = 0


215  //


216  // Find the z for which our particle has the same distance from the zaxis


217  // as the mirror surface.


218  //


219  // substitute r^2


220  //


221  //


222  // Equation to solve:


223  // 


224  //


225  // 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


226  //


227  // z = (b + sqrt(b*b  4ac))/(2*a)


228  //


229  // a = 1+v^2  a = v^2


230  // b =  2*(2*F+q*v)  b =  2*(2/F+q*v)


231  // c = q^2  c = q^2


232  // 


233  //


234  // substitute b := 2*b'


235  //


236  // z = (2*b' + 2*sqrt(b'*b'  ac))/(2*a)


237  // z = ( b' + sqrt(b'*b'  ac))/a


238  // z = (b'/a + sqrt(b'*b'  ac))/a


239  //


240  // substitute f := b'/a


241  //


242  // z = (f + sqrt(f^2  c/a)


243  //


244  // =======================================================================================


245  //


246  // After z of the incident point has been determined the position p is


247  // propagated along u to the plane with z=z. Now it is checked if the


248  // mirror was really hit (this is implemented in HasHit).


249  // From the position on the surface and the mirrors curvature we can


250  // now calculate the normal vector at the incident point.


251  // This normal vector is smeared out with MMirror::PSF (basically a


252  // random gaussian) and then the trajectory is reflected on the


253  // resulting normal vector.


254  //


255  Bool_t MMirror::ExecuteReflection(MQuaternion &p, MQuaternion &u) const


256  {


257  // If the zcomponenet of the direction vector is normalized to 1


258  // the calculation of the incident points becomes very simple and


259  // the resulting z is just the zcoordinate of the incident point.


260  const TVector2 v(u.XYvector()/u.Z());


261  const TVector2 q(p.XYvector());


262 


263  // Radius of curvature


264  const Double_t G = 2*fFocalLength;


265 


266  // Find the incident point of the vector to the mirror


267  // u corresponds to downwaqrd going particles, thus we use u here


268  const Double_t b = G  q*v;


269  const Double_t a = v.Mod2();


270  const Double_t c = q.Mod2();


271 


272  // Solution for q spherical (a+1) (parabolic mirror (a) instead of (a+1))


273  const Double_t f = b/(a+1);


274  const Double_t g = c/(a+1);


275 


276  // Solution of second order polynomial (transformed: a>0)


277  // (The second solution can be omitted, it is the intersection


278  // with the upper part of the sphere)


279  // const Double_t dz = a==0 ? c/(2*b) : f  TMath::Sqrt(f*f  g);


280  const Double_t z = f  TMath::Sqrt(f*f  g);


281 


282  // Move the photon along its trajectory to the x/y plane of the


283  // mirror's coordinate frame. Therefor stretch the vector


284  // until its zcomponent is the distance from the vector origin


285  // until the vector hits the mirror surface.


286  // p += z/u.Z()*u;


287  // p is at the mirror plane and we want to propagate back to the mirror surface


288  p.PropagateZ(u, z);


289 


290  // MirrorShape: Now check if the photon really hit the mirror or just missed it


291  if (!HasHit(p))


292  return kFALSE;


293 


294  // Get normal vector for reflection by calculating the derivatives


295  // of a spherical mirror along x and y


296  const Double_t d = TMath::Sqrt(G*G  p.R2());


297 


298  // This is a normal vector at the incident point


299  TVector3 n(p.X(), p.Y(), d);


300  // This is the obvious solution for the normal vector


301  // TVector3 n(p.X()/d, p.Y()/d, 1));


302 


303  if (fSigmaPSF>0)


304  n += SimPSF(n);


305 


306  // Changes also the sign of the zdirection of flight


307  // This is faster giving identical results


308  u *= MReflection(n);


309  //u *= MReflection(p.X(), p.Y(), d);


310 


311  return kTRUE;


312  }


313 


314  // 


315  //


316  // Converts the coordinates into the coordinate frame of the mirror.


317  // Executes the reflection calling ExecuteReflection and converts


318  // the coordinates back.


319  // Depending on whether the mirror was hit kTRUE or kFALSE is returned.


320  // It the mirror was not hit the result coordinates are wrong.


321  //


322  Bool_t MMirror::ExecuteMirror(MQuaternion &p, MQuaternion &u) const


323  {


324  // Move the mirror to the point of origin and rotate the position into


325  // the individual mirrors coordinate frame.


326  // Rotate the direction vector into the mirror's coordinate frame


327  p = fPos;


328  p *= fTilt;


329  u *= fTilt;


330 


331  // Move the photon along its trajectory to the x/y plane of the


332  // mirror's coordinate frame. Therefor stretch the vector


333  // until its zcomponent vanishes.


334  //p = p.Z()/u.Z()*u;


335 


336  // p is at the reflector plane and we want to propagate back to the mirror plane


337  p.PropagateZ0(u);


338 


339  // Now try to propagate the photon from the plane to the mirror


340  // and reflect its direction vector on the mirror.


341  if (!ExecuteReflection(p, u))


342  return kFALSE;


343 


344  // Derotate from mirror coordinates and shift the photon back to


345  // reflector coordinates.


346  // Derotate the direction vector


347  u *= fTilt.Inverse();


348  p *= fTilt.Inverse();


349  p += fPos;


350 


351  return kTRUE;


352  }


353 


354  // Jeder Spiegel sollte eine Liste aller andern Spiegel in der


355  // reihenfolge Ihrer Entfernung enthalten. Wir starten mit der Suche


356  // immer beim zuletzt getroffenen Spiegel!


357  //


358  // 


359  //


360  // Loops over all mirrors of the reflector. After doing a rough check


361  // whether the mirror can be hit at all the reflection is executed


362  // calling the ExecuteMirror function of the mirrors.


363  //


364  // If a mirror was hit its index is retuened, 1 otherwise.


365  //


366  // FIXME: Do to lopping over all mirrors for all photons this is the


367  // most time consuming function in teh reflector simulation. By a more


368  // intelligent way of finding the right mirror then just testing all


369  // this could be accelerated a lot.


370  //


371  Int_t MReflector::ExecuteReflector(MQuaternion &p, MQuaternion &u) const


372  {


373  //static const TObjArray *arr = &((MMirror*)fMirrors[0])>fNeighbors;


374 


375  // This way of access is somuch faster than the program is


376  // a few percent slower if accessed by UncheckedAt


377  const MMirror **s = GetFirstPtr();


378  const MMirror **e = s+GetNumMirrors();


379  //const MMirror **s = (const MMirror**)fMirrors.GetObjectRef(0);


380  //const MMirror **e = s+fMirrors.GetEntriesFast();


381  //const MMirror **s = (const MMirror**)arr>GetObjectRef(0);


382  //const MMirror **e = s+arr>GetEntriesFast();


383 


384  // Loop over all mirrors


385  for (const MMirror **m=s; m<e; m++)


386  {


387  const MMirror &mirror = **m;


388 


389  // FIXME: Can we speed up using lookup tables or


390  // indexed tables?


391 


392  // MirrorShape: Check if this mirror can be hit at all


393  // This is to avoid time consuming calculation if there is no


394  // chance of a coincidence.


395  // FIXME: Inmprove search algorithm (2D Binary search?)


396  if (!mirror.CanHit(p))


397  continue;


398 


399  // Make a local copy of position and direction which can be


400  // changed by ExecuteMirror.


401  MQuaternion q(p);


402  MQuaternion v(u);


403 


404  // Check if this mirror is hit, and if it is hit return


405  // the reflected position and direction vector.


406  // If the mirror is missed we go on with the next mirror.


407  if (!mirror.ExecuteMirror(q, v))


408  continue;


409 


410  // We hit a mirror. Restore the local copy of position and


411  // direction back into p und u.


412  p = q;


413  u = v;


414 


415  //arr = &mirror>fNeighbors;


416 


417  return ms;


418  }


419 


420  return 1;


421  }


422 


423  // 


424  //


425  // Converts the photons into the telscope coordinate frame using the


426  // pointing position from MPointingPos.


427  //


428  // Reflects all photons on all mirrors and stores the final photons on


429  // the focal plane. Also intermediate photons are stored for debugging.


430  //


431  Int_t MSimReflector::Process()


432  {


433  // Get arrays from event container


434  TClonesArray &arr = fEvt>GetArray();


435 


436  // Because we knwo in advance what the maximum storage space could


437  // be we allocated it in advance (or shrink it if it was extremely


438  // huge before)


439  // Note, that the drawback is that an extremly large event


440  // will take about five times its storage space


441  // for a moment even if a lot from it is unused.


442  // It will be freed in the next step.


443  fMirror0>Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData


444  fMirror2>Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData


445  fMirror3>Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData


446  fMirror4>Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData


447 


448  // Initialize mirror properties


449  const Double_t F = fGeomCam>GetCameraDist()*100; // Focal length [cm]


450 


451  // Local sky coordinates (direction of telescope axis)


452  const Double_t zd = fPointing>GetZdRad(); // x==north


453  const Double_t az = fPointing>GetAzRad();


454 


455  // Rotation matrix to derotate sky


456  // For the new coordinate system see the Wiki


457  TRotation rot; // The signs are positive because we align the incident point on ground to the telescope axis


458  rot.RotateZ( az); // Rotate point on ground to align it with the telescope axis


459  rot.RotateX(zd); // tilt the point from ground to make it parallel to the mirror plane


460 


461  // Now get the impact point from Corsikas output


462  const TVector3 impact(fEvtHeader>GetX(), fEvtHeader>GetY(), 0);


463 


464  // Counter for number of total and final events


465  UInt_t cnt[6] = { 0, 0, 0, 0, 0, 0 };


466 


467  const Int_t num = arr.GetEntriesFast();


468  for (Int_t idx=0; idx<num; idx++)


469  {


470  MPhotonData *dat = static_cast<MPhotonData*>(arr.UncheckedAt(idx));


471 


472  // w is pointing away from the direction the photon comes from


473  // CORSIKAorig: x(north), y(west), z(up), t(time)


474  // NOW: x(east), y(north), z(up), t(time)


475  MQuaternion p(dat>GetPosQ()); // z=0


476  MQuaternion w(dat>GetDirQ()); // z<0


477 


478  // Shift the coordinate system to the telescope. Corsika's


479  // coordinate system is always w.r.t. to the particle axis


480  p = impact;


481 


482  // Rotate the coordinates into the reflector's coordinate system.


483  // It is assumed that the zplane is parallel to the focal plane.


484  // (The reflector coordinate system is defined by the telescope orientation)


485  p *= rot;


486  w *= rot;


487 


488  // > Simulate starlight!


489  // w.fVectorPart.SetXYZ(0.2/17, 0.2/17, (1TMath::Hypot(0.3, 0.2)/17));


490 


491  // Now propagate the photon to the zplane in the new coordinate system


492  p.PropagateZ0(w);


493 


494  // Store new position and direction in the reflector's coordinate frame


495  dat>SetPosition(p);


496  dat>SetDirection(w);


497 


498  (*fMirror0)[cnt[0]++] = *dat;


499  //*static_cast<MPhotonData*>(cpy0.UncheckedAt(cnt[0]++)) = *dat;


500 


501  // Check if the photon has hit the camera housing and holding


502  if (fGeomCam>HitFrame(p, w))


503  continue;


504 


505  // FIXME: Do we really need this one??


506  //(*fMirror1)[cnt[1]++] = *dat;


507  //*static_cast<MPhotonData*>(cpy1.UncheckedAt(cnt[1]++)) = *dat;


508 


509  // Check if the reflector can be hit at all


510  if (!fReflector>CanHit(p))


511  continue;


512 


513  (*fMirror2)[cnt[2]++] = *dat;


514  //*static_cast<MPhotonData*>(cpy2.UncheckedAt(cnt[2]++)) = *dat;


515 


516  // Now execute the reflection of the photon on the mirrors' surfaces


517  const Int_t num = fReflector>ExecuteReflector(p, w);


518  if (num<0)


519  continue;


520 


521  // Set new position and direction (w.r.t. to the reflector's coordinate system)


522  // Set also the index of the mirror which was hit as tag.


523  dat>SetTag(num);


524  dat>SetPosition(p);


525  dat>SetDirection(w);


526 


527  (*fMirror3)[cnt[3]++] = *dat;


528  //*static_cast<MPhotonData*>(cpy3.UncheckedAt(cnt[3]++)) = *dat;


529 


530  // Propagate the photon along its trajectory to the focal plane z=F


531  p.PropagateZ(w, F);


532 


533  // Store new position


534  dat>SetPosition(p);


535 


536  (*fMirror4)[cnt[4]++] = *dat;


537  //*static_cast<MPhotonData*>(cpy4.UncheckedAt(cnt[4]++)) = *dat;


538 


539  // FIXME: It make make sense to move this out of this class


540  // It is detector specific not reflector specific


541  // Discard all photons which definitly can not hit the detector surface


542  if (!fGeomCam>HitDetector(p, fDetectorMargin))


543  continue;


544 


545  // Copy this event to the next 'new' in the list


546  *static_cast<MPhotonData*>(arr.UncheckedAt(cnt[5]++)) = *dat;


547  }


548 


549  // Now we shrink the array to a storable size (for details see


550  // MPhotonEvent::Shrink).


551  fMirror0>Shrink(cnt[0]);


552  //fMirror1>Shrink(cnt[1]);


553  fMirror2>Shrink(cnt[2]);


554  fMirror3>Shrink(cnt[3]);


555  fMirror4>Shrink(cnt[4]);


556  fEvt>Shrink(cnt[5]);


557 


558  // Doesn't seem to be too time consuming. But we could also sort later!


559  // (after cones, inside the camera)


560  fEvt>Sort(kTRUE);


561 


562  // FIXME FIXME FIXME: Set maxindex, first and last time.


563  // SetMaxIndex(fReflector>GetNumMirrors()1)


564  // if (fEvt>GetNumPhotons())


565  // {


566  // SetTime(fEvt>GetFirst()>GetTime(), fEvt>GetLast()>GetTime());


567  // }


568 


569  return kTRUE;


570  }


571 


572  // 


573  //


574  // DetectorMargin: 0


575  //


576  Int_t MSimReflector::ReadEnv(const TEnv &env, TString prefix, Bool_t print)


577  {


578  Bool_t rc = kFALSE;


579  if (IsEnvDefined(env, prefix, "DetectorMargin", print))


580  {


581  rc = kTRUE;


582  fDetectorMargin = GetEnvValue(env, prefix, "DetectorMargin", 0);


583  }


584 


585  return rc;


586  }

