FillEnergies.cxx

Go to the documentation of this file.

#include "CAFMaker/FillEnergies.h"

#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"

#include "RecoBase/Cluster.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Shower.h"

#include "LSTME/prediction/LSTMEnergy.h"
#include "NumuEnergy/NumuE.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "StandardRecord/StandardRecord.h"

namespace caf
{
  //......................................................................
  void FillNumuEnergyVars(const numue::NumuE& e,
                          caf::SRNumuEnergy& sre)
  {
    sre.E             =  e.E();
    sre.calccE        =  e.CalCCE();
    sre.trkqeE        =  e.TrkQEE();
    sre.trknonqeE     =  e.TrkNonQEE();
    sre.trkccE        =  e.TrkCCE();
    sre.angleE        =  e.AngleQEE();
    sre.angleerror    =  e.AngleQEError();
    sre.recomuonE     =  e.RecoMuonE();
    sre.ucrecomuonE1trk   =  e.UCMuonESingle();
    sre.ucrecomuonE   =  e.UCMuonENonSingle();
    sre.recotrkcchadE =  e.RecoTrkCCHadE();
    sre.hadcalE       =  e.HadCalE();
    sre.hadtrkE       =  e.HadTrkE();
    sre.ndtrklenact   =  e.NDTrkLenAct();
    sre.ndtrklencat   =  e.NDTrkLenCat();
    sre.ndtrkcalactE  =  e.NDTrkCalAct();
    sre.ndtrkcaltranE =  e.NDTrkCalTran();
    sre.ndtrkcalcatE  =  e.NDTrkCalCat();
    sre.ndhadcalactE  =  e.NDHadCalAct();
    sre.ndhadcaltranE =  e.NDHadCalTran();
    sre.ndhadcalcatE  =  e.NDHadCalCat();
    sre.ndhadtrkactE  =  e.NDHadTrkAct();
    sre.ndhadtrktranE =  e.NDHadTrkTran();
    sre.ndhadtrkcatE  =  e.NDHadTrkCat();
    sre.ndtrktranx    =  e.NDTrkTranX();
    sre.ndtrktrany    =  e.NDTrkTranY();

    FillHadClustVars(e.HadCluster(), sre.hadclust);

    sre.mc.truemuonE         =  e.MCTrueMuonE();
    sre.mc.truemuoncatcherE  =  e.MCTrueMuonCatcherE();
    sre.mc.truegoodmuon      =  e.MCGoodTrueMuon();
  }

  void FillNumuLSTMEnergyVars(const LSTME::LSTMEnergy& artE, float& lstmmuon, float& lstmnu)
  {
    lstmmuon = artE.primaryE;
    lstmnu   = artE.totalE;
  }

  //......................................................................
  void FillNueEnergyVars(const rb::Cluster& slice,const SRVertexBranch& srVertexBranch
                         , caf::SRNueEnergy& sre)
  {
      sre.rawgevslice = slice.TotalGeV();
//    sre.rawgevshw = shw.TotalGeV();
//    sre.rawgevhad = slice.TotalGeV()-shw.TotalGeV();

       if(srVertexBranch.elastic.IsValid == 0  ) return;
       if(srVertexBranch.elastic.fuzzyk.png.size() < 1) return;

       double CVNem_CalE = 0.0;
       double TotCalE = 0.0;
       for(unsigned int fuzzkProngIdx = 0; fuzzkProngIdx < srVertexBranch.elastic.fuzzyk.png.size(); ++fuzzkProngIdx){
        const SRFuzzyKProng &srPng = srVertexBranch.elastic.fuzzyk.png[fuzzkProngIdx];       
        double png_CalE = srVertexBranch.elastic.fuzzyk.png[fuzzkProngIdx].shwlid.calE;

        TotCalE += png_CalE;
        double emPID = srPng.cvnpart.electronid + srPng.cvnpart.pizeroid + srPng.cvnpart.photonid; 
        double hadPID = srPng.cvnpart.muonid + srPng.cvnpart.protonid + srPng.cvnpart.neutronid + srPng.cvnpart.pionid + srPng.cvnpart.otherid;
        if ( emPID < 0 ) continue;
        if ( emPID >= hadPID   ) CVNem_CalE += png_CalE;
        else continue;
        }//end iPng
        if( CVNem_CalE == 0.0) CVNem_CalE = srVertexBranch.elastic.fuzzyk.png[0].shwlid.calE;
        sre.rawgevshw = CVNem_CalE;

        if(srVertexBranch.elastic.fuzzyk.png2d.size() < 1) TotCalE +=0.0;
        else{
            for(unsigned int iPng = 0; iPng < srVertexBranch.elastic.fuzzyk.png2d.size(); ++iPng){
            TotCalE += srVertexBranch.elastic.fuzzyk.png2d[iPng].calE;
                }
        }
        TotCalE += srVertexBranch.elastic.fuzzyk.orphCalE;
        sre.rawgevhad = std::max(TotCalE - CVNem_CalE, 0.0);
  }

