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

Legend:

Unmodified
Added
Removed

  • 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:
Note: See TracChangeset for help on using the changeset viewer.