Changeset 1604 for trunk/MagicSoft/Mars/mhist

Timestamp:

11/14/02 10:51:54 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHFlux.cc (modified) (18 diffs)
• MHFlux.h (modified) (1 diff)

trunk/MagicSoft/Mars/mhist/MHFlux.cc

@@ -63 +63 @@
 MHFlux::MHFlux(const TH2D &h2d,  const Bool_t Draw,
                const TString varname, const TString unit)
- : fHOrig(), fHUnfold(), fHFlux() 
+    : fHOrig(), fHUnfold(), fHFlux()
 {
-  //    if (varname == NULL  ||  unit == NULL)
-  if (varname == ""  ||  unit == "")
-    {
-      *fLog << "MHFlux : varname or unit not defined" << endl;
-    }
+    if (varname.IsNull() || unit.IsNull())
+        *fLog << warn << dbginf << "varname or unit not defined" << endl;
 
     fVarname = varname;
@@ -92 +89 @@
 
     fHOrig.SetDirectory(NULL);
-    fHOrig.SetXTitle("E-est [GeV]         ");
+    fHOrig.SetXTitle("E-est [GeV]");
     fHOrig.SetYTitle(strg);
     fHOrig.Sumw2();
@@ -109 +106 @@
 
     fHUnfold.SetDirectory(NULL);
-    fHUnfold.SetXTitle("E-unfold [GeV]         ");
+    fHUnfold.SetXTitle("E-unfold [GeV]");
     fHUnfold.SetYTitle(strg);
     fHUnfold.Sumw2();
@@ -127 +124 @@
 
     fHFlux.SetDirectory(NULL);
-    fHFlux.SetXTitle("E-unfold [GeV]         ");
+    fHFlux.SetXTitle("E-unfold [GeV]");
     fHFlux.SetYTitle(strg);
     fHFlux.Sumw2();
@@ -161 +158 @@
       TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");
 
-      char txt0[100];
       strg0  = fVarname;
-      strg0 += "-bin %d";
-      sprintf(txt0, strg0, n);
+      strg0 += "-bin ";
+      strg0 += n;
 
       TString strg1("No.of photons vs. E-est for ");
-      strg1 += txt0;
-      new TCanvas(txt0,strg1);
+      strg1 += strg0;
+
+      new TCanvas(strg0, strg1);
       // TCanvas &c = *MakeDefCanvas(txt0, strg1);
       // gROOT->SetSelectedPad(NULL);
@@ -174 +171 @@
       gPad->SetLogx();
 
-      h.SetName(txt0);
+      h.SetName(strg0);
       h.SetTitle(strg1);
-      h.SetXTitle("E-est [GeV]            ");
+      h.SetXTitle("E-est [GeV]");
       h.SetYTitle("No.of photons");
       h.DrawCopy();
@@ -226 +223 @@
       strg0 += fVarname;
       TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
-
-      char txt0[100];
       strg0  = fVarname;
-      strg0 += "-bin %d";
-      sprintf(txt0, strg0, n);
+      strg0 += "-bin ";
+      strg0 += n;
 
       TString strg1("No.of photons vs. E-unfold for ");
-      strg1 += txt0;
-      new TCanvas(txt0,strg1);
+      strg1 += strg0;
+
+      new TCanvas(strg0, strg1);
 
       // TCanvas &c = *MakeDefCanvas(txt0, strg1);
@@ -241 +237 @@
       gPad->SetLogx();
 
-      h.SetName(txt0);
+      h.SetName(strg0);
       h.SetTitle(strg1);
-      h.SetXTitle("E-unfold [GeV]            ");
+      h.SetXTitle("E-unfold [GeV]");
       h.SetYTitle("No.of photons");
       h.DrawCopy();
@@ -276 +272 @@
   //              for the individual bins of the variable Var
 
+    const TAxis &axex = *((TH2*)aeff)->GetXaxis();
+    const TAxis &axey = *((TH2*)aeff)->GetYaxis();
+
   //....................................
   // define dummy histogram *aeff
@@ -287 +286 @@
   SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
 
-  const Int_t netru    = aeff->GetNbinsX();
-  const Int_t ntheta   = aeff->GetNbinsY();
+  const Int_t netru  = aeff->GetNbinsX();
+  const Int_t ntheta = aeff->GetNbinsY();
 
   for (int j=1; j<=netru; j++)
   {
-    for (int k=1; k<=ntheta; k++)
-    {
-      Double_t cont = 10000.0;;
-      ((TH1*)aeff)->SetBinContent(j,k,cont);
-
-      Double_t dcont = 100.0;
-      ((TH1*)aeff)->SetBinError(j,k,dcont);
-    }
+      for (int k=1; k<=ntheta; k++)
+      {
+          Double_t cont = 10000.0;
+          ((TH1*)aeff)->SetBinContent(j, k, cont);
+
+          Double_t dcont = 100.0;
+          ((TH1*)aeff)->SetBinError(j, k, dcont);
+      }
   }
   // *fLog << "Dummy aeff : netru =" << netru << ",  ntheta = " << ntheta << endl;
@@ -327 +326 @@
   SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
 
-  Int_t    errflag;
-  Double_t c0, c1, c2;  
-  Double_t t1, t2, t3;
-  Double_t a1, a2, a3;
-
-  Double_t *aeffbar;
-  aeffbar = new Double_t[nEtru];
-  Double_t *daeffbar;
-  daeffbar = new Double_t[nEtru];
+  Double_t *aeffbar  = new Double_t[nEtru];
+  Double_t *daeffbar = new Double_t[nEtru];
 
   Double_t aeffEunfVar;
   Double_t daeffEunfVar;
-
 
   //------   start n loop   ------
@@ -349 +340 @@
     // determine Theta bins (k1, k2, k3) for interpolation in Theta
     // k0 denotes the Theta bin from whicvh the error is copied
-    Int_t k0=0, k1=0, k2=0, k3=0;
+    Int_t k0=0;
+    Int_t k1=0;
+    Int_t k2=0;
+    Int_t k3=0;
+
     for (int k=3; k<=nTheta; k++)
     {
-      Double_t Thetalow = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(k);
-      if (Thetabar < Thetalow)
-      {
-        k1 = k-2;
-        k2 = k-1;
-        k3 = k;
+        Double_t Thetalow = axey.GetBinLowEdge(k);
+        if (Thetabar < Thetalow)
+        {
+            k1 = k-2;
+            k2 = k-1;
+            k3 = k;
+            k0 = k2;
+            break;
+        }
+    }
+
+    if (k3 == 0)
+    {
+        k1 = nTheta-2;
+        k2 = nTheta-1;
+        k3 = nTheta;
         k0 = k2;
-        break;
-      }
-    }  
-
-    if (k3 == 0)
-    {
-      k1 = nTheta-2;
-      k2 = nTheta-1;
-      k3 = nTheta;
-      k0 = k2;
     }
 
-    if (Thetabar <  ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2))
+    if (Thetabar <  axey.GetBinLowEdge(2))
       k0 = 1;
