Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7402) +++ trunk/MagicSoft/Mars/Changelog (revision 7403) @@ -27,4 +27,8 @@ * mtools/MTFillMatrix.[h,cc]: - added the possibility to use PreCuts like in MJOptim + + * mjobs/MJSpectrum.[h,cc]: + - implemented long awaiting E^2dN/dE plot + - updated displaying result (formula) Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7402) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7403) @@ -76,5 +76,5 @@ #include "MTaskEnv.h" #include "MFillH.h" -#include "MHillasCalc.h" +//#include "MHillasCalc.h" //#include "MSrcPosCalc.h" #include "MContinue.h" @@ -87,5 +87,5 @@ : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0), fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE), - fNoThetaWeights(kFALSE) + fNoThetaWeights(kFALSE), fWeighting(kTRUE) { fName = name ? name : "MJSpectrum"; @@ -551,4 +551,88 @@ } +TString MJSpectrum::FormString(const TF1 &f, Byte_t type) +{ + const Double_t p0 = f.GetParameter(0); + const Double_t p1 = f.GetParameter(1); + + const Double_t e0 = f.GetParError(0); + const Double_t e1 = f.GetParError(1); + + const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1))); + const Double_t exp = TMath::Power(10, np); + + const Float_t fe0 = TMath::Log10(e0); + const Float_t fe1 = TMath::Log10(e1); + + // calc number of (gueltige ziffern) digits to be displayed + const char *f0 = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0"; + const char *f1 = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0"; + + TString str; + switch (type) + { + case 0: + { + const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1); + const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0); + + str = Form(fmt0.Data(), p1/exp, e1/exp, np); + str += Form(fmt1.Data(), p0, e0); + str += "\\frac{ph}{TeVm^{2}s}"; + } + break; + case 1: + str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF()); + break; + case 2: + str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF())); + break; + } + return str; +} + +TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const +{ + TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4); + f.SetParameter(0, -2.5); + f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000))); + f.SetLineColor(kBlue); + spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped + f.DrawCopy("same"); + + /* + const Double_t p0 = f.GetParameter(0); + const Double_t p1 = f.GetParameter(1); + + const Double_t e0 = f.GetParError(0); + const Double_t e1 = f.GetParError(1); + + const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1))); + const Double_t exp = TMath::Power(10, np); + */ + TString str = FormString(f); + /* + str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np); + str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0); + str += "\\frac{ph}{TeVm^{2}s}"; + */ + TLatex tex; + tex.SetTextSize(0.045); + tex.SetBit(TLatex::kTextNDC); + tex.SetTextAlign(31); + tex.DrawLatex(0.89, 0.935, str); + + str = FormString(f, 1); + tex.DrawLatex(0.89, 0.83, str); + + str = FormString(f, 2); + tex.DrawLatex(0.89, 0.735, str); + + TArrayD res(2); + res[0] = f.GetParameter(0); + res[1] = f.GetParameter(1); + return res; +} + // -------------------------------------------------------------------------- // @@ -568,4 +652,7 @@ cout << "Effective On time: " << ontime << "s" << endl; + +// for (int i=0; i2) + { + //Float_t E = spectrum.GetXaxis()->GetBinWidth(i+1); + //Float_t c = log10(E)*3.012e-1+1.905e-1; + //spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*c); + //spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*c); + spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1)); + spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1)); + } spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000); - spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000); + spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) /spectrum.GetBinWidth(i+1)*1000); } @@ -623,11 +717,33 @@ spectrum.SetMinimum(1e-12); spectrum.SetXTitle("E [GeV]"); - spectrum.SetYTitle("N/TeVsm^{2}"); + spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]"); + TH1D *spec = (TH1D*)spectrum.DrawCopy(); + + // Calculate Spectrum * E^2 + for (int i=0; iAddTab("Check"); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + spectrum.SetName("FluxStd"); + spectrum.SetMarkerStyle(kFullDotMedium); + spectrum.SetTitle("Differential flux times E^{2}"); + spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]"); + spectrum.SetDirectory(0); spectrum.DrawCopy(); + // Fit Spectrum + c1.cd(4); + return FitSpectrum(*spec); + /* TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4); f.SetParameter(0, -2.87); f.SetParameter(1, 1.9e-6); - f.SetLineColor(kGreen); + f.SetLineColor(kBlue); spectrum.Fit(&f, "NIM", "", 100, 5000); f.DrawCopy("same"); @@ -657,5 +773,5 @@ TArrayD res(2); res[0] = f.GetParameter(0); - res[1] = f.GetParameter(1); + res[1] = f.GetParameter(1);*/ /* @@ -691,6 +807,6 @@ f.SetLineStyle(kDashed); f.DrawCopy("same"); + return res; */ - return res; } @@ -754,5 +870,6 @@ // Calls PlotSame // -Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const +//#include +Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale, TH1D &result) const { *fLog << inf << "Reading from file: " << fPathIn << endl; @@ -773,5 +890,5 @@ } - TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); + TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); if (!excess) return kFALSE; @@ -789,5 +906,5 @@ excess->SetTitle("Excess events vs Size (data/black, mc/blue)"); - excess = excess->DrawCopy("E2"); + excess = (TH1D*)excess->DrawCopy("E2"); // Don't do this on the original object! excess->SetStats(kFALSE); @@ -801,7 +918,72 @@ if ((o=plist.FindObject("ExcessMC"))) { - TH1 *histsel = (TH1F*)o->FindObject(""); + TH1 *histsel = (TH1D*)o->FindObject(""); if (histsel) { + // Calculate the weights such that the number of total + // MC events after cuts will stay constant + + /* + Int_t first = -1; + for (int i=0; iGetNbinsX(); i++) + if (excess->GetBinContent(i+1)>0) + { + first = i+1; + break; + } + + Int_t last = -1; + for (int i=excess->GetNbinsX(); i>0; i--) + if (excess->GetBinContent(i)>0) + { + last = i; + break; + }*/ + /* + Int_t first = 1; + Int_t last = excess->GetNbinsX(); + + const Double_t xmin = TMath::Log10(excess->GetBinLowEdge(first)); + const Double_t xmax = TMath::Log10(excess->GetBinLowEdge(last+1)); + const Int_t xnum = last-first+1; + + // Double_t xmin = TMath::Log10(excess->GetXaxis()->GetXmin()); + // Double_t xmax = TMath::Log10(excess->GetXaxis()->GetXmax()); + // Int_t xnum = excess->GetNbinsX(); + + TArrayD arr1(xnum); + TArrayD arr2(xnum); + for (int i=0; iGetBinContent(i+first); + Double_t v2 = histsel->GetBinContent(i+first); + arr1[i] = v1<1 ? -1 : TMath::Log10(v1); + arr2[i] = v2<1 ? -1 : TMath::Log10(v2); + } + cout << endl << endl; + + TSpline5 sp1("Excess", xmin, xmax, arr1.GetArray(), xnum+1); + TSpline5 sp2("MC", xmin, xmax, arr2.GetArray(), xnum+1); + + MBinning bins(xnum*5, TMath::Power(10, xmin), TMath::Power(10, xmax), "", "log"); + bins.Apply(result); + + for (int i=0; iGetBinLowEdge(i+1); + Double_t up = result.GetXaxis()->GetBinUpEdge(i+1); + + Double_t x = TMath::Log10(lo*up)/2; + Double_t v1 = sp1.Eval(x); + Double_t v2 = sp2.Eval(x); + Double_t rc = v1<-0.9 || v2<-0.9 ? 0 : TMath::Power(10, v1-v2); + result.SetBinContent(i+1, rc); + } + */ + excess->Copy(result); + result.Divide(histsel); + result.SetDirectory(0); + + // Apply additional absolute scale factor if (scale<0) scale = excess->Integral()/histsel->Integral(); @@ -814,4 +996,5 @@ histsel->SetStats(kFALSE); + // Do some Chi^2 and Kolmogorov tests fLog->Separator("Kolmogorov Test"); histsel->KolmogorovTest(excess, "DX"); @@ -819,4 +1002,5 @@ const Double_t p = histsel->Chi2Test(excess, "P"); + // Plot result TLatex tex; tex.SetBit(TLatex::kTextNDC); @@ -840,4 +1024,179 @@ c.cd(6); PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale); + + return kTRUE; +} + +Bool_t MJSpectrum::DisplaySpectrumW(TH1D &energy0, TH1 &energy, TH1 &hsize, TH1D &weights, TH2D &hall, Double_t scale, TArrayD &res) const +{ + *fLog << inf << "Reading from file: " << fPathIn << endl; + + TFile file(fPathIn, "READ"); + if (!file.IsOpen()) + { + *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl; + return kFALSE; + } + + file.cd(); + MStatusArray arr; + if (arr.Read()<=0) + { + *fLog << "MStatusDisplay not found in file... abort." << endl; + return kFALSE; + } + + TH1D *exc = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); + if (!exc) + return kFALSE; + + // energy: energy distribution after weighting + // energy0: energy distribution before weighting + // weights: weights to weight events depending on their size + + TCanvas &c = fDisplay->AddTab("SpectrumW"); + c.Divide(3,2); + c.cd(1); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + TH1D *e = (TH1D*)exc->DrawCopy("E2"); + // Don't do this on the original object! + e->SetName("ExcessData"); + e->SetTitle("Excess events vs Size (data/black, mc/blue)"); + e->SetDirectory(0); + e->SetStats(kFALSE); + e->SetBit(TH1::kNoStats); + e->ResetBit(TH1::kNoTitle); + e->SetMarkerStyle(kFullDotMedium); + e->SetFillColor(kBlack); + e->SetFillStyle(0); + e->SetName("Excess "); + e->SetDirectory(0); + + e = (TH1D*)hsize.DrawCopy("E1 same"); + // Don't do this on the original object! + e->SetName("ExcessMC"); + e->SetDirectory(0); + e->SetBit(kCanDelete); + e->SetBit(TH1::kNoStats); + e->ResetBit(TH1::kNoTitle); + e->SetLineColor(kBlue); + e->SetStats(kFALSE); + + c.cd(2); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetGrid(); + + weights.SetName("WeightsSize"); + weights.SetDirectory(0); + weights.SetBit(TH1::kNoStats); + weights.ResetBit(TH1::kNoTitle); + weights.SetTitle("Size weights: size(data)/size(mc)"); + weights.DrawCopy(); + + c.cd(4); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + energy.SetBit(TH1::kNoStats); + energy.ResetBit(TH1::kNoTitle); + energy.SetTitle("Normalized energy distribution before (black) and after (blue) weighting"); + energy.SetLineColor(kBlue); + TH1 *he = energy.DrawCopy(); + TH1 *hf = energy0.DrawCopy("same"); + he->Scale(1./he->Integral()); + hf->Scale(1./hf->Integral()); + he->SetName("EnergyOrig"); + hf->SetName("EnergyWeighted"); + he->SetDirectory(0); + hf->SetDirectory(0); + + c.cd(5); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetGrid(); + energy0.SetBit(TH1::kNoStats); + energy0.ResetBit(TH1::kNoTitle); + energy0.SetName("WeightsEnergy"); + energy0.SetTitle("Energy weights: E(before weighting)/E(after weighting)"); + energy0.SetDirectory(0); + energy0.Divide(&energy); + energy0.DrawCopy(); + + c.cd(3); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + TH1D *h1 = (TH1D*)hall.ProjectionY("ProjAllY", -1, 9999, "E"); + h1->SetTitle("Normalized spectrum before and after Energy weights applied"); + h1->SetDirectory(NULL); + h1->SetXTitle("E [GeV]"); + h1->SetBit(kCanDelete); + h1->SetBit(TH1::kNoStats); + he = h1->DrawCopy(); + h1->SetLineColor(kBlue); + // Convert the mc spectrum before weights to the mc spectrum after weights + // Therefor weights after/before are applied + h1->Divide(&energy0); + hf = h1->DrawCopy("same"); + he->Scale(1./he->GetBinContent(he->GetNbinsX()/2)); + hf->Scale(1./hf->GetBinContent(hf->GetNbinsX()/2)); + he->SetName("SpectrumOrig"); + hf->SetName("SpectrumWeighted"); + he->SetDirectory(0); + hf->SetDirectory(0); + + c.cd(6); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + h1->SetName("Spectrum"); + h1->SetTitle("Spectrum"); + h1->SetXTitle("E [GeV]"); + h1->SetYTitle("dN/dE [N/TeVsm^{2}]"); + h1->SetDirectory(0); + h1->Scale(1./scale); + + for (int i=0; iGetNbinsX(); i++) + { + h1->SetBinContent(i+1, h1->GetBinContent(i+1)/h1->GetBinWidth(i+1)*1000); + h1->SetBinError( i+1, h1->GetBinError(i+1) /h1->GetBinWidth(i+1)*1000); + } + + TH1D *h2 = (TH1D*)h1->DrawCopy(); + res = FitSpectrum(*h2); + + // Calculate and Display Spectrum * E^2 + TCanvas *c1 = fDisplay->GetCanvas("Check"); + if (c1) + { + c1->cd(); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + + for (int i=0; iGetNbinsX(); i++) + { + h1->SetBinContent(i+1, h1->GetBinContent(i+1)*h1->GetBinWidth(i+1)*h1->GetBinWidth(i+1)); + h1->SetBinError( i+1, h1->GetBinError(i+1) *h1->GetBinWidth(i+1)*h1->GetBinWidth(i+1)); + } + h1->SetName("FluxW"); + h1->SetTitle("Differential flux times E^{2}"); + h1->SetDirectory(0); + h1->SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]"); + h1->SetLineColor(kBlue); + h1->SetMarkerColor(kBlue); + h1->SetMarkerStyle(kFullDotMedium); + h1->Draw("same"); + } return kTRUE; @@ -946,168 +1305,314 @@ // ------------------------- Final loop -------------------------- - *fLog << endl; - fLog->Separator("Calculate efficiencies"); - *fLog << endl; - - MTaskList tlist2; - plist.Replace(&tlist2); - - MReadMarsFile read("Events"); - read.DisableAutoScheme(); - set.AddFilesOn(read); - - // Selector to get correct (final) theta-distribution - temp1.SetXTitle("MPointingPos.fZd"); - - MH3 mh3(temp1); - - MFEventSelector2 sel2(mh3); - sel2.SetHistIsProbability(); - - MContinue contsel(&sel2); - contsel.SetInverted(); - - // Get correct source position - //MSrcPosCalc calc; - - // Calculate corresponding Hillas parameters - /* - MHillasCalc hcalc1; - MHillasCalc hcalc2("MHillasCalcAnti"); - hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc); - hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); - hcalc2.SetNameHillasSrc("MHillasSrcAnti"); - hcalc2.SetNameSrcPosCam("MSrcPosAnti"); - */ - - // Fill collection area and energy estimator (unfolding) - // Make sure to use the same binning for MHCollectionArea and MHEnergyEst - MHCollectionArea area; - area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); - MHEnergyEst hest; - - MFillH fill3(&area, "", "FillCollectionArea"); - MFillH fill4(&hest, "", "FillEnergyEst"); - MFillH fill5("MHThreshold", "", "FillThreshold"); - fill3.SetWeight(); - fill4.SetWeight(); - fill5.SetWeight(); - fill3.SetNameTab("ColArea"); - fill4.SetNameTab("E-Est"); - fill5.SetNameTab("Threshold"); - - MH3 hsize("MHillas.fSize"); - hsize.SetName("ExcessMC"); - hsize.Sumw2(); - - MBinning bins(size, "BinningExcessMC"); - plist.AddToList(&hsize); - plist.AddToList(&bins); - - MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre"); - MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost"); - MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); - MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); - MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); - MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost"); - MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); - MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC"); - fill1a.SetNameTab("PreCut"); - fill2a.SetNameTab("PostCut"); - fill3a.SetNameTab("VsSize"); - fill4a.SetNameTab("HilExt"); - fill5a.SetNameTab("HilSrc"); - fill6a.SetNameTab("ImgPar"); - fill7a.SetNameTab("NewPar"); - fill8a.SetBit(MFillH::kDoNotDisplay); - fill1a.SetWeight(); - fill2a.SetWeight(); - fill3a.SetWeight(); - fill4a.SetWeight(); - fill5a.SetWeight(); - fill6a.SetWeight(); - fill7a.SetWeight(); - fill8a.SetWeight(); - - MEnergyEstimate est; - MTaskEnv taskenv1("EstimateEnergy"); - taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); - - tlist2.AddToList(&read); - if (!fRawMc && fNoThetaWeights) - tlist2.AddToList(&contsel); - //tlist2.AddToList(&calc); - //tlist2.AddToList(&hcalc1); - //tlist2.AddToList(&hcalc2); - tlist2.AddToList(&weight); - tlist2.AddToList(&fill1a); - tlist2.AddToList(fCut0); - tlist2.AddToList(fCut1); - tlist2.AddToList(fCut2); - tlist2.AddToList(fCut3); - tlist2.AddToList(&taskenv1); - tlist2.AddToList(&fill3); - tlist2.AddToList(&fill4); - tlist2.AddToList(&fill5); - tlist2.AddToList(&fill2a); - tlist2.AddToList(&fill3a); - tlist2.AddToList(&fill4a); - tlist2.AddToList(&fill5a); - tlist2.AddToList(&fill6a); - tlist2.AddToList(&fill7a); - tlist2.AddToList(&fill8a); - //tlist2.AddToList(&fill9a); - - MEvtLoop loop2("FillMonteCarlo"); // ***** fName ***** - loop2.SetParList(&plist); - loop2.SetDisplay(fDisplay); - loop2.SetLogStream(fLog); - - if (!SetupEnv(loop2)) - return kFALSE; - - if (!loop2.Eventloop(fMaxEvents)) - { - *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl; - return kFALSE; - } - - if (!loop2.GetDisplay()) - { - *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; - return kFALSE; - } - - gLog.Separator("Energy Estimator"); - if (plist.FindObject("EstimateEnergy")) - plist.FindObject("EstimateEnergy")->Print(); - - gLog.Separator("Spectrum"); - - // -------------------------- Spectrum ---------------------------- - - // Calculate and display spectrum (N/TeVsm^2 at 1TeV) - TArrayD res(DisplaySpectrum(area, excess, hest, ontime)); - - // Spectrum fitted (convert res[1] from TeV to GeV) - TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0])); - - // Number of events this spectrum would produce per s and m^2 - Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax()); - - // scale with effective collection area to get the event rate (N/s) - // scale with the effective observation time to absolute observed events - n *= area.GetCollectionAreaAbs()*ontime; // N - - // Now calculate the scale factor from the number of events - // produced and the number of events which should have been - // observed with our telescope in the time ontime - const Double_t scale = n/area.GetEntries(); - - // Print normalization constant - cout << "MC normalization factor: " << scale << endl; - - // Overlay normalized plots - DisplaySize(plist, scale); + TH1D weights; + TH1D energy0; + weights.SetDirectory(0); + + const Int_t n = fWeighting ? 2 : 1; + for (int i=0; iSeparator("Calculate efficiencies"); + *fLog << endl; + + if (i==1) + weight.SetWeightsSize(&weights); + + + MTaskList tlist2; + plist.Replace(&tlist2); + + MReadMarsFile read("Events"); + read.DisableAutoScheme(); + set.AddFilesOn(read); + + // Selector to get correct (final) theta-distribution + temp1.SetXTitle("MPointingPos.fZd"); + + MH3 mh3(temp1); + + MFEventSelector2 sel2(mh3); + sel2.SetHistIsProbability(); + + MContinue contsel(&sel2); + contsel.SetInverted(); + +//*** // Get correct source position +//*** MSrcPosCalc calc; +//*** +//*** // Calculate corresponding Hillas parameters +//*** MHillasCalc hcalc1; +//*** MHillasCalc hcalc2("MHillasCalcAnti"); +//*** hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc); +//*** hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); +//*** hcalc2.SetNameHillasSrc("MHillasSrcAnti"); +//*** hcalc2.SetNameSrcPosCam("MSrcPosAnti"); + + // Fill collection area and energy estimator (unfolding) + // Make sure to use the same binning for MHCollectionArea and MHEnergyEst + MHCollectionArea area; + area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); + MHEnergyEst hest; + + MFillH fill3(&area, "", "FillCollectionArea"); + MFillH fill4(&hest, "", "FillEnergyEst"); + MFillH fill5("MHThreshold", "", "FillThreshold"); + fill3.SetWeight(); + fill4.SetWeight(); + fill5.SetWeight(); + fill3.SetNameTab("ColArea"); + fill4.SetNameTab("E-Est"); + fill5.SetNameTab("Threshold"); + + MH3 hsize("MHillas.fSize"); + hsize.SetName("ExcessMC"); + hsize.Sumw2(); + + MH3 energy("MMcEvt.fEnergy"); + energy.SetName("EnergyEst"); // FIXME: Should be EnergyMC, but binning has name EnergyEst + energy.SetLogy(); + energy.Sumw2(); + + MBinning bins(size, "BinningExcessMC"); + plist.AddToList(&hsize); + plist.AddToList(&energy); + plist.AddToList(&bins); + + MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre"); + MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost"); + MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); + MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); + MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); + MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost"); + MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); + MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC"); + MFillH fill9a("EnergyEst [MH3]", "", "FillEnergy"); + fill1a.SetNameTab("PreCut"); + fill2a.SetNameTab("PostCut"); + fill3a.SetNameTab("VsSize"); + fill4a.SetNameTab("HilExt"); + fill5a.SetNameTab("HilSrc"); + fill6a.SetNameTab("ImgPar"); + fill7a.SetNameTab("NewPar"); + fill9a.SetNameTab("Energy"); + fill8a.SetBit(MFillH::kDoNotDisplay); + fill9a.SetBit(MFillH::kDoNotDisplay); + fill1a.SetWeight(); + fill2a.SetWeight(); + fill3a.SetWeight(); + fill4a.SetWeight(); + fill5a.SetWeight(); + fill6a.SetWeight(); + fill7a.SetWeight(); + fill8a.SetWeight(); + fill9a.SetWeight(); + + MEnergyEstimate est; + MTaskEnv taskenv1("EstimateEnergy"); + taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); + + if (i==1) + { + est.SetName("SetEnergyMC"); + est.SetRule("MMcEvt.fEnergy"); + } + + tlist2.AddToList(&read); + if (!fRawMc && fNoThetaWeights) + tlist2.AddToList(&contsel); +//*** tlist2.AddToList(&calc); +//*** tlist2.AddToList(&hcalc1); +//*** tlist2.AddToList(&hcalc2); + tlist2.AddToList(&weight); + if (i==n-1) + tlist2.AddToList(&fill1a); + tlist2.AddToList(fCut0); + tlist2.AddToList(fCut1); + tlist2.AddToList(fCut2); + tlist2.AddToList(fCut3); + tlist2.AddToList(i==1 ? (MTask*)&est : (MTask*)&taskenv1); + tlist2.AddToList(&fill3); + if (i==0) + tlist2.AddToList(&fill4); + tlist2.AddToList(&fill5); + if (i==n-1) + { + tlist2.AddToList(&fill2a); + tlist2.AddToList(&fill3a); + tlist2.AddToList(&fill4a); + tlist2.AddToList(&fill5a); + tlist2.AddToList(&fill6a); + tlist2.AddToList(&fill7a); + } + tlist2.AddToList(&fill8a); + tlist2.AddToList(&fill9a); + + MEvtLoop loop2("FillMonteCarlo"); // ***** fName ***** + loop2.SetParList(&plist); + loop2.SetDisplay(fDisplay); + loop2.SetLogStream(fLog); + +//***<<<<<<< MJSpectrum.cc + if (!SetupEnv(loop2)) + return kFALSE; +//***======= +//*** MFEventSelector2 sel2(mh3); +//*** sel2.SetHistIsProbability(); +//*** +//*** MContinue contsel(&sel2); +//*** contsel.SetInverted(); +//*** +//*** // Get correct source position +//*** //MSrcPosCalc calc; +//*** +//*** // Calculate corresponding Hillas parameters +//*** /* +//*** MHillasCalc hcalc1; +//*** MHillasCalc hcalc2("MHillasCalcAnti"); +//*** hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc); +//*** hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc); +//*** hcalc2.SetNameHillasSrc("MHillasSrcAnti"); +//*** hcalc2.SetNameSrcPosCam("MSrcPosAnti"); +//*** */ +//*** +//*** // Fill collection area and energy estimator (unfolding) +//*** // Make sure to use the same binning for MHCollectionArea and MHEnergyEst +//*** MHCollectionArea area; +//*** area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); +//*** MHEnergyEst hest; +//*** +//*** MFillH fill3(&area, "", "FillCollectionArea"); +//*** MFillH fill4(&hest, "", "FillEnergyEst"); +//*** MFillH fill5("MHThreshold", "", "FillThreshold"); +//*** fill3.SetWeight(); +//*** fill4.SetWeight(); +//*** fill5.SetWeight(); +//*** fill3.SetNameTab("ColArea"); +//*** fill4.SetNameTab("E-Est"); +//*** fill5.SetNameTab("Threshold"); +//*** +//*** MH3 hsize("MHillas.fSize"); +//*** hsize.SetName("ExcessMC"); +//*** hsize.Sumw2(); +//*** +//*** MBinning bins(size, "BinningExcessMC"); +//*** plist.AddToList(&hsize); +//*** plist.AddToList(&bins); +//*** +//*** MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre"); +//*** MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost"); +//*** MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost"); +//*** MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost"); +//*** MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost"); +//*** MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost"); +//*** MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost"); +//*** MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC"); +//*** fill1a.SetNameTab("PreCut"); +//*** fill2a.SetNameTab("PostCut"); +//*** fill3a.SetNameTab("VsSize"); +//*** fill4a.SetNameTab("HilExt"); +//*** fill5a.SetNameTab("HilSrc"); +//*** fill6a.SetNameTab("ImgPar"); +//*** fill7a.SetNameTab("NewPar"); +//*** fill8a.SetBit(MFillH::kDoNotDisplay); +//*** fill1a.SetWeight(); +//*** fill2a.SetWeight(); +//*** fill3a.SetWeight(); +//*** fill4a.SetWeight(); +//*** fill5a.SetWeight(); +//*** fill6a.SetWeight(); +//*** fill7a.SetWeight(); +//*** fill8a.SetWeight(); +//***>>>>>>> 1.17 + + if (!loop2.Eventloop(fMaxEvents)) + { + *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl; + return kFALSE; + } + +//***<<<<<<< MJSpectrum.cc + if (!loop2.GetDisplay()) + { + *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; + return kFALSE; + } +//***======= +//*** tlist2.AddToList(&read); +//*** if (!fRawMc && fNoThetaWeights) +//*** tlist2.AddToList(&contsel); +//*** //tlist2.AddToList(&calc); +//*** //tlist2.AddToList(&hcalc1); +//*** //tlist2.AddToList(&hcalc2); +//*** tlist2.AddToList(&weight); +//*** tlist2.AddToList(&fill1a); +//*** tlist2.AddToList(fCut0); +//*** tlist2.AddToList(fCut1); +//*** tlist2.AddToList(fCut2); +//*** tlist2.AddToList(fCut3); +//*** tlist2.AddToList(&taskenv1); +//*** tlist2.AddToList(&fill3); +//*** tlist2.AddToList(&fill4); +//*** tlist2.AddToList(&fill5); +//*** tlist2.AddToList(&fill2a); +//*** tlist2.AddToList(&fill3a); +//*** tlist2.AddToList(&fill4a); +//*** tlist2.AddToList(&fill5a); +//*** tlist2.AddToList(&fill6a); +//*** tlist2.AddToList(&fill7a); +//*** tlist2.AddToList(&fill8a); +//*** //tlist2.AddToList(&fill9a); +//*** +//*** MEvtLoop loop2("FillMonteCarlo"); // ***** fName ***** +//*** loop2.SetParList(&plist); +//*** loop2.SetDisplay(fDisplay); +//*** loop2.SetLogStream(fLog); +//***>>>>>>> 1.17 + + gLog.Separator("Energy Estimator"); + if (plist.FindObject("EstimateEnergy")) + plist.FindObject("EstimateEnergy")->Print(); + + gLog.Separator("Spectrum"); + + // -------------------------- Spectrum ---------------------------- + + + TArrayD res; + if (i==0) + { + // Calculate and display spectrum (N/TeVsm^2 at 1TeV) + res = DisplaySpectrum(area, excess, hest, ontime); + + ((TH1D&)energy.GetHist()).Copy(energy0); + } + else + { + const Double_t scale = ontime*area.GetCollectionAreaAbs(); + if (!DisplaySpectrumW(energy0, energy.GetHist(), hsize.GetHist(), weights, fSimpleMode ? hist : (TH2D&)mh1.GetHist(), scale, res)) + return kFALSE; + } + + // Spectrum fitted (convert res[1] from TeV to GeV) + TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0])); + + // Number of events this spectrum would produce per s and m^2 + Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax()); + + // scale with effective collection area to get the event rate (N/s) + // scale with the effective observation time to absolute observed events + n *= area.GetCollectionAreaAbs()*ontime; // N + + // Now calculate the scale factor from the number of events + // produced and the number of events which should have been + // observed with our telescope in the time ontime + const Double_t scale = i==0 ? n/area.GetEntries() : 1; + + // Print normalization constant + cout << "MC normalization factor: " << scale << endl; + + // Overlay normalized plots, and calculate weights for second loop + DisplaySize(plist, scale, weights); + } // check if output should be written Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 7402) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 7403) @@ -6,4 +6,5 @@ #endif +class TF1; class TH1; class TH1D; @@ -33,4 +34,5 @@ Bool_t fRawMc; Bool_t fNoThetaWeights; + Bool_t fWeighting; // Read Input @@ -46,6 +48,8 @@ void DisplayResult(const TH2D &mh1) const; Bool_t IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const; + TArrayD FitSpectrum(TH1D &spectrum) const; TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const; - Bool_t DisplaySize(MParList &plist, Double_t scale) const; + Bool_t DisplaySpectrumW(TH1D &energy0, TH1 &energy, TH1 &hsize, TH1D &weights, TH2D &hall, Double_t scale, TArrayD &res) const; + Bool_t DisplaySize(MParList &plist, Double_t scale, TH1D &result) const; Bool_t PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const; @@ -59,6 +63,9 @@ void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; } void EnableRawMc(Bool_t b=kTRUE) { fRawMc=b; } + void EnableWeighting(Bool_t b=kTRUE) { fWeighting=b; } void SetEnergyEstimator(const MTask *task); + + static TString FormString(const TF1 &f, Byte_t type=0); ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum