SpaceChargeProtoDUNEdp.cxx File Reference

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include "math.h"
#include "stdio.h"
#include "dunesim/SpaceCharge/SpaceChargeProtoDUNEdp.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "larcore/CoreUtils/ServiceUtil.h"
#include "larcore/Geometry/Geometry.h"
#include "cetlib_except/exception.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "TFile.h"
#include "TH3.h"
#include "TTree.h"
#include "TLeaf.h"

Go to the source code of this file.

Functions
bool	isWithinHist (geo::Point_t const &point, TH3F *hist)
bool	isWithinHistOuter (geo::Point_t const &point, TH3F *hist)
double	ExtrapolateAtEdge (TH3F *hist, double x, double y, double z)
bool	IsBorder (TAxis *ax, double v)
double	GetSymmetric (double x, double xc)

Function Documentation

double ExtrapolateAtEdge	(	TH3F *	hist,
		double	x,
		double	y,
		double	z
	)

Definition at line 318 of file SpaceChargeProtoDUNEdp.cxx.

                                                                  {
  /*
     Want to get the value at x, close to the wall
    |----+----+--+--+---->
         x    c  a  b
     c is the center of the voxel (interpolation limit)
     b is the symetric of x wrt to c
     a is between c and b
     ->make linear interpolation from B, A to X

    */

     int binx = hist->GetXaxis()->FindBin(x);
     int biny = hist->GetYaxis()->FindBin(y);
     int binz = hist->GetZaxis()->FindBin(z);


     //get coordinate of the center of the voxel
     double xcenter = IsBorder(hist->GetXaxis(),x) ? hist->GetXaxis()->GetBinCenter(binx) : x;
     double ycenter = IsBorder(hist->GetYaxis(),y) ? hist->GetYaxis()->GetBinCenter(biny) : y;
     double zcenter = IsBorder(hist->GetZaxis(),z) ? hist->GetZaxis()->GetBinCenter(binz) : z;

     int Nx = hist->GetXaxis()->GetNbins();
     int Ny = hist->GetYaxis()->GetNbins();
     int Nz = hist->GetZaxis()->GetNbins();

     /*     if (Nz==1){ // when the GAr map is used (only one Z bin at the moment)
        return InterpolateAtEdge(x,y,z);// LZ
        }*/

      //special case when the starting point is at the center of the voxel
     if(x==xcenter && IsBorder(hist->GetXaxis(),x)){
         int binx2 = (binx==Nx) ? binx-1 : 2;
         xcenter =  0.5*(hist->GetXaxis()->GetBinCenter(binx2)+x);
        }
     if(y==ycenter&& IsBorder(hist->GetYaxis(),y)){
       int biny2 = (biny==Ny) ? biny-1 : 2;
       ycenter =  0.5*(hist->GetYaxis()->GetBinCenter(biny2)+y);
     }
     if(z==zcenter&& IsBorder(hist->GetZaxis(), z)){
            int binz2 = (binz==Nz) ? binz-1 : 2;
         zcenter =  0.5*(hist->GetZaxis()->GetBinCenter(binz2)+z);
       }


      //Get the symmetric of edgy point wrt to center
     double xb = IsBorder(hist->GetXaxis(),x) ? GetSymmetric(x,xcenter) : x;
     double yb = IsBorder(hist->GetYaxis(),y) ? GetSymmetric(y,ycenter) : y;
     double zb = IsBorder(hist->GetZaxis(),z) ? GetSymmetric(z,zcenter) : z;

       //Get middle point of [symmetric, center]
     double xa = IsBorder(hist->GetXaxis(),x) ?  0.5*(xcenter+xb) : x;
     double ya = IsBorder(hist->GetYaxis(),y) ?  0.5*(ycenter+yb) : y;
     double za = IsBorder(hist->GetZaxis(),z) ?  0.5*(zcenter+zb) : z;

     double valb = hist->Interpolate(xb, yb, zb); //symmetric point value
     double vala =  hist->Interpolate(xa, ya, za);
     double val = 0.;

     double Nedges = 0.;
     if(IsBorder(hist->GetXaxis(),x)) Nedges = Nedges+1.;
     if(IsBorder(hist->GetYaxis(),y)) Nedges = Nedges+1.;
     if(IsBorder(hist->GetZaxis(),z)) Nedges = Nedges+1.;



     double delta = vala - valb;//

     double dx = IsBorder(hist->GetXaxis(),x) ? xa-xb : 1;
     double dy = IsBorder(hist->GetYaxis(),y) ? ya-yb : 1;
     double dz = IsBorder(hist->GetZaxis(),z) ? za-zb : 1;

     val = x*delta/dx + (xa*valb - xb*vala)/dx
         +  y*delta/dy + (ya*valb - yb*vala)/dy
         +  z*delta/dz + (za*valb - zb*vala)/dz;


     val = val/Nedges;
     return val;

  }

double GetSymmetric	(	double	x,
		double	xc
	)

Definition at line 314 of file SpaceChargeProtoDUNEdp.cxx.

                                        {
     return xc + (xc-x);
   }

bool IsBorder	(	TAxis *	ax,
		double	v
	)

Definition at line 302 of file SpaceChargeProtoDUNEdp.cxx.

                                   {

       /*
              In the external half of the first/last bin or outside the range
               */
    if( ax->GetNbins()==1) return true;
    if( v<= ax->GetBinCenter(1) ) return true;
    //a->GetNbins() or GetNbins()-1..?
    if( v>= ax->GetBinCenter(ax->GetNbins()) ) return true;

    return false;
  }

bool isWithinHist	(	geo::Point_t const &	point,
		TH3F *	hist
	)

Definition at line 414 of file SpaceChargeProtoDUNEdp.cxx.

{
  if(point.X()> hist->GetXaxis()->GetBinCenter(1)  &&
     point.X()< hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) &&
     point.Y()> hist->GetYaxis()->GetBinCenter(1) &&
     point.Y()< hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) &&
     point.Z()> hist->GetZaxis()->GetBinCenter(1) &&
     point.Z()< hist->GetZaxis()->GetBinCenter(hist->GetZaxis()->GetNbins()) ){

    return true;}
  return false;

}

bool isWithinHistOuter	(	geo::Point_t const &	point,
		TH3F *	hist
	)

Definition at line 400 of file SpaceChargeProtoDUNEdp.cxx.

{
  if(point.X()> hist->GetXaxis()->GetBinLowEdge(0) &&
     point.X()< hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins()+1) &&
     point.Y()> hist->GetYaxis()->GetBinLowEdge(0) &&
     point.Y()< hist->GetYaxis()->GetBinLowEdge(hist->GetYaxis()->GetNbins()+1) &&
     point.Z()> hist->GetZaxis()->GetBinLowEdge(0) &&
     point.Z()< hist->GetZaxis()->GetBinLowEdge(hist->GetZaxis()->GetNbins()+1)  ){

    return true;}
  return false;

}