Changeset 2170

Timestamp:

06/11/03 13:44:37 (22 years ago)

Author:

rwagner

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mhist/MHOnSubtraction.cc (modified) (37 diffs)
• mhist/MHOnSubtraction.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+ 2003/06/11: Robert Wagner
+
+   * mhist/MHOnSubtraction.[h,cc]
+     - Some bugfixes, e.g. concerning binning of result histograms
+     - Improvements in output
 
 

trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc

@@ -26 +26 @@
 //////////////////////////////////////////////////////////////////////////////
 //                                                                          //
-//  MHOnSubtraction                                                            //
+//  MHOnSubtraction                                                         //
+//                                                                          //
 //                                                                          //
 //  extracts the gamma signal from a pure ON-signal given in an             //
 //  ALPHA-Energy-Theta histogram. The class will take this histogram from   //
-//  the parameter list and will provide a MHArray of MH3 histograms in      //
-//  energy, one per theta bin.                                              //
+//  the parameter list and will provide result histograms in the            //
+//  parameter list.                                                         //
+//                                                                          //
+//  This class still is work in progress.                                   //
+//                                                                          //
 //                                                                          //
 //////////////////////////////////////////////////////////////////////////////
 
-// This part of MARS is code still evolving. Please do not remove commentary
-// parts and clearly state where and how you changed the code, if you did so.
+// This part of MARS is code still evolving. Please do not change the code 
+// without prior feedback by the author.
 
 #include "MHOnSubtraction.h"
@@ -64 +68 @@
 // Default Constructor.
 //
-MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0),
-                                                            fMaxRedChiSq(0), 
-                                                            fSlope(20.0)
+MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0),  fMaxRedChiSq(0), fSlope(20.0)
 { 
   fName  = name  ? name  : "MHOnSubtraction";
@@ -96 +98 @@
 // measurements have to be provided.
 //
+// This function underestimates the true significance for it does not take
+// into account errors on the event numbers. A more accurate variation wil
+// be provided soon
+//
 Double_t MHOnSubtraction::CalcSignificance(Double_t nOn, Double_t nOff, Double_t theta)
 { 
-
-  if (nOn<=0) {
+  if (nOn<=0)
     *fLog << warn << "Got " << nOn << " total events, " << flush;
-  }
-
-  if (nOff<=0) {
-    *fLog << warn << "Got " << nOff << " background events, " << flush;    
-  }
-
+  if (nOff<=0)
+    *fLog << warn << "Got " << nOff << " background events, " << flush;
   if (nOn<=0 || nOff<=0.0) {
     *fLog << warn << "returning significance of 0.0" << endl;
     return 0.0;
   }
-
   return sqrt(2*(nOn*log((1+theta)*nOn/(theta*(nOn+nOff)))
                  +nOff*log((1+theta)*(nOff/(nOn+nOff)))));
@@ -117 +117 @@
 
 // -----------------------------------------------------------------------
+//
+// This function takes a first look at a given ALPHA distribution
+// and determines if a fit is applicable.
 //
 // Fits a given TH1 containing an ALPHA distribution with a combined
@@ -177 +180 @@
 
    *fLog << dbg << "Plot contains " << 
-     outerEvents << " outer ev (" << outerEvents/outerBins << "/bin), " << flush;
-   *fLog << dbg <<
+     outerEvents << " outer ev (" << outerEvents/outerBins << "/bin), " <<
      innerEvents << " inner ev (" << innerEvents/innerBins << "/bin) " << endl;
 
   if ((outerBinsZero!=0) || (innerBinsZero != 0))
     *fLog << warn << "There are ";
-
   if (outerBinsZero != 0)
     *fLog << dbg << outerBinsZero << " empty outer bins ";
-
   if (innerBinsZero != 0)
     *fLog << dbg <<innerBinsZero << " empty inner bins ";
-
   if ((outerBinsZero!=0) || (innerBinsZero != 0))
     *fLog << endl;
 
   if (outerEvents == 0.0) {
-     *fLog << warn << "No events in OFF region. Assuming background = 0" << endl;
-
+     *fLog << warn << "No events in OFF region. Assuming background = 0" 
+           << endl;
      TF1 *po = new TF1("pol0"+funcName,"pol0(0)",-89.5,89.5);
      po->SetLineWidth(2);
@@ -203 +202 @@
      alphaHisto.GetListOfFunctions()->Add(po);
      alphaHisto.SetLineColor(97);
-     alphaHisto.DrawCopy(); //rwagner
-     po->Draw("SAME");
-
+     if (draw) {
+       alphaHisto.DrawCopy(); //rwagner
+       po->Draw("SAME");
+     }
      sigLiMa = 0;
      lowerBin = 1;
@@ -213 +213 @@
      return kFALSE; //No gaus fit applied
   }
+
   if (outerEvents/outerBins < 2.5) {
     *fLog << warn << "Only " << /*setprecision(2) <<*/ outerEvents/outerBins 
@@ -224 +225 @@
      po->SetParameter(0,outerEvents/outerBins);
      alphaHisto.GetListOfFunctions()->Add(po);    
-     alphaHisto.DrawCopy(); //rwagner
-     po->Draw("SAME");
+     if (draw) {
+       alphaHisto.DrawCopy(); //rwagner
+       po->Draw("SAME");
+     }
     
      Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0);
@@ -231 +234 @@
        alphaHisto.GetNbinsX()- centerBin : centerBin;
      
-
      Int_t lowerSignalEdge = centerBin;
      for (Int_t i=3; i < maxBin; i++) {
-       if ((double)alphaHisto.GetBinContent(centerBin-i) < outerEvents/outerBins) {
+       if ((double)alphaHisto.GetBinContent(centerBin-i) 
+           < outerEvents/outerBins) {
          lowerSignalEdge = centerBin-i+1;
          break;
@@ -243 +246 @@
      Int_t upperSignalEdge = centerBin;
      for (Int_t i=3; i < maxBin; i++) {
-       if ((double)alphaHisto.GetBinContent(centerBin+i) <= outerEvents/outerBins) {
+       if ((double)alphaHisto.GetBinContent(centerBin+i) 
+           <= outerEvents/outerBins) {
          upperSignalEdge=centerBin+i-1;
          break;
@@ -251 +255 @@
 
      double nOnInt = 0;
-     for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++) {
+     for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++)
        nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins;
-     }
      
-     double nOffInt = (upperSignalEdge - lowerSignalEdge + 1) * (outerEvents/outerBins);
+     double nOffInt = (upperSignalEdge - lowerSignalEdge + 1) 
+       * (outerEvents/outerBins);
      
      sigLiMa = CalcSignificance(nOnInt, nOffInt, 1);
-    
-     *fLog << inf << "Estimated significance is " << sigLiMa << " sigma " << endl;   
+
+     if (sigLiMa < 3) 
+       alphaHisto.SetLineColor(2);   
+
+   
+     *fLog << inf << "Estimated significance is " << sigLiMa 
+           << " sigma " << endl;   
      *fLog << inf << "Signal region is " 
-           << lowerBin << " < ALPHA < " << upperBin << " (Most likely you wanna ignore this)" << endl;
+           << lowerBin << " < ALPHA < " << upperBin 
+           << " (Most likely you wanna ignore this)" << endl;
 
      return kFALSE; //No gaus fit applied    
-
   }
   
@@ -325 +334 @@
   
   sigLiMa = CalcSignificance(nOnInt, nOffInt, 1);
-
-  *fLog << inf << "Fit estimates significance to be " << sigLiMa << " sigma " << endl; 
+//   if (sigLiMa < 3)
+//     alphaHisto.SetLineColor(2);   
+
+
+  *fLog << inf << "Fit estimates significance to be " 
+        << sigLiMa << " sigma " << endl; 
 
   *fLog << inf << "Fit yields signal region to be " 
@@ -336 +349 @@
 }
 
+
+
+
+
 // -----------------------------------------------------------------------
 //
 // Does the actual extraction of the gamma signal. For performance 
-// reasons, the fit from MHOnSubtraction::FitHistogram is used and not redone.
+// reasons, fits already done by MHOnSubtraction::FitHistogram are used 
+// and not redone.
 // From it, a polynomial function for the background is evaluated and the
 // gamma signal is extracted in the region given by lowerBin < ALPHA < 
@@ -347 +365 @@
 //
 Bool_t MHOnSubtraction::ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa,
-                              Double_t &lowerBin, Double_t &upperBin,
-                              Double_t &gammaSignal, Double_t &errorGammaSignal,
-                              Double_t &off, Double_t &errorOff,
-                              Float_t signalRegionFactor, const Bool_t draw, TString funcName)
+   Double_t &lowerBin, Double_t &upperBin,
+   Double_t &gammaSignal, Double_t &errorGammaSignal,
+   Double_t &off, Double_t &errorOff,
+   Float_t signalRegionFactor, const Bool_t draw, TString funcName,
+   TPad *drawPad, Int_t drawBase)
 {
   TF1 *gausPol = alphaHisto.GetFunction("gauspol3"+funcName);
@@ -358 +377 @@
     *fLog << err << "Fatal: ALPHA histogram has no gauspol or pol "  
           << " fit attached to it." << endl;
-//     if (draw) {
-      TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no fit)");
-      lab->SetFillStyle(3000);
-      lab->SetBorderSize(0);
-      lab->DrawClone();  
-      lab->SetBit(kCanDelete);
-
-//     }
+    TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no fit)");
+    lab->SetFillStyle(3000);
+    lab->SetBorderSize(0);
+    lab->DrawClone();  
+    lab->SetBit(kCanDelete);
     return kFALSE;
   } 
