d3sigma_calc.C
Go to the documentation of this file.
1 // *********************************************************************
2 // Xsec (3D) of single kaon production - Martti Nirkko (28th Nov 2014)
3 // Compile and run in terminal: root -l -q d3sigma_calc.C+g
4 // *********************************************************************
5 
6 #include "code/singlekaon_xsec.cxx"
7 #include <TMath.h>
8 #include <TString.h>
9 
10 // Compile using: root -l -b -q d3sigma_calc.C+g
11 void d3sigma_calc() {
12 
13  // ***************************
14  // * STEERING PARAMETERS *
15  // ***************************
16  const int verbose = 0; // debugging mode on/off
17  const int COMP = 10; // time complexity (runs with O(n^3))
18  const int pLeptonIsUsed = 1; // use momentum instead of kinetic energy
19  const int thetaIsUsed = 1; // use angle instead of cosine
20  const int varPhiKqIsUsed = 1; // use phi_kq as function of E_nu
21  const double pi = TMath::Pi();
22 
23  // NOTE: All results for a changing "pleptonIsUsed" and "thetaIsUsed" give the same results within
24  // about 1%. Here, both are set to 1 because this is closest to what the Fortran code does.
25  // To use this in comparison, it is important to use the same parameters as d4sigma_calc.C !
26 
27  const int type = 2; // lepton: 1=electron, 2=muon, 3=tau
28  const int reac = 3; // reaction: 1=NN, 2=NP, 3=PP
29 
30  // Phi_kq will run with energy from phikqmin to phikqmax (can be constant)
31  double phikqmin, phikqmax, FRAC; // used to set value of phi_kq [pi]
32  if (varPhiKqIsUsed) {
33  phikqmin = 0.50*pi;
34  phikqmax = 0.55*pi;
35  } else {
36  printf("Please enter phi_kq as a fraction of pi: ");
37  scanf("%lf", &FRAC);
38  printf("Phi_kq set to fixed value %.2lf*pi...\n", FRAC);
39  phikqmin = FRAC*pi;
40  phikqmax = FRAC*pi;
41  }
42 
43  double Mlep = 0.; // lepton mass [GeV]
44  if (type==1) Mlep = 0.510998928e-3;
45  else if (type==2) Mlep = 0.1056583715;
46  else if (type==3) Mlep = 1.77682;
47  else {std::cout<<"ERROR: Invalid type!"<<std::endl; return;}
48 
49  double Mkaon = 0.; // kaon mass [GeV]
50  if (reac==1) Mkaon = 0.493677;
51  else if (reac==2) Mkaon = 0.497614;
52  else if (reac==3) Mkaon = 0.493677;
53  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
54 
55  double Enu = 0.; // neutrino energy [GeV]
56  singlekaon_xsec *inst = new singlekaon_xsec();
57  inst->init(Enu, type, reac);
58 
59  double Emin = inst->threshold;
60  double Emax = 3.1;
61  int Esteps = (int)(20.0*(Emax-Emin));
62 
63  // Initialisation of variables
64  int i,j,k,itenu;
65  const int nsteps2 = 10*COMP, nsteps1 = 10*COMP, nsteps3 = 10*COMP;
66  double varmin1, varmin2, varmin3, varmax1, varmax2, varmax3, temp;
67  double binsvar2[nsteps2+1], binsvar1[nsteps1+1], binsvar3[nsteps3+1];
68  double varvals2[nsteps2], varvals1[nsteps1], varvals3[nsteps3];
69 
70  double Tlep, Elep; // LEPTON ENERGY [GeV]
71  double Tkaon, Ekaon; // KAON ENERGY [GeV]
72  double theta; // LEPTON ANGLE [rad] - FIXED VALUE IN 3D CASE !!!
73  double phikq; // KAON ANGLE [rad]
74 
75  double diff3sigma;
76  double diff2sigma;
77  double diff1sigma;
78  double sigma; // CROSS SECTION [cm^2]
79 
81  std::cout << " E [GeV]\txsec [cm^2]" << std::endl;
82  std::cout << "-----------------------------" << std::endl;
83 
84  std::string fpath = "./data/";
86  if (varPhiKqIsUsed) fname = "d3sigma_phi=var.out";
87  else fname = Form("d3sigma_phi=%.2lfpi.out",FRAC);
88  FILE *outfile = fopen((fpath+fname).c_str(),"w");
89  fprintf(outfile, "# E [GeV]\txsec [cm^2]\n");
90  fprintf(outfile, "#-----------------------------\n");
91 
92  for (itenu=0; itenu<=Esteps; itenu++) {
93 
94  // Initialisation of energy
95  Enu = Emin + itenu*0.05;
96  inst->init(Enu, type, reac);
97 
98  // Value of phi_kq is a linear function of energy (3D case)
99  phikq = phikqmin + (double)itenu/Esteps*(phikqmax-phikqmin);
100 
101  // Integration ranges
102  varmin1 = 0.; varmax1 = Enu-Mkaon-Mlep; // Tkaon
103  // varmin2, varmax2 filled in loop below
104  if (thetaIsUsed) { // theta
105  varmin3 = 0.;
106  varmax3 = pi;
107  } else { // cos(theta)
108  varmin3 = -1.0;
109  varmax3 = 1.0;
110  }
111 
112  // Calculate edges of bins
113  for (i=0; i<=nsteps1; i++) {
114  binsvar1[i] = varmin1 + i*(varmax1-varmin1)/nsteps1;
115  }
116  // binsvar2 filled in loop below
117  for (k=0; k<=nsteps3; k++) {
118  if (thetaIsUsed) {
119  temp = varmax3 - k*(varmax3-varmin3)/nsteps3; // theta
120  binsvar3[k] = cos(temp); // cos(theta)
121  } else {
122  binsvar3[k] = varmin3 + k*(varmax3-varmin3)/nsteps3;
123  }
124  if (verbose && !itenu && k)
125  std::cout << "Bin #" << k << " has cos(theta) = " << 0.5*(binsvar3[k]+binsvar3[k-1])
126  << ",\t width = " << binsvar3[k]-binsvar3[k-1] << std::endl;
127  }
128 
129  // Calculate central bin values (average of left and right edge of bin)
130  for (i=0; i<nsteps1; i++) varvals1[i] = 0.5*(binsvar1[i+1]+binsvar1[i]);
131  // binsvar2 filled in loop below
132  for (k=0; k<nsteps3; k++) varvals3[k] = 0.5*(binsvar3[k+1]+binsvar3[k]);
133 
134  sigma = 0.; // reset integral over Tkaon
135 
136  // **********************
137  // Integrate over Tkaon
138  // **********************
139  for (i=0; i<nsteps1; i++) {
140  if (binsvar1[i+1]==binsvar1[i]) continue;
141  Tkaon = varvals1[i];
142  Ekaon = Tkaon+Mkaon;
143 
144  // ----------------------------------------------
145  // Define binning for second integration variable
146  // ----------------------------------------------
147  if (pLeptonIsUsed) { // plep
148  varmin2 = 0.;
149  varmax2 = sqrt((Enu-Ekaon)*(Enu-Ekaon)-Mlep*Mlep);
150  } else { // Tlep
151  varmin2 = 0.;
152  varmax2 = Enu-Ekaon-Mlep;
153  }
154  for (j=0; j<=nsteps2; j++) {
155  temp = varmin2 + j*(varmax2-varmin2)/nsteps2; // plep OR Tlep
156  if (pLeptonIsUsed) {
157  binsvar2[j] = sqrt(temp*temp+Mlep*Mlep)-Mlep; // plep -> Tlep
158  } else {
159  binsvar2[j] = temp; // Tlep
160  }
161  if (verbose && !pLeptonIsUsed && !itenu && j)
162  std::cout << "Bin edge #" << j << " has Tlep = " << binsvar2[j]
163  << ",\t width = " << binsvar2[j]-binsvar2[j-1] << std::endl;
164  }
165  for (j=0; j<nsteps2; j++) varvals2[j] = 0.5*(binsvar2[j+1]+binsvar2[j]);
166  // ----------------------------------------------
167 
168  diff1sigma = 0.; // reset integral over plep
169 
170  // *******************************
171  // Integrate over plep (or Tlep)
172  // *******************************
173  for (j=0; j<nsteps2; j++) {
174  if (binsvar2[j+1]==binsvar2[j]) continue;
175  Tlep = varvals2[j];
176  Elep = Tlep+Mlep;
177 
178  diff2sigma = 0.; // reset integral over theta
179 
180  // ************************************
181  // Integrate over theta (or costheta)
182  // ************************************
183  for (k=0; k<nsteps3; k++) {
184  if (binsvar3[k+1]==binsvar3[k]) continue;
185  theta = acos(varvals3[k]);
186 
187  // Calculate 3D differential xsec
188  diff3sigma = inst->diffxsec(Tlep,Tkaon,theta,phikq);
189  diff3sigma *= pow(2*pi,2);
190 
191  // This Jacobian no longer needed, since binning is adjusted to cos(theta)
192  // diff3sigma *= sin(theta);
193 
194  diff2sigma += diff3sigma*(binsvar3[k+1]-binsvar3[k]);
195  }
196 
197  // Multiplication with Jacobian for transformation plep->Tlep (it is E/p)
198  diff2sigma *= Elep/sqrt(Elep*Elep-Mlep*Mlep);
199 
200  diff1sigma += diff2sigma*(binsvar2[j+1]-binsvar2[j]);
201  }
202 
203  sigma += diff1sigma*(binsvar1[i+1]-binsvar1[i]);
204  }
205 
206  // Total cross section for this energy
207  std::cout << Form("%12.6lf",Enu) << "\t" << sigma << std::endl;
208  fprintf(outfile, "%lf\t%e\n",Enu,sigma);
209  }
210 
211  fclose(outfile);
212  std::cout << std::endl << "Output written to file: " << fpath << fname << std::endl << std::endl;
213 }
214 
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
T acos(T number)
Definition: d0nt_math.hpp:54
void d3sigma_calc()
Definition: d3sigma_calc.C:11
fclose(fg1)
printf("%d Experimental points found\n", nlines)
double diffxsec(double Tlep, double Tkaon, double theta, double phikq)
const double j
Definition: BetheBloch.cxx:29
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
T cos(T number)
Definition: d0nt_math.hpp:78
double Emax
void init(double Etot, int type, int reac)
FILE * outfile
Definition: dump_event.C:13
enum BeamMode string