RegPixelMap.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
///// \file    RegPixelMap.h
///// \brief   RegPixelMap for Regression modified from Dominick's PixelMap
///// \author  Jianming Bian - bianjm@uci.edu
//////////////////////////////////////////////////////////////////////////

#include <cassert>
#include <iostream>
#include "CVN/func/RegPixelMap.h"

namespace cvn
{

  RegPixelMap::RegPixelMap(unsigned int nPlane, unsigned int nCell,
                     const Boundary& bound):
  fNPlane(nPlane),
  fNCell(nCell),
  fPE(nPlane*nCell),
  fPEX(nPlane*nCell),
  fPEY(nPlane*nCell),
  fPur(nPlane*nCell),
  fPurX(nPlane*nCell),
  fPurY(nPlane*nCell),
  fLab(nPlane*nCell),
  fLabX(nPlane*nCell),
  fLabY(nPlane*nCell),
  fBound(bound)
  {}


  RegPixelMap::RegPixelMap(unsigned int nPlane, unsigned int nCell,
                     unsigned int vtxPlane, unsigned int vtxXCell, unsigned int vtxYCell,
                     const Boundary& bound):
  fNPlane(nPlane),
  fNCell(nCell),
  fVtxPlane(vtxPlane),
  fVtxXCell(vtxXCell),
  fVtxYCell(vtxYCell),
  fPE(nPlane*nCell),
  fPEX(nPlane*nCell),
  fPEY(nPlane*nCell),
  fPur(nPlane*nCell),
  fPurX(nPlane*nCell),
  fPurY(nPlane*nCell),
  fLab(nPlane*nCell),
  fLabX(nPlane*nCell),
  fLabY(nPlane*nCell),
  fBound(bound)
  {}



  void RegPixelMap::FillInputVector(float* input) const
  {
    unsigned int i = 0;

    for(const auto& pe:fPE){
      input[i] = pe;
      ++i;
    }

  }


  void RegPixelMap::Add(const unsigned int& plane, const unsigned int& cell,
                     const float& pe, const HitType label = kEmptyHit, const double purity=0.0 )
  {
    unsigned int view = plane % 2;
    if(fBound.IsWithin(plane, cell)){
      fPE[GlobalToIndex(plane,cell)] += pe/2100.;
      fLab[GlobalToIndex(plane,cell)] = label;
      fPur[GlobalToIndexSingle(plane,cell)] = purity;
      if(view==0){
        fPEX[GlobalToIndexSingle(plane,cell)] += pe/2100.;//Why +=?
        fLabX[GlobalToIndexSingle(plane,cell)] = label;
        fPurX[GlobalToIndexSingle(plane,cell)] = purity;
      }
      else{
        fPEY[GlobalToIndexSingle(plane,cell)] += pe/2100.;
        fLabY[GlobalToIndexSingle(plane,cell)] = label;
        fPurY[GlobalToIndexSingle(plane,cell)] = purity;
      }
    }
  }

  void RegPixelMap::Add(const unsigned int plane, const unsigned int cell, const int view,
                     const float pe, const HitType label = kEmptyHit, const double purity=0.0)
  { 
    if(fBound.IsWithin(plane, cell)){
      fPE[GlobalToIndex(plane,cell)] += pe/2100.;
      fLab[GlobalToIndex(plane,cell)] = label;
      fPur[GlobalToIndexSingle(plane,cell)] = purity;
    }
    int localplaneid = plane - (int)fVtxPlane;
    if(view==0){
      int localcellid = cell - (int)fVtxXCell;
      if(localplaneid>=-30&&localplaneid<=120&&localcellid>=-70&&localcellid<=70){
        int index = 141*(localplaneid+30) + (localcellid+70);
        if(index>=0&&index<(int)fPEX.size()){
//    std::cout<<"in PixelMap, view, localplane, localcell, plane, cell, index,  pe*50 "<<view<<" "<<localplaneid<<" "<<localcellid<<" "<<(localplaneid+30)/2<<" "<<localcellid+70<<" "<<index<<" "<<pe*50<<std::endl;
          fPEX[index] = pe;
          fLabX[index] = label;
          fPurX[index] = purity;
        }
      }
    }
    if(view==1) {
      int localcellid = cell - (int)fVtxYCell;
      if(localplaneid>=-30&&localplaneid<=120&&localcellid>=-70&&localcellid<=70){
        int index = 141*(localplaneid+30) + (localcellid+70);
        if(index>=0&&index<(int)fPEY.size()){
//    std::cout<<"in PixelMap, view, plane, cell, index, pe*50 "<<view<<" "<<localplaneid<<" "<<localcellid<<" "<<(localplaneid+30-1)/2<<" "<<localcellid+70<<" "<<index<<" "<<pe*50<<std::endl;
          fPEY[index] = pe;
          fLabY[index] = label;
          fPurY[index] = purity;
        }
      }
    }
  }

  unsigned int  RegPixelMap::GlobalToIndex(const unsigned int& plane,
                               const unsigned int& cell)
  {
    unsigned int view = plane % 2;

    unsigned int internalPlane =  plane - fBound.FirstPlane();
    unsigned int internalCell  =  cell - fBound.FirstCell(view);

    unsigned int index = internalPlane * fNCell + internalCell % fNCell;

    assert(index < fPE.size());

    return index;
  }

  unsigned int  RegPixelMap::LocalToIndex(const unsigned int& plane,
                               const unsigned int& cell) const
  {
    unsigned int index = plane * fNCell + cell % fNCell;

    assert(index < fPE.size());
    return index;
  }

  unsigned int  RegPixelMap::GlobalToIndexSingle(const unsigned int& plane,
                               const unsigned int& cell)
  {
    unsigned int view = plane % 2;

//    unsigned int internalPlane =  ((plane - fBound.FirstPlane())/2);//from standard Pixelmap
    unsigned int internalPlane =  ((plane - fBound.FirstPlane()));

    unsigned int internalCell  =  cell - fBound.FirstCell(view);

    unsigned int index = internalPlane * fNCell + internalCell % fNCell;

    assert(index < fPEX.size());

    return index;
  }



  void RegPixelMap::Print()
  {

    // Start by doing even planes
    for(unsigned int iCell = 0; iCell < fNCell; ++iCell)
    {
      for(unsigned int iPlane = 0; iPlane < fNPlane; iPlane += 2)
      {
        unsigned int index = LocalToIndex(iPlane, iCell);
        if( fPE[index] > 0)
        {
          std::cout << "*";
        }
        else
        {
          std::cout << " ";
        }

      }
      std::cout << std::endl;
    }
    // Then do odd planes
    for(unsigned int iCell = 0; iCell < fNCell; ++iCell)
    {
      for(unsigned int iPlane = 1; iPlane < fNPlane; iPlane += 2)
      {
        unsigned int index = LocalToIndex(iPlane, iCell);
        if( fPE[index] > 0)
        {
          std::cout << "*";
        }
        else
        {
          std::cout << " ";
        }

      }
      std::cout << std::endl;
    }

  }

  TH2F* RegPixelMap::ToTH2(TString name) const
  {

    // Create a histogram, use twice as many cells to distinguish views
    TH2F* hist = new TH2F(name, ";Plane;Cell", fNPlane, 0, fNPlane,
                          fNCell*2, 0, fNCell*2);
    // TH2F* hist = new TH2F("RegPixelMap", ";Plane;Cell", fNPlane, 0, fNPlane,
    //                                                 fNCell*2, 0, fNCell*2);

    for(unsigned int iPlane = 0; iPlane < fNPlane; ++iPlane)
    {
      for(unsigned int iCell = 0; iCell < fNCell; ++iCell)
      {
        // Add 1 to in each bin to skip underflow
        hist->SetBinContent(iPlane+1, iCell + fNCell*(iPlane%2) + 1,
                            fPE[LocalToIndex(iPlane, iCell)]);

      }
    }
    return hist;
  }

  TH2F* RegPixelMap::ToLabTH2(TString name) const
  {

    // Create a histogram, use twice as many cells to distinguish views
    TH2F* hist = new TH2F(name, ";Plane;Cell", fNPlane, 0, fNPlane,
                                                    fNCell*2, 0, fNCell*2);

    for(unsigned int iPlane = 0; iPlane < fNPlane; ++iPlane)
    {
      for(unsigned int iCell = 0; iCell < fNCell; ++iCell)
      {

        float pe  =  fPE[LocalToIndex(iPlane, iCell)];
        float lab = (float)fLab[LocalToIndex(iPlane, iCell)];
        // Draw label 0 as non-empty bin
        if ( pe > 0 && lab < 1 ) lab = 0.1;

        // Add 1 to in each bin to skip underflow
        hist->SetBinContent(iPlane+1, iCell + fNCell*(iPlane%2) + 1, lab);

      }
    }
    return hist;
  }

  TH2F* RegPixelMap::SingleViewToTH2(const unsigned int& view, TString name) const
  {

    TH2F* hist = new TH2F(name, ";Plane;Cell", (fNPlane/2)+(fNPlane%2),
                          0, (fNPlane/2)+(fNPlane%2), fNCell, 0, fNCell);

    for(unsigned int iPlane = 0; iPlane < ((fNPlane/2)+(fNPlane%2)); ++iPlane)
    {
      for(unsigned int iCell = 0; iCell < fNCell; ++iCell)
      {

        double pe = (view == 0) ? fPEX[LocalToIndex(iPlane, iCell)] : fPEY[LocalToIndex(iPlane, iCell)];
        
        // Add 1 to in each bin to skip underflow
        hist->SetBinContent(iPlane+1, iCell+1, pe);

      }
    }
    return hist;
  }

  TH2F* RegPixelMap::SingleViewToLabTH2(const unsigned int& view, TString name) const
  {

    TH2F* hist = new TH2F(name, ";Plane;Cell", (fNPlane/2)+(fNPlane%2),
                          0, (fNPlane/2)+(fNPlane%2), fNCell, 0, fNCell);

    for(unsigned int iPlane = 0; iPlane < ((fNPlane/2)+(fNPlane%2)); ++iPlane)
    {
      for(unsigned int iCell = 0; iCell < fNCell; ++iCell)
      {
        double pe  = (view == 0) ? fPEX[LocalToIndex(iPlane, iCell)] : fPEY[LocalToIndex(iPlane, iCell)];
        double lab = (view == 0) ? fLabX[LocalToIndex(iPlane, iCell)] : fLabY[LocalToIndex(iPlane, iCell)];

        // Draw label 0 as non-empty bin
        if ( pe > 0 && lab < 1 ) lab = 0.1;

        // Add 1 to in each bin to skip underflow
        hist->SetBinContent(iPlane+1, iCell+1, lab);

      }
    }
    return hist;
  }

  std::ostream& operator<<(std::ostream& os, const RegPixelMap& m)
  {
    os << "RegPixelMap with " << m.NPixel() << " pixels, "
                           << m.NPlane() << " planes"
                << " by "  << m.NCell()  << " cells" ;
    return os;
  }
}