Index: trunk/MagicSoft/Cosy/tpoint/plot.C =================================================================== --- trunk/MagicSoft/Cosy/tpoint/plot.C (revision 7146) +++ trunk/MagicSoft/Cosy/tpoint/plot.C (revision 7147) @@ -121,5 +121,5 @@ }; -const Int_t counts = 12-4; +const Int_t counts = 13-4; Description_t desc[counts] = { @@ -139,5 +139,7 @@ {"05051", "TPoints Residuals 5/2004-1" , "tpoint/tpoint0505-1.txt"}, // Mirror alignment has been fixed - {"05052", "TPoints Residuals 5/2004-2" , "tpoint/tpoint0505-2.txt"} + {"05052", "TPoints Residuals 5/2004-2" , "tpoint/tpoint0505-2.txt"}, + // Fixes to pointing model due to fixing a screw + {"0506", "TPoints Residuals 6/2004" , "tpoint/tpoint0506.txt"} }; Index: trunk/MagicSoft/Cosy/tpoint/tpoint0506.txt =================================================================== --- trunk/MagicSoft/Cosy/tpoint/tpoint0506.txt (revision 7147) +++ trunk/MagicSoft/Cosy/tpoint/tpoint0506.txt (revision 7147) @@ -0,0 +1,28 @@ +# ------------- After fixing the Shaftencoder screw ------------------ +# ------------- Done 7.6.2005 ------------------ + +# --- 08.06.2005 --- 22:09:17.757 +# Her Zeta 40, 2.81 mag, Zd 45, Az 75 +75.33279 56.64848 217.8299 47.64721 16.68806 31.60306 0.07855271 -0.0245924 53529.923122 203.8 3.292 +# Rasteban, Zd 43 deg, Az 43 deg, 2.79 mag. +42.3295 47.82316 184.8677 38.85592 17.50722 52.30139 0.1119815 -0.04028769 53529.932071 189.2 3.191 +# Cyg Delta-18, Zd 63, Az 50, 2.87 mag, +51.08745 26.89075 193.625 17.93327 19.74944 45.13083 0.1227235 -0.01832511 53529.936881 205.9 3.397 +# Eltanin, Zd 44 deg, Az 44 deg, 2.23 mag. +44.12169 46.2001 186.6513 37.22763 17.94333 51.48889 0.1068942 -0.0286855 53529.941836 198.9 2.802 +# Sulafat, Zd 55, Az 68, 3.24 mag, +68.52952 35.33758 211.0607 26.37314 18.9825 32.68944 0.1130563 -0.02093151 53529.947116 216 3.649 +# Algenubi, Zd 68, Az 285, 2.98 mag, +-74.31378 22.35625 68.19518 13.51402 9.764167 23.77417 0.1017909 -0.01241812 53529.954137 211.3 3.421 +# Uma Theta 25, Zd 64, Az -43, 3.17 mag, +-43.33788 26.1297 99.1489 17.27408 9.5475 51.67722 0.09226138 0.01776548 53529.96392 197.4 3.951 +# Uma Psi 52, Zd 50, Az -55, 3.01 mag, +-54.35609 39.60993 88.152 30.77825 11.16111 44.49861 0.1141847 -0.018317 53529.969753 218.8 3.759 +# Cor Carioli, Zd 31, Az -63, 2.9 mag, +-62.88338 58.14462 79.59492 49.27889 12.93389 38.31833 0.08449887 -0.02779159 53529.975535 217.2 3.269 +# Zosma, Zd 59, Az 278, 2.56 mag, +-82.27258 31.09557 60.22991 22.25793 11.235 20.52361 0.1089809 -0.01533695 53529.982977 203.7 3.243 + +# ------------- New IE value in Pointing Model ------------------ +# ------------- Since 9./10.6.2005 (incl)------------------ + Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7146) +++ trunk/MagicSoft/Mars/Changelog (revision 7147) @@ -21,4 +21,19 @@ -*-*- END OF LINE -*-*- + 2005/06/13 Thomas Bretz + + * mcalib/MCalibrationHiLoCam.h: + - added Print to //*MENU* + + * mhflux/MHThetaSq.[h,cc]: + - added resources for fNumBinsSignal and fNumBinsTotal + + * mjobs/MJSpectrum.[h,cc]: + - implemented weighting in theta, so we get better statistics + - improved output + - added plotting other spectras + + + 2005/06/10 Daniela Dorner Index: trunk/MagicSoft/Mars/mcalib/MCalibrationHiLoCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationHiLoCam.h (revision 7146) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationHiLoCam.h (revision 7147) @@ -19,5 +19,5 @@ // Prints - void Print(Option_t *o="") const; + void Print(Option_t *o="") const; //*MENU* // Others Index: trunk/MagicSoft/Mars/mhflux/MHThetaSq.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHThetaSq.cc (revision 7146) +++ trunk/MagicSoft/Mars/mhflux/MHThetaSq.cc (revision 7147) @@ -34,4 +34,9 @@ // For more detailes see MHAlpha. // +// Version 2: +// --------- +// + UInt_t fNumBinsSignal; +// + UInt_t fNumBinsTotal; +// ////////////////////////////////////////////////////////////////////////////// #include "MHThetaSq.h" @@ -57,5 +62,5 @@ // MHThetaSq::MHThetaSq(const char *name, const char *title) - : MHAlpha(name, title), fThetaSq(0) + : MHAlpha(name, title), fThetaSq(0), fNumBinsSignal(3), fNumBinsTotal(45) { // @@ -120,6 +125,6 @@ // Calculate bining which fits alpha-cut const Double_t intmax = fit->GetSignalIntegralMax(); - const UInt_t nbins = 75; - const UInt_t nsig = 5; + const UInt_t nbins = fNumBinsTotal; + const UInt_t nsig = fNumBinsSignal; MBinning binsa(nbins, 0, nbins*intmax/nsig); @@ -201,2 +206,21 @@ // fMap[4] = fMatrix->AddColumn("MTime.GetAxisTime"); } + +Int_t MHThetaSq::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Int_t rc = MHAlpha::ReadEnv(env, prefix, print); + if (rc==kERROR) + return kERROR; + + if (IsEnvDefined(env, prefix, "NumBinsSignal", print)) + { + SetNumBinsSignal(GetEnvValue(env, prefix, "NumBinsSignal", (Int_t)fNumBinsSignal)); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "NumBinsTotal", print)) + { + SetNumBinsTotal(GetEnvValue(env, prefix, "NumBinsTotal", (Int_t)fNumBinsTotal)); + rc = kTRUE; + } + return rc; +} Index: trunk/MagicSoft/Mars/mhflux/MHThetaSq.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHThetaSq.h (revision 7146) +++ trunk/MagicSoft/Mars/mhflux/MHThetaSq.h (revision 7147) @@ -13,4 +13,7 @@ MParameterD *fThetaSq; //! + UInt_t fNumBinsSignal; + UInt_t fNumBinsTotal; + Bool_t GetParameter(const MParList &pl); Double_t GetVal() const; @@ -23,8 +26,14 @@ void InitMapping(MHMatrix *mat, Int_t type=0); + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + public: MHThetaSq(const char *name=NULL, const char *title=NULL); - ClassDef(MHThetaSq, 1) // Theta-Plot which is fitted online + void SetNumBinsSignal(UInt_t n) { fNumBinsSignal=TMath::Max(n, 1U); } + void SetNumBinsTotal(UInt_t n) { fNumBinsTotal =TMath::Max(n, 1U); } + + ClassDef(MHThetaSq, 2) // Theta-Plot which is fitted online }; Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7146) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7147) @@ -38,4 +38,5 @@ #include #include +#include #include #include @@ -84,5 +85,6 @@ MJSpectrum::MJSpectrum(const char *name, const char *title) : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0), - fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE) + fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE), + fNoThetaWeights(kFALSE) { fName = name ? name : "MJSpectrum"; @@ -336,5 +338,5 @@ // Calculate the Probability temp1.Divide(&temp2); - temp1.Scale(1./temp1.GetMaximum()); + temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral()); // Some cosmetics: Name, Axis, etc. @@ -348,4 +350,8 @@ } +// -------------------------------------------------------------------------- +// +// Display the final theta distribution. +// void MJSpectrum::DisplayResult(const TH2D &h2) const { @@ -376,4 +382,11 @@ } +// -------------------------------------------------------------------------- +// +// Fills the excess histogram (vs E-est) from the events stored in the +// ganymed result file and therefor estimates the energy. +// +// The resulting histogram excess-vs-energy ist copied into h2. +// Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const { @@ -531,4 +544,13 @@ } +// -------------------------------------------------------------------------- +// +// Calculate the final spectrum from: +// - collection area +// - excess +// - correction coefficients +// - ontime +// and display it +// TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const { @@ -573,10 +595,13 @@ weights.DrawCopy(); - //spectrum.Divide(&weights); - spectrum.Multiply(&weights); - spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)"); + //spectrum.Multiply(&weights); + spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)"); + spectrum.SetBit(TH1::kNoStats); for (int i=0; iSetGridx(); gPad->SetGridy(); + spectrum.SetMinimum(1e-12); spectrum.SetXTitle("E [GeV]"); spectrum.SetYTitle("N/TeVsm^{2}"); spectrum.DrawCopy(); - TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4); + TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4); f.SetParameter(0, -2.87); f.SetParameter(1, 1.9e-6); f.SetLineColor(kGreen); - spectrum.Fit(&f, "NIM", "", 55, 2e4); + spectrum.Fit(&f, "NIM", "", 100, 4000); f.DrawCopy("same"); - /* - TString str; - str += "(1.68#pm0.15)10^{-7}"; - str += "(\\frac{E}{TeV})^{-2.59#pm0.06}"; - str += "\\frac{ph}{TeVm^{2}s}"; - - TLatex tex; - tex.DrawLatex(2e2, 7e-5, str); - */ + const Double_t p0 = f.GetParameter(0); + const Double_t p1 = f.GetParameter(1); + + const Double_t e0 = f.GetParError(0); + const Double_t e1 = f.GetParError(1); + + const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1))); + const Double_t exp = TMath::Power(10, np); + + TString str; + str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np); + str += Form("(\\frac{E}{TeV})^{-%.2f#pm%.2f}", p0, e0); + str += "\\frac{ph}{TeVm^{2}s}"; + + TLatex tex; + tex.SetTextSize(0.045); + tex.SetBit(TLatex::kTextNDC); + tex.DrawLatex(0.45, 0.935, str); + + str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF()); + tex.DrawLatex(0.70, 0.83, str); TArrayD res(2); res[0] = f.GetParameter(0); res[1] = f.GetParameter(1); + + + // Plot other spectra from Whipple + f.SetParameter(0, -2.45); + f.SetParameter(1, 3.3e-7); + f.SetRange(300, 8000); + f.SetLineColor(kBlack); + f.SetLineStyle(kDashed); + f.DrawCopy("same"); + + // Plot other spectra from Cangaroo + f.SetParameter(0, -2.53); + f.SetParameter(1, 2.0e-7); + f.SetRange(7000, 50000); + f.SetLineColor(kBlack); + f.SetLineStyle(kDashed); + f.DrawCopy("same"); + + // Plot other spectra from Robert + f.SetParameter(0, -2.59); + f.SetParameter(1, 2.58e-7); + f.SetRange(150, 1500); + f.SetLineColor(kBlack); + f.SetLineStyle(kDashed); + f.DrawCopy("same"); + + // Plot other spectra from HEGRA + f.SetParameter(0, -2.61); + f.SetParameter(1, 2.7e-7); + f.SetRange(1000, 20000); + f.SetLineColor(kBlack); + f.SetLineStyle(kDashed); + f.DrawCopy("same"); + return res; } +// -------------------------------------------------------------------------- +// +// Scale some image parameter plots using the scale factor and plot them +// together with the corresponding MC histograms. +// Called from DisplaySize +// Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const { @@ -664,4 +742,9 @@ } +// -------------------------------------------------------------------------- +// +// Take a lot of histograms and plot them together in one plot. +// Calls PlotSame +// Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const { @@ -698,5 +781,5 @@ gPad->SetGridy(); - excess->SetTitle("Number of excess events vs Size (data, mc/blue)"); + excess->SetTitle("Excess events vs Size (data/black, mc/blue)"); excess = excess->DrawCopy("E2"); // Don't do this on the original object! @@ -727,5 +810,9 @@ histsel->KolmogorovTest(excess, "DX"); fLog->Separator("Chi^2 Test"); - histsel->Chi2Test(excess, "P"); + const Double_t p = histsel->Chi2Test(excess, "P"); + + TLatex tex; + tex.SetBit(TLatex::kTextNDC); + tex.DrawLatex(0.7, 0.93, Form("P(\\chi^{2})=%.0f", p*100)); } } @@ -814,4 +901,7 @@ if (!GetThetaDistribution(temp1, temp2)) return kFALSE; + + if (!fNoThetaWeights) + weight.SetZdWeights(&temp1); TH1D excess; @@ -930,5 +1020,5 @@ tlist2.AddToList(&read); - if (!fRawMc) + if (!fRawMc && fNoThetaWeights) tlist2.AddToList(&contsel); tlist2.AddToList(&calc);