Changeset 4836


Ignore:
Timestamp:
09/02/04 17:47:06 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r4829 r4836  
    4040
    4141   * mhist/MHAlpha.cc:
    42      - paint significance
     42     - paint significance and othe rinformations
     43     - unified fit in Finalize and Paint
     44     - replaced significance calculation by Li/Ma
    4345
    4446   * mhvstime/MHVsTime.[h,cc]:
    … …  
    5355   * msignal/MArrivalTime.[h,cc]:
    5456     - added Print()
     57
     58   * manalysis/MEventRateCalc.[h,cc]:
     59     - added the difference in time between two events into the output
     60     - made setup more flexible
     61
     62   * mbase/MContinue.cc:
     63     - fixed a bug which caused a problem if MContinue was not in the
     64       main tasklist
     65
     66   * mimage/MHImagePar.[h,cc], mimage/MHNewImagePar.[h,cc]:
     67     - added Paint function to support logarithmic y-axis scles
    5568
    5669

  • trunk/MagicSoft/Mars/mhist/MHAlpha.cc

    r4828 r4836  
    101101}
    102102
    103 // --------------------------------------------------------------------------
    104 //
    105 // Calculate Significance as
    106 // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
    107 //
    108 Double_t MHAlpha::Significance(Double_t s, Double_t b)
    109 {
    110     const Double_t k = b==0 ? 0 : s/b;
    111     const Double_t f = s+k*k*b;
    112 
    113     return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
    114 }
    115 
    116103Bool_t MHAlpha::SetupFill(const MParList *pl)
    117104{
    … …  
    158145void MHAlpha::Paint(Option_t *opt)
    159146{
     147    DoFit();
     148}
     149
     150// --------------------------------------------------------------------------
     151//
     152// Draw the histogram
     153//
     154void MHAlpha::Draw(Option_t *opt)
     155{
     156    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
     157    pad->SetBorderMode(0);
     158
     159    fHist.Draw();
     160
     161    AppendPad("");
     162}
     163
     164// --------------------------------------------------------------------------
     165//
     166// This is a preliminary implementation of a alpha-fit procedure for
     167// all possible source positions. It will be moved into its own
     168// more powerfull class soon.
     169//
     170// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
     171//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
     172// or
     173//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
     174//
     175// Parameter [1] is fixed to 0 while the alpha peak should be
     176// symmetric around alpha=0.
     177//
     178// Parameter [4] is fixed to 0 because the first derivative at
     179// alpha=0 should be 0, too.
     180//
     181// In a first step the background is fitted between bgmin and bgmax,
     182// while the parameters [0]=0 and [2]=1 are fixed.
     183//
     184// In a second step the signal region (alpha<sigmax) is fittet using
     185// the whole function with parameters [1], [3], [4] and [5] fixed.
     186//
     187// The number of excess and background events are calculated as
     188//   s = int(0, sigint, gaus(0)+pol2(3))
     189//   b = int(0, sigint,         pol2(3))
     190//
     191// The Significance is calculated using the Significance() member
     192// function.
     193//
     194Bool_t MHAlpha::DoFit(Double_t &sig, Bool_t paint)
     195{
    160196    Float_t sigint=fSigInt;
    161197    Float_t sigmax=fSigMax;
    162     Float_t bgmin =fBgMin;
    163     Float_t bgmax =fBgMax;
     198    Float_t bgmin=fBgMin;
     199    Float_t bgmax=fBgMax;
    164200    Byte_t polynom=fPolynom;
    165201
    … …  
    175211
    176212    TH1 *h=&fHist;
    177     //const Double_t alpha0 = h->GetBinContent(1);
    178     //const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);
     213    const Double_t alpha0 = h->GetBinContent(1);
     214    const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);
    179215
    180216    // Check for the regios which is not filled...
    181     const Double_t alpha0 = h->GetBinContent(1);
    182217    if (alpha0==0)
    183         return;
     218        return kFALSE; //*fLog << warn << "Histogram empty." << endl;
    184219
    185220    // First fit a polynom in the off region
    … …  
    191226        func.ReleaseParameter(i);
    192227
    193     //h->Fit(&func, "N0Q", "", bgmin, bgmax);
    194228    h->Fit(&func, "N0Q", "", bgmin, bgmax);
    195     func.SetRange(0, 90);
    196     func.SetLineColor(kRed);
    197     func.Paint("same");
    198 
    199     //h4a.Fill(func.GetChisquare());
    200     //h5a.Fill(func.GetProb());
     229
     230    const Double_t chisq1 = func.GetChisquare()/func.GetNDF();
     231
     232    // ------------------------------------
     233    if (paint)
     234    {
     235        func.SetRange(0, 90);
     236        func.SetLineColor(kRed);
     237        func.Paint("same");
     238    }
     239    // ------------------------------------
    201240
    202241    func.ReleaseParameter(0);
    … …  
    217256    func.SetParameter(2, sigmax*0.75);
    218257
    219     //h->Fit(&func, "N0Q", "", 0, sigmax);
    220258    h->Fit(&func, "N0Q", "", 0, sigmax);
    221     func.SetLineColor(kGreen);
    222     func.Paint("same");
    223 
    224     const Double_t w=fHist.GetXaxis()->GetBinWidth(1);
    225 
    226     const Double_t s   = func.Integral(0, sigint)/w;
     259
     260    const Double_t chisq2    = func.GetChisquare()/func.GetNDF();
     261    const Double_t sigmagaus = func.GetParameter(2);
     262
     263    // ------------------------------------
     264    if (paint)
     265    {
     266        func.SetLineColor(kGreen);
     267        func.Paint("same");
     268    }
     269    // ------------------------------------
     270
     271    const Double_t s   = func.Integral(0, sigint)/alphaw;
    227272    func.SetParameter(0, 0);
    228273    func.SetParameter(2, 1);
    229     const Double_t b   = func.Integral(0, sigint)/w;
    230     const Double_t sig = MMath::SignificanceLiMaSigned(s, b);
    231 
    232     TLatex text;
    233     text.PaintLatex(35, fHist.GetMaximum()*1.1, 0, 0.05, Form("\\sigma=%.1f E=%d", sig, (int)(s-b)));
    234 }
    235 
    236 // --------------------------------------------------------------------------
    237 //
    238 // Draw the histogram
    239 //
    240 void MHAlpha::Draw(Option_t *opt)
    241 {
    242     TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    243     pad->SetBorderMode(0);
    244 
    245     fHist.Draw();
    246 
    247     AppendPad("");
    248 }
    249 
    250 // --------------------------------------------------------------------------
    251 //
    252 // This is a preliminary implementation of a alpha-fit procedure for
    253 // all possible source positions. It will be moved into its own
    254 // more powerfull class soon.
    255 //
    256 // The fit function is "gaus(0)+pol2(3)" which is equivalent to:
    257 //   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
    258 // or
    259 //   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
    260 //
    261 // Parameter [1] is fixed to 0 while the alpha peak should be
    262 // symmetric around alpha=0.
    263 //
    264 // Parameter [4] is fixed to 0 because the first derivative at
    265 // alpha=0 should be 0, too.
    266 //
    267 // In a first step the background is fitted between bgmin and bgmax,
    268 // while the parameters [0]=0 and [2]=1 are fixed.
    269 //
    270 // In a second step the signal region (alpha<sigmax) is fittet using
    271 // the whole function with parameters [1], [3], [4] and [5] fixed.
    272 //
    273 // The number of excess and background events are calculated as
    274 //   s = int(0, sigint, gaus(0)+pol2(3))
    275 //   b = int(0, sigint,         pol2(3))
    276 //
    277 // The Significance is calculated using the Significance() member
    278 // function.
    279 //
     274    const Double_t b   = func.Integral(0, sigint)/alphaw;
     275
     276    sig = MMath::SignificanceLiMaSigned(s, b);
     277
     278    // ------------------------------------
     279    if (paint)
     280    {
     281        TLatex text;
     282        text.PaintLatex(9, fHist.GetMaximum()*1.11, 0, 0.04,
     283                        Form("\\sigma_{Li/Ma}=%.1f  \\omega=%.1f\\circ  E=%d  (\\chi_{b}^{2}/N=%.1f  \\chi_{s}^{2}/N=%.1f)",
     284                             sig, sigmagaus, (int)(s-b), chisq1, chisq2));
     285    }
     286    // ------------------------------------
     287
     288    return kTRUE;
     289}
     290
    280291Bool_t MHAlpha::Finalize()
    281292{
    282     //*fLog << dbg << "Fitting..." << endl;
    283 
    284     Float_t sigint=fSigInt;
    285     Float_t sigmax=fSigMax;
    286     Float_t bgmin=fBgMin;
    287     Float_t bgmax=fBgMax;
    288     Byte_t polynom=fPolynom;
    289 
    290     //                      xmin, xmax, npar
    291     //TF1 func("MyFunc", fcn, 0, 90, 6);
    292     // Implementing the function yourself is only about 5% faster
    293     TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    294     TArrayD maxpar(func.GetNpar());
    295 
    296     //*fLog << "Setup pars..." << endl;
    297 
    298     func.FixParameter(1, 0);
    299     func.FixParameter(4, 0);
    300     func.SetParLimits(2, 0, 90);
    301     func.SetParLimits(3, -1, 1);
    302 
    303     //    TStopwatch clk;
    304     //    clk.Start();
    305    /*
    306     *fLog << inf;
    307     *fLog << "Signal fit:     alpha < " << sigmax << endl;
    308     *fLog << "Integration:    alpha < " << sigint << endl;
    309     *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
    310     *fLog << "Polynom order:  " << (int)polynom << endl;
    311     *fLog << "Fitting False Source Plot..." << flush;
    312     */
    313     TH1 *h=&fHist;
    314     const Double_t alpha0 = h->GetBinContent(1);
    315     const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);
    316 
    317     // Check for the regios which is not filled...
    318     if (alpha0==0)
     293    Double_t sig = 0;
     294
     295    if (!DoFit(sig))
    319296    {
    320297        *fLog << warn << "Histogram empty." << endl;
    321         return kTRUE;
    322     }
    323 
    324     // First fit a polynom in the off region
    325     func.FixParameter(0, 0);
    326     func.FixParameter(2, 1);
    327     func.ReleaseParameter(3);
    328 
    329     //*fLog << "Release pars..." << endl;
    330 
    331     for (int i=5; i<func.GetNpar(); i++)
    332         func.ReleaseParameter(i);
    333 
    334     //*fLog << dbg << "Fit 1" << endl;
    335     //h->Fit(&func, "N0Q", "", bgmin, bgmax);
    336     h->Fit(&func, "N0Q", "", bgmin, bgmax);
    337     //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
    338 
    339     //h4a.Fill(func.GetChisquare());
    340     //h5a.Fill(func.GetProb());
    341 
    342     //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
    343     //func.SetParLimits(2, 5, 90);
    344 
    345     //*fLog << dbg << "Change params..." << endl;
    346 
    347     func.ReleaseParameter(0);
    348     //func.ReleaseParameter(1);
    349     func.ReleaseParameter(2);
    350     func.FixParameter(3, func.GetParameter(3));
    351     for (int i=5; i<func.GetNpar(); i++)
    352         func.FixParameter(i, func.GetParameter(i));
    353 
    354     // Do not allow signals smaller than the background
    355     const Double_t A  = alpha0-func.GetParameter(3);
    356     const Double_t dA = TMath::Abs(A);
    357     func.SetParLimits(0, -dA*4, dA*4);
    358     func.SetParLimits(2, 0, 90);
    359 
    360     // Now fit a gaus in the on region on top of the polynom
    361     func.SetParameter(0, A);
    362     func.SetParameter(2, sigmax*0.75);
    363 
    364     //*fLog << dbg << "Fit 2..." << endl;
    365     //h->Fit(&func, "N0Q", "", 0, sigmax);
    366     h->Fit(&func, "N0Q", "", 0, sigmax);
    367     //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
    368 
    369     //*fLog << dbg << "Copy par..." << endl;
    370     TArrayD p(func.GetNpar(), func.GetParameters());
    371 
    372     // Fill results into some histograms
    373     /*
    374     h0a.Fill(p[0]);
    375     h0b.Fill(p[3]);
    376     h1.Fill(p[1]);
    377     h2.Fill(p[2]);
    378     if (polynom>1)
    379         h3.Fill(p[5]);
    380     h4b.Fill(func.GetChisquare());
    381     //h5b.Fill(func.GetProb());
    382     */
    383     // Implementing the integral as analytical function
    384     // gives the same result in the order of 10e-5
    385     // and it is not faster at all...
    386 
    387     //*fLog << dbg << "Int1 par..." << endl;
    388     const Double_t s = func.Integral(0, sigint)/alphaw;
    389 
    390     func.SetParameter(0, 0);
    391     //func.SetParameter(2, 1);
    392 
    393     //*fLog << dbg << "Int2 par..." << endl;
    394     const Double_t b   = func.Integral(0, sigint)/alphaw;
    395     const Double_t sig = Significance(s, b);
    396 
    397     //*fLog << dbg << "Set Val "<<sig<<" ..." << fResult << endl;
     298        return kFALSE;
     299    }
     300
    398301    if (fResult)
    399302        fResult->SetVal(sig);
    400303
    401     *fLog << inf << "Found Significance: " << sig << "  (width=";
    402     *fLog << func.GetParameter(2) << "°, excess=" << s-b << ")" << endl;
    403 
    404     //*fLog << dbg << "Done." << endl;
    405304    return kTRUE;
    406     /*
    407     if (sig!=0)
    408         h6.Fill(sig);
    409 
    410     if (sig<maxs)
    411     {
    412         maxs = sig;
    413         maxx = ix;
    414         maxy = iy;
    415         maxpar = p;
    416         }
    417        */
    418 /*
    419 
    420     *fLog << inf << "done." << endl;
    421 
    422     h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    423     h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    424 
    425     //hists->SetMinimum(GetMinimumGT(*hists));
    426     histb->SetMinimum(GetMinimumGT(*histb));
    427 
    428 //    clk.Stop();
    429 //    clk.Print();
    430     TCanvas *c=new TCanvas;
    431 
    432     c->Divide(3,2, 0, 0);
    433     c->cd(1);
    434     gPad->SetBorderMode(0);
    435     hists->Draw("colz");
    436     hists->SetBit(kCanDelete);
    437     c->cd(2);
    438     gPad->SetBorderMode(0);
    439     hist->Draw("colz");
    440     hist->SetBit(kCanDelete);
    441     c->cd(3);
    442     gPad->SetBorderMode(0);
    443     histb->Draw("colz");
    444     histb->SetBit(kCanDelete);
    445     c->cd(4);
    446     gPad->Divide(1,3, 0, 0);
    447     TVirtualPad *p=gPad;
    448     p->SetBorderMode(0);
    449     p->cd(1);
    450     gPad->SetBorderMode(0);
    451     h0b.DrawCopy();
    452     h0a.DrawCopy("same");
    453     p->cd(2);
    454     gPad->SetBorderMode(0);
    455     h3.DrawCopy();
    456     p->cd(3);
    457     gPad->SetBorderMode(0);
    458     h2.DrawCopy();
    459     c->cd(6);
    460     gPad->Divide(1,2, 0, 0);
    461     TVirtualPad *q=gPad;
    462     q->SetBorderMode(0);
    463     q->cd(1);
    464     gPad->SetBorderMode(0);
    465     gPad->SetBorderMode(0);
    466     h4b.DrawCopy();
    467     h4a.DrawCopy("same");
    468     h6.DrawCopy("same");
    469     q->cd(2);
    470     gPad->SetBorderMode(0);
    471     //h5b.DrawCopy();
    472     //h5a.DrawCopy("same");
    473     c->cd(5);
    474     gPad->SetBorderMode(0);
    475     if (maxx>0 && maxy>0)
    476     {
    477         const char *title = Form(" \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f) ",
    478                                  fHist.GetXaxis()->GetBinCenter(maxx),
    479                                  fHist.GetYaxis()->GetBinCenter(maxy), maxs);
    480 
    481         TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
    482         result->SetDirectory(NULL);
    483         result->SetNameTitle("Result \\alpha", title);
    484         result->SetBit(kCanDelete);
    485         result->SetXTitle("\\alpha [\\circ]");
    486         result->SetYTitle("Counts");
    487         result->Draw();
    488 
    489         TF1 f1("", func.GetExpFormula(), 0, 90);
    490         TF1 f2("", func.GetExpFormula(), 0, 90);
    491         f1.SetParameters(maxpar.GetArray());
    492         f2.SetParameters(maxpar.GetArray());
    493         f2.FixParameter(0, 0);
    494         f2.FixParameter(1, 0);
    495         f2.FixParameter(2, 1);
    496         f1.SetLineColor(kGreen);
    497         f2.SetLineColor(kRed);
    498 
    499         f2.DrawCopy("same");
    500         f1.DrawCopy("same");
    501 
    502         TPaveText *leg = new TPaveText(0.35, 0.10, 0.90, 0.35, "brNDC");
    503         leg->SetBorderSize(1);
    504         leg->SetTextSize(0.04);
    505         leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
    506         //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
    507         leg->AddLine(0, 0.65, 0, 0.65);
    508         leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
    509         leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
    510         leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
    511         leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
    512         leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
    513         if (polynom>1)
    514             leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
    515         leg->SetBit(kCanDelete);
    516         leg->Draw();
    517 
    518         q->cd(2);
    519 
    520         TGraph *g = new TGraph;
    521         g->SetPoint(0, 0, 0);
    522 
    523         Int_t max=0;
    524         Float_t maxsig=0;
    525         for (int i=1; i<89; i++)
    526         {
    527             const Double_t s = f1.Integral(0, (float)i);
    528             const Double_t b = f2.Integral(0, (float)i);
    529 
    530             const Double_t sig = Significance(s, b);
    531 
    532             g->SetPoint(g->GetN(), i, sig);
    533 
    534             if (sig>maxsig)
    535             {
    536                 max = i;
    537                 maxsig = sig;
    538             }
    539         }
    540 
    541         g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
    542         g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
    543         g->GetHistogram()->SetYTitle("Significance");
    544         g->SetBit(kCanDelete);
    545         g->Draw("AP");
    546 
    547         leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
    548         leg->SetBorderSize(1);
    549         leg->SetTextSize(0.1);
    550         leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max));
    551         leg->SetBit(kCanDelete);
    552         leg->Draw();
    553     }*/
    554305}
    555306

  • trunk/MagicSoft/Mars/mhist/MHAlpha.h

    r3779 r4836  
    3333    Int_t fMap;               //!
    3434
     35    Bool_t DoFit(Double_t &, Bool_t);
     36    Bool_t DoFit(Double_t &d) { return DoFit(d, kFALSE); }
     37    void   DoFit() { Double_t d; DoFit(d, kTRUE); }
     38
    3539public:
    3640    MHAlpha(const char *name=NULL, const char *title=NULL);
    … …  
    4650    void StopMapping();
    4751
    48     static Double_t Significance(Double_t s, Double_t b);
    49 
    5052    ClassDef(MHAlpha, 1) // Alpha-Plot which is fitted online
    5153};
Note: See TracChangeset for help on using the changeset viewer.