checkBeam.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 "TCut.h"
11 #include "TString.h"
12 #include "TMath.h"
13 
14 /**
15  * Macro for checking beam FWHM
16  *
17  * Usage:
18  * root -l RootScripts/iaeaBenchmark/fragmentEnergy.C++
19  */
20 void checkBeam() {
21 ////////////////////////////////////////
22 ////// Importing data /////////
23 ////////////////////////////////////////
24  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
25  dir.ReplaceAll("fragmentEnergy.C","");
26  dir.ReplaceAll("/./","/");
27 
28 TString macroPath(gROOT->GetMacroPath());
29 gROOT->SetMacroPath(macroPath + ":RootScripts/iaeaBenchmark");
30 //gROOT->LoadMacro("rootlogon.C");
31 //gROOT->SetStyle("clearRetro"); //For stylesheet
32 
33  ifstream in;
34  in.open(Form("experimentalData/iaeaBenchmark/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
35  Float_t f1,f2,f3, f4,f5,f6;
36  Int_t nlines = 0;
37  TFile *f = new TFile("fragmentEnergy.root","RECREATE");
38  TNtuple *ntuple = new TNtuple("ntuple","Data from ascii file","Energy:He:B:H:Li:Be");
39 
40  Char_t DATAFLAG[4];
41  Int_t NDATA;
42  Char_t n1[6], n2[2], n3[2], n4[2], n5[2], n6[2];
43  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
44  in >> n1 >> n2 >> n3 >> n4 >> n5 >> n6; // Read column titles: 'Energy He B [...]'
45 
46  cout <<n1<<" "<<n2<<" "<<n3<<" "<<n4<<" "<<n5<<" "<<n6<<"\n";
47  while (1) {
48  in >> f1 >> f2 >> f3 >>f4 >> f5 >> f6;
49  if (!in.good()) break;
50  if (nlines < 500 ) printf("%f %0.2f %0.2f %0.2f %0.2f %0.2f \n",f1,f2,f3,f4,f5,f6);
51  ntuple->Fill(f1,f2,f3,f4,f5,f6);
52  nlines++;
53  }
54 
55  printf(" found %d points\n",nlines);
56 
57  //Let's pull in the monte carlo simulation results
58  TFile *MCData = TFile::Open("IAEA_15.9.root");
59  TH1F* MC_helium = (TH1F*)MCData->Get("heliumEnergyAfterPhantom");
60  TH1F* MC_hydrogen = (TH1F*)MCData->Get("hydrogenEnergyAfterPhantom");
61 //scale and plot
62  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
63 
64  //Block bellow pulls out the simulation's metadata from the metadata ntuple.
65  TNtuple *metadata = (TNtuple*) MCData->Get("metaData");
66  Float_t events, detectorDistance,waterThickness,beamEnergy,energyError,phantomCenterDistance;
67  metadata->SetBranchAddress("events",&events);
68  metadata->SetBranchAddress("waterThickness",&waterThickness);
69  metadata->SetBranchAddress("detectorDistance",&detectorDistance);
70  metadata->SetBranchAddress("beamEnergy",&beamEnergy);
71  metadata->SetBranchAddress("energyError",&energyError);
72  metadata->SetBranchAddress("phantomCenterDistance",&phantomCenterDistance);
73  metadata->GetEntry(0); //there is just one row to consider.
74 
75 //analysis numbers based on metadata
76  Double_t scatteringDistance = detectorDistance - phantomCenterDistance;
77  Double_t detectorSideLength = 4; //hardcoded, we have a zero angle square detector
78 //good to keep for ref. G4 might give weird units due to change.
79 metadata->Scan();
80 
81  //TCanvas* c3 = new TCanvas();
82  fragments->SetLineColor(kRed);
83  fragments->SetMarkerStyle(22);
84 
85  //fragments->Draw("posY:posZ", "abs(posZ) < 2 && abs(posY) < 2");
86  TH1F* posYHisto = new TH1F("vertical","Vertical distribution of beam",200,-2.,2.);
87  TH1F* posZHisto = new TH1F("horizontal","horizontal distribution of beam",200,-2.,2.);
88  TH2F* posHisto = new TH2F("posHisto","Distribution of beam",200,-2.,2.,200,-2.,2.);
89  //fragments->Project("posHisto","posY:posZ","abs(posZ) < 2 && abs(posY) < 2");
90  fragments->Project("vertical","posZ","abs(posZ) < 2 && abs(posY) < 2 && Z == 6");
91  fragments->Project("horizontal","posY","abs(posZ) < 2 && abs(posY) < 2 && Z == 6");
92  Float_t maxVal = posYHisto->GetMaximum();
93  std::cout << "maximum: " << maxVal << endl;
94  //fragments->Scan("posY:posZ","abs(posZ) < 2 && abs(posY) < 2");
95  //posYHisto->Draw();
96 // posZHisto->Draw("same");
97 
98  Float_t fwhm = 0.0, middle = 0.0, curVal;
99  int fwhmBin;
100  //So these dots have been binned and a fwhm is calculated.
101  for(int i = 0; i < posYHisto->GetNbinsX(); i++){
102 
103  curVal = posYHisto->GetBinContent(i);
104  if(pow(maxVal/2 - middle, 2.0) > pow(maxVal/2 - curVal, 2.0)){
105  fwhm = 2*TMath::Abs(posYHisto->GetBinCenter(i));
106  middle = curVal;
107  fwhmBin = i;
108  }else{
109  }
110  }
111  //posYHisto->SetBinContent(fwhmBin, 0);
112  /*
113  TNtuple* where = new TNtuple("where","where","x:y");
114  where->Fill(fwhmBin*posYHisto->GetBinWidth(0), 0);
115  where->Fill((200-fwhmBin)*posYHisto->GetBinWidth(0), 0);
116  where->->SetMarkerStyle(22); //triangle
117  where->SetMarkerColor(kRed);
118  */
119  TF1* fitgaus = new TF1("fitgaus","gaus");
120  fitgaus->SetLineColor(2);
121  posYHisto->Fit(fitgaus,"");
122  posYHisto->Draw();
123  std::cout << "fitted FWHM from normal distribution is: " << fitgaus->GetParameter(2)*10*2.35482 << endl;
124  //where->Draw("x:y","same");
125  //posYHisto->Smooth(150);
126  //posYHisto->Draw();
127  //posHisto->Draw();
128  std::cout << "Calculated (closest point) FWHM of Monte-Carlo simulation to be: " << fwhm*10 << " mm" << endl;
129  std::cout << "beam contained " << fragments->GetEntries("Z == 6") << " carbon nuclei." << endl;
130 
131 
132  //c3->SaveAs("checkBeam.png");
133 
134 }
135 
Float_t f4
enum BeamMode kRed
TTree * ntuple
Definition: macro.C:6
constexpr T pow(T x)
Definition: pow.h:75
void checkBeam()
Definition: checkBeam.C:20
Float_t f2
nlines
Float_t f3
printf("%d Experimental points found\n", nlines)
Float_t f1
OStream cout
Definition: OStream.cxx:6
ifstream in
Definition: comparison.C:7
TDirectory * dir
Definition: macro.C:5
void events(int which)
Definition: Cana.C:52