GeomAsym.cc
Go to the documentation of this file.
1 
19 #include "GeomAsym.h"
20 
21 #include <Math/PdfFuncMathCore.h>
22 #include <Math/ProbFuncMathCore.h>
23 #include <Math/QuantFuncMathCore.h>
24 
25 using namespace GeomAsymNS;
26 
27 
28 GeomAsym::GeomAsym(int DetectorType)
29 {
30  if (DetectorType < 0 || DetectorType > 2) {
31  printf("GeomAsym: Wrong Detector Type %d \n", DetectorType);
32  exit(1);
33  }
34  fDetectorType = DetectorType;
35  pzMeanLimit = 0.999;
36  pzRMSLimit = 1. - pzMeanLimit;
37  pzIntegralLimit = 1. - 1.e-6;
38 }
39 
40 
41 //------------------------------------------------------------------
42 //----Gamma function ( GSL implementation/ ROOT implementation/...)
43 //-------------------------------------------------------------------
44 double
45 GeomAsym::gamma_cdf(double x, double alpha, double beta)
46 {
47  return ROOT::Math::gamma_cdf(x, alpha, beta, 0.);
48 }
49 
50 
51 double
52 GeomAsym::gamma_quantile(double q, double alpha, double beta)
53 {
54  return ROOT::Math::gamma_quantile(q, alpha, beta);
55 }
56 
57 
58 //------------------------------------------------------------------
59 //----Parameterization of mean Pz
60 //-------------------------------------------------------------------
61 
62 double
63 GeomAsym::GetPzMeanPar(double r, int icomp, int ipar)
64 {
65  if (ipar < 0 || ipar > nPar)
66  return 0;
67 
68  const double* par = 0;
69  if (fDetectorType == 0)
70  par = parPzMean_WCD[ipar][icomp];
71  else if (fDetectorType == 1)
72  par = parPzMean_Scin[ipar][icomp];
73  else if (fDetectorType == 2)
74  par = parPzMean_MD[ipar][icomp];
75  else {
76  printf("unknown combination");
77  exit(1);
78  }
79 
80  const double r_km = r / 1e5;
81  const double val0km = par[0];
82  const double val1km = par[1];
83  const double B = par[2];
84 
85  const double C = (val1km - val0km) / (exp(-B) - 1);
86  const double A = val0km - C;
87 
88  double val = A + C * exp(-B * r_km);
89  if (val >= pzMeanLimit)
90  val = pzMeanLimit;
91  return val;
92 }
93 
94 
95 double
96 GeomAsym::GetPzMean(double DX, double r, int icomp) // DX (g/cm^2) r (cm)
97 {
98  if (icomp < 0 || icomp > 3)
99  return 0;
100  if (fDetectorType == 2 && icomp > 0)
101  return 1;
102 
103  double valDX0 = GetPzMeanPar(r, icomp, 0);
104  double valDX1 = GetPzMeanPar(r, icomp, 1);
105  if (valDX1 < valDX0)
106  return (valDX1 + valDX0) / 2.;
107 
108  double DX0 = DX_Eval[icomp][0];
109  double DX1 = DX_Eval[icomp][1];
110 
111  double X = (DX - DX0) / (DX1 - DX0);
112  double C = (valDX1 - valDX0) / (exp(-1.) - 1.0);
113  double A = valDX0 - C;
114  double pz_mean = A + C * exp(-X);
115 
116  // <pz> is calculated in pz=(0.,1.)
117  if (pz_mean < 0)
118  return 0;
119  if (pz_mean > pzMeanLimit)
120  return pzMeanLimit;
121  return pz_mean;
122 }
123 
124 
125 //------------------------------------------------------------------
126 //----Parameterization of RMS Pz
127 //-------------------------------------------------------------------
128 
129 double
130 GeomAsym::GetPzRMSPar(double r, int icomp, int ipar)
131 {
132  if (ipar < 0 || ipar > nPar)
133  return 0;
134 
135  const double* par = 0;
136  if (fDetectorType == 0)
137  par = parPzRMS_WCD[ipar][icomp];
138  else if (fDetectorType == 1)
139  par = parPzRMS_Scin[ipar][icomp];
140  else
141  par = parPzRMS_MD[ipar][icomp];
142 
143  double r_km = r / 1.e5;
144  double val0km = par[0];
145  double val1km = par[1];
146  double B = par[2];
147 
148  double C = (val1km - val0km) / (exp(-B) - 1.0);
149  double A = val0km - C;
150 
151  double val = A + C * exp(-B * r_km);
152  if (val <= pzRMSLimit)
153  val = pzRMSLimit;
154  return val;
155 }
156 
157 
158 double
159 GeomAsym::GetPzRMS(double DX, double r, int icomp) // DX (g/cm^2) r (cm)
160 {
161  if (icomp < 0 || icomp > 3)
162  return 0;
163  if (fDetectorType == 2 && icomp > 0)
164  return pzRMSLimit;
165 
166  double valDX0 = GetPzRMSPar(r, icomp, 0);
167  double valDX1 = GetPzRMSPar(r, icomp, 1);
168 
169  double DX0 = DX_Eval[icomp][0];
170  double DX1 = DX_Eval[icomp][1];
171 
172  double X = (DX - DX0) / (DX1 - DX0);
173  double C = (valDX1 - valDX0) / (exp(-1.) - 1.0);
174  double A = valDX0 - C;
175  double pz_rms = A + C * exp(-X);
176 
177  // RMS pz is calculated in pz=(0.,1.)
178  if (pz_rms < pzRMSLimit)
179  return pzRMSLimit;
180  if (pz_rms > 1)
181  return 1;
182  return pz_rms;
183 }
184 
185 
186 //--------------------------------------
187 // Amod x Tmod (cos(theta_p),pz,r)
188 //--------------------------------------
189 
190 double
191 GeomAsym::GetTmodAmod_Par_i(double r, int itheta, int icomp, int ipar)
192 {
193  const double r_km = r / 1e5;
194  //-----
195  const double* par = 0;
196  if (ipar == 0) {
197  if (fDetectorType == 0)
198  par = parR0_WCD[icomp][itheta];
199  else if (fDetectorType == 1)
200  par = parR0_Scin[icomp][itheta];
201  else if (fDetectorType == 2)
202  par = parR0_MD[icomp][itheta];
203  else {
204  printf("unknow combination\n");
205  exit(1);
206  }
207  } else if (ipar == 1) {
208  if (fDetectorType == 0)
209  par = parR1_WCD[icomp][itheta];
210  else if (fDetectorType == 1)
211  par = parR1_Scin[icomp][itheta];
212  else if (fDetectorType == 2)
213  par = parR1_MD[icomp][itheta];
214  else {
215  printf("unknown combination\n");
216  exit(1);
217  }
218  } else if (ipar == 2) {
219  if (fDetectorType == 0)
220  par = parR2_WCD[icomp][itheta];
221  else if (fDetectorType == 1)
222  par = parR2_Scin[icomp][itheta];
223  else if (fDetectorType == 2)
224  par = parR2_MD[icomp][itheta];
225  else {
226  printf("unknown combination\n");
227  exit(1);
228  }
229  }
230 
231  //-----
232 
233  double val = 0;
234  const bool UseCte = ((fDetectorType == 0 && ipar == 2 && icomp > 0) ||
235  (fDetectorType == 2 && ipar == 1 && icomp == 0) ||
236  (fDetectorType == 1 && ipar == 0 && !(icomp > 0 && itheta == 0)) ||
237  (fDetectorType == 1 && ipar == 1 && icomp == 0) ||
238  (fDetectorType == 1 && ipar == 2 && !(icomp > 0 && itheta == 0)));
239 
240  if (UseCte)
241  val = par[0];
242  else {
243  double y0 = par[0];
244  double y1 = par[1];
245  double B = par[2];
246 
247  double C = (y1 - y0) / (exp(-B) - 1.0);
248  double A = y0 - C;
249 
250  val = A + C * exp(-B * r_km);
251  }
252  return val;
253 }
254 
255 
256 double
257 GeomAsym::GetTmodAmod(double ctheta_p, double* par, int icomp)
258 {
259  double cthetaEval[2] = { 0.3, 0.6 };
260 
261  double x = 1 - ctheta_p;
262  double x0 = 1 - cthetaEval[0];
263  double x1 = 1 - cthetaEval[1];
264 
265  double y0 = par[0];
266  double val;
267 
268  if (icomp == 0 && fDetectorType != 2) {
269  val = 1 + x * (y0 - 1) / x0;
270  if (val < 0)
271  val = 0;
272  } else {
273  double y1 = par[1];
274 
275  double b = ((y1 - 1) / x1 - (y0 - 1) / x0) / (x1 - x0);
276  double a = (y0 - 1 - x0 * x0 * b) / x0;
277 
278  val = 1 + x * a + x * x * b;
279  if (val < 0)
280  val = 0;
281  }
282  return val;
283 }
284 
285 
286 double
287 GeomAsym::GetTmodAmod(double ctheta_p, double pz, double r, int icomp)
288 {
289  if (fDetectorType == 2 && icomp > 0)
290  return 1.; // TmodAmod not calculated
291  if (pz > 1.000)
292  pz = 1.000;
293  if (ctheta_p < 0.)
294  return 0.;
295 
296  double par[2];
297  for (int itheta = 0; itheta < 2; itheta++) {
298  if (icomp == 0 && itheta == 1 && fDetectorType != 2) {
299  par[itheta] = 0;
300  continue;
301  }
302 
303  double x = (PzEval[icomp][0] - pz) / (PzEval[icomp][0] - PzEval[icomp][1]);
304 
305  double p0 = GetTmodAmod_Par_i(r, itheta, icomp, 0);
306  double p1 = GetTmodAmod_Par_i(r, itheta, icomp, 1);
307  double p2 = GetTmodAmod_Par_i(r, itheta, icomp, 2);
308  if (p2 < 0.01) {
309  p0 = (p0 + p1) / 2.; //cte
310  p1 = p0;
311  p2 = 0.01;
312  }
313 
314  double B = p2;
315  double C = (p1 - p0) / (exp(-B) - 1.0);
316  double A = p0 - C;
317 
318  par[itheta] = A + C * exp(-B * x);
319  }
320 
321  double TmodAmod = GetTmodAmod(ctheta_p, par, icomp);
322  if (fDetectorType == 2)
323  TmodAmod *= ctheta_p;
324  return TmodAmod;
325 }
326 
327 
328 //------------------------------------------------------------------------------------------------------------
329 // Utilities to calculate the truncated <pz> and the truncated dS0dpz integral
330 //------------------------------------------------------------------------------------------------------------
331 
332 double
333 GeomAsym::GetARad(double r, int icomp)
334 {
335  if (fDetectorType == 2 && icomp > 0)
336  return 1; // Not calculated
337 
338  double r_km = r / 1e5;
339 
340  const double* par = 0;
341  if (fDetectorType == 0)
342  par = ARadPar_WCD[icomp];
343  else if (fDetectorType == 1)
344  par = ARadPar_Scin[icomp];
345  else if (fDetectorType == 2)
346  par = ARadPar_MD[icomp];
347  else {
348  printf("unknown combination\n");
349  exit(1);
350  }
351 
352  double y0 = par[0];
353  double y1 = par[1];
354  double B = par[2];
355 
356  double C = (y1 - y0) / (exp(-B) - 1);
357  double A = y0 - C;
358 
359  double val = A + C * exp(-B * r_km);
360  return val;
361 }
362 
363 
364 double
365 GeomAsym::GetCosTheta_p(double pz, double r, double theta, double psi, int icomp)
366 {
367  if (pz > 1 || pz < -1)
368  return 0;
369  const double ARad = GetARad(r, icomp);
370  const double cosTheta_p = pz * cos(theta) + ARad * sqrt(1 - pz * pz) * sin(theta) * cos(psi);
371  if (cosTheta_p < -1)
372  return -1;
373  if (cosTheta_p > 1)
374  return 1;
375  return cosTheta_p;
376 }
377 
378 
379 double
380 GeomAsym::GetPzCut(double r, double theta, double psi, int icomp)
381 {
382  double cut = ctheta_pCut[fDetectorType];
383  if (fabs(psi - M_PI / 2) < 1e-5 * 180 / M_PI)
384  return cut / cos(theta);
385 
386  const double ARad = GetARad(r, icomp);
387  const double alpha = sin(theta) * cos(psi) * ARad;
388  const double a = cos(theta) * cos(theta) + alpha * alpha;
389  const double b = -2*cut * cos(theta);
390  const double c = (cut * cut - alpha * alpha);
391 
392  double d = b * b - 4*a * c;
393  if (d < 0)
394  d = 0;
395  double pz_cut;
396  if (cos(psi) > 0)
397  pz_cut = (-b - sqrt(d)) / 2 / a;
398  else
399  pz_cut = (-b + sqrt(d)) / 2 / a;
400 
401  if (pz_cut < -1)
402  return -1;
403  if (pz_cut > 1)
404  return 1;
405  return pz_cut;
406 }
407 
408 
409 //------------------------------------------------------------------------------------------------------------
410 // Utilities to calculate the truncated <pz> and the truncated dS0dpz integral
411 //------------------------------------------------------------------------------------------------------------
412 
413 double
414 GeomAsym::GetPzMeanRange(double alpha, double beta, double pz0, double pz1)
415 {
416  const double B = gamma_cdf(1 - pz0, alpha, beta) - gamma_cdf(1 - pz1, alpha, beta);
417  const double A = gamma_cdf(1 - pz0, alpha + 1, beta) - gamma_cdf(1 - pz1, alpha + 1, beta);
418  const double pz = 1 - alpha * beta * A / B;
419  return pz;
420 }
421 
422 
423 double
424 GeomAsym::GetShadow(double DX, double r, double theta, double psi, int icomp) // DX (g/cm^2) r (cm) theta/psi (rad)
425 {
426  if (fDetectorType == 2 && icomp > 0)
427  return 1; // TmodAmod not calculated
428 
429  // <pz>, Variance pz
430  const double pz_mean = GetPzMean(DX, r, icomp);
431  const double pz_variance = pow(GetPzRMS(DX, r, icomp), 2.);
432 
433  // Parameters of Gamma distribution
434  // Gamma dist: x=(0,inf) mean=alpha*beta variance=alpha*beta^2 | but mean/RMS estimated x=(0,1) differences are negligible
435  const double beta = pz_variance / (1 - pz_mean);
436  const double alpha = (1 - pz_mean) * (1 - pz_mean) / pz_variance;
437  const double I0 = gamma_cdf(2, alpha, beta);
438 
439  //The cut in the p_z distribution
440  double pz_cut = GetPzCut(r, theta, psi, icomp);
441 
442  //The shadowing factor
443  const double fshadow = gamma_cdf(1 - pz_cut, alpha, beta) / I0;
444 
445  return fshadow;
446 }
447 
448 
449 //------------------------------------------------------------------------------------------------------------
450 // Multiplying factor to convert the signal in a spherical tank (S0) to the signal in an Auger tank (S)
451 // Corrections due to shadowing effects/Tmod/Amod
452 //------------------------------------------------------------------------------------------------------------
453 double
454 GeomAsym::GetGeomAsym(double DX, double r, double theta, double psi, int icomp) // DX (g/cm^2) r (cm) theta/psi (rad)
455 {
456  int nSteps = 4;
457  return GetGeomAsym(DX, r, theta, psi, icomp, nSteps);
458 }
459 
460 
461 double
462 GeomAsym::GetGeomAsym(double DX, double r, double theta, double psi, int icomp, int nSteps) // DX (g/cm^2) r (cm) theta/psi (rad)
463 {
464  if (fDetectorType == 2 && icomp > 0)
465  return 1; // TmodAmod not calculated
466 
467  // const double pzMin = -1., pzMax = 1.;
468 
469  // <pz>, Variance pz
470  const double pz_mean = GetPzMean(DX, r, icomp);
471  const double pz_variance = pow(GetPzRMS(DX, r, icomp), 2.);
472 
473  // Parameters of Gamma distribution
474  // Gamma dist: x=(0,inf) mean=alpha*beta variance=alpha*beta^2 | but mean/RMS estimated in x=(0,1) differences are negligible
475  const double beta = pz_variance / (1 - pz_mean);
476  const double alpha = (1 - pz_mean) * (1 - pz_mean) / pz_variance;
477  const double I0 = gamma_cdf(2, alpha, beta);
478 
479  //The cut in the p_z distribution
480  const double pz_cut = GetPzCut(r, theta, psi, icomp);
481  //
482  double cdf[2], pz[2];
483  pz[0] = pz_cut;
484  cdf[0] = gamma_cdf(1 - pz[0], alpha, beta) / I0;
485 
486  double cdfStep = cdf[0] / double(nSteps);
487  if (pz[0] < pzIntegralLimit) {
488  while (true) {
489  cdf[1] = cdf[0] - cdfStep;
490  pz[1] = 1 - gamma_quantile(cdf[1] * I0, alpha, beta);
491  if (pz[1] < pzIntegralLimit)
492  break;
493  nSteps += 1;
494  cdfStep = cdf[0] / double(nSteps);
495  }
496  }
497 
498  //printf("nSteps=%d DX=%5.2lf cdfStep=%lf <pz>=%e RMS pz %e pz_cut %e\n",nSteps,DX,cdfStep,pz_mean,sqrt(pz_variance),pz_cut);
499 
500  double val = 0;
501  for (int istep = 0; istep < nSteps; ++istep) {
502  bool stop = false;
503  if (pz[0] > pzIntegralLimit) {
504  pz[1] = 1;
505  cdf[1] = 0;
506  stop = true;
507  } else {
508  cdf[1] = cdf[0] - cdfStep;
509  pz[1] = 1 - gamma_quantile(cdf[1] * I0, alpha, beta);
510  if (pz[1] > pzIntegralLimit) {
511  pz[1] = 1; //Fixes bug in gamma_quantile
512  cdf[1] = 0;
513  stop = true;
514  }
515  }
516  const double pzBin = GetPzMeanRange(alpha, beta, pz[0], pz[1]);
517  const double ctheta_p = GetCosTheta_p(pzBin, r, theta, psi, icomp);
518  const double AmodTmod = GetTmodAmod(ctheta_p, pzBin , r, icomp);
519  val += ((cdf[0] - cdf[1]) * AmodTmod);
520 
521  //printf("pz=(%e,%e) (%5.4lf -> %5.4lf) cdf=(%e,%e) %e AmodTmod=%e ctheta_p=%e | Bin integral=%e \n",pz[0],pz[1],
522  // pzBin,pz[1]-pz[0],cdf[0],cdf[1],cdf[0]-cdf[1],AmodTmod,ctheta_p,(cdf[0]-cdf[1])*AmodTmod );
523 
524  if (stop)
525  break;
526  cdf[0] = cdf[1];
527  pz[0] = pz[1];
528  }
529  //printf("Integral %lf \n\n\n",val);
530 
531  return val;
532 }
GeomAsymNS::parR2_MD
const double parR2_MD[4][2][3]
Definition: GeomAsym.h:156
GeomAsymNS::GeomAsym::GetCosTheta_p
double GetCosTheta_p(double pz, double r, double theta, double psi, int icomp)
Definition: GeomAsym.cc:365
GeomAsymNS::parPzRMS_MD
const double parPzRMS_MD[2][4][3]
Definition: GeomAsym.h:91
GeomAsymNS::ARadPar_MD
const double ARadPar_MD[4][3]
Definition: GeomAsym.h:58
GeomAsymNS::GeomAsym::GetGeomAsym
double GetGeomAsym(double DX, double r, double theta, double psi, int icomp)
Definition: GeomAsym.cc:454
GeomAsymNS::PzEval
const double PzEval[4][2]
Definition: GeomAsym.h:39
GeomAsymNS::GeomAsym::GetTmodAmod_Par_i
double GetTmodAmod_Par_i(double r, int itheta, int icomp, int ipar)
Definition: GeomAsym.cc:191
GeomAsymNS::GeomAsym::GeomAsym
GeomAsym(int DetectorType)
Definition: GeomAsym.cc:28
GeomAsymNS::ARadPar_Scin
const double ARadPar_Scin[4][3]
Definition: GeomAsym.h:51
std::pow
double pow(const double x, const unsigned int i)
Definition: MuonArrivalTime.h:17
exit
int exit
Definition: dump1090.h:237
GeomAsymNS::parPzMean_WCD
const double parPzMean_WCD[2][4][3]
Definition: GeomAsym.h:69
GeomAsymNS::GeomAsym::GetPzCut
double GetPzCut(double r, double theta, double psi, int icomp)
Definition: GeomAsym.cc:380
GeomAsymNS::GeomAsym::GetTmodAmod
double GetTmodAmod(double ctheta_p, double *par, int icomp)
Definition: GeomAsym.cc:257
GeomAsymNS::GeomAsym::GetPzMean
double GetPzMean(double DX, double r, int icomp)
Definition: GeomAsym.cc:96
GeomAsymNS::GeomAsym::gamma_cdf
double gamma_cdf(double x, double alpha, double beta)
Definition: GeomAsym.cc:45
GeomAsymNS::parR1_WCD
const double parR1_WCD[4][2][3]
Definition: GeomAsym.h:106
GeomAsymNS::parR0_WCD
const double parR0_WCD[4][2][3]
Definition: GeomAsym.h:99
GeomAsymNS::GeomAsym::GetPzRMS
double GetPzRMS(double DX, double r, int icomp)
Definition: GeomAsym.cc:159
GeomAsymNS::GeomAsym::pzIntegralLimit
double pzIntegralLimit
Definition: GeomAsym.h:171
GeomAsymNS::parR1_Scin
const double parR1_Scin[4][2][3]
Definition: GeomAsym.h:127
GeomAsymNS::parR0_MD
const double parR0_MD[4][2][3]
Definition: GeomAsym.h:142
GeomAsymNS::GeomAsym::GetARad
double GetARad(double r, int icomp)
Definition: GeomAsym.cc:333
GeomAsymNS::ctheta_pCut
const double ctheta_pCut[3]
Definition: GeomAsym.h:40
GeomAsymNS::parPzMean_MD
const double parPzMean_MD[2][4][3]
Definition: GeomAsym.h:87
GeomAsymNS::parR2_Scin
const double parR2_Scin[4][2][3]
Definition: GeomAsym.h:134
GeomAsymNS::parR1_MD
const double parR1_MD[4][2][3]
Definition: GeomAsym.h:149
GeomAsymNS::parPzRMS_Scin
const double parPzRMS_Scin[2][4][3]
Definition: GeomAsym.h:82
GeomAsymNS::DX_Eval
const double DX_Eval[4][2]
Definition: GeomAsym.h:38
GeomAsymNS::parPzMean_Scin
const double parPzMean_Scin[2][4][3]
Definition: GeomAsym.h:78
GeomAsymNS::parPzRMS_WCD
const double parPzRMS_WCD[2][4][3]
Definition: GeomAsym.h:73
GeomAsymNS::GeomAsym::GetPzMeanRange
double GetPzMeanRange(double alpha, double beta, double pz0, double pz1)
Definition: GeomAsym.cc:414
GeomAsymNS::GeomAsym::gamma_quantile
double gamma_quantile(double x, double alpha, double beta)
Definition: GeomAsym.cc:52
GeomAsymNS::nPar
const int nPar
Definition: GeomAsym.h:37
GeomAsymNS::GeomAsym::GetShadow
double GetShadow(double DX, double r, double theta, double psi, int icomp)
Definition: GeomAsym.cc:424
GeomAsymNS::GeomAsym::GetPzRMSPar
double GetPzRMSPar(double r, int icomp, int ipar)
Definition: GeomAsym.cc:130
GeomAsymNS::GeomAsym::GetPzMeanPar
double GetPzMeanPar(double r, int icomp, int ipar)
Definition: GeomAsym.cc:63
GeomAsymNS::ARadPar_WCD
const double ARadPar_WCD[4][3]
Definition: GeomAsym.h:44
GeomAsymNS::parR0_Scin
const double parR0_Scin[4][2][3]
Definition: GeomAsym.h:120
GeomAsymNS::parR2_WCD
const double parR2_WCD[4][2][3]
Definition: GeomAsym.h:113
, generated on Tue Sep 26 2023.