Changeset 5776

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• mhflux/MAlphaFitter.cc (modified) (3 diffs)
• mhflux/MAlphaFitter.h (modified) (1 diff)
• mhflux/MHAlpha.cc (modified) (5 diffs)
• mhflux/MHFalseSource.cc (modified) (28 diffs)
• mhflux/MHFalseSource.h (modified) (3 diffs)
• mraw/MRawFileRead.cc (modified) (5 diffs)
• mraw/MRawFileRead.h (modified) (1 diff)
• mraw/MRawRunHeader.cc (modified) (3 diffs)

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc

@@ -227 +227 @@
         return kFALSE;
 
-    fChiSqSignal = fit.GetChiSqSignal();
-    fChiSqBg     = fit.GetChiSqBg();
+    fChiSqSignal  = fit.GetChiSqSignal();
+    fChiSqBg      = fit.GetChiSqBg();
     fCoefficients = fit.GetCoefficients();
 
@@ -393 +393 @@
     const TString name1(Form("TempAlpha%06d_on",  gRandom->Integer(1000000)));
     const TString name0(Form("TempAlpha%06d_off", gRandom->Integer(1000000)));
+
     TH1D *h1 = hon.ProjectionZ(name1, -1, 9999, bin, bin, "E");
     TH1D *h0 = hof.ProjectionZ(name0, -1, 9999, bin, bin, "E");
@@ -486 +487 @@
     if (scaleof!=0)
         of.Scale(scaleon/scaleof);
-}
+    else
+        of.Reset();
+}

trunk/MagicSoft/Mars/mhflux/MAlphaFitter.h

@@ -55 +55 @@
 public:
     // Implementing the function yourself is only about 5% faster
-    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)
+    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)
     {
         fName  = name  ? name  : "MAlphaFitter";

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -185 +185 @@
     const Int_t n = fHAlpha.GetNbinsY();
 
+    fHEnergy.SetEntries(0);
+
     for (int i=1; i<=n; i++)
     {
@@ -201 +203 @@
 
     const Int_t n = fHAlpha.GetNbinsX();
+
+    fHTheta.SetEntries(0);
 
     for (int i=1; i<=n; i++)
@@ -218 +222 @@
     fHTheta.Reset();
 
-    if (fName!=(TString)"MHAlphaOff")
-    {
-        MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff");
+    if (fName!=(TString)"MHAlphaOff" && fOffData==NULL)
+    {
+        MHAlpha *hoff = (MHAlpha*)pl->FindObject("MHAlphaOff", "MHAlpha");
         if (!hoff)
             *fLog << inf << "No MHAlphaOff [MHAlpha] found... using current data only!" << endl;
@@ -400 +404 @@
     // Reset Histogram
     fHAlphaTime.Reset();
+    fHAlphaTime.SetEntries(0);
 
     // Get number of bins
@@ -646 +651 @@
 
         h = fOffData->ProjectionZ("ProjAlphaOff", -1, 9999, -1, 9999, "E");
-        h->SetBit(TH1::kNoTitle|TH1::kNoStats);
+        h->SetBit(TH1::kNoTitle);
         h->SetXTitle("\\alpha [\\circ]");
         h->SetYTitle("Counts");

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();

trunk/MagicSoft/Mars/mhflux/MHFalseSource.h

@@ -39 +39 @@
     TH3D    fHist;               // Alpha vs. x and y
 
+    const TH3D *fHistOff;
+
     Double_t fRa;
     Double_t fDec;
@@ -46 +48 @@
 
     void ProjectAll(TH2D *h);
-    void ProjectOff(TH2D *h, TH2D *all);
-    void ProjectOn(TH2D *h, TH2D *all);
+    void ProjectOff(const TH3D &src, TH2D *h, TH2D *all);
+    void ProjectOn(const TH3D &src, TH2D *h, TH2D *all);
+    void ProjectOnOff(TH2D *h, TH2D *all);
 
     TObject *GetCatalog();
@@ -77 +80 @@
     void SetBgMeanMinus5() { SetBgMean(fBgMean-5); } //*MENU*
 
+    void SetOffData(const MHFalseSource &fs) {
+        fHistOff  = &fs.fHist;
+        fMinDist  =  fs.fMinDist;
+        fMaxDist  =  fs.fMaxDist;
+        fMaxDW    =  fs.fMaxDW;
+        fMinDW    =  fs.fMinDW;
+        fAlphaCut =  fs.fAlphaCut;
+        fBgMean   =  fs.fBgMean;
+    }
+
     void Paint(Option_t *opt="");
     void Draw(Option_t *option="");

trunk/MagicSoft/Mars/mraw/MRawFileRead.cc

@@ -165 +165 @@
 // This opens the next file in the list and deletes its name from the list.
 //
-Bool_t MRawFileRead::OpenNextFile()
+Int_t MRawFileRead::OpenNextFile()
 {
     //
@@ -209 +209 @@
     //
     if (!fRawRunHeader->ReadEvt(*fIn))
-        return kFALSE;
+        return kERROR;
 
     if (!(*fIn))
     {
         *fLog << err << "Error: Accessing file '" << name << "'" << endl;
-        return kFALSE;
+        return kERROR;
     }
 
@@ -242 +242 @@
     {
         *fLog << err << dbginf << "MTaskList not found... abort." << endl;
-        return kFALSE;
+        return kERROR;
     }
 
@@ -270 +270 @@
     fNumFile=0;
     fNumEvents=0;
-    return OpenNextFile();
+    return OpenNextFile()==kTRUE;
 }
 
@@ -348 +348 @@
         // If an event could not be read from file try to open new file
         //
-        if (!OpenNextFile())
-            return kFALSE;
+        const Int_t rc = OpenNextFile();
+        if (rc!=kTRUE)
+            return rc;
     }
     return kTRUE; 

trunk/MagicSoft/Mars/mraw/MRawFileRead.h

@@ -22 +22 @@
     UInt_t     fInterleave;
 
-    Bool_t OpenNextFile();
+    Int_t OpenNextFile();
 
     Int_t PreProcess(MParList *pList);

trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc

@@ -83 +83 @@
 
 const UShort_t MRawRunHeader::kMagicNumber      = 0xc0c0;
-const Byte_t   MRawRunHeader::kMaxFormatVersion =      5;
+const Byte_t   MRawRunHeader::kMaxFormatVersion =      2;
 
 // --------------------------------------------------------------------------
@@ -139 +139 @@
     if (fMagicNumber != kMagicNumber && fMagicNumber != kMagicNumber+1)
     {
-        *fLog << err << "Error: Wrong Magic Number (0x" << hex << fMagicNumber << "): Not a Magic File!" << endl;
+        *fLog << err << "ERROR - Wrong Magic Number (0x" << hex << fMagicNumber << "): Not a Magic File!" << endl;
         return kFALSE;
     }
@@ -151 +151 @@
     if (fFormatVersion>kMaxFormatVersion)
     {
-        *fLog << err << "WARNING - File vormat V" << fFormatVersion << " not implemented!" << endl;
+        *fLog << err << "ERROR - File vormat V" << fFormatVersion << " not implemented!" << endl;
         return kFALSE;
     }