Changeset 4889 for trunk


Ignore:
Timestamp:
09/08/04 18:49:00 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r4887 r4889  
    1919
    2020                                                 -*-*- END OF LINE -*-*-
     21
     22 2004/09/08: Thomas Bretz
     23 
     24   * mbadpixels/MBadPixelsCam.[h,cc], mbase/MParContainer.[h,cc],
     25     mhvstime/MHPixVsTime.cc, mhvstime/MHSectorVsTime.cc:
     26     - replaces ifstream by istream in AsciiRead
     27
     28   * mbase/MTime.[h,cc]:
     29     - fixed comment about SetTimeFormat
     30     - added AsciiRead
     31     - added AsciiWrite
     32     - added Minus1ns
     33
     34   * mfileio/MWriteAsciiFile.cc:
     35     - write all containers if one has its SetReadyToSaveFlag set
     36
     37   * mhist/MHEffectiveOnTime.[h,cc]:
     38     - for MEffectiveOnTime fit the whole projection instead
     39       of using the sum of the theta-bins
     40
     41   * mhvstime/MHVsTime.[h,cc]:
     42     - replaces ifstream by istream in AsciiRead
     43     - fixed to support MStatusDisplay
     44     - do not fill the same time twice
     45     - added support for error bars
     46
     47   * mjobs/MJStar.cc:
     48     - replaced MReadMarsFile by MReadReports
     49     - added MEventRateCalc and corresponding histogram
     50     - added MHEffectiveOnTime
     51
     52
    2153
    2254 2004/09/07: Thomas Bretz

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.cc

    r4793 r4889  
    3333#include "MBadPixelsCam.h"
    3434
    35 #include <fstream>
     35#include <iostream>
    3636
    3737#include <TClonesArray.h>
    … …  
    530530//   1235: 17 193 292 293
    531531//
    532 void MBadPixelsCam::AsciiRead(ifstream &fin, UInt_t run=0)
     532void MBadPixelsCam::AsciiRead(istream &fin, UInt_t run=0)
    533533{
    534534    Int_t len;

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCam.h

    r4679 r4889  
    4444    Short_t GetNumMaxCluster(const MGeomCam &geom, Int_t aidx=-1) { return GetNumMaxCluster(MBadPixelsPix::kUnsuitableRun, geom, aidx); }
    4545
    46     void   AsciiRead(ifstream &fin, UInt_t run);
    47     void   AsciiRead(ifstream &fin) { AsciiRead(fin, 0); }
     46    void   AsciiRead(istream &fin, UInt_t run);
     47    void   AsciiRead(istream &fin) { AsciiRead(fin, 0); }
    4848    Bool_t AsciiWrite(ostream &out, UInt_t run) const;
    4949    Bool_t AsciiWrite(ostream &out) const { return AsciiWrite(out, 0); }

  • trunk/MagicSoft/Mars/mbase/MParContainer.cc

    r4601 r4889  
    5656
    5757#include <ctype.h>        // isdigit
    58 #include <fstream>        // ofstream, AsciiWrite
     58#include <fstream>        // ofstream
    5959
    6060#include <TEnv.h>         // Env::Lookup
    … …  
    296296//  container, overload this function.
    297297//
    298 void MParContainer::AsciiRead(ifstream &fin)
     298void MParContainer::AsciiRead(istream &fin)
    299299{
    300300    *fLog << warn << "To use the the ascii input of " << GetName();

  • trunk/MagicSoft/Mars/mbase/MParContainer.h

    r4828 r4889  
    112112    Bool_t WriteDataMember(ostream &out, const TList *list) const;
    113113
    114     virtual void AsciiRead(ifstream &fin);
     114    virtual void AsciiRead(istream &fin);
    115115    virtual Bool_t AsciiWrite(ostream &out) const;
    116116

  • trunk/MagicSoft/Mars/mbase/MTime.cc

    r4887 r4889  
    170170// local time (while here we return UTC) such, that you may encounter
    171171// strange offsets. You can get rid of this by calling:
    172 //    TAxis::SetTimeFormat("[your-format] %F1995-01-01 00:00:00");
    173 //
    174 // Be carefull: It seems that root takes sommer and winter time into account!
    175 //              In some circumstances you may need
    176 //    TAxis::SetTimeFormat("[your-format] %F1995-01-00 23:00:00");
     172//    TAxis::SetTimeFormat("[your-format] %F1995-01-01 00:00:00 GMT");
    177173//
    178174Double_t MTime::GetAxisTime() const
    … …  
    573569    fTime -= 11*kHour;
    574570}
     571
     572void MTime::Minus1ns()
     573{
     574    if (fNanoSec>0)
     575    {
     576        fNanoSec--;
     577        return;
     578    }
     579
     580    fTime -= 1;
     581    fNanoSec = 999999;
     582
     583    if ((Long_t)fTime>=-(Long_t)kDay*11)
     584        return;
     585
     586    fTime = 13*kDay-1;
     587    fMjd--;
     588}   
     589
    575590/*
    576591MTime MTime::operator-(const MTime &tm1)
    … …  
    688703    SetMjd(mean);
    689704}
     705
     706void MTime::AsciiRead(istream &fin)
     707{
     708    fin >> *this;
     709}
     710
     711Bool_t MTime::AsciiWrite(ostream &out) const
     712{
     713    out << *this;
     714    return out;
     715}

  • trunk/MagicSoft/Mars/mbase/MTime.h

    r4887 r4889  
    118118    istream &ReadBinary(istream &fin);
    119119
     120    void AsciiRead(istream &fin);
     121    Bool_t AsciiWrite(ostream &out) const;
     122
    120123    // Conversion functions
    121124    operator double() const;   //[s]
    … …  
    124127    // Calculation functions
    125128    void AddMilliSeconds(UInt_t ms);
     129    void Minus1ns();
    126130    void SetMean(const MTime &t0, const MTime &t1);
    127131    void SetMean(Double_t t0, Double_t t1);

  • trunk/MagicSoft/Mars/mfileio/MWriteAsciiFile.cc

    r4694 r4889  
    2424
    2525/////////////////////////////////////////////////////////////////////////////
    26 //                                                                         //
    27 // MWriteAsciiFile                                                         //
    28 //                                                                         //
    29 // If you want to store a single container into an Ascii file you have     //
    30 // to use this class. You must know the name of the file you wanne write   //
    31 // (you should know it) and the name of the container you want to write.   //
    32 // This can be the name of the class or a given name, which identifies     //
    33 // the container in a parameter container list (MParList).                 //
    34 // The container is written to the ascii file if its ReadyToSave flag is   //
    35 // set (MParContainer)                                                     //
    36 //                                                                         //
    37 // You can write more than one container in one line of the file, see      //
    38 // AddContainer.                                                           //
    39 //                                                                         //
    40 // You can also write single data members of a container (like fWidth      //
    41 // of MHillas). For more details see AddContainer. Make sure, that a       //
    42 // getter method for the data member exist. The name of the method         //
    43 // must be the same than the data member itself, but the f must be         //
    44 // replaced by a Get.                                                      //
    45 //                                                                         //
     26//
     27// MWriteAsciiFile
     28//
     29// If you want to store a single container into an Ascii file you have
     30// to use this class. You must know the name of the file you wanne write
     31// (you should know it) and the name of the container you want to write.
     32// This can be the name of the class or a given name, which identifies
     33// the container in a parameter container list (MParList).
     34// The container is written to the ascii file if its ReadyToSave flag is
     35// set (MParContainer)
     36//
     37// You can write more than one container in one line of the file, see
     38// AddContainer.
     39//
     40// You can also write single data members of a container (like fWidth
     41// of MHillas). For more details see AddContainer. Make sure, that a
     42// getter method for the data member exist. The name of the method
     43// must be the same than the data member itself, but the f must be
     44// replaced by a Get.
     45//
    4646/////////////////////////////////////////////////////////////////////////////
    47 
    4847#include "MWriteAsciiFile.h"
    4948
    … …  
    5150
    5251#include <TMethodCall.h> // TMethodCall, AsciiWrite
     52
     53#include "MIter.h"
    5354
    5455#include "MDataList.h"   // MDataList
    … …  
    190191Bool_t MWriteAsciiFile::CheckAndWrite()
    191192{
     193    MParContainer *obj = NULL;
     194
     195    //
     196    // Check for the Write flag
     197    //
     198    Bool_t write = kFALSE;
     199    MIter Next(&fList);
     200    while ((obj=Next()))
     201        if (obj->IsReadyToSave())
     202        {
     203            write = kTRUE;
     204            break;
     205        }
     206
     207    //
     208    // Do not write if not at least one ReadyToSave-flag set
     209    //
     210    if (!write)
     211        return kTRUE;
     212
    192213    Bool_t written = kFALSE;
    193 
    194     MParContainer *obj = NULL;
    195 
    196214    Int_t num = fList.GetEntries();
    197215
    198     TIter Next(&fList);
    199     while ((obj=(MParContainer*)Next()))
     216    Next.Reset();
     217    while ((obj=Next()))
    200218    {
    201         //
    202         // Check for the Write flag
    203         //
    204         if (!obj->IsReadyToSave())
    205             continue;
     219        if (written)
     220            *fOut << " ";
    206221
    207222        //
    … …  
    224239
    225240        if (num!=0)
    226             *fLog << warn << "Warning - given number of objects doesn't fit number of written objects." << endl;
     241            *fLog << warn << "WARNING - Number of objects written mismatch!" << endl;
    227242    }
    228243    return kTRUE;

  • trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.cc

    r4887 r4889  
    4242#include <TCanvas.h>
    4343#include <TMinuit.h>
     44#include <TRandom.h>
    4445
    4546#include "MTime.h"
    … …  
    6263//
    6364MHEffectiveOnTime::MHEffectiveOnTime(const char *name, const char *title)
    64     : /*fLastTime(0), fFirstTime(-1),*/ fIsFinalized(kFALSE), fInterval(60)
     65    : fPointPos(0), fTime(0), fParam(0), fIsFinalized(kFALSE), fInterval(60)
    6566{
    6667    //
    … …  
    141142}
    142143
    143 Double_t testval = 0;
    144 Double_t testerr = 0;
     144// FIXME: Just for a preliminary check
     145static Double_t testval = 0;
     146static Double_t testerr = 0;
    145147
    146148// --------------------------------------------------------------------------
    … …  
    185187Bool_t MHEffectiveOnTime::Finalize()
    186188{
    187     Fit();
     189    FitThetaBins();
     190    Calc();
     191
    188192    fIsFinalized = kTRUE;
     193
    189194    return kTRUE;
    190195}
    191196
    192 void MHEffectiveOnTime::Fit()
    193 {
     197void MHEffectiveOnTime::FitThetaBins()
     198{
     199    const TString name = Form("CalcTheta%d", (UInt_t)gRandom->Uniform(999999999));
     200
    194201    // nbins = number of Theta bins
    195202    const Int_t nbins = fHTimeDiff.GetNbinsY();
    196203
     204    TH1D *h=0;
    197205    for (int i=1; i<=nbins; i++)
    198206    {
    199207        //        TH1D &h = *hist->ProjectionX("Calc-theta", i, i);
    200         TH1D *h = fHTimeDiff.ProjectionX("CalcTheta", i, i, "E");
    201 
    202         const Double_t Nm = h->Integral();
    203 
    204         if (Nm<=0)
     208        h = fHTimeDiff.ProjectionX(name, i, i, "E");
     209
     210        Double_t res[7];
     211        if (!FitH(h, res))
    205212            continue;
    206213
    207         // determine range (yq[0], yq[1]) of time differences
    208         // where fit should be performed;
    209         // require a fraction >=xq[0] of all entries to lie below yq[0]
    210         //     and a fraction <=xq[1] of all entries to lie below yq[1];
    211         // within the range (yq[0], yq[1]) there must be no empty bin;
    212         // choose pedestrian approach as long as GetQuantiles is not available
    213         Double_t xq[2] = { 0.15, 0.95 };
    214         Double_t yq[2];
    215 
    216         // GetQuantiles doesn't seem to be available in root 3.01/06
    217         h->GetQuantiles(2, yq, xq);
    218 
    219         // Nmdel = Nm * binwidth,  with Nm = number of observed events
    220         const Double_t Nmdel = h->Integral("width");
    221 
    222         //
    223         // Setup Poisson function for the fit:
    224         // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
    225         //
    226         // parameter 0 = lambda
    227         // parameter 1 = N0*del     
    228         //
    229         TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
    230         func.SetParNames("lambda", "N0", "del");
    231 
    232         func.SetParameter(0, 100);       // Hz
    233         func.SetParameter(1, Nm);
    234         func.FixParameter(2, Nmdel/Nm);
    235 
    236         // options : 0  do not plot the function
    237         //           I  use integral of function in bin rather than value at bin center
    238         //           R  use the range specified in the function range
    239         //           Q  quiet mode
    240         h->Fit(&func, "0IRQ");
    241 
    242         const Double_t chi2   = func.GetChisquare();
    243         const Int_t    NDF    = func.GetNDF();
    244 
    245         // was fit successful ?
    246         if (NDF>0 && chi2<2.5*NDF)
    247         {
    248             const Double_t lambda = func.GetParameter(0);
    249             const Double_t N0     = func.GetParameter(1);
    250             const Double_t prob   = func.GetProb();
    251 
    252             /*
    253              *fLog << all << "Nm/lambda=" << Nm/lambda << "  chi2/NDF=";
    254              *fLog << (NDF ? chi2/NDF : 0.0) << "  lambda=";
    255              *fLog << lambda << "  N0=" << N0 << endl;
    256              */
    257 
    258             Double_t emat[2][2];
    259             gMinuit->mnemat((Double_t*)emat, 2);
    260 
    261             const Double_t dldl   = emat[0][0];
    262             //const Double_t dN0dN0 = emat[1][1];
    263 
    264             const Double_t teff   = Nm/lambda;
    265             const Double_t dteff  = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
    266 
    267             const Double_t dl     = TMath::Sqrt(dldl);
    268 
    269             //const Double_t kappa  = Nm/N0;
    270             //const Double_t Rdead  = 1.0 - kappa;
    271             //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
    272 
    273             // the effective on time is Nm/lambda
    274             fHEffOn.SetBinContent(i, teff);
    275             fHEffOn.SetBinError  (i, dteff);
    276 
    277             // plot chi2-probability of fit
    278             fHProb.SetBinContent(i, prob*100);
    279 
    280             // plot chi2/NDF of fit
    281             fHChi2.SetBinContent(i, NDF ? chi2/NDF : 0.0);
    282 
    283             // lambda of fit
    284             fHLambda.SetBinContent(i, lambda);
    285             fHLambda.SetBinError  (i,     dl);
    286 
    287             // N0 of fit
    288             fHN0.SetBinContent(i, N0);
    289 
    290             // Rdead (from fit) is the fraction from real time lost by the dead time
    291             //fHRdead.SetBinContent(i, Rdead);
    292             //fHRdead.SetBinError  (i,dRdead);
    293         }
    294 
     214        // the effective on time is Nm/lambda
     215        fHEffOn.SetBinContent(i, res[0]);
     216        fHEffOn.SetBinError  (i, res[1]);
     217
     218        // plot chi2-probability of fit
     219        fHProb.SetBinContent(i, res[2]);
     220
     221        // plot chi2/NDF of fit
     222        fHChi2.SetBinContent(i, res[3]);
     223
     224        // lambda of fit
     225        fHLambda.SetBinContent(i, res[4]);
     226        fHLambda.SetBinError  (i, res[5]);
     227
     228        // N0 of fit
     229        fHN0.SetBinContent(i, res[6]);
     230
     231        // Rdead (from fit) is the fraction from real time lost by the dead time
     232        //fHRdead.SetBinContent(i, Rdead);
     233        //fHRdead.SetBinError  (i,dRdead);
     234    }
     235
     236    // Histogram is reused via gROOT->FindObject()
     237    // Need to be deleted only once
     238    if (h)
    295239        delete h;
    296     }
    297 }
    298 
     240}
     241
     242Bool_t MHEffectiveOnTime::FitH(TH1D *h, Double_t *res) const
     243{
     244    const Double_t Nm = h->Integral();
     245
     246    // FIXME: Do fit only if contents of bin has changed
     247    if (Nm<=0)
     248        return kFALSE;
     249
     250    // determine range (yq[0], yq[1]) of time differences
     251    // where fit should be performed;
     252    // require a fraction >=xq[0] of all entries to lie below yq[0]
     253    //     and a fraction <=xq[1] of all entries to lie below yq[1];
     254    // within the range (yq[0], yq[1]) there must be no empty bin;
     255    // choose pedestrian approach as long as GetQuantiles is not available
     256    Double_t xq[2] = { 0.05, 0.95 };
     257    Double_t yq[2];
     258    h->GetQuantiles(2, yq, xq);
     259
     260    // Nmdel = Nm * binwidth,  with Nm = number of observed events
     261    const Double_t Nmdel = h->Integral("width");
     262
     263    //
     264    // Setup Poisson function for the fit:
     265    // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
     266    //
     267    // parameter 0 = lambda
     268    // parameter 1 = N0*del
     269    //
     270    TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
     271    func.SetParNames("lambda", "N0", "del");
     272
     273    func.SetParameter(0, 100);       // Hz
     274    func.SetParameter(1, Nm);
     275    func.FixParameter(2, Nmdel/Nm);
     276
     277    // options : 0  do not plot the function
     278    //           I  use integral of function in bin rather than value at bin center
     279    //           R  use the range specified in the function range
     280    //           Q  quiet mode
     281    h->Fit(&func, "0IRQ");
     282
     283    const Double_t chi2   = func.GetChisquare();
     284    const Int_t    NDF    = func.GetNDF();
     285
     286    // was fit successful ?
     287    if (NDF<=0 || chi2>=2.5*NDF)
     288        return kFALSE;
     289
     290    const Double_t lambda = func.GetParameter(0);
     291    const Double_t N0     = func.GetParameter(1);
     292    const Double_t prob   = func.GetProb();
     293
     294    /*
     295     *fLog << all << "Nm/lambda=" << Nm/lambda << "  chi2/NDF=";
     296     *fLog << (NDF ? chi2/NDF : 0.0) << "  lambda=";
     297     *fLog << lambda << "  N0=" << N0 << endl;
     298     */
     299
     300    Double_t emat[2][2];
     301    gMinuit->mnemat((Double_t*)emat, 2);
     302
     303    const Double_t dldl   = emat[0][0];
     304    //const Double_t dN0dN0 = emat[1][1];
     305
     306    const Double_t teff   = Nm/lambda;
     307    const Double_t dteff  = teff * TMath::Sqrt(dldl/(lambda*lambda) + 1.0/Nm);
     308
     309    const Double_t dl     = TMath::Sqrt(dldl);
     310
     311    //const Double_t kappa  = Nm/N0;
     312    //const Double_t Rdead  = 1.0 - kappa;
     313    //const Double_t dRdead = kappa * TMath::Sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
     314
     315    // the effective on time is Nm/lambda
     316    res[0] = teff;
     317    res[1] = dteff;
     318
     319    // plot chi2-probability of fit
     320    res[2] = prob*100;
     321
     322    // plot chi2/NDF of fit
     323    res[3] = NDF ? chi2/NDF : 0.0;
     324
     325    // lambda of fit
     326    res[4] = lambda;
     327    res[5] = dl;
     328
     329    // N0 of fit
     330    res[6] = N0;
     331
     332    // Rdead (from fit) is the fraction from real time lost by the dead time
     333    //fHRdead.SetBinContent(i, Rdead);
     334    //fHRdead.SetBinError  (i,dRdead);
     335
     336    return kTRUE;
     337}
    299338
    300339void MHEffectiveOnTime::Paint(Option_t *opt)
    … …  
    320359
    321360        if (!fIsFinalized)
    322             Fit();
     361            FitThetaBins();
    323362    }
    324363    if (o==(TString)"paint")
    … …  
    386425void MHEffectiveOnTime::Calc()
    387426{
    388     const Double_t val = fHEffOn.Integral();
    389 
    390     Double_t error = 0;
    391     for (int i=0; i<fHEffOn.GetXaxis()->GetNbins(); i++)
    392         error += fHEffOn.GetBinError(i);
     427    TH1D *h = fHTimeDiff.ProjectionX("", -1, 99999, "E");
     428    h->SetDirectory(0);
     429
     430    Double_t res[7];
     431    Bool_t rc = FitH(h, res);
     432    delete h;
     433
     434    if (!rc)
     435        return;
     436
     437    const Double_t val   = res[0];
     438    const Double_t error = res[1];
    393439
    394440    fParam->SetVal(val-fEffOnTime0, error-fEffOnErr0);
    … …  
    403449    MTime now(*fTime);
    404450    now.AddMilliSeconds(fInterval*1000);
     451
    405452    *fLog <<all << now << " - ";// << fLastTime-fTime;
    406453    *fLog << Form("T_{eff} = %.1fs \\pm %.1fs",
    … …  
    435482
    436483        *fTime = *time;
    437         // fTime->MinusNull()
     484
     485        // Make this just a ns before the first event
     486        fTime->Minus1ns();
    438487    }
    439488
    440489    if (*fTime+fInterval<*time)
    441490    {
    442         Fit();
     491        FitThetaBins();
    443492        Calc();
    444493    }

  • trunk/MagicSoft/Mars/mhist/MHEffectiveOnTime.h

    r4887 r4889  
    2424{
    2525private:
    26     MPointingPos *fPointPos;   //!
    27     MTime         fLastTime;   //!
     26    MPointingPos   *fPointPos;   //!
     27    MTime           fLastTime;   //!
    2828
    29     Double_t      fEffOnTime0; //!
    30     Double_t      fEffOnErr0;  //!
     29    Double_t        fEffOnTime0; //!
     30    Double_t        fEffOnErr0;  //!
    3131
    32     MTime          *fTime;
    33     MParameterDerr *fParam;
     32    MTime          *fTime;       //!
     33    MParameterDerr *fParam;      //!
    3434
    3535    TH2D fHTimeDiff;
    … …  
    4444    Int_t fInterval;
    4545
    46     void Fit();
     46    Bool_t FitH(TH1D *h, Double_t *res) const;
     47    void FitThetaBins();
    4748    void Calc();
    4849

  • trunk/MagicSoft/Mars/mhvstime/MHPixVsTime.cc

    r3539 r4889  
    194194    {
    195195        TAxis *axe = h->GetXaxis();
    196         axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00");
     196        axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
    197197        axe->SetTimeDisplay(1);
    198198        axe->SetLabelSize(0.033);

  • trunk/MagicSoft/Mars/mhvstime/MHSectorVsTime.cc

    r3394 r4889  
    247247    {
    248248        TAxis *axe = h->GetXaxis();
    249         axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00");
     249        axe->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
    250250        axe->SetTimeDisplay(1);
    251251        axe->SetLabelSize(0.025);

  • trunk/MagicSoft/Mars/mhvstime/MHVsTime.cc

    r4828 r4889  
    4949#include <TCanvas.h>
    5050
    51 #include <TGraph.h>
     51#include <TGraphErrors.h>
    5252
    5353#include "MLog.h"
    … …  
    7171// see MDataChain.
    7272//
    73 MHVsTime::MHVsTime(const char *rule)
    74     : fGraph(NULL), fData(NULL), fScale(1), fMaxPts(-1), fUseEventNumber(0)
    75 { 
     73MHVsTime::MHVsTime(const char *rule, const char *error)
     74    : fGraph(NULL), fData(NULL), fError(0), fScale(1), fMaxPts(-1),
     75    fNumEvents(1), fUseEventNumber(0)
     76{
    7677    fName  = gsDefName;
    7778    fTitle = gsDefTitle;
    … …  
    8081        return;
    8182
    82     fGraph = new TGraph;
    8383    fData = new MDataChain(rule);
    8484
    85     fGraph->SetMarkerStyle(kFullDotMedium);
     85    if (error)
     86        fError = new MDataChain(error);
     87
     88    fGraph = error ? new TGraphErrors : new TGraph;
     89    fGraph->SetPoint(0, 0, 0); // Dummy point!
     90    fGraph->SetEditable();     // Used as flag: First point? yes/no
    8691}
    8792
    … …  
    116121Bool_t MHVsTime::SetupFill(const MParList *plist)
    117122{
    118     // reset histogram (necessary if the same eventloop is run more than once)
    119     //fGraph->Reset();
    120 
    121     if (fData && !fData->PreProcess(plist))
     123    if (!fGraph || !fData)
     124    {
     125        *fLog << err << "ERROR - MHVsTime cannot be used with its default constructor!" << endl;
    122126        return kFALSE;
    123 
    124     if (fGraph)
    125     {
    126         delete fGraph;
    127         fGraph = new TGraph;
    128     }
     127    }
     128
     129    if (!fData->PreProcess(plist))
     130        return kFALSE;
     131
     132    if (fError && !fError->PreProcess(plist))
     133        return kFALSE;
     134
     135    fGraph->Set(1);
     136    fGraph->SetPoint(0, 0, 0); // Dummy point!
     137    fGraph->SetEditable();     // Used as flag: First point? yes/no
    129138
    130139    TString title(fData ? GetRule() : (TString)"Histogram");
    … …  
    183192        if (!*tm)
    184193            return kTRUE;
     194
     195        // Do not fill events with equal time
     196        if (*tm==fLast || *tm==MTime())
     197            return kTRUE;
     198
     199        fLast = *tm;
     200
    185201        t = tm->GetAxisTime();
    186202    }
    187203
    188     const Double_t v = fData->GetValue()*fScale;
    189 
    190     if (fMaxPts>0 && fGraph->GetN()>fMaxPts)
    191         fGraph->RemovePoint(0);
    192 
    193     fMean += v;
     204    const Double_t v = fData->GetValue();
     205    const Double_t e = fError ? fError->GetValue() : 0;
     206
     207    //*fLog << all << "ADD " << v << " " << e << endl;
     208
     209    fMean    += v;
     210    fMeanErr += e;
    194211    fN++;
    195212
    196213    if (fN==fNumEvents)
    197214    {
    198         fGraph->SetPoint(fGraph->GetN(), t, fMean/fN);
     215        if (fMaxPts>0 && fGraph->GetN()>fMaxPts || fGraph->IsEditable())
     216        {
     217            fGraph->RemovePoint(0);
     218            fGraph->SetEditable(kFALSE);
     219        }
     220
     221        fGraph->SetPoint(fGraph->GetN(), t, fMean/fN*fScale);
     222
     223        if (fError)
     224            static_cast<TGraphErrors*>(fGraph)->SetPointError(fGraph->GetN()-1, 0, fMeanErr/fN*fScale);
     225
    199226        fMean = 0;
     227        fMeanErr = 0;
    200228        fN = 0;
    201229    }
    202230
    203231    return kTRUE;
     232}
     233
     234void MHVsTime::Paint(Option_t *opt)
     235{
     236    // SetPoint deletes the histogram!
     237    if (fUseEventNumber)
     238        fGraph->GetHistogram()->SetXTitle("Event Number");
     239    else
     240    {
     241        fGraph->GetHistogram()->SetXTitle("Time");
     242        fGraph->GetHistogram()->GetXaxis()->SetLabelSize(0.033);
     243        fGraph->GetHistogram()->GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00 GMT");
     244        fGraph->GetHistogram()->GetXaxis()->SetTimeDisplay(1);
     245    }
     246
     247    fGraph->GetHistogram()->SetMarkerStyle(kFullDotMedium);
     248    fGraph->GetHistogram()->SetYTitle(fAxisTitle.IsNull() ? GetRule() : fAxisTitle);
     249    if (fTitle!=gsDefTitle)
     250        fGraph->GetHistogram()->SetTitle(fTitle);
     251
     252    if (TestBit(kIsLogy))
     253        gPad->SetLogy();
     254
     255    // This is a workaround if the TGraph has only one point.
     256    // Otherwise MStatusDisplay::Update hangs.
     257    gPad->GetListOfPrimitives()->Remove(fGraph);
     258    gPad->GetListOfPrimitives()->Add(fGraph, fGraph->GetN()<2 ? "A" : opt);
    204259}
    205260
    … …  
    212267void MHVsTime::Draw(Option_t *opt)
    213268{
     269    if (!fGraph)
     270        return;
     271
    214272    if (fGraph->GetN()==0)
    215273        return;
    … …  
    218276    pad->SetBorderMode(0);
    219277
    220     AppendPad("");
    221 
    222278    TString str(opt);
    223 
    224     if (fUseEventNumber)
    225         fGraph->GetHistogram()->SetXTitle("Event Number");
    226     else
    227     {
    228         fGraph->GetHistogram()->SetXTitle("Time");
    229         fGraph->GetHistogram()->GetXaxis()->SetTimeFormat("%H:%M:%S %F1995-01-01 00:00:00");
    230         fGraph->GetHistogram()->GetXaxis()->SetTimeDisplay(1);
    231         fGraph->GetHistogram()->GetXaxis()->SetLabelSize(0.033);
    232     }
    233     fGraph->GetHistogram()->SetYTitle(GetRule());
    234279
    235280    if (!str.Contains("A"))
    … …  
    244289    }
    245290
     291    AppendPad(str);
    246292    fGraph->Draw(str);
    247     if (fGraph->TestBit(kIsLogy))
    248         pad->SetLogy();
    249 
    250     pad->Modified();
    251     pad->Update();
    252293}
    253294

  • trunk/MagicSoft/Mars/mhvstime/MHVsTime.h

    r4828 r4889  
    44#ifndef MARS_MH
    55#include "MH.h"
     6#endif
     7
     8#ifndef MARS_MTime
     9#include "MTime.h"
    610#endif
    711
    … …  
    1519    TGraph     *fGraph;     // Histogram to fill
    1620    MDataChain *fData;      // Object from which the data is filled
     21    MDataChain *fError;     // Object from which the error is filled
    1722    Double_t    fScale;     // Scale for axis (eg unit)
    1823    Int_t       fMaxPts;    // Maximum number of data points
    1924
    2025    Int_t       fNumEvents; // Number of events to average
     26
    2127    Double_t    fMean;      //! Mean value
    22     Int_t       fN;
     28    Double_t    fMeanErr;   //! Mean error
     29    Int_t       fN;         //! Number of entries in fMean
     30    MTime       fLast;      //! For checks
    2331
     32    TString     fAxisTitle;
    2433
    2534    enum {
    26         kIsLogy         = BIT(18),
    27         kUseEventNumber = BIT(20)
     35        kIsLogy = BIT(18)
    2836    };
    2937
    … …  
    3139
    3240public:
    33     MHVsTime(const char *rule=NULL);
     41    MHVsTime(const char *rule=NULL, const char *ruleerr=NULL);
    3442    ~MHVsTime();
    3543
    … …  
    4048    void SetName(const char *name);
    4149    void SetTitle(const char *title);
     50
     51    void SetLogy(Bool_t b=kTRUE) { b ? SetBit(kIsLogy) : ResetBit(kIsLogy); }
     52    void SetAxisTitle(const char *y) { fAxisTitle=y; }
    4253
    4354    Bool_t SetupFill(const MParList *pList);
    … …  
    5970
    6071    void Draw(Option_t *opt=NULL);
     72    void Paint(Option_t *opt=NULL);
    6173
    6274    MParContainer *New() const;

  • trunk/MagicSoft/Mars/mjobs/MJStar.cc

    r4760 r4889  
    4343#include "MStatusDisplay.h"
    4444
     45#include "MH3.h"
     46#include "MHVsTime.h"
    4547#include "MHCamEvent.h"
    46 
    47 #include "MReadMarsFile.h"
     48#include "MBinning.h"
     49
     50#include "MReadReports.h"
    4851#include "MGeomApply.h"
     52#include "MEventRateCalc.h"
    4953#include "MImgCleanStd.h"
    5054#include "MHillasCalc.h"
    5155#include "MFillH.h"
    5256#include "MWriteRootFile.h"
     57
     58#include "MPointingPosCalc.h"
     59//#include "MSrcPosFromModel.h"
    5360
    5461ClassImp(MJStar);
    … …  
    138145    plist.AddToList(&tlist);
    139146
    140     MReadMarsFile read("Events");
    141     read.DisableAutoScheme();
     147    MReadReports read;
     148    read.AddTree("Events", "MTime.", kTRUE);
     149    read.AddTree("Drive");
     150    //read.AddTree("Trigger");
     151    //read.AddTree("Camera");
     152    //read.AddTree("CC");
     153    //read.AddTree("Currents");
    142154    read.AddFiles(iter);
    143     //read.AddFiles(fnamein);
     155
     156    // ------------------ Setup general tasks ----------------
    144157
    145158    MGeomApply             apply; // Only necessary to craete geometry
     159    MEventRateCalc         rate;
     160/*
     161    MEventRateCalc         rate10000;
     162    rate10000.SetNameEventRate("MEventRate10000");
     163    rate10000.SetNumEvents(10000);
     164 */
    146165    //MBadPixelsMerge        merge(&badpix);
    147166    MImgCleanStd           clean;
    148167    MHillasCalc            hcalc;
    149168
     169    // ------------------ Setup histograms and fill tasks ----------------
    150170    MHCamEvent evt0("Cleaned");
    151171    evt0.SetType(0);
     172
     173    MH3 h1("MEventRate.fRate");
     174    h1.SetName("MHEventRate");
     175    h1.SetLogy();
     176/*
     177    MH3 h12("MEventRate10000.fRate");
     178    h12.SetName("MHEventRate");
     179    h12.SetLogy();
     180 */
     181    MBinning b1("BinningMHEventRate");
     182    b1.SetEdges(150, 0, 1500);
     183    plist.AddToList(&b1);
     184
     185    MHVsTime h2("MEffectiveOnTime.fVal", "MEffectiveOnTime.fErr");
     186    h2.SetAxisTitle("T_{eff}");
     187    h2.SetTitle("Effective On-Time T_{eff} vs. Time");
     188
    152189    MFillH fill0(&evt0, "MCerPhotEvt",            "FillCerPhotEvt");
    153190    MFillH fill1("MHHillas",      "MHillas",      "FillHillas");
    154     MFillH fill2("MHHillasExt",   "MHillasExt",   "FillHillasExt");
     191    MFillH fill2("MHHillasExt",   "",             "FillHillasExt");
    155192    MFillH fill3("MHHillasSrc",   "MHillasSrc",   "FillHillasSrc");
    156193    MFillH fill4("MHImagePar",    "MImagePar",    "FillImagePar");
    157194    MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
    158     MFillH fill6("MHCerPhot");
    159 
    160     MWriteRootFile write(2, "images/{s/_Y_/_I_}");
    161     write.AddContainer("MMcEvt",        "Events", kFALSE);
     195    MFillH fill6("MHImageParTime","MImageParTime","FillImageParTime");
     196    MFillH fill7("MHNewImagePar2","MNewImagePar2","FillNewImagePar2");
     197    MFillH fill8(&h1,             "",             "FillEventRate");
     198    MFillH fill9("MHEffectiveOnTime", "MTime",    "FillEffOnTime");
     199    MFillH filla(&h2,             "MTimeEffectiveOnTime", "FillEffOnTimeVsTime");
     200    //MFillH fillb(&h12, "", "FillEvtRate2");
     201    //MFillH fill9("MHCerPhot");
     202
     203    fill8.SetNameTab("EvtRate");
     204    fill9.SetNameTab("EffOnTime");
     205    fill9.SetNameTab("EffOnVsTime");
     206
     207    // ------------------ Setup write task ----------------
     208
     209    MWriteRootFile write(2, Form("%s{s/_Y_/_I_}", fPathOut.Data()), fOverwrite);
     210    // Data
    162211    write.AddContainer("MHillas",       "Events");
    163212    write.AddContainer("MHillasExt",    "Events");
    … …  
    165214    write.AddContainer("MImagePar",     "Events");
    166215    write.AddContainer("MNewImagePar",  "Events");
     216    write.AddContainer("MNewImagePar2", "Events");
     217    write.AddContainer("MImageParTime", "Events");
    167218    write.AddContainer("MTime",         "Events");
    168219    write.AddContainer("MRawEvtHeader", "Events");
    169     write.AddContainer("MRawRunHeader", "RunHeaders");
    170     write.AddContainer("MBadPixelsCam", "RunHeaders");
    171     write.AddContainer("MGeomCam",      "RunHeaders");
     220    // Monte Carlo
     221    write.AddContainer("MMcEvt",              "Events", kFALSE);
     222    write.AddContainer("MMcTrig",             "Events", kFALSE);
     223    // Run Header
     224    write.AddContainer("MRawRunHeader",       "RunHeaders");
     225    write.AddContainer("MBadPixelsCam",       "RunHeaders");
     226    write.AddContainer("MGeomCam",            "RunHeaders");
    172227    //write.AddContainer("MObservatory", "RunHeaders");
     228    // Monte Carlo Headers
     229    write.AddContainer("MMcTrigHeader",       "RunHeaders", kFALSE);
     230    write.AddContainer("MMcConfigRunHeader",  "RunHeaders", kFALSE);
     231    write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
     232    // Drive
     233    //write.AddContainer("MSrcPosCam",   "Drive");
     234    write.AddContainer("MPointingPos", "Drive");
     235    write.AddContainer("MReportDrive", "Drive");
     236    write.AddContainer("MTimeDrive",   "Drive");
     237    // Effective On Time
     238    write.AddContainer("MEffectiveOnTime",     "EffectiveOnTime");
     239    write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime");
     240
     241    MTaskList tlist2;
     242    tlist2.AddToList(&apply);
     243    tlist2.AddToList(&rate);
     244    //tlist2.AddToList(&rate10000);
     245    tlist2.AddToList(&fill8);
     246    tlist2.AddToList(&fill9);
     247    tlist2.AddToList(&filla);
     248    //tlist2.AddToList(&fillb);
     249    tlist2.AddToList(&clean);
     250    tlist2.AddToList(&fill0);
     251    tlist2.AddToList(&hcalc);
     252    tlist2.AddToList(&fill1);
     253    tlist2.AddToList(&fill2);
     254    tlist2.AddToList(&fill3);
     255    tlist2.AddToList(&fill4);
     256    tlist2.AddToList(&fill5);
     257    tlist2.AddToList(&fill6);
     258    tlist2.AddToList(&fill7);
     259    //tlist2.AddToList(&fill9);
     260
     261    MPointingPosCalc pcalc;
     262    //MSrcPosFromModel srcpos;
     263
     264    MTaskList tlist3;
     265    tlist3.AddToList(&pcalc);
     266    //tlist3.AddToList(&srcpos);
    173267
    174268    tlist.AddToList(&read);
    175     tlist.AddToList(&apply);
    176     tlist.AddToList(&clean);
    177     tlist.AddToList(&fill0);
    178     tlist.AddToList(&hcalc);
    179     tlist.AddToList(&fill1);
    180     tlist.AddToList(&fill2);
    181     tlist.AddToList(&fill3);
    182     tlist.AddToList(&fill4);
    183     tlist.AddToList(&fill5);
    184     //tlist.AddToList(&fill6);
     269    tlist.AddToList(&tlist3, "Drive");
     270    tlist.AddToList(&tlist2, "Events");
    185271    tlist.AddToList(&write);
    186272
Note: See TracChangeset for help on using the changeset viewer.