source: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc@ 1510

Last change on this file since 1510 was 1449, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 5.6 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
19! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de)
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27// //
28// //
29//////////////////////////////////////////////////////////////////////////////
30#include "MPhoton.h"
31
32#include <iostream.h>
33
34#include <TF1.h>
35#include <TH1.h>
36#include <TPad.h>
37#include <TCanvas.h>
38
39ClassImp(MPhoton);
40
41Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k)
42{
43 const Double_t m2 = x[0]; // m2: (E0/sqrt(s))^2
44
45 const Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2
46
47 const Double_t beta2 = 1.-m2;
48 const Double_t beta = sqrt(beta2);
49
50 const Double_t p1 = r0*r0*TMath::Pi()/2;
51
52 // ----- Extreme Relativistic -------
53 // return p1*2 * m*m*m* (log(2./m)-1);
54 // ----------------------------------
55
56 const Double_t p2 = 3.-beta2*beta2;
57 const Double_t p3 = log((1.+beta)/(1.-beta));
58 const Double_t p4 = beta*2*(1.+m2);
59
60 const Double_t sigma = p1*m2*(p2*p3-p4); // [m^2]
61
62 return sigma;
63}
64
65Double_t MPhoton::Int1(Double_t *x, Double_t *k)
66{
67 const Double_t costheta = x[0];
68
69 const Double_t Eg = k[0];
70 const Double_t Ep = k[1];
71
72 const Double_t E0 = 511e-6; // [GeV]
73
74 Double_t s = E0/Eg*E0/Ep/(1.-costheta)/2;
75 if (s>1) // Why is this necessary???
76 return 0;
77
78 const Double_t sigma = Sigma_gg(&s); // [m^2]
79
80 return sigma/2 * (1.-costheta); // [m^2]
81}
82
83Double_t MPhoton::Int2(Double_t *x, Double_t *k)
84{
85 const Double_t E0 = 511e-6; // [GeV]
86
87 Double_t Ep = x[0];
88 Double_t z = k[1];
89
90 const Double_t Eg = k[0];
91
92 Double_t val[2] = { Eg, Ep };
93
94 static TF1 f("int1", Int1, 0, 0, 2);
95
96 const Double_t from = -1.0;
97 const Double_t to = 1.-E0*E0/(2.*Eg*Ep); // Originally Was: 1.
98 const Double_t int1 = f.Integral(from, to, val, 1e-2); // [m^2]
99 const Double_t planck = MParticle::Planck(&Ep, &z); // [GeV^2]
100
101 const Double_t res = planck * int1;
102
103 return res; // [GeV^2 m^2]
104}
105
106// --------------------------------------------------------------------------
107//
108// Returns 0 in case IL becomes (numerically) infinite.
109//
110Double_t MPhoton::InteractionLength(Double_t *x, Double_t *k)
111{
112 Double_t E0 = 511e-6; // [GeV]
113 Double_t c = 299792458; // [m/s]
114 Double_t e = 1.602176462e-19; // [C]
115 Double_t h = 1e-9/e*6.62606876e-34; // [GeVs]
116 Double_t hc = h*c; // [GeVm]
117 Double_t pc = 1./3.258; // [pc/ly]
118 Double_t ly = 3600.*24.*365.*c; // [m/ly]
119
120 Double_t Eg = x[0];
121 Double_t z = k ? k[0] : 0;
122
123 if (Eg<100)
124 return 1e50;
125
126 Double_t val[2] = { Eg, z };
127
128 static TF1 f("int2", Int2, 0, 0, 2);
129
130 Double_t lolim = E0*E0/Eg;
131 Double_t inf = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
132
133 if (Eg<5e4)
134 //inf = 3e-11*(z+1)*pow(10, 4.7*0.5-log10(Eg)*0.5);
135 inf = 3e-11*(z+1)*pow(10, 4.7-log10(Eg));
136 //inf = 3e-11*(z+1)*pow(10, 7.0-log10(Eg)*1.5);
137 //inf = 3e-11*(z+1)*pow(10, 8.2-log10(Eg)*1.75);
138 //inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
139
140 Double_t int2 = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2]
141
142 if (int2==0)
143 {
144 //cout << "---> Int2==0 <---" << endl;
145 return 0;
146 }
147
148 /* Planck constants: konst */
149 const Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
150 int2 *= konst; // [1 / m]
151
152 Double_t res = 1./ int2; // [m]
153 res *= pc/ly * 1e-3; // [kpc]
154
155 if (res > 1e50) return 1e50;
156 if (res < 0) return 1e35;
157
158 return res; //[kpc]
159}
160
161Double_t MPhoton::GetInteractionLength(Double_t energy, Double_t z)
162{
163 return InteractionLength(&energy, &z);
164}
165
166Double_t MPhoton::GetInteractionLength() const
167{
168 return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ);
169}
170
171void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to, Option_t *opt)
172{
173 if (!gPad)
174 new TCanvas("ILPhoton", "Mean Interaction Length Photon");
175 else
176 gPad->GetVirtCanvas()->cd(4);
177
178 TF1 f1("length", InteractionLength, from, to, 1);
179 f1.SetParameter(0, z);
180
181 gPad->SetLogx();
182 gPad->SetLogy();
183 gPad->SetGrid();
184 f1.SetMinimum(1);
185 f1.SetMaximum(1e9);
186 f1.SetLineWidth(1);
187
188 TH1 &h=*f1.DrawCopy(opt)->GetHistogram();
189
190 h.SetTitle("Mean Interaction Length (Photon)");
191 h.SetXTitle("E [GeV]");
192 h.SetYTitle("x [kpc]");
193
194 gPad->Modified();
195 gPad->Update();
196}
197
198void MPhoton::DrawInteractionLength() const
199{
200 DrawInteractionLength(fZ);
201}
202
203void MPhoton::Fill(TH1 &h, Double_t idx, Double_t w) const
204{
205 h.Fill(fEnergy, pow(fEnergy, idx)*w);
206}
Note: See TracBrowser for help on using the repository browser.