source: trunk/Mars/hawc/atmprof.C@ 19932

Last change on this file since 19932 was 19770, checked in by tbretz, 6 years ago
Fit an atmospheric profile from an atmprof file to get the coefficients used in MAtmosphere to calculate the extinction coefficients.
File size: 6.4 KB
Line 
1Double_t Log10(Double_t lg)
2{
3 // Avoid numerical limits
4 return lg<=1e-6 ? -6 : log10(lg);
5}
6
7// The mass overburden in each layer is given by
8// T = AATM + BATM * exp(-h/CATM)
9// Note that the AATM are correlated with the layer border
10Double_t fcn(Double_t *xx, Double_t *par)
11{
12 Double_t *A = par;
13 Double_t *B = par+5;
14 Double_t *C = par+10;
15 Double_t *H = par+15;
16
17 Double_t x = xx[0]*1e5;
18
19 double Ai = A[0];
20
21 for (int i=0; i<4; i++)
22 {
23 if (x>=H[i] && x<H[i+1])
24 return Log10(Ai + B[i]*exp(-x/C[i]));
25
26 Ai += B[i] *exp(-H[i+1]/C[i]); // y-value of current bin function [i] at transition border H[i+1] to next layer
27 Ai -= B[i+1]*exp(-H[i+1]/C[i+1]); // y-value of next bin function [i+1] at transition border H[i+1] to next layer
28 }
29
30 return Log10(Ai + B[4]*exp(-x/C[4]));
31}
32
33void atmprof()
34{
35 // As in CORSIKA the atmosphere for atmabs.dat is required up to 50 km
36 // the best is to do a fit which is bext up to at least 50 km
37 double MAX = 50.5;
38
39 // MAGIC Winter (MAGIC Summer: atmprof9)
40 TGraph atmprof("/home/tbretz/SW/corsika-76900/bernlohr/atmprof8.dat", "%lg %*lg %lg, %*lg");
41
42 // The last point has a thickness of 0 (log!)
43 atmprof.RemovePoint(atmprof.GetN()-1);
44
45 // Calculate the logarithm (better fit)
46 double miny = log10(atmprof.GetY()[0]);
47 for (int i=0; i<atmprof.GetN(); i++)
48 {
49 atmprof.GetY()[i] = log10(atmprof.GetY()[i]);
50 if (atmprof.GetX()[i]<MAX)
51 miny = TMath::Min(miny, atmprof.GetY()[i]);
52 }
53
54 // =============================================================================
55
56 TCanvas *c = new TCanvas;
57 c->Divide(1, 2);
58 c->cd(1);
59
60 atmprof.SetMarkerStyle(kFullDotMedium);
61 atmprof.GetHistogram()->GetXaxis()->SetRangeUser(0, MAX);
62 atmprof.GetHistogram()->GetYaxis()->SetRangeUser(miny, atmprof.GetY()[0]*1.05);
63 atmprof.DrawClone("ALP");
64
65 TMarker m;
66 m.SetMarkerStyle(kFullDotLarge);
67
68 // =============================================================================
69
70 // Define the fit function betwene 0 and 50km with 15 free parameters
71 TF1 f("thickness", fcn, 0, MAX, 4*5);
72 f.SetNpx(10000);
73
74 // This is the Keilhauer U.S. Standard Atmosphere (as starting values)
75 Double_t params[] =
76 {
77 -149.801663, -57.932486, 0.63631894, 4.35453690e-4, 0.01128292, // A
78 1183.6071, 1143.0425, 1322.9748, 655.67307, 1, // B
79 954248.34, 800005.34, 629568.93, 737521.77, 1e9, // C
80 0, 7.0e5, 11.4e5, 37.0e5, 100e5, // H
81
82 };
83 f.SetParameters(params);
84
85 // =============================================================================
86 // Calc residuals
87
88 TGraph g1;
89 for (int i=0; i<atmprof.GetN(); i++)
90 {
91 double x = atmprof.GetX()[i];
92 if (x>MAX)
93 break;
94
95 double y = atmprof.GetY()[i];
96 g1.SetPoint(g1.GetN(), x, 1-pow(10, f.Eval(x))/pow(10, y));
97 }
98
99
100 // =============================================================================
101 // Plot and print
102
103 m.SetMarkerColor(kRed);
104 f.SetLineColor(kRed);
105
106 f.DrawCopy("same");
107
108 cout << endl;
109 for (int i=0; i<5; i++)
110 {
111 double A[5] = { f.GetParameters()[0], 0, 0, 0, 0 };
112 double *B = f.GetParameters()+5;
113 double *C = f.GetParameters()+10;
114 double *H = f.GetParameters()+15;
115
116 A[1] = A[0] + B[0]*exp(-H[1]/C[0]) - B[1]*exp(-H[1]/C[1]);
117 A[2] = A[1] + B[1]*exp(-H[2]/C[1]) - B[2]*exp(-H[2]/C[2]);
118 A[3] = A[2] + B[2]*exp(-H[3]/C[2]) - B[3]*exp(-H[3]/C[3]);
119 A[4] = A[3] + B[3]*exp(-H[4]/C[3]) - B[4]*exp(-H[4]/C[4]);
120
121 cout << Form("A[%d]=%13f B[%d]=%13f C[%d]=%13f", i, A[i], i, B[i], i, C[i]);
122 if (i<4)
123 {
124 cout << Form(" H=%13f", H[i+1]/1e5);
125 m.DrawMarker(H[i+1]/1e5, f.Eval(H[i+1]/1e5));
126 }
127 cout << endl;
128 }
129
130 // =============================================================================
131
132 // The parameters above 50km should not be left free
133 f.FixParameter( 1, params[1]);
134 f.FixParameter( 2, params[2]);
135 f.FixParameter( 3, params[3]);
136 f.FixParameter( 4, params[4]);
137 f.FixParameter( 9, params[9]);
138 f.FixParameter(14, params[14]);
139 f.FixParameter(15, 0);
140
141 // Fit
142 float max = 5;
143 while (max<MAX)
144 {
145 atmprof.Fit(&f, "", "", 0, max);
146 max += 5;
147 }
148
149 // =============================================================================
150 // Plot and print
151
152 m.SetMarkerColor(kBlue);
153 f.SetLineColor(kBlue);
154
155 f.DrawCopy("same");
156
157 cout << endl;
158 for (int i=0; i<5; i++)
159 {
160 double A[5] = { f.GetParameters()[0], 0, 0, 0, 0 };
161 double *B = f.GetParameters()+5;
162 double *C = f.GetParameters()+10;
163 double *H = f.GetParameters()+15;
164
165 A[1] = A[0] + B[0]*exp(-H[1]/C[0]) - B[1]*exp(-H[1]/C[1]);
166 A[2] = A[1] + B[1]*exp(-H[2]/C[1]) - B[2]*exp(-H[2]/C[2]);
167 A[3] = A[2] + B[2]*exp(-H[3]/C[2]) - B[3]*exp(-H[3]/C[3]);
168 A[4] = A[3] + B[3]*exp(-H[4]/C[3]) - B[4]*exp(-H[4]/C[4]);
169
170 cout << Form("A[%d]=%13f B[%d]=%13f C[%d]=%13f", i, A[i], i, B[i], i, C[i]);
171 if (i<4)
172 {
173 cout << Form(" H=%13f", H[i+1]/1e5);
174 m.DrawMarker(H[i+1]/1e5, f.Eval(H[i+1]/1e5));
175 }
176 cout << endl;
177 }
178
179 // =============================================================================
180 // Plot residuals
181
182 c->cd(2);
183
184 TGraph g2;
185 for (int i=0; i<atmprof.GetN(); i++)
186 {
187 double x = atmprof.GetX()[i];
188 if (x>MAX)
189 break;
190
191 double y = atmprof.GetY()[i];
192 g2.SetPoint(g2.GetN(), x, 1-pow(10, f.Eval(x))/pow(10, y));
193 }
194
195 g1.SetLineColor(kRed);
196 g1.SetMarkerColor(kRed);
197 g1.SetMarkerStyle(kFullDotMedium);
198
199 g2.SetLineColor(kBlue);
200 g2.SetMarkerColor(kBlue);
201 g2.SetMarkerStyle(kFullDotMedium);
202
203 TH1C h("", "Residuals", 100, 0, MAX);
204 Double_t xmin, ymin, xmax, ymax;
205 g1.ComputeRange(xmin, ymin, xmax, ymax);
206 Double_t maxy = TMath::Max(fabs(ymin)*1.05, fabs(ymax)*1.05);
207 h.SetMinimum(-maxy);
208 h.SetMaximum( maxy);
209 h.SetStats(kFALSE);
210 h.DrawCopy();
211
212 g1.GetHistogram()->GetXaxis()->SetRangeUser(0, MAX);
213 g1.DrawClone("LP");
214 g2.DrawClone("LP");
215
216 g2.ComputeRange(xmin, ymin, xmax, ymax);
217 maxy = TMath::Max(fabs(ymin)*1.05, fabs(ymax)*1.05);
218
219 cout << "\nMax residual after fit: " << maxy*100 << " %\n" << endl;
220}
Note: See TracBrowser for help on using the repository browser.