Changeset 6978


Ignore:
Timestamp:
04/25/05 15:37:35 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r6976 r6978  
    2121
    2222                                                 -*-*- END OF LINE -*-*-
     23
     24 2005/04/25 Thomas Bretz
     25
     26   * ganymed.cc:
     27     - changed policy of writing the resulting events to the result file
     28
     29   * sponde.cc:
     30     - added commandline option to use all monte carlos
     31     - added command line option to read the MCs more accurate
     32
     33   * sponde.rc:
     34     - added
     35
     36   * mbase/MStatusDisplay.[h,cc]:
     37     - added some code to get Tab by name
     38     - fixed a typo in a status line output
     39
     40   * mhbase/MH.[h,cc], mhbase/MH3.[h,cc], mhflux/MHFalseSource.h,
     41     mhist/MHCamEvent.[h,cc], mhist/MHCamEventRot.h,
     42     mhist/MHEvent.h, mhist/MHStarMap.h, mhist/MHTriggerLvl0.[h,cc],
     43     mhistmc/MHMcTriggerLvl2.[h,cc], mhvstime/MHPixVsTime.[h,cc],
     44     mhvstime/MHSectorVsTime.[h,cc], mimage/MHHillas.[h,cc],
     45     mimage/MHHillasExt.[h,cc], mimage/MHHillasSrc.[h,cc],
     46     mimage/MHImagePar.[h,cc], mimage/MHNewImagePar.[h,cc]:
     47     - changed GetHistByName to be const-qualified to be compatible
     48       with FindObject
     49     - added some FindObject function to call GetHistByName
     50
     51   * mhflux/MHAlpha.[h,cc]:
     52     - changed such, that it can be forced to display the excess
     53       events versus size
     54
     55   * mjobs/MJCut.[h,cc]:
     56     - display number of excess events versus Size per default
     57     - removed energy estimator
     58
     59   * mjobs/MJOptimize.cc:
     60     - display number of excess events verss size after optimization
     61
     62   * mjobs/MJSpectrum.[h,cc]:
     63     - implemented setting up energy estimator
     64     - replaced some gLog by fLog
     65     - display comparison of image parameters
     66
     67
    2368
    2469 2005/04/22 Thomas Bretz
  • trunk/MagicSoft/Mars/NEWS

    r6964 r6978  
    2121   - Changed the default paths for calibrated data and image files.
    2222     (The implemented access to these files doesn't yet exist)
     23
     24   - ganymed (MJCut and MJOptimize) now displayes the number of
     25     excess events versus size. The energy estimation is done in
     26     MJSpectrum (sponde)
    2327
    2428
  • trunk/MagicSoft/Mars/mbase/MStatusDisplay.cc

    r6954 r6978  
    178178// --------------------------------------------------------------------------
    179179
     180TGCompositeFrame *MStatusDisplay::GetTabContainer(const char *name) const
     181{
     182#if ROOT_VERSION_CODE < ROOT_VERSION(4,03,05)
     183    if (!fTab)
     184        return 0;
     185
     186   TGFrameElement *el;
     187   TGTabElement *tab = 0;
     188   TGCompositeFrame *comp = 0;
     189
     190   TIter next(fTab->GetList());
     191   next();           // skip first container
     192
     193   while ((el = (TGFrameElement *) next())) {
     194      el = (TGFrameElement *) next();
     195      comp = (TGCompositeFrame *) el->fFrame;
     196      next();
     197      tab = (TGTabElement *)el->fFrame;
     198      if (name == tab->GetText()->GetString()) {
     199         return comp;
     200      }
     201   }
     202
     203   return 0;
     204#else
     205   return fTab ? fTab->GetTabContainer(name) : 0;
     206#endif
     207}
     208
     209TGTabElement *MStatusDisplay::GetTabTab(const char *name) const
     210{
     211#if ROOT_VERSION_CODE < ROOT_VERSION(4,03,05)
     212    if (!fTab)
     213        return 0;
     214
     215   TGFrameElement *el;
     216   TGTabElement *tab = 0;
     217
     218   TIter next(fTab->GetList());
     219   next();           // skip first container
     220
     221   while ((el = (TGFrameElement *) next())) {
     222      next();
     223      tab = (TGTabElement *)el->fFrame;
     224      if (name == tab->GetText()->GetString()) {
     225         return tab;
     226      }
     227   }
     228
     229   return 0;
     230#else
     231   return fTab ? fTab->GetTabTab(name) : 0;
     232#endif
     233}
     234// --------------------------------------------------------------------------
     235
     236
    180237// --------------------------------------------------------------------------
    181238//
     
    13411398
    13421399    case kLogAppend:
    1343         SetStatusLine1("Appending logg...");
     1400        SetStatusLine1("Appending log...");
    13441401        SetStatusLine2("");
    13451402        *fLog << inf << "Appending log... " << flush;
  • trunk/MagicSoft/Mars/mbase/MStatusDisplay.h

    r6938 r6978  
    3030class TGTextView;
    3131class TGStatusBar;
     32class TGTabElement;
    3233class TGProgressBar;
    3334class TGHProgressBar;
     
    122123    Bool_t HandleEvent(Event_t *event);
    123124
     125    TGCompositeFrame *GetTabContainer(const char *name) const;
     126    TGTabElement     *GetTabTab(const char *name) const;
     127
    124128    Bool_t HandleTimer(TTimer *timer=NULL);
    125129    void UpdateTab(TGCompositeFrame *f);
  • trunk/MagicSoft/Mars/mhbase/MFillH.cc

    r6947 r6978  
    543543     */
    544544
    545     TVirtualPad *save = gPad;
    546     if (fCanvas)
    547         fCanvas->cd();
     545//    TVirtualPad *save = gPad;
     546//    if (fCanvas)
     547//        fCanvas->cd();
    548548
    549549    const Bool_t rc = fH->Fill(fParContainer, fWeight?fWeight->GetVal():1);
    550550    fH->SetNumExecutions(GetNumExecutions()+1);
    551551
    552     if (save && fCanvas)
    553         save->cd();
     552//    if (save && fCanvas)
     553//        save->cd();
    554554    return rc;
    555555}
  • trunk/MagicSoft/Mars/mhbase/MH.cc

    r6949 r6978  
    117117// a pointer to a root histogram from the MH-derived class.
    118118//
    119 TH1 *MH::GetHistByName(const TString name)
     119TH1 *MH::GetHistByName(const TString name) const
    120120{
    121121    *fLog << warn << GetDescriptor() << ": GetHistByName not overloaded! Can't be used!" << endl;
  • trunk/MagicSoft/Mars/mhbase/MH.h

    r6948 r6978  
    4949    virtual TString GetDataMember() const { return ""; }
    5050
    51     virtual TH1 *GetHistByName(const TString name);
     51    virtual TH1 *GetHistByName(const TString name) const;
    5252
    5353    static TCanvas *MakeDefCanvas(TString name="", const char *title="",
  • trunk/MagicSoft/Mars/mhbase/MH3.h

    r5620 r6978  
    7272    const TH1 &GetHist() const { return *fHist; }
    7373
    74     TH1 *GetHistByName(const TString name="") { return fHist; }
     74    TH1 *GetHistByName(const TString name="") const { return fHist; }
     75    TObject *FindObject(const TObject *obj) const { return 0; }
     76    TObject *FindObject(const char *name) const
     77    {
     78        return (TObject*)GetHistByName(name);
     79    }
    7580
    7681    void SetColors() const;
  • trunk/MagicSoft/Mars/mhflux/MHFalseSource.h

    r5957 r6978  
    5959    MHFalseSource(const char *name=NULL, const char *title=NULL);
    6060
    61     TH1 *GetHistByName(const TString name) { return &fHist; }
    62 
    6361    void FitSignificance(Float_t sigint=10, Float_t sigmax=75, Float_t bgmin=45, Float_t bgmax=85, Byte_t polynom=2); //*MENU*
    6462    void FitSignificanceStd() { FitSignificance(); } //*MENU*
  • trunk/MagicSoft/Mars/mjobs/MJCut.cc

    r6958 r6978  
    4848#include "MPrint.h"
    4949#include "MContinue.h"
    50 #include "MEnergyEstimate.h"
     50//#include "MEnergyEstimate.h"
    5151#include "MTaskEnv.h"
    5252#include "MSrcPosCalc.h"
     
    5656
    5757#include "../mhflux/MAlphaFitter.h"
     58#include "../mhflux/MHAlpha.h"
     59
    5860#include "MH3.h"
    5961#include "MBinning.h"
     
    7173// Default constructor.  Set defaults for fStoreSummary, fStoreresult,
    7274// fWriteOnly, fIsWobble and fFullDisplay to kFALSE and initialize
    73 // fEstimateEnergy and fCalcHadronness with NULL.
     75// /*fEstimateEnergy and*/ fCalcHadronness with NULL.
    7476//
    7577MJCut::MJCut(const char *name, const char *title)
    76     : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE),
     78    : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
    7779    fIsWobble(kFALSE), fIsMonteCarlo(kFALSE),  fFullDisplay(kFALSE), /*fSubstraction(kFALSE),*/
    78     fEstimateEnergy(0), fCalcHadronness(0)
     80    /*fEstimateEnergy(0),*/ fCalcHadronness(0)
    7981{
    8082    fName  = name  ? name  : "MJCut";
     
    8890MJCut::~MJCut()
    8991{
    90     if (fEstimateEnergy)
    91         delete fEstimateEnergy;
     92    //if (fEstimateEnergy)
     93    //    delete fEstimateEnergy;
    9294    if (fCalcHadronness)
    9395        delete fCalcHadronness;
     
    132134// Setup a task estimating the energy. The given task is cloned.
    133135//
     136/*
    134137void MJCut::SetEnergyEstimator(const MTask *task)
    135138{
     
    138141    fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
    139142}
     143*/
    140144
    141145// --------------------------------------------------------------------------
     
    246250    TObjArray arr;
    247251
     252    // Save all MBinnings
    248253    TIter Next(plist);
    249254    TObject *o=0;
     
    252257            arr.Add(o);
    253258
     259    // Save also the result, not only the setup
     260    const MHAlpha *halpha = (MHAlpha*)plist.FindObject("MHAlpha");
     261    if (halpha)
     262        arr.Add((TObject*)(&halpha->GetAlphaFitter()));
     263
    254264    const TString fname(GetOutputFile(num));
    255265
     266    // If requested, write to already open output file
    256267    if (fNameResult.IsNull() && fStoreResult)
    257268    {
     
    413424        fit.SetScaleMode(MAlphaFitter::kNone);
    414425
     426    MHAlpha halphaoff("MHAlphaOff");
     427    halphaoff.ForceUsingSize();
     428    plist.AddToList(&halphaoff);
     429
    415430    MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
    416431    MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
     
    452467    SetupWriter(write1, "WriteAfterCut3");
    453468
    454 
     469/*
    455470    MEnergyEstimate est;
    456471
    457472    MTaskEnv taskenv1("EstimateEnergy");
    458473    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
    459 
     474  */
    460475    MTaskEnv taskenv2("CalcHadronness");
    461476    taskenv2.SetDefault(fCalcHadronness);
     
    495510    if (fIsWobble)
    496511        tlist2.AddToList(&hcalc2);
    497     tlist2.AddToList(&taskenv1);
     512    //tlist2.AddToList(&taskenv1);
    498513    tlist2.AddToList(&taskenv2);
    499514    tlist2.AddToList(&cont0);
     
    539554
    540555    TObjArray cont;
    541     cont.Add(&fit);
    542556    cont.Add(&cont0);
    543557    cont.Add(&cont1);
    544558    cont.Add(&cont2);
    545559    cont.Add(&cont3);
    546     if (taskenv1.GetTask())
    547         cont.Add(taskenv1.GetTask());
     560    //if (taskenv1.GetTask())
     561    //    cont.Add(taskenv1.GetTask());
    548562    if (taskenv2.GetTask())
    549563        cont.Add(taskenv2.GetTask());
     
    641655    }
    642656    */
     657    MHAlpha halphaon("MHAlpha");
     658    halphaon.ForceUsingSize();
     659    plist.AddToList(&halphaon);
     660
    643661    MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
    644662    MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
  • trunk/MagicSoft/Mars/mjobs/MJCut.h

    r6949 r6978  
    2626    TString fNameOutput;
    2727
    28     MTask *fEstimateEnergy;
     28    //MTask *fEstimateEnergy;
    2929    MTask *fCalcHadronness;
    3030
     
    5757    void SetNameOutFile(const char *name="")      { fNameOutput=name; }
    5858
    59     void SetEnergyEstimator(const MTask *task=0);
     59    //void SetEnergyEstimator(const MTask *task=0);
    6060    void SetHadronnessCalculator(const MTask *task=0);
    6161
  • trunk/MagicSoft/Mars/mjobs/MJOptimize.cc

    r6965 r6978  
    914914    histof.SkipHistTheta();
    915915    //histof.SkipHistEnergy();
    916     histon.InitMapping(&m, 0);
    917     histof.InitMapping(&m, 0);
     916    histon.ForceUsingSize();
     917    histof.ForceUsingSize();
     918    histon.InitMapping(&m, 1);
     919    histof.InitMapping(&m, 1);
    918920
    919921    if (filter && filter->InheritsFrom(MFMagicCuts::Class()))
  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

    r6966 r6978  
    6868#include "MFEventSelector2.h"
    6969#include "MFDataMember.h"
     70#include "MEnergyEstimate.h"
     71#include "MTaskEnv.h"
    7072#include "MFillH.h"
    7173#include "MHillasCalc.h"
     
    7981MJSpectrum::MJSpectrum(const char *name, const char *title)
    8082    : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
    81     fRefill(kFALSE), fSimpleMode(kTRUE)
     83    fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
    8284{
    8385    fName  = name  ? name  : "MJSpectrum";
     
    99101}
    100102
     103// --------------------------------------------------------------------------
     104//
     105// Setup a task estimating the energy. The given task is cloned.
     106//
     107void MJSpectrum::SetEnergyEstimator(const MTask *task)
     108{
     109    if (fEstimateEnergy)
     110        delete fEstimateEnergy;
     111    fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
     112}
     113
    101114Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
    102115{
     
    125138}
    126139
    127 Bool_t MJSpectrum::ReadInput(const MParList &plist)
    128 {
    129     const TString fname = fPathIn;
    130 
    131     *fLog << inf << "Reading from file: " << fname << endl;
    132 
    133     TFile file(fname, "READ");
     140void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
     141{
     142    fLog->Separator("Alpha Fitter");
     143    *fLog << all;
     144    fit.Print();
     145
     146    fLog->Separator("Used Cuts");
     147    fCut0->Print();
     148    fCut1->Print();
     149    fCut2->Print();
     150    fCut3->Print();
     151
     152    //gLog.Separator("Energy Estimator");
     153    //fEstimateEnergy->Print();
     154}
     155
     156Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
     157{
     158    *fLog << inf << "Reading from file: " << fPathIn << endl;
     159
     160    TFile file(fPathIn, "READ");
    134161    if (!file.IsOpen())
    135162    {
    136         *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
    137         return kFALSE;
    138     }
     163        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
     164        return -1;
     165    }
     166
     167    MStatusArray arr;
     168    if (arr.Read()<=0)
     169    {
     170        *fLog << "MStatusDisplay not found in file... abort." << endl;
     171        return -1;
     172    }
     173
     174    TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",  "TH1D", "OnTime");
     175    TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "MHAlpha");
     176    if (!vstime || !size)
     177        return -1;
     178
     179    vstime->Copy(h1);
     180    size->Copy(h2);
     181    h1.SetDirectory(0);
     182    h2.SetDirectory(0);
     183
     184    if (fDisplay)
     185        arr.DisplayIn(*fDisplay, "MHAlpha");
    139186
    140187    if (!ReadTask(fCut0, "Cut0"))
    141         return kFALSE;
     188        return -1;
    142189    if (!ReadTask(fCut1, "Cut1"))
    143         return kFALSE;
     190        return -1;
    144191    if (!ReadTask(fCut2, "Cut2"))
    145         return kFALSE;
     192        return -1;
    146193    if (!ReadTask(fCut3, "Cut3"))
    147         return kFALSE;
    148 
    149     if (!ReadTask(fEstimateEnergy, "EstimateEnergy"))
    150         return kFALSE;
    151 
    152     TObjArray arr;
     194        return -1;
     195
     196    TObjArray arrread;
    153197
    154198    TIter Next(plist);
     
    156200    while ((o=Next()))
    157201        if (o->InheritsFrom(MBinning::Class()))
    158             arr.Add(o);
    159 
    160     arr.Add(plist.FindObject("MAlphaFitter"));
    161 
    162     return ReadContainer(arr);
    163 }
    164 
    165 void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
    166 {
    167     gLog.Separator("Alpha Fitter");
    168     *fLog << all;
    169     fit.Print();
    170 
    171     gLog.Separator("Used Cuts");
    172     fCut0->Print();
    173     fCut1->Print();
    174     fCut2->Print();
    175     fCut3->Print();
    176 
    177     gLog.Separator("Energy Estimator");
    178     fEstimateEnergy->Print();
    179 }
    180 
    181 Float_t MJSpectrum::ReadHistograms(TH1D &h1, TH1D &h2) const
    182 {
    183     TFile file(fPathIn, "READ");
    184     if (!file.IsOpen())
    185     {
    186         *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
     202            arrread.Add(o);
     203
     204    arrread.Add(plist.FindObject("MAlphaFitter"));
     205
     206    if (!ReadContainer(arrread))
    187207        return -1;
    188     }
    189 
    190     MStatusArray *arr = (MStatusArray*)file.Get("MStatusDisplay");
    191     if (!arr)
    192     {
    193         gLog << "MStatusDisplay not found in file... abort." << endl;
    194         return -1;
    195     }
    196 
    197     TH1D *vstime = (TH1D*)arr->FindObjectInCanvas("Theta",        "TH1D", "OnTime");
    198     TH1D *excen  = (TH1D*)arr->FindObjectInCanvas("ExcessEnergy", "TH1D", "MHAlpha");
    199     if (!vstime || !excen)
    200         return -1;
    201 
    202     vstime->Copy(h1);
    203     excen->Copy(h2);
    204     h1.SetDirectory(0);
    205     h2.SetDirectory(0);
    206 
    207     if (fDisplay)
    208         arr->DisplayIn(*fDisplay, "MHAlpha");
    209 
    210     delete arr;
    211208
    212209    return vstime->Integral();
     
    303300void MJSpectrum::DisplayResult(const TH2D &h2) const
    304301{
    305     if (!fDisplay->CdCanvas("ZdDist"))
     302    if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
    306303        return;
    307304
     
    342339    read.AddFile(fPathIn);
    343340
    344     MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
    345     MFillH fill2("MHAlpha", "MHillasSrc", "FillAlpha");
     341    MEnergyEstimate est;
     342    MTaskEnv taskenv1("EstimateEnergy");
     343    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
     344
     345    MFillH fill1("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlphaOff");
     346    MFillH fill2("MHAlpha",              "MHillasSrc", "FillAlphaOn");
    346347
    347348    MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
     
    352353
    353354    tlist.AddToList(&read);
    354     tlist.AddToList(fEstimateEnergy);
     355    tlist.AddToList(&taskenv1);
    355356    tlist.AddToList(&f0);
    356357    tlist.AddToList(&f1);
     
    389390    halpha->GetHEnergy().Copy(h2);
    390391    h2.SetDirectory(0);
     392
     393    return kTRUE;
     394}
     395
     396Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
     397{
     398    MTaskList tlist1;
     399    plist.Replace(&tlist1);
     400
     401    MReadMarsFile readmc("OriginalMC");
     402    //readmc.DisableAutoScheme();
     403    set.AddFilesOn(readmc);
     404    readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
     405    readmc.EnableBranch("MMcEvtBasic.fEnergy");
     406
     407    mh1.SetLogy();
     408    mh1.SetLogz();
     409    mh1.SetName("ThetaE");
     410
     411    MFillH fill0(&mh1);
     412    //fill0.SetDrawOption("projx only");
     413
     414    MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
     415    MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
     416    if (bins2 && bins3)
     417    {
     418        bins2->SetName("BinningThetaEY");
     419        bins3->SetName("BinningThetaEX");
     420    }
     421    tlist1.AddToList(&readmc);
     422
     423    temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
     424    MH3 mh3mc(temp1);
     425
     426    MFEventSelector2 sel1(mh3mc);
     427    sel1.SetHistIsProbability();
     428
     429    fill0.SetFilter(&sel1);
     430
     431    if (!fRawMc)
     432        tlist1.AddToList(&sel1);
     433    tlist1.AddToList(&fill0);
     434
     435    MEvtLoop loop1(fName);
     436    loop1.SetParList(&plist);
     437    loop1.SetLogStream(fLog);
     438    loop1.SetDisplay(fDisplay);
     439
     440    if (!SetupEnv(loop1))
     441        return kFALSE;
     442
     443    if (!loop1.Eventloop(fMaxEvents))
     444    {
     445        *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
     446        return kFALSE;
     447    }
     448
     449    tlist1.PrintStatistics();
     450
     451    if (!loop1.GetDisplay())
     452    {
     453        *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
     454        return kFALSE;
     455    }
     456
     457    if (bins2 && bins3)
     458    {
     459        bins2->SetName("BinningEnergyEst");
     460        bins3->SetName("BinningTheta");
     461    }
     462    return kTRUE;
     463}
     464
     465void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
     466{
     467    TH1D collarea(area.GetHEnergy());
     468    TH1D spectrum(excess);
     469    TH1D weights;
     470    hest.GetWeights(weights);
     471
     472    cout << "Effective On time: " << ontime << "s" << endl;
     473
     474    spectrum.SetDirectory(NULL);
     475    spectrum.SetBit(kCanDelete);
     476    spectrum.Scale(1./ontime);
     477    spectrum.Divide(&collarea);
     478    spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
     479    spectrum.SetYTitle("N/sm^{2}");
     480
     481    TCanvas &c1 = fDisplay->AddTab("Spectrum");
     482    c1.Divide(2,2);
     483    c1.cd(1);
     484    gPad->SetBorderMode(0);
     485    gPad->SetLogx();
     486    gPad->SetLogy();
     487    gPad->SetGridx();
     488    gPad->SetGridy();
     489    collarea.DrawCopy();
     490
     491    c1.cd(2);
     492    gPad->SetBorderMode(0);
     493    gPad->SetLogx();
     494    gPad->SetLogy();
     495    gPad->SetGridx();
     496    gPad->SetGridy();
     497    spectrum.DrawCopy();
     498
     499    c1.cd(3);
     500    gPad->SetBorderMode(0);
     501    gPad->SetLogx();
     502    gPad->SetLogy();
     503    gPad->SetGridx();
     504    gPad->SetGridy();
     505    weights.DrawCopy();
     506
     507    //spectrum.Divide(&weights);
     508    spectrum.Multiply(&weights);
     509    spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
     510
     511    for (int i=0; i<excess.GetNbinsX(); i++)
     512    {
     513        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
     514        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
     515    }
     516
     517    c1.cd(4);
     518    gPad->SetBorderMode(0);
     519    gPad->SetLogx();
     520    gPad->SetLogy();
     521    gPad->SetGridx();
     522    gPad->SetGridy();
     523    spectrum.SetXTitle("E [GeV]");
     524    spectrum.SetYTitle("N/TeVsm^{2}");
     525    spectrum.DrawCopy();
     526
     527    TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
     528    f.SetParameter(0, -2.87);
     529    f.SetParameter(1, 1.9e-6);
     530    f.SetLineColor(kGreen);
     531    spectrum.Fit(&f, "NI", "", 55, 2e4);
     532    f.DrawCopy("same");
     533
     534    /*
     535     TString str;
     536     str += "(1.68#pm0.15)10^{-7}";
     537     str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
     538     str += "\\frac{ph}{TeVm^{2}s}";
     539
     540     TLatex tex;
     541     tex.DrawLatex(2e2, 7e-5, str);
     542     */
     543}
     544
     545Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
     546{
     547    cout << name << endl;
     548    TString same(name);
     549    same += "Same";
     550
     551    TH1 *h1  = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
     552    TH1 *h2  = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
     553    if (!h1 || !h2)
     554        return kFALSE;
     555
     556    TObject *obj = plist.FindObject(plot);
     557    if (!obj)
     558    {
     559        *fLog << warn << plot << " not in parameter list... skipping." << endl;
     560        return kFALSE;
     561    }
     562
     563    TH1 *h3  = (TH1*)obj->FindObject(name);
     564    if (!h3)
     565    {
     566        *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
     567        return kFALSE;
     568    }
     569
     570
     571    const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
     572    const Double_t scale = fit ? fit->GetScaleFactor() : 1;
     573
     574    gPad->SetBorderMode(0);
     575    h2->SetLineColor(kBlack);
     576    h3->SetLineColor(kBlue);
     577    h2->Add(h1, -scale);
     578
     579    h2->Scale(1./h2->Integral());
     580    h3->Scale(1./h3->Integral());
     581
     582    h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
     583
     584    h2 = h2->DrawCopy();
     585    h3 = h3->DrawCopy("same");
     586
     587    // Don't do this on the original object!
     588    h2->SetStats(kFALSE);
     589    h3->SetStats(kFALSE);
     590
     591    return kTRUE;
     592}
     593
     594Bool_t MJSpectrum::DisplaySize(MParList &plist) const
     595{
     596    *fLog << inf << "Reading from file: " << fPathIn << endl;
     597
     598    cout << "Opening..." << endl;
     599    TFile file(fPathIn, "READ");
     600    if (!file.IsOpen())
     601    {
     602        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
     603        return kFALSE;
     604    }
     605
     606    cout << "Reading..." << endl;
     607    file.cd();
     608    MStatusArray arr;
     609    if (arr.Read()<=0)
     610    {
     611        *fLog << "MStatusDisplay not found in file... abort." << endl;
     612        return kFALSE;
     613    }
     614
     615    cout << "Searching..." << endl;
     616
     617    TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "MHAlpha");
     618    if (!excess)
     619        return kFALSE;
     620
     621    cout << "Displaying..." << endl;
     622
     623    // ------------------- Plot excess versus size -------------------
     624
     625    TCanvas &c = fDisplay->AddTab("Excess");
     626    c.Divide(3,2);
     627    c.cd(1);
     628    gPad->SetBorderMode(0);
     629    gPad->SetLogx();
     630    gPad->SetLogy();
     631    gPad->SetGridx();
     632    gPad->SetGridy();
     633
     634    excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
     635    excess->SetStats(kFALSE);
     636    excess->Scale(1./excess->Integral());
     637    excess->DrawCopy();
     638
     639    TObject *o=0;
     640    if ((o=plist.FindObject("ExcessSize")))
     641    {
     642        TH1F *histsel = (TH1F*)o->FindObject("");
     643        histsel->SetStats(kFALSE);
     644        histsel->Scale(1./histsel->Integral());
     645        histsel->SetLineColor(kBlue);
     646        histsel->SetBit(kCanDelete);
     647        histsel->DrawCopy("same");
     648    }
     649
     650    // -------------- Comparison of Image Parameters --------------
     651    c.cd(2);
     652    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost");
     653
     654    c.cd(3);
     655    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
     656
     657    c.cd(4);
     658    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost");
     659
     660    c.cd(5);
     661    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost");
     662
     663    c.cd(6);
     664    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost");
    391665
    392666    return kTRUE;
     
    432706    plist.AddToList(&fit);
    433707
    434     if (!ReadInput(plist))
    435         return kFALSE;
     708    TH1D temp1, size;
     709    const Float_t ontime = ReadInput(plist, temp1, size);
     710    if (ontime<0)
     711    {
     712        *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
     713        return kFALSE;
     714    }
    436715
    437716    PrintSetup(fit);
    438 
    439     TH1D temp1, excess;
    440     const Float_t ontime = ReadHistograms(temp1, excess);
    441     if (ontime<0)
    442     {
    443         *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
    444         return kFALSE;
    445     }
     717    bins3.SetEdges(temp1, 'x');
    446718
    447719    TH1D temp2(temp1);
     
    452724        return kFALSE;
    453725
    454     if (fRefill && !Refill(plist, temp2))
     726    TH1D excess;
     727    if (!Refill(plist, excess))
    455728        return kFALSE;
    456729
     
    460733    {
    461734        hist.UseCurrentStyle();
    462         MH::SetBinning(&hist, temp1.GetXaxis(), excess.GetXaxis());
     735        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
    463736        if (!ReadOrigMCDistribution(set, hist))
    464737            return kFALSE;
    465738
    466         for (int y=0; y<hist.GetNbinsY(); y++)
    467             for (int x=0; x<hist.GetNbinsX(); x++)
    468                 hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
     739        if (!fRawMc)
     740        {
     741            for (int y=0; y<hist.GetNbinsY(); y++)
     742                for (int x=0; x<hist.GetNbinsX(); x++)
     743                    hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
     744        }
    469745    }
    470746    else
    471     {
    472         MTaskList tlist1;
    473         plist.Replace(&tlist1);
    474 
    475         MReadMarsFile readmc("OriginalMC");
    476         //readmc.DisableAutoScheme();
    477         set.AddFilesOn(readmc);
    478         readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
    479         readmc.EnableBranch("MMcEvtBasic.fEnergy");
    480 
    481         mh1.SetLogy();
    482         mh1.SetLogz();
    483         mh1.SetName("ThetaE");
    484 
    485         MFillH fill0(&mh1);
    486         //fill0.SetDrawOption("projx only");
    487 
    488         MBinning binsx(bins3, "BinningThetaEX");
    489         MBinning binsy(bins2, "BinningThetaEY");
    490         plist.AddToList(&binsx);
    491         plist.AddToList(&binsy);
    492         tlist1.AddToList(&readmc);
    493 
    494         temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
    495         MH3 mh3mc(temp1);
    496 
    497         MFEventSelector2 sel1(mh3mc);
    498         sel1.SetHistIsProbability();
    499 
    500         fill0.SetFilter(&sel1);
    501 
    502         tlist1.AddToList(&sel1);
    503         tlist1.AddToList(&fill0);
    504 
    505         MEvtLoop loop1(fName);
    506         loop1.SetParList(&plist);
    507         loop1.SetLogStream(fLog);
    508         loop1.SetDisplay(fDisplay);
    509 
    510         if (!SetupEnv(loop1))
     747        if (!IntermediateLoop(plist, mh1, temp1, set))
    511748            return kFALSE;
    512 
    513         if (!loop1.Eventloop(fMaxEvents))
    514         {
    515             *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
    516             return kFALSE;
    517         }
    518 
    519         tlist1.PrintStatistics();
    520 
    521         if (!loop1.GetDisplay())
    522         {
    523             *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
    524             return kFALSE;
    525         }
    526     }
    527749
    528750    DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
     
    572794    MFillH fill4(&hest, "", "FillEnergyEst");
    573795
     796    MH3 hsize("MHillas.fSize");
     797    //MH3 henergy("MEnergyEst.fVal");
     798    hsize.SetName("ExcessSize");
     799    //henergy.SetName("EnergyEst");
     800    MBinning bins(size, "BinningExcessSize");
     801    plist.AddToList(&hsize);
     802    //plist.AddToList(&henergy);
     803    plist.AddToList(&bins);
     804
    574805    MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
    575806    MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
     
    579810    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
    580811    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
     812    MFillH fill8a("ExcessSize     [MH3]",           "",             "FillExcessSize");
     813    //MFillH fill9a("EnergyEst      [MH3]",           "",             "FillExcessEEst");
    581814    fill1a.SetNameTab("PreCut");
    582815    fill2a.SetNameTab("PostCut");
     
    586819    fill6a.SetNameTab("ImgPar");
    587820    fill7a.SetNameTab("NewPar");
     821    fill8a.SetBit(MFillH::kDoNotDisplay);
     822    //fill9a.SetBit(MFillH::kDoNotDisplay);
     823
     824    MEnergyEstimate est;
     825    MTaskEnv taskenv1("EstimateEnergy");
     826    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
    588827
    589828    tlist2.AddToList(&read);
    590     tlist2.AddToList(&contsel);
     829    if (!fRawMc)
     830        tlist2.AddToList(&contsel);
    591831    tlist2.AddToList(&calc);
    592832    tlist2.AddToList(&hcalc1);
     
    597837    tlist2.AddToList(fCut2);
    598838    tlist2.AddToList(fCut3);
    599     tlist2.AddToList(fEstimateEnergy);
     839    tlist2.AddToList(&taskenv1);
    600840    tlist2.AddToList(&fill3);
    601841    tlist2.AddToList(&fill4);
     
    606846    tlist2.AddToList(&fill6a);
    607847    tlist2.AddToList(&fill7a);
     848    tlist2.AddToList(&fill8a);
     849    //tlist2.AddToList(&fill9a);
    608850
    609851    MEvtLoop loop2(fName);
     
    631873    // -------------------------- Spectrum ----------------------------
    632874
    633     TH1D collarea(area.GetHEnergy());
    634     TH1D weights;
    635     hest.GetWeights(weights);
    636 
    637     cout << "Effective On time: " << ontime << "s" << endl;
    638 
    639     excess.SetDirectory(NULL);
    640     excess.SetBit(kCanDelete);
    641     excess.Scale(1./ontime);
    642     excess.Divide(&collarea);
    643     excess.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
    644     excess.SetYTitle("N/sm^{2}");
    645 
    646     TCanvas &c = fDisplay->AddTab("Spectrum");
    647     c.Divide(2,2);
    648     c.cd(1);
    649     gPad->SetBorderMode(0);
    650     gPad->SetLogx();
    651     gPad->SetLogy();
    652     gPad->SetGridx();
    653     gPad->SetGridy();
    654     collarea.DrawCopy();
    655 
    656     c.cd(2);
    657     gPad->SetBorderMode(0);
    658     gPad->SetLogx();
    659     gPad->SetLogy();
    660     gPad->SetGridx();
    661     gPad->SetGridy();
    662     excess.DrawCopy();
    663 
    664     c.cd(3);
    665     gPad->SetBorderMode(0);
    666     gPad->SetLogx();
    667     gPad->SetLogy();
    668     gPad->SetGridx();
    669     gPad->SetGridy();
    670     weights.DrawCopy();
    671 
    672     excess.Divide(&weights);
    673     //excess.Multiply(&weights);
    674     excess.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
    675 
    676     for (int i=0; i<excess.GetNbinsX(); i++)
    677     {
    678         excess.SetBinContent(i+1, excess.GetBinContent(i+1)/excess.GetBinWidth(i+1)*1000);
    679         excess.SetBinError(i+1,   excess.GetBinError(i+1)/  excess.GetBinWidth(i+1)*1000);
    680     }
    681 
    682     excess.SetYTitle("N/TeVsm^{2}");
    683 
    684     c.cd(4);
    685     gPad->SetBorderMode(0);
    686     gPad->SetLogx();
    687     gPad->SetLogy();
    688     gPad->SetGridx();
    689     gPad->SetGridy();
    690     excess.DrawCopy();
    691 
    692     TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
    693     f.SetParameter(0, -2.87);
    694     f.SetParameter(1, 1.9e-6);
    695     f.SetLineColor(kGreen);
    696     excess.Fit(&f, "NI", "", 55, 2e4);
    697     f.DrawCopy("same");
     875    DisplaySpectrum(area, excess, hest, ontime);
     876    DisplaySize(plist);
    698877
    699878    if (!fPathOut.IsNull())
    700879        fDisplay->SaveAsRoot(fPathOut);
    701880
    702     /*
    703      TString str;
    704      str += "(1.68#pm0.15)10^{-7}";
    705      str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
    706      str += "\\frac{ph}{TeVm^{2}s}";
    707 
    708      TLatex tex;
    709      tex.DrawLatex(2e2, 7e-5, str);
    710      */
    711 
    712881    return kTRUE;
    713882}
  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.h

    r6966 r6978  
    1414class MParList;
    1515class MDataSet;
     16class MHEnergyEst;
    1617class MAlphaFitter;
     18class MStatusArray;
     19class MHCollectionArea;
    1720
    1821class MJSpectrum : public MJob
     
    2730    Bool_t fRefill;
    2831    Bool_t fSimpleMode;
     32    Bool_t fRawMc;
    2933
     34    // Read Input
    3035    Bool_t  ReadTask(MTask* &task, const char *name) const;
    31     Bool_t  ReadInput(const MParList &plist);
    32     void    PrintSetup(const MAlphaFitter &fit) const;
    33     Float_t ReadHistograms(TH1D &h1, TH1D &h2) const;
     36    Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size);
    3437    Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const;
    3538    Bool_t  GetThetaDistribution(TH1D &temp1, TH1D &temp2) const;
     39    Bool_t  Refill(MParList &plist, TH1D &h) const;
     40
     41    // Display Output
     42    void    PrintSetup(const MAlphaFitter &fit) const;
    3643    void    DisplayResult(const TH2D &mh1) const;
    37     Bool_t  Refill(MParList &plist, TH1D &h) const;
     44    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set) const;
     45    void    DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
     46    Bool_t  DisplaySize(MParList &plist) const;
     47    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const;
    3848
    3949public:
     
    4555    void EnableRefilling(Bool_t b=kTRUE)  { fRefill=b; }
    4656    void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; }
     57    void EnableRawMc(Bool_t b=kTRUE)      { fRawMc=b; }
     58
     59    void SetEnergyEstimator(const MTask *task);
    4760
    4861    ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum
Note: See TracChangeset for help on using the changeset viewer.