Difference between revisions of "SPIM Brem Lab Instructions"

From New IAC Wiki
Jump to navigation Jump to search
 
(36 intermediate revisions by 2 users not shown)
Line 1: Line 1:
The objective of this lab is to alter the BetheBloch program to
+
The objective of this lab is to alter the ExN02PhysicsList.cc program to
include the Bremsstrahlung Physics process in order to determine
+
ONLY include the Bremsstrahlung Physics process in order to determine
 
the photon energy distribution due to bremsstrahlung and compare
 
the photon energy distribution due to bremsstrahlung and compare
 
that distribution with experiment.
 
that distribution with experiment.
Line 6: Line 6:
 
Figure #2 in  
 
Figure #2 in  
  
http://physics.isu.edu/~tforest/Classes/NucSim/Day8/Mondelaers_XXInt.Linac_Conf._Brem_E-spectrum.pdf
+
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
 
shows the photon energy distribution when 15 MeV electrons
 
impinge on a 4mm thick target of Graphite (C12) and Tantalum.   
 
impinge on a 4mm thick target of Graphite (C12) and Tantalum.   
  
Let's alter the BetheBloch program to output the photon kinetic
+
Let's alter the GEANT4 program to output the photon kinetic
 
energy, position, and momentum.
 
energy, position, and momentum.
  
Line 19: Line 22:
  
  
 
+
<pre>
 
   if( fTrack->GetDefinition()->GetPDGEncoding()==22 &&  fTrack->GetCurrentStepNumber()==1 && fTrack->GetKineticEnergy()>0)
 
   if( fTrack->GetDefinition()->GetPDGEncoding()==22 &&  fTrack->GetCurrentStepNumber()==1 && fTrack->GetKineticEnergy()>0)
 
     {
 
     {
 
 
       //      G4cout  << "  Photon Energy "   
 
       //      G4cout  << "  Photon Energy "   
 
       outfile <<  fTrack->GetKineticEnergy() << "    "
 
       outfile <<  fTrack->GetKineticEnergy() << "    "
      << fTrack->GetPosition().x()<< "    "
+
      << fTrack->GetPosition().x()<< "    "
      << fTrack->GetPosition().y()<< "    "
+
      << fTrack->GetPosition().y()<< "    "
      << fTrack->GetPosition().z()<< "    "
+
      << fTrack->GetPosition().z()<< "    "
      << fTrack->GetMomentum().x() << "    "
+
      << fTrack->GetMomentum().x() << "    "
      << fTrack->GetMomentum().y() << "    "
+
      << fTrack->GetMomentum().y() << "    "
      << fTrack->GetMomentum().z() << "    "
+
      << fTrack->GetMomentum().z() << "    "
      << G4endl;
+
      << G4endl;
 
   }
 
   }
  
 
+
</pre>
 
make sure the above code isn't embedded in another if statement
 
make sure the above code isn't embedded in another if statement
  
Line 50: Line 52:
 
Change the target matter variable to tantalum (Tnt)
 
Change the target matter variable to tantalum (Tnt)
  
3.) After you check that things are working right  
+
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
 
Run 10000 events at 15 MeV
  
create the file run1.mac it has the following commands
+
create the file run1.mac with the following commands
  
 +
<pre>
 
/gun/particle e-
 
/gun/particle e-
 
/gun/energy 15 MeV
 
/gun/energy 15 MeV
Line 61: Line 65:
 
/run/beamOn 10000
 
/run/beamOn 10000
 
exit
 
exit
 +
</pre>
  
now you can run the simulation in "batch mode" ie without visualization and re-direct teh output to a file
+
now you can run the simulation in "batch mode" ie without visualization and re-direct the output to a file
  
ie: > exampleN02 run1.mac > /dev/null
+
ie: > ./exampleN02 run1.mac > /dev/null &
  
4.) You should have a file called SigmaR.txt which has entries
+
4.) You may have a file called Range.date which has entries
 
that look like
 
that look like
  
 
+
<pre>
 
0.082831    -0.102156    -0.155741    -800    -0.00347244    -0.00405495    0.0826588     
 
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.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     
 
0.332485    -172.26    -254.795    800    -0.0621649    -0.0604049    0.320987     
 +
</pre>
  
 
+
create a file for C12 and Tantalum (make may need to rename
create a file for C12 and Tantalum (make sure to rename
 
 
SigmaR.txt to something else or it will be written over when you
 
SigmaR.txt to something else or it will be written over when you
 
run the program again)
 
run the program again)
Line 83: Line 88:
 
root tree
 
root tree
  
