View source for SPIM Brem Lab Instructions

The objective of this lab is to alter the ExN02PhysicsList.cc program to
ONLY 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/Classes/NucSim/Day8/Mondelaers_XXInt.Linac_Conf._Brem_E-spectrum.pdf

[[File:Mondelaers_XXIntLinacConf.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 GEANT4 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()


<pre>
  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;
  }

</pre>
make sure the above code isn't embedded in another if statement

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


Change the target matter variable to tantalum (Tnt)

3.) After you check that things are working right  (check ExN02PhysicsList.cc to be sure that only the Bremsstrahlung process is turned on for the electrons).

Run 10000 events at 15 MeV

create the file run1.mac it has the following commands

/gun/particle e-
/gun/energy 15 MeV
/event/verbose 0
/tracking/verbose 1
/run/beamOn 10000
exit

now you can run the simulation in "batch mode" ie without visualization and re-direct teh output to a file

ie: > exampleN02 run1.mac > /dev/null

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;
 }


Copy the above program into a file called Brem.C and save it into the same subdirectory as the output file you created with the simulation.

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 don't want things deleted).

to run this program run root from the same subdirectory as Brem.C and just type ".x Brem.C" at the root prompt

(to run root set the environmental variable ROOTSYS to point to the root subdirectory and then type $ROOTSYS/bin/root)

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:

Run root and give it the filename of the root file (Brem.root) on the command line (or you could type new TBrowser(): and load it from the GUI).

ie: root Brem.root

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.14","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 little tricky to have two root files open at the same 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 with 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.


commands to generate a pdf file from the template are

 latex filename.tex

 dvipdf filename

you can look at the text using 

  xdvi filename.dvi

you can convert the dvi file to a PDF file using the commands below

  dvips filename.dvi

  ps2pdf filename.ps

You can look at the pdf file using the command

 xpdf filename.pdf

or 

 evince filename.pdf




---
  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
referenced experiment, a description of the simulation, and
analysis of GEANT4 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.

---

[http://wiki.iac.isu.edu/index.php/HomeWork_Simulations_of_Particle_Interactions_with_Matter Back to Homework Assigments]