Changeset 1668 for trunk

Timestamp:

11/25/02 09:12:47 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mbase/MEvtLoop.cc (modified) (1 diff)
• mbase/MFilter.h (modified) (1 diff)
• meventdisp/MGCamDisplay.cc (modified) (1 diff)
• meventdisp/MGCamDisplay.h (modified) (1 diff)
• meventdisp/MGEvtDisplay.cc (modified) (4 diffs)
• meventdisp/MGEvtDisplay.h (modified) (2 diffs)
• meventdisp/MGFadcDisp.cc (modified) (1 diff)
• meventdisp/MGFadcDisp.h (modified) (1 diff)
• mfileio/MCT1ReadAscii.cc (modified) (1 diff)
• mfileio/MCT1ReadPreProc.cc (modified) (10 diffs)
• mfileio/MCT1ReadPreProc.h (modified) (4 diffs)
• mfileio/MRead.cc (modified) (1 diff)
• mfileio/MRead.h (modified) (1 diff)
• mfileio/MReadMarsFile.cc (modified) (2 diffs)
• mfileio/MReadMarsFile.h (modified) (1 diff)
• mfileio/MReadTree.cc (modified) (10 diffs)
• mfileio/MReadTree.h (modified) (1 diff)
• mfilter/MF.cc (modified) (4 diffs)
• mfilter/MFDataChain.cc (modified) (1 diff)
• mfilter/MFEventSelector.cc (modified) (1 diff)
• mfilter/MFEventSelector.h (modified) (2 diffs)
• mfilter/MFilterList.cc (modified) (3 diffs)
• mhist/MBinning.cc (modified) (4 diffs)
• mhist/MBinning.h (modified) (3 diffs)
• mhist/MFillH.cc (modified) (1 diff)
• mhist/MH.cc (modified) (1 diff)
• mhist/MH.h (modified) (1 diff)
• mhist/MH3.h (modified) (1 diff)
• mhist/MHAlphaEnergyTheta.cc (modified) (10 diffs)
• mhist/MHAlphaEnergyTheta.h (modified) (1 diff)
• mhist/MHAlphaEnergyTime.h (modified) (1 diff)
• mhist/MHEffOnTime.cc (modified) (10 diffs)
• mhist/MHEffOnTime.h (modified) (3 diffs)
• mhist/MHEffOnTimeTheta.h (modified) (1 diff)
• mhist/MHEffOnTimeTime.h (modified) (1 diff)
• mhist/MHFlux.cc (modified) (14 diffs)
• mhist/MHFlux.h (modified) (3 diffs)
• mhist/MHGamma.cc (modified) (6 diffs)
• mhist/MHGamma.h (modified) (1 diff)
• mhist/MHHadronness.cc (modified) (12 diffs)
• mhist/MHHadronness.h (modified) (1 diff)
• mhist/MHMcCollectionArea.cc (modified) (2 diffs)
• mhist/MHMcCollectionArea.h (modified) (1 diff)
• mhist/MHMcEnergyMigration.h (modified) (1 diff)
• mhist/MHThetabarTheta.cc (modified) (2 diffs)
• mhist/MHThetabarTheta.h (modified) (2 diffs)
• mhist/MHThetabarTime.h (modified) (1 diff)
• mhist/MHTimeDiffTheta.cc (modified) (2 diffs)
• mhist/MHTimeDiffTheta.h (modified) (2 diffs)
• mhist/MHTimeDiffTime.cc (modified) (2 diffs)
• mhist/MHTimeDiffTime.h (modified) (2 diffs)
• mmain/MAnalysis.cc (modified) (1 diff)
• mmain/MAnalysis.h (modified) (1 diff)
• mmain/MBrowser.cc (modified) (2 diffs)
• mmain/MBrowser.h (modified) (2 diffs)
• mmain/MCameraDisplay.cc (modified) (2 diffs)
• mmain/MCameraDisplay.h (modified) (1 diff)
• mmain/MDataCheck.cc (modified) (1 diff)
• mmain/MDataCheck.h (modified) (1 diff)
• mmain/MEvtDisp.cc (modified) (2 diffs)
• mmain/MEvtDisp.h (modified) (1 diff)
• mmain/MMars.cc (modified) (1 diff)
• mmain/MMonteCarlo.cc (modified) (1 diff)
• mmain/MMonteCarlo.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                                -*-*- END -*-*-
+ 2002/11/25: Thomas Bretz
+
+   * macros/multidimdist2.C:
+     - renamed eventloops (instances had same names)
+     - fixed a type in PrintStatistics (the gamma statistics
+       were printed two times)
+
+   * mbase/MEvtLoop.cc:
+     - take the lowest value (entries from MRead or user input)
+       for the progress bar
+     - reset the progress bar
+
+   * mbase/MFilter.h:
+     - added 'private'
+
+   * meventdisp/MGCamDisplay.[h,cc], meventdisp/MGEvtDisplay.[h,cc],
+     meventdisp/MGFadcDisp.[h,cc], mmain/MMonteCarlo.[h,cc],
+     mmain/MAnalysis.[h,cc], mmain/MBrowser.[h,cc], 
+     mmain/MCameraDisplay.[h,cc], mmain/MDataCheck.[h,cc],
+     mmain/MEvtDisp.[h,cc], mmain/MMars.cc:
+     - changed from TTransientFrame to TMainFrame (with this I
+       get decorations, eg. Close Button)
+
+   * meventdisp/MGEvtDisplay.cc:
+     - Update the layout each time the fEvtInfo has changed
+
+   * mfileio/MCT1ReadAscii.cc, mfileio/MCT1ReadPreProc.cc:
+     - delete return of gSystem->ExpandPathName
+
+   * mfileio/MCT1ReadPreProc.[h,cc]:
+     - added output of Time
+     - added usage of Selector
+     - changed MTask basics to be private
+
+   * mfileio/MRead.[h,cc]:
+     - added comment about selector
+     - added Selector-stuff
+
+   * mfileio/MReadMarsFile.[h,cc], mfileio/MReadTree.[h,cc]:
+     - added 'entries' argument to AddFile
+
+   * mfileio/MReadTree.[h,cc]:
+     - added workaround for a root bug when a file doesn't exist
+     - changed AddFiles to use Add(TChain*)
+     - changed to use Selector
+
+   * mfilter/MF.cc:
+     - Set debug level to suppress output when MFDataChain is created
+
+   * mfilter/MFEventSelector.h:
+     - changed Pre//PostProcess to private
+     
+   * mfilter/MF.cc, mfilter/MFilterList.cc:
+     - changed the use of Pre//PostProcess to CallPre//PostProcess
+   
+   * mhist/MBinning.[h,cc]:
+     - changed comments
+     - added SetEdgesCos
+
+   * mhist/MFillH.[h,cc]:
+     - added GetBinCenterLog
+
+   * mhist/MH3.h:
+     - added default argument to GetHistByName
+
+   * mhist/MHAlphaEnergyTheta.[h,cc], mhist/MHAlphaEnergyTime.h,
+     mhist/MHEffOnTime.[h,cc], mhist/MHEffOnTimeTheta.h,
+     mhist/MHEffOnTimeTime.h, mhist/MHFlux.[h,cc], mhist/MHGamma.[h,cc],
+     mhist/MHMcEnergyMigration.h, mhist/MHThetabarTheta.[h,cc],
+     mhist/MHThetabarTime.h:
+     - changed the output
+     - changed the algorithms to be more modular (more usage of member
+       function)
+     - changed ClassDef to 0
+     - fixed some small bugs (access of TArray[n])
+
+   * mhist/MHHadronness.[h,cc]:
+     - removed shortest distance to (0,1) stuff
+
+   * mhist/MHMcCollectionArea.h:
+     - changed Fill to Double_t
+
+   * mhist/MHTimeDiffTheta.[h,cc], mhist/MHTimeDiffTime.[h,cc]:
+     - in a first draft changed to use 200ns timing of CT1
+     - changed ClassDef to 0
+
+
 
  2002/11/22: Thomas Bretz

trunk/MagicSoft/Mars/mbase/MEvtLoop.cc

@@ -208 +208 @@
     if (fProgress)
     {
+        fProgress->Reset();
+
+        Int_t entries = INT_MAX;
+
+#ifdef __MARS__
+        // limits.h
+        MRead *read = (MRead*)fTaskList->FindObject("MRead");
+        if (read && read->GetEntries()>0)
+            entries = read->GetEntries();
+#endif
+
         if (maxcnt>0)
-            fProgress->SetRange(0, maxcnt);
-#ifdef __MARS__
+            fProgress->SetRange(0, TMath::Min(maxcnt, entries));
         else
-        {
-            MRead *read = (MRead*)fTaskList->FindObject("MRead");
-            if (read && read->GetEntries()>0)
-                fProgress->SetRange(0, read->GetEntries());
-        }
-#endif
+            if (entries!=INT_MAX)
+                fProgress->SetRange(0, entries);
     }
 

trunk/MagicSoft/Mars/mbase/MFilter.h

@@ -10 +10 @@
 class MFilter : public MTask
 {
+private:
+    virtual Bool_t PreProcess(MParList *pList);
+    virtual Bool_t Process();
+    virtual Bool_t PostProcess();
+
 public:
     MFilter(const char *name=NULL, const char *title=NULL);
 
     virtual Bool_t IsExpressionTrue() const = 0;
-
-    virtual Bool_t PreProcess(MParList *pList);
-    virtual Bool_t Process();
-    virtual Bool_t PostProcess();
-
     virtual TString GetRule() const;
 

trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.cc

@@ -148 +148 @@
 //
 MGCamDisplay::MGCamDisplay(const char *filename,
-                           const TGWindow *p, const TGWindow *main,
+                           const TGWindow *p, /*const TGWindow *main,*/
                            UInt_t w, UInt_t h)
-: MGEvtDisplay(filename, "Events", p, main, w, h), fDisplay(NULL)
+: MGEvtDisplay(filename, "Events", p, /*main,*/ w, h), fDisplay(NULL)
 {
     //

trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.h

@@ -32 +32 @@
 public:
     MGCamDisplay(const char *filename,
-                 const TGWindow *p, const TGWindow *main,
+                 const TGWindow *p, /*const TGWindow *main,*/
                  UInt_t w, UInt_t h);
 

trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.cc

@@ -141 +141 @@
     fEvtInfo = new TGLabel(top2, new TGString(""));
     fList->Add(fEvtInfo);
-
     top2->AddFrame(fEvtInfo, laystd);
 
@@ -393 +392 @@
 
 MGEvtDisplay::MGEvtDisplay(const char *fname, const char *tname,
-                           const TGWindow *p, const TGWindow *main,
+                           const TGWindow *p, /*const TGWindow *main,*/
                            UInt_t w, UInt_t h)
-    : TGTransientFrame(p, main, w, h), fInitOk(kFALSE)
+//    : TGTransientFrame(p, main, w, h), fInitOk(kFALSE)
+: TGMainFrame(p, w, h), fInitOk(kFALSE)
 {
     //
@@ -472 +472 @@
         break;
     default:
-        txt += "Unknown Particle Id";
+        txt += "Unknown Particle Id#";
+        txt += evt->GetPartId();
     }
 
@@ -493 +494 @@
 
     fEvtInfo->SetText(txt);
+
+    //
+    // Seems to be necessary to newly layout the upper part to display
+    // the whole line of text
+    //
+    TGFrame &f = *(TGFrame*)fEvtInfo->GetParent()->GetParent();
+    f.Layout();
+    f.MapSubwindows();
 }
 

trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.h

@@ -21 +21 @@
 class MReadTree;
 
-class MGEvtDisplay : public TGTransientFrame
+class MGEvtDisplay : public TGMainFrame/*TGTransientFrame*/
 {
 private:
@@ -74 +74 @@
 public:
     MGEvtDisplay(const char *fname, const char *tname,
-                 const TGWindow *p, const TGWindow *main,
+                 const TGWindow *p, /*const TGWindow *main,*/
                  UInt_t w, UInt_t h);
 

trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.cc

@@ -113 +113 @@
 //
 MGFadcDisp::MGFadcDisp(const char *filename, const char *treename,
-                       const TGWindow *p, const TGWindow *main,
+                       const TGWindow *p, /*const TGWindow *main,*/
                        UInt_t w, UInt_t h)
-    : MGEvtDisplay(filename, treename, p, main, w, h)
+    : MGEvtDisplay(filename, treename, p, /*main,*/ w, h)
 {
     //

trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.h

@@ -23 +23 @@
 
     MGFadcDisp(const char *filename, const char *treename,
-               const TGWindow *p, const TGWindow *main, UInt_t w, UInt_t h);
+               const TGWindow *p, /*const TGWindow *main,*/ UInt_t w, UInt_t h);
 
     virtual Bool_t ProcessMessage(Long_t msg, Long_t parm1, Long_t parm2);

trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc

@@ -127 +127 @@
     const char *name = file->GetName();
 
-    fIn = new ifstream(gSystem->ExpandPathName(name));
+    const char *expname = gSystem->ExpandPathName(name);
+    fIn = new ifstream(expname);
+    delete [] expname;
 
     const Bool_t noexist = !(*fIn);

trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
-!   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
+!   Author(s): Thomas Bretz 11/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2001
+!   Copyright: MAGIC Software Development, 2000-2002
 !
 !
@@ -29 +28 @@
 //
 // Reads a output file of the CT1 preproc.
+//
+// Implements usage of a selector (see MRead) Use such a filter to skip
+// events before reading! But never use a filter which needs read data
+// as input...
 //
 //  Input Containers:
@@ -58 +61 @@
 #include "MLogManip.h"
 
+#include "MTime.h"
+#include "MFilter.h"
+
 #include "MParList.h"
 #include "MCerPhotEvt.h"
@@ -112 +118 @@
 void MCT1ReadPreProc::AddFile(const char *txt)
 {
-    const TString fname = gSystem->ExpandPathName(txt);
+    const char *name = gSystem->ExpandPathName(txt);
+    TString fname(name);
+    delete [] name;
 
     if (!CheckHeader(fname))
@@ -502 +510 @@
     // open the file which is the first one in the chain
     //
-    const TString name  = file->GetName();
-    const TString fname = gSystem->ExpandPathName(name);
+    const TString name = file->GetName();
+
+    const char *expname = gSystem->ExpandPathName(name);
+    const TString fname(expname);
+    delete [] expname;
 
     //
@@ -648 +659 @@
 
     //
+    //  look for the time class in the plist
+    //
+    fTime = (MTime*)pList->FindCreateObj("MTime");
+    if (!fTime)
+        return kFALSE;
+
+    //
     //  look for the pedestal class in the plist
     //
@@ -685 +703 @@
     fNumRuns       = 0;
 
-    return kTRUE;
+    return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE;
 }
 
@@ -695 +713 @@
 void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event)
 {
+    //  int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
+    //  int   isecfrac_200ns;     // fractional part of isecs_since_midday
+    fTime->SetTime(event.isecfrac_200ns, event.isecs_since_midday);
+    fTime->SetReadyToSave();
+
     /*
      --- USEFULL? NEEDED? ---
-     int   isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
-     int   isecfrac_200ns;     // fractional part of isecs_since_midday
      short snot_ok_flags;      // the bits in these two bytes are flags for additional information on the event: Everything OK =: all Bits = 0
 
@@ -862 +883 @@
             return kFALSE;
 
+    //
+    // Check for a selector. If one is given and returns kFALSE
+    // skip this event.
+    //
+    if (GetSelector())
+    {
+        //
+        // Make sure selector is processed
+        //
+        if (!GetSelector()->CallProcess())
+        {
+            *fLog << err << dbginf << "Processing Selector failed." << endl;
+            return kFALSE;
+        }
+
+        //
+        // Skip Event
+        //
+        if (!GetSelector()->IsExpressionTrue())
+        {
+            fIn->seekg(sizeof(struct eventrecord), ios::cur);
+
+            fNumEvents++;
+            fNumEventsInRun++;
+
+            return kCONTINUE;
+        }
+    }
+
     // event data to be read from the file
     struct eventrecord event;
@@ -890 +940 @@
     }
 
-    return kTRUE;
-}
+    return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE;
+}

trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h

@@ -7 +7 @@
 
 class TList;
+class MTime;
 class MMcEvt;
 class MMcTrig;
@@ -27 +28 @@
     MCerPhotEvt  *fNphot;   // the data container for all data.
     MPedestalCam *fPedest;  // ct1 pedestals
+    MTime        *fTime;    // event time
     MMcEvt       *fMcEvt;   // monte carlo data container for MC files
     MMcTrig      *fMcTrig;  // mc data container for trigger information
@@ -52 +54 @@
     void   ProcessEvent(const struct eventrecord &event);
 
+    Bool_t PreProcess(MParList *pList);
+    Bool_t Process();
+    Bool_t PostProcess();
+
 public:
     MCT1ReadPreProc(const char *filename=NULL,
@@ -61 +67 @@
     void AddFile(const char *fname);
 
-    Bool_t PreProcess(MParList *pList);
-    Bool_t Process();
-    Bool_t PostProcess();
-
     UInt_t GetEntries() { return fEntries; }
 

trunk/MagicSoft/Mars/mfileio/MRead.cc

@@ -29 +29 @@
 // Base class for all reading tasks                                        //
 //                                                                         //
+// You can set a selector. Depending on the impelementation in the derived //
+// class it can be used to skip events, if the filter return kFALSE.       //
+// Make sure that the selector (filter) doesn't need information which     //
+// doesn't exist before reading an event!                                  //
+//                                                                         //
 /////////////////////////////////////////////////////////////////////////////
 #include "MRead.h"

trunk/MagicSoft/Mars/mfileio/MRead.h

@@ -6 +6 @@
 #endif
 
+class MFilter;
+
 class MRead : public MTask
 {
+private:
+    MFilter *fSelector;
+
 public:
+    MRead() : fSelector(NULL) {}
 
     virtual UInt_t GetEntries() = 0;
+
+    void SetSelector(MFilter *f) { fSelector = f; }
+    MFilter *GetSelector() { return fSelector; }
 
     ClassDef(MRead, 0)  // Base class for a reading task

trunk/MagicSoft/Mars/mfileio/MReadMarsFile.cc

@@ -98 +98 @@
 //  chains -1 is returned, otherwise the number of files which were added.
 //
-Int_t MReadMarsFile::AddFile(const char *fname)
+Int_t MReadMarsFile::AddFile(const char *fname, Int_t entries)
 {
     //
@@ -107 +107 @@
     //
     Int_t n1 = fRun->AddFile(fname);
-    Int_t n2 = MReadTree::AddFile(fname);
+    Int_t n2 = MReadTree::AddFile(fname, entries);
 
     return n1 != n2 ? -1 : n1;

trunk/MagicSoft/Mars/mfileio/MReadMarsFile.h

@@ -24 +24 @@
     ~MReadMarsFile();
 
-    Int_t AddFile(const char *fname);
+    Int_t AddFile(const char *fname, Int_t entries=-1);
 
     ClassDef(MReadMarsFile, 1)  // Reads a tree from file(s) and the information from the 'RunHeader'-tree

trunk/MagicSoft/Mars/mfileio/MReadTree.cc

@@ -191 +191 @@
 //  AddFile returns the number of files added to the chain.
 //
-Int_t MReadTree::AddFile(const char *fname)
-{
+//  For more information see TChain::Add
+//
+Int_t MReadTree::AddFile(const char *fname, Int_t entries)
+{
+#if ROOT_VERSION_CODE < ROOT_VERSION(3,03,01)
+    //
+    // This is a workaround to get rid of crashed if the file doesn't
+    // exist due to non initialized TFile::fProcessIDs
+    //
+    //  (Code taken from TFile::TFile
+    //
+    const char *name;
+
+    TString newname;
+
+    if ((name = gSystem->ExpandPathName(fname)))
+    {
+        newname = name;
+        delete [] name;
+    }
+
+    if (newname.IsNull())
+    {
+        *fLog << err << dbginf << "Error expanding path " << fname << "." << endl;
+        return 0;
+    }
+
+    if (gSystem->AccessPathName(newname, kFileExists))
+    {
+        *fLog << err << "ERROR - File '" << fname << "' does not exist." << endl;
+        return 0;
+    }
+
+    fname = newname.Data();
+#endif
+
     //
     // FIXME! A check is missing whether the file already exists or not.
     //
-    //
-    // returns the number of file which were added
-    //
-
-    const Int_t numfiles = fChain->Add(fname);
+    const Int_t numfiles = fChain->Add(fname, entries);
 
     if (numfiles>0)
@@ -208 +238 @@
 }
 
+/*
+ // --------------------------------------------------------------------------
+ //
+ //
+ Int_t MReadTree::AddFile(const TChainElement &obj)
+ {
+     return AddFile(obj.GetTitle(), obj.GetEntries());
+ }
+*/
+
 // --------------------------------------------------------------------------
 //
@@ -216 +256 @@
 Int_t MReadTree::AddFiles(const MReadTree &read)
 {
-    Int_t rc = 0;
-
-    TIter Next(read.fChain->GetListOfFiles());
-    TObject *obj = NULL;
-    while ((obj=Next()))
-        rc += AddFile(obj->GetTitle());
+    const Int_t rc = fChain->Add(read.fChain);
 
     if (rc>0)
         SetBit(kChainWasChanged);
+
+    /*
+     Int_t rc = 0;
+
+     TIter Next(read.fChain->GetListOfFiles());
+     TObject *obj = NULL;
+     while ((obj=Next()))
+         rc += AddFile(*(TChainElement*)obj);
+    */
 
     return rc;
@@ -241 +285 @@
         return;
 
-    *fLog << inf << GetDescriptor() << ": Branch choosing method enabled (only enabled branches are read)." << endl;
+    *fLog << inf << GetDescriptor() << ": Branch choosing enabled (only enabled branches are read)." << endl;
     fChain->SetBranchStatus("*", kFALSE);
     fBranchChoosing = kTRUE;
@@ -255 +299 @@
 void MReadTree::EnableBranch(const char *name)
 {
+    if (fChain->GetListOfFiles()->GetEntries()==0)
+    {
+        *fLog << err << "Chain contains no file... Branch '";
+        *fLog << name << "' ignored." << endl;
+        return;
+    }
+
     EnableBranchChoosing();
 
@@ -449 +500 @@
 //  root environment) the branch is skipped and an error message is printed
 //  out.
+//  If a selector is specified it is preprocessed after the
+//  MReadTree::PreProcess
 //
 Bool_t MReadTree::PreProcess(MParList *pList)
@@ -576 +629 @@
     fChain->SetNotify(this);
 
-    return kTRUE;
+    return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE;
 }
 
@@ -626 +679 @@
 //  (Remark: The position can also be set by some member functions
 //  If the end of the file is reached the Eventloop is stopped.
+//  In case an event selector is given its value is checked before
+//  reading the event. If it returns kAFLSE the event is skipped.
 //
 #if ROOT_VERSION_CODE < ROOT_VERSION(3,02,06)
@@ -658 +713 @@
 #endif
 
-    Bool_t rc = fChain->GetEntry(fNumEntry++) != 0;
+    if (GetSelector())
+    {
+        //
+        // Make sure selector is processed
+        //
+        if (!GetSelector()->CallProcess())
+        {
+            *fLog << err << dbginf << "Processing Selector failed." << endl;
+            return kFALSE;
+        }
+
+        //
+        // Skip Event
+        //
+        if (!GetSelector()->IsExpressionTrue())
+        {
+            fNumEntry++;
+            return kCONTINUE;
+        }
+    }
+
+    const Bool_t rc = fChain->GetEntry(fNumEntry++) != 0;
 
     if (rc)
@@ -664 +740 @@
 
     return rc;
+}
+
+// --------------------------------------------------------------------------
+//
+//  If a selector is given the selector is post processed
+//
+Bool_t MReadTree::PostProcess()
+{
+    return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE;
 }
 

trunk/MagicSoft/Mars/mfileio/MReadTree.h

@@ -64 +64 @@
     virtual void   Print(Option_t *opt="") const;
 
-    virtual Int_t  AddFile(const char *fname);
+    virtual Int_t  AddFile(const char *fname, Int_t entries=-1);
     virtual Int_t  AddFiles(const MReadTree &read);
 
     virtual Bool_t PreProcess(MParList *pList);
     virtual Bool_t Process();
+    virtual Bool_t PostProcess();
 
     virtual Bool_t Notify();

trunk/MagicSoft/Mars/mfilter/MF.cc

@@ -244 +244 @@
     if (isrule)
     {
+        Int_t lvl = gLog.GetDebugLevel();
+        gLog.SetDebugLevel(1);
         newfilter = new MFDataChain(text.Data(), c, num);
         newfilter->SetName(Form("Chain%02d%c%f", level, c, num));
+        gLog.SetDebugLevel(lvl);
     }
     else
@@ -407 +410 @@
     }
 
-    if (!fFilter->PreProcess(plist))
+    if (!fFilter->CallPreProcess(plist))
     {
         *fLog << err << dbginf << "PreProcessing filters in ";
@@ -423 +426 @@
 Bool_t MF::Process()
 {
-    return fFilter->Process();
+    return fFilter->CallProcess();
 }
 
@@ -432 +435 @@
 Bool_t MF::PostProcess()
 {
-    return fFilter->PostProcess();
+    return fFilter->CallPostProcess();
 }
 

trunk/MagicSoft/Mars/mfilter/MFDataChain.cc

@@ -58 +58 @@
 //
 MFDataChain::MFDataChain(const char *member, const char type, const Double_t val,
-                           const char *name, const char *title)
+                         const char *name, const char *title)
     : fData(member), fValue(val)
 {

trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc

@@ -147 +147 @@
 Bool_t MFEventSelector::Process()
 {
-    Float_t evt = gRandom->Uniform();
+    const Float_t evt = gRandom->Uniform();
 
     if (fNumSelectEvts>0)

trunk/MagicSoft/Mars/mfilter/MFEventSelector.h

@@ -27 +27 @@
     void StreamPrimitive(ofstream &out) const;
 
+    Bool_t PreProcess(MParList *pList);
+    Bool_t Process();
+    Bool_t PostProcess();
+
     enum { kNumTotalFromFile = BIT(14) };
     /*
@@ -47 +51 @@
      */
 
-    Bool_t PreProcess(MParList *pList);
-    Bool_t Process();
-    Bool_t PostProcess();
-
     ClassDef(MFEventSelector, 0) // A Filter to select events from files
 };

trunk/MagicSoft/Mars/mfilter/MFilterList.cc

@@ -188 +188 @@
     //
     while ((filter=(MFilter*)Next()))
-        if (!filter->PreProcess(pList))
+        if (!filter->CallPreProcess(pList))
         {
             *fLog << err << "Error - Preprocessing Filter ";
@@ -212 +212 @@
     //
     while ((filter=(MFilter*)Next()))
-        if (!filter->Process())
+        if (!filter->CallProcess())
             return kFALSE;
 
@@ -232 +232 @@
     //
     while ((filter=(MFilter*)Next()))
-        if (!filter->PostProcess())
+        if (!filter->CallPostProcess())
             return kFALSE;
 

trunk/MagicSoft/Mars/mhist/MBinning.cc

@@ -64 +64 @@
 // --------------------------------------------------------------------------
 //
-// Specify the number of bins (not the number of edges), the lowest and
-// highest Edge.
+// Specify the number of bins <nbins> (not the number of edges), the
+// lowest <lo> and highest <up> Edge (of your histogram)
 //
 void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up)
@@ -79 +79 @@
 // --------------------------------------------------------------------------
 //
-// Specify the number of bins (not the number of edges), the lowest and
-// highest Edge.
+// Specify the number of bins <nbins> (not the number of edges), the
+// lowest <lo> and highest <up> Edge (of your histogram)
 //
 void MBinning::SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up)
@@ -92 +92 @@
 
     fType = kIsLogarithmic;
+}
+
+// --------------------------------------------------------------------------
+//
+// Specify the number of bins <nbins> (not the number of edges), the
+// lowest [deg] <lo> and highest [deg] <up> Edge (of your histogram)
+//
+void MBinning::SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up)
+{
+    // if (lo==0) ...
+    const Axis_t ld = lo/kRad2Deg;
+    const Axis_t ud = up/kRad2Deg;
+
+    const Double_t binsize = (cos(ld)-cos(ud))/nbins;
+    fEdges.Set(nbins+1);
+    for (int i=0; i<=nbins; i++)
+        fEdges[i] = (acos(binsize*i) + ld)*kRad2Deg;
+
+    fType = kIsCosinic;
 }
 
@@ -132 +151 @@
         return;
 
-    if (IsLinear() || IsLogarithmic())
+    if (IsLinear() || IsLogarithmic() || IsCosinic())
     {
         out << "   " << GetUniqueName() << ".SetEdges";
         if (IsLogarithmic())
             out << "Log";
+        if (IsCosinic())
+            out << "Cos";
         out << "(" << GetNumBins() << ", " << GetEdgeLo() << ", " << GetEdgeHi() << ");" << endl;
         return;

trunk/MagicSoft/Mars/mhist/MBinning.h

@@ -25 +25 @@
         kIsLinear,
         kIsLogarithmic,
+        kIsCosinic,
         kIsUserArray
     };
@@ -39 +40 @@
     void SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up);
     void SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up);
+    void SetEdgesCos(const Int_t nbins, const Axis_t lo, Axis_t up);
 
     // FIXME: ROOT workaround: "operator[] const" missing
@@ -51 +53 @@
     Bool_t IsLinear() const { return fType==kIsLinear; }
     Bool_t IsLogarithmic() const { return fType==kIsLogarithmic; }
+    Bool_t IsCosinic() const { return fType==kIsCosinic; }
     Bool_t IsDefault() const { return fType==kIsDefault; }
     Bool_t IsUserArray() const { return fType==kIsUserArray; }

trunk/MagicSoft/Mars/mhist/MFillH.cc

@@ -354 +354 @@
         {
             if (cls==name)
-                *fLog << inf << "Object '" << fHName << "' not found in parlist... trying to create one." << endl;
+                *fLog << inf << "Object '" << fHName << "' not found in parlist... creating." << endl;
             obj = pList->FindCreateObj(cls, name);
         }

trunk/MagicSoft/Mars/mhist/MH.cc

@@ -546 +546 @@
 // --------------------------------------------------------------------------
 //
+//  Returns the bin center in a logarithmic scale. If the given bin
+//  number is <1 it is set to 1. If it is =GetNbins() it is set to
+//  GetNbins()
+//
+Double_t MH::GetBinCenterLog(const TAxis &axe, Int_t nbin)
+{
+    if (nbin>axe.GetNbins())
+        nbin = axe.GetNbins();
+
+    if (nbin<1)
+        nbin = 1;
+
+    const Double_t lo = axe.GetBinLowEdge(nbin);
+    const Double_t hi = axe.GetBinUpEdge(nbin);
+
+    const Double_t val = log10(lo) + log10(hi);
+
+    return pow(10, val/2);
+}
+
+// --------------------------------------------------------------------------
+//
 // Draws a copy of the two given histograms. Uses title as the pad title.
 // Also layout the two statistic boxes and a legend.

trunk/MagicSoft/Mars/mhist/MH.h

@@ -54 +54 @@
     static void    ScaleAxis(TH1 *bins, Double_t fx=1, Double_t fy=1, Double_t fz=1);
 
+    static Double_t GetBinCenterLog(const TAxis &axe, Int_t nbin);
+
     static void DrawCopy(const TH1 &hist1, const TH1 &hist2, const TString title);
     static void Draw(TH1 &hist1, TH1 &hist2, const TString title);

trunk/MagicSoft/Mars/mhist/MH3.h

@@ -56 +56 @@
     const TH1 &GetHist() const { return *fHist; }
 
-    TH1 *GetHistByName(const TString name) { return fHist; }
+    TH1 *GetHistByName(const TString name="") { return fHist; }
 
     void SetColors() const;

trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc

@@ -64 +64 @@
     fHist.SetDirectory(NULL);
 
-    fHist.SetTitle("3D-plot of alpha, E-est, Theta");
+    fHist.SetTitle("3D-plot of alpha, E_{est}, Theta");
     fHist.SetXTitle("\\alpha [\\circ]");
-    fHist.SetYTitle("E-est [GeV]            ");
+    fHist.SetYTitle("E_{est} [GeV]");
     fHist.SetZTitle("\\Theta [\\circ]");
 }
@@ -79 +79 @@
    if (!fEnergy)
    {
-       *fLog << err << dbginf << "MHAlphaEnergyTheta : MEnergyEst not found... aborting." << endl;
+       *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl;
        return kFALSE;
    }
@@ -86 +86 @@
    if (!fMcEvt)
    {
-       *fLog << err << dbginf << "MHAlphaEnergyTheta : MMcEvt not found... aborting." << endl;
+       *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
        return kFALSE;
    }
@@ -95 +95 @@
    if (!binsenergy || !binsalphaflux || !binstheta)
    {
-       *fLog << err << dbginf << "MHAlphaEnergyTheta : At least one MBinning not found... aborting." << endl;
+       *fLog << err << dbginf << "At least one MBinning not found... aborting." << endl;
        return kFALSE;      
    }
@@ -103 +103 @@
    fHist.Sumw2(); 
 
-   fHist.Sumw2(); 
-
    return kTRUE;
 }
@@ -116 +114 @@
     MHillasSrc &hil = *(MHillasSrc*)par;
 
-    fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg);
+    fHist.Fill(fabs(hil.GetAlpha()), fEnergy->GetEnergy(),
+               fMcEvt->GetTelescopeTheta()*kRad2Deg);
 
     return kTRUE;
@@ -148 +147 @@
 
     h->SetTitle("Distribution of E-est [GeV]");
-    h->SetXTitle("E-est [GeV]            ");
+    h->SetXTitle("E_{est} [GeV]");
     h->SetYTitle("Counts");
 
@@ -186 +185 @@
 
     c.cd(1);
-    h = ((TH3*)(&fHist))->Project3D("expro");
+    h = ((TH3*)(&fHist))->Project3D(fName+"_expro");
 
     h->SetTitle("Distribution of \\alpha [\\circ]");
@@ -196 +195 @@
 
     c.cd(2);
-    h = ((TH3*)(&fHist))->Project3D("eypro");
+    h = ((TH3*)(&fHist))->Project3D(fName+"_eypro");
 
     h->SetTitle("Distribution of E-est [GeV]");
-    h->SetXTitle("E-est [GeV]            ");
+    h->SetXTitle("E_{est} [GeV]");
     h->SetYTitle("Counts");
 
@@ -207 +206 @@
 
     c.cd(3);
-    h = ((TH3*)(&fHist))->Project3D("ezpro");
+    h = ((TH3*)(&fHist))->Project3D(fName+"_ezpro");
 
     h->SetTitle("Distribution of \\Theta [\\circ]");

trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h

@@ -44 +44 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHAlphaEnergyTheta, 1) //3D-histogram in alpha, Energy and theta
+    ClassDef(MHAlphaEnergyTheta, 0) //3D-histogram in alpha, Energy and theta
 };
 

trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h

@@ -47 +47 @@
     TH1D *IntegrateEestTime(const char *title, Bool_t Draw);
    
-    ClassDef(MHAlphaEnergyTime, 1) //3D-histogram in alpha, Energy and time
+    ClassDef(MHAlphaEnergyTime, 0) //3D-histogram in alpha, Energy and time
 };
 

trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc

@@ -54 +54 @@
 #include "MLogManip.h"
 
+#include "MHTimeDiffTheta.h"
+
 ClassImp(MHEffOnTime);
 
@@ -68 +70 @@
     fUnit     = unit;
 
-    TString strg(fVarname);
-    strg += "  ";
-    strg += fUnit;
+    TString strg = fVarname + " " + fUnit;
 
     //
     //   set the name and title of this object
     //
-    TString strg1("MHEffOnTime-");
-    strg1 += fVarname;
-    fName  = strg1;
-    fTitle =  "1-D histogram of Eff On Time";
+    fName  = TString("MHEffOnTime-")+fVarname;
+    fTitle = "1-D histogram of Eff On Time";
 
     // effective on time versus Var
     fHEffOn.SetName("EffOn");
-    TString strg2("effective On Time vs. ");
-    strg2 += fVarname;
-    fHEffOn.SetTitle(strg2);
+    fHEffOn.SetTitle(TString("T_{on, eff} vs. ")+fVarname);
 
     fHEffOn.SetDirectory(NULL);
 
     fHEffOn.SetXTitle(strg);
-    fHEffOn.SetYTitle("effective ontime [s]");
+    fHEffOn.SetYTitle("T_{on, eff} [s]");
 
     // chi2-probability versus Var
     fHProb.SetName("Chi2-prob");
-    TString strg3("Chi2-prob of OnTimeFit vs. ");
+    TString strg3("\\chi^{2}-prob of OnTimeFit vs. ");
     strg3 += fVarname;
     fHProb.SetTitle(strg3);
@@ -100 +96 @@
 
     fHProb.SetXTitle(strg);
-    fHProb.SetYTitle("chi2-probability");
+    fHProb.SetYTitle("\\chi^{2}-probability");
 
     // lambda versus Var
     fHLambda.SetName("lambda");
-    TString strg4("lambda from OnTimeFit vs. ");
-    strg4 += fVarname;
-    fHLambda.SetTitle(strg4);
+    fHLambda.SetTitle(TString("\\lambda from OnTimeFit vs. ")+fVarname);
 
     fHLambda.SetDirectory(NULL);
@@ -116 +110 @@
     // Rdead is the fraction from real time lost by the dead time
     fHRdead.SetName("Rdead");
-    TString strg5("Rdead of  OnTimeFit vs. ");
-    strg5 += fVarname;
-    fHRdead.SetTitle(strg5);
+
+    fHRdead.SetTitle(TString("Rdead of OnTimeFit vs. ")+fVarname);
 
     fHRdead.SetDirectory(NULL);
@@ -127 +120 @@
 }
 
+// --------------------------------------------------------------------
+//
+//  determine range (yq[0], yq[1]) of time differences
+//  where fit should be performed;
+//  require a fraction >=min of all entries to lie below yq[0]
+//      and a fraction <=max of all entries to lie below yq[1];
+//  within the range (yq[0], yq[1]) there must be no empty bin;
+//
+void MHEffOnTime::GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const
+{
+#if ROOT_VERSION_CODE < ROOT_VERSION(3,02,07)
+    // WOrkaround for missing const qualifier of TH1::Integral
+    TH1 &h1 = (TH1&)h;
+
+    // choose pedestrian approach as long as GetQuantiles is
+    // not available
+    const Int_t    jbins = h1.GetNbinsX();
+    const Double_t Nm    = h1.Integral();
+
+    const Double_t xq[2] = { 0.15*Nm, 0.98*Nm };
+
+    yq[0] = yq[1] = h1.GetBinLowEdge(jbins+1);
+
+    for (int j=1; j<=jbins; j++)
+        if (h1.Integral(2, j) >= xq[0])
+        {
+            yq[0] = h1.GetBinLowEdge(j);
+            break;
+        }
+
+    for (int j=1; j<=jbins; j++)
+        if (h1.Integral(1, j) >= xq[1] || h1.GetBinContent(j)==0)
+        {
+            yq[1] = h1.GetBinLowEdge(j);
+            break;
+        }
+#else
+    // GetQuantiles doesn't seem to be available in root 3.01/06
+    Double_t xq[2] = { min, max };
+    h.GetQuantiles(2, yq, xq);
+#endif
+}
+
+void MHEffOnTime::DrawBin(TH1 &h, Int_t i) const
+{
+    TString strg1 = fVarname+"-bin #";
+    strg1 += i;
+
+    new TCanvas(strg1, strg1);
+
+    gPad->SetLogy();
+    gStyle->SetOptFit(1011);
+
+    TString name="Bin_";
+    name += i;
+
+    h.SetName(name);
+    h.SetXTitle("\\Delta t [s]");
+    h.SetYTitle("Counts");
+    h.DrawCopy();
+}
+
+Bool_t MHEffOnTime::CalcResults(const TF1 &func, Double_t Nm, Int_t i)
+{
+    const Double_t lambda = func.GetParameter(0);
+    const Double_t N0     = func.GetParameter(1);
+    const Double_t prob   = func.GetProb();
+    const Int_t    NDF    = func.GetNDF();
+
+    Double_t xmin, xmax;
+    ((TF1&)func).GetRange(xmin, xmax);
+
+    *fLog << inf;
+    *fLog << "Fitted bin #" << i << " from " << xmin << " to " << xmax;
+    *fLog << ",  got: lambda=" << lambda << "Hz N0=" << N0 << endl;
+
+    if (prob<=0.01)
+    {
+        *fLog << warn << "WARN - Fit bin#" << i << " gives:";
+        *fLog << " Chi^2-Probab(" << prob << ")<0.01";
+        *fLog << " NoFitPts=" << func.GetNumberFitPoints();
+        *fLog << " Chi^2=" << func.GetChisquare();
+    }
+
+    // was fit successful ?
+    if (NDF<=0 || /*prob<=0.001 ||*/ lambda<=0 || N0<=0)
+    {
+        *fLog << warn << dbginf << "Fit failed bin #" << i << ": ";
+        if (NDF<=0)
+            *fLog << " NDF(Number of Degrees of Freedom)=0";
+        if (lambda<=0)
+            *fLog << " Parameter#0(lambda)=0";
+        if (N0<=0)
+            *fLog << " Parameter#1(N0)=0";
+        *fLog << endl;
+
+        return kFALSE;
+    }
+
+    //
+    // -------------- start error calculation ----------------
+    //
+    Double_t emat[2][2];
+    gMinuit->mnemat(&emat[0][0], 2);
+
+    //
+    // Rdead : fraction of real time lost by the dead time
+    // kappa = number of observed events (Nm) divided by
+    //         the number of genuine events (N0)
+    // Teff  : effective on-time
+    //
+    const Double_t Teff  = Nm/lambda;
+    const Double_t kappa = Nm/N0;
+    const Double_t Rdead = 1.0 - kappa;
+
+    const Double_t dldl   = emat[0][0];
+    const Double_t dN0dN0 = emat[1][1];
+
+    const Double_t dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm);
+    const Double_t dl = sqrt(dldl);
+    const Double_t dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
+    //
+    // -------------- end error calculation ----------------
+    //
+
+    // the effective on time is Nm/lambda
+    fHEffOn.SetBinContent(i,  Teff);
+    fHEffOn.SetBinError  (i, dTeff);
+
+    // plot chi2-probability of fit
+    fHProb.SetBinContent(i, prob);
+
+    // lambda from fit
+    fHLambda.SetBinContent(i, lambda);
+    fHLambda.SetBinError  (i,     dl);
+
+    // Rdead from fit
+    fHRdead.SetBinContent(i, Rdead);
+    fHRdead.SetBinError  (i,dRdead);
+
+    return kTRUE;
+}
+
+void MHEffOnTime::ResetBin(Int_t i)
+{
+    fHEffOn.SetBinContent (i, 1.e-20);
+    fHProb.SetBinContent  (i, 1.e-20);
+    fHLambda.SetBinContent(i, 1.e-20);
+    fHRdead.SetBinContent (i, 1.e-20);
+}
+
 // -----------------------------------------------------------------------
 //
@@ -132 +276 @@
 // time differences
 //
-void MHEffOnTime::Calc(TH2D *hist, const Bool_t Draw)
-{
-     // Draw != 0   means the distribution of time differences should be drawn
-
+void MHEffOnTime::Calc(const TH2D *hist, const Bool_t draw)
+{
     // nbins = number of Var bins
     const Int_t nbins = hist->GetNbinsY();
 
-    // start  of 'for' loop ---------------------------
     for (int i=1; i<=nbins; i++)
     {
-      TString strg0("Calc-");
-      strg0 += fVarname;
-      TH1D &h = *hist->ProjectionX(strg0, i, i, "E");
+        TH1D &h = *((TH2D*)hist)->ProjectionX(TString("Calc-")+fVarname,
+                                              i, i, "E");
+        if (draw)
+            DrawBin(h, i);
+
+        ResetBin(i);
 
         // Nmdel = Nm * binwidth,  with Nm = number of observed events
+        const Double_t Nm    = h.Integral();
         const Double_t Nmdel = h.Integral("width");
-        const Double_t Nm    = h.Integral();
-
-        //...................................................
-        // determine range (yq[0], yq[1]) of time differences 
-        // where fit should be performed;
-        // require a fraction >=xq[0] of all entries to lie below yq[0]
-        //     and a fraction <=xq[1] of all entries to lie below yq[1];  
-        // within the range (yq[0], yq[1]) there must be no empty bin;
-        // choose pedestrian approach as long as GetQuantiles is not available
-
-        Double_t xq[2] = { 0.15, 0.95 };
+        if (Nm <= 0)
+        {
+            *fLog << warn << dbginf << "Nm<0 for bin #" << i << endl;
+            delete &h;
+            continue;
+        }
+
         Double_t yq[2];
-
-        // GetQuantiles doesn't seem to be available in root 3.01/06
-        // h.GetQuantiles(2,yq,xq);
-
-        const Int_t    jbins  = h.GetNbinsX();
-
-        fHEffOn.SetBinContent (i, 1.e-20);
-        fHProb.SetBinContent  (i, 1.e-20);
-        fHLambda.SetBinContent(i, 1.e-20);
-        fHRdead.SetBinContent (i, 1.e-20);
-
-        // start  of 'if' loop ---------------------------
-        if (Nm > 0.0)
-        {
-            // del is bin width in time difference
-            const Double_t del = Nmdel/Nm;
-
-            Double_t sum1 = 0.0;
-            yq[0]  = h.GetBinLowEdge(jbins+1);
-            for (int j=1; j<=jbins; j++)
-            {
-                if (sum1 >= xq[0]*Nm)
-                {
-                    yq[0] = h.GetBinLowEdge(j);
-                    break;
-                }
-                sum1 += h.GetBinContent(j);
-            }
-        
-            Double_t sum2 = 0.0;
-            yq[1] = h.GetBinLowEdge(jbins+1);
-            for (int j=1; j<=jbins; j++)
-            {
-                const Double_t content = h.GetBinContent(j);
-                sum2 += content;
-                if (sum2 >= xq[1]*Nm || content == 0.0)
-                {
-                    yq[1] = h.GetBinLowEdge(j);
-                    break;
-                }
-            }
-
-            //...................................................
-
-            // parameter 0 = lambda
-            // parameter 1 = N0           with N0 = ideal number of events
-            // parameter 2 = del (fixed)  with del = bin width of time difference
-
-            TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
-            func.FixParameter(2, del); 
-
-            func.SetParameter(0, 100); // [Hz]
-            func.SetParameter(1, Nm);
-
-            func.SetParLimits(0,   0,  1000);    // [Hz]
-            func.SetParLimits(1,   0, 10*Nm);
-
-            func.SetParName(0, "lambda");
-            func.SetParName(1, "N0");
-            func.SetParName(2, "del");
-
-            // options : 0  (=zero) do not plot the function
-            //           I  use integral of function in bin rather than value at bin center
-            //           R  use the range specified in the function range
-            //           Q  quiet mode
-            h.Fit("Poisson", "0IRQ");
-
-            const Double_t lambda = func.GetParameter(0);
-            const Double_t N0     = func.GetParameter(1);
-            // const Double_t chi2   = func.GetChisquare();
-            const Double_t Prob   = func.GetProb();
-            const Int_t    NDF    = func.GetNDF();
-
-            // was fit successful ?
-            if (NDF>0  &&  Prob>0.01  &&  lambda>0.0  &&  N0>0.0)
-            {
-                // start error calculation .....................
-
-                Double_t emat[2][2];
-                // Double_t eplus, eminus, eparab, gcc;
-                gMinuit->mnemat(&emat[0][0], 2);
-
-                // for (int m=0; m<2; m++)
-                // {
-                //   *fLog << "emat[" << m << "][n] = "; 
-                //   for (int n=0; n<2; n++)
-                //   {
-                //     *fLog << emat[m][n] << "  ";
-                //   }
-                //   *fLog << endl;
-                // }
-
-                // *fLog << "eplus, eminus, eparab, gcc :" << endl; 
-                // for (int k=0; k<2; k++)
-                // {
-                //   gMinuit->mnerrs(k, eplus, eminus, eparab, gcc);
-                //   *fLog << eplus << " " << eminus << " " << eparab << " "
-                //         << gcc << endl;
-                // }
-
-                // Rdead : fraction of real time lost by the dead time
-                // kappa = number of observed events (Nm) divided by
-                //         the number of genuine events (N0)
-                // Teff  : effective on-time
-
-                Double_t Rdead, kappa, Teff, dTeff, dldl, dl, dRdead;
-                Double_t dN0dN0;
-                //Double_t dldN0,lN0corr;
-                Teff  = Nm/lambda;
-                kappa = Nm/N0;
-                Rdead = 1.0 - kappa;
-
-                dldl   = emat[0][0];
-                dN0dN0 = emat[1][1];
-                // dldN0  = emat[0][1];
-
-                dTeff = Teff * sqrt(dldl/(lambda*lambda) + 1.0/Nm);
-                dl = sqrt(dldl);
-                dRdead = kappa * sqrt(dN0dN0/(N0*N0) + 1.0/Nm);
-
-                // lN0corr = dldN0 / sqrt(dldl * dN0dN0);
-                // *fLog << "MHEffOnTime : correlation between lambda and N0 = "
-                //       << lN0corr << endl;
-
-                // end error calculation .....................
-
-                // the effective on time is Nm/lambda
-                fHEffOn.SetBinContent(i,  Teff);
-                fHEffOn.SetBinError  (i, dTeff);
-
-                // plot chi2-probability of fit
-                fHProb.SetBinContent(i, Prob);
-
-                // lambda from fit
-                fHLambda.SetBinContent(i, lambda);
-                fHLambda.SetBinError  (i,     dl);
-
-                // Rdead from fit
-                fHRdead.SetBinContent(i, Rdead);
-                fHRdead.SetBinError  (i,dRdead);
-            }
-
-            //........................
-            // draw the distribution of time differences if requested
-            //
-            if (Draw == kTRUE)
-              {
-                char txt[100];
-                TString strg1(fVarname);
-                strg1 += "-bin %d";
-                sprintf(txt, strg1, i);
-                new TCanvas(txt, txt);
-                // TCanvas &c = *MakeDefCanvas(txt, txt);
-                // gROOT->SetSelectedPad(NULL);
-
-                gPad->SetLogy();
-                gStyle->SetOptFit(1011);
-                // gStyle->SetOptStat(1);
-
-                h.SetName(txt);
-                h.SetXTitle("time difference [s]");
-                h.SetYTitle("Counts");
-                h.DrawCopy();
-
-                func.SetRange(yq[0], yq[1]); // Range of Drawing
-                func.DrawCopy("same");
-
-                // c.Modified();
-                // c.Update();
-              }
-            //........................
-        }
-        // end  of 'if' loop ---------------------------
+        GetQuantiles(h, 0.15, 0.95, yq);
+
+        //
+        // Setup Poisson function for the fit:
+        // lambda [Hz], N0 = ideal no of evts, del = bin width of dt
+        //
+        TF1 func("Poisson", " [1]*[2] * [0] * exp(-[0] *x)", yq[0], yq[1]);
+        func.SetParNames("lambda", "N0", "del");
+
+        func.SetParameter(0, 1);
+        func.SetParameter(1, Nm);
+        func.FixParameter(2, Nmdel/Nm);
+
+        // For me (TB) it seems that setting parameter limits isn't necessary
+
+        // For a description of the options see TH1::Fit
+        h.Fit("Poisson", "0IRQ");
+
+        // Calc and fill results of fit into the histograms
+        const Bool_t rc = CalcResults(func, Nm, i);
+
+        // draw the distribution of time differences if requested
+        if (rc && draw)
+            func.DrawCopy("same");
 
         delete &h;
     }
-    // end  of 'for' loop ---------------------------
-    //    gStyle->SetOptStat(1111);
-
 }
 
@@ -347 +336 @@
 Bool_t MHEffOnTime::SetupFill(const MParList *plist)
 {
-    TString strg0("Binning");
-    strg0 += fVarname;
-
-    const MBinning* binsVar = (MBinning*)plist->FindObject(strg0);
+    TString bn = "Binning"+fVarname;
+
+    const MBinning* binsVar = (MBinning*)plist->FindObject(bn);
 
     if (!binsVar)
     {
-        *fLog << err << dbginf << strg0 
-              << " [MBinning] not found... aborting." << endl;
+        *fLog << err << dbginf << bn << " [MBinning] not found... aborting." << endl;
         return kFALSE;
     }
@@ -387 +374 @@
 TObject *MHEffOnTime::DrawClone(Option_t *opt) const
 {
-    TString strg0("EffOnTime");
-    strg0 += fVarname;    
-    TString strg1("Results from on time fit vs. ");
-    strg1 += fVarname;    
-
-    TCanvas &c = *MakeDefCanvas(strg0, strg1);
+    TCanvas &c = *MakeDefCanvas(TString("EffOnTime")+fVarname,
+                                TString("Results from on time fit vs. ")+fVarname);
     c.Divide(2, 2);
 
@@ -421 +404 @@
 void MHEffOnTime::Draw(Option_t *opt) 
 {
-    TString strg0("EffOnTime");
-    strg0 += fVarname;    
-    TString strg1("Results from on time fit vs. ");
-    strg1 += fVarname;    
-
     if (!gPad)
-        MakeDefCanvas(strg0, strg1);
+        MakeDefCanvas(TString("EffOnTime")+fVarname,
+                      TString("Results from on time fit vs. ")+fVarname);
 
     gPad->Divide(2,2);
@@ -447 +426 @@
 }
 
-
-
-
-
-
+void MHEffOnTime::Calc(const MHTimeDiffTheta &hist, const Bool_t Draw)
+{
+    Calc(hist.GetHist(), Draw);
+}
+

trunk/MagicSoft/Mars/mhist/MHEffOnTime.h

@@ -10 +10 @@
 
 class TH2D;
+class MHTimeDiffTheta;
 class MParList;
 
@@ -20 +21 @@
     TH1D fHRdead;
 
-    const char *fVarname;
-    const char *fUnit;
+    TString fVarname;
+    TString fUnit;
+
+    void ResetBin(Int_t i);
+    Bool_t CalcResults(const TF1 &func, Double_t Nm, Int_t i);
+    void DrawBin(TH1 &h, Int_t i) const;
+    void GetQuantiles(const TH1 &h, Double_t min, Double_t max, Double_t yq[2]) const;
 
 public:
@@ -32 +38 @@
     const TH1D *GetHist() const { return &fHEffOn; }
 
-    void Calc(TH2D *hist, const Bool_t Draw);
+    void Calc(const TH2D *hist, const Bool_t Draw);
+    void Calc(const MHTimeDiffTheta &hist, const Bool_t Draw);
 
     void Draw(Option_t *option="");
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHEffOnTime, 1) //1D-plot of Delta t vs. Var
+    ClassDef(MHEffOnTime, 0) //1D-plot of Delta t vs. Var
 };
 

trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h

@@ -35 +35 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHEffOnTimeTheta, 1) //1D-plot of Delta t vs. Theta
+    ClassDef(MHEffOnTimeTheta, 0) //1D-plot of Delta t vs. Theta
 };
 

trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h

@@ -35 +35 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHEffOnTimeTime, 1) //1D-plot of Delta t vs. time
+    ClassDef(MHEffOnTimeTime, 0) //1D-plot of Delta t vs. time
 };
 

trunk/MagicSoft/Mars/mhist/MHFlux.cc

@@ -42 +42 @@
 #include <TProfile.h>
 
-
 #include <TCanvas.h>
 
@@ -53 +52 @@
 #include "MLogManip.h"
 
+#include "MHThetabarTheta.h"
+#include "MHEffOnTime.h"
+#include "MHGamma.h"
+
 ClassImp(MHFlux);
 
+MHFlux::MHFlux(const MHGamma &hist, const TString varname, const TString unit)
+    : fHOrig(), fHUnfold(), fHFlux()
+{
+    const TH2D &h2d = *hist.GetProject();
+
+    if (varname.IsNull() || unit.IsNull())
+        *fLog << warn << dbginf << "varname or unit not defined" << endl;
+
+    fVarname = varname;
+    fUnit    = unit;
+
+    // char txt[100];
+
+    // original distribution of E-est for different bins 
+    //                       of the variable (Theta or time)
+    // sprintf(txt, "gammas vs. E-est and %s",varname);
+
+    ((TH2D&)h2d).Copy(fHOrig);
+
+    fHOrig.SetName("E_est");
+    fHOrig.SetTitle(TString("No.of gammas vs. E-est and ")+fVarname);
+
+    fHOrig.SetDirectory(NULL);
+    fHOrig.SetXTitle("E_{est} [GeV]");
+    fHOrig.SetYTitle(fVarname+fUnit);
+    //fHOrig.Sumw2();
+
+    SetBinning((TH2*)&fHOrig, (TH2*)&h2d);
+
+    fHOrig.Copy(fHUnfold);
+
+    // unfolded distribution of E-unfold for different bins 
+    //                       of the variable (Theta or time)
+    // sprintf(txt, "gammas vs. E-unfold and %s",varname);
+    fHUnfold.SetName("E-unfolded");
+    fHUnfold.SetTitle(TString("No.of gammas vs. E-unfold and ")+fVarname);
+
+    fHUnfold.SetDirectory(NULL);
+    fHUnfold.SetXTitle("E_{unfold} [GeV]");
+    fHUnfold.SetYTitle(fVarname+fUnit);
+    //fHUnfold.Sumw2();
+    
+    SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);
+
+
+    // absolute photon flux vs. E-unfold
+    //          for different bins of the variable (Theta or time)
+    //
+    // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
+    fHFlux.SetName("photon flux");
+    fHFlux.SetTitle(TString("Gamma flux [1/(s m^2 GeV) vs. E-unfold and ")+fVarname);
+
+    fHFlux.SetDirectory(NULL);
+    fHFlux.SetXTitle("E_{unfold} [GeV]");
+    fHFlux.SetYTitle(fVarname+fUnit);
+    fHFlux.Sumw2();
+
+    SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);
+}
 
 // --------------------------------------------------------------------------
@@ -61 +123 @@
 //                      and the units for the variable
 // 
-MHFlux::MHFlux(const TH2D &h2d,  const Bool_t Draw,
-               const TString varname, const TString unit)
+MHFlux::MHFlux(const TH2D &h2d, const TString varname, const TString unit)
     : fHOrig(), fHUnfold(), fHFlux()
 {
@@ -70 +131 @@
     fVarname = varname;
     fUnit    = unit;
-
-    TString strg(varname);
-    strg += unit;
 
     // char txt[100];
@@ -80 +138 @@
     // sprintf(txt, "gammas vs. E-est and %s",varname);
 
-    TString strg1 = "no.of gammas vs. E-est and ";
-    strg1 += varname;
-    
     ((TH2D&)h2d).Copy(fHOrig);
 
-    fHOrig.SetName("E-est");
-    fHOrig.SetTitle(strg1);
+    fHOrig.SetName("E_est");
+    fHOrig.SetTitle(TString("No.of gammas vs. E-est and ")+fVarname);
 
     fHOrig.SetDirectory(NULL);
-    fHOrig.SetXTitle("E-est [GeV]");
-    fHOrig.SetYTitle(strg);
-    fHOrig.Sumw2();
+    fHOrig.SetXTitle("E_{est} [GeV]");
+    fHOrig.SetYTitle(fVarname+fUnit);
+    //fHOrig.Sumw2();
+
+    // copy fHOrig into fHUnfold in case no unfolding is done
+    fHOrig.Copy(fHUnfold);
 
     SetBinning((TH2*)&fHOrig, (TH2*)&h2d);
@@ -99 +157 @@
     //                       of the variable (Theta or time)
     // sprintf(txt, "gammas vs. E-unfold and %s",varname);
-    TString strg2 = "no.of gammas vs. E-unfold and ";
-    strg2 += varname;
-
     fHUnfold.SetName("E-unfolded");
-    fHUnfold.SetTitle(strg2);
+    fHUnfold.SetTitle(TString("No.of gammas vs. E-unfold and ")+fVarname);
 
     fHUnfold.SetDirectory(NULL);
-    fHUnfold.SetXTitle("E-unfold [GeV]");
-    fHUnfold.SetYTitle(strg);
-    fHUnfold.Sumw2();
+    fHUnfold.SetXTitle("E_{unfold} [GeV]");
+    fHUnfold.SetYTitle(fVarname+fUnit);
+    //fHUnfold.Sumw2();
     
     SetBinning((TH2*)&fHUnfold, (TH2*)&fHOrig);
@@ -117 +172 @@
     //
     // sprintf(txt, "gamma flux [1/(s m2 GeV) vs. E-unfold and %s",varname);
-    TString strg3 = "gamma flux [1/(s m2 GeV) vs. E-unfold and ";
-    strg3 += varname;
-
     fHFlux.SetName("photon flux");
-    fHFlux.SetTitle(strg3);
+    fHFlux.SetTitle(TString("Gamma flux [1/(s m^{2} GeV)] vs. E-unfold and ")+fVarname);
 
     fHFlux.SetDirectory(NULL);
-    fHFlux.SetXTitle("E-unfold [GeV]");
-    fHFlux.SetYTitle(strg);
+    fHFlux.SetXTitle("E_{unfold} [GeV]");
+    fHFlux.SetYTitle(fVarname+fUnit);
     fHFlux.Sumw2();
 
     SetBinning((TH2*)&fHFlux, (TH2*)&fHUnfold);
-
-
-    // copy fHOrig into fHUnfold in case no unfolding is done
-    const Int_t nEunf   = fHUnfold.GetNbinsX();
-    const Int_t nVar    = fHUnfold.GetNbinsY();
-    for (int m=1; m<=nEunf; m++)
-    {
-      for (int n=1; n<=nVar; n++)
-      {
-        Double_t cont  = fHOrig.GetBinContent(m,n);
-        Double_t dcont = fHOrig.GetBinError(m,n);
-        fHUnfold.SetBinContent(m,n,cont);
-        fHUnfold.SetBinError(m,n,dcont);
-      }
-    }
-  //..............................................
-  // draw the No.of photons vs. E-est 
-  // for the individual bins of the variable Var
-
-  if (Draw == kTRUE)
-  {
-    const Int_t nVar = fHOrig.GetNbinsY();
-
-    for (int n=1; n<=nVar; n++)
-    {
-      TString strg0("Orig-");
-      strg0 += fVarname;
-      TH1D &h = *fHOrig.ProjectionX(strg0, n, n, "E");
-
-      strg0  = fVarname;
-      strg0 += "-bin ";
-      strg0 += n;
-
-      TString strg1("No.of photons vs. E-est for ");
-      strg1 += strg0;
-
-      new TCanvas(strg0, strg1);
-      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
-      // gROOT->SetSelectedPad(NULL);
-
-      gPad->SetLogx();
-
-      h.SetName(strg0);
-      h.SetTitle(strg1);
-      h.SetXTitle("E-est [GeV]");
-      h.SetYTitle("No.of photons");
-      h.DrawCopy();
-
-      // c.Modified();
-      // c.Update();
-    }
-  }
-  //........................
 }
 
@@ -200 +199 @@
 Bool_t MHFlux::Fill(const MParContainer *par)
 {
-   return kTRUE;
+    return kTRUE;
 }
 
@@ -208 +207 @@
 // Unfold the distribution in E-est
 //
-void MHFlux::Unfold(const Bool_t Draw)
-{
-  //..............................................
-  // draw the No.of photons vs. E-unfold 
-  // for the individual bins of the variable Var
-
-  if (Draw == kTRUE)
-  {
-    const Int_t nVar = fHUnfold.GetNbinsY();
-
-    for (int n=1; n<=nVar; n++)
-    {
-      TString strg0("Unfold-");
-      strg0 += fVarname;
-      TH1D &h = *fHUnfold.ProjectionX(strg0, n, n, "E");
-      strg0  = fVarname;
-      strg0 += "-bin ";
-      strg0 += n;
-
-      TString strg1("No.of photons vs. E-unfold for ");
-      strg1 += strg0;
-
-      new TCanvas(strg0, strg1);
-
-      // TCanvas &c = *MakeDefCanvas(txt0, strg1);
-      // gROOT->SetSelectedPad(NULL);
-
-      gPad->SetLogx();
-
-      h.SetName(strg0);
-      h.SetTitle(strg1);
-      h.SetXTitle("E-unfold [GeV]");
-      h.SetYTitle("No.of photons");
-      h.DrawCopy();
-
-      // c.Modified();
-      // c.Update();
-    }
-  }
-  //........................
-}
-
+void MHFlux::Unfold()
+{
+}
+
+void MHFlux::CalcFlux(const MHEffOnTime &teff, const MHThetabarTheta &thetabar,
+                      const TH2D *aeff)
+{
+    CalcFlux(teff.GetHist(), thetabar.GetHist(), aeff);
+}
+
+Double_t MHFlux::ParabInterpolLog(const TAxis &axe, Int_t j,
+                                  Double_t y[], Double_t Ebar) const
+{
+    const double t1 = log10(axe.GetBinLowEdge(j-1)) + log10(axe.GetBinUpEdge(j-1));
+    const double t2 = log10(axe.GetBinLowEdge(j))   + log10(axe.GetBinUpEdge(j));
+    const double t3 = log10(axe.GetBinLowEdge(j+1)) + log10(axe.GetBinUpEdge(j+1));
+
+    const Double_t lebar = log10(Ebar);
+
+    return Parab(t1/2, t2/2, t3/2, y[j-2], y[j-1], y[j], lebar);
+}
+
+// --------------------------------------------------------------------
+//
+//  determine bins for interpolation (k3 is the middle one) in bar.
+//  k0 denotes the bin from which the error is copied
+//
+void MHFlux::FindBins(const TAxis &axe, const Double_t bar, Int_t &k3, Int_t &k0) const
+{
+    const Int_t n = axe.GetNbins();
+
+    k3 = axe.FindFixBin(bar);
+    k0 = k3;
+
+    if (k3<2)
+    {
+        k3 = 2;
+        if (bar<axe.GetBinLowEdge(2))
+            k0 = 1;
+    }
+
+    if (k3>n-1)
+    {
+        k3 = n-1;
+        if (bar>axe.GetBinLowEdge(n))
+            k0 = n;
+    }
+
+    if (bar>=axe.GetBinLowEdge(1) && bar<=axe.GetBinUpEdge(n))
+        return;
+
+    *fLog << dbginf << "extrapolation: bar = " << bar;
+    *fLog << ", min =" << axe.GetBinLowEdge(1);
+    *fLog << ", max =" << axe.GetBinUpEdge(n) << endl;
+}
+
+Double_t MHFlux::ParabInterpolCos(const TAxis &axe, const TH2D *aeff, Int_t j, Int_t k3, Double_t val) const
+{
+    const double t1 = cos( axe.GetBinCenter (k3-1) );
+    const double t2 = cos( axe.GetBinCenter (k3)   );
+    const double t3 = cos( axe.GetBinCenter (k3+1) );
+
+    const double a1 = aeff->GetBinContent(j, k3-1);
+    const double a2 = aeff->GetBinContent(j, k3);
+    const double a3 = aeff->GetBinContent(j, k3+1);
+
+    return Parab(t1, t2, t3, a1, a2, a3, val);
+}
 
 // -------------------------------------------------------------------------
@@ -259 +284 @@
 //
 void MHFlux::CalcFlux(const TH1D *teff, const TProfile *thetabar,
-                      const TH2D *aeff, const Bool_t Draw)
-{
-  // Note that fHUnfold  has bins in Eunf and Var 
-  //           *teff     has bins in Var  (the same bins in Var as fHUnfold)
-  //           *thetabar has bins in Var  (the same bins in Var as fHUnfold)
-  //           *aeff     has bins in Etru and Theta
-  //                     (where in general the binning in Etru is different
-  //                      from the binning in Eunf)
-  // The variable Var may be 'time' or 'Theta'
-
-  // Draw = kTRUE means the differential flux vs E-unf should be drawn
-  //              for the individual bins of the variable Var
+                      const TH2D *aeff)
+{
+    //
+    // Note that fHUnfold  has bins in Eunf and Var
+    //           *teff     has bins in Var  (the same bins in Var as fHUnfold)
+    //           *thetabar has bins in Var  (the same bins in Var as fHUnfold)
+    //           *aeff     has bins in Etru and Theta
+    //                     (where in general the binning in Etru is different
+    //                      from the binning in Eunf)
+    // The variable Var may be 'time' or 'Theta'
 
     const TAxis &axex = *((TH2*)aeff)->GetXaxis();
     const TAxis &axey = *((TH2*)aeff)->GetYaxis();
 
-  //....................................
-  // define dummy histogram *aeff
-  ((TH1*)aeff)->Sumw2();
-  MBinning binsetru("BinningEtru");
-  binsetru.SetEdgesLog(10, 10, 1e3);
-
-  MBinning binsthetatru("BinningThetatru");
-  binsthetatru.SetEdges(7, -2.5, 32.5);
-  //SetBinning((TH1*)aeff, &binsetru, &binsthetatru);
-  SetBinning((TH2*)aeff, &binsetru, &binsthetatru);
-
-  const Int_t netru  = aeff->GetNbinsX();
-  const Int_t ntheta = aeff->GetNbinsY();
-
-  for (int j=1; j<=netru; j++)
-  {
-      for (int k=1; k<=ntheta; k++)
-      {
-          Double_t cont = 10000.0;
-          ((TH1*)aeff)->SetBinContent(j, k, cont);
-
-          Double_t dcont = 100.0;
-          ((TH1*)aeff)->SetBinError(j, k, dcont);
-      }
-  }
-  // *fLog << "Dummy aeff : netru =" << netru << ",  ntheta = " << ntheta << endl;
-  //....................................
-
-  // number of Eunf and Var bins   (histograms : fHUnfold, fHFlux)
-  const Int_t nEunf    = fHFlux.GetNbinsX();
-  const Int_t nVar     = fHFlux.GetNbinsY();
-
-  // number of Etru and Theta bins (histogram *aeff of collection area)
-
-  const Int_t nEtru    = aeff->GetNbinsX();
-  const Int_t nTheta   = aeff->GetNbinsY();
-
-  //*fLog << "nEunf =" << nEunf << ",  nVar = "   << nVar   << endl;
-  //*fLog << "nEtru =" << nEtru << ",  nTheta = " << nTheta << endl;
-
-  //...................................................................
-  // calculate effective collection area 
-  //    for the Eunf and Var bins of the histogram fHUnfold
-  //    from the histogram *aeff, which has bins in Etru and Theta
-  // the result is the histogram fHAeff
-  //
-
-  TH2D fHAeff;
-  fHAeff.Sumw2();
-  SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
-
-  Double_t *aeffbar  = new Double_t[nEtru];
-  Double_t *daeffbar = new Double_t[nEtru];
-
-  Double_t aeffEunfVar;
-  Double_t daeffEunfVar;
-
-  //------   start n loop   ------
-  for (int n=1; n<=nVar; n++)
-  {
-    Double_t Thetabar = thetabar->GetBinContent(n);
-    Double_t cosThetabar = cos(Thetabar);
-
-    // determine Theta bins (k1, k2, k3) for interpolation in Theta
-    // k0 denotes the Theta bin from whicvh the error is copied
-    Int_t k0=0;
-    Int_t k1=0;
-    Int_t k2=0;
-    Int_t k3=0;
-
-    for (int k=3; k<=nTheta; k++)
-    {
-        Double_t Thetalow = axey.GetBinLowEdge(k);
-        if (Thetabar < Thetalow)
+    if (axex.GetNbins()<3)
+    {
+        *fLog << err << "ERROR - Number of Energy bins <3 not implemented!" << endl;
+        return;
+    }
+
+    if (axey.GetNbins()<3)
+        *fLog << warn << "WARNING - Less than 3 theta-bins not supported very well!" << endl;
+
+    //
+    // calculate effective collection area
+    //    for the Eunf and Var bins of the histogram fHUnfold
+    //    from the histogram *aeff, which has bins in Etru and Theta
+    // the result is the histogram fHAeff
+    //
+    TH2D fHAeff;
+    SetBinning((TH2*)&fHAeff, (TH2*)&fHUnfold);
+    fHAeff.Sumw2();
+
+    //
+    // ------ start loops ------
+    //
+    const Int_t nEtru = aeff->GetNbinsX();
+
+    Double_t *aeffbar  = new Double_t[nEtru];
+    Double_t *daeffbar = new Double_t[nEtru];
+
+    const Int_t nVar = fHFlux.GetNbinsY();
+    for (int n=1; n<=nVar; n++)
+    {
+        const Double_t tbar = thetabar->GetBinContent(n);
+        const Double_t costbar = cos(tbar/kRad2Deg);
+
+        // determine bins for interpolation (k3, k0)
+        Int_t kv, ke;
+        FindBins(axey, tbar, kv, ke);
+
+        //
+        // calculate effective collection area at Theta = Thetabar
+        // by quadratic interpolation in cos(Theta);
+        // do this for each bin of Etru
+        //
+        for (int j=1; j<=nEtru; j++)
         {
-            k1 = k-2;
-            k2 = k-1;
-            k3 = k;
-            k0 = k2;
-            break;
-        }
-    }
-
-    if (k3 == 0)
-    {
-        k1 = nTheta-2;
-        k2 = nTheta-1;
-        k3 = nTheta;
-        k0 = k2;
-    }
-
-    if (Thetabar <  axey.GetBinLowEdge(2))
-      k0 = 1;
-    else
-        if (Thetabar >  axey.GetBinLowEdge(nTheta))
-            k0 = nTheta;
-
-    Double_t Thetamin = axey.GetBinLowEdge(1);
-    Double_t Thetamax = axey.GetBinLowEdge(nTheta+1);
-    if (Thetabar < Thetamin  ||  Thetabar > Thetamax)
-    {
-      *fLog << "MHFlux.cc : extrapolation in Theta; Thetabar = " << Thetabar
-            << ",  Thetamin =" << Thetamin
-            << ",  Thetamax =" << Thetamax << endl;
-    } 
-
-    //*fLog << "Var bin "   << n  << ":"  <<  endl;
-    //*fLog << "Thetabar= " << Thetabar 
-    //      << ",  k0= "    << k0
-    //      << ",  k1= "    << k1
-    //      << ",  k2= "    << k2
-    //      << ",  k3= "    << k3         <<  endl;
- 
-
-    // calculate effective collection area at Theta = Thetabar
-    // by quadratic interpolation in cos(Theta);
-    // do this for each bin of Etru 
-    //
-
-    for (int j=1; j<=nEtru; j++)
-    {
-        double c0 = 0;
-        double c1 = 0;
-        double c2 = 0;
-
-        const double t1 = cos( axey.GetBinCenter (k1) );
-        const double t2 = cos( axey.GetBinCenter (k2) );
-        const double t3 = cos( axey.GetBinCenter (k3) );
-
-        const double a1 = aeff->GetBinContent(j, k1);
-        const double a2 = aeff->GetBinContent(j, k2);
-        const double a3 = aeff->GetBinContent(j, k3);
-
-        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
-        aeffbar[j]  = c0 + c1*cosThetabar + c2*cosThetabar*cosThetabar;
-        daeffbar[j] = aeff->GetBinError(j,k0);
-
-        //*fLog << "Etru bin " << j <<  ":  tbar= " << Thetabar
-        //      << ",  abar= "  << aeffbar[j]
-        //      << ",  dabar= " << daeffbar[j] << endl;
-    }
-
-    //---   start m loop ---
-    // calculate effective collection area at (E = Ebar, Theta = Thetabar)
-    // by quadratic interpolation in log10(Etru)
-    // do this for each bin of Eunf
-    //
-    for (int m=1; m<=nEunf; m++)
-    {
-        Double_t log10Ebar = 0.5 * ( log10( fHUnfold.GetXaxis()->GetBinLowEdge(m)  )+
-                                     log10( fHUnfold.GetXaxis()->GetBinLowEdge(m+1)) );
-        Double_t Ebar = pow(10.0, log10Ebar);
-
-        // determine Etru bins (j1, j2, j3) for interpolation in E
-        // j0 denotes the Etru bin from which the error is copied
-        Int_t j0=0;
-        Int_t j1=0;
-        Int_t j2=0;
-        Int_t j3=0;
-
-        for (int j=3; j<=nEtru; j++)
-        {
-            Double_t Elow = axex.GetBinLowEdge(j);
-            if (Ebar < Elow)
+            if (axey.GetNbins()<3)
             {
-                j1 = j-2;
-                j2 = j-1;
-                j3 = j;
-                j0 = j2;
-                break;
+                // FIXME: Other interpolation?
+                aeffbar[j-1]  = aeff->GetBinContent(j, n);
+                daeffbar[j-1] = aeff->GetBinError(j, n);
+            }
+            else
+            {
+                aeffbar[j-1]  = ParabInterpolCos(axey, aeff, j, kv, costbar);
+                daeffbar[j-1] = aeff->GetBinError(j, ke);
             }
         }
 
-        if (j3 == 0)
+        //
+        // calculate effective collection area at (E = Ebar, Theta = Thetabar)
+        // by quadratic interpolation in log10(Etru)
+        // do this for each bin of Eunf
+        //
+        CalcEffCol(axex, fHAeff, n, aeffbar, daeffbar);
+    }
+
+    delete aeffbar;
+    delete daeffbar;
+
+    //
+    // now calculate the absolute gamma flux
+    //
+    CalcAbsGammaFlux(*teff, fHAeff);
+}
+
+// --------------------------------------------------------------------
+//
+//  calculate effective collection area at (E = Ebar, Theta = Thetabar)
+//  by quadratic interpolation in log10(Etru)
+//  do this for each bin of Eunf
+//
+void MHFlux::CalcEffCol(const TAxis &axex, TH2D &fHAeff, Int_t n, Double_t aeffbar[], Double_t daeffbar[])
+{
+    const Int_t nEunf = fHFlux.GetNbinsX();
+
+    const TAxis &unfx = *fHUnfold.GetXaxis();
+
+    for (int m=1; m<=nEunf; m++)
+    {
+        const Double_t Ebar = GetBinCenterLog(unfx, m);
+
+        Int_t j0, j3;
+        FindBins(axex, Ebar, j3, j0);
+
+        const Double_t v = ParabInterpolLog(axex, j3, aeffbar, Ebar);
+
+        fHAeff.SetBinContent(m,n, v);
+        fHAeff.SetBinError(m,n, daeffbar[j0-1]);
+    }
+}
+
+// --------------------------------------------------------------------
+//
+//  calculate the absolute gamma flux
+//
+void MHFlux::CalcAbsGammaFlux(const TH1D &teff, const TH2D &fHAeff)
+{
+    const Int_t nEunf = fHFlux.GetNbinsX();
+    const Int_t nVar  = fHFlux.GetNbinsY();
+
+    for (int m=1; m<=nEunf; m++)
+    {
+        const Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
+
+        for (int n=1; n<=nVar; n++)
         {
-            j1 = nEtru-2;
-            j2 = nEtru-1;
-            j3 = nEtru;
-            j0 = j2;
+            const Double_t Ngam  = fHUnfold.GetBinContent(m,n);
+            const Double_t Aeff  = fHAeff.GetBinContent(m,n);
+            const Double_t Effon = teff.GetBinContent(n);
+
+            const Double_t c1 = fHUnfold.GetBinError(m,n)/Ngam;
+            const Double_t c2 = teff.GetBinError(n)      /Effon;
+            const Double_t c3 = fHAeff.GetBinError(m,n)  /Aeff;
+
+            const Double_t cont  = Ngam / (DeltaE * Effon * Aeff);
+            const Double_t dcont = sqrt(c1*c1 + c2*c2 + c3*c3);
+
+            //
+            // Out of Range
+            //
+            const Bool_t oor = Ngam<=0 || DeltaE<=0 || Effon<=0 || Aeff<=0;
+
+            if (oor)
+                *fLog << warn << "MHFlux::CalcAbsGammaFlux(" << m << "," << n << ") ";
+
+            if (Ngam<=0)
+                *fLog << " Ngam=0";
+            if (DeltaE<=0)
+                *fLog << " DeltaE=0";
+            if (Effon<=0)
+                *fLog << " Effon=0";
+            if (Aeff<=0)
+                *fLog << " Aeff=0";
+
+            if (oor)
+                *fLog << endl;
+
+            fHFlux.SetBinContent(m,n, oor ? 1e-20 : cont);
+            fHFlux.SetBinError(m,n,   oor ? 1e-20 : dcont*cont);
         }
-
-        if (Ebar <  axex.GetBinLowEdge(2))
-            j0 = 1;
-        else
-            if (Ebar >  axex.GetBinLowEdge(nEtru))
-                j0 = nEtru;
-
-        Double_t Etrumin = axex.GetBinLowEdge(1);
-        Double_t Etrumax = axex.GetBinLowEdge(nEtru+1);
-        if (Ebar < Etrumin  ||  Ebar > Etrumax)
-        {
-            *fLog << "MHFlux.cc : extrapolation in Energy; Ebar = " << Ebar
-                << ",  Etrumin =" << Etrumin
-                << ",  Etrumax =" << Etrumax << endl;
-        }
-
-        //*fLog << "Var bin "   << n  << ":"  <<  endl;
-        //*fLog << "Ebar= " << Ebar
-        //      << ",  j1= "    << j1
-        //      << ",  j2= "    << j2
-        //      << ",  j3= "    << j3         <<  endl;
-
-
-        double c0=0.0;
-        double c1=0.0;
-        double c2=0.0;
-
-        const double t1 = 0.5 * ( log10( axex.GetBinLowEdge (j1)  )+
-                                  log10( axex.GetBinLowEdge (j1+1)) );
-        const double t2 = 0.5 * ( log10( axex.GetBinLowEdge (j2)  )+
-                                  log10( axex.GetBinLowEdge (j2+1)) );
-        const double t3 = 0.5 * ( log10( axex.GetBinLowEdge (j3)  )+
-                                  log10( axex.GetBinLowEdge (j3+1)) );
-
-        const double a1 = aeffbar[j1];
-        const double a2 = aeffbar[j2];
-        const double a3 = aeffbar[j3];
-
-        Parab(t1, t2, t3, a1, a2, a3, &c0, &c1, &c2);
-        aeffEunfVar  = c0 + c1*log10(Ebar) + c2*log10(Ebar)*log10(Ebar);
-        daeffEunfVar = daeffbar[j0];
-
-        //*fLog << "Eunf bin " << m     <<  ":  Ebar= " << Ebar
-        //      << ",  aeffEunfVar = "  << aeffEunfVar
-        //      << ",  daeffEunfVar = " << daeffEunfVar << endl;
-
-        fHAeff.SetBinContent(m,n,aeffEunfVar);
-        fHAeff.SetBinError(m,n,daeffEunfVar);
-    }
-    //---   end m loop ---
-  }
-  //------   end n loop   ------
-  delete aeffbar;
-
-  //...................................................................
-  // now calculate the absolute gamma flux
-  //
-  for (int m=1; m<=nEunf; m++)
-  {
-      Double_t DeltaE = fHFlux.GetXaxis()->GetBinWidth(m);
-
-      for (int n=1; n<=nVar; n++)
-      {
-          Double_t Ngam   = fHUnfold.GetBinContent(m,n);
-          Double_t dNgam  = fHUnfold.GetBinError(m,n);
-
-          Double_t Aeff   = fHAeff.GetBinContent(m,n);
-          Double_t dAeff  = fHAeff.GetBinError(m,n);
-
-          Double_t Effon  = teff->GetBinContent(n);
-          Double_t dEffon = teff->GetBinError(n);
-
-          Double_t Cont, dCont;
-          if (Ngam>0 && DeltaE>0 && Effon>0 && Aeff>0)
-          {
-              Cont  = Ngam / (DeltaE * Effon * Aeff);
-              dCont = Cont * sqrt( dNgam *dNgam  / (Ngam*Ngam) +
-                                   dEffon*dEffon / (Effon*Effon) +
-                                   dAeff *dAeff  / (Aeff*Aeff)  );
-          }
-          else
-          {
-              Cont  = 1.e-20;
-              dCont = 1.e-20;
-          }
-
-          fHFlux.SetBinContent(m,n,Cont);
-          fHFlux.SetBinError(m,n,dCont);
-
-          //*fLog << "Eunf bin "    << m      << ",  Var bin " << n
-          //      << ":  Ngam = "   << Ngam   << ",  Flux = "
-          //      << Cont  << ", dFlux = " << dCont << endl;
-          //*fLog << ",  DeltaE = " << DeltaE << ",  Effon = " << Effon
-          //      << ",  Aeff = "   << Aeff   << endl;
-      }
-  }
-
-  //..............................................
-  // draw the differential photon flux vs. E-unf 
-  // for the individual bins of the variable Var
-
-  if (Draw == kTRUE)
-  {
-      for (int n=1; n<=nVar; n++)
-      {
-          TString strg0("Flux-");
-          strg0 += fVarname;
-
-          TH1D &h = *fHFlux.ProjectionX(strg0, n, n, "E");
-
-          TString strg1("Photon flux vs. E-unfold for ");
-          TString strg2 = fVarname;
-
-          strg2 += "-bin ";
-          strg2 += n;
-
-          TString strg3 = strg1 + strg2;
-          new TCanvas(strg2, strg3);
-          // TCanvas &c = *MakeDefCanvas(txt, txt);
-          // gROOT->SetSelectedPad(NULL);
-
-          gPad->SetLogx();
-
-          h.SetName(strg2);
-          h.SetTitle(strg3);
-          h.SetXTitle("E-unfold [GeV]            ");
-          h.SetYTitle("photons / (s m2 GeV)");
-          h.DrawCopy();
-
-          // c.Modified();
-          // c.Update();
-      }
-  }
-  //........................
+    }
+}
+
+// --------------------------------------------------------------------
+//
+// draw the differential photon flux vs. E-unf
+// for the individual bins of the variable Var
+//
+void MHFlux::DrawFluxProjectionX(Option_t *opt) const
+{
+    const Int_t nVar = fHFlux.GetNbinsY();
+
+    for (int n=1; n<=nVar; n++)
+    {
+        TString strg0("Flux-");
+
+        TH1D &h = *((TH2D)fHFlux).ProjectionX(strg0+fVarname, n, n, "E");
+
+        TString strg1 = "Photon flux vs. E_{unfold} for ";
+        TString strg2 = fVarname+"-bin #";
+        strg2 += n;
+
+        new TCanvas(strg2, strg1+strg2);
+        gPad->SetLogx();
+        gPad->SetLogy();
+
+        TString name = fVarname+"bin_";
+        name += n;
+
+        h.SetName(name);
+        h.SetTitle(strg1+strg2);
+        h.SetXTitle("E_{unfold} [GeV]");
+        h.SetYTitle("photons / (s m^{2} GeV)");
+        h.GetXaxis()->SetLabelOffset(-0.025);
+        h.GetXaxis()->SetTitleOffset(1.1);
+        h.GetXaxis()->SetNdivisions(1);
+        h.GetYaxis()->SetTitleOffset(1.25);
+        h.DrawCopy();
+     }
+}
+
+void MHFlux::DrawOrigProjectionX(Option_t *opt) const
+{
+    const Int_t nVar = fHOrig.GetNbinsY();
+
+    for (int n=1; n<=nVar; n++)
+    {
+        TString strg0 = "Orig-";
+        strg0 += fVarname;
+        strg0 += "_";
+        strg0 += n;
+
+        TH1D &h = *((TH2D)fHOrig).ProjectionX(strg0, n, n, "E");
+
+        TString strg1("No.of photons vs. E-est for ");
+        strg1 += fVarname+"-bin ";
+        strg1 += n;
+
+        new TCanvas(strg0, strg1);
+
+        gPad->SetLogx();
+        gPad->SetLogy();
+
+        h.SetName(strg0);
+        h.SetTitle(strg1);
+        h.SetXTitle("E_{est} [GeV]");
+        h.GetXaxis()->SetLabelOffset(-0.025);
+        h.GetXaxis()->SetTitleOffset(1.1);
+        h.GetXaxis()->SetNdivisions(1);
+        h.SetYTitle("No.of photons");
+        h.DrawCopy();
+    }
 }
 
 // -------------------------------------------------------------------------
 //
