// DCAProc.cc
// Contact person: Kevin Nuckolls <knuckolls@berkeley.edu>
// See DCAProc.hh for more details.
//-------------------------------------------------------//

#include <RAT/DCAProc.hh>
#include <CLHEP/Units/SystemOfUnits.h>
#include <RAT/Log.hh>

using namespace std;

namespace RAT
{
  DCAProc::DCAProc() : Processor( "dcaProc" ) { }

  void DCAProc::SetI( const std::string& param, const int value )
  {
    fFlag = 0x0;
    if (param == "bitmask") {
      fFlag = fFlag | value;
      if (fFlag & 0x1) {
        info << "DCAProc: Using plotting parameter slot" << "\n";
      } // Plots slot occupancy for each crate according to nhits
      if (fFlag & 0x2) {
        info << "DCAProc: Using plotting parameter crate" << "\n";
      } // Plots crate occupancy according to nhits
      if (fFlag & 0x4) {
        info << "DCAProc: Using plotting parameter channel" << "\n";
      } // Plots channel occupancy according to electronic space
        // polar angle, and azimuthal angle
      if (fFlag & 0x8) {
        info << "DCAProc: Using plotting parameter flagged" << "\n";
      } // Plots total flagged events and nhit distribution for each cut
      if (fFlag & 0x10) {
        info << "DCAProc: Using plotting parameter timediff" << "\n";
      } // Plots time difference between like-cut events
    } else {
      throw ParamUnknown(param);
    }
  }


  void DCAProc::SetS( const std::string& param, const std::string& value  )
  {
    if (param == "type") {
      if (value == "slot") {
        fFlag_Vector.push_back(0x1);
        info << "DCAProc: Using plotting parameter slot" << "\n";
      } else if (value == "crate") {
        fFlag_Vector.push_back(0x2);
        info << "DCAProc: Using plotting parameter crate" << "\n";
      } else if (value == "channel") {
        fFlag_Vector.push_back(0x4);
        info << "DCAProc: Using plotting parameter channel" << "\n";
      } else if (value == "flagged") {
        fFlag_Vector.push_back(0x8);
        info << "DCAProc: Using plotting parameter flagged" << "\n";
      } else if (value == "timediff") {
        fFlag_Vector.push_back(0x10);
        info << "DCAProc: Using plotting parameter timediff" << "\n";
      } else {
        throw ParamUnknown(param);
      }
    } else if (param == "file") {
      fFilename = value;
      info << "DCAProc: Using output filename " << fFilename.c_str() << "\n";
    } else {
      throw ParamUnknown(param);
    }
  }

  void DCAProc::BeginOfRun( DS::Run& )
  {
    const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
    RAT::DU::DataCleaningBits rDataCleaningBits = RAT::DU::Utility::Get()->GetDataCleaningBits();
    fPMTCount = pmtInfo.GetCount();
    fNumberOfEvents = 0;

    for (size_t j=0; j<fFlag_Vector.size(); j++) {
      fFlag = fFlag | fFlag_Vector[j];
    }

    // Instantiation for slot occupancy
    if (fFlag & 0x1) {
      fSNHITS.resize(4);
      for ( int j=0; j<4; j++) {
        fSNHITS[j].resize(304);
      }
    }

    // Instantiation for crate occupancy
    if (fFlag & 0x2) {
      fCRNHITS.resize(4);
      for ( int j=0; j<4; j++) {
        fCRNHITS[j].resize(19);
      }
    }

    // Instantiation for channel occupancy
    if (fFlag & 0x4) {
      fCHNHITS.resize(12);
      for ( int j=0; j<4; j++) {
        fCHNHITS[j].resize(fPMTCount);
        fCHNHITS[4+j].resize(200);
        fCHNHITS[8+j].resize(360);
      }
      fTOTALCOSINE.resize(200);
      fTOTALPHI.resize(360);
      TVector2 PMTPosFlat;
    }

    // Instantiation for flagged vs nhits
    if (fFlag & 0x8) {
      size_t numberOfCuts = rDataCleaningBits.GetInverseMapLast()->first+1;
      fNHITSCut.resize(numberOfCuts);
      fMaxNhits.resize(numberOfCuts);
      fCutMatrix.resize(numberOfCuts);
      fCutsPerEvent.resize(numberOfCuts);
    for( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        fMaxNhits[iDCBits->first] = 120; // Initial setting, can be expanded while running
        fNHITSCut[iDCBits->first].resize(fMaxNhits[iDCBits->first]);
        fCutMatrix[iDCBits->first].resize(numberOfCuts);
      }
      fTotalCutCount.resize(numberOfCuts);
    }

    // Instantiation for flagged event time difference
    if (fFlag & 0x10) {
      size_t numberOfCuts = rDataCleaningBits.GetInverseMapLast()->first+1;
      fTIMECut.resize(numberOfCuts);
      fMaxTimeDiff.resize(numberOfCuts);
    for( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        fMaxTimeDiff[iDCBits->first] = 120; // Initial setting, can be expanded while running
        fTIMECut[iDCBits->first].resize(fMaxTimeDiff[iDCBits->first]);
      }
      fLastEventTime.resize(numberOfCuts);
    }
  }

