Changeset 2176 for trunk/MagicSoft/Mars

Timestamp:

06/15/03 10:31:36 (21 years ago)

Author:

rwagner

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• NEWS (modified) (1 diff)
• mhist/MHOnSubtraction.cc (modified) (13 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+ 2003/06/13: Robert Wagner
+   * mhist/MHOnSubtraction.cc
+     - removed casts from double to Double_t found by gcc 3.3
+     - added MHOnSubtraction::CalcLightCurve, a methods towards a 
+       lightcurve 
+
 
  2003/06/13: Thomas Bretz (making Mars work with gcc 3.3 on Suse 8.2)
@@ -112 +119 @@
    * mreflector/MRflEvtData.cc, mreflector/MRflSinglePhoton.cc:
      - fixed definition of Print
+
 
 

trunk/MagicSoft/Mars/NEWS

@@ -3 +3 @@
  *** Version 0.9
    
+   - added signal subtraction for pure on data by means of fitting
+     the background in the off region or by performing a combined
+     signal/background fit. Provides necessary histograms for
+     obtaining energy spectra and a light curve.
+
    - added classes to perform and study the selection of the 
      2nd Level Trigger on MC data (example in triglvl2.C macro)

trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc

@@ -237 +237 @@
      Int_t lowerSignalEdge = centerBin;
      for (Int_t i=3; i < maxBin; i++) {
-       if ((double)alphaHisto.GetBinContent(centerBin-i) 
+       if ((Double_t)alphaHisto.GetBinContent(centerBin-i) 
            < outerEvents/outerBins) {
          lowerSignalEdge = centerBin-i+1;
@@ -247 +247 @@
      Int_t upperSignalEdge = centerBin;
      for (Int_t i=3; i < maxBin; i++) {
-       if ((double)alphaHisto.GetBinContent(centerBin+i) 
+       if ((Double_t)alphaHisto.GetBinContent(centerBin+i) 
            <= outerEvents/outerBins) {
          upperSignalEdge=centerBin+i-1;
@@ -255 +255 @@
      if (centerBin>upperSignalEdge) upperSignalEdge=centerBin;
 
-     double nOnInt = 0;
+     Double_t nOnInt = 0;
      for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++)
        nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins;
      
-     double nOffInt = (upperSignalEdge - lowerSignalEdge + 1) 
+     Double_t nOffInt = (upperSignalEdge - lowerSignalEdge + 1) 
        * (outerEvents/outerBins);
      
@@ -292 +292 @@
   alphaHisto.SetDirectory(NULL); //rwagner
   
-  double gausMean  = gp->GetParameter(1);
-  double gausSigma = gp->GetParameter(2);
+  Double_t gausMean  = gp->GetParameter(1);
+  Double_t gausSigma = gp->GetParameter(2);
 
 //   TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-signalRegionFactor*gausSigma+gausMean,
@@ -315 +315 @@
   // this allows us to
   // (1) determine n_{ON}, n_{OFF}
-  double scalingFactor = 
+  Double_t scalingFactor = 
     (alphaHisto.GetXaxis()->GetBinLowEdge(alphaHisto.GetNbinsX()+1)
      - alphaHisto.GetXaxis()->GetBinLowEdge(1)) / alphaHisto.GetNbinsX();
   
-  double nOnInt  = (gp->Integral(-signalRegionFactor*gausSigma+gausMean, 
+  Double_t nOnInt  = (gp->Integral(-signalRegionFactor*gausSigma+gausMean, 
                                  signalRegionFactor*gausSigma+gausMean) / scalingFactor);
-  double nOffInt = (p3->Integral(-signalRegionFactor*gausSigma+gausMean, 
+  Double_t nOffInt = (p3->Integral(-signalRegionFactor*gausSigma+gausMean, 
                                  signalRegionFactor*gausSigma+gausMean) / scalingFactor);
 
@@ -390 +390 @@
 
   if (gausPol) {
-    double gausMean  = gausPol->GetParameter(1);
-    double gausSigma = gausPol->GetParameter(2);
+    Double_t gausMean  = gausPol->GetParameter(1);
+    Double_t gausSigma = gausPol->GetParameter(2);
     po = new TF1("po"+funcName,"pol3(0)",
                   -signalRegionFactor*gausSigma+gausMean,
@@ -455 +455 @@
   } // pol
     
-  double nOn = 0;
-  double eNOn = 0;
+  Double_t nOn = 0;
+  Double_t eNOn = 0;
 
   Int_t binNumberLow = alphaHisto.GetXaxis()->FindBin(lowerBin);
@@ -469 +469 @@
   // Evaluate background
   
-  double nOff = 0;
-  double eNOff = 0;
+  Double_t nOff = 0;
+  Double_t eNOff = 0;
   for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
     Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
@@ -500 +500 @@
   *fLog << inf << "Significance: "<<sigLiMa<<" sigma "<<endl;
 
-  if (draw) {  
+  if (draw) { 
     TPaveText *pt = new TPaveText(0.11,.74,.57,.88,"NDC ");
     char tx[60];
@@ -514 +514 @@
     pt->Draw("");
   }  
-
   if (draw||drawPad) {
+
     // Plot significance vs alpha_cut
     
@@ -528 +528 @@
     
     for (Int_t i=0; i < maxBin; i++) {
-      double nOn = 0;  double eNOn = 0;
-      double nOff = 0; double eNOff = 0;
+      Double_t nOn = 0;  Double_t eNOn = 0;
+      Double_t nOff = 0; Double_t eNOff = 0;
       for (Int_t bin=centerBin-i; bin<centerBin+i+1; bin++) {
         nOn += alphaHisto.GetBinContent(bin);
@@ -802 +802 @@
     thetaHisto.Fit("expF");
     
-    double expoSlope = e->GetParameter(1);
+    Double_t expoSlope = e->GetParameter(1);
     
     *fLog << inf << "Determined slope for energy distribution is " 
           << expoSlope << endl;
 
-    double integralEvts =
+    Double_t integralEvts =
       e->Integral(aetHisto->GetYaxis()->GetBinLowEdge(1), 
                   aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
@@ -1066 +1066 @@
   return kTRUE;
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+// -----------------------------------------------------------------------
+//
+// Extraction of Gamma signal is performed on a TH3D histogram in
+// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
+// in energy and Theta.
+//
+Bool_t MHOnSubtraction::CalcLightCurve(TH3 *aetHisto, MParList *parlist, const Bool_t draw)
+{
+  // Analyze aetHisto
+  // -------------------------------------
+  Int_t alphaBins = aetHisto->GetNbinsX(); 
+  Int_t energyBins = aetHisto->GetNbinsY();  
+  Int_t timeBins =  aetHisto->GetNbinsZ();  
+
+  *fLog << "Total events:      " << aetHisto->GetEntries() << endl;
+  *fLog << "Total energy bins: " << energyBins << endl;
+  *fLog << "Total time bins:   " << timeBins << endl;
+
+  // We output results in an array of histograms to the parameter list.
+  // Create a template for such a histogram, needed by MH3
+  // -------------------------------------
+  MBinning *binsmh3x = new MBinning("BinningMH3X");
+  // dummy binning, real one follows some lines down
+  binsmh3x->SetEdges(timeBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
+                     aetHisto->GetZaxis()->GetBinLowEdge(timeBins+1));
+  parlist->AddToList(binsmh3x); 
+  
+  MH3 *signalHisto = new MH3(1); // Signal(t)
+  signalHisto->SetupFill(parlist);
+  parlist->AddToList(signalHisto);
+  signalHisto->SetName("MHOnSubtractionSignal");
+  signalHisto->SetTitle("Gamma Events");
+
+  MH3 *offHisto = new MH3(1); // Off(t)
+  offHisto->SetupFill(parlist);
+  parlist->AddToList(offHisto);
+  offHisto->SetName("MHOnSubtractionOff");
+  offHisto->SetTitle("Background Events");
+
+  MH3 *signifHisto = new MH3(1); // Significance(t)
+  signifHisto->SetupFill(parlist);
+  parlist->AddToList(signifHisto);
+  signifHisto->SetName("MHOnSubtractionSignif");
+  signifHisto->SetTitle("Significance");
+
+  TH1D& signalTH1DHist = (TH1D&)signalHisto->GetHist();
+  TH1D& offTH1DHist =    (TH1D&)offHisto->GetHist();
+  TH1D& signifTH1DHist = (TH1D&)signifHisto->GetHist();
+  
+  signalTH1DHist.Reset();
+  signalTH1DHist.SetTitle("Gamma Signal");
+  signalTH1DHist.SetXTitle("Time [MJD]");
+  signalTH1DHist.SetYTitle("Events");
+  signalHisto->SetBinning(&signalTH1DHist, aetHisto->GetZaxis());
+  signalTH1DHist.Sumw2();
+  
+  offTH1DHist.Reset();
+  offTH1DHist.SetTitle("Background Signal");
+  offTH1DHist.SetXTitle("Time [MJD]");
+  offTH1DHist.SetYTitle("Events");
+  offHisto->SetBinning(&offTH1DHist,  aetHisto->GetZaxis());
+  offTH1DHist.Sumw2();  
+
+  signifTH1DHist.Reset();
+  signifTH1DHist.SetTitle("Significance");
+  signifTH1DHist.SetXTitle("Time [MJD]");
+  signifTH1DHist.SetYTitle("Significance #sigma");
+  signifHisto->SetBinning(&signifTH1DHist, aetHisto->GetZaxis());
+  signifTH1DHist.Sumw2();
+ 
+  //The following histogram is an additional histogram needed for
+  //the light curve
+
+  TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
+  c4->Divide(8,7,0,0,0);
+
+ 
+
+  // The following loop fixes a common integration region for each
+  // energy bin and alerts when no significant excess is found.
+  
+  Double_t minLowerBin = -10; // minimum integration range
+  Double_t maxUpperBin =  10; // minimum integration range
+  Double_t maxAlpha =  70; // maximum integration range
+  Double_t sumSigLiMa = 0;
+
+  TH1F* alphaHisto[timeBins];
+
+  for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {  
+
+    *fLog << dbg << "--------------------------------------" << endl;
+    *fLog << dbg << "Processing ALPHA plots for time bin " << timeBin << endl;           
+    *fLog << dbg << "--------------------------------------" << endl;
+   
+    aetHisto->GetZaxis()->SetRange(timeBin, timeBin);
+
+    char hname[80];
+    sprintf(hname, "MHOnSubtractionTimeBin%d", timeBin);    
+    alphaHisto[timeBin-1] =
+      new TH1F(hname, "ALPHA histogram",
+               alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
+               aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
+    alphaHisto[timeBin-1]->SetDirectory(NULL);
+    alphaHisto[timeBin-1]->Sumw2();
+    alphaHisto[timeBin-1] = (TH1F*)aetHisto->Project3D("xe");
+    alphaHisto[timeBin-1]->SetName(hname);
+    alphaHisto[timeBin-1]->SetDirectory(NULL);
+
+    sprintf(hname, "%6.0f < t < %6.0f", 
+            aetHisto->GetZaxis()->GetBinLowEdge(timeBin),
+            aetHisto->GetZaxis()->GetBinLowEdge(timeBin+1));
+    *fLog << inf << "There are " << alphaHisto[timeBin-1]->GetEntries() 
+          << " entries in " << hname << endl;
+    alphaHisto[timeBin-1]->SetTitle(hname);
+      
+    // Find signal region and significance in histogram
+    Double_t lowerBin, upperBin, sSigLiMa;
+    char funcName[20];
+    sprintf(funcName, "%2d", timeBin);
+
+    Bool_t drawFit = kTRUE;
+
+    c4->cd(timeBin);
+//     alphaHisto[timeBin-1]->SetMarkerColor(3);
+    alphaHisto[timeBin-1]->DrawCopy();
+
+    c4->Update();
+
+    fSigniPlotColor = 0;
+ ;
+    Bool_t fit = FitHistogram(*alphaHisto[timeBin-1], sSigLiMa, 
+                              lowerBin, upperBin, (Float_t)3.0, drawFit, 
+                              funcName);
+
+    if (fit) { 
+      if (sSigLiMa < 3)
+        *fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
+      sumSigLiMa+=sSigLiMa;
+   
+      // Evaluate integration range
+      lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
+      minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;    
+      upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
+      maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
+    }
+
+
+    
+
+
+
+  } //timeBin
+  
+  *fLog << inf << "=> Integration range will be " << minLowerBin << " < ALPHA < "
+        << maxUpperBin << "," << endl;    
+  *fLog << inf << "   Summed estimated significance is " << sumSigLiMa << endl;
+      
+  // we have an idea of what is going on in the ALPHA plot...
+  // So now we can indeed extract the signal.
+         
+  sumSigLiMa = 0;
+  for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {
+    *fLog << inf << "Processing ALPHA distribution for time bin " << timeBin << endl;
+    if (alphaHisto[timeBin-1]->GetEntries() == 0) {
+       *fLog << warn << "No attempt to extract a signal from 0 events." << endl;       
+       if (draw) {
+         TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
+         lab->SetBorderSize(0);
+         lab->Draw();  
+       }
+    } else {
+      char funcName[20];
+      sprintf(funcName, "%2d", timeBin);
+      
+      Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;    
+
+      Bool_t drawFit = kFALSE;     
+
+      ExtractSignal(*alphaHisto[timeBin-1],
+                    sigLiMa, minLowerBin, maxUpperBin, 
+                    gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
+                    funcName, NULL, 1/*drawBase*/);
+
+      sumSigLiMa += sigLiMa;      
+
+      fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;
+
+      // Fill Signal 
+      TH1D &h = (TH1D&)(signalHisto->GetHist());
+      h.SetBinContent(timeBin, gammaSignal);
+      h.SetBinError(timeBin, errorGammaSignal);
+      h.SetEntries(h.GetEntries()+gammaSignal);
+      
+      // Fill Off Events
+      TH1D &ho = (TH1D&)(offHisto->GetHist());
+      ho.SetBinContent(timeBin, off);
+      ho.SetBinError(timeBin, errorOff);
+      ho.SetEntries(ho.GetEntries()+off);
+      
+      // Fill Significance
+      TH1D &hg = (TH1D&)(signifHisto->GetHist());
+      hg.SetBinContent(timeBin, sigLiMa);
+      hg.SetEntries(hg.GetEntries()+sigLiMa);           
+    }
+    
+  }//timeBin
+
+  *fLog << inf << "Summed significance is " << sumSigLiMa << endl;
+
+
+  if (draw) {
+    Double_t sigLiMa, lowerBin, upperBin;    
+    FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
+    fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");
+
+    *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
+    // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
+    
+    // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
+    // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
+    
+    fChiSquareHisto->DrawClone();
+    fSignificanceHisto->DrawClone();
+    fSummedAlphaPlots->DrawClone();
+  }
+  
+  return kTRUE;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+