-// Draw copies of the histograms
+//  Draw copies of the histograms
 //
 TObject *MHFlux::DrawClone(Option_t *opt) const
@@ -622 +552 @@
 // -------------------------------------------------------------------------
 //
-// Draw the histograms
+//  Draw the histograms
 //
 void MHFlux::Draw(Option_t *opt)
@@ -642 +572 @@
     gPad->Modified();
     gPad->Update();
+}
+
+Double_t MHFlux::Parab(Double_t x1, Double_t x2, Double_t x3,
+                       Double_t y1, Double_t y2, Double_t y3,
+                       Double_t val)
+{
+    Double_t c0, c1, c2;
+    Parab(x1, x2, x3, y1, y2, y3, &c0, &c1, &c2);
+    return c0 + c1*val + c2*val*val;
 }
 
@@ -660 +599 @@
         - x2*x1*x1 - x3*x2*x2 - x1*x3*x3;
 
-    if (det == 0.0)
+    if (det==0)
     {
         *a = 0;
@@ -686 +625 @@
     *c = (ai31*y1 + ai32*y2 + ai33*y3) * det1;
 
-    //yt1 = *a + *b * x1 + *c * x1*x1;
-    //yt2 = *a + *b * x2 + *c * x2*x2;
-    //yt3 = *a + *b * x3 + *c * x3*x3;
-
-    //*fLog << "x1 = " << x1 << ",  x2 = " << x2 << ",  x3 = " << x3 << endl;
-    //*fLog << "y1 = " << y1 << ",  y2 = " << y2 << ",  y3 = " << y3 << endl;
-    //*fLog << "yt1 = " << yt1 << ",  yt2 = " << yt2
-    //    << ",  yt3 = " << yt3 << endl;
-    //*fLog << "*a = " << *a << ",  *b = " << *b << ",  *c= " << *c
-    //      << ",  *errflag = " << *errflag << endl;
-
     return kTRUE;
 }

trunk/MagicSoft/Mars/mhist/MHFlux.h

@@ -20 +20 @@
 class TH2D;
 
+class MHGamma;
+class MHEffOnTime;
+class MHThetabarTheta;
+
 // base class MH is used because it defines "MakeDefCanvas"
 // if base class MH is used one has to define the member function Fill
@@ -35 +39 @@
     // all these plots for different bins of the variable (Theta or time)
 
+    void CalcEffCol(const TAxis &axex, TH2D &aeff, Int_t n, Double_t aeffbar[], Double_t daeffbar[]);
+    Double_t ParabInterpolCos(const TAxis &axe, const TH2D *aeff, Int_t j, Int_t k3, Double_t val) const;
+    void FindBins(const TAxis &axe, const Double_t bar, Int_t &k3, Int_t &k0) const;
+    void CalcAbsGammaFlux(const TH1D &teff, const TH2D &fHAeff);
+    Double_t ParabInterpolLog(const TAxis &axe, Int_t j,
+                              Double_t y[], Double_t Ebar) const;
+
 public:
-    MHFlux(const TH2D &h2d, const Bool_t Draw,
-           const TString varname, const TString unit);
+    MHFlux(const TH2D &h2d, const TString varname, const TString unit);
+    MHFlux(const MHGamma &h2d, const TString varname, const TString unit);
 
     Bool_t Fill(const MParContainer *par);
 
-    void Unfold(const Bool_t Draw);
+    void Unfold();
     void CalcFlux(const TH1D *teff, const TProfile *thetabar,
-                  const TH2D *aeff, const Bool_t Draw);
+                  const TH2D *aeff);
+
+    void CalcFlux(const MHEffOnTime &teff,
+                  const MHThetabarTheta &thetabar,
+                  const TH2D *aeff);
 
     void Draw(Option_t *option="");
     TObject *DrawClone(Option_t *option="") const;
+    void DrawFluxProjectionX(Option_t *opt="") const;
+    void DrawOrigProjectionX(Option_t *opt="") const;
 
     const TH2D *GetHOrig()       { return &fHOrig; }
@@ -56 +73 @@
                         double *a, double *b, double *c);
 
-    ClassDef(MHFlux, 1) //2D-plots (original, unfolded, flux)
+    static Double_t Parab(double x1, double x2, double x3,
+                          double y1, double y2, double y3,
+                          double val);
+
+    ClassDef(MHFlux, 0) //2D-plots (original, unfolded, flux)
 };
 

trunk/MagicSoft/Mars/mhist/MHGamma.cc

@@ -38 +38 @@
 #include <math.h>
 
-#include <TH1.h>
 #include <TH2.h>
 #include <TH3.h>
 
+#include "MHAlphaEnergyTheta.h"
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -52 +53 @@
 // Default Constructor. 
 //
-MHGamma::MHGamma()
-{
+MHGamma::MHGamma(const TString &name, const TString &title)
+    : fHist(NULL), fProject(NULL)
+{
+    fName  = name.IsNull()  ? (TString)"MHGamma" : name;
+    fTitle = title.IsNull() ? (TString)"3D-histogram of Alpha, E_est, Theta (Gamma sample)" : title;
 }
 
@@ -64 +68 @@
   return kTRUE;
 }
