source: trunk/MagicSoft/Mars/mhist/MHFlux.cc@ 1575

Last change on this file since 1575 was 1573, checked in by wittek, 24 years ago
change type of fVarname and fUnit from *char to TString
File size: 20.0 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): Wolfgang Wittek 5/2002 <mailto:wittek@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26// //
27// MHFlux //
28// //
29// calculates absolute photon fluxes //
30// from the distributions of the estimated energy //
31// for the different bins in some variable 'Var' //
32// (Var = Theta or time) //
33// //
34//////////////////////////////////////////////////////////////////////////////
35
36#include "MHFlux.h"
37
38#include <TStyle.h>
39
40#include <TF1.h>
41#include <TH2.h>
42#include <TProfile.h>
43
44
45#include <TCanvas.h>
46
47#include "MTime.h"
48
49#include "MBinning.h"
50#include "MParList.h"
51
52#include "MLog.h"
53#include "MLogManip.h"
54
55ClassImp(MHFlux);
56
57
58// --------------------------------------------------------------------------
59//
60// Default Constructor. It sets the variable name (Theta or time)
61// and the units for the variable
62//
63MHFlux::MHFlux(const TH2D &h2d, const Bool_t Draw,
64 // const char *varname, const char *unit)
65 const TString varname, const TString unit)
66 : fHOrig(), fHUnfold(), fHFlux(), fVarname(), fUnit()
67{
68 // if (varname == NULL || unit == NULL)
69 if (varname == "" || unit == "")
70 {
71 *fLog << "MHFlux : varname or unit not defined" << endl;
72 }
73
74 fVarname = varname;
75 fUnit = unit;
76
77 TString strg(varname);
78 strg += unit;
79
80 // char txt[100];
81
82 // original distribution of E-est for different bins
83 // of the variable (Theta or time)
84 // sprintf(txt, "gammas vs. E-est and %s",varname);
85
86 TString strg1 = "no.of gammas vs. E-est and ";
87 strg1 += varname;
88
89 ((TH2D&)h2d).Copy(fHOrig);
90
91 fHOrig.SetName("E-est");
92 fHOrig.SetTitle(strg1);
93
94 fHOrig.SetDirectory(NULL);
95 fHOrig.SetXTitle("E-est [GeV] ");
96 fHOrig.SetYTitle(strg);
97 fHOrig.Sumw2();
98
99 SetBinning((TH2*)&fHOrig, (TH2*)&h2d);
100
101
102 // unfolded distribution of E-unfold for different bins
103 // of the variable (Theta or time)
104 // sprintf(txt, "gammas vs. E-unfold and %s",varname);
105 TString strg2 = "no.of gammas vs. E-unfold and ";
106 strg2 += varname;
107
108 fHUnfold.SetName("E-unfolded");
109 fHUnfold.SetTitle(strg2);
110
111 fHUnfold.SetDirectory(NULL);
112 fHUnfold.SetXTitle("E-unfold [GeV] ");
113 fHUnfold.SetYTitle(strg);
114 fHUnfold.Sumw2();
115
116 SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);
117
118
119 // absolute photon flux vs. E-unfold
120 // for different bins of the variable (Theta or time)
121 //
122 // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
123 TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and ";
124 strg3 += varname;
125
126 fHFlux.SetName("photon flux");
127 fHFlux.SetTitle(strg3);
128
129 fHFlux.SetDirectory(NULL);
130 fHFlux.SetXTitle("E-unfold [GeV] ");
131 fHFlux.SetYTitle(strg);
132 fHFlux.Sumw2();
133
134 SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);
135
136
137 // copy fHOrig into fHUnfold in case no unfolding is done
138 const Int_t nEunf = fHUnfold.GetNbinsX();
139 const Int_t nVar = fHUnfold.GetNbinsY();
140 for (int m=1; m<=nEunf; m++)
141 {
142 for (int n=1; n<=nVar; n++)
143 {
144 Double_t cont = fHOrig.GetBinContent(m,n);
145 Double_t dcont = fHOrig.GetBinError(m,n);
146 fHUnfold.SetBinContent(m,n,cont);
147 fHUnfold.SetBinError(m,n,dcont);
148 }
149 }
150 //..............................................
151 // draw the No.of photons vs. E-est
152 // for the individual bins of the variable Var
153
154 if (Draw == kTRUE)
155 {
156 const Int_t nVar = fHOrig.GetNbinsY();
157
158 for (int n=1; n<=nVar; n++)
159 {
160 TString strg0("Orig-");
161 strg0 += fVarname;
162 TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");
163
164 char txt0[100];
165 strg0 = fVarname;
166 strg0 += "-bin %d";
167 sprintf(txt0, strg0, n);
168
169 TString strg1("No.of photons vs. E-est for ");
170 strg1 += txt0;
171 new TCanvas(txt0,strg1);
172 // TCanvas &c = *MakeDefCanvas(txt0, strg1);
173 // gROOT->SetSelectedPad(NULL);
174
175 gPad->SetLogx();
176
177 h.SetName(txt0);
178 h.SetTitle(strg1);
179 h.SetXTitle("E-est [GeV] ");
180 h.SetYTitle("No.of photons");
181 h.DrawCopy();
182
183 // c.Modified();
184 // c.Update();
185 }
186 }
187 //........................
188}
189
190// -------------------------------------------------------------------------
191//
192// Dummy Fill (has to be included because in base class MH Fill is set to 0);
193// without the dummy Fill one gets the error message :
194//
195// Error: Can't call MHFlux::MHFlux(evttime,"time","[s]") in current scope
196// FILE:macros/flux.C LINE:465
197// Possible candidates are...
198// filename line:size busy function type and name (in MHFlux)
199// filename line:size busy function type and name (in MH)
200// filename line:size busy function type and name (in MParContainer)
201// filename line:size busy function type and name (in TObject)
202//
203Bool_t MHFlux::Fill(const MParContainer *par)
204{
205 return kTRUE;
206}
207
208
209// -------------------------------------------------------------------------
210//
211// Unfold the distribution in E-est
212//
213void MHFlux::Unfold(const Bool_t Draw)
214{
215 //..............................................
216 // draw the No.of photons vs. E-unfold
217 // for the individual bins of the variable Var
218
219 if (Draw == kTRUE)
220 {
221 const Int_t nVar = fHUnfold.GetNbinsY();
222
223 for (int n=1; n<=nVar; n++)
224 {
225 TString strg0("Unfold-");
226 strg0 += fVarname;
227 TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
228
229 char txt0[100];
230 strg0 = fVarname;
231 strg0 += "-bin %d";
232 sprintf(txt0, strg0, n);
233
234 TString strg1("No.of photons vs. E-unfold for ");
235 strg1 += txt0;
236 new TCanvas(txt0,strg1);
237
238 // TCanvas &c = *MakeDefCanvas(txt0, strg1);
239 // gROOT->SetSelectedPad(NULL);
240
241 gPad->SetLogx();
242
243 h.SetName(txt0);
244 h.SetTitle(strg1);
245 h.SetXTitle("E-unfold [GeV] ");
246 h.SetYTitle("No.of photons");
247 h.DrawCopy();
248
249 // c.Modified();
250 // c.Update();
251 }
252 }
253 //........................
254}
255
256
257// -------------------------------------------------------------------------
258//
259// Calculate photon flux by dividing the distribution in Eunf (fHUnfold) by
260// the width of the energy interval (deltaE)
261// the effective ontime (*teff)
262// and the effective collection area (*aeff)
263//
264void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar,
265 const TH2D *aeff, const Bool_t Draw)
266{
267 // Note that fHUnfold has bins in Eunf and Var
268 // *teff has bins in Var (the same bins in Var as fHUnfold)
269 // *thetabar has bins in Var (the same bins in Var as fHUnfold)
270 // *aeff has bins in Etru and Theta
271 // (where in general the binning in Etru is different
272 // from the binning in Eunf)
273 // The variable Var may be 'time' or 'Theta'
274
275 // Draw = kTRUE means the differential flux vs E-unf should be drawn
276 // for the individual bins of the variable Var
277
278 //....................................
279 // define dummy histogram *aeff
280 ((TH1*)aeff)->Sumw2();
281 MBinning binsetru("BinningEtru");
282 binsetru.SetEdgesLog(10, 10, 1e3);
283
284 MBinning binsthetatru("BinningThetatru");
285 binsthetatru.SetEdges(7, -2.5, 32.5);
286 //SetBinning((TH1*)aeff, &binsetru, &binsthetatru);
287 SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
288
289 const Int_t netru = aeff->GetNbinsX();
290 const Int_t ntheta = aeff->GetNbinsY();
291
292 for (int j=1; j<=netru; j++)
293 {
294 for (int k=1; k<=ntheta; k++)
295 {
296 Double_t cont = 10000.0;;
297 ((TH1*)aeff)->SetBinContent(j,k,cont);
298
299 Double_t dcont = 100.0;
300 ((TH1*)aeff)->SetBinError(j,k,dcont);
301 }
302 }
303 // *fLog << "Dummy aeff : netru =" << netru << ", ntheta = " << ntheta << endl;
304 //....................................
305
306 // number of Eunf and Var bins (histograms : fHUnfold, fHFlux)
307 const Int_t nEunf = fHFlux.GetNbinsX();
308 const Int_t nVar = fHFlux.GetNbinsY();
309
310 // number of Etru and Theta bins (histogram *aeff of collection area)
311
312 const Int_t nEtru = aeff->GetNbinsX();
313 const Int_t nTheta = aeff->GetNbinsY();
314
315 //*fLog << "nEunf =" << nEunf << ", nVar = " << nVar << endl;
316 //*fLog << "nEtru =" << nEtru << ", nTheta = " << nTheta << endl;
317
318 //...................................................................
319 // calculate effective collection area
320 // for the Eunf and Var bins of the histogram fHUnfold
321 // from the histogram *aeff, which has bins in Etru and Theta
322 // the result is the histogram fHAeff
323 //
324
325 TH2D fHAeff;
326 fHAeff.Sumw2();
327 SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
328
329 Int_t errflag;
330 Double_t c0, c1, c2;
331 Double_t t1, t2, t3;
332 Double_t a1, a2, a3;
333
334 Double_t *aeffbar;
335 aeffbar = new Double_t[nEtru];
336 Double_t *daeffbar;
337 daeffbar = new Double_t[nEtru];
338
339 Double_t aeffEunfVar;
340 Double_t daeffEunfVar;
341
342
343 //------ start n loop ------
344 for (int n=1; n<=nVar; n++)
345 {
346 Double_t Thetabar = thetabar->GetBinContent(n);
347 Double_t cosThetabar = cos(Thetabar);
348
349 // determine Theta bins (k1, k2, k3) for interpolation in Theta
350 // k0 denotes the Theta bin from whicvh the error is copied
351 Int_t k0=0, k1=0, k2=0, k3=0;
352 for (int k=3; k<=nTheta; k++)
353 {
354 Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(k);
355 if (Thetabar < Thetalow)
356 {
357 k1 = k-2;
358 k2 = k-1;
359 k3 = k;
360 k0 = k2;
361 break;
362 }
363 }
364
365 if (k3 == 0)
366 {
367 k1 = nTheta-2;
368 k2 = nTheta-1;
369 k3 = nTheta;
370 k0 = k2;
371 }
372
373 if (Thetabar < ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2))
374 k0 = 1;
375 else if (Thetabar > ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta))
376 k0 = nTheta;
377
378 Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1);
379 Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta+1);
380 if (Thetabar < Thetamin || Thetabar > Thetamax)
381 {
382 *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar
383 << ", Thetamin =" << Thetamin
384 << ", Thetamax =" << Thetamax << endl;
385 }
386
387 //*fLog << "Var bin " << n << ":" << endl;
388 //*fLog << "Thetabar= " << Thetabar
389 // << ", k0= " << k0
390 // << ", k1= " << k1
391 // << ", k2= " << k2
392 // << ", k3= " << k3 << endl;
393
394
395 // calculate effective collection area at Theta = Thetabar
396 // by quadratic interpolation in cos(Theta);
397 // do this for each bin of Etru
398 //
399
400 for (int j=1; j<=nEtru; j++)
401 {
402 c0 = 0.0;
403 c1 = 0.0;
404 c2 = 0.0;
405
406 t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) );
407 t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) );
408 t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) );
409
410 a1 = aeff->GetBinContent(j,k1);
411 a2 = aeff->GetBinContent(j,k2);
412 a3 = aeff->GetBinContent(j,k3);
413
414 Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
415 aeffbar[j] = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
416 daeffbar[j] = aeff->GetBinError(j,k0);
417
418 //*fLog << "Etru bin " << j << ": tbar= " << Thetabar
419 // << ", abar= " << aeffbar[j]
420 // << ", dabar= " << daeffbar[j] << endl;
421 }
422
423 //--- start m loop ---
424 // calculate effective collection area at (E = Ebar, Theta = Thetabar)
425 // by quadratic interpolation in log10(Etru)
426 // do this for each bin of Eunf
427 //
428 for (int m=1; m<=nEunf; m++)
429 {
430 Double_t log10Ebar = 0.5 *
431 ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m) )
432 +log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
433 Double_t Ebar = pow(10.0, log10Ebar);
434
435 // determine Etru bins (j1, j2, j3) for interpolation in E
436 // j0 denotes the Etru bin from which the error is copied
437 Int_t j0=0, j1=0, j2=0, j3=0;
438
439 for (int j=3; j<=nEtru; j++)
440 {
441 Double_t Elow = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(j);
442 if (Ebar < Elow)
443 {
444 j1 = j-2;
445 j2 = j-1;
446 j3 = j;
447 j0 = j2;
448 break;
449 }
450 }
451
452 if (j3 == 0)
453 {
454 j1 = nEtru-2;
455 j2 = nEtru-1;
456 j3 = nEtru;
457 j0 = j2;
458 }
459
460 if (Ebar < ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2))
461 j0 = 1;
462 else if (Ebar > ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru))
463 j0 = nEtru;
464
465 Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1);
466 Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru+1);
467 if (Ebar < Etrumin || Ebar > Etrumax)
468 {
469 *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
470 << ", Etrumin =" << Etrumin
471 << ", Etrumax =" << Etrumax << endl;
472 }
473
474 //*fLog << "Var bin " << n << ":" << endl;
475 //*fLog << "Ebar= " << Ebar
476 // << ", j1= " << j1
477 // << ", j2= " << j2
478 // << ", j3= " << j3 << endl;
479
480
481 c0=0.0;
482 c1=0.0;
483 c2=0.0;
484
485 t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1) )
486 +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) );
487
488 t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2) )
489 +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) );
490
491 t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3) )
492 +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) );
493
494
495 a1 = aeffbar[j1];
496 a2 = aeffbar[j2];
497 a3 = aeffbar[j3];
498
499 Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
500 aeffEunfVar = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
501 daeffEunfVar = daeffbar[j0];
502
503 //*fLog << "Eunf bin " << m << ": Ebar= " << Ebar
504 // << ", aeffEunfVar = " << aeffEunfVar
505 // << ", daeffEunfVar = " << daeffEunfVar << endl;
506
507 fHAeff.SetBinContent(m,n,aeffEunfVar);
508 fHAeff.SetBinError(m,n,daeffEunfVar);
509 }
510 //--- end m loop ---
511 }
512 //------ end n loop ------
513 delete aeffbar;
514
515 //...................................................................
516 // now calculate the absolute gamma flux
517 //
518 for (int m=1; m<=nEunf; m++)
519 {
520 Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
521
522 for (int n=1; n<=nVar; n++)
523 {
524 Double_t Ngam = fHUnfold.GetBinContent(m,n);
525 Double_t dNgam = fHUnfold.GetBinError(m,n);
526
527 Double_t Aeff = fHAeff.GetBinContent(m,n);
528 Double_t dAeff = fHAeff.GetBinError(m,n);
529
530 Double_t Effon = teff->GetBinContent(n);
531 Double_t dEffon = teff->GetBinError(n);
532
533 Double_t Cont, dCont;
534 if (Ngam > 0.0 && DeltaE > 0.0 && Effon > 0.0 && Aeff > 0.0)
535 {
536 Cont = Ngam / (DeltaE * Effon * Aeff);
537 dCont = Cont * sqrt( dNgam*dNgam / (Ngam*Ngam)
538 + dEffon*dEffon / (Effon*Effon)
539 + dAeff*dAeff / (Aeff*Aeff) );
540 }
541 else
542 {
543 Cont = 1.e-20;
544 dCont = 1.e-20;
545 }
546
547 fHFlux.SetBinContent(m,n,Cont);
548 fHFlux.SetBinError(m,n,dCont);
549
550 //*fLog << "Eunf bin " << m << ", Var bin " << n
551 // << ": Ngam = " << Ngam << ", Flux = "
552 // << Cont << ", dFlux = " << dCont << endl;
553 //*fLog << ", DeltaE = " << DeltaE << ", Effon = " << Effon
554 // << ", Aeff = " << Aeff << endl;
555 }
556 }
557
558 //..............................................
559 // draw the differential photon flux vs. E-unf
560 // for the individual bins of the variable Var
561
562 if (Draw == kTRUE)
563 {
564 for (int n=1; n<=nVar; n++)
565 {
566 TString strg0("Flux-");
567 strg0 += fVarname;
568 TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
569
570 char txt[100];
571 TString strg1("Photon flux vs. E-unfold for ");
572 TString strg2 = fVarname;
573 strg2 += "-bin %d";
574 sprintf(txt, strg2, n);
575 TString strg3 = strg1 + txt;
576
577 new TCanvas(txt, strg3);
578 // TCanvas &c = *MakeDefCanvas(txt, txt);
579 // gROOT->SetSelectedPad(NULL);
580
581 gPad->SetLogx();
582
583 h.SetName(txt);
584 h.SetTitle(strg3);
585 h.SetXTitle("E-unfold [GeV] ");
586 h.SetYTitle("photons / (s m2 GeV)");
587 h.DrawCopy();
588
589 // c.Modified();
590 // c.Update();
591 }
592 }
593 //........................
594}
595
596// -------------------------------------------------------------------------
597//
598// Draw copies of the histograms
599//
600TObject *MHFlux::DrawClone(Option_t *opt) const
601{
602 TCanvas &c = *MakeDefCanvas("flux", "Orig - Unfold - Flux plots");
603 c.Divide(2, 2);
604
605 gROOT->SetSelectedPad(NULL);
606
607 c.cd(1);
608 ((TH2*)&fHOrig)->DrawCopy("");
609 gPad->SetLogx();
610
611 c.cd(2);
612 ((TH2*)&fHUnfold)->DrawCopy("");
613 gPad->SetLogx();
614
615 c.cd(3);
616 ((TH2*)&fHFlux)->DrawCopy("");
617 gPad->SetLogx();
618
619 c.Modified();
620 c.Update();
621
622 return &c;
623}
624
625// -------------------------------------------------------------------------
626//
627// Draw the histograms
628//
629void MHFlux::Draw(Option_t *opt)
630{
631 if (!gPad)
632 MakeDefCanvas("flux", "orig-unfold-flux plots");
633
634 gPad->Divide(2,2);
635
636 gPad->cd(1);
637 fHOrig.Draw(opt);
638
639 gPad->cd(2);
640 fHUnfold.Draw(opt);
641
642 gPad->cd(3);
643 fHFlux.Draw(opt);
644
645 gPad->Modified();
646 gPad->Update();
647}
648
649// -------------------------------------------------------------------------
650//
651// Quadratic interpolation
652//
653// *** calculate the parameters of a parabula
654// y = a + b*x + c*x**2 = F(x)
655// such that yi = F(xi) for (i=1,3)
656//
657void MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
658 Double_t y1, Double_t y2, Double_t y3,
659 Double_t *a, Double_t *b, Double_t *c, Int_t *errflag)
660{
661 double ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33;
662 double det,det1;
663 //double yt1,yt2,yt3;
664
665 det = x2*x3*x3 + x1*x2*x2 + x3*x1*x1
666 - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
667
668 if (det != 0.0)
669 {
670 *errflag = 0;
671 det1 = 1.0/det;
672
673 ai11 = x2*x3*x3 - x3*x2*x2;
674 ai12 = x3*x1*x1 - x1*x3*x3;
675 ai13 = x1*x2*x2 - x2*x1*x1;
676
677 ai21 = x2*x2 - x3*x3;
678 ai22 = x3*x3 - x1*x1;
679 ai23 = x1*x1 - x2*x2;
680
681 ai31 = x3 - x2;
682 ai32 = x1 - x3;
683 ai33 = x2 - x1;
684
685 *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1;
686 *b = (ai21*y1 + ai22*y2 + ai23*y3) * det1;
687 *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1;
688
689 //yt1 = *a + *b * x1 + *c * x1*x1;
690 //yt2 = *a + *b * x2 + *c * x2*x2;
691 //yt3 = *a + *b * x3 + *c * x3*x3;
692
693 //*fLog << "x1 = " << x1 << ", x2 = " << x2 << ", x3 = " << x3 << endl;
694 //*fLog << "y1 = " << y1 << ", y2 = " << y2 << ", y3 = " << y3 << endl;
695 //*fLog << "yt1 = " << yt1 << ", yt2 = " << yt2
696 // << ", yt3 = " << yt3 << endl;
697 //*fLog << "*a = " << *a << ", *b = " << *b << ", *c= " << *c
698 // << ", *errflag = " << *errflag << endl;
699
700 return;
701 }
702
703 *errflag = 1;
704 *a = 0.0;
705 *b = 0.0;
706 *c = 0.0;
707 return;
708}
Note: See TracBrowser for help on using the repository browser.