Geant4_10
fragmentYieldsPlot.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 the yield calculations for fragments
17  * below the angle of 10 degrees.
18  * Then plots the data oin comparison to Haettner
19  *
20  * There is no fit here yet.
21  *
22  * @author Gillis Danielsen
23  * **************************/
24 
25 
27  gStyle->SetOptStat(0000000000); //remove the for this graphs totally redundant statbox
28  int ZNumGiven;
29  cout << "Enter fragment Z-number (eg. 1): ";
30  cin >> ZNumGiven;
31 
32  TCanvas *c1 = new TCanvas("fragmentYieldsPlot", "Total yield of fragments zero to ten degrees as function of depth");
33 
34  TString fragmentNameChoices[6] = {"H","He","Li","Be","B","C"};
35 
36  TString fragmentName = fragmentNameChoices[ZNumGiven-1];
37 
38  std::cout << fragmentName << endl;
39 
40  TH1F* dummyHisto = new TH1F("dummyHisto", fragmentName + " yields 0-10 degrees" ,100, 0.0,40); //Dummyhisto fix for missing TNtuple methods.
41  dummyHisto->SetXTitle("Depth (cm)");
42  dummyHisto->SetYTitle("N/N0");
43 
44  ifstream in;
45  TString experimentalDataPath = "experimentalData/iaeaBenchmark/yields/TDK" + fragmentName + ".dat";
46 
47  ifstream in;
48 
49  //Pull in ascii/exfor-style data
50  in.open(experimentalDataPath);
51 
52  Float_t f1,f2;
53  Int_t nlines = 0;
54  TFile *f = new TFile("fragmentAngularDistribution.root","RECREATE");
55  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","x:y");
56 
57  Char_t DATAFLAG[4];
58  Int_t NDATA;
59  Char_t n1[15], n2[15];
60  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
61  in >> n1 >> n2; // Read column titles: 'Energy He B [...]'
62 
63  cout <<n1<<" "<<n2<<"\n";
64  while (1) {
65  in >> f1 >> f2;
66  if (!in.good()) break;
67  if (nlines < 500 ) printf("%f %f\n",f1,f2);
68  ntuple->Fill(f1,f2);
69  nlines++;
70  }
71  std::cout << "Imported " << nlines << " lines from data-file" << endl;
72 
73 
74 
75 
76  TNtuple *simData = new TNtuple("ntuple","Data from ascii file","depth:H:He:Li:Be:B:C");
77 
78 // gROOT->SetStyle("clearRetro");
79  //this will be used as base for pulling the experimental data
80  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
81  dir.ReplaceAll("fragmentAngularDistribution.C","");
82  dir.ReplaceAll("/./","/");
83  ifstream in;
84  simData->Fill(0.0,0.0,0.0,0.0,0.0,0.0,1.0);
85 for(int j = 1; j <= 40;j=j=j+1){
86  TString pDepth, fragment, Znum, normToOneAtZeroAngle;
87  pDepth = Form("%i",j);
88 /*
89  cout << "Enter phantom depth (eg. 27.9, see experimentalData directory for choices): ";
90  cin >> pDepth;
91 */
92  TString simulationDataPath = "IAEA_" + pDepth + ".root";
93 
94  //Let's pull in the simulation-data
95  //TFile *MCData = TFile::Open("IAEA_" + pDepth + ".root");
96  TFile *MCData = TFile::Open(simulationDataPath);
97  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
98 
99  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
100  TNtuple *metadata = (TNtuple*) MCData->Get("metaData");
101  Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
102  metadata->SetBranchAddress("events",&events);
103  metadata->SetBranchAddress("waterThickness",&waterThickness);
104  metadata->SetBranchAddress("detectorDistance",&detectorDistance);
105  metadata->SetBranchAddress("beamEnergy",&beamEnergy);
106  metadata->SetBranchAddress("energyError",&energyError);
107  metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
108  metadata->GetEntry(0); //there is just one row to consider.
109  //ALL UNITS ARE cm!
110  Double_t scatteringDistance = detectorDistance - phantomCenterDistance; //temporarily hard-coded, should be distance from target-center to detector
111 
112  Double_t degrees = 10.0;
113  Double_t r, rMin, rMax, graphMaximum = 0.0;
114  Double_t norming = events*.999;
115  TString rMinString;
116  TString rMaxString;
117 
118  rMinString = "0.00";
119  rMaxString = Form("%f", scatteringDistance*TMath::ATan(degrees*TMath::DegToRad()));
120 
121  Double_t H = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",1) + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
122  Double_t He = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",2) + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
123  Double_t Li = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",3) + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
124  Double_t Be = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",4) + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
125  Double_t B = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",5) + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
126  Double_t C = ((Double_t*) fragments->GetEntries("(Z == " + TString::Format("%i",6) + " && sqrt(posY^2 + posZ^2) < " + rMaxString + "&& sqrt(posY*posY + posZ*posZ) > " + rMinString + ")")) / norming;
127  simData->Fill(waterThickness,H,He,Li,Be,B,C);
128 
129  }
130  simData->Scan();
131  simData->SetMarkerStyle(2); //filled dot
132  simData->SetMarkerColor(kBlue);
133  graphMaximum = TMath::Max(graphMaximum, simData->GetMaximum(fragmentName));
134  graphMaximum = TMath::Max(graphMaximum, ntuple->GetMaximum("y"));
135  dummyHisto->SetMaximum(graphMaximum + .05*graphMaximum);
136  dummyHisto->Draw();
137 
138  simData->Draw(fragmentName + ":depth","","p,same");
139 
140  ntuple->SetMarkerStyle(22); //triangle
141  ntuple->SetMarkerColor(kRed);
142  ntuple->Draw("y:x","","p,same");
143  c1->SaveAs("fragmentYieldsFor" + fragmentName + ".png");
144 }
void fragmentYieldsPlot()
ifstream in
Definition: comparison.C:7
double ATan(double)
TCanvas * c1
Definition: plotHisto.C:7
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