[Gate-users] Track-lenght fluence spectrum
David Boersma
david.boersma at igp.uu.se
Thu Nov 30 16:34:00 CET 2017
Hi Tony,
I don't know of an actor that does this, but a while ago I wrote a crude
ROOT macro (see attachment) to compute something similar to what you
describe (track energy weighted with track-length-in-volume). It assumes
that your volume is a rectangular box and that you have your tracks
collected in a phase space file in ROOT format, as collected with the
GATE phase space actor.
A box has 6 surfaces. First we extend those surfaces into infinite
planes. Then we find for a given track the 6 "punch through points" of
that track through those planes. Then we find which of those ptps are
actually on a box surface. There should be at most 2 such points. If
indeed 2 such points are found, then the length inside the box can be
computed; otherwise (0 points on a box surface) the track apparently
missed the box.
For arbitrary non-rectangular volumes the algorithm needs to be a lot
smarter, especially if the volumes are allowed to be non-convex. :-)
If anyone decides that this crude ROOT macro should be turned into a
GATE actor, then the code needs a lot of polishing and documentation...
HTH,
David
Den 20/11/2017 kl. 12:20, skrev tony youness:
> Dear Gate users,
>
> I am actually working on external radiotherapy applications with 6MV
> photon beam. I need to generate electron fluence spectrum differential
> in energy in a volume of water in terms of *Path length per unit volume
> *(ΦE=∑path length of particle i in V / V) and *not* as particle counting
> (dN/da). Is there an already implemented actor that can do that in GATE?
>
> Another question regarding Filters. Is it possible to apply two filters
> at the same time, for exemple when applying particle filter is it
> possible to filter on particle type (e-) and parentparticle(gamma)?
>
> Best Regards,
>
> Tony Younes
>
>
>
> _______________________________________________
> Gate-users mailing list
> Gate-users at lists.opengatecollaboration.org
> http://lists.opengatecollaboration.org/mailman/listinfo/gate-users
>
-------------- next part --------------
#include <TString.h>
#include <TChain.h>
#include <TMath.h>
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
TString xyzprint(Float_t X,Float_t Y,Float_t Z,Float_t dX,Float_t dY,Float_t dZ){
TString bzzz = "(X=";
bzzz += X;
bzzz += ", Y=";
bzzz += Y;
bzzz += ", Z=";
bzzz += Z;
bzzz += ", dX=";
bzzz += dX;
bzzz += ", dY=";
bzzz += dY;
bzzz += ", dZ=";
bzzz += dZ;
bzzz += ")";
return bzzz;
}
void tlcalc(int nmax=0,TString treename="PhaseSpace",TString rootfilename="PhaseSpace.root",double boxDX=8.,double boxDY=8.,double boxDZ=1.){
TChain *chain = new TChain(treename);
chain->Add(rootfilename);
double volume = 8 * boxDX * boxDY * boxDZ;
double LEsum = 0.;
// Char_t lastVolumeName[256];
char lastVolumeName[256];
Float_t X, Y, Z;
Float_t dX, dY, dZ;
Float_t Ekine;
chain->SetBranchAddress("ProductionVolume", &lastVolumeName);
chain->SetBranchAddress("Ekine",&Ekine);
chain->SetBranchAddress("X",&X);
chain->SetBranchAddress("Y",&Y);
chain->SetBranchAddress("Z",&Z);
chain->SetBranchAddress("dX",&dX);
chain->SetBranchAddress("dY",&dY);
chain->SetBranchAddress("dZ",&dZ);
Long64_t n = chain->GetEntries();
cout << "chain has " << n << " entries" << endl;
Long64_t i;
Long64_t nDirect = 0;
Long64_t nScattered = 0;
if (nmax == 0 || nmax>n) nmax = n;
for (i=0; i<nmax; i++){
if ( (i%100000) == 0 ){
int ileft = nmax - i;
cerr << "i=" << i << ", " << ileft << " still to go... " << endl;
}
chain->GetEntry(i);
if (X>=boxDX) continue;
if (Y>=boxDY) continue;
if (Z>=boxDZ) continue;
if (X<=-boxDX) continue;
if (Y<=-boxDY) continue;
if (Z<=-boxDZ) continue;
/**/
TString ee = lastVolumeName;
if (ee.Contains("target_syn") || ee.Contains("target_support_top") || ee.Contains("target_support_bottom")){
++nDirect;
} else {
++nScattered;
continue;
}
/**/
// ptp = punch through point (intersection of track with a wall of the box)
double ptp1[3];
double ptp2[3];
int nfound = 0;
double Zx1 = X + dX * (boxDZ-Z)/dZ;
double Zy1 = Y + dY * (boxDZ-Z)/dZ;
double Zz1 = boxDZ;
if ( (fabs(Zx1)<boxDX) && (fabs(Zy1)<boxDY) ){
nfound += 1;
if (nfound == 1){
ptp1[0] = Zx1;
ptp1[1] = Zy1;
ptp1[2] = Zz1;
// cerr << "Z1 is ptp1" << endl;
} else if (nfound == 2){
ptp2[0] = Zx1;
ptp2[1] = Zy1;
ptp2[2] = Zz1;
// cerr << "Z1 is ptp2" << endl;
} else {
cerr << "Z1 ERROR! xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
cerr << "Z1 ERROR! ptp1=(" << ptp1[0] << "," << ptp1[1] << "," << ptp1[2] << ")" << endl;
cerr << "Z1 ERROR! ptp2=(" << ptp2[0] << "," << ptp2[1] << "," << ptp2[2] << ")" << endl;
}
}
double Zx2 = X - dX * (Z+boxDZ)/dZ;
double Zy2 = Y - dY * (Z+boxDZ)/dZ;
double Zz2 = -boxDZ;
if ( (fabs(Zx2)<boxDX) && (fabs(Zy2)<boxDY) ){
nfound += 1;
if (nfound == 1){
ptp1[0] = Zx2;
ptp1[1] = Zy2;
ptp1[2] = Zz2;
// cerr << "Z2 is ptp1" << endl;
} else if (nfound == 2){
ptp2[0] = Zx2;
ptp2[1] = Zy2;
ptp2[2] = Zz2;
// cerr << "Z2 is ptp2" << endl;
} else {
cerr << "Z2 ERROR! xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
cerr << "Z2 ERROR! ptp1=(" << ptp1[0] << "," << ptp1[1] << "," << ptp1[2] << ")" << endl;
cerr << "Z2 ERROR! ptp2=(" << ptp2[0] << "," << ptp2[1] << "," << ptp2[2] << ")" << endl;
}
}
double Xx1 = boxDX;
double Xy1 = Y + dY * (boxDX-X)/dX;
double Xz1 = Z + dZ * (boxDX-X)/dX;
if ( (fabs(Xy1)<boxDY) && (fabs(Xz1)<boxDZ) ){
nfound += 1;
if (nfound == 1){
ptp1[0] = Xx1;
ptp1[1] = Xy1;
ptp1[2] = Xz1;
// cerr << "X1 is ptp1" << endl;
} else if (nfound == 2){
ptp2[0] = Xx1;
ptp2[1] = Xy1;
ptp2[2] = Xz1;
// cerr << "X1 is ptp2" << endl;
} else {
cerr << "X1 ERROR! i=" << i << " xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
cerr << "X1 ERROR! ptp1=(" << ptp1[0] << "," << ptp1[1] << "," << ptp1[2] << ")" << endl;
cerr << "X1 ERROR! ptp2=(" << ptp2[0] << "," << ptp2[1] << "," << ptp2[2] << ")" << endl;
}
}
double Xx2 = -boxDX;
double Xy2 = Y - dY * (X+boxDX)/dX;
double Xz2 = Z - dZ * (X+boxDX)/dX;
if ( (fabs(Xy2)<boxDY) && (fabs(Xz2)<boxDZ) ){
nfound += 1;
if (nfound == 1){
ptp1[0] = Xx2;
ptp1[1] = Xy2;
ptp1[2] = Xz2;
// cerr << "X2 is ptp1" << endl;
} else if (nfound == 2){
ptp2[0] = Xx2;
ptp2[1] = Xy2;
ptp2[2] = Xz2;
// cerr << "X2 is ptp2" << endl;
} else {
cerr << "X2 ERROR! i=" << i << " xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
cerr << "X2 ERROR! ptp1=(" << ptp1[0] << "," << ptp1[1] << "," << ptp1[2] << ")" << endl;
cerr << "X2 ERROR! ptp2=(" << ptp2[0] << "," << ptp2[1] << "," << ptp2[2] << ")" << endl;
}
}
double Yx1 = X + dX * (boxDY-Y)/dY;
double Yy1 = boxDY;
double Yz1 = Z + dZ * (boxDY-Y)/dY;
if ( (fabs(Yx1)<boxDX) && (fabs(Yz1)<boxDZ) ){
nfound += 1;
if (nfound == 1){
ptp1[0] = Yx1;
ptp1[1] = Yy1;
ptp1[2] = Yz1;
// cerr << "Y1 is ptp1" << endl;
} else if (nfound == 2){
ptp2[0] = Yx1;
ptp2[1] = Yy1;
ptp2[2] = Yz1;
// cerr << "Y1 is ptp2" << endl;
} else {
cerr << "Y1 ERROR! i=" << i << " xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
cerr << "Y1 ERROR! ptp1=(" << ptp1[0] << "," << ptp1[1] << "," << ptp1[2] << ")" << endl;
cerr << "Y1 ERROR! ptp2=(" << ptp2[0] << "," << ptp2[1] << "," << ptp2[2] << ")" << endl;
}
}
double Yx2 = X - dX * (Y+boxDY)/dY;
double Yy2 = -boxDY;
double Yz2 = Z - dZ * (Y+boxDY)/dY;
if ( (fabs(Yx2)<boxDX) && (fabs(Yz2)<boxDZ) ){
nfound += 1;
if (nfound == 1){
ptp1[0] = Yx2;
ptp1[1] = Yy2;
ptp1[2] = Yz2;
// cerr << "Y2 is ptp1" << endl;
} else if (nfound == 2){
ptp2[0] = Yx2;
ptp2[1] = Yy2;
ptp2[2] = Yz2;
// cerr << "Y2 is ptp2" << endl;
} else {
cerr << "Y2 ERROR! i=" << i << " xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
cerr << "Y2 ERROR! ptp1=(" << ptp1[0] << "," << ptp1[1] << "," << ptp1[2] << ")" << endl;
cerr << "Y2 ERROR! ptp2=(" << ptp2[0] << "," << ptp2[1] << "," << ptp2[2] << ")" << endl;
}
}
if (nfound != 2){
// cerr << "ERROR! xyz=" << xyzprint(X,Y,Z,dX,dY,dZ) << " ! nfound=" << nfound << endl;
continue;
}
double length2 = (ptp2[0]-ptp1[0])*(ptp2[0]-ptp1[0]);
length2 += (ptp2[1]-ptp1[1])*(ptp2[1]-ptp1[1]);
length2 += (ptp2[2]-ptp1[2])*(ptp2[2]-ptp1[2]);
double length = sqrt(length2);
LEsum += length * Ekine; // cm * MeV
}
cout << "nDirect=" << nDirect << " and nScattered=" << nScattered << endl;
cout << "energy fluence is " << LEsum << " / volume = " << LEsum/volume << " MeV/cm2" << endl;
}
More information about the Gate-users
mailing list