  TVector2 DCAProc::TransformCoord(
       const TVector3& V1,
       const TVector3& V2,
       const TVector3& V3,
       const TVector2& A1,
       const TVector2& A2,
       const TVector2& A3,
       const TVector3& P )
  {
    // Transform the coord P relative to V1,2,3 to the
    // returned value relative to A1,2,3
    TVector3 xV = V2 - V1;
    TVector3 yV = ( ( V3 - V1 ) + ( V3 - V2 ) ) * 0.5;
    TVector3 zV = xV.Cross( yV ).Unit();

    double planeD = V1.Dot( zV );
    double t = planeD / P.Dot( zV );

    TVector3 localP = t*P - V1;
    TVector2 xA = A2 - A1;
    TVector2 yA = ( ( A3 - A1 ) +( A3 - A2 ) )  * 0.5;

    double convUnits = xA.Mod() / xV.Mag();

    TVector2 result;
    result = localP.Dot( xV.Unit() ) * xA.Unit() * convUnits;
    result += localP.Dot( yV.Unit() ) * yA.Unit() * convUnits + A1;
    return result;
  }

  void DCAProc::SphereToIcosahedron(
      TVector3& pointOnSphere, // The spherical Position
      TVector2& resultPosition, // The returned position
      const double totalLength, // Total size of the 2d grid (default to 1)
      const double rotation) // Rotation: 0 for "Phil" style, 2.12 for "SNO" style
  {
    pointOnSphere.RotateX(rotation);
    TGraph* gA = new TGraph();
    TGraph2D* gV = new TGraph2D();
    int gAC = 0, gVC = 0;

    const double t = ( 1.0 + sqrt( 5.0 ) ) / 2.0;
    const TVector3 V2 = TVector3( t * t, 0.0, t * t * t ).Unit();
    const TVector3 V6 = TVector3( -t * t, 0.0, t * t * t ).Unit(); // V6
    const TVector3 V12 = TVector3( 0.0, t * t * t, t * t ).Unit(); // V12
    const TVector3 V17 = TVector3( 0.0, -t * t * t, t * t ).Unit(); // V17
    const TVector3 V27 = TVector3( t * t * t, t * t, 0.0 ).Unit(); // V27
    const TVector3 V31 = TVector3( -t * t * t, t * t, 0.0 ).Unit(); // V31
    const TVector3 V33 = TVector3( -t * t * t, -t * t, 0.0 ).Unit(); // V33
    const TVector3 V37 = TVector3( t * t * t, -t * t, 0.0 ).Unit(); // V37
    const TVector3 V46 = TVector3( 0.0, t * t * t, -t * t ).Unit(); // V46
    const TVector3 V51 = TVector3( 0.0, -t * t * t, -t * t ).Unit(); // V51
    const TVector3 V54 = TVector3( t * t, 0.0, -t * t * t ).Unit(); // V54
    const TVector3 V58 = TVector3( -t * t, 0.0, -t * t * t ).Unit(); // V58

    gV->SetPoint( gVC, V2.X(), V2.Y(), V2.Z() ); gVC++;
    gV->SetPoint( gVC, V6.X(), V6.Y(), V6.Z() ); gVC++;
    gV->SetPoint( gVC, V12.X(), V12.Y(), V12.Z() ); gVC++;
    gV->SetPoint( gVC, V17.X(), V17.Y(), V17.Z() ); gVC++;
    gV->SetPoint( gVC, V27.X(), V27.Y(), V27.Z() ); gVC++;
    gV->SetPoint( gVC, V31.X(), V31.Y(), V31.Z() ); gVC++;
    gV->SetPoint( gVC, V33.X(), V33.Y(), V33.Z() ); gVC++;
    gV->SetPoint( gVC, V37.X(), V37.Y(), V37.Z() ); gVC++;
    gV->SetPoint( gVC, V46.X(), V46.Y(), V46.Z() ); gVC++;
    gV->SetPoint( gVC, V51.X(), V51.Y(), V51.Z() ); gVC++;
    gV->SetPoint( gVC, V54.X(), V54.Y(), V54.Z() ); gVC++;
    gV->SetPoint( gVC, V58.X(), V58.Y(), V58.Z() ); gVC++;

    // Faces {{ 2, 6, 17}, { 2, 12, 6}, { 2, 17, 37}, { 2, 37, 27}, { 2, 27, 12}, {37, 54, 27},
    // {27, 54, 46}, {27, 46, 12}, {12, 46, 31}, {12, 31, 6}, { 6, 31, 33}, { 6, 33, 17},
    // {17, 33, 51}, {17, 51, 37}, {37, 51, 54}, {58, 54, 51}, {58, 46, 54}, {58, 31, 46},
    // {58, 33, 31}, {58, 51, 33}}
    vector<TVector3> IcosahedronCentres;

    IcosahedronCentres.push_back( ( V2 + V6 + V17 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V12 + V6 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V17 + V37 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V37 + V27 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V27 + V12 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V37 + V54 + V27 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V27 + V54 + V46 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V27 + V46 + V12 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V12 + V46 + V31 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V12 + V31 + V6 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V6 + V31 + V33 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V6 + V33 + V17 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V17 + V33 + V51 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V17 + V51 + V37 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V37 + V51 + V54 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V54 + V51 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V46 + V54 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V31 + V46 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V33 + V31 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V51 + V33 ) * ( 1.0 / 3.0 ) );

    vector<double> distFromCentre;
    unsigned int uLoop;
    for( uLoop = 0; uLoop < IcosahedronCentres.size(); uLoop++ )
      distFromCentre.push_back( ( IcosahedronCentres[uLoop] - pointOnSphere ).Mag() );

    int face = min_element( distFromCentre.begin(), distFromCentre.end() ) - distFromCentre.begin() + 1;
    double a = totalLength / 5.5;
    double b = a * sqrt( 3.0 ) / 2.0;
    TVector2 A12a = TVector2( a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12a.X(), A12a.Y() ); gAC++;
    TVector2 A12b = TVector2( 3.0 * a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12b.X(), A12b.Y() ); gAC++;
    TVector2 A12c = TVector2( 5.0 * a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12c.X(), A12c.Y() ); gAC++;
    TVector2 A12d = TVector2( 7.0 *a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12d.X(), A12d.Y() ); gAC++;
    TVector2 A12e = TVector2( 9.0 * a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12e.X(), A12e.Y() ); gAC++;
    TVector2 A2a = TVector2( 0.0, b );
    gA->SetPoint( gAC, A2a.X(), A2a.Y() ); gAC++;
    TVector2 A2b = TVector2( 5.0 * a, b );
    gA->SetPoint( gAC, A2b.X(), A2b.Y() ); gAC++;
    TVector2 A17a = TVector2( a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A17a.X(), A17a.Y() ); gAC++;
    TVector2 A17b = TVector2( 11.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A17b.X(), A17b.Y() ); gAC++;
    TVector2 A51a = TVector2( a, 3.0 * b );
    gA->SetPoint( gAC, A51a.X(), A51a.Y() ); gAC++;
    TVector2 A51b = TVector2( 2.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51b.X(), A51b.Y() ); gAC++;
    TVector2 A51c = TVector2( 3.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51c.X(), A51c.Y() ); gAC++;
    TVector2 A51d = TVector2( 4.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51d.X(), A51d.Y() ); gAC++;
    TVector2 A51e = TVector2( 5.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51e.X(), A51e.Y() ); gAC++;
    TVector2 A27 = TVector2( 4.0 * a, b );
    gA->SetPoint( gAC, A27.X(), A27.Y() ); gAC++;
    TVector2 A46 = TVector2( 3.0 * a, b );
    gA->SetPoint( gAC, A46.X(), A46.Y() ); gAC++;
    TVector2 A31 = TVector2( 2.0 * a, b );
    gA->SetPoint( gAC, A31.X(), A31.Y() ); gAC++;
    TVector2 A6 = TVector2( a, b );
    gA->SetPoint( gAC, A6.X(), A6.Y() ); gAC++;
    TVector2 A37 = TVector2( 9.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A37.X(), A37.Y() ); gAC++;
    TVector2 A33 = TVector2( 3.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A33.X(), A33.Y() ); gAC++;
    TVector2 A58 = TVector2( 5.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A58.X(), A58.Y() ); gAC++;
    TVector2 A54 = TVector2( 7.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A54.X(), A54.Y() ); gAC++;

    TVector3 xV, yV;
    TVector2 xA, yA;
    switch(face)
      {
      case 1:// {2, 6, 17}
        resultPosition = DCAProc::TransformCoord( V2, V6, V17, A2a, A6, A17a, pointOnSphere );
        break;
      case 2:// {2, 12, 6}
        resultPosition = DCAProc::TransformCoord( V2, V12, V6, A2a, A12a, A6, pointOnSphere );
        break;
      case 3:// {2, 17, 37}
        resultPosition = DCAProc::TransformCoord( V2, V17, V37, A2b, A17b, A37, pointOnSphere );
        break;
      case 4:// {2, 37, 27}
        resultPosition = DCAProc::TransformCoord( V2, V37, V27, A2b, A37, A27, pointOnSphere );
        break;
      case 5:// {2, 27, 12}
        resultPosition = DCAProc::TransformCoord( V2, V27, V12, A2b, A27, A12e, pointOnSphere );
        break;
      case 6:// {37, 54, 27}
        resultPosition = DCAProc::TransformCoord( V37, V54, V27, A37, A54, A27, pointOnSphere );
        break;
      case 7:// {27, 54, 46}
        resultPosition = DCAProc::TransformCoord( V27, V54, V46, A27, A54, A46, pointOnSphere );
        break;
      case 8:// {27, 46, 12}
        resultPosition = DCAProc::TransformCoord( V27, V46, V12, A27, A46, A12d, pointOnSphere );
        break;
      case 9:// {12, 46, 31}
        resultPosition = DCAProc::TransformCoord( V12, V46, V31, A12c, A46, A31, pointOnSphere );
        break;
      case 10:// {12, 31, 6}
        resultPosition = DCAProc::TransformCoord( V12, V31, V6, A12b, A31, A6, pointOnSphere );
        break;
      case 11:// {6, 31, 33}
        resultPosition = DCAProc::TransformCoord( V6, V31, V33, A6, A31, A33, pointOnSphere );
        break;
      case 12:// { 6, 33, 17}
        resultPosition = DCAProc::TransformCoord( V6, V33, V17, A6, A33, A17a, pointOnSphere );
        break;
      case 13:// {17, 33, 51}
        resultPosition = DCAProc::TransformCoord( V17, V33, V51, A17a, A33, A51a, pointOnSphere );
        break;
      case 14:// {17, 51, 37}
        resultPosition = DCAProc::TransformCoord( V17, V51, V37, A17b, A51e, A37, pointOnSphere );
        break;
      case 15:// {37, 51, 54}
        resultPosition = DCAProc::TransformCoord( V37, V51, V54, A37, A51d, A54, pointOnSphere );
        break;
      case 16:// {58, 54, 51}
        resultPosition = DCAProc::TransformCoord( V58, V54, V51, A58, A54, A51c, pointOnSphere );
        break;
      case 17:// {58, 46, 54}
        resultPosition = DCAProc::TransformCoord( V58, V46, V54, A58, A46, A54, pointOnSphere );
        break;
      case 18:// {58, 31, 46}
        resultPosition = DCAProc::TransformCoord( V58, V31, V46, A58, A31, A46, pointOnSphere );
        break;
      case 19:// {58, 33, 31}
        resultPosition = DCAProc::TransformCoord( V58, V33, V31, A58, A33, A31, pointOnSphere );
        break;
      case 20:// {58, 51, 33}
        resultPosition = DCAProc::TransformCoord( V58, V51, V33, A58, A51b, A33, pointOnSphere );
        break;
      }
    return;
  }

