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

Last change on this file since 1578 was 1578, checked in by wittek, 22 years ago
*** empty log message ***
File size: 20.1 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// (abstract member function));
194// without the dummy Fill one gets the error message :
195//
196// Error: Can't call MHFlux::MHFlux(evttime,"time","[s]") in current scope
197// FILE:macros/flux.C LINE:465
198// Possible candidates are...
199// filename line:size busy function type and name (in MHFlux)
200// filename line:size busy function type and name (in MH)
201// filename line:size busy function type and name (in MParContainer)
202// filename line:size busy function type and name (in TObject)
203//
204Bool_t MHFlux::Fill(const MParContainer *par)
205{
206 return kTRUE;
207}
208
209
210// -------------------------------------------------------------------------
211//
212// Unfold the distribution in E-est
213//
214void MHFlux::Unfold(const Bool_t Draw)
215{
216 //..............................................
217 // draw the No.of photons vs. E-unfold
218 // for the individual bins of the variable Var
219
220 if (Draw == kTRUE)
221 {
222 const Int_t nVar = fHUnfold.GetNbinsY();
223
224 for (int n=1; n<=nVar; n++)
225 {
226 TString strg0("Unfold-");
227 strg0 += fVarname;
228 TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
229
230 char txt0[100];
231 strg0 = fVarname;
232 strg0 += "-bin %d";
233 sprintf(txt0, strg0, n);
234
235 TString strg1("No.of photons vs. E-unfold for ");
236 strg1 += txt0;
237 new TCanvas(txt0,strg1);
238
239 // TCanvas &c = *MakeDefCanvas(txt0, strg1);
240 // gROOT->SetSelectedPad(NULL);
241
242 gPad->SetLogx();
243
244 h.SetName(txt0);
245 h.SetTitle(strg1);
246 h.SetXTitle("E-unfold [GeV] ");
247 h.SetYTitle("No.of photons");
248 h.DrawCopy();
249
250 // c.Modified();
251 // c.Update();
252 }
253 }
254 //........................
255}
256
257
258// -------------------------------------------------------------------------
259//
260// Calculate photon flux by dividing the distribution in Eunf (fHUnfold) by
261// the width of the energy interval (deltaE)
262// the effective ontime (*teff)
263// and the effective collection area (*aeff)
264//
265void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar,
266 const TH2D *aeff, const Bool_t Draw)
267{
268 // Note that fHUnfold has bins in Eunf and Var
269 // *teff has bins in Var (the same bins in Var as fHUnfold)
270 // *thetabar has bins in Var (the same bins in Var as fHUnfold)
271 // *aeff has bins in Etru and Theta
272 // (where in general the binning in Etru is different
273 // from the binning in Eunf)
274 // The variable Var may be 'time' or 'Theta'
275
276 // Draw = kTRUE means the differential flux vs E-unf should be drawn
277 // for the individual bins of the variable Var
278
279 //....................................
280 // define dummy histogram *aeff
281 ((TH1*)aeff)->Sumw2();
282 MBinning binsetru("BinningEtru");
283 binsetru.SetEdgesLog(10, 10, 1e3);
284
285 MBinning binsthetatru("BinningThetatru");
286 binsthetatru.SetEdges(7, -2.5, 32.5);
287 //SetBinning((TH1*)aeff, &binsetru, &binsthetatru);
288 SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
289
290 const Int_t netru = aeff->GetNbinsX();
291 const Int_t ntheta = aeff->GetNbinsY();
292
293 for (int j=1; j<=netru; j++)
294 {
295 for (int k=1; k<=ntheta; k++)
296 {
297 Double_t cont = 10000.0;;
298 ((TH1*)aeff)->SetBinContent(j,k,cont);
299
300 Double_t dcont = 100.0;
301 ((TH1*)aeff)->SetBinError(j,k,dcont);
302 }
303 }
304 // *fLog << "Dummy aeff : netru =" << netru << ", ntheta = " << ntheta << endl;
305 //....................................
306
307 // number of Eunf and Var bins (histograms : fHUnfold, fHFlux)
308 const Int_t nEunf = fHFlux.GetNbinsX();
309 const Int_t nVar = fHFlux.GetNbinsY();
310
311 // number of Etru and Theta bins (histogram *aeff of collection area)
312
313 const Int_t nEtru = aeff->GetNbinsX();
314 const Int_t nTheta = aeff->GetNbinsY();
315
316 //*fLog << "nEunf =" << nEunf << ", nVar = " << nVar << endl;
317 //*fLog << "nEtru =" << nEtru << ", nTheta = " << nTheta << endl;
318
319 //...................................................................
320 // calculate effective collection area
321 // for the Eunf and Var bins of the histogram fHUnfold
322 // from the histogram *aeff, which has bins in Etru and Theta
323 // the result is the histogram fHAeff
324 //
325
326 TH2D fHAeff;
327 fHAeff.Sumw2();
328 SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
329
330 Int_t errflag;
331 Double_t c0, c1, c2;
332 Double_t t1, t2, t3;
333 Double_t a1, a2, a3;
334
335 Double_t *aeffbar;
336 aeffbar = new Double_t[nEtru];
337 Double_t *daeffbar;
338 daeffbar = new Double_t[nEtru];
339
340 Double_t aeffEunfVar;
341 Double_t daeffEunfVar;
342
343
344 //------ start n loop ------
345 for (int n=1; n<=nVar; n++)
346 {
347 Double_t Thetabar = thetabar->GetBinContent(n);
348 Double_t cosThetabar = cos(Thetabar);
349
350 // determine Theta bins (k1, k2, k3) for interpolation in Theta
351 // k0 denotes the Theta bin from whicvh the error is copied
352 Int_t k0=0, k1=0, k2=0, k3=0;
353 for (int k=3; k<=nTheta; k++)
354 {
355 Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(k);
356 if (Thetabar < Thetalow)
357 {
358 k1 = k-2;
359 k2 = k-1;
360 k3 = k;
361 k0 = k2;
362 break;
363 }
364 }
365
366 if (k3 == 0)
367 {
368 k1 = nTheta-2;
369 k2 = nTheta-1;
370 k3 = nTheta;
371 k0 = k2;
372 }
373
374 if (Thetabar < ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2))
375 k0 = 1;
376 else if (Thetabar > ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta))
377 k0 = nTheta;
378
379 Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1);
380 Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta+1);
381 if (Thetabar < Thetamin || Thetabar > Thetamax)
382 {
383 *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar
384 << ", Thetamin =" << Thetamin
385 << ", Thetamax =" << Thetamax << endl;
386 }
387
388 //*fLog << "Var bin " << n << ":" << endl;
389 //*fLog << "Thetabar= " << Thetabar
390 // << ", k0= " << k0
391 // << ", k1= " << k1
392 // << ", k2= " << k2
393 // << ", k3= " << k3 << endl;
394
395
396 // calculate effective collection area at Theta = Thetabar
397 // by quadratic interpolation in cos(Theta);
398 // do this for each bin of Etru
399 //
400
401 for (int j=1; j<=nEtru; j++)
402 {
403 c0 = 0.0;
404 c1 = 0.0;
405 c2 = 0.0;
406
407 t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) );
408 t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) );
409 t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) );
410
411 a1 = aeff->GetBinContent(j,k1);
412 a2 = aeff->GetBinContent(j,k2);
413 a3 = aeff->GetBinContent(j,k3);
414
415 Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
416 aeffbar[j] = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
417 daeffbar[j] = aeff->GetBinError(j,k0);
418
419 //*fLog << "Etru bin " << j << ": tbar= " << Thetabar
420 // << ", abar= " << aeffbar[j]
421 // << ", dabar= " << daeffbar[j] << endl;
422 }
423
424 //--- start m loop ---
425 // calculate effective collection area at (E = Ebar, Theta = Thetabar)
426 // by quadratic interpolation in log10(Etru)
427 // do this for each bin of Eunf
428 //
429 for (int m=1; m<=nEunf; m++)
430 {
431 Double_t log10Ebar = 0.5 *
432 ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m) )
433 +log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
434 Double_t Ebar = pow(10.0, log10Ebar);
435
436 // determine Etru bins (j1, j2, j3) for interpolation in E
437 // j0 denotes the Etru bin from which the error is copied
438 Int_t j0=0, j1=0, j2=0, j3=0;
439
440 for (int j=3; j<=nEtru; j++)
441 {
442 Double_t Elow = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(j);
443 if (Ebar < Elow)
444 {
445 j1 = j-2;
446 j2 = j-1;
447 j3 = j;
448 j0 = j2;
449 break;
450 }
451 }
452
453 if (j3 == 0)
454 {
455 j1 = nEtru-2;
456 j2 = nEtru-1;
457 j3 = nEtru;
458 j0 = j2;
459 }
460
461 if (Ebar < ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2))
462 j0 = 1;
463 else if (Ebar > ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru))
464 j0 = nEtru;
465
466 Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1);
467 Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru+1);
468 if (Ebar < Etrumin || Ebar > Etrumax)
469 {
470 *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
471 << ", Etrumin =" << Etrumin
472 << ", Etrumax =" << Etrumax << endl;
473 }
474
475 //*fLog << "Var bin " << n << ":" << endl;
476 //*fLog << "Ebar= " << Ebar
477 // << ", j1= " << j1
478 // << ", j2= " << j2
479 // << ", j3= " << j3 << endl;
480
481
482 c0=0.0;
483 c1=0.0;
484 c2=0.0;
485
486 t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1) )
487 +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) );
488
489 t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2) )
490 +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) );
491
492 t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3) )
493 +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) );
494
495
496 a1 = aeffbar[j1];
497 a2 = aeffbar[j2];
498 a3 = aeffbar[j3];
499
500 Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
501 aeffEunfVar = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
502 daeffEunfVar = daeffbar[j0];
503
504 //*fLog << "Eunf bin " << m << ": Ebar= " << Ebar
505 // << ", aeffEunfVar = " << aeffEunfVar
506 // << ", daeffEunfVar = " << daeffEunfVar << endl;
507
508 fHAeff.SetBinContent(m,n,aeffEunfVar);
509 fHAeff.SetBinError(m,n,daeffEunfVar);
510 }
511 //--- end m loop ---
512 }
513 //------ end n loop ------
514 delete aeffbar;
515
516 //...................................................................
517 // now calculate the absolute gamma flux
518 //
519 for (int m=1; m<=nEunf; m++)
520 {
521 Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
522
523 for (int n=1; n<=nVar; n++)
524 {
525 Double_t Ngam = fHUnfold.GetBinContent(m,n);
526 Double_t dNgam = fHUnfold.GetBinError(m,n);
527
528 Double_t Aeff = fHAeff.GetBinContent(m,n);
529 Double_t dAeff = fHAeff.GetBinError(m,n);
530
531 Double_t Effon = teff->GetBinContent(n);
532 Double_t dEffon = teff->GetBinError(n);
533
534 Double_t Cont, dCont;
535 if (Ngam > 0.0 && DeltaE > 0.0 && Effon > 0.0 && Aeff > 0.0)
536 {
537 Cont = Ngam / (DeltaE * Effon * Aeff);
538 dCont = Cont * sqrt( dNgam*dNgam / (Ngam*Ngam)
539 + dEffon*dEffon / (Effon*Effon)
540 + dAeff*dAeff / (Aeff*Aeff) );
541 }
542 else
543 {
544 Cont = 1.e-20;
545 dCont = 1.e-20;
546 }
547
548 fHFlux.SetBinContent(m,n,Cont);
549 fHFlux.SetBinError(m,n,dCont);
550
551 //*fLog << "Eunf bin " << m << ", Var bin " << n
552 // << ": Ngam = " << Ngam << ", Flux = "
553 // << Cont << ", dFlux = " << dCont << endl;
554 //*fLog << ", DeltaE = " << DeltaE << ", Effon = " << Effon
555 // << ", Aeff = " << Aeff << endl;
556 }
557 }
558
559 //..............................................
560 // draw the differential photon flux vs. E-unf
561 // for the individual bins of the variable Var
562
563 if (Draw == kTRUE)
564 {
565 for (int n=1; n<=nVar; n++)
566 {
567 TString strg0("Flux-");
568 strg0 += fVarname;
569 TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
570
571 char txt[100];
572 TString strg1("Photon flux vs. E-unfold for ");
573 TString strg2 = fVarname;
574 strg2 += "-bin %d";
575 sprintf(txt, strg2, n);
576 TString strg3 = strg1 + txt;
577
578 new TCanvas(txt, strg3);
579 // TCanvas &c = *MakeDefCanvas(txt, txt);
580 // gROOT->SetSelectedPad(NULL);
581
582 gPad->SetLogx();
583
584 h.SetName(txt);
585 h.SetTitle(strg3);
586 h.SetXTitle("E-unfold [GeV] ");
587 h.SetYTitle("photons / (s m2 GeV)");
588 h.DrawCopy();
589
590 // c.Modified();
591 // c.Update();
592 }
593 }
594 //........................
595}
596
597// -------------------------------------------------------------------------
598//
599// Draw copies of the histograms
600//
601TObject *MHFlux::DrawClone(Option_t *opt) const
602{
603 TCanvas &c = *MakeDefCanvas("flux", "Orig - Unfold - Flux plots");
604 c.Divide(2, 2);
605
606 gROOT->SetSelectedPad(NULL);
607
608 c.cd(1);
609 ((TH2*)&fHOrig)->DrawCopy("");
610 gPad->SetLogx();
611
612 c.cd(2);
613 ((TH2*)&fHUnfold)->DrawCopy("");
614 gPad->SetLogx();
615
616 c.cd(3);
617 ((TH2*)&fHFlux)->DrawCopy("");
618 gPad->SetLogx();
619
620 c.Modified();
621 c.Update();
622
623 return &c;
624}
625
626// -------------------------------------------------------------------------
627//
628// Draw the histograms
629//
630void MHFlux::Draw(Option_t *opt)
631{
632 if (!gPad)
633 MakeDefCanvas("flux", "orig-unfold-flux plots");
634
635 gPad->Divide(2,2);
636
637 gPad->cd(1);
638 fHOrig.Draw(opt);
639
640 gPad->cd(2);
641 fHUnfold.Draw(opt);
642
643 gPad->cd(3);
644 fHFlux.Draw(opt);
645
646 gPad->Modified();
647 gPad->Update();
648}
649
650// -------------------------------------------------------------------------
651//
652// Quadratic interpolation
653//
654// *** calculate the parameters of a parabula
655// y = a + b*x + c*x**2 = F(x)
656// such that yi = F(xi) for (i=1,3)
657//
658void MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
659 Double_t y1, Double_t y2, Double_t y3,
660 Double_t *a, Double_t *b, Double_t *c, Int_t *errflag)
661{
662 double ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33;
663 double det,det1;
664 //double yt1,yt2,yt3;
665
666 det = x2*x3*x3 + x1*x2*x2 + x3*x1*x1
667 - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
668
669 if (det != 0.0)
670 {
671 *errflag = 0;
672 det1 = 1.0/det;
673
674 ai11 = x2*x3*x3 - x3*x2*x2;
675 ai12 = x3*x1*x1 - x1*x3*x3;
676 ai13 = x1*x2*x2 - x2*x1*x1;
677
678 ai21 = x2*x2 - x3*x3;
679 ai22 = x3*x3 - x1*x1;
680 ai23 = x1*x1 - x2*x2;
681
682 ai31 = x3 - x2;
683 ai32 = x1 - x3;
684 ai33 = x2 - x1;
685
686 *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1;
687 *b = (ai21*y1 + ai22*y2 + ai23*y3) * det1;
688 *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1;
689
690 //yt1 = *a + *b * x1 + *c * x1*x1;
691 //yt2 = *a + *b * x2 + *c * x2*x2;
692 //yt3 = *a + *b * x3 + *c * x3*x3;
693
694 //*fLog << "x1 = " << x1 << ", x2 = " << x2 << ", x3 = " << x3 << endl;
695 //*fLog << "y1 = " << y1 << ", y2 = " << y2 << ", y3 = " << y3 << endl;
696 //*fLog << "yt1 = " << yt1 << ", yt2 = " << yt2
697 // << ", yt3 = " << yt3 << endl;
698 //*fLog << "*a = " << *a << ", *b = " << *b << ", *c= " << *c
699 // << ", *errflag = " << *errflag << endl;
700
701 return;
702 }
703
704 *errflag = 1;
705 *a = 0.0;
706 *b = 0.0;
707 *c = 0.0;
708 return;
709}
Note: See TracBrowser for help on using the repository browser.