#include <PndRadMapBoxMesh.h>

Public Member Functions
	PndRadMapBoxMesh ()

	PndRadMapBoxMesh (PndRadMapBoxMesh &m)

	PndRadMapBoxMesh (const char *Name, int Xbins, Double_t Xlow, Double_t Xhigh, int Ybins, Double_t Ylow, Double_t Yhigh, int Zbins, Double_t Zlow, Double_t Zhigh)

	PndRadMapBoxMesh (const char *name, int xbins, Double_t xlow, Double_t xhigh, int ybins, Double_t ylow, Double_t yhigh, Double_t zlow, Double_t zhigh, orientation plane=ZX, quantity Quantity=Edep)

	~PndRadMapBoxMesh ()

void	SetFilter (const char *filter)

void	SetQuantity (quantity Quantity=Edep)

void	SetOrientation (orientation plane, Double_t rotate=99999, axis Ax=Xx)

void	SetOrientation (Double_t rotate=99999, axis Ax=Xx)

void	SetVerbosityLevel (int verbose=0)

void	Fill (FairRadMapPoint *p)

void	Transform (Double_t X, Double_t Y, Double_t Z)

void	Transform (Double_t X, Double_t Y, Double_t Z, Double_t &X0, Double_t &Y0, Double_t &Z0)

void	Transform (TVector3 InV, TVector3 &OutV)

void	Scale (Double_t sca)

void	Save (TFile *fout)

void	Save ()

TH2D *	GetHisto ()

PndRadMapPlane *	GetPlane ()

Double_t	CalcFluence (FairRadMapPoint *p)

Protected Member Functions
bool	IsInside (Double_t X, Double_t Y, Double_t Z)

bool	IsInside (FairRadMapPoint *p)

bool	IsInside ()

void	Fill (Int_t gBin, Double_t val)

void	Fill (Double_t X, Double_t Y, Double_t Z, Double_t we=1)

void	makeHisto (const char *Orient, Double_t rotate, int Hbins, Double_t Hlow, Double_t Hhigh, int Vbins, Double_t Vlow, Double_t Vhigh, Double_t dlow)

Protected Attributes
TH2D *	_MeshHisto

TH2I *	_StatHisto

TH2D *	_EnergyHisto

TString	_Name

orientation	_orientation

Double_t	_rotate

axis	_axis

Double_t	_volume

quantity	_quantity

int	_Xbins

Double_t	_Xlow

Double_t	_Xhigh

int	_Ybins

Double_t	_Ylow

Double_t	_Yhigh

int	_Zbins

Double_t	_Zlow

Double_t	_Zhigh

Double_t	_X

Double_t	_Y

Double_t	_Z

Double_t	_tX

Double_t	_tY

Double_t	_tZ

TFormula	_filter

int	_verbose

PndRadMapPlane *	_plane

TVector3	InterSection

Private Attributes
bool	_isSurfaceQuantity

TDatabasePDG *	pdg

TParticlePDG *	pdgpart

Detailed Description

Definition at line 50 of file PndRadMapBoxMesh.h.

Constructor & Destructor Documentation

PndRadMapBoxMesh::PndRadMapBoxMesh ( )

PndRadMapBoxMesh::PndRadMapBoxMesh ( PndRadMapBoxMesh & m )

Definition at line 87 of file PndRadMapBoxMesh.cxx.

References _isSurfaceQuantity, _MeshHisto, _Name, _orientation, _plane, _rotate, _StatHisto, _volume, and pdg.

                                                      {
   _MeshHisto = (TH2D*)m._MeshHisto->Clone();
   _StatHisto = (TH2I*)m._StatHisto->Clone();
   _orientation = m._orientation;
   _rotate = m._rotate;
   _volume = m._volume;
   _Name = m._Name;
   _plane = NULL;
   _isSurfaceQuantity = false;
   pdg = TDatabasePDG::Instance();
   pdg -> ReadPDGTable();
 }

PndRadMapBoxMesh::PndRadMapBoxMesh	(	const char *	Name,
		int	Xbins,
		Double_t	Xlow,
		Double_t	Xhigh,
		int	Ybins,
		Double_t	Ylow,
		Double_t	Yhigh,
		int	Zbins,
		Double_t	Zlow,
		Double_t	Zhigh
	)

Definition at line 101 of file PndRadMapBoxMesh.cxx.

References _filter, _isSurfaceQuantity, _Name, _verbose, _Xbins, _Xhigh, _Xlow, _Ybins, _Yhigh, _Ylow, _Zbins, _Zhigh, _Zlow, pdg, and TString.

                                                                             {//in cm!!!!
   _verbose = 0;
   _Name = TString(name);
   _Xbins = Xbins; _Xlow = Xlow; _Xhigh = Xhigh;
   _Ybins = Ybins; _Ylow = Ylow; _Yhigh = Yhigh;
   _Zbins = Zbins; _Zlow = Zlow; _Zhigh = Zhigh;
   _filter = TFormula("formula", "1");
 
   if((_Xlow == _Xhigh) ||
      (_Ylow == _Yhigh) ||
      (_Zlow == _Zhigh))
     _isSurfaceQuantity = true;
   pdg = TDatabasePDG::Instance();
   pdg -> ReadPDGTable();
 
   // if(_verbose)
   std::cout << "PndRadMapBoxMesh::PndRadMapBoxMesh()\n"
             << _Xlow << ' ' << _Xhigh << ' '
             << _Ylow << ' ' << _Yhigh << ' '
             << _Zlow << ' ' << _Zhigh << ' '
             << _isSurfaceQuantity << std::endl;
 
 }

PndRadMapBoxMesh::PndRadMapBoxMesh	(	const char *	name,
		int	xbins,
		Double_t	xlow,
		Double_t	xhigh,
		int	ybins,
		Double_t	ylow,
		Double_t	yhigh,
		Double_t	zlow,
		Double_t	zhigh,
		orientation	plane = `ZX`,
		quantity	Quantity = `Edep`
	)

Definition at line 128 of file PndRadMapBoxMesh.cxx.

