UtilMisc.cxx

Go to the documentation of this file.

#include <cassert>
#include <string>
#include <sstream>

#include "TMath.h"
#include <cmath>

#include "MessageService/MsgService.h"

#include "AtNuEvent/AtmosEvent.h"
#include "AtNuEvent/AtmosScintHit.h"
#include "AtNuEvent/AtmosStrip.h"
#include "AtNuEvent/AtmosMC.h"


#include "AtNuUtils/UtilMisc.h"

CVSID("$Id: UtilMisc.cxx,v 1.5 2008/04/04 21:50:46 bspeak Exp $");

const double C45 = TMath::Sqrt(2.) / 2.; //Cosine of 45 degrees

const double SinSq2Th = 0.95;
const double SinSq2ThHi = 1.0;
const double SinSq2ThLo = 0.9;

const double DMSq = 0.0021;
const double DMSqHi = 0.003;
const double DMSqLo = 0.0015;

const double REarth = 6371.0;
const double RMine = 6360.0;

double UtilMisc::OscProb(const AtmosMC &MCInfo, double *sigma) {

  //To start, only concern with numu CC events
  if (MCInfo.IDact!=1 || TMath::Abs(MCInfo.IDnu)!=14) {
    if(sigma) *sigma = 0.0;
    return(0.0);
  }

  double NuEnergy = TMath::Abs(MCInfo.Enu);

  //Use cosy as the zenith angle
  double NuMomentum = TMath::Sqrt((MCInfo.PnuX * MCInfo.PnuX) + 
    (MCInfo.PnuY * MCInfo.PnuY) + (MCInfo.PnuZ * MCInfo.PnuZ));

  double NuCosY = MCInfo.PnuY / NuMomentum;

  return UtilMisc::OscProb(NuCosY,NuEnergy,sigma);
}

double UtilMisc::OscProb(double CosY, double NuE, double *sigma) {
  MSG("Util",Msg::kVerbose) << "UtilMisc::OscProb" << endl;
  double NuL = (REarth*REarth) - (RMine*RMine)*(1.0-(CosY*CosY));
  NuL = -1.0*RMine*CosY + TMath::Sqrt(NuL);

  double OscFact = TMath::Sin(1.27 * DMSq * NuL / NuE);
  OscFact = OscFact * OscFact * SinSq2Th;

  if (sigma) {
    double OscHi = TMath::Sin(1.27 * DMSqHi * NuL / NuE);
    OscHi = OscHi * OscHi;
    double OscLo = TMath::Sin(1.27 * DMSqLo * NuL / NuE);
    OscLo = OscLo * OscLo;

    double OscFactHi = SinSq2ThHi * TMath::Max(OscHi,OscLo);
    double OscFactLo = SinSq2ThLo * TMath::Min(OscHi,OscLo);

    *sigma = TMath::Abs(OscFactHi - OscFactLo) / 2.0;
    OscFact = (OscFactHi + OscFactLo) / 2.0;
  }

  return(OscFact);
}

bool UtilMisc::IdIsChargedLepton(int Id) {
  if(TMath::Abs(Id) == 11) return true;
  if(TMath::Abs(Id) == 13) return true;
  if(TMath::Abs(Id) == 15) return true;
  return false;
}

bool UtilMisc::IdIsNeutralLepton(int Id) {
  if(TMath::Abs(Id) == 12) return true;
  if(TMath::Abs(Id) == 14) return true;
  if(TMath::Abs(Id) == 16) return true;
  return false;
}

bool UtilMisc::IdIsNeutrino(int Id) {
  return UtilMisc::IdIsNeutralLepton(Id);
}

bool UtilMisc::IdIsLepton(int Id) {
  return UtilMisc::IdIsChargedLepton(Id) ||
         UtilMisc::IdIsNeutralLepton(Id);
}

/* IType defines an index type used through out where:
 * IType = 0(kNuECC) : nu_e CC
 *       = 1(kNuMuCC) : nu_mu CC
 *       = 2(kNuTauCC) : nu_tau CC
 *       = 3(kNuNC) : nu NC
 *       = 4(kCRMu) : CR muon
 *       = 5(kCRNeutron) : CR muon induced neutron
 *       = 6(kMCOther) : MC other
 *       = 7(kData) : Data
 */
short UtilMisc::IType(const AtmosEvent *event) {
  if(event->NScintHits == 0) return(kData);//No ScintHits recorded, must be data
  if(event->MCInfo.IDnu == 0) {
    if(event->MCInfo.Emu == 0.) return kCRNeutron;
    return(kCRMu);
  }
  if(event->MCInfo.IDact == 0) return(kNuNC);//nu NC
  if(TMath::Abs(event->MCInfo.IDnu)==16) return(kNuTauCC);//nu_tau CC, call it NC
  if(TMath::Abs(event->MCInfo.IDnu)==14) return(kNuMuCC);//nu_mu CC
  if(TMath::Abs(event->MCInfo.IDnu)==12) return(kNuECC);//nue CC
  return(kMCOther);

  MSG("UtilMisc",Msg::kSynopsis) << "What is this with IDnu = " << event->MCInfo.IDnu
         << " and IDact = " << event->MCInfo.IDact << endl;
  MSG("UtilMisc",Msg::kWarning) << "This event fails classifaction: " << event->Run
         << ", " << event->SubRun << ", " << event->Snarl << endl;
  return (-1);
}

