Changeset 1604 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
11/14/02 10:51:54 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r1603 r1604  
    55   * mmain/MAnalysis.cc, mmain/MMonteCarlo.cc, mmain/MDataCheck.cc:
    66     - removed SetProgressBar of reader
     7
     8   * mhist/MHFlux.cc:
     9     - localized some variables
     10     - get rid of old c-style sprintf
     11     - return errorflag in Parab as return value!
    712
    813
  • trunk/MagicSoft/Mars/mgeom/MGeomMirror.cc

    r1600 r1604  
    3939// Initializes a Mirror geometry with 0 values, except for fMirrorID.
    4040//
    41 MGeomMirror::MGeomMirror(Int_t mir, const char *name=NULL, const char *title=NULL)
     41MGeomMirror::MGeomMirror(Int_t mir, const char *name, const char *title)
    4242{
    4343    fName  = name  ? name  : "MGeomMirror";
     
    7474//
    7575void MGeomMirror::SetMirrorContent(Int_t mir, Float_t focal, Float_t curv_x,
    76                                    Float_t curv_y, Float_t lin_x, Float_t lin_y,
     76                                   Float_t curv_y, Float_t lin_x, Float_t lin_y,
    7777                                   Float_t lin_z, Float_t theta, Float_t phi,
    7878                                   Float_t x_n, Float_t y_n, Float_t z_n,
  • trunk/MagicSoft/Mars/mhist/MHFlux.cc

    r1587 r1604  
    6363MHFlux::MHFlux(const TH2D &h2d,  const Bool_t Draw,
    6464               const TString varname, const TString unit)
    65  : fHOrig(), fHUnfold(), fHFlux()
     65    : fHOrig(), fHUnfold(), fHFlux()
    6666{
    67   //    if (varname == NULL  ||  unit == NULL)
    68   if (varname == ""  ||  unit == "")
    69     {
    70       *fLog << "MHFlux : varname or unit not defined" << endl;
    71     }
     67    if (varname.IsNull() || unit.IsNull())
     68        *fLog << warn << dbginf << "varname or unit not defined" << endl;
    7269
    7370    fVarname = varname;
     
    9289
    9390    fHOrig.SetDirectory(NULL);
    94     fHOrig.SetXTitle("E-est [GeV]         ");
     91    fHOrig.SetXTitle("E-est [GeV]");
    9592    fHOrig.SetYTitle(strg);
    9693    fHOrig.Sumw2();
     
    109106
    110107    fHUnfold.SetDirectory(NULL);
    111     fHUnfold.SetXTitle("E-unfold [GeV]         ");
     108    fHUnfold.SetXTitle("E-unfold [GeV]");
    112109    fHUnfold.SetYTitle(strg);
    113110    fHUnfold.Sumw2();
     
    127124
    128125    fHFlux.SetDirectory(NULL);
    129     fHFlux.SetXTitle("E-unfold [GeV]         ");
     126    fHFlux.SetXTitle("E-unfold [GeV]");
    130127    fHFlux.SetYTitle(strg);
    131128    fHFlux.Sumw2();
     
    161158      TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");
    162159
    163       char txt0[100];
    164160      strg0  = fVarname;
    165       strg0 += "-bin %d";
    166       sprintf(txt0, strg0, n);
     161      strg0 += "-bin ";
     162      strg0 += n;
    167163
    168164      TString strg1("No.of photons vs. E-est for ");
    169       strg1 += txt0;
    170       new TCanvas(txt0,strg1);
     165      strg1 += strg0;
     166
     167      new TCanvas(strg0, strg1);
    171168      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
    172169      // gROOT->SetSelectedPad(NULL);
     
    174171      gPad->SetLogx();
    175172
    176       h.SetName(txt0);
     173      h.SetName(strg0);
    177174      h.SetTitle(strg1);
    178       h.SetXTitle("E-est [GeV]            ");
     175      h.SetXTitle("E-est [GeV]");
    179176      h.SetYTitle("No.of photons");
    180177      h.DrawCopy();
     
    226223      strg0 += fVarname;
    227224      TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
    228 
    229       char txt0[100];
    230225      strg0  = fVarname;
    231       strg0 += "-bin %d";
    232       sprintf(txt0, strg0, n);
     226      strg0 += "-bin ";
     227      strg0 += n;
    233228
    234229      TString strg1("No.of photons vs. E-unfold for ");
    235       strg1 += txt0;
    236       new TCanvas(txt0,strg1);
     230      strg1 += strg0;
     231
     232      new TCanvas(strg0, strg1);
    237233
    238234      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
     
    241237      gPad->SetLogx();
    242238
    243       h.SetName(txt0);
     239      h.SetName(strg0);
    244240      h.SetTitle(strg1);
    245       h.SetXTitle("E-unfold [GeV]            ");
     241      h.SetXTitle("E-unfold [GeV]");
    246242      h.SetYTitle("No.of photons");
    247243      h.DrawCopy();
     
    276272  //              for the individual bins of the variable Var
    277273
     274    const TAxis &axex = *((TH2*)aeff)->GetXaxis();
     275    const TAxis &axey = *((TH2*)aeff)->GetYaxis();
     276
    278277  //....................................
    279278  // define dummy histogram *aeff
     
    287286  SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
    288287
    289   const Int_t netru    = aeff->GetNbinsX();
    290   const Int_t ntheta   = aeff->GetNbinsY();
     288  const Int_t netru  = aeff->GetNbinsX();
     289  const Int_t ntheta = aeff->GetNbinsY();
    291290
    292291  for (int j=1; j<=netru; j++)
    293292  {
    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     }
     293      for (int k=1; k<=ntheta; k++)
     294      {
     295          Double_t cont = 10000.0;
     296          ((TH1*)aeff)->SetBinContent(j, k, cont);
     297
     298          Double_t dcont = 100.0;
     299          ((TH1*)aeff)->SetBinError(j, k, dcont);
     300      }
    302301  }
    303302  // *fLog << "Dummy aeff : netru =" << netru << ",  ntheta = " << ntheta << endl;
     
    327326  SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
    328327
    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];
     328  Double_t *aeffbar  = new Double_t[nEtru];
     329  Double_t *daeffbar = new Double_t[nEtru];
    338330
    339331  Double_t aeffEunfVar;
    340332  Double_t daeffEunfVar;
    341 
    342333
    343334  //------   start n loop   ------
     
    349340    // determine Theta bins (k1, k2, k3) for interpolation in Theta
    350341    // k0 denotes the Theta bin from whicvh the error is copied
    351     Int_t k0=0, k1=0, k2=0, k3=0;
     342    Int_t k0=0;
     343    Int_t k1=0;
     344    Int_t k2=0;
     345    Int_t k3=0;
     346
    352347    for (int k=3; k<=nTheta; k++)
    353348    {
    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;
     349        Double_t Thetalow = axey.GetBinLowEdge(k);
     350        if (Thetabar < Thetalow)
     351        {
     352            k1 = k-2;
     353            k2 = k-1;
     354            k3 = k;
     355            k0 = k2;
     356            break;
     357        }
     358    }
     359
     360    if (k3 == 0)
     361    {
     362        k1 = nTheta-2;
     363        k2 = nTheta-1;
     364        k3 = nTheta;
    360365        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;
    371366    }
    372367
    373     if (Thetabar <  ((TH1*)aeff)->GetYaxis()->GetBinLowEdge(2))
     368    if (Thetabar <  axey.GetBinLowEdge(2))
    374369      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);   
     370    else
     371        if (Thetabar >  axey.GetBinLowEdge(nTheta))
     372            k0 = nTheta;
     373
     374    Double_t Thetamin = axey.GetBinLowEdge(1);
     375    Double_t Thetamax = axey.GetBinLowEdge(nTheta+1);
    380376    if (Thetabar < Thetamin  ||  Thetabar > Thetamax)
    381377    {
     
    400396    for (int j=1; j<=nEtru; j++)
    401397    {
    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;
     398        double c0 = 0;
     399        double c1 = 0;
     400        double c2 = 0;
     401
     402        const double t1 = cos( axey.GetBinCenter (k1) );
     403        const double t2 = cos( axey.GetBinCenter (k2) );
     404        const double t3 = cos( axey.GetBinCenter (k3) );
     405
     406        const double a1 = aeff->GetBinContent(j, k1);
     407        const double a2 = aeff->GetBinContent(j, k2);
     408        const double a3 = aeff->GetBinContent(j, k3);
     409
     410        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
     411        aeffbar[j]  = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
     412        daeffbar[j] = aeff->GetBinError(j,k0);
     413
     414        //*fLog << "Etru bin " << j <<  ":  tbar= " << Thetabar
     415        //      << ",  abar= "  << aeffbar[j]
     416        //      << ",  dabar= " << daeffbar[j] << endl;
    421417    }
    422418
     
    428424    for (int m=1; m<=nEunf; m++)
    429425    {
    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)
     426        Double_t log10Ebar = 0.5 * ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m)  )+
     427                                     log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
     428        Double_t Ebar = pow(10.0, log10Ebar);
     429
     430        // determine Etru bins (j1, j2, j3) for interpolation in E
     431        // j0 denotes the Etru bin from which the error is copied
     432        Int_t j0=0;
     433        Int_t j1=0;
     434        Int_t j2=0;
     435        Int_t j3=0;
     436
     437        for (int j=3; j<=nEtru; j++)
    443438        {
    444           j1 = j-2;
    445           j2 = j-1;
    446           j3 = j;
    447           j0 = j2;
    448           break;
     439            Double_t Elow = axex.GetBinLowEdge(j);
     440            if (Ebar < Elow)
     441            {
     442                j1 = j-2;
     443                j2 = j-1;
     444                j3 = j;
     445                j0 = j2;
     446                break;
     447            }
    449448        }
    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);
     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 <  axex.GetBinLowEdge(2))
     459            j0 = 1;
     460        else
     461            if (Ebar >  axex.GetBinLowEdge(nEtru))
     462                j0 = nEtru;
     463
     464        Double_t Etrumin = axex.GetBinLowEdge(1);
     465        Double_t Etrumax = axex.GetBinLowEdge(nEtru+1);
     466        if (Ebar < Etrumin  ||  Ebar > Etrumax)
     467        {
     468            *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
     469                << ",  Etrumin =" << Etrumin
     470                << ",  Etrumax =" << Etrumax << endl;
     471        }
     472
     473        //*fLog << "Var bin "   << n  << ":"  <<  endl;
     474        //*fLog << "Ebar= " << Ebar
     475        //      << ",  j1= "    << j1
     476        //      << ",  j2= "    << j2
     477        //      << ",  j3= "    << j3         <<  endl;
     478
     479
     480        double c0=0.0;
     481        double c1=0.0;
     482        double c2=0.0;
     483
     484        const double t1 = 0.5 * ( log10( axex.GetBinLowEdge (j1)  )+
     485                                  log10( axex.GetBinLowEdge (j1+1)) );
     486        const double t2 = 0.5 * ( log10( axex.GetBinLowEdge (j2)  )+
     487                                  log10( axex.GetBinLowEdge (j2+1)) );
     488        const double t3 = 0.5 * ( log10( axex.GetBinLowEdge (j3)  )+
     489                                  log10( axex.GetBinLowEdge (j3+1)) );
     490
     491        const double a1 = aeffbar[j1];
     492        const double a2 = aeffbar[j2];
     493        const double a3 = aeffbar[j3];
     494
     495        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
     496        aeffEunfVar  = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
     497        daeffEunfVar = daeffbar[j0];
     498
     499        //*fLog << "Eunf bin " << m     <<  ":  Ebar= " << Ebar
     500        //      << ",  aeffEunfVar = "  << aeffEunfVar
     501        //      << ",  daeffEunfVar = " << daeffEunfVar << endl;
     502
     503        fHAeff.SetBinContent(m,n,aeffEunfVar);
     504        fHAeff.SetBinError(m,n,daeffEunfVar);
    509505    }
    510506    //---   end m loop ---
     
    518514  for (int m=1; m<=nEunf; m++)
    519515  {
    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)
     516      Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
     517
     518      for (int n=1; n<=nVar; n++)
    535519      {
    536         Cont  = Ngam / (DeltaE * Effon * Aeff);
    537         dCont = Cont * sqrt(   dNgam*dNgam   / (Ngam*Ngam)
    538                              + dEffon*dEffon / (Effon*Effon)
    539                              + dAeff*dAeff   / (Aeff*Aeff)  ); 
     520          Double_t Ngam   = fHUnfold.GetBinContent(m,n);
     521          Double_t dNgam  = fHUnfold.GetBinError(m,n);
     522
     523          Double_t Aeff   = fHAeff.GetBinContent(m,n);
     524          Double_t dAeff  = fHAeff.GetBinError(m,n);
     525
     526          Double_t Effon  = teff->GetBinContent(n);
     527          Double_t dEffon = teff->GetBinError(n);
     528
     529          Double_t Cont, dCont;
     530          if (Ngam>0 && DeltaE>0 && Effon>0 && Aeff>0)
     531          {
     532              Cont  = Ngam / (DeltaE * Effon * Aeff);
     533              dCont = Cont * sqrt( dNgam *dNgam  / (Ngam*Ngam) +
     534                                   dEffon*dEffon / (Effon*Effon) +
     535                                   dAeff *dAeff  / (Aeff*Aeff)  );
     536          }
     537          else
     538          {
     539              Cont  = 1.e-20;
     540              dCont = 1.e-20;
     541          }
     542
     543          fHFlux.SetBinContent(m,n,Cont);
     544          fHFlux.SetBinError(m,n,dCont);
     545
     546          //*fLog << "Eunf bin "    << m      << ",  Var bin " << n
     547          //      << ":  Ngam = "   << Ngam   << ",  Flux = "
     548          //      << Cont  << ", dFlux = " << dCont << endl;
     549          //*fLog << ",  DeltaE = " << DeltaE << ",  Effon = " << Effon
     550          //      << ",  Aeff = "   << Aeff   << endl;
    540551      }
    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     }
    556552  }
    557553
     
    562558  if (Draw == kTRUE)
    563559  {
    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     }
     560      for (int n=1; n<=nVar; n++)
     561      {
     562          TString strg0("Flux-");
     563          strg0 += fVarname;
     564
     565          TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
     566
     567          TString strg1("Photon flux vs. E-unfold for ");
     568          TString strg2 = fVarname;
     569
     570          strg2 += "-bin ";
     571          strg2 += n;
     572
     573          TString strg3 = strg1 + strg2;
     574          new TCanvas(strg2, strg3);
     575          // TCanvas &c = *MakeDefCanvas(txt, txt);
     576          // gROOT->SetSelectedPad(NULL);
     577
     578          gPad->SetLogx();
     579
     580          h.SetName(strg2);
     581          h.SetTitle(strg3);
     582          h.SetXTitle("E-unfold [GeV]            ");
     583          h.SetYTitle("photons / (s m2 GeV)");
     584          h.DrawCopy();
     585
     586          // c.Modified();
     587          // c.Update();
     588      }
    592589  }
    593590  //........................
     
    655652//     such that       yi = F(xi)       for (i=1,3)
    656653//
    657 void 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)
     654Bool_t MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
     655                     Double_t y1, Double_t y2, Double_t y3,
     656                     Double_t *a, Double_t *b, Double_t *c)
    660657{
    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
     658    const double det =
     659        + x2*x3*x3 + x1*x2*x2 + x3*x1*x1
    666660        - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
    667661
    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;
     662    if (det == 0.0)
     663    {
     664        *a = 0;
     665        *b = 0;
     666        *c = 0;
     667        return kFALSE;
     668    }
     669
     670    const double det1 = 1.0/det;
     671
     672    const double ai11 = x2*x3*x3 - x3*x2*x2;
     673    const double ai12 = x3*x1*x1 - x1*x3*x3;
     674    const double ai13 = x1*x2*x2 - x2*x1*x1;
     675
     676    const double ai21 = x2*x2 - x3*x3;
     677    const double ai22 = x3*x3 - x1*x1;
     678    const double ai23 = x1*x1 - x2*x2;
     679
     680    const double ai31 = x3 - x2;
     681    const double ai32 = x1 - x3;
     682    const double ai33 = x2 - x1;
    684683
    685684    *a = (ai11*y1 + ai12*y2 + ai13*y3) * det1;
     
    690689    //yt2 = *a + *b * x2 + *c * x2*x2;
    691690    //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;
     691
     692    //*fLog << "x1 = " << x1 << ",  x2 = " << x2 << ",  x3 = " << x3 << endl;
     693    //*fLog << "y1 = " << y1 << ",  y2 = " << y2 << ",  y3 = " << y3 << endl;
     694    //*fLog << "yt1 = " << yt1 << ",  yt2 = " << yt2
     695    //    << ",  yt3 = " << yt3 << endl;
     696    //*fLog << "*a = " << *a << ",  *b = " << *b << ",  *c= " << *c
     697    //      << ",  *errflag = " << *errflag << endl;
     698
     699    return kTRUE;
    708700}
  • trunk/MagicSoft/Mars/mhist/MHFlux.h

    r1587 r1604  
    5252    const TH2D *GetHFlux()       { return &fHFlux; }
    5353
    54     void Parab(double x1, double x2, double x3,
    55                double y1, double y2, double y3,
    56                double *a, double *b, double *c, int *errflag);
     54    static Bool_t Parab(double x1, double x2, double x3,
     55                        double y1, double y2, double y3,
     56                        double *a, double *b, double *c);
    5757
    5858    ClassDef(MHFlux, 1) //2D-plots (original, unfolded, flux)
Note: See TracChangeset for help on using the changeset viewer.