@@ -414 +430 @@
 
   if (pol) {
-
     po = pol;
     
@@ -449 +464 @@
   } //for bin
   eNOn = sqrt(eNOn);
-  
-  
+     
   // Evaluate background
   
@@ -463 +477 @@
   eNOff = sqrt(fabs(nOff));
   
-  if (nOn==0) { // there should not be a negative number of signal events
+  if (nOn==0)   // there should not be a negative number of signal events
     nOff=0;
-  }
-
+  
   if (nOff<0) { // there should not be a negative number of off events
     nOff=0;
@@ -472 +485 @@
   }
     
-
   *fLog << inf << "nEvts = " << nOn << "+-" << eNOn << ", ";
 
@@ -498 +510 @@
     pt->SetBorderSize(0);
     pt->SetTextAlign(12);
-    pt->Draw();
+    pt->Draw("");
   }  
 
-  if (draw) {
+  if (draw||drawPad) {
     // Plot significance vs alpha_cut
     
@@ -507 +519 @@
     Int_t maxBin    = centerBin > alphaHisto.GetNbinsX()- centerBin ? 
       alphaHisto.GetNbinsX()- centerBin : centerBin;
-    
-    
+         
     const Int_t totalBins = centerBin;
     Float_t alpha[totalBins];
     Float_t signi[totalBins];
-
+    Float_t maxSigni = 0;
+    
     for (Int_t i=0; i < maxBin; i++) {
       double nOn = 0;  double eNOn = 0;
@@ -529 +541 @@
          alphaHisto.GetXaxis()->GetBinLowEdge(centerBin-i))/2;
       signi[i] = CalcSignificance(nOn, nOff, 1);  
+      maxSigni = maxSigni > signi[i] ? maxSigni : signi[i];
     }
-        
-    TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600);
-    c3->cd();
+
+    if (!drawPad) {
+      TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600);
+      c3->cd();
+    }
+
+    if (drawPad) 
+      drawPad->cd(drawBase-1);
+    
     TGraph* gr = new TGraph(totalBins-2,alpha,signi);
-    gr->Draw("ACP");
-    gr->GetXaxis()->SetTitle("|#alpha_{max}|");
-    gr->GetYaxis()->SetTitle("Significance");
+    TString drawOpt = "L";
+    
+    if (draw || (drawPad && fSigniPlotIndex == 0))
+      drawOpt += "A";
+    
+    gr->Draw(drawOpt);
+    if (drawPad && fSigniPlotIndex == 0) {
+      gr->GetXaxis()->SetTitle("|#alpha_{max}|");
+      gr->GetYaxis()->SetTitle("Significance");
+    }
     gr->SetMarkerStyle(2);
     gr->SetMarkerSize(1);