-    else if (Thetabar >  ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta))
-      k0 = nTheta;
-
-    Double_t Thetamin = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(1);   
-    Double_t Thetamax = ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(nTheta+1);   
+    else
+        if (Thetabar >  axey.GetBinLowEdge(nTheta))
+            k0 = nTheta;
+
+    Double_t Thetamin = axey.GetBinLowEdge(1);
+    Double_t Thetamax = axey.GetBinLowEdge(nTheta+1);
     if (Thetabar < Thetamin  ||  Thetabar > Thetamax)
     {
@@ -400 +396 @@
     for (int j=1; j<=nEtru; j++)
     {
-      c0 = 0.0;
-      c1 = 0.0;
-      c2 = 0.0;
-
-      t1 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k1) );
-      t2 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k2) );
-      t3 = cos( ((TH1*)aeff)->GetYaxis()->GetBinCenter (k3) );
-
-      a1 = aeff->GetBinContent(j,k1);
-      a2 = aeff->GetBinContent(j,k2);
-      a3 = aeff->GetBinContent(j,k3);
-
-      Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
-      aeffbar[j]  = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
-      daeffbar[j] = aeff->GetBinError(j,k0);
-
-      //*fLog << "Etru bin " << j <<  ":  tbar= " << Thetabar 
-      //      << ",  abar= "  << aeffbar[j] 
-      //      << ",  dabar= " << daeffbar[j] << endl;
+        double c0 = 0;
+        double c1 = 0;
+        double c2 = 0;
+
+        const double t1 = cos( axey.GetBinCenter (k1) );
+        const double t2 = cos( axey.GetBinCenter (k2) );
+        const double t3 = cos( axey.GetBinCenter (k3) );
+
+        const double a1 = aeff->GetBinContent(j, k1);
+        const double a2 = aeff->GetBinContent(j, k2);
+        const double a3 = aeff->GetBinContent(j, k3);
+
+        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
+        aeffbar[j]  = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
+        daeffbar[j] = aeff->GetBinError(j,k0);
+
+        //*fLog << "Etru bin " << j <<  ":  tbar= " << Thetabar
+        //      << ",  abar= "  << aeffbar[j]
+        //      << ",  dabar= " << daeffbar[j] << endl;
     }
 
