Ignore:
Timestamp:
11/15/05 17:06:04 (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

    r7403 r7404  
    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), fWeighting(kTRUE)
     89    fNoThetaWeights(kFALSE)
    9090{
    9191    fName  = name  ? name  : "MJSpectrum";
    … …  
    601601    f.DrawCopy("same");
    602602
    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     */
    613603    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     */
     604
    619605    TLatex tex;
    620606    tex.SetTextSize(0.045);
    … …  
    652638
    653639    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);
    657640
    658641    spectrum.SetDirectory(NULL);
    … …  
    689672    weights.DrawCopy();
    690673
    691     //spectrum.Multiply(&weights); // Currentyl we skip the error distribution!
     674    //spectrum.Multiply(&weights);
    692675    spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
    693676    spectrum.SetBit(TH1::kNoStats);
    … …  
    695678    for (int i=0; i<excess.GetNbinsX(); i++)
    696679    {
    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         }
     680        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
     681        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
    706682
    707683        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
    708         spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  /spectrum.GetBinWidth(i+1)*1000);
     684        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
    709685    }
    710686
    … …  
    726702        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
    727703    }
     704
     705    // Display dN/dE*E^2 for conveinience
    728706    fDisplay->AddTab("Check");
    729707    gPad->SetLogx();
    … …  
    741719    c1.cd(4);
    742720    return FitSpectrum(*spec);
    743     /*
     721
     722/*
    744723    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
    745724    f.SetParameter(0, -2.87);
    746725    f.SetParameter(1, 1.9e-6);
    747     f.SetLineColor(kBlue);
     726    f.SetLineColor(kGreen);
    748727    spectrum.Fit(&f, "NIM", "", 100, 5000);
    749728    f.DrawCopy("same");
    … …  
    773752    TArrayD res(2);
    774753    res[0] = f.GetParameter(0);
    775     res[1] = f.GetParameter(1);*/
    776 
     754    res[1] = f.GetParameter(1);
     755
     756    return res;
     757  */
    777758/*
    778759    // Plot other spectra  from Whipple
    … …  
    807788    f.SetLineStyle(kDashed);
    808789    f.DrawCopy("same");
    809     return res;
    810790 */
    811791}
    … …  
    870850// Calls PlotSame
    871851//
    872 //#include <TSpline.h>
    873 Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale, TH1D &result) const
     852Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
    874853{
    875854    *fLog << inf << "Reading from file: " << fPathIn << endl;
    … …  
    890869    }
    891870
    892     TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
     871    TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
    893872    if (!excess)
    894873        return kFALSE;
    … …  
    906885
    907886    excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
    908     excess = (TH1D*)excess->DrawCopy("E2");
     887    excess = excess->DrawCopy("E2");
    909888    // Don't do this on the original object!
    910889    excess->SetStats(kFALSE);
    … …  
    918897    if ((o=plist.FindObject("ExcessMC")))
    919898    {
    920         TH1 *histsel = (TH1D*)o->FindObject("");
     899        TH1 *histsel = (TH1F*)o->FindObject("");
    921900        if (histsel)
    922901        {
    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
    988902            if (scale<0)
    989903                scale = excess->Integral()/histsel->Integral();
    … …  
    996910            histsel->SetStats(kFALSE);
    997911
    998             // Do some Chi^2 and Kolmogorov tests
    999912            fLog->Separator("Kolmogorov Test");
    1000913            histsel->KolmogorovTest(excess, "DX");
    … …  
    1002915            const Double_t p = histsel->Chi2Test(excess, "P");
    1003916
    1004             // Plot result
    1005917            TLatex tex;
    1006918            tex.SetBit(TLatex::kTextNDC);
    … …  
    1024936    c.cd(6);
    1025937    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);
    1026 
    1027     return kTRUE;
    1028 }
    1029 
    1030 Bool_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     }
    1201938
    1202939    return kTRUE;
    … …  
    13051042    // ------------------------- Final loop --------------------------
    13061043
    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     }
     1044    *fLog << endl;
     1045    fLog->Separator("Calculate efficiencies");
     1046    *fLog << endl;
     1047
     1048    MTaskList tlist2;
     1049    plist.Replace(&tlist2);
     1050
     1051    MReadMarsFile read("Events");
     1052    read.DisableAutoScheme();
     1053    set.AddFilesOn(read);
     1054
     1055    // Selector to get correct (final) theta-distribution
     1056    temp1.SetXTitle("MPointingPos.fZd");
     1057
     1058    MH3 mh3(temp1);
     1059
     1060    MFEventSelector2 sel2(mh3);
     1061    sel2.SetHistIsProbability();
     1062
     1063    MContinue contsel(&sel2);
     1064    contsel.SetInverted();
     1065
     1066    // Get correct source position
     1067    //MSrcPosCalc calc;
     1068
     1069    // Calculate corresponding Hillas parameters
     1070    /*
     1071     MHillasCalc hcalc1;
     1072     MHillasCalc hcalc2("MHillasCalcAnti");
     1073     hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
     1074     hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
     1075     hcalc2.SetNameHillasSrc("MHillasSrcAnti");
     1076     hcalc2.SetNameSrcPosCam("MSrcPosAnti");
     1077     */
     1078
     1079    // Fill collection area and energy estimator (unfolding)
     1080    // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
     1081    MHCollectionArea area;
     1082    area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
     1083    MHEnergyEst      hest;
     1084
     1085    MFillH fill3(&area, "", "FillCollectionArea");
     1086    MFillH fill4(&hest, "", "FillEnergyEst");
     1087    MFillH fill5("MHThreshold", "", "FillThreshold");
     1088    fill3.SetWeight();
     1089    fill4.SetWeight();
     1090    fill5.SetWeight();
     1091    fill3.SetNameTab("ColArea");
     1092    fill4.SetNameTab("E-Est");
     1093    fill5.SetNameTab("Threshold");
     1094
     1095    MH3 hsize("MHillas.fSize");
     1096    hsize.SetName("ExcessMC");
     1097    hsize.Sumw2();
     1098
     1099    MBinning bins(size, "BinningExcessMC");
     1100    plist.AddToList(&hsize);
     1101    plist.AddToList(&bins);
     1102
     1103    MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
     1104    MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
     1105    MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
     1106    MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
     1107    MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
     1108    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
     1109    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
     1110    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
     1111    fill1a.SetNameTab("PreCut");
     1112    fill2a.SetNameTab("PostCut");
     1113    fill3a.SetNameTab("VsSize");
     1114    fill4a.SetNameTab("HilExt");
     1115    fill5a.SetNameTab("HilSrc");
     1116    fill6a.SetNameTab("ImgPar");
     1117    fill7a.SetNameTab("NewPar");
     1118    fill8a.SetBit(MFillH::kDoNotDisplay);
     1119    fill1a.SetWeight();
     1120    fill2a.SetWeight();
     1121    fill3a.SetWeight();
     1122    fill4a.SetWeight();
     1123    fill5a.SetWeight();
     1124    fill6a.SetWeight();
     1125    fill7a.SetWeight();
     1126    fill8a.SetWeight();
     1127
     1128    MEnergyEstimate est;
     1129    MTaskEnv taskenv1("EstimateEnergy");
     1130    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
     1131
     1132    tlist2.AddToList(&read);
     1133    if (!fRawMc && fNoThetaWeights)
     1134        tlist2.AddToList(&contsel);
     1135    //tlist2.AddToList(&calc);
     1136    //tlist2.AddToList(&hcalc1);
     1137    //tlist2.AddToList(&hcalc2);
     1138    tlist2.AddToList(&weight);
     1139    tlist2.AddToList(&fill1a);
     1140    tlist2.AddToList(fCut0);
     1141    tlist2.AddToList(fCut1);
     1142    tlist2.AddToList(fCut2);
     1143    tlist2.AddToList(fCut3);
     1144    tlist2.AddToList(&taskenv1);
     1145    tlist2.AddToList(&fill3);
     1146    tlist2.AddToList(&fill4);
     1147    tlist2.AddToList(&fill5);
     1148    tlist2.AddToList(&fill2a);
     1149    tlist2.AddToList(&fill3a);
     1150    tlist2.AddToList(&fill4a);
     1151    tlist2.AddToList(&fill5a);
     1152    tlist2.AddToList(&fill6a);
     1153    tlist2.AddToList(&fill7a);
     1154    tlist2.AddToList(&fill8a);
     1155    //tlist2.AddToList(&fill9a);
     1156
     1157    MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
     1158    loop2.SetParList(&plist);
     1159    loop2.SetDisplay(fDisplay);
     1160    loop2.SetLogStream(fLog);
     1161
     1162    if (!SetupEnv(loop2))
     1163        return kFALSE;
     1164
     1165    if (!loop2.Eventloop(fMaxEvents))
     1166    {
     1167        *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
     1168        return kFALSE;
     1169    }
     1170
     1171    if (!loop2.GetDisplay())
     1172    {
     1173        *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
     1174        return kFALSE;
     1175    }
     1176
     1177    gLog.Separator("Energy Estimator");
     1178    if (plist.FindObject("EstimateEnergy"))
     1179        plist.FindObject("EstimateEnergy")->Print();
     1180
     1181    gLog.Separator("Spectrum");
     1182
     1183    // -------------------------- Spectrum ----------------------------
     1184
     1185    // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
     1186    TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
     1187
     1188    // Spectrum fitted (convert res[1] from TeV to GeV)
     1189    TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
     1190
     1191    // Number of events this spectrum would produce per s and m^2
     1192    Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
     1193
     1194    // scale with effective collection area to get the event rate (N/s)
     1195    // scale with the effective observation time to absolute observed events
     1196    n *= area.GetCollectionAreaAbs()*ontime; // N
     1197
     1198    // Now calculate the scale factor from the number of events
     1199    // produced and the number of events which should have been
     1200    // observed with our telescope in the time ontime
     1201    const Double_t scale = n/area.GetEntries();
     1202
     1203    // Print normalization constant
     1204    cout << "MC normalization factor:  " << scale << endl;
     1205
     1206    // Overlay normalized plots
     1207    DisplaySize(plist, scale);
    16171208
    16181209    // check if output should be written

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

    r7403 r7404  
    3434    Bool_t fRawMc;
    3535    Bool_t fNoThetaWeights;
    36     Bool_t fWeighting;
    3736
    3837    // Read Input
    … …  
    5049    TArrayD FitSpectrum(TH1D &spectrum) const;
    5150    TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) 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;
     51    Bool_t  DisplaySize(MParList &plist, Double_t scale) const;
    5452    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const;
    5553
    … …  
    6361    void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; }
    6462    void EnableRawMc(Bool_t b=kTRUE)      { fRawMc=b; }
    65     void EnableWeighting(Bool_t b=kTRUE)  { fWeighting=b; }
    6663
    6764    void SetEnergyEstimator(const MTask *task);
Note: See TracChangeset for help on using the changeset viewer.