-    gr->SetMarkerColor(4);
-    gr->SetLineColor(4);
-    gr->SetLineWidth(2);
+    gr->SetMarkerColor(4+fSigniPlotIndex);
+    gr->SetLineColor(4+fSigniPlotIndex);
+    gr->SetLineWidth(1);
     gr->SetTitle("Significance vs ALPHA integration range");
     gr->Draw("LP");
+
+    fSigniPlotIndex++;      
+    
   } //draw
 
@@ -724 +753 @@
     // So now we can indeed extract the signal.
          
-    cout << "STARTIGN REAL CALC "<< endl;
+//     cout << "STARTIGN REAL CALC "<< endl;
 
     sumSigLiMa = 0;
@@ -757 +786 @@
 
     //fit exponential function to data
-    TF1 *e = new TF1("expF","expo",0,5);
+    TF1* e = new TF1("expF","expo",0,5);
     e->SetLineWidth(3);
     e->SetLineColor(thetaBin);
@@ -847 +876 @@
   // -------------------------------------
   MBinning *binsmh3x = new MBinning("BinningMH3X");
+  // dummy binning, real one follows some lines down
   binsmh3x->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
-                     aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
+                        aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
   parlist->AddToList(binsmh3x);  
 
   MBinning *binsmh3y = new MBinning("BinningMH3Y");
+  // dummy binning, real one follows some lines down
   binsmh3y->SetEdges(thetaBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
                      aetHisto->GetZaxis()->GetBinLowEdge(thetaBins+1));
@@ -894 +925 @@
   signalTH2DHist.SetXTitle("E [GeV]");
   signalTH2DHist.SetYTitle("theta [deg]");
+  signalHisto->SetBinning(&signalTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
   signalTH2DHist.Sumw2();
   
@@ -900 +932 @@
   offTH2DHist.SetXTitle("E [GeV]");
   offTH2DHist.SetYTitle("theta [deg]");
+  offHisto->SetBinning(&offTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
   offTH2DHist.Sumw2();  
 
@@ -906 +939 @@
   signifTH2DHist.SetXTitle("E [GeV]");
   signifTH2DHist.SetYTitle("theta [deg]");
+  signifHisto->SetBinning(&signifTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
   signifTH2DHist.Sumw2();
  
@@ -1005 +1039 @@
   
     c4->cd((energyBins+3)*(thetaBin-1)+2);
-    
     sprintf(hname,"c4_%d",(energyBins+3)*(thetaBin-1)+2);     
     TPad* tp = (TPad*)(gROOT->FindObject(hname));
-    tp->SetLogx();
-    
+    tp->SetLogx();     
     signalTH1DHist.SetLineColor(2);
-    signalTH1DHist.DrawCopy();
-    
+    signalTH1DHist.DrawCopy();     
     offTH1DHist.SetLineColor(4);
-    offTH1DHist.DrawCopy("SAME");
-    
-    
+    offTH1DHist.DrawCopy("SAME");          
     c4->Update();
-    
-    
-    
+           
     signalTH2DHist.SetEntries(signalTH2DHist.GetEntries()+signalTH1DHist.GetEntries());
     offTH2DHist.SetEntries(offTH2DHist.GetEntries()+offTH1DHist.GetEntries());
@@ -1032 +1059 @@
   }
 
+
   c4->Update();
-  
+
   return kTRUE;
 }
@@ -1084 +1112 @@
 Bool_t MHOnSubtraction::TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t draw, TPad *drawPad, Int_t drawBase)
 {
+
+  fSigniPlotColor = 0;
+
   // Analyze aeHisto
   // -------------------------------------
@@ -1113 +1144 @@
     signalHisto->SetTitle("Extracted Signal");
     parlist->AddToList(signalHisto);
- 
+    signalHisto->SetBinning(&((TH1D&)signalHisto->GetHist()), 
+                            aeHisto->GetYaxis());
+
     offHisto = new MH3(1); // Off(E)
     offHisto->SetName("MHOnSubtractionOff");
     offHisto->SetTitle("Off Signal");
     parlist->AddToList(offHisto);
+    offHisto->SetBinning(&((TH1D&)offHisto->GetHist()), 
+                            aeHisto->GetYaxis());
     
     signifHisto = new MH3(1); // Significance(E)
@@ -1123 +1158 @@
     signifHisto->SetTitle("Significance");
     parlist->AddToList(signifHisto);
+    signifHisto->SetBinning(&((TH1D&)signifHisto->GetHist()), 
+                            aeHisto->GetYaxis());
+
   }
   if (fChiSquareHisto==0x0)
