Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 4888) +++ trunk/MagicSoft/Mars/Changelog (revision 4889) @@ -19,4 +19,36 @@ -*-*- END OF LINE -*-*- + + 2004/09/08: Thomas Bretz + + * mbadpixels/MBadPixelsCam.[h,cc], mbase/MParContainer.[h,cc], + mhvstime/MHPixVsTime.cc, mhvstime/MHSectorVsTime.cc: + - replaces ifstream by istream in AsciiRead + + * mbase/MTime.[h,cc]: + - fixed comment about SetTimeFormat + - added AsciiRead + - added AsciiWrite + - added Minus1ns + + * mfileio/MWriteAsciiFile.cc: + - write all containers if one has its SetReadyToSaveFlag set + + * mhist/MHEffectiveOnTime.[h,cc]: + - for MEffectiveOnTime fit the whole projection instead + of using the sum of the theta-bins + + * mhvstime/MHVsTime.[h,cc]: + - replaces ifstream by istream in AsciiRead + - fixed to support MStatusDisplay + - do not fill the same time twice + - added support for error bars + + * mjobs/MJStar.cc: + - replaced MReadMarsFile by MReadReports + - added MEventRateCalc and corresponding histogram + - added MHEffectiveOnTime + + 2004/09/07: Thomas Bretz Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc =================================================================== --- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc (revision 4888) +++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc (revision 4889) @@ -33,5 +33,5 @@ #include "MBadPixelsCam.h" -#include +#include #include @@ -530,5 +530,5 @@ // 1235: 17 193 292 293 // -void MBadPixelsCam::AsciiRead(ifstream &fin, UInt_t run=0) +void MBadPixelsCam::AsciiRead(istream &fin, UInt_t run=0) { Int_t len; Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.h =================================================================== --- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.h (revision 4888) +++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.h (revision 4889) @@ -44,6 +44,6 @@ Short_t GetNumMaxCluster(const MGeomCam &geom, Int_t aidx=-1) { return GetNumMaxCluster(MBadPixelsPix::kUnsuitableRun, geom, aidx); } - void AsciiRead(ifstream &fin, UInt_t run); - void AsciiRead(ifstream &fin) { AsciiRead(fin, 0); } + void AsciiRead(istream &fin, UInt_t run); + void AsciiRead(istream &fin) { AsciiRead(fin, 0); } Bool_t AsciiWrite(ostream &out, UInt_t run) const; Bool_t AsciiWrite(ostream &out) const { return AsciiWrite(out, 0); } Index: trunk/MagicSoft/Mars/mbase/MParContainer.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MParContainer.cc (revision 4888) +++ trunk/MagicSoft/Mars/mbase/MParContainer.cc (revision 4889) @@ -56,5 +56,5 @@ #include // isdigit -#include // ofstream, AsciiWrite +#include // ofstream #include // Env::Lookup @@ -296,5 +296,5 @@ // container, overload this function. // -void MParContainer::AsciiRead(ifstream &fin) +void MParContainer::AsciiRead(istream &fin) { *fLog << warn << "To use the the ascii input of " << GetName(); Index: trunk/MagicSoft/Mars/mbase/MParContainer.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MParContainer.h (revision 4888) +++ trunk/MagicSoft/Mars/mbase/MParContainer.h (revision 4889) @@ -112,5 +112,5 @@ Bool_t WriteDataMember(ostream &out, const TList *list) const; - virtual void AsciiRead(ifstream &fin); + virtual void AsciiRead(istream &fin); virtual Bool_t AsciiWrite(ostream &out) const; Index: trunk/MagicSoft/Mars/mbase/MTime.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MTime.cc (revision 4888) +++ trunk/MagicSoft/Mars/mbase/MTime.cc (revision 4889) @@ -170,9 +170,5 @@ // local time (while here we return UTC) such, that you may encounter // strange offsets. You can get rid of this by calling: -// TAxis::SetTimeFormat("[your-format] %F1995-01-01 00:00:00"); -// -// Be carefull: It seems that root takes sommer and winter time into account! -// In some circumstances you may need -// TAxis::SetTimeFormat("[your-format] %F1995-01-00 23:00:00"); +// TAxis::SetTimeFormat("[your-format] %F1995-01-01 00:00:00 GMT"); // Double_t MTime::GetAxisTime() const @@ -573,4 +569,23 @@ fTime -= 11*kHour; } + +void MTime::Minus1ns() +{ + if (fNanoSec>0) + { + fNanoSec--; + return; + } + + fTime -= 1; + fNanoSec = 999999; + + if ((Long_t)fTime>=-(Long_t)kDay*11) + return; + + fTime = 13*kDay-1; + fMjd--; +} + /* MTime MTime::operator-(const MTime &tm1) @@ -688,2 +703,13 @@ SetMjd(mean); } + +void MTime::AsciiRead(istream &fin) +{ + fin >> *this; +} + +Bool_t MTime::AsciiWrite(ostream &out) const +{ + out << *this; + return out; +} Index: trunk/MagicSoft/Mars/mbase/MTime.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MTime.h (revision 4888) +++ trunk/MagicSoft/Mars/mbase/MTime.h (revision 4889) @@ -118,4 +118,7 @@ istream &ReadBinary(istream &fin); + void AsciiRead(istream &fin); + Bool_t AsciiWrite(ostream &out) const; + // Conversion functions operator double() const; //[s] @@ -124,4 +127,5 @@ // Calculation functions void AddMilliSeconds(UInt_t ms); + void Minus1ns(); void SetMean(const MTime &t0, const MTime &t1); void SetMean(Double_t t0, Double_t t1); Index: trunk/MagicSoft/Mars/mfileio/MWriteAsciiFile.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MWriteAsciiFile.cc (revision 4888) +++ trunk/MagicSoft/Mars/mfileio/MWriteAsciiFile.cc (revision 4889) @@ -24,26 +24,25 @@ ///////////////////////////////////////////////////////////////////////////// -// // -// MWriteAsciiFile // -// // -// If you want to store a single container into an Ascii file you have // -// to use this class. You must know the name of the file you wanne write // -// (you should know it) and the name of the container you want to write. // -// This can be the name of the class or a given name, which identifies // -// the container in a parameter container list (MParList). // -// The container is written to the ascii file if its ReadyToSave flag is // -// set (MParContainer) // -// // -// You can write more than one container in one line of the file, see // -// AddContainer. // -// // -// You can also write single data members of a container (like fWidth // -// of MHillas). For more details see AddContainer. Make sure, that a // -// getter method for the data member exist. The name of the method // -// must be the same than the data member itself, but the f must be // -// replaced by a Get. // -// // +// +// MWriteAsciiFile +// +// If you want to store a single container into an Ascii file you have +// to use this class. You must know the name of the file you wanne write +// (you should know it) and the name of the container you want to write. +// This can be the name of the class or a given name, which identifies +// the container in a parameter container list (MParList). +// The container is written to the ascii file if its ReadyToSave flag is +// set (MParContainer) +// +// You can write more than one container in one line of the file, see +// AddContainer. +// +// You can also write single data members of a container (like fWidth +// of MHillas). For more details see AddContainer. Make sure, that a +// getter method for the data member exist. The name of the method +// must be the same than the data member itself, but the f must be +// replaced by a Get. +// ///////////////////////////////////////////////////////////////////////////// - #include "MWriteAsciiFile.h" @@ -51,4 +50,6 @@ #include // TMethodCall, AsciiWrite + +#include "MIter.h" #include "MDataList.h" // MDataList @@ -190,18 +191,32 @@ Bool_t MWriteAsciiFile::CheckAndWrite() { + MParContainer *obj = NULL; + + // + // Check for the Write flag + // + Bool_t write = kFALSE; + MIter Next(&fList); + while ((obj=Next())) + if (obj->IsReadyToSave()) + { + write = kTRUE; + break; + } + + // + // Do not write if not at least one ReadyToSave-flag set + // + if (!write) + return kTRUE; + Bool_t written = kFALSE; - - MParContainer *obj = NULL; - Int_t num = fList.GetEntries(); - TIter Next(&fList); - while ((obj=(MParContainer*)Next())) + Next.Reset(); + while ((obj=Next())) { - // - // Check for the Write flag - // - if (!obj->IsReadyToSave()) - continue; + if (written) + *fOut << " "; // @@ -224,5 +239,5 @@ if (num!=0) - *fLog << warn << "Warning - given number of objects doesn't fit number of written objects." << endl; + *fLog << warn << "WARNING - Number of objects written mismatch!" << endl; } return kTRUE; Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4888) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc (revision 4889) @@ -42,4 +42,5 @@ #include #include +#include #include "MTime.h" @@ -62,5 +63,5 @@ // MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title) - : /*fLastTime(0), fFirstTime(-1),*/ fIsFinalized(kFALSE), fInterval(60) + : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), fInterval(60) { // @@ -141,6 +142,7 @@ } -Double_t testval = 0; -Double_t testerr = 0; +// FIXME: Just for a preliminary check +static Double_t testval = 0; +static Double_t testerr = 0; // -------------------------------------------------------------------------- @@ -185,116 +187,153 @@ Bool_t MHEffectiveOnTime::Finalize() { - Fit(); + FitThetaBins(); + Calc(); + fIsFinalized = kTRUE; + return kTRUE; } -void MHEffectiveOnTime::Fit() -{ +void MHEffectiveOnTime::FitThetaBins() +{ + const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999)); + // nbins = number of Theta bins const Int_t nbins = fHTimeDiff.GetNbinsY(); + TH1D *h=0; 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) + h = fHTimeDiff.ProjectionX(name, i, i, "E"); + + Double_t res[7]; + if (!FitH(h, res)) 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); - } - + // the effective on time is Nm/lambda + fHEffOn.SetBinContent(i, res[0]); + fHEffOn.SetBinError (i, res[1]); + + // plot chi2-probability of fit + fHProb.SetBinContent(i, res[2]); + + // plot chi2/NDF of fit + fHChi2.SetBinContent(i, res[3]); + + // lambda of fit + fHLambda.SetBinContent(i, res[4]); + fHLambda.SetBinError (i, res[5]); + + // N0 of fit + fHN0.SetBinContent(i, res[6]); + + // Rdead (from fit) is the fraction from real time lost by the dead time + //fHRdead.SetBinContent(i, Rdead); + //fHRdead.SetBinError (i,dRdead); + } + + // Histogram is reused via gROOT->FindObject() + // Need to be deleted only once + if (h) delete h; - } -} - +} + +Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res) const +{ + const Double_t Nm = h->Integral(); + + // FIXME: Do fit only if contents of bin has changed + if (Nm<=0) + return kFALSE; + + // 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.05, 0.95 }; + Double_t yq[2]; + 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) + return kFALSE; + + 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 + res[0] = teff; + res[1] = dteff; + + // plot chi2-probability of fit + res[2] = prob*100; + + // plot chi2/NDF of fit + res[3] = NDF ? chi2/NDF : 0.0; + + // lambda of fit + res[4] = lambda; + res[5] = dl; + + // N0 of fit + res[6] = N0; + + // Rdead (from fit) is the fraction from real time lost by the dead time + //fHRdead.SetBinContent(i, Rdead); + //fHRdead.SetBinError (i,dRdead); + + return kTRUE; +} void MHEffectiveOnTime::Paint(Option_t *opt) @@ -320,5 +359,5 @@ if (!fIsFinalized) - Fit(); + FitThetaBins(); } if (o==(TString)"paint") @@ -386,9 +425,16 @@ void MHEffectiveOnTime::Calc() { - const Double_t val = fHEffOn.Integral(); - - Double_t error = 0; - for (int i=0; iGetNbins(); i++) - error += fHEffOn.GetBinError(i); + TH1D *h = fHTimeDiff.ProjectionX("", -1, 99999, "E"); + h->SetDirectory(0); + + Double_t res[7]; + Bool_t rc = FitH(h, res); + delete h; + + if (!rc) + return; + + const Double_t val = res[0]; + const Double_t error = res[1]; fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0); @@ -403,4 +449,5 @@ MTime now(*fTime); now.AddMilliSeconds(fInterval*1000); + *fLog <MinusNull() + + // Make this just a ns before the first event + fTime->Minus1ns(); } if (*fTime+fInterval<*time) { - Fit(); + FitThetaBins(); Calc(); } Index: trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4888) +++ trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h (revision 4889) @@ -24,12 +24,12 @@ { private: - MPointingPos *fPointPos; //! - MTime fLastTime; //! + MPointingPos *fPointPos; //! + MTime fLastTime; //! - Double_t fEffOnTime0; //! - Double_t fEffOnErr0; //! + Double_t fEffOnTime0; //! + Double_t fEffOnErr0; //! - MTime *fTime; - MParameterDerr *fParam; + MTime *fTime; //! + MParameterDerr *fParam; //! TH2D fHTimeDiff; @@ -44,5 +44,6 @@ Int_t fInterval; - void Fit(); + Bool_t FitH(TH1D *h, Double_t *res) const; + void FitThetaBins(); void Calc(); Index: trunk/MagicSoft/Mars/mhvstime/MHPixVsTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhvstime/MHPixVsTime.cc (revision 4888) +++ trunk/MagicSoft/Mars/mhvstime/MHPixVsTime.cc (revision 4889) @@ -194,5 +194,5 @@ { TAxis *axe = h->GetXaxis(); - axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00"); + axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); axe->SetTimeDisplay(1); axe->SetLabelSize(0.033); Index: trunk/MagicSoft/Mars/mhvstime/MHSectorVsTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhvstime/MHSectorVsTime.cc (revision 4888) +++ trunk/MagicSoft/Mars/mhvstime/MHSectorVsTime.cc (revision 4889) @@ -247,5 +247,5 @@ { TAxis *axe = h->GetXaxis(); - axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00"); + axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); axe->SetTimeDisplay(1); axe->SetLabelSize(0.025); Index: trunk/MagicSoft/Mars/mhvstime/MHVsTime.cc =================================================================== --- trunk/MagicSoft/Mars/mhvstime/MHVsTime.cc (revision 4888) +++ trunk/MagicSoft/Mars/mhvstime/MHVsTime.cc (revision 4889) @@ -49,5 +49,5 @@ #include -#include +#include #include "MLog.h" @@ -71,7 +71,8 @@ // see MDataChain. // -MHVsTime::MHVsTime(const char *rule) - : fGraph(NULL), fData(NULL), fScale(1), fMaxPts(-1), fUseEventNumber(0) -{ +MHVsTime::MHVsTime(const char *rule, const char *error) + : fGraph(NULL), fData(NULL), fError(0), fScale(1), fMaxPts(-1), + fNumEvents(1), fUseEventNumber(0) +{ fName = gsDefName; fTitle = gsDefTitle; @@ -80,8 +81,12 @@ return; - fGraph = new TGraph; fData = new MDataChain(rule); - fGraph->SetMarkerStyle(kFullDotMedium); + if (error) + fError = new MDataChain(error); + + fGraph = error ? new TGraphErrors : new TGraph; + fGraph->SetPoint(0, 0, 0); // Dummy point! + fGraph->SetEditable(); // Used as flag: First point? yes/no } @@ -116,15 +121,19 @@ Bool_t MHVsTime::SetupFill(const MParList *plist) { - // reset histogram (necessary if the same eventloop is run more than once) - //fGraph->Reset(); - - if (fData && !fData->PreProcess(plist)) + if (!fGraph || !fData) + { + *fLog << err << "ERROR - MHVsTime cannot be used with its default constructor!" << endl; return kFALSE; - - if (fGraph) - { - delete fGraph; - fGraph = new TGraph; - } + } + + if (!fData->PreProcess(plist)) + return kFALSE; + + if (fError && !fError->PreProcess(plist)) + return kFALSE; + + fGraph->Set(1); + fGraph->SetPoint(0, 0, 0); // Dummy point! + fGraph->SetEditable(); // Used as flag: First point? yes/no TString title(fData ? GetRule() : (TString)"Histogram"); @@ -183,23 +192,69 @@ if (!*tm) return kTRUE; + + // Do not fill events with equal time + if (*tm==fLast || *tm==MTime()) + return kTRUE; + + fLast = *tm; + t = tm->GetAxisTime(); } - const Double_t v = fData->GetValue()*fScale; - - if (fMaxPts>0 && fGraph->GetN()>fMaxPts) - fGraph->RemovePoint(0); - - fMean += v; + const Double_t v = fData->GetValue(); + const Double_t e = fError ? fError->GetValue() : 0; + + //*fLog << all << "ADD " << v << " " << e << endl; + + fMean += v; + fMeanErr += e; fN++; if (fN==fNumEvents) { - fGraph->SetPoint(fGraph->GetN(), t, fMean/fN); + if (fMaxPts>0 && fGraph->GetN()>fMaxPts || fGraph->IsEditable()) + { + fGraph->RemovePoint(0); + fGraph->SetEditable(kFALSE); + } + + fGraph->SetPoint(fGraph->GetN(), t, fMean/fN*fScale); + + if (fError) + static_cast(fGraph)->SetPointError(fGraph->GetN()-1, 0, fMeanErr/fN*fScale); + fMean = 0; + fMeanErr = 0; fN = 0; } return kTRUE; +} + +void MHVsTime::Paint(Option_t *opt) +{ + // SetPoint deletes the histogram! + if (fUseEventNumber) + fGraph->GetHistogram()->SetXTitle("Event Number"); + else + { + fGraph->GetHistogram()->SetXTitle("Time"); + fGraph->GetHistogram()->GetXaxis()->SetLabelSize(0.033); + fGraph->GetHistogram()->GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT"); + fGraph->GetHistogram()->GetXaxis()->SetTimeDisplay(1); + } + + fGraph->GetHistogram()->SetMarkerStyle(kFullDotMedium); + fGraph->GetHistogram()->SetYTitle(fAxisTitle.IsNull() ? GetRule() : fAxisTitle); + if (fTitle!=gsDefTitle) + fGraph->GetHistogram()->SetTitle(fTitle); + + if (TestBit(kIsLogy)) + gPad->SetLogy(); + + // This is a workaround if the TGraph has only one point. + // Otherwise MStatusDisplay::Update hangs. + gPad->GetListOfPrimitives()->Remove(fGraph); + gPad->GetListOfPrimitives()->Add(fGraph, fGraph->GetN()<2 ? "A" : opt); } @@ -212,4 +267,7 @@ void MHVsTime::Draw(Option_t *opt) { + if (!fGraph) + return; + if (fGraph->GetN()==0) return; @@ -218,18 +276,5 @@ pad->SetBorderMode(0); - AppendPad(""); - TString str(opt); - - if (fUseEventNumber) - fGraph->GetHistogram()->SetXTitle("Event Number"); - else - { - fGraph->GetHistogram()->SetXTitle("Time"); - fGraph->GetHistogram()->GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00"); - fGraph->GetHistogram()->GetXaxis()->SetTimeDisplay(1); - fGraph->GetHistogram()->GetXaxis()->SetLabelSize(0.033); - } - fGraph->GetHistogram()->SetYTitle(GetRule()); if (!str.Contains("A")) @@ -244,10 +289,6 @@ } + AppendPad(str); fGraph->Draw(str); - if (fGraph->TestBit(kIsLogy)) - pad->SetLogy(); - - pad->Modified(); - pad->Update(); } Index: trunk/MagicSoft/Mars/mhvstime/MHVsTime.h =================================================================== --- trunk/MagicSoft/Mars/mhvstime/MHVsTime.h (revision 4888) +++ trunk/MagicSoft/Mars/mhvstime/MHVsTime.h (revision 4889) @@ -4,4 +4,8 @@ #ifndef MARS_MH #include "MH.h" +#endif + +#ifndef MARS_MTime +#include "MTime.h" #endif @@ -15,15 +19,19 @@ TGraph *fGraph; // Histogram to fill MDataChain *fData; // Object from which the data is filled + MDataChain *fError; // Object from which the error is filled Double_t fScale; // Scale for axis (eg unit) Int_t fMaxPts; // Maximum number of data points Int_t fNumEvents; // Number of events to average + Double_t fMean; //! Mean value - Int_t fN; + Double_t fMeanErr; //! Mean error + Int_t fN; //! Number of entries in fMean + MTime fLast; //! For checks + TString fAxisTitle; enum { - kIsLogy = BIT(18), - kUseEventNumber = BIT(20) + kIsLogy = BIT(18) }; @@ -31,5 +39,5 @@ public: - MHVsTime(const char *rule=NULL); + MHVsTime(const char *rule=NULL, const char *ruleerr=NULL); ~MHVsTime(); @@ -40,4 +48,7 @@ void SetName(const char *name); void SetTitle(const char *title); + + void SetLogy(Bool_t b=kTRUE) { b ? SetBit(kIsLogy) : ResetBit(kIsLogy); } + void SetAxisTitle(const char *y) { fAxisTitle=y; } Bool_t SetupFill(const MParList *pList); @@ -59,4 +70,5 @@ void Draw(Option_t *opt=NULL); + void Paint(Option_t *opt=NULL); MParContainer *New() const; Index: trunk/MagicSoft/Mars/mjobs/MJStar.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJStar.cc (revision 4888) +++ trunk/MagicSoft/Mars/mjobs/MJStar.cc (revision 4889) @@ -43,12 +43,19 @@ #include "MStatusDisplay.h" +#include "MH3.h" +#include "MHVsTime.h" #include "MHCamEvent.h" - -#include "MReadMarsFile.h" +#include "MBinning.h" + +#include "MReadReports.h" #include "MGeomApply.h" +#include "MEventRateCalc.h" #include "MImgCleanStd.h" #include "MHillasCalc.h" #include "MFillH.h" #include "MWriteRootFile.h" + +#include "MPointingPosCalc.h" +//#include "MSrcPosFromModel.h" ClassImp(MJStar); @@ -138,26 +145,68 @@ plist.AddToList(&tlist); - MReadMarsFile read("Events"); - read.DisableAutoScheme(); + MReadReports read; + read.AddTree("Events", "MTime.", kTRUE); + read.AddTree("Drive"); + //read.AddTree("Trigger"); + //read.AddTree("Camera"); + //read.AddTree("CC"); + //read.AddTree("Currents"); read.AddFiles(iter); - //read.AddFiles(fnamein); + + // ------------------ Setup general tasks ---------------- MGeomApply apply; // Only necessary to craete geometry + MEventRateCalc rate; +/* + MEventRateCalc rate10000; + rate10000.SetNameEventRate("MEventRate10000"); + rate10000.SetNumEvents(10000); + */ //MBadPixelsMerge merge(&badpix); MImgCleanStd clean; MHillasCalc hcalc; + // ------------------ Setup histograms and fill tasks ---------------- MHCamEvent evt0("Cleaned"); evt0.SetType(0); + + MH3 h1("MEventRate.fRate"); + h1.SetName("MHEventRate"); + h1.SetLogy(); +/* + MH3 h12("MEventRate10000.fRate"); + h12.SetName("MHEventRate"); + h12.SetLogy(); + */ + MBinning b1("BinningMHEventRate"); + b1.SetEdges(150, 0, 1500); + plist.AddToList(&b1); + + MHVsTime h2("MEffectiveOnTime.fVal", "MEffectiveOnTime.fErr"); + h2.SetAxisTitle("T_{eff}"); + h2.SetTitle("Effective On-Time T_{eff} vs. Time"); + MFillH fill0(&evt0, "MCerPhotEvt", "FillCerPhotEvt"); MFillH fill1("MHHillas", "MHillas", "FillHillas"); - MFillH fill2("MHHillasExt", "MHillasExt", "FillHillasExt"); + MFillH fill2("MHHillasExt", "", "FillHillasExt"); MFillH fill3("MHHillasSrc", "MHillasSrc", "FillHillasSrc"); MFillH fill4("MHImagePar", "MImagePar", "FillImagePar"); MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar"); - MFillH fill6("MHCerPhot"); - - MWriteRootFile write(2, "images/{s/_Y_/_I_}"); - write.AddContainer("MMcEvt", "Events", kFALSE); + MFillH fill6("MHImageParTime","MImageParTime","FillImageParTime"); + MFillH fill7("MHNewImagePar2","MNewImagePar2","FillNewImagePar2"); + MFillH fill8(&h1, "", "FillEventRate"); + MFillH fill9("MHEffectiveOnTime", "MTime", "FillEffOnTime"); + MFillH filla(&h2, "MTimeEffectiveOnTime", "FillEffOnTimeVsTime"); + //MFillH fillb(&h12, "", "FillEvtRate2"); + //MFillH fill9("MHCerPhot"); + + fill8.SetNameTab("EvtRate"); + fill9.SetNameTab("EffOnTime"); + fill9.SetNameTab("EffOnVsTime"); + + // ------------------ Setup write task ---------------- + + MWriteRootFile write(2, Form("%s{s/_Y_/_I_}", fPathOut.Data()), fOverwrite); + // Data write.AddContainer("MHillas", "Events"); write.AddContainer("MHillasExt", "Events"); @@ -165,22 +214,59 @@ write.AddContainer("MImagePar", "Events"); write.AddContainer("MNewImagePar", "Events"); + write.AddContainer("MNewImagePar2", "Events"); + write.AddContainer("MImageParTime", "Events"); write.AddContainer("MTime", "Events"); write.AddContainer("MRawEvtHeader", "Events"); - write.AddContainer("MRawRunHeader", "RunHeaders"); - write.AddContainer("MBadPixelsCam", "RunHeaders"); - write.AddContainer("MGeomCam", "RunHeaders"); + // Monte Carlo + write.AddContainer("MMcEvt", "Events", kFALSE); + write.AddContainer("MMcTrig", "Events", kFALSE); + // Run Header + write.AddContainer("MRawRunHeader", "RunHeaders"); + write.AddContainer("MBadPixelsCam", "RunHeaders"); + write.AddContainer("MGeomCam", "RunHeaders"); //write.AddContainer("MObservatory", "RunHeaders"); + // Monte Carlo Headers + write.AddContainer("MMcTrigHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcConfigRunHeader", "RunHeaders", kFALSE); + write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE); + // Drive + //write.AddContainer("MSrcPosCam", "Drive"); + write.AddContainer("MPointingPos", "Drive"); + write.AddContainer("MReportDrive", "Drive"); + write.AddContainer("MTimeDrive", "Drive"); + // Effective On Time + write.AddContainer("MEffectiveOnTime", "EffectiveOnTime"); + write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime"); + + MTaskList tlist2; + tlist2.AddToList(&apply); + tlist2.AddToList(&rate); + //tlist2.AddToList(&rate10000); + tlist2.AddToList(&fill8); + tlist2.AddToList(&fill9); + tlist2.AddToList(&filla); + //tlist2.AddToList(&fillb); + tlist2.AddToList(&clean); + tlist2.AddToList(&fill0); + tlist2.AddToList(&hcalc); + tlist2.AddToList(&fill1); + tlist2.AddToList(&fill2); + tlist2.AddToList(&fill3); + tlist2.AddToList(&fill4); + tlist2.AddToList(&fill5); + tlist2.AddToList(&fill6); + tlist2.AddToList(&fill7); + //tlist2.AddToList(&fill9); + + MPointingPosCalc pcalc; + //MSrcPosFromModel srcpos; + + MTaskList tlist3; + tlist3.AddToList(&pcalc); + //tlist3.AddToList(&srcpos); tlist.AddToList(&read); - tlist.AddToList(&apply); - tlist.AddToList(&clean); - tlist.AddToList(&fill0); - tlist.AddToList(&hcalc); - tlist.AddToList(&fill1); - tlist.AddToList(&fill2); - tlist.AddToList(&fill3); - tlist.AddToList(&fill4); - tlist.AddToList(&fill5); - //tlist.AddToList(&fill6); + tlist.AddToList(&tlist3, "Drive"); + tlist.AddToList(&tlist2, "Events"); tlist.AddToList(&write);