void UtilMisc::ReportEvent(const AtmosEvent *event) {
  MSG("Util",Msg::kSynopsis) << event->Run << "-" << event->SubRun << "-" << event->Snarl
         << endl;
}

std::string UtilMisc::EventAsString(const AtmosEvent *event) {
  std::ostringstream oss;
  oss << event->Run << "-" << event->SubRun << "-" << event->Snarl;
  return oss.str();
}

static const double SM1EdgesZ[2] = { 0.00, 14.75};
static const double SM2EdgesZ[2] = {15.96, 29.94};
static const double HalfSideL = 4. / (1. + sqrt(2.));

short UtilMisc::DetectorWall(double *Vtx, double *Cos, double *DetVtx) {
  if(!Vtx || !Cos || !DetVtx) return(-1);
  const double VtxX = Vtx[0];
  const double VtxY = Vtx[1];
  const double VtxU = C45 * (Vtx[1] + Vtx[0]);
  const double VtxV = C45 * (Vtx[1] - Vtx[0]);
  const double VtxL[4] = {VtxV, VtxY, VtxU, VtxX};

  const double VtxZ = Vtx[2];
  const int VtxSM = 1 + (int)(VtxZ/15.);

  //For uncontained vertex, set Wall to Vtx and return -1;
  if (TMath::Abs(VtxU) > 4.0 || TMath::Abs(VtxV) > 4.0 ||
      TMath::Abs(VtxX) > 4.0 || TMath::Abs(VtxY) > 4.0 ||
      VtxZ < SM1EdgesZ[0] || VtxZ > SM2EdgesZ[1] ||
      (SM1EdgesZ[1] < VtxZ && VtxZ < SM2EdgesZ[0])) {
    for (int i=0; i<3; i++)  DetVtx[i] = Vtx[i];
    return -1;
  }

  const double CosX = Cos[0];
  const double CosY = Cos[1];
  const double CosZ = Cos[2];
  const double CosU = C45 * (Cos[1] + Cos[0]);
  const double CosV = C45 * (Cos[1] - Cos[0]);
  const double CosL[4] = {CosV,CosY,CosU,CosX};

  DetVtx[0] = 0.; DetVtx[1] = 0.; DetVtx[2] = 0.;

  double ZProj(0.), OProj(0.), Wall(0.);
  int Side(-1);
  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(CosL[i]) < 1.e-6) continue;//DirCos can't be small
    Wall = CosL[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 = VtxL[j] + (CosL[j]/CosL[i])*(Wall-VtxL[i]);
    if(fabs(OProj) > HalfSideL) continue;

    //Calculated the ZProj for the case that gets to this point
    ZProj = VtxZ + (CosZ/CosL[i])*(Wall-VtxL[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] && VtxZ < SM1EdgesZ[1]) ||
      (ZProj<SM2EdgesZ[0] && VtxZ > SM2EdgesZ[0])) {
    //Low Z edge of SM1
    if(CosZ<0.0 && VtxSM==1) {
      DetVtx[2] = SM1EdgesZ[0];
      Side = 8;
    }
    //High Z edge of SM1
    else if(CosZ>0.0 && VtxSM==1) {
      DetVtx[2] = SM1EdgesZ[1];
      Side = 9;
    }
    //Low Z edge of SM2
    else if(CosZ<0.0 && VtxSM==2) {
      DetVtx[2] = SM2EdgesZ[0];
      Side = 10;
    }
    //High Z edge of SM2
    else if(CosZ>0.0 && VtxSM==2) {
      DetVtx[2] = SM2EdgesZ[1];
      Side = 11;
    }
    //WTF edge of WTF, should crash and make a core dump
    else {
      MSG("UtilMisc", Msg::kError) << "WTF, puking on shoes" << endl;
      assert(false);
    }
    DetVtx[0] = VtxX + (CosX/CosZ)*(DetVtx[2]-VtxZ);
    DetVtx[1] = VtxY + (CosY/CosZ)*(DetVtx[2]-VtxZ);
  }
  else {
    DetVtx[2] = ZProj;
    Side = CosL[i] > 0.0 ? i : i+4;
    switch(i) {
     //Detector Vertex on V edge
     case 0:
      DetVtx[0] = C45 * (OProj - Wall);
      DetVtx[1] = C45 * (OProj + Wall);
      break;
     //Detector Vertex on Y edge
     case 1:
      DetVtx[1] = Wall;
      DetVtx[0] = OProj;
      break;
     //Detector Vertex on U edge
     case 2:
      DetVtx[0] = C45 * (Wall - OProj);
      DetVtx[1] = C45 * (Wall + OProj);
      break;
     //Detector Vertex on X edge
     case 3:
      DetVtx[1] = OProj;
      DetVtx[0] = Wall;
      break;
     //Detector Vertex on WTF edge, should crash and make a core dump
     default:
      MSG("UtilMisc", Msg::kError) << "WTF, puking on shoes" << endl;
      assert(false);
      break;
    }
  }
  return Side;
}

---