VetoShower.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/AtmosShower.h"
#include "AtNuEvent/AtmosShieldPlank.h"
#include "MessageService/MsgService.h"

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

CVSID("$Id: VetoShower.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;

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

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

VetoShower::~VetoShower() {
  //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;
}

VetoShower::VetoShower() {
  InTimePlanks = new TClonesArray("PlankToShw");
  Zero();
}

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

VetoShower::VetoShower(const AtmosShower *shower,
                       const TClonesArray *ShieldPlankList) {
  InTimePlanks = new TClonesArray("PlankToShw");
  Zero();
  Fill(shower, ShieldPlankList);
}

void VetoShower::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 VetoShower::Fill(const AtmosEvent *event) {
  //This assumes a single shower over which to find the veto correlation
  const AtmosShower* shower =
    dynamic_cast<const AtmosShower*>(event->ShowerList->At(0));
  if(!shower) return;

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

  Fill(shower, event->ShieldPlankList);
}

void VetoShower::Fill(const AtmosShower *shower,
                      const TClonesArray *ShieldPlankList) {
  assert(shower);
  assert(ShieldPlankList);
  Zero();

  //vector<int> itp;//Temporary for InTimePlanes
  //vector<double> stp;//Temporary for SteelToPlanks
  //vector<double> ztp;//Temporary for ZToPlanks
  //vector<double> pep;//Temporary for PEPlanks

  const double VtxX = shower->VtxX;
  const double VtxY = shower->VtxY;
  const double VtxZ = shower->VtxZ;
  const double VtxTime = shower->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)) {
      SectionQInTime[plank->Section-1] += plank->QPE[0]+plank->QPE[1];

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

      //itp.push_back(plank->Plane);
      //pep.push_back(plank->QPE[0]+plank->QPE[1]);

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

      double DZ(0.);
      if(!ShwUnderPlank) 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 (!ShwUnderPlank) {
        RP = TMath::Sqrt((DZ*DZ) + 
               ((plank->X-VtxX)*(plank->X-VtxX)) +
               ((plank->X-VtxX)*(plank->X-VtxX)) );
        double ShwVtx[3] = {VtxX, VtxY, VtxZ};
        double ShwCos[3];
        ShwCos[0] = (plank->X-VtxX) / RP;
        ShwCos[1] = (plank->Y-VtxY) / RP;
        ShwCos[2] = (DZ) / RP;
        double DetVtx[3];
        int Side = UtilMisc::DetectorWall(ShwVtx, ShwCos, DetVtx);
        if (Side==9 || Side==10) {
          MSG("VetoShower",Msg::kDebug) << "Cross on SM Gap" << endl;
        }
        SteelToPlank = TMath::Abs((ShwVtx[2]-DetVtx[2])/ShwCos[2]);
        SteelToPlank = SteelToPlank * 2.54 / 5.99;
      }
      //stp.push_back(SteelToPlank);
      new((*InTimePlanks)[NInTimePlanks++]) PlankToShw(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];
        }
      }
      */
    }
    else {
      SectionQOutTime[plank->Section-1] += plank->QPE[0]+plank->QPE[1];
    }
  }

  //Put the itp vector in the C-style array InTimePlanes
  /*
  NPlanksInTime = itp.size();
  if (NPlanksInTime > 0) {
    InTimePlanes = new int[NPlanksInTime];
    SteelToPlanks = new double[NPlanksInTime];
    ZToPlanks = new double[NPlanksInTime];
    PEPlanks = new double[NPlanksInTime];
    for (int i=0; i<NPlanksInTime; i++) {
      InTimePlanes[i] = itp[i];
      SteelToPlanks[i] = stp[i];
      ZToPlanks[i] = ztp[i];
      PEPlanks[i] = pep[i];
    }
  }
  */
}

---