Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7403) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7404) @@ -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), fWeighting(kTRUE) + fNoThetaWeights(kFALSE) { fName = name ? name : "MJSpectrum"; @@ -601,20 +601,6 @@ 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); @@ -652,7 +638,4 @@ 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)*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); } @@ -726,4 +702,6 @@ spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) *spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1)); } + + // Display dN/dE*E^2 for conveinience fDisplay->AddTab("Check"); gPad->SetLogx(); @@ -741,9 +719,10 @@ 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(kBlue); + f.SetLineColor(kGreen); spectrum.Fit(&f, "NIM", "", 100, 5000); f.DrawCopy("same"); @@ -773,6 +752,8 @@ TArrayD res(2); res[0] = f.GetParameter(0); - res[1] = f.GetParameter(1);*/ - + res[1] = f.GetParameter(1); + + return res; + */ /* // Plot other spectra from Whipple @@ -807,5 +788,4 @@ f.SetLineStyle(kDashed); f.DrawCopy("same"); - return res; */ } @@ -870,6 +850,5 @@ // Calls PlotSame // -//#include -Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale, TH1D &result) const +Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const { *fLog << inf << "Reading from file: " << fPathIn << endl; @@ -890,5 +869,5 @@ } - TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); + TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); if (!excess) return kFALSE; @@ -906,5 +885,5 @@ excess->SetTitle("Excess events vs Size (data/black, mc/blue)"); - excess = (TH1D*)excess->DrawCopy("E2"); + excess = excess->DrawCopy("E2"); // Don't do this on the original object! excess->SetStats(kFALSE); @@ -918,72 +897,7 @@ if ((o=plist.FindObject("ExcessMC"))) { - TH1 *histsel = (TH1D*)o->FindObject(""); + TH1 *histsel = (TH1F*)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(); @@ -996,5 +910,4 @@ histsel->SetStats(kFALSE); - // Do some Chi^2 and Kolmogorov tests fLog->Separator("Kolmogorov Test"); histsel->KolmogorovTest(excess, "DX"); @@ -1002,5 +915,4 @@ const Double_t p = histsel->Chi2Test(excess, "P"); - // Plot result TLatex tex; tex.SetBit(TLatex::kTextNDC); @@ -1024,179 +936,4 @@ 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; @@ -1305,314 +1042,168 @@ // ------------------------- Final loop -------------------------- - 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); - } + *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); // check if output should be written Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 7403) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 7404) @@ -34,5 +34,4 @@ Bool_t fRawMc; Bool_t fNoThetaWeights; - Bool_t fWeighting; // Read Input @@ -50,6 +49,5 @@ TArrayD FitSpectrum(TH1D &spectrum) const; TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) 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 DisplaySize(MParList &plist, Double_t scale) const; Bool_t PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const; @@ -63,5 +61,4 @@ 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); Index: trunk/MagicSoft/Mars/msql/MSQLServer.cc =================================================================== --- trunk/MagicSoft/Mars/msql/MSQLServer.cc (revision 7403) +++ trunk/MagicSoft/Mars/msql/MSQLServer.cc (revision 7404) @@ -444,4 +444,5 @@ Int_t MSQLServer::SelectDataBase(const char *dbname) /*FOLD00*/ { + fDataBase = dbname; return fType==kIsServer ? fServ->SelectDataBase(dbname) : 0; } @@ -496,5 +497,5 @@ switch (fType) { - case kIsServer: return Form("%s://%s:%d", fServ->GetDBMS(), fServ->GetHost(), fServ->GetPort()); + case kIsServer: return Form("%s://%s:%d/%s", fServ->GetDBMS(), fServ->GetHost(), fServ->GetPort(), fDataBase.Data()); case kIsDataBase: return GetNameDataBase(); case kIsTable: return GetNameTable();