@@ -428 +424 @@
     for (int m=1; m<=nEunf; m++)
     {
-      Double_t log10Ebar = 0.5 * 
-               ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m)  )
-                +log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
-      Double_t Ebar = pow(10.0, log10Ebar);
-
-      // determine Etru bins (j1, j2, j3) for interpolation in E
-      // j0 denotes the Etru bin from which the error is copied
-      Int_t j0=0, j1=0, j2=0, j3=0;
-
-      for (int j=3; j<=nEtru; j++)
-      {
-        Double_t Elow = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(j);
-        if (Ebar < Elow)
+        Double_t log10Ebar = 0.5 * ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m)  )+
+                                     log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
+        Double_t Ebar = pow(10.0, log10Ebar);
+
+        // determine Etru bins (j1, j2, j3) for interpolation in E
+        // j0 denotes the Etru bin from which the error is copied
+        Int_t j0=0;
+        Int_t j1=0;
+        Int_t j2=0;
+        Int_t j3=0;
+
+        for (int j=3; j<=nEtru; j++)
         {
-          j1 = j-2;
-          j2 = j-1;
-          j3 = j;
-          j0 = j2;
-          break;
+            Double_t Elow = axex.GetBinLowEdge(j);
+            if (Ebar < Elow)
+            {
+                j1 = j-2;
+                j2 = j-1;
+                j3 = j;
+                j0 = j2;
+                break;
+            }
         }
-      }  
-
-      if (j3 == 0)
-      {
-        j1 = nEtru-2;
-        j2 = nEtru-1;
-        j3 = nEtru;
-        j0 = j2;
-      }
-
-      if (Ebar <  ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(2))
-        j0 = 1;
-      else if (Ebar >  ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru))
-        j0 = nEtru;
-
-      Double_t Etrumin = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(1);   
-      Double_t Etrumax = ((TH1*)aeff)->GetXaxis()->GetBinLowEdge(nEtru+1);   
-      if (Ebar < Etrumin  ||  Ebar > Etrumax)
-      {
-        *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
-              << ",  Etrumin =" << Etrumin
-              << ",  Etrumax =" << Etrumax << endl;
-      } 
-
-      //*fLog << "Var bin "   << n  << ":"  <<  endl;
-      //*fLog << "Ebar= " << Ebar 
-      //      << ",  j1= "    << j1
-      //      << ",  j2= "    << j2
-      //      << ",  j3= "    << j3         <<  endl;
-
-
-      c0=0.0; 
-      c1=0.0;
-      c2=0.0;  
-
-      t1 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1)  )
-                  +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j1+1)) );
-
-      t2 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2)  )
-                  +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j2+1)) );
-
-      t3 = 0.5 * ( log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3)  )
-                  +log10( ((TH1*)aeff)->GetXaxis()->GetBinLowEdge (j3+1)) );
-
-
-      a1 = aeffbar[j1];
-      a2 = aeffbar[j2];
-      a3 = aeffbar[j3];
-
-      Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2, &errflag);
-      aeffEunfVar  = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
-      daeffEunfVar = daeffbar[j0];
-
-      //*fLog << "Eunf bin " << m     <<  ":  Ebar= " << Ebar 
-      //      << ",  aeffEunfVar = "  << aeffEunfVar 
-      //      << ",  daeffEunfVar = " << daeffEunfVar << endl;
-
-      fHAeff.SetBinContent(m,n,aeffEunfVar);
-      fHAeff.SetBinError(m,n,daeffEunfVar);
+
+        if (j3 == 0)
+        {
+            j1 = nEtru-2;
+            j2 = nEtru-1;
+            j3 = nEtru;
+            j0 = j2;
+        }
+
+        if (Ebar <  axex.GetBinLowEdge(2))
+            j0 = 1;
+        else
+            if (Ebar >  axex.GetBinLowEdge(nEtru))
+                j0 = nEtru;
+
+        Double_t Etrumin = axex.GetBinLowEdge(1);
+        Double_t Etrumax = axex.GetBinLowEdge(nEtru+1);
+        if (Ebar < Etrumin  ||  Ebar > Etrumax)
+        {
+            *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
+                << ",  Etrumin =" << Etrumin
+                << ",  Etrumax =" << Etrumax << endl;
+        }
+
+        //*fLog << "Var bin "   << n  << ":"  <<  endl;
+        //*fLog << "Ebar= " << Ebar
+        //      << ",  j1= "    << j1
+        //      << ",  j2= "    << j2
+        //      << ",  j3= "    << j3         <<  endl;
+
+
+        double c0=0.0;
+        double c1=0.0;
+        double c2=0.0;
+
+        const double t1 = 0.5 * ( log10( axex.GetBinLowEdge (j1)  )+
+                                  log10( axex.GetBinLowEdge (j1+1)) );
+        const double t2 = 0.5 * ( log10( axex.GetBinLowEdge (j2)  )+
+                                  log10( axex.GetBinLowEdge (j2+1)) );
+        const double t3 = 0.5 * ( log10( axex.GetBinLowEdge (j3)  )+
+                                  log10( axex.GetBinLowEdge (j3+1)) );
+
+        const double a1 = aeffbar[j1];
+        const double a2 = aeffbar[j2];
+        const double a3 = aeffbar[j3];
+
+        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
+        aeffEunfVar  = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
+        daeffEunfVar = daeffbar[j0];
+
+        //*fLog << "Eunf bin " << m     <<  ":  Ebar= " << Ebar
+        //      << ",  aeffEunfVar = "  << aeffEunfVar
+        //      << ",  daeffEunfVar = " << daeffEunfVar << endl;
+
+        fHAeff.SetBinContent(m,n,aeffEunfVar);
+        fHAeff.SetBinError(m,n,daeffEunfVar);
     }
     //---   end m loop ---
