Changeset 6890 for trunk/MagicSoft/Mars/mhflux

Timestamp:

03/29/05 09:56:16 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhflux

Files:

• FluxLinkDef.h (modified) (1 diff)
• MHAlpha.cc (modified) (5 diffs)
• MHEnergyEst.cc (modified) (11 diffs)
• MHEnergyEst.h (modified) (1 diff)
• MHFalseSource.cc (modified) (11 diffs)

trunk/MagicSoft/Mars/mhflux/FluxLinkDef.h

@@ -10 +10 @@
 #pragma link C++ class MHEnergyEst+;
 #pragma link C++ class MHFalseSource+;
+#pragma link C++ class MHEnergyEst+;
 #pragma link C++ class MHEffectiveOnTime+;
 

trunk/MagicSoft/Mars/mhflux/MHAlpha.cc

@@ -138 +138 @@
     MBinning binsa, binse, binst;
     binsa.SetEdges(18, 0, 90);
-    binse.SetEdgesLog(25, 10, 100000);
+    binse.SetEdgesLog(15, 10, 100000);
     binst.SetEdgesCos(50, 0, 60);
     binse.Apply(fHEnergy);
@@ -467 +467 @@
 
         alpha  = hil->GetAlpha();
+
         if (fHillas)
             size = fHillas->GetSize();
@@ -472 +473 @@
         theta  = fPointPos ? fPointPos->GetZd()   : 0;
     }
+
+    //if (size>0)
+    //    alpha /= (2.4 + 1.13*(log10((energy-14)/0.37)-5)*(log10((energy-14)/0.37)-5))/15;
 
     // enhance histogram if necessary
@@ -563 +567 @@
             }
         }
+        else
+            hon->SetMinimum(0);
         FitEnergyBins();
         FitThetaBins();
@@ -748 +754 @@
     Int_t n = fHAlpha.GetNbinsY();
     Int_t nc = (Int_t)(TMath::Sqrt((Float_t)n-1)+1);
-    c->Divide(nc, nc, 0, 0);
+    c->Divide(nc, nc, 1e-10, 1e-10);
 
     // Do not store the 'final' result in fFit

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -147 +147 @@
 
     fChisq = 0;
+    fHEnergy.Reset();
+    fHImpact.Reset();
+    fHResolution.Reset();
 
     return kTRUE;
@@ -159 +162 @@
     const Double_t eest  = fEnergy->GetEnergy();
     const Double_t etru  = fMatrix ? GetVal(0) : fMcEvt->GetEnergy();
+    const Double_t imp   = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
     const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
-    const Double_t imp   = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
     const Double_t res   = (eest-etru)/etru;
 
     fHEnergy.Fill(eest, etru, theta, w);
-    fHResolution.Fill(eest, etru, res, w);
-    fHImpact.Fill(imp, res, w);
-
-    fChisq += res*res;
+    fHResolution.Fill(eest, etru, TMath::Abs(res), w);
+    fHImpact.Fill(imp, TMath::Abs(res), w);
+
+    fChisq += TMath::Abs(res);//*res;
+    fBias  += res;
 
     return kTRUE;
@@ -175 +179 @@
 {
     fChisq /= GetNumExecutions();
-
-    fResult->SetVal(fChisq);
-
-    *fLog << all << "Mean Energy Resoltuion: " << Form("%.1f%%", TMath::Sqrt(fChisq)*100) << endl;
+    fBias  /= GetNumExecutions();
+
+    Double_t res = fChisq; //TMath::Sqrt(fChisq - fBias*fBias);
+
+    fResult->SetVal(TMath::IsNaN(res)?0:res);/// GetNumExecutions());
+
+    *fLog << all << "Mean Energy Resoltuion: " << Form("%.1f%%", fResult->GetVal()*100) << endl;
+    *fLog << all << "Energy Bias at:         " << Form("%.1f%%", fBias*100) << endl;
 
     return kTRUE;
@@ -223 +231 @@
         {
             pad->GetPad(1)->GetPad(2)->cd(1);
-            /*h =*/ fHImpact.ProjectionX("Impact", -1, 9999, "e");
-
-            pad->GetPad(1)->GetPad(2)->cd(2);
-            if ((hx=(TH1D*)gPad->FindObject("EnergyEst_z")))
+            if ((hx=(TH1D*)gPad->FindObject("EnergyEst_ez")))
             {
-                TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_z");
+                TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum_ez");
                 hx->Reset();
                 hx->Add(h2);
                 delete h2;
             }
+
+            //pad->GetPad(1)->GetPad(2)->cd(2);
+            ///*h =*/ fHImpact.ProjectionX("Impact", -1, 9999, "e");
         }
     }
@@ -273 +281 @@
 
     TH1D *hx = 0;
