EMComponentIVR.cc
Go to the documentation of this file.
1 
9 #include "EMComponentIVR.h"
10 
11 #include <utl/ErrorLogger.h>
12 #include <utl/AugerUnits.h>
13 
14 #include <sstream>
15 #include <cmath>
16 #include <TMath.h>
17 
18 using namespace tls;
19 using namespace utl;
20 using namespace std;
21 
22 EMComponentIVR::EMComponentIVR(utl::Branch branch)
23 {
24  branch.GetChild("asymmetry").GetData(fAsymmetry);
25 
26  ostringstream msg;
27  msg << "\n Option:\n"
28  " Azimuthal asymmetry: " << (fAsymmetry ? "On" : "Off");
29  INFO(msg);
30 }
31 
32 EMComponentIVR::~EMComponentIVR()
33 {}
34 
35 double
36 EMComponentIVR::SignalRatio(double x, double y, double /* rmu */, double theta, double phi)
37 {
38 
39  // I.Valino parameterization. 30/05/2008
40 
41 
42 
43  // Provisional correction for theta < 60deg. // valid zenith range: th = [60,90] deg
44  // This provides an underestimate value.
45 
46  if (theta < GetThetaMin()) theta = GetThetaMin();
47  else if (theta > GetThetaMax()) theta = GetThetaMax();
48 
49 
50 
51  // x,y in Auger coord system and on ground plane in meters
52  // theta and phi in rad
53 
54  double xp, yp, zp;
55  double xShower, yShower;
56 
57  double th = theta/degree; // in deg
58 
59  // Shower coord. system: rotation around the z-axis
60  //setting x axis along the arrival direction
61 
62  xp = cos(phi)*x + sin(phi)*y;
63  yp =-sin(phi)*x + cos(phi)*y;
64  zp = 0.;
65 
66  // Shower plane: rotate around the xp axis to put the zp axis along the arrival direction.
67  xShower = cos(theta)*xp - sin(theta)*zp;
68  yShower = yp;
69 
70 
71  // Internal azimuth angle on the shower plane :
72 
73  double azi = atan2(yShower,xShower);
74  double azideg = azi*TMath::RadToDeg(); // in deg
75 
76 
77  double r = sqrt(xShower*xShower + yShower*yShower); // distance from the core on shower plane and in meters
78 
79  double rlog = log10(r);
80  const double rlogup=3.7;
81  double meanratio = 0;
82 
83  if (rlog<2)
84  meanratio=ffitnearEm(rlog,theta)/ffitnearMu(rlog,theta);
85  else if (2<=rlog && rlog<rlogup)
86  meanratio=ffitEm(rlog,theta)/ffitMu(rlog,theta);
87  else
88  meanratio=ffitEm(rlogup,theta)/ffitMu(rlogup,theta);
89 
90  double asym = 0;
91  double thup=68.; // Asymmetry in zenith range: th = [60,68] deg
92  if (fAsymmetry)
93  {
94  if (th<=thup)
95  {
96  if (rlog<=rlogup)
97  asym = ratioasym(rlog,th,azideg); // th and ph in deg
98  else
99  asym = ratioasym(rlogup,th,azideg);
100  }
101  }
102 
103  double ratio = meanratio*(1+asym);
104 
105  return ratio;
106 }
107 
108 double
109 EMComponentIVR::ffitMu(double r,double theta)
110  const
111 {
112 
113  double costh=cos(theta);
114  // here r is log10(r)
115 
116  double A1 = 2.398568e+00;
117  double A2 = 6.620959e+01;
118  double A3 = 1.371171e+02;
119  double A4 = -1.482223e+02;
120  double A = A1 + (A2*costh) + (A3*pow(costh,2.)) + (A4*pow(costh,3.));
121 
122  double B = 2000.;
123 
124  double C1 = 6.364131e-01;
125  double C2 = 5.880061e-01;
126  double C3 = 4.025384e-03;
127  double C = C1 + (C2*costh) + (C3*pow(costh,2.));
128 
129  double D1 = 4.299209e+00;
130  double D2 = 5.668097e+00;
131  double D3 = -1.345011e+01;
132  double D4 = 1.195753e+01;
133  double D = D1 + (D2*costh) + (D3*pow(costh,2.)) + (D4*pow(costh,3.));
134 
135 
136 
137  double E1 = -6.214005e-03;
138  double E2 = 1.735574e-02;
139  double E3 = 3.047811e-02;
140  double E4 = -3.364831e-01;
141  double E5 = 3.505095e-01;
142  double E = E1 + (E2*costh) + (E3*pow(costh,2.))+ (E4*pow(costh,3.))+ (E5*pow(costh,4.));
143 
144 
145  double xx=pow(10,r);
146 
147  double t1 = xx/B;
148  double t2= 1+(xx/B);
149 
150  double fmu= A *pow(t1,-C)*pow(t2,-D) + E;
151  return fmu;
152 
153 }
154 
155 double EMComponentIVR::ffitEm(double r,double theta)
156  const
157 {
158 
159  double costh=cos(theta);
160  // here r is log10(r)
161 
162  double thetathr1 = 70.*TMath::DegToRad();
163  double thetathr2 = 72.*TMath::DegToRad();
164 
165  double A1,A2,A3,A4,A;
166  double C1,C2,C3,C;
167  double D1,D2,D3,D4,D;
168  double E1,E2,E;
169  double B=2000.;
170 
171 
172  if (theta<thetathr1)
173  {
174 
175  D1 = 9.689558e+00;
176  D2 = -1.361415e+01;
177  D = D1 + (D2*costh);
178 
179  E1 = 8.498621e-03;
180  E2 = -2.592573e-02;
181  }
182  else
183  {
184  D1 = 4.028661e+00 ;
185  D2 = 3.251934e+00;
186  D3 = 1.407006e+01;
187  D4 = -4.465588e+01;
188  D = D1 + (D2*costh)+ (D3*pow(costh,2.)) +(D4*pow(costh,3.));
189 
190  E1 = -7.711403e-04 ;
191  E2 = 6.489431e-04;
192  }
193 
194  E = E1 + (E2*costh);
195 
196  if (theta<thetathr2)
197  {
198 
199  A1 = 1.38090e-01;
200  A2 = 8.50140e+01;
201  A3 = -2.86508e+02;
202  A4 = 2.52676e+02;
203 
204  A = A1 + (A2*costh)+ (A3*pow(costh,2.)) + (A4*pow(costh,3.));
205 
206 
207  C1 = 2.74957e+00 ;
208  C2 = -1.30478e+01;
209  C3 = 2.33996e+01;
210 
211 
212  C = C1 + (C2*costh)+ (C3*pow(costh,2.));
213  }
214  else
215  {
216 
217  A1 = 3.836692e-01;
218  A2 = 9.103229e+00;
219  A3 = 2.780868e+01;
220  A = A1 + (A2*costh)+ (A3*pow(costh,2.));
221 
222  C1 = 1.023005e+00;
223  C2 = -1.936833e-01;
224  C = C1 + (C2*costh);
225 
226  }
227 
228  double xx=pow(10,r);
229  double t1 = xx/B;
230  double t2= 1 +(xx/B);
231 
232  double fem=A *pow(t1,-C)*pow(t2,-D) + E;
233  return fem;
234 }
235 
236 
237 double
238 EMComponentIVR::ffitnearMu(double r,double theta)
239  const
240 {
241  double costh=cos(theta);
242  // here r is log10(r)
243 
244 
245  double A1 = 2.658940e+02;
246  double A2 = 7.586829e+03;
247  double A3 = 3.467820e+04;
248  double A4 = 5.973019e+03;
249  double A5 = 5.735970e+04;
250  double A = A1 + (A2*costh) + (A3*pow(costh,2.)) + (A4*pow(costh,3.))+ (A5*pow(costh,4.));
251 
252  double B = 0.6;
253 
254  double ld=2.;
255  double d=pow(10.,ld);
256  double fgeneral=ffitMu(ld,theta);
257  double fact=pow(d,-B);
258  double denom=fgeneral/(A*fact);
259  double C = -d/log(denom);
260 
261 
262  double xx=pow(10,r);
263  double t1 = xx;
264  double t2 = -xx/C;
265 
266  double fmu2 = A*pow(t1,-B)*exp(t2);
267  return fmu2;
268 
269 }
270 
271 double
272 EMComponentIVR::ffitnearEm(double r,double theta)
273  const
274 {
275  // here r is log10(r)
276 
277  double costh=cos(theta);
278  double thetathr = 72.*TMath::DegToRad();
279 
280  double C = 300.;
281  double B1,B2,B3,B4,B;
282 
283 
284  if (theta < thetathr) // theta in rad
285  {
286  B1 = 1.902538e+01;
287  B2 = -1.455794e+02;
288  B3 = 3.724824e+02;
289  B4 = -3.001939e+02;
290  B = B1 + (B2*costh) + (B3*pow(costh,2.)) + (B4*pow(costh,3.));
291  }
292  else
293  {
294 
295  B1 = 1.072279e+00;
296  B2 = -1.003438e+00;
297  B = B1 + (B2*costh);
298  }
299 
300  double ld=2.;
301  double d=pow(10.,ld);
302  double fgeneral=ffitEm(ld,theta);
303  double fact1=pow(d,-B);
304  double fact2=exp(-d/C);
305 
306  double A = fgeneral/(fact1*fact2);
307 
308  double xx=pow(10,r);
309  double t1 = xx;
310  double t2 = -xx/C;
311 
312  double fem2 = A*pow(t1,-B)*exp(t2);
313  return fem2;
314 }
315 
316 
317 
318 double
319 EMComponentIVR::ratioasym(double r,double theta, double phi)
320  const
321 {
322 
323  // here r is log10(r), theta & phi in deg
324 
325  double xx=pow(10.,r);
326  double xxc = 1000.;
327 
328  double phirad = phi*TMath::DegToRad(); // in rad
329 
330  if (phi>180.) phi = phi-360.;
331 
332 
333  double p0a = 2053.5;
334  double p0b = -30.1068;
335 
336  double p0 = p0a + (p0b*theta);
337 
338 
339  double p1a = 80.258;
340  double p1b = 21.3549;
341  double p1c = 11.1303;
342  double var = (90.-theta-p1b)/sqrt(p1c);
343  double p1 =p1a*(1-TMath::Erfc(var));
344 
345 
346  double p2a =0.0526881;
347  double p2b =-0.00175893;
348  double p2c =1.52561e-05;
349  double p2 = p2a + (p2b*theta)+ (p2c*theta*theta);
350 
351 
352 
353  double p3 = 0.492845;
354  double p4 = -0.281903;
355 
356 
357  double Norm = p0*((1./TMath::Pi())*(p1/(pow(phi,2)+pow(p1,2)))-p2);
358 
359  double B = p3 *cos(p4*phirad);
360 
361  double A = Norm /pow(xxc,B);
362 
363 
364  double asy = A*pow(xx,B);
365 
366  return asy;
367 }
368 
369 
370 // Configure (x)emacs for this file ...
371 // Local Variables:
372 // mode:c++
373 // compile-command: "make -C .. -k"
374 // End:
degree
const double degree
Definition: PhysicsConstants.h:3
tls::EMComponentIVR::SignalRatio
virtual double SignalRatio(double x, double y, double rmu, double theta, double phi)
Definition: EMComponentIVR.cc:36
tls::EMComponentIVR::ffitnearMu
double ffitnearMu(double r, double theta) const
Definition: EMComponentIVR.cc:238
INFO
#define INFO(message)
Macro for logging informational messages.
Definition: ErrorLogger.h:161
tls::EMComponentIVR::ffitMu
double ffitMu(double r, double theta) const
Definition: EMComponentIVR.cc:109
utl::Branch::GetChild
Branch GetChild(const std::string &childName) const
Get child of this Branch by child name.
Definition: Branch.cc:211
std::pow
double pow(const double x, const unsigned int i)
Definition: MuonArrivalTime.h:17
utl::Branch
Class representing a document branch.
Definition: Branch.h:107
tls::EMComponentIVR::ffitEm
double ffitEm(double r, double theta) const
Definition: EMComponentIVR.cc:155
utl::Branch::GetData
void GetData(bool &b) const
Overloads of the GetData member template function.
Definition: Branch.cc:644
t2
a second level trigger
Definition: XbT2.h:8
tls::EMComponentIVR::ffitnearEm
double ffitnearEm(double r, double theta) const
Definition: EMComponentIVR.cc:272
tls::EMComponentIVR::ratioasym
double ratioasym(double r, double theta, double phi) const
Definition: EMComponentIVR.cc:319
tls::EMComponentIVR::EMComponentIVR
EMComponentIVR(utl::Branch branch)
Definition: EMComponentIVR.cc:22
, generated on Tue Sep 26 2023.