Changeset 7094 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
05/27/05 09:46:22 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
2 added
6 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r7091 r7094  
    2929     - implement usage of extractor resolution (only if UseExtractorRes
    3030       is set).
     31
     32
     33 2005/05/27 Thomas Bretz
     34
     35   * Makefile:
     36     - removed mmontecarlo directory
     37
     38   * mmontecarlo/MMcEnergyEst.[h,cc],
     39     mmontecarlo/MMcTimeGenerate.[h,cc],
     40     mmontecarlo/MMcWeightEnergySpecCalc.[h,cc]:
     41     - removed
     42
     43   * sponde.rc:
     44     - added new line for weighted spectral index
     45
     46   * mbadpixels/MBadPixelsCalc.[h,cc]:
     47     - added an option to perform the checks also in PostProcess
     48
     49   * mhbase/MFillH.h:
     50     - added default argument to SetWeight
     51
     52   * mhbase/MH3.h:
     53     - added Sumw2() member function
     54
     55   * mhflux/MHCollectionArea.[h,cc]:
     56     - added TLatex output to plots
     57     - added some Getter
     58
     59   * mjobs/MJSpectrum.[h,cc]:
     60     - implemented the possibility to weight the monte carlo spectrum
     61       to a new index or function. More details can be found
     62       in MMcSpectrumWeight
     63     - slightly changed the plot comparing the size distributions
     64     - scale the comparsison plots by the resulting spectrum
     65
     66   * mjobs/MJob.[h,cc]:
     67     - added a member function to check ReadEnv of a single
     68       container
     69
     70   * mjobs/Makefile:
     71     - added -I../mmc
     72
     73   * mmc/MMcEvt.[hxx,cxx], mmc/McEvtBasic.[hxx,cxx]:
     74     - changed the inheritance: MMcEvt now derives from MMcEvtBasic
     75       so that both classes are interchangable
     76     - increased both class versions
     77     - chaged the default partictle in MMcEvtBasic from
     78       kGAMMA to kUNDEFINED
     79     - added new particle type: kUNDEFINED
     80
     81   * mhflux/MMcSpectrumWeight.[h,cc]:
     82     - added
    3183
    3284

  • trunk/MagicSoft/Mars/Makefile

    r6979 r7094  
    5656          msql \
    5757          mimage \
    58           mmontecarlo \
    5958          mhft \
    6059          mmc \

  • trunk/MagicSoft/Mars/NEWS

    r7082 r7094  
    22
    33 *** Version <cvs>
     4
     5   - general: MMcEvt now derived from MMcEvtBasic which should
     6     have no influence on compatibility with older camera files
     7
     8   - sponde: The input MC spectrum can now be weighted to fake a
     9     different spectrum. This is done via MMcSpectrumWeight. For
     10     more details see the class description and sponde.rc
     11
     12   - sponde: The paremetr comparsion plots are not scaled by
     13     their entries anymore. Instead the MC plot is scaled by using
     14     the result spectrum of the analysis. If the input MC spectrum
     15     and the result spectrum has different slopes the absolut
     16     normalization is normally wrong.
    417
    518

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc

    r7035 r7094  
    8686//
    8787MBadPixelsCalc::MBadPixelsCalc(const char *name, const char *title)
    88     : fPedestalLevel(3), fPedestalLevelVariance(5), fNamePedPhotCam("MPedPhotCam")
     88    : fPedestalLevel(3), fPedestalLevelVariance(5), fNamePedPhotCam("MPedPhotCam"),
     89    fCheckInProcess(kTRUE), fCheckInPostProcess(kFALSE)
    8990{
    9091    fName  = name  ? name  : gsDefName.Data();
    … …  
    129130// --------------------------------------------------------------------------
    130131//
    131 // Check the pedestal RMS of every pixel with respect to the mean pedestal RMS
    132 // of the camera (Sigmabar).
    133 //
    134 // The pixels can be set as blind if the pedestalRMS is too big or 0.
    135 //
    136 // If you don't want to use this option set the PedestalLevel<=0;
    137 //
    138 //     MBadPixelsCalc calc;
    139 //     calc.SetPedestalLevel(-1);
    140 /*
    141 void MBadPixelsCalc::CheckPedestalRMS() const
    142 {
    143     const Int_t entries = fPedPhotCam->GetSize();
    144 
    145     const Float_t meanPedRMS = fSigmabar->GetSigmabar();
    146 
    147     for (Int_t i=0; i<entries; i++)
    148     {
    149         //
    150         // get pixel as entry from list
    151         //
    152         const Double_t nratio    = fGeomCam->GetPixRatio(i);
    153         const Double_t pixPedRms = (*fPedPhotCam)[i].GetRms();
    154 
    155         if (pixPedRms*nratio > fPedestalLevel * meanPedRMS || pixPedRms == 0)
    156             (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
    157     }
    158 }
    159 */
    160 
    161 // --------------------------------------------------------------------------
    162 //
    163132// Check the pedestal Rms of the pixels: compute with 2 iterations the mean
    164133// for inner and outer pixels. Set as blind the pixels with too small or
    165134// too high pedestal Rms with respect to the mean.
    166135//
    167 Bool_t MBadPixelsCalc::CheckPedestalRms() const
    168 {
     136Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsPix::UnsuitableType_t type) const
     137{
     138    if (fPedestalLevel<=0 && fPedestalLevelVariance<=0)
     139        return kTRUE;
     140
    169141    const Int_t entries = fPedPhotCam->GetSize();
    170142
    … …  
    266238            continue;
    267239
    268         (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
     240        (*fBadPixels)[i].SetUnsuitable(type);
     241        if (type==MBadPixelsPix::kUnsuitableRun)
     242            (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeadPedestalRms);
     243
    269244        bads++;
    270245    }
    … …  
    281256}
    282257
    283 
    284258// --------------------------------------------------------------------------
    285259//
    … …  
    287261Int_t MBadPixelsCalc::Process()
    288262{
    289     if (fPedestalLevel>0 || fPedestalLevelVariance>0)
    290         CheckPedestalRms();
    291 
    292     return kTRUE;
     263    return fCheckInProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableEvt) : kTRUE;
     264}
     265
     266// --------------------------------------------------------------------------
     267//
     268//
     269Int_t MBadPixelsCalc::PostProcess()
     270{
     271    return fCheckInPostProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableRun) : kTRUE;
    293272}
    294273

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h

    r5777 r7094  
    44#ifndef MARS_MTask
    55#include "MTask.h"
     6#endif
     7#ifndef MARS_MBadPixelsPix
     8#include "MBadPixelsPix.h"
    69#endif
    710
    … …  
    2124
    2225    TString fNamePedPhotCam; // name of the 'MPedPhotCam' container
    23    
    24     //    void CheckPedestalRMS() const;
    25     Bool_t CheckPedestalRms() const;
    2626
     27    Bool_t fCheckInProcess;
     28    Bool_t fCheckInPostProcess;
     29
     30    // MBadPixelsCalc
     31    Bool_t CheckPedestalRms(MBadPixelsPix::UnsuitableType_t t) const;
     32
     33    // MTask
    2734    Int_t PreProcess(MParList *pList);
    2835    Int_t Process();
     36    Int_t PostProcess();
     37
     38    // MParContainer
    2939    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    3040
    … …  
    3242    MBadPixelsCalc(const char *name=NULL, const char *title=NULL);
    3343
     44    // Setter
    3445    void SetPedestalLevel(Float_t f)         { fPedestalLevel=f; }
    3546    void SetPedestalLevelVariance(Float_t f) { fPedestalLevelVariance=f; }
    3647    void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; }
     48
     49    void EnableCheckInProcess(Bool_t b=kTRUE)     { fCheckInProcess = b; }
     50    void EnableCheckInPostProcess(Bool_t b=kTRUE) { fCheckInPostProcess = b; }
    3751
    3852    ClassDef(MBadPixelsCalc, 1) // Task to find bad pixels (star, broken pixels, etc)

  • trunk/MagicSoft/Mars/mhbase/MFillH.h

    r6915 r7094  
    6363
    6464    void SetWeight(MParameterD *w)   { fWeight = w; }
    65     void SetWeight(const char *name) { fWeightName = name; }
     65    void SetWeight(const char *name="MWeight") { fWeightName = name; }
    6666
    6767    void      SetDrawOption(Option_t *option="");

  • trunk/MagicSoft/Mars/mhbase/MH3.h

    r6978 r7094  
    2626    Double_t    fScale[3];       // Scale for the three axis (eg unit)
    2727
    28     // TString     fNameProfX;      //! This should make sure, that gROOT doen't confuse the profile with something else
    29     // TString     fNameProfY;      //! This should make sure, that gROOT doen't confuse the profile with something else
    30 
    3128    void StreamPrimitive(ofstream &out) const;
    3229
    … …  
    3633        kIsLogz = BIT(19)
    3734    };
    38 
    39 //    Int_t DistancetoPrimitive(Int_t px, Int_t py);
    40 //    void  ExecuteEvent(Int_t event, Int_t px, Int_t py);
    4135
    4236public:
    … …  
    4842    ~MH3();
    4943
     44    // Setter
    5045    void SetScaleX(Double_t scale) { fScale[0] = scale; }
    5146    void SetScaleY(Double_t scale) { fScale[1] = scale; }
    … …  
    5651    void SetLogz(Bool_t b=kTRUE) { b ? fHist->SetBit(kIsLogz) : fHist->ResetBit(kIsLogz); }
    5752
     53    void Sumw2() const { if (fHist) fHist->Sumw2(); }
     54
     55    // Getter
    5856    Int_t GetDimension() const { return fDimension; }
    5957    Int_t GetNbins() const;
    6058    Int_t FindFixBin(Double_t x, Double_t y=0, Double_t z=0) const;
    6159
    62     void SetName(const char *name);
    63     void SetTitle(const char *title);
    64 
    65     Bool_t SetupFill(const MParList *pList);
    66     Bool_t Fill(const MParContainer *par, const Stat_t w=1);
     60    TH1 &GetHist() { return *fHist; }
     61    const TH1 &GetHist() const { return *fHist; }
    6762
    6863    TString GetDataMember() const;
    6964    TString GetRule(const Char_t axis='x') const;
    7065
    71     TH1 &GetHist() { return *fHist; }
    72     const TH1 &GetHist() const { return *fHist; }
     66    // MH
     67    Bool_t SetupFill(const MParList *pList);
     68    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
    7369
    7470    TH1 *GetHistByName(const TString name="") const { return fHist; }
    … …  
    7975    }
    8076
     77    // MParContainer
     78    MParContainer *New() const;
     79
     80    // TObject
     81    void SetName(const char *name);
     82    void SetTitle(const char *title);
     83
    8184    void SetColors() const;
    8285    void Draw(Option_t *opt=NULL);
    8386    void Paint(Option_t *opt="");
    84 
    85     MParContainer *New() const;
    8687
    8788    ClassDef(MH3, 1) // Generalized 1/2/3D-histogram for Mars variables

  • trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc

    r6976 r7094  
    3131#include "MHCollectionArea.h"
    3232
     33#include <TLatex.h>
    3334#include <TCanvas.h>
    3435#include <TPaveStats.h>
    … …  
    8182    fHistAll.SetTitle("Number of events produced");
    8283    fHistAll.SetXTitle("\\Theta [deg]");
    83     fHistAll.SetYTitle("E [GeV]");
     84    fHistAll.SetYTitle("E_{mc} [GeV]");
    8485    fHistAll.SetDirectory(NULL);
    8586    fHistAll.UseCurrentStyle();
    8687
    8788    fHEnergy.SetName("CollEnergy");
    88     fHEnergy.SetTitle("Collection Area vs Energy");
     89    fHEnergy.SetTitle("Collection Area");
    8990    fHEnergy.SetXTitle("E [GeV]");
    90     fHEnergy.SetYTitle("A [m^{2}]");
     91    fHEnergy.SetYTitle("A_{eff} [m^{2}]");
    9192    fHEnergy.SetDirectory(NULL);
    9293    fHEnergy.UseCurrentStyle();
    … …  
    119120    // read from run header, but it is not yet in!!
    120121    //
    121     const Float_t r1 = 0;
    122     const Float_t r2 = fMcAreaRadius;
    123 
    124     //*fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <<endl<<endl;
    125     const Float_t total_area = TMath::Pi() * (r2*r2 - r1*r1);
     122    const Float_t totalarea = GetCollectionAreaAbs();//TMath::Pi() * (r2*r2 - r1*r1);
    126123
    127124    for (Int_t ix=1; ix<=nbinx; ix++)
    … …  
    141138            continue;
    142139
    143         const Double_t eff = Ns/Na;
    144 
    145         const Double_t efferr = TMath::Sqrt((1.-eff)*Ns)/Na;
     140        const Double_t eff         = Ns/Na;
     141        const Double_t efferr      = TMath::Sqrt((1.-eff)*Ns)/Na;
    146142       
    147         const Float_t col_area  =  eff * total_area;
    148         const Float_t col_area_error =  efferr * total_area;
    149 
    150         fHEnergy.SetBinContent(ix, col_area);
    151         fHEnergy.SetBinError(ix, col_area_error);
     143        const Float_t colarea      =  eff    * totalarea;
     144        const Float_t colareaerror =  efferr * totalarea;
     145
     146        fHEnergy.SetBinContent(ix, colarea);
     147        fHEnergy.SetBinError(ix,   colareaerror);
    152148    }
    153149
    … …  
    275271void MHCollectionArea::Paint(Option_t *option)
    276272{
     273    if (TString(option)=="paint3")
     274    {
     275        /*
     276        TH1 *h = dynamic_cast<TH1*>(gPad->FindObject("Efficiency"));
     277        if (h)
     278        {
     279            const TString txt = Form("N/N_{0}=%.2f",
     280                                 GetCollectionAreaEff(),
     281                                 GetCollectionAreaAbs(), fMcAreaRadius);
     282
     283        TLatex text(0.31, 0.95, txt);
     284        text.SetBit(TLatex::kTextNDC);
     285        text.SetTextSize(0.04);
     286        text.Paint();*/
     287        return;
     288    }
     289    if (TString(option)=="paint4")
     290    {
     291        const TString txt = Form("A_{eff}=%.0fm^{2}  A_{abs}=%.0fm^{2}  r=%.0fm",
     292                                 GetCollectionAreaEff(),
     293                                 GetCollectionAreaAbs(), fMcAreaRadius);
     294
     295        TLatex text(0.31, 0.95, txt);
     296        text.SetBit(TLatex::kTextNDC);
     297        text.SetTextSize(0.04);
     298        text.Paint();
     299        return;
     300    }
     301
    277302    TVirtualPad *pad;
    278303
    … …  
    405430        h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency");
    406431        h->SetDirectory(NULL);
     432        AppendPad("paint3");
    407433    }
    408434    else
    … …  
    416442        gPad->SetGridy();
    417443        fHEnergy.Draw();
     444        AppendPad("paint4");
    418445    }
    419446    else
    … …  
    428455    //*fLog << energy << " " << theta << endl;
    429456
    430     fHistSel.Fill(theta, energy);
     457    fHistSel.Fill(theta, energy, weight);
    431458
    432459    return kTRUE;

  • trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h

    r6976 r7094  
    5656    const TH1D &GetHEnergy() const { return fHEnergy; }
    5757
     58    Double_t GetEntries() const { return fHistAll.Integral(); }
     59    Double_t GetCollectionAreaEff() const { return fHEnergy.Integral(); }
     60    Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; }
     61
    5862/*
    5963    void DrawAll(Option_t *option="");

  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

    r7071 r7094  
    5252#include "MBinning.h"
    5353#include "MDataSet.h"
     54#include "MMcCorsikaRunHeader.h"
    5455
    5556// Spectrum
    … …  
    5859#include "../mhflux/MHCollectionArea.h"
    5960#include "../mhflux/MHEnergyEst.h"
     61#include "../mhflux/MMcSpectrumWeight.h"
    6062
    6163// Eventloop
    … …  
    155157}
    156158
     159// --------------------------------------------------------------------------
     160//
     161// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
     162// the dataset.
     163// The simulated energy range and spectral slope is initialized from
     164// there.
     165// In the following eventloops the forced check in MMcSpectrumWeight
     166// ensures, that the spectral slope and energy range doesn't change.
     167//
     168Bool_t  MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
     169{
     170    fLog->Separator("Initialize energy weighting");
     171
     172    MParList l;
     173    l.AddToList(&w);
     174    if (!CheckEnv(l))
     175    {
     176        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
     177        return kFALSE;
     178    }
     179
     180    TChain chain("RunHeaders");
     181    set.AddFilesOn(chain);
     182
     183    MMcCorsikaRunHeader *h=0;
     184    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
     185    chain.GetEntry(1);
     186
     187    if (!h)
     188    {
     189        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
     190        return kFALSE;
     191    }
     192
     193    if (!w.Set(*h))
     194    {
     195        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
     196        return kFALSE;
     197    }
     198
     199    w.Print();
     200    return kTRUE;
     201}
     202
    157203Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
    158204{
    … …  
    211257}
    212258
    213 Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
     259Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
    214260{
    215261    // Some debug output
    216262    fLog->Separator("Compiling original MC distribution");
    217263
    218     *fLog << inf << "Please stand by, this may take a while..." << flush;
     264    weight.SetNameMcEvt("MMcEvtBasic");
     265    const TString w(weight.GetFormulaWeights());
     266    weight.SetNameMcEvt();
     267
     268    *fLog << inf << "Using weights: " << w << endl;
     269    *fLog << "Please stand by, this may take a while..." << flush;
    219270
    220271    if (fDisplay)
    … …  
    236287        h.SetYTitle("E [GeV]");
    237288        h.SetZTitle("Counts");
    238         chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
     289        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
    239290    }
    240291    else
    … …  
    243294        h.SetXTitle("\\Theta [\\circ]");
    244295        h.SetYTitle("Counts");
    245         chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
     296        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
    246297    }
    247298    h.SetDirectory(0);
    … …  
    413464}
    414465
    415 Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
     466Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
    416467{
    417468    MTaskList tlist1;
    … …  
    439490    }
    440491    tlist1.AddToList(&readmc);
     492    tlist1.AddToList(&weight);
    441493
    442494    temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
    … …  
    477529        bins3->SetName("BinningTheta");
    478530    }
     531
    479532    return kTRUE;
    480533}
    481534
    482 void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
     535TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
    483536{
    484537    TH1D collarea(area.GetHEnergy());
    … …  
    546599    f.SetParameter(1, 1.9e-6);
    547600    f.SetLineColor(kGreen);
    548     spectrum.Fit(&f, "NI", "", 55, 2e4);
     601    spectrum.Fit(&f, "NIM", "", 55, 2e4);
    549602    f.DrawCopy("same");
    550603
    … …  
    558611     tex.DrawLatex(2e2, 7e-5, str);
    559612     */
    560 }
    561 
    562 Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
    563 {
    564     cout << name << endl;
     613
     614    TArrayD res(2);
     615    res[0] = f.GetParameter(0);
     616    res[1] = f.GetParameter(1);
     617    return res;
     618}
     619
     620Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
     621{
    565622    TString same(name);
    566623    same += "Same";
    … …  
    587644
    588645    const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
    589     const Double_t scale = fit ? fit->GetScaleFactor() : 1;
     646    const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
    590647
    591648    gPad->SetBorderMode(0);
    592649    h2->SetLineColor(kBlack);
    593650    h3->SetLineColor(kBlue);
    594     h2->Add(h1, -scale);
    595 
    596     h2->Scale(1./h2->Integral());
    597     h3->Scale(1./h3->Integral());
     651    h2->Add(h1, -ascale);
     652
     653    //h2->Scale(1./ontime);   //h2->Integral());
     654    h3->Scale(scale);    //h3->Integral());
    598655
    599656    h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
    … …  
    609666}
    610667
    611 Bool_t MJSpectrum::DisplaySize(MParList &plist) const
     668Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
    612669{
    613670    *fLog << inf << "Reading from file: " << fPathIn << endl;
    … …  
    644701
    645702    excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
    646     excess->Scale(1./excess->Integral());
    647     excess = excess->DrawCopy();
     703    excess = excess->DrawCopy("E2");
    648704    // Don't do this on the original object!
    649705    excess->SetStats(kFALSE);
     706    excess->SetMarkerStyle(kFullDotMedium);
     707    excess->SetFillColor(kBlack);
     708    excess->SetFillStyle(0);
     709    excess->SetName("Excess  ");
     710    excess->SetDirectory(0);
    650711
    651712    TObject *o=0;
    652     if ((o=plist.FindObject("ExcessSize")))
     713    if ((o=plist.FindObject("ExcessMC")))
    653714    {
    654715        TH1 *histsel = (TH1F*)o->FindObject("");
    655716        if (histsel)
    656717        {
    657             histsel->Scale(1./histsel->Integral());
     718            if (scale<0)
     719                scale = excess->Integral()/histsel->Integral();
     720
     721            histsel->Scale(scale);
    658722            histsel->SetLineColor(kBlue);
    659723            histsel->SetBit(kCanDelete);
    660             histsel = histsel->DrawCopy("same");
     724            histsel = histsel->DrawCopy("E1 same");
    661725            // Don't do this on the original object!
    662726            histsel->SetStats(kFALSE);
     727
     728            fLog->Separator("Kolmogorov Test");
     729            histsel->KolmogorovTest(excess, "DX");
     730            fLog->Separator("Chi^2 Test");
     731            histsel->Chi2Test(excess, "P");
    663732        }
    664733    }
    … …  
    666735    // -------------- Comparison of Image Parameters --------------
    667736    c.cd(2);
    668     PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost");
     737    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost", scale);
    669738
    670739    c.cd(3);
    671     PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
     740    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
    672741
    673742    c.cd(4);
    674     PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost");
     743    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost", scale);
    675744
    676745    c.cd(5);
    677     PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost");
     746    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost", scale);
    678747
    679748    c.cd(6);
    680     PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost");
     749    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);
    681750
    682751    return kTRUE;
    … …  
    734803    }
    735804
     805    MMcSpectrumWeight weight;
     806    if (!InitWeighting(set, weight))
     807        return kFALSE;
     808
    736809    PrintSetup(fit);
    737810    bins3.SetEdges(temp1, 'x');
    738811
    739812    TH1D temp2(temp1);
    740     if (!ReadOrigMCDistribution(set, temp2))
     813    if (!ReadOrigMCDistribution(set, temp2, weight))
    741814        return kFALSE;
    742815
    … …  
    754827        hist.UseCurrentStyle();
    755828        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
    756         if (!ReadOrigMCDistribution(set, hist))
     829        if (!ReadOrigMCDistribution(set, hist, weight))
    757830            return kFALSE;
    758831
    … …  
    762835                for (int x=0; x<hist.GetNbinsX(); x++)
    763836                    hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
     837            //hist.SetEntries(hist.Integral());
    764838        }
    765839    }
    766840    else
    767         if (!IntermediateLoop(plist, mh1, temp1, set))
     841    {
     842        weight.SetNameMcEvt("MMcEvtBasic");
     843        if (!IntermediateLoop(plist, mh1, temp1, set, weight))
    768844            return kFALSE;
     845        weight.SetNameMcEvt();
     846    }
    769847
    770848    DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
    … …  
    813891    MFillH fill3(&area, "", "FillCollectionArea");
    814892    MFillH fill4(&hest, "", "FillEnergyEst");
     893    fill3.SetWeight();
     894    fill4.SetWeight();
    815895
    816896    MH3 hsize("MHillas.fSize");
    817     //MH3 henergy("MEnergyEst.fVal");
    818     hsize.SetName("ExcessSize");
    819     //henergy.SetName("EnergyEst");
    820     MBinning bins(size, "BinningExcessSize");
     897    hsize.SetName("ExcessMC");
     898    hsize.Sumw2();
     899
     900    MBinning bins(size, "BinningExcessMC");
    821901    plist.AddToList(&hsize);
    822     //plist.AddToList(&henergy);
    823902    plist.AddToList(&bins);
    824903
    … …  
    830909    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
    831910    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
    832     MFillH fill8a("ExcessSize     [MH3]",           "",             "FillExcessSize");
    833     //MFillH fill9a("EnergyEst      [MH3]",           "",             "FillExcessEEst");
     911    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
    834912    fill1a.SetNameTab("PreCut");
    835913    fill2a.SetNameTab("PostCut");
    … …  
    840918    fill7a.SetNameTab("NewPar");
    841919    fill8a.SetBit(MFillH::kDoNotDisplay);
    842     //fill9a.SetBit(MFillH::kDoNotDisplay);
     920    fill1a.SetWeight();
     921    fill2a.SetWeight();
     922    fill3a.SetWeight();
     923    fill4a.SetWeight();
     924    fill5a.SetWeight();
     925    fill6a.SetWeight();
     926    fill7a.SetWeight();
     927    fill8a.SetWeight();
    843928
    844929    MEnergyEstimate est;
    … …  
    852937    tlist2.AddToList(&hcalc1);
    853938    tlist2.AddToList(&hcalc2);
     939    tlist2.AddToList(&weight);
    854940    tlist2.AddToList(&fill1a);
    855941    tlist2.AddToList(fCut0);
    … …  
    891977    // -------------------------- Spectrum ----------------------------
    892978
    893     DisplaySpectrum(area, excess, hest, ontime);
    894     DisplaySize(plist);
    895 
     979    // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
     980    TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
     981
     982    // Spectrum fitted (convert res[1] from TeV to GeV)
     983    TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
     984
     985    // Number of events this spectrum would produce per s and m^2
     986    Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
     987
     988    // scale with effective collection area to get the event rate (N/s)
     989    // scale with the effective observation time to absolute observed events
     990    n *= area.GetCollectionAreaAbs()*ontime; // N
     991
     992    // Now calculate the scale factor from the number of events
     993    // produced and the number of events which should have been
     994    // observed with our telescope in the time ontime
     995    const Double_t scale = n/area.GetEntries();
     996
     997    // Print normalization constant
     998    cout << "MC normalization factor:  " << scale << endl;
     999
     1000    // Overlay normalized plots
     1001    DisplaySize(plist, scale);
     1002
     1003    // check if output should be written
    8961004    if (!fPathOut.IsNull())
    8971005        fDisplay->SaveAsRoot(fPathOut);

  • trunk/MagicSoft/Mars/mjobs/MJSpectrum.h

    r6978 r7094  
    1818class MStatusArray;
    1919class MHCollectionArea;
     20class MMcSpectrumWeight;
    2021
    2122class MJSpectrum : public MJob
    … …  
    3536    Bool_t  ReadTask(MTask* &task, const char *name) const;
    3637    Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size);
    37     Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const;
     38    Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &w) const;
    3839    Bool_t  GetThetaDistribution(TH1D &temp1, TH1D &temp2) const;
    3940    Bool_t  Refill(MParList &plist, TH1D &h) const;
     41    Bool_t  InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
    4042
    4143    // Display Output
    4244    void    PrintSetup(const MAlphaFitter &fit) const;
    4345    void    DisplayResult(const TH2D &mh1) 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;
     46    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const;
     47    TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
     48    Bool_t  DisplaySize(MParList &plist, Double_t scale) const;
     49    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const;
    4850
    4951public:

  • trunk/MagicSoft/Mars/sponde.rc

    r6977 r7094  
    11EstimateEnergy.Rule: (0.380075+(MPointingPos.fZd*MPointingPos.fZd*0.00109028))*pow(MHillas.fSize,0.892462)
     2
     3#MMcSpectrumWeight.NewSlope: -2.26
Note: See TracChangeset for help on using the changeset viewer.