3 #include "TInterpreter.h"
5 #include "TApplication.h"
28 gStyle->SetOptStat(0000000000);
29 gROOT->SetStyle(
"clearRetro");
31 TString
dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
32 dir.ReplaceAll(
"basic.C",
"");
33 dir.ReplaceAll(
"/./",
"/");
36 TString pDepth, fragment, Znum;
37 cout <<
"Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
39 cout <<
"Enter fragment Z-number (eg. 1): ";
41 cout <<
"Enter fragment name (Znum 1 -> H,Znum 2->He...): ";
44 TString experimentalDataPath =
"experimentalData/iaeaBenchmark/angularDistributions/" + pDepth +
"/" + fragment +
"" + pDepth +
".dat";
45 TString simulationDataPath =
"IAEA_" + pDepth +
".root";
47 in.open(experimentalDataPath);
50 TFile *
f =
new TFile(
"fragmentAngularDistribution.root",
"RECREATE");
51 TNtuple *
ntuple =
new TNtuple(
"ntuple",
"Data from ascii file",
"x:y");
55 Char_t n1[15], n2[15];
56 in >> DATAFLAG >> NDATA ;
59 cout <<n1<<
" "<<n2<<
"\n";
62 if (!in.good())
break;
63 if (nlines < 500 )
printf(
"%f %f\n",f1,f2);
67 printf(
" found %d points\n",nlines);
70 TFile *simulation = TFile::Open(simulationDataPath);
71 TNtuple *fragments = (TNtuple*) simulation->Get(
"fragmentNtuple");
74 TNtuple *metadata = (TNtuple*) simulation->Get(
"metaData");
75 Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
76 metadata->SetBranchAddress(
"events",&events);
77 metadata->SetBranchAddress(
"waterThickness",&waterThickness);
78 metadata->SetBranchAddress(
"detectorDistance",&detectorDistance);
79 metadata->SetBranchAddress(
"beamEnergy",&beamEnergy);
80 metadata->SetBranchAddress(
"energyError",&energyError);
81 metadata->SetBranchAddress(
"phantomCenterDistance",&phantomCenterDistance);
82 metadata->GetEntry(0);
88 TString Znum = Form(
"%i", 1);
103 TH1F *
hist1 =
new TH1F(
"hist1",
"Fragment angular distribution", binAmount, 0, maxDegrees);
105 TH1F *symmetricHist =
new TH1F(
"symmetricHist",
"Fragment angular distr.", 2*binAmount, -1*maxDegrees, maxDegrees);
110 Double_t scatteringDistance = detectorDistance - phantomCenterDistance;
111 TString sdstring = Form(
"%.4f", scatteringDistance);
113 TCanvas *
c1 =
new TCanvas(
"histograms",
"Angular distribution at certain phantom thickness, sd" + sdstring);
120 fragments->Project(
"hist1",
"57.29577*atan((posZ^2+posY^2)/" + sdstring +
")",
"(Z == " + Znum +
")");
126 std::cout <<
"bin-number" <<
"\t" <<
"value" <<
"\t" <<
"solid angle segment" << endl;
127 for(
int bin = 0;
bin <= hist1->GetNbinsX();
bin++){
128 value = hist1->GetBinContent(
bin);
129 width = hist1->GetBinWidth(
bin);
130 degrees = hist1->GetBinCenter(
bin);
132 binNormalization = 2*TMath::Pi()*(TMath::Cos(TMath::DegToRad()*(degrees-deltaPhi)) - TMath::Cos(TMath::DegToRad()*(degrees+deltaPhi)));
133 symmetricHist->SetBinContent(
bin+hist1->GetNbinsX(), value/(binNormalization*events));
134 symmetricHist->SetBinContent(hist1->GetNbinsX()-
bin+1, value/(binNormalization*events));
136 symmetricHist->SetAxisRange(-2,14);
137 symmetricHist->SetXTitle(
"Angle (degrees)");
138 symmetricHist->SetYTitle(
"(N/N0) / [sr]");
141 TF1* fitgaus =
new TF1(
"fitgaus",
"gaus");
142 TF2* fitexpo =
new TF1(
"fitexpo",
"expo");
143 fitgaus->SetLineColor(2);
144 fitexpo->SetLineColor(4);
145 symmetricHist->Fit(fitgaus,
"");
146 symmetricHist->Fit(fitexpo,
"+");
147 symmetricHist->Draw();
156 c1->SaveAs(
"angularDistributionHistogramWithFits.png");
void fragmentAngularDistributionHistogram()
const XML_Char int const XML_Char * value
printf("%d Experimental points found\n", nlines)