LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
braggPeak.C
Go to the documentation of this file.
1 #include "Riostream.h"
2 #include "TSystem.h"
3 #include "TInterpreter.h"
4 #include "TROOT.h"
5 #include "TApplication.h"
6 #include "TFile.h"
7 #include "TNtuple.h"
8 #include "TCanvas.h"
9 #include "TH1F.h"
10 #include "THStack.h"
11 #include "TCut.h"
12 #include "TString.h"
13 #include "TMath.h"
14 
15 /***************************
16  * Root script that produces comparison between hadrontherapy
17  * mc-simulation data and E.Haettner's experimental data.
18  *
19  * prompts user if to normalize to zero depth=1
20  *
21  * @author Gillis Danielsen
22  * **************************/
23 
24 
25 void braggPeak() {
26  TString normToZeroPos;
27  cout << "Normalize to first bin? (Y/N):";
28  cin >> normToZeroPos;
29 
30  TCanvas *c1 = new TCanvas("Bragg curve", "Bragg curve comparison");
31 
32  //TCanvas *c1 = new TCanvas("BraggPeaks", "Energy depositions along x-axis in phantom");
33  gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox
34 
35  gROOT->SetStyle("clearRetro");
36  //this will be used as base for pulling the experimental data
37  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
38  dir.ReplaceAll("basic.C","");
39  dir.ReplaceAll("/./","/");
40  ifstream in;
41  in.open(Form("experimentalData/iaeaBenchmark/400tdk.dat",dir.Data())); //also set Znum!
42  TString Znum = "2"; //set here what will be put in the selection for z, does not automatically change imoprted data.
43  Float_t f1,f2,f3,f4;
44  Int_t nlines = 0;
45  TFile *f = new TFile("fragmentAngularDistribution.root","RECREATE");
46  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","d:i:err1:err2");
47 
48  Char_t DATAFLAG[4];
49  Int_t NDATA;
50  Char_t n1[15], n2[15], n3[15], n4[15];
51  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 4'
52  in >> n1 >> n2 >> n3 >> n4; // Read column titles: 'Energy He B [...]'
53 
54  cout <<n1<<" "<<n2<<" "<<n3<<" "<<n4<<"\n";
55  while (1) {
56  in >> f1 >> f2 >> f3 >> f4;
57  if (!in.good()){
58  break;
59  }
60  if (nlines < 500 ) printf("%f %f %f %f\n",f1,f2,f3,f4);
61  ntuple->Fill(f1,f2,f3,f4);
62  nlines++;
63  }
64 
65  //Let's pull in the simulation data
66  TFile* simulation = TFile::Open("IAEA_braggPeak.root");
67  TH1F* simBragg = (TH1F*) simulation->Get("braggPeak");
68 
69  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
70  TNtuple *metadata = (TNtuple*) simulation->Get("metaData");
71  Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
72  metadata->SetBranchAddress("events",&events);
73  metadata->SetBranchAddress("waterThickness",&waterThickness);
74  metadata->SetBranchAddress("detectorDistance",&detectorDistance);
75  metadata->SetBranchAddress("beamEnergy",&beamEnergy);
76  metadata->SetBranchAddress("energyError",&energyError);
77  metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
78  metadata->GetEntry(0); //there is just one row to consider.
79 
80  //good to keep for ref. G4 might give weird units due to change.
81  metadata->Scan();
82 
83  simBragg->GetXaxis()->SetLimits(0.0, waterThickness);
84  simBragg->SetLineColor(kBlue);
85  simBragg->SetXTitle("Depth in phantom (cm)");
86  //simBragg->GetXaxis()->SetTitleOffset(1);
87  simBragg->GetYaxis()->SetTitleOffset(1.5);
88  std::cout << "Maximum (Bragg peak) for simulation data is at: " << simBragg->GetBinCenter(simBragg->GetMaximumBin()) + 0.478/2 + 0.027 + 0.073 << endl;
89  std::cout << "Bin width is " << simBragg->GetBinWidth(simBragg->GetMaximumBin()) << endl;
90 
91  TString scaleTuple;
92  if(normToZeroPos == "Y"){
93  Float_t normElement;
94  ntuple->SetBranchAddress("i",&normElement);
95  ntuple->GetEntry(0);
96  scaleTuple = Form("/%f", normElement);
97  simBragg->Scale(1.0/simBragg->GetBinContent(0));
98  simBragg->SetYTitle("Relative ionization");
99  }else{
100  simBragg->Scale(1.0/(100*events*simBragg->GetBinWidth(0))); // 100 for converting to MeV/m
101  simBragg->SetYTitle("Ionization (MeV/m)");
102  scaleTuple = "";
103  }
104 
105  //simBragg->ShowPeaks(); //Can be used to mark out bragg peaks specificly.
106  simBragg->Draw();
107  ntuple->SetMarkerColor(kRed);
108  ntuple->SetMarkerStyle(22);
109  std::cout << ntuple->GetEntries() << endl;
110  ntuple->Draw("i" + scaleTuple + ":d-(0.478+0.027+0.073)","","p,same"); // the minuses are the WE's, here only real water depth is plotted.
111  c1->SaveAs("braggPeakComparisonToData_norm_" + normToZeroPos + ".png");
112 }
Float_t f4
TTree * ntuple
Definition: macro.C:6
DecayPhysicsFactory f3
SynchrotronAndGN f2
TFile f
Definition: plotHisto.C:6
nlines
TCanvas * c1
Definition: plotHisto.C:7
void braggPeak()
Definition: braggPeak.C:25
EmPhysicsFactory f1
printf("%d Experimental points found\n", nlines)
ifstream in
Definition: comparison.C:7
TDirectory * dir
Definition: macro.C:5