/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz, 4/2005 ! ! Copyright: MAGIC Software Development, 2000-2005 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MJSpectrum // // Program to calculate spectrum // ///////////////////////////////////////////////////////////////////////////// #include "MJSpectrum.h" // Root #include #include #include #include #include #include #include #include #include // Environment #include "MLog.h" #include "MLogManip.h" #include "MStatusArray.h" #include "MStatusDisplay.h" // Container #include "MH3.h" #include "MBinning.h" #include "MDataSet.h" #include "MMcCorsikaRunHeader.h" // Spectrum #include "../mhflux/MAlphaFitter.h" #include "../mhflux/MHAlpha.h" #include "../mhflux/MHCollectionArea.h" //#include "../mhflux/MHThreshold.h" #include "../mhflux/MHEnergyEst.h" #include "../mhflux/MMcSpectrumWeight.h" // Eventloop #include "MEvtLoop.h" #include "MTaskList.h" #include "MParList.h" // Tasks/Filter #include "MReadMarsFile.h" #include "MReadMarsFile.h" #include "MFEventSelector2.h" #include "MFDataMember.h" #include "MEnergyEstimate.h" #include "MTaskEnv.h" #include "MFillH.h" #include "MHillasCalc.h" //#include "MSrcPosCalc.h" #include "MContinue.h" ClassImp(MJSpectrum); using namespace std; MJSpectrum::MJSpectrum(const char *name, const char *title) : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0), fEstimateEnergy(0), fCalcHadronness(0), fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE), fNoThetaWeights(kFALSE) { fName = name ? name : "MJSpectrum"; fTitle = title ? title : "Standard program to calculate spectrum"; } MJSpectrum::~MJSpectrum() { if (fCutQ) delete fCutQ; if (fCut0) delete fCut0; if (fCut1) delete fCut1; if (fCut2) delete fCut2; if (fCut3) delete fCut3; if (fCutS) delete fCutS; if (fEstimateEnergy) delete fEstimateEnergy; if (fCalcHadronness) delete fCalcHadronness; } // -------------------------------------------------------------------------- // // Setup a task estimating the energy. The given task is cloned. // void MJSpectrum::SetEnergyEstimator(const MTask *task) { if (fEstimateEnergy) delete fEstimateEnergy; fEstimateEnergy = task ? (MTask*)task->Clone() : 0; } Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const { if (task) { delete task; task = 0; } task = (MTask*)gFile->Get(name); if (!task) { if (!mustexist) return kTRUE; *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl; return kFALSE; } if (!task->InheritsFrom(MTask::Class())) { *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl; delete task; return kFALSE; } task->SetName(name); if (dynamic_cast(task)) dynamic_cast(task)->SetAllowEmpty(); return kTRUE; } void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const { fLog->Separator("Alpha Fitter"); *fLog << all; fit.Print("result"); fLog->Separator("Used Cuts"); fCutQ->Print(); fCut0->Print(); fCut1->Print(); fCutS->Print(); fCut2->Print(); fCut3->Print(); } // -------------------------------------------------------------------------- // // Read the first MMcCorsikaRunHeader from the RunHeaders tree in // the dataset. // The simulated energy range and spectral slope is initialized from // there. // In the following eventloops the forced check in MMcSpectrumWeight // ensures, that the spectral slope and energy range doesn't change. // Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const { fLog->Separator("Initialize energy weighting"); // Read Resources if (!CheckEnv(w)) { *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl; return kFALSE; } // Read the first corsika RunHeader from the MC files TChain chain("RunHeaders"); if (!set.AddFilesOn(chain)) return kFALSE; MMcCorsikaRunHeader *h=0; chain.SetBranchAddress("MMcCorsikaRunHeader.", &h); chain.GetEntry(1); if (!h) { *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl; return kFALSE; } // Propagate the run header to MMcSpectrumWeight if (!w.Set(*h)) { *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl; return kFALSE; } // Print new setup of MMcSpectrumWeight w.Print(); return kTRUE; } Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2) { *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 -1; } MStatusArray arr; if (arr.Read()<=0) { *fLog << "MStatusDisplay not found in file... abort." << endl; return -1; } TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime"); TH1D *size = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); if (!vstime || !size) return -1; vstime->Copy(h1); size->Copy(h2); h1.SetDirectory(0); h2.SetDirectory(0); if (fDisplay) arr.DisplayIn(*fDisplay, "Hist"); if (!ReadTask(fCutQ, "CutQ")) return -1; if (!ReadTask(fCut0, "Cut0")) return -1; if (!ReadTask(fCut1, "Cut1")) return -1; if (!ReadTask(fCut2, "Cut2")) return -1; if (!ReadTask(fCut3, "Cut3")) return -1; if (!ReadTask(fCalcHadronness, "CalcHadronness", kFALSE)) return -1; TObjArray arrread; TIter Next(plist); TObject *o=0; while ((o=Next())) if (o->InheritsFrom(MBinning::Class())) arrread.Add(o); arrread.Add(plist.FindObject("MAlphaFitter")); if (!ReadContainer(arrread)) return -1; return vstime->Integral(); } Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const { // Some debug output fLog->Separator("Compiling original MC distribution"); // The name of the input container is MMcEvtBasic weight.SetNameMcEvt("MMcEvtBasic"); // Get the corresponding rule for the weights const TString w(weight.GetFormulaWeights()); // Set back to MMcEvt weight.SetNameMcEvt(); *fLog << inf << "Using weights: " << w << endl; *fLog << "Please stand by, this may take a while..." << flush; if (fDisplay) fDisplay->SetStatusLine1("Compiling MC distribution..."); // Create chain TChain chain("OriginalMC"); if (!set.AddFilesOn(chain)) return kFALSE; // Prepare histogram h.Reset(); h.Sumw2(); // Fill histogram from chain h.SetDirectory(gROOT); if (h.InheritsFrom(TH2::Class())) { h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)"); h.SetXTitle("\\Theta [\\circ]"); h.SetYTitle("E [GeV]"); h.SetZTitle("Counts"); chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff"); } else { h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)"); h.SetXTitle("\\Theta [\\circ]"); h.SetYTitle("Counts"); chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff"); } h.SetDirectory(0); *fLog << "done." << endl; if (fDisplay) fDisplay->SetStatusLine2("done."); if (h.GetEntries()>0) return kTRUE; *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl; return h.GetEntries()>0; } Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const { // Display some stuff if (fDisplay) { TCanvas &c = fDisplay->AddTab("ZdDist"); c.Divide(2,2); // On-Time vs. Theta c.cd(1); gPad->SetBorderMode(0); temp1.DrawCopy(); // Number of MC events (produced) vs Theta c.cd(2); gPad->SetBorderMode(0); temp2.SetName("NVsTheta"); temp2.DrawCopy("hist"); c.cd(4); gPad->SetBorderMode(0); c.cd(3); gPad->SetBorderMode(0); } // Calculate the Probability temp1.Divide(&temp2); temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral()); // Some cosmetics: Name, Axis, etc. temp1.SetName("ProbVsTheta"); temp1.SetTitle("Probability vs. Zenith Angle to choose MC events"); temp1.SetYTitle("Probability"); if (fDisplay) temp1.DrawCopy("hist"); return kTRUE; } // -------------------------------------------------------------------------- // // Display the final theta distribution. // Bool_t MJSpectrum::DisplayResult(const TH2D &h2) const { if (!fDisplay || !fDisplay->CdCanvas("ZdDist")) { *fLog << err << "ERROR - Display or tab ZdDist vanished... abort." << endl; return kFALSE; } TH1D &proj = *h2.ProjectionX(); proj.SetNameTitle("ThetaFinal", "Final Theta Distribution"); proj.SetXTitle("\\Theta [\\circ]"); proj.SetYTitle("Counts"); proj.SetLineColor(kBlue); proj.SetDirectory(0); proj.SetBit(kCanDelete); TVirtualPad *pad = gPad; pad->cd(4); proj.DrawCopy(); pad->cd(1); TH1D *theta = (TH1D*)gPad->FindObject("Theta"); if (!theta) { *fLog << err << "ERROR - Theta-Histogram vanished... cannot proceed." << endl; return kFALSE; } if (proj.GetMaximum()==0) { *fLog << err; *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap"; *fLog << " with the the Zenith Angle distribution of your observation."; theta->SetLineColor(kRed); return kFALSE;; } proj.Scale(theta->GetMaximum()/proj.GetMaximum()); theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum())); proj.Draw("same"); // Compare both histograms const Double_t prob = proj.Chi2Test(theta, ""); if (prob==1) return kTRUE; if (prob>0.99) { *fLog << inf; *fLog << "The Zenith Angle distribution of your Monte Carlos fits well"; *fLog << "with the the Zenith Angle distribution of your observation."; return kTRUE; } // <0.01: passen gar nicht // ==1: passen perfekt // >0.99: passer sehr gut if (prob<0.01) { *fLog << err; *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit"; *fLog << " with the the Zenith Angle distribution of your observation."; theta->SetLineColor(kRed); return kFALSE; } *fLog << warn; *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well"; *fLog << " with the the Zenith Angle distribution of your observation."; return kTRUE; } // -------------------------------------------------------------------------- // // Fills the excess histogram (vs E-est) from the events stored in the // ganymed result file and therefor estimates the energy. // // The resulting histogram excess-vs-energy ist copied into h2. // Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2)/*const*/ { // Try to find the class used to determin the signal! TString cls("MHAlpha"); if (fDisplay) { TCanvas *c = fDisplay->GetCanvas("Hist"); if (c) { TIter Next(c->GetListOfPrimitives()); TObject *obj=0; while ((obj=Next())) if (obj->InheritsFrom(MHAlpha::Class())) break; if (obj) cls = obj->ClassName(); } } // Now fill the histogram *fLog << endl; fLog->Separator("Refill Excess"); *fLog << endl; MTaskList tlist; plist.AddToList(&tlist); MReadTree read("Events"); read.DisableAutoScheme(); read.AddFile(fPathIn); MTaskEnv taskenv0("CalcHadronness"); taskenv0.SetDefault(fCalcHadronness); MEnergyEstimate est; MTaskEnv taskenv1("EstimateEnergy"); taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); MContinue *cont = new MContinue("", "CutS"); cont->SetAllowEmpty(); if (fCutS) delete fCutS; fCutS = cont; // FIXME: Create HistE and HistEOff to be able to modify it from // the resource file. MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff"); MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE"); MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData"); MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData"); fill1.SetFilter(&f0); fill2.SetFilter(&f1); tlist.AddToList(&read); //tlist.AddToList(&taskenv0); // not necessary, stored in file! tlist.AddToList(fCutS); tlist.AddToList(&taskenv1); tlist.AddToList(&f0); tlist.AddToList(&f1); tlist.AddToList(&fill1); tlist.AddToList(&fill2); // by setting it here it is distributed to all consecutive tasks tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); MEvtLoop loop("RefillExcess"); // ***** fName ***** loop.SetParList(&plist); loop.SetDisplay(fDisplay); loop.SetLogStream(fLog); if (!SetupEnv(loop)) return kFALSE; if (!loop.Eventloop()) { *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl; return kFALSE; } if (!loop.GetDisplay()) { *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; return kFALSE; } const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE"); if (!halpha) { *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl; return kFALSE; } halpha->GetHEnergy().Copy(h2); h2.SetDirectory(0); return kTRUE; } Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const { MTaskList tlist1; plist.Replace(&tlist1); MReadMarsFile readmc("OriginalMC"); //readmc.DisableAutoScheme(); if (!set.AddFilesOn(readmc)) return kFALSE; readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta"); readmc.EnableBranch("MMcEvtBasic.fEnergy"); mh1.SetLogy(); mh1.SetLogz(); mh1.SetName("ThetaE"); MFillH fill0(&mh1); //fill0.SetDrawOption("projx only"); MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst"); MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta"); if (bins2 && bins3) { bins2->SetName("BinningThetaEY"); bins3->SetName("BinningThetaEX"); } tlist1.AddToList(&readmc); tlist1.AddToList(&weight); temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg"); MH3 mh3mc(temp1); MFEventSelector2 sel1(mh3mc); sel1.SetHistIsProbability(); fill0.SetFilter(&sel1); if (!fRawMc) tlist1.AddToList(&sel1); tlist1.AddToList(&fill0); // by setting it here it is distributed to all consecutive tasks tlist1.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); MEvtLoop loop1("IntermediateLoop"); // ***** fName ***** loop1.SetParList(&plist); loop1.SetLogStream(fLog); loop1.SetDisplay(fDisplay); if (!SetupEnv(loop1)) return kFALSE; if (!loop1.Eventloop(fMaxEvents)) { *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl; return kFALSE; } if (!loop1.GetDisplay()) { *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; return kFALSE; } if (bins2 && bins3) { bins2->SetName("BinningEnergyEst"); bins3->SetName("BinningTheta"); } return kTRUE; } 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"); TString str = FormString(f); 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; } // -------------------------------------------------------------------------- // // Calculate the final spectrum from: // - collection area // - excess // - correction coefficients // - ontime // and display it // TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const { // Create copies of the histograms TH1D collarea(area.GetHEnergy()); TH1D spectrum(excess); TH1D weights; // Get spill-over corrections from energy estimation hest.GetWeights(weights); // Print effective on-time cout << "Effective On time: " << ontime << "s" << endl; // Setup spectrum plot spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)"); spectrum.SetYTitle("N/sm^{2}"); spectrum.SetDirectory(NULL); spectrum.SetBit(kCanDelete); // Divide by collection are and on-time spectrum.Scale(1./ontime); spectrum.Divide(&collarea); // Draw spectrum before applying spill-over corrections TCanvas &c1 = fDisplay->AddTab("Spectrum"); c1.Divide(2,2); c1.cd(1); gPad->SetBorderMode(0); gPad->SetLogx(); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); collarea.DrawCopy(); c1.cd(2); gPad->SetBorderMode(0); gPad->SetLogx(); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); TH1D *spec=(TH1D*)spectrum.DrawCopy(); FitSpectrum(*spec); c1.cd(3); gPad->SetBorderMode(0); gPad->SetLogx(); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); weights.DrawCopy(); // Apply spill-over correction (done't take the errors into account) // They are supposed to be identical with the errors of the // collection area and cancel out. //spectrum.Multiply(&weights); spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)"); spectrum.SetBit(TH1::kNoStats); for (int i=0; iSetBorderMode(0); gPad->SetLogx(); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); spectrum.SetMinimum(1e-12); spectrum.SetXTitle("E [GeV]"); spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]"); 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); spectrum.Fit(&f, "NIM", "", 100, 5000); 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; 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.DrawLatex(0.45, 0.935, str); str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF()); tex.DrawLatex(0.70, 0.83, str); TArrayD res(2); res[0] = f.GetParameter(0); res[1] = f.GetParameter(1); return res; */ /* // Plot other spectra from Whipple f.SetParameter(0, -2.45); f.SetParameter(1, 3.3e-7); f.SetRange(300, 8000); f.SetLineColor(kBlack); f.SetLineStyle(kDashed); f.DrawCopy("same"); // Plot other spectra from Cangaroo f.SetParameter(0, -2.53); f.SetParameter(1, 2.0e-7); f.SetRange(7000, 50000); f.SetLineColor(kBlack); f.SetLineStyle(kDashed); f.DrawCopy("same"); // Plot other spectra from Robert f.SetParameter(0, -2.59); f.SetParameter(1, 2.58e-7); f.SetRange(150, 1500); f.SetLineColor(kBlack); f.SetLineStyle(kDashed); f.DrawCopy("same"); // Plot other spectra from HEGRA f.SetParameter(0, -2.61); f.SetParameter(1, 2.7e-7); f.SetRange(1000, 20000); f.SetLineColor(kBlack); f.SetLineStyle(kDashed); f.DrawCopy("same"); */ } // -------------------------------------------------------------------------- // // Scale some image parameter plots using the scale factor and plot them // together with the corresponding MC histograms. // Called from DisplaySize // Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const { TString same(name); same += "Same"; TH1 *h1 = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab); TH1 *h2 = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab); if (!h1 || !h2) return kFALSE; TObject *obj = plist.FindObject(plot); if (!obj) { *fLog << warn << plot << " not in parameter list... skipping." << endl; return kFALSE; } TH1 *h3 = (TH1*)obj->FindObject(name); if (!h3) { *fLog << warn << name << " not found in " << plot << "... skipping." << endl; return kFALSE; } const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter"); const Double_t ascale = fit ? fit->GetScaleFactor() : 1; gPad->SetBorderMode(0); h2->SetLineColor(kBlack); h3->SetLineColor(kBlue); h2->Add(h1, -ascale); //h2->Scale(1./ontime); //h2->Integral()); h3->Scale(scale); //h3->Integral()); h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum())); h2 = h2->DrawCopy(); h3 = h3->DrawCopy("same"); // Don't do this on the original object! h2->SetStats(kFALSE); h3->SetStats(kFALSE); return kTRUE; } // -------------------------------------------------------------------------- // // Take a lot of histograms and plot them together in one plot. // Calls PlotSame // Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) 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; } TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist"); if (!excess) return kFALSE; // ------------------- Plot excess versus size ------------------- TCanvas &c = fDisplay->AddTab("Excess"); c.Divide(3,2); c.cd(1); gPad->SetBorderMode(0); gPad->SetLogx(); gPad->SetLogy(); gPad->SetGridx(); gPad->SetGridy(); excess->SetTitle("Excess events vs Size (data/black, mc/blue)"); excess = excess->DrawCopy("E2"); // Don't do this on the original object! excess->SetStats(kFALSE); excess->SetMarkerStyle(kFullDotMedium); excess->SetFillColor(kBlack); excess->SetFillStyle(0); excess->SetName("Excess "); excess->SetDirectory(0); TObject *o=0; if ((o=plist.FindObject("ExcessMC"))) { TH1 *histsel = (TH1F*)o->FindObject(""); if (histsel) { if (scale<0) scale = excess->Integral()/histsel->Integral(); histsel->Scale(scale); histsel->SetLineColor(kBlue); histsel->SetBit(kCanDelete); histsel = histsel->DrawCopy("E1 same"); // Don't do this on the original object! histsel->SetStats(kFALSE); fLog->Separator("Kolmogorov Test"); histsel->KolmogorovTest(excess, "DX"); fLog->Separator("Chi^2 Test"); const Double_t p = histsel->Chi2Test(excess, "P"); TLatex tex; tex.SetBit(TLatex::kTextNDC); tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100)); } } // -------------- Comparison of Image Parameters -------------- c.cd(2); PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale); c.cd(3); PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale); c.cd(4); PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale); c.cd(5); PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale); c.cd(6); PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale); return kTRUE; } // -------------------------------------------------------------------------- // void MJSpectrum::DisplayCutEfficiency(const MHCollectionArea &area0, const MHCollectionArea &area1) const { if (!fDisplay) return; const TH1D &trig = area0.GetHEnergy(); TH1D &cut = (TH1D&)*area1.GetHEnergy().Clone(); fDisplay->AddTab("CutEff"); gPad->SetBorderMode(0); gPad->SetFrameBorderMode(0); gPad->SetLogx(); gPad->SetGridx(); gPad->SetGridy(); cut.Divide(&trig); cut.Scale(100); cut.SetNameTitle("CutEff", "Background Supression: Cut efficiency (after star)"); cut.SetYTitle("\\eta [%]"); cut.SetDirectory(0); cut.SetMinimum(0); cut.SetMaximum(100); cut.SetBit(kCanDelete); cut.Draw(); TLine line; line.SetLineColor(kBlue); line.SetLineWidth(2); line.SetLineStyle(kDashed); line.DrawLine(cut.GetBinLowEdge(1), 50, cut.GetBinLowEdge(cut.GetNbinsX()+1), 50); } Bool_t MJSpectrum::Process(const MDataSet &set) { if (!set.IsValid()) { *fLog << err << "ERROR - DataSet invalid!" << endl; return kFALSE; } CheckEnv(); // -------------------------------------------------------------------------------- *fLog << inf; fLog->Separator(GetDescriptor()); *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl; *fLog << endl; MBinning bins1("BinningAlpha"); MBinning bins2("BinningEnergyEst"); MBinning bins3("BinningTheta"); MBinning bins4("BinningFalseSource"); MBinning bins5("BinningWidth"); MBinning bins6("BinningLength"); MBinning bins7("BinningDist"); MBinning bins8("BinningM3Long"); MBinning bins9("BinningM3Trans"); MBinning bins0("BinningSlope"); MBinning binsa("BinningAsym"); MBinning binsb("BinningConc1"); MAlphaFitter fit; MParList plist; plist.AddToList(&bins0); plist.AddToList(&bins1); plist.AddToList(&bins3); plist.AddToList(&bins4); plist.AddToList(&bins5); plist.AddToList(&bins6); plist.AddToList(&bins7); plist.AddToList(&bins8); plist.AddToList(&bins9); plist.AddToList(&binsa); plist.AddToList(&binsb); plist.AddToList(&fit); TH1D temp1, size; const Float_t ontime = ReadInput(plist, temp1, size); if (ontime<0) { *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl; return kFALSE; } plist.AddToList(&bins2); // Initialize weighting to a new spectru // m as defined in the resource file MMcSpectrumWeight weight; if (!InitWeighting(set, weight)) return kFALSE; bins3.SetEdges(temp1, 'x'); // Read the MC distribution as produced by corsika into // temp2 (1D) and apply the weights previously determined TH1D temp2(temp1); if (!ReadOrigMCDistribution(set, temp2, weight)) return kFALSE; // Calculate the weights according to the zenith angle distribution // of the raw-data which have to be applied to the MC events if (!GetThetaDistribution(temp1, temp2)) return kFALSE; // Tell the weighting task about the ZA-weights if (!fNoThetaWeights) weight.SetWeightsZd(&temp1); // Refill excess histogram to determin the excess events TH1D excess; if (!Refill(plist, excess)) return kFALSE; // Print the setup and result of the MAlphaFitter // Print used cuts PrintSetup(fit); TH2D hist; MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy"); if (fSimpleMode) { // Now we read the MC distribution as produced by corsika again // with the spectral weights applied into a histogram vs. // zenith angle and energy hist.UseCurrentStyle(); MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/); if (!ReadOrigMCDistribution(set, hist, weight)) return kFALSE; // No we apply the the zenith-angle-weights to the corsika produced // MC distribution. Unfortunately this must be done manually // because we are multiplying column by column if (!fRawMc) { for (int y=0; ySeparator("Calculate efficiencies"); *fLog << endl; MTaskList tlist2; plist.Replace(&tlist2); MReadMarsFile read("Events"); read.DisableAutoScheme(); if (!set.AddFilesOn(read)) return kFALSE; // 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 area0("TriggerArea"); MHCollectionArea area1; area0.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); area1.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist()); MHEnergyEst hest; MFillH fill30(&area0, "", "FillTriggerArea"); MFillH fill31(&area1, "", "FillCollectionArea"); MFillH fill4(&hest, "", "FillEnergyEst"); MFillH fill5("MHThreshold", "", "FillThreshold"); fill30.SetWeight(); fill31.SetWeight(); fill4.SetWeight(); fill5.SetWeight(); fill30.SetNameTab("TrigArea"); fill31.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(); MTaskEnv taskenv0("CalcHadronness"); taskenv0.SetDefault(fCalcHadronness); MEnergyEstimate est; MTaskEnv taskenv1("EstimateEnergy"); taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); tlist2.AddToList(&read); // If no weighting should be applied but the events should // be thrown away according to the theta distribution // it is enabled here if (!fRawMc && fNoThetaWeights) tlist2.AddToList(&contsel); //tlist2.AddToList(&calc); //tlist2.AddToList(&hcalc1); //tlist2.AddToList(&hcalc2); tlist2.AddToList(&weight); tlist2.AddToList(&fill1a); tlist2.AddToList(&fill30); tlist2.AddToList(fCutQ); tlist2.AddToList(fCut0); tlist2.AddToList(&taskenv0); tlist2.AddToList(fCutS); tlist2.AddToList(fCut1); tlist2.AddToList(fCut2); tlist2.AddToList(fCut3); tlist2.AddToList(&taskenv1); tlist2.AddToList(&fill31); 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); // by setting it here it is distributed to all consecutive tasks tlist2.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime); 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(area1, 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 *= area1.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/area1.GetEntries(); // Print normalization constant cout << "MC normalization factor: " << scale << endl; // Display cut efficiency DisplayCutEfficiency(area0, area1); // Overlay normalized plots DisplaySize(plist, scale); // check if output should be written if (fPathOut.IsNull()) return kTRUE; // Write the output TObjArray cont; cont.Add((TObject*)GetEnv()); if (fDisplay) cont.Add(fDisplay); return WriteContainer(cont, "", "RECREATE"); }