References _isSurfaceQuantity, _Name, _orientation, _verbose, _Xbins, _Xhigh, _Xlow, _Ybins, _Yhigh, _Ylow, _Zbins, _Zhigh, _Zlow, pdg, SetQuantity(), TString, XY, XZ, YX, YZ, ZX, and ZY.

                                                      {//in cm!!!!
 
   _verbose = 0;
 
   _Name = TString(name);
 
   _orientation = plane;
   SetQuantity(Quantity);
 
   switch (_orientation){
   case XY:
     _Xbins = xbins;
     _Xlow  = xlow;
     _Xhigh = xhigh;
 
     _Ybins = ybins;
     _Ylow  = ylow;
     _Yhigh = yhigh;
 
     _Zlow  = zlow;
     _Zhigh = zhigh;
     break;
   case YX:
     _Xbins = ybins;
     _Xlow = ylow;
     _Xhigh = yhigh;
 
     _Ybins = xbins;
     _Ylow = xlow;
     _Yhigh = xhigh;
 
     _Zlow  = zlow;
     _Zhigh = zhigh;
     break;
 
   case XZ:
     _Xbins = xbins;
     _Xlow = xlow;
     _Xhigh = xhigh;
 
     _Ylow  = zlow;
     _Yhigh = zhigh;
 
     _Zbins = ybins;
     _Zlow = ylow;
     _Zhigh = yhigh;
     break;
   case ZX:
     _Xbins = ybins;
     _Xlow = ylow;
     _Xhigh = yhigh;
 
     _Zbins = xbins;
     _Zlow = xlow;
     _Zhigh = xhigh;
 
     _Ylow  = zlow;
     _Yhigh = zhigh;
     break;
 
   case YZ:
     _Xlow  = zlow;
     _Xhigh = zhigh;
 
     _Ybins = xbins;
     _Ylow = xlow;
     _Yhigh = xhigh;
 
     _Zbins = ybins;
     _Zlow = ylow;
     _Zhigh = yhigh;
     break;
   case ZY:
     _Xlow  = zlow;
     _Xhigh = zhigh;
 
     _Ybins = ybins;
     _Ylow = ylow;
     _Yhigh = yhigh;
 
     _Zbins = xbins;
     _Zlow = xlow;
     _Zhigh = xhigh;
     break;
   };
 
   if((_Xlow == _Xhigh) ||
      (_Ylow == _Yhigh) ||
      (_Zlow == _Zhigh))
     _isSurfaceQuantity = true;
 
   pdg = TDatabasePDG::Instance();
   pdg -> ReadPDGTable();
 }

PndRadMapBoxMesh::~PndRadMapBoxMesh ( )

Definition at line 228 of file PndRadMapBoxMesh.cxx.

References _MeshHisto, _StatHisto, and _verbose.

                                    {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::~PndRadMapBoxMesh()\n";
   delete _StatHisto;
   delete _MeshHisto;
 }

Member Function Documentation

Double_t PndRadMapBoxMesh::CalcFluence ( FairRadMapPoint * p )

Definition at line 834 of file PndRadMapBoxMesh.cxx.

References _MeshHisto, _orientation, _StatHisto, _verbose, _Xhigh, _Xlow, _Yhigh, _Ylow, _Zhigh, _Zlow, counter, Double_t, dx, dy, dz, IsInside(), PndRadMapPlane::LineIntersection(), momentumfilter(), Transform(), verbose, X, XY, XZ, Y, YX, YZ, Z, ZX, and ZY.

                                                         {
 
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::CalcFluence(FairRadMapPoint *p)" << std::endl;
 
   Double_t xx, yy, zz, dx = 0, dy = 0, dz =0;
   // float eff;
   std::vector<Corner> CV;
   std::vector<PndRadMapPlane> PV;
   //The Z should be between zlow and zhigh!!
   Transform(p->GetXOut(), p->GetYOut(), p->GetZOut(),
             xx, yy, zz);
   std::pair<int, double> minp;
   bool OK= false, OK2= false;
 
   if(IsInside(p)){
     if(IsInside(xx, yy, zz)) OK = true;
     switch (_orientation){
     case XY: case YX:
       dz = _Zlow;
       break;
     case XZ: case ZX:
       dy = _Ylow;
       break;
     case ZY: case YZ:
       dx = _Xlow;
       break;
     }
 
     TVector3 poststepv   = TVector3(p->GetXOut(), p->GetYOut(), p->GetZOut());
     TVector3 prestepv    = TVector3(p->GetX(), p->GetY(), p->GetZ());
     unsigned int prebin  = _MeshHisto->FindBin(prestepv.X(),
                                                prestepv.Y());
     unsigned int postbin = _MeshHisto->FindBin(poststepv.X(),
                                                poststepv.Y());
     PndRadMapPlane plane;
     int prebinX  = 0, prebinY  = 0, prebinZ  = 0;
     int postbinX = 0, postbinY = 0, postbinZ = 0;
     _MeshHisto->GetBinXYZ(postbin, postbinX, postbinY, postbinZ);
     _MeshHisto->GetBinXYZ(prebin, prebinX, prebinY, prebinZ);
     TVector3 tn(0, 0, 0);
 
     if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 1" << std::endl;
 
     // 1, only two planes: the narrowest bin lower and higher plane
     // 2, if the resulting two points in the range of the other two axis-limits
     // 3, start the stepping from one of those points
 
     // -------- 1 -------
 
     TVector3 isl, lbeg, lend;
     double min_x, min_y, min_z, max_x, max_y, max_z;
 
     if(!OK){//if the 'out' is not in the given range
       //the 'in' and the 'out' is outside of the given range, so it can be that it will cross it
 
       if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 2" << std::endl;
       OK2 = false;
       min_x = _Xlow;
       min_y = _Ylow;
       min_z = _Zlow;
 
       max_x = _Xhigh;
       max_y = _Yhigh;
       max_z = _Zhigh;
       if(verbose) std::cout << "min_z: " << min_z << ' ' << max_z << ' ' << dz << std::endl;
 
       //the 'front' and 'back' XY planes
       plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz),
                              TVector3(max_x+dx, min_y+dy, min_z+dz),
                              TVector3(min_x+dx, max_y+dy, min_z+dz));//XY @ z=z_min
       isl = plane.LineIntersection(prestepv, poststepv);
       if(((min_x < isl.X()) && (isl.X() < max_x)) &&
          ((min_y < isl.Y()) && (isl.Y() < max_y))){
         OK2 = true;
       }
       if(verbose){
         std::cout << "Xcross1: " << min_x << ' ' << isl.X() << ' ' << max_x << std::endl;
         std::cout << "Ycross1: " << min_y << ' ' << isl.Y() << ' ' << max_y << std::endl;
         std::cout << "=> " << OK2 << std::endl;
       }
       lbeg = isl;
 
       plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, max_z+dz),
                              TVector3(max_x+dx, min_y+dy, max_z+dz),
                              TVector3(min_x+dx, max_y+dy, max_z+dz));//XY @ z-z_max
       isl = plane.LineIntersection(prestepv, poststepv);
       if(((min_x < isl.X()) && (isl.X() < max_x)) &&
          ((min_y < isl.Y()) && (isl.Y() < max_y))){
         OK2 = OK2 && true;
       }
       if(verbose){
         std::cout << "Xcross2: " << min_x << ' ' << isl.X() << ' ' << max_x << std::endl;
         std::cout << "Ycross2: " << min_y << ' ' << isl.Y() << ' ' << max_y << std::endl;
         std::cout << "=> " << OK2 << std::endl;
       }
       lend = isl;
 
       if(OK2){//the intersections are on the XY planes withing boarders at both z_min and z_max
         if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 3" << std::endl;
         prestepv  = lbeg;
         poststepv = lend;
         prebin    = _MeshHisto->FindBin(prestepv.X(),
                                         prestepv.Y());
         postbin   = _MeshHisto->FindBin(poststepv.X(),
                                         poststepv.Y());
         _MeshHisto->GetBinXYZ(prebin, prebinX, prebinY, prebinZ);
         _MeshHisto->GetBinXYZ(postbin, postbinX, postbinY, postbinZ);
       }else{
         if(verbose) std::cout << "PndRadMapBoxMesh::CalcFluence - 4" << std::endl;
       }
     }else{
       // if(prebin != postbin){
       int counter = 0;
       while((prebin != postbin) && !tn.Z()){
         counter++;
         PV.clear();
 
         //       max
         //  /---/      y|
         // /---/|       |   /z
         // ||  ||       |  /
         // |/--|/       | /
         // -----        |/
         // min          +----------
         //                        x
 
 
         // .              corner1          corner2          corner3
         //plane1   [(min, min, min), (min, min, max), (max, min, min)] ZX
         //plane2   [(min, min, min), (max, min, min), (min, max, min)] XY
         //plane3   [(min, min, min), (min, min, max), (min, max, min)] ZY
 
         //plane4   [(min, max, min), (min, max, max), (max, max, min)] ZX
         //plane5   [(min, min, max), (max, min, max), (min, max, max)] XY
         //plane6   [(max, min, min), (max, min, max), (max, max, min)] ZY
 
 
         // //plane4   [(max, max, max), (max, max, min), (min, max, max)] XZ
         // //plane5   [(max, max, max), (min, max, max), (max, min, max)]
         // //plane6   [(max, max, max), (max, max, min), (max, min, max)]
 
 
         min_x = (_MeshHisto->GetXaxis()->GetBinLowEdge(prebinX) < 0) ? 0 : _MeshHisto->GetXaxis()->GetBinLowEdge(prebinX);
         min_y = (_MeshHisto->GetYaxis()->GetBinLowEdge(prebinY) < 0) ? 0 : _MeshHisto->GetYaxis()->GetBinLowEdge(prebinY);
         min_z = (_MeshHisto->GetZaxis()->GetBinLowEdge(prebinZ) < 0) ? 0 : _MeshHisto->GetZaxis()->GetBinLowEdge(prebinZ);
 
         max_x = (_MeshHisto->GetXaxis()->GetBinUpEdge(prebinX) < 0) ? 0 : _MeshHisto->GetXaxis()->GetBinUpEdge(prebinX);
         max_y = (_MeshHisto->GetYaxis()->GetBinUpEdge(prebinY) < 0) ? 0 : _MeshHisto->GetYaxis()->GetBinUpEdge(prebinY);
         max_z = (_MeshHisto->GetZaxis()->GetBinUpEdge(prebinZ) < 0) ? 0 : _MeshHisto->GetZaxis()->GetBinUpEdge(prebinZ);
         max_z += _MeshHisto->GetZaxis()->GetBinWidth(1);
 
         if(verbose){
           std::cout << counter << " - Corners [[" << min_x << ", " << min_y << ", " << min_z << "] - ";
           std::cout <<                        "[" << max_x << ", " << max_y << ", " << max_z << "]]\n";
         }
 
         //plane1
         plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz),
                                TVector3(min_x+dx, min_y+dy, max_z+dz),
                                TVector3(max_x+dx, min_y+dy, min_z+dz));
         PV.push_back(plane);
         //plane2
         plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz),
                                TVector3(max_x+dx, min_y+dy, min_z+dz),
                                TVector3(min_x+dx, max_y+dy, min_z+dz));
         PV.push_back(plane);
         //plane3
         plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, min_z+dz),
                                TVector3(min_x+dx, min_y+dy, max_z+dz),
                                TVector3(min_x+dx, max_y+dy, min_z+dz));
         PV.push_back(plane);
         //plane4
         plane = PndRadMapPlane(TVector3(min_x+dx, max_y+dy, min_z+dz),
                                TVector3(min_x+dx, max_y+dy, max_z+dz),
                                TVector3(max_x+dx, max_y+dy, min_z+dz));
         PV.push_back(plane);
         //plane5
         plane = PndRadMapPlane(TVector3(min_x+dx, min_y+dy, max_z+dz),
                                TVector3(max_x+dx, min_y+dy, max_z+dz),
                                TVector3(min_x+dx, max_y+dy, max_z+dz));
         PV.push_back(plane);
         //plane6
         plane = PndRadMapPlane(TVector3(max_x+dx, min_y+dy, min_z+dz),
                                TVector3(max_x+dx, min_y+dy, max_z+dz),
                                TVector3(max_x+dx, max_y+dy, min_z+dz));
         PV.push_back(plane);
 
         if(verbose)
           std::cout << "Start\n["
                     << prebin << " {" << prestepv.X() << ' ' << prestepv.Y() << ' ' << prestepv.Z() << "}("
                     << prebinX << ' ' << prebinY << ' ' << prebinZ << ")\n"
                     << ' ' << postbin << " {" << poststepv.X() << ' ' << poststepv.Y() << ' ' << poststepv.Z() << "}("
                     << postbinX << ' ' << postbinY << ' ' << postbinZ << ")]\n\n";
 
         minp = std::make_pair(-1, 111111);//(which plane, what is the disrance)
 
         //calculate the intersection with every each planes
         isl = PV.at(0).LineIntersection(prestepv, poststepv);
         if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)))
           if(minp.second > (isl-prestepv).Mag())
             minp = std::make_pair(0, (isl-prestepv).Mag());
         if(verbose)
           std::cout << "XZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") "
                     << (poststepv-prestepv).Mag() << ", "
                     << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") "
                     << (isl-prestepv).Mag() << ' '
                     << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' '
                     << (minp.second > (isl-prestepv).Mag()) << ' '
                     << "[" << minp.first << ", " << minp.second << std::endl;
 
         if(!OK2){ //FIXME: OK2 was not set before!
 
           isl = PV.at(1).LineIntersection(prestepv, poststepv);
           if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)))
             if(minp.second > (isl-prestepv).Mag())
               minp = std::make_pair(1, (isl-prestepv).Mag());
           if(verbose)
             std::cout << "XY (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") "
                       << (poststepv-prestepv).Mag() << ", "
                       << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") "
                       << (isl-prestepv).Mag() << ' '
                       << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' '
                       << (minp.second > (isl-prestepv).Mag()) << ' '
                       << "[" << minp.first << ", " << minp.second << "]" << std::endl;
         }
 
         isl = PV.at(2).LineIntersection(prestepv, poststepv);
         if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)))
           if(minp.second > (isl-prestepv).Mag())
             minp = std::make_pair(2, (isl-prestepv).Mag());
         if(verbose)
           std::cout << "YZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") "
                     << (poststepv-prestepv).Mag() << ", "
                     << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") "
                     << (isl-prestepv).Mag() << ' '
                     << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' '
                     << (minp.second > (isl-prestepv).Mag()) << ' '
                     << "[" << minp.first << ", " << minp.second << "]" << std::endl;
 
         isl = PV.at(3).LineIntersection(prestepv, poststepv);
         if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)))
           if(minp.second > (isl-prestepv).Mag())
             minp = std::make_pair(3, (isl-prestepv).Mag());
         if(verbose)
           std::cout << "XZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") "
                     << (poststepv-prestepv).Mag() << ", "
                     << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") "
                     << (isl-prestepv).Mag() << ' '
                     << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' '
                     << (minp.second > (isl-prestepv).Mag()) << ' '
                     << "[" << minp.first << ", " << minp.second  << "]" << std::endl;
 
         if(!OK2){ //FIXME: OK2 was not set before!
           isl = PV.at(4).LineIntersection(prestepv, poststepv);
           if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)))
             if(minp.second > (isl-prestepv).Mag())
               minp = std::make_pair(4, (isl-prestepv).Mag());
           if(verbose)
             std::cout << "XY (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") "
                       << (poststepv-prestepv).Mag() << ", "
                       << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") "
                       << (isl-prestepv).Mag() << ' '
                       << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' '
                       << (minp.second > (isl-prestepv).Mag()) << ' '
                       << "[" << minp.first << ", " << minp.second  << "]" << std::endl;
         }
 
         isl = PV.at(5).LineIntersection(prestepv, poststepv);
         if(momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)))
           if(minp.second > (isl-prestepv).Mag())
             minp = std::make_pair(5, (isl-prestepv).Mag());
         if(verbose)
           std::cout << "YZ (" << (isl).X() << ", " << (isl).Y() << ", " << (isl).Z() << ") "
                     << (poststepv-prestepv).Mag() << ", "
                     << "(" << (isl-prestepv).X() << ", " << (isl-prestepv).Y() << ", " << (isl-prestepv).Z() << ") "
                     << (isl-prestepv).Mag() << ' '
                     << momentumfilter(TVector3(p->GetPx(), p->GetPy(), p->GetPz()), (isl-prestepv)) << ' '
                     << (minp.second > (isl-prestepv).Mag()) << ' '
                     << "[" << minp.first << ", " << minp.second  << "]" << std::endl;
 
         if(verbose)
           std::cout << "\nmomentum: "
                     << p->GetPx() << ' '
                     << p->GetPy() << ' '
                     << p->GetPz() << std::endl;
 
         //minp stores the smallest possible step
         if(0 <= minp.first){
           if(verbose){
             std::cout << "\nmomentum: "
                       << p->GetPx() << ' '
                       << p->GetPy() << ' '
                       << p->GetPz() << std::endl;
             std::cout << "[(----)] " << minp.first << ' ' << minp.second << std::endl;
             std::cout << "normal: "
                       << PV.at(minp.first).Normal().X() << ' '
                       << PV.at(minp.first).Normal().Y() << ' '
                       << PV.at(minp.first).Normal().Z() << std::endl;
           }
 
           isl = PV.at(minp.first).LineIntersection(prestepv, poststepv);
           Double_t step = (isl-prestepv).Mag();
           prestepv = isl;
 
           _MeshHisto->Fill(prestepv.X(), prestepv.Y(), step);
           _StatHisto->Fill(prestepv.X(), prestepv.Y());
 
           tn = PV.at(minp.first).Normal();
           float ds = 0;
           if(tn.X()){
             ds = _MeshHisto->GetXaxis()->GetBinWidth(1)/10.;
             if(verbose) std::cout << "Xprebin: " << prebin << " ---> ";
             if(p->GetPx() > 0){
               prebin = _MeshHisto->FindBin(isl.X() + ds,
                                            isl.Y());
               if(verbose) std::cout << "(" << isl.X() << " ---> " << isl.X() + ds << ' '  << p->GetPx() << " ) ";
             }else{
               prebin = _MeshHisto->FindBin(isl.X() - ds,
                                            isl.Y());
               if(verbose) std::cout << "(" << isl.X() << " ---> " << isl.X() - ds << ' '  << p->GetPx() << " ) ";
             }
             if(verbose) std::cout << prebin << std::endl;
           }
           if(tn.Y()){
             ds = _MeshHisto->GetYaxis()->GetBinWidth(1)/10.;
             if(verbose) std::cout << "Yprebin: " << prebin << " ---> ";
             if(p->GetPy() > 0){
               prebin = _MeshHisto->FindBin(isl.X(),
                                            isl.Y() + ds);
               if(verbose) std::cout << "1(" << isl.Y() << " ---> " << isl.Y() + ds << " ) ";
             }else{
               prebin = _MeshHisto->FindBin(isl.X(),
                                            isl.Y()-ds);
               if(verbose) std::cout << "1(" << isl.Y() << " ---> " << isl.Y() - ds << " ) ";
             }
             if(verbose) std::cout << prebin << std::endl;
           }
 
           _MeshHisto->GetBinXYZ(prebin, prebinX, prebinY, prebinZ);
           if(verbose) std::cout << "Stop\n["
                                 << prebin << " {" << prestepv.X() << ' ' << prestepv.Y() << ' ' << prestepv.Z() << "}("
                                 << prebinX << ' ' << prebinY << ' ' << prebinZ << ")\n"
                                 << ' ' << postbin << " {" << poststepv.X() << ' ' << poststepv.Y() << ' ' << poststepv.Z() << "}("
                                 << postbinX << ' ' << postbinY << ' ' << postbinZ << ")]\n\n";
 
           if(IsInside(p->GetX(), p->GetY(), p->GetZ())){
             if(verbose) std::cout << " Bingo\n";
           }
         }
       }
     }
   }
   return minp.second;
 }