@@ -518 +514 @@
   for (int m=1; m<=nEunf; m++)
   {
-    Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
-
-    for (int n=1; n<=nVar; n++)
-    {
-      Double_t Ngam   = fHUnfold.GetBinContent(m,n);
-      Double_t dNgam  = fHUnfold.GetBinError(m,n);
-
-      Double_t Aeff   = fHAeff.GetBinContent(m,n);
-      Double_t dAeff  = fHAeff.GetBinError(m,n);
-
-      Double_t Effon  = teff->GetBinContent(n);
-      Double_t dEffon = teff->GetBinError(n);
-
-      Double_t Cont, dCont;
-      if (Ngam > 0.0  &&  DeltaE > 0.0  &&  Effon > 0.0  &&  Aeff > 0.0)
+      Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
+
+      for (int n=1; n<=nVar; n++)
       {
-        Cont  = Ngam / (DeltaE * Effon * Aeff);
-        dCont = Cont * sqrt(   dNgam*dNgam   / (Ngam*Ngam)
-                             + dEffon*dEffon / (Effon*Effon)
-                             + dAeff*dAeff   / (Aeff*Aeff)  );  
+          Double_t Ngam   = fHUnfold.GetBinContent(m,n);
+          Double_t dNgam  = fHUnfold.GetBinError(m,n);
+
+          Double_t Aeff   = fHAeff.GetBinContent(m,n);
+          Double_t dAeff  = fHAeff.GetBinError(m,n);
+
+          Double_t Effon  = teff->GetBinContent(n);
+          Double_t dEffon = teff->GetBinError(n);
+
+          Double_t Cont, dCont;
+          if (Ngam>0 && DeltaE>0 && Effon>0 && Aeff>0)
+          {
+              Cont  = Ngam / (DeltaE * Effon * Aeff);
+              dCont = Cont * sqrt( dNgam *dNgam  / (Ngam*Ngam) +
+                                   dEffon*dEffon / (Effon*Effon) +
+                                   dAeff *dAeff  / (Aeff*Aeff)  );
+          }
+          else
+          {
+              Cont  = 1.e-20;
+              dCont = 1.e-20;
+          }
+
+          fHFlux.SetBinContent(m,n,Cont);
+          fHFlux.SetBinError(m,n,dCont);
+
+          //*fLog << "Eunf bin "    << m      << ",  Var bin " << n
+          //      << ":  Ngam = "   << Ngam   << ",  Flux = "
+          //      << Cont  << ", dFlux = " << dCont << endl;
+          //*fLog << ",  DeltaE = " << DeltaE << ",  Effon = " << Effon
+          //      << ",  Aeff = "   << Aeff   << endl;
       }
-      else
-      {
-        Cont  = 1.e-20;
-        dCont = 1.e-20;
-      }
-
-      fHFlux.SetBinContent(m,n,Cont);
-      fHFlux.SetBinError(m,n,dCont);
-
-      //*fLog << "Eunf bin "    << m      << ",  Var bin " << n 
-      //      << ":  Ngam = "   << Ngam   << ",  Flux = "  
-      //      << Cont  << ", dFlux = " << dCont << endl; 
-      //*fLog << ",  DeltaE = " << DeltaE << ",  Effon = " << Effon
-      //      << ",  Aeff = "   << Aeff   << endl;
-    }
   }
 
@@ -562 +558 @@
   if (Draw == kTRUE)
   {
-    for (int n=1; n<=nVar; n++)
-    {
-      TString strg0("Flux-");
-      strg0 += fVarname;
-      TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
-
-      char txt[100];
-      TString strg1("Photon flux vs. E-unfold for ");
-      TString strg2 = fVarname;
-      strg2 += "-bin %d";
-      sprintf(txt, strg2, n);
-      TString strg3 = strg1 + txt;
-  
-      new TCanvas(txt, strg3);
-      // TCanvas &c = *MakeDefCanvas(txt, txt);
-      // gROOT->SetSelectedPad(NULL);
-
-      gPad->SetLogx();
-
-      h.SetName(txt);
-      h.SetTitle(strg3);
-      h.SetXTitle("E-unfold [GeV]            ");
-      h.SetYTitle("photons / (s m2 GeV)");
-      h.DrawCopy();
-
-      // c.Modified();
-      // c.Update();
-    }
+      for (int n=1; n<=nVar; n++)
+      {
+          TString strg0("Flux-");
+          strg0 += fVarname;
+
+          TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
+
+          TString strg1("Photon flux vs. E-unfold for ");
+          TString strg2 = fVarname;
+
+          strg2 += "-bin ";
+          strg2 += n;
+
+          TString strg3 = strg1 + strg2;
+          new TCanvas(strg2, strg3);
+          // TCanvas &c = *MakeDefCanvas(txt, txt);
+          // gROOT->SetSelectedPad(NULL);
+
+          gPad->SetLogx();
+
+          h.SetName(strg2);
+          h.SetTitle(strg3);
+          h.SetXTitle("E-unfold [GeV]            ");
+          h.SetYTitle("photons / (s m2 GeV)");
+          h.DrawCopy();
+
+          // c.Modified();
+          // c.Update();
+      }
   }
   //........................
@@ -655 +652 @@
 //     such that       yi = F(xi)       for (i=1,3)
 //