-    if ((hx=(TH1D*)gPad->FindObject("Prof")))
+    if ((hx=(TH1D*)gPad->FindObject(Form("Prof%s", h.GetName()))))
     {
-        hx = hyx->ProfileX("Prof", -1, 9999, "s");
+        hx = hyx->ProfileX(Form("Prof%s", h.GetName()), -1, 9999, "s");
 
         if (logy && hx->GetMaximum()>0)
@@ -287 +295 @@
     gPad->SetLogx();
 
-    gROOT->GetListOfCleanups()->Add(gPad); // WHY?
+    //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
 
     TH2D *h2 = (TH2D*)h.Project3D(how);
-    TH1D *h1 = h2->ProfileX("Prof", -1, 9999, "s");
-
-    h1->SetDirectory(NULL);
-    //h1->SetBit(kCanDelete);
-    h1->SetLineWidth(2);
-    h1->SetLineColor(kRed); // PROBLEM!
-    h1->SetStats(kFALSE);
-
     h2->SetDirectory(NULL);
     h2->SetBit(kCanDelete);
     h2->SetFillColor(kBlue);
 
-    h1->Draw("E3");
-    h2->Draw("boxsame");
+    TH1D *h1 = h2->ProfileX(Form("Prof%s", h.GetName()), -1, 9999, "s");
+    h1->SetDirectory(NULL);
+    h1->SetBit(kCanDelete);
+    h1->SetLineWidth(2);
+    h1->SetLineColor(kRed);
+    h1->SetFillStyle(4000);
+    h1->SetStats(kFALSE);
+
+ 
+    //h1->Draw("E3");
+    h2->Draw();
     h1->Draw("Chistsame");
 
@@ -323 +332 @@
 
     pad->SetBorderMode(0);
-    pad->Divide(2, 1, 0, 0);
+    pad->Divide(2, 1, 1e-10, 1e-10);
 
     TH1 *h;
@@ -330 +339 @@
     gPad->SetBorderMode(0);
 
-    gPad->Divide(1, 2, 0, 0);
+    gPad->Divide(1, 2, 1e-10, 1e-10);
 
     TVirtualPad *pad2 = gPad;
@@ -365 +374 @@
 
     TVirtualPad *pad3 = gPad;
-    pad3->Divide(2, 1, 0, 0);
+    pad3->Divide(2, 1, 1e-10, 1e-10);
     pad3->cd(1);
-    gPad->SetBorderMode(0);
+    gPad->SetBorderMode(0);/*
     h = fHImpact.ProjectionX("Impact", -1, 9999, "e");
     h->SetBit(TH1::kNoStats);
@@ -374 +383 @@
     h->SetXTitle("Impact [m]");
     h->SetBit(kCanDelete);
-    h->Draw();
-
-    pad3->cd(2);
-    gPad->SetBorderMode(0);
+    h->Draw();*/
     h = fHEnergy.Project3D("ez");
     h->SetTitle("Distribution of Theta");
@@ -386 +392 @@
     h->Draw();
 
+    pad3->cd(2);
+    gPad->SetBorderMode(0);
+
     pad->cd(2);
     gPad->SetBorderMode(0);
 
-    gPad->Divide(1, 3, 0, 0);
+    gPad->Divide(1, 3, 1e-10, 1e-10);
     pad2 = gPad;
 

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.h

@@ -35 +35 @@
 
     Double_t fChisq;
+    Double_t fBias;
 
     TH1 *MakePlot(TH3 &h, const char *how);

trunk/MagicSoft/Mars/mhflux/MHFalseSource.cc

@@ -410 +410 @@
 
     // Get projection for range
-    TH2D *p = (TH2D*)src.Project3D("yx_off");
+    TH2D *p = (TH2D*)src.Project3D("yx_off_NULL");
 
     // Reset range
     axe.SetRange(0,9999);
 
+//#if ROOT_VERSION_CODE < ROOT_VERSION(4,02,00)
     // Move contents from projection to h2
     h2->Reset();
@@ -422 +423 @@
     // Delete p
     delete p;
+/*#else
+    p->Scale(all->GetMaximum());
+    p->Divide(all);
+#endif*/
 
     // Set Minimum as minimum value Greater Than 0
@@ -442 +447 @@
 
     // Get projection for range
-    TH2D *p = (TH2D*)src.Project3D("yx_on");
+    TH2D *p = (TH2D*)src.Project3D("yx_on_NULL");
 
     // Reset range
     axe.SetRange(0,9999);
 
+//#if ROOT_VERSION_CODE < ROOT_VERSION(4,02,00)
     // Move contents from projection to h3
     h3->Reset();
@@ -454 +460 @@
     // Delete p
     delete p;
+/*#else
+    p->Scale(all->GetMaximum());
+    p->Divide(all);
+#endif*/
 
     // Set Minimum as minimum value Greater Than 0
@@ -469 +479 @@
 
     // Get projection for range
-    TH2D *p = (TH2D*)fHist.Project3D(Form("yx_%d", gRandom->Uniform(999999)));
-    p->SetDirectory(0);
+#if ROOT_VERSION_CODE < ROOT_VERSION(4,02,00)
+    TH2D *p = (TH2D*)fHist.Project3D("yx_all");
 
     // Move contents from projection to h3
@@ -476 +486 @@
     h3->Add(p);
     delete p;
+#else
+    fHist.Project3D("yx_all");
+#endif
 
     // Set Minimum as minimum value Greater Than 0
@@ -518 +531 @@
     TH2D* h5;
 
-    /*
-     fHistProjAll  = Form("All_%p",  this);
-     fHistProjOn   = Form("On_%p",   this);
-     fHistProjOff  = Form("Off_%p",  this);
-     fHistProjDiff = Form("Diff_%p", this);
-     fHistProjAll  = Form("All_%p",  this);
-     */
-
     // Update projection of all-events
     padsave->GetPad(2)->cd(3);
@@ -559 +564 @@
         const Int_t nx = h4->GetXaxis()->GetNbins();
         const Int_t ny = h4->GetYaxis()->GetNbins();
-        //const Int_t nr = nx*nx + ny*ny;
 
         Int_t maxx=nx/2;
@@ -566 +570 @@
         Int_t max = h4->GetBin(nx, ny);
 
+        h4->SetEntries(0);
         for (int ix=1; ix<=nx; ix++)
             for (int iy=1; iy<=ny; iy++)
@@ -651 +656 @@
     AppendPad("");
 
-    pad->Divide(1, 2, 0, 0.03);
+    pad->Divide(1, 2, 1e-10, 0.03);
 
 //    TObject *catalog = GetCatalog();
@@ -1048 +1053 @@
     gStyle->SetPalette(1, 0);
 
-    c->Divide(3,2, 0, 0);
+    c->Divide(3,2, 1e-10, 1e-10);
     c->cd(1);
     gPad->SetBorderMode(0);