d4sigma_hist.C
Go to the documentation of this file.
1 // *********************************************************************
2 // Generate 4D diff-xsec histograms - Martti Nirkko (28th Nov 2014)
3 // Compile and run in terminal: root -l -q d4sigma_hist.C+g
4 // *********************************************************************
5 
6 #include "code/singlekaon_xsec.cxx"
7 #include <TMath.h>
8 #include <TFile.h>
9 #include <TH3D.h>
10 
11 // Compile using: root -l -b -q d4sigma_hist.C+g
12 void d4sigma_hist() {
13 
14  // ***************************
15  // * STEERING PARAMETERS *
16  // ***************************
17  const int COMP = 10; // time complexity (runs with O(n^4)...)
18  const int pLeptonIsUsed = 0; // use momentum instead of kinetic energy
19  const int thetaIsUsed = 0; // use angle instead of cosine
20  const double pi = TMath::Pi();
21 
22  // NOTE: All results for a changing "pleptonIsUsed" and "thetaIsUsed" give the same results within
23  // about 1% regarding the total cross-section. Here, for making histograms of the differential
24  // 1D cross-section, I use the variables I want to actually plot (Tlep and costheta)
25 
26  int type = 2; // lepton: 1=electron, 2=muon, 3=tau
27  int reac = 3; // reaction: 1=NN, 2=NP, 3=PP
28 
29  double Enu; // neutrino energy [GeV]
30  printf("Please enter neutrino energy: ");
31  scanf("%lf", &Enu);
32  printf("E_nu set to %.3lf GeV...\n", Enu);
33  singlekaon_xsec *inst = new singlekaon_xsec();
34  inst->init(Enu, type, reac);
35 
36  double Mlep = 0.; // lepton mass [GeV]
37  if (type==1) Mlep = 0.510998928e-3;
38  else if (type==2) Mlep = 0.1056583715;
39  else if (type==3) Mlep = 1.77682;
40  else {std::cout<<"ERROR: Invalid type!"<<std::endl; return;}
41 
42  double Mkaon = 0.; // kaon mass [GeV]
43  if (reac==1) Mkaon = 0.493677;
44  else if (reac==2) Mkaon = 0.497614;
45  else if (reac==3) Mkaon = 0.493677;
46  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
47 
48  // Initialisation of variables
49  int i,j,k,l;
50  const int nsteps1 = 10*COMP, nsteps2 = 10*COMP, nsteps3 = 10*COMP, nsteps4 = 2*COMP;
51  double varmin1, varmin2, varmin3, varmin4, varmax1, varmax2, varmax3, varmax4, temp;
52  double binsvar1[nsteps1+1], binsvar2[nsteps2+1], binsvar3[nsteps3+1], binsvar4[nsteps4+1];
53  double varvals1[nsteps1], varvals2[nsteps2], varvals3[nsteps3], varvals4[nsteps4];
54 
55  double Tlep, Elep; // LEPTON ENERGY [GeV]
56  double Tkaon; // KAON ENERGY [GeV]
57  double theta; // LEPTON ANGLE [rad]
58  double phikq; // KAON ANGLE [rad]
59  double diff4sigma; // DIFFERENTIAL CROSS SECTION [cm^2/GeV^2/rad^2]
60 
61  // Integration ranges
62  varmin1 = 0.; varmax1 = Enu-Mkaon-Mlep; // Tkaon
63  if (pLeptonIsUsed) { // plep
64  varmin2 = 0.;
65  varmax2 = sqrt((Enu-Mkaon)*(Enu-Mkaon)-Mlep*Mlep);
66  } else { // Tlep
67  varmin2 = 0.;
68  varmax2 = Enu-Mkaon-Mlep;
69  }
70  if (thetaIsUsed) { // theta
71  varmin3 = 0.;
72  varmax3 = pi;
73  } else { // cos(theta)
74  varmin3 = -1.0;
75  varmax3 = 1.0;
76  }
77  varmin4 = -pi; varmax4 = pi; // phi_kq
78 
79  // Calculate edges of bins
80  for (i=0; i<=nsteps1; i++) {
81  binsvar1[i] = varmin1 + i*(varmax1-varmin1)/nsteps1;
82  }
83  for (j=0; j<=nsteps2; j++) {
84  temp = varmin2 + j*(varmax2-varmin2)/nsteps2; // plep OR Tlep
85  if (pLeptonIsUsed) {
86  binsvar2[j] = sqrt(temp*temp+Mlep*Mlep)-Mlep; // plep -> Tlep
87  } else {
88  binsvar2[j] = temp; // Tlep
89  }
90  }
91  for (k=0; k<=nsteps3; k++) {
92  if (thetaIsUsed) {
93  temp = varmax3 - k*(varmax3-varmin3)/nsteps3; // theta
94  binsvar3[k] = cos(temp); // cos(theta)
95  } else {
96  binsvar3[k] = varmin3 + k*(varmax3-varmin3)/nsteps3;
97  }
98  }
99  for (l=0; l<=nsteps4; l++) {
100  binsvar4[l] = varmin4 + l*(varmax4-varmin4)/nsteps4;
101  }
102 
103  // Calculate edges of bins (NEW / TEST)
104  /*for (i=0; i<=nsteps1; i++) binsvar1[i] = varmin1 + i*(varmax1-varmin1)/nsteps1;
105  for (j=0; j<=nsteps2; j++) binsvar2[j] = varmin2 + j*(varmax2-varmin2)/nsteps2; // plep OR Tlep
106  for (k=0; k<=nsteps3; k++) binsvar3[k] = varmin3 + k*(varmax3-varmin3)/nsteps3; // theta OR costh
107  for (l=0; l<=nsteps4; l++) binsvar4[l] = varmin4 + l*(varmax4-varmin4)/nsteps4;*/
108 
109  // Calculate central bin values (average of left and right edge of bin)
110  for (i=0; i<nsteps1; i++) varvals1[i] = 0.5*(binsvar1[i+1]+binsvar1[i]);
111  for (j=0; j<nsteps2; j++) varvals2[j] = 0.5*(binsvar2[j+1]+binsvar2[j]);
112  for (k=0; k<nsteps3; k++) varvals3[k] = 0.5*(binsvar3[k+1]+binsvar3[k]);
113  for (l=0; l<nsteps4; l++) varvals4[l] = 0.5*(binsvar4[l+1]+binsvar4[l]);
114 
115  TH3D *hist[nsteps4];
116  for (l=0; l<nsteps4; l++) {
117  hist[l] = new TH3D(Form("hist_%2d",l), "diff4sigma", nsteps1, binsvar1, nsteps2, binsvar2, nsteps3, binsvar3);
118  }
119 
120  // CALCULATE CROSS-SECTION
121  // -----------------------
122 
123  for (i=0; i<nsteps1; i++) {
124  if (binsvar1[i+1]==binsvar1[i]) continue;
125  Tkaon = varvals1[i];
126 
127  for (j=0; j<nsteps2; j++) {
128  if (binsvar2[j+1]==binsvar2[j]) continue;
129  Tlep = varvals2[j];
130  Elep = Tlep+Mlep;
131 
132  for (k=0; k<nsteps3; k++) {
133  if (binsvar3[k+1]==binsvar3[k]) continue;
134  theta = acos(varvals3[k]);
135 
136  for (l=0; l<nsteps4; l++) {
137  if (binsvar4[l+1]==binsvar4[l]) continue;
138  phikq = varvals4[l];
139 
140  // Calculate 4D differential xsec
141  diff4sigma = inst->diffxsec(Tlep,Tkaon,theta,phikq);
142  diff4sigma *= 2*pi;
143 
144  // This Jacobian no longer needed, since binning is adjusted to cos(theta)
145  // diff4sigma *= sin(theta);
146 
147  // Multiplication with Jacobian for transformation plep->Tlep (it is E/p)
148  diff4sigma *= Elep/sqrt(Elep*Elep-Mlep*Mlep);
149 
150  hist[l]->SetBinContent(i+1,j+1,k+1,diff4sigma);
151  }
152 
153  }
154 
155  }
156 
157  if ((i+1)%5==0)
158  std::cout << "Processing... (" << 100.*(i+1.)/nsteps1 << "%)" << std::endl;
159 
160  }
161 
162  std::string fname = Form("data/d4sigma_hist_%3.1lfGeV.root", Enu);
163  TFile* outfile = new TFile(fname.c_str(), "RECREATE");
164  for (l=0; l<nsteps4; l++) hist[l]->Write(Form("d4sigma_hist_%d",l));
165  outfile->Close();
166  std::cout << std::endl << "Output written to file: " << fname << std::endl << std::endl;
167 
168 }
169 
void d4sigma_hist()
Definition: d4sigma_hist.C:12
T sqrt(T number)
Definition: d0nt_math.hpp:156
T acos(T number)
Definition: d0nt_math.hpp:54
printf("%d Experimental points found\n", nlines)
double diffxsec(double Tlep, double Tkaon, double theta, double phikq)
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
T cos(T number)
Definition: d0nt_math.hpp:78
void init(double Etot, int type, int reac)
FILE * outfile
Definition: dump_event.C:13
gm Write()
enum BeamMode string