SPIM Brem Lab Instructions
The objective of this lab is to alter the BetheBloch program to include the Bremsstrahlung Physics process in order to determine the photon energy distribution due to bremsstrahlung and compare that distribution with experiment.
Figure #2 in
http://physics.isu.edu/~tforest/NucSim/Day8/Mondelaers_XXInt.Linac_Conf._Brem_E-spectrum.pdf shows the photon energy distribution when 15 MeV electrons impinge on a 4mm thick target of Graphite (C12) and Tantalum.
Let's alter the BetheBloch program to output the photon kinetic energy, position, and momentum.
1.) In the file src/ExN02SteppingVerbose.cc add the code below to the function
void ExN02SteppingVerbose::StepInfo()
if( fTrack->GetDefinition()->GetPDGEncoding()==22 && fTrack->GetCurrentStepNumber()==1 && fTrack->GetKineticEnergy()>0) {
// G4cout << " Photon Energy " outfile << fTrack->GetKineticEnergy() << " "
<< fTrack->GetPosition().x()<< " " << fTrack->GetPosition().y()<< " " << fTrack->GetPosition().z()<< " " << fTrack->GetMomentum().x() << " " << fTrack->GetMomentum().y() << " " << fTrack->GetMomentum().z() << " " << G4endl;
}
2.) Now edit src/ExN02DetectorConstruction.cc
and add the material Tantalum to the list
G4Material* Tnt = new G4Material("Tantalum", z=73., a= 180.9479*g/mole, density= 16.654*g/cm3);
set target length to 4 mm
fTargetLength = 0.4 * cm; // Full length of Target
3.) After you chack that things are working right
Run 10000 events at 15 MeV
the end of vis.mac should look like:
/gun/particle e- /gun/energy 15 MeV /event/verbose 0 /tracking/verbose 1 /run/beamOn 10000 exit
4.) You should have a file called SigmaR.txt which has entries that look like
0.082831 -0.102156 -0.155741 -800 -0.00347244 -0.00405495 0.0826588
0.346969 0.499903 -0.474851 -800 0.0683589 -0.0569303 0.335371
0.332485 -172.26 -254.795 800 -0.0621649 -0.0604049 0.320987
create a file for C12 and Tantalum (make sure to rename
SigmaR.txt to something else or it will be written over when you
run the program again)
5.) I wrote the following ROOT macro to read in the data into a
root tree
void Brem() {
struct evt_t { Int_t event; Float_t KE, pos[3],mom[3]; };
ifstream in; in.open("Brem.txt");
evt_t evt; Int_t nlines=0; TFile *f = new TFile("Brem.root","RECREATE"); TTree *tree = new TTree("Brem","Brem data from ascii file"); tree->Branch("evt",&evt.event,"event/I:ke/F:posx:posy:posz:px:py:pz"); while(in.good()){ evt.event=nlines; in >> evt.KE >> evt.pos[0] >> evt.pos[1] >> evt.pos[2] >> evt.mom[0] >> evt.mom[1] >> evt.mom[2]; /*
printf( " %d %f %f %f %f %f %f %f\n",
evt.event, evt.KE, evt.pos[0], evt.pos[1], evt.pos[2], evt.mom[0], evt.mom[1], evt.mom[2] );
*/ nlines++; tree->Fill();
} tree->Print(); tree->Write(); in.close(); delete tree; delete f;
}
Notice it expects the data file to be called "Brem.txt" and it
creates the root file "Brem.root" (you will need to rename files
if you dont want things deleted).
to run this program just type ".x Brem.C" at the root prompt
A root file is created called "Brem.root". ".q" root and rename the file so you won't write over it the next time you run the root program.
6.) Now analyze the root file:
a.) first, from within root, define a variable to point to the root tree and assign it to the data.
root[0] TTree *tree=Brem;
Now try to plot the scattering angle of the outgoing photon with respect to the Z-axis
root[1]Brem->Draw("atan((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.18","pz>0");
Notice "pz>0" above is a cut that only looks at forward going
photons (no backward scattered photons are plotted).
The experiment reports that the lower Z target creates photons which have a range of angles that are smaller than the high Z target.
You can check this by comparing the above plot for the C12 root tree and the Tantalum root Tree.
It is a litle tricky to have two root files open at the sam time for plotting but it can be done.
Use the Browser to open both files ("new TBrowser" opens the browser window and clicking on file names opens the file)
When you click on the file name listed under "ROOT files" ROOT will direct all commands to that file. So if you click on file1.root and then do
root[0] TTree *tree=Brem;
You will now analyze the entries in file1.root and histograms will be save there.
If you click on "file2.root" and again execute
root[0] TTree *tree=Brem;
you will redefine the tree to be associated ith file2.root.
You can move between the two files by clicking on the name and redefining th tree.
Suppose you want to create a histogram now of the photon scattering angles for C12 and Tantalum.
click on the C12.root file and execute the command
root[0] TTree *tree=Brem;
then define a 1-D histogram which spans the angular range.
root[1]TH1F *C12hist=new TH1F("BremAngle_C12","BremAngle_C12",90,0,90);
now fill the histgram with the data
root[2]Brem->Draw("atan((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.18 >> BremAngle_C12","pz>0");
you will see a histogram save under the file "C12.root"
Do the same thing for Tantalum
first click on the "Tantalum.root" file name then
root[0] TTree *tree=Brem;
then define a 1-D histogram which spans the angular range.
root[1]TH1F *Tanthist=new TH1F("BremAngle_Tant","BremAngle_Tant",90,0,90);
now fill the histgram with the data
root[2]Brem->Draw("atan((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.18 >> BremAngle_Tant","pz>0");
you will see a histogram save under the file "Tantalum.root"
You can overlay the two histograms using the "same" switch when
drawing the histograms
click on "C12.rot" file and then do
BremAngle_C12->Draw();
now click on "Tantalum.root" and do
BremAngle_Tant->Draw("same");
The two distributions are now on the same plot for comparison. Notice that the angular distribution for the tantalum target is wider (goes to higher angles) than the C12 distribution, as suggested by the article.
b.) Now we want to create the Histograms for Figure 2 in the article. I was unable to determine what photon angular range was subtended by the detector in the experiment. You can clearly see that if you cut on the photon angle, the energy distribution of the Tantalum changes.
click on the "Tantalum.root" file and redefine th tree pointer
root[0]TTree *tree=Brem
now create a histogram to store the energy distribution:
root[1]TH1F *TantKEhist=new TH1F("BremKE_Tant","BremKE_Tant",500,0,1);
now fill the histogram with the KE using different angle cuts and watch how the distribution changes.
root[2]Brem->Draw("ke >> BremKE_Tant","pz>0 && (px*px+py*py)/(px*px+py*py+pz*pz) < 0.5");
or
root[2]Brem->Draw("ke >> BremKE_Tant","pz>0 && (px*px+py*py)/(px*px+py*py+pz*pz) < 0.1");
Your goal for this lab is to change the angle cut until you get
something close to Figure 2 in the paper (this won't be a precise
method for determining the angle only a qualitative one).
You will use the TeX template I gave you to write up your result
with a graph included.
Note: My grading scheme is as follows
8 points/10 if you supply a pdf file with an energy specta from GEANT4 for C12 and Tantalum
9 points /10 if you include a spectrum of photon scattering angles
10 points /10 if you write a complete description summarizing the references experiment and a description of the simulation and analysis of GENT4 output. The grade is not based on document length but document completeness. You can provide enough details about the experiment in 2 paragraphs to be complete. The GEANT simulation may need a few more paragraphs in which you focus on detector geometry, physicslists, and the writing of tracking variables to a file. I imagine 4 paragraphs could adequately describe the analysis of the GEANT4 output.