Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6977) +++ trunk/MagicSoft/Mars/Changelog (revision 6978) @@ -21,4 +21,49 @@ -*-*- END OF LINE -*-*- + + 2005/04/25 Thomas Bretz + + * ganymed.cc: + - changed policy of writing the resulting events to the result file + + * sponde.cc: + - added commandline option to use all monte carlos + - added command line option to read the MCs more accurate + + * sponde.rc: + - added + + * mbase/MStatusDisplay.[h,cc]: + - added some code to get Tab by name + - fixed a typo in a status line output + + * mhbase/MH.[h,cc], mhbase/MH3.[h,cc], mhflux/MHFalseSource.h, + mhist/MHCamEvent.[h,cc], mhist/MHCamEventRot.h, + mhist/MHEvent.h, mhist/MHStarMap.h, mhist/MHTriggerLvl0.[h,cc], + mhistmc/MHMcTriggerLvl2.[h,cc], mhvstime/MHPixVsTime.[h,cc], + mhvstime/MHSectorVsTime.[h,cc], mimage/MHHillas.[h,cc], + mimage/MHHillasExt.[h,cc], mimage/MHHillasSrc.[h,cc], + mimage/MHImagePar.[h,cc], mimage/MHNewImagePar.[h,cc]: + - changed GetHistByName to be const-qualified to be compatible + with FindObject + - added some FindObject function to call GetHistByName + + * mhflux/MHAlpha.[h,cc]: + - changed such, that it can be forced to display the excess + events versus size + + * mjobs/MJCut.[h,cc]: + - display number of excess events versus Size per default + - removed energy estimator + + * mjobs/MJOptimize.cc: + - display number of excess events verss size after optimization + + * mjobs/MJSpectrum.[h,cc]: + - implemented setting up energy estimator + - replaced some gLog by fLog + - display comparison of image parameters + + 2005/04/22 Thomas Bretz Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 6977) +++ trunk/MagicSoft/Mars/NEWS (revision 6978) @@ -21,4 +21,8 @@ - Changed the default paths for calibrated data and image files. (The implemented access to these files doesn't yet exist) + + - ganymed (MJCut and MJOptimize) now displayes the number of + excess events versus size. The energy estimation is done in + MJSpectrum (sponde) Index: trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc (revision 6977) +++ trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc (revision 6978) @@ -178,4 +178,61 @@ // -------------------------------------------------------------------------- +TGCompositeFrame *MStatusDisplay::GetTabContainer(const char *name) const +{ +#if ROOT_VERSION_CODE < ROOT_VERSION(4,03,05) + if (!fTab) + return 0; + + TGFrameElement *el; + TGTabElement *tab = 0; + TGCompositeFrame *comp = 0; + + TIter next(fTab->GetList()); + next(); // skip first container + + while ((el = (TGFrameElement *) next())) { + el = (TGFrameElement *) next(); + comp = (TGCompositeFrame *) el->fFrame; + next(); + tab = (TGTabElement *)el->fFrame; + if (name == tab->GetText()->GetString()) { + return comp; + } + } + + return 0; +#else + return fTab ? fTab->GetTabContainer(name) : 0; +#endif +} + +TGTabElement *MStatusDisplay::GetTabTab(const char *name) const +{ +#if ROOT_VERSION_CODE < ROOT_VERSION(4,03,05) + if (!fTab) + return 0; + + TGFrameElement *el; + TGTabElement *tab = 0; + + TIter next(fTab->GetList()); + next(); // skip first container + + while ((el = (TGFrameElement *) next())) { + next(); + tab = (TGTabElement *)el->fFrame; + if (name == tab->GetText()->GetString()) { + return tab; + } + } + + return 0; +#else + return fTab ? fTab->GetTabTab(name) : 0; +#endif +} +// -------------------------------------------------------------------------- + + // -------------------------------------------------------------------------- // @@ -1341,5 +1398,5 @@ case kLogAppend: - SetStatusLine1("Appending logg..."); + SetStatusLine1("Appending log..."); SetStatusLine2(""); *fLog << inf << "Appending log... " << flush; Index: trunk/MagicSoft/Mars/mbase/MStatusDisplay.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MStatusDisplay.h (revision 6977) +++ trunk/MagicSoft/Mars/mbase/MStatusDisplay.h (revision 6978) @@ -30,4 +30,5 @@ class TGTextView; class TGStatusBar; +class TGTabElement; class TGProgressBar; class TGHProgressBar; @@ -122,4 +123,7 @@ Bool_t HandleEvent(Event_t *event); + TGCompositeFrame *GetTabContainer(const char *name) const; + TGTabElement *GetTabTab(const char *name) const; + Bool_t HandleTimer(TTimer *timer=NULL); void UpdateTab(TGCompositeFrame *f); Index: trunk/MagicSoft/Mars/mhbase/MFillH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MFillH.cc (revision 6977) +++ trunk/MagicSoft/Mars/mhbase/MFillH.cc (revision 6978) @@ -543,13 +543,13 @@ */ - TVirtualPad *save = gPad; - if (fCanvas) - fCanvas->cd(); +// TVirtualPad *save = gPad; +// if (fCanvas) +// fCanvas->cd(); const Bool_t rc = fH->Fill(fParContainer, fWeight?fWeight->GetVal():1); fH->SetNumExecutions(GetNumExecutions()+1); - if (save && fCanvas) - save->cd(); +// if (save && fCanvas) +// save->cd(); return rc; } Index: trunk/MagicSoft/Mars/mhbase/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.cc (revision 6977) +++ trunk/MagicSoft/Mars/mhbase/MH.cc (revision 6978) @@ -117,5 +117,5 @@ // a pointer to a root histogram from the MH-derived class. // -TH1 *MH::GetHistByName(const TString name) +TH1 *MH::GetHistByName(const TString name) const { *fLog << warn << GetDescriptor() << ": GetHistByName not overloaded! Can't be used!" << endl; Index: trunk/MagicSoft/Mars/mhbase/MH.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.h (revision 6977) +++ trunk/MagicSoft/Mars/mhbase/MH.h (revision 6978) @@ -49,5 +49,5 @@ virtual TString GetDataMember() const { return ""; } - virtual TH1 *GetHistByName(const TString name); + virtual TH1 *GetHistByName(const TString name) const; static TCanvas *MakeDefCanvas(TString name="", const char *title="", Index: trunk/MagicSoft/Mars/mhbase/MH3.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH3.h (revision 6977) +++ trunk/MagicSoft/Mars/mhbase/MH3.h (revision 6978) @@ -72,5 +72,10 @@ const TH1 &GetHist() const { return *fHist; } - TH1 *GetHistByName(const TString name="") { return fHist; } + TH1 *GetHistByName(const TString name="") const { return fHist; } + TObject *FindObject(const TObject *obj) const { return 0; } + TObject *FindObject(const char *name) const + { + return (TObject*)GetHistByName(name); + } void SetColors() const; Index: trunk/MagicSoft/Mars/mhflux/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 6977) +++ trunk/MagicSoft/Mars/mhflux/MHFalseSource.h (revision 6978) @@ -59,6 +59,4 @@ MHFalseSource(const char *name=NULL, const char *title=NULL); - TH1 *GetHistByName(const TString name) { return &fHist; } - void FitSignificance(Float_t sigint=10, Float_t sigmax=75, Float_t bgmin=45, Float_t bgmax=85, Byte_t polynom=2); //*MENU* void FitSignificanceStd() { FitSignificance(); } //*MENU* Index: trunk/MagicSoft/Mars/mjobs/MJCut.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6977) +++ trunk/MagicSoft/Mars/mjobs/MJCut.cc (revision 6978) @@ -48,5 +48,5 @@ #include "MPrint.h" #include "MContinue.h" -#include "MEnergyEstimate.h" +//#include "MEnergyEstimate.h" #include "MTaskEnv.h" #include "MSrcPosCalc.h" @@ -56,4 +56,6 @@ #include "../mhflux/MAlphaFitter.h" +#include "../mhflux/MHAlpha.h" + #include "MH3.h" #include "MBinning.h" @@ -71,10 +73,10 @@ // Default constructor. Set defaults for fStoreSummary, fStoreresult, // fWriteOnly, fIsWobble and fFullDisplay to kFALSE and initialize -// fEstimateEnergy and fCalcHadronness with NULL. +// /*fEstimateEnergy and*/ fCalcHadronness with NULL. // MJCut::MJCut(const char *name, const char *title) - : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE), + : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE), fIsWobble(kFALSE), fIsMonteCarlo(kFALSE), fFullDisplay(kFALSE), /*fSubstraction(kFALSE),*/ - fEstimateEnergy(0), fCalcHadronness(0) + /*fEstimateEnergy(0),*/ fCalcHadronness(0) { fName = name ? name : "MJCut"; @@ -88,6 +90,6 @@ MJCut::~MJCut() { - if (fEstimateEnergy) - delete fEstimateEnergy; + //if (fEstimateEnergy) + // delete fEstimateEnergy; if (fCalcHadronness) delete fCalcHadronness; @@ -132,4 +134,5 @@ // Setup a task estimating the energy. The given task is cloned. // +/* void MJCut::SetEnergyEstimator(const MTask *task) { @@ -138,4 +141,5 @@ fEstimateEnergy = task ? (MTask*)task->Clone() : 0; } +*/ // -------------------------------------------------------------------------- @@ -246,4 +250,5 @@ TObjArray arr; + // Save all MBinnings TIter Next(plist); TObject *o=0; @@ -252,6 +257,12 @@ arr.Add(o); + // Save also the result, not only the setup + const MHAlpha *halpha = (MHAlpha*)plist.FindObject("MHAlpha"); + if (halpha) + arr.Add((TObject*)(&halpha->GetAlphaFitter())); + const TString fname(GetOutputFile(num)); + // If requested, write to already open output file if (fNameResult.IsNull() && fStoreResult) { @@ -413,4 +424,8 @@ fit.SetScaleMode(MAlphaFitter::kNone); + MHAlpha halphaoff("MHAlphaOff"); + halphaoff.ForceUsingSize(); + plist.AddToList(&halphaoff); + MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha"); MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS"); @@ -452,10 +467,10 @@ SetupWriter(write1, "WriteAfterCut3"); - +/* MEnergyEstimate est; MTaskEnv taskenv1("EstimateEnergy"); taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); - + */ MTaskEnv taskenv2("CalcHadronness"); taskenv2.SetDefault(fCalcHadronness); @@ -495,5 +510,5 @@ if (fIsWobble) tlist2.AddToList(&hcalc2); - tlist2.AddToList(&taskenv1); + //tlist2.AddToList(&taskenv1); tlist2.AddToList(&taskenv2); tlist2.AddToList(&cont0); @@ -539,11 +554,10 @@ TObjArray cont; - cont.Add(&fit); cont.Add(&cont0); cont.Add(&cont1); cont.Add(&cont2); cont.Add(&cont3); - if (taskenv1.GetTask()) - cont.Add(taskenv1.GetTask()); + //if (taskenv1.GetTask()) + // cont.Add(taskenv1.GetTask()); if (taskenv2.GetTask()) cont.Add(taskenv2.GetTask()); @@ -641,4 +655,8 @@ } */ + MHAlpha halphaon("MHAlpha"); + halphaon.ForceUsingSize(); + plist.AddToList(&halphaon); + MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha"); MFillH ffs2("MHFalseSource", "MHillas", "FillFS"); Index: trunk/MagicSoft/Mars/mjobs/MJCut.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCut.h (revision 6977) +++ trunk/MagicSoft/Mars/mjobs/MJCut.h (revision 6978) @@ -26,5 +26,5 @@ TString fNameOutput; - MTask *fEstimateEnergy; + //MTask *fEstimateEnergy; MTask *fCalcHadronness; @@ -57,5 +57,5 @@ void SetNameOutFile(const char *name="") { fNameOutput=name; } - void SetEnergyEstimator(const MTask *task=0); + //void SetEnergyEstimator(const MTask *task=0); void SetHadronnessCalculator(const MTask *task=0); Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 6977) +++ trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 6978) @@ -914,6 +914,8 @@ histof.SkipHistTheta(); //histof.SkipHistEnergy(); - histon.InitMapping(&m, 0); - histof.InitMapping(&m, 0); + histon.ForceUsingSize(); + histof.ForceUsingSize(); + histon.InitMapping(&m, 1); + histof.InitMapping(&m, 1); if (filter && filter->InheritsFrom(MFMagicCuts::Class())) Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 6977) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 6978) @@ -68,4 +68,6 @@ #include "MFEventSelector2.h" #include "MFDataMember.h" +#include "MEnergyEstimate.h" +#include "MTaskEnv.h" #include "MFillH.h" #include "MHillasCalc.h" @@ -79,5 +81,5 @@ MJSpectrum::MJSpectrum(const char *name, const char *title) : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0), - fRefill(kFALSE), fSimpleMode(kTRUE) + fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE) { fName = name ? name : "MJSpectrum"; @@ -99,4 +101,15 @@ } +// -------------------------------------------------------------------------- +// +// 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) const { @@ -125,30 +138,61 @@ } -Bool_t MJSpectrum::ReadInput(const MParList &plist) -{ - const TString fname = fPathIn; - - *fLog << inf << "Reading from file: " << fname << endl; - - TFile file(fname, "READ"); +void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const +{ + fLog->Separator("Alpha Fitter"); + *fLog << all; + fit.Print(); + + fLog->Separator("Used Cuts"); + fCut0->Print(); + fCut1->Print(); + fCut2->Print(); + fCut3->Print(); + + //gLog.Separator("Energy Estimator"); + //fEstimateEnergy->Print(); +} + +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 " << fname << endl; - return kFALSE; - } + *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", "MHAlpha"); + if (!vstime || !size) + return -1; + + vstime->Copy(h1); + size->Copy(h2); + h1.SetDirectory(0); + h2.SetDirectory(0); + + if (fDisplay) + arr.DisplayIn(*fDisplay, "MHAlpha"); if (!ReadTask(fCut0, "Cut0")) - return kFALSE; + return -1; if (!ReadTask(fCut1, "Cut1")) - return kFALSE; + return -1; if (!ReadTask(fCut2, "Cut2")) - return kFALSE; + return -1; if (!ReadTask(fCut3, "Cut3")) - return kFALSE; - - if (!ReadTask(fEstimateEnergy, "EstimateEnergy")) - return kFALSE; - - TObjArray arr; + return -1; + + TObjArray arrread; TIter Next(plist); @@ -156,57 +200,10 @@ while ((o=Next())) if (o->InheritsFrom(MBinning::Class())) - arr.Add(o); - - arr.Add(plist.FindObject("MAlphaFitter")); - - return ReadContainer(arr); -} - -void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const -{ - gLog.Separator("Alpha Fitter"); - *fLog << all; - fit.Print(); - - gLog.Separator("Used Cuts"); - fCut0->Print(); - fCut1->Print(); - fCut2->Print(); - fCut3->Print(); - - gLog.Separator("Energy Estimator"); - fEstimateEnergy->Print(); -} - -Float_t MJSpectrum::ReadHistograms(TH1D &h1, TH1D &h2) const -{ - TFile file(fPathIn, "READ"); - if (!file.IsOpen()) - { - *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl; + arrread.Add(o); + + arrread.Add(plist.FindObject("MAlphaFitter")); + + if (!ReadContainer(arrread)) return -1; - } - - MStatusArray *arr = (MStatusArray*)file.Get("MStatusDisplay"); - if (!arr) - { - gLog << "MStatusDisplay not found in file... abort." << endl; - return -1; - } - - TH1D *vstime = (TH1D*)arr->FindObjectInCanvas("Theta", "TH1D", "OnTime"); - TH1D *excen = (TH1D*)arr->FindObjectInCanvas("ExcessEnergy", "TH1D", "MHAlpha"); - if (!vstime || !excen) - return -1; - - vstime->Copy(h1); - excen->Copy(h2); - h1.SetDirectory(0); - h2.SetDirectory(0); - - if (fDisplay) - arr->DisplayIn(*fDisplay, "MHAlpha"); - - delete arr; return vstime->Integral(); @@ -303,5 +300,5 @@ void MJSpectrum::DisplayResult(const TH2D &h2) const { - if (!fDisplay->CdCanvas("ZdDist")) + if (!fDisplay || !fDisplay->CdCanvas("ZdDist")) return; @@ -342,6 +339,10 @@ read.AddFile(fPathIn); - MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha"); - MFillH fill2("MHAlpha", "MHillasSrc", "FillAlpha"); + MEnergyEstimate est; + MTaskEnv taskenv1("EstimateEnergy"); + taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est); + + MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlphaOff"); + MFillH fill2("MHAlpha", "MHillasSrc", "FillAlphaOn"); MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData"); @@ -352,5 +353,5 @@ tlist.AddToList(&read); - tlist.AddToList(fEstimateEnergy); + tlist.AddToList(&taskenv1); tlist.AddToList(&f0); tlist.AddToList(&f1); @@ -389,4 +390,277 @@ halpha->GetHEnergy().Copy(h2); h2.SetDirectory(0); + + return kTRUE; +} + +Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const +{ + MTaskList tlist1; + plist.Replace(&tlist1); + + MReadMarsFile readmc("OriginalMC"); + //readmc.DisableAutoScheme(); + set.AddFilesOn(readmc); + 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); + + temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg"); + MH3 mh3mc(temp1); + + MFEventSelector2 sel1(mh3mc); + sel1.SetHistIsProbability(); + + fill0.SetFilter(&sel1); + + if (!fRawMc) + tlist1.AddToList(&sel1); + tlist1.AddToList(&fill0); + + MEvtLoop loop1(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; + } + + tlist1.PrintStatistics(); + + if (!loop1.GetDisplay()) + { + *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl; + return kFALSE; + } + + if (bins2 && bins3) + { + bins2->SetName("BinningEnergyEst"); + bins3->SetName("BinningTheta"); + } + return kTRUE; +} + +void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const +{ + TH1D collarea(area.GetHEnergy()); + TH1D spectrum(excess); + TH1D weights; + hest.GetWeights(weights); + + cout << "Effective On time: " << ontime << "s" << endl; + + spectrum.SetDirectory(NULL); + spectrum.SetBit(kCanDelete); + spectrum.Scale(1./ontime); + spectrum.Divide(&collarea); + spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)"); + spectrum.SetYTitle("N/sm^{2}"); + + 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(); + spectrum.DrawCopy(); + + c1.cd(3); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + weights.DrawCopy(); + + //spectrum.Divide(&weights); + spectrum.Multiply(&weights); + spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)"); + + for (int i=0; iSetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + spectrum.SetXTitle("E [GeV]"); + spectrum.SetYTitle("N/TeVsm^{2}"); + spectrum.DrawCopy(); + + TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4); + f.SetParameter(0, -2.87); + f.SetParameter(1, 1.9e-6); + f.SetLineColor(kGreen); + spectrum.Fit(&f, "NI", "", 55, 2e4); + f.DrawCopy("same"); + + /* + TString str; + str += "(1.68#pm0.15)10^{-7}"; + str += "(\\frac{E}{TeV})^{-2.59#pm0.06}"; + str += "\\frac{ph}{TeVm^{2}s}"; + + TLatex tex; + tex.DrawLatex(2e2, 7e-5, str); + */ +} + +Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const +{ + cout << name << endl; + 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 scale = fit ? fit->GetScaleFactor() : 1; + + gPad->SetBorderMode(0); + h2->SetLineColor(kBlack); + h3->SetLineColor(kBlue); + h2->Add(h1, -scale); + + h2->Scale(1./h2->Integral()); + h3->Scale(1./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; +} + +Bool_t MJSpectrum::DisplaySize(MParList &plist) const +{ + *fLog << inf << "Reading from file: " << fPathIn << endl; + + cout << "Opening..." << endl; + TFile file(fPathIn, "READ"); + if (!file.IsOpen()) + { + *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl; + return kFALSE; + } + + cout << "Reading..." << endl; + file.cd(); + MStatusArray arr; + if (arr.Read()<=0) + { + *fLog << "MStatusDisplay not found in file... abort." << endl; + return kFALSE; + } + + cout << "Searching..." << endl; + + TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "MHAlpha"); + if (!excess) + return kFALSE; + + cout << "Displaying..." << endl; + + // ------------------- 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("Number of excess events vs Size (data, mc/blue)"); + excess->SetStats(kFALSE); + excess->Scale(1./excess->Integral()); + excess->DrawCopy(); + + TObject *o=0; + if ((o=plist.FindObject("ExcessSize"))) + { + TH1F *histsel = (TH1F*)o->FindObject(""); + histsel->SetStats(kFALSE); + histsel->Scale(1./histsel->Integral()); + histsel->SetLineColor(kBlue); + histsel->SetBit(kCanDelete); + histsel->DrawCopy("same"); + } + + // -------------- Comparison of Image Parameters -------------- + c.cd(2); + PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost"); + + c.cd(3); + PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost"); + + c.cd(4); + PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost"); + + c.cd(5); + PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost"); + + c.cd(6); + PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost"); return kTRUE; @@ -432,16 +706,14 @@ plist.AddToList(&fit); - if (!ReadInput(plist)) - return kFALSE; + 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; + } PrintSetup(fit); - - TH1D temp1, excess; - const Float_t ontime = ReadHistograms(temp1, excess); - if (ontime<0) - { - *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl; - return kFALSE; - } + bins3.SetEdges(temp1, 'x'); TH1D temp2(temp1); @@ -452,5 +724,6 @@ return kFALSE; - if (fRefill && !Refill(plist, temp2)) + TH1D excess; + if (!Refill(plist, excess)) return kFALSE; @@ -460,69 +733,18 @@ { hist.UseCurrentStyle(); - MH::SetBinning(&hist, temp1.GetXaxis(), excess.GetXaxis()); + MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/); if (!ReadOrigMCDistribution(set, hist)) return kFALSE; - for (int y=0; yAddTab("Spectrum"); - c.Divide(2,2); - c.cd(1); - gPad->SetBorderMode(0); - gPad->SetLogx(); - gPad->SetLogy(); - gPad->SetGridx(); - gPad->SetGridy(); - collarea.DrawCopy(); - - c.cd(2); - gPad->SetBorderMode(0); - gPad->SetLogx(); - gPad->SetLogy(); - gPad->SetGridx(); - gPad->SetGridy(); - excess.DrawCopy(); - - c.cd(3); - gPad->SetBorderMode(0); - gPad->SetLogx(); - gPad->SetLogy(); - gPad->SetGridx(); - gPad->SetGridy(); - weights.DrawCopy(); - - excess.Divide(&weights); - //excess.Multiply(&weights); - excess.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)"); - - for (int i=0; iSetBorderMode(0); - gPad->SetLogx(); - gPad->SetLogy(); - gPad->SetGridx(); - gPad->SetGridy(); - excess.DrawCopy(); - - TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4); - f.SetParameter(0, -2.87); - f.SetParameter(1, 1.9e-6); - f.SetLineColor(kGreen); - excess.Fit(&f, "NI", "", 55, 2e4); - f.DrawCopy("same"); + DisplaySpectrum(area, excess, hest, ontime); + DisplaySize(plist); if (!fPathOut.IsNull()) fDisplay->SaveAsRoot(fPathOut); - /* - TString str; - str += "(1.68#pm0.15)10^{-7}"; - str += "(\\frac{E}{TeV})^{-2.59#pm0.06}"; - str += "\\frac{ph}{TeVm^{2}s}"; - - TLatex tex; - tex.DrawLatex(2e2, 7e-5, str); - */ - return kTRUE; } Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 6977) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 6978) @@ -14,5 +14,8 @@ class MParList; class MDataSet; +class MHEnergyEst; class MAlphaFitter; +class MStatusArray; +class MHCollectionArea; class MJSpectrum : public MJob @@ -27,13 +30,20 @@ Bool_t fRefill; Bool_t fSimpleMode; + Bool_t fRawMc; + // Read Input Bool_t ReadTask(MTask* &task, const char *name) const; - Bool_t ReadInput(const MParList &plist); - void PrintSetup(const MAlphaFitter &fit) const; - Float_t ReadHistograms(TH1D &h1, TH1D &h2) const; + Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size); Bool_t ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const; Bool_t GetThetaDistribution(TH1D &temp1, TH1D &temp2) const; + Bool_t Refill(MParList &plist, TH1D &h) const; + + // Display Output + void PrintSetup(const MAlphaFitter &fit) const; void DisplayResult(const TH2D &mh1) const; - Bool_t Refill(MParList &plist, TH1D &h) const; + Bool_t IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set) const; + void DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const; + Bool_t DisplaySize(MParList &plist) const; + Bool_t PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const; public: @@ -45,4 +55,7 @@ void EnableRefilling(Bool_t b=kTRUE) { fRefill=b; } void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; } + void EnableRawMc(Bool_t b=kTRUE) { fRawMc=b; } + + void SetEnergyEstimator(const MTask *task); ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum