Ignore:
Timestamp:
11/18/02 11:38:31 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/macros/estimate.C

    r1633 r1634  
    6262    eest.Add("MHillasSrc");
    6363
     64    //
     65    // Use this to change the binnign of the histograms to CT1-style
     66    //
     67    Bool_t usect1 = kTRUE;
     68
     69    //
     70    // Here you set up the coefficients of the parametrization
     71    // (MEnergyEstParam)
     72    //
    6473    TArrayD fA(5);
    6574    fA[0] =  -2539; // [cm]
     
    7786    fB[5] =  4.14e-8; // [GeV]
    7887    fB[6] =    -0.06;
    79     /* fb5=size*width, ir 29.8% (RMS=23%)
    80      fB[0] =    -8.64; // [GeV]
    81      fB[1] =   -0.069; // [GeV]
    82      fB[2] =  0.00066; // [GeV]
    83      fB[3] =    0.033; // [GeV]
    84      fB[4] = 0.000226; // [GeV]
    85      fB[5] =  4.14e-8; // [GeV]
    86      fB[6] =    -0.06;*/
    87     /* size*dist, ir 30%
    88      fB[0] =   -32; // [GeV]
    89      fB[1] =  -0.0295; // [GeV]
    90      fB[2] =   0.00112; // [GeV]
    91      fB[3] =   0.0644; // [GeV]
    92      fB[4] = 1e-5; // [GeV]
    93      fB[5] =  6.3e-6; // [GeV]
    94      fB[6] =   -0.06;
    95      */
    96     /* size*width, ir 30%
    97      fB[0] =   -17; // [GeV]
    98      fB[1] =  -0.06; // [GeV]
    99      fB[2] =   0.0013; // [GeV]
    100      fB[3] =   0.039; // [GeV]
    101      fB[4] = 0.00036; // [GeV]
    102      fB[5] =  5.6e-6; // [GeV]
    103      fB[6] =   -0.06;
    104      */
    105     /* fb5=width, ir 31%
    106      fB[0] =   0.14; // [GeV]
    107      fB[1] =  -0.078; // [GeV]
    108      fB[2] =   0.0003; // [GeV]
    109      fB[3] =  -0.171; // [GeV]
    110      fB[4] = 0.00038; // [GeV]
    111      fB[5] =  3e-5; // [GeV]
    112      fB[6] =   -0.06;
    113      */
    114     /* size*length, ir  32%
    115      fB[0] =   -1.14; // [GeV]
    116      fB[1] =  -0.075; // [GeV]
    117      fB[2] =   0.0008; // [GeV]
    118      fB[3] =  -0.148; // [GeV]
    119      fB[4] = 0.00037; // [GeV]
    120      fB[5] =  8e-7; // [GeV]
    121      fB[6] =   -0.06;
    122      */
    123     /* ir/dist 35%
    124      fB[0] =     7.7; // [GeV]
    125      fB[1] =  -0.076; // [GeV]
    126      fB[2] =   0.011; // [GeV]
    127      fB[3] =  -0.175; // [GeV]
    128      fB[4] = 0.00038; // [GeV]
    129      fB[5] = -0.0001; // [GeV]
    130      fB[6] =   -0.06;
    131      */
    132 
    13388
    13489    eest.SetCoeffA(fA);
    13590    eest.SetCoeffB(fB);
    13691
    137     MH3 mh3e("MMcEvt.fEnergy",     "abs(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
    138     MH3 mh3i("MMcEvt.fImpact/100", "abs(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
     92    MH3 mh3e("MMcEvt.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
     93    MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
     94    MH3 mh3eo("MMcEvt.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
     95    MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
     96
     97    MH3 mh3e2("MEnergyEst.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
     98    MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
     99    MH3 mh3eo2("MEnergyEst.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
     100    MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
     101
     102    MH3 mhe("MMcEvt.fEnergy",     "MEnergyEst.fEnergy");
     103    MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100");
    139104
    140105    mh3e.SetName("HistEnergy");
    141106    mh3i.SetName("HistImpact");
     107    mh3eo.SetName("HistEnergyOffset");
     108    mh3io.SetName("HistImpactOffset");
     109
     110    mh3e2.SetName("HistEnergy");
     111    mh3i2.SetName("HistImpact");
     112    mh3eo2.SetName("HistEnergyOffset");
     113    mh3io2.SetName("HistImpactOffset");
     114
     115    mhe.SetName("HistEE");
     116    mhi.SetName("HistII");
    142117
    143118    MFillH hfille(&mh3e);
    144119    MFillH hfilli(&mh3i);
     120    MFillH hfilleo(&mh3eo);
     121    MFillH hfillio(&mh3io);
     122
     123    MFillH hfille2(&mh3e2);
     124    MFillH hfilli2(&mh3i2);
     125    MFillH hfilleo2(&mh3eo2);
     126    MFillH hfillio2(&mh3io2);
     127
     128    MFillH hfillee(&mhe);
     129    MFillH hfillii(&mhi);
    145130
    146131    MBinning binsex("BinningHistEnergyX");
     
    148133    MBinning binsix("BinningHistImpactX");
    149134    MBinning binsiy("BinningHistImpactY");
    150 
    151     binsex.SetEdgesLog(50, 10, 1e4);
    152     binsey.SetEdges(50, 0, 7);
    153 
    154     binsix.SetEdges(50, 0, 300);
    155     binsiy.SetEdges(50, 0, 3);
     135    MBinning binseox("BinningHistEnergyOffsetX");
     136    MBinning binseoy("BinningHistEnergyOffsetY");
     137    MBinning binsiox("BinningHistImpactOffsetX");
     138    MBinning binsioy("BinningHistImpactOffsetY");
     139    MBinning binseex("BinningHistEEX");
     140    MBinning binsiix("BinningHistIIX");
     141    MBinning binseey("BinningHistEEY");
     142    MBinning binsiiy("BinningHistIIY");
     143
     144    binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4);
     145    binsey.SetEdges(50, 0, 1.75);
     146    binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
     147    binseoy.SetEdges(50, -1.75, 1.75);
     148
     149    binsix.SetEdges(50, 0, usect1 ? 450 : 300);
     150    binsiy.SetEdges(50, 0, 1.75);
     151    binsiox.SetEdges(50, 0, usect1 ? 450 : 300);
     152    binsioy.SetEdges(50, -1.75, 1.75);
     153
     154    binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3);
     155    binseey.SetEdgesLog(50, usect1 ? 300 : 1,  usect1 ? 50000 : 2e3);
     156    binsiix.SetEdges(50, 0, usect1 ? 450 : 300);
     157    binsiiy.SetEdges(50, 0, usect1 ? 450 : 150);
    156158
    157159    plist.AddToList(&binsex);
     
    159161    plist.AddToList(&binsix);
    160162    plist.AddToList(&binsiy);
     163    plist.AddToList(&binseox);
     164    plist.AddToList(&binseoy);
     165    plist.AddToList(&binsiox);
     166    plist.AddToList(&binsioy);
     167    plist.AddToList(&binseex);
     168    plist.AddToList(&binseey);
     169    plist.AddToList(&binsiix);
     170    plist.AddToList(&binsiiy);
    161171
    162172    //
     
    170180    tlist2.AddToList(&hfille);
    171181    tlist2.AddToList(&hfilli);
     182    tlist2.AddToList(&hfilleo);
     183    tlist2.AddToList(&hfillio);
     184
     185    tlist2.AddToList(&hfille2);
     186    tlist2.AddToList(&hfilli2);
     187    tlist2.AddToList(&hfilleo2);
     188    tlist2.AddToList(&hfillio2);
     189
     190    tlist2.AddToList(&hfillee);
     191    tlist2.AddToList(&hfillii);
    172192
    173193    /*
     
    179199    // Create and setup the eventloop
    180200    //
     201    MProgressBar bar;
     202
    181203    MEvtLoop evtloop;
     204    evtloop.SetProgressBar(&bar);
    182205    evtloop.SetParList(&plist);
    183206
     
    190213    tlist.PrintStatistics();
    191214
    192     mh3i.DrawClone("PROFX");
    193     mh3e.DrawClone("PROFX");
     215    const TString text = "\\sqrt{<y>}=%.0f%%";
     216
     217    TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
     218    c->Divide(2,2);
     219    c->cd(1);
     220    mh3i.DrawClone("PROFXnonew");
     221    TLatex *t = new TLatex(180, 1.6, Form(text, sqrt(mh3i.GetHist().GetMean(2))*100));
     222    t->Draw();
     223    c->cd(2);
     224    mh3e.DrawClone("PROFXnonew");
     225    t = new TLatex(2.7, 1.6, Form(text, sqrt(mh3e.GetHist().GetMean(2))*100));
     226    t->Draw();
     227    c->cd(3);
     228    mh3io.DrawClone("PROFXnonew");
     229    c->cd(4);
     230    mh3eo.DrawClone("PROFXnonew");
     231
     232    c=new TCanvas("Est2", "Estimates vs. E_{est}");
     233    c->Divide(2,2);
     234    c->cd(1);
     235    mh3i2.DrawClone("PROFXnonew");
     236    c->cd(2);
     237    mh3e2.DrawClone("PROFXnonew");
     238    c->cd(3);
     239    mh3io2.DrawClone("PROFXnonew");
     240    c->cd(4);
     241    mh3eo2.DrawClone("PROFXnonew");
     242
     243    mhe.DrawClone("PROFX");
     244    mhi.DrawClone("PROFX");
    194245}
    195 
Note: See TracChangeset for help on using the changeset viewer.