showerreco::ShowerRecoAlg Class Reference

#include "ShowerRecoAlg.h"

Inheritance diagram for showerreco::ShowerRecoAlg:

Public Member Functions
	ShowerRecoAlg ()
	Default constructor. More...
virtual	~ShowerRecoAlg ()
	Default destructor. More...
virtual void	Reset ()
	Function to reset algorithm, to be called @ beginning of each event. More...
void	SetUseArea (bool on)
	Function to decide if to use Area or Pulse Amplitude for calculations. More...
void	setEcorrection (bool on)
	Function to decide if to use ModBox (True) or Birks (False) for calorimetry. More...
virtual void	AppendInputClusters (const std::vector< cluster::ClusterParamsAlg > &cpan_v)
	Setter for a matched combination of clusters. More...
std::vector< recob::Shower >	Reconstruct ()
	Execute reconstruction. More...
virtual void	Verbose (bool on=true)
	Verbosity switch. More...
void	CaloAlgo (::calo::CalorimetryAlg *alg)
	Calorimetry algorithm setter. More...

Protected Member Functions
virtual ::recob::Shower	RecoOneShower (const std::vector< ::showerreco::ShowerCluster_t > &)
	Function to reconstruct a shower. More...
virtual void	ProcessInputClusters ()
	Function to reorganize input cluster information. More...
virtual ::recob::Shower	RecoOneShower (const std::vector< showerreco::ShowerCluster_t > &clusters)=0
	Function to reconstruct one shower. More...

Protected Attributes
util::GeometryUtilities *	fGSer
bool	fVerbosity
	Verbosity flag. More...
::calo::CalorimetryAlg *	fCaloAlg
	Calorimetry algorithm. More...
std::vector< std::vector< showerreco::ShowerCluster_t > >	fInputClusters
	Input clusters. More...

Private Attributes
bool	_Ecorrection
	Boolean -> decide if to use energy correction or not. More...
double	fcalodEdxlength
double	fdEdxlength
bool	fUseArea

Detailed Description

User defined class ShowerRecoAlg ... these comments are used to generate doxygen documentation!

Definition at line 28 of file ShowerRecoAlg.h.

Constructor & Destructor Documentation

showerreco::ShowerRecoAlg::ShowerRecoAlg

(

)

Default constructor.

Definition at line 8 of file ShowerRecoAlg.cxx.

References _Ecorrection, fcalodEdxlength, fdEdxlength, fGSer, fUseArea, showerreco::ShowerRecoAlgBase::fVerbosity, and util::GeometryUtilities::GetME().

                               : ShowerRecoAlgBase(), fGSer(nullptr)
  {
    
    if(!fGSer) fGSer = (util::GeometryUtilities*)(util::GeometryUtilities::GetME());
    
    fcalodEdxlength=1000;
    fdEdxlength=2.4;
    fUseArea = true;
    fVerbosity = true;
    _Ecorrection = true;
  }

virtual showerreco::ShowerRecoAlg::~ShowerRecoAlg ( )

inlinevirtual

Default destructor.

Definition at line 36 of file ShowerRecoAlg.h.

{}

Member Function Documentation

void showerreco::ShowerRecoAlgBase::AppendInputClusters ( const std::vector< cluster::ClusterParamsAlg > &  cpan_v )

virtualinherited

Setter for a matched combination of clusters.

Definition at line 18 of file ShowerRecoAlgBase.cxx.

References showerreco::ShowerRecoAlgBase::fInputClusters.

Referenced by showerreco::ShowerRecoManager::Process().

  {
    std::vector< ::showerreco::ShowerCluster_t> clusters;
    clusters.reserve(cpan_v.size());

    for(auto const& cpan : cpan_v) {

      clusters.push_back( ::showerreco::ShowerCluster_t() );

      (*clusters.rbegin()).start_point = cpan.GetParams().start_point;
      (*clusters.rbegin()).end_point   = cpan.GetParams().end_point;
      (*clusters.rbegin()).angle_2d    = cpan.GetParams().angle_2d;
      (*clusters.rbegin()).plane_id    = cpan.Plane();
      (*clusters.rbegin()).hit_vector  = cpan.GetHitVector();

    }

    fInputClusters.push_back(clusters);

  }

void showerreco::ShowerRecoAlgBase::CaloAlgo ( ::calo::CalorimetryAlg *  alg )

inlineinherited

Calorimetry algorithm setter.

Definition at line 64 of file ShowerRecoAlgBase.h.

Referenced by ShowerReco3D::ShowerReco3D().

{ fCaloAlg = alg; }

virtual void showerreco::ShowerRecoAlgBase::ProcessInputClusters ( )

inlineprotectedvirtualinherited

Function to reorganize input cluster information.

Definition at line 69 of file ShowerRecoAlgBase.h.

Referenced by showerreco::ShowerRecoAlgBase::Reconstruct().

    { return; }

std::vector<::recob::Shower > showerreco::ShowerRecoAlgBase::Reconstruct ( )

inherited

Execute reconstruction.

Definition at line 39 of file ShowerRecoAlgBase.cxx.

References showerreco::ShowerRecoAlgBase::fCaloAlg, showerreco::ShowerRecoAlgBase::fInputClusters, showerreco::ShowerRecoAlgBase::ProcessInputClusters(), and showerreco::ShowerRecoAlgBase::RecoOneShower().

Referenced by showerreco::ShowerRecoManager::Process().

  {

    if(!fCaloAlg) throw ShowerRecoException("Calorimetry algorithm must be provided!");

    ProcessInputClusters();

    std::vector< ::recob::Shower> output;

    output.reserve(fInputClusters.size());

    for(auto const& clusters : fInputClusters)

      output.push_back( RecoOneShower(clusters) );

    return output;

  }

recob::Shower showerreco::ShowerRecoAlg::RecoOneShower ( const std::vector< ::showerreco::ShowerCluster_t > &  clusters )

protected

Function to reconstruct a shower.

third loop to get only points inside of 1RMS of value.

Definition at line 21 of file ShowerRecoAlg.cxx.

References _Ecorrection, calo::CalorimetryAlg::dEdx_AMP(), calo::CalorimetryAlg::dEdx_AREA(), showerreco::energy::DEFAULT_ECorr, calo::CalorimetryAlg::ElectronsFromADCArea(), calo::CalorimetryAlg::ElectronsFromADCPeak(), showerreco::ShowerRecoAlgBase::fCaloAlg, fcalodEdxlength, fdEdxlength, fGSer, fUseArea, showerreco::ShowerRecoAlgBase::fVerbosity, util::GeometryUtilities::Get2DDistance(), util::GeometryUtilities::Get3DaxisN(), util::GeometryUtilities::GetDirectionCosines(), util::GeometryUtilities::GetPointOnLine(), util::kGeVToElectrons, calo::CalorimetryAlg::LifetimeCorrection(), util::GeometryUtilities::Nplanes(), util::GeometryUtilities::PitchInView(), recob::Shower::set_dedx(), recob::Shower::set_direction(), recob::Shower::set_length(), recob::Shower::set_start_point(), recob::Shower::set_total_best_plane(), recob::Shower::set_total_energy(), recob::Shower::set_total_MIPenergy(), util::GeometryUtilities::TimeToCm(), w, and util::GeometryUtilities::WireToCm().

