Changeset 1668


Ignore:
Timestamp:
11/25/02 09:12:47 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
65 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r1667 r1668  
    11                                                               -*-*- END -*-*-
     2 2002/11/25: Thomas Bretz
     3
     4   * macros/multidimdist2.C:
     5     - renamed eventloops (instances had same names)
     6     - fixed a type in PrintStatistics (the gamma statistics
     7       were printed two times)
     8
     9   * mbase/MEvtLoop.cc:
     10     - take the lowest value (entries from MRead or user input)
     11       for the progress bar
     12     - reset the progress bar
     13
     14   * mbase/MFilter.h:
     15     - added 'private'
     16
     17   * meventdisp/MGCamDisplay.[h,cc], meventdisp/MGEvtDisplay.[h,cc],
     18     meventdisp/MGFadcDisp.[h,cc], mmain/MMonteCarlo.[h,cc],
     19     mmain/MAnalysis.[h,cc], mmain/MBrowser.[h,cc],
     20     mmain/MCameraDisplay.[h,cc], mmain/MDataCheck.[h,cc],
     21     mmain/MEvtDisp.[h,cc], mmain/MMars.cc:
     22     - changed from TTransientFrame to TMainFrame (with this I
     23       get decorations, eg. Close Button)
     24
     25   * meventdisp/MGEvtDisplay.cc:
     26     - Update the layout each time the fEvtInfo has changed
     27
     28   * mfileio/MCT1ReadAscii.cc, mfileio/MCT1ReadPreProc.cc:
     29     - delete return of gSystem->ExpandPathName
     30
     31   * mfileio/MCT1ReadPreProc.[h,cc]:
     32     - added output of Time
     33     - added usage of Selector
     34     - changed MTask basics to be private
     35
     36   * mfileio/MRead.[h,cc]:
     37     - added comment about selector
     38     - added Selector-stuff
     39
     40   * mfileio/MReadMarsFile.[h,cc], mfileio/MReadTree.[h,cc]:
     41     - added 'entries' argument to AddFile
     42
     43   * mfileio/MReadTree.[h,cc]:
     44     - added workaround for a root bug when a file doesn't exist
     45     - changed AddFiles to use Add(TChain*)
     46     - changed to use Selector
     47
     48   * mfilter/MF.cc:
     49     - Set debug level to suppress output when MFDataChain is created
     50
     51   * mfilter/MFEventSelector.h:
     52     - changed Pre//PostProcess to private
     53     
     54   * mfilter/MF.cc, mfilter/MFilterList.cc:
     55     - changed the use of Pre//PostProcess to CallPre//PostProcess
     56   
     57   * mhist/MBinning.[h,cc]:
     58     - changed comments
     59     - added SetEdgesCos
     60
     61   * mhist/MFillH.[h,cc]:
     62     - added GetBinCenterLog
     63
     64   * mhist/MH3.h:
     65     - added default argument to GetHistByName
     66
     67   * mhist/MHAlphaEnergyTheta.[h,cc], mhist/MHAlphaEnergyTime.h,
     68     mhist/MHEffOnTime.[h,cc], mhist/MHEffOnTimeTheta.h,
     69     mhist/MHEffOnTimeTime.h, mhist/MHFlux.[h,cc], mhist/MHGamma.[h,cc],
     70     mhist/MHMcEnergyMigration.h, mhist/MHThetabarTheta.[h,cc],
     71     mhist/MHThetabarTime.h:
     72     - changed the output
     73     - changed the algorithms to be more modular (more usage of member
     74       function)
     75     - changed ClassDef to 0
     76     - fixed some small bugs (access of TArray[n])
     77
     78   * mhist/MHHadronness.[h,cc]:
     79     - removed shortest distance to (0,1) stuff
     80
     81   * mhist/MHMcCollectionArea.h:
     82     - changed Fill to Double_t
     83
     84   * mhist/MHTimeDiffTheta.[h,cc], mhist/MHTimeDiffTime.[h,cc]:
     85     - in a first draft changed to use 200ns timing of CT1
     86     - changed ClassDef to 0
     87
     88
    289
    390 2002/11/22: Thomas Bretz
  • trunk/MagicSoft/Mars/mbase/MEvtLoop.cc

    r1657 r1668  
    208208    if (fProgress)
    209209    {
     210        fProgress->Reset();
     211
     212        Int_t entries = INT_MAX;
     213
     214#ifdef __MARS__
     215        // limits.h
     216        MRead *read = (MRead*)fTaskList->FindObject("MRead");
     217        if (read && read->GetEntries()>0)
     218            entries = read->GetEntries();
     219#endif
     220
    210221        if (maxcnt>0)
    211             fProgress->SetRange(0, maxcnt);
    212 #ifdef __MARS__
     222            fProgress->SetRange(0, TMath::Min(maxcnt, entries));
    213223        else
    214         {
    215             MRead *read = (MRead*)fTaskList->FindObject("MRead");
    216             if (read && read->GetEntries()>0)
    217                 fProgress->SetRange(0, read->GetEntries());
    218         }
    219 #endif
     224            if (entries!=INT_MAX)
     225                fProgress->SetRange(0, entries);
    220226    }
    221227
  • trunk/MagicSoft/Mars/mbase/MFilter.h

    r1481 r1668  
    1010class MFilter : public MTask
    1111{
     12private:
     13    virtual Bool_t PreProcess(MParList *pList);
     14    virtual Bool_t Process();
     15    virtual Bool_t PostProcess();
     16
    1217public:
    1318    MFilter(const char *name=NULL, const char *title=NULL);
    1419
    1520    virtual Bool_t IsExpressionTrue() const = 0;
    16 
    17     virtual Bool_t PreProcess(MParList *pList);
    18     virtual Bool_t Process();
    19     virtual Bool_t PostProcess();
    20 
    2121    virtual TString GetRule() const;
    2222
  • trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.cc

    r1540 r1668  
    148148//
    149149MGCamDisplay::MGCamDisplay(const char *filename,
    150                            const TGWindow *p, const TGWindow *main,
     150                           const TGWindow *p, /*const TGWindow *main,*/
    151151                           UInt_t w, UInt_t h)
    152 : MGEvtDisplay(filename, "Events", p, main, w, h), fDisplay(NULL)
     152: MGEvtDisplay(filename, "Events", p, /*main,*/ w, h), fDisplay(NULL)
    153153{
    154154    //
  • trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.h

    r1540 r1668  
    3232public:
    3333    MGCamDisplay(const char *filename,
    34                  const TGWindow *p, const TGWindow *main,
     34                 const TGWindow *p, /*const TGWindow *main,*/
    3535                 UInt_t w, UInt_t h);
    3636
  • trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.cc

    r1600 r1668  
    141141    fEvtInfo = new TGLabel(top2, new TGString(""));
    142142    fList->Add(fEvtInfo);
    143 
    144143    top2->AddFrame(fEvtInfo, laystd);
    145144
     
    393392
    394393MGEvtDisplay::MGEvtDisplay(const char *fname, const char *tname,
    395                            const TGWindow *p, const TGWindow *main,
     394                           const TGWindow *p, /*const TGWindow *main,*/
    396395                           UInt_t w, UInt_t h)
    397     : TGTransientFrame(p, main, w, h), fInitOk(kFALSE)
     396//    : TGTransientFrame(p, main, w, h), fInitOk(kFALSE)
     397: TGMainFrame(p, w, h), fInitOk(kFALSE)
    398398{
    399399    //
     
    472472        break;
    473473    default:
    474         txt += "Unknown Particle Id";
     474        txt += "Unknown Particle Id#";
     475        txt += evt->GetPartId();
    475476    }
    476477
     
    493494
    494495    fEvtInfo->SetText(txt);
     496
     497    //
     498    // Seems to be necessary to newly layout the upper part to display
     499    // the whole line of text
     500    //
     501    TGFrame &f = *(TGFrame*)fEvtInfo->GetParent()->GetParent();
     502    f.Layout();
     503    f.MapSubwindows();
    495504}
    496505
  • trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.h

    r1540 r1668  
    2121class MReadTree;
    2222
    23 class MGEvtDisplay : public TGTransientFrame
     23class MGEvtDisplay : public TGMainFrame/*TGTransientFrame*/
    2424{
    2525private:
     
    7474public:
    7575    MGEvtDisplay(const char *fname, const char *tname,
    76                  const TGWindow *p, const TGWindow *main,
     76                 const TGWindow *p, /*const TGWindow *main,*/
    7777                 UInt_t w, UInt_t h);
    7878
  • trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.cc

    r1172 r1668  
    113113//
    114114MGFadcDisp::MGFadcDisp(const char *filename, const char *treename,
    115                        const TGWindow *p, const TGWindow *main,
     115                       const TGWindow *p, /*const TGWindow *main,*/
    116116                       UInt_t w, UInt_t h)
    117     : MGEvtDisplay(filename, treename, p, main, w, h)
     117    : MGEvtDisplay(filename, treename, p, /*main,*/ w, h)
    118118{
    119119    //
  • trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.h

    r1172 r1668  
    2323
    2424    MGFadcDisp(const char *filename, const char *treename,
    25                const TGWindow *p, const TGWindow *main, UInt_t w, UInt_t h);
     25               const TGWindow *p, /*const TGWindow *main,*/ UInt_t w, UInt_t h);
    2626
    2727    virtual Bool_t ProcessMessage(Long_t msg, Long_t parm1, Long_t parm2);
  • trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc

    r1583 r1668  
    127127    const char *name = file->GetName();
    128128
    129     fIn = new ifstream(gSystem->ExpandPathName(name));
     129    const char *expname = gSystem->ExpandPathName(name);
     130    fIn = new ifstream(expname);
     131    delete [] expname;
    130132
    131133    const Bool_t noexist = !(*fIn);
  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

    r1666 r1668  
    1616!
    1717!
    18 !   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
    19 !   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
     18!   Author(s): Thomas Bretz 11/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
    2019!
    21 !   Copyright: MAGIC Software Development, 2000-2001
     20!   Copyright: MAGIC Software Development, 2000-2002
    2221!
    2322!
     
    2928//
    3029// Reads a output file of the CT1 preproc.
     30//
     31// Implements usage of a selector (see MRead) Use such a filter to skip
     32// events before reading! But never use a filter which needs read data
     33// as input...
    3134//
    3235//  Input Containers:
     
    5861#include "MLogManip.h"
    5962
     63#include "MTime.h"
     64#include "MFilter.h"
     65
    6066#include "MParList.h"
    6167#include "MCerPhotEvt.h"
     
    112118void MCT1ReadPreProc::AddFile(const char *txt)
    113119{
    114     const TString fname = gSystem->ExpandPathName(txt);
     120    const char *name = gSystem->ExpandPathName(txt);
     121    TString fname(name);
     122    delete [] name;
    115123
    116124    if (!CheckHeader(fname))
     
    502510    // open the file which is the first one in the chain
    503511    //
    504     const TString name  = file->GetName();
    505     const TString fname = gSystem->ExpandPathName(name);
     512    const TString name = file->GetName();
     513
     514    const char *expname = gSystem->ExpandPathName(name);
     515    const TString fname(expname);
     516    delete [] expname;
    506517
    507518    //
     
    648659
    649660    //
     661    //  look for the time class in the plist
     662    //
     663    fTime = (MTime*)pList->FindCreateObj("MTime");
     664    if (!fTime)
     665        return kFALSE;
     666
     667    //
    650668    //  look for the pedestal class in the plist
    651669    //
     
    685703    fNumRuns       = 0;
    686704
    687     return kTRUE;
     705    return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE;
    688706}
    689707
     
    695713void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event)
    696714{
     715    //  int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
     716    //  int   isecfrac_200ns;     // fractional part of isecs_since_midday
     717    fTime->SetTime(event.isecfrac_200ns, event.isecs_since_midday);
     718    fTime->SetReadyToSave();
     719
    697720    /*
    698721     --- USEFULL? NEEDED? ---
    699      int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
    700      int   isecfrac_200ns;     // fractional part of isecs_since_midday
    701722     short snot_ok_flags;      // the bits in these two bytes are flags for additional information on the event: Everything OK =: all Bits = 0
    702723
     
    862883            return kFALSE;
    863884
     885    //
     886    // Check for a selector. If one is given and returns kFALSE
     887    // skip this event.
     888    //
     889    if (GetSelector())
     890    {
     891        //
     892        // Make sure selector is processed
     893        //
     894        if (!GetSelector()->CallProcess())
     895        {
     896            *fLog << err << dbginf << "Processing Selector failed." << endl;
     897            return kFALSE;
     898        }
     899
     900        //
     901        // Skip Event
     902        //
     903        if (!GetSelector()->IsExpressionTrue())
     904        {
     905            fIn->seekg(sizeof(struct eventrecord), ios::cur);
     906
     907            fNumEvents++;
     908            fNumEventsInRun++;
     909
     910            return kCONTINUE;
     911        }
     912    }
     913
    864914    // event data to be read from the file
    865915    struct eventrecord event;
     
    890940    }
    891941
    892     return kTRUE;
    893 }
     942    return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE;
     943}
  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h

    r1664 r1668  
    77
    88class TList;
     9class MTime;
    910class MMcEvt;
    1011class MMcTrig;
     
    2728    MCerPhotEvt  *fNphot;   // the data container for all data.
    2829    MPedestalCam *fPedest;  // ct1 pedestals
     30    MTime        *fTime;    // event time
    2931    MMcEvt       *fMcEvt;   // monte carlo data container for MC files
    3032    MMcTrig      *fMcTrig;  // mc data container for trigger information
     
    5254    void   ProcessEvent(const struct eventrecord &event);
    5355
     56    Bool_t PreProcess(MParList *pList);
     57    Bool_t Process();
     58    Bool_t PostProcess();
     59
    5460public:
    5561    MCT1ReadPreProc(const char *filename=NULL,
     
    6167    void AddFile(const char *fname);
    6268
    63     Bool_t PreProcess(MParList *pList);
    64     Bool_t Process();
    65     Bool_t PostProcess();
    66 
    6769    UInt_t GetEntries() { return fEntries; }
    6870
  • trunk/MagicSoft/Mars/mfileio/MRead.cc

    r1600 r1668  
    2929// Base class for all reading tasks                                        //
    3030//                                                                         //
     31// You can set a selector. Depending on the impelementation in the derived //
     32// class it can be used to skip events, if the filter return kFALSE.       //
     33// Make sure that the selector (filter) doesn't need information which     //
     34// doesn't exist before reading an event!                                  //
     35//                                                                         //
    3136/////////////////////////////////////////////////////////////////////////////
    3237#include "MRead.h"
  • trunk/MagicSoft/Mars/mfileio/MRead.h

    r1600 r1668  
    66#endif
    77
     8class MFilter;
     9
    810class MRead : public MTask
    911{
     12private:
     13    MFilter *fSelector;
     14
    1015public:
     16    MRead() : fSelector(NULL) {}
    1117
    1218    virtual UInt_t GetEntries() = 0;
     19
     20    void SetSelector(MFilter *f) { fSelector = f; }
     21    MFilter *GetSelector() { return fSelector; }
    1322
    1423    ClassDef(MRead, 0)  // Base class for a reading task
  • trunk/MagicSoft/Mars/mfileio/MReadMarsFile.cc

    r1667 r1668  
    9898//  chains -1 is returned, otherwise the number of files which were added.
    9999//
    100 Int_t MReadMarsFile::AddFile(const char *fname)
     100Int_t MReadMarsFile::AddFile(const char *fname, Int_t entries)
    101101{
    102102    //
     
    107107    //
    108108    Int_t n1 = fRun->AddFile(fname);
    109     Int_t n2 = MReadTree::AddFile(fname);
     109    Int_t n2 = MReadTree::AddFile(fname, entries);
    110110
    111111    return n1 != n2 ? -1 : n1;
  • trunk/MagicSoft/Mars/mfileio/MReadMarsFile.h

    r1664 r1668  
    2424    ~MReadMarsFile();
    2525
    26     Int_t AddFile(const char *fname);
     26    Int_t AddFile(const char *fname, Int_t entries=-1);
    2727
    2828    ClassDef(MReadMarsFile, 1)  // Reads a tree from file(s) and the information from the 'RunHeader'-tree
  • trunk/MagicSoft/Mars/mfileio/MReadTree.cc

    r1664 r1668  
    191191//  AddFile returns the number of files added to the chain.
    192192//
    193 Int_t MReadTree::AddFile(const char *fname)
    194 {
     193//  For more information see TChain::Add
     194//
     195Int_t MReadTree::AddFile(const char *fname, Int_t entries)
     196{
     197#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,01)
     198    //
     199    // This is a workaround to get rid of crashed if the file doesn't
     200    // exist due to non initialized TFile::fProcessIDs
     201    //
     202    //  (Code taken from TFile::TFile
     203    //
     204    const char *name;
     205
     206    TString newname;
     207
     208    if ((name = gSystem->ExpandPathName(fname)))
     209    {
     210        newname = name;
     211        delete [] name;
     212    }
     213
     214    if (newname.IsNull())
     215    {
     216        *fLog << err << dbginf << "Error expanding path " << fname << "." << endl;
     217        return 0;
     218    }
     219
     220    if (gSystem->AccessPathName(newname, kFileExists))
     221    {
     222        *fLog << err << "ERROR - File '" << fname << "' does not exist." << endl;
     223        return 0;
     224    }
     225
     226    fname = newname.Data();
     227#endif
     228
    195229    //
    196230    // FIXME! A check is missing whether the file already exists or not.
    197231    //
    198     //
    199     // returns the number of file which were added
    200     //
    201 
    202     const Int_t numfiles = fChain->Add(fname);
     232    const Int_t numfiles = fChain->Add(fname, entries);
    203233
    204234    if (numfiles>0)
     
    208238}
    209239
     240/*
     241 // --------------------------------------------------------------------------
     242 //
     243 //
     244 Int_t MReadTree::AddFile(const TChainElement &obj)
     245 {
     246     return AddFile(obj.GetTitle(), obj.GetEntries());
     247 }
     248*/
     249
    210250// --------------------------------------------------------------------------
    211251//
     
    216256Int_t MReadTree::AddFiles(const MReadTree &read)
    217257{
    218     Int_t rc = 0;
    219 
    220     TIter Next(read.fChain->GetListOfFiles());
    221     TObject *obj = NULL;
    222     while ((obj=Next()))
    223         rc += AddFile(obj->GetTitle());
     258    const Int_t rc = fChain->Add(read.fChain);
    224259
    225260    if (rc>0)
    226261        SetBit(kChainWasChanged);
     262
     263    /*
     264     Int_t rc = 0;
     265
     266     TIter Next(read.fChain->GetListOfFiles());
     267     TObject *obj = NULL;
     268     while ((obj=Next()))
     269         rc += AddFile(*(TChainElement*)obj);
     270    */
    227271
    228272    return rc;
     
    241285        return;
    242286
    243     *fLog << inf << GetDescriptor() << ": Branch choosing method enabled (only enabled branches are read)." << endl;
     287    *fLog << inf << GetDescriptor() << ": Branch choosing enabled (only enabled branches are read)." << endl;
    244288    fChain->SetBranchStatus("*", kFALSE);
    245289    fBranchChoosing = kTRUE;
     
    255299void MReadTree::EnableBranch(const char *name)
    256300{
     301    if (fChain->GetListOfFiles()->GetEntries()==0)
     302    {
     303        *fLog << err << "Chain contains no file... Branch '";
     304        *fLog << name << "' ignored." << endl;
     305        return;
     306    }
     307
    257308    EnableBranchChoosing();
    258309
     
    449500//  root environment) the branch is skipped and an error message is printed
    450501//  out.
     502//  If a selector is specified it is preprocessed after the
     503//  MReadTree::PreProcess
    451504//
    452505Bool_t MReadTree::PreProcess(MParList *pList)
     
    576629    fChain->SetNotify(this);
    577630
    578     return kTRUE;
     631    return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE;
    579632}
    580633
     
    626679//  (Remark: The position can also be set by some member functions
    627680//  If the end of the file is reached the Eventloop is stopped.
     681//  In case an event selector is given its value is checked before
     682//  reading the event. If it returns kAFLSE the event is skipped.
    628683//
    629684#if ROOT_VERSION_CODE < ROOT_VERSION(3,02,06)
     
    658713#endif
    659714
    660     Bool_t rc = fChain->GetEntry(fNumEntry++) != 0;
     715    if (GetSelector())
     716    {
     717        //
     718        // Make sure selector is processed
     719        //
     720        if (!GetSelector()->CallProcess())
     721        {
     722            *fLog << err << dbginf << "Processing Selector failed." << endl;
     723            return kFALSE;
     724        }
     725
     726        //
     727        // Skip Event
     728        //
     729        if (!GetSelector()->IsExpressionTrue())
     730        {
     731            fNumEntry++;
     732            return kCONTINUE;
     733        }
     734    }
     735
     736    const Bool_t rc = fChain->GetEntry(fNumEntry++) != 0;
    661737
    662738    if (rc)
     
    664740
    665741    return rc;
     742}
     743
     744// --------------------------------------------------------------------------
     745//
     746//  If a selector is given the selector is post processed
     747//
     748Bool_t MReadTree::PostProcess()
     749{
     750    return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE;
    666751}
    667752
  • trunk/MagicSoft/Mars/mfileio/MReadTree.h

    r1600 r1668  
    6464    virtual void   Print(Option_t *opt="") const;
    6565
    66     virtual Int_t  AddFile(const char *fname);
     66    virtual Int_t  AddFile(const char *fname, Int_t entries=-1);
    6767    virtual Int_t  AddFiles(const MReadTree &read);
    6868
    6969    virtual Bool_t PreProcess(MParList *pList);
    7070    virtual Bool_t Process();
     71    virtual Bool_t PostProcess();
    7172
    7273    virtual Bool_t Notify();
  • trunk/MagicSoft/Mars/mfilter/MF.cc

    r1666 r1668  
    244244    if (isrule)
    245245    {
     246        Int_t lvl = gLog.GetDebugLevel();
     247        gLog.SetDebugLevel(1);
    246248        newfilter = new MFDataChain(text.Data(), c, num);
    247249        newfilter->SetName(Form("Chain%02d%c%f", level, c, num));
     250        gLog.SetDebugLevel(lvl);
    248251    }
    249252    else
     
    407410    }
    408411
    409     if (!fFilter->PreProcess(plist))
     412    if (!fFilter->CallPreProcess(plist))
    410413    {
    411414        *fLog << err << dbginf << "PreProcessing filters in ";
     
    423426Bool_t MF::Process()
    424427{
    425     return fFilter->Process();
     428    return fFilter->CallProcess();
    426429}
    427430
     
    432435Bool_t MF::PostProcess()
    433436{
    434     return fFilter->PostProcess();
     437    return fFilter->CallPostProcess();
    435438}
    436439
  • trunk/MagicSoft/Mars/mfilter/MFDataChain.cc

    r1666 r1668  
    5858//
    5959MFDataChain::MFDataChain(const char *member, const char type, const Double_t val,
    60                            const char *name, const char *title)
     60                         const char *name, const char *title)
    6161    : fData(member), fValue(val)
    6262{
  • trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc

    r1600 r1668  
    147147Bool_t MFEventSelector::Process()
    148148{
    149     Float_t evt = gRandom->Uniform();
     149    const Float_t evt = gRandom->Uniform();
    150150
    151151    if (fNumSelectEvts>0)
  • trunk/MagicSoft/Mars/mfilter/MFEventSelector.h

    r1600 r1668  
    2727    void StreamPrimitive(ofstream &out) const;
    2828
     29    Bool_t PreProcess(MParList *pList);
     30    Bool_t Process();
     31    Bool_t PostProcess();
     32
    2933    enum { kNumTotalFromFile = BIT(14) };
    3034    /*
     
    4751     */
    4852
    49     Bool_t PreProcess(MParList *pList);
    50     Bool_t Process();
    51     Bool_t PostProcess();
    52 
    5353    ClassDef(MFEventSelector, 0) // A Filter to select events from files
    5454};
  • trunk/MagicSoft/Mars/mfilter/MFilterList.cc

    r1493 r1668  
    188188    //
    189189    while ((filter=(MFilter*)Next()))
    190         if (!filter->PreProcess(pList))
     190        if (!filter->CallPreProcess(pList))
    191191        {
    192192            *fLog << err << "Error - Preprocessing Filter ";
     
    212212    //
    213213    while ((filter=(MFilter*)Next()))
    214         if (!filter->Process())
     214        if (!filter->CallProcess())
    215215            return kFALSE;
    216216
     
    232232    //
    233233    while ((filter=(MFilter*)Next()))
    234         if (!filter->PostProcess())
     234        if (!filter->CallPostProcess())
    235235            return kFALSE;
    236236
  • trunk/MagicSoft/Mars/mhist/MBinning.cc

    r1487 r1668  
    6464// --------------------------------------------------------------------------
    6565//
    66 // Specify the number of bins (not the number of edges), the lowest and
    67 // highest Edge.
     66// Specify the number of bins <nbins> (not the number of edges), the
     67// lowest <lo> and highest <up> Edge (of your histogram)
    6868//
    6969void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up)
     
    7979// --------------------------------------------------------------------------
    8080//
    81 // Specify the number of bins (not the number of edges), the lowest and
    82 // highest Edge.
     81// Specify the number of bins <nbins> (not the number of edges), the
     82// lowest <lo> and highest <up> Edge (of your histogram)
    8383//
    8484void MBinning::SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up)
     
    9292
    9393    fType = kIsLogarithmic;
     94}
     95
     96// --------------------------------------------------------------------------
     97//
     98// Specify the number of bins <nbins> (not the number of edges), the
     99// lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
     100//
     101void MBinning::SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up)
     102{
     103    // if (lo==0) ...
     104    const Axis_t ld = lo/kRad2Deg;
     105    const Axis_t ud = up/kRad2Deg;
     106
     107    const Double_t binsize = (cos(ld)-cos(ud))/nbins;
     108    fEdges.Set(nbins+1);
     109    for (int i=0; i<=nbins; i++)
     110        fEdges[i] = (acos(binsize*i) + ld)*kRad2Deg;
     111
     112    fType = kIsCosinic;
    94113}
    95114
     
    132151        return;
    133152
    134     if (IsLinear() || IsLogarithmic())
     153    if (IsLinear() || IsLogarithmic() || IsCosinic())
    135154    {
    136155        out << "   " << GetUniqueName() << ".SetEdges";
    137156        if (IsLogarithmic())
    138157            out << "Log";
     158        if (IsCosinic())
     159            out << "Cos";
    139160        out << "(" << GetNumBins() << ", " << GetEdgeLo() << ", " << GetEdgeHi() << ");" << endl;
    140161        return;
  • trunk/MagicSoft/Mars/mhist/MBinning.h

    r1572 r1668  
    2525        kIsLinear,
    2626        kIsLogarithmic,
     27        kIsCosinic,
    2728        kIsUserArray
    2829    };
     
    3940    void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up);
    4041    void SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up);
     42    void SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up);
    4143
    4244    // FIXME: ROOT workaround: "operator[] const" missing
     
    5153    Bool_t IsLinear() const { return fType==kIsLinear; }
    5254    Bool_t IsLogarithmic() const { return fType==kIsLogarithmic; }
     55    Bool_t IsCosinic() const { return fType==kIsCosinic; }
    5356    Bool_t IsDefault() const { return fType==kIsDefault; }
    5457    Bool_t IsUserArray() const { return fType==kIsUserArray; }
  • trunk/MagicSoft/Mars/mhist/MFillH.cc

    r1629 r1668  
    354354        {
    355355            if (cls==name)
    356                 *fLog << inf << "Object '" << fHName << "' not found in parlist... trying to create one." << endl;
     356                *fLog << inf << "Object '" << fHName << "' not found in parlist... creating." << endl;
    357357            obj = pList->FindCreateObj(cls, name);
    358358        }
  • trunk/MagicSoft/Mars/mhist/MH.cc

    r1648 r1668  
    546546// --------------------------------------------------------------------------
    547547//
     548//  Returns the bin center in a logarithmic scale. If the given bin
     549//  number is <1 it is set to 1. If it is =GetNbins() it is set to
     550//  GetNbins()
     551//
     552Double_t MH::GetBinCenterLog(const TAxis &axe, Int_t nbin)
     553{
     554    if (nbin>axe.GetNbins())
     555        nbin = axe.GetNbins();
     556
     557    if (nbin<1)
     558        nbin = 1;
     559
     560    const Double_t lo = axe.GetBinLowEdge(nbin);
     561    const Double_t hi = axe.GetBinUpEdge(nbin);
     562
     563    const Double_t val = log10(lo) + log10(hi);
     564
     565    return pow(10, val/2);
     566}
     567
     568// --------------------------------------------------------------------------
     569//
    548570// Draws a copy of the two given histograms. Uses title as the pad title.
    549571// Also layout the two statistic boxes and a legend.
  • trunk/MagicSoft/Mars/mhist/MH.h

    r1574 r1668  
    5454    static void    ScaleAxis(TH1 *bins, Double_t fx=1, Double_t fy=1, Double_t fz=1);
    5555
     56    static Double_t GetBinCenterLog(const TAxis &axe, Int_t nbin);
     57
    5658    static void DrawCopy(const TH1 &hist1, const TH1 &hist2, const TString title);
    5759    static void Draw(TH1 &hist1, TH1 &hist2, const TString title);
  • trunk/MagicSoft/Mars/mhist/MH3.h

    r1574 r1668  
    5656    const TH1 &GetHist() const { return *fHist; }
    5757
    58     TH1 *GetHistByName(const TString name) { return fHist; }
     58    TH1 *GetHistByName(const TString name="") { return fHist; }
    5959
    6060    void SetColors() const;
  • trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc

    r1415 r1668  
    6464    fHist.SetDirectory(NULL);
    6565
    66     fHist.SetTitle("3D-plot of alpha, E-est, Theta");
     66    fHist.SetTitle("3D-plot of alpha, E_{est}, Theta");
    6767    fHist.SetXTitle("\\alpha [\\circ]");
    68     fHist.SetYTitle("E-est [GeV]            ");
     68    fHist.SetYTitle("E_{est} [GeV]");
    6969    fHist.SetZTitle("\\Theta [\\circ]");
    7070}
     
    7979   if (!fEnergy)
    8080   {
    81        *fLog << err << dbginf << "MHAlphaEnergyTheta : MEnergyEst not found... aborting." << endl;
     81       *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
    8282       return kFALSE;
    8383   }
     
    8686   if (!fMcEvt)
    8787   {
    88        *fLog << err << dbginf << "MHAlphaEnergyTheta : MMcEvt not found... aborting." << endl;
     88       *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
    8989       return kFALSE;
    9090   }
     
    9595   if (!binsenergy || !binsalphaflux || !binstheta)
    9696   {
    97        *fLog << err << dbginf << "MHAlphaEnergyTheta : At least one MBinning not found... aborting." << endl;
     97       *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
    9898       return kFALSE;     
    9999   }
     
    103103   fHist.Sumw2();
    104104
    105    fHist.Sumw2();
    106 
    107105   return kTRUE;
    108106}
     
    116114    MHillasSrc &hil = *(MHillasSrc*)par;
    117115
    118     fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg);
     116    fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(),
     117               fMcEvt->GetTelescopeTheta()*kRad2Deg);
    119118
    120119    return kTRUE;
     
    148147
    149148    h->SetTitle("Distribution of E-est [GeV]");
    150     h->SetXTitle("E-est [GeV]            ");
     149    h->SetXTitle("E_{est} [GeV]");
    151150    h->SetYTitle("Counts");
    152151
     
    186185
    187186    c.cd(1);
    188     h = ((TH3*)(&fHist))->Project3D("expro");
     187    h = ((TH3*)(&fHist))->Project3D(fName+"_expro");
    189188
    190189    h->SetTitle("Distribution of \\alpha [\\circ]");
     
    196195
    197196    c.cd(2);
    198     h = ((TH3*)(&fHist))->Project3D("eypro");
     197    h = ((TH3*)(&fHist))->Project3D(fName+"_eypro");
    199198
    200199    h->SetTitle("Distribution of E-est [GeV]");
    201     h->SetXTitle("E-est [GeV]            ");
     200    h->SetXTitle("E_{est} [GeV]");
    202201    h->SetYTitle("Counts");
    203202
     
    207206
    208207    c.cd(3);
    209     h = ((TH3*)(&fHist))->Project3D("ezpro");
     208    h = ((TH3*)(&fHist))->Project3D(fName+"_ezpro");
    210209
    211210    h->SetTitle("Distribution of \\Theta [\\circ]");
  • trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h

    r1574 r1668  
    4444    TObject *DrawClone(Option_t *option="") const;
    4545
    46     ClassDef(MHAlphaEnergyTheta, 1) //3D-histogram in alpha, Energy and theta
     46    ClassDef(MHAlphaEnergyTheta, 0) //3D-histogram in alpha, Energy and theta
    4747};
    4848
  • trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h

    r1574 r1668  
    4747    TH1D *IntegrateEestTime(const char *title, Bool_t Draw);
    4848   
    49     ClassDef(MHAlphaEnergyTime, 1) //3D-histogram in alpha, Energy and time
     49    ClassDef(MHAlphaEnergyTime, 0) //3D-histogram in alpha, Energy and time
    5050};
    5151
  • trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc

    r1411 r1668  
    5454#include "MLogManip.h"
    5555
     56#include "MHTimeDiffTheta.h"
     57
    5658ClassImp(MHEffOnTime);
    5759
     
    6870    fUnit     = unit;
    6971
    70     TString strg(fVarname);
    71     strg += "  ";
    72     strg += fUnit;
     72    TString strg = fVarname + " " + fUnit;
    7373
    7474    //
    7575    //   set the name and title of this object
    7676    //
    77     TString strg1("MHEffOnTime-");
    78     strg1 += fVarname;
    79     fName  = strg1;
    80     fTitle =  "1-D histogram of Eff On Time";
     77    fName  = TString("MHEffOnTime-")+fVarname;
     78    fTitle = "1-D histogram of Eff On Time";
    8179
    8280    // effective on time versus Var
    8381    fHEffOn.SetName("EffOn");
    84     TString strg2("effective On Time vs. ");
    85     strg2 += fVarname;
    86     fHEffOn.SetTitle(strg2);
     82    fHEffOn.SetTitle(TString("T_{on, eff} vs. ")+fVarname);
    8783
    8884    fHEffOn.SetDirectory(NULL);
    8985
    9086    fHEffOn.SetXTitle(strg);
    91     fHEffOn.SetYTitle("effective ontime [s]");
     87    fHEffOn.SetYTitle("T_{on, eff} [s]");
    9288
    9389    // chi2-probability versus Var
    9490    fHProb.SetName("Chi2-prob");
    95     TString strg3("Chi2-prob of OnTimeFit vs. ");
     91    TString strg3("\\chi^{2}-prob of OnTimeFit vs. ");
    9692    strg3 += fVarname;
    9793    fHProb.SetTitle(strg3);
     
    10096
    10197    fHProb.SetXTitle(strg);
    102     fHProb.SetYTitle("chi2-probability");
     98    fHProb.SetYTitle("\\chi^{2}-probability");
    10399
    104100    // lambda versus Var
    105101    fHLambda.SetName("lambda");
    106     TString strg4("lambda from OnTimeFit vs. ");
    107     strg4 += fVarname;
    108     fHLambda.SetTitle(strg4);
     102    fHLambda.SetTitle(TString("\\lambda from OnTimeFit vs. ")+fVarname);
    109103
    110104    fHLambda.SetDirectory(NULL);
     
    116110    // Rdead is the fraction from real time lost by the dead time
    117111    fHRdead.SetName("Rdead");
    118     TString strg5("Rdead of  OnTimeFit vs. ");
    119     strg5 += fVarname;
    120     fHRdead.SetTitle(strg5);
     112
     113    fHRdead.SetTitle(TString("Rdead of OnTimeFit vs. ")+fVarname);
    121114
    122115    fHRdead.SetDirectory(NULL);
     
    127120}
    128121
     122// --------------------------------------------------------------------
     123//
     124//  determine range (yq[0], yq[1]) of time differences
     125//  where fit should be performed;
     126//  require a fraction >=min of all entries to lie below yq[0]
     127//      and a fraction <=max of all entries to lie below yq[1];
     128//  within the range (yq[0], yq[1]) there must be no empty bin;
     129//
     130void MHEffOnTime::GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const
     131{
     132#if ROOT_VERSION_CODE < ROOT_VERSION(3,02,07)
     133    // WOrkaround for missing const qualifier of TH1::Integral
     134    TH1 &h1 = (TH1&)h;
     135
     136    // choose pedestrian approach as long as GetQuantiles is
     137    // not available
     138    const Int_t    jbins = h1.GetNbinsX();
     139    const Double_t Nm    = h1.Integral();
     140
     141    const Double_t xq[2] = { 0.15*Nm, 0.98*Nm };
     142
     143    yq[0] = yq[1] = h1.GetBinLowEdge(jbins+1);
     144
     145    for (int j=1; j<=jbins; j++)
     146        if (h1.Integral(2, j) >= xq[0])
     147        {
     148            yq[0] = h1.GetBinLowEdge(j);
     149            break;
     150        }
     151
     152    for (int j=1; j<=jbins; j++)
     153        if (h1.Integral(1, j) >= xq[1] || h1.GetBinContent(j)==0)
     154        {
     155            yq[1] = h1.GetBinLowEdge(j);
     156            break;
     157        }
     158#else
     159    // GetQuantiles doesn't seem to be available in root 3.01/06
     160    Double_t xq[2] = { min, max };
     161    h.GetQuantiles(2, yq, xq);
     162#endif
     163}
     164
     165void MHEffOnTime::DrawBin(TH1 &h, Int_t i) const
     166{
     167    TString strg1 = fVarname+"-bin #";
     168    strg1 += i;
     169
     170    new TCanvas(strg1, strg1);
     171
     172    gPad->SetLogy();
     173    gStyle->SetOptFit(1011);
     174
     175    TString name="Bin_";
     176    name += i;
     177
     178    h.SetName(name);
     179    h.SetXTitle("\\Delta t [s]");
     180    h.SetYTitle("Counts");
     181    h.DrawCopy();
     182}
     183
     184Bool_t MHEffOnTime::CalcResults(const TF1 &func, Double_t Nm, Int_t i)
     185{
     186    const Double_t lambda = func.GetParameter(0);
     187    const Double_t N0     = func.GetParameter(1);
     188    const Double_t prob   = func.GetProb();
     189    const Int_t    NDF    = func.GetNDF();
     190
     191    Double_t xmin, xmax;
     192    ((TF1&)func).GetRange(xmin, xmax);
     193
     194    *fLog << inf;
     195    *fLog << "Fitted bin #" << i << " from " << xmin << " to " << xmax;
     196    *fLog << ",  got: lambda=" << lambda << "Hz N0=" << N0 << endl;
     197
     198    if (prob<=0.01)
     199    {
     200        *fLog << warn << "WARN - Fit bin#" << i << " gives:";
     201        *fLog << " Chi^2-Probab(" << prob << ")<0.01";
     202        *fLog << " NoFitPts=" << func.GetNumberFitPoints();
     203        *fLog << " Chi^2=" << func.GetChisquare();
     204    }
     205
     206    // was fit successful ?
     207    if (NDF<=0 || /*prob<=0.001 ||*/ lambda<=0 || N0<=0)
     208    {
     209        *fLog << warn << dbginf << "Fit failed bin #" << i << ": ";
     210        if (NDF<=0)
     211            *fLog << " NDF(Number of Degrees of Freedom)=0";
     212        if (lambda<=0)
     213            *fLog << " Parameter#0(lambda)=0";
     214        if (N0<=0)
     215            *fLog << " Parameter#1(N0)=0";
     216        *fLog << endl;
     217
     218        return kFALSE;
     219    }
     220
     221    //
     222    // -------------- start error calculation ----------------
     223    //
     224    Double_t emat[2][2];
     225    gMinuit->mnemat(&emat[0][0], 2);
     226
     227    //
     228    // Rdead : fraction of real time lost by the dead time
     229    // kappa = number of observed events (Nm) divided by
     230    //         the number of genuine events (N0)
     231    // Teff  : effective on-time
     232    //
     233    const Double_t Teff  = Nm/lambda;
     234    const Double_t kappa = Nm/N0;
     235    const Double_t Rdead = 1.0 - kappa;
     236
     237    const Double_t dldl   = emat[0][0];
     238    const Double_t dN0dN0 = emat[1][1];
     239
     240    const Double_t dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm);
     241    const Double_t dl = sqrt(dldl);
     242    const Double_t dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
     243    //
     244    // -------------- end error calculation ----------------
     245    //
     246
     247    // the effective on time is Nm/lambda
     248    fHEffOn.SetBinContent(i,  Teff);
     249    fHEffOn.SetBinError  (i, dTeff);
     250
     251    // plot chi2-probability of fit
     252    fHProb.SetBinContent(i, prob);
     253
     254    // lambda from fit
     255    fHLambda.SetBinContent(i, lambda);
     256    fHLambda.SetBinError  (i,     dl);
     257
     258    // Rdead from fit
     259    fHRdead.SetBinContent(i, Rdead);
     260    fHRdead.SetBinError  (i,dRdead);
     261
     262    return kTRUE;
     263}
     264
     265void MHEffOnTime::ResetBin(Int_t i)
     266{
     267    fHEffOn.SetBinContent (i, 1.e-20);
     268    fHProb.SetBinContent  (i, 1.e-20);
     269    fHLambda.SetBinContent(i, 1.e-20);
     270    fHRdead.SetBinContent (i, 1.e-20);
     271}
     272
    129273// -----------------------------------------------------------------------
    130274//
     
    132276// time differences
    133277//
    134 void MHEffOnTime::Calc(TH2D *hist, const Bool_t Draw)
    135 {
    136      // Draw != 0   means the distribution of time differences should be drawn
    137 
     278void MHEffOnTime::Calc(const TH2D *hist, const Bool_t draw)
     279{
    138280    // nbins = number of Var bins
    139281    const Int_t nbins = hist->GetNbinsY();
    140282
    141     // start  of 'for' loop ---------------------------
    142283    for (int i=1; i<=nbins; i++)
    143284    {
    144       TString strg0("Calc-");
    145       strg0 += fVarname;
    146       TH1D &h = *hist->ProjectionX(strg0, i, i, "E");
     285        TH1D &h = *((TH2D*)hist)->ProjectionX(TString("Calc-")+fVarname,
     286                                              i, i, "E");
     287        if (draw)
     288            DrawBin(h, i);
     289
     290        ResetBin(i);
    147291
    148292        // Nmdel = Nm * binwidth,  with Nm = number of observed events
     293        const Double_t Nm    = h.Integral();
    149294        const Double_t Nmdel = h.Integral("width");
    150         const Double_t Nm    = h.Integral();
    151 
    152         //...................................................
    153         // determine range (yq[0], yq[1]) of time differences
    154         // where fit should be performed;
    155         // require a fraction >=xq[0] of all entries to lie below yq[0]
    156         //     and a fraction <=xq[1] of all entries to lie below yq[1]; 
    157         // within the range (yq[0], yq[1]) there must be no empty bin;
    158         // choose pedestrian approach as long as GetQuantiles is not available
    159 
    160         Double_t xq[2] = { 0.15, 0.95 };
     295        if (Nm <= 0)
     296        {
     297            *fLog << warn << dbginf << "Nm<0 for bin #" << i << endl;
     298            delete &h;
     299            continue;
     300        }
     301
    161302        Double_t yq[2];
    162 
    163         // GetQuantiles doesn't seem to be available in root 3.01/06
    164         // h.GetQuantiles(2,yq,xq);
    165 
    166         const Int_t    jbins  = h.GetNbinsX();
    167 
    168         fHEffOn.SetBinContent (i, 1.e-20);
    169         fHProb.SetBinContent  (i, 1.e-20);
    170         fHLambda.SetBinContent(i, 1.e-20);
    171         fHRdead.SetBinContent (i, 1.e-20);
    172 
    173         // start  of 'if' loop ---------------------------
    174         if (Nm > 0.0)
    175         {
    176             // del is bin width in time difference
    177             const Double_t del = Nmdel/Nm;
    178 
    179             Double_t sum1 = 0.0;
    180             yq[0]  = h.GetBinLowEdge(jbins+1);
    181             for (int j=1; j<=jbins; j++)
    182             {
    183                 if (sum1 >= xq[0]*Nm)
    184                 {
    185                     yq[0] = h.GetBinLowEdge(j);
    186                     break;
    187                 }
    188                 sum1 += h.GetBinContent(j);
    189             }
    190        
    191             Double_t sum2 = 0.0;
    192             yq[1] = h.GetBinLowEdge(jbins+1);
    193             for (int j=1; j<=jbins; j++)
    194             {
    195                 const Double_t content = h.GetBinContent(j);
    196                 sum2 += content;
    197                 if (sum2 >= xq[1]*Nm || content == 0.0)
    198                 {
    199                     yq[1] = h.GetBinLowEdge(j);
    200                     break;
    201                 }
    202             }
    203 
    204             //...................................................
    205 
    206             // parameter 0 = lambda
    207             // parameter 1 = N0           with N0 = ideal number of events
    208             // parameter 2 = del (fixed)  with del = bin width of time difference
    209 
    210             TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
    211             func.FixParameter(2, del);
    212 
    213             func.SetParameter(0, 100); // [Hz]
    214             func.SetParameter(1, Nm);
    215 
    216             func.SetParLimits(0,   0,  1000);    // [Hz]
    217             func.SetParLimits(1,   0, 10*Nm);
    218 
    219             func.SetParName(0, "lambda");
    220             func.SetParName(1, "N0");
    221             func.SetParName(2, "del");
    222 
    223             // options : 0  (=zero) do not plot the function
    224             //           I  use integral of function in bin rather than value at bin center
    225             //           R  use the range specified in the function range
    226             //           Q  quiet mode
    227             h.Fit("Poisson", "0IRQ");
    228 
    229             const Double_t lambda = func.GetParameter(0);
    230             const Double_t N0     = func.GetParameter(1);
    231             // const Double_t chi2   = func.GetChisquare();
    232             const Double_t Prob   = func.GetProb();
    233             const Int_t    NDF    = func.GetNDF();
    234 
    235             // was fit successful ?
    236             if (NDF>0  &&  Prob>0.01  &&  lambda>0.0  &&  N0>0.0)
    237             {
    238                 // start error calculation .....................
    239 
    240                 Double_t emat[2][2];
    241                 // Double_t eplus, eminus, eparab, gcc;
    242                 gMinuit->mnemat(&emat[0][0], 2);
    243 
    244                 // for (int m=0; m<2; m++)
    245                 // {
    246                 //   *fLog << "emat[" << m << "][n] = ";
    247                 //   for (int n=0; n<2; n++)
    248                 //   {
    249                 //     *fLog << emat[m][n] << "  ";
    250                 //   }
    251                 //   *fLog << endl;
    252                 // }
    253 
    254                 // *fLog << "eplus, eminus, eparab, gcc :" << endl;
    255                 // for (int k=0; k<2; k++)
    256                 // {
    257                 //   gMinuit->mnerrs(k, eplus, eminus, eparab, gcc);
    258                 //   *fLog << eplus << " " << eminus << " " << eparab << " "
    259                 //         << gcc << endl;
    260                 // }
    261 
    262                 // Rdead : fraction of real time lost by the dead time
    263                 // kappa = number of observed events (Nm) divided by
    264                 //         the number of genuine events (N0)
    265                 // Teff  : effective on-time
    266 
    267                 Double_t Rdead, kappa, Teff, dTeff, dldl, dl, dRdead;
    268                 Double_t dN0dN0;
    269                 //Double_t dldN0,lN0corr;
    270                 Teff  = Nm/lambda;
    271                 kappa = Nm/N0;
    272                 Rdead = 1.0 - kappa;
    273 
    274                 dldl   = emat[0][0];
    275                 dN0dN0 = emat[1][1];
    276                 // dldN0  = emat[0][1];
    277 
    278                 dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm);
    279                 dl = sqrt(dldl);
    280                 dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
    281 
    282                 // lN0corr = dldN0 / sqrt(dldl * dN0dN0);
    283                 // *fLog << "MHEffOnTime : correlation between lambda and N0 = "
    284                 //       << lN0corr << endl;
    285 
    286                 // end error calculation .....................
    287 
    288                 // the effective on time is Nm/lambda
    289                 fHEffOn.SetBinContent(i,  Teff);
    290                 fHEffOn.SetBinError  (i, dTeff);
    291 
    292                 // plot chi2-probability of fit
    293                 fHProb.SetBinContent(i, Prob);
    294 
    295                 // lambda from fit
    296                 fHLambda.SetBinContent(i, lambda);
    297                 fHLambda.SetBinError  (i,     dl);
    298 
    299                 // Rdead from fit
    300                 fHRdead.SetBinContent(i, Rdead);
    301                 fHRdead.SetBinError  (i,dRdead);
    302             }
    303 
    304             //........................
    305             // draw the distribution of time differences if requested
    306             //
    307             if (Draw == kTRUE)
    308               {
    309                 char txt[100];
    310                 TString strg1(fVarname);
    311                 strg1 += "-bin %d";
    312                 sprintf(txt, strg1, i);
    313                 new TCanvas(txt, txt);
    314                 // TCanvas &c = *MakeDefCanvas(txt, txt);
    315                 // gROOT->SetSelectedPad(NULL);
    316 
    317                 gPad->SetLogy();
    318                 gStyle->SetOptFit(1011);
    319                 // gStyle->SetOptStat(1);
    320 
    321                 h.SetName(txt);
    322                 h.SetXTitle("time difference [s]");
    323                 h.SetYTitle("Counts");
    324                 h.DrawCopy();
    325 
    326                 func.SetRange(yq[0], yq[1]); // Range of Drawing
    327                 func.DrawCopy("same");
    328 
    329                 // c.Modified();
    330                 // c.Update();
    331               }
    332             //........................
    333         }
    334         // end  of 'if' loop ---------------------------
     303        GetQuantiles(h, 0.15, 0.95, yq);
     304
     305        //
     306        // Setup Poisson function for the fit:
     307        // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
     308        //
     309        TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
     310        func.SetParNames("lambda", "N0", "del");
     311
     312        func.SetParameter(0, 1);
     313        func.SetParameter(1, Nm);
     314        func.FixParameter(2, Nmdel/Nm);
     315
     316        // For me (TB) it seems that setting parameter limits isn't necessary
     317
     318        // For a description of the options see TH1::Fit
     319        h.Fit("Poisson", "0IRQ");
     320
     321        // Calc and fill results of fit into the histograms
     322        const Bool_t rc = CalcResults(func, Nm, i);
     323
     324        // draw the distribution of time differences if requested
     325        if (rc && draw)
     326            func.DrawCopy("same");
    335327
    336328        delete &h;
    337329    }
    338     // end  of 'for' loop ---------------------------
    339     //    gStyle->SetOptStat(1111);
    340 
    341330}
    342331
     
    347336Bool_t MHEffOnTime::SetupFill(const MParList *plist)
    348337{
    349     TString strg0("Binning");
    350     strg0 += fVarname;
    351 
    352     const MBinning* binsVar = (MBinning*)plist->FindObject(strg0);
     338    TString bn = "Binning"+fVarname;
     339
     340    const MBinning* binsVar = (MBinning*)plist->FindObject(bn);
    353341
    354342    if (!binsVar)
    355343    {
    356         *fLog << err << dbginf << strg0
    357               << " [MBinning] not found... aborting." << endl;
     344        *fLog << err << dbginf << bn << " [MBinning] not found... aborting." << endl;
    358345        return kFALSE;
    359346    }
     
    387374TObject *MHEffOnTime::DrawClone(Option_t *opt) const
    388375{
    389     TString strg0("EffOnTime");
    390     strg0 += fVarname;   
    391     TString strg1("Results from on time fit vs. ");
    392     strg1 += fVarname;   
    393 
    394     TCanvas &c = *MakeDefCanvas(strg0, strg1);
     376    TCanvas &c = *MakeDefCanvas(TString("EffOnTime")+fVarname,
     377                                TString("Results from on time fit vs. ")+fVarname);
    395378    c.Divide(2, 2);
    396379
     
    421404void MHEffOnTime::Draw(Option_t *opt)
    422405{
    423     TString strg0("EffOnTime");
    424     strg0 += fVarname;   
    425     TString strg1("Results from on time fit vs. ");
    426     strg1 += fVarname;   
    427 
    428406    if (!gPad)
    429         MakeDefCanvas(strg0, strg1);
     407        MakeDefCanvas(TString("EffOnTime")+fVarname,
     408                      TString("Results from on time fit vs. ")+fVarname);
    430409
    431410    gPad->Divide(2,2);
     
    447426}
    448427
    449 
    450 
    451 
    452 
    453 
     428void MHEffOnTime::Calc(const MHTimeDiffTheta &hist, const Bool_t Draw)
     429{
     430    Calc(hist.GetHist(), Draw);
     431}
     432
  • trunk/MagicSoft/Mars/mhist/MHEffOnTime.h

    r1411 r1668  
    1010
    1111class TH2D;
     12class MHTimeDiffTheta;
    1213class MParList;
    1314
     
    2021    TH1D fHRdead;
    2122
    22     const char *fVarname;
    23     const char *fUnit;
     23    TString fVarname;
     24    TString fUnit;
     25
     26    void ResetBin(Int_t i);
     27    Bool_t CalcResults(const TF1 &func, Double_t Nm, Int_t i);
     28    void DrawBin(TH1 &h, Int_t i) const;
     29    void GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const;
    2430
    2531public:
     
    3238    const TH1D *GetHist() const { return &fHEffOn; }
    3339
    34     void Calc(TH2D *hist, const Bool_t Draw);
     40    void Calc(const TH2D *hist, const Bool_t Draw);
     41    void Calc(const MHTimeDiffTheta &hist, const Bool_t Draw);
    3542
    3643    void Draw(Option_t *option="");
    3744    TObject *DrawClone(Option_t *option="") const;
    3845
    39     ClassDef(MHEffOnTime, 1) //1D-plot of Delta t vs. Var
     46    ClassDef(MHEffOnTime, 0) //1D-plot of Delta t vs. Var
    4047};
    4148
  • trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h

    r1292 r1668  
    3535    TObject *DrawClone(Option_t *option="") const;
    3636
    37     ClassDef(MHEffOnTimeTheta, 1) //1D-plot of Delta t vs. Theta
     37    ClassDef(MHEffOnTimeTheta, 0) //1D-plot of Delta t vs. Theta
    3838};
    3939
  • trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h

    r1294 r1668  
    3535    TObject *DrawClone(Option_t *option="") const;
    3636
    37     ClassDef(MHEffOnTimeTime, 1) //1D-plot of Delta t vs. time
     37    ClassDef(MHEffOnTimeTime, 0) //1D-plot of Delta t vs. time
    3838};
    3939
  • trunk/MagicSoft/Mars/mhist/MHFlux.cc

    r1604 r1668  
    4242#include <TProfile.h>
    4343
    44 
    4544#include <TCanvas.h>
    4645
     
    5352#include "MLogManip.h"
    5453
     54#include "MHThetabarTheta.h"
     55#include "MHEffOnTime.h"
     56#include "MHGamma.h"
     57
    5558ClassImp(MHFlux);
    5659
     60MHFlux::MHFlux(const MHGamma &hist, const TString varname, const TString unit)
     61    : fHOrig(), fHUnfold(), fHFlux()
     62{
     63    const TH2D &h2d = *hist.GetProject();
     64
     65    if (varname.IsNull() || unit.IsNull())
     66        *fLog << warn << dbginf << "varname or unit not defined" << endl;
     67
     68    fVarname = varname;
     69    fUnit    = unit;
     70
     71    // char txt[100];
     72
     73    // original distribution of E-est for different bins
     74    //                       of the variable (Theta or time)
     75    // sprintf(txt, "gammas vs. E-est and %s",varname);
     76
     77    ((TH2D&)h2d).Copy(fHOrig);
     78
     79    fHOrig.SetName("E_est");
     80    fHOrig.SetTitle(TString("No.of gammas vs. E-est and ")+fVarname);
     81
     82    fHOrig.SetDirectory(NULL);
     83    fHOrig.SetXTitle("E_{est} [GeV]");
     84    fHOrig.SetYTitle(fVarname+fUnit);
     85    //fHOrig.Sumw2();
     86
     87    SetBinning((TH2*)&fHOrig, (TH2*)&h2d);
     88
     89    fHOrig.Copy(fHUnfold);
     90
     91    // unfolded distribution of E-unfold for different bins
     92    //                       of the variable (Theta or time)
     93    // sprintf(txt, "gammas vs. E-unfold and %s",varname);
     94    fHUnfold.SetName("E-unfolded");
     95    fHUnfold.SetTitle(TString("No.of gammas vs. E-unfold and ")+fVarname);
     96
     97    fHUnfold.SetDirectory(NULL);
     98    fHUnfold.SetXTitle("E_{unfold} [GeV]");
     99    fHUnfold.SetYTitle(fVarname+fUnit);
     100    //fHUnfold.Sumw2();
     101   
     102    SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);
     103
     104
     105    // absolute photon flux vs. E-unfold
     106    //          for different bins of the variable (Theta or time)
     107    //
     108    // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
     109    fHFlux.SetName("photon flux");
     110    fHFlux.SetTitle(TString("Gamma flux [1/(s m^2 GeV) vs. E-unfold and ")+fVarname);
     111
     112    fHFlux.SetDirectory(NULL);
     113    fHFlux.SetXTitle("E_{unfold} [GeV]");
     114    fHFlux.SetYTitle(fVarname+fUnit);
     115    fHFlux.Sumw2();
     116
     117    SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);
     118}
    57119
    58120// --------------------------------------------------------------------------
     
    61123//                      and the units for the variable
    62124//
    63 MHFlux::MHFlux(const TH2D &h2d,  const Bool_t Draw,
    64                const TString varname, const TString unit)
     125MHFlux::MHFlux(const TH2D &h2d, const TString varname, const TString unit)
    65126    : fHOrig(), fHUnfold(), fHFlux()
    66127{
     
    70131    fVarname = varname;
    71132    fUnit    = unit;
    72 
    73     TString strg(varname);
    74     strg += unit;
    75133
    76134    // char txt[100];
     
    80138    // sprintf(txt, "gammas vs. E-est and %s",varname);
    81139
    82     TString strg1 = "no.of gammas vs. E-est and ";
    83     strg1 += varname;
    84    
    85140    ((TH2D&)h2d).Copy(fHOrig);
    86141
    87     fHOrig.SetName("E-est");
    88     fHOrig.SetTitle(strg1);
     142    fHOrig.SetName("E_est");
     143    fHOrig.SetTitle(TString("No.of gammas vs. E-est and ")+fVarname);
    89144
    90145    fHOrig.SetDirectory(NULL);
    91     fHOrig.SetXTitle("E-est [GeV]");
    92     fHOrig.SetYTitle(strg);
    93     fHOrig.Sumw2();
     146    fHOrig.SetXTitle("E_{est} [GeV]");
     147    fHOrig.SetYTitle(fVarname+fUnit);
     148    //fHOrig.Sumw2();
     149
     150    // copy fHOrig into fHUnfold in case no unfolding is done
     151    fHOrig.Copy(fHUnfold);
    94152
    95153    SetBinning((TH2*)&fHOrig, (TH2*)&h2d);
     
    99157    //                       of the variable (Theta or time)
    100158    // sprintf(txt, "gammas vs. E-unfold and %s",varname);
    101     TString strg2 = "no.of gammas vs. E-unfold and ";
    102     strg2 += varname;
    103 
    104159    fHUnfold.SetName("E-unfolded");
    105     fHUnfold.SetTitle(strg2);
     160    fHUnfold.SetTitle(TString("No.of gammas vs. E-unfold and ")+fVarname);
    106161
    107162    fHUnfold.SetDirectory(NULL);
    108     fHUnfold.SetXTitle("E-unfold [GeV]");
    109     fHUnfold.SetYTitle(strg);
    110     fHUnfold.Sumw2();
     163    fHUnfold.SetXTitle("E_{unfold} [GeV]");
     164    fHUnfold.SetYTitle(fVarname+fUnit);
     165    //fHUnfold.Sumw2();
    111166   
    112167    SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);
     
    117172    //
    118173    // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
    119     TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and ";
    120     strg3 += varname;
    121 
    122174    fHFlux.SetName("photon flux");
    123     fHFlux.SetTitle(strg3);
     175    fHFlux.SetTitle(TString("Gamma flux [1/(s m^{2} GeV)] vs. E-unfold and ")+fVarname);
    124176
    125177    fHFlux.SetDirectory(NULL);
    126     fHFlux.SetXTitle("E-unfold [GeV]");
    127     fHFlux.SetYTitle(strg);
     178    fHFlux.SetXTitle("E_{unfold} [GeV]");
     179    fHFlux.SetYTitle(fVarname+fUnit);
    128180    fHFlux.Sumw2();
    129181
    130182    SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);
    131 
    132 
    133     // copy fHOrig into fHUnfold in case no unfolding is done
    134     const Int_t nEunf   = fHUnfold.GetNbinsX();
    135     const Int_t nVar    = fHUnfold.GetNbinsY();
    136     for (int m=1; m<=nEunf; m++)
    137     {
    138       for (int n=1; n<=nVar; n++)
    139       {
    140         Double_t cont  = fHOrig.GetBinContent(m,n);
    141         Double_t dcont = fHOrig.GetBinError(m,n);
    142         fHUnfold.SetBinContent(m,n,cont);
    143         fHUnfold.SetBinError(m,n,dcont);
    144       }
    145     }
    146   //..............................................
    147   // draw the No.of photons vs. E-est
    148   // for the individual bins of the variable Var
    149 
    150   if (Draw == kTRUE)
    151   {
    152     const Int_t nVar = fHOrig.GetNbinsY();
    153 
    154     for (int n=1; n<=nVar; n++)
    155     {
    156       TString strg0("Orig-");
    157       strg0 += fVarname;
    158       TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");
    159 
    160       strg0  = fVarname;
    161       strg0 += "-bin ";
    162       strg0 += n;
    163 
    164       TString strg1("No.of photons vs. E-est for ");
    165       strg1 += strg0;
    166 
    167       new TCanvas(strg0, strg1);
    168       // TCanvas &c = *MakeDefCanvas(txt0, strg1);
    169       // gROOT->SetSelectedPad(NULL);
    170 
    171       gPad->SetLogx();
    172 
    173       h.SetName(strg0);
    174       h.SetTitle(strg1);
    175       h.SetXTitle("E-est [GeV]");
    176       h.SetYTitle("No.of photons");
    177       h.DrawCopy();
    178 
    179       // c.Modified();
    180       // c.Update();
    181     }
    182   }
    183   //........................
    184183}
    185184
     
    200199Bool_t MHFlux::Fill(const MParContainer *par)
    201200{
    202    return kTRUE;
     201    return kTRUE;
    203202}
    204203
     
    208207// Unfold the distribution in E-est
    209208//
    210 void MHFlux::Unfold(const Bool_t Draw)
    211 {
    212   //..............................................
    213   // draw the No.of photons vs. E-unfold
    214   // for the individual bins of the variable Var
    215 
    216   if (Draw == kTRUE)
    217   {
    218     const Int_t nVar = fHUnfold.GetNbinsY();
    219 
    220     for (int n=1; n<=nVar; n++)
    221     {
    222       TString strg0("Unfold-");
    223       strg0 += fVarname;
    224       TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
    225       strg0  = fVarname;
    226       strg0 += "-bin ";
    227       strg0 += n;
    228 
    229       TString strg1("No.of photons vs. E-unfold for ");
    230       strg1 += strg0;
    231 
    232       new TCanvas(strg0, strg1);
    233 
    234       // TCanvas &c = *MakeDefCanvas(txt0, strg1);
    235       // gROOT->SetSelectedPad(NULL);
    236 
    237       gPad->SetLogx();
    238 
    239       h.SetName(strg0);
    240       h.SetTitle(strg1);
    241       h.SetXTitle("E-unfold [GeV]");
    242       h.SetYTitle("No.of photons");
    243       h.DrawCopy();
    244 
    245       // c.Modified();
    246       // c.Update();
    247     }
    248   }
    249   //........................
    250 }
    251 
     209void MHFlux::Unfold()
     210{
     211}
     212
     213void MHFlux::CalcFlux(const MHEffOnTime &teff, const MHThetabarTheta &thetabar,
     214                      const TH2D *aeff)
     215{
     216    CalcFlux(teff.GetHist(), thetabar.GetHist(), aeff);
     217}
     218
     219Double_t MHFlux::ParabInterpolLog(const TAxis &axe, Int_t j,
     220                                  Double_t y[], Double_t Ebar) const
     221{
     222    const double t1 = log10(axe.GetBinLowEdge(j-1)) + log10(axe.GetBinUpEdge(j-1));
     223    const double t2 = log10(axe.GetBinLowEdge(j))   + log10(axe.GetBinUpEdge(j));
     224    const double t3 = log10(axe.GetBinLowEdge(j+1)) + log10(axe.GetBinUpEdge(j+1));
     225
     226    const Double_t lebar = log10(Ebar);
     227
     228    return Parab(t1/2, t2/2, t3/2, y[j-2], y[j-1], y[j], lebar);
     229}
     230
     231// --------------------------------------------------------------------
     232//
     233//  determine bins for interpolation (k3 is the middle one) in bar.
     234//  k0 denotes the bin from which the error is copied
     235//
     236void MHFlux::FindBins(const TAxis &axe, const Double_t bar, Int_t &k3, Int_t &k0) const
     237{
     238    const Int_t n = axe.GetNbins();
     239
     240    k3 = axe.FindFixBin(bar);
     241    k0 = k3;
     242
     243    if (k3<2)
     244    {
     245        k3 = 2;
     246        if (bar<axe.GetBinLowEdge(2))
     247            k0 = 1;
     248    }
     249
     250    if (k3>n-1)
     251    {
     252        k3 = n-1;
     253        if (bar>axe.GetBinLowEdge(n))
     254            k0 = n;
     255    }
     256
     257    if (bar>=axe.GetBinLowEdge(1) && bar<=axe.GetBinUpEdge(n))
     258        return;
     259
     260    *fLog << dbginf << "extrapolation: bar = " << bar;
     261    *fLog << ", min =" << axe.GetBinLowEdge(1);
     262    *fLog << ", max =" << axe.GetBinUpEdge(n) << endl;
     263}
     264
     265Double_t MHFlux::ParabInterpolCos(const TAxis &axe, const TH2D *aeff, Int_t j, Int_t k3, Double_t val) const
     266{
     267    const double t1 = cos( axe.GetBinCenter (k3-1) );
     268    const double t2 = cos( axe.GetBinCenter (k3)   );
     269    const double t3 = cos( axe.GetBinCenter (k3+1) );
     270
     271    const double a1 = aeff->GetBinContent(j, k3-1);
     272    const double a2 = aeff->GetBinContent(j, k3);
     273    const double a3 = aeff->GetBinContent(j, k3+1);
     274
     275    return Parab(t1, t2, t3, a1, a2, a3, val);
     276}
    252277
    253278// -------------------------------------------------------------------------
     
    259284//
    260285void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar,
    261                       const TH2D *aeff, const Bool_t Draw)
    262 {
    263   // Note that fHUnfold  has bins in Eunf and Var
    264   //           *teff     has bins in Var  (the same bins in Var as fHUnfold)
    265   //           *thetabar has bins in Var  (the same bins in Var as fHUnfold)
    266   //           *aeff     has bins in Etru and Theta
    267   //                     (where in general the binning in Etru is different
    268   //                      from the binning in Eunf)
    269   // The variable Var may be 'time' or 'Theta'
    270 
    271   // Draw = kTRUE means the differential flux vs E-unf should be drawn
    272   //              for the individual bins of the variable Var
     286                      const TH2D *aeff)
     287{
     288    //
     289    // Note that fHUnfold  has bins in Eunf and Var
     290    //           *teff     has bins in Var  (the same bins in Var as fHUnfold)
     291    //           *thetabar has bins in Var  (the same bins in Var as fHUnfold)
     292    //           *aeff     has bins in Etru and Theta
     293    //                     (where in general the binning in Etru is different
     294    //                      from the binning in Eunf)
     295    // The variable Var may be 'time' or 'Theta'
    273296
    274297    const TAxis &axex = *((TH2*)aeff)->GetXaxis();
    275298    const TAxis &axey = *((TH2*)aeff)->GetYaxis();
    276299
    277   //....................................
    278   // define dummy histogram *aeff
    279   ((TH1*)aeff)->Sumw2();
    280   MBinning binsetru("BinningEtru");
    281   binsetru.SetEdgesLog(10, 10, 1e3);
    282 
    283   MBinning binsthetatru("BinningThetatru");
    284   binsthetatru.SetEdges(7, -2.5, 32.5);
    285   //SetBinning((TH1*)aeff, &binsetru, &binsthetatru);
    286   SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
    287 
    288   const Int_t netru  = aeff->GetNbinsX();
    289   const Int_t ntheta = aeff->GetNbinsY();
    290 
    291   for (int j=1; j<=netru; j++)
    292   {
    293       for (int k=1; k<=ntheta; k++)
    294       {
    295           Double_t cont = 10000.0;
    296           ((TH1*)aeff)->SetBinContent(j, k, cont);
    297 
    298           Double_t dcont = 100.0;
    299           ((TH1*)aeff)->SetBinError(j, k, dcont);
    300       }
    301   }
    302   // *fLog << "Dummy aeff : netru =" << netru << ",  ntheta = " << ntheta << endl;
    303   //....................................
    304 
    305   // number of Eunf and Var bins   (histograms : fHUnfold, fHFlux)
    306   const Int_t nEunf    = fHFlux.GetNbinsX();
    307   const Int_t nVar     = fHFlux.GetNbinsY();
    308 
    309   // number of Etru and Theta bins (histogram *aeff of collection area)
    310 
    311   const Int_t nEtru    = aeff->GetNbinsX();
    312   const Int_t nTheta   = aeff->GetNbinsY();
    313 
    314   //*fLog << "nEunf =" << nEunf << ",  nVar = "   << nVar   << endl;
    315   //*fLog << "nEtru =" << nEtru << ",  nTheta = " << nTheta << endl;
    316 
    317   //...................................................................
    318   // calculate effective collection area
    319   //    for the Eunf and Var bins of the histogram fHUnfold
    320   //    from the histogram *aeff, which has bins in Etru and Theta
    321   // the result is the histogram fHAeff
    322   //
    323 
    324   TH2D fHAeff;
    325   fHAeff.Sumw2();
    326   SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
    327 
    328   Double_t *aeffbar  = new Double_t[nEtru];
    329   Double_t *daeffbar = new Double_t[nEtru];
    330 
    331   Double_t aeffEunfVar;
    332   Double_t daeffEunfVar;
    333 
    334   //------   start n loop   ------
    335   for (int n=1; n<=nVar; n++)
    336   {
    337     Double_t Thetabar = thetabar->GetBinContent(n);
    338     Double_t cosThetabar = cos(Thetabar);
    339 
    340     // determine Theta bins (k1, k2, k3) for interpolation in Theta
    341     // k0 denotes the Theta bin from whicvh the error is copied
    342     Int_t k0=0;
    343     Int_t k1=0;
    344     Int_t k2=0;
    345     Int_t k3=0;
    346 
    347     for (int k=3; k<=nTheta; k++)
    348     {
    349         Double_t Thetalow = axey.GetBinLowEdge(k);
    350         if (Thetabar < Thetalow)
     300    if (axex.GetNbins()<3)
     301    {
     302        *fLog << err << "ERROR - Number of Energy bins <3 not implemented!" << endl;
     303        return;
     304    }
     305
     306    if (axey.GetNbins()<3)
     307        *fLog << warn << "WARNING - Less than 3 theta-bins not supported very well!" << endl;
     308
     309    //
     310    // calculate effective collection area
     311    //    for the Eunf and Var bins of the histogram fHUnfold
     312    //    from the histogram *aeff, which has bins in Etru and Theta
     313    // the result is the histogram fHAeff
     314    //
     315    TH2D fHAeff;
     316    SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
     317    fHAeff.Sumw2();
     318
     319    //
     320    // ------ start loops ------
     321    //
     322    const Int_t nEtru = aeff->GetNbinsX();
     323
     324    Double_t *aeffbar  = new Double_t[nEtru];
     325    Double_t *daeffbar = new Double_t[nEtru];
     326
     327    const Int_t nVar = fHFlux.GetNbinsY();
     328    for (int n=1; n<=nVar; n++)
     329    {
     330        const Double_t tbar = thetabar->GetBinContent(n);
     331        const Double_t costbar = cos(tbar/kRad2Deg);
     332
     333        // determine bins for interpolation (k3, k0)
     334        Int_t kv, ke;
     335        FindBins(axey, tbar, kv, ke);
     336
     337        //
     338        // calculate effective collection area at Theta = Thetabar
     339        // by quadratic interpolation in cos(Theta);
     340        // do this for each bin of Etru
     341        //
     342        for (int j=1; j<=nEtru; j++)
    351343        {
    352             k1 = k-2;
    353             k2 = k-1;
    354             k3 = k;
    355             k0 = k2;
    356             break;
    357         }
    358     }
    359 
    360     if (k3 == 0)
    361     {
    362         k1 = nTheta-2;
    363         k2 = nTheta-1;
    364         k3 = nTheta;
    365         k0 = k2;
    366     }
    367 
    368     if (Thetabar <  axey.GetBinLowEdge(2))
    369       k0 = 1;
    370     else
    371         if (Thetabar >  axey.GetBinLowEdge(nTheta))
    372             k0 = nTheta;
    373 
    374     Double_t Thetamin = axey.GetBinLowEdge(1);
    375     Double_t Thetamax = axey.GetBinLowEdge(nTheta+1);
    376     if (Thetabar < Thetamin  ||  Thetabar > Thetamax)
    377     {
    378       *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar
    379             << ",  Thetamin =" << Thetamin
    380             << ",  Thetamax =" << Thetamax << endl;
    381     }
    382 
    383     //*fLog << "Var bin "   << n  << ":"  <<  endl;
    384     //*fLog << "Thetabar= " << Thetabar
    385     //      << ",  k0= "    << k0
    386     //      << ",  k1= "    << k1
    387     //      << ",  k2= "    << k2
    388     //      << ",  k3= "    << k3         <<  endl;
    389  
    390 
    391     // calculate effective collection area at Theta = Thetabar
    392     // by quadratic interpolation in cos(Theta);
    393     // do this for each bin of Etru
    394     //
    395 
    396     for (int j=1; j<=nEtru; j++)
    397     {
    398         double c0 = 0;
    399         double c1 = 0;
    400         double c2 = 0;
    401 
    402         const double t1 = cos( axey.GetBinCenter (k1) );
    403         const double t2 = cos( axey.GetBinCenter (k2) );
    404         const double t3 = cos( axey.GetBinCenter (k3) );
    405 
    406         const double a1 = aeff->GetBinContent(j, k1);
    407         const double a2 = aeff->GetBinContent(j, k2);
    408         const double a3 = aeff->GetBinContent(j, k3);
    409 
    410         Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
    411         aeffbar[j]  = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
    412         daeffbar[j] = aeff->GetBinError(j,k0);
    413 
    414         //*fLog << "Etru bin " << j <<  ":  tbar= " << Thetabar
    415         //      << ",  abar= "  << aeffbar[j]
    416         //      << ",  dabar= " << daeffbar[j] << endl;
    417     }
    418 
    419     //---   start m loop ---
    420     // calculate effective collection area at (E = Ebar, Theta = Thetabar)
    421     // by quadratic interpolation in log10(Etru)
    422     // do this for each bin of Eunf
    423     //
    424     for (int m=1; m<=nEunf; m++)
    425     {
    426         Double_t log10Ebar = 0.5 * ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m)  )+
    427                                      log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
    428         Double_t Ebar = pow(10.0, log10Ebar);
    429 
    430         // determine Etru bins (j1, j2, j3) for interpolation in E
    431         // j0 denotes the Etru bin from which the error is copied
    432         Int_t j0=0;
    433         Int_t j1=0;
    434         Int_t j2=0;
    435         Int_t j3=0;
    436 
    437         for (int j=3; j<=nEtru; j++)
    438         {
    439             Double_t Elow = axex.GetBinLowEdge(j);
    440             if (Ebar < Elow)
     344            if (axey.GetNbins()<3)
    441345            {
    442                 j1 = j-2;
    443                 j2 = j-1;
    444                 j3 = j;
    445                 j0 = j2;
    446                 break;
     346                // FIXME: Other interpolation?
     347                aeffbar[j-1]  = aeff->GetBinContent(j, n);
     348                daeffbar[j-1] = aeff->GetBinError(j, n);
     349            }
     350            else
     351            {
     352                aeffbar[j-1]  = ParabInterpolCos(axey, aeff, j, kv, costbar);
     353                daeffbar[j-1] = aeff->GetBinError(j, ke);
    447354            }
    448355        }
    449356
    450         if (j3 == 0)
     357        //
     358        // calculate effective collection area at (E = Ebar, Theta = Thetabar)
     359        // by quadratic interpolation in log10(Etru)
     360        // do this for each bin of Eunf
     361        //
     362        CalcEffCol(axex, fHAeff, n, aeffbar, daeffbar);
     363    }
     364
     365    delete aeffbar;
     366    delete daeffbar;
     367
     368    //
     369    // now calculate the absolute gamma flux
     370    //
     371    CalcAbsGammaFlux(*teff, fHAeff);
     372}
     373
     374// --------------------------------------------------------------------
     375//
     376//  calculate effective collection area at (E = Ebar, Theta = Thetabar)
     377//  by quadratic interpolation in log10(Etru)
     378//  do this for each bin of Eunf
     379//
     380void MHFlux::CalcEffCol(const TAxis &axex, TH2D &fHAeff, Int_t n, Double_t aeffbar[], Double_t daeffbar[])
     381{
     382    const Int_t nEunf = fHFlux.GetNbinsX();
     383
     384    const TAxis &unfx = *fHUnfold.GetXaxis();
     385
     386    for (int m=1; m<=nEunf; m++)
     387    {
     388        const Double_t Ebar = GetBinCenterLog(unfx, m);
     389
     390        Int_t j0, j3;
     391        FindBins(axex, Ebar, j3, j0);
     392
     393        const Double_t v = ParabInterpolLog(axex, j3, aeffbar, Ebar);
     394
     395        fHAeff.SetBinContent(m,n, v);
     396        fHAeff.SetBinError(m,n, daeffbar[j0-1]);
     397    }
     398}
     399
     400// --------------------------------------------------------------------
     401//
     402//  calculate the absolute gamma flux
     403//
     404void MHFlux::CalcAbsGammaFlux(const TH1D &teff, const TH2D &fHAeff)
     405{
     406    const Int_t nEunf = fHFlux.GetNbinsX();
     407    const Int_t nVar  = fHFlux.GetNbinsY();
     408
     409    for (int m=1; m<=nEunf; m++)
     410    {
     411        const Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
     412
     413        for (int n=1; n<=nVar; n++)
    451414        {
    452             j1 = nEtru-2;
    453             j2 = nEtru-1;
    454             j3 = nEtru;
    455             j0 = j2;
     415            const Double_t Ngam  = fHUnfold.GetBinContent(m,n);
     416            const Double_t Aeff  = fHAeff.GetBinContent(m,n);
     417            const Double_t Effon = teff.GetBinContent(n);
     418
     419            const Double_t c1 = fHUnfold.GetBinError(m,n)/Ngam;
     420            const Double_t c2 = teff.GetBinError(n)      /Effon;
     421            const Double_t c3 = fHAeff.GetBinError(m,n)  /Aeff;
     422
     423            const Double_t cont  = Ngam / (DeltaE * Effon * Aeff);
     424            const Double_t dcont = sqrt(c1*c1 + c2*c2 + c3*c3);
     425
     426            //
     427            // Out of Range
     428            //
     429            const Bool_t oor = Ngam<=0 || DeltaE<=0 || Effon<=0 || Aeff<=0;
     430
     431            if (oor)
     432                *fLog << warn << "MHFlux::CalcAbsGammaFlux(" << m << "," << n << ") ";
     433
     434            if (Ngam<=0)
     435                *fLog << " Ngam=0";
     436            if (DeltaE<=0)
     437                *fLog << " DeltaE=0";
     438            if (Effon<=0)
     439                *fLog << " Effon=0";
     440            if (Aeff<=0)
     441                *fLog << " Aeff=0";
     442
     443            if (oor)
     444                *fLog << endl;
     445
     446            fHFlux.SetBinContent(m,n, oor ? 1e-20 : cont);
     447            fHFlux.SetBinError(m,n,   oor ? 1e-20 : dcont*cont);
    456448        }
    457 
    458         if (Ebar <  axex.GetBinLowEdge(2))
    459             j0 = 1;
    460         else
    461             if (Ebar >  axex.GetBinLowEdge(nEtru))
    462                 j0 = nEtru;
    463 
    464         Double_t Etrumin = axex.GetBinLowEdge(1);
    465         Double_t Etrumax = axex.GetBinLowEdge(nEtru+1);
    466         if (Ebar < Etrumin  ||  Ebar > Etrumax)
    467         {
    468             *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
    469                 << ",  Etrumin =" << Etrumin
    470                 << ",  Etrumax =" << Etrumax << endl;
    471         }
    472 
    473         //*fLog << "Var bin "   << n  << ":"  <<  endl;
    474         //*fLog << "Ebar= " << Ebar
    475         //      << ",  j1= "    << j1
    476         //      << ",  j2= "    << j2
    477         //      << ",  j3= "    << j3         <<  endl;
    478 
    479 
    480         double c0=0.0;
    481         double c1=0.0;
    482         double c2=0.0;
    483 
    484         const double t1 = 0.5 * ( log10( axex.GetBinLowEdge (j1)  )+
    485                                   log10( axex.GetBinLowEdge (j1+1)) );
    486         const double t2 = 0.5 * ( log10( axex.GetBinLowEdge (j2)  )+
    487                                   log10( axex.GetBinLowEdge (j2+1)) );
    488         const double t3 = 0.5 * ( log10( axex.GetBinLowEdge (j3)  )+
    489                                   log10( axex.GetBinLowEdge (j3+1)) );
    490 
    491         const double a1 = aeffbar[j1];
    492         const double a2 = aeffbar[j2];
    493         const double a3 = aeffbar[j3];
    494 
    495         Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
    496         aeffEunfVar  = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
    497         daeffEunfVar = daeffbar[j0];
    498 
    499         //*fLog << "Eunf bin " << m     <<  ":  Ebar= " << Ebar
    500         //      << ",  aeffEunfVar = "  << aeffEunfVar
    501         //      << ",  daeffEunfVar = " << daeffEunfVar << endl;
    502 
    503         fHAeff.SetBinContent(m,n,aeffEunfVar);
    504         fHAeff.SetBinError(m,n,daeffEunfVar);
    505     }
    506     //---   end m loop ---
    507   }
    508   //------   end n loop   ------
    509   delete aeffbar;
    510 
    511   //...................................................................
    512   // now calculate the absolute gamma flux
    513   //
    514   for (int m=1; m<=nEunf; m++)
    515   {
    516       Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
    517 
    518       for (int n=1; n<=nVar; n++)
    519       {
    520           Double_t Ngam   = fHUnfold.GetBinContent(m,n);
    521           Double_t dNgam  = fHUnfold.GetBinError(m,n);
    522 
    523           Double_t Aeff   = fHAeff.GetBinContent(m,n);
    524           Double_t dAeff  = fHAeff.GetBinError(m,n);
    525 
    526           Double_t Effon  = teff->GetBinContent(n);
    527           Double_t dEffon = teff->GetBinError(n);
    528 
    529           Double_t Cont, dCont;
    530           if (Ngam>0 && DeltaE>0 && Effon>0 && Aeff>0)
    531           {
    532               Cont  = Ngam / (DeltaE * Effon * Aeff);
    533               dCont = Cont * sqrt( dNgam *dNgam  / (Ngam*Ngam) +
    534                                    dEffon*dEffon / (Effon*Effon) +
    535                                    dAeff *dAeff  / (Aeff*Aeff)  );
    536           }
    537           else
    538           {
    539               Cont  = 1.e-20;
    540               dCont = 1.e-20;
    541           }
    542 
    543           fHFlux.SetBinContent(m,n,Cont);
    544           fHFlux.SetBinError(m,n,dCont);
    545 
    546           //*fLog << "Eunf bin "    << m      << ",  Var bin " << n
    547           //      << ":  Ngam = "   << Ngam   << ",  Flux = "
    548           //      << Cont  << ", dFlux = " << dCont << endl;
    549           //*fLog << ",  DeltaE = " << DeltaE << ",  Effon = " << Effon
    550           //      << ",  Aeff = "   << Aeff   << endl;
    551       }
    552   }
    553 
    554   //..............................................
    555   // draw the differential photon flux vs. E-unf
    556   // for the individual bins of the variable Var
    557 
    558   if (Draw == kTRUE)
    559   {
    560       for (int n=1; n<=nVar; n++)
    561       {
    562           TString strg0("Flux-");
    563           strg0 += fVarname;
    564 
    565           TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
    566 
    567           TString strg1("Photon flux vs. E-unfold for ");
    568           TString strg2 = fVarname;
    569 
    570           strg2 += "-bin ";
    571           strg2 += n;
    572 
    573           TString strg3 = strg1 + strg2;
    574           new TCanvas(strg2, strg3);
    575           // TCanvas &c = *MakeDefCanvas(txt, txt);
    576           // gROOT->SetSelectedPad(NULL);
    577 
    578           gPad->SetLogx();
    579 
    580           h.SetName(strg2);
    581           h.SetTitle(strg3);
    582           h.SetXTitle("E-unfold [GeV]            ");
    583           h.SetYTitle("photons / (s m2 GeV)");
    584           h.DrawCopy();
    585 
    586           // c.Modified();
    587           // c.Update();
    588       }
    589   }
    590   //........................
     449    }
     450}
     451
     452// --------------------------------------------------------------------
     453//
     454// draw the differential photon flux vs. E-unf
     455// for the individual bins of the variable Var
     456//
     457void MHFlux::DrawFluxProjectionX(Option_t *opt) const
     458{
     459    const Int_t nVar = fHFlux.GetNbinsY();
     460
     461    for (int n=1; n<=nVar; n++)
     462    {
     463        TString strg0("Flux-");
     464
     465        TH1D &h = *((TH2D)fHFlux).ProjectionX(strg0+fVarname, n, n, "E");
     466
     467        TString strg1 = "Photon flux vs. E_{unfold} for ";
     468        TString strg2 = fVarname+"-bin #";
     469        strg2 += n;
     470
     471        new TCanvas(strg2, strg1+strg2);
     472        gPad->SetLogx();
     473        gPad->SetLogy();
     474
     475        TString name = fVarname+"bin_";
     476        name += n;
     477
     478        h.SetName(name);
     479        h.SetTitle(strg1+strg2);
     480        h.SetXTitle("E_{unfold} [GeV]");
     481        h.SetYTitle("photons / (s m^{2} GeV)");
     482        h.GetXaxis()->SetLabelOffset(-0.025);
     483        h.GetXaxis()->SetTitleOffset(1.1);
     484        h.GetXaxis()->SetNdivisions(1);
     485        h.GetYaxis()->SetTitleOffset(1.25);
     486        h.DrawCopy();
     487     }
     488}
     489
     490void MHFlux::DrawOrigProjectionX(Option_t *opt) const
     491{
     492    const Int_t nVar = fHOrig.GetNbinsY();
     493
     494    for (int n=1; n<=nVar; n++)
     495    {
     496        TString strg0 = "Orig-";
     497        strg0 += fVarname;
     498        strg0 += "_";
     499        strg0 += n;
     500
     501        TH1D &h = *((TH2D)fHOrig).ProjectionX(strg0, n, n, "E");
     502
     503        TString strg1("No.of photons vs. E-est for ");
     504        strg1 += fVarname+"-bin ";
     505        strg1 += n;
     506
     507        new TCanvas(strg0, strg1);
     508
     509        gPad->SetLogx();
     510        gPad->SetLogy();
     511
     512        h.SetName(strg0);
     513        h.SetTitle(strg1);
     514        h.SetXTitle("E_{est} [GeV]");
     515        h.GetXaxis()->SetLabelOffset(-0.025);
     516        h.GetXaxis()->SetTitleOffset(1.1);
     517        h.GetXaxis()->SetNdivisions(1);
     518        h.SetYTitle("No.of photons");
     519        h.DrawCopy();
     520    }
    591521}
    592522
    593523// -------------------------------------------------------------------------
    594524//
    595 // Draw copies of the histograms
     525//  Draw copies of the histograms
    596526//
    597527TObject *MHFlux::DrawClone(Option_t *opt) const
     
    622552// -------------------------------------------------------------------------
    623553//
    624 // Draw the histograms
     554//  Draw the histograms
    625555//
    626556void MHFlux::Draw(Option_t *opt)
     
    642572    gPad->Modified();
    643573    gPad->Update();
     574}
     575
     576Double_t MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
     577                       Double_t y1, Double_t y2, Double_t y3,
     578                       Double_t val)
     579{
     580    Double_t c0, c1, c2;
     581    Parab(x1, x2, x3, y1, y2, y3, &c0, &c1, &c2);
     582    return c0 + c1*val + c2*val*val;
    644583}
    645584
     
    660599        - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
    661600
    662     if (det == 0.0)
     601    if (det==0)
    663602    {
    664603        *a = 0;
     
    686625    *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1;
    687626
    688     //yt1 = *a + *b * x1 + *c * x1*x1;
    689     //yt2 = *a + *b * x2 + *c * x2*x2;
    690     //yt3 = *a + *b * x3 + *c * x3*x3;
    691 
    692     //*fLog << "x1 = " << x1 << ",  x2 = " << x2 << ",  x3 = " << x3 << endl;
    693     //*fLog << "y1 = " << y1 << ",  y2 = " << y2 << ",  y3 = " << y3 << endl;
    694     //*fLog << "yt1 = " << yt1 << ",  yt2 = " << yt2
    695     //    << ",  yt3 = " << yt3 << endl;
    696     //*fLog << "*a = " << *a << ",  *b = " << *b << ",  *c= " << *c
    697     //      << ",  *errflag = " << *errflag << endl;
    698 
    699627    return kTRUE;
    700628}
  • trunk/MagicSoft/Mars/mhist/MHFlux.h

    r1604 r1668  
    2020class TH2D;
    2121
     22class MHGamma;
     23class MHEffOnTime;
     24class MHThetabarTheta;
     25
    2226// base class MH is used because it defines "MakeDefCanvas"
    2327// if base class MH is used one has to define the member function Fill
     
    3539    // all these plots for different bins of the variable (Theta or time)
    3640
     41    void CalcEffCol(const TAxis &axex, TH2D &aeff, Int_t n, Double_t aeffbar[], Double_t daeffbar[]);
     42    Double_t ParabInterpolCos(const TAxis &axe, const TH2D *aeff, Int_t j, Int_t k3, Double_t val) const;
     43    void FindBins(const TAxis &axe, const Double_t bar, Int_t &k3, Int_t &k0) const;
     44    void CalcAbsGammaFlux(const TH1D &teff, const TH2D &fHAeff);
     45    Double_t ParabInterpolLog(const TAxis &axe, Int_t j,
     46                              Double_t y[], Double_t Ebar) const;
     47
    3748public:
    38     MHFlux(const TH2D &h2d, const Bool_t Draw,
    39           const TString varname, const TString unit);
     49    MHFlux(const TH2D &h2d, const TString varname, const TString unit);
     50    MHFlux(const MHGamma &h2d, const TString varname, const TString unit);
    4051
    4152    Bool_t Fill(const MParContainer *par);
    4253
    43     void Unfold(const Bool_t Draw);
     54    void Unfold();
    4455    void CalcFlux(const TH1D *teff, const TProfile *thetabar,
    45                   const TH2D *aeff, const Bool_t Draw);
     56                  const TH2D *aeff);
     57
     58    void CalcFlux(const MHEffOnTime &teff,
     59                  const MHThetabarTheta &thetabar,
     60                  const TH2D *aeff);
    4661
    4762    void Draw(Option_t *option="");
    4863    TObject *DrawClone(Option_t *option="") const;
     64    void DrawFluxProjectionX(Option_t *opt="") const;
     65    void DrawOrigProjectionX(Option_t *opt="") const;
    4966
    5067    const TH2D *GetHOrig()       { return &fHOrig; }
     
    5673                        double *a, double *b, double *c);
    5774
    58     ClassDef(MHFlux, 1) //2D-plots (original, unfolded, flux)
     75    static Double_t Parab(double x1, double x2, double x3,
     76                          double y1, double y2, double y3,
     77                          double val);
     78
     79    ClassDef(MHFlux, 0) //2D-plots (original, unfolded, flux)
    5980};
    6081
  • trunk/MagicSoft/Mars/mhist/MHGamma.cc

    r1422 r1668  
    3838#include <math.h>
    3939
    40 #include <TH1.h>
    4140#include <TH2.h>
    4241#include <TH3.h>
    4342
     43#include "MHAlphaEnergyTheta.h"
     44
    4445#include "MLog.h"
    4546#include "MLogManip.h"
     
    5253// Default Constructor.
    5354//
    54 MHGamma::MHGamma()
    55 {
     55MHGamma::MHGamma(const TString &name, const TString &title)
     56    : fHist(NULL), fProject(NULL)
     57{
     58    fName  = name.IsNull()  ? (TString)"MHGamma" : name;
     59    fTitle = title.IsNull() ? (TString)"3D-histogram of Alpha, E_est, Theta (Gamma sample)" : title;
    5660}
    5761
     
    6468  return kTRUE;
    6569}
     70
     71TH3D *MHGamma::Subtract(const MHAlphaEnergyTheta &h1, const MHAlphaEnergyTheta &h2)
     72{
     73    return Subtract(h1.GetHist(), h2.GetHist());
     74}
     75
     76TObject *MHGamma::DrawClone(Option_t *opt="") const
     77{
     78    DrawClone1();
     79    DrawClone2();
     80
     81    return NULL;
     82}
     83
     84void MHGamma::DrawClone1() const
     85{
     86    if (!fHist)
     87        return;
     88
     89    //
     90    // ------------- Part 1 ---------------------
     91    //
     92    TString titlex = fHist->GetXaxis()->GetTitle();
     93    TString titley = fHist->GetYaxis()->GetTitle();
     94    TString titlez = fHist->GetYaxis()->GetTitle();
     95
     96    TString canvtitle = "3D-plot "; //of ";
     97    /*
     98     canvtitle += titlex;
     99     canvtitle += ", ";
     100     canvtitle += titley;
     101     canvtitle += ", ";
     102     canvtitle += titlez+" ";
     103     */
     104    canvtitle += "for \'gamma\' sample";
     105
     106    TCanvas &c = *MakeDefCanvas("Alpha", canvtitle);
     107
     108    c.Divide(2, 2);
     109
     110    gROOT->SetSelectedPad(NULL);
     111
     112    TH1 *h;
     113
     114    c.cd(1);
     115    h = ((TH3D*)(fHist))->Project3D(fName+"_ex");
     116
     117    TString title= "Source-Antisource: ";
     118    h->SetTitle(title + titlex);
     119    h->SetXTitle(titlex);
     120    h->SetYTitle("Counts");
     121
     122    h->Draw();
     123    h->SetBit(kCanDelete);
     124
     125    c.cd(2);
     126    h = ((TH3D*)(fHist))->Project3D(fName+"_ey");
     127
     128    h->SetTitle(title + titley);
     129    h->SetXTitle(titley);
     130    h->SetYTitle("Counts");
     131
     132    h->Draw();
     133    h->SetBit(kCanDelete);
     134    gPad->SetLogx();
     135
     136    c.cd(3);
     137    h = ((TH3D*)(fHist))->Project3D(fName+"_ez");
     138
     139    h->SetTitle(title + titlez);
     140    h->SetXTitle(titlez);
     141    h->SetYTitle("Counts");
     142
     143    h->Draw();
     144    h->SetBit(kCanDelete);
     145
     146    c.cd(4);
     147    ((TH3D*)fHist)->DrawCopy();
     148
     149    c.Modified();
     150    c.Update();
     151}
     152
    66153
    67154// --------------------------------------------------------------------------
     
    70157//          fHist(gamma) = h1(source) - h2(antisource)
    71158//
    72 TH3D *MHGamma::Subtract(const TH3D *h1, const TH3D *h2,
    73                         const char *name, const char *title,
    74                         Bool_t Draw)
    75 {
    76     TH3D *fHist;
    77     fHist = new TH3D();
    78     fHist->SetName(name);
    79     fHist->SetTitle(title);
    80 
     159TH3D *MHGamma::Subtract(const TH3D *h1, const TH3D *h2)
     160{
     161    if (fHist)
     162        delete fHist;
     163
     164    fHist = new TH3D;
     165    fHist->SetName(fName);
     166    fHist->SetTitle(fTitle);
    81167    fHist->SetDirectory(NULL);
    82168
    83     // SetBinning((TH3D*)fHist, (TH3D*)h1);
    84169    SetBinning((TH1*)fHist, (TH1*)h1);
    85170
    86     TString strg1 =   (((TH1*)h1)->GetXaxis())->GetTitle();
    87     TString strg2 =   (((TH1*)h1)->GetYaxis())->GetTitle();
    88     TString strg3 =   (((TH1*)h1)->GetZaxis())->GetTitle();
    89     fHist->SetXTitle(strg1);
    90     fHist->SetYTitle(strg2);
    91     fHist->SetZTitle(strg3);
     171    fHist->SetXTitle((((TH1*)h1)->GetXaxis())->GetTitle());
     172    fHist->SetYTitle((((TH1*)h1)->GetYaxis())->GetTitle());
     173    fHist->SetZTitle((((TH1*)h1)->GetZaxis())->GetTitle());
    92174
    93175    fHist->Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME!
    94 
    95     //...........................................................
    96     // draw histogram
    97     if (Draw == kTRUE)
    98     {
    99       TString strg7 = "3D-plot of ";
    100       strg7 += strg1;
    101       strg7 += ",";
    102       strg7 += strg2;
    103       strg7 += ",";
    104       strg7 += strg3;
    105       strg7 += "  for \'gamma\' sample";
    106 
    107       TCanvas &c = *MakeDefCanvas("Alpha", strg7);
    108 
    109       c.Divide(2, 2);
    110 
    111       gROOT->SetSelectedPad(NULL);
    112 
    113       TH1 *h;
    114 
    115       c.cd(1);
    116       h = ((TH3D*)(fHist))->Project3D("ex");
    117 
    118       TString strg0 = "SRC-ASRC :    ";
    119       TString strg4 = strg0 + strg1;
    120       h->SetTitle(strg4);
    121       h->SetXTitle(strg1);
    122       h->SetYTitle("Counts");
    123 
    124       h->Draw();
    125       h->SetBit(kCanDelete);
    126 
    127       c.cd(2);
    128       h = ((TH3D*)(fHist))->Project3D("ey");
    129 
    130       TString strg5 = strg0 + strg2;
    131       h->SetTitle(strg5);
    132       h->SetXTitle(strg2);
    133       h->SetYTitle("Counts");
    134 
    135       h->Draw();
    136       h->SetBit(kCanDelete);
    137       gPad->SetLogx();
    138 
    139       c.cd(3);
    140       h = ((TH3D*)(fHist))->Project3D("ez");
    141 
    142       TString strg6 = strg0 + strg3;
    143       h->SetTitle(strg6);
    144       h->SetXTitle(strg3);
    145       h->SetYTitle("Counts");
    146 
    147       h->Draw();
    148       h->SetBit(kCanDelete);
    149 
    150       c.cd(4);
    151       ((TH3D*)fHist)->DrawCopy();
    152 
    153       c.Modified();
    154       c.Update();
    155     }
    156176
    157177    return fHist;
     
    162182// Integrate fHist(gamma) in the alpha range (lo, up)
    163183//
    164 TH2D *MHGamma::GetAlphaProjection(TH3D *fHist, Axis_t lo, Axis_t up,
    165                                   Bool_t Drawp)
     184TH2D *MHGamma::GetAlphaProjection(Axis_t lo, Axis_t up)
    166185{
    167186    if (up < lo)
     
    204223    axe.SetRange(ilo, iup);
    205224
    206     TH2D &h2D = *(TH2D *)fHist->Project3D("ezy");
    207 
    208     TString strg0 = "2D-plot of ";
    209     TString strg1 = (fHist->GetYaxis())->GetTitle();
    210     TString strg2 = (fHist->GetZaxis())->GetTitle();
    211     strg0 += strg2;
    212     strg0 += " vs. ";
    213     strg0 += strg1;
    214     h2D.SetTitle(strg0);
    215     h2D.SetXTitle(strg1);
    216     h2D.SetYTitle(strg2);
    217 
    218 
    219     //...........................................................
    220     // draw histogram
    221     if (Drawp == kTRUE)
    222     {
    223       char txt[100];
    224       TString strg3 = "No.of Gammas vs. ";
    225       strg3 += strg1;
    226       strg3 += " and ";
    227       strg3 += strg2;
    228       sprintf(txt, "   (%.1f < alpha < %.1f deg)", lo, up);
    229       strg3 += txt;
    230 
    231       TCanvas &c = *MakeDefCanvas("Gamma", strg3);
    232 
    233       c.Divide(2, 2);
    234 
    235       gROOT->SetSelectedPad(NULL);
    236 
    237       TH1 *h;
    238 
    239       c.cd(1);
    240       h = h2D.ProjectionX("xpro", -1, 9999, "E");
    241       TString strg0 = "No.of gammas : ";
    242       TString strg7 = strg0 + strg1;
    243       h->SetTitle(strg7);
    244       h->SetXTitle(strg1);
    245       h->SetYTitle("Counts");
    246 
    247       h->Draw();
    248       h->SetBit(kCanDelete);
    249       gPad->SetLogx();
    250 
    251       c.cd(2);
    252       h = h2D.ProjectionY("ypro", -1, 9999, "E");
    253       TString strg8 = strg0 + strg2;
    254       h->SetTitle(strg8);
    255       h->SetXTitle(strg2);
    256       h->SetYTitle("Counts");
    257 
    258       h->Draw();
    259       h->SetBit(kCanDelete);
    260 
    261       c.cd(3);
    262 
    263       h2D.DrawCopy();
    264       gPad->SetLogx();
    265 
    266       c.Modified();
    267       c.Update();
    268     }
    269     //...........................................................
    270 
    271     return &h2D;
    272 }
     225    fLo = lo;
     226    fHi = up;
     227
     228    if (fProject)
     229        delete fProject;
     230    fProject = (TH2D*)fHist->Project3D(fName+"_ezy");
     231
     232    const TString title = "2D-plot of ";
     233    const TString titley = fHist->GetYaxis()->GetTitle();
     234    const TString titlez = fHist->GetZaxis()->GetTitle();
     235
     236    fProject->SetTitle(title + titley + " vs. " + titlez);
     237    fProject->SetXTitle(titley);
     238    fProject->SetYTitle(titlez);
     239
     240    return fProject;
     241}
     242
     243void MHGamma::DrawClone2() const
     244{
     245    if (!fProject)
     246        return;
     247
     248    const TString titley = fHist->GetYaxis()->GetTitle();
     249    const TString titlez = fHist->GetZaxis()->GetTitle();
     250
     251    TString canvtitle = "No.of Gammas ";//vs. ";
     252    /*
     253     canvtitle += titley;
     254     canvtitle += " and ";
     255     canvtitle += titlez;
     256     */
     257    canvtitle += Form("(%.1f < alpha < %.1f deg)", fLo, fHi);
     258
     259    TCanvas &c = *MakeDefCanvas("Gamma", canvtitle);
     260
     261    c.Divide(2, 2);
     262
     263    gROOT->SetSelectedPad(NULL);
     264
     265    TH1 *h;
     266
     267    c.cd(1);
     268    h = fProject->ProjectionX(fName+"_xpro", -1, 9999, "E");
     269    TString title = "No.of gammas: ";
     270    h->SetTitle(title+titley);
     271    h->SetXTitle(titley);
     272    h->SetYTitle("Counts");
     273
     274    h->Draw();
     275    h->SetBit(kCanDelete);
     276    gPad->SetLogx();
     277
     278    c.cd(2);
     279    h = fProject->ProjectionY(fName+"_ypro", -1, 9999, "E");
     280    h->SetTitle(title+titlez);
     281    h->SetXTitle(titlez);
     282    h->SetYTitle("Counts");
     283
     284    h->Draw();
     285    h->SetBit(kCanDelete);
     286
     287    c.cd(3);
     288
     289    fProject->DrawCopy();
     290    gPad->SetLogx();
     291
     292    c.Modified();
     293    c.Update();
     294}
  • trunk/MagicSoft/Mars/mhist/MHGamma.h

    r1413 r1668  
    66#endif
    77
    8 #ifndef ROOT_TH3
    9 #include "TH3.h"
    10 #endif
    11 
    12 #ifndef ROOT_TH2
    13 #include "TH2.h"
    14 #endif
    15 
    168class TH2D;
    179class TH3D;
    1810
     11class MHAlphaEnergyTheta;
     12
    1913class MHGamma : public MH
    2014{
     15private:
     16    TH3D *fHist;    //!
     17    TH2D *fProject; //!
     18
     19    Axis_t fLo; //!
     20    Axis_t fHi; //!
     21
    2122public:
    22     MHGamma();
     23    MHGamma(const TString &name="", const TString &title="");
    2324
    2425    Bool_t Fill(const MParContainer *par);
    2526
    26     TH3D *Subtract(const TH3D *h1, const TH3D *h2,
    27                    const char *name, const char *title, Bool_t Draw=kFALSE);
     27    TH3D *Subtract(const TH3D *h1, const TH3D *h2);
    2828
    29     TH2D *GetAlphaProjection(TH3D *fHist, Axis_t lo, Axis_t up,
    30                              Bool_t Drawp=kFALSE);
     29    TH3D *Subtract(const MHAlphaEnergyTheta &h1, const MHAlphaEnergyTheta &h2);
    3130
     31    TH2D *GetAlphaProjection(Axis_t lo, Axis_t up);
    3232
    33     ClassDef(MHGamma, 1) // manipulation of alpha distributions
     33    TObject *DrawClone(Option_t *opt="") const;
     34    void DrawClone1() const;
     35    void DrawClone2() const;
     36
     37    const TH2D *GetProject() const { return fProject; }
     38
     39    ClassDef(MHGamma, 0) // manipulation of alpha distributions
    3440};
    3541
  • trunk/MagicSoft/Mars/mhist/MHHadronness.cc

    r1610 r1668  
    4545//    - black: acceptance of gammas (Ag) vs. the hadroness
    4646//    - red:   acceptance of non gammas (Ah) vs. the hadroness
    47 //    - blue:  2D distance of (acceptance_hadrons, acceptances_gammas)
    48 //             to optimum (0, 1)
    49 //             1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
    5047//  * Bottom Left Corner:
    5148//    Naive quality factor: Ag/sqrt(Ah)
     
    5653////////////////////////////////////////////////////////////////////////////
    5754#include "MHHadronness.h"
     55
     56/*
     57 // - blue:  2D distance of (acceptance_hadrons, acceptances_gammas)
     58 //          to optimum (0, 1):  1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
     59 */
    5860
    5961#include <TPad.h>
     
    113115    fQfac->SetTitle(" Naive Quality factor ");
    114116
    115     fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
    116     fMinDist->SetXTitle("Hadronness");
    117     fMinDist->SetYTitle("Distance");
     117    /*
     118     fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
     119     fMinDist->SetXTitle("Hadronness");
     120     fMinDist->SetYTitle("Distance");
     121     */
    118122
    119123    //fQfac->SetDirectory(NULL);
    120124    fGhness->SetDirectory(NULL);
    121125    fPhness->SetDirectory(NULL);
    122     fMinDist->SetDirectory(NULL);
     126    //fMinDist->SetDirectory(NULL);
    123127    fIntGhness->SetDirectory(NULL);
    124128    fIntPhness->SetDirectory(NULL);
     
    136140    delete fIntPhness;
    137141    delete fQfac;
    138     delete fMinDist;
     142    //    delete fMinDist;
    139143    delete fGraph;
    140144}
     
    266270            continue;
    267271
    268         fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
    269 
    270         Double_t val = ig/sqrt(ip);
     272        //fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
     273
     274        const Double_t val = ig/sqrt(ip);
    271275        fQfac->SetPoint(i, ig, val);
    272276
     
    360364
    361365    *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
    362     *fLog << " acc(hadron)  acc(gamma)" << endl <<endl;
    363 
    364     *fLog << "    0.005   "  << Form("%6.3f", GetGammaAcceptance(0.005)) << endl;
    365     *fLog << "    0.02    "  << Form("%6.3f", GetGammaAcceptance(0.02)) << endl;
    366     *fLog << "    0.05    "  << Form("%6.3f", GetGammaAcceptance(0.05)) << endl;
    367     *fLog << "    0.1     "  << Form("%6.3f", GetGammaAcceptance(0.1 )) << endl;
    368     *fLog << "    0.2     "  << Form("%6.3f", GetGammaAcceptance(0.2 )) << endl;
    369     *fLog << "    0.3     "  << Form("%6.3f", GetGammaAcceptance(0.3 )) << endl;
    370     *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << "       0.1  " << endl;
    371     *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << "       0.2  " << endl;
    372     *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << "       0.3  " << endl;
    373     *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << "       0.4  " << endl;
    374     *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << "       0.5  " << endl;
    375     *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << "       0.6  " << endl;
    376     *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << "       0.7  " << endl;
    377     *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << "       0.8  " << endl;
     366    *fLog << "acc(hadron)  acc(gamma)  acc(g)/acc(h)" << endl <<endl;
     367
     368    *fLog << "    0.005    " << Form("%6.3f", GetGammaAcceptance(0.005)) << "    " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
     369    *fLog << "    0.02     " << Form("%6.3f", GetGammaAcceptance(0.02))  << "    " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
     370    *fLog << "    0.05     " << Form("%6.3f", GetGammaAcceptance(0.05))  << "    " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
     371    *fLog << "    0.1      " << Form("%6.3f", GetGammaAcceptance(0.1 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
     372    *fLog << "    0.2      " << Form("%6.3f", GetGammaAcceptance(0.2 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
     373    *fLog << "    0.3      " << Form("%6.3f", GetGammaAcceptance(0.3 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
     374    *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << "        0.1  " << endl;
     375    *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << "        0.2  " << endl;
     376    *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << "        0.3  " << endl;
     377    *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << "        0.4  " << endl;
     378    *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << "        0.5  " << endl;
     379    *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << "        0.6  " << endl;
     380    *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << "        0.7  " << endl;
     381    *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << "        0.8  " << endl;
    378382    *fLog << endl;
    379383
    380     const Int_t h = fMinDist->GetMaximumBin();
    381     *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
    382     *fLog << "  Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
    383     *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
     384    /*
     385     const Int_t h = fMinDist->GetMaximumBin();
     386     *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
     387     *fLog << "  Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
     388     *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
     389     */
    384390
    385391    *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
     
    412418    Getiphness()->SetLineColor(kRed);
    413419    Getiphness()->DrawCopy("same");
    414     fMinDist->SetLineColor(kBlue);
    415     fMinDist->DrawCopy("same");
     420    /*
     421     fMinDist->SetLineColor(kBlue);
     422     fMinDist->DrawCopy("same");
     423     */
    416424    c.cd(3);
    417425    //GetQfac()->DrawCopy();
     
    449457        gPad->Update();
    450458    }
     459    /*
    451460    const Int_t h = fMinDist->GetMaximumBin();
    452461    TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
     
    455464    m->SetBit(kCanDelete);
    456465    m->Draw();
     466    */
    457467    return &c;
    458468}
     
    482492    Getiphness()->SetLineColor(kRed);
    483493    Getiphness()->Draw("same");
    484     fMinDist->SetLineColor(kBlue);
    485     fMinDist->Draw("same");
     494    /*
     495     fMinDist->SetLineColor(kBlue);
     496     fMinDist->Draw("same");
     497     */
    486498    gPad->cd(3);
    487499    //GetQfac()->Draw();
     
    516528        gPad->Update();
    517529    }
     530    /*
    518531    const Int_t h = fMinDist->GetMaximumBin();
    519532    TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
     
    522535    m->SetBit(kCanDelete);
    523536    m->Draw();
    524 }
     537    */
     538}
  • trunk/MagicSoft/Mars/mhist/MHHadronness.h

    r1574 r1668  
    1515{
    1616private:
    17     const MMcEvt *fMcEvt;            //!
    18     const MHadronness *fHadronness;    //!
     17    const MMcEvt *fMcEvt;           //!
     18    const MHadronness *fHadronness; //!
    1919
    20     TH1D* fPhness;        //-> Hadrons Hadronness
    21     TH1D* fGhness;        //-> Gammas Hadronness
    22     TH1D* fIntPhness;     //-> Hadrons Acceptance
    23     TH1D* fIntGhness;     //-> Gammas Acceptance
    24     TH1D* fMinDist;      //-> Minimum Distance to optimum acceptance
     20    TH1D* fPhness;    //-> Hadrons Hadronness
     21    TH1D* fGhness;    //-> Gammas Hadronness
     22    TH1D* fIntPhness; //-> Hadrons Acceptance
     23    TH1D* fIntGhness; //-> Gammas Acceptance
     24    //TH1D* fMinDist; //-> Minimum Distance to optimum acceptance
    2525
    26     TGraph *fQfac;        //-> Quality factor
    27     TGraph *fGraph;       //-> gamma acceptance vs. hadron acceptance
     26    TGraph *fQfac;    //-> Quality factor
     27    TGraph *fGraph;   //-> gamma acceptance vs. hadron acceptance
    2828
    2929public:
  • trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc

    r1642 r1668  
    112112// Fill data into the histogram which contains all showers
    113113//
    114 void MHMcCollectionArea::FillAll(Float_t energy, Float_t radius)
     114void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius)
    115115{
    116116    fHistAll->Fill(energy, radius);
     
    121121// Fill data into the histogram which contains the selected showers
    122122//
    123 void MHMcCollectionArea::FillSel(Float_t energy, Float_t radius)
     123void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius)
    124124{
    125125    fHistSel->Fill(energy, radius);
  • trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h

    r1610 r1668  
    2626    ~MHMcCollectionArea();
    2727
    28     void FillAll(Float_t energy, Float_t radius);
    29     void FillSel(Float_t energy, Float_t radius);
     28    void FillAll(Double_t energy, Double_t radius);
     29    void FillSel(Double_t energy, Double_t radius);
    3030
    3131    void DrawAll(Option_t *option="");
  • trunk/MagicSoft/Mars/mhist/MHMcEnergyMigration.h

    r1574 r1668  
    4040    TObject *DrawClone(Option_t *option="") const;
    4141
    42     ClassDef(MHMcEnergyMigration, 1) //3D-histogram   E-true E-est Theta
     42    ClassDef(MHMcEnergyMigration, 0) //3D-histogram   E-true E-est Theta
    4343
    4444};
  • trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc

    r1412 r1668  
    9999    gROOT->SetSelectedPad(NULL);
    100100
    101     ((TProfile*)&fHist)->DrawCopy(opt);
     101    ((TProfile)fHist).DrawCopy(opt);
    102102
    103103    c.Modified();
     
    128128Bool_t MHThetabarTheta::Fill(const MParContainer *par)
    129129{
    130     fHist.Fill(fMcEvt->GetTheta()*kRad2Deg, fMcEvt->GetTheta()*kRad2Deg);
     130    const Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
     131
     132    fHist.Fill(theta, theta);
    131133
    132134    return kTRUE;
    133135}
    134 
    135 
    136 
    137 
    138 
    139 
    140 
    141 
    142 
    143 
    144 
    145 
    146 
  • trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h

    r1574 r1668  
    1717{
    1818private:
    19     MTime  *fTime;   //!
    20     MMcEvt *fMcEvt;  //!
     19    MTime   *fTime;   //!
     20    MMcEvt  *fMcEvt;  //!
    2121
    22     TProfile    fHist;
     22    TProfile fHist;
    2323
    2424public:
     
    3636    TObject *DrawClone(Option_t *option="") const;
    3737
    38     ClassDef(MHThetabarTheta, 1) //Profile histogram Thetabar vs. time
     38    ClassDef(MHThetabarTheta, 0) //Profile histogram Thetabar vs. time
    3939
    4040};
  • trunk/MagicSoft/Mars/mhist/MHThetabarTime.h

    r1574 r1668  
    3939    TObject *DrawClone(Option_t *option="") const;
    4040
    41     ClassDef(MHThetabarTime, 1) //Profile histogram Thetabar vs. time
     41    ClassDef(MHThetabarTime, 0) //Profile histogram Thetabar vs. time
    4242
    4343};
  • trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.cc

    r1295 r1668  
    187187    gPad->Modified();
    188188    gPad->Update();
    189 
    190189}
    191190
     
    196195Bool_t MHTimeDiffTheta::Fill(const MParContainer *par)
    197196{
    198     const Int_t time = fTime->GetTimeLo();
    199 
    200     fHist.Fill(0.0001*(time-fLastTime), fMcEvt->GetTheta()*kRad2Deg);
     197    const Double_t time = 200e-9*fTime->GetTimeLo() + fTime->GetTimeHi();
     198
     199    fHist.Fill(time-fLastTime, fMcEvt->GetTelescopeTheta()*kRad2Deg);
    201200    fLastTime = time;
    202201
  • trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.h

    r1574 r1668  
    1919    MTime  *fTime;   //!
    2020    MMcEvt *fMcEvt;  //!
    21     Int_t  fLastTime;
     21    Double_t fLastTime;
    2222
    2323    TH2D    fHist;
     
    3737    TObject *DrawClone(Option_t *option="") const;
    3838
    39     ClassDef(MHTimeDiffTheta, 1) //2D-histogram  time-diff vs. Theta
     39    ClassDef(MHTimeDiffTheta, 0) //2D-histogram  time-diff vs. Theta
    4040};
    4141
  • trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.cc

    r1295 r1668  
    100100{
    101101
    102     TCanvas &c = *MakeDefCanvas("DiffTimeTime", "Distrib of \\Delta t, time");
     102    TCanvas &c = *MakeDefCanvas("DiffTimeTime", "Distrib of dt and t");
    103103
    104104    c.Divide(2, 2);
     
    186186Bool_t MHTimeDiffTime::Fill(const MParContainer *par)
    187187{
    188     const Int_t time = fTime->GetTimeLo();
     188    const Double_t time = 200e-9*fTime->GetTimeLo() + fTime->GetTimeHi();
    189189
    190     fHist.Fill(0.0001*(time-fLastTime), 0.0001*time);
    191 
     190    fHist.Fill(time-fLastTime, time);
    192191    fLastTime = time;
    193192
  • trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.h

    r1574 r1668  
    1717private:
    1818    MTime *fTime;   //!
    19     Int_t  fLastTime;
     19    Double_t  fLastTime;
    2020
    2121    TH2D   fHist;
     
    3535    TObject *DrawClone(Option_t *option="") const;
    3636
    37     ClassDef(MHTimeDiffTime, 1) //2D-histogram  time-diff vs. time
     37    ClassDef(MHTimeDiffTime, 0) //2D-histogram  time-diff vs. time
    3838};
    3939
  • trunk/MagicSoft/Mars/mmain/MAnalysis.cc

    r1656 r1668  
    100100}
    101101
    102 MAnalysis::MAnalysis(const TGWindow *main, const TGWindow *p,
     102MAnalysis::MAnalysis(/*const TGWindow *main,*/ const TGWindow *p,
    103103                     const UInt_t w, const UInt_t h)
    104 : MBrowser(main, p, w, h)
     104: MBrowser(/*main,*/ p, w, h)
    105105{
    106106    AddButtons();
  • trunk/MagicSoft/Mars/mmain/MAnalysis.h

    r1076 r1668  
    2828
    2929public:
    30     MAnalysis(const TGWindow *main=NULL, const TGWindow *p=NULL,
     30    MAnalysis(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
    3131              const UInt_t w=500, const UInt_t h=500) ;
    3232
  • trunk/MagicSoft/Mars/mmain/MBrowser.cc

    r1172 r1668  
    275275}
    276276
    277 MBrowser::MBrowser(const TGWindow *main, const TGWindow *p,
    278                      const UInt_t w, const UInt_t h)
    279 : TGTransientFrame(p?p:gClient->GetRoot(),
    280                    main?main:gClient->GetRoot(), w, h)
     277MBrowser::MBrowser(/*const TGWindow *main,*/ const TGWindow *p,
     278                   const UInt_t w, const UInt_t h)
     279//: TGTransientFrame(p?p:gClient->GetRoot(),
     280//                   main?main:gClient->GetRoot(), w, h)
     281: TGMainFrame(p?p:gClient->GetRoot(), w, h)
    281282{
    282283    fList = new MGList;
     
    291292    fList->Add(linesep);
    292293    AddFrame(linesep, laylinesep);
    293 
    294294
    295295    //
  • trunk/MagicSoft/Mars/mmain/MBrowser.h

    r1172 r1668  
    1919class TGFileContainer;
    2020
    21 class MBrowser : public TGTransientFrame
     21class MBrowser : public TGMainFrame/*TGTransientFrame*/
    2222{
    2323private:
     
    6262
    6363 public:
    64      MBrowser(const TGWindow *main=NULL, const TGWindow *p=NULL,
    65                const UInt_t w=500, const UInt_t h=500) ;
     64     MBrowser(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
     65              const UInt_t w=500, const UInt_t h=500) ;
    6666
    6767     virtual ~MBrowser();
  • trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc

    r1086 r1668  
    4343//  'new MCameraDisplay()'
    4444//
    45 MCameraDisplay::MCameraDisplay(const TGWindow *main, const TGWindow *p,
    46                        const UInt_t w, const UInt_t h)
    47 : MBrowser(main, p, w, h)
     45MCameraDisplay::MCameraDisplay(/*const TGWindow *main,*/ const TGWindow *p,
     46                               const UInt_t w, const UInt_t h)
     47: MBrowser(/*main,*/ p, w, h)
    4848{
    4949    TGTextButton *disp = new TGTextButton(fTop2, "Display Events", kButDisplay);
     
    9191        case kButDisplay:
    9292            new MGCamDisplay(fInputFile,
    93                              fClient->GetRoot(), this, 600, 500);
     93                             fClient->GetRoot(), /*this,*/ 600, 500);
    9494            return kTRUE;
    9595        }
  • trunk/MagicSoft/Mars/mmain/MCameraDisplay.h

    r1016 r1668  
    99{
    1010public:
    11     MCameraDisplay(const TGWindow *main=NULL, const TGWindow *p=NULL,
     11    MCameraDisplay(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
    1212                   const UInt_t w=500, const UInt_t h=500) ;
    1313
  • trunk/MagicSoft/Mars/mmain/MDataCheck.cc

    r1652 r1668  
    5454//  'new MDataCheck()'
    5555//
    56 MDataCheck::MDataCheck(const TGWindow *main, const TGWindow *p,
     56MDataCheck::MDataCheck(/*const TGWindow *main,*/ const TGWindow *p,
    5757                       const UInt_t w, const UInt_t h)
    58 : MBrowser(main, p, w, h)
     58: MBrowser(/*main,*/ p, w, h)
    5959{
    6060    TGTextButton *pedadc = new TGTextButton(fTop1, "ADC Spectra of Pedestals", kButPedAdc);
  • trunk/MagicSoft/Mars/mmain/MDataCheck.h

    r1652 r1668  
    1616
    1717public:
    18     MDataCheck(const TGWindow *main=NULL, const TGWindow *p=NULL,
     18    MDataCheck(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
    1919               const UInt_t w=500, const UInt_t h=500) ;
    2020
  • trunk/MagicSoft/Mars/mmain/MEvtDisp.cc

    r1086 r1668  
    3939};
    4040
    41 MEvtDisp::MEvtDisp(const TGWindow *main, const TGWindow *p,
     41MEvtDisp::MEvtDisp(/*const TGWindow *main,*/ const TGWindow *p,
    4242                   const UInt_t w, const UInt_t h)
    43 : MBrowser(main, p, w, h)
     43: MBrowser(/*main,*/ p, w, h)
    4444{
    4545    TGTextButton *fadcevt = new TGTextButton(fTop1, "FADC Disp for Events",  kButDispEvt);
     
    9696        case kButDispEvt:
    9797            new MGFadcDisp(fInputFile, "Events",
    98                            fClient->GetRoot(), this, 600, 500);
     98                           fClient->GetRoot(), /*this,*/ 600, 500);
    9999            return kTRUE;
    100100
    101101        case kButDispPedestal:
    102102            new MGFadcDisp(fInputFile , "PedEvts",
    103                            fClient->GetRoot(), this, 600, 500);
     103                           fClient->GetRoot(), /*this,*/ 600, 500);
    104104            return kTRUE;
    105105
    106106        case kButDispCalibration:
    107107            new MGFadcDisp(fInputFile , "CalEvts",
    108                            fClient->GetRoot(), this, 600, 500);
     108                           fClient->GetRoot(), /*this,*/ 600, 500);
    109109            return kTRUE;
    110110        }
  • trunk/MagicSoft/Mars/mmain/MEvtDisp.h

    r1016 r1668  
    99{
    1010public:
    11     MEvtDisp(const TGWindow *main=NULL, const TGWindow *p=NULL,
     11    MEvtDisp(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
    1212             const UInt_t w=500, const UInt_t h=500) ;
    1313
  • trunk/MagicSoft/Mars/mmain/MMars.cc

    r1644 r1668  
    277277
    278278            case kButEvtDisplay:
    279                 new MEvtDisp(this);
     279                new MEvtDisp(/*this*/);
    280280                return kTRUE;
    281281
    282282            case kButDataCheck:
    283                 new MDataCheck(this);
     283                new MDataCheck(/*this*/);
    284284                return kTRUE;
    285285
    286286            case kButAnalysis:
    287                 new MAnalysis(this);
     287                new MAnalysis(/*this*/);
    288288                return kTRUE;
    289289
    290290            case kButMonteCarlo:
    291                 new MMonteCarlo(this);
     291                new MMonteCarlo(/*this*/);
    292292                return kTRUE;
    293293
    294294            case kButCameraDisplay:
    295                 new MCameraDisplay(this);
     295                new MCameraDisplay(/*this*/);
    296296                return kTRUE;
    297297
  • trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc

    r1646 r1668  
    162162}
    163163
    164 MMonteCarlo::MMonteCarlo(const TGWindow *main, const TGWindow *p,
     164MMonteCarlo::MMonteCarlo(/*const TGWindow *main,*/ const TGWindow *p,
    165165                         const UInt_t w, const UInt_t h)
    166 : MBrowser(main, p, w, h)
     166: MBrowser(/*main,*/ p, w, h)
    167167{
    168168    AddButtons();
  • trunk/MagicSoft/Mars/mmain/MMonteCarlo.h

    r1030 r1668  
    2828
    2929public:
    30     MMonteCarlo(const TGWindow *main=NULL, const TGWindow *p=NULL,
     30    MMonteCarlo(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
    3131                const UInt_t w=500, const UInt_t h=500);
    3232
Note: See TracChangeset for help on using the changeset viewer.