DeconvolveAlg.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file    DeconvolveAlg.cxx
/// \brief   Separates the hadronic energy and electromagnetic energy in 
///          a cell
///
/// \author  Kanika Sachdev (ksachdev@physics.umn.edu)
////////////////////////////////////////////////////////////////////////

#include "messagefacility/MessageLogger/MessageLogger.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"

#include "ShowerLID/DeconvolveAlg.h"

#include "Calibrator/Calibrator.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "Utilities/func/MathUtil.h"

#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TMatrixD.h"

#include <string.h>

namespace slid
{

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////

  DeconvolveAlg::DeconvolveAlg(fhicl::ParameterSet const& pset)
  {
    this->reconfigure(pset);
    mf::LogInfo("DeconvolveAlg") << " DeconvolveAlg::DeconvolveAlg()\n";
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////

  DeconvolveAlg::~DeconvolveAlg()
  {

  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void DeconvolveAlg::reconfigure(const fhicl::ParameterSet& pset)
  {
    fExpConst = pset.get< float >("ExpConst");
  }

  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void DeconvolveAlg::GetWeights( std::vector<rb::Prong> & prongs,
                                  const std::vector< float > & totpe)
  {
    std::vector<rb::WeightedHit>  wHits;
    const int nPro = prongs.size();
    
    // This map holds info about the cells that are shared by 
    // multiple prongs. The key is the channel, and the value 
    // is another map. In the inner map, the key is the index 
    // of the prong and the value is the index of the cell in
    // that prong. Therefore, proIdx[chan][prgIdx] = [cellIdx]
    std::map< uint32_t , std::map< int,  int > > proIdx ;

    for( int iPro = 0; iPro < nPro; ++iPro){

      // loop over cells in ith prong
      for( unsigned int iCell = 0; iCell < prongs[iPro].NCell(); ++iCell ){
        const art::Ptr<rb::CellHit>& ihit = prongs[iPro].Cell(iCell);

        // In the following loop, we will find which other
        // prongs the ihits[iCell] belongs to.

        for( int jPro = iPro+1; jPro < nPro; ++jPro){

          // loop over cells in jth prong
          for( unsigned int jCell = 0 ; jCell < prongs[jPro].NCell(); ++jCell){
            const art::Ptr<rb::CellHit>& jhit = prongs[jPro].Cell(jCell);

            if( *jhit == *ihit){
              proIdx[ ihit->Channel() ][ iPro ] = iCell;
              proIdx[ ihit->Channel() ][ jPro ] = jCell;
            }

          }// end loop over jprong cells
        }// end j loop over prongs 
      }// end loop over iprong cells
    }// end i loop over prongs
    
    // We now have all the info we need to set weights:
    for( auto & cmap : proIdx){
      double weights[ cmap.second.size()];
      double total = 0;
      int    count = 0;
      for( auto & pmap : cmap.second){

        rb::Prong iProng  = prongs[ pmap.first ];
        int iCell    = pmap.second;
        double dist  = DistanceToCore( iProng.Cell( iCell), iProng );
        if (dist == 9999)
          weights[count] = -1;
        else{
            
          // Shower lateral profile constant from a fit to 
          // (Transverse dE/dx)/(shower Energy) as a function 
          // of distance to shower core for Monte Carlo electron.

          weights[count] = totpe[ pmap.first] * exp( -fExpConst*dist);
          total  += weights[count];
        }
        count = count +1;
      }// end loop over proIdx.second
      
      count=0;
      for( auto & pmap : cmap.second){

        if( weights[count] >= 0 )
          prongs[pmap.first].SetWeight( pmap.second, weights[count]/total);
        else
          prongs[pmap.first].SetWeight( pmap.second, 0.0 );
        count = count +1;
      }
    }// end loop over proIdx
    
    return;
    
  }// end of GetWeights


  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                             
  
  float DeconvolveAlg::DistanceToCore( art::Ptr<rb::CellHit> cHit,
                                       rb::Prong &  prong)
  {

    double dist = 9999, di;

    // If the prong is 2D, reclustering a cell from the other view
    // makes no sense. Bail.
    if( prong.Is2D() && ( prong.View() != cHit->View() ) )
      return dist;

    TVector3 start  = prong.Start();
    TVector3 dir    = prong.Dir();

    TVector3 xyz;
    fGeom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
    const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);

    // geo::Closest approach takes in 4 3D points-                                                                      
    // first two are start and stop points of a                                                                       
    // line and the nxt two, of the other.                                                                           
    // It returns the square of the distance of                                                                     
    // closest approach between two line.                                                                             
    
    dist = geo::ClosestApproach(start, start+dir,
                                xyz, xyz+dxyz, &di);
    
    return sqrt(dist);
    
  } // end of DistanceToCore                                                                  


  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                             



  /*********************************************************************************************/
  ////////////////////////////////////////////////////////////////////////                                             

}// End of namespace jmprong