Referenced by setEcorrection().

  {
    
    ::recob::Shower result;
    //
    // Reconstruct and store
    //
    std::vector < util::PxPoint > fStartPoint;    // for each plane
    std::vector < util::PxPoint > fEndPoint;    // for each plane
    std::vector < double > fOmega2D;    // for each plane
    
    std::vector < double > fEnergy   (fGSer->Nplanes(),-1);    // for each plane
    std::vector < double > fMIPEnergy(fGSer->Nplanes(),-1);    // for each plane
    std::vector < double > fdEdx     (fGSer->Nplanes(),-1);      
    std::vector <unsigned char> fPlaneID;
    
    // First Get Start Points
    for(auto const & cl : clusters)
      {
        fStartPoint.push_back(cl.start_point);    // for each plane
        fEndPoint.push_back(cl.end_point);    // for each plane
        fOmega2D.push_back(cl.angle_2d);
        fPlaneID.push_back(cl.plane_id);
        if(fVerbosity) {
          std::cout << " planes : " << cl.plane_id
                    << " " << cl.start_point.t 
                    << " " << cl.start_point.w 
                    << " " << cl.end_point.t 
                    << " " << cl.end_point.w 
                    <<" angle2d "<<  cl.angle_2d  
                    << std::endl;   
        }
      }
    
    
    //decide best two planes. for now, using length in wires - flatter is better.
    int index_to_use[2]={0,1}; 
    int best_plane=-1;
    //int good_plane=-1;
    double best_length=0;
    double good_length=0;
    for (size_t ip=0;ip<fPlaneID.size();ip++)
      {
        double dist = fabs( fEndPoint[ip].w - fStartPoint[ip].w );
        if(dist > best_length )
          {
            good_length = best_length;
            //good_plane  = best_plane;
            index_to_use[1] = index_to_use[0];

            best_length = dist;
            best_plane = fPlaneID.at(ip);
            index_to_use[0] = ip;
          }
        else if(dist > good_length)
          {
            good_length = dist;
            //good_plane  = fPlaneID.at(ip);
            index_to_use[1] = ip;
          }
      }

    // Second Calculate 3D angle and effective pitch and start point 
    double xphi=0,xtheta=0;
    
    
    fGSer->Get3DaxisN(fPlaneID[index_to_use[0]],
                      fPlaneID[index_to_use[1]],
                      fOmega2D[index_to_use[0]]*TMath::Pi()/180.,
                      fOmega2D[index_to_use[1]]*TMath::Pi()/180.,
                      xphi,
                      xtheta);
    
    if(fVerbosity)
      std::cout << " new angles: " << xphi << " " << xtheta << std::endl; 
    
    
    
    double xyz[3];
    // calculate start point here?
    fGSer-> GetXYZ(&(fStartPoint[index_to_use[0]]),&(fStartPoint[index_to_use[1]]),xyz);
    
    if(fVerbosity)
      std::cout << " XYZ:  " << xyz[0] << " " << xyz[1] << " " << xyz[2] << std::endl;
    
    // Third calculate dE/dx and total energy for all planes, because why not?
    //for(auto const &clustit : clusters)
    for(size_t cl_index=0; cl_index < fPlaneID.size(); ++cl_index) 
      {
        int plane = fPlaneID.at(cl_index);
        double newpitch=fGSer->PitchInView(plane,xphi,xtheta);
        if(plane == best_plane) best_length *= newpitch / fGSer->WireToCm();

        if(fVerbosity)
          std::cout << std::endl << " Plane: " << plane << std::endl;
        
        double totEnergy=0;
        double totLowEnergy=0;
        double totHighEnergy=0;
        double totMIPEnergy=0;
        int direction=-1;
        //double RMS_dedx=0;
        double dEdx_av=0,dedx_final=0;
        int npoints_first=0, npoints_sec=0;
        
        //override direction if phi (XZ angle) is less than 90 degrees 
        // this is shady, and needs to be rethought.
        if(fabs(fOmega2D.at(cl_index)) < 90)
          direction=1;
        
        auto const& hitlist = clusters.at(cl_index).hit_vector;
        std::vector<util::PxHit> local_hitlist;
        local_hitlist.reserve(hitlist.size());

        for(const auto& theHit : hitlist){
          
          double dEdx_new=0;
          double hitElectrons = 0;
          //double Bcorr_half;
          //double dEdx_sub;
          // double dEdx_MIP;
          
          //    dEdx_new = fCaloAlg->dEdx_AREA(theHit, newpitch );
          //Bcorr_half = 2.*fCaloAlg->dEdx_AREA(theHit->Charge()/2.,theHit->PeakTime(), newpitch, plane); ; 
          //dEdx_sub = fCaloAlg->dEdx_AREA(theHit->Charge()-PION_CORR,theHit->PeakTime(), newpitch, plane); ; 
          // dEdx_MIP = fCaloAlg->dEdx_AREA_forceMIP(theHit, newpitch ); 
          if(!fUseArea)
            {
              dEdx_new = fCaloAlg->dEdx_AMP(theHit.peak / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
              hitElectrons = fCaloAlg->ElectronsFromADCPeak(theHit.peak, plane);
            }
          else
            {
              dEdx_new = fCaloAlg->dEdx_AREA(theHit.charge / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
              hitElectrons = fCaloAlg->ElectronsFromADCArea(theHit.charge, plane);
            }

          hitElectrons *= fCaloAlg->LifetimeCorrection(theHit.t / fGSer->TimeToCm());
          
          totEnergy += hitElectrons * 1.e3 / (::util::kGeVToElectrons);

          //totNewCnrg+=dEdx_MIP;
        //  if(dEdx_new < 3.5 && dEdx_new >0 )
        //    {
        //      totLowEnergy +=dEdx_new*newpitch;
        //      totMIPEnergy += dEdx_new*newpitch;
        //    }
        //  else
         //   {
         //     totHighEnergy += dEdx_new*newpitch;
              int multiplier=1;
              if(plane < 2) multiplier =2 ;
              if(!fUseArea){
                totMIPEnergy += theHit.peak *  0.0061*multiplier;
              }
              else{
                totMIPEnergy += theHit.charge *  0.00336*multiplier;
              }
          //  }
          
           if(fVerbosity && dEdx_new >1.9 && dEdx_new <2.1)
            std::cout << "dEdx_new " << dEdx_new << " " <<dEdx_new/theHit.charge*newpitch << " "<< theHit.charge *0.0033*multiplier/newpitch << std::endl;
          
          util::PxPoint OnlinePoint; 
          // calculate the wire,time coordinates of the hit projection on to the 2D shower axis
          fGSer->GetPointOnLine(fOmega2D.at(cl_index),
                                &(fStartPoint.at(cl_index)),
                                &theHit,
                                OnlinePoint);
          
          double ortdist=fGSer->Get2DDistance(&OnlinePoint,&theHit);
          double linedist=fGSer->Get2DDistance(&OnlinePoint,&(fStartPoint.at(cl_index)));
          
          
          //calculate the distance from the vertex using the effective pitch metric 
          double wdist=((theHit.w-fStartPoint.at(cl_index).w)*newpitch)*direction;  //wdist is always positive
          
          if(fVerbosity && dEdx_new < 3.5 )
            std::cout << " CALORIMETRY:" 
                      << " Pitch " <<newpitch 
                      << " dist: " << wdist 
                      << " dE/dx: " << dEdx_new << "MeV/cm "  
                      << " average: " <<  totEnergy  
                      << "hit: wire, time " << theHit.w/fGSer->WireToCm() << " " << theHit.t/fGSer->TimeToCm() 
                      << "total energy" << totEnergy 
                      << std::endl;
          
          //   if( (fabs(wdist)<fcalodEdxlength)&&(fabs(wdist)>0.2)){  
          if( (wdist<fcalodEdxlength)&&(wdist>0.2)){    
            
            //vdEdx.push_back(dEdx_new);
            //vresRange.push_back(fabs(wdist));
            //vdQdx.push_back(theHit->Charge(true)/newpitch);
            //Trk_Length=wdist;
            
            //fTrkPitchC=fNPitch[set][plane];
            //Kin_En+=dEdx_new*newpitch;
            //npoints_calo++;
            //sum+=dEdx_new;
            
            
            
            //fDistribChargeADC[set].push_back(ch_adc);  //vector with the first De/Dx points
            //     fDistribChargeMeV[set].push_back(dEdx_new);  //vector with the first De/Dx points
            //     fDistribHalfChargeMeV[set].push_back(Bcorr_half);
            //     fDistribChargeposition[set].push_back(wdist);  //vector with the first De/Dx points' positions 
            //     fDistribChargeMeVMIPsub[set].push_back(dEdx_sub);
            //  
            
            
            //first pass at average dE/dx
            if(wdist<fdEdxlength 
               //take no hits before vertex (depending on direction)
               && ((direction==1 && theHit.w > fStartPoint.at(cl_index).w) || (direction==-1 && theHit.w < fStartPoint.at(cl_index).w) )
               && ortdist<3.5 
               && linedist < fdEdxlength ){
              
              dEdx_av+= dEdx_new;
              local_hitlist.push_back(theHit);
              
              npoints_first++;
              if(fVerbosity)
                std::cout << " CALORIMETRY:" << " Pitch " <<newpitch 
                          << " dist: " << wdist 
                          << " dE/dx: " << dEdx_new << "MeV/cm "  
                          << " average: " <<  dEdx_av/npoints_first  
                          << " hit: wire, time " << theHit.w << " " << theHit.t 
                          << " line,ort " << linedist << " " << ortdist
                          << " direction " << direction 
                          << std::endl;
            }
            
          }//end inside range if statement
          
        }// end first loop on hits.
        
        
        
        double mevav2cm=0;
        double fRMS_corr=0;
        
        //calculate average dE/dx
        if(npoints_first>0)     {
          mevav2cm=dEdx_av/npoints_first;
          
        }
        
        //loop only on subset of hits
        for(auto const & theHit : local_hitlist){
          double dEdx=0;
          if(fUseArea)
            { 
              dEdx = fCaloAlg->dEdx_AREA(theHit.charge / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
            }
          else  //this will hopefully go away, once all of the calibration factors are calculated.
            {
              dEdx = fCaloAlg->dEdx_AMP(theHit.peak / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
            }
          //fDistribAfterMin[set].push_back(MinBefore);
          //fDistribBeforeMin[set].push_back(MinAfter);
          //  auto position = hitIter - local_hitlist.begin() ; 
          
          fRMS_corr+=(dEdx-mevav2cm)*(dEdx-mevav2cm); 
          
        }
        
        
        
        if(npoints_first>0)     
          {
            fRMS_corr=TMath::Sqrt(fRMS_corr/npoints_first);
          }
        
        //loop only on subset of hits

        for(auto const & theHit : local_hitlist){
          double dEdx=0;
          if(fUseArea)
            { 
              dEdx = fCaloAlg->dEdx_AREA(theHit.charge / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
            }
          else  //this will hopefully go away, once all of the calibration factors are calculated.
            {
              dEdx = fCaloAlg->dEdx_AMP(theHit.peak / newpitch, theHit.t / fGSer->TimeToCm(), theHit.plane);
            }
          //fDistribAfterMin[set].push_back(MinBefore);
          //fDistribBeforeMin[set].push_back(MinAfter);
          
          if( ((dEdx > (mevav2cm-fRMS_corr) ) && (dEdx < (mevav2cm+fRMS_corr) )) 
              || 
              (newpitch > 0.3*fdEdxlength ) 
              ) {
            dedx_final+= dEdx; 
            npoints_sec++;
          }
        }
        
        if(npoints_sec>0){
          dedx_final/=npoints_sec;
        }
        
        if(fVerbosity){ 
          std::cout << " total ENERGY, birks: " << totEnergy << " MeV "  
                    << " |average:  " << dedx_final  << std::endl
                    << " Energy:  lo:" << totLowEnergy << " hi: " << totHighEnergy 
                    << " MIP corrected " << totMIPEnergy << std::endl;
        }
        // if Energy correction factor to be used
        // then get it from ShowerCalo.h
        if (_Ecorrection) { totEnergy *= showerreco::energy::DEFAULT_ECorr; }
        fEnergy[plane]    = totEnergy;    // for each plane
        fMIPEnergy[plane] = totMIPEnergy;
        fdEdx[plane]      = dedx_final; 
    
        // break;
      } // end loop on clusters

    result.set_total_MIPenergy(fMIPEnergy);
    result.set_total_best_plane(best_plane);
    result.set_length(best_length);
    result.set_total_energy(fEnergy);
    //result.set_total_energy_err  (std::vector< Double_t > q)            { fSigmaTotalEnergy = q;        }
    
    double dirs[3]={0};
    fGSer->GetDirectionCosines(xphi,xtheta,dirs);
    TVector3 vdirs(dirs[0],dirs[1],dirs[2]);
    
    TVector3 vxyz(xyz[0], xyz[1], xyz[2]);
    
    result.set_direction(vdirs); 
    //result.set_direction_err (TVector3 dir_e)      { fSigmaDCosStart = dir_e; }
    result.set_start_point(vxyz); 
    //result.set_start_point_err (TVector3 xyz_e)      { fSigmaXYZstart = xyz_e; }
    result.set_dedx  (fdEdx);
    //result.set_dedx_err  (std::vector< Double_t > q)            { fSigmadEdx = q;        }

    // done
    return result;
    }

virtual ::recob::Shower showerreco::ShowerRecoAlgBase::RecoOneShower ( const std::vector< showerreco::ShowerCluster_t > &  clusters )

protectedpure virtualinherited

Function to reconstruct one shower.

Referenced by showerreco::ShowerRecoAlgBase::Reconstruct().

virtual void showerreco::ShowerRecoAlg::Reset ( )

inlinevirtual

Function to reset algorithm, to be called @ beginning of each event.

Implements showerreco::ShowerRecoAlgBase.

Definition at line 39 of file ShowerRecoAlg.h.

References showerreco::ShowerRecoAlgBase::Reset().

{ ShowerRecoAlgBase::Reset(); }

void showerreco::ShowerRecoAlg::setEcorrection ( bool  on )

inline

Function to decide if to use ModBox (True) or Birks (False) for calorimetry.

Function to set whether to use E correction

Definition at line 48 of file ShowerRecoAlg.h.

References _Ecorrection, and RecoOneShower().

Referenced by ShowerReco3D::ShowerReco3D().

{ _Ecorrection = on; }

void showerreco::ShowerRecoAlg::SetUseArea ( bool  on )

inline

Function to decide if to use Area or Pulse Amplitude for calculations.

Definition at line 42 of file ShowerRecoAlg.h.

References fUseArea.

Referenced by ShowerReco3D::ShowerReco3D().

{ fUseArea = on; }

virtual void showerreco::ShowerRecoAlgBase::Verbose ( bool  on = true )

inlinevirtualinherited

Verbosity switch.

Definition at line 61 of file ShowerRecoAlgBase.h.

Referenced by ShowerReco3D::ShowerReco3D().

{ fVerbosity=on; }

Member Data Documentation

bool showerreco::ShowerRecoAlg::_Ecorrection

private

Boolean -> decide if to use energy correction or not.

Definition at line 62 of file ShowerRecoAlg.h.

Referenced by RecoOneShower(), setEcorrection(), and ShowerRecoAlg().

::calo::CalorimetryAlg* showerreco::ShowerRecoAlgBase::fCaloAlg

protectedinherited

Calorimetry algorithm.

Definition at line 81 of file ShowerRecoAlgBase.h.

Referenced by showerreco::ShowerRecoAlgBase::Reconstruct(), and RecoOneShower().

double showerreco::ShowerRecoAlg::fcalodEdxlength

private

Definition at line 64 of file ShowerRecoAlg.h.

Referenced by RecoOneShower(), and ShowerRecoAlg().

double showerreco::ShowerRecoAlg::fdEdxlength

private

Definition at line 65 of file ShowerRecoAlg.h.

Referenced by RecoOneShower(), and ShowerRecoAlg().

util::GeometryUtilities* showerreco::ShowerRecoAlg::fGSer

protected

Definition at line 57 of file ShowerRecoAlg.h.

Referenced by RecoOneShower(), and ShowerRecoAlg().

std::vector<std::vector<showerreco::ShowerCluster_t> > showerreco::ShowerRecoAlgBase::fInputClusters

protectedinherited

Input clusters.

Definition at line 84 of file ShowerRecoAlgBase.h.

Referenced by showerreco::ShowerRecoAlgBase::AppendInputClusters(), showerreco::ShowerRecoAlgBase::Reconstruct(), and showerreco::ShowerRecoAlgBase::Reset().

bool showerreco::ShowerRecoAlg::fUseArea

private

Definition at line 66 of file ShowerRecoAlg.h.

Referenced by RecoOneShower(), SetUseArea(), and ShowerRecoAlg().

bool showerreco::ShowerRecoAlgBase::fVerbosity

protectedinherited

Verbosity flag.

Definition at line 78 of file ShowerRecoAlgBase.h.

Referenced by RecoOneShower(), ShowerRecoAlg(), and showerreco::ShowerRecoAlgBase::ShowerRecoAlgBase().

---

The documentation for this class was generated from the following files:

• larreco/v06_64_02/source/larreco/ShowerFinder/ShowerReco3D/ShowerRecoAlg.h
• larreco/v06_64_02/source/larreco/ShowerFinder/ShowerReco3D/ShowerRecoAlg.cxx