3 #include "TInterpreter.h"
5 #include "TApplication.h"
17 gStyle->SetOptStat(0000000000);
21 TString pDepth, fragment, Znum, normToOneAtZeroAngle;
22 cout <<
"Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
24 TString simulationDataPath =
"IAEA_" + pDepth +
".root";
25 TString
dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
26 dir.ReplaceAll(
"basic.C",
"");
27 dir.ReplaceAll(
"/./",
"/");
29 in.open(Form(
"experimentalData/iaeaBenchmark/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
32 TFile *
f =
new TFile(
"fragmentEnergyWithAngularDistribution.root",
"RECREATE");
33 TNtuple *
ntuple =
new TNtuple(
"ntuple",
"Data from ascii file",
"Energy:He:B:H:Li:Be");
37 Char_t n1[6], n2[2], n3[2], n4[2], n5[2], n6[2];
38 in >> DATAFLAG >> NDATA ;
39 in >> n1 >> n2 >> n3 >> n4 >> n5 >> n6;
41 cout <<n1<<
" "<<n2<<
" "<<n3<<
" "<<n4<<
" "<<n5<<
" "<<n6<<
"\n";
43 in >> f1 >> f2 >> f3 >>f4 >> f5 >> f6;
44 if (!in.good())
break;
45 if (nlines < 500 )
printf(
"%f %0.2f %0.2f %0.2f %0.2f %0.2f \n",f1,f2,f3,f4,f5,f6);
46 ntuple->Fill(f1,f2,f3,f4,f5,f6);
52 TFile *MCData = TFile::Open(
"IAEA_200000.root");
53 TNtuple *fragments = (TNtuple*) MCData->Get(
"fragmentNtuple");
56 TNtuple *metadata = (TNtuple*) MCData->Get(
"metaData");
57 Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
58 metadata->SetBranchAddress(
"events",&events);
59 metadata->SetBranchAddress(
"waterThickness",&waterThickness);
60 metadata->SetBranchAddress(
"detectorDistance",&detectorDistance);
61 metadata->SetBranchAddress(
"beamEnergy",&beamEnergy);
62 metadata->SetBranchAddress(
"energyError",&energyError);
63 metadata->SetBranchAddress(
"phantomCenterDistance",&phantomCenterDistance);
64 metadata->GetEntry(0);
68 std::cout <<
"Recieved metadata-row: " << events <<
" " << detectorDistance <<
" " << waterThickness <<
" " << beamEnergy <<
" " << energyError <<
" " << phantomCenterDistance;
74 Double_t binWidth = maxEnergy / binAmount;
75 TH1F *
hist1 =
new TH1F(
"hist1",
"", binAmount, 0.0, maxEnergy);
76 TH1F *
hist2 =
new TH1F(
"hist2",
"", binAmount, 0.0, maxEnergy);
77 TH1F *
hist3 =
new TH1F(
"hist3",
"", binAmount, 0.0, maxEnergy);
78 TH1F *
hist4 =
new TH1F(
"hist4",
"", binAmount, 0.0, maxEnergy);
79 TH1F *
hist5 =
new TH1F(
"hist5",
"", binAmount, 0.0, maxEnergy);
80 TH1F *
hist6 =
new TH1F(
"hist6",
"", binAmount, 0.0, maxEnergy);
81 TH1F *
hist7 =
new TH1F(
"hist7",
"", binAmount, 0.0, maxEnergy);
82 TH1F *
hist8 =
new TH1F(
"hist8",
"", binAmount, 0.0, maxEnergy);
83 TH1F *
hist9 =
new TH1F(
"hist9",
"", binAmount, 0.0, maxEnergy);
84 for(
int k = 1; k <= 6; k++){
85 TString Znum = Form(
"%i", k);
86 hist1->SetTitle(
"Z=" + Znum);
91 Double_t scatteringDistance = detectorDistance - phantomCenterDistance;
94 Double_t r, rMin, rMax, deltaOmega, normFloat;
95 TString rMinString, rMaxString, normString;
98 TCanvas *c3 =
new TCanvas(
"histograms",
"Distribution (at different angles)");
101 std::cout <<
"The following numbers also make it possible to make number of fragments comparison to the graph in A1 of E.Haettner\n";
102 for(
Double_t j = 0.0; j <= 8.0; j=j+1.0){
104 degrees = j * TMath::DegToRad();
107 r = scatteringDistance * TMath::Tan(degrees);
111 Double_t deltaPhi =
TMath::ATan((TMath::Cos(degrees)*detectorSideLength)/(2*scatteringDistance));
112 rMin = TMath::Max(0.0,r - (detectorSideLength/(2*TMath::Cos(degrees))));
113 rMax = rMin + ((detectorSideLength*TMath::Sin(degrees))/TMath::Tan((TMath::Pi()/2) - degrees - deltaPhi)) + (detectorSideLength*TMath::Cos(degrees));
114 rMinString = Form(
"%f", rMin);
115 rMaxString = Form(
"%f", rMax);
117 deltaPhi = degrees -
TMath::ATan(TMath::Tan(degrees) - detectorSideLength/(2*scatteringDistance));
119 deltaOmega = 2*TMath::Pi()*(TMath::Cos(TMath::Max(0.0,degrees-deltaPhi)) - TMath::Cos(degrees+deltaPhi));
121 deltaOmega = 4 *
TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
123 normFloat = deltaOmega * events * binWidth;
124 normString = Form(
"/%f", normFloat);
127 histName = Form(
"hist%i", i);
129 fragments->Project(histName,
"energy",
"(Z == " + Znum +
" && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString +
"&& sqrt(posY*posY + posZ*posZ) > " + rMinString +
")" + normString);
131 fragments->Project(histName,
"energy",
"(Z == " + Znum +
" && energy > 0 && posZ < " + rMaxString +
"&& posY < " + rMaxString +
" && posY > 0 && posZ > 0)" + normString);
133 int numEntries = fragments->GetEntries(
"(Z == " + Znum +
" && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString +
"&& sqrt(posY*posY + posZ*posZ) > " + rMinString +
")");
134 std::cout <<
"\nj: "<< numEntries/(deltaOmega * events) <<
" entries for " << j;
140 hist1->SetLineColor(kBlue);
143 hist2->SetLineColor(kGreen);
146 hist3->SetLineColor(kRed);
152 hist5->SetLineColor(kGreen + 3);
158 hist7->SetLineColor(kRed);
168 leg =
new TLegend(0.9,0.7,1,1);
169 leg->SetHeader(
"Angles");
170 leg->AddEntry(hist1,
"0",
"l");
171 leg->AddEntry(hist2,
"1",
"l");
172 leg->AddEntry(hist3,
"2",
"l");
173 leg->AddEntry(hist5,
"4",
"l");
174 leg->AddEntry(hist7,
"6",
"l");
177 c3->SaveAs(
"AEDistrib" + Znum +
".png");
void fragmentEnergyDistributionDifferentAngles()
printf("%d Experimental points found\n", nlines)