Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4886) +++ /trunk/MagicSoft/Mars/Changelog (revision 4887) @@ -44,4 +44,24 @@ and maybe to secure everything from compiler bugs a bit more + * callisto.cc, star.cc: + - do not allow implicit batch-mode + + * showlog.cc: + - removed an unused variable + + * manalysis/MEventRateCalc.[h,cc]: + - started implementing a corrsponding time + + * manalysis/MParameters.[h,cc], manalysis/AnalysisLinkDef.h: + - added MParameterDerr + + * mbadpixels/MHBadPixels.[h,cc]: + - removed obsolete fPedPhotCam + + * mbase/MTime.[h,cc]: + - added AddMilliSeconds + - added SetMean + + @@ -68,4 +88,5 @@ - included functions AddCamContent(MArrayD...) and SetCamContent(MArrayD... ) + Index: /trunk/MagicSoft/Mars/callisto.cc =================================================================== --- /trunk/MagicSoft/Mars/callisto.cc (revision 4886) +++ /trunk/MagicSoft/Mars/callisto.cc (revision 4887) @@ -262,5 +262,5 @@ TApplication app("Callisto", &argc, argv); - if (!gROOT->IsBatch() && !gClient) + if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch) { gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl; Index: /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 4886) +++ /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 4887) @@ -40,4 +40,5 @@ #pragma link C++ class MParameterI+; #pragma link C++ class MParameterD+; +#pragma link C++ class MParameterDerr+; //#pragma link C++ class MParameters+; Index: /trunk/MagicSoft/Mars/manalysis/MEventRateCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEventRateCalc.cc (revision 4886) +++ /trunk/MagicSoft/Mars/manalysis/MEventRateCalc.cc (revision 4887) @@ -88,4 +88,5 @@ const TString MEventRateCalc::gsNameEventRate = "MEventRate"; const TString MEventRateCalc::gsNameTimeDiff = "MTimeDiff"; +const TString MEventRateCalc::gsNameTimeRate = "MTimeRate"; const Int_t MEventRateCalc::gsNumEvents = 1000; @@ -96,6 +97,7 @@ // MEventRateCalc::MEventRateCalc(const char *name, const char *title) - : fNameTime(gsNameTime), fNameEventRate(gsNameEventRate), - fNameTimeDiff(gsNameTimeDiff), fTimes(gsNumEvents) + : fNameEventRate(gsNameEventRate), fNameTime(gsNameTime), + fNameTimeRate(gsNameTimeRate), fNameTimeDiff(gsNameTimeDiff), + fTimes(gsNumEvents) { fName = name ? name : gsDefName.Data(); @@ -118,5 +120,5 @@ if (!fTime) { - *fLog << err << "MTime not found... aborting." << endl; + *fLog << err << AddSerialNumber(fNameTime) << " [MTime] not found... aborting." << endl; return kFALSE; } @@ -124,4 +126,8 @@ fRate = (MEventRate*)pList->FindCreateObj("MEventRate", AddSerialNumber(fNameEventRate)); if (!fRate) + return kFALSE; + + fTimeRate = (MTime*)pList->FindCreateObj("MTime", AddSerialNumber(fNameTimeRate)); + if (!fTimeRate) return kFALSE; @@ -145,20 +151,30 @@ const ULong_t exec = GetNumExecutions()-1; + // Calculate the rate const UInt_t n = fTimes.GetSize(); + // Get the positon in the ring-buffer const UInt_t n1 = exec; const UInt_t n2 = exec>=n ? exec+1 : 0; + // Store the current event time fTimes[n1%n] = *fTime; + // Get the number of events const UInt_t cnt = n1SetVal(exec==0 ? -1 : fTimes[n1%n] - fTimes[(n1+n-1)%n]); + fTimeDiff->SetReadyToSave(); + + // Store the rate fRate->SetRate(exec>1?rate:0, cnt); fRate->SetReadyToSave(); - fTimeDiff->SetVal(exec==0 ? -1 : fTimes[n1%n] - fTimes[(n1+n-1)%n]); - fTimeDiff->SetReadyToSave(); + // Store the corresponding time + fTimeRate->SetMean(fTimes[n1%n], fTimes[n2%n]); return kTRUE; Index: /trunk/MagicSoft/Mars/manalysis/MEventRateCalc.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEventRateCalc.h (revision 4886) +++ /trunk/MagicSoft/Mars/manalysis/MEventRateCalc.h (revision 4887) @@ -19,7 +19,8 @@ static const TString gsDefTitle; //! Default title of container + static const TString gsNameEventRate; //! default name of rate container static const TString gsNameTime; //! Default name of time container - static const TString gsNameEventRate; //! default name of rate container static const TString gsNameTimeDiff; //! default name of time-diff container + static const TString gsNameTimeRate; //! default name of time-rate container static const Int_t gsNumEvents; //! Default number of events @@ -27,13 +28,14 @@ MTime *fTime; //! pointer to event time + MTime *fTimeRate; //! pointer to time of event rate MEventRate *fRate; //! pointer to rate storage container MParameterD *fTimeDiff; //! Difference of time between two consecutive events + TString fNameEventRate; // name of event rate container TString fNameTime; // name of time container - TString fNameEventRate; // name of event rate container + TString fNameTimeRate; // name of event rate time container TString fNameTimeDiff; // name of time-diff container TArrayD fTimes; //! internal array to store the last n event times - Int_t PreProcess(MParList *pList); @@ -48,7 +50,8 @@ void SetNumEvents(ULong_t num) { fTimes.Set(num); } + void SetNameEventRate(const char *name) { fNameEventRate = name; } void SetNameTime(const char *name) { fNameTime = name; } - void SetNameEventRate(const char *name) { fNameEventRate = name; } void SetNameTimeDiff(const char *name) { fNameTimeDiff = name; } + void SetNameTimeRate(const char *name) { fNameTimeRate = name; } ClassDef(MEventRateCalc, 1)// Task to calculate event rates Index: /trunk/MagicSoft/Mars/manalysis/MParameters.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MParameters.cc (revision 4886) +++ /trunk/MagicSoft/Mars/manalysis/MParameters.cc (revision 4887) @@ -18,5 +18,5 @@ ! Author(s): Thomas Bretz 03/2003 ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -42,4 +42,5 @@ ClassImp(MParameterD); ClassImp(MParameterI); +ClassImp(MParameterDerr); //ClassImp(MParameter); @@ -60,4 +61,14 @@ // Default constructor. // +MParameterDerr::MParameterDerr(const char *name, const char *title) +{ + fName = name ? name : "MParameterDerr"; + fTitle = title ? title : "Storgare container for general parameters (double) and its error"; +} + +// -------------------------------------------------------------------------- +// +// Default constructor. +// MParameterI::MParameterI(const char *name, const char *title) { Index: /trunk/MagicSoft/Mars/manalysis/MParameters.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MParameters.h (revision 4886) +++ /trunk/MagicSoft/Mars/manalysis/MParameters.h (revision 4887) @@ -18,4 +18,20 @@ ClassDef(MParameterD, 1) // Container to hold a generalized parameters (double) +}; + +class MParameterDerr : public MParContainer +{ +private: + Double_t fVal; + Double_t fErr; + +public: + MParameterDerr(const char *name=NULL, const char *title=NULL); + + void SetVal(Double_t v, Double_t e) { fVal = v; fErr = e; } + Double_t GetVal() const { return fVal; } + Double_t GetErr() const { return fErr; } + + ClassDef(MParameterDerr, 1) // Container to hold a generalized parameters (double) and its Error }; Index: /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.cc =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.cc (revision 4886) +++ /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.cc (revision 4887) @@ -35,5 +35,4 @@ #include "MBadPixelsCam.h" #include "MGeomCam.h" -#include "MPedPhotCam.h" #include "MParList.h" #include "MBinning.h" @@ -51,4 +50,5 @@ // MHBadPixels::MHBadPixels(const char *name, const char *title) + : fNamePedPhotCam("MPedPhotCam") { fName = name ? name : "MHBadPixels"; @@ -67,5 +67,4 @@ fBadN.SetYTitle("number of bad pixels"); - fNamePedPhotCam = "MPedPhotCamFromData"; } @@ -79,9 +78,7 @@ if (!fCam) { - *fLog << err << "MHBadPixels::SetupFill; MGeomCam not found... aborting." << endl; + *fLog << err << "MGeomCam not found... aborting." << endl; return kFALSE; } - *fLog << "MHBadPixels::SetupFill; fCam = " << fCam << endl; - fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); if (!fPointPos) @@ -91,14 +88,4 @@ } - - fPedPhot = (MPedPhotCam*)plist->FindObject(AddSerialNumber(fNamePedPhotCam), "MPedPhotCam"); - if (!fPedPhot) - { - *fLog << err << AddSerialNumber(fNamePedPhotCam) - << "[MPedPhotCam] not found... aborting." << endl; - return kFALSE; - } - fPedPhot->InitSize(fCam->GetNumPixels()); - //---------------------------------------------------- // Get Theta Binning @@ -106,12 +93,10 @@ if (!binstheta) { - *fLog << err << "Object 'BinningTheta' [MBinning] not found... aborting" << endl; + *fLog << err << "BinningTheta [MBinning] not found... aborting" << endl; return kFALSE; } // Get binning for pixel number - const UInt_t npix1 = fPedPhot->GetSize()+1; - - *fLog << "MHBadPixels::SetupFill; npix1 = " << npix1 << endl; + const UInt_t npix1 = fCam->GetNumPixels()+1; MBinning binspix("BinningPixel"); @@ -121,8 +106,4 @@ SetBinning(&fBadId, binstheta, &binspix); SetBinning(&fBadN, binstheta, &binspix); - - //---------------------------------------------------- - *fLog << inf << "Name of MPedPhotCam container : " << fNamePedPhotCam - << endl; return kTRUE; @@ -141,19 +122,16 @@ Double_t theta = fPointPos->GetZd(); - const MBadPixelsCam *fBadPixels = (MBadPixelsCam*)par; + const MBadPixelsCam *badpixels = (MBadPixelsCam*)par; - const UInt_t entries = fPedPhot->GetSize(); + const UInt_t entries = badpixels->GetSize(); + UInt_t nb = 0; for (UInt_t i=0; iUpdate(); } -//========================================================================== - - - - - - - - - - - Index: /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.h =================================================================== --- /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.h (revision 4886) +++ /trunk/MagicSoft/Mars/mbadpixels/MHBadPixels.h (revision 4887) @@ -18,5 +18,4 @@ private: MGeomCam *fCam; //! - MPedPhotCam *fPedPhot; //! MPointingPos *fPointPos; //! @@ -37,7 +36,4 @@ const TH2D *GetBadN() const { return &fBadN; } - TH2 *GetBadIdByName(const TString name) { return &fBadId; } - TH2 *GetBadNByName(const TString name) { return &fBadN; } - void Draw(Option_t* option = ""); Bool_t SetupFill(const MParList *plist); Index: /trunk/MagicSoft/Mars/mbase/MTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTime.cc (revision 4886) +++ /trunk/MagicSoft/Mars/mbase/MTime.cc (revision 4887) @@ -563,2 +563,127 @@ return fin; } + +void MTime::AddMilliSeconds(UInt_t ms) +{ + fTime += ms; + + fTime += 11*kHour; + fMjd += (Long_t)fTime/kDay; + fTime = (Long_t)fTime%kDay; + fTime -= 11*kHour; +} +/* +MTime MTime::operator-(const MTime &tm1) +{ + const MTime &tm0 = *this; + + MTime t0 = tm0>tm1 ? tm0 : tm1; + const MTime &t1 = tm0>tm1 ? tm1 : tm0; + + if (t0.fNanoSec999999) + { + t0.fNanoSec -= 1000000; + t0.fTime += kDay; + } + + t0.fTime += t1.fTime; + + if ((Long_t)t0.fTime>=(Long_t)kHour*13) + { + t0.fTime -= kDay; + t0.fMjd++; + } + + t0.fMjd += t1.fMjd; + + return t0; +} + +void MTime::operator+=(const MTime &t) +{ + *this = *this+t; +} +*/ + +void MTime::SetMean(const MTime &t0, const MTime &t1) +{ + // This could be an operator+ + *this = t0; + + fNanoSec += t1.fNanoSec; + + if (fNanoSec>999999) + { + fNanoSec -= 1000000; + fTime += kDay; + } + + fTime += t1.fTime; + + if ((Long_t)fTime>=(Long_t)kHour*13) + { + fTime -= kDay; + fMjd++; + } + + fMjd += t1.fMjd; + + // This could be an operator/ + if ((Long_t)fTime<0) + { + fTime += kDay; + fMjd--; + } + + Int_t reminder = fMjd%2; + fMjd /= 2; + + fTime += reminder*kDay; + reminder = (Long_t)fTime%2; + fTime /= 2; + + fNanoSec += reminder*1000000; + fNanoSec /= 2; + + fTime += 11*kHour; + fMjd += (Long_t)fTime/kDay; + fTime = (Long_t)fTime%kDay; + fTime -= 11*kHour; +} + +void MTime::SetMean(Double_t t0, Double_t t1) +{ + const Double_t mean = (t0+t1)*(500./kDay); + SetMjd(mean); +} Index: /trunk/MagicSoft/Mars/mbase/MTime.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MTime.h (revision 4886) +++ /trunk/MagicSoft/Mars/mbase/MTime.h (revision 4887) @@ -35,4 +35,5 @@ UInt_t fNanoSec; // [ns] NanoSec part of TimeOfDay (<1000000) + ULong_t GetTime24() const { @@ -75,4 +76,5 @@ void Now(); + // Setter functions Bool_t Set(UShort_t y, Byte_t m, Byte_t d, Byte_t h=13, Byte_t min=0, Byte_t s=0, UShort_t ms=0, UInt_t ns=0); void Set(const struct timeval &tv); @@ -85,4 +87,6 @@ Bool_t SetMjd(UInt_t mjd, ULong_t ms, UInt_t ns=0); void SetMjd(Double_t m); + + // Getter functions Double_t GetMjd() const; Double_t GetGmst() const; @@ -111,60 +115,32 @@ UInt_t Sec() const { Byte_t h, m, s; GetTime(h,m,s); return s; } + // I/O functions istream &ReadBinary(istream &fin); - operator double() const //[s] - { - return ((Double_t)fMjd*kDay+(Long_t)fTime+fNanoSec/1e6)/1000; - } - double operator()() const //[s] - { - return operator double(); - } - - bool operator<(const MTime &t) const - { - if (fMjd(const MTime &t) const - { - if (fMjd>t.fMjd) - return true; - if (fMjd==t.fMjd && fTime>t.fTime) - return true; - if (fMjd==t.fMjd && fTime==t.fTime && fNanoSec>t.fNanoSec) - return true; - return false; - } - - bool operator<=(const MTime &t) const - { - return !operator>(t); - } - - bool operator>=(const MTime &t) const - { - return !operator<(t); - } - - bool operator==(const MTime &t) const - { - return fNanoSec==t.fNanoSec && fTime==t.fTime && fMjd==t.fMjd; - } - - bool operator!=(const MTime &t) const - { - return fNanoSec!=t.fNanoSec || fTime!=t.fTime || fMjd!=t.fMjd; - } - - bool operator!() const - { - return fNanoSec==0 && (ULong_t)fTime==0 && fMjd==0; - } + // Conversion functions + operator double() const; //[s] + double operator()() const; //[s] + + // Calculation functions + void AddMilliSeconds(UInt_t ms); + void SetMean(const MTime &t0, const MTime &t1); + void SetMean(Double_t t0, Double_t t1); + /* + MTime operator-(const MTime &tm1); + void operator-=(const MTime &t); + MTime operator+(const MTime &t1); + void operator+=(const MTime &t); + */ + + // Base comparison operators + bool operator<(const MTime &t) const; + bool operator>(const MTime &t) const; + + // Derived comparison operators + bool operator<=(const MTime &t) const; + bool operator>=(const MTime &t) const; + bool operator==(const MTime &t) const; + bool operator!=(const MTime &t) const; + bool operator!() const; ClassDef(MTime, 3) //A generalized MARS time stamp @@ -186,3 +162,60 @@ } -#endif +inline MTime::operator double() const //[s] +{ + return ((Double_t)fMjd*kDay+(Long_t)fTime+fNanoSec/1e6)/1000; +} + +inline double MTime::operator()() const //[s] +{ + return operator double(); +} + +inline bool MTime::operator<(const MTime &t) const +{ + if (fMjd(const MTime &t) const +{ + if (fMjd>t.fMjd) + return true; + if (fMjd==t.fMjd && fTime>t.fTime) + return true; + if (fMjd==t.fMjd && fTime==t.fTime && fNanoSec>t.fNanoSec) + return true; + return false; +} + +inline bool MTime::operator<=(const MTime &t) const +{ + return !operator>(t); +} + +inline bool MTime::operator>=(const MTime &t) const +{ + return !operator<(t); +} + +inline bool MTime::operator==(const MTime &t) const +{ + return fNanoSec==t.fNanoSec && fTime==t.fTime && fMjd==t.fMjd; +} + +inline bool MTime::operator!=(const MTime &t) const +{ + return fNanoSec!=t.fNanoSec || fTime!=t.fTime || fMjd!=t.fMjd; +} + +inline bool MTime::operator!() const +{ + return fNanoSec==0 && (ULong_t)fTime==0 && fMjd==0; +} + +#endif Index: /trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4887) +++ /trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4887) @@ -0,0 +1,451 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 8/2002 +! Author(s): Wolfgang Wittek, 1/2002 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHEffectiveOnTime +// +// Fills a 2-D histogram Time-difference vs. Theta +// +// From this histogram the effective on-time is determined by a fit. +// The result of the fit (see Fit()) and the fit-parameters (like chi^2) +// are stored in corresponding histograms +// +////////////////////////////////////////////////////////////////////////////// + +#include "MHEffectiveOnTime.h" + +#include +#include +#include +#include + +#include "MTime.h" +#include "MParameters.h" +#include "MPointingPos.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHEffectiveOnTime); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title) + : /*fLastTime(0), fFirstTime(-1),*/ fIsFinalized(kFALSE), fInterval(60) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHEffectiveOnTime"; + fTitle = title ? title : "2-D histogram in Theta and time difference"; + + // Main histogram + fHTimeDiff.SetXTitle("\\Delta t [s]"); + fHTimeDiff.SetYTitle("\\Theta [\\circ]"); + fHTimeDiff.UseCurrentStyle(); + fHTimeDiff.SetDirectory(NULL); + + // effective on time versus theta + fHEffOn.SetName("EffOnTime"); + fHEffOn.SetTitle("Effective On Time T_{eff}"); + fHEffOn.SetXTitle("\\Theta [\\circ]"); + fHEffOn.SetYTitle("T_{eff} [s]"); + fHEffOn.UseCurrentStyle(); + fHEffOn.SetDirectory(NULL); + fHEffOn.GetYaxis()->SetTitleOffset(1.2); + //fHEffOn.Sumw2(); + + // chi2/NDF versus theta + fHChi2.SetName("Chi2/NDF"); + fHChi2.SetTitle("\\chi^{2}/NDF of Effective On Time Fit"); + fHChi2.SetXTitle("\\Theta [\\circ]"); + fHChi2.SetYTitle("\\chi^{2}/NDF"); + fHChi2.UseCurrentStyle(); + fHChi2.SetDirectory(NULL); + fHChi2.GetYaxis()->SetTitleOffset(1.2); + //fHChi2.Sumw2(); + + // chi2 probability + fHProb.SetName("Prob"); + fHProb.SetTitle("\\chi^{2} Probability of Fit"); + fHProb.SetXTitle("\\Theta [\\circ]"); + fHProb.SetYTitle("p [%]"); + fHProb.UseCurrentStyle(); + fHProb.SetDirectory(NULL); + fHProb.GetYaxis()->SetTitleOffset(1.2); + fHProb.SetMaximum(101); + //fHChi2.Sumw2(); + + // lambda versus theta + fHLambda.SetName("lambda"); + fHLambda.SetTitle("lambda of Effective On Time Fit"); + fHLambda.SetXTitle("\\Theta [\\circ]"); + fHLambda.SetYTitle("\\lambda [Hz]"); + fHLambda.UseCurrentStyle(); + fHLambda.SetDirectory(NULL); + // fHLambda.Sumw2(); + + // N0 versus theta + fHN0.SetName("N0"); + fHN0.SetTitle("Ideal number of events N_{0}"); + fHN0.SetXTitle("\\Theta [\\circ]"); + fHN0.SetYTitle("N_{0}"); + fHN0.UseCurrentStyle(); + fHN0.SetDirectory(NULL); + //fHN0del.Sumw2(); + + // setup binning + MBinning btheta("BinningTheta"); + btheta.SetEdgesCos(51, 0, 60); + + MBinning btime("BinningTimeDiff"); + btime.SetEdges(50, 0, 0.1); + + MH::SetBinning(&fHTimeDiff, &btime, &btheta); + + btheta.Apply(fHEffOn); + btheta.Apply(fHChi2); + btheta.Apply(fHLambda); + btheta.Apply(fHN0); + btheta.Apply(fHProb); +} + +Double_t testval = 0; +Double_t testerr = 0; + +// -------------------------------------------------------------------------- +// +// Set the binnings and prepare the filling of the histogram +// +Bool_t MHEffectiveOnTime::SetupFill(const MParList *plist) +{ + fPointPos = (MPointingPos*)plist->FindObject("MPointingPos"); + if (!fPointPos) + { + *fLog << err << dbginf << "MPointingPos not found... aborting." << endl; + return kFALSE; + } + + // FIXME: Remove const-qualifier from base-class! + fTime = (MTime*)const_cast(plist)->FindCreateObj("MTime", "MTimeEffectiveOnTime"); + if (!fTime) + return kFALSE; + fParam = (MParameterDerr*)const_cast(plist)->FindCreateObj("MParameterDerr", "MEffectiveOnTime"); + if (!fParam) + return kFALSE; + + const MBinning* binsdtime = (MBinning*)plist->FindObject("BinningTimeDiff"); + const MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); + if (!binstheta || !binsdtime) + *fLog << warn << dbginf << "At least one MBinning not found... ignored." << endl; + else + { + SetBinning(&fHTimeDiff, binsdtime, binstheta); + + binstheta->Apply(fHEffOn); + binstheta->Apply(fHChi2); + binstheta->Apply(fHLambda); + binstheta->Apply(fHN0); + binstheta->Apply(fHProb); + } + + return kTRUE; +} + +Bool_t MHEffectiveOnTime::Finalize() +{ + Fit(); + fIsFinalized = kTRUE; + return kTRUE; +} + +void MHEffectiveOnTime::Fit() +{ + // nbins = number of Theta bins + const Int_t nbins = fHTimeDiff.GetNbinsY(); + + for (int i=1; i<=nbins; i++) + { + // TH1D &h = *hist->ProjectionX("Calc-theta", i, i); + TH1D *h = fHTimeDiff.ProjectionX("CalcTheta", i, i, "E"); + + const Double_t Nm = h->Integral(); + + if (Nm<=0) + continue; + + // determine range (yq[0], yq[1]) of time differences + // where fit should be performed; + // require a fraction >=xq[0] of all entries to lie below yq[0] + // and a fraction <=xq[1] of all entries to lie below yq[1]; + // within the range (yq[0], yq[1]) there must be no empty bin; + // choose pedestrian approach as long as GetQuantiles is not available + Double_t xq[2] = { 0.15, 0.95 }; + Double_t yq[2]; + + // GetQuantiles doesn't seem to be available in root 3.01/06 + h->GetQuantiles(2, yq, xq); + + // Nmdel = Nm * binwidth, with Nm = number of observed events + const Double_t Nmdel = h->Integral("width"); + + // + // Setup Poisson function for the fit: + // lambda [Hz], N0 = ideal no of evts, del = bin width of dt + // + // parameter 0 = lambda + // parameter 1 = N0*del + // + TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]); + func.SetParNames("lambda", "N0", "del"); + + func.SetParameter(0, 100); // Hz + func.SetParameter(1, Nm); + func.FixParameter(2, Nmdel/Nm); + + // options : 0 do not plot the function + // I use integral of function in bin rather than value at bin center + // R use the range specified in the function range + // Q quiet mode + h->Fit(&func, "0IRQ"); + + const Double_t chi2 = func.GetChisquare(); + const Int_t NDF = func.GetNDF(); + + // was fit successful ? + if (NDF>0 && chi2<2.5*NDF) + { + const Double_t lambda = func.GetParameter(0); + const Double_t N0 = func.GetParameter(1); + const Double_t prob = func.GetProb(); + + /* + *fLog << all << "Nm/lambda=" << Nm/lambda << " chi2/NDF="; + *fLog << (NDF ? chi2/NDF : 0.0) << " lambda="; + *fLog << lambda << " N0=" << N0 << endl; + */ + + Double_t emat[2][2]; + gMinuit->mnemat((Double_t*)emat, 2); + + const Double_t dldl = emat[0][0]; + //const Double_t dN0dN0 = emat[1][1]; + + const Double_t teff = Nm/lambda; + const Double_t dteff = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm); + + const Double_t dl = TMath::Sqrt(dldl); + + //const Double_t kappa = Nm/N0; + //const Double_t Rdead = 1.0 - kappa; + //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm); + + // the effective on time is Nm/lambda + fHEffOn.SetBinContent(i, teff); + fHEffOn.SetBinError (i, dteff); + + // plot chi2-probability of fit + fHProb.SetBinContent(i, prob*100); + + // plot chi2/NDF of fit + fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0); + + // lambda of fit + fHLambda.SetBinContent(i, lambda); + fHLambda.SetBinError (i, dl); + + // N0 of fit + fHN0.SetBinContent(i, N0); + + // Rdead (from fit) is the fraction from real time lost by the dead time + //fHRdead.SetBinContent(i, Rdead); + //fHRdead.SetBinError (i,dRdead); + } + + delete h; + } +} + + +void MHEffectiveOnTime::Paint(Option_t *opt) +{ + TVirtualPad *padsave = gPad; + + TString o(opt); + if (o==(TString)"fit") + { + TH1D *h0=0; + + padsave->cd(1); + if ((h0 = (TH1D*)gPad->FindObject("TimeDiff"))) + { + h0 = fHTimeDiff.ProjectionX("TimeDiff", -1, 9999, "E"); + if (h0->GetEntries()>0) + gPad->SetLogy(); + } + + padsave->cd(2); + if ((h0 = (TH1D*)gPad->FindObject("Theta"))) + fHTimeDiff.ProjectionY("Theta", -1, 9999, "E"); + + if (!fIsFinalized) + Fit(); + } + if (o==(TString)"paint") + { + Double_t error = 0; + for (int i=0; iGetNbins(); i++) + error += fHEffOn.GetBinError(i); + + TLatex text(0.45, 0.94, Form("T_{eff} = %.1fs \\pm %.1fs", fHEffOn.Integral(), error)); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint(); + } + gPad = padsave; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHEffectiveOnTime::Draw(Option_t *opt) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad("fit"); + + pad->Divide(2,2); + + TH1 *h; + + pad->cd(1); + gPad->SetBorderMode(0); + h = fHTimeDiff.ProjectionX("TimeDiff", -1, 9999, "E"); + h->SetTitle("Distribution of \\Delta t [s]"); + h->SetXTitle("\\Delta t [s]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->Draw(); + + pad->cd(2); + gPad->SetBorderMode(0); + h = fHTimeDiff.ProjectionY("Theta", -1, 9999, "E"); + h->SetTitle("Distribution of \\Theta [\\circ]"); + h->SetXTitle("\\Theta [\\circ]"); + h->SetYTitle("Counts"); + h->SetDirectory(NULL); + h->SetMarkerStyle(kFullDotMedium); + h->SetBit(kCanDelete); + h->GetYaxis()->SetTitleOffset(1.1); + h->Draw(); + + pad->cd(3); + gPad->SetBorderMode(0); + fHEffOn.Draw(); + AppendPad("paint"); + + pad->cd(4); + gPad->SetBorderMode(0); + fHProb.Draw(); +} + +void MHEffectiveOnTime::Calc() +{ + const Double_t val = fHEffOn.Integral(); + + Double_t error = 0; + for (int i=0; iGetNbins(); i++) + error += fHEffOn.GetBinError(i); + + fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0); + fParam->SetReadyToSave(); + + testval += fParam->GetVal(); + testerr += fParam->GetErr(); + + fEffOnTime0 = val; + fEffOnErr0 = error; + + MTime now(*fTime); + now.AddMilliSeconds(fInterval*1000); + *fLog <GetVal(), fParam->GetErr()); + *fLog << Form(" %.1f %.1f %.1f %.1f", + val, testval, error, testerr) << endl; + + fTime->AddMilliSeconds(fInterval*1000); +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHEffectiveOnTime::Fill(const MParContainer *par, const Stat_t w) +{ + const MTime *time = dynamic_cast(par); + if (!time) + { + *fLog << err << "ERROR - MHEffectiveOnTimeTime::Fill without argument or container doesn't inherit from MTime... abort." << endl; + return kFALSE; + } + + if (*fTime==MTime()) + { + *fLog << all << *time << " - " << *fTime << " " << TMath::Floor(*time) << endl; + fEffOnTime0 = 0; + fEffOnErr0 = 0; + + fParam->SetVal(0, 0); + fParam->SetReadyToSave(); + + *fTime = *time; + // fTime->MinusNull() + } + + if (*fTime+fInterval<*time) + { + Fit(); + Calc(); + } + + fHTimeDiff.Fill(*time-fLastTime, fPointPos->GetZd(), w); + fLastTime = *time; + + return kTRUE; +} + Index: /trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4887) +++ /trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4887) @@ -0,0 +1,70 @@ +#ifndef MARS_MHEffectiveOnTime +#define MARS_MHEffectiveOnTime + +#ifndef MARS_MH +#include "MH.h" +#endif +#ifndef MARS_MTime +#include "MTime.h" +#endif +#ifndef ROOT_TH1 +#include +#endif +#ifndef ROOT_TH2 +#include +#endif + +class MTime; +class MPointingPos; +class MParameterDerr; + +class MParList; + +class MHEffectiveOnTime : public MH +{ +private: + MPointingPos *fPointPos; //! + MTime fLastTime; //! + + Double_t fEffOnTime0; //! + Double_t fEffOnErr0; //! + + MTime *fTime; + MParameterDerr *fParam; + + TH2D fHTimeDiff; + TH1D fHEffOn; + TH1D fHChi2; + TH1D fHProb; + TH1D fHN0; + TH1D fHLambda; + + Bool_t fIsFinalized; + + Int_t fInterval; + + void Fit(); + void Calc(); + +public: + MHEffectiveOnTime(const char *name=NULL, const char *title=NULL); + + void SetInterval(Int_t i) { fInterval=i; } + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + const TH2D *GetHist() { return &fHTimeDiff; } + const TH2D *GetHist() const { return &fHTimeDiff; } + + TH1 *GetHistByName(const TString name) { return &fHTimeDiff; } + + void Draw(Option_t *option=""); + void Paint(Option_t *opt=""); + + ClassDef(MHEffectiveOnTime, 1) // 2D-histogram to determin effective on-time vs. theta +}; + +#endif + Index: /trunk/MagicSoft/Mars/showlog.cc =================================================================== --- /trunk/MagicSoft/Mars/showlog.cc (revision 4886) +++ /trunk/MagicSoft/Mars/showlog.cc (revision 4887) @@ -32,6 +32,4 @@ { static const TRegexp regexp("[][[][0-9]+[m]"); - - int i=0; while (1) Index: /trunk/MagicSoft/Mars/star.cc =================================================================== --- /trunk/MagicSoft/Mars/star.cc (revision 4886) +++ /trunk/MagicSoft/Mars/star.cc (revision 4887) @@ -182,5 +182,5 @@ TApplication app("Star", &argc, argv); - if (!gROOT->IsBatch() && !gClient) + if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch) { gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl;