Changeset 5776 for trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc

Timestamp:

01/10/05 17:55:11 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc (modified) (28 diffs)

trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc

@@ -147 +147 @@
 using namespace std;
 
-class MHillasExt;
+//class MHillasExt;
 
 // --------------------------------------------------------------------------
@@ -155 +155 @@
 MHFalseSource::MHFalseSource(const char *name, const char *title)
     : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5),
-    fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1)
+    fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinDW(-1), fMaxDW(-1),
+    fHistOff(0)
 {
     //
@@ -174 +175 @@
 void MHFalseSource::MakeSymmetric(TH1 *h)
 {
+    h->SetMinimum();
+    h->SetMaximum();
+
     const Float_t min = TMath::Abs(h->GetMinimum());
     const Float_t max = TMath::Abs(h->GetMaximum());
@@ -248 +252 @@
 // Also search for MTime, MObservatory and MPointingPos
 //
-MHillasExt *ext=0;
+//MHillasExt *ext=0;
 Bool_t MHFalseSource::SetupFill(const MParList *plist)
 {
@@ -257 +261 @@
         return kFALSE;
     }
+    /*
     ext = (MHillasExt*)plist->FindObject("MHillasExt");
     if (!ext)
@@ -263 +268 @@
         return kFALSE;
     }
-
+    */
     fMm2Deg = geom->GetConvMm2Deg();
 
-    MBinning binsa;
-    binsa.SetEdges(18, 0, 90);
-
-    const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
-    if (!bins)
+    if (fName!=(TString)"MHFalseSourceOff" && fHistOff==NULL)
     {
-        const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
-
-        MBinning b;
-        b.SetEdges(20, -r, r);
-        SetBinning(&fHist, &b, &b, &binsa);
+        MHFalseSource *hoff = (MHFalseSource*)plist->FindObject("MHFalseSourceOff", "MHFalseSource");
+        if (!hoff)
+            *fLog << inf << "No MHFalseSourceOff [MHFalseSource] found... using current data only!" << endl;
+        else
+        {
+            *fLog << inf << "MHFalseSource [MHFalseSource] found... using on-off mode!" << endl;
+            SetOffData(*hoff);
+        }
     }
+
+    if (fHistOff)
+        MH::SetBinning(&fHist, fHistOff);
     else
-        SetBinning(&fHist, bins, bins, &binsa);
+    {
+        MBinning binsa;
+        binsa.SetEdges(18, 0, 90);
+
+        const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
+        if (!bins)
+        {
+            const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
+
+            MBinning b;
+            b.SetEdges(20, -r, r);
+            SetBinning(&fHist, &b, &b, &binsa);
+        }
+        else
+            SetBinning(&fHist, bins, bins, &binsa);
+    }
 
     fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
@@ -399 +421 @@
 // the same number of bins than for on-data
 //
-void MHFalseSource::ProjectOff(TH2D *h2, TH2D *all)
-{
-    TAxis &axe = *fHist.GetZaxis();
+void MHFalseSource::ProjectOff(const TH3D &src, TH2D *h2, TH2D *all)
+{
+    TAxis &axe = *src.GetZaxis();
 
     // Find range to cut (left edge and width)
@@ -413 +435 @@
 
     // Get projection for range
-    TH2D *p = (TH2D*)fHist.Project3D("yx_off");
+    TH2D *p = (TH2D*)src.Project3D("yx_off");
 
     // Reset range
@@ -435 +457 @@
 // range (0, fAlphaCut)
 //
-void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all)
-{
-    TAxis &axe = *fHist.GetZaxis();
+void MHFalseSource::ProjectOn(const TH3D &src, TH2D *h3, TH2D *all)
+{
+    TAxis &axe = *src.GetZaxis();
 
     // Find range to cut
@@ -445 +467 @@
 
     // Get projection for range
-    TH2D *p = (TH2D*)fHist.Project3D("yx_on");
+    TH2D *p = (TH2D*)src.Project3D("yx_on");
 
     // Reset range
@@ -482 +504 @@
     // Set Minimum as minimum value Greater Than 0
     h3->SetMinimum(GetMinimumGT(*h3));
+}
+
+void MHFalseSource::ProjectOnOff(TH2D *h2, TH2D *h0)
+{
+    ProjectOn(*fHistOff, h2, h0);
+
+    TH2D h;
+    MH::SetBinning(&h, h2);
+
+    // Divide by number of entries in off region (of off-data)
+    ProjectOff(*fHistOff, &h, h0);
+    h2->Divide(&h);
+
+    // Multiply by number of entries in off region (of on-data)
+    ProjectOff(fHist, &h, h0);
+    h2->Multiply(&h);
+
+    // Recalculate Minimum
+    h2->SetMinimum(GetMinimumGT(*h2));
 }
 
@@ -518 +559 @@
     padsave->GetPad(1)->cd(1);
     if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
-        ProjectOn(h3, h0);
+        ProjectOn(fHist, h3, h0);
 
     // Update projection of off-events
     padsave->GetPad(1)->cd(3);
     if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
-        ProjectOff(h2, h0);
+    {
+        if (!fHistOff)
+            ProjectOff(fHist, h2, h0);
+        else
+            ProjectOnOff(h2, h0);
+    }
 
     padsave->GetPad(2)->cd(2);
@@ -579 +625 @@
             TH1 *h = fHist.ProjectionZ("Alpha_z", maxx, maxx, maxy, maxy);
             h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma_{max}=%.1f)", x, y, s));
+
+            TH1D *h0=0;
+            if ((h0 = (TH1D*)gPad->FindObject("AlphaOff_z")))
+            {
+                fHistOff->ProjectionZ("AlphaOff_z", maxx, maxx, maxy, maxy);
+
+                const Int_t f = h0->GetXaxis()->FindFixBin(fBgMean-fAlphaCut/2);
+                const Int_t l = h0->GetXaxis()->FindFixBin(fAlphaCut)+f-1;
+                h0->Scale(h1->Integral(f, l)/h0->Integral(f, l));
+                //h0->Scale(h1->GetEntries()/h0->GetEntries());
+
+            }
         }
     }
@@ -620 +678 @@
     pad->Divide(1, 2, 0, 0.03);
 
-    TObject *catalog = GetCatalog();
+//    TObject *catalog = GetCatalog();
 
     // Initialize upper part
     pad->cd(1);
     // Make sure that the catalog is deleted only once
-    gROOT->GetListOfCleanups()->Add(gPad);
+    // Normally this is not done by root, because it is not necessary...
+    // but in this case it is necessary, because the catalog is
+    // deleted by the first pad and the second one tries to do the same!
+//    gROOT->GetListOfCleanups()->Add(gPad);
     gPad->SetBorderMode(0);
     gPad->Divide(3, 1);
@@ -640 +701 @@
     h3->Draw("colz");
     h3->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
 
     // PAD #2
@@ -654 +715 @@
     h4->Draw("colz");
     h4->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
 
     // PAD #3
     pad->GetPad(1)->cd(3);
     gPad->SetBorderMode(0);
-    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
-    TH1 *h2 = fHist.Project3D("yx_off");
+    TH1 *h2 = 0;
+    if (fHistOff)
+    {
+        fHistOff->GetZaxis()->SetRangeUser(0,fAlphaCut);
+        h2 = fHistOff->Project3D("yx_off");
+        fHistOff->GetZaxis()->SetRange(0,9999);
+    }
+    else
+    {
+        fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
+        h2 = fHist.Project3D("yx_off");
+        fHist.GetZaxis()->SetRange(0,9999);
+    }
     h2->SetDirectory(NULL);
     h2->SetXTitle(fHist.GetXaxis()->GetTitle());
@@ -666 +738 @@
     h2->Draw("colz");
     h2->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
 
     // Initialize lower part
     pad->cd(2);
     // Make sure that the catalog is deleted only once
-    gROOT->GetListOfCleanups()->Add(gPad);
+//    gROOT->GetListOfCleanups()->Add(gPad);
     gPad->SetBorderMode(0);
     gPad->Divide(3, 1);
@@ -685 +757 @@
     h1->Draw();
     h1->SetBit(kCanDelete);
+    if (fHistOff)
+    {
+        h1->SetLineColor(kGreen);
+
+        h1 = fHistOff->ProjectionZ("AlphaOff_z");
+        h1->SetDirectory(NULL);
+        h1->SetTitle("Distribution of \\alpha");
+        h1->SetXTitle(fHistOff->GetZaxis()->GetTitle());
+        h1->SetYTitle("Counts");
+        h1->Draw("same");
+        h1->SetBit(kCanDelete);
+        h1->SetLineColor(kRed);
+    }
 
     // PAD #5
@@ -695 +780 @@
     h5->Draw("colz");
     h5->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
 
     // PAD #6
@@ -706 +791 @@
     h0->Draw("colz");
     h0->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
 }
 
@@ -785 +870 @@
 void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
 {
-    TObject *catalog = GetCatalog();
+//    TObject *catalog = GetCatalog();
 
     TH1D h0a("A",          "", 50,   0, 4000);
@@ -868 +953 @@
     fit.SetPolynomOrder(polynom);
 
-    TH1D *h=0;
+    TH1D *h=0, *hoff=0;
     for (int ix=1; ix<=nx; ix++)
         for (int iy=1; iy<=ny; iy++)
@@ -879 +964 @@
                 continue;
 
-
             h->Scale(entries/h->GetEntries());
 
-            if (!fit.Fit(*h))
+            if (fHistOff)
+            {
+                hoff = fHistOff->ProjectionZ("AlphaFitOff", ix, ix, iy, iy);
+                hoff->Scale(entries/hoff->GetEntries());
+            }
+
+            if (!fit.Fit(*h, hoff))
                 continue;
 
@@ -894 +984 @@
             h1.Fill(fit.GetGausMu());  // mu
             h2.Fill(fit.GetGausSigma());  // sigma-gaus
-            if (polynom>1)
+            if (polynom>1 && !fHistOff)
                 h3.Fill(fit.GetCoefficient(5));
             h4b.Fill(fit.GetChiSqSignal());
@@ -945 +1035 @@
     hists->Draw("colz");
     hists->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
     c->cd(2);
     gPad->SetBorderMode(0);
@@ -976 +1066 @@
 
 
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
     c->cd(3);
     gPad->SetBorderMode(0);
     histb->Draw("colz");
     histb->SetBit(kCanDelete);
-    catalog->Draw("mirror same *");
+//    catalog->Draw("mirror same *");
     c->cd(4);
     gPad->Divide(1,3, 0, 0);
@@ -1018 +1108 @@
                                  hist->GetYaxis()->GetBinCenter(maxy), maxs);
 
-        TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
-        result->Scale(entries/h->GetEntries());
-
-        result->SetDirectory(NULL);
-        result->SetNameTitle("Result \\alpha", title);
-        result->SetBit(kCanDelete);
-        result->SetXTitle("\\alpha [\\circ]");
-        result->SetYTitle("Counts");
-        result->Draw();
-
-        TF1 f1("f1", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
-        TF1 f2("f2", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
+        h = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
+        h->Scale(entries/h->GetEntries());
+
+        h->SetDirectory(NULL);
+        h->SetNameTitle("Result \\alpha", title);
+        h->SetBit(kCanDelete);
+        h->SetXTitle("\\alpha [\\circ]");
+        h->SetYTitle("Counts");
+        h->Draw();
+
+        if (fHistOff)
+        {
+            h->SetLineColor(kGreen);
+
+            TH1D *hof=fHistOff->ProjectionZ("AlphaFitOff", maxx, maxx, maxy, maxy);
+            hof->Scale(entries/hof->GetEntries());
+
+            fit.Scale(*(TH1D*)hof, *(TH1D*)h);
+
+            hof->SetLineColor(kRed);
+            hof->SetDirectory(NULL);
+            hof->SetNameTitle("Result \\alpha", title);
+            hof->SetBit(kCanDelete);
+            hof->SetXTitle("\\alpha [\\circ]");
+            hof->SetYTitle("Counts");
+            hof->Draw("same");
+
+            TH1D *diff = (TH1D*)h->Clone("AlphaFitOnOff");
+            diff->Add(hof, -1);
+            diff->SetLineColor(kBlue);
+            diff->SetNameTitle("Result On-Off \\alpha", title);
+            diff->SetBit(kCanDelete);
+            diff->SetXTitle("\\alpha [\\circ]");
+            diff->SetYTitle("Counts");
+            diff->Draw("same");
+
+            h->SetMinimum(diff->GetMinimum()<0 ? diff->GetMinimum()*1.2 : 0);
+
+            TLine l;
+            l.DrawLine(0, 0, 90, 0);
+        }
+
+        TF1 f1("f1", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90);
+        TF1 f2("f2", Form("gaus(0) + pol%d(3)", fHistOff ? 0 : polynom), 0, 90);
         f1.SetParameters(maxpar.GetArray());
         f2.SetParameters(maxpar.GetArray());
@@ -1053 +1175 @@
         leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
         if (polynom>1)
-            leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
+            leg->AddText(0.60, 0.14, Form("c=%.2f", !fHistOff?maxpar[5]:0))->SetTextAlign(11);
         leg->SetBit(kCanDelete);
         leg->Draw();