Changeset 4975


Ignore:
Timestamp:
09/13/04 13:28:58 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhist
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r4966 r4975  
    4848//
    4949//  For each "time-bin" the histogram is fitted and the resulting effective
    50 //  on-time is stored in the fHEffOnTime histogram. Each entry in this
     50//  on-time is stored in the fHTimeEffOn histogram. Each entry in this
    5151//  histogram is the effective observation time between the upper and
    5252//  lower edges of the bins.
    … …  
    8383
    8484#include <TLatex.h>
     85#include <TGaxis.h>
    8586#include <TCanvas.h>
    8687#include <TPaveStats.h>
    … …  
    131132
    132133    // effective on time versus theta
    133     fHEffOnTheta.SetName("EffOnTheta");
    134     fHEffOnTheta.SetTitle("Effective On Time T_{eff}");
    135     fHEffOnTheta.SetXTitle("\\Theta [\\circ]");
    136     fHEffOnTheta.SetYTitle("T_{eff} [s]");
    137     fHEffOnTheta.UseCurrentStyle();
    138     fHEffOnTheta.SetDirectory(NULL);
    139     fHEffOnTheta.GetYaxis()->SetTitleOffset(1.2);
     134    fHThetaEffOn.SetName("EffOnTheta");
     135    fHThetaEffOn.SetTitle("Effective On Time T_{eff}");
     136    fHThetaEffOn.SetXTitle("\\Theta [\\circ]");
     137    fHThetaEffOn.SetYTitle("T_{eff} [s]");
     138    fHThetaEffOn.UseCurrentStyle();
     139    fHThetaEffOn.SetDirectory(NULL);
     140    fHThetaEffOn.GetYaxis()->SetTitleOffset(1.2);
     141    fHThetaEffOn.GetYaxis()->SetTitleColor(kBlue);
     142    fHThetaEffOn.SetLineColor(kBlue);
    140143    //fHEffOn.Sumw2();
    141144
    142145    // effective on time versus time
    143     fHEffOnTime.SetName("EffOnTime");
    144     fHEffOnTime.SetTitle("Effective On Time T_{eff}");
    145     fHEffOnTime.SetXTitle("Time");
    146     fHEffOnTime.SetYTitle("T_{eff} [s]");
    147     fHEffOnTime.UseCurrentStyle();
    148     fHEffOnTime.SetDirectory(NULL);
    149     fHEffOnTime.GetYaxis()->SetTitleOffset(1.2);
    150     fHEffOnTime.GetXaxis()->SetLabelSize(0.033);
    151     fHEffOnTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
    152     fHEffOnTime.GetXaxis()->SetTimeDisplay(1);
    153     fHEffOnTime.Sumw2();
     146    fHTimeEffOn.SetName("EffOnTime");
     147    fHTimeEffOn.SetTitle("Effective On Time T_{eff}");
     148    fHTimeEffOn.SetXTitle("Time");
     149    fHTimeEffOn.SetYTitle("T_{eff} [s]");
     150    fHTimeEffOn.UseCurrentStyle();
     151    fHTimeEffOn.SetDirectory(NULL);
     152    fHTimeEffOn.GetYaxis()->SetTitleOffset(1.2);
     153    fHTimeEffOn.GetXaxis()->SetLabelSize(0.033);
     154    fHTimeEffOn.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
     155    fHTimeEffOn.GetXaxis()->SetTimeDisplay(1);
     156    fHTimeEffOn.GetYaxis()->SetTitleColor(kBlue);
     157    fHTimeEffOn.SetLineColor(kBlue);
    154158
    155159    // chi2 probability
    156     fHProbTheta.SetName("ProbTheta");
    157     fHProbTheta.SetTitle("\\chi^{2} Probability of Fit");
    158     fHProbTheta.SetXTitle("\\Theta [\\circ]");
    159     fHProbTheta.SetYTitle("p [%]");
    160     fHProbTheta.UseCurrentStyle();
    161     fHProbTheta.SetDirectory(NULL);
    162     fHProbTheta.GetYaxis()->SetTitleOffset(1.2);
    163     fHProbTheta.SetMaximum(101);
     160    fHThetaProb.SetName("ProbTheta");
     161    fHThetaProb.SetTitle("\\chi^{2} Probability of Fit");
     162    fHThetaProb.SetXTitle("\\Theta [\\circ]");
     163    fHThetaProb.SetYTitle("p [%]");
     164    fHThetaProb.UseCurrentStyle();
     165    fHThetaProb.SetDirectory(NULL);
     166    fHThetaProb.GetYaxis()->SetTitleOffset(1.2);
     167    fHThetaProb.SetMaximum(101);
     168    fHThetaProb.GetYaxis()->SetTitleColor(kBlue);
     169    fHThetaProb.SetLineColor(kBlue);
    164170
    165171    // chi2 probability
    166     fHProbTime.SetName("ProbTime");
    167     fHProbTime.SetTitle("\\chi^{2} Probability of Fit");
    168     fHProbTime.SetXTitle("Time");
    169     fHProbTime.SetYTitle("p [%]");
    170     fHProbTime.UseCurrentStyle();
    171     fHProbTime.SetDirectory(NULL);
    172     fHProbTime.GetYaxis()->SetTitleOffset(1.2);
    173     fHProbTime.GetXaxis()->SetLabelSize(0.033);
    174     fHProbTime.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
    175     fHProbTime.GetXaxis()->SetTimeDisplay(1);
    176     fHProbTime.SetMaximum(101);
     172    fHTimeProb.SetName("ProbTime");
     173    fHTimeProb.SetTitle("\\chi^{2} Probability of Fit");
     174    fHTimeProb.SetXTitle("Time");
     175    fHTimeProb.SetYTitle("p [%]");
     176    fHTimeProb.UseCurrentStyle();
     177    fHTimeProb.SetDirectory(NULL);
     178    fHTimeProb.GetYaxis()->SetTitleOffset(1.2);
     179    fHTimeProb.GetXaxis()->SetLabelSize(0.033);
     180    fHTimeProb.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
     181    fHTimeProb.GetXaxis()->SetTimeDisplay(1);
     182    fHTimeProb.SetMaximum(101);
     183    fHTimeProb.GetYaxis()->SetTitleColor(kBlue);
     184    fHTimeProb.SetLineColor(kBlue);
    177185
    178186    // lambda versus theta
    179     fHLambda.SetName("lambda");
    180     fHLambda.SetTitle("lambda of Effective On Time Fit");
    181     fHLambda.SetXTitle("\\Theta [\\circ]");
    182     fHLambda.SetYTitle("\\lambda [Hz]");
    183     fHLambda.UseCurrentStyle();
    184     fHLambda.SetDirectory(NULL);
    185     // fHLambda.Sumw2();
    186 
    187     // N0 versus theta
    188     fHN0.SetName("N0");
    189     fHN0.SetTitle("Ideal number of events N_{0}");
    190     fHN0.SetXTitle("\\Theta [\\circ]");
    191     fHN0.SetYTitle("N_{0}");
    192     fHN0.UseCurrentStyle();
    193     fHN0.SetDirectory(NULL);
    194     //fHN0del.Sumw2();
    195 
    196     // chi2/NDF versus theta
     187    fHThetaLambda.SetName("LambdaTheta");
     188    fHThetaLambda.SetTitle("Slope (Rate) vs Theta");
     189    fHThetaLambda.SetXTitle("\\Theta [\\circ]");
     190    fHThetaLambda.SetYTitle("S");
     191    fHThetaLambda.UseCurrentStyle();
     192    fHThetaLambda.SetDirectory(NULL);
     193    fHThetaLambda.SetLineColor(kGreen);
     194
     195    // lambda versus time
     196    fHTimeLambda.SetName("LambdaTime");
     197    fHTimeLambda.SetTitle("Slope (Rate) vs Time");
     198    fHTimeLambda.SetXTitle("\\Time [\\circ]");
     199    fHTimeLambda.SetYTitle("S");
     200    fHTimeLambda.UseCurrentStyle();
     201    fHTimeLambda.SetDirectory(NULL);
     202    fHTimeLambda.GetYaxis()->SetTitleOffset(1.2);
     203    fHTimeLambda.GetXaxis()->SetLabelSize(0.033);
     204    fHTimeLambda.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
     205    fHTimeLambda.GetXaxis()->SetTimeDisplay(1);
     206    fHTimeLambda.SetLineColor(kGreen);
     207
     208    // NDoF versus theta
     209    fHThetaNDF.SetName("NDofTheta");
     210    fHThetaNDF.SetTitle("Number of Degrees of freedom vs Theta");
     211    fHThetaNDF.SetXTitle("\\Theta [\\circ]");
     212    fHThetaNDF.SetYTitle("NDoF [#]");
     213    fHThetaNDF.UseCurrentStyle();
     214    fHThetaNDF.SetDirectory(NULL);
     215    fHThetaNDF.SetLineColor(kGreen);
     216
     217    // NDoF versus time
    197218    /*
    198      fHChi2.SetName("Chi2/NDF");
    199      fHChi2.SetTitle("\\chi^{2}/NDF of Effective On Time Fit");
    200      fHChi2.SetXTitle("\\Theta [\\circ]");
    201      fHChi2.SetYTitle("\\chi^{2}/NDF");
    202      fHChi2.UseCurrentStyle();
    203      fHChi2.SetDirectory(NULL);
    204      fHChi2.GetYaxis()->SetTitleOffset(1.2);
    205      //fHChi2.Sumw2();
    206      */
    207 
     219    fHTimeNDF.SetName("NDofTime");
     220    fHTimeNDF.SetTitle("Number of Degrees of freedom vs Time");
     221    fHTimeNDF.SetXTitle("Time");
     222    fHTimeNDF.SetYTitle("NDoF [#]");
     223    fHTimeNDF.UseCurrentStyle();
     224    fHTimeNDF.SetDirectory(NULL);
     225    fHTimeNDF.GetYaxis()->SetTitleOffset(1.2);
     226    fHTimeNDF.GetXaxis()->SetLabelSize(0.033);
     227    fHTimeNDF.GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
     228    fHTimeNDF.GetXaxis()->SetTimeDisplay(1);
     229    fHTimeNDF.SetLineColor(kBlue);
     230    */
    208231    // setup binning
    209232    MBinning btheta("BinningTheta");
    … …  
    217240    btime.Apply(fH1DeltaT);
    218241
    219     btheta.Apply(fHEffOnTheta);
    220     btheta.Apply(fHLambda);
    221     btheta.Apply(fHN0);
    222     btheta.Apply(fHProbTheta);
     242    btheta.Apply(fHThetaEffOn);
     243    btheta.Apply(fHThetaLambda);
     244    btheta.Apply(fHThetaNDF);
     245    btheta.Apply(fHThetaProb);
    223246    //btheta.Apply(fHChi2);
    224247}
    … …  
    251274   if (binstheta)
    252275   {
    253        binstheta->Apply(fHEffOnTheta);
    254        binstheta->Apply(fHLambda);
    255        binstheta->Apply(fHN0);
    256        binstheta->Apply(fHProbTheta);
     276       binstheta->Apply(fHThetaEffOn);
     277       binstheta->Apply(fHThetaLambda);
     278       binstheta->Apply(fHThetaNDF);
     279       binstheta->Apply(fHThetaProb);
    257280       //binstheta->Apply(fHChi2);
    258281   }
    … …  
    295318    // parameter 1 = N0*del
    296319    //
    297     TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
    298     func.SetParNames("lambda", "N0", "del");
     320    static TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)");
     321    //func.SetParNames("lambda", "N0", "del");
    299322
    300323    func.SetParameter(0, 100);       // Hz
    … …  
    302325    func.FixParameter(2, Nmdel/Nm);
    303326
    304     // options : 0  do not plot the function
     327    // options : N  do not store the function, do not draw
    305328    //           I  use integral of function in bin rather than value at bin center
    306329    //           R  use the range specified in the function range
    307330    //           Q  quiet mode
    308     h->Fit(&func, "0IRQ");
     331    h->Fit(&func, "NIQ", "", yq[0], yq[1]);
    309332
    310333    const Double_t chi2 = func.GetChisquare();
    … …  
    325348
    326349    const Double_t lambda = func.GetParameter(0);
    327     const Double_t N0     = func.GetParameter(1);
     350    //const Double_t N0     = func.GetParameter(1);
    328351    const Double_t prob   = func.GetProb();
    329352
    … …  
    362385    res[4] = dl;
    363386
    364     // N0 of fit
    365     res[5] = N0;
     387    // NDoF of fit
     388    res[5] = NDF;
    366389
    367390    // Rdead (from fit) is the fraction from real time lost by the dead time
    … …  
    379402void MHEffectiveOnTime::FitThetaBins()
    380403{
    381     fHEffOnTheta.Reset();
    382     fHProbTheta.Reset();
    383     fHLambda.Reset();
    384     fHN0.Reset();
     404    fHThetaEffOn.Reset();
     405    fHThetaProb.Reset();
     406    fHThetaLambda.Reset();
     407    fHThetaNDF.Reset();
    385408
    386409    const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
    … …  
    400423
    401424        // the effective on time is Nm/lambda
    402         fHEffOnTheta.SetBinContent(i, res[0]);
    403         fHEffOnTheta.SetBinError  (i, res[1]);
     425        fHThetaEffOn.SetBinContent(i, res[0]);
     426        fHThetaEffOn.SetBinError  (i, res[1]);
    404427
    405428        // plot chi2-probability of fit
    406         fHProbTheta.SetBinContent(i, res[2]);
     429        fHThetaProb.SetBinContent(i, res[2]);
    407430
    408431        // plot chi2/NDF of fit
    … …  
    410433
    411434        // lambda of fit
    412         fHLambda.SetBinContent(i, res[3]);
    413         fHLambda.SetBinError  (i, res[4]);
    414 
    415         // N0 of fit
    416         fHN0.SetBinContent(i, res[5]);
     435        fHThetaLambda.SetBinContent(i, res[3]);
     436        fHThetaLambda.SetBinError  (i, res[4]);
     437
     438        // NDoF of fit
     439        fHThetaNDF.SetBinContent(i, res[5]);
    417440
    418441        // Rdead (from fit) is the fraction from real time lost by the dead time
    … …  
    449472
    450473    // Get number of bins
    451     const Int_t n = fHEffOnTime.GetNbinsX();
     474    const Int_t n = fHTimeEffOn.GetNbinsX();
    452475
    453476    // Enhance binning
    454477    MBinning bins;
    455     bins.SetEdges(fHEffOnTime, 'x');
     478    bins.SetEdges(fHTimeEffOn, 'x');
    456479    bins.AddEdge(fLastTime.GetAxisTime());
    457     bins.Apply(fHEffOnTime);
    458     bins.Apply(fHProbTime);
     480    bins.Apply(fHTimeEffOn);
     481    bins.Apply(fHTimeProb);
     482    bins.Apply(fHTimeLambda);
     483    //bins.Apply(fHTimeNDF);
    459484
    460485    //
    461486    // Fill histogram
    462487    //
    463     fHEffOnTime.SetBinContent(n, res[0]);
    464     fHEffOnTime.SetBinError(n, res[1]);
    465 
    466     fHProbTime.SetBinContent(n, res[2]);
     488    fHTimeEffOn.SetBinContent(n, res[0]);
     489    fHTimeEffOn.SetBinError(n, res[1]);
     490
     491    fHTimeProb.SetBinContent(n, res[2]);
     492
     493    fHTimeLambda.SetBinContent(n, res[3]);
     494    fHTimeLambda.SetBinError(n, res[4]);
     495
     496    //fHTimeNDF.SetBinContent(n, res[5]);
    467497
    468498    //
    … …  
    500530        MBinning bins;
    501531        bins.SetEdges(1, time->GetAxisTime()-fNumEvents/200, time->GetAxisTime());
    502         bins.Apply(fHEffOnTime);
    503         bins.Apply(fHProbTime);
     532        bins.Apply(fHTimeEffOn);
     533        bins.Apply(fHTimeProb);
     534        bins.Apply(fHTimeLambda);
     535        //bins.Apply(fHTimeNDF);
    504536
    505537        fParam->SetVal(0, 0);
    … …  
    540572    FitThetaBins();
    541573    FitTimeBin();
    542     MH::RemoveFirstBin(fHEffOnTime);
    543     MH::RemoveFirstBin(fHProbTime);
     574    MH::RemoveFirstBin(fHTimeEffOn);
     575    MH::RemoveFirstBin(fHTimeProb);
     576    MH::RemoveFirstBin(fHTimeLambda);
     577    //MH::RemoveFirstBin(fHTimeNDF);
    544578
    545579    fIsFinalized = kTRUE;
    … …  
    558592    text.SetTextSize(0.04);
    559593    text.Paint();
     594}
     595
     596void MHEffectiveOnTime::PaintProb(TH1 &h) const
     597{
     598    Double_t sum = 0;
     599    Int_t    n = 0;
     600    for (int i=0; i<h.GetNbinsX(); i++)
     601        if (h.GetBinContent(i+1)>0)
     602        {
     603            sum += h.GetBinContent(i+1);
     604            n++;
     605        }
     606
     607    if (n==0)
     608        return;
     609
     610    TLatex text(0.45, 0.94, Form("\\bar{p} = %.1f%%  (n=%d)", sum/n, n));
     611    text.SetBit(TLatex::kTextNDC);
     612    text.SetTextSize(0.04);
     613    text.Paint();
     614}
     615
     616void MHEffectiveOnTime::UpdateRightAxis(TH1 &h)
     617{
     618    cout << "Update... " << flush;
     619    const Double_t max = h.GetMaximum()*1.1;
     620    if (max==0)
     621        return;
     622
     623    h.SetNormFactor(h.Integral()*gPad->GetUymax()/max);
     624
     625    TGaxis *axis = (TGaxis*)gPad->FindObject("RightAxis");
     626    if (!axis)
     627        return;
     628
     629    cout << axis << " ";
     630    cout << "X: " << gPad->GetUxmax() << "  ";
     631    cout << "Y: " << gPad->GetUymin() << "/" << gPad->GetUymax() << "  ";
     632    cout << "max=" << max << endl;
     633
     634    axis->SetX1(gPad->GetUxmax());
     635    axis->SetX2(gPad->GetUxmax());
     636    axis->SetY1(gPad->GetUymin());
     637    axis->SetY2(gPad->GetUymax());
     638    axis->SetWmax(max);
    560639}
    561640
    … …  
    621700    }
    622701
     702    if (o==(TString)"timendf")
     703    {
     704        //    UpdateRightAxis(fHTimeNDF);
     705        // FIXME: first bin?
     706        PaintProb(fHTimeProb);
     707    }
     708
     709    if (o==(TString)"thetandf")
     710    {
     711        UpdateRightAxis(fHThetaNDF);
     712        // FIXME: first bin?
     713        PaintProb(fHThetaProb);
     714    }
     715
    623716    h=0;
    624717    if (o==(TString)"theta")
    625         h = &fHEffOnTheta;
     718    {
     719        h = &fHThetaEffOn;
     720        UpdateRightAxis(fHThetaLambda);
     721    }
    626722    if (o==(TString)"time")
    627         h = &fHEffOnTime;
     723    {
     724        h = &fHTimeEffOn;
     725        UpdateRightAxis(fHTimeLambda);
     726    }
    628727
    629728    if (!h)
    … …  
    635734
    636735    PaintText(h->Integral(), error);
     736}
     737
     738void MHEffectiveOnTime::DrawRightAxis(const char *title)
     739{
     740    TGaxis *axis = new TGaxis(gPad->GetUxmax(), gPad->GetUymin(),
     741                              gPad->GetUxmax(), gPad->GetUymax(),
     742                              0, 1, 510, "+L");
     743    axis->SetName("RightAxis");
     744    axis->SetTitle(title);
     745    axis->SetTitleOffset(0.9);
     746    axis->SetTextColor(kGreen);
     747    axis->CenterTitle();
     748    axis->SetBit(kCanDelete);
     749    axis->Draw();
    637750}
    638751
    … …  
    669782    pad->GetPad(1)->cd(2);
    670783    gPad->SetBorderMode(0);
    671     fHProbTime.Draw();
    672     AppendPad("paint");
     784    fHTimeProb.Draw();
     785    AppendPad("timendf");
     786    //fHTimeNDF.Draw("same");
     787    //DrawRightAxis("NDF");
    673788
    674789    pad->GetPad(1)->cd(3);
    675790    gPad->SetBorderMode(0);
    676     fHEffOnTime.Draw();
     791    fHTimeEffOn.Draw();
    677792    AppendPad("time");
     793    fHTimeLambda.Draw("same");
     794    DrawRightAxis("Slope");
    678795
    679796    pad->cd(2);
    … …  
    695812    pad->GetPad(2)->cd(2);
    696813    gPad->SetBorderMode(0);
    697     fHProbTheta.Draw();
     814    fHThetaProb.Draw();
     815    AppendPad("thetandf");
     816    fHThetaNDF.Draw("same");
     817    DrawRightAxis("NDF");
    698818
    699819    pad->GetPad(2)->cd(3);
    700820    gPad->SetBorderMode(0);
    701     fHEffOnTheta.Draw();
     821    fHThetaEffOn.Draw();
    702822    AppendPad("theta");
    703 }
     823    fHThetaLambda.Draw("same");
     824    DrawRightAxis("Slope");
     825}

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

    r4925 r4975  
    1313#ifndef ROOT_TH2
    1414#include <TH2.h>
     15#endif
     16#ifndef ROOT_TF1
     17#include <TF1.h>
    1518#endif
    1619
    … …  
    3336    TH1D fH1DeltaT;      //! Counts vs Delta T (for a time interval)
    3437
    35     TH1D fHEffOnTheta;   // Effective On time versus Theta
    36     TH1D fHEffOnTime;    // Effective On time versus Time
     38    TH1D fHThetaEffOn;   // Effective On time versus Theta
     39    TH1D fHThetaProb;    // Chisq prob fit of Effective On time versus Theta
     40    TH1D fHThetaNDF;     // NDF vs Theta
     41    TH1D fHThetaLambda;  // Slope (rate) vs Theta
    3742
    38     TH1D fHProbTheta;    // Chisq prob fit of Effective On time versus Theta
    39     TH1D fHProbTime;     // Chisq prob fit of Effective On time versus Time
     43    TH1D fHTimeEffOn;    // Effective On time versus Time
     44    TH1D fHTimeProb;     // Chisq prob fit of Effective On time versus Time
     45    //TH1D fHTimeNDF;      // NDF vs Time
     46    TH1D fHTimeLambda;   // Slope (rate) vs Time
    4047
    41     TH1D fHN0;
    42     TH1D fHLambda;
    43 
    44     Bool_t fIsFinalized; // Flag telling you whether fHEffOnTheta is the final result
     48    Bool_t fIsFinalized; // Flag telling you whether fHThetaEffOn is the final result
    4549
    4650    Int_t fNumEvents;    // Number of events to be used for a bin in time
    … …  
    5256    void FitThetaBins();
    5357    void FitTimeBin();
     58    void PaintProb(TH1 &h) const;
    5459    void PaintText(Double_t val, Double_t error) const;
     60    void DrawRightAxis(const char *title);
     61    void UpdateRightAxis(TH1 &h);
    5562
    5663public:
    … …  
    6370    Bool_t Finalize();
    6471
    65     const TH1D &GetHEffOnTheta() const { return fHEffOnTheta; }
    66     const TH1D &GetHEffOnTime() const { return fHEffOnTime; }
     72    const TH1D &GetHEffOnTheta() const { return fHThetaEffOn; }
     73    const TH1D &GetHEffOnTime() const { return fHTimeEffOn; }
    6774
    6875    void Draw(Option_t *option="");
Note: See TracChangeset for help on using the changeset viewer.