+
+TH3D *MHGamma::Subtract(const MHAlphaEnergyTheta &h1, const MHAlphaEnergyTheta &h2)
+{
+    return Subtract(h1.GetHist(), h2.GetHist());
+}
+
+TObject *MHGamma::DrawClone(Option_t *opt="") const
+{
+    DrawClone1();
+    DrawClone2();
+
+    return NULL;
+}
+
+void MHGamma::DrawClone1() const
+{
+    if (!fHist)
+        return;
+
+    //
+    // ------------- Part 1 ---------------------
+    //
+    TString titlex = fHist->GetXaxis()->GetTitle();
+    TString titley = fHist->GetYaxis()->GetTitle();
+    TString titlez = fHist->GetYaxis()->GetTitle();
+
+    TString canvtitle = "3D-plot "; //of ";
+    /*
+     canvtitle += titlex;
+     canvtitle += ", ";
+     canvtitle += titley;
+     canvtitle += ", ";
+     canvtitle += titlez+" ";
+     */
+    canvtitle += "for \'gamma\' sample";
+
+    TCanvas &c = *MakeDefCanvas("Alpha", canvtitle);
+
+    c.Divide(2, 2);
+
+    gROOT->SetSelectedPad(NULL);
+
+    TH1 *h;
+
+    c.cd(1);
+    h = ((TH3D*)(fHist))->Project3D(fName+"_ex");
+
+    TString title= "Source-Antisource: ";
+    h->SetTitle(title + titlex);
+    h->SetXTitle(titlex);
+    h->SetYTitle("Counts");
+
+    h->Draw();
+    h->SetBit(kCanDelete);
+
+    c.cd(2);
+    h = ((TH3D*)(fHist))->Project3D(fName+"_ey");
+
+    h->SetTitle(title + titley);
+    h->SetXTitle(titley);
+    h->SetYTitle("Counts");
+
+    h->Draw();
+    h->SetBit(kCanDelete);
+    gPad->SetLogx();
+
+    c.cd(3);
+    h = ((TH3D*)(fHist))->Project3D(fName+"_ez");
+
+    h->SetTitle(title + titlez);
+    h->SetXTitle(titlez);
+    h->SetYTitle("Counts");
+
+    h->Draw();
+    h->SetBit(kCanDelete);
+
+    c.cd(4);
+    ((TH3D*)fHist)->DrawCopy();
+
+    c.Modified();
+    c.Update();
+}
+
 
 // --------------------------------------------------------------------------
@@ -70 +157 @@
 //          fHist(gamma) = h1(source) - h2(antisource)
 // 
-TH3D *MHGamma::Subtract(const TH3D *h1, const TH3D *h2,
-                        const char *name, const char *title,
-                        Bool_t Draw)
-{
-    TH3D *fHist;
-    fHist = new TH3D();
-    fHist->SetName(name);
-    fHist->SetTitle(title);
-
+TH3D *MHGamma::Subtract(const TH3D *h1, const TH3D *h2)
+{
+    if (fHist)
+        delete fHist;
+
+    fHist = new TH3D;
+    fHist->SetName(fName);
+    fHist->SetTitle(fTitle);
     fHist->SetDirectory(NULL);
 
-    // SetBinning((TH3D*)fHist, (TH3D*)h1);
     SetBinning((TH1*)fHist, (TH1*)h1);
 
-    TString strg1 =   (((TH1*)h1)->GetXaxis())->GetTitle();
-    TString strg2 =   (((TH1*)h1)->GetYaxis())->GetTitle();
-    TString strg3 =   (((TH1*)h1)->GetZaxis())->GetTitle();
-    fHist->SetXTitle(strg1);
-    fHist->SetYTitle(strg2);
-    fHist->SetZTitle(strg3);
+    fHist->SetXTitle((((TH1*)h1)->GetXaxis())->GetTitle());
+    fHist->SetYTitle((((TH1*)h1)->GetYaxis())->GetTitle());
+    fHist->SetZTitle((((TH1*)h1)->GetZaxis())->GetTitle());
 
     fHist->Add((TH1*)h1, (TH1*)h2, 1, -1); // ROOT: FIXME!
-
-    //...........................................................
-    // draw histogram
-    if (Draw == kTRUE)
-    {
-      TString strg7 = "3D-plot of ";
-      strg7 += strg1;
-      strg7 += ",";
-      strg7 += strg2;
-      strg7 += ",";
-      strg7 += strg3;
-      strg7 += "  for \'gamma\' sample";
-
-      TCanvas &c = *MakeDefCanvas("Alpha", strg7); 
-
-      c.Divide(2, 2);
-
-      gROOT->SetSelectedPad(NULL);
-
-      TH1 *h;
-
-      c.cd(1);
-      h = ((TH3D*)(fHist))->Project3D("ex");
-
-      TString strg0 = "SRC-ASRC :    ";
-      TString strg4 = strg0 + strg1;
-      h->SetTitle(strg4);
-      h->SetXTitle(strg1);
-      h->SetYTitle("Counts");
-
-      h->Draw();
-      h->SetBit(kCanDelete);
-
-      c.cd(2);
-      h = ((TH3D*)(fHist))->Project3D("ey");
-
-      TString strg5 = strg0 + strg2;
-      h->SetTitle(strg5);
-      h->SetXTitle(strg2);
-      h->SetYTitle("Counts");
-
-      h->Draw();
-      h->SetBit(kCanDelete);
-      gPad->SetLogx();
-
-      c.cd(3);
-      h = ((TH3D*)(fHist))->Project3D("ez");
-
-      TString strg6 = strg0 + strg3;
-      h->SetTitle(strg6);
-      h->SetXTitle(strg3);
-      h->SetYTitle("Counts");
-
-      h->Draw();
-      h->SetBit(kCanDelete);
-
-      c.cd(4);
-      ((TH3D*)fHist)->DrawCopy();
-
-      c.Modified();
-      c.Update();
-    }
 
     return fHist;
@@ -162 +182 @@
 // Integrate fHist(gamma) in the alpha range (lo, up)
 // 
-TH2D *MHGamma::GetAlphaProjection(TH3D *fHist, Axis_t lo, Axis_t up,
-                                  Bool_t Drawp)
+TH2D *MHGamma::GetAlphaProjection(Axis_t lo, Axis_t up)
 {
     if (up < lo)
@@ -204 +223 @@
     axe.SetRange(ilo, iup);
 
-    TH2D &h2D = *(TH2D *)fHist->Project3D("ezy");
-
-    TString strg0 = "2D-plot of ";
-    TString strg1 = (fHist->GetYaxis())->GetTitle();
-    TString strg2 = (fHist->GetZaxis())->GetTitle();
-    strg0 += strg2;
-    strg0 += " vs. ";
-    strg0 += strg1;
-    h2D.SetTitle(strg0);
-    h2D.SetXTitle(strg1);
-    h2D.SetYTitle(strg2);
-
-
-    //...........................................................
-    // draw histogram
-    if (Drawp == kTRUE)
-    {
-      char txt[100];
-      TString strg3 = "No.of Gammas vs. ";
-      strg3 += strg1;
-      strg3 += " and ";
-      strg3 += strg2;
-      sprintf(txt, "   (%.1f < alpha < %.1f deg)", lo, up);
-      strg3 += txt;
-
-      TCanvas &c = *MakeDefCanvas("Gamma", strg3);
-
-      c.Divide(2, 2);
-
-      gROOT->SetSelectedPad(NULL);
-
-      TH1 *h;
-
-      c.cd(1);
-      h = h2D.ProjectionX("xpro", -1, 9999, "E");
-      TString strg0 = "No.of gammas : ";
-      TString strg7 = strg0 + strg1;
-      h->SetTitle(strg7);
-      h->SetXTitle(strg1);
-      h->SetYTitle("Counts");
-
-      h->Draw();
-      h->SetBit(kCanDelete);
-      gPad->SetLogx();
-
-      c.cd(2);
-      h = h2D.ProjectionY("ypro", -1, 9999, "E");
-      TString strg8 = strg0 + strg2;
-      h->SetTitle(strg8);
-      h->SetXTitle(strg2);
-      h->SetYTitle("Counts");
-
-      h->Draw();
-      h->SetBit(kCanDelete);
-
-      c.cd(3);
-
-      h2D.DrawCopy();
-      gPad->SetLogx();
-
-      c.Modified();
-      c.Update();
-    }
-    //...........................................................
-
-    return &h2D;
-}
+    fLo = lo;
+    fHi = up;
+
+    if (fProject)
+        delete fProject;
+    fProject = (TH2D*)fHist->Project3D(fName+"_ezy");
+
+    const TString title = "2D-plot of ";
+    const TString titley = fHist->GetYaxis()->GetTitle();
+    const TString titlez = fHist->GetZaxis()->GetTitle();
+
+    fProject->SetTitle(title + titley + " vs. " + titlez);
+    fProject->SetXTitle(titley);
+    fProject->SetYTitle(titlez);
+
+    return fProject;
+}
+
+void MHGamma::DrawClone2() const
+{
+    if (!fProject)
+        return;
+
+    const TString titley = fHist->GetYaxis()->GetTitle();
+    const TString titlez = fHist->GetZaxis()->GetTitle();
+
+    TString canvtitle = "No.of Gammas ";//vs. ";
+    /*
+     canvtitle += titley;
+     canvtitle += " and ";
+     canvtitle += titlez;
+     */
+    canvtitle += Form("(%.1f < alpha < %.1f deg)", fLo, fHi);
+
+    TCanvas &c = *MakeDefCanvas("Gamma", canvtitle);
+
+    c.Divide(2, 2);
+
+    gROOT->SetSelectedPad(NULL);
+
+    TH1 *h;
+
+    c.cd(1);
+    h = fProject->ProjectionX(fName+"_xpro", -1, 9999, "E");
+    TString title = "No.of gammas: ";
+    h->SetTitle(title+titley);
+    h->SetXTitle(titley);
+    h->SetYTitle("Counts");
+
+    h->Draw();
+    h->SetBit(kCanDelete);
+    gPad->SetLogx();
+
+    c.cd(2);
+    h = fProject->ProjectionY(fName+"_ypro", -1, 9999, "E");
+    h->SetTitle(title+titlez);
+    h->SetXTitle(titlez);
+    h->SetYTitle("Counts");
+
+    h->Draw();
+    h->SetBit(kCanDelete);
+
+    c.cd(3);
+
+    fProject->DrawCopy();
+    gPad->SetLogx();
+
+    c.Modified();
+    c.Update();
+}

trunk/MagicSoft/Mars/mhist/MHGamma.h

@@ -6 +6 @@
 #endif
 
-#ifndef ROOT_TH3
-#include "TH3.h"
-#endif
-
-#ifndef ROOT_TH2
-#include "TH2.h"
-#endif
-
 class TH2D;
 class TH3D;
 
+class MHAlphaEnergyTheta;
+
 class MHGamma : public MH 
 {
+private:
+    TH3D *fHist;    //!
+    TH2D *fProject; //!
+
+    Axis_t fLo; //!
+    Axis_t fHi; //!
+
 public:
-    MHGamma();
+    MHGamma(const TString &name="", const TString &title="");
 
     Bool_t Fill(const MParContainer *par);
 
-    TH3D *Subtract(const TH3D *h1, const TH3D *h2,
-                   const char *name, const char *title, Bool_t Draw=kFALSE);
+    TH3D *Subtract(const TH3D *h1, const TH3D *h2);
 
-    TH2D *GetAlphaProjection(TH3D *fHist, Axis_t lo, Axis_t up, 
-                             Bool_t Drawp=kFALSE);
+    TH3D *Subtract(const MHAlphaEnergyTheta &h1, const MHAlphaEnergyTheta &h2);
 
+    TH2D *GetAlphaProjection(Axis_t lo, Axis_t up);
 
-    ClassDef(MHGamma, 1) // manipulation of alpha distributions
+    TObject *DrawClone(Option_t *opt="") const;
+    void DrawClone1() const;
+    void DrawClone2() const;
+
+    const TH2D *GetProject() const { return fProject; }
+
+    ClassDef(MHGamma, 0) // manipulation of alpha distributions
 };
 

trunk/MagicSoft/Mars/mhist/MHHadronness.cc

@@ -45 +45 @@
 //    - black: acceptance of gammas (Ag) vs. the hadroness
 //    - red:   acceptance of non gammas (Ah) vs. the hadroness
-//    - blue:  2D distance of (acceptance_hadrons, acceptances_gammas)
-//             to optimum (0, 1)
-//             1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
 //  * Bottom Left Corner:
 //    Naive quality factor: Ag/sqrt(Ah)
@@ -56 +53 @@
 ////////////////////////////////////////////////////////////////////////////
 #include "MHHadronness.h"
+
+/*
+ // - blue:  2D distance of (acceptance_hadrons, acceptances_gammas)
+ //          to optimum (0, 1):  1-sqrt(Ag*Ag + (1-Ah)*(1-Ah))
+ */
 
 #include <TPad.h>
@@ -113 +115 @@
     fQfac->SetTitle(" Naive Quality factor ");
 
-    fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
-    fMinDist->SetXTitle("Hadronness");
-    fMinDist->SetYTitle("Distance");
+    /*
+     fMinDist = new TH1D("MinDist", "Minimum Dist to (0, 1)", nbins, 0, 1);
+     fMinDist->SetXTitle("Hadronness");
+     fMinDist->SetYTitle("Distance");
+     */
 
     //fQfac->SetDirectory(NULL);
     fGhness->SetDirectory(NULL);
     fPhness->SetDirectory(NULL);
-    fMinDist->SetDirectory(NULL);
+    //fMinDist->SetDirectory(NULL);
     fIntGhness->SetDirectory(NULL);
     fIntPhness->SetDirectory(NULL);
@@ -136 +140 @@
     delete fIntPhness;
     delete fQfac;
-    delete fMinDist;
+    //    delete fMinDist;
     delete fGraph;
 }
@@ -266 +270 @@
             continue;
 
-        fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
-
-        Double_t val = ig/sqrt(ip);
+        //fMinDist->SetBinContent(i, 1-sqrt(ip*ip + (ig-1)*(ig-1)));
+
+        const Double_t val = ig/sqrt(ip);
         fQfac->SetPoint(i, ig, val);
 
@@ -360 +364 @@
 
     *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl;
-    *fLog << " acc(hadron)  acc(gamma)" << endl <<endl;
-
-    *fLog << "    0.005   "  << Form("%6.3f", GetGammaAcceptance(0.005)) << endl;
-    *fLog << "    0.02    "  << Form("%6.3f", GetGammaAcceptance(0.02))  << endl;
-    *fLog << "    0.05    "  << Form("%6.3f", GetGammaAcceptance(0.05))  << endl;
-    *fLog << "    0.1     "  << Form("%6.3f", GetGammaAcceptance(0.1 ))  << endl;
-    *fLog << "    0.2     "  << Form("%6.3f", GetGammaAcceptance(0.2 ))  << endl;
-    *fLog << "    0.3     "  << Form("%6.3f", GetGammaAcceptance(0.3 ))  << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << "       0.1  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << "       0.2  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << "       0.3  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << "       0.4  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << "       0.5  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << "       0.6  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << "       0.7  " << endl;
-    *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << "       0.8  " << endl;
+    *fLog << "acc(hadron)  acc(gamma)  acc(g)/acc(h)" << endl <<endl;
+
+    *fLog << "    0.005    " << Form("%6.3f", GetGammaAcceptance(0.005)) << "    " << Form("%6.3f", GetGammaAcceptance(0.005)/0.005) << endl;
+    *fLog << "    0.02     " << Form("%6.3f", GetGammaAcceptance(0.02))  << "    " << Form("%6.3f", GetGammaAcceptance(0.02)/0.02) << endl;
+    *fLog << "    0.05     " << Form("%6.3f", GetGammaAcceptance(0.05))  << "    " << Form("%6.3f", GetGammaAcceptance(0.05)/0.05) << endl;
+    *fLog << "    0.1      " << Form("%6.3f", GetGammaAcceptance(0.1 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.1)/0.1) << endl;
+    *fLog << "    0.2      " << Form("%6.3f", GetGammaAcceptance(0.2 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.2)/0.2) << endl;
+    *fLog << "    0.3      " << Form("%6.3f", GetGammaAcceptance(0.3 ))  << "    " << Form("%6.3f", GetGammaAcceptance(0.3)/0.3) << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.1)) << "        0.1  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.2)) << "        0.2  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.3)) << "        0.3  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.4)) << "        0.4  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.5)) << "        0.5  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.6)) << "        0.6  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.7)) << "        0.7  " << endl;
+    *fLog << Form("%6.3f", GetHadronAcceptance(0.8)) << "        0.8  " << endl;
     *fLog << endl;
 
