SPIM Brem Lab Instructions

From New IAC Wiki
Revision as of 18:42, 22 October 2009 by Oborn (talk | contribs)
Jump to navigation Jump to search

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

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 BetheBlochPhysicsLists.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

the above 2 commands should first create a device independent file for your document and the second command converts it to a pdf file. 



 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.

Back to Homework Assigments

Retrieved from "https://wiki.iac.isu.edu/index.php?title=SPIM_Brem_Lab_Instructions&oldid=23327"