d4sigma_plot.C
Go to the documentation of this file.
1 // *********************************************************************
2 // Plot differential cross-sections - Martti Nirkko (28th Nov 2014)
3 // Compile and run in terminal: root -l -b -q d4sigma_plot.C+g
4 // *********************************************************************
5 
6 #include <TMath.h>
7 #include <TH3D.h>
8 #include <TFile.h>
9 #include <TCanvas.h>
10 #include <TPad.h>
11 #include <TGraph.h>
12 #include <TString.h>
13 #include <TStyle.h>
14 #include <iostream>
15 
16 // Compile using: root -l -b -q d4sigma_plot.C+g
17 void d4sigma_plot() {
18 
19  // Should be the same as when input rootfile was generated
20  const int COMP = 10;
21 
22  // Initialisation of variables
23  const double pi = TMath::Pi();
24  const int nsteps1 = 10*COMP, nsteps2 = 10*COMP, nsteps3 = 10*COMP, nsteps4 = 2*COMP;
25  double varmin1, varmin2, varmin3, varmin4, varmax1, varmax2, varmax3, varmax4;
26  int i,j,k,l;
27 
28  double Enu; // neutrino energy [GeV]
29  printf("Please enter neutrino energy: ");
30  scanf("%lf", &Enu);
31  printf("Trying to find input file for E_nu = %3.1lf GeV...\n", Enu);
32 
33  std::string fname = Form("data/d4sigma_hist_%3.1lfGeV.root", Enu);
34  TFile *file = new TFile(fname.c_str());
35 
36  if (!file) {
37  printf("ERROR: File not found!");
38  return;
39  }
40 
41  TH3D *hist[nsteps4];
42  for (l=0; l<nsteps4; l++) {
43  hist[l] = (TH3D*)file->Get(Form("d4sigma_hist_%d",l));
44  }
45 
46  // Get minimum/maximum of variables
47  varmin1 = hist[0]->GetXaxis()->GetXmin(); varmax1 = hist[0]->GetXaxis()->GetXmax();
48  varmin2 = hist[0]->GetYaxis()->GetXmin(); varmax2 = hist[0]->GetYaxis()->GetXmax();
49  varmin3 = hist[0]->GetZaxis()->GetXmin(); varmax3 = hist[0]->GetZaxis()->GetXmax();
50  varmin4 = -pi; varmax4 = pi;
51 
52  // Get bin edges for each variable
53  double binedge1[nsteps1+1], binedge2[nsteps2+1], binedge3[nsteps3+1];
54  binedge1[0] = varmin1; binedge2[0] = varmin2; binedge3[0] = varmin3;
55  for (i=1; i<nsteps1; i++) binedge1[i] = hist[0]->GetXaxis()->GetBinLowEdge(i+1);
56  for (j=1; j<nsteps2; j++) binedge2[j] = hist[0]->GetYaxis()->GetBinLowEdge(j+1);
57  for (k=1; k<nsteps3; k++) binedge3[k] = hist[0]->GetZaxis()->GetBinLowEdge(k+1);
58  binedge1[nsteps1] = varmax1; binedge2[nsteps2] = varmax2; binedge3[nsteps3] = varmax3;
59 
60  // Get bin widths for each variable
61  double binwidth1[nsteps1], binwidth2[nsteps2], binwidth3[nsteps3], binwidth4[nsteps4];
62  for (i=0; i<nsteps1; i++) binwidth1[i] = hist[0]->GetXaxis()->GetBinWidth(i+1);
63  for (j=0; j<nsteps2; j++) binwidth2[j] = hist[0]->GetYaxis()->GetBinWidth(j+1);
64  for (k=0; k<nsteps3; k++) binwidth3[k] = hist[0]->GetZaxis()->GetBinWidth(k+1);
65  for (l=0; l<nsteps4; l++) binwidth4[l] = (varmax4-varmin4)/nsteps4;
66 
67  // Initialise histograms with supposedly clever bins
68  TH1D* dTk = new TH1D("dTk", "d#sigma/dT_{kaon}", nsteps1, binedge1);
69  TH1D* dTl = new TH1D("dTl", "d#sigma/dT_{lepton}", nsteps2, binedge2);
70  TH1D* dct = new TH1D("dct", "d#sigma/dcos(#theta_{#nul})", nsteps3, binedge3);
71  TH1D* dph = new TH1D("dph", "d#sigma/d#phi_{kq}", nsteps4, varmin4, varmax4);
72 
73  // Get differential cross-section over T_kaon
74  double diff1Tk, diff2Tk, diff3Tk;
75  for (i=0; i<nsteps1; i++) {
76  diff1Tk = 0.;
77  for (j=0; j<nsteps2; j++) {
78  diff2Tk = 0.;
79  for (k=0; k<nsteps3; k++) {
80  diff3Tk = 0.;
81  for (l=0; l<nsteps4; l++) {
82  diff3Tk += hist[l]->GetBinContent(i+1,j+1,k+1)*binwidth4[l];
83  }
84  diff2Tk += diff3Tk*binwidth3[k];
85  }
86  diff1Tk += diff2Tk*binwidth2[j];
87  }
88  dTk->SetBinContent(i+1, diff1Tk*binwidth1[i]);
89  }
90 
91  // Get differential cross-section over T_lepton
92  double diff1Tl, diff2Tl, diff3Tl;
93  for (j=0; j<nsteps2; j++) {
94  diff1Tl = 0.;
95  for (i=0; i<nsteps1; i++) {
96  diff2Tl = 0.;
97  for (k=0; k<nsteps3; k++) {
98  diff3Tl = 0.;
99  for (l=0; l<nsteps4; l++) {
100  diff3Tl += hist[l]->GetBinContent(i+1,j+1,k+1)*binwidth4[l];
101  }
102  diff2Tl += diff3Tl*binwidth3[k];
103  }
104  diff1Tl += diff2Tl*binwidth1[i];
105  }
106  dTl->SetBinContent(j+1, diff1Tl*binwidth2[j]);
107  }
108 
109  // Get differential cross-section over cos(theta)
110  double diff1ct, diff2ct, diff3ct;
111  for (k=0; k<nsteps3; k++) {
112  diff1ct = 0.;
113  for (i=0; i<nsteps1; i++) {
114  diff2ct = 0.;
115  for (j=0; j<nsteps2; j++) {
116  diff3ct = 0.;
117  for (l=0; l<nsteps4; l++) {
118  diff3ct += hist[l]->GetBinContent(i+1,j+1,k+1)*binwidth4[l];
119  }
120  diff2ct += diff3ct*binwidth2[j];
121  }
122  diff1ct += diff2ct*binwidth1[i];
123  }
124  dct->SetBinContent(k+1, diff1ct*binwidth3[k]);
125  }
126 
127  // Get differential cross-section over phi_kq
128  double diff1ph, diff2ph, diff3ph;
129  for (l=0; l<nsteps4; l++) {
130  diff1ph = 0.;
131  for (i=0; i<nsteps1; i++) {
132  diff2ph = 0.;
133  for (j=0; j<nsteps2; j++) {
134  diff3ph = 0.;
135  for (k=0; k<nsteps3; k++) {
136  diff3ph += hist[l]->GetBinContent(i+1,j+1,k+1)*binwidth3[k];
137  }
138  diff2ph += diff3ph*binwidth2[j];
139  }
140  diff1ph += diff2ph*binwidth1[i];
141  }
142  dph->SetBinContent(l+1, diff1ph*binwidth4[l]);
143  }
144 
145  double sum1 = dTk->Integral();
146  double sum2 = dTl->Integral();
147  double sum3 = dct->Integral();
148  double sum4 = dph->Integral();
149 
150  printf("Integrals are:\n %.6e\n %.6e\n %.6e\n %.6e\n",sum1,sum2,sum3,sum4);
151 
152  dTk->SetMinimum(0);
153  dTl->SetMinimum(0);
154  dct->SetMinimum(0);
155  dph->SetMinimum(0);
156 
157  gStyle->SetOptStat(0);
158 
159  std::string outname = Form("d4sigma_plot_%3.1lfGeV", Enu);
160 
161  dTl->SetXTitle("T_{lepton} [GeV]");
162  dTk->SetXTitle("T_{kaon} [GeV]");
163  dct->SetXTitle("cos(#theta_{#nul}) [ ]");
164  dph->SetXTitle("#phi_{Kq} [rad]");
165 
166  TCanvas *c1 = new TCanvas("c1", "Differential cross-sections", 800, 600);
167  c1->Divide(2,2);
168  c1->cd(1); dTl->Draw();
169  c1->cd(2); dTk->Draw();
170  c1->cd(3); dct->Draw();
171  c1->cd(4); dph->Draw();
172  c1->Print(("images/"+outname+".png").c_str()); c1->Close();
173 
174  TFile* outfile = new TFile(("data/"+outname+".root").c_str(), "RECREATE");
175  dTl->Write("dsigma_dTlepton");
176  dTk->Write("dsigma_dTkaon");
177  dct->Write("dsigma_dcostheta");
178  dph->Write("dsigma_dphikq");
179  outfile->Close();
180  std::cout << std::endl << "Output written to file: " << outname << ".root" << std::endl << std::endl;
181 
182 }
183 
void d4sigma_plot()
Definition: d4sigma_plot.C:17
correl_xv GetYaxis() -> SetDecimals()
correl_xv GetXaxis() -> SetDecimals()
printf("%d Experimental points found\n", nlines)
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24
ratio_hxv GetZaxis() -> SetRangeUser(0, 2)
FILE * outfile
Definition: dump_event.C:13
enum BeamMode string