void Brem() {
+
<pre>
 
+
 
  struct evt_t {
+
void Brem() {
 +
    struct evt_t {
 
     Int_t event;
 
     Int_t event;
 
     Float_t KE, pos[3],mom[3];
 
     Float_t KE, pos[3],mom[3];
 
   };
 
   };
 
 
   ifstream in;
 
   ifstream in;
 
   in.open("Brem.txt");
 
   in.open("Brem.txt");
 
 
   evt_t evt;
 
   evt_t evt;
 
   Int_t nlines=0;
 
   Int_t nlines=0;
Line 100: Line 104:
 
   while(in.good()){
 
   while(in.good()){
 
     evt.event=nlines;
 
     evt.event=nlines;
     in >> evt.KE >> evt.pos[0] >> evt.pos[1] >> evt.pos[2]
+
     in >> evt.KE >> evt.pos[0] >> evt.pos[1] >> evt.pos[2]   >> evt.mom[0]  
      >> evt.mom[0] >> evt.mom[1] >> evt.mom[2];
+
>> evt.mom[1] >> evt.mom[2];
 
     /*
 
     /*
 
+
       printf( " %d  %f %f %f %f %f %f %f\n", evt.event, evt.KE, evt.pos[0],  
       printf( " %d  %f %f %f %f %f %f %f\n",
+
evt.pos[1], evt.pos[2], evt.mom[0], evt.mom[1], evt.mom[2] );
    evt.event, evt.KE, evt.pos[0], evt.pos[1], evt.pos[2],  
 
  evt.mom[0], evt.mom[1], evt.mom[2] );
 
 
 
 
     */
 
     */
 
     nlines++;
 
     nlines++;
 
     tree->Fill();
 
     tree->Fill();
 
 
   }
 
   }
 
   tree->Print();
 
   tree->Print();
   tree->Write();
+
   tree->Write();
 
 
 
   in.close();
 
   in.close();
 
   delete tree;
 
   delete tree;
 
   delete f;
 
   delete f;
}
+
}
  
 +
</pre>
  
 
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.
 
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.
Line 126: Line 126:
 
Notice it expects the data file to be called "Brem.txt" and it
 
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
 
creates the root file "Brem.root" (you will need to rename files
if you dont want things deleted).
+
if you don't want things overwritten).
  
 
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 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
 
A root file is created called "Brem.root".  ".q" root and rename
Line 136: Line 138:
 
6.)  Now analyze the root file:
 
6.)  Now analyze the root file:
  
a.) first, from within root, define a variable to point to the
+
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).
root tree and assign it to the data.
 
  
 +
ie: root Brem.root
  
root[0] TTree *tree=Brem;
+
a.) first, from within root, clock on the Brem.root file name under the "ROOT files" subdirectory
  
 
Now try to plot the scattering angle of the outgoing photon with
 
Now try to plot the scattering angle of the outgoing photon with
 
respect to the Z-axis
 
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");
+
root[0]Brem->Draw("asin((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.14","pz>0");
  
  
Line 158: Line 160:
 
tree and the Tantalum root Tree.
 
tree and the Tantalum root Tree.
  
It is a litle tricky to have two root files open at the sam time
+
It is a little tricky to have two root files open at the same time
 
for plotting but it can be done.
 
for plotting but it can be done.
  
Line 177: Line 179:
 
root[0] TTree *tree=Brem;  
 
root[0] TTree *tree=Brem;  
  
you will redefine the tree to be associated ith file2.root.
+
you will redefine the tree to be associated with file2.root.
  
 
You can move between the two files by clicking on the name and
 
You can move between the two files by clicking on the name and
Line 195: Line 197:
 
now fill the histgram with the data
 
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");
+
root[2]Brem->Draw("asin((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.14 >> BremAngle_C12","pz>0");
  
 
you will see a histogram save under the file "C12.root"
 
you will see a histogram save under the file "C12.root"
Line 212: Line 214:
 
now fill the histgram with the data
 
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");
+
root[2]Brem->Draw("asin((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.14 >> BremAngle_Tant","pz>0");
  
 
you will see a histogram save under the file "Tantalum.root"
 
you will see a histogram save under the file "Tantalum.root"
Line 220: Line 222:
 
drawing the histograms
 
drawing the histograms
  
click on "C12.rot" file and then do  
+
click on "C12.root" file and then do  
  
 
  BremAngle_C12->Draw();
 
  BremAngle_C12->Draw();
Line 268: Line 270:
 
with a graph included.
 
with a graph included.
  
******
+
 
 +
commands to generate a pdf file from the template are
 +
 
 +
pdflatex filename.tex
 +
 
 +
evince filename.pdf
 +
 
 +
 
 +
 
 +
 
 +
---
 
   Note: My grading scheme is as follows
 
   Note: My grading scheme is as follows
  
Line 277: Line 289:
  
 
10 points /10 if you write a complete description summarizing the
 
10 points /10 if you write a complete description summarizing the
references experiment and a description of the simulation and
+
referenced experiment, a description of the simulation, and
analysis of GENT4 output.  The grade is not based on document
+
analysis of GEANT4 output.  The grade is not based on document
 
length but document completeness.  You can provide enough details
 
length but document completeness.  You can provide enough details
 
about the experiment in 2 paragraphs to be complete.  The GEANT
 
about the experiment in 2 paragraphs to be complete.  The GEANT
Line 286: Line 298:
 
describe the analysis of the GEANT4 output.
 
describe the analysis of the GEANT4 output.
  
*********
+
---
  
[http://www.iac.isu.edu/mediawiki/index.php/HomeWork_Simulations_of_Particle_Interactions_with_Matter Back to Homework Assigments]
+
[http://wiki.iac.isu.edu/index.php/HomeWork_Simulations_of_Particle_Interactions_with_Matter Back to Homework Assigments]

Latest revision as of 21:24, 30 September 2015

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()


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

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 with 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 the output to a file

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

4.) You may have a file called Range.date 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 may need 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 overwritten).

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, clock on the Brem.root file name under the "ROOT files" subdirectory

Now try to plot the scattering angle of the outgoing photon with respect to the Z-axis

root[0]Brem->Draw("asin((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("asin((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.14 >> 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("asin((px*px+py*py)/(px*px+py*py+pz*pz))*180/3.14 >> 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.root" 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

pdflatex filename.tex
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.

---

Back to Homework Assigments