Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 8678) +++ trunk/MagicSoft/Mars/Changelog (revision 8679) @@ -25,4 +25,30 @@ the ListOfCleanups instead of removing the newly created canvas from the global list of canvases + + * mimage/MHHillas.cc, mimage/MHHillasExt.cc, mimage/MHHillasSrc.cc, + mimage/MHImagePar.cc, mimage/MHNewImagePar.cc: + - enables a lot of grids + + * mhbase/MBinning.[h,cc]: + - addec constructor with a TArrayD + + * mhbase/MH.[h,cc]: + - added two new member functions GetRange and GetRangeUser + + * mjoptim/MJOptimizeDisp.[h,cc]: + - added new option UseThetaSq + + * mjoptim/MJOptimizeEnergy.[h,cc]: + - added possibility to use also logarithmic fits + + * mmc/MMcCorsikaRunHeader.cc: + - improved Print function + + * msignal/MExtractTimeAndChargeSpline.cc: + - initialize rise and fall time as suggested by valgrind + + * mtools/MChisqEval.[h,cc]: + - implemented the possibility to use weights + Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 8678) +++ trunk/MagicSoft/Mars/NEWS (revision 8679) @@ -132,4 +132,17 @@ + PostCleanType: 3 + - star/ganymed: The old spark-cuts have been replaced by new ones. + These new spark cuts have been cross checked with the sequences + 84720 (before splitter), 101041 (after splitter) and 223539 (new + FADCs). In all cases they seperate pretty well, but of course not + perfect. Monte Carlos have suggested to tighten the cuts a little + bit further at low sizes. This has been done. + + - ganymed: IMPORTANT - The spark cuts have been changed to match + all available data more or less well. This might mean that for + YOUR data there are still sparks visible. + Whenever you do an analysis you have to make sure that NO sparks + survive your cuts! + - ganymed: The old hadronness and size cuts have gotten new indices (10, 11). The old indices 8 and 9 now mean an additional condition @@ -188,17 +201,4 @@ paremtrization set parameter 8 to 0. - - star/ganymed: The old spark-cuts have been replaced by new ones. - These new spark cuts have been cross checked with the sequences - 84720 (before splitter), 101041 (after splitter) and 223539 (new - FADCs). In all cases they seperate pretty well, but of course not - perfect. Monte Carlos have suggested to tighten the cuts a little - bit further at low sizes. This has been done. - - - ganymed: IMPORTANT - The spark cuts have been changed to match - all available data more or less well. This might mean that for - YOUR data there are still sparks visible. - Whenever you do an analysis you have to make sure that NO sparks - survive your cuts! - - ganymed: The distribution of observation time versus zenith angle is now shown for off- AND on-data. @@ -209,5 +209,4 @@ all off-regions from this file. Please do not forget to set the scaling in optimwobble.C properly: - fit.SetScaleMode(MAlphaFitter::kNone); fit.SetScaleUser(1./3); it is not yet read in automatically (will follow soon). @@ -244,9 +243,9 @@ the information is lost. - - sponde: added a new tab "EventDist" showing the unweighted event - distribution of your sample after cuts. This tells you how many - events with this energy you had in your MC files. The same information - you get from the error bars of the weighted histograms, but this is - less intuitive. + - sponde: added a new tab "EventDist" showing the unweighted real + absolute number of events of your sample after cuts. This tells + you how many events with this energy you had in your MC files. + The same information you get from the error bars of the weighted + histograms, but this is less intuitive. Index: trunk/MagicSoft/Mars/mhbase/MBinning.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MBinning.cc (revision 8678) +++ trunk/MagicSoft/Mars/mhbase/MBinning.cc (revision 8679) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MBinning.cc,v 1.17 2006-10-17 17:16:00 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MBinning.cc,v 1.18 2007-08-19 21:35:00 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -134,4 +134,16 @@ // -------------------------------------------------------------------------- // +// Initialize Binning from TArrayD. +// +MBinning::MBinning(const TArrayD &axis, const char *name, const char *title) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + SetEdges(axis); +} + +// -------------------------------------------------------------------------- +// // Search in the parameter list for the binning with name "name". If found, // set the edges and title accordingly. Default is name of object. Index: trunk/MagicSoft/Mars/mhbase/MBinning.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 8678) +++ trunk/MagicSoft/Mars/mhbase/MBinning.h (revision 8679) @@ -43,4 +43,5 @@ MBinning(const TH1 &h, const Char_t axis='x', const char *name=0, const char *title=0); MBinning(const TAxis &a, const char *name=0, const char *title=0); + MBinning(const TArrayD &a, const char *name=0, const char *title=0); void Copy(TObject &named) const Index: trunk/MagicSoft/Mars/mhbase/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.cc (revision 8678) +++ trunk/MagicSoft/Mars/mhbase/MH.cc (revision 8679) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.33 2007-08-06 14:44:13 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MH.cc,v 1.34 2007-08-19 21:36:00 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -20,5 +20,5 @@ ! Author(s): Thomas Bretz 07/2001 ! -! Copyright: MAGIC Software Development, 2000-2002 +! Copyright: MAGIC Software Development, 2000-2007 ! ! @@ -1409,4 +1409,37 @@ // -------------------------------------------------------------------------- // +// Return the first and last bin of the histogram which is not 0 +// +void MH::GetRange(const TH1 &h, Int_t &lo, Int_t &hi) +{ + lo = 0; + hi = 1; + + for (int i=2; i<=h.GetNbinsX(); i++) + { + if (lo==0 && h.GetBinContent(i)>0) + lo = i; + + if (h.GetBinContent(i)>0) + hi = i; + } +} + +// -------------------------------------------------------------------------- +// +// Return the lower edge of the first and the upper edge of the last bin +// of the histogram which is not 0 +// +void MH::GetRangeUser(const TH1 &h, Axis_t &lo, Axis_t &hi) +{ + Int_t f, l; + GetRange(h, f, l); + + lo = h.GetBinLowEdge(f); + hi = h.GetBinLowEdge(l+1); +} + +// -------------------------------------------------------------------------- +// // See MTask::PrintSkipped // Index: trunk/MagicSoft/Mars/mhbase/MH.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.h (revision 8678) +++ trunk/MagicSoft/Mars/mhbase/MH.h (revision 8679) @@ -119,4 +119,7 @@ static void ProjectionY(TH1D &dest, const TH2 &src, Int_t firstxbin=-1, Int_t lastxbin=9999); + static void GetRangeUser(const TH1 &hist, Axis_t &lo, Axis_t &hi); + static void GetRange(const TH1 &hist, Int_t &lo, Int_t &hi); + static TObject *FindObjectInPad(const char *name, TVirtualPad *pad=NULL); Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc (revision 8678) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc (revision 8679) @@ -178,5 +178,7 @@ MChisqEval eval; - eval.SetY1("sqrt(ThetaSquared.fVal)"); + eval.SetY1(fUseThetaSq?"ThetaSquared.fVal":"sqrt(ThetaSquared.fVal)"); + if (weights) + eval.SetNameWeight(); MMatrixLoop loop(&m); Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.h (revision 8678) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.h (revision 8679) @@ -10,6 +10,11 @@ class MJOptimizeDisp : public MJOptimize { +private: + Bool_t fUseThetaSq; + public: - MJOptimizeDisp() : MJOptimize() { } + MJOptimizeDisp() : MJOptimize(), fUseThetaSq(kFALSE) { } + + void EnableThetaSq(Bool_t b=kTRUE) { fUseThetaSq=b; } // Special optimizing routines Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.cc =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.cc (revision 8678) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.cc (revision 8679) @@ -66,4 +66,5 @@ // histograms #include "../mhflux/MHEnergyEst.h" +#include "../mtools/MChisqEval.h" // tasks @@ -108,4 +109,5 @@ MHMatrix m("M"); AddRulesToMatrix(m); + const Int_t map = m.AddColumn("MMcEvt.fEnergy"); parlist.AddToList(&m); @@ -136,4 +138,9 @@ fill.SetWeight(); + MChisqEval eval; + eval.SetY1(fOptimLog?Form("log10(MEnergyEst.fVal/M[%d])", map):Form("MEnergyEst.fVal-M[%d]", map)); + if (weights) + eval.SetNameWeight(); + MMatrixLoop loop(&m); @@ -143,4 +150,5 @@ tasklist.AddToList(weights); tasklist.AddToList(&fill); + tasklist.AddToList(&eval); // Optimize with the tasklist in this parameterlist Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.h (revision 8678) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.h (revision 8679) @@ -10,6 +10,11 @@ class MJOptimizeEnergy : public MJOptimize { +private: + Bool_t fOptimLog; + public: - MJOptimizeEnergy() : MJOptimize() { } + MJOptimizeEnergy() : MJOptimize(), fOptimLog(kFALSE) { } + + void EnableOptimLog(Bool_t b=kTRUE) { fOptimLog=b; } // Special optimizing routines Index: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 8678) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (revision 8679) @@ -1,4 +1,4 @@ /* ======================================================================== *\ -! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.68 2007-05-09 12:15:53 tbretz Exp $ +! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.69 2007-08-19 21:40:03 tbretz Exp $ ! -------------------------------------------------------------------------- ! @@ -179,5 +179,5 @@ // MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title) - : fHeightTm(0.5), fExtractionType(MExtralgoSpline::kIntegralRel) + : fRiseTimeHiGain(0), fFallTimeHiGain(0), fHeightTm(0.5), fExtractionType(MExtralgoSpline::kIntegralRel) { Index: trunk/MagicSoft/Mars/mtools/MChisqEval.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MChisqEval.cc (revision 8678) +++ trunk/MagicSoft/Mars/mtools/MChisqEval.cc (revision 8679) @@ -49,5 +49,6 @@ const TString MChisqEval::gsDefTitle = "Evaluate a chisq"; -MChisqEval::MChisqEval(const char *name, const char *title) : fData0(NULL), fData1(NULL), fNameResult("MinimizationValue") +MChisqEval::MChisqEval(const char *name, const char *title) + : fData0(NULL), fData1(NULL), fWeight(NULL), fNameResult("MinimizationValue") { fName = name ? name : gsDefName.Data(); @@ -55,5 +56,6 @@ } -MChisqEval::MChisqEval(MData *y1, const char *name, const char *title) : fData0(NULL), fData1(NULL), fNameResult("MinimizationValue") +MChisqEval::MChisqEval(MData *y1, const char *name, const char *title) +: fData0(NULL), fData1(NULL), fWeight(NULL), fNameResult("MinimizationValue") { fName = name ? name : gsDefName.Data(); @@ -62,5 +64,6 @@ } -MChisqEval::MChisqEval(MData *y1, MData *y2, const char *name, const char *title) : fData0(NULL), fData1(NULL), fNameResult("MinimizationValue") +MChisqEval::MChisqEval(MData *y1, MData *y2, const char *name, const char *title) +: fData0(NULL), fData1(NULL), fWeight(NULL), fNameResult("MinimizationValue") { fName = name ? name : gsDefName.Data(); @@ -112,4 +115,5 @@ { fChisq = 0; + fSumW = 0; if (!fData0) @@ -122,4 +126,11 @@ if (!fData1->PreProcess(plist)) return kFALSE; + + if (!fNameWeight.IsNull()) + { + fWeight = (MParameterD*)plist->FindObject(fNameWeight, "MParameterD"); + if (!fWeight) + return kFALSE; + } fResult = (MParameterD*)plist->FindCreateObj("MParameterD", fNameResult); @@ -138,5 +149,9 @@ const Double_t err = fData1 ? y1*y2 : 1; - fChisq += dy*dy/err; + const Double_t w = fWeight ? fWeight->GetVal() : 1; + + fChisq += w*dy*dy/err; + fSumW += w; + return kTRUE; } @@ -145,5 +160,5 @@ { if (GetNumExecutions()>0) - fChisq /= GetNumExecutions(); + fChisq /= fSumW; fResult->SetVal(fChisq); Index: trunk/MagicSoft/Mars/mtools/MChisqEval.h =================================================================== --- trunk/MagicSoft/Mars/mtools/MChisqEval.h (revision 8678) +++ trunk/MagicSoft/Mars/mtools/MChisqEval.h (revision 8679) @@ -15,11 +15,15 @@ static const TString gsDefTitle; - Double_t fChisq; //! Evaluated chi square - MParameterD *fResult; //! Storage for result - MData *fData0; // Data Member one (monte carlo data or chisq function) MData *fData1; // Data Member two (measured data) + MParameterD *fWeight; //! Storage for weight + MParameterD *fResult; //! Storage for result + TString fNameResult; + TString fNameWeight; + + Double_t fChisq; //! Evaluated chi square + Double_t fSumW; //! Sum of weights void StreamPrimitive(ostream &out) const; @@ -43,4 +47,5 @@ void SetOwner(Bool_t o=kTRUE) { o ? SetBit(kIsOwner) : ResetBit(kIsOwner); } + void SetNameWeight(TString w="MWeight") { fNameWeight=w; } Double_t GetChisq() const { return fChisq; } Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.cc =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.cc (revision 8678) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.cc (revision 8679) @@ -206,6 +206,7 @@ // Print Information about the Geometry of the camera // - *fLog << all << GetTitle() <<". Number of Telescopes: " << fNumCT << endl; - - fTelescopes.Print(); + *fLog << all << GetTitle() <<":" << endl; + *fLog << " Spectral Slope: " << fSlopeSpec << " from " << fELowLim << "GeV to " << fEUppLim << "GeV" << endl; + *fLog << " Number of Telescopes: " << fNumCT << endl; + fTelescopes.Print(); }