void PndRadMapBoxMesh::Fill ( FairRadMapPoint * p )

Definition at line 459 of file PndRadMapBoxMesh.cxx.

References _EnergyHisto, _filter, _orientation, _plane, _quantity, _verbose, _volume, Density, Dose, Double_t, Edep, EnergyFluence, Fluence, Flux, InterSection, Kerma, PndRadMapPlane::LineIntersection(), Mass, mom, pdg, pdgpart, pid(), SimpleFluence, CAMath::Sqrt(), theta, Twos, val, XY, XZ, YX, YZ, ZX, and ZY.

Referenced by Fill().

                                              {
   if(_verbose)
     std::cout << "void PndRadMapBoxMesh::Fill(FairRadMapPoint *p)\n";
 
   Int_t pid = p->GetPdg();
   Double_t mom = TMath::Sqrt((p->GetPx()*p->GetPx()) +
                              (p->GetPy()*p->GetPy()) +
                              (p->GetPz()*p->GetPz()));
   pdgpart = pdg->GetParticle(pid);
   Int_t cha = -1000;// to be implemented
   Double_t pdgmass = -1;
   Double_t val = 0;
 
   Double_t edep = 0;
   Double_t dens = 0;
   Double_t mass = p->GetMass();
   TVector3 poststepv, prestepv;
   Double_t lX = 0, lY = 0, lZ = 0;
 
   if(pdgpart){
     pdgmass = pdgpart->Mass();
     cha = pdgpart->Charge();
   }
 
   if(_filter.Eval(pid, cha, mom)){
 
     val = 0;
 
     edep = p->GetEnergyLoss()*1.60217657*1e-10;//GeV->Joule
     dens = p->GetDensity();
     dens *= 0.001;//g/cm3 -> kg/cm3
     mass = p->GetMass();
     poststepv   = TVector3(p->GetXOut(), p->GetYOut(), p->GetZOut());
     prestepv    = TVector3(p->GetX(), p->GetY(), p->GetZ());
     lX = 0;
     lY = 0;
     lZ = 0;
 
     if(_verbose)
       std::cout << "Density: " << dens
                 << ", mass: " << mass
                 << " (" << (dens*_volume)
                 << ", " << _volume << ") "
                 << dens*mass << ' '
                 << "charge: " << cha << "\n";
 
     lX = p->GetXOut();
     lY = p->GetYOut();
     lZ = p->GetZOut();
 
     bool OK = false;
     switch (_quantity){
 
     case Edep:
       val = edep;
       if(!std::isnan(val) && val !=0)
         // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
         Fill(lX, lY, lZ, val);
       break;
 
     case Dose:
       val = edep/(dens*_volume);//mass=density*vol,dens[kg/cm^3],vol[cm^3]!
       if(!std::isnan(val) && val !=0){
         // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
         Fill(lX, lY, lZ, val);
       }
       break;
 
     case Fluence:
       InterSection = _plane->LineIntersection(poststepv, prestepv);
       switch (_orientation){
       case XY: case YX:
         if(((prestepv.Z()  <= InterSection.Z()) && (InterSection.Z() <= poststepv.Z())) ||
            ((poststepv.Z() <= InterSection.Z()) && (InterSection.Z() <= prestepv.Z())))
           OK = true;
         break;
 
       case XZ: case ZX:
         if(((prestepv.Y()  <= InterSection.Y()) && (InterSection.Y() <= poststepv.Y())) ||
            ((poststepv.Y() <= InterSection.Y()) && (InterSection.Y() <= prestepv.Y())))
           OK = true;
         break;
 
       case YZ: case ZY:
         if(((prestepv.X()  <= InterSection.X()) && (InterSection.X() <= poststepv.X())) ||
            ((poststepv.X() <= InterSection.X()) && (InterSection.X() <= prestepv.X())))
           OK = true;
         break;
       };
       // std::cout << "test0\n"
       //           << "[" << prestepv.X()  << "<=" << InterSection.X() << "<=" << poststepv.X() << "] - " << OK << "\n"
       //           << "[" << poststepv.X() << "<=" << InterSection.X() << "<=" << prestepv.X()  << "] - " << OK << "\n"
       //           << "[" << prestepv.Y()  << "<=" << InterSection.Y() << "<=" << poststepv.Y() << "] - " << OK << "\n"
       //           << "[" << poststepv.Y() << "<=" << InterSection.Y() << "<=" << prestepv.Y()  << "] - " << OK << "\n"
       //           << "[" << prestepv.Z()  << "<=" << InterSection.Z() << "<=" << poststepv.Z() << "] - " << OK << "\n"
       //           << "[" << poststepv.Z() << "<=" << InterSection.Z() << "<=" << prestepv.Z()  << "] - " << OK << "\n\n";
 
       // // if(((prestepv.Z() <= InterSection.Z()) && (InterSection.Z() <= poststepv.Z())) ||
       // //    ((poststepv.Z() <= InterSection.Z()) && (InterSection.Z() <= prestepv.Z()))){
       if(OK){
         // std::cout << "testing 0...\n";
         // std::cout << "testing 0... "
         //   << _Xlow << "<=" << InterSection.X() << "<=" << _Xhigh << " -> "
         //   << ((_Xlow <= InterSection.X()) && (InterSection.X() <= _Xhigh)) << ' ' << (InterSection.X() <= _Xhigh) << ' ' << (_Xlow <= InterSection.X()) << ' '
         //   << _Ylow << "<=" << InterSection.Y() << "<=" << _Yhigh << " -> "
         //   << ((_Ylow <= InterSection.Y()) && (InterSection.Y() <= _Yhigh)) << ' '
         //   << _Zlow << "<=" << InterSection.Z() << "<=" << _Zhigh << " -> "
         //   << ((_Zlow <= InterSection.Z()) && (InterSection.Z() <= _Zhigh)) << std::endl;
 
         Fill(InterSection.X(),
              InterSection.Y(),
              InterSection.Z());  //normalisation to the area ... now it is OK
 
 
         Double_t theta = InterSection.Theta()*TMath::RadToDeg();
         Double_t K = TMath::Sqrt(mom*mom + pdgmass*pdgmass) - pdgmass;
 
         _EnergyHisto->Fill(theta, K);           //uncommented by MP, Sep 8
         // val = CalcFluence(p);
         // if(!isnan(val) && val !=0){
         //   Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
       }
       break;
 
     case Density:
       val = dens;
       if(!std::isnan(val) && val !=0){
         Fill(lX, lY, lZ, val);
         // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
       }
       break;
 
     case Mass:
       val = (dens*_volume);
       if(!std::isnan(val) && val !=0){
         Fill(lX, lY, lZ, val);
         // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
       }
       break;
 
     case SimpleFluence:{
       poststepv   = TVector3(p->GetXOut(), p->GetYOut(), p->GetZOut());
       prestepv    = TVector3(p->GetX(), p->GetY(), p->GetZ());
       val = (poststepv-prestepv).Mag()/1e2;// cm->m
       val /= (_volume/1e6);//cm^3->m^3 -> 1/m^2
       if(!std::isnan(val) && val!=0){
         if(_verbose > 1){
           std::cout << "Pre: " << prestepv.X() << ", " << prestepv.Y() << ", " << prestepv.Z() << "\n";
           std::cout << "Pre: " << poststepv.X() << ", " << poststepv.Y() << ", " << poststepv.Z() << "\n";
           std::cout << "diff: " << (poststepv-prestepv).Mag() << "; val: " << val << std::endl;
           std::cout << "volume: " << _volume << "\n\n";
         }
         Fill(lX, lY, lZ, val);
         // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
       }
     }break;
 
     case Twos:
       val = 2;
       if(!std::isnan(val) && val !=0){
         Fill(lX, lY, lZ, val);
         // Fill(p->GetXOut(), p->GetYOut(), p->GetZOut(), val);
       }
       break;
     case EnergyFluence:
       break;
     case Flux:
       break;
     case Kerma:
       break;
     }
   }
 }

void PndRadMapBoxMesh::Fill	(	Int_t	gBin,
		Double_t	val
	)

protected

Definition at line 634 of file PndRadMapBoxMesh.cxx.

References _MeshHisto, _verbose, and Fill().

                                                    {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::Fill(" << gBin << ", " << val << ")\n";
 
   Int_t binX, binY, binZ;
   _MeshHisto->GetBinXYZ(gBin, binX, binY, binZ);
   Fill(binX, binY, binZ, val);
 }

void PndRadMapBoxMesh::Fill	(	Double_t	X,
		Double_t	Y,
		Double_t	Z,
		Double_t	we = `1`
	)

protected

Definition at line 643 of file PndRadMapBoxMesh.cxx.

References _MeshHisto, _orientation, _StatHisto, _verbose, _X, _Y, _Z, IsInside(), Transform(), XY, XZ, YX, YZ, ZX, and ZY.

                                                                           {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::Fill(" << X << ", " << Y << ", " << Z << ", " << we << ")\n";
 
   Transform(X, Y, Z);
 
   // std::cout << "0 - PndRadMapBoxMesh::Fill(" << X << ", " << Y << ", " << Z << ", " << we << ")\n";
   // std::cout << "1 - PndRadMapBoxMesh::Fill(" << _X << ", " << _Y << ", " << _Z << ", " << we << ")\n";
 
   // if(IsInside(_X, _Y, _Z)){
   if(IsInside()){
     // std::cout << "2 - PndRadMapBoxMesh::Fill(" << _X << ", " << _Y << ", " << _Z << ", " << we << ")\n";
 
     switch (_orientation){
 
     case XY:
       _MeshHisto->Fill(_X, _Y, we);
       _StatHisto->Fill(_X, _Y);
       break;
     case YX:
       _MeshHisto->Fill(_Y, _X, we);
       _StatHisto->Fill(_Y, _X);
       break;
 
     case XZ:
       _MeshHisto->Fill(_X, _Z, we);
       _StatHisto->Fill(_X, _Z);
       break;
     case ZX:
       _MeshHisto->Fill(_Z, _X, we);
       _StatHisto->Fill(_Z, _X);
       break;
 
     case YZ:
       _MeshHisto->Fill(_Y, _Z, we);
       _StatHisto->Fill(_Y, _Z);
       break;
     case ZY:
       _MeshHisto->Fill(_Z, _Y, we);
       _StatHisto->Fill(_Z, _Y);
       break;
     };
   }
 }

TH2D * PndRadMapBoxMesh::GetHisto ( )

Definition at line 694 of file PndRadMapBoxMesh.cxx.

References _MeshHisto, and _verbose.

                                 {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::GetHisto()\n";
   return _MeshHisto;
 }

PndRadMapPlane * PndRadMapBoxMesh::GetPlane ( )

Definition at line 700 of file PndRadMapBoxMesh.cxx.

References _plane, and _verbose.

                                           {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::GetPlane()\n";
   return _plane;
 }

bool PndRadMapBoxMesh::IsInside	(	Double_t	X,
		Double_t	Y,
		Double_t	Z
	)

protected

Definition at line 736 of file PndRadMapBoxMesh.cxx.

References _verbose, _Xhigh, _Xlow, _Yhigh, _Ylow, _Zhigh, _Zlow, and ok.

                                                                  {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::IsInside(" << X << ", " << Y << ", " << Z << ")\n";
 
   bool ok = false;
 
   if(((_Xlow <= X) && (X <= _Xhigh)) &&
      ((_Ylow <= Y) && (Y <= _Yhigh)) &&
      ((_Zlow <= Z) && (Z <= _Zhigh)))
     ok = true;
   return ok;
 }

bool PndRadMapBoxMesh::IsInside ( FairRadMapPoint * p )

protected

Definition at line 780 of file PndRadMapBoxMesh.cxx.

References _verbose, _Xhigh, _Xlow, _Yhigh, _Ylow, _Zhigh, _Zlow, ok, and Transform().

                                                  {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::IsInside(FairRadMapPoint *p)\n";
 
   bool ok = false;
   double xin, yin, zin, xout, yout, zout;
   Transform(p->GetX(), p->GetY(), p->GetZ(),
             xin, yin, zin);
   Transform(p->GetXOut(), p->GetYOut(), p->GetZOut(),
             xout, yout, zout);
 
   if((((xin   <= _Xlow) && (_Xhigh <= xout))   ||
       ((xout  <= _Xlow) && (_Xhigh <= xin ))   ||
       ((_Xlow <= xout)  && (xout   <= _Xhigh))) &&
      (((yin   <= _Ylow) && (_Yhigh <= yout))   ||
       ((yout  <= _Ylow) && (_Yhigh <= yin ))   ||
       ((_Ylow <= yout)  && (yout   <= _Yhigh))) &&
      (((zin   <= _Zlow) && (_Zhigh <= zout))   ||
       ((zout  <= _Zlow) && (_Zhigh <= zin ))   ||
       ((_Zlow <= zout)  && (zout   <= _Zhigh))))
     ok = true;
 
   return ok;
 }

bool PndRadMapBoxMesh::IsInside ( )

protected

Definition at line 749 of file PndRadMapBoxMesh.cxx.

References _orientation, _verbose, _X, _Xhigh, _Xlow, _Y, _Yhigh, _Ylow, _Z, _Zhigh, _Zlow, ok, XY, XZ, YX, YZ, ZX, and ZY.

