validation.C
Go to the documentation of this file.
1 // *********************************************************************
2 // Validation of single kaon kinematics - Martti Nirkko (28th Nov 2014)
3 // Compile and run in terminal: root -l -b -q validation.C+g
4 // *********************************************************************
5 
6 #include "code/singlekaon_xsec.cxx"
7 #include <TH1D.h>
8 #include <TCanvas.h>
9 #include <TFile.h>
10 #include <TLorentzVector.h>
11 #include <TVector3.h>
12 #include <TMath.h>
13 #include <TStyle.h>
14 
15 void validation() {
16  // Validate kinematical distributions of single kaon production
17  // Have: E_nu, T_l, T_K, cos(theta_l), phi_Kq (generated randomly)
18  // Want: Generate differential xsec histograms for these variables
19 
20  // ***************************
21  // * STEERING PARAMETERS *
22  // ***************************
23  const int COMP = 10; // time complexity (runs with O(n^4)...)
24  const int pLeptonIsUsed = 0; // use momentum instead of kinetic energy
25  const int thetaIsUsed = 0; // use angle instead of cosine
26  const double pi = TMath::Pi();
27 
28  // NOTE: To compare with the output of d4sigma_hist.C, use the same input variables above. Here,
29  // for plotting the differential 1D cross-section, I want to use Tlep and costheta.
30 
31  const int type = 2; // lepton: 1=electron, 2=muon, 3=tau
32  const int reac = 3; // reaction: 1=NN, 2=NP, 3=PP
33 
34  const double NEvents = 1e6; // Number of events to generate
35 
36  // Initialise random number seed
37  srand (time(NULL));
38 
39  double Enu; // neutrino energy [GeV]
40  printf("Please enter neutrino energy: ");
41  scanf("%lf", &Enu);
42  printf("Trying to find input file for E_nu = %3.1lf GeV...\n", Enu);
43  singlekaon_xsec *inst = new singlekaon_xsec();
44  inst->init(Enu, type, reac);
45 
46  // Get output from d4sigma_plot.C to compare histograms
47  std::string fname = Form("data/d4sigma_plot_%3.1lfGeV.root", Enu);
48  TFile *infile = new TFile(fname.c_str());
49  if (!infile) {
50  printf("ERROR: File not found!");
51  return;
52  }
53  TH1D *dTl = (TH1D*)infile->Get("dsigma_dTlepton");
54  TH1D *dTk = (TH1D*)infile->Get("dsigma_dTkaon");
55  TH1D *dct = (TH1D*)infile->Get("dsigma_dcostheta");
56  TH1D *dph = (TH1D*)infile->Get("dsigma_dphikq");
57 
58  std::string strl;
59  double ml = 0.; // lepton mass [GeV]
60  if (type==1) { ml = 0.510998928e-3; strl = "e"; }
61  else if (type==2) { ml = 0.1056583715; strl = "#mu"; }
62  else if (type==3) { ml = 1.77682; strl = "#tau"; }
63  else {std::cout<<"ERROR: Invalid lepton type!"<<std::endl; return;}
64 
65  std::string strK;
66  double mK = 0.; // kaon mass [GeV]
67  if (reac==1) { mK = 0.493677; strK = "K^{+}"; }
68  else if (reac==2) { mK = 0.497614; strK = "K^{0}"; }
69  else if (reac==3) { mK = 0.493677; strK = "K^{+}"; }
70  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
71 
72  std::string strN0, strN1;
73  double mN = 0.; // nucleon mass [GeV]
74  if (reac==1) { mN = 0.939565379; strN0 = "n"; strN1 = "n"; }
75  else if (reac==2) { mN = 0.939565379; strN0 = "n"; strN1 = "p"; }
76  else if (reac==3) { mN = 0.938272046; strN0 = "p"; strN1 = "p"; }
77  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
78 
79  std::string reaction = Form( "%.2lf GeV #nu_{%s} %s #rightarrow %s^{-} %s %s", Enu, strl.c_str(),
80  strN0.c_str(), strl.c_str(), strN1.c_str(), strK.c_str() );
81 
82  // Number of bins
83  const int NBINS = 10*COMP;
84 
85  // Initialize histograms
86  TH1D *hTK = new TH1D("hTK", "Kaon kinetic energy", NBINS, 0., Enu-mK-ml);
87  TH1D *hTl = new TH1D("hTl", "Lepton kinetic energy", NBINS, 0., Enu-mK-ml);
88  TH1D *hct = new TH1D("hct", "Lepton polar angle", NBINS, -1., 1.);
89  TH1D *hph = new TH1D("hph", "Kaon azimuthal angle", NBINS/5, -pi, pi);
90 
91  TH1D *hQ2 = new TH1D("hQ2", "Momentum transfer", NBINS, 0., Enu-mK-ml);
92  TH1D *hTN = new TH1D("hTN", "Nucleon kinetic energy", NBINS, 0., Enu-mK-ml);
93  TH1D *hcl = new TH1D("hcl", "Lepton polar angle", NBINS, -1., 1.);
94  TH1D *hcK = new TH1D("hcK", "Kaon polar angle", NBINS, -1., 1.);
95  TH1D *hcN = new TH1D("hcN", "Nucleon polar angle", NBINS, -1., 1.);
96  TH1D *hpl = new TH1D("hpl", "Lepton azimuthal angle", NBINS, -pi, pi);
97  TH1D *hpK = new TH1D("hpK", "Kaon azimuthal angle", NBINS, -pi, pi);
98  TH1D *hpN = new TH1D("hpN", "Nucleon azimuthal angle", NBINS, -pi, pi);
99 
100  // Histograms for momenta (measured by detector)
101  double maxmom = sqrt(pow(Enu-mK-ml + mN, 2) - mN*mN);
102  TH1D *hmoml = new TH1D("hmoml", "Lepton momentum", NBINS, 0., maxmom);
103  TH1D *hmomK = new TH1D("hmomK", "Kaon momentum", NBINS, 0., maxmom);
104  TH1D *hmomN = new TH1D("hmomN", "Nucleon momentum", NBINS, 0., maxmom);
105  TH1D *hmomlx = new TH1D("hmomlx", "Lepton momentum", NBINS, 0., maxmom);
106  TH1D *hmomKx = new TH1D("hmomKx", "Kaon momentum", NBINS, 0., maxmom);
107  TH1D *hmomNx = new TH1D("hmomNx", "Nucleon momentum", NBINS, 0., maxmom);
108  TH1D *hmomly = new TH1D("hmomly", "Lepton momentum", NBINS, 0., maxmom);
109  TH1D *hmomKy = new TH1D("hmomKy", "Kaon momentum", NBINS, 0., maxmom);
110  TH1D *hmomNy = new TH1D("hmomNy", "Nucleon momentum", NBINS, 0., maxmom);
111  TH1D *hmomlz = new TH1D("hmomlz", "Lepton momentum", NBINS, 0., maxmom);
112  TH1D *hmomKz = new TH1D("hmomKz", "Kaon momentum", NBINS, 0., maxmom);
113  TH1D *hmomNz = new TH1D("hmomNz", "Nucleon momentum", NBINS, 0., maxmom);
114 
115  // Initialise all needed parameters
116  double TK_max, TK, Tl_max, Tl, pl_max, pl, costh_l, phi_l, phi_Kq, xsec;
117  double El, theta, Q2;
118  TLorentzVector P4_nu, P4_lep, q;
119  TVector3 lep3, p_l, p_Q, p_K, p_N, p_Kq, p_Nq;
120  double TN, EK, pK, pN;
121  double pQ, costh_Kq, theta_Kq, rot_Z, rot_X;
122  double costh_K, costh_N, phi_K, phi_N;
123  int sign;
124 
125  // Unit 3-vectors and neutrino momentum
126  const TVector3 e1(1,0,0);
127  const TVector3 e2(0,1,0);
128  const TVector3 e3(0,0,1);
129  const TVector3 p_nu = Enu*e3;
130 
131  // Generate large number of events
132  for (int i=0; i<(int)NEvents; i++) {
133  xsec = -999.;
134  while(xsec <= 0) { // retry until valid kinematics are found
135  TK_max = Enu - mK - ml; // maximal allowed kaon kinetic energy
136  TK = TK_max*(rand()/(double)RAND_MAX); // kaon kinetic energy [0, TK_max]
137  Tl_max = Enu - mK - ml; // maximal allowed lepton kinetic energy (uncorrelated!)
138  pl_max = sqrt((Tl_max+ml)*(Tl_max+ml)-ml*ml); // maximal allowed lepton momentum
139  if (pLeptonIsUsed) {
140  pl = pl_max*(rand()/(double)RAND_MAX); // lepton momentum [0, pl_max]
141  Tl = sqrt(pl*pl+ml*ml)-ml;
142  } else {
143  Tl = Tl_max*(rand()/(double)RAND_MAX); // lepton kinetic energy [0, Tl_max]
144  }
145  if (thetaIsUsed) {
146  theta = pi*(rand()/(double)RAND_MAX); // lepton polar angle [0, pi]
147  costh_l = cos(theta);
148  } else {
149  costh_l = 2.*(rand()/(double)RAND_MAX)-1.; // cosine of lepton polar angle [-1, 1]
150  theta = acos(costh_l);
151  }
152  phi_l = pi*(2.*(rand()/(double)RAND_MAX)-1.); // lepton azimuthal angle [-pi, pi]
153  phi_Kq = pi*(2.*(rand()/(double)RAND_MAX)-1.); // kaon azimuthal angle [-pi, pi]
154  xsec = inst->diffxsec(Tl,TK,theta,phi_Kq); // 4D-differential cross-section
155  }
156 
157  // INFO: In dsig_dQ2.f, the value for cos(theta) is changed to the following:
158  // cos(theta)=(2.0*Enu*El-aml*aml+aq2)/(2.0*Enu*alepvec);
159  // aq2 is the (negative) integration variable for Q^2, going from 0 to -Enu
160 
161  // Calculate total energy and momentum of particles
162  El = Tl + ml; pl = sqrt(El*El - ml*ml);
163  EK = TK + mK; pK = sqrt(EK*EK - mK*mK);
164 
165  // Calculate Q^2 from lepton kinematics
166  P4_nu = TLorentzVector(0.,0.,Enu,Enu);
167  lep3 = TVector3(0.,0.,0.);
168  lep3.SetMagThetaPhi(pl,theta,0.);
169  P4_lep = TLorentzVector(lep3,El);
170  q = P4_nu - P4_lep;
171  Q2 = -q.Mag2(); // --> Q^2 histogram
172 
173  // Calculate 3-momentum of lepton
174  p_l[0] = pl*sin(theta)*cos(phi_l);
175  p_l[1] = pl*sin(theta)*sin(phi_l);
176  p_l[2] = pl*costh_l;
177 
178  // Calculate 3-momentum of momentum transfer to hadronic system
179  p_Q = -p_l;
180  p_Q[2] += Enu;
181  pQ = p_Q.Mag();
182  TVector3 dirQ = p_Q.Unit();
183 
184  // In the momentum transfer plane, new angles... (see eq.17 of notes)
185  costh_Kq = (pQ*pQ + pK*pK + mN*mN - (Enu-El-EK+mN)*(Enu-El-EK+mN))/(2.*pQ*pK);
186  theta_Kq = acos(costh_Kq);
187 
188  // Get rotation angle from z-axis to Q^2-direction
189  rot_Z = acos(e3*(p_Q.Unit())); // rotate counter-clockwise, towards x-axis by [0,pi]
190 
191  // Get rotation angle from x-axis to Q^2-direction (projected to transverse plane)
192  sign = 0;
193  if (p_Q[1] > 0) sign = +1;
194  else sign = -1;
195  TVector3 p_Qt(p_Q[0],p_Q[1],0);
196  rot_X = sign*acos(e1*(p_Qt.Unit())); // rotate around z-axis by [-pi,pi]
197 
198  // Calculate 3-momentum of kaon (in Q-frame)
199  p_Kq[0] = pK*sin(theta_Kq)*cos(phi_Kq);
200  p_Kq[1] = pK*sin(theta_Kq)*sin(phi_Kq);
201  p_Kq[2] = pK*costh_Kq;
202 
203  // Calculate 3-momentum of nucleon (in Q-frame)
204  p_Nq = -p_Kq;
205  p_Nq[2] += pQ;
206 
207  // Rotate particles into correct frame
208  p_K = p_Kq;
209  p_K.RotateY(rot_Z);
210  p_K.RotateZ(rot_X);
211  p_N = p_Nq;
212  p_N.RotateY(rot_Z);
213  p_N.RotateZ(rot_X);
214 
215  // Kinetic energy of nucleon
216  pN = p_N.Mag();
217  TN = sqrt(pN*pN+mN*mN) - mN; // --> TN histogram
218 
219  // Momenta of particles
220  hmoml->Fill(pl, xsec);
221  hmomK->Fill(pK, xsec);
222  hmomN->Fill(pN, xsec);
223  hmomlx->Fill(p_l[0], xsec);
224  hmomKx->Fill(p_K[0], xsec);
225  hmomNx->Fill(p_N[0], xsec);
226  hmomly->Fill(p_l[1], xsec);
227  hmomKy->Fill(p_K[1], xsec);
228  hmomNy->Fill(p_N[1], xsec);
229  hmomlz->Fill(p_l[2], xsec);
230  hmomKz->Fill(p_K[2], xsec);
231  hmomNz->Fill(p_N[2], xsec);
232 
233  // Calculate polar angles w.r.t. neutrino
234  costh_l = p_l.CosTheta();
235  costh_K = p_K.CosTheta();
236  costh_N = p_N.CosTheta();
237 
238  // CAUTION: TESTING!!!
239  // Rotate to lepton frame
240  p_l.RotateUz(p_l.Unit());
241  p_K.RotateUz(p_l.Unit());
242  p_N.RotateUz(p_l.Unit());
243 
244  // Calculate azimuthal angles w.r.t. neutrino (z-axis)
245  phi_l = p_l.Phi() - p_l.Phi();
246  phi_K = p_K.Phi() - p_l.Phi();
247  phi_N = p_N.Phi() - p_l.Phi();
248 
249  // Multiply Jacobians, if needed
250  //if (thetaIsUsed) xsec *= sin(theta); // costh -> theta
251  if (!pLeptonIsUsed) xsec *= El/pl; // pl -> Tlep
252 
253  // Weight entries by xsec
254  hTK->Fill(TK, xsec);
255  hTl->Fill(Tl, xsec);
256  hct->Fill(costh_l, xsec);
257  hph->Fill(phi_Kq, xsec);
258 
259  hQ2->Fill(Q2, xsec);
260  hTN->Fill(TN, xsec);
261  hcl->Fill(costh_l, xsec);
262  hcK->Fill(costh_K, xsec);
263  hcN->Fill(costh_N, xsec);
264  hpl->Fill(phi_l, xsec); // should always be zero
265  hpK->Fill(phi_K, xsec);
266  hpN->Fill(phi_N, xsec);
267 
268  if (i%((int)NEvents/20)==0) cout << "Processing... (" << i/((int)NEvents/100) << "%)" << endl;
269  }
270 
271  // Scale histograms to cross-section
272  double totalxsec = dTk->Integral();
273  hTK->Scale(totalxsec/hTK->Integral());
274  hTl->Scale(totalxsec/hTl->Integral());
275  hct->Scale(totalxsec/hct->Integral());
276  hph->Scale(totalxsec/hph->Integral());
277  hQ2->Scale(totalxsec/hQ2->Integral());
278  hTN->Scale(totalxsec/hTN->Integral());
279  hcl->Scale(totalxsec/hcl->Integral());
280  hcK->Scale(totalxsec/hcK->Integral());
281  hcN->Scale(totalxsec/hcN->Integral());
282  hpl->Scale(totalxsec/hpl->Integral());
283  hpK->Scale(totalxsec/hpK->Integral());
284  hpN->Scale(totalxsec/hpN->Integral());
285 
286  // Make lines thicker & colourful
287  hTK->SetLineWidth(2); hTK->SetLineColor(kRed);
288  hTl->SetLineWidth(2); hTl->SetLineColor(kRed);
289  hct->SetLineWidth(2); hct->SetLineColor(kRed);
290  hph->SetLineWidth(2); hph->SetLineColor(kRed);
291  hQ2->SetLineWidth(2); //hQ2->SetLineColor(kBlue);
292  hTN->SetLineWidth(2); hTN->SetLineColor(kOrange);
293  hcl->SetLineWidth(2); hcl->SetLineColor(kRed);
294  hcK->SetLineWidth(2); hcK->SetLineColor(kBlack);
295  hcN->SetLineWidth(2); hcN->SetLineColor(kOrange);
296  hpl->SetLineWidth(2); hpl->SetLineColor(kRed);
297  hpK->SetLineWidth(2); hpK->SetLineColor(kBlack);
298  hpN->SetLineWidth(2); hpN->SetLineColor(kOrange);
299  hmoml->SetLineWidth(2); hmoml->SetLineColor(kRed);
300  hmomK->SetLineWidth(2); hmomK->SetLineColor(kBlack);
301  hmomN->SetLineWidth(2); hmomN->SetLineColor(kOrange);
302 
303  hmomlx->SetLineWidth(2); hmomlx->SetLineColor(kRed);
304  hmomKx->SetLineWidth(2); hmomKx->SetLineColor(kBlack);
305  hmomNx->SetLineWidth(2); hmomNx->SetLineColor(kOrange);
306  hmomly->SetLineWidth(2); hmomly->SetLineColor(kRed);
307  hmomKy->SetLineWidth(2); hmomKy->SetLineColor(kBlack);
308  hmomNy->SetLineWidth(2); hmomNy->SetLineColor(kOrange);
309  hmomlz->SetLineWidth(2); hmomlz->SetLineColor(kRed);
310  hmomKz->SetLineWidth(2); hmomKz->SetLineColor(kBlack);
311  hmomNz->SetLineWidth(2); hmomNz->SetLineColor(kOrange);
312 
313  // Play with stat box
314  gStyle->SetOptStat(0);
315  //gStyle->SetStatX(0.5);
316  //gStyle->SetStatY(0.88);
317 
318  // Set logarithmic scale on/off
319  int logscale = 0;
320  if(!logscale) {
321  hTK->SetMinimum(0);
322  hTl->SetMinimum(0);
323  hct->SetMinimum(0);
324  hph->SetMinimum(0);
325 
326  hQ2->SetMinimum(0);
327  hTN->SetMinimum(0);
328  hcl->SetMinimum(0); hcK->SetMinimum(0); hcN->SetMinimum(0);
329  hpl->SetMinimum(0); hpK->SetMinimum(0); hpN->SetMinimum(0);
330 
331  hmoml->SetMinimum(0); hmomK->SetMinimum(0); hmomN->SetMinimum(0);
332  }
333 
334  hTK->SetMaximum(1.2*max(hTK->GetMaximum(),dTk->GetMaximum()));
335  hTl->SetMaximum(1.2*max(hTl->GetMaximum(),dTl->GetMaximum()));
336  hct->SetMaximum(1.2*max(hct->GetMaximum(),dct->GetMaximum()));
337  hph->SetMaximum(1.2*max(hph->GetMaximum(),dph->GetMaximum()));
338 
339  std::string outname = Form("%3.1lfGeV_", Enu);
340 
341  // Plot distributions of primary variables
342  TCanvas *c1 = new TCanvas("c1", "Kinematics of single kaon production", 1200, 900);
343  c1->SetTitle(reaction.c_str());
344  c1->Divide(2,2);
345  TPad* p[4];
346  for (int k=0; k<4; k++) {
347  p[k] = (TPad*)c1->GetPad(k+1);
348  if (logscale) p[k]->SetLogy();
349  }
350  c1->cd(1);
351  hTl->SetXTitle("T_{#mu} [GeV]");
352  hTl->Draw(); dTl->Draw("same");
353  c1->cd(2);
354  hTK->SetXTitle("T_{K} [GeV]");
355  hTK->Draw(); dTk->Draw("same");
356  c1->cd(3);
357  hct->SetXTitle("cos(#theta_{#mu}) [ ]");
358  hct->Draw(); dct->Draw("same");
359  c1->cd(4);
360  hph->SetXTitle("#phi_{Kq} [rad]");
361  hph->Draw(); dph->Draw("same");
362  c1->Print(("images/kinematics_"+outname+"01_input.png").c_str()); c1->Close();
363 
364  // Plot distributions of secondary variables (e.g. Q^2, cos(theta_kaon), ...)
365  TCanvas *c2 = new TCanvas("c2", "Q^2", 1200, 900);
366  c2->SetTitle(reaction.c_str());
367  c2->Divide(2,2);
368  //TPad* p[4];
369  for (int k=0; k<4; k++) {
370  p[k] = (TPad*)c2->GetPad(k+1);
371  if (logscale) p[k]->SetLogy();
372  }
373  c2->cd(1);
374  hQ2->SetTitle("Momentum transfer to nucleon"); hQ2->SetXTitle("Q^{2} [GeV^{2}]");
375  hQ2->SetMaximum(1.2*hQ2->GetMaximum());
376  hQ2->Draw();
377  c2->cd(2);
378  hTl->SetTitle("Kinetic energy of particle"); hTl->SetXTitle("T_{x} [GeV]");
379  hTl->SetLineColor(kRed); hTK->SetLineColor(kBlack);
380  hTl->SetMaximum(max(1.2*hTN->GetMaximum(),max(hTl->GetMaximum(),hTK->GetMaximum())));
381  hTl->Draw(); hTK->Draw("same"); hTN->Draw("same");
382  c2->cd(3);
383  hcl->SetTitle("Polar angle w.r.t. neutrino"); hcl->SetXTitle("cos(#theta_{#nux}) [ ]");
384  hcl->SetMaximum(1.2*max(hcN->GetMaximum(),max(hcl->GetMaximum(),hcK->GetMaximum())));
385  hcl->Draw(); hcK->Draw("same"); hcN->Draw("same");
386  c2->cd(4);
387  hpl->SetTitle("Azimuthal angle w.r.t. lepton"); hpl->SetXTitle("#phi_{lx} [rad]");
388  hpl->SetMaximum(1.2*max(hpK->GetMaximum(),hpN->GetMaximum()));
389  hpl->Draw(); hpK->Draw("same"); hpN->Draw("same");
390  c2->Print(("images/kinematics_"+outname+"02_output.png").c_str()); c2->Close();
391 
392  // Plot momenta of final state particles (measured by the detectors)
393  TCanvas *c3 = new TCanvas("c3", "Momenta in lab frame", 1200, 900);
394  c3->SetTitle(reaction.c_str());
395  c3->Divide(2,2);
396  for (int k=0; k<4; k++) {
397  p[k] = (TPad*)c3->GetPad(k+1);
398  if (logscale) p[k]->SetLogy();
399  }
400  c3->cd(1);
401  hmomlx->SetTitle("X-momentum"); hmomlx->SetXTitle("p_{x} [GeV/c]");
402  hmomlx->SetMaximum(1.2*max(hmomKx->GetMaximum(),hmomlx->GetMaximum()));
403  hmomlx->Draw(); hmomKx->Draw("same"); hmomNx->Draw("same");
404  c3->cd(2);
405  hmomly->SetTitle("Y-momentum"); hmomly->SetXTitle("p_{y} [GeV/c]");
406  hmomly->SetMaximum(1.2*max(hmomKy->GetMaximum(),hmomly->GetMaximum()));
407  hmomly->Draw(); hmomKy->Draw("same"); hmomNy->Draw("same");
408  c3->cd(3);
409  hmomlz->SetTitle("Z-momentum"); hmomlz->SetXTitle("p_{z} [GeV/c]");
410  hmomlz->SetMaximum(1.2*max(hmomKz->GetMaximum(),hmomlz->GetMaximum()));
411  hmomlz->Draw(); hmomKz->Draw("same"); hmomNz->Draw("same");
412  c3->cd(4);
413  hmoml->SetTitle("Total momentum"); hmoml->SetXTitle("|p| [GeV/c]");
414  hmoml->SetMaximum(1.2*max(hmomK->GetMaximum(),hmoml->GetMaximum()));
415  hmoml->Draw(); hmomK->Draw("same"); hmomN->Draw("same");
416  c3->Print(("images/kinematics_"+outname+"03_momenta.png").c_str()); c3->Close();
417 
418  // Generate rootfile for GENIE comparison
419  TFile* outfile = new TFile(("data/kinematics_"+outname+"validation.root").c_str(), "RECREATE");
420  // Kinetic energies of final state particles
421  hTl->Write("T_lepton");
422  hTK->Write("T_kaon");
423  hTN->Write("T_nucleon");
424  // Azimuthal angles w.r.t. neutrino
425  hct->Write("cos(theta_lepton)");
426  hcK->Write("cos(theta_kaon)");
427  hcN->Write("cos(theta_nucleon)");
428  // Radial angles w.r.t. neutrino (z-axis) and lepton (x-axis)
429  hpl->Write("phi_l");
430  hpK->Write("phi_kaon");
431  hpN->Write("phi_nucleon");
432  // Other variables related to momentum transfer to nucleus
433  hph->Write("phi_Kq");
434  hQ2->Write("Q_squared");
435  outfile->Close();
436 
437  if (hTl) delete hTl;
438  if (hTK) delete hTK;
439  if (hTN) delete hTN;
440  if (hct) delete hct;
441  if (hcK) delete hcK;
442  if (hcN) delete hcN;
443  if (hph) delete hph;
444  if (hpl) delete hpl;
445  if (hpK) delete hpK;
446  if (hpN) delete hpN;
447  if (hQ2) delete hQ2;
448  if (hmoml) delete hmoml;
449  if (hmomK) delete hmomK;
450  if (hmomN) delete hmomN;
451  if (hmomlx) delete hmomlx;
452  if (hmomKx) delete hmomKx;
453  if (hmomNx) delete hmomNx;
454  if (hmomly) delete hmomly;
455  if (hmomKy) delete hmomKy;
456  if (hmomNy) delete hmomNy;
457  if (hmomlz) delete hmomlz;
458  if (hmomKz) delete hmomKz;
459  if (hmomNz) delete hmomNz;
460  if (inst) delete inst;
461 }
462 
enum BeamMode kOrange
enum BeamMode kRed
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
const char * p
Definition: xmltok.h:285
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
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
c2
Definition: demo5.py:33
string infile
printf("%d Experimental points found\n", nlines)
double diffxsec(double Tlep, double Tkaon, double theta, double phikq)
Double_t xsec[nknots]
Definition: testXsec.C:47
OStream cout
Definition: OStream.cxx:6
T sin(T number)
Definition: d0nt_math.hpp:132
void validation()
Definition: validation.C:15
T cos(T number)
Definition: d0nt_math.hpp:78
c1
Definition: demo5.py:24
void init(double Etot, int type, int reac)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
FILE * outfile
Definition: dump_event.C:13
def sign(x)
Definition: canMan.py:197
enum BeamMode string