{{lastupdated_at}} by {{lastupdated_by}}

h1. TRA file format

The @TRA@ file format is used by MEGALib to store event information. Check the MEGALib methods @src/revan/src/MRERawEvent::ParseLine@ to better understand how to interpret the format. Note the inheritance of this class:
<pre><code class="cpp">
class MRERawEvent : public MRESE, public MRotationInterface
</code></pre>

h2. TI

Parameters:
# [float] Defines the event time in seconds. Time zero corresponds to @1970-01-01T01:00:00@.

Here is the relevant code in @src/revan/src/MRERawEvent::ParseLine@:
<pre><code class="cpp">
m_EventTime.Set(Line)
bool MTime::Set(const char* Line)
m_Seconds
m_NanoSeconds
}
</code></pre>Here a GammaLib code snippet for @TI=1465776080.278382500@ that gives the event time in UTC:
<pre><code class="python">
time  = gammalib.GTime('1970-01-01T01:00:00')
time += 1465776080.278382500
print(time.utc())
2016-06-13T01:00:54
</code></pre>

h2. RX

Parameters:
# [float] @m_DetectorRotationXAxis[0]@ [cm] (default: 1.0)
# [float] @m_DetectorRotationXAxis[1]@ [cm] (default: 0.0)
# [float] @m_DetectorRotationXAxis[2]@ [cm] (default: 0.0)

h2. RZ

Parameters:
# [float] @m_DetectorRotationZAxis[0]@ [cm] (default: 0.0)
# [float] @m_DetectorRotationZAxis[1]@ [cm] (default: 0.0)
# [float] @m_DetectorRotationZAxis[2]@ [cm] (default: 1.0)

h2. GX

Galactic longitude of detector X axis.
Parameters:
# [float] @Longitude@ [deg]
# [float] @Latitude@ [deg]
<pre><code class="cpp">
SetGalacticPointingXAxis(Longitude, Latitude)
</code></pre>

h2. GZ

Galactic longitude of detector Z axis.
Parameters:
# [float] @Longitude@ [deg]
# [float] @Latitude@ [deg]
<pre><code class="cpp">
SetGalacticPointingZAxis(Longitude, Latitude)
</code></pre>

h2. HX

Horizon of detector X axis.
Parameters:
# [float] @Azimuth@ [deg]
# [float] @Altitude@ [deg]

The information is set by
<pre><code class="cpp">
m_HorizonPointingXAxis.SetMagThetaPhi(1.0, (90-Altitude)*c_Rad, Azimuth*c_Rad)
</code></pre>which uses the @MVector::SetMagThetaPhi@ method
<pre><code class="cpp">
void MVector::SetMagThetaPhi(double Magnitude, double Theta, double Phi) 
{
  // Set in spherical coordinates

  double A = fabs(Magnitude);
  m_X = A*sin(Theta)*cos(Phi);
  m_Y = A*sin(Theta)*sin(Phi);
  m_Z = A*cos(Theta);
}
</code></pre>

h2. HZ

Horizon of detector Z axis.
Parameters:
# [float] @Azimuth@ [deg]
# [float] @Altitude@ [deg]

The same method
<pre><code class="cpp">
m_HorizonPointingZAxis.SetMagThetaPhi(1.0, (90-Altitude)*c_Rad, Azimuth*c_Rad)
</code></pre>is used (see HX)

h2. CE

Energy of Compton event, sets members of the @MComptonEvent@ class.
Parameters:
# [float] @m_Eg@ [keV]
# [float] @m_dEg@ [keV]
# [float] @m_Ee@ [keV]
# [float] @m_dEe@ [keV]

h2. CD

Direction of Compton events, sets members of the @MComptonEvent@ class.
Parameters:
# [float] @m_C1[0]@ [cm]
# [float] @m_C1[1]@ [cm]
# [float] @m_C1[2]@ [dm]
# [float] @m_dC1[0]@ [cm]
# [float] @m_dC1[1]@ [cm]
# [float] @m_dC1[2]@ [dm]
# [float] @m_C2[0]@ [cm]
# [float] @m_C2[1]@ [cm]
# [float] @m_C2[2]@ [dm]
# [float] @m_dC2[0]@ [cm]
# [float] @m_dC2[1]@ [cm]
# [float] @m_dC2[2]@ [dm]
# [float] @m_Ce[0]@ [cm]
# [float] @m_Ce[1]@ [cm]
# [float] @m_Ce[2]@ [dm]
# [float] @m_dCe[0]@ [cm]
# [float] @m_dCe[1]@ [cm]
# [float] @m_dCe[2]@ [dm]