  DCAProc::~DCAProc()
  {
    TFile *fCreateFile = TFile::Open(fFilename.c_str(),"RECREATE");
    fCreateFile->cd();
    const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
    fPMTCount = pmtInfo.GetCount();
    static const int arr[] = {10,20,40,60};
    vector<int> nhitcuts (arr, arr + sizeof(arr) / sizeof(arr[0]) );

    if (fFlag & 0x1) { // Slot Occupancy, Type:"slot"
      vector<TH1F*> sHistList;

      // Set up slot occupancy histograms
      for ( int j=0; j<4; j++ ) {
        for ( int k=0; k<19; k++) {
          std::ostringstream _crate;
          _crate << k;
          std::string crate = _crate.str();
          std::ostringstream _nhitcut;
          _nhitcut << nhitcuts[j];
          std::string nhitcut = _nhitcut.str();
          std::ostringstream _totalEvents;
          _totalEvents << fNumberOfEvents;
          std::string totalEvents = _totalEvents.str();
          std::string name = "Nhit_" + nhitcut + ",Crate_" + crate;
          std::string title = "Slot Occupancy Nhit>" + nhitcut + " Crate "
                              + crate + ", Total Events Analyzed: " + totalEvents;
          sHistList.push_back(new TH1F(name.c_str(),title.c_str(),16,0,16));
        }
      }

      // Fill "Slot Occupancy" histograms
      for (int j=0; j<4; j++) { // Each Crate
        for (int k=0; k<19; k++) { // Each NHit Threshold
          for (int l=0; l<16; l++) { // Each Slot in Crate
            sHistList[19*j+k]->Fill(l,fSNHITS[j][(16*k)+l]);
          }
        }
      }

      for (int j=0; j<76; j++) {
        sHistList[j]->GetXaxis()->SetTitle("Slot Number within Crate");
        sHistList[j]->GetYaxis()->SetTitle("Number of Hits within Slot");
        sHistList[j]->SetStats(0);
        sHistList[j]->Write();
      }
    }

    if (fFlag & 0x2) { // Crate Occupancy, Type:"slot"
      std::vector<TH1F*> crHistList;

      // Set up crate occupancy histograms
      for ( int j=0; j<4; j++ ) {
        std::ostringstream _nhitcut;
        _nhitcut << nhitcuts[j];
        std::string nhitcut = _nhitcut.str();
        std::ostringstream _totalEvents;
        _totalEvents << fNumberOfEvents;
        std::string totalEvents = _totalEvents.str();
        std::string name = "Nhit_" + nhitcut + ",All_Crate";
        std::string title = "Crate Occupancy Nhit>" + nhitcut
                             + ", Total Events Analyzed: " + totalEvents;
        crHistList.push_back(new TH1F(name.c_str(),title.c_str(),19,0,19));
      }

      // Fill "Crate Occupancy" histograms
      for (int j=0; j<4; j++) { // Each NHit Threshold
        for (int k=0; k<19; k++) { // Each Crate
        crHistList[j]->Fill(k,fCRNHITS[j][k]);
        }
      }

      for (int j=0; j<4; j++) {
        crHistList[j]->GetXaxis()->SetTitle("Crate Number");
        crHistList[j]->GetYaxis()->SetTitle("Number of Hits within Crate");
        crHistList[j]->SetStats(0);
        crHistList[j]->Write();
      }
    }

    if (fFlag & 0x4) { // Channel Occupancy, Type:"channel"
      std::vector<TH1F*> chHistList;
      static const string arr2[] = {"Electronic_Space","Cos_Theta","Phi","Percent_Cos","Percent_Phi"};
      vector<string> chanoccnames (arr2, arr2 + sizeof(arr2) / sizeof(arr2[0]) );
      static const int arr3[] = {304,0,fPMTCount,40,-1,1,72,0,360,40,-1,1,72,0,360};
      vector<int> bins_bounds (arr3, arr3 + sizeof(arr3) / sizeof(arr3[0]) );

      // Set up histograms for channel occupancy
      for ( int j=0; j<5; j++ ) {
        for ( int k=0; k<4; k++) {
          std::ostringstream _channame;
          _channame << chanoccnames[j];
          std::string channame = _channame.str();
          std::ostringstream _nhitcut;
          _nhitcut << nhitcuts[k];
          std::string nhitcut = _nhitcut.str();
          std::ostringstream _totalEvents;
          _totalEvents << fNumberOfEvents;
          std::string totalEvents = _totalEvents.str();
          std::string name = "Nhit_" + nhitcut + channame;
          std::string title = "Channel Occupancy Nhit>" + nhitcut + " by "
                              + channame + ", Total Events Analyzed: " + totalEvents;
          chHistList.push_back(new TH1F(name.c_str(),title.c_str(),bins_bounds[3*j],
                               bins_bounds[3*j+1],bins_bounds[3*j+2]));
        }
      }

      // Assemble vectors of where all PMTs are in order to normalize later
      for (int ip=0; ip<fPMTCount; ip++) {
        TVector3 totposvect = pmtInfo.GetPosition(ip);
        int totcosine = int((100*(totposvect[2])/(sqrt(totposvect[0]*totposvect[0]
                        +totposvect[1]*totposvect[1]+totposvect[2]*totposvect[2])))+100);
        if (totcosine < 200 and totcosine >=0 ) {
          fTOTALCOSINE[totcosine]++;
        }
        int totphi = int(acos(totposvect[0]/(sqrt(totposvect[0]*totposvect[0]
                    +totposvect[1]*totposvect[1])))*57.325);
        if (totposvect[1] < 0) {
          totphi = 360 - totphi;
        }
        if (totphi < 360 and totphi >= 0 ) {
          fTOTALPHI[totphi]++;
        }
      }

      // Fill "Channel Occupancy: Electronic Space" histograms
      for (int j=0; j<4; j++) { // Each Nhit Threshold
        for (int k=0; k<fPMTCount; k++) { // Each PMT ID
          chHistList[j]->Fill(k,fCHNHITS[j][k]);
        }
      }
      // Fill "Channel Occupancy: Cosine Theta" histograms
      for (int j=0; j<4; j++) { // Each Nhit Threshold
        for (int k=0; k<200; k++) { // Each Bin of Cos(theta)
          chHistList[4+j]->Fill(0.01*(k-100),fCHNHITS[4+j][k]);
        }
      }
      // Fill "Channel Occupancy: Phi" histograms
      for (int j=0; j<4; j++) { // Each Nhit Threshold
        for (int k=0; k<360; k++) { // Each Degree of Phi
          chHistList[8+j]->Fill(k,fCHNHITS[8+j][k]);
        }
      }
      // Fill "Channel Occupancy: Normalized Cosine Theta" histograms
      for (int j=0; j<4; j++) { // Each Nhit Threshold
        for (int k=0; k<200; k++) { // Each Bin of Cos(theta)
          if (fTOTALCOSINE[k] > 0) {
            chHistList[12+j]->Fill(0.01*(k-100),fCHNHITS[4+j][k]/fTOTALCOSINE[k]);
          } else {
            chHistList[12+j]->Fill(0.01*(k-100),0);
          }
        }
      }
      // Fill "Channel Occupancy: Normalized Phi" histograms
      for (int j=0; j<4; j++) { // Each Nhit Threshold
        for (int k=0; k<360; k++) { // Each Degree of Phi
          if (fTOTALPHI[k] > 0) {
            chHistList[16+j]->Fill(k,fCHNHITS[8+j][k]/fTOTALPHI[k]);
          } else {
            chHistList[16+j]->Fill(k,0);
          }
        }
      }

      for (int j=0; j<4; j++) {
        chHistList[j]->GetXaxis()->SetTitle("PMT ID");
        chHistList[4+j]->GetXaxis()->SetTitle("Cosine of Polar Angle");
        chHistList[8+j]->GetXaxis()->SetTitle("Azimuthal Angle (in degrees)");
        chHistList[12+j]->GetXaxis()->SetTitle("Cosine of Polar Angle");
        chHistList[16+j]->GetXaxis()->SetTitle("Azimuthal Angle (in degrees)");
        chHistList[j]->GetYaxis()->SetTitle("Number of Hits per Channel");
        chHistList[4+j]->GetYaxis()->SetTitle("Number of Hits per Bin");
        chHistList[8+j]->GetYaxis()->SetTitle("Number of Hits per Bin");
        chHistList[12+j]->GetYaxis()->SetTitle("Normalized Number of Hits per Bin");
        chHistList[16+j]->GetYaxis()->SetTitle("Normalized Number of Hits per Bin");
      }
      for (int j=0; j<20; j++) {
        chHistList[j]->SetStats(0);
        chHistList[j]->Write();
      }

      // Setting up "Channel Occupancy: Flat Map"
      TCanvas c2("Flat_Map","Flat Map Hit Distribution",10,10,800,600);
      TGraph2D *coverage = new TGraph2D();
      TView3D *view = new TView3D();
      c2.cd();
      coverage->Draw("zcolpcol");
      view->RotateView(270,0);
      TVector2 PMTPosFlat;

      // Transform the 3D position of the pmt (in the TVector3 PMTPos) to a 2D position (fill PMTPosFlat)
      // Fill the 2D graph, in my case, the z-axis is the number of hits
      for (int j=0; j<fPMTCount; j++) {
        TVector3 PMTPos = pmtInfo.GetPosition(j);
        DCAProc::SphereToIcosahedron(PMTPos,PMTPosFlat,1.0,2.12);
        coverage->SetPoint(j,PMTPosFlat.X(),PMTPosFlat.Y(),fCHNHITS[0][j]);
      }

      // Adjust aesthetics of flat map
      coverage->SetMarkerStyle(8);
      coverage->SetMarkerSize(0.5);
      gStyle->SetLabelSize(0.02);
      coverage->GetXaxis()->SetLabelSize(0.03);
      coverage->GetYaxis()->SetLabelSize(0.03);
      std::ostringstream _totalEvents;
      _totalEvents << fNumberOfEvents;
      std::string totalEvents = _totalEvents.str();
      std::string FlatMapTitle = "Flat Map Hit Distribution, Total Events Analyzed: " + totalEvents;
      coverage->SetTitle(FlatMapTitle.c_str());
      coverage->Draw("zcolpcol");
      c2.Update();
      c2.Write();
    }

    if (fFlag & 0x8) { // Flagged Events vs. Nhit Analysis, Type:"flagged"
      std::vector<TH1F*> flHistList;
      RAT::DU::DataCleaningBits rDataCleaningBits = RAT::DU::Utility::Get()->GetDataCleaningBits();
      size_t numberOfCuts = rDataCleaningBits.GetInverseMapLast()->first+1;
      size_t cutsWithoutWB = 17;

      // Set up histograms for flagged vs.nhit
      for( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        std::ostringstream _totalEvents;
        _totalEvents << fNumberOfEvents;
        std::string totalEvents = _totalEvents.str();
        std::string name = rDataCleaningBits.GetBitName(iDCBits->first) + "_nhit";
        std::string title = "Flagged Events:" + rDataCleaningBits.GetBitName(iDCBits->first) + " vs. Nhits"
                              + ", Total Events Analyzed: " + totalEvents;
        flHistList.push_back(new TH1F(name.c_str(),title.c_str(),
          fMaxNhits[iDCBits->first]/2,0,fMaxNhits[iDCBits->first]));
      }
      std::ostringstream _totalEvents;
      _totalEvents << fNumberOfEvents;
      std::string totalEvents = _totalEvents.str();
      std::string TotalFlaggedTitle = "Total Flagged Events for Each Cut, Total Events Analyzed: "
                                     + totalEvents;
      flHistList.push_back(new TH1F("TotalFlagged",TotalFlaggedTitle.c_str(),numberOfCuts,0,numberOfCuts));
      std::string CutsPerEventTitle = "Number of Cuts per Event, Total Events Analyzed: "
                                     + totalEvents;
      flHistList.push_back(new TH1F("CutsPerEvent",CutsPerEventTitle.c_str(),cutsWithoutWB,1,cutsWithoutWB));
      std::string DC_Cut_Matrix_Title = "Data Cleaning Cut Matrix, Total Events Analyzed: " + totalEvents;
      TH2F* DC_Cut_Matrix = new TH2F("Cut_Matrix",DC_Cut_Matrix_Title.c_str(),numberOfCuts,0,numberOfCuts,
                                    numberOfCuts,0,numberOfCuts);

      // Fill "Flagged vs. Nhit" histograms
      for( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        for (int k=0; k<fMaxNhits[iDCBits->first]; k++) {
          flHistList[iDCBits->first]->Fill(k,fNHITSCut[iDCBits->first][k]);
        }
      }

      // Label and fill "Total Flagged" histogram, label other histograms, fill CutsPerEvent histogram
      for ( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        flHistList[numberOfCuts]->GetXaxis()->SetBinLabel(iDCBits->first+1,rDataCleaningBits.GetBitName(iDCBits->first).c_str());
        flHistList[numberOfCuts]->Fill(iDCBits->first,fTotalCutCount[iDCBits->first]);
        DC_Cut_Matrix->GetXaxis()->SetBinLabel(iDCBits->first+1,rDataCleaningBits.GetBitName(iDCBits->first).c_str());
        DC_Cut_Matrix->GetYaxis()->SetBinLabel(iDCBits->first+1,rDataCleaningBits.GetBitName(iDCBits->first).c_str());
        flHistList[iDCBits->first]->GetXaxis()->SetTitle("NHit Value");
        flHistList[iDCBits->first]->GetYaxis()->SetTitle("Number of Events per Bin");
        flHistList[iDCBits->first]->SetStats(0);
        flHistList[iDCBits->first]->Write();
        flHistList[numberOfCuts+1]->Fill(iDCBits->first,fCutsPerEvent[iDCBits->first]);
      }
      flHistList[numberOfCuts]->GetYaxis()->SetTitle("Number of Events per Cut");
      flHistList[numberOfCuts]->SetStats(0);
      for ( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        for ( map<size_t, string>::iterator jDCBits = rDataCleaningBits.GetInverseMapBegin();
             jDCBits != rDataCleaningBits.GetInverseMapEnd(); jDCBits++ ) {
          DC_Cut_Matrix->Fill(iDCBits->first,jDCBits->first,
                                            fCutMatrix[iDCBits->first][jDCBits->first]);
        }
      }

      std::string c1Title = "Cut Matrix, Total Events Analyzed: " + totalEvents;
      TCanvas c1("Data_Cleaning_Cut_Matrix",c1Title.c_str(),10,10,800,600);
      std::string c3Title = "Total Flagged Events for Each Cut, Total Events Analyzed: " + totalEvents;
      TCanvas c3("TotalFlaggedEvents",c3Title.c_str(),10,10,800,600);
      c1.cd();
      c1.SetGridx(1);
      c1.SetGridy(1);
      DC_Cut_Matrix->Draw("text zcol");
      DC_Cut_Matrix->SetStats(0);
      DC_Cut_Matrix->GetXaxis()->SetRangeUser(0,cutsWithoutWB);
      DC_Cut_Matrix->GetYaxis()->SetRangeUser(0,cutsWithoutWB);
      c1.Update();
      c1.Write();
      c3.cd();
      flHistList[numberOfCuts]->GetYaxis()->SetRangeUser(0,1.02*fNumberOfEvents);
      flHistList[numberOfCuts]->Draw();
      TGaxis *axis = new TGaxis(numberOfCuts*gPad->GetUxmax(),gPad->GetUymin(), // Create new percentage axis
        numberOfCuts*gPad->GetUxmax(),1.02*fNumberOfEvents*gPad->GetUymax(),0,1.02,20,"+L");
      axis->SetTitle("Fraction of Events per Cut");
      axis->Draw();
      c3.Update();
      c3.Write();
      flHistList[numberOfCuts+1]->GetXaxis()->SetTitle("Number of Cuts per Event");
      flHistList[numberOfCuts+1]->GetYaxis()->SetTitle("Number of Events");
      flHistList[numberOfCuts+1]->SetStats(0);
      flHistList[numberOfCuts+1]->Write();
    }

    if (fFlag & 0x10) { // Time Difference Between Flagged Events, Type:"timediff"
      std::vector<TH1F*> tHistList;
      RAT::DU::DataCleaningBits rDataCleaningBits = RAT::DU::Utility::Get()->GetDataCleaningBits();

      // Set up histograms for time difference
      for( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        std::ostringstream _totalEvents;
        _totalEvents << fNumberOfEvents;
        std::string totalEvents = _totalEvents.str();
        std::string name = rDataCleaningBits.GetBitName(iDCBits->first) + "_time";
        std::string title = "Time Difference Between Events Flagged by "
                             + rDataCleaningBits.GetBitName(iDCBits->first)
                             + ", Total Events Analyzed: " + totalEvents;
        tHistList.push_back(new TH1F(name.c_str(),title.c_str(),
          fMaxTimeDiff[iDCBits->first]/2,0,fMaxTimeDiff[iDCBits->first]));
      }

      // Fill and label "Time Difference Between Flagged Events" histograms
      for( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
           iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
        for (int k=0; k<fMaxTimeDiff[iDCBits->first]; k++) {
          tHistList[iDCBits->first]->Fill(k,fTIMECut[iDCBits->first][k]);
        }
        tHistList[iDCBits->first]->GetXaxis()->SetTitle("Time Difference (in tenths of seconds)");
        tHistList[iDCBits->first]->GetYaxis()->SetTitle("Number of Events per Bin");
        tHistList[iDCBits->first]->SetStats(0);
        tHistList[iDCBits->first]->Write();
      }
    }
    // Save and close file
    fCreateFile->Close();
  }

  Processor::Result DCAProc::DSEvent(DS::Run&, DS::Entry& ds)
  {
    int nevC = ds.GetEVCount();
    if (nevC > 0) {
      fNumberOfEvents++;
      RAT::DS::EV& ev = ds.GetEV(0);
      int npmtCal = ev.GetCalPMTs().GetCount();
      static const int arr[] = {10,20,40,60};
      vector<int> nhitcuts (arr, arr + sizeof(arr) / sizeof(arr[0]) );

      if (fFlag & 0x3) { // Slot or Crate Occupancy
        for (int ip=0; ip<npmtCal; ip++) {
          RAT::DS::PMTCal& pmt = ev.GetCalPMTs().GetPMT(ip);
          int cratenumber = fBits->GetCrate(pmt.GetID());
          int slotnumber = fBits->GetCard(pmt.GetID());

          // Fill vectors for slot and crate occupancy
          for (int j=0; j<4; j++) {
            if (npmtCal >= nhitcuts[j]) {
              if (fFlag & 0x1) { // Slot occupancy
                fSNHITS[j][cratenumber*16 + slotnumber]++;
              } else if (fFlag & 0x2) { // Crate occupancy
                fCRNHITS[j][cratenumber]++;
              }
            }
          }
        }
      }

      if (fFlag & 0x4) { // Channel Occupancy
        const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
        for (int ip=0; ip<npmtCal; ++ip) {
          RAT::DS::PMTCal& pmt = ev.GetCalPMTs().GetPMT(ip);
          TVector3 posvect = pmtInfo.GetPosition(pmt.GetID());

          // Calculate cosine of polar angle
          int cosine = int((100*(posvect[2])/(sqrt(posvect[0]*posvect[0]
                        +posvect[1]*posvect[1]+posvect[2]*posvect[2])))+100);
          // Calculate azimuthal angle
          int phi = int(acos(posvect[0]/(sqrt(posvect[0]*posvect[0]
                    +posvect[1]*posvect[1])))*57.325);
          // Arc cosine returns elements of [0,Pi], so extend range to [0,2*Pi]
          if (posvect[1] < 0) {
            phi = 360 - phi;
          }

          // Fill vectors for channel occupancy
          for (int j=0; j<4; j++) {
            if (npmtCal >= nhitcuts[j]) {
              fCHNHITS[j][pmt.GetID()]++;
              fCHNHITS[4+j][cosine]++;
              fCHNHITS[8+j][phi]++;
            }
          }
        }
      }

      if (fFlag & 0x18) { // Flagged vs. Nhits and Flagged Event Time Difference
        double scalingfactor = 0.1; // Used to neatly visualize time difference histogram
        RAT::DU::DataCleaningBits rDataCleaningBits = RAT::DU::Utility::Get()->GetDataCleaningBits();
        size_t numberOfCuts = rDataCleaningBits.GetInverseMapLast()->first+1;
        const RAT::DS::DataQCFlags& rDataQCFlags = ev.GetDataCleaningFlags();
        // Log current time for possible replacement of the last event time for cuts
        double currentEventTime = (double)(ev.GetUniversalTime().GetNanoSeconds()/1000000000)
                                  + (double)(ev.GetUniversalTime().GetSeconds());
        vector<int> numberEventCuts;
        numberEventCuts.resize(numberOfCuts);
        if ( rDataQCFlags.ExistFlags(0) ) { // Checks flags exist
          const RAT::DS::BitMask& rBitMaskFlags = rDataQCFlags.GetFlags(0);
          const RAT::DS::BitMask& rBitMaskApplied = rDataQCFlags.GetApplied(0);

          // Loop through data cleaning bits
          for ( map<size_t, string>::iterator iDCBits = rDataCleaningBits.GetInverseMapBegin();
                iDCBits != rDataCleaningBits.GetInverseMapEnd(); iDCBits++ ) {
            Bool_t rApplied = rBitMaskApplied.Get( iDCBits->first );
            if ( rApplied ) { // Checks cut has been applied
              Bool_t rFlag = rBitMaskFlags.Get( iDCBits->first );
              if( !rFlag ) { // Fill if event fails the cut
                if (fFlag & 0x8) { // Flagged vs. Nhits
                  if (npmtCal > fMaxNhits[iDCBits->first]) { // Resizing range if necessary
                    fMaxNhits[iDCBits->first] = npmtCal + 1;
                    fNHITSCut[iDCBits->first].resize(fMaxNhits[iDCBits->first]);
                  }
                  fNHITSCut[iDCBits->first][npmtCal]++; // Keep track of nhit value for cut
                  fTotalCutCount[iDCBits->first]++; // Keep track of which cut from which this event came
                  numberEventCuts[iDCBits->first]++;
                }
                if (fFlag & 0x10) { // Flagged Event Time Difference
                  if (fLastEventTime[iDCBits->first] != 0) { // Check if this is the first event cut
                                                             // by this particular cut
                    int timeDifference = (int)((currentEventTime - fLastEventTime[iDCBits->first])/scalingfactor);
                    if (timeDifference > fMaxTimeDiff[iDCBits->first]) { // Resizing range if necessary
                      fMaxTimeDiff[iDCBits->first] = (timeDifference + 1);
                      fTIMECut[iDCBits->first].resize(fMaxTimeDiff[iDCBits->first]);
                    }
                    fTIMECut[iDCBits->first][timeDifference]++; // Keep track of time value for cut
                  }
                  fLastEventTime[iDCBits->first] = currentEventTime; // Log new last seen time
                }
              }
            }
          }
          for (size_t j=0; j<(numberOfCuts-1); j++) { // Checking each pair of cuts to find correlations
            if (numberEventCuts[j] == 1) {
              for (size_t k=1; k<(numberOfCuts - j); k++) {
                if (numberEventCuts[numberOfCuts - k] == 1) {
                  fCutMatrix[j][numberOfCuts - k]++;
                  fCutMatrix[numberOfCuts - k][j]++;
                }
              }
            }
          }
          numberEventCuts.resize(17); // Pairs only matter between non-Water Blind cuts
          fCutsPerEvent[std::accumulate(numberEventCuts.begin(),numberEventCuts.end(),0)]++;
        }
      }
    }
  return OK;
  }
} //namespace RAT