Referenced by CalcFluence(), and Fill().

                                {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::IsInside()\n";
 
   bool ok = false;
 
   switch (_orientation){
   case XY: case YX:
   if(((_Xlow <= _X) && (_X <= _Xhigh)) &&
      ((_Ylow <= _Y) && (_Y <= _Yhigh)))
     ok = true;
   break;
   case XZ: case ZX:
   if(((_Xlow <= _X) && (_X <= _Xhigh)) &&
      ((_Zlow <= _Z) && (_Z <= _Zhigh)))
     ok = true;
   break;
   case YZ: case ZY:
   if(((_Ylow <= _Y) && (_Y <= _Yhigh)) &&
      ((_Zlow <= _Z) && (_Z <= _Zhigh)))
     ok = true;
   break;
   };
 
   // if(((_Xlow <= _X) && (_X <= _Xhigh)) &&
   //    ((_Ylow <= _Y) && (_Y <= _Yhigh)) &&
   //    ((_Zlow <= _Z) && (_Z <= _Zhigh)))
   //   ok = true;
   return ok;
 }

void PndRadMapBoxMesh::makeHisto	(	const char *	Orient,
		Double_t	rotate,
		int	Hbins,
		Double_t	Hlow,
		Double_t	Hhigh,
		int	Vbins,
		Double_t	Vlow,
		Double_t	Vhigh,
		Double_t	dlow
	)

protected

Definition at line 805 of file PndRadMapBoxMesh.cxx.

References _axis, _EnergyHisto, _isSurfaceQuantity, _MeshHisto, _Name, _orientation, _plane, _rotate, _StatHisto, _verbose, and C().

Referenced by SetOrientation().

                                                 {
   if(_verbose){
     std::cout << "PndRadMapBoxMesh::makeHisto(" << Orient << ", " << rotate << ",\n";
     std::cout << "                   " << Hbins << ", " << Hlow << ", " << Hhigh << ",\n";
     std::cout << "                   " << Vbins << ", " << Vlow << ", " << Vhigh << ")\n";
   }
   char C[32];
   sprintf(C, "%s_Histo_%s_%g", _Name.Data(), Orient, _rotate);
   _MeshHisto = new TH2D(C, C,
                         Hbins, Hlow, Hhigh,
                         Vbins, Vlow, Vhigh);
   sprintf(C, "%s_StatHisto_%s_%g", _Name.Data(), Orient, _rotate);
   _StatHisto = new TH2I(C, C,
                         Hbins, Hlow, Hhigh,
                         Vbins, Vlow, Vhigh);
 
   if(_isSurfaceQuantity){
     sprintf(C, "%s_EnergyHisto_%s_%g", _Name.Data(), Orient, _rotate);
     _EnergyHisto = new TH2D(C, C,
                             180,  0, 180,
                             10000, 0, 10);
     _plane = new PndRadMapPlane(dlow, _rotate, _orientation, _axis);
 
   }
 }

void PndRadMapBoxMesh::Save ( TFile * fout )

Definition at line 706 of file PndRadMapBoxMesh.cxx.

References _EnergyHisto, _isSurfaceQuantity, _MeshHisto, _StatHisto, _verbose, and TString.

                                       {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::Save(TFile* fout)\n";
 
   TObjArray* toa = TString(_MeshHisto->GetName()).Tokenize('_');
   TString dirname = ((TObjString*)toa->At(0))->GetString();
   dirname += '_';
   dirname += ((TObjString*)toa->At(2))->GetString().Data();
   dirname += '_';
   dirname += ((TObjString*)toa->At(3))->GetString().Data();
   std::cout << dirname.Data() << std::endl;
 
   // fout->mkdir(_MeshHisto->GetName());
   // fout->cd   (_MeshHisto->GetName());
   fout->mkdir(dirname.Data());
   fout->cd   (dirname.Data());
   _MeshHisto->Clone()->Write();
   _StatHisto->Clone()->Write();
   if(_isSurfaceQuantity) _EnergyHisto->Clone()->Write();
   fout->cd();
 }

void PndRadMapBoxMesh::Save ( )

Definition at line 728 of file PndRadMapBoxMesh.cxx.

References _EnergyHisto, _isSurfaceQuantity, _MeshHisto, _StatHisto, and _verbose.

                            {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::Save()\n";
   _MeshHisto->Write();
   _StatHisto->Write();
   if(_isSurfaceQuantity) _EnergyHisto->Clone()->Write();
 }

void PndRadMapBoxMesh::Scale ( Double_t sca )

Definition at line 688 of file PndRadMapBoxMesh.cxx.

References _MeshHisto, and _verbose.

                                         {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::Scale(" << sca << ")" << std::endl;
   _MeshHisto->Scale(sca);
 }

void PndRadMapBoxMesh::SetFilter ( const char * filter )

Definition at line 235 of file PndRadMapBoxMesh.cxx.

References _filter, _verbose, and TString.

                                                   {
   //Pid=x, Ch=y, Mom=z,
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::SetFilter(" << filter << ")\n";
   TString temps(filter);
   temps.ReplaceAll("Pid", "x");
   temps.ReplaceAll("Ch" , "y");
   temps.ReplaceAll("Mom", "z");
 
   _filter = TFormula("formula", temps.Data());
 }

void PndRadMapBoxMesh::SetOrientation	(	orientation	plane,
		Double_t	rotate = `99999`,
		axis	Ax = `Xx`
	)

Definition at line 253 of file PndRadMapBoxMesh.cxx.

References _orientation, and _verbose.

                                               {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::SetOrientation(" << plane << ", " << rotate << ", " << Ax << ")\n";
   _orientation = plane;
   SetOrientation(rotate,Ax);
 }

void PndRadMapBoxMesh::SetOrientation	(	Double_t	rotate = `99999`,
		axis	Ax = `Xx`
	)

Definition at line 262 of file PndRadMapBoxMesh.cxx.

