Ignore:
Timestamp:
04/25/05 15:37:35 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

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