h2. TL

Track length, sets member of @MComptonEvent@ class.
Parameters:
# [float] @m_TrackLength@ [cm]

h2. TE

Track initial deposit, sets member of @MComptonEvent@ class.
Parameters:
# [float] @m_TrackInitialDeposit@ [keV]

h2. LA

Lever arm, sets member of @MComptonEvent@ class.
Parameters:
# [float] @m_LeverArm@ [cm]

h2. SQ

Sequence length, sets member of @MComptonEvent@ class.
Parameters:
# [int] @m_SequenceLength@

h2. PQ

Clustering quality factor, sets member of @MComptonEvent@ class.
Parameters:
# [float] @m_ClusteringQualityFactor@

h2. CT

Compton quality factor, sets member of @MComptonEvent@ class.
Parameters:
# [float] @m_ComptonQualityFactor1@
# [float] @m_ComptonQualityFactor2@ (optional, default: 1.0)

h2. Computation of Chi and Psi

In the detector system, the event scatter direction is defined by Chi and Psi, where Chi is an azimuth angle and Psi is a zenith angle. The computation of Chi and Psi is done in MEGAlib in @src/response/src/MResponseImagingBinnedMode::Analyze@:
<pre><code class="cpp">
  // Get the data space information
  MRotation Rotation = Compton->GetDetectorRotationMatrix();

  double Phi = Compton->Phi()*c_Deg;
  MVector Dg = -Compton->Dg();
  Dg = Rotation*Dg;
  double Chi = Dg.Phi()*c_Deg;
  while (Chi < -180) Chi += 360.0;
  while (Chi > +180) Chi -= 360.0;
  double Psi = Dg.Theta()*c_Deg;
</code></pre>where
<pre><code class="cpp">
  // Direction of the second gamma-ray:
  m_Dg = (m_C2 - m_C1).Unit();
</code></pre>is the direction of the second gamma-ray, defined in @src/misc/src/MComptonEvent::Validate@.

h2. Computation of Phi

The Compton scatter angle is computed in MEGAlib in @src/global/misc/src/MComptonEvent::CalculatePhi@
<pre><code class="cpp">
bool MComptonEvent::CalculatePhi()
{
  // Compute the compton scatter angle due to the standard equation
  // i.e neglect the movement of the electron,
  // which would lead to a Doppler-broadening
  //
  // Attention, make sure to test before you call this method,
  //            that the Compton-kinematics is correct

  double Value = 1 - c_E0 * (1/m_Eg - 1/(m_Ee + m_Eg));

  if (Value <= -1 || Value >= 1) {
    return false;
  }

  m_Phi = acos(Value);

  return true;
}
</code></pre>

h2. Computation of the Earth horizon angle

