Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 1667) +++ /trunk/MagicSoft/Mars/Changelog (revision 1668) @@ -1,3 +1,90 @@ -*-*- 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 Index: /trunk/MagicSoft/Mars/mbase/MEvtLoop.cc =================================================================== --- /trunk/MagicSoft/Mars/mbase/MEvtLoop.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mbase/MEvtLoop.cc (revision 1668) @@ -208,14 +208,20 @@ 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); } Index: /trunk/MagicSoft/Mars/mbase/MFilter.h =================================================================== --- /trunk/MagicSoft/Mars/mbase/MFilter.h (revision 1667) +++ /trunk/MagicSoft/Mars/mbase/MFilter.h (revision 1668) @@ -10,13 +10,13 @@ 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; Index: /trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.cc =================================================================== --- /trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.cc (revision 1667) +++ /trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.cc (revision 1668) @@ -148,7 +148,7 @@ // 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) { // Index: /trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.h =================================================================== --- /trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.h (revision 1667) +++ /trunk/MagicSoft/Mars/meventdisp/MGCamDisplay.h (revision 1668) @@ -32,5 +32,5 @@ public: MGCamDisplay(const char *filename, - const TGWindow *p, const TGWindow *main, + const TGWindow *p, /*const TGWindow *main,*/ UInt_t w, UInt_t h); Index: /trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.cc =================================================================== --- /trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.cc (revision 1667) +++ /trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.cc (revision 1668) @@ -141,5 +141,4 @@ fEvtInfo = new TGLabel(top2, new TGString("")); fList->Add(fEvtInfo); - top2->AddFrame(fEvtInfo, laystd); @@ -393,7 +392,8 @@ 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,5 +472,6 @@ break; default: - txt += "Unknown Particle Id"; + txt += "Unknown Particle Id#"; + txt += evt->GetPartId(); } @@ -493,4 +494,12 @@ 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(); } Index: /trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.h =================================================================== --- /trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.h (revision 1667) +++ /trunk/MagicSoft/Mars/meventdisp/MGEvtDisplay.h (revision 1668) @@ -21,5 +21,5 @@ class MReadTree; -class MGEvtDisplay : public TGTransientFrame +class MGEvtDisplay : public TGMainFrame/*TGTransientFrame*/ { private: @@ -74,5 +74,5 @@ 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); Index: /trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.cc =================================================================== --- /trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.cc (revision 1667) +++ /trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.cc (revision 1668) @@ -113,7 +113,7 @@ // 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) { // Index: /trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.h =================================================================== --- /trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.h (revision 1667) +++ /trunk/MagicSoft/Mars/meventdisp/MGFadcDisp.h (revision 1668) @@ -23,5 +23,5 @@ 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); Index: /trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc (revision 1668) @@ -127,5 +127,7 @@ 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); Index: /trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 1668) @@ -16,8 +16,7 @@ ! ! -! Author(s): Thomas Bretz 12/2000 -! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 11/2002 ! -! Copyright: MAGIC Software Development, 2000-2001 +! Copyright: MAGIC Software Development, 2000-2002 ! ! @@ -29,4 +28,8 @@ // // 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,4 +61,7 @@ #include "MLogManip.h" +#include "MTime.h" +#include "MFilter.h" + #include "MParList.h" #include "MCerPhotEvt.h" @@ -112,5 +118,7 @@ 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,6 +510,9 @@ // 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,4 +659,11 @@ // + // 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,5 +703,5 @@ fNumRuns = 0; - return kTRUE; + return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE; } @@ -695,8 +713,11 @@ 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,4 +883,33 @@ 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,4 +940,4 @@ } - return kTRUE; -} + return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE; +} Index: /trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 1668) @@ -7,4 +7,5 @@ class TList; +class MTime; class MMcEvt; class MMcTrig; @@ -27,4 +28,5 @@ 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,4 +54,8 @@ void ProcessEvent(const struct eventrecord &event); + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + public: MCT1ReadPreProc(const char *filename=NULL, @@ -61,8 +67,4 @@ void AddFile(const char *fname); - Bool_t PreProcess(MParList *pList); - Bool_t Process(); - Bool_t PostProcess(); - UInt_t GetEntries() { return fEntries; } Index: /trunk/MagicSoft/Mars/mfileio/MRead.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MRead.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MRead.cc (revision 1668) @@ -29,4 +29,9 @@ // 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" Index: /trunk/MagicSoft/Mars/mfileio/MRead.h =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MRead.h (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MRead.h (revision 1668) @@ -6,9 +6,18 @@ #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 Index: /trunk/MagicSoft/Mars/mfileio/MReadMarsFile.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MReadMarsFile.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MReadMarsFile.cc (revision 1668) @@ -98,5 +98,5 @@ // 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,5 +107,5 @@ // Int_t n1 = fRun->AddFile(fname); - Int_t n2 = MReadTree::AddFile(fname); + Int_t n2 = MReadTree::AddFile(fname, entries); return n1 != n2 ? -1 : n1; Index: /trunk/MagicSoft/Mars/mfileio/MReadMarsFile.h =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MReadMarsFile.h (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MReadMarsFile.h (revision 1668) @@ -24,5 +24,5 @@ ~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 Index: /trunk/MagicSoft/Mars/mfileio/MReadTree.cc =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 1668) @@ -191,14 +191,44 @@ // 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,4 +238,14 @@ } +/* + // -------------------------------------------------------------------------- + // + // + Int_t MReadTree::AddFile(const TChainElement &obj) + { + return AddFile(obj.GetTitle(), obj.GetEntries()); + } +*/ + // -------------------------------------------------------------------------- // @@ -216,13 +256,17 @@ 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,5 +285,5 @@ 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,4 +299,11 @@ 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,4 +500,6 @@ // 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,5 +629,5 @@ fChain->SetNotify(this); - return kTRUE; + return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE; } @@ -626,4 +679,6 @@ // (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,5 +713,26 @@ #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,4 +740,13 @@ return rc; +} + +// -------------------------------------------------------------------------- +// +// If a selector is given the selector is post processed +// +Bool_t MReadTree::PostProcess() +{ + return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE; } Index: /trunk/MagicSoft/Mars/mfileio/MReadTree.h =================================================================== --- /trunk/MagicSoft/Mars/mfileio/MReadTree.h (revision 1667) +++ /trunk/MagicSoft/Mars/mfileio/MReadTree.h (revision 1668) @@ -64,9 +64,10 @@ 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(); Index: /trunk/MagicSoft/Mars/mfilter/MF.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MF.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfilter/MF.cc (revision 1668) @@ -244,6 +244,9 @@ 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,5 +410,5 @@ } - if (!fFilter->PreProcess(plist)) + if (!fFilter->CallPreProcess(plist)) { *fLog << err << dbginf << "PreProcessing filters in "; @@ -423,5 +426,5 @@ Bool_t MF::Process() { - return fFilter->Process(); + return fFilter->CallProcess(); } @@ -432,5 +435,5 @@ Bool_t MF::PostProcess() { - return fFilter->PostProcess(); + return fFilter->CallPostProcess(); } Index: /trunk/MagicSoft/Mars/mfilter/MFDataChain.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFDataChain.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfilter/MFDataChain.cc (revision 1668) @@ -58,5 +58,5 @@ // 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) { Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector.cc (revision 1668) @@ -147,5 +147,5 @@ Bool_t MFEventSelector::Process() { - Float_t evt = gRandom->Uniform(); + const Float_t evt = gRandom->Uniform(); if (fNumSelectEvts>0) Index: /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 1667) +++ /trunk/MagicSoft/Mars/mfilter/MFEventSelector.h (revision 1668) @@ -27,4 +27,8 @@ void StreamPrimitive(ofstream &out) const; + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + enum { kNumTotalFromFile = BIT(14) }; /* @@ -47,8 +51,4 @@ */ - Bool_t PreProcess(MParList *pList); - Bool_t Process(); - Bool_t PostProcess(); - ClassDef(MFEventSelector, 0) // A Filter to select events from files }; Index: /trunk/MagicSoft/Mars/mfilter/MFilterList.cc =================================================================== --- /trunk/MagicSoft/Mars/mfilter/MFilterList.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mfilter/MFilterList.cc (revision 1668) @@ -188,5 +188,5 @@ // while ((filter=(MFilter*)Next())) - if (!filter->PreProcess(pList)) + if (!filter->CallPreProcess(pList)) { *fLog << err << "Error - Preprocessing Filter "; @@ -212,5 +212,5 @@ // while ((filter=(MFilter*)Next())) - if (!filter->Process()) + if (!filter->CallProcess()) return kFALSE; @@ -232,5 +232,5 @@ // while ((filter=(MFilter*)Next())) - if (!filter->PostProcess()) + if (!filter->CallPostProcess()) return kFALSE; Index: /trunk/MagicSoft/Mars/mhist/MBinning.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MBinning.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MBinning.cc (revision 1668) @@ -64,6 +64,6 @@ // -------------------------------------------------------------------------- // -// Specify the number of bins (not the number of edges), the lowest and -// highest Edge. +// Specify the number of bins (not the number of edges), the +// lowest and highest Edge (of your histogram) // void MBinning::SetEdges(const Int_t nbins, const Axis_t lo, Axis_t up) @@ -79,6 +79,6 @@ // -------------------------------------------------------------------------- // -// Specify the number of bins (not the number of edges), the lowest and -// highest Edge. +// Specify the number of bins (not the number of edges), the +// lowest and highest Edge (of your histogram) // void MBinning::SetEdgesLog(const Int_t nbins, const Axis_t lo, Axis_t up) @@ -92,4 +92,23 @@ fType = kIsLogarithmic; +} + +// -------------------------------------------------------------------------- +// +// Specify the number of bins (not the number of edges), the +// lowest [deg] and highest [deg] 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,9 +151,11 @@ 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; Index: /trunk/MagicSoft/Mars/mhist/MBinning.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MBinning.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MBinning.h (revision 1668) @@ -25,4 +25,5 @@ kIsLinear, kIsLogarithmic, + kIsCosinic, kIsUserArray }; @@ -39,4 +40,5 @@ 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,4 +53,5 @@ 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; } Index: /trunk/MagicSoft/Mars/mhist/MFillH.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MFillH.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MFillH.cc (revision 1668) @@ -354,5 +354,5 @@ { 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); } Index: /trunk/MagicSoft/Mars/mhist/MH.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MH.cc (revision 1668) @@ -546,4 +546,26 @@ // -------------------------------------------------------------------------- // +// 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. Index: /trunk/MagicSoft/Mars/mhist/MH.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MH.h (revision 1668) @@ -54,4 +54,6 @@ 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); Index: /trunk/MagicSoft/Mars/mhist/MH3.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MH3.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MH3.h (revision 1668) @@ -56,5 +56,5 @@ const TH1 &GetHist() const { return *fHist; } - TH1 *GetHistByName(const TString name) { return fHist; } + TH1 *GetHistByName(const TString name="") { return fHist; } void SetColors() const; Index: /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.cc (revision 1668) @@ -64,7 +64,7 @@ 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,5 +79,5 @@ if (!fEnergy) { - *fLog << err << dbginf << "MHAlphaEnergyTheta : MEnergyEst not found... aborting." << endl; + *fLog << err << dbginf << "MEnergyEst not found... aborting." << endl; return kFALSE; } @@ -86,5 +86,5 @@ if (!fMcEvt) { - *fLog << err << dbginf << "MHAlphaEnergyTheta : MMcEvt not found... aborting." << endl; + *fLog << err << dbginf << "MMcEvt not found... aborting." << endl; return kFALSE; } @@ -95,5 +95,5 @@ 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,6 +103,4 @@ fHist.Sumw2(); - fHist.Sumw2(); - return kTRUE; } @@ -116,5 +114,6 @@ 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,5 +147,5 @@ h->SetTitle("Distribution of E-est [GeV]"); - h->SetXTitle("E-est [GeV] "); + h->SetXTitle("E_{est} [GeV]"); h->SetYTitle("Counts"); @@ -186,5 +185,5 @@ c.cd(1); - h = ((TH3*)(&fHist))->Project3D("expro"); + h = ((TH3*)(&fHist))->Project3D(fName+"_expro"); h->SetTitle("Distribution of \\alpha [\\circ]"); @@ -196,8 +195,8 @@ 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,5 +206,5 @@ c.cd(3); - h = ((TH3*)(&fHist))->Project3D("ezpro"); + h = ((TH3*)(&fHist))->Project3D(fName+"_ezpro"); h->SetTitle("Distribution of \\Theta [\\circ]"); Index: /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTheta.h (revision 1668) @@ -44,5 +44,5 @@ 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 }; Index: /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHAlphaEnergyTime.h (revision 1668) @@ -47,5 +47,5 @@ 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 }; Index: /trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHEffOnTime.cc (revision 1668) @@ -54,4 +54,6 @@ #include "MLogManip.h" +#include "MHTimeDiffTheta.h" + ClassImp(MHEffOnTime); @@ -68,30 +70,24 @@ 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,11 +96,9 @@ 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,7 +110,6 @@ // 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,4 +120,155 @@ } +// -------------------------------------------------------------------- +// +// 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,211 +276,56 @@ // 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,13 +336,11 @@ 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,10 +374,6 @@ 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,11 +404,7 @@ 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,7 +426,7 @@ } - - - - - +void MHEffOnTime::Calc(const MHTimeDiffTheta &hist, const Bool_t Draw) +{ + Calc(hist.GetHist(), Draw); +} + Index: /trunk/MagicSoft/Mars/mhist/MHEffOnTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffOnTime.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHEffOnTime.h (revision 1668) @@ -10,4 +10,5 @@ class TH2D; +class MHTimeDiffTheta; class MParList; @@ -20,6 +21,11 @@ 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,10 +38,11 @@ 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 }; Index: /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTheta.h (revision 1668) @@ -35,5 +35,5 @@ 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 }; Index: /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h (revision 1668) @@ -35,5 +35,5 @@ 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 }; Index: /trunk/MagicSoft/Mars/mhist/MHFlux.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHFlux.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHFlux.cc (revision 1668) @@ -42,5 +42,4 @@ #include - #include @@ -53,6 +52,69 @@ #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,6 +123,5 @@ // 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,7 +131,4 @@ fVarname = varname; fUnit = unit; - - TString strg(varname); - strg += unit; // char txt[100]; @@ -80,16 +138,16 @@ // 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,14 +157,11 @@ // 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,69 +172,13 @@ // // 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,5 +199,5 @@ Bool_t MHFlux::Fill(const MParContainer *par) { - return kTRUE; + return kTRUE; } @@ -208,46 +207,72 @@ // 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 (barn-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,339 +284,244 @@ // 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,5 +552,5 @@ // ------------------------------------------------------------------------- // -// Draw the histograms +// Draw the histograms // void MHFlux::Draw(Option_t *opt) @@ -642,4 +572,13 @@ 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,5 +599,5 @@ - x2*x1*x1 - x3*x2*x2 - x1*x3*x3; - if (det == 0.0) + if (det==0) { *a = 0; @@ -686,15 +625,4 @@ *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; } Index: /trunk/MagicSoft/Mars/mhist/MHFlux.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHFlux.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHFlux.h (revision 1668) @@ -20,4 +20,8 @@ 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,16 +39,29 @@ // 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,5 +73,9 @@ 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) }; Index: /trunk/MagicSoft/Mars/mhist/MHGamma.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHGamma.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHGamma.cc (revision 1668) @@ -38,8 +38,9 @@ #include -#include #include #include +#include "MHAlphaEnergyTheta.h" + #include "MLog.h" #include "MLogManip.h" @@ -52,6 +53,9 @@ // 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,4 +68,87 @@ 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,88 +157,21 @@ // 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,6 +182,5 @@ // 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,69 +223,72 @@ 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(); +} Index: /trunk/MagicSoft/Mars/mhist/MHGamma.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHGamma.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHGamma.h (revision 1668) @@ -6,30 +6,36 @@ #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 }; Index: /trunk/MagicSoft/Mars/mhist/MHHadronness.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHHadronness.cc (revision 1668) @@ -45,7 +45,4 @@ // - 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,4 +53,9 @@ //////////////////////////////////////////////////////////////////////////// #include "MHHadronness.h" + +/* + // - blue: 2D distance of (acceptance_hadrons, acceptances_gammas) + // to optimum (0, 1): 1-sqrt(Ag*Ag + (1-Ah)*(1-Ah)) + */ #include @@ -113,12 +115,14 @@ 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,5 +140,5 @@ delete fIntPhness; delete fQfac; - delete fMinDist; + // delete fMinDist; delete fGraph; } @@ -266,7 +270,7 @@ 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,26 +364,28 @@ *fLog << "Used " << fGhness->GetEntries() << " Gammas and " << fPhness->GetEntries() << " Hadrons." << endl; - *fLog << " acc(hadron) acc(gamma)" << endl <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,6 +418,8 @@ Getiphness()->SetLineColor(kRed); Getiphness()->DrawCopy("same"); - fMinDist->SetLineColor(kBlue); - fMinDist->DrawCopy("same"); + /* + fMinDist->SetLineColor(kBlue); + fMinDist->DrawCopy("same"); + */ c.cd(3); //GetQfac()->DrawCopy(); @@ -449,4 +457,5 @@ gPad->Update(); } + /* const Int_t h = fMinDist->GetMaximumBin(); TMarker *m = new TMarker(fIntPhness->GetBinContent(h), @@ -455,4 +464,5 @@ m->SetBit(kCanDelete); m->Draw(); + */ return &c; } @@ -482,6 +492,8 @@ Getiphness()->SetLineColor(kRed); Getiphness()->Draw("same"); - fMinDist->SetLineColor(kBlue); - fMinDist->Draw("same"); + /* + fMinDist->SetLineColor(kBlue); + fMinDist->Draw("same"); + */ gPad->cd(3); //GetQfac()->Draw(); @@ -516,4 +528,5 @@ gPad->Update(); } + /* const Int_t h = fMinDist->GetMaximumBin(); TMarker *m = new TMarker(fIntPhness->GetBinContent(h), @@ -522,3 +535,4 @@ m->SetBit(kCanDelete); m->Draw(); -} + */ +} Index: /trunk/MagicSoft/Mars/mhist/MHHadronness.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHHadronness.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHHadronness.h (revision 1668) @@ -15,15 +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: Index: /trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc (revision 1668) @@ -112,5 +112,5 @@ // 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,5 +121,5 @@ // 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); Index: /trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.h (revision 1668) @@ -26,6 +26,6 @@ ~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=""); Index: /trunk/MagicSoft/Mars/mhist/MHMcEnergyMigration.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHMcEnergyMigration.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHMcEnergyMigration.h (revision 1668) @@ -40,5 +40,5 @@ 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 }; Index: /trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHThetabarTheta.cc (revision 1668) @@ -99,5 +99,5 @@ gROOT->SetSelectedPad(NULL); - ((TProfile*)&fHist)->DrawCopy(opt); + ((TProfile)fHist).DrawCopy(opt); c.Modified(); @@ -128,19 +128,8 @@ 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; } - - - - - - - - - - - - - Index: /trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHThetabarTheta.h (revision 1668) @@ -17,8 +17,8 @@ { private: - MTime *fTime; //! - MMcEvt *fMcEvt; //! + MTime *fTime; //! + MMcEvt *fMcEvt; //! - TProfile fHist; + TProfile fHist; public: @@ -36,5 +36,5 @@ TObject *DrawClone(Option_t *option="") const; - ClassDef(MHThetabarTheta, 1) //Profile histogram Thetabar vs. time + ClassDef(MHThetabarTheta, 0) //Profile histogram Thetabar vs. time }; Index: /trunk/MagicSoft/Mars/mhist/MHThetabarTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHThetabarTime.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHThetabarTime.h (revision 1668) @@ -39,5 +39,5 @@ TObject *DrawClone(Option_t *option="") const; - ClassDef(MHThetabarTime, 1) //Profile histogram Thetabar vs. time + ClassDef(MHThetabarTime, 0) //Profile histogram Thetabar vs. time }; Index: /trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.cc (revision 1668) @@ -187,5 +187,4 @@ gPad->Modified(); gPad->Update(); - } @@ -196,7 +195,7 @@ 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; Index: /trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHTimeDiffTheta.h (revision 1668) @@ -19,5 +19,5 @@ MTime *fTime; //! MMcEvt *fMcEvt; //! - Int_t fLastTime; + Double_t fLastTime; TH2D fHist; @@ -37,5 +37,5 @@ TObject *DrawClone(Option_t *option="") const; - ClassDef(MHTimeDiffTheta, 1) //2D-histogram time-diff vs. Theta + ClassDef(MHTimeDiffTheta, 0) //2D-histogram time-diff vs. Theta }; Index: /trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.cc (revision 1668) @@ -100,5 +100,5 @@ { - TCanvas &c = *MakeDefCanvas("DiffTimeTime", "Distrib of \\Delta t, time"); + TCanvas &c = *MakeDefCanvas("DiffTimeTime", "Distrib of dt and t"); c.Divide(2, 2); @@ -186,8 +186,7 @@ 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; Index: /trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.h (revision 1667) +++ /trunk/MagicSoft/Mars/mhist/MHTimeDiffTime.h (revision 1668) @@ -17,5 +17,5 @@ private: MTime *fTime; //! - Int_t fLastTime; + Double_t fLastTime; TH2D fHist; @@ -35,5 +35,5 @@ TObject *DrawClone(Option_t *option="") const; - ClassDef(MHTimeDiffTime, 1) //2D-histogram time-diff vs. time + ClassDef(MHTimeDiffTime, 0) //2D-histogram time-diff vs. time }; Index: /trunk/MagicSoft/Mars/mmain/MAnalysis.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MAnalysis.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MAnalysis.cc (revision 1668) @@ -100,7 +100,7 @@ } -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(); Index: /trunk/MagicSoft/Mars/mmain/MAnalysis.h =================================================================== --- /trunk/MagicSoft/Mars/mmain/MAnalysis.h (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MAnalysis.h (revision 1668) @@ -28,5 +28,5 @@ 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) ; Index: /trunk/MagicSoft/Mars/mmain/MBrowser.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MBrowser.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MBrowser.cc (revision 1668) @@ -275,8 +275,9 @@ } -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,5 +292,4 @@ fList->Add(linesep); AddFrame(linesep, laylinesep); - // Index: /trunk/MagicSoft/Mars/mmain/MBrowser.h =================================================================== --- /trunk/MagicSoft/Mars/mmain/MBrowser.h (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MBrowser.h (revision 1668) @@ -19,5 +19,5 @@ class TGFileContainer; -class MBrowser : public TGTransientFrame +class MBrowser : public TGMainFrame/*TGTransientFrame*/ { private: @@ -62,6 +62,6 @@ 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(); Index: /trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MCameraDisplay.cc (revision 1668) @@ -43,7 +43,7 @@ // '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,5 +91,5 @@ case kButDisplay: new MGCamDisplay(fInputFile, - fClient->GetRoot(), this, 600, 500); + fClient->GetRoot(), /*this,*/ 600, 500); return kTRUE; } Index: /trunk/MagicSoft/Mars/mmain/MCameraDisplay.h =================================================================== --- /trunk/MagicSoft/Mars/mmain/MCameraDisplay.h (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MCameraDisplay.h (revision 1668) @@ -9,5 +9,5 @@ { 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) ; Index: /trunk/MagicSoft/Mars/mmain/MDataCheck.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MDataCheck.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MDataCheck.cc (revision 1668) @@ -54,7 +54,7 @@ // '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); Index: /trunk/MagicSoft/Mars/mmain/MDataCheck.h =================================================================== --- /trunk/MagicSoft/Mars/mmain/MDataCheck.h (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MDataCheck.h (revision 1668) @@ -16,5 +16,5 @@ 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) ; Index: /trunk/MagicSoft/Mars/mmain/MEvtDisp.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MEvtDisp.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MEvtDisp.cc (revision 1668) @@ -39,7 +39,7 @@ }; -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,15 +96,15 @@ 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; } Index: /trunk/MagicSoft/Mars/mmain/MEvtDisp.h =================================================================== --- /trunk/MagicSoft/Mars/mmain/MEvtDisp.h (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MEvtDisp.h (revision 1668) @@ -9,5 +9,5 @@ { 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) ; Index: /trunk/MagicSoft/Mars/mmain/MMars.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MMars.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MMars.cc (revision 1668) @@ -277,21 +277,21 @@ 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; Index: /trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc =================================================================== --- /trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MMonteCarlo.cc (revision 1668) @@ -162,7 +162,7 @@ } -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(); Index: /trunk/MagicSoft/Mars/mmain/MMonteCarlo.h =================================================================== --- /trunk/MagicSoft/Mars/mmain/MMonteCarlo.h (revision 1667) +++ /trunk/MagicSoft/Mars/mmain/MMonteCarlo.h (revision 1668) @@ -28,5 +28,5 @@ 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);