Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 8046) +++ trunk/MagicSoft/Mars/Changelog (revision 8047) @@ -18,4 +18,38 @@ -*-*- END OF LINE -*-*- + 2006/10/11 Thomas Bretz + + * mhbase/MH.[h,cc]: + - added a function to calculate binomial errors including weights + (this was added in root 5.13/04, but necessary for older versions) + + * mhflux/MHCollectionArea.[h,cc]: + - added Sumw2() to the constructor so that the weights array gets + correctly initialize + - replaced the calculation of the binomial errors by the + corresponding root-function and the new MH function + - made sure that in all histogram operations the errors are + properly propagated + - let ReInit determine fMcRadius from MMcRunHeader + + * mhflux/MHEnergyEst.cc, mhflux/MHThreshold.cc + - fixed the order in the constructor such that the Sumw2() does + correctly initialize the weights array + + * mhflux/MMcSpectrumWeight.cc: + - a minor code reordering + + * mjobs/MJSpectrum.cc: + - made sure that the histogram with the corsika spectrum has + the errors initialized and thus takes the weights correctly + into account + - corresponding to this changed some draw option to get the + same plots (hist) as before + - added a lot of comments to the code + - when the zenith angle weights are applied to the MC distribution + make sure that also the errors are correctly treated. + + + 2006/10/10 Thomas Bretz Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 8046) +++ trunk/MagicSoft/Mars/NEWS (revision 8047) @@ -3,4 +3,8 @@ *** Version + - sponde: In the calculation of the collection area(s) and the + distribution for MOnte Carlo and estimated energy the error + calculation was wrong because root didn't take the errors + properly into account... fixed. Index: trunk/MagicSoft/Mars/mhbase/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.cc (revision 8046) +++ trunk/MagicSoft/Mars/mhbase/MH.cc (revision 8047) @@ -1316,4 +1316,29 @@ // -------------------------------------------------------------------------- // +// This is a workaround for rrot-version <5.13/04 to get correct +// binomial errors even if weights have been used, do: +// h->Divide(h1, h2, 5, 2, "b"); +// MH::SetBinomialErrors(*h, *h1, *h2, 5, 2); +// +// see http://root.cern.ch/phpBB2/viewtopic.php?p=14818 +// +void MH::SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1, Double_t c2) +{ + for (Int_t binx=0; binx<=hres.GetNbinsX()+1; binx++) + { + const Double_t b1 = h1.GetBinContent(binx); + const Double_t b2 = h2.GetBinContent(binx); + const Double_t w = c2*b2 ? (c1*b1)/(c2*b2) : 0; + const Double_t e1 = h2.GetBinError(binx); + const Double_t e2 = h1.GetBinError(binx); + + const Double_t rc = ((1-2*w)*e1*e1+w*w*e2*e2)/(b2*b2); + + hres.SetBinError(binx, TMath::Sqrt(TMath::Abs(rc))); + } +} + +// -------------------------------------------------------------------------- +// // See MTask::PrintSkipped // Index: trunk/MagicSoft/Mars/mhbase/MH.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.h (revision 8046) +++ trunk/MagicSoft/Mars/mhbase/MH.h (revision 8047) @@ -76,4 +76,6 @@ static void SetBinning(TH1 *h, const TH1 *x); + static void SetBinomialErrors(TH1 &hres, const TH1 &h1, const TH1 &h2, Double_t c1=1, Double_t c2=1); + static void RemoveFirstBin(TH1 &h); Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 8046) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 8047) @@ -64,5 +64,5 @@ // MHCollectionArea::MHCollectionArea(const char *name, const char *title) - : fMcEvt(0), /*fEnergy(0),*/ fMcAreaRadius(300.), fIsExtern(kFALSE) + : fMcEvt(0), fMcAreaRadius(-1), fIsExtern(kFALSE) { // initialize the histogram for the distribution r vs E @@ -82,4 +82,5 @@ fHistSel.SetDirectory(NULL); fHistSel.UseCurrentStyle(); + fHistSel.SetLineColor(kBlue); fHistAll.SetName("AllEvts"); @@ -105,4 +106,10 @@ MH::SetBinning(&fHistSel, &binst, &binse); MH::SetBinning(&fHistAll, &binst, &binse); + + // For some unknown reasons this must be called after + // the binning has been initialized at least once + fHistSel.Sumw2(); + fHistAll.Sumw2(); + fHEnergy.Sumw2(); } @@ -114,45 +121,22 @@ void MHCollectionArea::CalcEfficiency() { - TH1D *hsel = fHistSel.ProjectionY(); - TH1D *hall = fHistAll.ProjectionY(); - - const Int_t nbinx = hsel->GetNbinsX(); - - // ----------------------------------------------------------- - // - // Impact parameter range: TO BE FIXED! Impact parameter shoud be - // read from run header, but it is not yet in!! + TH1D *hsel = fHistSel.ProjectionY("Spy", -1, 9999, "E");; + TH1D *hall = fHistAll.ProjectionY("Apy", -1, 9999, "E"); + + // + // Impact parameter range. // const Float_t totalarea = GetCollectionAreaAbs();//TMath::Pi() * (r2*r2 - r1*r1); - for (Int_t ix=1; ix<=nbinx; ix++) - { - const Float_t Na = hall->GetBinContent(ix); - - if (Na <= 0) - continue; - - const Float_t Ns = hsel->GetBinContent(ix); - - // Since Na is an estimate of the total number of showers generated - // in the energy bin, it may happen that Ns (triggered showers) is - // larger than Na. In that case, the bin is skipped: - - if (Na < Ns) - continue; - - const Double_t eff = Ns/Na; - const Double_t efferr = TMath::Sqrt((1.-eff)*Ns)/Na; - - const Float_t colarea = eff * totalarea; - const Float_t colareaerror = efferr * totalarea; - - fHEnergy.SetBinContent(ix, colarea); - fHEnergy.SetBinError(ix, colareaerror); - } + // "b" option: calculate binomial errors + fHEnergy.Divide(hsel, hall, totalarea, 1, "b"); +#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04) + MH::SetBinomialErrors(fHEnergy, *hsel, *hall, totalarea, 1); +#endif delete hsel; delete hall; } + Bool_t MHCollectionArea::SetupFill(const MParList *pl) @@ -196,4 +180,6 @@ MH::SetBinning(&fHistAll, &binst, &binse); + fMcAreaRadius = -1; + return kTRUE; } @@ -201,25 +187,24 @@ Bool_t MHCollectionArea::ReInit(MParList *plist) { + MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader"); + if (!runheader) + { + *fLog << err << "MMcRunHeader not found... abort." << endl; + return kFALSE; + } + + // FIXME: Does this need some weighting with the number of produced events? + if (runheader->GetImpactMax()>fMcAreaRadius*100) + { + fMcAreaRadius = 0.01*runheader->GetImpactMax(); // cm->m + *fLog << inf << "Maximum simulated impact: " << fMcAreaRadius << "m" << endl; + } + if (fIsExtern) return kTRUE; - MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader"); - if (!runheader) - { - *fLog << err << "MMcRunHeader not found... abort." << endl; - return kFALSE; - } - fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers(); *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl; - /* - if (fTheta>=0 && fTheta!=runheader->GetTelesTheta()) - { - *fLog << warn << "Warning - Read files have different TelesTheta ("; - *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl; - } - fTheta = runheader->GetTelesTheta(); - */ if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion()) *fLog << warn << "Warning - Read files have different Corsika versions..." << endl; @@ -352,6 +337,9 @@ if (h1 && h2 && h) { - h->Divide(h2, h1); - h->SetMinimum(0); + h->Divide(h2, h1, 1, 1, "b"); +#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04) + MH::SetBinomialErrors(*h, *h2, *h1); +#endif + h->SetMinimum(0); } @@ -433,5 +421,8 @@ gPad->SetLogx(); h = h2->DrawCopy(); - h->Divide(h1); + h->Divide(h2, h1, 1, 1, "b"); +#if ROOT_VERSION_CODE < ROOT_VERSION(5,13,04) + MH::SetBinomialErrors(*h, *h2, *h1); +#endif h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency"); h->SetDirectory(NULL); @@ -456,9 +447,7 @@ Bool_t MHCollectionArea::Fill(const MParContainer *par, const Stat_t weight) { - const Double_t energy = fMcEvt->GetEnergy()/*fEnergy->GetVal()*/; + const Double_t energy = fMcEvt->GetEnergy(); const Double_t theta = fMcEvt->GetTelescopeTheta()*TMath::RadToDeg(); - //*fLog << energy << " " << theta << endl; - fHistSel.Fill(theta, energy, weight); Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 8046) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 8047) @@ -60,19 +60,4 @@ Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; } -/* - void DrawAll(Option_t *option=""); - void DrawSel(Option_t *option=""); - - const TH1D *GetHist() { return fHistCol; } - const TH1D *GetHist() const { return fHistCol; } - - - - void CalcEfficiency(); - void CalcEfficiency2(); - - void Calc(const MHMcEnergyImpact &mcsel, const MHMcEnergyImpact &mcall); - void Calc(const MHMcEfficiency &heff); - */ void SetMcAreaRadius(Float_t x) { fMcAreaRadius = x; } Index: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 8046) +++ trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 8047) @@ -94,8 +94,4 @@ fHImpact.SetYTitle("E_{est}/E_{mc}-1"); - fHEnergy.Sumw2(); - fHImpact.Sumw2(); - fHResolution.Sumw2(); - MBinning binsi, binse, binst, binsr; binse.SetEdgesLog(25, 10, 1000000); @@ -107,4 +103,10 @@ SetBinning(&fHImpact, &binsi, &binsr); SetBinning(&fHResolution, &binse, &binse, &binsr); + + // For some unknown reasons this must be called after + // the binning has been initialized at least once + fHEnergy.Sumw2(); + fHResolution.Sumw2(); + fHImpact.Sumw2(); } @@ -243,4 +245,5 @@ { // Project into EnergyEst_ey + // the "e" ensures that errors are calculated TH1D *h1 = (TH1D*)fHEnergy.Project3D("rtlprmft_ex"); // E_Est TH1D *h2 = (TH1D*)fHEnergy.Project3D("rtlprmft_ey"); // E_mc Index: trunk/MagicSoft/Mars/mhflux/MHThreshold.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHThreshold.cc (revision 8046) +++ trunk/MagicSoft/Mars/mhflux/MHThreshold.cc (revision 8047) @@ -79,8 +79,9 @@ fHEnergy.SetDirectory(NULL); fHEnergy.UseCurrentStyle(); - fHEnergy.Sumw2(); MBinning binse(50, 20, 50000, "", "log"); binse.Apply(fHEnergy); + + fHEnergy.Sumw2(); } Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 8046) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 8047) @@ -443,6 +443,4 @@ Int_t MMcSpectrumWeight::Process() { - const Double_t e = fMcEvt->GetEnergy(); - Double_t w = 1; @@ -459,4 +457,5 @@ } + const Double_t e = fMcEvt->GetEnergy(); fWeight->SetVal(fFunc->Eval(e)*w); Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8046) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 8047) @@ -176,4 +176,5 @@ fLog->Separator("Initialize energy weighting"); + // Read Resources if (!CheckEnv(w)) { @@ -182,4 +183,5 @@ } + // Read the first corsika RunHeader from the MC files TChain chain("RunHeaders"); set.AddFilesOn(chain); @@ -195,4 +197,5 @@ } + // Propagate the run header to MMcSpectrumWeight if (!w.Set(*h)) { @@ -201,4 +204,5 @@ } + // Print new setup of MMcSpectrumWeight w.Print(); return kTRUE; @@ -266,6 +270,9 @@ fLog->Separator("Compiling original MC distribution"); + // The name of the input container is MMcEvtBasic weight.SetNameMcEvt("MMcEvtBasic"); + // Get the corresponding rule for the weights const TString w(weight.GetFormulaWeights()); + // Set back to MMcEvt weight.SetNameMcEvt(); @@ -282,4 +289,5 @@ // Prepare histogram h.Reset(); + h.Sumw2(); // Fill histogram from chain @@ -331,5 +339,5 @@ gPad->SetBorderMode(0); temp2.SetName("NVsTheta"); - temp2.DrawCopy(); + temp2.DrawCopy("hist"); c.cd(4); @@ -349,5 +357,5 @@ temp1.SetYTitle("Probability"); if (fDisplay) - temp1.DrawCopy(); + temp1.DrawCopy("hist"); return kTRUE; @@ -580,5 +588,5 @@ case 0: { - const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1); + const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}", f1, f1); const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0); @@ -638,18 +646,26 @@ TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const { + // Create copies of the histograms TH1D collarea(area.GetHEnergy()); TH1D spectrum(excess); TH1D weights; + + // Get spill-over corrections from energy estimation hest.GetWeights(weights); + // Print effective on-time cout << "Effective On time: " << ontime << "s" << endl; + // Setup spectrum plot + spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)"); + spectrum.SetYTitle("N/sm^{2}"); spectrum.SetDirectory(NULL); spectrum.SetBit(kCanDelete); + + // Divide by collection are and on-time spectrum.Scale(1./ontime); spectrum.Divide(&collarea); - spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)"); - spectrum.SetYTitle("N/sm^{2}"); - + + // Draw spectrum before applying spill-over corrections TCanvas &c1 = fDisplay->AddTab("Spectrum"); c1.Divide(2,2); @@ -679,4 +695,7 @@ weights.DrawCopy(); + // Apply spill-over correction (done't take the errors into account) + // They are supposed to be identical with the errors of the + // collection area and cancel out. //spectrum.Multiply(&weights); spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)"); @@ -840,5 +859,5 @@ //h2->Scale(1./ontime); //h2->Integral()); - h3->Scale(scale); //h3->Integral()); + h3->Scale(scale); //h3->Integral()); h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum())); @@ -1001,21 +1020,30 @@ plist.AddToList(&bins2); + // Initialize weighting to a new spectru + // m as defined in the resource file MMcSpectrumWeight weight; if (!InitWeighting(set, weight)) return kFALSE; + // Print the setup of the MAlphaFitter PrintSetup(fit); bins3.SetEdges(temp1, 'x'); + // Read the MC distribution as produced by corsika into + // temp2 (1D) and apply the weights previously determined TH1D temp2(temp1); if (!ReadOrigMCDistribution(set, temp2, weight)) return kFALSE; + // Calculate the weights according to the zenith angle distribution + // of the raw-data which have to be applied to the MC events if (!GetThetaDistribution(temp1, temp2)) return kFALSE; + // Tell the weighting task about the ZA-weights if (!fNoThetaWeights) weight.SetWeightsZd(&temp1); + // Refill excess histogram to determin the excess events TH1D excess; if (!Refill(plist, excess)) @@ -1026,4 +1054,7 @@ if (fSimpleMode) { + // Now we read the MC distribution as produced by corsika again + // with the spectral weights applied into a histogram vs. + // zenith angle and energy hist.UseCurrentStyle(); MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/); @@ -1031,14 +1062,21 @@ return kFALSE; + // No we apply the the zenith-angle-weights to the corsika produced + // MC distribution. Unfortunately this must be done manually + // because we are multiplying column by column if (!fRawMc) { for (int y=0; y