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, 20002010


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. To switch off this check set detector margin to 1.


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 = r^2/4F


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  r^2/4F = 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*z*(2*F+q*v) + 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*b' with b' = 2*F+q*v  b =  2*b' with b' = 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 downward 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 A = fShape ? a : a+1;


274 


275  const Double_t f = b/A;


276  const Double_t g = c/A;


277 


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


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


280  // with the upper part of the sphere)


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


282 


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


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


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


286  // until the vector hits the mirror surface.


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


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


289  p.PropagateZ(u, z);


290 


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


292  if (!HasHit(p))


293  return kFALSE;


294 


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


296  // of a the mirror's surface along x and y


297  const Double_t d = fShape ? G : TMath::Sqrt(G*G  p.R2());


298 


299  // The solution for the normal vector is


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


301  // Since the normal vector doesn't need to be of normal


302  // length we can avoid an obsolete division


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


304 


305  if (fSigmaPSF>0)


306  n += SimPSF(n);


307 


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


309  // This is faster giving identical results


310  u *= MReflection(n);


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


312 


313  return kTRUE;


314  }


315 


316  // 


317  //


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


319  // Executes the reflection calling ExecuteReflection and converts


320  // the coordinates back.


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


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


323  //


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


325  {


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


327  // the individual mirrors coordinate frame.


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


329  p = fPos;


330  p *= fTilt;


331  u *= fTilt;


332 


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


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


335  // until its zcomponent vanishes.


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


337 


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


339  p.PropagateZ0(u);


340 


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


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


343  if (!ExecuteReflection(p, u))


344  return kFALSE;


345 


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


347  // reflector coordinates.


348  // Derotate the direction vector


349  u *= fTilt.Inverse();


350  p *= fTilt.Inverse();


351  p += fPos;


352 


353  return kTRUE;


354  }


355 


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


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


358  // immer beim zuletzt getroffenen Spiegel!


359  //


360  // 


361  //


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


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


364  // calling the ExecuteMirror function of the mirrors.


365  //


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


367  //


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


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


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


371  // this could be accelerated a lot.


372  //


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


374  {


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


376 


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


378  // a few percent slower if accessed by UncheckedAt


379  const MMirror **s = GetFirstPtr();


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


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


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


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


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


385 


386  // Loop over all mirrors


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


388  {


389  const MMirror &mirror = **m;


390 


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


392  // indexed tables?


393 


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


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


396  // chance of a coincidence.


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


398  if (!mirror.CanHit(p))


399  continue;


400 


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


402  // changed by ExecuteMirror.


403  MQuaternion q(p);


404  MQuaternion v(u);


405 


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


407  // the reflected position and direction vector.


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


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


410  continue;


411 


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


413  // direction back into p und u.


414  p = q;


415  u = v;


416 


417  //arr = &mirror>fNeighbors;


418 


419  return ms;


420  }


421 


422  return 1;


423  }


424 


425  // 


426  //


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


428  // pointing position from MPointingPos.


429  //


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


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


432  //


433  Int_t MSimReflector::Process()


434  {


435  // Get arrays from event container


436  TClonesArray &arr = fEvt>GetArray();


437 


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


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


440  // huge before)


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


442  // will take about five times its storage space


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


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


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


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


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


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


449 


450  // Initialize mirror properties


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


452 


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


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


455  const Double_t az = fPointing>GetAzRad();


456 


457  // Rotation matrix to derotate sky


458  // For the new coordinate system see the Wiki


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


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


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


462 


463  // Now get the impact point from Corsikas output


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


465 


466  // Counter for number of total and final events


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


468 


469  const Int_t num = arr.GetEntriesFast();


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


471  {


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


473 


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


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


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


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


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


479 


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


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


482  p = impact;


483 


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


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


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


487  p *= rot;


488  w *= rot;


489 


490  // > Simulate starlight!


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


492 


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


494  p.PropagateZ0(w);


495 


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


497  dat>SetPosition(p);


498  dat>SetDirection(w);


499 


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


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


502 


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


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


505  continue;


506 


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


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


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


510 


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


512  if (!fReflector>CanHit(p))


513  continue;


514 


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


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


517 


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


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


520  if (num<0)


521  continue;


522 


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


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


525  dat>SetTag(num);


526  dat>SetPosition(p);


527  dat>SetDirection(w);


528 


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


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


531 


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


533  p.PropagateZ(w, F);


534 


535  // Store new position


536  dat>SetPosition(p);


537 


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


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


540 


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


542  // It is detector specific not reflector specific


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


544  if (fDetectorMargin>=0 && !fGeomCam>HitDetector(p, fDetectorMargin))


545  continue;


546 


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


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


549  }


550 


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


552  // MPhotonEvent::Shrink).


553  fMirror0>Shrink(cnt[0]);


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


555  fMirror2>Shrink(cnt[2]);


556  fMirror3>Shrink(cnt[3]);


557  fMirror4>Shrink(cnt[4]);


558  fEvt>Shrink(cnt[5]);


559 


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


561  // (after cones, inside the camera)


562  fEvt>Sort(kTRUE);


563 


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


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


566  // if (fEvt>GetNumPhotons())


567  // {


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


569  // }


570 


571  return kTRUE;


572  }


573 


574  // 


575  //


576  // DetectorMargin: 0


577  //


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


579  {


580  Bool_t rc = kFALSE;


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


582  {


583  rc = kTRUE;


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


585  }


586 


587  return rc;


588  }

