piN_asymmetric_decay.C
Go to the documentation of this file.
1 //
2 // test asymmetric decays of low multiplicity (2-body pi+N) hadronic systems
3 // using data (spherical harmonic expansion) from Nucl.Phys.B343(1990) 285-309
4 //
5 // To run:
6 // root[0] .L piN_asymmetric_decay.C++
7 // root[1] plot_input_data(); //or
8 // root[1] run_evgen_test(W);
9 //
10 // Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
11 // University of Liverpool & STFC Rutherford Appleton Lab
12 //
13 
14 //
15 // Refresher on spherical harmonics:
16 //
17 // Y^{m}_{l} := Yml = N * exp(i*m*phi) * Pml(costh)
18 //
19 // Y01 = (1/2) * sqrt(3/pi) * costh
20 // Y02 = (1/4) * sqrt(5/pi) * (3*cos^2th - 1)
21 // Y11 = -(1/2) * sqrt(3/2pi) * sinth * exp(i*phi)
22 // Y12 = -(1/2) * sqrt(15/2pi) * sinth*costh * exp(i*phi)
23 // Y22 = (1/4) * sqrt(15/2pi) * sin^2th * exp(i*2*phi)
24 //
25 
26 #include <TMath.h>
27 #include <TH1D.h>
28 #include <TH2D.h>
29 #include <TGraphAsymmErrors.h>
30 #include <TFile.h>
31 #include <TNtuple.h>
32 #include <TSpline.h>
33 #include <TRandom3.h>
34 #include <TGenPhaseSpace.h>
35 #include <TLorentzVector.h>
36 
37 // p pi+
38 const int kNWBins = 6;
39 double Wmin [kNWBins] = { 1.100, 1.180, 1.230, 1.280, 1.500, 1.700 };
40 double Wmax [kNWBins] = { 1.180, 1.230, 1.280, 1.500, 1.700, 2.000 };
41 double Wc [kNWBins] = { 1.140, 1.205, 1.255, 1.390, 1.600, 1.850 };
42 
43 double CoeffReY01 [kNWBins] = { -0.22, -0.01, 0.29, 0.32, 0.61, 0.84 };
44 double CoeffReY02 [kNWBins] = { -0.06, -0.13, -0.19, -0.08, 0.49, 0.80 };
45 double CoeffReY11 [kNWBins] = { 0.21, -0.02, -0.20, -0.05, 0.12, -0.03 };
46 double CoeffImY11 [kNWBins] = { -0.09, -0.06, -0.06, 0.01, -0.08, -0.17 };
47 double CoeffReY12 [kNWBins] = { 0.26, 0.05, 0.02, 0.00, 0.09, -0.05 };
48 double CoeffImY12 [kNWBins] = { 0.07, 0.00, 0.16, -0.01, -0.04, -0.03 };
49 double CoeffReY22 [kNWBins] = { -0.04, 0.08, 0.08, 0.04, 0.09, -0.08 };
50 double CoeffImY22 [kNWBins] = { -0.07, -0.07, -0.01, 0.06, -0.15, -0.03 };
51 
52 double dCoeffReY01 [kNWBins] = { 0.11, 0.07, 0.09, 0.09, 0.13, 0.13 };
53 double dCoeffReY02 [kNWBins] = { 0.12, 0.07, 0.10, 0.09, 0.15, 0.14 };
54 double dCoeffReY11 [kNWBins] = { 0.09, 0.06, 0.08, 0.06, 0.08, 0.07 };
55 double dCoeffImY11 [kNWBins] = { 0.08, 0.06, 0.08, 0.06, 0.09, 0.09 };
56 double dCoeffReY12 [kNWBins] = { 0.08, 0.06, 0.08, 0.06, 0.07, 0.06 };
57 double dCoeffImY12 [kNWBins] = { 0.07, 0.05, 0.08, 0.06, 0.08, 0.07 };
58 double dCoeffReY22 [kNWBins] = { 0.07, 0.05, 0.07, 0.06, 0.09, 0.09 };
59 double dCoeffImY22 [kNWBins] = { 0.08, 0.05, 0.08, 0.06, 0.08, 0.07 };
60 
61 // cubic splines (built from data above)
62 //
63 TSpline3 splCoeffReY01("splCoeffReY01",Wc,CoeffReY01,kNWBins,"",0,0);
64 TSpline3 splCoeffReY02("splCoeffReY02",Wc,CoeffReY02,kNWBins,"",1,1);
65 TSpline3 splCoeffReY11("splCoeffReY11",Wc,CoeffReY11,kNWBins,"",1,1);
66 TSpline3 splCoeffImY11("splCoeffImY11",Wc,CoeffImY11,kNWBins,"",1,1);
67 TSpline3 splCoeffReY12("splCoeffReY12",Wc,CoeffReY12,kNWBins,"",1,1);
68 TSpline3 splCoeffImY12("splCoeffImY12",Wc,CoeffImY12,kNWBins,"",1,1);
69 TSpline3 splCoeffReY22("splCoeffReY22",Wc,CoeffReY22,kNWBins,"",1,1);
70 TSpline3 splCoeffImY22("splCoeffImY22",Wc,CoeffImY22,kNWBins,"",1,1);
71 
72 // consts
73 //
74 const double kPi = 3.14;
75 
76 // prototypes
77 //
78 double angular_distribution_pdf(double W, double th, double phi);
79 double coeff (double W, int m, int l, bool real);
80 double Y (double th, double phi, int m, int l, bool real);
81 double Y00 (double th, double phi, bool real);
82 double Y01 (double th, double phi, bool real);
83 double Y02 (double th, double phi, bool real);
84 double Y11 (double th, double phi, bool real);
85 double Y12 (double th, double phi, bool real);
86 double Y22 (double th, double phi, bool real);
87 
88 //_______________________________________________________________________________
89 void run_evgen_test(double W)
90 {
91 // generate events with both the old & new scheme for comparison
92 //
93  TFile f("./events.root","recreate");
94  TNtuple nt("nt","","m:i:th:phi");
95 
96  TRandom rndgen(18928928);
97 
98  const int nev = 10000000;
99 
100  double th_min = 0;
101  double th_max = kPi;
102  double phi_min = 0;
103  double phi_max = 2*kPi;
104 
105 // scan for max
106 //
107  const int nth = 180;
108  const int nphi = 360;
109  double dth = (th_max - th_min ) / (nth -1);
110  double dphi = (phi_max - phi_min) / (nphi-1);
111 
112  double pdf_max = -1;
113  for(int i=0; i<nth; i++)
114  for(int j=0; j<nphi; j++)
115  pdf_max = TMath::Max(pdf_max,
116  angular_distribution_pdf(W,th_min+i*dth,phi_min+j*dphi));
117  pdf_max *= 1.2;
118 
119 // generate N theta,phi pairs using the new method
120 //
121  for(int iev=0; iev<nev; iev++) {
122  bool generated=false;
123  while(!generated) {
124  double th = rndgen.Rndm() * kPi;
125  double phi = rndgen.Rndm() * 2*kPi;
126  double t = rndgen.Rndm() * pdf_max;
127  double pdf = angular_distribution_pdf(W,th,phi);
128  generated = (t<pdf);
129  if(generated) { nt.Fill(1,iev,th,phi); }
130  }
131  }
132 
133 // generate N theta,phi pairs using the old method (phase space decay)
134 //
135  const int Np = 2;
136  const double MN = 0.938;
137  const double Mpi = 0.139;
138  TLorentzVector p4(0.0, 0.0, 0.0, W);
139  Double_t masses[Np] = { MN, Mpi } ;
140 
141  TGenPhaseSpace phspg;
142  Bool_t is_allowed = phspg.SetDecay(p4, Np, masses);
143  assert(is_allowed);
144 
145  double max_wght = -1;
146  for (Int_t iev=0; iev<200; iev++) {
147  double weight = phspg.Generate();
148  max_wght = TMath::Max(max_wght, weight);
149  }
150  max_wght *= 1.2;
151 
152  for(int iev=0; iev<nev; iev++) {
153  bool generated=false;
154  while(!generated) {
155  double weight = phspg.Generate();
156  double t = rndgen.Rndm() * max_wght;
157  generated = (t<weight);
158  if(generated) {
159  TLorentzVector * p4pi = phspg.GetDecay(1);
160  double th = p4pi->Theta();
161  double phi = p4pi->Phi() + kPi;
162  nt.Fill(-1,iev,th,phi);
163  }
164  }
165  }
166 
167  nt.Write();
168  f.Close();
169 }
170 //_______________________________________________________________________________
171 void plot_input_data(void)
172 {
173 // save all distributions that are input to the new 2-body decay scheme
174 //
175  TFile f("./inputs.root","recreate");
176  TNtuple nt("nt","","th:phi:W:sum:r01:r02:r11:i11:r12:i12:r22:i22");
177 
178  const int nw = 11;
179  const int nth = 180;
180  const int nphi = 360;
181 
182  double th_min = 0;
183  double th_max = kPi;
184  double phi_min = 0;
185  double phi_max = 2*kPi;
186  double dth = (th_max - th_min ) / (nth -1);
187  double dphi = (phi_max - phi_min) / (nphi-1);
188 
189  double W[nw] = { 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0 };
190 
191  for(int i=0; i<nth; i++) {
192  double th = th_min + i*dth;
193  for(int j=0; j<nphi; j++) {
194  double phi = phi_min + j*dphi;
195  for(int k=0; k<nw; k++) {
196  double w = W[k];
197 
198  double sum = angular_distribution_pdf(w, th, phi);
199  double r01 = coeff(w, 0,1,true ) * Y(th,phi, 0,1,true );
200  double r02 = coeff(w, 0,2,true ) * Y(th,phi, 0,2,true );
201  double r11 = coeff(w, 1,1,true ) * Y(th,phi, 1,1,true );
202  double i11 = coeff(w, 1,1,false) * Y(th,phi, 1,1,false);
203  double r12 = coeff(w, 1,2,true ) * Y(th,phi, 1,2,true );
204  double i12 = coeff(w, 1,2,false) * Y(th,phi, 1,2,false);
205  double r22 = coeff(w, 2,2,true ) * Y(th,phi, 2,2,true );
206  double i22 = coeff(w, 2,2,false) * Y(th,phi, 2,2,false);
207 
208  nt.Fill(th,phi,w,sum,r01,r02,r11,i11,r12,i12,r22,i22);
209  }//w
210  }//phi
211  }//theta
212 
213  double wmin = 1.1;
214  double wmax = 2.0;
215  double nwb = 100;
216  TH1D cr01("cr01","coefficients of R{Y01} vs W",nwb,wmin,wmax);
217  TH1D cr02("cr02","coefficients of R{Y02} vs W",nwb,wmin,wmax);
218  TH1D cr11("cr11","coefficients of R{Y11} vs W",nwb,wmin,wmax);
219  TH1D ci11("ci11","coefficients of I{Y11} vs W",nwb,wmin,wmax);
220  TH1D cr12("cr12","coefficients of R{Y12} vs W",nwb,wmin,wmax);
221  TH1D ci12("ci12","coefficients of I{Y12} vs W",nwb,wmin,wmax);
222  TH1D cr22("cr22","coefficients of R{Y22} vs W",nwb,wmin,wmax);
223  TH1D ci22("ci22","coefficients of I{Y22} vs W",nwb,wmin,wmax);
224 
225  for(int i=1; i<cr01.GetNbinsX(); i++) {
226  double w = cr01.GetBinCenter(i);
227  cr01.Fill( w, coeff(w,0,1,true) );
228  cr02.Fill( w, coeff(w,0,2,true) );
229  cr11.Fill( w, coeff(w,1,1,true) );
230  ci11.Fill( w, coeff(w,1,1,false) );
231  cr12.Fill( w, coeff(w,1,2,true) );
232  ci12.Fill( w, coeff(w,1,2,false) );
233  cr22.Fill( w, coeff(w,2,2,true) );
234  ci22.Fill( w, coeff(w,2,2,false) );
235  }
236 
237  double dW[kNWBins];
238  for(int i=0; i<kNWBins; i++) {dW[i] = Wmax[i]-Wc[i];}
239 
240  TGraphAsymmErrors dr01(kNWBins,Wc,CoeffReY01,dW,dW,dCoeffReY01,dCoeffReY01);
241  TGraphAsymmErrors dr02(kNWBins,Wc,CoeffReY02,dW,dW,dCoeffReY02,dCoeffReY02);
242  TGraphAsymmErrors dr11(kNWBins,Wc,CoeffReY11,dW,dW,dCoeffReY11,dCoeffReY11);
243  TGraphAsymmErrors di11(kNWBins,Wc,CoeffImY11,dW,dW,dCoeffImY11,dCoeffImY11);
244  TGraphAsymmErrors dr12(kNWBins,Wc,CoeffReY12,dW,dW,dCoeffReY12,dCoeffReY12);
245  TGraphAsymmErrors di12(kNWBins,Wc,CoeffImY12,dW,dW,dCoeffImY12,dCoeffImY12);
246  TGraphAsymmErrors dr22(kNWBins,Wc,CoeffReY22,dW,dW,dCoeffReY22,dCoeffReY22);
247  TGraphAsymmErrors di22(kNWBins,Wc,CoeffImY22,dW,dW,dCoeffImY22,dCoeffImY22);
248 
249  dr01.Write("dr01");
250  dr02.Write("dr02");
251  dr11.Write("dr11");
252  di11.Write("di11");
253  dr12.Write("dr12");
254  di12.Write("di12");
255  dr22.Write("dr22");
256  di22.Write("di22");
257  cr01.Write();
258  cr02.Write();
259  cr11.Write();
260  ci11.Write();
261  cr12.Write();
262  ci12.Write();
263  cr22.Write();
264  ci22.Write();
265  nt.Write();
266  f.Close();
267 }
268 //_______________________________________________________________________________
269 //_______________________________________________________________________________
270 // Code implementing the new 2-body decays
271 //_______________________________________________________________________________
272 //_______________________________________________________________________________
273 double angular_distribution_pdf(double W, double th, double phi)
274 {
275 // return Sum{m,l} { (interpolated coeff) x (spherical harmonic)}
276 //
277  double fr = 1 +
278  coeff(W, 0,1,true ) * Y(th,phi, 0,1,true ) +
279  coeff(W, 0,2,true ) * Y(th,phi, 0,2,true ) +
280  coeff(W, 1,1,true ) * Y(th,phi, 1,1,true ) +
281  coeff(W, 1,2,true ) * Y(th,phi, 1,2,true ) +
282  coeff(W, 2,2,true ) * Y(th,phi, 2,2,true );
283 
284  double fi = coeff(W, 1,1,false) * Y(th,phi, 1,1,false) +
285  coeff(W, 1,2,false) * Y(th,phi, 1,2,false) +
286  coeff(W, 2,2,false) * Y(th,phi, 2,2,false);
287 
288  double f = TMath::Sqrt(fr*fr+fi*fi);
289  return f;
290 }
291 //_______________________________________________________________________________
292 double coeff(double W, int m, int l, bool real)
293 {
294  if (m==0 && l==0 && real) return 1;
295  else if (m==0 && l==1 && real) return splCoeffReY01.Eval(W);
296  else if (m==0 && l==2 && real) return splCoeffReY02.Eval(W);
297  else if (m==1 && l==1 && real) return splCoeffReY11.Eval(W);
298  else if (m==1 && l==1 && !real) return splCoeffImY11.Eval(W);
299  else if (m==1 && l==2 && real) return splCoeffReY12.Eval(W);
300  else if (m==1 && l==2 && !real) return splCoeffImY12.Eval(W);
301  else if (m==2 && l==2 && real) return splCoeffReY22.Eval(W);
302  else if (m==2 && l==2 && !real) return splCoeffImY22.Eval(W);
303  else
304  return 0;
305 }
306 //_______________________________________________________________________________
307 double Y(double th, double phi, int m, int l, bool real)
308 {
309  if (m==0 && l==0) return Y00(th,phi,real);
310  else if (m==0 && l==1) return Y01(th,phi,real);
311  else if (m==0 && l==2) return Y02(th,phi,real);
312  else if (m==1 && l==1) return Y11(th,phi,real);
313  else if (m==1 && l==2) return Y12(th,phi,real);
314  else if (m==2 && l==2) return Y22(th,phi,real);
315  else
316  return 0;
317 }
318 //_______________________________________________________________________________
319 double Y00(double /*th*/, double /*phi*/, bool real)
320 {
321  if(!real) return 0;
322  return 0.5/TMath::Sqrt(kPi);
323 }
324 //_______________________________________________________________________________
325 double Y01(double th, double /*phi*/, bool real)
326 {
327  if(!real) return 0;
328  return 0.5 * TMath::Sqrt(3./kPi) * TMath::Cos(th);
329 }
330 //_______________________________________________________________________________
331 double Y02(double th, double /*phi*/, bool real)
332 {
333  if(!real) return 0;
334  return 0.25 * TMath::Sqrt(5./kPi) * (3.* TMath::Power(TMath::Cos(th),2.) - 1);
335 }
336 //_______________________________________________________________________________
337 double Y11(double th, double phi, bool real)
338 {
339  double a = -0.5 * TMath::Sqrt(1.5/kPi) * TMath::Sin(th);
340  if(!real) return a * TMath::Sin(phi);
341  else return a * TMath::Cos(phi);
342 }
343 //_______________________________________________________________________________
344 double Y12(double th, double phi, bool real)
345 {
346  double a = -0.5 * TMath::Sqrt(7.5/kPi) * TMath::Cos(th) * TMath::Sin(th);
347  if(!real) return a * TMath::Sin(phi);
348  else return a * TMath::Cos(phi);
349 }
350 //_______________________________________________________________________________
351 double Y22(double th, double phi, bool real)
352 {
353  double a = 0.25 * TMath::Sqrt(7.5/kPi) * TMath::Power( TMath::Sin(th), 2.);
354  if(!real) return a * TMath::Sin(2*phi);
355  else return a * TMath::Cos(2*phi);
356 }
357 //_______________________________________________________________________________
int iev
Definition: runWimpSim.h:118
void run_evgen_test(double W)
const double kPi
void plot_input_data(void)
double Y12(double th, double phi, bool real)
const Var weight
TSpline3 splCoeffImY22("splCoeffImY22", Wc, CoeffImY22, kNWBins,"", 1, 1)
double Wc[kNWBins]
double dCoeffImY22[kNWBins]
double Wmax[kNWBins]
double CoeffImY12[kNWBins]
double CoeffReY02[kNWBins]
double dCoeffReY02[kNWBins]
TSpline3 splCoeffImY12("splCoeffImY12", Wc, CoeffImY12, kNWBins,"", 1, 1)
double dCoeffReY22[kNWBins]
double Y22(double th, double phi, bool real)
double dCoeffReY11[kNWBins]
double Y00(double th, double phi, bool real)
double angular_distribution_pdf(double W, double th, double phi)
TSpline3 splCoeffReY11("splCoeffReY11", Wc, CoeffReY11, kNWBins,"", 1, 1)
TSpline3 splCoeffReY12("splCoeffReY12", Wc, CoeffReY12, kNWBins,"", 1, 1)
const double a
double Y01(double th, double phi, bool real)
const int kNWBins
const double j
Definition: BetheBloch.cxx:29
TSpline3 splCoeffReY01("splCoeffReY01", Wc, CoeffReY01, kNWBins,"", 0, 0)
double coeff(double W, int m, int l, bool real)
double CoeffImY11[kNWBins]
double Y(double th, double phi, int m, int l, bool real)
double CoeffReY22[kNWBins]
double dCoeffReY12[kNWBins]
double Y11(double th, double phi, bool real)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double dCoeffImY11[kNWBins]
TSpline3 splCoeffImY11("splCoeffImY11", Wc, CoeffImY11, kNWBins,"", 1, 1)
assert(nhit_max >=nhit_nbins)
double Wmin[kNWBins]
double Y02(double th, double phi, bool real)
TSpline3 splCoeffReY02("splCoeffReY02", Wc, CoeffReY02, kNWBins,"", 1, 1)
double CoeffReY12[kNWBins]
double CoeffReY01[kNWBins]
Double_t sum
Definition: plot.C:31
double CoeffImY22[kNWBins]
TNtuple * nt
Definition: drawXsec.C:2
double dCoeffImY12[kNWBins]
Float_t w
Definition: plot.C:20
#define W(x)
TSpline3 splCoeffReY22("splCoeffReY22", Wc, CoeffReY22, kNWBins,"", 1, 1)
double CoeffReY11[kNWBins]
double dCoeffReY01[kNWBins]