In MEGAlib, an intersection test with the Earth horizon is done in @src/mimrec/src/MEarthHorizon::IsEventFromEarthByIntersectionTest@
<pre><code class="cpp">
bool MEarthHorizon::IsEventFromEarthByIntersectionTest(MPhysicalEvent* Event, bool DumpOutput) const
{
  if (Event->GetType() == MPhysicalEvent::c_Compton) {
    MComptonEvent* C = dynamic_cast<MComptonEvent*>(Event);

    double Phi = C->Phi();

    MVector ConeAxis = -C->Dg();
    
    // Rotate the ConeAxis into the Earth system:
    ConeAxis = C->GetHorizonPointingRotationMatrix()*ConeAxis;

    
    // Distance between the Earth cone axis and the Compton cone axis:
    double AxisDist = m_PositionEarth.Angle(ConeAxis);
     
    if (fabs(AxisDist - Phi) > m_HorizonAngle) {
      return false;
    } else {
      if (DumpOutput == true) {
        mout<<"ID "<<Event->GetId()<<": Cone intersects Earth: "<<AxisDist + Phi<<" < "<<m_HorizonAngle<<endl;
      }
      return true;
    }
  } else if (Event->GetType() == MPhysicalEvent::c_Pair) {
    MPairEvent* P = dynamic_cast<MPairEvent*>(Event); 
    double AxisDist = m_PositionEarth.Angle(P->GetOrigin());
    
    if (AxisDist < m_HorizonAngle) {
      if (DumpOutput == true) {
        mout<<"ID "<<Event->GetId()<<": Origin inside Earth: "<<AxisDist<<" < "<<m_HorizonAngle<<endl;
      }
      return true;
    }
  }

  return false;
}
</code></pre>Note that
<pre><code class="cpp">
MRotation MRotationInterface::GetHorizonPointingRotationMatrix() const
{
  // Return the rotation matrix of this event

  // Verify that x and z axis are at right angle:
  if (fabs(m_HorizonPointingXAxis.Angle(m_HorizonPointingZAxis) - c_Pi/2.0)*c_Deg > 0.1) {
    cout<<"Event "<<m_Id<<": Horizon axes are not at right angle, but: "<<m_HorizonPointingXAxis.Angle(m_HorizonPointingZAxis)*c_Deg<<" deg"<<endl;
  }

  // First compute the y-Axis vector:
  MVector m_HorizonPointingYAxis = m_HorizonPointingZAxis.Cross(m_HorizonPointingXAxis);

  return MRotation(m_HorizonPointingXAxis.X(), m_HorizonPointingYAxis.X(), m_HorizonPointingZAxis.X(),
                   m_HorizonPointingXAxis.Y(), m_HorizonPointingYAxis.Y(), m_HorizonPointingZAxis.Y(),
                   m_HorizonPointingXAxis.Z(), m_HorizonPointingYAxis.Z(), m_HorizonPointingZAxis.Z());
}
</code></pre>and
<pre><code class="cpp">
double MVector::Angle(const MVector& V) const 
{
  // Calculate the angle between two vectors:
  // cos Angle = (v dot w) / (|v| x |w|)

  // Protect against division by zero:
  double Denom = Mag()*V.Mag();
  if (Denom == 0) {
    return 0.0;
  } else {
    double Value = Dot(V)/Denom;
    if (Value >  1.0) Value =  1.0;
    if (Value < -1.0) Value = -1.0;
    return acos(Value);
  }
}
</code></pre>and
<pre><code class="cpp">
m_PositionEarth = MVector(0, 0, -1);
</code></pre>which is the position of Earth (center) in detector coordinates and
<pre><code class="cpp">
m_HorizonAngle = 90*c_Rad;
</code></pre>is the (azimuth-) angle from the Earth position to Earth horizon (the values indicate the default values of the @MEarthHorizon@ constructor). The values are set using
<pre><code class="cpp">
bool MEarthHorizon::SetEarthHorizon(const MVector& PositionEarth, 
                                    const double HorizonAngle)
{
  m_PositionEarth = PositionEarth.Unit();
  m_HorizonAngle = HorizonAngle;

  return true;
}
</code></pre>which is called in @src/mimrec/src/MEventSelector.cxx@.

Alternatively, a probabilistic evaluation is done using
<pre><code class="cpp">
bool MEarthHorizon::IsEventFromEarthByProbabilityTest(MPhysicalEvent* Event, bool DumpOutput) const
{
  massert(Event != 0);

  if (Event->GetType() == MPhysicalEvent::c_Compton) {
    MComptonEvent* C = dynamic_cast<MComptonEvent*>(Event);

    // Take care of scatter angles larger than 90 deg:
    double Phi = C->Phi();
    MVector ConeAxis = C->Dg();
    // Rotate the ConeAxis into the Earth system:
    ConeAxis = C->GetHorizonPointingRotationMatrix()*ConeAxis;

    MVector Origin = C->DiOnCone();

    // That's the trick, but I don't remember what it means...
    if (Phi > c_Pi/2.0) {
      Phi = c_Pi - Phi;
    } else {
      ConeAxis *= -1;
    }

    // Distance between the Earth cone axis and the Compton cone axis:
    double EarthConeaxisDist = m_PositionEarth.Angle(ConeAxis);

    // Now determine both angles between the cone and Earth (simpler now since phi always <= 90)
    // First the one towards the Earth -- if it is smaller than 0 just use fabs
    double AngleEarthCone1 = EarthConeaxisDist - Phi;
    if (AngleEarthCone1 < 0) AngleEarthCone1 = fabs(AngleEarthCone1); // please don't simplify
    // Then the one away from Earth -- if it is > 180 deg then use the smaller angle
    double AngleEarthCone2 = EarthConeaxisDist + Phi;
    if (AngleEarthCone2 > c_Pi) AngleEarthCone2 -= 2*(AngleEarthCone2 - c_Pi);
    
    // Case a: Cone is completely inside Earth
    if (AngleEarthCone1 < m_HorizonAngle && AngleEarthCone2 < m_HorizonAngle) {
      mdebug<<"EHC: Cone is completely inside Earth"<<endl;
      if (m_MaxProbability < 1.0) {
        if (DumpOutput == true) {
          mout<<"ID "<<Event->GetId()<<": Cone inside Earth"<<endl;
        }
        return true;
      } else {
        // If the probability is 1.0, we are OK with events from Earth
        return false;
      }
    }
    // Case b: Cone is completely outside Earth
    else if (AngleEarthCone1 > m_HorizonAngle && AngleEarthCone2 > m_HorizonAngle) {
      mdebug<<"EHC: Cone is completely outside Earth"<<endl;
      return false;
    }
    // Case c: Cone intersects horizon
    else {
      // Determine the intersection points on the Compton cone:
      mdebug<<"EHC: Cone intersects horizon"<<endl;
      
      if (sin(EarthConeaxisDist) == 0 || sin(Phi) == 0) {
        merr<<"Numerical boundary: Scattered gamma-ray flew in direction of "
            <<"the Earth axis and the Compton scatter angle is identical with "
            <<"horizon angle (horizon angle="<<m_HorizonAngle*c_Deg
            <<", Earth-coneaxis-distance="<<EarthConeaxisDist*c_Deg
            <<") - or we have a Compton backscattering (180 deg) " 
            <<"or no scattering at all (phi="<<Phi*c_Deg<<")... "
            <<"Nevertheless, I am not rejecting event "<<Event->GetId()<<show;
        return false;
      }
      // Law of cosines for the sides of spherical triangles
      double EarthConeaxisIntersectionAngle = 
        acos((cos(m_HorizonAngle)-cos(EarthConeaxisDist)*cos(Phi))/(sin(EarthConeaxisDist)*sin(Phi)));

      if (C->HasTrack() == true && m_ValidComptonResponse == true) {
        // Now we have to determine the distance to the origin on the cone:

        double ConeaxisOriginDist = ConeAxis.Angle(Origin);
        double EarthOriginDist = m_PositionEarth.Angle(Origin);

        if (sin(ConeaxisOriginDist) == 0) {
          merr<<"Numerical boundary: The photon's origin is identical with the cone axis...!"
              <<" Nevertheless, I am not rejecting this event "<<Event->GetId()<<show;
          return false;
        }

        // Law of cosines for the sides of spherical triangles
        double EarthConeaxisOriginAngle = 
          acos((cos(EarthOriginDist) - cos(EarthConeaxisDist)*cos(ConeaxisOriginDist))/
               (sin(EarthConeaxisDist)*sin(ConeaxisOriginDist)));

        // The intersections on the cone are at EarthConeaxisOriginAngle +-EarthConeaxisIntersectionAngle
        // Let's figure out the probabilities:

        double Probability = 0;

        // Case A: "Maximum" on cone (not necessarily origin) is from *within* Earth:
        if (EarthOriginDist < m_HorizonAngle) { 
          // If the intersection are in different Origin-Coneaxis-Hemispheres:
          if (EarthConeaxisOriginAngle + EarthConeaxisIntersectionAngle < c_Pi) {
            Probability = m_ComptonResponse.GetInterpolated((EarthConeaxisIntersectionAngle - EarthConeaxisOriginAngle)*c_Deg, C->Ee()) +
              m_ComptonResponse.GetInterpolated((EarthConeaxisIntersectionAngle + EarthConeaxisOriginAngle)*c_Deg, C->Ee());
            Probability = 0.5*Probability;

            massert(EarthConeaxisIntersectionAngle - EarthConeaxisOriginAngle >= 0);
            massert(EarthConeaxisIntersectionAngle - EarthConeaxisOriginAngle <= c_Pi);
            massert(EarthConeaxisIntersectionAngle + EarthConeaxisOriginAngle >= 0);
            massert(EarthConeaxisIntersectionAngle + EarthConeaxisOriginAngle <= c_Pi);

          } 
          // otherwise:
          else {
            Probability = m_ComptonResponse.GetInterpolated((2*c_Pi - EarthConeaxisIntersectionAngle - EarthConeaxisOriginAngle)*c_Deg, C->Ee()) -
              m_ComptonResponse.GetInterpolated((EarthConeaxisIntersectionAngle - EarthConeaxisOriginAngle)*c_Deg, C->Ee());
            Probability = 1-0.5*Probability;
          }
        } 
        // Case B: "Maximum" on cone (not necessarily origin) is from *outside* Earth:
        else {
          // If the intersection are in different hemispheres Origin-Coneaxis-Hemispheres:
          if (EarthConeaxisOriginAngle + EarthConeaxisIntersectionAngle > c_Pi) {
            Probability = m_ComptonResponse.GetInterpolated((EarthConeaxisOriginAngle - EarthConeaxisIntersectionAngle)*c_Deg, C->Ee()) +
              m_ComptonResponse.GetInterpolated((c_Pi + EarthConeaxisIntersectionAngle - EarthConeaxisOriginAngle)*c_Deg, C->Ee());
            Probability = 1-0.5*Probability;
          } 
          // otherwise:
          else {
            Probability = m_ComptonResponse.GetInterpolated((EarthConeaxisOriginAngle + EarthConeaxisIntersectionAngle)*c_Deg, C->Ee()) -
              m_ComptonResponse.GetInterpolated((EarthConeaxisOriginAngle - EarthConeaxisIntersectionAngle)*c_Deg, C->Ee());
            Probability = 0.5*Probability;

          } 
        }

        if (Probability > m_MaxProbability) {
          if (DumpOutput == true) {
            mout<<"ID "<<Event->GetId()<<": Probability higher max probability: "<<Probability<<" > "<<m_MaxProbability<<endl;
          }
          return true;
        }
      } else {
        // The probability is simply determined by the length of the segment within Earth
        if (EarthConeaxisIntersectionAngle/c_Pi > m_MaxProbability) {
          if (DumpOutput == true) {
            mout<<"ID "<<Event->GetId()<<": Probability higher max probability: "<<EarthConeaxisIntersectionAngle/c_Pi<<" > "<<m_MaxProbability<<endl;
          }
          return true;
        }
      }
    } 
  } else if (Event->GetType() == MPhysicalEvent::c_Pair) {
    MPairEvent* P = dynamic_cast<MPairEvent*>(Event); 
    double AxisDist = m_PositionEarth.Angle(P->GetOrigin());
    
    if (AxisDist < m_HorizonAngle) {
      if (DumpOutput == true) {
        mout<<"ID "<<Event->GetId()<<": Origin inside Earth: "<<AxisDist<<" < "<<m_HorizonAngle<<endl;
      }
      return true;
    }
  }

  return false;
}
</code></pre>

h2. MEGAlib constants

The following constants are used in MEGAlib:
<pre><code class="cpp">
const double c_Pi = 3.14159265358979311600;
const double c_TwoPi = 2*c_Pi;
const double c_Sqrt2Pi = 2.506628274631;
const double c_Rad = c_Pi / 180.0;
const double c_Deg = 180.0 / c_Pi;
const double c_SpeedOfLight = 29.9792458E+9; // cm/s
const double c_E0 = 510.999; // keV
const double c_FarAway = 1E20; // cm
const double c_LargestEnergy = 0.999*numeric_limits<float>::max();
const MVector c_NullVector(0.0, 0.0, 0.0);
</code></pre>