Ignore:
Timestamp:
11/15/05 17:01:07 (19 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

    r7361 r7403  
    7676#include "MTaskEnv.h"
    7777#include "MFillH.h"
    78 #include "MHillasCalc.h"
     78//#include "MHillasCalc.h"
    7979//#include "MSrcPosCalc.h"
    8080#include "MContinue.h"
    … …  
    8787    : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
    8888    fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE),
    89     fNoThetaWeights(kFALSE)
     89    fNoThetaWeights(kFALSE), fWeighting(kTRUE)
    9090{
    9191    fName  = name  ? name  : "MJSpectrum";
    … …  
    551551}
    552552
     553TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
     554{
     555    const Double_t p0 = f.GetParameter(0);
     556    const Double_t p1 = f.GetParameter(1);
     557
     558    const Double_t e0 = f.GetParError(0);
     559    const Double_t e1 = f.GetParError(1);
     560
     561    const Int_t    np  = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
     562    const Double_t exp = TMath::Power(10, np);
     563
     564    const Float_t fe0 = TMath::Log10(e0);
     565    const Float_t fe1 = TMath::Log10(e1);
     566
     567    // calc number of (gueltige ziffern) digits to be displayed
     568    const char *f0  = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0";
     569    const char *f1  = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0";
     570
     571    TString str;
     572    switch (type)
     573    {
     574    case 0:
     575        {
     576            const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",   f1, f1);
     577            const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
     578
     579            str  = Form(fmt0.Data(), p1/exp, e1/exp, np);
     580            str += Form(fmt1.Data(), p0, e0);
     581            str += "\\frac{ph}{TeVm^{2}s}";
     582        }
     583        break;
     584    case 1:
     585        str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
     586        break;
     587    case 2:
     588        str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(),  f.GetNDF()));
     589        break;
     590    }
     591    return str;
     592}
     593
     594TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
     595{
     596    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
     597    f.SetParameter(0, -2.5);
     598    f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000)));
     599    f.SetLineColor(kBlue);
     600    spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped
     601    f.DrawCopy("same");
     602
     603    /*
     604    const Double_t p0 = f.GetParameter(0);
     605    const Double_t p1 = f.GetParameter(1);
     606
     607    const Double_t e0 = f.GetParError(0);
     608    const Double_t e1 = f.GetParError(1);
     609
     610    const Int_t    np  = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
     611    const Double_t exp = TMath::Power(10, np);
     612    */
     613    TString str = FormString(f);
     614    /*
     615    str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
     616    str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
     617    str += "\\frac{ph}{TeVm^{2}s}";
     618    */
     619    TLatex tex;
     620    tex.SetTextSize(0.045);
     621    tex.SetBit(TLatex::kTextNDC);
     622    tex.SetTextAlign(31);
     623    tex.DrawLatex(0.89, 0.935, str);
     624
     625    str = FormString(f, 1);
     626    tex.DrawLatex(0.89, 0.83, str);
     627
     628    str = FormString(f, 2);
     629    tex.DrawLatex(0.89, 0.735, str);
     630
     631    TArrayD res(2);
     632    res[0] = f.GetParameter(0);
     633    res[1] = f.GetParameter(1);
     634    return res;
     635}
     636
    553637// --------------------------------------------------------------------------
    554638//
    … …  
    568652
    569653    cout << "Effective On time: " << ontime << "s" << endl;
     654
     655//    for (int i=0; i<collarea.GetNbinsX(); i++)
     656//        collarea.SetBinError(i+1, 0);//collarea.GetBinError(i+1)/2);
    570657
    571658    spectrum.SetDirectory(NULL);
    … …  
    602689    weights.DrawCopy();
    603690
    604     //spectrum.Multiply(&weights);
     691    //spectrum.Multiply(&weights); // Currentyl we skip the error distribution!
    605692    spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
    606693    spectrum.SetBit(TH1::kNoStats);
    … …  
    608695    for (int i=0; i<excess.GetNbinsX(); i++)
    609696    {
    610         spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
    611         spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
     697        if (i+1>2)
     698        {
     699            //Float_t E = spectrum.GetXaxis()->GetBinWidth(i+1);
     700            //Float_t c = log10(E)*3.012e-1+1.905e-1;
     701            //spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*c);
     702            //spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*c);
     703            spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
     704            spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
     705        }
    612706
    613707        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
    614         spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
     708        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  /spectrum.GetBinWidth(i+1)*1000);
    615709    }
    616710
    … …  
    623717    spectrum.SetMinimum(1e-12);
    624718    spectrum.SetXTitle("E [GeV]");
    625     spectrum.SetYTitle("N/TeVsm^{2}");
     719    spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
     720    TH1D *spec = (TH1D*)spectrum.DrawCopy();
     721
     722    // Calculate Spectrum * E^2
     723    for (int i=0; i<spectrum.GetNbinsX(); i++)
     724    {
     725        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
     726        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
     727    }
     728    fDisplay->AddTab("Check");
     729    gPad->SetLogx();
     730    gPad->SetLogy();
     731    gPad->SetGrid();
     732
     733    spectrum.SetName("FluxStd");
     734    spectrum.SetMarkerStyle(kFullDotMedium);
     735    spectrum.SetTitle("Differential flux times E^{2}");
     736    spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
     737    spectrum.SetDirectory(0);
    626738    spectrum.DrawCopy();
    627739
     740    // Fit Spectrum
     741    c1.cd(4);
     742    return FitSpectrum(*spec);
     743    /*
    628744    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
    629745    f.SetParameter(0, -2.87);
    630746    f.SetParameter(1, 1.9e-6);
    631     f.SetLineColor(kGreen);
     747    f.SetLineColor(kBlue);
    632748    spectrum.Fit(&f, "NIM", "", 100, 5000);
    633749    f.DrawCopy("same");
    … …  
    657773    TArrayD res(2);
    658774    res[0] = f.GetParameter(0);
    659     res[1] = f.GetParameter(1);
     775    res[1] = f.GetParameter(1);*/
    660776
    661777/*
    … …  
    691807    f.SetLineStyle(kDashed);
    692808    f.DrawCopy("same");
     809    return res;
    693810 */
    694     return res;
    695811}
    696812
    … …  
    754870// Calls PlotSame
    755871//
    756 Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
     872//#include <TSpline.h>
     873Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale, TH1D &result) const
    757874{
    758875    *fLog << inf << "Reading from file: " << fPathIn << endl;
    … …  
    773890    }
    774891
    775     TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
     892    TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
    776893    if (!excess)
    777894        return kFALSE;
    … …  
    789906
    790907    excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
    791     excess = excess->DrawCopy("E2");
     908    excess = (TH1D*)excess->DrawCopy("E2");
    792909    // Don't do this on the original object!
    793910    excess->SetStats(kFALSE);
    … …  
    801918    if ((o=plist.FindObject("ExcessMC")))
    802919    {
    803         TH1 *histsel = (TH1F*)o->FindObject("");
     920        TH1 *histsel = (TH1D*)o->FindObject("");
    804921        if (histsel)
    805922        {
     923            // Calculate the weights such that the number of total
     924            // MC events after cuts will stay constant
     925
     926            /*
     927            Int_t first = -1;
     928            for (int i=0; i<excess->GetNbinsX(); i++)
     929                if (excess->GetBinContent(i+1)>0)
     930                {
     931                    first = i+1;
     932                    break;
     933                }
     934
     935            Int_t last = -1;
     936            for (int i=excess->GetNbinsX(); i>0; i--)
     937                if (excess->GetBinContent(i)>0)
     938                {
     939                    last = i;
     940                    break;
     941                    }*/
     942            /*
     943            Int_t first = 1;
     944            Int_t last  = excess->GetNbinsX();
     945
     946            const Double_t xmin = TMath::Log10(excess->GetBinLowEdge(first));
     947            const Double_t xmax = TMath::Log10(excess->GetBinLowEdge(last+1));
     948            const Int_t    xnum = last-first+1;
     949
     950            //            Double_t xmin = TMath::Log10(excess->GetXaxis()->GetXmin());
     951            //            Double_t xmax = TMath::Log10(excess->GetXaxis()->GetXmax());
     952            //            Int_t    xnum = excess->GetNbinsX();
     953
     954            TArrayD arr1(xnum);
     955            TArrayD arr2(xnum);
     956            for (int i=0; i<xnum; i++)
     957            {
     958                Double_t v1 = excess->GetBinContent(i+first);
     959                Double_t v2 = histsel->GetBinContent(i+first);
     960                arr1[i] = v1<1 ? -1 : TMath::Log10(v1);
     961                arr2[i] = v2<1 ? -1 : TMath::Log10(v2);
     962            }
     963            cout << endl << endl;
     964
     965            TSpline5 sp1("Excess", xmin, xmax, arr1.GetArray(), xnum+1);
     966            TSpline5 sp2("MC",     xmin, xmax, arr2.GetArray(), xnum+1);
     967
     968            MBinning bins(xnum*5, TMath::Power(10, xmin), TMath::Power(10, xmax), "", "log");
     969            bins.Apply(result);
     970
     971            for (int i=0; i<result.GetNbinsX(); i++)
     972            {
     973                Double_t lo = result.GetXaxis()->GetBinLowEdge(i+1);
     974                Double_t up = result.GetXaxis()->GetBinUpEdge(i+1);
     975
     976                Double_t x  = TMath::Log10(lo*up)/2;
     977                Double_t v1 = sp1.Eval(x);
     978                Double_t v2 = sp2.Eval(x);
     979                Double_t rc = v1<-0.9 || v2<-0.9 ? 0 : TMath::Power(10, v1-v2);
     980                result.SetBinContent(i+1, rc);
     981            }
     982            */
     983            excess->Copy(result);
     984            result.Divide(histsel);
     985            result.SetDirectory(0);
     986
     987            // Apply additional absolute scale factor
    806988            if (scale<0)
    807989                scale = excess->Integral()/histsel->Integral();
    … …  
    814996            histsel->SetStats(kFALSE);
    815997
     998            // Do some Chi^2 and Kolmogorov tests
    816999            fLog->Separator("Kolmogorov Test");
    8171000            histsel->KolmogorovTest(excess, "DX");
    … …  
    8191002            const Double_t p = histsel->Chi2Test(excess, "P");
    8201003
     1004            // Plot result
    8211005            TLatex tex;
    8221006            tex.SetBit(TLatex::kTextNDC);
    … …  
    8401024    c.cd(6);
    8411025    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);
     1026
     1027    return kTRUE;
     1028}
     1029
     1030Bool_t MJSpectrum::DisplaySpectrumW(TH1D &energy0, TH1 &energy, TH1 &hsize, TH1D &weights, TH2D &hall, Double_t scale, TArrayD &res) const
     1031{
     1032    *fLog << inf << "Reading from file: " << fPathIn << endl;
     1033
     1034    TFile file(fPathIn, "READ");
     1035    if (!file.IsOpen())
     1036    {
     1037        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
     1038        return kFALSE;
     1039    }
     1040
     1041    file.cd();
     1042    MStatusArray arr;
     1043    if (arr.Read()<=0)
     1044    {
     1045        *fLog << "MStatusDisplay not found in file... abort." << endl;
     1046        return kFALSE;
     1047    }
     1048
     1049    TH1D *exc = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
     1050    if (!exc)
     1051        return kFALSE;
     1052
     1053    // energy: energy distribution after weighting
     1054    // energy0: energy distribution before weighting
     1055    // weights: weights to weight events depending on their size
     1056
     1057    TCanvas &c = fDisplay->AddTab("SpectrumW");
     1058    c.Divide(3,2);
     1059    c.cd(1);
     1060    gPad->SetBorderMode(0);
     1061    gPad->SetLogx();
     1062    gPad->SetLogy();
     1063    gPad->SetGrid();
     1064
     1065    TH1D *e = (TH1D*)exc->DrawCopy("E2");
     1066    // Don't do this on the original object!
     1067    e->SetName("ExcessData");
     1068    e->SetTitle("Excess events vs Size (data/black, mc/blue)");
     1069    e->SetDirectory(0);
     1070    e->SetStats(kFALSE);
     1071    e->SetBit(TH1::kNoStats);
     1072    e->ResetBit(TH1::kNoTitle);
     1073    e->SetMarkerStyle(kFullDotMedium);
     1074    e->SetFillColor(kBlack);
     1075    e->SetFillStyle(0);
     1076    e->SetName("Excess  ");
     1077    e->SetDirectory(0);
     1078
     1079    e = (TH1D*)hsize.DrawCopy("E1 same");
     1080    // Don't do this on the original object!
     1081    e->SetName("ExcessMC");
     1082    e->SetDirectory(0);
     1083    e->SetBit(kCanDelete);
     1084    e->SetBit(TH1::kNoStats);
     1085    e->ResetBit(TH1::kNoTitle);
     1086    e->SetLineColor(kBlue);
     1087    e->SetStats(kFALSE);
     1088
     1089    c.cd(2);
     1090    gPad->SetBorderMode(0);
     1091    gPad->SetLogx();
     1092    gPad->SetGrid();
     1093
     1094    weights.SetName("WeightsSize");
     1095    weights.SetDirectory(0);
     1096    weights.SetBit(TH1::kNoStats);
     1097    weights.ResetBit(TH1::kNoTitle);
     1098    weights.SetTitle("Size weights: size(data)/size(mc)");
     1099    weights.DrawCopy();
     1100
     1101    c.cd(4);
     1102    gPad->SetBorderMode(0);
     1103    gPad->SetLogx();
     1104    gPad->SetLogy();
     1105    gPad->SetGrid();
     1106    energy.SetBit(TH1::kNoStats);
     1107    energy.ResetBit(TH1::kNoTitle);
     1108    energy.SetTitle("Normalized energy distribution before (black) and after (blue) weighting");
     1109    energy.SetLineColor(kBlue);
     1110    TH1 *he = energy.DrawCopy();
     1111    TH1 *hf = energy0.DrawCopy("same");
     1112    he->Scale(1./he->Integral());
     1113    hf->Scale(1./hf->Integral());
     1114    he->SetName("EnergyOrig");
     1115    hf->SetName("EnergyWeighted");
     1116    he->SetDirectory(0);
     1117    hf->SetDirectory(0);
     1118
     1119    c.cd(5);
     1120    gPad->SetBorderMode(0);
     1121    gPad->SetLogx();
     1122    gPad->SetGrid();
     1123    energy0.SetBit(TH1::kNoStats);
     1124    energy0.ResetBit(TH1::kNoTitle);
     1125    energy0.SetName("WeightsEnergy");
     1126    energy0.SetTitle("Energy weights: E(before weighting)/E(after weighting)");
     1127    energy0.SetDirectory(0);
     1128    energy0.Divide(&energy);
     1129    energy0.DrawCopy();
     1130
     1131    c.cd(3);
     1132    gPad->SetBorderMode(0);
     1133    gPad->SetLogx();
     1134    gPad->SetLogy();
     1135    gPad->SetGrid();
     1136
     1137    TH1D *h1 = (TH1D*)hall.ProjectionY("ProjAllY", -1, 9999, "E");
     1138    h1->SetTitle("Normalized spectrum before and after Energy weights applied");
     1139    h1->SetDirectory(NULL);
     1140    h1->SetXTitle("E [GeV]");
     1141    h1->SetBit(kCanDelete);
     1142    h1->SetBit(TH1::kNoStats);
     1143    he = h1->DrawCopy();
     1144    h1->SetLineColor(kBlue);
     1145    // Convert the mc spectrum before weights to the mc spectrum after weights
     1146    // Therefor weights after/before are applied
     1147    h1->Divide(&energy0);
     1148    hf = h1->DrawCopy("same");
     1149    he->Scale(1./he->GetBinContent(he->GetNbinsX()/2));
     1150    hf->Scale(1./hf->GetBinContent(hf->GetNbinsX()/2));
     1151    he->SetName("SpectrumOrig");
     1152    hf->SetName("SpectrumWeighted");
     1153    he->SetDirectory(0);
     1154    hf->SetDirectory(0);
     1155
     1156    c.cd(6);
     1157    gPad->SetBorderMode(0);
     1158    gPad->SetLogx();
     1159    gPad->SetLogy();
     1160    gPad->SetGrid();
     1161
     1162    h1->SetName("Spectrum");
     1163    h1->SetTitle("Spectrum");
     1164    h1->SetXTitle("E [GeV]");
     1165    h1->SetYTitle("dN/dE [N/TeVsm^{2}]");
     1166    h1->SetDirectory(0);
     1167    h1->Scale(1./scale);
     1168
     1169    for (int i=0; i<h1->GetNbinsX(); i++)
     1170    {
     1171        h1->SetBinContent(i+1, h1->GetBinContent(i+1)/h1->GetBinWidth(i+1)*1000);
     1172        h1->SetBinError(  i+1, h1->GetBinError(i+1)  /h1->GetBinWidth(i+1)*1000);
     1173    }
     1174
     1175    TH1D *h2 = (TH1D*)h1->DrawCopy();
     1176    res = FitSpectrum(*h2);
     1177
     1178    // Calculate and Display Spectrum * E^2
     1179    TCanvas *c1 = fDisplay->GetCanvas("Check");
     1180    if (c1)
     1181    {
     1182        c1->cd();
     1183        gPad->SetLogx();
     1184        gPad->SetLogy();
     1185        gPad->SetGrid();
     1186
     1187        for (int i=0; i<h1->GetNbinsX(); i++)
     1188        {
     1189            h1->SetBinContent(i+1, h1->GetBinContent(i+1)*h1->GetBinWidth(i+1)*h1->GetBinWidth(i+1));
     1190            h1->SetBinError(  i+1, h1->GetBinError(i+1)  *h1->GetBinWidth(i+1)*h1->GetBinWidth(i+1));
     1191        }
     1192        h1->SetName("FluxW");
     1193        h1->SetTitle("Differential flux times E^{2}");
     1194        h1->SetDirectory(0);
     1195        h1->SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
     1196        h1->SetLineColor(kBlue);
     1197        h1->SetMarkerColor(kBlue);
     1198        h1->SetMarkerStyle(kFullDotMedium);
     1199        h1->Draw("same");
     1200    }
    8421201
    8431202    return kTRUE;
    … …  
    9461305    // ------------------------- Final loop --------------------------
    9471306
    948     *fLog << endl;
    949     fLog->Separator("Calculate efficiencies");
    950     *fLog << endl;
    951 
    952     MTaskList tlist2;
    953     plist.Replace(&tlist2);
    954 
    955     MReadMarsFile read("Events");
    956     read.DisableAutoScheme();
    957     set.AddFilesOn(read);
    958 
    959     // Selector to get correct (final) theta-distribution
    960     temp1.SetXTitle("MPointingPos.fZd");
    961 
    962     MH3 mh3(temp1);
    963 
    964     MFEventSelector2 sel2(mh3);
    965     sel2.SetHistIsProbability();
    966 
    967     MContinue contsel(&sel2);
    968     contsel.SetInverted();
    969 
    970     // Get correct source position
    971     //MSrcPosCalc calc;
    972 
    973     // Calculate corresponding Hillas parameters
    974     /*
    975      MHillasCalc hcalc1;
    976      MHillasCalc hcalc2("MHillasCalcAnti");
    977      hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
    978      hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
    979      hcalc2.SetNameHillasSrc("MHillasSrcAnti");
    980      hcalc2.SetNameSrcPosCam("MSrcPosAnti");
    981      */
    982 
    983     // Fill collection area and energy estimator (unfolding)
    984     // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
    985     MHCollectionArea area;
    986     area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
    987     MHEnergyEst      hest;
    988 
    989     MFillH fill3(&area, "", "FillCollectionArea");
    990     MFillH fill4(&hest, "", "FillEnergyEst");
    991     MFillH fill5("MHThreshold", "", "FillThreshold");
    992     fill3.SetWeight();
    993     fill4.SetWeight();
    994     fill5.SetWeight();
    995     fill3.SetNameTab("ColArea");
    996     fill4.SetNameTab("E-Est");
    997     fill5.SetNameTab("Threshold");
    998 
    999     MH3 hsize("MHillas.fSize");
    1000     hsize.SetName("ExcessMC");
    1001     hsize.Sumw2();
    1002 
    1003     MBinning bins(size, "BinningExcessMC");
    1004     plist.AddToList(&hsize);
    1005     plist.AddToList(&bins);
    1006 
    1007     MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
    1008     MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
    1009     MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
    1010     MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
    1011     MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
    1012     MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
    1013     MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
    1014     MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
    1015     fill1a.SetNameTab("PreCut");
    1016     fill2a.SetNameTab("PostCut");
    1017     fill3a.SetNameTab("VsSize");
    1018     fill4a.SetNameTab("HilExt");
    1019     fill5a.SetNameTab("HilSrc");
    1020     fill6a.SetNameTab("ImgPar");
    1021     fill7a.SetNameTab("NewPar");
    1022     fill8a.SetBit(MFillH::kDoNotDisplay);
    1023     fill1a.SetWeight();
    1024     fill2a.SetWeight();
    1025     fill3a.SetWeight();
    1026     fill4a.SetWeight();
    1027     fill5a.SetWeight();
    1028     fill6a.SetWeight();
    1029     fill7a.SetWeight();
    1030     fill8a.SetWeight();
    1031 
    1032     MEnergyEstimate est;
    1033     MTaskEnv taskenv1("EstimateEnergy");
    1034     taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
    1035 
    1036     tlist2.AddToList(&read);
    1037     if (!fRawMc && fNoThetaWeights)
    1038         tlist2.AddToList(&contsel);
    1039     //tlist2.AddToList(&calc);
    1040     //tlist2.AddToList(&hcalc1);
    1041     //tlist2.AddToList(&hcalc2);
    1042     tlist2.AddToList(&weight);
    1043     tlist2.AddToList(&fill1a);
    1044     tlist2.AddToList(fCut0);
    1045     tlist2.AddToList(fCut1);
    1046     tlist2.AddToList(fCut2);
    1047     tlist2.AddToList(fCut3);
    1048     tlist2.AddToList(&taskenv1);
    1049     tlist2.AddToList(&fill3);
    1050     tlist2.AddToList(&fill4);
    1051     tlist2.AddToList(&fill5);
    1052     tlist2.AddToList(&fill2a);
    1053     tlist2.AddToList(&fill3a);
    1054     tlist2.AddToList(&fill4a);
    1055     tlist2.AddToList(&fill5a);
    1056     tlist2.AddToList(&fill6a);
    1057     tlist2.AddToList(&fill7a);
    1058     tlist2.AddToList(&fill8a);
    1059     //tlist2.AddToList(&fill9a);
    1060 
    1061     MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
    1062     loop2.SetParList(&plist);
    1063     loop2.SetDisplay(fDisplay);
    1064     loop2.SetLogStream(fLog);
    1065 
    1066     if (!SetupEnv(loop2))
    1067         return kFALSE;
    1068 
    1069     if (!loop2.Eventloop(fMaxEvents))
    1070     {
    1071         *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
    1072         return kFALSE;
    1073     }
    1074 
    1075     if (!loop2.GetDisplay())
    1076     {
    1077         *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
    1078         return kFALSE;
    1079     }
    1080 
    1081     gLog.Separator("Energy Estimator");
    1082     if (plist.FindObject("EstimateEnergy"))
    1083         plist.FindObject("EstimateEnergy")->Print();
    1084 
    1085     gLog.Separator("Spectrum");
    1086 
    1087     // -------------------------- Spectrum ----------------------------
    1088 
    1089     // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
    1090     TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
    1091 
    1092     // Spectrum fitted (convert res[1] from TeV to GeV)
    1093     TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
    1094 
    1095     // Number of events this spectrum would produce per s and m^2
    1096     Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
    1097 
    1098     // scale with effective collection area to get the event rate (N/s)
    1099     // scale with the effective observation time to absolute observed events
    1100     n *= area.GetCollectionAreaAbs()*ontime; // N
    1101 
    1102     // Now calculate the scale factor from the number of events
    1103     // produced and the number of events which should have been
    1104     // observed with our telescope in the time ontime
    1105     const Double_t scale = n/area.GetEntries();
    1106 
    1107     // Print normalization constant
    1108     cout << "MC normalization factor:  " << scale << endl;
    1109 
    1110     // Overlay normalized plots
    1111     DisplaySize(plist, scale);
     1307    TH1D weights;
     1308    TH1D energy0;
     1309    weights.SetDirectory(0);
     1310
     1311    const Int_t n = fWeighting ? 2 : 1;
     1312    for (int i=0; i<n; i++)
     1313    {
     1314        *fLog << endl;
     1315        fLog->Separator("Calculate efficiencies");
     1316        *fLog << endl;
     1317
     1318        if (i==1)
     1319            weight.SetWeightsSize(&weights);
     1320
     1321
     1322        MTaskList tlist2;
     1323        plist.Replace(&tlist2);
     1324
     1325        MReadMarsFile read("Events");
     1326        read.DisableAutoScheme();
     1327        set.AddFilesOn(read);
     1328
     1329        // Selector to get correct (final) theta-distribution
     1330        temp1.SetXTitle("MPointingPos.fZd");
     1331
     1332        MH3 mh3(temp1);
     1333
     1334        MFEventSelector2 sel2(mh3);
     1335        sel2.SetHistIsProbability();
     1336
     1337        MContinue contsel(&sel2);
     1338        contsel.SetInverted();
     1339
     1340//***        // Get correct source position
     1341//***        MSrcPosCalc calc;
     1342//***
     1343//***        // Calculate corresponding Hillas parameters
     1344//***        MHillasCalc hcalc1;
     1345//***        MHillasCalc hcalc2("MHillasCalcAnti");
     1346//***        hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
     1347//***        hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
     1348//***        hcalc2.SetNameHillasSrc("MHillasSrcAnti");
     1349//***        hcalc2.SetNameSrcPosCam("MSrcPosAnti");
     1350
     1351        // Fill collection area and energy estimator (unfolding)
     1352        // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
     1353        MHCollectionArea area;
     1354        area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
     1355        MHEnergyEst      hest;
     1356
     1357        MFillH fill3(&area, "", "FillCollectionArea");
     1358        MFillH fill4(&hest, "", "FillEnergyEst");
     1359        MFillH fill5("MHThreshold", "", "FillThreshold");
     1360        fill3.SetWeight();
     1361        fill4.SetWeight();
     1362        fill5.SetWeight();
     1363        fill3.SetNameTab("ColArea");
     1364        fill4.SetNameTab("E-Est");
     1365        fill5.SetNameTab("Threshold");
     1366
     1367        MH3 hsize("MHillas.fSize");
     1368        hsize.SetName("ExcessMC");
     1369        hsize.Sumw2();
     1370
     1371        MH3 energy("MMcEvt.fEnergy");
     1372        energy.SetName("EnergyEst"); // FIXME: Should be EnergyMC, but binning has name EnergyEst
     1373        energy.SetLogy();
     1374        energy.Sumw2();
     1375
     1376        MBinning bins(size, "BinningExcessMC");
     1377        plist.AddToList(&hsize);
     1378        plist.AddToList(&energy);
     1379        plist.AddToList(&bins);
     1380
     1381        MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
     1382        MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
     1383        MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
     1384        MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
     1385        MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
     1386        MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
     1387        MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
     1388        MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
     1389        MFillH fill9a("EnergyEst      [MH3]",           "",             "FillEnergy");
     1390        fill1a.SetNameTab("PreCut");
     1391        fill2a.SetNameTab("PostCut");
     1392        fill3a.SetNameTab("VsSize");
     1393        fill4a.SetNameTab("HilExt");
     1394        fill5a.SetNameTab("HilSrc");
     1395        fill6a.SetNameTab("ImgPar");
     1396        fill7a.SetNameTab("NewPar");
     1397        fill9a.SetNameTab("Energy");
     1398        fill8a.SetBit(MFillH::kDoNotDisplay);
     1399        fill9a.SetBit(MFillH::kDoNotDisplay);
     1400        fill1a.SetWeight();
     1401        fill2a.SetWeight();
     1402        fill3a.SetWeight();
     1403        fill4a.SetWeight();
     1404        fill5a.SetWeight();
     1405        fill6a.SetWeight();
     1406        fill7a.SetWeight();
     1407        fill8a.SetWeight();
     1408        fill9a.SetWeight();
     1409
     1410        MEnergyEstimate est;
     1411        MTaskEnv taskenv1("EstimateEnergy");
     1412        taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
     1413
     1414        if (i==1)
     1415        {
     1416            est.SetName("SetEnergyMC");
     1417            est.SetRule("MMcEvt.fEnergy");
     1418        }
     1419
     1420        tlist2.AddToList(&read);
     1421        if (!fRawMc && fNoThetaWeights)
     1422            tlist2.AddToList(&contsel);
     1423//***        tlist2.AddToList(&calc);
     1424//***        tlist2.AddToList(&hcalc1);
     1425//***        tlist2.AddToList(&hcalc2);
     1426        tlist2.AddToList(&weight);
     1427        if (i==n-1)
     1428            tlist2.AddToList(&fill1a);
     1429        tlist2.AddToList(fCut0);
     1430        tlist2.AddToList(fCut1);
     1431        tlist2.AddToList(fCut2);
     1432        tlist2.AddToList(fCut3);
     1433        tlist2.AddToList(i==1 ? (MTask*)&est : (MTask*)&taskenv1);
     1434        tlist2.AddToList(&fill3);
     1435        if (i==0)
     1436            tlist2.AddToList(&fill4);
     1437        tlist2.AddToList(&fill5);
     1438        if (i==n-1)
     1439        {
     1440            tlist2.AddToList(&fill2a);
     1441            tlist2.AddToList(&fill3a);
     1442            tlist2.AddToList(&fill4a);
     1443            tlist2.AddToList(&fill5a);
     1444            tlist2.AddToList(&fill6a);
     1445            tlist2.AddToList(&fill7a);
     1446        }
     1447        tlist2.AddToList(&fill8a);
     1448        tlist2.AddToList(&fill9a);
     1449
     1450        MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
     1451        loop2.SetParList(&plist);
     1452        loop2.SetDisplay(fDisplay);
     1453        loop2.SetLogStream(fLog);
     1454
     1455//***<<<<<<< MJSpectrum.cc
     1456        if (!SetupEnv(loop2))
     1457            return kFALSE;
     1458//***=======
     1459//***    MFEventSelector2 sel2(mh3);
     1460//***    sel2.SetHistIsProbability();
     1461//***
     1462//***    MContinue contsel(&sel2);
     1463//***    contsel.SetInverted();
     1464//***
     1465//***    // Get correct source position
     1466//***    //MSrcPosCalc calc;
     1467//***
     1468//***    // Calculate corresponding Hillas parameters
     1469//***    /*
     1470//***     MHillasCalc hcalc1;
     1471//***     MHillasCalc hcalc2("MHillasCalcAnti");
     1472//***     hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
     1473//***     hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
     1474//***     hcalc2.SetNameHillasSrc("MHillasSrcAnti");
     1475//***     hcalc2.SetNameSrcPosCam("MSrcPosAnti");
     1476//***     */
     1477//***
     1478//***    // Fill collection area and energy estimator (unfolding)
     1479//***    // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
     1480//***    MHCollectionArea area;
     1481//***    area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
     1482//***    MHEnergyEst      hest;
     1483//***
     1484//***    MFillH fill3(&area, "", "FillCollectionArea");
     1485//***    MFillH fill4(&hest, "", "FillEnergyEst");
     1486//***    MFillH fill5("MHThreshold", "", "FillThreshold");
     1487//***    fill3.SetWeight();
     1488//***    fill4.SetWeight();
     1489//***    fill5.SetWeight();
     1490//***    fill3.SetNameTab("ColArea");
     1491//***    fill4.SetNameTab("E-Est");
     1492//***    fill5.SetNameTab("Threshold");
     1493//***
     1494//***    MH3 hsize("MHillas.fSize");
     1495//***    hsize.SetName("ExcessMC");
     1496//***    hsize.Sumw2();
     1497//***
     1498//***    MBinning bins(size, "BinningExcessMC");
     1499//***    plist.AddToList(&hsize);
     1500//***    plist.AddToList(&bins);
     1501//***
     1502//***    MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
     1503//***    MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
     1504//***    MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
     1505//***    MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
     1506//***    MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
     1507//***    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
     1508//***    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
     1509//***    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
     1510//***    fill1a.SetNameTab("PreCut");
     1511//***    fill2a.SetNameTab("PostCut");
     1512//***    fill3a.SetNameTab("VsSize");
     1513//***    fill4a.SetNameTab("HilExt");
     1514//***    fill5a.SetNameTab("HilSrc");
     1515//***    fill6a.SetNameTab("ImgPar");
     1516//***    fill7a.SetNameTab("NewPar");
     1517//***    fill8a.SetBit(MFillH::kDoNotDisplay);
     1518//***    fill1a.SetWeight();
     1519//***    fill2a.SetWeight();
     1520//***    fill3a.SetWeight();
     1521//***    fill4a.SetWeight();
     1522//***    fill5a.SetWeight();
     1523//***    fill6a.SetWeight();
     1524//***    fill7a.SetWeight();
     1525//***    fill8a.SetWeight();
     1526//***>>>>>>> 1.17
     1527
     1528        if (!loop2.Eventloop(fMaxEvents))
     1529        {
     1530            *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
     1531            return kFALSE;
     1532        }
     1533
     1534//***<<<<<<< MJSpectrum.cc
     1535        if (!loop2.GetDisplay())
     1536        {
     1537            *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
     1538            return kFALSE;
     1539        }
     1540//***=======
     1541//***    tlist2.AddToList(&read);
     1542//***    if (!fRawMc && fNoThetaWeights)
     1543//***        tlist2.AddToList(&contsel);
     1544//***    //tlist2.AddToList(&calc);
     1545//***    //tlist2.AddToList(&hcalc1);
     1546//***    //tlist2.AddToList(&hcalc2);
     1547//***    tlist2.AddToList(&weight);
     1548//***    tlist2.AddToList(&fill1a);
     1549//***    tlist2.AddToList(fCut0);
     1550//***    tlist2.AddToList(fCut1);
     1551//***    tlist2.AddToList(fCut2);
     1552//***    tlist2.AddToList(fCut3);
     1553//***    tlist2.AddToList(&taskenv1);
     1554//***    tlist2.AddToList(&fill3);
     1555//***    tlist2.AddToList(&fill4);
     1556//***    tlist2.AddToList(&fill5);
     1557//***    tlist2.AddToList(&fill2a);
     1558//***    tlist2.AddToList(&fill3a);
     1559//***    tlist2.AddToList(&fill4a);
     1560//***    tlist2.AddToList(&fill5a);
     1561//***    tlist2.AddToList(&fill6a);
     1562//***    tlist2.AddToList(&fill7a);
     1563//***    tlist2.AddToList(&fill8a);
     1564//***    //tlist2.AddToList(&fill9a);
     1565//***
     1566//***    MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
     1567//***    loop2.SetParList(&plist);
     1568//***    loop2.SetDisplay(fDisplay);
     1569//***    loop2.SetLogStream(fLog);
     1570//***>>>>>>> 1.17
     1571
     1572        gLog.Separator("Energy Estimator");
     1573        if (plist.FindObject("EstimateEnergy"))
     1574            plist.FindObject("EstimateEnergy")->Print();
     1575
     1576        gLog.Separator("Spectrum");
     1577
     1578        // -------------------------- Spectrum ----------------------------
     1579
     1580
     1581        TArrayD res;
     1582        if (i==0)
     1583        {
     1584            // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
     1585            res = DisplaySpectrum(area, excess, hest, ontime);
     1586
     1587            ((TH1D&)energy.GetHist()).Copy(energy0);
     1588        }
     1589        else
     1590        {
     1591            const Double_t scale = ontime*area.GetCollectionAreaAbs();
     1592            if (!DisplaySpectrumW(energy0, energy.GetHist(), hsize.GetHist(), weights, fSimpleMode ? hist : (TH2D&)mh1.GetHist(), scale, res))
     1593                return kFALSE;
     1594        }
     1595
     1596        // Spectrum fitted (convert res[1] from TeV to GeV)
     1597        TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
     1598
     1599        // Number of events this spectrum would produce per s and m^2
     1600        Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
     1601
     1602        // scale with effective collection area to get the event rate (N/s)
     1603        // scale with the effective observation time to absolute observed events
     1604        n *= area.GetCollectionAreaAbs()*ontime; // N
     1605
     1606        // Now calculate the scale factor from the number of events
     1607        // produced and the number of events which should have been
     1608        // observed with our telescope in the time ontime
     1609        const Double_t scale = i==0 ? n/area.GetEntries() : 1;
     1610
     1611        // Print normalization constant
     1612        cout << "MC normalization factor:  " << scale << endl;
     1613
     1614        // Overlay normalized plots, and calculate weights for second loop
     1615        DisplaySize(plist, scale, weights);
     1616    }
    11121617
    11131618    // check if output should be written

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

    r7148 r7403  
    66#endif
    77
     8class TF1;
    89class TH1;
    910class TH1D;
    … …  
    3334    Bool_t fRawMc;
    3435    Bool_t fNoThetaWeights;
     36    Bool_t fWeighting;
    3537
    3638    // Read Input
    … …  
    4648    void    DisplayResult(const TH2D &mh1) const;
    4749    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const;
     50    TArrayD FitSpectrum(TH1D &spectrum) const;
    4851    TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
    49     Bool_t  DisplaySize(MParList &plist, Double_t scale) const;
     52    Bool_t  DisplaySpectrumW(TH1D &energy0, TH1 &energy, TH1 &hsize, TH1D &weights, TH2D &hall, Double_t scale, TArrayD &res) const;
     53    Bool_t  DisplaySize(MParList &plist, Double_t scale, TH1D &result) const;
    5054    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const;
    5155
    … …  
    5963    void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; }
    6064    void EnableRawMc(Bool_t b=kTRUE)      { fRawMc=b; }
     65    void EnableWeighting(Bool_t b=kTRUE)  { fWeighting=b; }
    6166
    6267    void SetEnergyEstimator(const MTask *task);
     68
     69    static TString FormString(const TF1 &f, Byte_t type=0);
    6370
    6471    ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum
Note: See TracChangeset for help on using the changeset viewer.