  //......................................................................
  void FillHadClustVars(const rb::Cluster& hadclust,
                        caf::SRHadClust& srhadclust)
  {
    static art::ServiceHandle<calib::Calibrator> cal;
    static art::ServiceHandle<geo::Geometry> geom;

    const double wx = (hadclust.NXCell() > 0) ? hadclust.W(hadclust.XCell(0).get()) : 0;
    const double wy = (hadclust.NYCell() > 0) ? hadclust.W(hadclust.YCell(0).get()) : 0;

    // Variables needed:
    // Holds contiguous planes info:
    unsigned int  cellsFromEdge = -1;

    // Hit counters:
    int nCalHit = 0;
    int nMipHit = 0;

    unsigned int minCellsFromEdge = 99999999;

    // Loop over hadclust in question here:
    for(unsigned int hitIdx = 0; hitIdx < hadclust.NCell(); ++hitIdx)
    {
      const art::Ptr<rb::CellHit>& chit = hadclust.Cell(hitIdx);
      // Get calbrated reco hit:
      const rb::RecoHit rhit(cal->MakeRecoHit(*chit,
                                              chit->View() == geo::kX ? wx : wy));

      const int planeNum = chit->Plane();
      cellsFromEdge = std::min((unsigned int)chit->Cell(), geom->Plane(planeNum)->Ncells() - 1 - chit->Cell());

      if(cellsFromEdge < minCellsFromEdge)
      {
        minCellsFromEdge = cellsFromEdge;
      }

      // Make sure things have been calibrated:
      if(!rhit.IsCalibrated())
      {
        //        MF_LOG_DEBUG("CAFMaker") << "Not calibrated?! "
        //        << chit->Plane() << " " << chit->Cell();
        continue;
      }

      // Count up both the total calibrated hits
      // and the MIP hits from 100-245 PECorr
      // Count up total hits:
      nCalHit++;
      // Count up MIP hits:
      if ( (rhit.PECorr()>100.0) && (rhit.PECorr()<245.0) ) nMipHit++;
    } // end of hadclust loop

    srhadclust.nhit           = hadclust.NCell();
    srhadclust.ncalhit        = nCalHit;
    srhadclust.nmiphit        = nMipHit;
    srhadclust.ncontplanes    = hadclust.MostContiguousPlanes(geo::kXorY);
    if (!(hadclust.NCell()==0)){

      srhadclust.firstplane     = hadclust.MinPlane();
      srhadclust.lastplane      = hadclust.MaxPlane();
      srhadclust.firstcell      = hadclust.MinCell(geo::kXorY);
      srhadclust.lastcell       = hadclust.MaxCell(geo::kXorY);
      srhadclust.ncellsfromedge = minCellsFromEdge;

      srhadclust.calE           = hadclust.CalorimetricEnergy();

      TVector3 boxmin, boxmax, meanpos;
      hadclust.MinMaxMeanXYZ(boxmin, boxmax, meanpos);
      srhadclust.boxmin = boxmin;
      srhadclust.boxmax = boxmax;
      srhadclust.meanpos = meanpos;

      // TBrowser can't draw DBL_MAX
      if (hadclust.NXCell()==0){
        srhadclust.boxmin.SetX(-5555);
        srhadclust.boxmax.SetX(5555);
      }
      if (hadclust.NYCell()==0){
        srhadclust.boxmin.SetY(-5555);
        srhadclust.boxmax.SetY(5555);
      }

    } // end of loop over variables to only fill for hadclust if it has hits
    else{
      srhadclust.firstplane     = 0;
      srhadclust.lastplane      = 0;
      srhadclust.firstcell      = 0;
      srhadclust.lastcell       = 0;
      srhadclust.ncellsfromedge = 0;

      srhadclust.calE           = -5;

      srhadclust.boxmin         = TVector3(-5555,-5555,-5555);
      srhadclust.boxmax         = TVector3(-5555,-5555,-5555);
      srhadclust.meanpos        = TVector3(-5555,-5555,-5555);
    }
  }

} // end namespace caf