Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2175) +++ trunk/MagicSoft/Mars/Changelog (revision 2176) @@ -1,3 +1,10 @@ -*-*- 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,4 +119,5 @@ * mreflector/MRflEvtData.cc, mreflector/MRflSinglePhoton.cc: - fixed definition of Print + Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 2175) +++ trunk/MagicSoft/Mars/NEWS (revision 2176) @@ -3,4 +3,9 @@ *** 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) Index: trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 2175) +++ trunk/MagicSoft/Mars/mhist/MHOnSubtraction.cc (revision 2176) @@ -237,5 +237,5 @@ 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,5 +247,5 @@ 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,9 +255,9 @@ 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,6 +292,6 @@ 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,11 +315,11 @@ // 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,6 +390,6 @@ 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,6 +455,6 @@ } // pol - double nOn = 0; - double eNOn = 0; + Double_t nOn = 0; + Double_t eNOn = 0; Int_t binNumberLow = alphaHisto.GetXaxis()->FindBin(lowerBin); @@ -469,6 +469,6 @@ // Evaluate background - double nOff = 0; - double eNOff = 0; + Double_t nOff = 0; + Double_t eNOff = 0; for (Int_t bin=binNumberLow; binDraw(""); } - if (draw||drawPad) { + // Plot significance vs alpha_cut @@ -528,6 +528,6 @@ 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; binGetParameter(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,4 +1066,327 @@ 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 <<")!"< 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; +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +