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

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