fragmentEnergyDistributionDifferentAngles.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 
16 
17 gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox
18 
19 // gROOT->SetStyle("clearRetro");
20 
21  TString pDepth, fragment, Znum, normToOneAtZeroAngle;
22  cout << "Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
23  cin >> pDepth;
24  TString simulationDataPath = "IAEA_" + pDepth + ".root";
25  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
26  dir.ReplaceAll("basic.C","");
27  dir.ReplaceAll("/./","/");
28  ifstream in;
29  in.open(Form("experimentalData/iaeaBenchmark/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
30  Float_t f1,f2,f3, f4,f5,f6;
31  Int_t nlines = 0;
32  TFile *f = new TFile("fragmentEnergyWithAngularDistribution.root","RECREATE");
33  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","Energy:He:B:H:Li:Be");
34 
35  Char_t DATAFLAG[4];
36  Int_t NDATA;
37  Char_t n1[6], n2[2], n3[2], n4[2], n5[2], n6[2];
38  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
39  in >> n1 >> n2 >> n3 >> n4 >> n5 >> n6; // Read column titles: 'Energy He B [...]'
40 
41  cout <<n1<<" "<<n2<<" "<<n3<<" "<<n4<<" "<<n5<<" "<<n6<<"\n";
42  while (1) {
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);
47  nlines++;
48  }
49 
50 
51  //Let's pull in the simulation-data
52  TFile *MCData = TFile::Open("IAEA_200000.root");
53  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
54 
55  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
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); //there is just one row to consider.
65 
66  //good to keep for ref. G4 might give weird units due to change.
67  metadata->Scan();
68  std::cout << "Recieved metadata-row: " << events << " " << detectorDistance << " " << waterThickness << " " << beamEnergy << " " << energyError << " " << phantomCenterDistance;
69 
70 
71 //A lot of hardcoded histograms, ugly
72  Double_t binAmount = 50.0; //casting from int failed somehow, so in float temporarily, fixme
73  Double_t maxEnergy = 450.0;
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);
87 
88 
89  //ALL UNITS ARE cm!
90  Double_t detectorSideLength = 4; //40mm, as e.haettner H1 detector
91  Double_t scatteringDistance = detectorDistance - phantomCenterDistance; //temporarily hard-coded, should be distance from target-center to detector
92  Double_t degrees; //< actually radians
93 
94  Double_t r, rMin, rMax, deltaOmega, normFloat;
95  TString rMinString, rMaxString, normString;
96  TString same = "";
97  TString histName;
98  TCanvas *c3 = new TCanvas("histograms", "Distribution (at different angles)");
99  int i = 0; //so that the degree steps can be varied to unevenly spaced values separate counter is used
100 
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){
103  i++;
104  degrees = j * TMath::DegToRad();
105  //std::cout << "plotting for Z = " << Znum << " at " << j << " degrees\n";
106  //Distance from straight beam at the requested angle
107  r = scatteringDistance * TMath::Tan(degrees);
108  //now the "detector is rotated around all possible perpendicularlynangle values to beamline".
109  //This forms an annulus with rMin and RMax as otuer and inner radiuses
110  //Notice this will give a bit of approximation at small angles where at 0 degrees this gives a round sensor.
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);
116  //normalization of the bins.
117  deltaPhi = degrees - TMath::ATan(TMath::Tan(degrees) - detectorSideLength/(2*scatteringDistance)); // this should be around arctan(detectorsidelength/sd)
118  if(j != 0.0){
119  deltaOmega = 2*TMath::Pi()*(TMath::Cos(TMath::Max(0.0,degrees-deltaPhi)) - TMath::Cos(degrees+deltaPhi));
120  }else{
121  deltaOmega = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
122  }
123  normFloat = deltaOmega * events * binWidth;
124  normString = Form("/%f", normFloat);
125 
126  // The following is veryvery ugly relies on a bunch of hardcoded histograms because other solutions did not work
127  histName = Form("hist%i", i);
128  if(j != 0.0){
129  fragments->Project(histName,"energy", "(Z == " + Znum + " && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")" + normString);
130  }else{
131  fragments->Project(histName,"energy", "(Z == " + Znum + " && energy > 0 && posZ < " + rMaxString + "&& posY < " + rMaxString + " && posY > 0 && posZ > 0)" + normString);
132  }
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;
135 
136  }
137 
138 //the ugly hardcoded histograms being plotted
139  //0 degrees
140  hist1->SetLineColor(kBlue);
141  hist1->Draw();
142  //1 degree
143  hist2->SetLineColor(kGreen);
144  hist2->Draw("same"); //add "same" when also plotting 0 degrees
145  //2 degrees
146  hist3->SetLineColor(kRed);
147  hist3->Draw("same");
148  //3 degrees
149  //hist4->SetLineColor(kGreen + 5);
150  //hist4->Draw("same");
151  //4 degrees
152  hist5->SetLineColor(kGreen + 3); //gives a darker shade of green
153  hist5->Draw("same");
154  //5 degrees
155  //hist6->SetLineColor(kRed);
156  //hist6->Draw("same");
157  //6 degrees
158  hist7->SetLineColor(kRed);
159  hist7->Draw("same");
160  //7 degrees
161  //hist8->SetLineColor(kRed);
162  //hist8->Draw("same");
163  //8 degrees
164  //hist9->SetLineColor(kRed);
165  //hist9->Draw("same");
166 
167  // Legends for the data
168  leg = new TLegend(0.9,0.7,1,1); //coordinates are fractions
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");
175  leg->Draw();
176 
177  c3->SaveAs("AEDistrib" + Znum + ".png");
178 }
179 
180  in.close();
181 
182  f->Write();
183 }
Float_t f4
enum BeamMode kRed
TH1D * hist2
Definition: plotHisto.C:12
TTree * ntuple
Definition: macro.C:6
TH1D * hist3
Definition: plotHisto.C:15
constexpr T pow(T x)
Definition: pow.h:75
Float_t f2
TH1D * hist6
Definition: plotHisto.C:27
nlines
TH1D * hist4
Definition: plotHisto.C:18
Float_t f3
TH1D * hist8
Definition: plotHisto.C:25
TH1D * hist5
Definition: plotHisto.C:21
TH1D * hist1
Definition: plotHisto.C:9
void fragmentEnergyDistributionDifferentAngles()
printf("%d Experimental points found\n", nlines)
Float_t f1
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
ifstream in
Definition: comparison.C:7
TDirectory * dir
Definition: macro.C:5
TRandom3 r(0)
TH1D * hist9
Definition: plotHisto.C:33
TH1D * hist7
Definition: plotHisto.C:31
enum BeamMode kGreen
enum BeamMode kBlue
void events(int which)
Definition: Cana.C:52