@@ -1142 +1180 @@
 
   TH1F* alphaHisto[energyBins];
+
+  fSigniPlotIndex = 0; // cf. ExtractSignal
  
   for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
@@ -1238 +1278 @@
                     sigLiMa, minLowerBin, maxUpperBin, 
                     gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
-                    funcName);
+                    funcName, drawPad, drawBase);
 
       sumSigLiMa += sigLiMa;      
@@ -1304 +1344 @@
 //
 // Extraction of Gamma signal is performed on a TH1 histogram in ALPHA
+// 
 //
 Bool_t MHOnSubtraction::Calc(TH1 *aHisto, MParList *parlist, 
-                             const Bool_t Draw, Float_t signalRegion)
+                             const Bool_t draw, Float_t signalRegion)
 { 
   // Find signal region and estimate significance
@@ -1323 +1364 @@
 
   ExtractSignal(*aHisto, sigLiMa, lowerBin, upperBin, 
-                gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3.0, kTRUE);
+                gammaSignal, errorGammaSignal, off, errorOff, 
+                (Float_t)3.0, draw);
   
   *fLog << inf << "Signal is " 

trunk/MagicSoft/Mars/mhist/MHOnSubtraction.h

@@ -17 +17 @@
 #include "TPaveLabel.h"
 #endif
+#ifndef ROOT_TFile
+#include "TFile.h"
+#endif
+
 
 
@@ -41 +45 @@
   Double_t fSlope;             // slope for exponential fit
   Int_t fThetaBin;
+  Int_t fSigniPlotIndex;
+  Int_t fSigniPlotColor;
+
+  Bool_t CalcAET(TH3D *histon, MParList *parlist, const Bool_t Draw);
+
+  Bool_t FitHistogram(TH1 &alphaHisto, Double_t &sigLiMa,
+                      Double_t &lowerBin, Double_t &upperBin,
+                      Float_t signalRegionFactor = 3, const Bool_t draw = kFALSE,
+                      TString funcName = "");
+
+  Bool_t ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa,
+                       Double_t &lowerBin, Double_t &upperBin,
+                       Double_t &gammaSignal, Double_t &errorGammaSignal,
+                       Double_t &off, Double_t &errorOff,
+                       Float_t signalRegionFactor = 3, const Bool_t draw = kFALSE,
+                       TString funcName = "", TPad *drawPad = 0, Int_t drawBase = 0);
+
+
 
 public:
@@ -55 +77 @@
 
   void SetExpoSlope(Double_t slope) { fSlope = slope; }
-
-  Bool_t FitHistogram(TH1 &alphaHisto, Double_t &sigLiMa,
-                      Double_t &lowerBin, Double_t &upperBin,
-                      Float_t signalRegionFactor = 3, const Bool_t draw = kFALSE,
-                      TString funcName = "");
-
-  Bool_t ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa,
-                       Double_t &lowerBin, Double_t &upperBin,
-                       Double_t &gammaSignal, Double_t &errorGammaSignal,
-                       Double_t &off, Double_t &errorOff,
-                       Float_t signalRegionFactor = 3, const Bool_t draw = kFALSE,
-                       TString funcName = "");
   
   Bool_t Calc(MParList *parlist, const Bool_t Draw);
-  Bool_t CalcAET(TH3D *histon, MParList *parlist, const Bool_t Draw);
 
   Bool_t Calc(TH3 *histon, MParList *parlist, const Bool_t Draw);
-  Bool_t TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t Draw, TPad *drawPad = 0, Int_t drawBase = 0);
+  Bool_t TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t Draw, 
+                 TPad *drawPad = 0, Int_t drawBase = 0);
   Bool_t Calc(TH1 *histon, MParList *parlist, const Bool_t Draw,
               Float_t signalRegion = 0);