VetoTrack.cxx

Go to the documentation of this file.

#include <cassert>
#include <math.h>
#include <vector>
#include <iostream>

#include "TString.h"
#include "TFile.h"
#include "TMath.h"

#include "AtNuEvent/AtmosEvent.h"
#include "AtNuEvent/AtmosTrack.h"
#include "AtNuEvent/AtmosShieldPlank.h"
#include "MessageService/MsgService.h"

#include "AtNuUtils/UtilMisc.h"
#include "AtNuUtils/VetoTrack.h"

CVSID("$Id: VetoTrack.cxx,v 1.4 2007/01/15 19:52:00 rhatcher Exp $");

const double WlsFibreN = 1.43;
const double ClearFibreN = 1.6;
const double TimeWin = 100e-9;

PlankToTrk::PlankToTrk() {
  Plane = 0;
  Section = 0;
  PlankPE = 0.;
  DZ = 0.;
  DSteel = 0.;
}

PlankToTrk::PlankToTrk(const AtmosShieldPlank *plank, float dz,
                       float dsteel) {
  Plane = plank->Plane;
  Section = plank->Section;
  PlankPE = plank->QPE[0] + plank->QPE[1];
  DZ = dz;
  DSteel = dsteel;
}

VetoTrack::~VetoTrack() {
  //if(InTimePlanes) delete [] InTimePlanes; InTimePlanes = 0;
  //if(SteelToPlanks) delete [] SteelToPlanks; SteelToPlanks = 0;
  //if(ZToPlanks) delete [] ZToPlanks; ZToPlanks = 0;
  //if(PEPlanks) delete [] PEPlanks; PEPlanks = 0;
  delete InTimePlanks;
}

VetoTrack::VetoTrack() {
  InTimePlanks = new TClonesArray("PlankToTrk");
  Zero();
}

VetoTrack::VetoTrack(const AtmosEvent *event) {
  InTimePlanks = new TClonesArray("PlankToTrk");
  Zero();
  Fill(event);
}

VetoTrack::VetoTrack(const AtmosTrack *track,
                     const TClonesArray *ShieldPlankList) {
  InTimePlanks = new TClonesArray("PlankToTrk");
  Zero();
  Fill(track, ShieldPlankList);
}

void VetoTrack::Zero() {
  memset(SectionQInTime, 0, 4*sizeof(double));
  memset(SectionQOutTime, 0, 4*sizeof(double));

  VetoQCorrSS = 0;

  //memset(VetoQPlank50cmZ, 0, 10*sizeof(double));
  //memset(VetoQPlank10cmSteel, 0, 10*sizeof(double));
  for (int i=0; i<1+InTimePlanks->GetLast(); i++){
    InTimePlanks->At(i)->Clear();
  }
  InTimePlanks->Delete();
  NInTimePlanks = 0;
}

void VetoTrack::Fill(const AtmosEvent *event) {
  //This assumes a single track over which to find the veto correlation
  const AtmosTrack* track =
    dynamic_cast<const AtmosTrack*>(event->TrackList->At(0));
  if(!track) return;

  if(event->NShieldPlanks == 0) return;
  if(!(event->ShieldPlankList)) return;

  Fill(track, event->ShieldPlankList);
}

void VetoTrack::Fill(const AtmosTrack *track,
                     const TClonesArray *ShieldPlankList) {
  assert(track);
  assert(ShieldPlankList);

  //project track vertex back to the shield
  const double VtxX = track->VtxX;
  const double VtxY = track->VtxY;
  const double VtxZ = track->VtxZ;
  const double VtxTime = track->VtxTime;

  //Which sections of the shield do we care about?
  int vtx_section(0);
  if (VtxZ>0.0 && VtxZ<15) {//sm1
    if(VtxZ<8.5) vtx_section|=1;
    if(VtxZ>6.9) vtx_section|=2;
  }
  else if (VtxZ>15.5 && VtxZ<31.0) {//sm2
    if(VtxZ<23.9) vtx_section|=4;
    if(VtxZ>21.8) vtx_section|=8;
  }

  const int nplanks = (int)(ShieldPlankList->GetEntries());
  double DTime[2] = {0., 0.};
  for (int iplnk = 0; iplnk< nplanks; ++iplnk) {
    const AtmosShieldPlank* plank =
      dynamic_cast<const AtmosShieldPlank*>
      (ShieldPlankList->At(iplnk));
    assert(plank);

    //Delta T Vtx
    for (unsigned int iend=0; iend<2; iend++) {
      DTime[iend] = plank->Tcal[iend];
      DTime[iend] -= plank->WlsPigtail[iend]*WlsFibreN/TMath::C();
      DTime[iend] -= plank->ClearFibre[iend]*ClearFibreN/TMath::C();
      DTime[iend] -= TMath::Abs(VtxZ-plank->Z[iend])*WlsFibreN/TMath::C();
    }

    //Add to veto counts if shield hit is in-time with event
    if ((plank->QPE[0]>0 && TMath::Abs(DTime[0]-VtxTime)<TimeWin) ||
        (plank->QPE[1]>0 && TMath::Abs(DTime[1]-VtxTime)<TimeWin)) {

      if ((1<<(plank->Section-1)&(vtx_section))) {
        VetoQCorrSS += plank->QPE[0]+plank->QPE[1];
      }

      bool TrkUnderPlank(false);
      if(plank->Z[0] < plank->Z[1]) 
        TrkUnderPlank =
          (plank->Z[0] < VtxZ && VtxZ < plank->Z[1]);
      else
        TrkUnderPlank =
          (plank->Z[1] < VtxZ && VtxZ < plank->Z[0]);

      double DZ(0.);
      if(!TrkUnderPlank) DZ = TMath::Min(TMath::Abs(plank->Z[0]-VtxZ),
                                         TMath::Abs(plank->Z[1]-VtxZ));
      //ztp.push_back(DZ);
      
      double SteelToPlank(0.);
      double RP(0.);
      if (!TrkUnderPlank) {
        RP = TMath::Sqrt((DZ*DZ) + 
               ((plank->X-VtxX)*(plank->X-VtxX)) +
               ((plank->X-VtxX)*(plank->X-VtxX)) );
        double TrkVtx[3] = {VtxX, VtxY, VtxZ};
        double TrkCos[3];
        TrkCos[0] = (plank->X-VtxX) / RP;
        TrkCos[1] = (plank->Y-VtxY) / RP;
        TrkCos[2] = (DZ) / RP;
        double DetVtx[3];
        int Side = UtilMisc::DetectorWall(TrkVtx, TrkCos, DetVtx);
        if (Side==9 || Side==10) {
          MSG("VetoTrack",Msg::kDebug) << "Cross on SM Gap" << endl;
        }
        SteelToPlank = TMath::Abs((TrkVtx[2]-DetVtx[2])/TrkCos[2]);
        SteelToPlank = SteelToPlank * 2.54 / 5.99;
      }
      //stp.push_back(SteelToPlank);
      new((*InTimePlanks)[NInTimePlanks++]) PlankToTrk(plank, DZ, SteelToPlank);

      /*
      for (int i=0; i<10; i++) {
        if (SteelToPlank <= i*0.1) {
          VetoQPlank10cmSteel[i] += plank->QPE[0]+plank->QPE[1];
        }
      }

      for (int i=0; i<10; i++) {
        if (DZ <= i*0.5) {
          VetoQPlank50cmZ[i] += plank->QPE[0]+plank->QPE[1];
        }
      }
      */
    }
  }
}

---