-    const Int_t h = fMinDist->GetMaximumBin();
-    *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
-    *fLog << "  Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
-    *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
+    /*
+     const Int_t h = fMinDist->GetMaximumBin();
+     *fLog << "Minimum Distance to (0, 1): " << Form("%.2f", 1-fMinDist->GetMaximum()) << " @ H=" << fMinDist->GetBinCenter(h) << endl;
+     *fLog << "  Acc Gammas = " << Form("%3.0f", fIntGhness->GetBinContent(h)*100) << "%, ";
+     *fLog << "Acc Hadrons = " << Form("%3.0f", fIntPhness->GetBinContent(h)*100) << "%" << endl;
+     */
 
     *fLog << "Q-Factor @ Acc Gammas=0.5: Q(0.5)=" << Form("%.1f", GetQ05()) << endl;
@@ -412 +418 @@
     Getiphness()->SetLineColor(kRed);
     Getiphness()->DrawCopy("same");
-    fMinDist->SetLineColor(kBlue);
-    fMinDist->DrawCopy("same");
+    /*
+     fMinDist->SetLineColor(kBlue);
+     fMinDist->DrawCopy("same");
+     */
     c.cd(3);
     //GetQfac()->DrawCopy();
@@ -449 +457 @@
         gPad->Update();
     }
+    /*
     const Int_t h = fMinDist->GetMaximumBin();
     TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
@@ -455 +464 @@
     m->SetBit(kCanDelete);
     m->Draw();
+    */
     return &c;
 }
@@ -482 +492 @@
     Getiphness()->SetLineColor(kRed);
     Getiphness()->Draw("same");
-    fMinDist->SetLineColor(kBlue);
-    fMinDist->Draw("same");
+    /*
+     fMinDist->SetLineColor(kBlue);
+     fMinDist->Draw("same");
+     */
     gPad->cd(3);
     //GetQfac()->Draw();
@@ -516 +528 @@
         gPad->Update();
     }
+    /*
     const Int_t h = fMinDist->GetMaximumBin();
     TMarker *m = new TMarker(fIntPhness->GetBinContent(h),
@@ -522 +535 @@
     m->SetBit(kCanDelete);
     m->Draw();
-}
+    */
+}

trunk/MagicSoft/Mars/mhist/MHHadronness.h

@@ -15 +15 @@
 {
 private:
-    const MMcEvt *fMcEvt;            //!
-    const MHadronness *fHadronness;    //!
+    const MMcEvt *fMcEvt;           //!
+    const MHadronness *fHadronness; //!
 
-    TH1D* fPhness;        //-> Hadrons Hadronness
-    TH1D* fGhness;        //-> Gammas Hadronness
-    TH1D* fIntPhness;     //-> Hadrons Acceptance
-    TH1D* fIntGhness;     //-> Gammas Acceptance
-    TH1D* fMinDist;       //-> Minimum Distance to optimum acceptance
+    TH1D* fPhness;    //-> Hadrons Hadronness
+    TH1D* fGhness;    //-> Gammas  Hadronness
+    TH1D* fIntPhness; //-> Hadrons Acceptance
+    TH1D* fIntGhness; //-> Gammas  Acceptance
+    //TH1D* fMinDist; //-> Minimum Distance to optimum acceptance
 
-    TGraph *fQfac;        //-> Quality factor
-    TGraph *fGraph;       //-> gamma acceptance vs. hadron acceptance
+    TGraph *fQfac;    //-> Quality factor
+    TGraph *fGraph;   //-> gamma acceptance vs. hadron acceptance
 
 public:

trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc

@@ -112 +112 @@
 // Fill data into the histogram which contains all showers
 //
-void MHMcCollectionArea::FillAll(Float_t energy, Float_t radius)
+void MHMcCollectionArea::FillAll(Double_t energy, Double_t radius)
 {
     fHistAll->Fill(energy, radius);
@@ -121 +121 @@
 // Fill data into the histogram which contains the selected showers
 //
-void MHMcCollectionArea::FillSel(Float_t energy, Float_t radius)
+void MHMcCollectionArea::FillSel(Double_t energy, Double_t radius)
 {
     fHistSel->Fill(energy, radius);

trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h

@@ -26 +26 @@
     ~MHMcCollectionArea();
 
-    void FillAll(Float_t energy, Float_t radius);
-    void FillSel(Float_t energy, Float_t radius);
+    void FillAll(Double_t energy, Double_t radius);
+    void FillSel(Double_t energy, Double_t radius);
 
     void DrawAll(Option_t *option="");

trunk/MagicSoft/Mars/mhist/MHMcEnergyMigration.h

@@ -40 +40 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHMcEnergyMigration, 1) //3D-histogram   E-true E-est Theta
+    ClassDef(MHMcEnergyMigration, 0) //3D-histogram   E-true E-est Theta
 
 };

trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc

@@ -99 +99 @@
     gROOT->SetSelectedPad(NULL);
 
-    ((TProfile*)&fHist)->DrawCopy(opt);
+    ((TProfile)fHist).DrawCopy(opt);
 
     c.Modified();
@@ -128 +128 @@
 Bool_t MHThetabarTheta::Fill(const MParContainer *par)
 {
-    fHist.Fill(fMcEvt->GetTheta()*kRad2Deg, fMcEvt->GetTheta()*kRad2Deg);
+    const Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
+
+    fHist.Fill(theta, theta);
 
     return kTRUE;
 }
-
-
-
-
-
-
-
-
-
-
-
-
-

trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h

@@ -17 +17 @@
 {
 private:
-    MTime  *fTime;   //!
-    MMcEvt *fMcEvt;  //!
+    MTime   *fTime;   //!
+    MMcEvt  *fMcEvt;  //!
 
-    TProfile    fHist;
+    TProfile fHist;
 
 public:
@@ -36 +36 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHThetabarTheta, 1) //Profile histogram Thetabar vs. time
+    ClassDef(MHThetabarTheta, 0) //Profile histogram Thetabar vs. time
 
 };

trunk/MagicSoft/Mars/mhist/MHThetabarTime.h

@@ -39 +39 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHThetabarTime, 1) //Profile histogram Thetabar vs. time
+    ClassDef(MHThetabarTime, 0) //Profile histogram Thetabar vs. time
 
 };

trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.cc

@@ -187 +187 @@
     gPad->Modified();
     gPad->Update();
-
 }
 
@@ -196 +195 @@
 Bool_t MHTimeDiffTheta::Fill(const MParContainer *par)
 {
-    const Int_t time = fTime->GetTimeLo();
-
-    fHist.Fill(0.0001*(time-fLastTime), fMcEvt->GetTheta()*kRad2Deg);
+    const Double_t time = 200e-9*fTime->GetTimeLo() + fTime->GetTimeHi();
+
+    fHist.Fill(time-fLastTime, fMcEvt->GetTelescopeTheta()*kRad2Deg);
     fLastTime = time;
 

trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.h

@@ -19 +19 @@
     MTime  *fTime;   //!
     MMcEvt *fMcEvt;  //!
-    Int_t   fLastTime;
+    Double_t fLastTime;
 
     TH2D    fHist;
@@ -37 +37 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHTimeDiffTheta, 1) //2D-histogram  time-diff vs. Theta
+    ClassDef(MHTimeDiffTheta, 0) //2D-histogram  time-diff vs. Theta
 };
 

trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.cc

@@ -100 +100 @@
 {
 
-    TCanvas &c = *MakeDefCanvas("DiffTimeTime", "Distrib of \\Delta t, time");
+    TCanvas &c = *MakeDefCanvas("DiffTimeTime", "Distrib of dt and t");
 
     c.Divide(2, 2);
@@ -186 +186 @@
 Bool_t MHTimeDiffTime::Fill(const MParContainer *par)
 {
-    const Int_t time = fTime->GetTimeLo();
+    const Double_t time = 200e-9*fTime->GetTimeLo() + fTime->GetTimeHi();
 
-    fHist.Fill(0.0001*(time-fLastTime), 0.0001*time);
-
+    fHist.Fill(time-fLastTime, time);
     fLastTime = time;
 

trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.h

@@ -17 +17 @@
 private:
     MTime *fTime;   //!
-    Int_t  fLastTime;
+    Double_t  fLastTime;
 
     TH2D   fHist;
@@ -35 +35 @@
     TObject *DrawClone(Option_t *option="") const;
 
-    ClassDef(MHTimeDiffTime, 1) //2D-histogram  time-diff vs. time
+    ClassDef(MHTimeDiffTime, 0) //2D-histogram  time-diff vs. time
 };
 

trunk/MagicSoft/Mars/mmain/MAnalysis.cc

@@ -100 +100 @@
 }
 
-MAnalysis::MAnalysis(const TGWindow *main, const TGWindow *p,
+MAnalysis::MAnalysis(/*const TGWindow *main,*/ const TGWindow *p,
                      const UInt_t w, const UInt_t h)
-: MBrowser(main, p, w, h)
+: MBrowser(/*main,*/ p, w, h)
 {
     AddButtons();

trunk/MagicSoft/Mars/mmain/MAnalysis.h

@@ -28 +28 @@
 
 public:
-    MAnalysis(const TGWindow *main=NULL, const TGWindow *p=NULL,
+    MAnalysis(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
               const UInt_t w=500, const UInt_t h=500) ;
 

trunk/MagicSoft/Mars/mmain/MBrowser.cc

@@ -275 +275 @@
 }
 
-MBrowser::MBrowser(const TGWindow *main, const TGWindow *p,
-                     const UInt_t w, const UInt_t h)
-: TGTransientFrame(p?p:gClient->GetRoot(),
-                   main?main:gClient->GetRoot(), w, h)
+MBrowser::MBrowser(/*const TGWindow *main,*/ const TGWindow *p,
+                   const UInt_t w, const UInt_t h)
+//: TGTransientFrame(p?p:gClient->GetRoot(),
+//                   main?main:gClient->GetRoot(), w, h)
+: TGMainFrame(p?p:gClient->GetRoot(), w, h)
 {
     fList = new MGList;
@@ -291 +292 @@
     fList->Add(linesep);
     AddFrame(linesep, laylinesep);
-
 
     //

trunk/MagicSoft/Mars/mmain/MBrowser.h

@@ -19 +19 @@
 class TGFileContainer;
 
-class MBrowser : public TGTransientFrame
+class MBrowser : public TGMainFrame/*TGTransientFrame*/
 {
 private:
@@ -62 +62 @@
 
  public: 
-     MBrowser(const TGWindow *main=NULL, const TGWindow *p=NULL,
-               const UInt_t w=500, const UInt_t h=500) ;
+     MBrowser(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
+              const UInt_t w=500, const UInt_t h=500) ;
 
      virtual ~MBrowser();

trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc

@@ -43 +43 @@
 //  'new MCameraDisplay()'
 //
-MCameraDisplay::MCameraDisplay(const TGWindow *main, const TGWindow *p,
-                       const UInt_t w, const UInt_t h)
-: MBrowser(main, p, w, h)
+MCameraDisplay::MCameraDisplay(/*const TGWindow *main,*/ const TGWindow *p,
+                               const UInt_t w, const UInt_t h)
+: MBrowser(/*main,*/ p, w, h)
 {
     TGTextButton *disp = new TGTextButton(fTop2, "Display Events", kButDisplay);
@@ -91 +91 @@
         case kButDisplay:
             new MGCamDisplay(fInputFile,
-                             fClient->GetRoot(), this, 600, 500);
+                             fClient->GetRoot(), /*this,*/ 600, 500);
             return kTRUE;
         }

trunk/MagicSoft/Mars/mmain/MCameraDisplay.h

@@ -9 +9 @@
 {
 public:
-    MCameraDisplay(const TGWindow *main=NULL, const TGWindow *p=NULL,
+    MCameraDisplay(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
                    const UInt_t w=500, const UInt_t h=500) ;
 

trunk/MagicSoft/Mars/mmain/MDataCheck.cc

@@ -54 +54 @@
 //  'new MDataCheck()'
 //
-MDataCheck::MDataCheck(const TGWindow *main, const TGWindow *p,
+MDataCheck::MDataCheck(/*const TGWindow *main,*/ const TGWindow *p,
                        const UInt_t w, const UInt_t h)
-: MBrowser(main, p, w, h)
+: MBrowser(/*main,*/ p, w, h)
 {
     TGTextButton *pedadc = new TGTextButton(fTop1, "ADC Spectra of Pedestals", kButPedAdc);

trunk/MagicSoft/Mars/mmain/MDataCheck.h

@@ -16 +16 @@
 
 public:
-    MDataCheck(const TGWindow *main=NULL, const TGWindow *p=NULL,
+    MDataCheck(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
                const UInt_t w=500, const UInt_t h=500) ;
 

trunk/MagicSoft/Mars/mmain/MEvtDisp.cc

@@ -39 +39 @@
 };
 
-MEvtDisp::MEvtDisp(const TGWindow *main, const TGWindow *p,
+MEvtDisp::MEvtDisp(/*const TGWindow *main,*/ const TGWindow *p,
                    const UInt_t w, const UInt_t h)
-: MBrowser(main, p, w, h)
+: MBrowser(/*main,*/ p, w, h)
 {
     TGTextButton *fadcevt = new TGTextButton(fTop1, "FADC Disp for Events",  kButDispEvt);
@@ -96 +96 @@
         case kButDispEvt:
             new MGFadcDisp(fInputFile, "Events",
-                           fClient->GetRoot(), this, 600, 500);
+                           fClient->GetRoot(), /*this,*/ 600, 500);
             return kTRUE;
 
         case kButDispPedestal:
             new MGFadcDisp(fInputFile , "PedEvts",
-                           fClient->GetRoot(), this, 600, 500);
+                           fClient->GetRoot(), /*this,*/ 600, 500);
             return kTRUE;
 
         case kButDispCalibration:
             new MGFadcDisp(fInputFile , "CalEvts",
-                           fClient->GetRoot(), this, 600, 500);
+                           fClient->GetRoot(), /*this,*/ 600, 500);
             return kTRUE;
         }

trunk/MagicSoft/Mars/mmain/MEvtDisp.h

@@ -9 +9 @@
 { 
 public:
-    MEvtDisp(const TGWindow *main=NULL, const TGWindow *p=NULL,
+    MEvtDisp(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
              const UInt_t w=500, const UInt_t h=500) ;
 

trunk/MagicSoft/Mars/mmain/MMars.cc

@@ -277 +277 @@
 
             case kButEvtDisplay:
-                new MEvtDisp(this);
+                new MEvtDisp(/*this*/);
                 return kTRUE;
 
             case kButDataCheck:
-                new MDataCheck(this);
+                new MDataCheck(/*this*/);
                 return kTRUE;
 
             case kButAnalysis:
-                new MAnalysis(this);
+                new MAnalysis(/*this*/);
                 return kTRUE;
 
             case kButMonteCarlo:
-                new MMonteCarlo(this);
+                new MMonteCarlo(/*this*/);
                 return kTRUE;
 
             case kButCameraDisplay:
-                new MCameraDisplay(this);
+                new MCameraDisplay(/*this*/);
                 return kTRUE;
 

trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc

@@ -162 +162 @@
 }
 
-MMonteCarlo::MMonteCarlo(const TGWindow *main, const TGWindow *p,
+MMonteCarlo::MMonteCarlo(/*const TGWindow *main,*/ const TGWindow *p,
                          const UInt_t w, const UInt_t h)
-: MBrowser(main, p, w, h)
+: MBrowser(/*main,*/ p, w, h)
 {
     AddButtons();

trunk/MagicSoft/Mars/mmain/MMonteCarlo.h

@@ -28 +28 @@
 
 public:
-    MMonteCarlo(const TGWindow *main=NULL, const TGWindow *p=NULL,
+    MMonteCarlo(/*const TGWindow *main=NULL,*/ const TGWindow *p=NULL,
                 const UInt_t w=500, const UInt_t h=500);