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