Geant4_10
fragmentEnergy.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 "TCut.h"
11 #include "TString.h"
12 #include "TMath.h"
13 
20 void fragmentEnergy() {
24  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
25  dir.ReplaceAll("fragmentEnergy.C","");
26  dir.ReplaceAll("/./","/");
27 
28 gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox
29 
30 
31 TString macroPath(gROOT->GetMacroPath());
32 gROOT->SetMacroPath(macroPath + ":RootScripts/iaeaBenchmark");
33 //gROOT->LoadMacro("rootlogon.C");
34 //gROOT->SetStyle("clearRetro"); //For stylesheet
35 
36  TString pDepth;
37  cout << "Enter phantom depth (eg. 27.9, see experimentalData directory for choices): " << endl;
38  cout << "Entering 27.9 will make script look for IAEA_27.9.root:";
39  cin >> pDepth;
40 
41  TString simulationDataPath = "IAEA_" + pDepth + ".root";
42  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","Energy:He:B:H:Li:Be");
43  if(pDepth == "27.9"){
44  //there is only hand.read data for this depth, fixme, is ugly
45  ifstream in;
46  in.open(Form("experimentalData/iaeaBenchmark/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
47  Float_t f1,f2,f3, f4,f5,f6;
48  Int_t nlines = 0;
49  TFile *f = new TFile("fragmentEnergy.root","RECREATE");
50 
51  Char_t DATAFLAG[4];
52  Int_t NDATA;
53  Char_t n1[6], n2[2], n3[2], n4[2], n5[2], n6[2];
54  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
55  in >> n1 >> n2 >> n3 >> n4 >> n5 >> n6; // Read column titles: 'Energy He B [...]'
56 
57  cout <<n1<<" "<<n2<<" "<<n3<<" "<<n4<<" "<<n5<<" "<<n6<<"\n";
58  while (1) {
59  in >> f1 >> f2 >> f3 >>f4 >> f5 >> f6;
60  if (!in.good()) break;
61  if (nlines < 500 ) printf("%f %0.2f %0.2f %0.2f %0.2f %0.2f \n",f1,f2,f3,f4,f5,f6);
62  ntuple->Fill(f1,f2,f3,f4,f5,f6);
63  nlines++;
64  }
65  ntuple->SetMarkerStyle(5);
66  ntuple->Draw("He:Energy","","l");
67  ntuple->Draw("B:Energy","","l,Same");
68  ntuple->Draw("H:Energy","","l,Same");
69  ntuple->Draw("Li:Energy","","l,Same");
70  ntuple->Draw("Be:Energy","","l,Same");
71  printf(" found %d points\n",nlines);
72  }
73 
74  //Let's pull in the monte carlo simulation results
75  TCanvas *mc = new TCanvas("mc", "Simulation");
76  TFile *MCData = TFile::Open(simulationDataPath);
77  TH1F* MC_helium = (TH1F*)MCData->Get("heliumEnergyAfterPhantom");
78  TH1F* MC_hydrogen = (TH1F*)MCData->Get("hydrogenEnergyAfterPhantom");
79 //scale and plot
80  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
81 
82  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
83  TNtuple *metadata = (TNtuple*) MCData->Get("metaData");
84  Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
85  metadata->SetBranchAddress("events",&events);
86  metadata->SetBranchAddress("waterThickness",&waterThickness);
87  metadata->SetBranchAddress("detectorDistance",&detectorDistance);
88  metadata->SetBranchAddress("beamEnergy",&beamEnergy);
89  metadata->SetBranchAddress("energyError",&energyError);
90  metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
91  metadata->GetEntry(0); //there is just one row to consider.
92 
93 //analysis numbers based on metadata
94  Double_t scatteringDistance = detectorDistance - phantomCenterDistance;
95  Double_t detectorSideLength = 4; //hardcoded, we have a zero angle square detector
96 //good to keep for ref. G4 might give weird units due to change.
97 metadata->Scan();
98 
99 
100  Double_t binAmount = 50.0; //casting from int failed somehow, so in float temporarily, fixme
101  Double_t maxEnergy = 450.0;
102  Double_t binWidth = maxEnergy / binAmount;
103  TH1F *histH = new TH1F("histH", "Hydrogen", binAmount, 0.0, maxEnergy);
104  //HistH will be black
105  TH1F *histHe = new TH1F("histHe", "Helium", binAmount, 0.0, maxEnergy);
106  histHe->SetLineColor(kRed);
107  TH1F *histLi = new TH1F("histLi", "Lithium", binAmount, 0.0, maxEnergy);
108  histLi->SetLineColor(kBlue);
109  TH1F *histBe = new TH1F("histBe", "Beryllium", binAmount, 0.0, maxEnergy);
110  histBe->SetLineColor(kGreen);
111  TH1F *histB = new TH1F("histB", "Boron", binAmount, 0.0, maxEnergy);
112  histB->SetLineColor(kYellow);
113  TH1F *histC = new TH1F("histC", "Carbon", binAmount, 0.0, maxEnergy);
114 
115 TH1F* histPos = new TH1F("histPos", "check position",100,-2000,2000);
116 //Solid angle according to \Omega = 4 \arcsin \frac {\alpha\beta} {\sqrt{(4d^2+\alpha^2)(4d^2+\beta^2)}}
117 Double_t steradians = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
118  std::cout << "Detector seen at solid angle: " << steradians << endl;
119  TString normalization(Form("/%f", steradians*events*binWidth));
120 
121  fragments->SetLineColor(kRed);
122  fragments->SetMarkerStyle(22);
123 
124 
125 TString halfSideLengthString(Form("%f", detectorSideLength/2));
126 
127 
128  fragments->SetLineColor(kGreen);
129  fragments->Project("histHe", "energy", "(Z == 2 && energy > 45 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
130  fragments->Project("histB", "energy", "(Z == 5 && energy > 45 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
131  fragments->Project("histH", "energy", "(Z == 1 && energy > 45 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
132  fragments->Project("histLi", "energy", "(Z == 3 && energy > 45 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
133  fragments->Project("histBe", "energy", "(Z == 4 && energy > 45 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
134  fragments->Project("histB", "energy", "(Z == 5 && energy > 45 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
135 
136  histH->SetMaximum(0.27);
137 
138 
139 /*
140 * This is cludgy but the previous plots can't contain <45MeV stuff but hte following calcualtions need them
141 * Thus there is another projection done with different selections, this is error-prone,
142 * se to that changes are made in both places.
143 * */
144  fragments->Project("histHe", "energy", "(Z == 2 && energy > 0 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization);
145  fragments->Project("histB", "energy", "(Z == 5 && energy > 0 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization, "same");
146  fragments->Project("histH", "energy", "(Z == 1 && energy > 0 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization, "same");
147  fragments->Project("histLi", "energy", "(Z == 3 && energy > 0 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization, "same");
148  fragments->Project("histBe", "energy", "(Z == 4 && energy > 0 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization, "same");
149  fragments->Project("histB", "energy", "(Z == 5 && energy > 0 && abs(posY) < " + halfSideLengthString + " && abs(posZ) < " + halfSideLengthString + " )" + normalization, "same");
150 
151  histH->SetXTitle("Energy per nucleon (MeV/u)");
152  histH->SetYTitle("(N/N0) [1/(MeV*sr)]");
153 
154  TCanvas *c3 = new TCanvas("histograms", "Energy distribution of charged fragments");
155  Int_t nEve = events; //redundant
156 
157  cout <<"H : " << histH->GetEntries() / nEve << endl;
158  histH->Draw("");
159 
160  cout <<"He : " << histHe->GetEntries() / nEve << endl;
161  histHe->Draw("same");
162 
163  cout <<"Li : " << histLi->GetEntries() / nEve << endl;
164  histLi->Draw("same");
165 
166 
167  cout <<"Be : " << histBe->GetEntries() / nEve << endl;
168  histBe->Draw("same");
169 
170  cout <<"B : " << histB->GetEntries() / nEve << endl;
171  histB->Draw("same");
172 
173  // Legends for the data
174  leg = new TLegend(0.8,0.5,1,1); //coordinates are fractions
175  leg->SetHeader("Fragments");
176  leg->AddEntry(histH,"Hydrogen","l");
177  leg->AddEntry(histHe,"Helium","l");
178  leg->AddEntry(histLi,"Lithium","l");
179  leg->AddEntry(histBe,"Beryllium","l");
180  leg->AddEntry(histB,"Boron","l");
181  leg->Draw();
182 
183  ntuple->SetMarkerStyle(22);
184  ntuple->Draw("H:Energy","","p,same");
185  ntuple->SetMarkerColor(kRed);
186  ntuple->Draw("He:Energy","","p,same");
187  ntuple->SetMarkerColor(kBlue);
188  ntuple->Draw("Li:Energy","","p,same");
189  ntuple->SetMarkerColor(kGreen);
190  ntuple->Draw("Be:Energy","","p,same");
191  ntuple->SetMarkerColor(kYellow);
192  ntuple->Draw("B:Energy","","p,same");
193 
194  c3->SaveAs("fragmentEnergyDistr.png");
195 
196  in.close();
197 
198  f->Write();
199 }
200 
in
ifstream in
Definition: comparison.C:7
f4
Float_t f4
Definition: comparison_ascii.C:56
UUtils::ASin
double ASin(double)
Definition: UUtils.hh:211
ntuple
TTree * ntuple
Definition: macro.C:6
f
TFile f
Definition: plotHisto.C:6
f1
Float_t f1
Definition: comparison_ascii.C:56
f2
Float_t f2
Definition: comparison_ascii.C:56
G4INCL::Int_t
G4int Int_t
Definition: G4INCLEventInfo.hh:59
f3
Float_t f3
Definition: comparison_ascii.C:56
fragmentEnergy
void fragmentEnergy()
Definition: fragmentEnergy.C:20
G4INCL::Float_t
G4float Float_t
Definition: G4INCLEventInfo.hh:61
printf
printf("%d Experimental points found\n", nlines)
G4INCL::Double_t
G4double Double_t
Definition: G4INCLEventInfo.hh:62
leg
leg
Definition: comparison.C:61
nlines
nlines
Definition: comparison_ascii.C:57
dir
TDirectory * dir
Definition: macro.C:5
Generated on Sat Dec 14 2013 14:32:59 for Geant4_10 by   doxygen 1.8.5