References _axis, _isSurfaceQuantity, _Name, _orientation, _rotate, _verbose, _volume, _Xbins, _Xhigh, _Xlow, _Ybins, _Yhigh, _Ylow, _Zbins, _Zhigh, _Zlow, makeHisto(), verbose, XY, XZ, YX, YZ, ZX, and ZY.

                                               {
 
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::SetOrientation(" << rotate << ", " << Ax << ")\n";
   _axis = Ax;
   _rotate = (rotate == 99999) ? 0 : rotate;
 
   if(_verbose)
     std::cout << "Volume[" << _Name.Data() << "][("
               << _Xhigh << ", " << _Xlow << ", " << _Xbins << "),("
               << _Yhigh << ", " << _Ylow << ", " << _Ybins << "),("
               << _Zhigh << ", " << _Zlow << ", " << _Zbins << ")]:\n";
 
   // char C[32];
   switch (_orientation){
   case XY:
     makeHisto("XY", _rotate,
               _Xbins, _Xlow, _Xhigh,
               _Ybins, _Ylow, _Yhigh,
               _Zlow/*, _Zhigh --unused parameter*/);
     break;
   case YX:
     makeHisto("YX", _rotate,
               _Ybins, _Ylow, _Yhigh,
               _Xbins, _Xlow, _Xhigh,
               _Zlow/*, _Zhigh --unused parameter*/);
     break;
 
   case XZ:
     makeHisto("XZ", _rotate,
               _Xbins, _Xlow, _Xhigh,
               _Zbins, _Zlow, _Zhigh,
               _Ylow/*, _Yhigh --unused parameter*/);
     break;
   case ZX:
     makeHisto("ZX", _rotate,
               _Zbins, _Zlow, _Zhigh,
               _Xbins, _Xlow, _Xhigh,
               _Ylow/*, _Yhigh --unused parameter*/);
     break;
 
   case YZ:
     makeHisto("YZ", _rotate,
               _Ybins, _Ylow, _Yhigh,
               _Zbins, _Zlow, _Zhigh,
               _Xlow/*, _Xhigh --unused parameter*/);
     break;
   case ZY:
     makeHisto("ZY", _rotate,
               _Zbins, _Zlow, _Zhigh,
               _Ybins, _Ylow, _Yhigh,
               _Xlow/*, _Xhigh --unused parameter*/);
     break;
   };
 
   if(!_isSurfaceQuantity)
     _volume = ((_Xhigh-_Xlow)/_Xbins) *
       ((_Yhigh-_Ylow)/_Ybins)*
       ((_Zhigh-_Zlow)/_Zbins);
   else
     _volume = 1;
 
   if(verbose)
     std::cout << "Volume: "
               << (_Xhigh-_Xlow)/_Xbins << ", "
               << (_Yhigh-_Ylow)/_Ybins << ", "
               << (_Zhigh-_Zlow)/_Zbins << " -> " << _volume << "\n";
 }

void PndRadMapBoxMesh::SetQuantity ( quantity Quantity = Edep )

Definition at line 247 of file PndRadMapBoxMesh.cxx.

References _quantity, and _verbose.

Referenced by PndRadMapBoxMesh().

                                                    {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::SetQuantity(" << Quantity << ")\n";
   _quantity = Quantity;
 }

void PndRadMapBoxMesh::SetVerbosityLevel ( int verbose = 0 )

Definition at line 333 of file PndRadMapBoxMesh.cxx.

References _verbose.

                                                    {
   _verbose = Verbose;
 }

void PndRadMapBoxMesh::Transform	(	Double_t	X,
		Double_t	Y,
		Double_t	Z
	)

Definition at line 339 of file PndRadMapBoxMesh.cxx.

References _verbose, _X, _Y, _Z, X, Y, and Z.

Referenced by CalcFluence(), Fill(), and IsInside().

                                                                   {
   if(_verbose)
     std::cout << "PndRadMapBoxMesh::Transform(" << X << ", " << Y << ", " << Z << ")\n";
 
   _X = X; _Y = Y; _Z = Z;
   Transform(   X,  Y,  Z,
                _X, _Y, _Z);
 }

void PndRadMapBoxMesh::Transform	(	Double_t	X,
		Double_t	Y,
		Double_t	Z,
		Double_t &	X0,
		Double_t &	Y0,
		Double_t &	Z0
	)

Definition at line 348 of file PndRadMapBoxMesh.cxx.

References _axis, _orientation, _rotate, _verbose, CAMath::Cos(), Double_t, R, CAMath::Sin(), CAMath::Sqrt(), X, Xx, XY, XZ, Y, YX, Yy, YZ, Z, ZX, ZY, and Zz.

                                                                           {
 
   if(_verbose){
     std::cout << "PndRadMapBoxMesh::Transform(" << X << ", " << Y << ", " << Z << ")\n";
     std::cout << "                            " << X0 << ", " << Y0 << ", " << Z0 << ")\n";
   }
   if(_rotate == 99999) _rotate = 0;
   if((_rotate == 0)){
     X0 = X; Y0 = Y; Z0 = Z;
   }else{
 
     Double_t V1 = 0;
     Double_t V2 = 0;
 
     X0 = X;
     Y0 = Y;
     Z0 = Z;
     switch (_axis){
     case Xx:
       V1 = Y;
       V2 = Z;
       X0 = X;
       break;
     case Yy:
       V1 = X;
       V2 = Z;
       Y0 = Y;
       break;
     case Zz:
       V1 = X;
       V2 = Y;
       Z0 = Z;
       break;
     };
 
     Double_t _Vh, _Vv;
     if(_rotate != 0){
       Double_t R = TMath::Sqrt(V1*V1 + V2*V2);
       Double_t rot0 = TMath::ATan(V1/V2)*TMath::RadToDeg();
       if(V2 < 0) rot0 += 180;
       Double_t drot = rot0-_rotate;
       _Vh = R*TMath::Sin(1*drot*TMath::DegToRad());
       _Vv = R*TMath::Cos(1*drot*TMath::DegToRad());
     }else{
       _Vh = V1;
       _Vv = V2;
     }
 
     switch (_axis){
     case Xx:
       switch (_orientation){
       case XY: case YX: case XZ: case ZX:
         Y0 = _Vv;
         Z0 = _Vh;
         break;
 
       case YZ:
         Y0 = _Vh;
         Z0 = _Vv;
         break;
       case ZY:
         Y0 = _Vv;
         Z0 = _Vh;
         break;
       };
       break;
 
     case Yy:
       switch (_orientation){
       case XY: case YX: case YZ: case ZY:
         X0 = _Vh;
         Z0 = _Vv;
         break;
 
       case XZ:
         X0 = _Vh;
         Z0 = _Vv;
         break;
       case ZX:
         X0 = _Vv;
         Z0 = _Vh;
         break;
       };
       break;
 
     case Zz:
       switch (_orientation){
       case XZ: case ZX: case YZ: case ZY:
         X0 = _Vh;
         Y0 = _Vv;
         break;
 
       case XY:
         X0 = _Vh;
         Y0 = _Vv;
         break;
       case YX:
         X0 = _Vv;
         Y0 = _Vh;
         break;
       };
       break;
     };
   }
   if(_verbose){
     std::cout << " -> PndRadMapBoxMesh::Transform(" << X << ", " << Y << ", " << Z << ")\n";
     std::cout << "                                " << X0 << ", " << Y0 << ", " << Z0 << ")\n";
   }
 }