LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
Functions
fragmentAngularDistribution.C File Reference
#include "Riostream.h"
#include "TSystem.h"
#include "TInterpreter.h"
#include "TROOT.h"
#include "TApplication.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "THStack.h"
#include "TCut.h"
#include "TString.h"
#include "TMath.h"

Go to the source code of this file.

Functions

void fragmentAngularDistribution ()
 

Function Documentation

void fragmentAngularDistribution ( )

Definition at line 31 of file fragmentAngularDistribution.C.

References c1, dir, f, larg4::f1, larg4::f2, in, nlines, ntuple, and printf().

31  {
32 
33  gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox
34  gROOT->SetStyle("clearRetro");
35  //this will be used as base for pulling the experimental data
36  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
37  dir.ReplaceAll("fragmentAngularDistribution.C","");
38  dir.ReplaceAll("/./","/");
39  ifstream in;
40 
41  int ZnumInt;
42  TString pDepth, fragment, Znum, normToOneAtZeroAngle;
43  cout << "Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
44  cin >> pDepth;
45  cout << "Enter fragment Z-number (eg. 1): ";
46  cin >> ZnumInt;
47  //cout << "Enter fragment name (Znum 1 -> H,Znum 2->He...): ";
48  //cin >> fragment;
49  TString fragmentNameChoices[6] = {"H","He","Li","Be","B","C"};
50  TString fragment = fragmentNameChoices[ZnumInt - 1];
51  Znum = Form("%i",ZnumInt);
52  cout << "Normalize to 1 at zero angle? (Y/N): ";
53  cin >> normToOneAtZeroAngle;
54 
55  TString experimentalDataPath = "experimentalData/iaeaBenchmark/angularDistributions/" + pDepth + "/" + fragment + "" + pDepth +".dat";
56  TString simulationDataPath = "IAEA_" + pDepth + ".root";
57 
58  TCanvas *c1 = new TCanvas("AngularDistribution", "Angular distribution with discrete measurement annuluses");
59 
60  //Pull in ascii/exfor-style data
61  in.open(experimentalDataPath);
62 
63  Float_t f1,f2;
64  Int_t nlines = 0;
65  TFile *f = new TFile("fragmentAngularDistribution.root","RECREATE");
66  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","x:y");
67 
68  Char_t DATAFLAG[4];
69  Int_t NDATA;
70  Char_t n1[15], n2[15];
71  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
72  in >> n1 >> n2; // Read column titles: 'Energy He B [...]'
73 
74  cout <<n1<<" "<<n2<<"\n";
75  while (1) {
76  in >> f1 >> f2;
77  if (!in.good()) break;
78  if (nlines < 500 ) printf("%f %f\n",f1,f2);
79  ntuple->Fill(f1,f2);
80  nlines++;
81  }
82  std::cout << "Imported " << nlines << " lines from data-file" << endl;
83 
84  //Let's pull in the simulation-data
85  //TFile *MCData = TFile::Open("IAEA_" + pDepth + ".root");
86  TFile *MCData = TFile::Open(simulationDataPath);
87  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
88 
89  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
90  TNtuple *metadata = (TNtuple*) MCData->Get("metaData");
91  Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
92  metadata->SetBranchAddress("events",&events);
93  metadata->SetBranchAddress("waterThickness",&waterThickness);
94  metadata->SetBranchAddress("detectorDistance",&detectorDistance);
95  metadata->SetBranchAddress("beamEnergy",&beamEnergy);
96  metadata->SetBranchAddress("energyError",&energyError);
97  metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
98  metadata->GetEntry(0); //there is just one row to consider.
99 
100  //good to keep for ref. G4 might give weird units due to change.
101  metadata->Scan();
102 
103  //ALL UNITS ARE cm!
104  Double_t detectorSideLength = 4; //40mm, as e.haettner H1 detector
105  Double_t scatteringDistance = detectorDistance - phantomCenterDistance; //temporarily hard-coded, should be distance from target-center to detector
106 
107  Double_t r;
108  Double_t degrees;
109  Double_t rMin;
110  Double_t rMax;
111  TString rMinString, rMaxString, experimentalNorm;
112  Double_t maxValue = 0.0; //When normalizing to 1 this will allways end up being 1)
113 
114  int i = 0; //so that the degree steps can be varied to unevenly spaced values separate counter is used
115  TNtuple* distrib = new TNtuple("angularDistrib","FragmentAngularDistrib","angle:particleAmount:normalized");
116  std::cout << "Fragments comparison to the graphs in appendices of E.Haettner\n";
117  std::cout << "Scattering distance: " << scatteringDistance << " cm" << endl ;
118  std::cout << "(scattering distance may vary with data-files too, see haettner A.1." << endl << endl;
119  //This will norm it to the zero degree entry to get rid of Emma's weird normalization
120  //fixme detectorsidelengthstring
121  rMinString = "0.00";
122  rMaxString = Form("%f", detectorSideLength/2);
123  //SA of a square Double_t zeroSAsquare = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
124  //SA with square approx squareApprox = (4*4) / (4*3.14*scatteringDistance*scatteringDistance);
125 
126  //First calculates the normalization from the zero position
127  //Normalization by events becomes redundant but is left in place for future needs
128 
129  Double_t deltaPhi = TMath::ATan((detectorSideLength/2)/scatteringDistance); //Angle where side of detector is found
130  /*
131  //Alternative normalization, here zero position is also done with annulus where rMin=0, the actual detector is a square though
132  //Difference with this approach and the other is very small
133  Double_t normEntries = fragments->GetEntries("(Z == " + Znum + " && energy > 0 && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")");
134  Double_t zeroSA = 2*TMath::Pi()*(TMath::Cos(0) - TMath::Cos(deltaPhi));
135  */
136  //fragments->Scan();
137  //Results are normalized by a square detector mimicing H1 with center at 0 degrees.
138  Double_t normEntries = fragments->GetEntries("(Z == " + Znum + " && posY < " + rMaxString + " && posY > -" + rMaxString + " && posZ > -" + rMaxString + " && posZ < " + rMaxString + ")");
139  Double_t zeroSA = 4 * TMath::ASin(pow(detectorSideLength,2.0) / (4*pow(scatteringDistance,2) + pow(detectorSideLength,2)) );
140  Double_t zeroYieldNormed = normEntries / (events * zeroSA);
141  if(normToOneAtZeroAngle == "Y"){
142  Double_t zeroNorm = zeroYieldNormed; //values normalized to one at zero
143  }else{
144  Double_t zeroNorm = 1.0; //non-zeronormalized values
145  }
146 
147  distrib->Fill(0,normEntries,zeroYieldNormed/zeroNorm); //< degrees, entyamount, normalized result for graph
148  //fragments->Scan(); //debug
149  std::cout << "Norming events: " << normEntries << endl;
150  //Loop through all other wanted angles, too large angles will fall outside reach of phantom window.
151  for(Double_t j = deltaPhi*TMath::RadToDeg(); j <= 15.0; j=j+.05){
152  i++;
153  degrees = j * TMath::DegToRad();
154  //Distance from straight beam at the requested angle
155  r = scatteringDistance * TMath::Tan(degrees);
156  //now the "detector is rotated around all possible perpendicularlynangle values to beamline".
157  //This forms an annulus with rMin and Rmax as outer and inner radiuses
158  //this will give a bit of approximation at small angles and at 0 degrees this gives a completely round sensor.
159 
160 
161  /*
162  * deltaPhi calculated so that phi+deltaphi points to one side of the detector and phi-deltaphi the other side
163  *
164  * Alternative 1: detector is moved but the normal is not pointed towards the scattering source.
165  * Alternative 2: detector is moved and pointed towards scattering source. (E.Haettner seems to use this)
166  *
167  * The difference in results is very minute though. (3% at largest angles)
168  */
169  //Alternative 1
170  //Double_t deltaPhi = degrees - TMath::ATan(TMath::Tan(degrees) - (detectorSideLength/(2*scatteringDistance)));
171  //rMin = TMath::Max(0.0,r - (detectorSideLength/2));
172  //rMax = r + (detectorSideLength/2);
173  //Alternative 2
174  Double_t deltaPhi = TMath::ATan((TMath::Cos(degrees)*detectorSideLength)/(2*scatteringDistance));
175  rMin = TMath::Max(0.0,r - (detectorSideLength/(2*TMath::Cos(degrees))));
176  rMax = rMin + ((detectorSideLength*TMath::Sin(degrees))/TMath::Tan((TMath::Pi()/2) - degrees - deltaPhi)) + (detectorSideLength*TMath::Cos(degrees));
177  rMinString = Form("%f", rMin);
178  rMaxString = Form("%f", rMax);
179  /*
180  * From Gunzert-marx. Solid angle of annulus with rmin trmax,
181  * a bit of an aproximation especially at small phi.
182  */
183  Double_t deltaOmega = 2*TMath::Pi()*(TMath::Cos(TMath::Max(0.0,degrees-deltaPhi)) - TMath::Cos(degrees+deltaPhi));
184  int numEntries = fragments->GetEntries("(Z == " + Znum + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")");
185  distrib->Fill(j,numEntries,numEntries/(deltaOmega * events * zeroNorm)); //< degrees, entyamount, normalized result for graph
186  distrib->Fill(-j,numEntries,numEntries/(deltaOmega * events * zeroNorm)); //< To get gaussian shape better visible
187  maxValue = TMath::Max(maxValue, numEntries/(deltaOmega * events * zeroNorm)); //< for calculation of FWHM
188  }
189  distrib->SetMarkerStyle(2); //filled dot
190  distrib->SetMarkerColor(kBlue);
191  ntuple->SetMarkerStyle(22); //triangle
192  ntuple->SetMarkerColor(kRed);
193 
194  TH1F* dummHisto = new TH1F("dummyHisto", fragment + ", " + Form("%.1f", waterThickness) + " cm",100, -3.0,14); //Dummyhisto fix for missing TNtuple methods.
195  dummyHisto->SetXTitle("Angle (degrees)");
196  dummyHisto->SetYTitle("(N/N0) [sr^-1]");
197  if(normToOneAtZeroAngle == "Y"){
198  dummyHisto->SetMaximum(1.1);
199  dummyHisto->SetYTitle("[sr^-1]");
200  Float_t zeroPosData; //This is where we store what we norm the experimental data with
201  Float_t zeroPosAngle; //okay, so this should be zero, but regrettably is not allways that
202  ntuple->SetBranchAddress("y",&zeroPosData);
203  ntuple->SetBranchAddress("x",&zeroPosAngle);
204  int row = 0;
205  ntuple->GetEntry(row); //Pull the first row, usually is the right one
206  while(zeroPosAngle*zeroPosAngle > .01){
207  row++;
208  ntuple->GetEntry(row);
209  if(row == ntuple->GetEntries()){
210  std::cerr << "Could not find zero angle data in imported experimental data. Change normalization or relax exactness of this check." << endl;
211  exit();
212  }
213  }
214 
215  std::cout << "For zero-position of experimental data using angle " << zeroPosAngle << " with amount " << zeroPosData << " on row " << row << endl;
216  experimentalNorm = Form("(1/%f)*", zeroPosData);
217  }else{
218  //nor normalization to 1 of data
219  dummyHisto->SetMaximum(ntuple->GetMaximum("y")+ ntuple->GetMaximum("y")*.1);
220  experimentalNorm = ""; //no norming of experimental resilts
221  }
222  dummyHisto->Draw();
223  ntuple->Draw(experimentalNorm + "y:x","","p,same");
224  distrib->Draw("normalized:angle","angle > -3 && angle < 14","p,same"); //similar axises to e.haettner
225 
226  //Calculate closest-point-FWHM.
227  Float_t fwhm = 0.0, middle = 0.0, currentX, currentY;
228  distrib->SetBranchAddress("normalized",&currentY);
229  distrib->SetBranchAddress("angle",&currentX);
230  for(int i = 0; i < distrib->GetEntries(); i++){
231  distrib->GetEntry(i);
232  if(pow(maxValue/2 - middle, 2.0) > pow(maxValue/2 - currentY, 2.0)){
233  fwhm = 2*currentX;
234  middle = currentY;
235  }else{
236  }
237  }
238  std::cout << "Calculated (closest point) FWHM of Monte-Carlo simulation to be: " << fwhm << " degrees" << endl;
239 
240  /*
241  * This code is left here because it allows to calcualte the values without using annuluses
242  * this of course is a bit more like the experimental data but statistically less precise.
243  for(Double_t p = 0.0;p < 14.0; p = p + 1.0){
244  rMinString = "0.00";
245  rMaxString = "2.00";
246  TString plusY = Form("(posY - %f)", p);
247  TString plusZ = Form("(posZ - %f)", p);
248  Double_t deltaPhi = TMath::ATan((detectorSideLength/2)/scatteringDistance);
249  Double_t normEntries = fragments->GetEntries("(Z == " + Znum + " && energy > 0 && sqrt(" + plusY + "^2 + " + plusZ + "^2) < " + rMaxString + "&& sqrt("+plusY + "^2 + " + plusZ + "^2) > " + rMinString + ")");
250  //std::cout << "(Z == " + Znum + " && energy > 0 && sqrt(" + plusY + "^2 + " + plusZ + "^2) < " + rMaxString + "&& sqrt("+plusY + "^2 + " + plusZ + "^2) > " + rMinString + ")";
251  Double_t zeroSA = 2*TMath::Pi()*(TMath::Cos(0) - TMath::Cos(deltaPhi));
252  std::cout << "with " << p << "cm the amount is " << normEntries << " / " << normEntries /(events*zeroSA) << endl;
253  }
254  */
255 
256 
257  //c1->SaveAs("angularDistrib_depth_" + pDepth + "_Z_" + Znum + "_normedToZero_" + normToOneAtZeroAngle + "_ComparedToEHaettner.png");
258  pDepth.ReplaceAll(".","");
259  c1->SaveAs("AD_" + pDepth + "_" + Znum + "_" + normToOneAtZeroAngle + ".png");
260  in.close();
261  f->Write();
262 }
ntuple
TTree * ntuple
Definition: macro.C:6
larg4::f2
SynchrotronAndGN f2
Definition: CustomPhysicsBuiltIns.cxx:11
f
TFile f
Definition: plotHisto.C:6
nlines
nlines
Definition: comparison_ascii.C:57
c1
TCanvas * c1
Definition: plotHisto.C:7
larg4::f1
EmPhysicsFactory f1
Definition: CustomPhysicsBuiltIns.cxx:10
printf
printf("%d Experimental points found\n", nlines)
in
ifstream in
Definition: comparison.C:7
dir
TDirectory * dir
Definition: macro.C:5