Geant4_10
fragmentAngularDistributionGM.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 a graph of the angular distribution of
17  * a certain type of charged fragment (default Z=1).
18  *
19  * Results are not stored in a histogram.
20  *
21  * This can be compared to measurements made with a square
22  * detector that is being moved around. Such as that of E.Haettner[1].
23  *
24  * Results are normalized to the 0-angle because documentation on E.Haettner's
25  * normalization is not found.
26  *
27  * @author Gillis Danielsen
28  * **************************/
29 
30 
32 
33  TCanvas *c1 = new TCanvas("AngularDistribution", "Angular distribution with discrete measurement annuluses");
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  ifstream in;
39 
40  //Settings for analysis
41  TString pDepth = "5.9"; //set ehre what depth to analyse, requries suitable root file and data file
42  TString fragment = "H"; //string of what fragment is being looked at H, He, Li ...
43  TString Anum = "1"; //set here what will be put in the selection for z, does not automatically change imoprted data.
44  //TString experimentalDataFile = "experimentalData/iaeaBenchmark/angularDistributions/" + pDepth + "/" + fragment + pDepth + ".dat";
45  //in.open(Form("experimentalData/iaeaBenchmark/angularDistributions/" + pDepth + "/" + fragment + pDepth + ".dat",dir.Data()));
46  in.open(Form("experimentalData/iaeaBenchmark/alternativeGMBeamAD.dat",dir.Data())); //fixme: redundant dir.data()
47  std::cout << "experimentalData/iaeaBenchmark/angularDistributions/" + pDepth + "/" + fragment + pDepth + ".dat" << endl;
48  std::cout << "didata" << dir.Data() << endl;
49  Float_t f1,f2;
50  Int_t nlines = 0;
51  TFile *f = new TFile("fragmentAngularDistribution.root","RECREATE");
52  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","x:y");
53 
54  Char_t DATAFLAG[4];
55  Int_t NDATA;
56  Char_t n1[15], n2[15];
57  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
58  in >> n1 >> n2; // Read column titles: 'Energy He B [...]'
59 
60  cout <<n1<<" "<<n2<<"\n";
61  while (1) {
62  in >> f1 >> f2;
63  if (!in.good()) break;
64  if (nlines < 500 ) printf("%f %f\n",f1,f2);
65  ntuple->Fill(f1,f2);
66  nlines++;
67  }
68  std::cout << "Imported " << nlines << " lines from data-file" << endl;
69 
70  //Let's pull in the simulation-data
71  //TFile *MCData = TFile::Open("IAEA_" + pDepth + ".root");
72  TFile *MCData = TFile::Open("IAEA_G-M.root");
73  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
74 
75  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
76  TNtuple *metadata = (TNtuple*) MCData->Get("metaData");
77  Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
78  metadata->SetBranchAddress("events",&events);
79  metadata->SetBranchAddress("waterThickness",&waterThickness);
80  metadata->SetBranchAddress("detectorDistance",&detectorDistance);
81  metadata->SetBranchAddress("beamEnergy",&beamEnergy);
82  metadata->SetBranchAddress("energyError",&energyError);
83  metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
84  metadata->GetEntry(0); //there is just one row to consider.
85 
86  //good to keep for ref. G4 might give weird units due to change.
87  metadata->Scan();
88 
89  //ALL UNITS ARE cm!
90  Double_t detectorSideLength = 4; //40mm, as e.haettner H1 detector, G-M uses 1.6cm sidelength (which is however to small for < 1 million events files
91  Double_t scatteringDistance = detectorDistance - phantomCenterDistance; //temporarily hard-coded, should be distance from target-center to detector
92 
93  Double_t r;
94  Double_t degrees;
95  Double_t rMin;
96  Double_t rMax;
97  TString rMinString;
98  TString rMaxString;
99  Double_t maxValue = 0.0; //When normalizing to 1 this will allways end up being 1)
100 
101  int i = 0; //so that the degree steps can be varied to unevenly spaced values separate counter is used
102  TNtuple* distrib = new TNtuple("angularDistrib","FragmentAngularDistrib","angle:particleAmount:normalized");
103  std::cout << "Fragments comparison to the graphs in appendices of E.Haettner\n";
104  std::cout << "Scattering distance: " << scatteringDistance << " cm" << endl ;
105  std::cout << "(scattering distance may vary with data-files too, see haettner A.1." << endl << endl;
106  //This will norm it to the zero degree entry to get rid of Emma's weird normalization
107  //fixme detectorsidelengthstring
108  rMinString = "0.00";
109  rMaxString = Form("%f", detectorSideLength/2);
110  //SA of a square Double_t zeroSAsquare = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
111  //SA with square approx squareApprox = (4*4) / (4*3.14*scatteringDistance*scatteringDistance);
112 
113  //First calculates the normalization from the zero position
114  //Normalization by events becomes redundant but is left in place for future needs
115 
116  Double_t deltaPhi = TMath::ATan((detectorSideLength/2)/scatteringDistance); //Angle where side of detector is found
117  /*
118  //Alternative normalization, here zero position is also done with annulus where rMin=0, the actual detector is a square though
119  //Difference with this approach and the other is very small
120  Double_t normEntries = fragments->GetEntries("(Z == " + Anum + " && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")");
121  Double_t zeroSA = 2*TMath::Pi()*(TMath::Cos(0) - TMath::Cos(deltaPhi));
122  */
123  //fragments->Scan();
124  //Results are normalized by a square detector mimicing H1 with center at 0 degrees.
125  Double_t normEntries = fragments->GetEntries("(A == " + Anum + " && posY < " + rMaxString + " && posY > -" + rMaxString + " && posZ > -" + rMaxString + " && posZ < " + rMaxString + ")");
126  Double_t zeroSA = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
127  Double_t zeroNorm = normEntries / (events * zeroSA);
128  distrib->Fill(0,normEntries,zeroNorm); //< degrees, entyamount, normalized result for graph
129  //fragments->Scan();
130  std::cout << "Norming events: " << normEntries << endl;
131  //Loop through all other wanted angles, too large angles will fall outside reach of phantom window.
132  for(Double_t j = deltaPhi*TMath::RadToDeg(); j <= 15.0; j=j+.05){
133  i++;
134  degrees = j * TMath::DegToRad();
135  //Distance from straight beam at the requested angle
136  r = scatteringDistance * TMath::Tan(degrees);
137  //now the "detector is rotated around all possible perpendicularlynangle values to beamline".
138  //This forms an annulus with rMin and Rmax as outer and inner radiuses
139  //this will give a bit of approximation at small angles and at 0 degrees this gives a completely round sensor.
140  Double_t deltaPhi = TMath::ATan((TMath::Cos(degrees)*detectorSideLength)/(2*scatteringDistance));
141  rMin = TMath::Max(0.0,r - (detectorSideLength/(2*TMath::Cos(degrees))));
142  rMax = rMin + ((detectorSideLength*TMath::Sin(degrees))/TMath::Tan((TMath::Pi()/2) - degrees - deltaPhi)) + (detectorSideLength*TMath::Cos(degrees));
143  rMinString = Form("%f", rMin);
144  rMaxString = Form("%f", rMax); /*
145  * From Gunzert-marx. Solid angle of annulus with rmin trmax,
146  * a bit of an aproximation especially at small phi.
147  */
148  Double_t deltaOmega = 2*TMath::Pi()*(TMath::Cos(TMath::Max(0.0,degrees-deltaPhi)) - TMath::Cos(degrees+deltaPhi));
149  int numEntries = fragments->GetEntries("(A == " + Anum + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")");
150  //distrib->Fill(j,numEntries,numEntries/(deltaOmega * events * zeroNorm)); //< degrees, entyamount, normalized result for graph
151  distrib->Fill(j,numEntries,numEntries/(deltaOmega * events));
152  //distrib->Fill(-j,numEntries,numEntries/(deltaOmega * events * zeroNorm)); //< To get gaussian shape better visible
153  distrib->Fill(-j,numEntries,numEntries/(deltaOmega * events));
154  maxValue = TMath::Max(maxValue, numEntries/(deltaOmega * events)); //< for calculation of FWHM
155  }
156  distrib->SetMarkerStyle(2); //filled dot
157  distrib->SetMarkerColor(kBlue);
158  distrib->Draw("normalized:angle","angle > -3 && angle < 14","p"); //similar axises to e.haettner
159  ntuple->SetMarkerStyle(22); //triangle
160  ntuple->SetMarkerColor(kRed);
161  Float_t zeroPosData; //This is where we store what we norm the experimental data with
162  Float_t zeroPosAngle; //okay, so this should be zero, but regrettably is not allways that
163  ntuple->SetBranchAddress("y",&zeroPosData);
164  ntuple->SetBranchAddress("x",&zeroPosAngle);
165  int row = 0;
166  ntuple->GetEntry(row); //Pull the first row, usually is the right one
167  while(zeroPosAngle*zeroPosAngle > .01){
168  row++;
169  ntuple->GetEntry(row);
170  if(row == ntuple->GetEntries()){
171  std::cerr << "Could not find zero angle data in imported experimental data. Change normalization or relax exactness of this check." << endl;
172  exit();
173  }
174  }
175  std::cout << "For zero-position of experimental data using angle " << zeroPosAngle << " with amount " << zeroPosData << " on row " << row << endl;
176  TString experimentalNorm = Form("(1/%f)*", zeroPosData);
177  ntuple->Draw(experimentalNorm + "y:x","","p,same");
178 
179  //Calculate FWHM, ineffective solution. (but works also without normalizing to 1)
180  Float_t fwhm = 0.0, middle = 0.0, currentX, currentY;
181  distrib->SetBranchAddress("normalized",&currentY);
182  distrib->SetBranchAddress("angle",&currentX);
183  for(int i = 0; i < distrib->GetEntries(); i++){
184  distrib->GetEntry(i);
185  if(pow(maxValue/2 - middle, 2.0) > pow(maxValue/2 - currentY, 2.0)){
186  fwhm = 2*currentX;
187  middle = currentY;
188  }else{
189  }
190  }
191  std::cout << "Calculated (closest point) FWHM of Monte-Carlo simulation to be: " << fwhm << " degrees" << endl;
192 
193  /*
194  * This code is left here because it allows to calcualte the values without using annuluses
195  * this of course is a bit more like the experimental data but statistically less precise.
196  for(Double_t p = 0.0;p < 14.0; p = p + 1.0){
197  rMinString = "0.00";
198  rMaxString = "2.00";
199  TString plusY = Form("(posY - %f)", p);
200  TString plusZ = Form("(posZ - %f)", p);
201  Double_t deltaPhi = TMath::ATan((detectorSideLength/2)/scatteringDistance);
202  Double_t normEntries = fragments->GetEntries("(Z == " + Anum + " && energy > 0 && sqrt(" + plusY + "^2 + " + plusZ + "^2) < " + rMaxString + "&& sqrt("+plusY + "^2 + " + plusZ + "^2) > " + rMinString + ")");
203  //std::cout << "(Z == " + Anum + " && energy > 0 && sqrt(" + plusY + "^2 + " + plusZ + "^2) < " + rMaxString + "&& sqrt("+plusY + "^2 + " + plusZ + "^2) > " + rMinString + ")";
204  Double_t zeroSA = 2*TMath::Pi()*(TMath::Cos(0) - TMath::Cos(deltaPhi));
205  std::cout << "with " << p << "cm the amount is " << normEntries << " / " << normEntries /(events*zeroSA) << endl;
206  }
207  */
208 
209 
210  c1->SaveAs("AD_for_Z_" + Anum + "_ComparedToGM.png");
211  in.close();
212  f->Write();
213 }
void fragmentAngularDistributionGM()
ifstream in
Definition: comparison.C:7
double ATan(double)
TCanvas * c1
Definition: plotHisto.C:7
double ASin(double)
Definition: UUtils.hh:211
TTree * ntuple
Definition: macro.C:6
TFile f
Definition: plotHisto.C:6
Float_t f1
Float_t f2
G4int Int_t
jump r
Definition: plot.C:36
G4float Float_t
printf("%d Experimental points found\n", nlines)
G4double Double_t
nlines
TDirectory * dir
Definition: macro.C:5