-void MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
-                   Double_t y1, Double_t y2, Double_t y3,
-                   Double_t *a, Double_t *b, Double_t *c, Int_t *errflag)
+Bool_t MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
+                     Double_t y1, Double_t y2, Double_t y3,
+                     Double_t *a, Double_t *b, Double_t *c)
 {
-  double  ai11,ai12,ai13,ai21,ai22,ai23,ai31,ai32,ai33;
-  double  det,det1;
-  //double  yt1,yt2,yt3;
-
-  det =   x2*x3*x3 + x1*x2*x2 + x3*x1*x1
+    const double det =
+        + x2*x3*x3 + x1*x2*x2 + x3*x1*x1
         - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
 
-  if (det != 0.0)
-  {
-    *errflag = 0;
-    det1 = 1.0/det;
-
-    ai11 = x2*x3*x3 - x3*x2*x2;
-    ai12 = x3*x1*x1 - x1*x3*x3;
-    ai13 = x1*x2*x2 - x2*x1*x1;
-
-    ai21 = x2*x2 - x3*x3;
-    ai22 = x3*x3 - x1*x1;
-    ai23 = x1*x1 - x2*x2;
-
-    ai31 = x3 - x2;
-    ai32 = x1 - x3;
-    ai33 = x2 - x1;
+    if (det == 0.0)
+    {
+        *a = 0;
+        *b = 0;
+        *c = 0;
+        return kFALSE;
+    }
+
+    const double det1 = 1.0/det;
+
+    const double ai11 = x2*x3*x3 - x3*x2*x2;
+    const double ai12 = x3*x1*x1 - x1*x3*x3;
+    const double ai13 = x1*x2*x2 - x2*x1*x1;
+
+    const double ai21 = x2*x2 - x3*x3;
+    const double ai22 = x3*x3 - x1*x1;
+    const double ai23 = x1*x1 - x2*x2;
+
+    const double ai31 = x3 - x2;
+    const double ai32 = x1 - x3;
+    const double ai33 = x2 - x1;
 
     *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1;
@@ -690 +689 @@
     //yt2 = *a + *b * x2 + *c * x2*x2;
     //yt3 = *a + *b * x3 + *c * x3*x3;
-    
-    //*fLog << "x1 = " << x1 << ",  x2 = " << x2 << ",  x3 = " << x3 << endl; 
-    //*fLog << "y1 = " << y1 << ",  y2 = " << y2 << ",  y3 = " << y3 << endl; 
-    //*fLog << "yt1 = " << yt1 << ",  yt2 = " << yt2  
-    //    << ",  yt3 = " << yt3 << endl; 
-    //*fLog << "*a = " << *a << ",  *b = " << *b << ",  *c= " << *c 
-    //      << ",  *errflag = " << *errflag << endl; 
-
-    return;
-  }
-
-  *errflag = 1;
-  *a = 0.0;
-  *b = 0.0;
-  *c = 0.0;
-  return;
+
+    //*fLog << "x1 = " << x1 << ",  x2 = " << x2 << ",  x3 = " << x3 << endl;
+    //*fLog << "y1 = " << y1 << ",  y2 = " << y2 << ",  y3 = " << y3 << endl;
+    //*fLog << "yt1 = " << yt1 << ",  yt2 = " << yt2
+    //    << ",  yt3 = " << yt3 << endl;
+    //*fLog << "*a = " << *a << ",  *b = " << *b << ",  *c= " << *c
+    //      << ",  *errflag = " << *errflag << endl;
+
+    return kTRUE;
 }

trunk/MagicSoft/Mars/mhist/MHFlux.h

@@ -52 +52 @@
     const TH2D *GetHFlux()       { return &fHFlux; }
 
-    void Parab(double x1, double x2, double x3,
-               double y1, double y2, double y3,
-               double *a, double *b, double *c, int *errflag);
+    static Bool_t Parab(double x1, double x2, double x3,
+                        double y1, double y2, double y3,
+                        double *a, double *b, double *c);
 
     ClassDef(MHFlux, 1) //2D-plots (original, unfolded, flux)