Ignore:
Timestamp:
01/10/05 17:55:11 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhflux
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

    r5692 r5776  
    227227        return kFALSE;
    228228
    229     fChiSqSignal = fit.GetChiSqSignal();
    230     fChiSqBg     = fit.GetChiSqBg();
     229    fChiSqSignal  = fit.GetChiSqSignal();
     230    fChiSqBg      = fit.GetChiSqBg();
    231231    fCoefficients = fit.GetCoefficients();
    232232
    … …  
    393393    const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
    394394    const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
     395
    395396    TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E");
    396397    TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E");
    … …  
    486487    if (scaleof!=0)
    487488        of.Scale(scaleon/scaleof);
    488 }
     489    else
     490        of.Reset();
     491}

  • trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h

    r5692 r5776  
    5555public:
    5656    // Implementing the function yourself is only about 5% faster
    57     MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kEntries), fScaleUser(1)
     57    MAlphaFitter(const char *name=0, const char *title=0) : fSigInt(15), fSigMax(75), fBgMin(45), fBgMax(85), fScaleMin(40), fScaleMax(80), fPolynomOrder(2), fCoefficients(3+fPolynomOrder+1), fFunc(new TF1("", Form("gaus(0) + pol%d(3)", fPolynomOrder), 0, 90)), fScaleMode(kOffRegion), fScaleUser(1)
    5858    {
    5959        fName  = name  ? name  : "MAlphaFitter";

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

    r5692 r5776  
    185185    const Int_t n = fHAlpha.GetNbinsY();
    186186
     187    fHEnergy.SetEntries(0);
     188
    187189    for (int i=1; i<=n; i++)
    188190    {
    … …  
    201203
    202204    const Int_t n = fHAlpha.GetNbinsX();
     205
     206    fHTheta.SetEntries(0);
    203207
    204208    for (int i=1; i<=n; i++)
    … …  
    218222    fHTheta.Reset();
    219223
    220     if (fName!=(TString)"MHAlphaOff")
    221     {
    222         MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff");
     224    if (fName!=(TString)"MHAlphaOff" && fOffData==NULL)
     225    {
     226        MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha");
    223227        if (!hoff)
    224228            *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
    … …  
    400404    // Reset Histogram
    401405    fHAlphaTime.Reset();
     406    fHAlphaTime.SetEntries(0);
    402407
    403408    // Get number of bins
    … …  
    646651
    647652        h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
    648         h->SetBit(TH1::kNoTitle|TH1::kNoStats);
     653        h->SetBit(TH1::kNoTitle);
    649654        h->SetXTitle("\\alpha [\\circ]");
    650655        h->SetYTitle("Counts");

  • trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc

    r5301 r5776  
    147147using namespace std;
    148148
    149 class MHillasExt;
     149//class MHillasExt;
    150150
    151151// --------------------------------------------------------------------------
    … …  
    155155MHFalseSource::MHFalseSource(const char *name, const char *title)
    156156    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5),
    157     fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1)
     157    fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1),
     158    fHistOff(0)
    158159{
    159160    //
    … …  
    174175void MHFalseSource::MakeSymmetric(TH1 *h)
    175176{
     177    h->SetMinimum();
     178    h->SetMaximum();
     179
    176180    const Float_t min = TMath::Abs(h->GetMinimum());
    177181    const Float_t max = TMath::Abs(h->GetMaximum());
    … …  
    248252// Also search for MTime, MObservatory and MPointingPos
    249253//
    250 MHillasExt *ext=0;
     254//MHillasExt *ext=0;
    251255Bool_t MHFalseSource::SetupFill(const MParList *plist)
    252256{
    … …  
    257261        return kFALSE;
    258262    }
     263    /*
    259264    ext = (MHillasExt*)plist->FindObject("MHillasExt");
    260265    if (!ext)
    … …  
    263268        return kFALSE;
    264269    }
    265 
     270    */
    266271    fMm2Deg = geom->GetConvMm2Deg();
    267272
    268     MBinning binsa;
    269     binsa.SetEdges(18, 0, 90);
    270 
    271     const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
    272     if (!bins)
     273    if (fName!=(TString)"MHFalseSourceOff" && fHistOff==NULL)
    273274    {
    274         const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
    275 
    276         MBinning b;
    277         b.SetEdges(20, -r, r);
    278         SetBinning(&fHist, &b, &b, &binsa);
     275        MHFalseSource *hoff = (MHFalseSource*)plist->FindObject("MHFalseSourceOff", "MHFalseSource");
     276        if (!hoff)
     277            *fLog << inf << "No MHFalseSourceOff [MHFalseSource] found... using current data only!" << endl;
     278        else
     279        {
     280            *fLog << inf << "MHFalseSource [MHFalseSource] found... using on-off mode!" << endl;
     281            SetOffData(*hoff);
     282        }
    279283    }
     284
     285    if (fHistOff)
     286        MH::SetBinning(&fHist, fHistOff);
    280287    else
    281         SetBinning(&fHist, bins, bins, &binsa);
     288    {
     289        MBinning binsa;
     290        binsa.SetEdges(18, 0, 90);
     291
     292        const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
     293        if (!bins)
     294        {
     295            const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
     296
     297            MBinning b;
     298            b.SetEdges(20, -r, r);
     299            SetBinning(&fHist, &b, &b, &binsa);
     300        }
     301        else
     302            SetBinning(&fHist, bins, bins, &binsa);
     303    }
    282304
    283305    fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
    … …  
    399421// the same number of bins than for on-data
    400422//
    401 void MHFalseSource::ProjectOff(TH2D *h2, TH2D *all)
    402 {
    403     TAxis &axe = *fHist.GetZaxis();
     423void MHFalseSource::ProjectOff(const TH3D &src, TH2D *h2, TH2D *all)
     424{
     425    TAxis &axe = *src.GetZaxis();
    404426
    405427    // Find range to cut (left edge and width)
    … …  
    413435
    414436    // Get projection for range
    415     TH2D *p = (TH2D*)fHist.Project3D("yx_off");
     437    TH2D *p = (TH2D*)src.Project3D("yx_off");
    416438
    417439    // Reset range
    … …  
    435457// range (0, fAlphaCut)
    436458//
    437 void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all)
    438 {
    439     TAxis &axe = *fHist.GetZaxis();
     459void MHFalseSource::ProjectOn(const TH3D &src, TH2D *h3, TH2D *all)
     460{
     461    TAxis &axe = *src.GetZaxis();
    440462
    441463    // Find range to cut
    … …  
    445467
    446468    // Get projection for range
    447     TH2D *p = (TH2D*)fHist.Project3D("yx_on");
     469    TH2D *p = (TH2D*)src.Project3D("yx_on");
    448470
    449471    // Reset range
    … …  
    482504    // Set Minimum as minimum value Greater Than 0
    483505    h3->SetMinimum(GetMinimumGT(*h3));
     506}
     507
     508void MHFalseSource::ProjectOnOff(TH2D *h2, TH2D *h0)
     509{
     510    ProjectOn(*fHistOff, h2, h0);
     511
     512    TH2D h;
     513    MH::SetBinning(&h, h2);
     514
     515    // Divide by number of entries in off region (of off-data)
     516    ProjectOff(*fHistOff, &h, h0);
     517    h2->Divide(&h);
     518
     519    // Multiply by number of entries in off region (of on-data)
     520    ProjectOff(fHist, &h, h0);
     521    h2->Multiply(&h);
     522
     523    // Recalculate Minimum
     524    h2->SetMinimum(GetMinimumGT(*h2));
    484525}
    485526
    … …  
    518559    padsave->GetPad(1)->cd(1);
    519560    if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
    520         ProjectOn(h3, h0);
     561        ProjectOn(fHist, h3, h0);
    521562
    522563    // Update projection of off-events
    523564    padsave->GetPad(1)->cd(3);
    524565    if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
    525         ProjectOff(h2, h0);
     566    {
     567        if (!fHistOff)
     568            ProjectOff(fHist, h2, h0);
     569        else
     570            ProjectOnOff(h2, h0);
     571    }
    526572
    527573    padsave->GetPad(2)->cd(2);
    … …  
    579625            TH1 *h = fHist.ProjectionZ("Alpha_z", maxx, maxx, maxy, maxy);
    580626            h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s));
     627
     628            TH1D *h0=0;
     629            if ((h0 = (TH1D*)gPad->FindObject("AlphaOff_z")))
     630            {
     631                fHistOff->ProjectionZ("AlphaOff_z", maxx, maxx, maxy, maxy);
     632
     633                const Int_t f = h0->GetXaxis()->FindFixBin(fBgMean-fAlphaCut/2);
     634                const Int_t l = h0->GetXaxis()->FindFixBin(fAlphaCut)+f-1;
     635                h0->Scale(h1->Integral(f, l)/h0->Integral(f, l));
     636                //h0->Scale(h1->GetEntries()/h0->GetEntries());
     637
     638            }
    581639        }
    582640    }
    … …  
    620678    pad->Divide(1, 2, 0, 0.03);
    621679
    622     TObject *catalog = GetCatalog();
     680//    TObject *catalog = GetCatalog();
    623681
    624682    // Initialize upper part
    625683    pad->cd(1);
    626684    // Make sure that the catalog is deleted only once
    627     gROOT->GetListOfCleanups()->Add(gPad);
     685    // Normally this is not done by root, because it is not necessary...
     686    // but in this case it is necessary, because the catalog is
     687    // deleted by the first pad and the second one tries to do the same!
     688//    gROOT->GetListOfCleanups()->Add(gPad);
    628689    gPad->SetBorderMode(0);
    629690    gPad->Divide(3, 1);
    … …  
    640701    h3->Draw("colz");
    641702    h3->SetBit(kCanDelete);
    642     catalog->Draw("mirror same *");
     703//    catalog->Draw("mirror same *");
    643704
    644705    // PAD #2
    … …  
    654715    h4->Draw("colz");
    655716    h4->SetBit(kCanDelete);
    656     catalog->Draw("mirror same *");
     717//    catalog->Draw("mirror same *");
    657718
    658719    // PAD #3
    659720    pad->GetPad(1)->cd(3);
    660721    gPad->SetBorderMode(0);
    661     fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
    662     TH1 *h2 = fHist.Project3D("yx_off");
     722    TH1 *h2 = 0;
     723    if (fHistOff)
     724    {
     725        fHistOff->GetZaxis()->SetRangeUser(0,fAlphaCut);
     726        h2 = fHistOff->Project3D("yx_off");
     727        fHistOff->GetZaxis()->SetRange(0,9999);
     728    }
     729    else
     730    {
     731        fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
     732        h2 = fHist.Project3D("yx_off");
     733        fHist.GetZaxis()->SetRange(0,9999);
     734    }
    663735    h2->SetDirectory(NULL);
    664736    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
    … …  
    666738    h2->Draw("colz");
    667739    h2->SetBit(kCanDelete);
    668     catalog->Draw("mirror same *");
     740//    catalog->Draw("mirror same *");
    669741
    670742    // Initialize lower part
    671743    pad->cd(2);
    672744    // Make sure that the catalog is deleted only once
    673     gROOT->GetListOfCleanups()->Add(gPad);
     745//    gROOT->GetListOfCleanups()->Add(gPad);
    674746    gPad->SetBorderMode(0);
    675747    gPad->Divide(3, 1);
    … …  
    685757    h1->Draw();
    686758    h1->SetBit(kCanDelete);
     759    if (fHistOff)
     760    {
     761        h1->SetLineColor(kGreen);
     762
     763        h1 = fHistOff->ProjectionZ("AlphaOff_z");
     764        h1->SetDirectory(NULL);
     765        h1->SetTitle("Distribution of \\alpha");
     766        h1->SetXTitle(fHistOff->GetZaxis()->GetTitle());
     767        h1->SetYTitle("Counts");
     768        h1->Draw("same");
     769        h1->SetBit(kCanDelete);
     770        h1->SetLineColor(kRed);
     771    }
    687772
    688773    // PAD #5
    … …  
    695780    h5->Draw("colz");
    696781    h5->SetBit(kCanDelete);
    697     catalog->Draw("mirror same *");
     782//    catalog->Draw("mirror same *");
    698783
    699784    // PAD #6
    … …  
    706791    h0->Draw("colz");
    707792    h0->SetBit(kCanDelete);
    708     catalog->Draw("mirror same *");
     793//    catalog->Draw("mirror same *");
    709794}
    710795
    … …  
    785870void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
    786871{
    787     TObject *catalog = GetCatalog();
     872//    TObject *catalog = GetCatalog();
    788873
    789874    TH1D h0a("A",          "", 50,   0, 4000);
    … …  
    868953    fit.SetPolynomOrder(polynom);
    869954
    870     TH1D *h=0;
     955    TH1D *h=0, *hoff=0;
    871956    for (int ix=1; ix<=nx; ix++)
    872957        for (int iy=1; iy<=ny; iy++)
    … …  
    879964                continue;
    880965
    881 
    882966            h->Scale(entries/h->GetEntries());
    883967
    884             if (!fit.Fit(*h))
     968            if (fHistOff)
     969            {
     970                hoff = fHistOff->ProjectionZ("AlphaFitOff", ix, ix, iy, iy);
     971                hoff->Scale(entries/hoff->GetEntries());
     972            }
     973
     974            if (!fit.Fit(*h, hoff))
    885975                continue;
    886976
    … …  
    894984            h1.Fill(fit.GetGausMu());  // mu
    895985            h2.Fill(fit.GetGausSigma());  // sigma-gaus
    896             if (polynom>1)
     986            if (polynom>1 && !fHistOff)
    897987                h3.Fill(fit.GetCoefficient(5));
    898988            h4b.Fill(fit.GetChiSqSignal());
    … …  
    9451035    hists->Draw("colz");
    9461036    hists->SetBit(kCanDelete);
    947     catalog->Draw("mirror same *");
     1037//    catalog->Draw("mirror same *");
    9481038    c->cd(2);
    9491039    gPad->SetBorderMode(0);
    … …  
    9761066
    9771067
    978     catalog->Draw("mirror same *");
     1068//    catalog->Draw("mirror same *");
    9791069    c->cd(3);
    9801070    gPad->SetBorderMode(0);
    9811071    histb->Draw("colz");
    9821072    histb->SetBit(kCanDelete);
    983     catalog->Draw("mirror same *");
     1073//    catalog->Draw("mirror same *");
    9841074    c->cd(4);
    9851075    gPad->Divide(1,3, 0, 0);
    … …  
    10181108                                 hist->GetYaxis()->GetBinCenter(maxy), maxs);
    10191109
    1020         TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
    1021         result->Scale(entries/h->GetEntries());
    1022 
    1023         result->SetDirectory(NULL);
    1024         result->SetNameTitle("Result \\alpha", title);
    1025         result->SetBit(kCanDelete);
    1026         result->SetXTitle("\\alpha [\\circ]");
    1027         result->SetYTitle("Counts");
    1028         result->Draw();
    1029 
    1030         TF1 f1("f1", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
    1031         TF1 f2("f2", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
     1110        h = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
     1111        h->Scale(entries/h->GetEntries());
     1112
     1113        h->SetDirectory(NULL);
     1114        h->SetNameTitle("Result \\alpha", title);
     1115        h->SetBit(kCanDelete);
     1116        h->SetXTitle("\\alpha [\\circ]");
     1117        h->SetYTitle("Counts");
     1118        h->Draw();
     1119
     1120        if (fHistOff)
     1121        {
     1122            h->SetLineColor(kGreen);
     1123
     1124            TH1D *hof=fHistOff->ProjectionZ("AlphaFitOff", maxx, maxx, maxy, maxy);
     1125            hof->Scale(entries/hof->GetEntries());
     1126
     1127            fit.Scale(*(TH1D*)hof, *(TH1D*)h);
     1128
     1129            hof->SetLineColor(kRed);
     1130            hof->SetDirectory(NULL);
     1131            hof->SetNameTitle("Result \\alpha", title);
     1132            hof->SetBit(kCanDelete);
     1133            hof->SetXTitle("\\alpha [\\circ]");
     1134            hof->SetYTitle("Counts");
     1135            hof->Draw("same");
     1136
     1137            TH1D *diff = (TH1D*)h->Clone("AlphaFitOnOff");
     1138            diff->Add(hof, -1);
     1139            diff->SetLineColor(kBlue);
     1140            diff->SetNameTitle("Result On-Off \\alpha", title);
     1141            diff->SetBit(kCanDelete);
     1142            diff->SetXTitle("\\alpha [\\circ]");
     1143            diff->SetYTitle("Counts");
     1144            diff->Draw("same");
     1145
     1146            h->SetMinimum(diff->GetMinimum()<0 ? diff->GetMinimum()*1.2 : 0);
     1147
     1148            TLine l;
     1149            l.DrawLine(0, 0, 90, 0);
     1150        }
     1151
     1152        TF1 f1("f1", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90);
     1153        TF1 f2("f2", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90);
    10321154        f1.SetParameters(maxpar.GetArray());
    10331155        f2.SetParameters(maxpar.GetArray());
    … …  
    10531175        leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
    10541176        if (polynom>1)
    1055             leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
     1177            leg->AddText(0.60, 0.14, Form("c=%.2f", !fHistOff?maxpar[5]:0))->SetTextAlign(11);
    10561178        leg->SetBit(kCanDelete);
    10571179        leg->Draw();

  • trunk/MagicSoft/Mars/mhflux/MHFalseSource.h

    r5301 r5776  
    3939    TH3D    fHist;               // Alpha vs. x and y
    4040
     41    const TH3D *fHistOff;
     42
    4143    Double_t fRa;
    4244    Double_t fDec;
    … …  
    4648
    4749    void ProjectAll(TH2D *h);
    48     void ProjectOff(TH2D *h, TH2D *all);
    49     void ProjectOn(TH2D *h, TH2D *all);
     50    void ProjectOff(const TH3D &src, TH2D *h, TH2D *all);
     51    void ProjectOn(const TH3D &src, TH2D *h, TH2D *all);
     52    void ProjectOnOff(TH2D *h, TH2D *all);
    5053
    5154    TObject *GetCatalog();
    … …  
    7780    void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU*
    7881
     82    void SetOffData(const MHFalseSource &fs) {
     83        fHistOff  = &fs.fHist;
     84        fMinDist  =  fs.fMinDist;
     85        fMaxDist  =  fs.fMaxDist;
     86        fMaxDW    =  fs.fMaxDW;
     87        fMinDW    =  fs.fMinDW;
     88        fAlphaCut =  fs.fAlphaCut;
     89        fBgMean   =  fs.fBgMean;
     90    }
     91
    7992    void Paint(Option_t *opt="");
    8093    void Draw(Option_t *option="");
Note: See TracChangeset for help on using the changeset viewer.