Changeset 8047 for trunk/MagicSoft


Ignore:
Timestamp:
10/11/06 08:55:36 (18 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r8046 r8047  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20 2006/10/11 Thomas Bretz
     21
     22   * mhbase/MH.[h,cc]:
     23     - added a function to calculate binomial errors including weights
     24       (this was added in root 5.13/04, but necessary for older versions)
     25
     26   * mhflux/MHCollectionArea.[h,cc]:
     27     - added Sumw2() to the constructor so that the weights array gets
     28       correctly initialize
     29     - replaced the calculation of the binomial errors by the
     30       corresponding root-function and the new MH function
     31     - made sure that in all histogram operations the errors are
     32       properly propagated
     33     - let ReInit determine fMcRadius from MMcRunHeader
     34
     35   * mhflux/MHEnergyEst.cc, mhflux/MHThreshold.cc
     36     - fixed the order in the constructor such that the Sumw2() does
     37       correctly initialize the weights array
     38
     39   * mhflux/MMcSpectrumWeight.cc:
     40     - a minor code reordering
     41
     42   * mjobs/MJSpectrum.cc:
     43     - made sure that the histogram with the corsika spectrum has
     44       the errors initialized and thus takes the weights correctly
     45       into account
     46     - corresponding to this changed some draw option to get the
     47       same plots (hist) as before
     48     - added a lot of comments to the code
     49     - when the zenith angle weights are applied to the MC distribution
     50       make sure that also the errors are correctly treated.
     51
     52
     53
    2054 2006/10/10 Thomas Bretz
    2155

  • trunk/MagicSoft/Mars/NEWS

    r8042 r8047  
    33 *** Version  <cvs>
    44
     5   - sponde: In the calculation of the collection area(s) and the
     6     distribution for MOnte Carlo and estimated energy the error
     7     calculation was wrong because root didn't take the errors
     8     properly into account... fixed.
    59
    610

  • trunk/MagicSoft/Mars/mhbase/MH.cc

    r7200 r8047  
    13161316// --------------------------------------------------------------------------
    13171317//
     1318// This is a workaround for rrot-version <5.13/04 to get correct
     1319// binomial errors even if weights have been used, do:
     1320//    h->Divide(h1, h2, 5, 2, "b");
     1321//    MH::SetBinomialErrors(*h, *h1, *h2, 5, 2);
     1322//
     1323// see http://root.cern.ch/phpBB2/viewtopic.php?p=14818
     1324//
     1325void MH::SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1, Double_t c2)
     1326{
     1327    for (Int_t binx=0; binx<=hres.GetNbinsX()+1; binx++)
     1328    {
     1329        const Double_t b1 = h1.GetBinContent(binx);
     1330        const Double_t b2 = h2.GetBinContent(binx);
     1331        const Double_t w  = c2*b2 ? (c1*b1)/(c2*b2) : 0;
     1332        const Double_t e1 = h2.GetBinError(binx);
     1333        const Double_t e2 = h1.GetBinError(binx);
     1334
     1335        const Double_t rc = ((1-2*w)*e1*e1+w*w*e2*e2)/(b2*b2);
     1336
     1337        hres.SetBinError(binx, TMath::Sqrt(TMath::Abs(rc)));
     1338    }
     1339}
     1340
     1341// --------------------------------------------------------------------------
     1342//
    13181343// See MTask::PrintSkipped
    13191344//

  • trunk/MagicSoft/Mars/mhbase/MH.h

    r7173 r8047  
    7676    static void SetBinning(TH1 *h, const TH1 *x);
    7777
     78    static void SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1=1, Double_t c2=1);
     79
    7880    static void RemoveFirstBin(TH1 &h);
    7981

  • trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc

    r7684 r8047  
    6464//
    6565MHCollectionArea::MHCollectionArea(const char *name, const char *title)
    66   : fMcEvt(0), /*fEnergy(0),*/ fMcAreaRadius(300.), fIsExtern(kFALSE)
     66  : fMcEvt(0), fMcAreaRadius(-1), fIsExtern(kFALSE)
    6767{
    6868    //   initialize the histogram for the distribution r vs E
    … …  
    8282    fHistSel.SetDirectory(NULL);
    8383    fHistSel.UseCurrentStyle();
     84    fHistSel.SetLineColor(kBlue);
    8485
    8586    fHistAll.SetName("AllEvts");
    … …  
    105106    MH::SetBinning(&fHistSel, &binst, &binse);
    106107    MH::SetBinning(&fHistAll, &binst, &binse);
     108
     109    // For some unknown reasons this must be called after
     110    // the binning has been initialized at least once
     111    fHistSel.Sumw2();
     112    fHistAll.Sumw2();
     113    fHEnergy.Sumw2();
    107114}
    108115
    … …  
    114121void MHCollectionArea::CalcEfficiency()
    115122{
    116     TH1D *hsel = fHistSel.ProjectionY();
    117     TH1D *hall = fHistAll.ProjectionY();
    118 
    119     const Int_t nbinx = hsel->GetNbinsX();
    120 
    121     // -----------------------------------------------------------
    122     //
    123     // Impact parameter range:  TO BE FIXED! Impact parameter shoud be
    124     // read from run header, but it is not yet in!!
     123    TH1D *hsel = fHistSel.ProjectionY("Spy", -1, 9999, "E");;
     124    TH1D *hall = fHistAll.ProjectionY("Apy", -1, 9999, "E");
     125
     126    //
     127    // Impact parameter range.
    125128    //
    126129    const Float_t totalarea = GetCollectionAreaAbs();//TMath::Pi() * (r2*r2 - r1*r1);
    127130
    128     for (Int_t ix=1; ix<=nbinx; ix++)
    129     {
    130         const Float_t Na = hall->GetBinContent(ix);
    131 
    132         if (Na <= 0)
    133             continue;
    134 
    135         const Float_t Ns = hsel->GetBinContent(ix);
    136 
    137         // Since Na is an estimate of the total number of showers generated
    138         // in the energy bin, it may happen that Ns (triggered showers) is
    139         // larger than Na. In that case, the bin is skipped:
    140 
    141         if (Na < Ns)
    142             continue;
    143 
    144         const Double_t eff         = Ns/Na;
    145         const Double_t efferr      = TMath::Sqrt((1.-eff)*Ns)/Na;
    146        
    147         const Float_t colarea      =  eff    * totalarea;
    148         const Float_t colareaerror =  efferr * totalarea;
    149 
    150         fHEnergy.SetBinContent(ix, colarea);
    151         fHEnergy.SetBinError(ix,   colareaerror);
    152     }
     131    // "b" option: calculate binomial errors
     132    fHEnergy.Divide(hsel, hall, totalarea, 1, "b");
     133#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
     134    MH::SetBinomialErrors(fHEnergy, *hsel, *hall, totalarea, 1);
     135#endif
    153136
    154137    delete hsel;
    155138    delete hall;
    156139}
     140
    157141
    158142Bool_t MHCollectionArea::SetupFill(const MParList *pl)
    … …  
    196180    MH::SetBinning(&fHistAll, &binst, &binse);
    197181
     182    fMcAreaRadius = -1;
     183
    198184    return kTRUE;
    199185}
    … …  
    201187Bool_t MHCollectionArea::ReInit(MParList *plist)
    202188{
     189    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
     190    if (!runheader)
     191    {
     192        *fLog << err << "MMcRunHeader not found... abort." << endl;
     193        return kFALSE;
     194    }
     195
     196    // FIXME: Does this need some weighting with the number of produced events?
     197    if (runheader->GetImpactMax()>fMcAreaRadius*100)
     198    {
     199        fMcAreaRadius = 0.01*runheader->GetImpactMax(); // cm->m
     200        *fLog << inf << "Maximum simulated impact: " << fMcAreaRadius << "m" << endl;
     201    }
     202
    203203    if (fIsExtern)
    204204        return kTRUE;
    205205
    206     MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
    207     if (!runheader)
    208     {
    209         *fLog << err << "MMcRunHeader not found... abort." << endl;
    210         return kFALSE;
    211     }
    212 
    213206    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
    214207    *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;
    215208
    216  /*
    217     if (fTheta>=0 && fTheta!=runheader->GetTelesTheta())
    218     {
    219         *fLog << warn << "Warning - Read files have different TelesTheta (";
    220         *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl;
    221     }
    222     fTheta = runheader->GetTelesTheta();
    223   */
    224209    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
    225210        *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
    … …  
    352337    if (h1 && h2 && h)
    353338    {
    354         h->Divide(h2, h1);
    355         h->SetMinimum(0);
     339            h->Divide(h2, h1, 1, 1, "b");
     340#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
     341            MH::SetBinomialErrors(*h, *h2, *h1);
     342#endif
     343            h->SetMinimum(0);
    356344    }
    357345
    … …  
    433421        gPad->SetLogx();
    434422        h = h2->DrawCopy();
    435         h->Divide(h1);
     423        h->Divide(h2, h1, 1, 1, "b");
     424#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04)
     425        MH::SetBinomialErrors(*h, *h2, *h1);
     426#endif
    436427        h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency");
    437428        h->SetDirectory(NULL);
    … …  
    456447Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight)
    457448{
    458     const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/;
     449    const Double_t energy = fMcEvt->GetEnergy();
    459450    const Double_t theta  = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
    460451
    461     //*fLog << energy << " " << theta << endl;
    462 
    463452    fHistSel.Fill(theta, energy, weight);
    464453

  • trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h

    r7683 r8047  
    6060    Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; }
    6161
    62 /*
    63     void DrawAll(Option_t *option="");
    64     void DrawSel(Option_t *option="");
    65 
    66     const TH1D *GetHist()       { return fHistCol; }
    67     const TH1D *GetHist() const { return fHistCol; }
    68 
    69 
    70 
    71     void CalcEfficiency();
    72     void CalcEfficiency2();
    73 
    74     void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall);
    75     void Calc(const MHMcEfficiency &heff);
    76   */
    7762    void SetMcAreaRadius(Float_t x) { fMcAreaRadius = x; }
    7863

  • trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

    r8021 r8047  
    9494    fHImpact.SetYTitle("E_{est}/E_{mc}-1");
    9595
    96     fHEnergy.Sumw2();
    97     fHImpact.Sumw2();
    98     fHResolution.Sumw2();
    99 
    10096    MBinning binsi, binse, binst, binsr;
    10197    binse.SetEdgesLog(25, 10, 1000000);
    … …  
    107103    SetBinning(&fHImpact,     &binsi, &binsr);
    108104    SetBinning(&fHResolution, &binse, &binse, &binsr);
     105
     106    // For some unknown reasons this must be called after
     107    // the binning has been initialized at least once
     108    fHEnergy.Sumw2();
     109    fHResolution.Sumw2();
     110    fHImpact.Sumw2();
    109111}
    110112
    … …  
    243245{
    244246    // Project into EnergyEst_ey
     247    // the "e" ensures that errors are calculated
    245248    TH1D *h1 = (TH1D*)fHEnergy.Project3D("rtlprmft_ex"); // E_Est
    246249    TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey"); // E_mc

  • trunk/MagicSoft/Mars/mhflux/MHThreshold.cc

    r7715 r8047  
    7979    fHEnergy.SetDirectory(NULL);
    8080    fHEnergy.UseCurrentStyle();
    81     fHEnergy.Sumw2();
    8281
    8382    MBinning binse(50, 20, 50000, "", "log");
    8483    binse.Apply(fHEnergy);
     84
     85    fHEnergy.Sumw2();
    8586}
    8687

  • trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc

    r8036 r8047  
    443443Int_t MMcSpectrumWeight::Process()
    444444{
    445     const Double_t e = fMcEvt->GetEnergy();
    446 
    447445    Double_t w = 1;
    448446
    … …  
    459457    }
    460458
     459    const Double_t e = fMcEvt->GetEnergy();
    461460    fWeight->SetVal(fFunc->Eval(e)*w);
    462461

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

    r7650 r8047  
    176176    fLog->Separator("Initialize energy weighting");
    177177
     178    // Read Resources
    178179    if (!CheckEnv(w))
    179180    {
    … …  
    182183    }
    183184
     185    // Read the first corsika RunHeader from the MC files
    184186    TChain chain("RunHeaders");
    185187    set.AddFilesOn(chain);
    … …  
    195197    }
    196198
     199    // Propagate the run header to MMcSpectrumWeight
    197200    if (!w.Set(*h))
    198201    {
    … …  
    201204    }
    202205
     206    // Print new setup of MMcSpectrumWeight
    203207    w.Print();
    204208    return kTRUE;
    … …  
    266270    fLog->Separator("Compiling original MC distribution");
    267271
     272    // The name of the input container is MMcEvtBasic
    268273    weight.SetNameMcEvt("MMcEvtBasic");
     274    // Get the corresponding rule for the weights
    269275    const TString w(weight.GetFormulaWeights());
     276    // Set back to MMcEvt
    270277    weight.SetNameMcEvt();
    271278
    … …  
    282289    // Prepare histogram
    283290    h.Reset();
     291    h.Sumw2();
    284292
    285293    // Fill histogram from chain
    … …  
    331339        gPad->SetBorderMode(0);
    332340        temp2.SetName("NVsTheta");
    333         temp2.DrawCopy();
     341        temp2.DrawCopy("hist");
    334342
    335343        c.cd(4);
    … …  
    349357    temp1.SetYTitle("Probability");
    350358    if (fDisplay)
    351         temp1.DrawCopy();
     359        temp1.DrawCopy("hist");
    352360
    353361    return kTRUE;
    … …  
    580588    case 0:
    581589        {
    582             const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",   f1, f1);
     590            const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",  f1, f1);
    583591            const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
    584592
    … …  
    638646TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
    639647{
     648    // Create copies of the histograms
    640649    TH1D collarea(area.GetHEnergy());
    641650    TH1D spectrum(excess);
    642651    TH1D weights;
     652
     653    // Get spill-over corrections from energy estimation
    643654    hest.GetWeights(weights);
    644655
     656    // Print effective on-time
    645657    cout << "Effective On time: " << ontime << "s" << endl;
    646658
     659    // Setup spectrum plot
     660    spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
     661    spectrum.SetYTitle("N/sm^{2}");
    647662    spectrum.SetDirectory(NULL);
    648663    spectrum.SetBit(kCanDelete);
     664
     665    // Divide by collection are and on-time
    649666    spectrum.Scale(1./ontime);
    650667    spectrum.Divide(&collarea);
    651     spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
    652     spectrum.SetYTitle("N/sm^{2}");
    653 
     668
     669    // Draw spectrum before applying spill-over corrections
    654670    TCanvas &c1 = fDisplay->AddTab("Spectrum");
    655671    c1.Divide(2,2);
    … …  
    679695    weights.DrawCopy();
    680696
     697    // Apply spill-over correction (done't take the errors into account)
     698    // They are supposed to be identical with the errors of the
     699    // collection area and cancel out.
    681700    //spectrum.Multiply(&weights);
    682701    spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
    … …  
    840859
    841860    //h2->Scale(1./ontime);   //h2->Integral());
    842     h3->Scale(scale);    //h3->Integral());
     861    h3->Scale(scale);         //h3->Integral());
    843862
    844863    h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
    … …  
    10011020    plist.AddToList(&bins2);
    10021021
     1022    // Initialize weighting to a new spectru
     1023    // m as defined in the resource file
    10031024    MMcSpectrumWeight weight;
    10041025    if (!InitWeighting(set, weight))
    10051026        return kFALSE;
    10061027
     1028    // Print the setup of the MAlphaFitter
    10071029    PrintSetup(fit);
    10081030    bins3.SetEdges(temp1, 'x');
    10091031
     1032    // Read the MC distribution as produced by corsika into
     1033    // temp2 (1D) and apply the weights previously determined
    10101034    TH1D temp2(temp1);
    10111035    if (!ReadOrigMCDistribution(set, temp2, weight))
    10121036        return kFALSE;
    10131037
     1038    // Calculate the weights according to the zenith angle distribution
     1039    // of the raw-data which have to be applied to the MC events
    10141040    if (!GetThetaDistribution(temp1, temp2))
    10151041        return kFALSE;
    10161042
     1043    // Tell the weighting task about the ZA-weights
    10171044    if (!fNoThetaWeights)
    10181045        weight.SetWeightsZd(&temp1);
    10191046
     1047    // Refill excess histogram to determin the excess events
    10201048    TH1D excess;
    10211049    if (!Refill(plist, excess))
    … …  
    10261054    if (fSimpleMode)
    10271055    {
     1056        // Now we read the MC distribution as produced by corsika again
     1057        // with the spectral weights applied into a histogram vs.
     1058        // zenith angle and energy
    10281059        hist.UseCurrentStyle();
    10291060        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
    … …  
    10311062            return kFALSE;
    10321063
     1064        // No we apply the the zenith-angle-weights to the corsika produced
     1065        // MC distribution. Unfortunately this must be done manually
     1066        // because we are multiplying column by column
    10331067        if (!fRawMc)
    10341068        {
    10351069            for (int y=0; y<hist.GetNbinsY(); y++)
    10361070                for (int x=0; x<hist.GetNbinsX(); x++)
     1071                {
    10371072                    hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
    1038             //hist.SetEntries(hist.Integral());
     1073                    hist.SetBinError(x, y,   hist.GetBinError(x, y)  *temp1.GetBinContent(x));
     1074                }
    10391075        }
    10401076    }
    10411077    else
    10421078    {
     1079        // This rereads the original MC spectrumand aaplies both
     1080        // weights spectral weights and ZA-weights.
    10431081        weight.SetNameMcEvt("MMcEvtBasic");
    10441082        if (!IntermediateLoop(plist, mh1, temp1, set, weight))
    … …  
    11461184
    11471185    tlist2.AddToList(&read);
     1186    // If no weighting should be applied but the events should
     1187    // be thrown away according to the theta distribution
     1188    // it is enabled here
    11481189    if (!fRawMc && fNoThetaWeights)
    11491190        tlist2.AddToList(&contsel);
Note: See TracChangeset for help on using the changeset viewer.