ShowerTrace.cxx

Go to the documentation of this file.

#include <cassert>
#include <iostream>

#include "ShowerTrace.h"

#include "MessageService/MsgService.h"
#include "AtNuEvent/AtmosShower.h"

#include <cmath>

static const double C45 = sqrt(2.) / 2.;
static const double HalfSideL = 4. / (1. + sqrt(2.));

static const double MaxPossibleTrace = 15.3;
static const double MaxPossibleTraceZ = 14.7;

//The below values needn't be overly accurate, since the vertex and
//direction find for a shower can't be the same for the far detector
static const double SM1EdgesZ[2] = { 0.00, 14.75};
static const double SM2EdgesZ[2] = {15.96, 29.94};

CVSID("$Id: ShowerTrace.cxx,v 1.3 2009/02/13 15:03:01 blake Exp $");

void ShowerTrace::Zero() {
  Side = -1;
  for(int i=0;i<3;i++) TraceVtx[i] = 0.;
 
  Trace = MaxPossibleTrace;
  TraceZ = MaxPossibleTraceZ;

  ReSide = -1;
  for(int i=0;i<3;i++) ReTraceVtx[i] = 0.;
  ReTrace = MaxPossibleTrace;
  ReTraceZ = MaxPossibleTraceZ;

  UpSide = -1;
  for(int i=0;i<3;i++) UpTraceVtx[i] = 0.;
  UpTrace = MaxPossibleTrace;
  UpTraceZ = MaxPossibleTraceZ;
}

void ShowerTrace::Fill(const AtmosShower *shower) {
  if(!shower) return;

  const double ShwVtxX = shower->VtxX;
  const double ShwVtxY = shower->VtxY;
  const double ShwVtxZ = shower->VtxZ;
  const double ShwVtxU = shower->VtxU;
  const double ShwVtxV = shower->VtxV;
  const double ShwVtxL[4] = {ShwVtxV,ShwVtxY,ShwVtxU,ShwVtxX};
  const int ShwVtxSM = 1 + (shower->VtxPlane/249);

  //For uncontained vertex, set TraceVtx to ShwVtx and return;
  if (fabs(ShwVtxU) > 4.0 ||
      fabs(ShwVtxV) > 4.0 ||
      fabs(ShwVtxX) > 4.0 ||
      fabs(ShwVtxY) > 4.0) {
    double ShwVtx[3] = {ShwVtxX, ShwVtxY, ShwVtxZ};
    for (int i=0;i<3;i++) {
      TraceVtx[i] = ShwVtx[i];
      ReTraceVtx[i] = ShwVtx[i];
      UpTraceVtx[i] = ShwVtx[i];
    }
    Trace = 0.; TraceZ = 0.;
    ReTrace = 0.; ReTraceZ = 0.;
    UpTrace = 0.; UpTraceZ = 0.;
    
    return;
  }

  const double ShwCosX = shower->VtxDirCosX;
  const double ShwCosY = shower->VtxDirCosY;
  const double ShwCosZ = shower->VtxDirCosZ;
  const double ShwCosU = shower->VtxDirCosU;
  const double ShwCosV = shower->VtxDirCosV;
  const double ShwCosL[4] = {ShwCosV,ShwCosY,ShwCosU,ShwCosX};

  double ZProj(0.), OProj(0.), Wall(0.);

  //Find the forward projection
  int i = 0;
  //loop over each octagonal edge pair (0=V, 1=Y, 2=U, 3=X)
  for (i=0; i<4; i++) {
    if(fabs(ShwCosL[i]) < 1.e-6) continue;//DirCos can't be small
    Wall = ShwCosL[i] > 0.0 ? 4.0 : -4.0;//Pointing up or down
    int j = i<2 ? i+2 : i-2;//Find the orthogonal coordinate index

    //OProj is the coordinate in the orthogonal 'j' coordinate where the
    //vertex points back to in the 'i' coordinate, if it sits outside of
    //the halflength of an octagonal side, move along
    OProj = ShwVtxL[j] + (ShwCosL[j]/ShwCosL[i])*(Wall-ShwVtxL[i]);
    if(fabs(OProj) > HalfSideL) continue;

    //Calculated the ZProj for the case that gets to this point
    ZProj = ShwVtxZ + (ShwCosZ/ShwCosL[i])*(Wall-ShwVtxL[i]);
    break;
  }

  //Firt examine the cases where the projection would land on the edge
  //of one of the super modules
  if ((i>3) || (ZProj<SM1EdgesZ[0]) || (ZProj>SM2EdgesZ[1]) ||
      (ZProj>SM1EdgesZ[1] && ShwVtxZ < SM1EdgesZ[1]) ||
      (ZProj<SM2EdgesZ[0] && ShwVtxZ > SM2EdgesZ[0])) {
    //Low Z edge of SM1
    if(ShwCosZ<0.0 && ShwVtxSM==1) {
      TraceVtx[2] = SM1EdgesZ[0];
      Side = 8;
    }
    //High Z edge of SM1
    else if(ShwCosZ>0.0 && ShwVtxSM==1) {
      TraceVtx[2] = SM1EdgesZ[1];
      Side = 9;
    }
    //Low Z edge of SM2
    else if(ShwCosZ<0.0 && ShwVtxSM==2) {
      TraceVtx[2] = SM2EdgesZ[0];
      Side = 10;
    }
    //High Z edge of SM2
    else if(ShwCosZ>0.0 && ShwVtxSM==2) {
      TraceVtx[2] = SM2EdgesZ[1];
      Side = 11;
    }
    //WTF edge of WTF, should crash and make a core dump
    else {
      MSG("ShowerTrace", Msg::kError) << "WTF, puking on shoes" << endl;
      assert(false);
    }
    TraceVtx[0] = ShwVtxX + (ShwCosX/ShwCosZ)*(TraceVtx[2]-ShwVtxZ);
    TraceVtx[1] = ShwVtxY + (ShwCosY/ShwCosZ)*(TraceVtx[2]-ShwVtxZ);
  }
  else {
    TraceVtx[2] = ZProj;
    Side = ShwCosL[i] > 0.0 ? i : i+4;
    switch(i) {
     //TraceVertex on V edge
     case 0:
      TraceVtx[0] = C45 * (OProj - Wall);
      TraceVtx[1] = C45 * (OProj + Wall);
      break;
     //TraceVertex on Y edge
     case 1:
      TraceVtx[1] = Wall;
      TraceVtx[0] = OProj;
      break;
     //TraceVertex on U edge
     case 2:
      TraceVtx[0] = C45 * (Wall - OProj);
      TraceVtx[1] = C45 * (Wall + OProj);
      break;
     //TraceVertex on X edge
     case 3:
      TraceVtx[1] = OProj;
      TraceVtx[0] = Wall;
      break;
     //TraceVertex on WTF edge, should crash and make a core dump
     default:
      MSG("ShowerTrace", Msg::kError) << "WTF, puking on shoes" << endl;
      assert(false);
      break;
    }
  }
  Trace = sqrt(
    (TraceVtx[0]-ShwVtxX)*(TraceVtx[0]-ShwVtxX) +
    (TraceVtx[1]-ShwVtxY)*(TraceVtx[1]-ShwVtxY) +
    (TraceVtx[2]-ShwVtxZ)*(TraceVtx[2]-ShwVtxZ));
  TraceZ = fabs(TraceVtx[2] - ShwVtxZ);

  //Find the reverse projection
  //loop over each octagonal edge pair (0=V, 1=Y, 2=U, 3=X)
  for (i=0; i<4; i++) {
    if(fabs(ShwCosL[i]) < 1.e-6) continue;//DirCos can't be small
    Wall = ShwCosL[i] > 0.0 ? -4.0 : 4.0;//Pointing up or down, pick opposite side
    int j = i<2 ? i+2 : i-2;
    //No need to negate the below angles since it would cancel in the
    //2-D projection, snap.
    OProj = ShwVtxL[j] + (ShwCosL[j]/ShwCosL[i])*(Wall-ShwVtxL[i]);
    if(fabs(OProj) > HalfSideL) continue;
    ZProj = ShwVtxZ + (ShwCosZ/ShwCosL[i])*(Wall-ShwVtxL[i]);
    break;
  }

  if ((i>3) || (ZProj<SM1EdgesZ[0]) || (ZProj>SM2EdgesZ[1]) ||
      (ZProj>SM1EdgesZ[1] && ShwVtxZ < SM1EdgesZ[1]) ||
      (ZProj<SM2EdgesZ[0] && ShwVtxZ > SM2EdgesZ[0])) {
    //Low Z edge of SM1
    if(ShwCosZ>0.0 && ShwVtxSM==1) {
      ReTraceVtx[2] = SM1EdgesZ[0];
      ReSide = 8;
    }
    //High Z edge of SM1
    else if(ShwCosZ<0.0 && ShwVtxSM==1) {
      ReTraceVtx[2] = SM1EdgesZ[1];
      ReSide = 9;
    }
    //Low Z edge of SM2
    else if(ShwCosZ>0.0 && ShwVtxSM==2) {
      ReTraceVtx[2] = SM2EdgesZ[0];
      ReSide = 10;
    }
    //High Z edge of SM2
    else if(ShwCosZ<0.0 && ShwVtxSM==2) {
      ReTraceVtx[2] = SM2EdgesZ[1];
      ReSide = 11;
    }
    //WTF edge of WTF, should crash and make a core dump
    else {
      MSG("ShowerTrace", Msg::kError) << "WTF, puking on shoes" << endl;
      assert(false);
    }
    ReTraceVtx[0] = ShwVtxX + (ShwCosX/ShwCosZ)*(TraceVtx[2]-ShwVtxZ);
    ReTraceVtx[1] = ShwVtxY + (ShwCosY/ShwCosZ)*(TraceVtx[2]-ShwVtxZ);
  }
  else {
    ReTraceVtx[2] = ZProj;
    ReSide = ShwCosL[i] < 0.0 ? i : i+4;
    switch(i) {
     //TraceVertex on V edge
     case 0:
      ReTraceVtx[0] = C45 * (OProj - Wall);
      ReTraceVtx[1] = C45 * (OProj + Wall);
      break;
     //TraceVertex on Y edge
     case 1:
      ReTraceVtx[1] = Wall;
      ReTraceVtx[0] = OProj;
      break;
     //TraceVertex on X edge
     case 2:
      ReTraceVtx[0] = C45 * (Wall - OProj);
      ReTraceVtx[1] = C45 * (Wall + OProj);
      break;
     //TraceVertex on U edge
     case 3:
      ReTraceVtx[1] = OProj;
      ReTraceVtx[0] = Wall;
      break;
     //TraceVertex on WTF edge, should crash and make a core dump
     default:
      MSG("ShowerTrace", Msg::kError) << "WTF, puking on shoes" << endl;
      assert(false);
      break;
    }
  }
  ReTrace = sqrt(
    (ReTraceVtx[0]-ShwVtxX)*(ReTraceVtx[0]-ShwVtxX) +
    (ReTraceVtx[1]-ShwVtxY)*(ReTraceVtx[1]-ShwVtxY) +
    (ReTraceVtx[2]-ShwVtxZ)*(ReTraceVtx[2]-ShwVtxZ));
  ReTraceZ = fabs(ReTraceVtx[2] - ShwVtxZ);

  //Find the up projection by the which ever trace vertex has a higher Y
  //coordinate
  if (TraceVtx[1] >= ReTraceVtx[1]) {
    for(int j=0; j<3; j++) UpTraceVtx[j] = TraceVtx[j];
    UpTrace = Trace;
    UpTraceZ = TraceZ;
    UpSide = Side;
  }
  else {
    for(int j=0; j<3; j++) UpTraceVtx[j] = ReTraceVtx[j];
    UpTrace = ReTrace;
    UpTraceZ = ReTraceZ;
    UpSide = ReSide;
  }
}

---